Hidden thermodynamics of mutant proteins: a molecular dynamics analysis

The physics-AI dialogue in drug design

A long path has led from the determination of the first protein structure in 1960 to the recent breakthroughs in protein science. Protein structure prediction and design methodologies based on machine learning (ML) have been recognized with the 2024 Nobel prize in Chemistry, but they would not have been possible without previous work and the input of many domain scientists. Challenges remain in the application of ML tools for the prediction of structural ensembles and their usage within the software pipelines for structure determination by crystallography or cryogenic electron microscopy. In the drug discovery workflow, ML techniques are being used in diverse areas such as scoring of docked poses, or the generation of molecular descriptors. As the ML techniques become more widespread, novel applications emerge which can profit from the large amounts of data available. Nevertheless, it is essential to balance the potential advantages against the environmental costs of ML deployment to decide if and when it is best to apply it. For hit to lead optimization ML tools can efficiently interpolate between compounds in large chemical series but free energy calculations by molecular dynamics simulations seem to be superior for designing novel derivatives. Importantly, the potential complementarity and/or synergism of physics-based methods (e.g., force field-based simulation models) and data-hungry ML techniques is growing strongly. Current ML methods have evolved from decades of research. It is now necessary for biologists, physicists, and computer scientists to fully understand advantages and limitations of ML techniques to ensure that the complementarity of physics-based methods and ML tools can be fully exploited for drug design.

The structural and mechanistic bases for the viral resistance to allosteric HIV-1 integrase inhibitor pirmitegravir

Allosteric HIV-1 integrase (IN) inhibitors (ALLINIs) are investigational antiretroviral agents that potently impair virion maturation by inducing hyper-multimerization of IN and inhibiting its interaction with viral genomic RNA. The pyrrolopyridine-based ALLINI pirmitegravir (PIR) has recently advanced into phase 2a clinical trials. Previous cell culture-based viral breakthrough assays identified the HIV-1(Y99H/A128T IN) variant that confers substantial resistance to this inhibitor. Here, we have elucidated the unexpected mechanism of viral resistance to PIR. Although both Tyr99 and Ala128 are positioned within the inhibitor binding V-shaped cavity at the IN catalytic core domain (CCD) dimer interface, the Y99H/A128T IN mutations did not substantially affect the direct binding of PIR to the CCD dimer or functional oligomerization of full-length IN. Instead, the drug-resistant mutations introduced a steric hindrance at the inhibitor-mediated interface between CCD and C-terminal domain (CTD) and compromised CTD binding to the CCDY99H/A128T + PIR complex. Consequently, full-length INY99H/A128T was substantially less susceptible to the PIR-induced hyper-multimerization than the WT protein, and HIV-1(Y99H/A128T IN) conferred >150-fold resistance to the inhibitor compared with the WT virus. By rationally modifying PIR, we have developed its analog EKC110, which readily induced hyper-multimerization of INY99H/A128T in vitro and was ~14-fold more potent against HIV-1(Y99H/A128T IN) than the parent inhibitor. These findings suggest a path for developing improved PIR chemotypes with a higher barrier to resistance for their potential clinical use.IMPORTANCEAntiretroviral therapies save the lives of millions of people living with HIV (PLWH). However, the evolution of multi-drug-resistant viral phenotypes is a major clinical problem, and there are limited or no treatment options for heavily treatment-experienced PLWH. Allosteric HIV-1 integrase inhibitors (ALLINIs) are a novel class of antiretroviral compounds that work by a unique mechanism of binding to the non-catalytic site on the viral protein and inducing aberrant integrase multimerization. Accordingly, ALLINIs potently inhibit both wild-type HIV-1 and all drug-resistant viral phenotypes that have so far emerged against currently used therapies. Pirmitegravir, a highly potent and safe investigational ALLINI, is currently advancing through clinical trials. Here, we have elucidated the structural and mechanistic bases behind the emergence of HIV-1 integrase mutations in infected cells that confer resistance to pirmitegravir. In turn, our findings allowed us to rationally develop an improved ALLINI with substantially enhanced potency against the pirmitegravir-resistant virus.

Alchemical Transformations and Beyond: Recent Advances and Real-World Applications of Free Energy Calculations in Drug Discovery

Computational methods constitute efficient strategies for screening and optimizing potential drug molecules. A critical factor in this process is the binding affinity between candidate molecules and targets, quantified as binding free energy. Among various estimation methods, alchemical transformation methods stand out for their theoretical rigor. Despite challenges in force field accuracy and sampling efficiency, advancements in algorithms, software, and hardware have increased the application of free energy perturbation (FEP) calculations in the pharmaceutical industry. Here, we review the practical applications of FEP in drug discovery projects since 2018, covering both ligand-centric and residue-centric transformations. We show that relative binding free energy calculations have steadily achieved chemical accuracy in real-world applications. In addition, we discuss alternative physics-based simulation methods and the incorporation of deep learning into free energy calculations.

GENESIS 2.1: High-Performance Molecular Dynamics Software for Enhanced Sampling and Free-Energy Calculations for Atomistic, Coarse-Grained, and Quantum Mechanics/Molecular Mechanics Models

GENeralized-Ensemble SImulation System (GENESIS) is a molecular dynamics (MD) software developed to simulate the conformational dynamics of a single biomolecule, as well as molecular interactions in large biomolecular assemblies and between multiple biomolecules in cellular environments. To achieve the latter purpose, the earlier versions of GENESIS emphasized high performance in atomistic MD simulations on massively parallel supercomputers, with or without graphics processing units (GPUs). Here, we implemented multiscale MD simulations that include atomistic, coarse-grained, and hybrid quantum mechanics/molecular mechanics (QM/MM) calculations. They demonstrate high performance and are integrated with enhanced conformational sampling algorithms and free-energy calculations without using external programs except for the QM programs. In this article, we review new functions, molecular models, and other essential features in GENESIS version 2.1 and discuss ongoing developments for future releases.

Binding characteristics of the doxepin E/Z-isomers to the histamine H1 receptor revealed by receptor-bound ligand analysis and molecular dynamics study

