Exploration of Protein Conformational Change with PELE and Meta-Dynamics

Tuning Selectivity in CalA Lipase: Beyond Tunnel Engineering

Engineering studies of Candida (Pseudozyma) antarctica lipase A (CalA) have demonstrated the potential of this enzyme in the selective hydrolysis of fatty acid esters of different chain lengths. CalA has been shown to bind substrates preferentially through an acyl-chain binding tunnel accessed via the hydrolytic active site; it has also been shown that selectivity for substrates of longer or shorter chain length can be tuned, for instance by modulating steric hindrance within the tunnel. Here we demonstrate that, whereas the tunnel region is certainly of paramount importance for substrate recognition, residues in distal regions of the enzyme can also modulate substrate selectivity. To this end, we investigate variants that carry one or more substitutions within the substrate tunnel as well as in distal regions. Combining experimental determination of the substrate selectivity using natural and synthetic substrates with computational characterization of protein dynamics and of tunnels, we deconvolute the effect of key substitutions and demonstrate that epistatic interactions contribute to procuring selectivity toward either long-chain or short/medium-chain fatty acid esters. We demonstrate that various mechanisms contribute to the diverse selectivity profiles, ranging from reshaping tunnel morphology and tunnel stabilization to obstructing the main substrate-binding tunnel, highlighting the dynamic nature of the substrate-binding region. This work provides important insights into the versatility of this robust lipase toward diverse applications.

Cosolvent and Dynamic Effects in Binding Pocket Search by Docking Simulations

The lack of conformational sampling in virtual screening projects can lead to inefficient results because many of the potential drugs may not be able to bind to the target protein during the static docking simulations. Here, we performed ensemble docking for around 2000 United States Food and Drug Administration (FDA)-approved drugs with the RNA-dependent RNA polymerase (RdRp) protein of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) as a target. The representative protein structures were generated by clustering classical molecular dynamics trajectories, which were evolved using three solvent scenarios, namely, pure water, benzene/water and phenol/water mixtures. The introduction of dynamic effects in the theoretical model showed improvement in docking results in terms of the number of strong binders and binding sites in the protein. Some of the discovered pockets were found only for the cosolvent simulations, where the nonpolar probes induced local conformational changes in the protein that lead to the opening of transient pockets. In addition, the selection of the ligands based on a combination of the binding free energy and binding free energy gap between the best two poses for each ligand provided more suitable binders than the selection of ligands based solely on one of the criteria. The application of cosolvent molecular dynamics to enhance the sampling of the configurational space is expected to improve the efficacy of virtual screening campaigns of future drug discovery projects.

aquaPELE: A Monte Carlo-Based Algorithm to Sample the Effects of Buried Water Molecules in Proteins

Water is frequently found inside proteins, carrying out important roles in catalytic reactions or molecular recognition tasks. Therefore, computational models that aim to study protein-ligand interactions usually have to include water effects through explicit or implicit approaches to obtain reliable results. While full explicit models might be too computationally daunting for some applications, implicit models are normally faster but omit some of the most important contributions of water. This is the case of our in-house software, called protein energy landscape exploration (PELE), which uses implicit models to speed up conformational explorations as much as possible; the lack of explicit water sampling, however, limits its model. In this work, we confront this problem with the development of aquaPELE. It is a new algorithm that extends the exploration capabilities while keeping efficiency as it employs a mixed implicit/explicit approach to also take into account the effects of buried water molecules. With an additional Monte Carlo (MC) routine, a set of explicit water molecules is perturbed inside protein cavities and their effects are dynamically adjusted to the current state of the system. As a result, this implementation can be used to predict the principal hydration sites or the rearrangement and displacement of conserved water molecules upon the binding of a ligand. We benchmarked this new tool focusing on estimating ligand binding modes and hydration sites in cavities with important interfacial water molecules, according to crystallographic structures. Results suggest that aquaPELE sets a fast and reliable alternative for molecular recognition studies in systems with a strong water-dependency.

Mapping enzyme-substrate interactions: its potential to study the mechanism of enzymes