Doxepin is an antihistamine and tricyclic antidepressant that binds to the histamine H1 receptor (H1R) with high affinity. Doxepin is an 85:15 mixture of the E- and Z-isomers. The Z-isomer is well known to be more effective than the E-isomer, whereas based on the crystal structure of the H1R/doxepin complex, the hydroxyl group of Thr1123.37 is close enough to form a hydrogen bond with the oxygen atom of the E-isomer. The detailed binding characteristics and reasons for the differences remain unclear. In this study, we analyzed doxepin isomers bound to the receptor following extraction from a purified H1R protein complexed with doxepin. The ratio of the E- and Z-isomers bound to wild-type (WT) H1R was 55:45, indicating that the Z-isomer was bound to WT H1R with an approximately 5.2-fold higher affinity than the E-isomer. For the T1123.37V mutant, the E/Z ratio was 89:11, indicating that both isomers have similar affinities. Free energy calculations using molecular dynamics (MD) simulations also reproduced the experimental results of the relative binding free energy differences between the isomers for WT and T1123.37V. Furthermore, MD simulations revealed that the hydroxyl group of T1123.37 did not form hydrogen bonds with the E-isomer, but with the adjacent residues in the binding pocket. Analysis of the receptor-bound doxepin and MD simulations suggested that the hydroxyl group of T1123.37 contributes to the formation of a chemical environment in the binding pocket, which is slightly more favorable for the Z-isomer without hydrogen bonding with doxepin.

Kartograf: A Geometrically Accurate Atom Mapper for Hybrid-Topology Relative Free Energy Calculations