With the increase of the need to use more sustainable processes for the industry in our society, the modeling of enzymes has become crucial to fully comprehend their mechanism of action and use this knowledge to enhance and design their properties. A lot of methods to study enzymes computationally exist and they have been classified on sequence-based, structure-based, and the more new artificial intelligence-based ones. Albeit the abundance of methods to help predict the function of an enzyme, molecular modeling is crucial when trying to understand the enzyme mechanism, as they aim to correlate atomistic information with experimental data. Among them, methods that simulate the system dynamics at a molecular mechanics level of theory (classical force fields) have shown to offer a comprehensive study. In this book chapter, we will analyze these techniques, emphasizing the importance of precise modeling of enzyme-substrate interactions. In the end, a brief explanation of the transference of the information from research studies to the industry is given accompanied with two examples of family enzymes where their modeling has helped their exploitation.

Discovery of processive catalysis by an exo-hydrolase with a pocket-shaped active site

Substrates associate and products dissociate from enzyme catalytic sites rapidly, which hampers investigations of their trajectories. The high-resolution structure of the native Hordeum exo-hydrolase HvExoI isolated from seedlings reveals that non-covalently trapped glucose forms a stable enzyme-product complex. Here, we report that the alkyl β-d-glucoside and methyl 6-thio-β-gentiobioside substrate analogues perfused in crystalline HvExoI bind across the catalytic site after they displace glucose, while methyl 2-thio-β-sophoroside attaches nearby. Structural analyses and multi-scale molecular modelling of nanoscale reactant movements in HvExoI reveal that upon productive binding of incoming substrates, the glucose product modifies its binding patterns and evokes the formation of a transient lateral cavity, which serves as a conduit for glucose departure to allow for the next catalytic round. This path enables substrate-product assisted processive catalysis through multiple hydrolytic events without HvExoI losing contact with oligo- or polymeric substrates. We anticipate that such enzyme plasticity could be prevalent among exo-hydrolases.

Enzyme substrates and products often diffuse too rapidly to assess the catalytic implications of these movements. Here, the authors characterise the structural basis of product and substrate diffusion for an exo-hydrolase and discover a substrate-product assisted processive catalytic mechanism.

Computational Design of Protein Function

The computational design of enzymes is a tremendous challenge for both chemistry and biochemistry. The ability to design stable and functional biocatalysts that could operate under different conditions to perform chemical reactions without precedent in nature, allowing the large-scale production of chemicals à la carte, would revolutionise both synthetic, pharmacologic and materials chemistry. Despite the great advances achieved, this highly multidisciplinary area of research is still in its infancy. This chapter describes the ‘inside-out’ protocol for computational enzyme design and both the achievements and limitations of the current technology are highlighted. Furthermore, molecular dynamics simulations have proved to be invaluable in the enzyme design process, constituting an important tool for discovering elusive catalytically relevant conformations of the engineered or designed enzyme. As a complement to the ‘inside-out’ design protocol, different examples where hybrid QM/MM approaches have been directly applied to discover beneficial mutations in rational computational enzyme design are described.

The ins and outs of vanillyl alcohol oxidase: Identification of ligand migration paths

Vanillyl alcohol oxidase (VAO) is a homo-octameric flavoenzyme belonging to the VAO/PCMH family. Each VAO subunit consists of two domains, the FAD-binding and the cap domain. VAO catalyses, among other reactions, the two-step conversion of p-creosol (2-methoxy-4-methylphenol) to vanillin (4-hydroxy-3-methoxybenzaldehyde). To elucidate how different ligands enter and exit the secluded active site, Monte Carlo based simulations have been performed. One entry/exit path via the subunit interface and two additional exit paths have been identified for phenolic ligands, all leading to the si side of FAD. We argue that the entry/exit path is the most probable route for these ligands. A fourth path leading to the re side of FAD has been found for the co-ligands dioxygen and hydrogen peroxide. Based on binding energies and on the behaviour of ligands in these four paths, we propose a sequence of events for ligand and co-ligand migration during catalysis. We have also identified two residues, His466 and Tyr503, which could act as concierges of the active site for phenolic ligands, as well as two other residues, Tyr51 and Tyr408, which could act as a gateway to the re side of FAD for dioxygen. Most of the residues in the four paths are also present in VAO’s closest relatives, eugenol oxidase and p-cresol methylhydroxylase. Key path residues show movements in our simulations that correspond well to conformations observed in crystal structures of these enzymes. Preservation of other path residues can be linked to the electron acceptor specificity and oligomerisation state of the three enzymes. This study is the first comprehensive overview of ligand and co-ligand migration in a member of the VAO/PCMH family, and provides a proof of concept for the use of an unbiased method to sample this process.