Relative binding free energy (RBFE) calculations have emerged as a powerful tool that supports ligand optimization in drug discovery. Despite many successes, the use of RBFEs can often be limited by automation problems, in particular, the setup of such calculations. Atom mapping algorithms are an essential component in setting up automatic large-scale hybrid-topology RBFE calculation campaigns. Traditional algorithms typically employ a 2D subgraph isomorphism solver (SIS) in order to estimate the maximum common substructure. SIS-based approaches can be limited by time-intensive operations and issues with capturing geometry-linked chemical properties, potentially leading to suboptimal solutions. To overcome these limitations, we have developed Kartograf, a geometric-graph-based algorithm that uses primarily the 3D coordinates of atoms to find a mapping between two ligands. In free energy approaches, the ligand conformations are usually derived from docking or other previous modeling approaches, giving the coordinates a certain importance. By considering the spatial relationships between atoms related to the molecule coordinates, our algorithm bypasses the computationally complex subgraph matching of SIS-based approaches and reduces the problem to a much simpler bipartite graph matching problem. Moreover, Kartograf effectively circumvents typical mapping issues induced by molecule symmetry and stereoisomerism, making it a more robust approach for atom mapping from a geometric perspective. To validate our method, we calculated mappings with our novel approach using a diverse set of small molecules and used the mappings in relative hydration and binding free energy calculations. The comparison with two SIS-based algorithms showed that Kartograf offers a fast alternative approach. The code for Kartograf is freely available on GitHub (https://github.com/OpenFreeEnergy/kartograf). While developed for the OpenFE ecosystem, Kartograf can also be utilized as a standalone Python package.

The structural and mechanistic bases for the viral resistance to allosteric HIV-1 integrase inhibitor pirmitegravir

Allosteric HIV-1 integrase (IN) inhibitors (ALLINIs) are investigational antiretroviral agents which potently impair virion maturation by inducing hyper-multimerization of IN and inhibiting its interaction with viral genomic RNA. The pyrrolopyridine-based ALLINI pirmitegravir (PIR) has recently advanced into Phase 2a clinical trials. Previous cell culture based viral breakthrough assays identified the HIV-1(Y99H/A128T IN) variant that confers substantial resistance to this inhibitor. Here, we have elucidated the unexpected mechanism of viral resistance to PIR. While both Tyr99 and Ala128 are positioned within the inhibitor binding V-shaped cavity at the IN catalytic core domain (CCD) dimer interface, the Y99H/A128T IN mutations did not substantially affect direct binding of PIR to the CCD dimer or functional oligomerization of full-length IN. Instead, the drug-resistant mutations introduced a steric hindrance at the inhibitor mediated interface between CCD and C-terminal domain (CTD) and compromised CTD binding to the CCDY99H/A128T + PIR complex. Consequently, full-length INY99H/A128T was substantially less susceptible to the PIR induced hyper-multimerization than the WT protein, and HIV-1(Y99H/A128T IN) conferred >150-fold resistance to the inhibitor compared to the WT virus. By rationally modifying PIR we have developed its analog EKC110, which readily induced hyper-multimerization of INY99H/A128T in vitro and was ~14-fold more potent against HIV-1(Y99H/A128T IN) than the parent inhibitor. These findings suggest a path for developing improved PIR chemotypes with a higher barrier to resistance for their potential clinical use.

Accurate and Efficient Multilevel Free Energy Simulations with Neural Network-Assisted Enhanced Sampling

Free energy differences (ΔF) are essential to quantitative characterization and understanding of chemical and biological processes. Their direct estimation with an accurate quantum mechanical potential is of great interest and yet impractical due to high computational cost and incompatibility with typical alchemical free energy protocols. One promising solution is the multilevel free energy simulation in which the estimate of ΔF at an inexpensive low level of theory is combined with the correction toward a higher level of theory. The poor configurational overlap generally expected between the two levels of theory, however, presents a major challenge. We overcome this challenge by using a deep neural network model and enhanced sampling simulations. An adversarial autoencoder is used to identify a low-dimensional (latent) space that compactly represents the degrees of freedom that encode the distinct distributions at the two levels of theory. Enhanced sampling in this latent space is then used to drive the sampling of configurations that predominantly contribute to the free energy correction. Results for both gas phase and condensed phase systems demonstrate that this data-driven approach offers high accuracy and efficiency with great potential for scalability to complex systems.

Can Free Energy Perturbation Simulations Coupled with Replica-Exchange Molecular Dynamics Study Ligands with Distributed Binding Sites?

Free energy perturbation coupled with replica exchange with solute tempering (FEP/REST) offers a rigorous approach to compute relative free energy changes for ligands. To determine the applicability of FEP/REST for the ligands with distributed binding poses, we considered two alchemical transformations involving three putative inhibitors I0, I1, and I2 of the Venezuelan equine encephalitis virus nuclear localization signal sequence binding to the importin-α (impα) transporter protein. I0 → I1 and I0 → I2 transformations, respectively, increase or decrease the polarity of the parent molecule. Our objective was three-fold─(i) to verify FEP/REST technical performance and convergence, (ii) to estimate changes in binding free energy ΔΔG, and (iii) to determine the utility of FEP/REST simulations for conformational binding analysis. Our results are as follows. First, our FEP/REST implementation properly follows FEP/REST formalism and produces converged ΔΔG estimates. Due to ligand inherent unbinding, the better FEP/REST strategy lies in performing multiple independent trajectories rather than extending their length. Second, I0 → I1 and I0 → I2 transformations result in overall minor changes in inhibitor binding free energy, slightly strengthening the affinity of I1 and weakening that of I2. Electrostatic interactions dominate binding interactions, determining the enthalpic changes. The two transformations cause opposite entropic changes, which ultimately govern binding affinities. Importantly, we confirm the validity of FEP/REST free energy estimates by comparing them with our previous REST simulations, directly probing binding of three ligands to impα. Third, we established that FEP/REST simulations can sample binding ensembles of ligands. Thus, FEP/REST can be applied (i) to study the energetics of the ligand binding without defined poses and showing minor differences in affinities |ΔΔG| ≲ 0.5 kcal/mol and (ii) to collect ligand binding conformational ensembles.

Exploring Routes to Enhance the Calculation of Free Energy Differences via Non-Equilibrium Work SQM/MM Switching Simulations Using Hybrid Charge Intermediates between MM and SQM Levels of Theory or Non-Linear Switching Schemes

Non-equilibrium work switching simulations and Jarzynski's equation are a reliable method for computing free energy differences, ΔAlow→high, between two levels of theory, such as a pure molecular mechanical (MM) and a quantum mechanical/molecular mechanical (QM/MM) description of a system of interest. Despite the inherent parallelism, the computational cost of this approach can quickly become very high. This is particularly true for systems where the core region, the part of the system to be described at different levels of theory, is embedded in an environment such as explicit solvent water. We find that even for relatively simple solute-water systems, switching lengths of at least 5 ps are necessary to compute ΔAlow→high reliably. In this study, we investigate two approaches towards an affordable protocol, with an emphasis on keeping the switching length well below 5 ps. Inserting a hybrid charge intermediate state with modified partial charges, which resembles the charge distribution of the desired high level, makes it possible to obtain reliable calculations with 2 ps switches. Attempts using step-wise linear switching paths, on the other hand, did not lead to improvement, i.e., a faster convergence for all systems. To understand these findings, we analyzed the solutes' properties as a function of the partial charges used and the number of water molecules in direct contact with the solute, and studied the time needed for water molecules to reorient themselves upon a change in the solute's charge distribution.

The Drug-Induced Interface That Drives HIV-1 Integrase Hypermultimerization and Loss of Function

Allosteric HIV-1 integrase (IN) inhibitors (ALLINIs) are an emerging class of small molecules that disrupt viral maturation by inducing the aberrant multimerization of IN. Here, we present cocrystal structures of HIV-1 IN with two potent ALLINIs, namely, BI-D and the drug candidate Pirmitegravir. The structures reveal atomistic details of the ALLINI-induced interface between the HIV-1 IN catalytic core and carboxyl-terminal domains (CCD and CTD). Projecting from their principal binding pocket on the IN CCD dimer, the compounds act as molecular glue by engaging a triad of invariant HIV-1 IN CTD residues, namely, Tyr226, Trp235, and Lys266, to nucleate the CTD-CCD interaction. The drug-induced interface involves the CTD SH3-like fold and extends to the beginning of the IN carboxyl-terminal tail region. We show that mutations of HIV-1 IN CTD residues that participate in the interface with the CCD greatly reduce the IN-aggregation properties of Pirmitegravir. Our results explain the mechanism of the ALLINI-induced condensation of HIV-1 IN and provide a reliable template for the rational development of this series of antiretrovirals through the optimization of their key contacts with the viral target. IMPORTANCE Despite the remarkable success of combination antiretroviral therapy, HIV-1 remains among the major causes of human suffering and loss of life in poor and developing nations. To prevail in this drawn-out battle with the pandemic, it is essential to continue developing advanced antiviral agents to fight drug resistant HIV-1 variants. Allosteric integrase inhibitors (ALLINIs) are an emerging class of HIV-1 antagonists that are orthogonal to the current antiretroviral drugs. These small molecules act as highly specific molecular glue, which triggers the aggregation of HIV-1 integrase. In this work, we present high-resolution crystal structures that reveal the crucial interactions made by two potent ALLINIs, namely, BI-D and Pirmitegravir, with HIV-1 integrase. Our results explain the mechanism of drug action and will inform the development of this promising class of small molecules for future use in antiretroviral regimens.

Leveraging the Sampling Efficiency of RE-EDS in OpenMM Using a Shifted Reaction-Field With an Atom-Based Cutoff

Replica-exchange enveloping distribution sampling (RE-EDS) is a pathway-independent multistate free-energy method, currently implemented in the GROMOS software package for molecular dynamics (MD) simulations. It has a high intrinsic sampling efficiency as the interactions between the unperturbed particles have to be calculated only once for multiple end-states. As a result, RE-EDS is an attractive method for the calculation of relative solvation and binding free energies. An essential requirement for reaching this high efficiency is the separability of the nonbonded interactions into solute-solute, solute-environment, and environment-environment contributions. Such a partitioning is trivial when using a Coulomb term with a reaction-field (RF) correction to model the electrostatic interactions, but not when using lattice- sum schemes. To avoid cutoff artifacts, the RF correction is typically used in combination with a charge-group based cutoff, which is not supported by most small-molecule force fields and other MD engines. To address this issue, we investigate the combination of RE-EDS simulations with a recently introduced RF scheme including a shifting function that enables the rigorous calculation of RF electrostatics with atom-based cutoffs. The resulting approach is validated by calculating solvation free energies with the generalized AMBER force field (GAFF) in water and chloroform using both the GROMOS software package and a proof-of-concept implementation in OpenMM.

Investigation of the Complexation between 4-Aminoazobenzene and Cucurbit[7]uril through a Combined Spectroscopic, Nuclear Magnetic Resonance, and Molecular Simulation Studies

Cucurbiturils are well known for their ability to form supramolecular systems with ultrahigh affinities binding. Inclusion complex between 4-aminoazobenzene and cucurbit[7]uril has been investigated in aqueous solution by ultraviolet (UV)-spectroscopy, ¹H NMR, and molecular simulations. 4-aminoazobenzene shows high affinity in acidic solutions while no association was detected in neutral solutions. The thermodynamic properties of complex formation are investigated using both UV spectroscopy and nuclear magnetic resonance (NMR) measurements. Our results highlight that the high binding constant between CB7 and 4AA (log K = 4.9) is the result of a large negative change in Δ_r H° (-19 kJ/mol) and a small positive change in TΔ_r S° (9 kJ/mol). The analysis of the experimental data lead to hypothesis on the structure of the complex. We have used molecular dynamics simulation to interpret experiments. Interestingly, the cis-trans isomerization of aminoazobenzene is considered. All the results are discussed and compared with those previously obtained with other host molecules.

Alchemical Free Energy Calculations to Investigate Protein-Protein Interactions: the Case of the CDC42/PAK1 Complex

Here, we show that alchemical free energy calculations can quantitatively compute the effect of mutations at the protein-protein interface. As a test case, we have used the protein complex formed by the small Rho-GTPase CDC42 and its downstream effector PAK1, a serine/threonine kinase. Notably, the CDC42/PAK1 complex offers a wealth of structural, mutagenesis, and binding affinity data because of its central role in cellular signaling and cancer progression. In this context, we have considered 16 mutations in the CDC42/PAK1 complex and obtained excellent agreement between computed and experimental data on binding affinity. Importantly, we also show that a careful analysis of the side-chain conformations in the mutated amino acids can considerably improve the computed estimates, solving issues related to sampling limitations. Overall, this study demonstrates that alchemical free energy calculations can conveniently be integrated into the design of experimental mutagenesis studies.

Collective Variables for Free Energy Surface Tailoring: Understanding and Modifying Functionality in Systems Dominated by Rare Events

We introduce a method for elucidating and modifying the functionality of systems dominated by rare events that relies on the semiautomated tuning of their underlying free energy surface. The proposed approach seeks to construct collective variables (CVs) that encode the essential information regarding the rare events of the system of interest. The appropriate CVs are identified using harmonic linear discriminant analysis (HLDA), a machine-learning-based method that is trained solely on data collected from short ordinary simulations in the relevant metastable states of the system. Utilizing the interpretable form of the resulting CVs, the critical interaction potentials that determine the system's rare transitions are identified and purposely modified to tailor the free energy surface in a manner that alters functionality as desired. The applicability of the method is illustrated in the context of three different systems, thereby demonstrating that thermodynamic and kinetic properties can be tractably modified with little to no prior knowledge or intuition.

Electronic Polarization Is Essential for the Stabilization and Dynamics of Buried Ion Pairs in Staphylococcal Nuclease Mutants

Buried charged residues play important roles in the modulation of protein stabilities and conformational dynamics and make crucial contributions to protein functions. Considering the generally nonpolar nature of protein interior, a key question concerns the contribution of electronic polarization to the stabilization and properties of buried charges. We answer this question by conducting free energy simulations using the latest polarizable CHARMM force field based on Drude oscillators for a series of Staphylococcal nuclease mutants that involve a buried Glu-Lys pair in different titration states and orientations. While a nonpolarizable model suggests that the ionized form of the buried Glu-Lys pair is more than 40 kcal/mol less stable than the charge-neutral form, the two titration states are comparable in stability when electronic polarization is included explicitly, a result better reconcilable with available experimental data. Analysis of free energy components suggests that additional stabilization of the ionized Glu-Lys pair has contributions from both the enhanced salt-bridge strength and stronger interaction between the ion-pair and surrounding protein residues and penetrated water. Despite the stronger direct interaction between Glu and Lys, the ion-pair exhibits considerably larger and faster structural fluctuations when polarization is included, due to compensation of interactions in the cavity. Collectively, observations from this work provide compelling evidence that electronic polarization is essential to the stability, hydration, dynamics, and therefore function of buried charges in proteins. Therefore, our study advocates for the explicit consideration of electronic polarization for mechanistic and engineering studies that implicate buried charged residues, such as enzymes and ion transporters.

RestraintMaker: a graph-based approach to select distance restraints in free-energy calculations with dual topology

The calculation of relative binding free energies (RBFE) involves the choice of the end-state/system representation, of a sampling approach, and of a free-energy estimator. System representations are usually termed "single topology" or "dual topology". As the terminology is often used ambiguously in the literature, a systematic categorization of the system representations is proposed here. In the dual-topology approach, the molecules are simulated as separate molecules. Such an approach is relatively easy to automate for high-throughput RBFE calculations compared to the single-topology approach. Distance restraints are commonly applied to prevent the molecules from drifting apart, thereby improving the sampling efficiency. In this study, we introduce the program RestraintMaker, which relies on a greedy algorithm to find (locally) optimal distance restraints between pairs of atoms based on geometric measures. The algorithm is further extended for multi-state methods such as enveloping distribution sampling (EDS) or multi-site [Formula: see text]-dynamics. The performance of RestraintMaker is demonstrated for toy models and for the calculation of relative hydration free energies. The Python program can be used in script form or through an interactive GUI within PyMol. The selected distance restraints can be written out in GROMOS or GROMACS file formats. Additionally, the program provides a human-readable JSON format that can easily be parsed and processed further. The code of RestraintMaker is freely available on GitHub https://github.com/rinikerlab/restraintmaker.

Virtual Double-System Single-Box for Absolute Dissociation Free Energy Calculations in GROMACS

We describe a step-by-step protocol for the computation of absolute dissociation free energy with GROMACS code and PLUMED library, which exploits a combination of advanced sampling techniques and nonequilibrium alchemical methodologies. The computational protocol has been automated through an open source Python middleware (HPC_Drug) which allows one to set up the GROMACS/PLUMED input files for execution on high performing computing facilities. The proposed protocol, by exploiting its inherent parallelism and the power of the GROMACS code on graphical processing units, has the potential to afford accurate and precise estimates of the dissociation constants in drug-receptor systems described at the atomistic level. The procedure has been applied to the calculation of the absolute dissociation free energy of PF-07321332, an oral antiviral proposed by Pfizer, with the main protease (3CL^pro) of SARS-CoV-2.

Reverse Protonation of Buried Ion-Pairs in Staphylococcal Nuclease Mutants

Although buried titratable residues in protein cavities are often of major functional importance, it is generally challenging to understand their properties such as the ionization state and factors of stabilization based on experimental studies alone. A specific set of examples involve buried Glu-Lys pairs in a series of variants of Staphylococcal nuclease, for which recent structural and thermodynamic studies appeared to suggest that both the stability and the ionization state of the buried Glu-Lys pair are sensitive to its orientation (i.e., Glu23-Lys36 vs Lys23-Glu36). To further clarify the situation, especially ionization states of the buried Glu-Lys pairs, we have conducted extensive molecular dynamics simulations and free energy computations. Microsecond molecular dynamics simulations show that the hydration level of the cavity depends on the orientation of the buried ion-pair therein as well as its ionization state; free energy simulations recapitulate the relative stability of Glu23-Lys36 (EK) vs Lys23-Glu36 (KE) mutants measured experimentally, although the difference is similar in magnitude regardless of the ionization state of the Glu-Lys pair. A complementary set of free energy simulations strongly suggests that, in contrast to the original suggestion in the experimental analysis, the Glu and Lys residues prefer to adopt their charge-neutral rather than the ionized states. This result is consistent with the low dielectric constant computed for water in the cavity, which makes it difficult for the protein cavity to stabilize a pair of charged Glu-Lys residues, even with water penetration. The current study highlights the role of free energy simulations in understanding the ionization state of buried titratable residues and the relevant energetic contributions, forming the basis for the rational design of buried charge networks in proteins.

Free-Energy Landscape of a pH-Modulated G·C Base Pair Transition from Watson-Crick to Hoogsteen State in Duplex DNA

In double-helical DNAs, the most stable Watson-Crick (WC) base pair (bp) can be in thermal equilibrium with much less abundant Hoogsteen (HG) bp by the spontaneous rotation of the glycosidic angle in purine bases. Previous experimental studies showed that in the case of a G·C bp, the population of the transient HG is enhanced as a protonated form (HG⁺) through the protonation of the cytosine base under weakly acidic conditions. Hence, pH is a key factor that can modulate this transition event from the WC to HG⁺ bp. In this study, to computationally probe the overall free-energy landscapes of this pH-modulated G·C HG breathing, a comprehensive classical molecular dynamics (MD) simulation protocol is proposed using an enhanced sampling MD in conjunction with the standard thermodynamic integration method. From this MD protocol proposed, the free-energy surfaces of the G·C bp transition from the WC to HG bp were constructed successfully at any pH range, producing pH-dependent free-energy quantities in close agreement with previously reported experimental results. The simulation protocol is expected to provide valuable atomistic insight into the DNA bp transition events coupled with protonation or tautomeric shift in a target bp.

Implementing and Assessing an Alchemical Method for Calculating Protein-Protein Binding Free Energy

Protein-protein binding is fundamental to most biological processes. It is important to be able to use computation to accurately estimate the change in protein-protein binding free energy due to mutations in order to answer biological questions that would be experimentally challenging, laborious, or time-consuming. Although nonrigorous free-energy methods are faster, rigorous alchemical molecular dynamics-based methods are considerably more accurate and are becoming more feasible with the advancement of computer hardware and molecular simulation software. Even with sufficient computational resources, there are still major challenges to using alchemical free-energy methods for protein-protein complexes, such as generating hybrid structures and topologies, maintaining a neutral net charge of the system when there is a charge-changing mutation, and setting up the simulation. In the current study, we have used the pmx package to generate hybrid structures and topologies, and a double-system/single-box approach to maintain the net charge of the system. To test the approach, we predicted relative binding affinities for two protein-protein complexes using a nonequilibrium alchemical method based on the Crooks fluctuation theorem and compared the results with experimental values. The method correctly identified stabilizing from destabilizing mutations for a small protein-protein complex, and a larger, more challenging antibody complex. Strong correlations were obtained between predicted and experimental relative binding affinities for both protein-protein systems.

Inhibitor binding influences the protonation states of histidines in SARS-CoV-2 main protease

The main protease (Mpro) of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is an attractive target for antiviral therapeutics. Recently, many high-resolution apo and inhibitor-bound structures of Mpro, a cysteine protease, have been determined, facilitating structure-based drug design. Mpro plays a central role in the viral life cycle by catalyzing the cleavage of SARS-CoV-2 polyproteins. In addition to the catalytic dyad His41-Cys145, Mpro contains multiple histidines including His163, His164, and His172. The protonation states of these histidines and the catalytic nucleophile Cys145 have been debated in previous studies of SARS-CoV Mpro, but have yet to be investigated for SARS-CoV-2. In this work we have used molecular dynamics simulations to determine the structural stability of SARS-CoV-2 Mpro as a function of the protonation assignments for these residues. We simulated both the apo and inhibitor-bound enzyme and found that the conformational stability of the binding site, bound inhibitors, and the hydrogen bond networks of Mpro are highly sensitive to these assignments. Additionally, the two inhibitors studied, the peptidomimetic N3 and an α-ketoamide, display distinct His41/His164 protonation-state-dependent stabilities. While the apo and the N3-bound systems favored N δ (HD) and N ϵ (HE) protonation of His41 and His164, respectively, the α-ketoamide was not stably bound in this state. Our results illustrate the importance of using appropriate histidine protonation states to accurately model the structure and dynamics of SARS-CoV-2 Mpro in both the apo and inhibitor-bound states, a necessary prerequisite for drug-design efforts.

Redox Potentials of Disulfide Bonds in LOXL2 Studied by Nonequilibrium Alchemical Simulation

Lysyl oxidase-like 2 (LOXL2) is a metalloenzyme that catalyzes the oxidative deamination ε-amino group of lysine. It is found that LOXL2 is a promotor for the metastasis and invasion of cancer cells. Disulfide bonds are important components in LOXL2, and they play a stabilizing role for protein structure or a functional role for regulating protein bioactivity. The redox potential of disulfide bond is one important property to determine the functional role of disulfide bond. In this study, we have calculated the reduction potential of all the disulfide bonds in LOXL2 by non-equilibrium alchemical simulations. Our results show that seven of seventeen disulfide bonds have high redox potentials between -182 and -298 mV and could have a functional role, viz., Cys573-Cys625, Cys579-Cys695, Cys657-Cys673, and Cys663-Cys685 in the catalytic domain, Cys351-Cys414, Cys464-Cys530, and Cys477-Cys543 in the scavenger receptor cysteine-rich (SRCR) domains. The disulfide bond of Cys351-Cys414 is predicted to play an allosteric function role, which could affect the metastasis and invasion of cancer cells. Other functional bonds have a catalytic role related to enzyme activity. The rest of disulfide bonds are predicted to play a structural role. Our study provides an important insight for the classification of disulfide bonds in LOXL2 and can be utilized for the drug design that targets the cysteine residues in LOXL2.

Molecular Simulation of αvβ6 Integrin Inhibitors

The urgent need for new treatments for the chronic lung disease idiopathic pulmonary fibrosis (IPF) motivates research into antagonists of the RGD binding integrin αvβ6, a protein linked to the initiation and progression of the disease. Molecular dynamics (MD) simulations of αvβ6 in complex with its natural ligand, pro-TGF-β1, show the persistence over time of a bidentate Arg-Asp ligand-receptor interaction and a metal chelate interaction between an aspartate on the ligand and an Mg²⁺ ion in the active site. This is typical of RGD binding ligands. Additional binding site interactions, which are not observed in the static crystal structure, are also identified. We investigate an RGD mimetic, which serves as a framework for a series of potential αvβ6 antagonists. The scaffold includes a derivative of the widely utilized 1,8-naphthyridine moiety, for which we present force field parameters, to enable MD and relative free energy perturbation (FEP) simulations. The MD simulations highlight the importance of hydrogen bonding and cation-π interactions. The FEP calculations predict relative binding affinities, within 1.5 kcal mol^-1, on average, of experiments.

Virtual Double-System Single-Box: A Nonequilibrium Alchemical Technique for Absolute Binding Free Energy Calculations: Application to Ligands of the SARS-CoV-2 Main Protease

In the context of drug-receptor binding affinity calculations using molecular dynamics techniques, we implemented a combination of Hamiltonian replica exchange (HREM) and a novel nonequilibrium alchemical methodology, called virtual double-system single-box, with increased accuracy, precision, and efficiency with respect to the standard nonequilibrium approaches. The method has been applied for the determination of absolute binding free energies of 16 newly designed noncovalent ligands of the main protease (3CL^pro) of SARS-CoV-2. The core structures of 3CL^pro ligands were previously identified using a multimodal structure-based ligand design in combination with docking techniques. The calculated binding free energies for four additional ligands with known activity (either for SARS-CoV or SARS-CoV-2 main protease) are also reported. The nature of binding in the 3CL^pro active site and the involved residues besides the CYS-HYS catalytic dyad have been thoroughly characterized by enhanced sampling simulations of the bound state. We have identified several noncongeneric compounds with predicted low micromolar activity for 3CL^pro inhibition, which may constitute possible lead compounds for the development of antiviral agents in Covid-19 treatment.

Inhibitor binding influences the protonation states of histidines in SARS-CoV-2 main protease

The main protease (M ^pro ) of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is an attractive target for antiviral therapeutics. Recently, many high-resolution apo and inhibitor-bound structures of M ^pro , a cysteine protease, have been determined, facilitating structure-based drug design. M ^pro plays a central role in the viral life cycle by catalyzing the cleavage of SARS-CoV-2 polyproteins. In addition to the catalytic dyad His41-Cys145, M ^pro contains multiple histidines including His163, His164, and His172. The protonation states of these histidines and the catalytic nu-cleophile Cys145 have been debated in previous studies of SARS-CoV M ^pro , but have yet to be investigated for SARS-CoV-2. In this work we have used molecular dynamics simulations to determine the structural stability of SARS-CoV-2 M ^pro as a function of the protonation assignments for these residues. We simulated both the apo and inhibitor-bound enzyme and found that the conformational stability of the binding site, bound inhibitors, and the hydrogen bond networks of M ^pro are highly sensitive to these assignments. Additionally, the two inhibitors studied, the peptidomimetic N3 and an α -ketoamide, display distinct His41/His164 protonation-state-dependent stabilities. While the apo and the N3-bound systems favored N _δ (HD) and N _ϵ (HE) protonation of His41 and His164, respectively, the α -ketoamide was not stably bound in this state. Our results illustrate the importance of using appropriate histidine protonation states to accurately model the structure and dynamics of SARS-CoV-2 M ^pro in both the apo and inhibitor-bound states, a necessary prerequisite for drug-design efforts.

Boosting Free-Energy Perturbation Calculations with GPU-Accelerated NAMD

Harnessing the power of graphics processing units (GPUs) to accelerate molecular dynamics (MD) simulations in the context of free-energy calculations has been a longstanding effort toward the development of versatile, high-performance MD engines. We report a new GPU-based implementation in NAMD of free-energy perturbation (FEP), one of the oldest, most popular importance-sampling approaches for the determination of free-energy differences that underlie alchemical transformations. Compared to the CPU implementation available since 2001 in NAMD, our benchmarks indicate that the new implementation of FEP in traditional GPU code is about four times faster, without any noticeable loss of accuracy, thereby paving the way toward more affordable free-energy calculations on large biological objects. Moreover, we have extended this new FEP implementation to a code path highly optimized for a single-GPU node, which proves to be up to nearly 30 times faster than the CPU implementation. Through optimized GPU performance, the present developments provide the community with a cost-effective solution for conducting FEP calculations. The new FEP-enabled code has been released with NAMD 3.0.

Scalable molecular dynamics on CPU and GPU architectures with NAMD

NAMDis a molecular dynamics program designed for high-performance simulations of very large biological objects on CPU- and GPU-based architectures. NAMD offers scalable performance on petascale parallel supercomputers consisting of hundreds of thousands of cores, as well as on inexpensive commodity clusters commonly found in academic environments. It is written in C++ and leans on Charm++ parallel objects for optimal performance on low-latency architectures. NAMD is a versatile, multipurpose code that gathers state-of-the-art algorithms to carry out simulations in apt thermodynamic ensembles, using the widely popular CHARMM, AMBER, OPLS, and GROMOS biomolecular force fields. Here, we review the main features of NAMD that allow both equilibrium and enhanced-sampling molecular dynamics simulations with numerical efficiency. We describe the underlying concepts utilized by NAMD and their implementation, most notably for handling long-range electrostatics; controlling the temperature, pressure, and pH; applying external potentials on tailored grids; leveraging massively parallel resources in multiple-copy simulations; and hybrid quantum-mechanical/molecular-mechanical descriptions. We detail the variety of options offered by NAMD for enhanced-sampling simulations aimed at determining free-energy differences of either alchemical or geometrical transformations and outline their applicability to specific problems. Last, we discuss the roadmap for the development of NAMD and our current efforts toward achieving optimal performance on GPU-based architectures, for pushing back the limitations that have prevented biologically realistic billion-atom objects to be fruitfully simulated, and for making large-scale simulations less expensive and easier to set up, run, and analyze. NAMD is distributed free of charge with its source code at www.ks.uiuc.edu.

Pseudo-Improper-Dihedral Model for Intrinsically Disordered Proteins

We present a new coarse-grained C_α-based protein model with a nonradial multibody pseudo-improper-dihedral potential that is transferable, time-independent, and suitable for molecular dynamics. It captures the nature of backbone and side-chain interactions between amino acid residues by adapting a simple improper dihedral term for a one-bead-per-residue model. It is parameterized for intrinsically disordered proteins and applicable to simulations of such proteins and their assemblies on millisecond time scales.

Absolute Binding Free Energy Calculation and Design of a Subnanomolar Inhibitor of Phosphodiesterase-10

Accurate prediction of absolute protein-ligand binding free energy could considerably enhance the success rate of structure-based drug design but is extremely challenging and time-consuming. Free energy perturbation (FEP) has been proven reliable but is limited to prediction of relative binding free energies of similar ligands (with only minor structural differences) in binding with a same drug target in practical drug design applications. Herein, a Gaussian algorithm-enhanced FEP (GA-FEP) protocol has been developed to enhance the FEP simulation performance, enabling to efficiently carry out the FEP simulations on vanishing the whole ligand and, thus, predict the absolute binding free energies (ABFEs). Using the GA-FEP protocol, the FEP simulations for the ABFE calculation (denoted as GA-FEP/ABFE) can achieve a satisfactory accuracy for both structurally similar and diverse ligands in a dataset of more than 100 receptor-ligand systems. Further, our GA-FEP/ABFE-guided lead optimization against phosphodiesterase-10 led to the discovery of a subnanomolar inhibitor (IC₅₀ = 0.87 nM, ∼2000-fold improvement in potency) with cocrystal confirmation.

Effect of the chemical environment of the DNA guanine quadruplex on the free energy of binding of Na and K ions

Guanine quadruplex (G-quadruplex) structures play a vital role in stabilizing the DNA genome and in protecting healthy cells from transforming into cancer cells. The structural stability of G-quadruplexes is greatly enhanced by the binding of monovalent cations such as Na⁺ or K⁺ into the interior axial channel. We computationally study the free energy of binding of Na⁺ and K⁺ ions to two intramolecular G-quadruplexes that differ considerably in their degree of rigidity and the presence or absence of terminal nucleotides. The goal of our study is two-fold. On the one hand, we study the free energy of binding every ion, which complements the experimental findings that report the average free energy for replacing Na⁺ with K⁺ ions. On the other hand, we examine the role of the G-quadruplex structure in the binding free energy. In the study, we employ all-atom molecular dynamics simulations and the alchemical transformation method for the computation of the free energies. To compare the cation-dependent contribution to the structural stability of G-quadruplexes, we use a two-step approach to calculate the individual free energy difference ΔG of binding two Na⁺ and two K⁺ to two G-quadruplexes: the unimolecular DNA d[T₂GT₂(G₃T)₃] (Protein Data Bank ID 2M4P) and the human telomeric DNA d[AGGG(TTAGGG)₃] (PDB ID 1KF1). In contrast to the experimental studies that estimate the average free energy of binding, we find a varying difference of approximately 2-9 kcal/mol between the free energy contribution of binding the first and second cation, Na⁺ or K⁺. Furthermore, we found that the free energy of binding K⁺ is not affected by the chemical nature of the two quadruplexes. By contrast, Na⁺ showed dependency on the G-quadruplex structure; the relatively small size allows Na⁺ to explore larger configurational space than K⁺. Numerical results presented here may offer reference values for future design of cationic drug-like ligands that replace the metal ions in G-quadruplexes.

Studying Functional Disulphide Bonds by Computer Simulations

Biochemical and structural data reveal important aspects of the properties and function of a protein disulphide bond. Molecular dynamics simulations can complement this experimental data and can yield valuable insights into the dynamical behavior of the disulphide bond within the protein environment. Due to the increasing accuracy of the underlying energetic description and the increasing computational power at hand, such simulations have now reached a level, at which they can also make quantitative and experimentally testable predictions. We here give an overview of the computational methods used to predict functional aspects of protein disulphides, including the prestress, protein allosteric effects upon thiol/disulphide exchange, and disulphide redox potentials. We then outline in detail the use of free-energy perturbation methods to calculate the redox potential of a protein disulphide bond of interest. In a step-by-step protocol, we describe the workflow within the MD suite Gromacs, including practical advice on the simulation setup and choice of parameters. For other disulphide-related simulation methods, we refer to resources available online.

Systematic Mutagenesis of Oncocin Reveals Enhanced Activity and Insights into the Mechanisms of Antimicrobial Activity

Oncocin is a proline-rich antimicrobial peptide that inhibits protein synthesis by binding to the bacterial ribosome. In this work, the antimicrobial activity of oncocin was improved by systematic peptide mutagenesis and activity evaluation. We found that a pair of cationic substitutions (P4K and L7K/R) improves the activity by 2-4 fold (p<0.05) against multiple Gram-negative bacteria. An in vitro transcription / translation assay indicated that the increased activity was not because of stronger ribosome binding. Rather a cellular internalization assay revealed a higher internalization rate for the optimized analogs thereby suggesting a mechanism to increase potency. In addition, we found that the optimized peptides' benefit is dependent upon nutrient-depleted media conditions. The molecular design and characterization strategies have broad potential for development of antimicrobial peptides.

Effect of Atomic Charges on Octanol-Water Partition Coefficient Using Alchemical Free Energy Calculation

The octanol-water partition coefficient (logPow) is an important index for measuring solubility, membrane permeability, and bioavailability in the drug discovery field. In this paper, the logPow values of 58 compounds were predicted by alchemical free energy calculation using molecular dynamics simulation. In free energy calculations, the atomic charges of the compounds are always fixed. However, they must be recalculated for each solvent. Therefore, three different sets of atomic charges were tested using quantum chemical calculations, taking into account vacuum, octanol, and water environments. The calculated atomic charges in the different environments do not necessarily influence the correlation between calculated and experimentally measured ∆Gwater values. The largest correlation coefficient values of the solvation free energy in water and octanol were 0.93 and 0.90, respectively. On the other hand, the correlation coefficient of logPow values calculated from free energies, the largest of which was 0.92, was sensitive to the combination of the solvation free energies calculated from the calculated atomic charges. These results reveal that the solvent assumed in the atomic charge calculation is an important factor determining the accuracy of predicted logPow values.

Constant-pH Molecular Dynamics Simulations for Large Biomolecular Systems

An increasingly important endeavor is to develop computational strategies that enable molecular dynamics (MD) simulations of biomolecular systems with spontaneous changes in protonation states under conditions of constant pH. The present work describes our efforts to implement the powerful constant-pH MD simulation method, based on a hybrid nonequilibrium MD/Monte Carlo (neMD/MC) technique within the highly scalable program NAMD. The constant-pH hybrid neMD/MC method has several appealing features; it samples the correct semigrand canonical ensemble rigorously, the computational cost increases linearly with the number of titratable sites, and it is applicable to explicit solvent simulations. The present implementation of the constant-pH hybrid neMD/MC in NAMD is designed to handle a wide range of biomolecular systems with no constraints on the choice of force field. Furthermore, the sampling efficiency can be adaptively improved on-the-fly by adjusting algorithmic parameters during the simulation. Illustrative examples emphasizing medium- and large-scale applications on next-generation supercomputing architectures are provided.

Free Energy Calculations of Microcin J25 Variants Binding to the FhuA Receptor

Computer simulations were performed to study the antimicrobial peptide microcin J25 (MJ25), a 21-mer peptide with an unusual lasso structure and high activity against Gram-negative bacteria. MJ25 has intracellular targets. The initial step of MJ25 acquisition in bacterial cells is binding to the outer-membrane receptor FhuA. Molecular dynamics simulations were implemented to study the binding mechanism of MJ25 to FhuA and to search for important binding residues. The absolute binding free energy calculated from combined free energy perturbation and thermodynamic integration methods agrees well with experimental data. In addition, computational mutation analysis revealed that His5 is the key residue responsible for MJ25 and FhuA association. We found that the number of hydrogen bonds is essential for binding of MJ25 to FhuA. This atomistic, quantitative insight sheds light on the mechanism of action of MJ25 and may pave a path for designing active MJ25 analogues.

Probing the stereospecificity of tyrosyl- and glutaminyl-tRNA synthetase with molecular dynamics

The stereospecificity of aminoacyl-tRNA synthetases helps exclude d-amino acids from protein synthesis and could perhaps be engineered to allow controlled d-amino acylation of tRNA. We use molecular dynamics simulations to probe the stereospecificity of the class I tyrosyl- and glutaminyl-tRNA synthetases (TyrRS, GlnRS), including wildtype enzymes and three point mutants suggested by three different protein design methods. l/d binding free energy differences are obtained by alchemically and reversibly transforming the ligand from L to D in simulations of the protein-ligand complex. The D81Q mutation in Escherichia coli TyrRS is homologous to the D81R mutant shown earlier to have inverted stereospecificity. D81Q is predicted to lead to a rotated ligand backbone and an increased, not a decreased l-Tyr preference. The E36Q mutation in Methanococcus jannaschii TyrRS has a predicted l/d binding free energy difference ΔΔG of just 0.5±0.9kcal/mol, compared to 3.1±0.8kcal/mol for the wildtype enzyme (favoring l-Tyr). The ligand ammonium position is preserved in the d-Tyr complex, while the carboxylate is shifted. Wildtype GlnRS has a similar preference for l-glutaminyl adenylate; the R260Q mutant has an increased preference, even though Arg260 makes a large contribution to the wildtype ΔΔG value.

Practically Efficient QM/MM Alchemical Free Energy Simulations: The Orthogonal Space Random Walk Strategy

The difference between free energy changes occurring at two chemical states can be rigorously estimated via alchemical free energy (AFE) simulations. Traditionally, most AFE simulations are carried out under the classical energy potential treatment; then, accuracy and applicability of AFE simulations are limited. In the present work, we integrate a recent second-order generalized ensemble strategy, the orthogonal space random walk (OSRW) method, into the combined quantum mechanical/molecular mechanical (QM/MM) potential based AFE simulation scheme. Thereby, within a commonly affordable simulation length, accurate QM/MM alchemical free energy simulations can be achieved. As revealed by the model study on the equilibrium of a tautomerization process of hydrated 3-hydroxypyrazole and by the model calculations of the redox potentials of two flavin derivatives, lumichrome (LC) and riboflavin (RF) in aqueous solution, the present OSRW-based scheme could be a viable path toward the realization of practically efficient QM/MM AFE simulations.

Accurate Estimates of Free Energy Changes in Charge Mutations

The ability to determine the effect of charge changes on the free energy is necessary for fundamental studies of the electrostatic contribution to protein binding and stability. Currently, calculations of differences in free energy for charge mutations carried out in systems with periodic boundary conditions must include an approximate self-energy correction that can be on the same order of magnitude as the calculated free energy change. Here, a new method for accurately calculating the free energy change associated with any alchemical mutation, regardless of charge, is presented. In this method, paired mutations of opposite charge exactly cancel the self-energy term because of its quadratic charge dependence. Since the self-energy term implicitly cancels within the method, a correction never needs to be applied, and the statistical uncertainty associated is thereby removed. An implementation procedure is described and applied to the free energy of ionic hydration and a charged amino acid mutation.

Interplay of Electrostatics and Hydrophobic Effects in the Metamorphic Protein Human Lymphotactin

The human lymphotactin (hLtn) is a protein that features two native states both of which are physiologically relevant: it is a monomer (hLtn10) at 10 °C with 200 mM salt and a dimer (hLtn40) at 40 °C and without salt. Here we focus on the networks of electrostatic and hydrophobic interactions that display substantial changes upon the conversion from hLtn10 to hLtn40 since they are expected to modulate the relative stability of the two folds. In addition to the Arg 23-Arg 43 interaction discussed in previous work, we find several other like-charge pairs that are likely important to the stability of hLtn10. Free energy perturbation calculations are carried out to explicitly evaluate the contribution of the Arg 23-Arg 43 interaction to the hLtn10 stability. hLtn40 features a larger number of salt bridges, and a set of hydrophobic residues undergo major changes in the solvent accessible surface area between hLtn10 and hLtn40, pointing to their importance to the relative stability of the two folds. We also discuss the use of explicit and implicit solvent simulations for characterizing the conformational ensembles under different solution conditions.