Enzymes are bionanomachines, which speed up chemical reactions in organisms. To understand how they achieve that, we need to study their mechanisms. Computational enzymology can show us what happens in the enzyme’s active site during a reaction. But molecules need first to reach the active site before a reaction can start. The process of substrate entry and product exit to the active site is often neglected when studying enzymes. However, these two events are of fundamental importance to the proper functioning of any enzyme. We are interested in these dynamic processes to complete our understanding of the mode of action of enzymes. In our work, we have studied substrate and product migration in vanillyl alcohol oxidase. This enzyme can produce the flavour vanillin and enantiopure alcohols, but also catalyses other reactions. The named products are of interest to the flavour- and fine-chemical industries.

Adaptive simulations, towards interactive protein-ligand modeling

Modeling the dynamic nature of protein-ligand binding with atomistic simulations is one of the main challenges in computational biophysics, with important implications in the drug design process. Although in the past few years hardware and software advances have significantly revamped the use of molecular simulations, we still lack a fast and accurate ab initio description of the binding mechanism in complex systems, available only for up-to-date techniques and requiring several hours or days of heavy computation. Such delay is one of the main limiting factors for a larger penetration of protein dynamics modeling in the pharmaceutical industry. Here we present a game-changing technology, opening up the way for fast reliable simulations of protein dynamics by combining an adaptive reinforcement learning procedure with Monte Carlo sampling in the frame of modern multi-core computational resources. We show remarkable performance in mapping the protein-ligand energy landscape, being able to reproduce the full binding mechanism in less than half an hour, or the active site induced fit in less than 5 minutes. We exemplify our method by studying diverse complex targets, including nuclear hormone receptors and GPCRs, demonstrating the potential of using the new adaptive technique in screening and lead optimization studies.

Simulating Substrate Recognition and Oxidation in Laccases: From Description to Design

To meet the very specific requirements demanded by industry, proteins must be appropriately tailored. Engineering laccases, to improve the oxidation of small molecules, with applications in multiple fields, is, however, a difficult task. Most efforts have concentrated on increasing the redox potential of the enzyme, but in recent work, we have pursued an alternate strategy to engineering these biocatalysts. In particular, we have found that redesigning substrate binding at the T1 pocket, guided by in silico methodologies, to be a more consistent option. In this work, we evaluate the robustness of our computational approach to estimate activity, emphasizing the importance of the binding event in laccase reactivity. Strengths and weaknesses of the protocol are discussed along with its potential for scoring large numbers of protein sequences and thus its significance in protein engineering.

Repurposing designed mutants: a valuable strategy for computer-aided laccase engineering – the case of POXA1b

The broad specificity of laccases, a direct consequence of their shallow binding site, makes this class of enzymes a suitable template to build specificity toward putative substrates.

Enhancing backbone sampling in Monte Carlo simulations using internal coordinates normal mode analysis

Normal mode methods are becoming a popular alternative to sample the conformational landscape of proteins. In this study, we describe the implementation of an internal coordinate normal mode analysis method and its application in exploring protein flexibility by using the Monte Carlo method PELE. This new method alternates two different stages, a perturbation of the backbone through the application of torsional normal modes, and a resampling of the side chains. We have evaluated the new approach using two test systems, ubiquitin and c-Src kinase, and the differences to the original ANM method are assessed by comparing both results to reference molecular dynamics simulations. The results suggest that the sampled phase space in the internal coordinate approach is closer to the molecular dynamics phase space than the one coming from a Cartesian coordinate anisotropic network model. In addition, the new method shows a great speedup (∼5-7×), making it a good candidate for future normal mode implementations in Monte Carlo methods.

The unravelling of the complex pattern of tyrosinase inhibition

Tyrosinases are responsible for melanin formation in all life domains. Tyrosinase inhibitors are used for the prevention of severe skin diseases, in skin-whitening creams and to avoid fruit browning, however continued use of many such inhibitors is considered unsafe. In this study we provide conclusive evidence of the inhibition mechanism of two well studied tyrosinase inhibitors, KA (kojic acid) and HQ (hydroquinone), which are extensively used in hyperpigmentation treatment. KA is reported in the literature with contradicting inhibition mechanisms, while HQ is described as both a tyrosinase inhibitor and a substrate. By visualization of KA and HQ in the active site of TyrBm crystals, together with molecular modeling, binding constant analysis and kinetic experiments, we have elucidated their mechanisms of inhibition, which was ambiguous for both inhibitors. We confirm that while KA acts as a mixed inhibitor, HQ can act both as a TyrBm substrate and as an inhibitor.

Molecular determinants for selective C25-hydroxylation of vitamins D2and D3by fungal peroxygenases

Regioselective hydroxylation of vitamin D byAgrocybe aegeritaandCoprinopsis cinereaperoxygenases was investigated in an experimental and computational study.

New Monte Carlo Based Technique To Study DNA-Ligand Interactions

We present a new all-atom Monte Carlo technique capable of performing quick and accurate DNA-ligand conformational sampling. In particular, and using the PELE software as a frame, we have introduced an additional force field, an implicit solvent, and an anisotropic network model to effectively map the DNA energy landscape. With these additions, we successfully generated DNA conformations for a test set composed of six DNA fragments of A-DNA and B-DNA. Moreover, trajectories generated for cisplatin and its hydrolysis products identified the best interacting compound and binding site, producing analogous results to microsecond molecular dynamics simulations. Furthermore, a combination of the Monte Carlo trajectories with Markov State Models produced noncovalent binding free energies in good agreement with the published molecular dynamics results, at a significantly lower computational cost. Overall our approach will allow a quick but accurate sampling of DNA-ligand interactions.

Unveiling prolyl oligopeptidase ligand migration by comprehensive computational techniques

Prolyl oligopeptidase (POP) is a large 80 kDa protease, which cleaves oligopeptides at the C-terminal side of proline residues and constitutes an important pharmaceutical target. Despite the existence of several crystallographic structures, there is an open debate about migration (entrance and exit) pathways for ligands, and their coupling with protein dynamics. Recent studies have shown the capabilities of molecular dynamics and classical force fields in describing spontaneous binding events and nonbiased ligand migration pathways. Due to POP's size and to the buried nature of its active site, an exhaustive sampling by means of conventional long enough molecular dynamics trajectories is still a nearly impossible task. Such a level of sampling, however, is possible with the breakthrough protein energy landscape exploration technique. Here, we present an exhaustive sampling of POP with a known inhibitor, Z-pro-prolinal. In >3000 trajectories Z-pro-prolinal explores all the accessible surface area, showing multiple entrance events into the large internal cavity through the pore in the β-propeller domain. Moreover, we modeled a natural substrate binding and product release by predicting the entrance of an undecapeptide substrate, followed by manual active site cleavage and nonbiased exit of one of the products (a dipeptide). The product exit shows preference from a flexible 18-amino acid residues loop, pointing to an overall mechanism where entrance and exit occur in different sites.

Study on the interaction of artificial and natural food colorants with human serum albumin: A computational point of view

Due to the high amount of artificial food colorants present in infants' diets, their adverse effects have been of major concern among the literature. Artificial food colorants have been suggested to affect children's behavior, being hyperactivity the most common disorder. In this study we compare binding affinities of a group of artificial colorants (sunset yellow, quinoline yellow, carmoisine, allura red and tartrazine) and their natural industrial equivalents (carminic acid, curcumin, peonidin-3-glucoside, cyanidin-3-glucoside) to human serum albumin (HSA) by a docking approach and further refinement through atomistic molecular dynamics simulations. Due to the protein-ligand conformational interface complexity, we used collective variable driven molecular dynamics to refine docking predictions and to score them according to a hydrogen-bond criterion. With this protocol, we were able to rank ligand affinities to HSA and to compare between the studied natural and artificial food additives. Our results show that the five artificial colorants studied bind better to HSA than their equivalent natural options, in terms of their H-bonding network, supporting the hypothesis of their potential risk to human health.

Insights into Laccase Engineering from Molecular Simulations: Toward a Binding-Focused Strategy

Understanding the molecular determinants of enzyme performance is of primary importance for the rational design of ad hoc mutants. A novel approach, which combines efficient conformational sampling and quick reactivity scoring, is used here to shed light on how substrate oxidation was improved during the directed evolution experiment of a fungal laccase (from Pycnoporus cinnabarinus), an industrially relevant class of oxidoreductases. It is found that the enhanced activity of the evolved enzyme is mainly the result of substrate arrangement in the active site, with no important change in the redox potential of the T1 copper. Mutations at the active site shift the binding mode into a more buried substrate position and provide a more favorable electrostatic environment for substrate oxidation. As a consequence, engineering the binding event seems to be a viable way to in silico evolution of oxidoreductases.

Steroid hydroxylation by basidiomycete peroxygenases: a combined experimental and computational study

The goal of this study is the selective oxyfunctionalization of steroids under mild and environmentally friendly conditions using fungal enzymes. With this purpose, peroxygenases from three basidiomycete species were tested for the hydroxylation of a variety of steroidal compounds, using H2O2 as the only cosubstrate. Two of them are wild-type enzymes from Agrocybe aegerita and Marasmius rotula, and the third one is a recombinant enzyme from Coprinopsis cinerea. The enzymatic reactions on free and esterified sterols, steroid hydrocarbons, and ketones were monitored by gas chromatography, and the products were identified by mass spectrometry. Hydroxylation at the side chain over the steroidal rings was preferred, with the 25-hydroxyderivatives predominating. Interestingly, antiviral and other biological activities of 25-hydroxycholesterol have been reported recently (M. Blanc et al., Immunity 38:106-118, 2013, http://dx.doi.org/10.1016/j.immuni.2012.11.004). However, hydroxylation in the ring moiety and terminal hydroxylation at the side chain also was observed in some steroids, the former favored by the absence of oxygenated groups at C-3 and by the presence of conjugated double bonds in the rings. To understand the yield and selectivity differences between the different steroids, a computational study was performed using Protein Energy Landscape Exploration (PELE) software for dynamic ligand diffusion. These simulations showed that the active-site geometry and hydrophobicity favors the entrance of the steroid side chain, while the entrance of the ring is energetically penalized. Also, a direct correlation between the conversion rate and the side chain entrance ratio could be established that explains the various reaction yields observed.

Differential control of heme reactivity in alpha and beta subunits of hemoglobin: a combined Raman spectroscopic and computational study

The use of hybrid hemoglobin (Hb), with mesoheme substituted for protoheme, allows separate monitoring of the α or β hemes along the allosteric pathway. Using resonance Raman (rR) spectroscopy in silica gel, which greatly slows protein motions, we have observed that the Fe-histidine stretching frequency, νFeHis, which is a monitor of heme reactivity, evolves between frequencies characteristic of the R and T states, for both α or β chains, prior to the quaternary R-T and T-R shifts. Computation of νFeHis, using QM/MM and the conformational search program PELE, produced remarkable agreement with experiment. Analysis of the PELE structures showed that the νFeHis shifts resulted from heme distortion and, in the α chain, Fe-His bond tilting. These results support the tertiary two-state model of ligand binding (Henry et al., Biophys. Chem. 2002, 98, 149). Experimentally, the νFeHis evolution is faster for β than for α chains, and pump-probe rR spectroscopy in solution reveals an inflection in the νFeHis time course at 3 μs for β but not for α hemes, an interval previously shown to be the first step in the R-T transition. In the α chain νFeHis dropped sharply at 20 μs, the final step in the R-T transition. The time courses are fully consistent with recent computational mapping of the R-T transition via conjugate peak refinement by Karplus and co-workers (Fischer et al., Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 5608). The effector molecule IHP was found to lower νFeHis selectively for α chains within the R state, and a binding site in the α1α2 cleft is suggested.

Atomic picture of ligand migration in toluene 4-monooxygenase

Computational modeling combined with mutational and activity assays was used to underline the substrate migration pathways in toluene 4-monooxygenase, a member of the important family of bacterial multicomponent monooxygenases (BMMs). In all structurally defined BMM hydroxylases, several hydrophobic cavities in the α-subunit map a preserved path from the protein surface to the diiron active site. Our results confirm the presence of two pathways by which different aromatic molecules can enter/escape the active site. While the substrate is observed to enter from both channels, the more hydrophilic product is withdrawn mainly from the shorter channel ending at residues D285 and E214. The long channel ends in the vicinity of S395, whose variants have been seen to affect activity and specificity. These mutational effects are clearly reproduced and rationalized by the in silico studies. Furthermore, the combined computational and experimental results highlight the importance of residue F269, which is located at the intersection of the two channels.

Monte Carlo Free Ligand Diffusion with Markov State Model Analysis and Absolute Binding Free Energy Calculations

Obtaining absolute binding free energies from unbiased ligand diffusion has attracted a significant amount of attention due to its implications in drug design. Several studies have used special purpose computers and software to achieve microsecond molecular dynamics which, combined with a Markov state model analysis, are capable of providing absolute binding free energies. We have recently developed a Monte Carlo based technique, PELE, capable of performing a dynamical exploration of the protein-ligand energy landscape including free ligand diffusion into the active site, at a fraction of the computational cost of molecular dynamics techniques. We demonstrate here the capabilities of our Monte Carlo technique in obtaining absolute binding free energies for a series of benzamidine like inhibitors into trypsin. Our results are in good agreement with experimental data and other molecular dynamics simulations, indicating that PELE can be a useful tool for quick estimates of binding free energies and mechanisms.

PELE web server: atomistic study of biomolecular systems at your fingertips

PELE, Protein Energy Landscape Exploration, our novel technology based on protein structure prediction algorithms and a Monte Carlo sampling, is capable of modelling the all-atom protein–ligand dynamical interactions in an efficient and fast manner, with two orders of magnitude reduced computational cost when compared with traditional molecular dynamics techniques. PELE’s heuristic approach generates trial moves based on protein and ligand perturbations followed by side chain sampling and global/local minimization. The collection of accepted steps forms a stochastic trajectory. Furthermore, several processors may be run in parallel towards a collective goal or defining several independent trajectories; the whole procedure has been parallelized using the Message Passing Interface. Here, we introduce the PELE web server, designed to make the whole process of running simulations easier and more practical by minimizing input file demand, providing user-friendly interface and producing abstract outputs (e.g. interactive graphs and tables). The web server has been implemented in C++ using Wt (http://www.webtoolkit.eu) and MySQL (http://www.mysql.com). The PELE web server, accessible at http://pele.bsc.es, is free and open to all users with no login requirement.

Molecular interactions of prodiginines with the BH3 domain of anti-apoptotic Bcl-2 family members

Prodigiosin and obatoclax, members of the prodiginines family, are small molecules with anti-cancer properties that are currently under preclinical and clinical trials. The molecular target(s) of these agents, however, is an open question. Combining experimental and computational techniques we find that prodigiosin binds to the BH3 domain in some BCL-2 protein families, which play an important role in the apoptotic programmed cell death. In particular, our results indicate a large affinity of prodigiosin for MCL-1, an anti-apoptotic member of the BCL-2 family. In melanoma cells, we demonstrate that prodigiosin activates the mitochondrial apoptotic pathway by disrupting MCL-1/BAK complexes. Computer simulations with the PELE software allow the description of the induced fit process, obtaining a detailed atomic view of the molecular interactions. These results provide new data to understand the mechanism of action of these molecules, and assist in the development of more specific inhibitors of anti-apoptotic BCL-2 proteins.

Single vs. multiple ligand pathways in globins: a computational view

Diatomic ligand migration in globins has been the subject of numerous studies. Still, a consensus picture for the ligand entrance is not clear, with a growing concern among experimental researchers that computational simulations always show multiple pathways for any globin. Modeling non-biased ligand entrance from conventional molecular dynamics techniques, however, has shown to be difficult (and expensive). Here we use our Monte Carlo methodology, capable of freely mapping ligand diffusion and the description of rare events, to two well-studied systems: myoglobin and the mini-hemoglobin from the sea worm Cerebratulus lacteus. Our results clearly show that the simulations are specific to the system providing a different trend in the entrance pathway, as expected from experiments. While Mb presents multiple entrance pathways, populating the well-known xenon cavities, in CerHb the ligand enters the protein only by one apolar channel. Most of the trajectories (64%) visiting myoglobin's active site though, are gated by the distal histidine. Such detailed information, accessible through the state of the art algorithms in PELE, is computationally inexpensive and available to all non-profit researchers. This article is part of a Special Issue entitled: Oxygen Binding and Sensing Proteins.