Self-guided Langevin dynamics via generalized Langevin equation

Intrinsic disorder and conformational coexistence in auxin coreceptors

AUXIN/INDOLE 3-ACETIC ACID (Aux/IAA) transcriptional repressor proteins and the TRANSPORT INHIBITOR RESISTANT 1/AUXIN SIGNALING F-BOX (TIR1/AFB) proteins to which they bind act as auxin coreceptors. While the structure of TIR1 has been solved, structural characterization of the regions of the Aux/IAA protein responsible for auxin perception has been complicated by their predicted disorder. Here, we use NMR, CD and molecular dynamics simulation to investigate the N-terminal domains of the Aux/IAA protein IAA17/AXR3. We show that despite the conformational flexibility of the region, a critical W-P bond in the core of the Aux/IAA degron motif occurs at a strikingly high (1:1) ratio of cis to trans isomers, consistent with the requirement of the cis conformer for the formation of the fully-docked receptor complex. We show that the N-terminal half of AXR3 is a mixture of multiple transiently structured conformations with a propensity for two predominant and distinct conformational subpopulations within the overall ensemble. These two states were modeled together with the C-terminal PB1 domain to provide the first complete simulation of an Aux/IAA. Using MD to recreate the assembly of each complex in the presence of auxin, both structural arrangements were shown to engage with the TIR1 receptor, and contact maps from the simulations match closely observations of NMR signal-decreases. Together, our results and approach provide a platform for exploring the functional significance of variation in the Aux/IAA coreceptor family and for understanding the role of intrinsic disorder in auxin signal transduction and other signaling systems.

A generic binding pocket for small molecule IKs activators at the extracellular inter-subunit interface of KCNQ1 and KCNE1 channel complexes

The cardiac IKs ion channel comprises KCNQ1, calmodulin, and KCNE1 in a dodecameric complex which provides a repolarizing current reserve at higher heart rates and protects from arrhythmia syndromes that cause fainting and sudden death. Pharmacological activators of IKs are therefore of interest both scientifically and therapeutically for treatment of IKs loss-of-function disorders. One group of chemical activators are only active in the presence of the accessory KCNE1 subunit and here we investigate this phenomenon using molecular modeling techniques and mutagenesis scanning in mammalian cells. A generalized activator binding pocket is formed extracellularly by KCNE1, the domain-swapped S1 helices of one KCNQ1 subunit and the pore/turret region made up of two other KCNQ1 subunits. A few residues, including K41, A44 and Y46 in KCNE1, W323 in the KCNQ1 pore, and Y148 in the KCNQ1 S1 domain, appear critical for the binding of structurally diverse molecules, but in addition, molecular modeling studies suggest that induced fit by structurally different molecules underlies the generalized nature of the binding pocket. Activation of IKs is enhanced by stabilization of the KCNQ1-S1/KCNE1/pore complex, which ultimately slows deactivation of the current, and promotes outward current summation at higher pulse rates. Our results provide a mechanistic explanation of enhanced IKs currents by these activator compounds and provide a map for future design of more potent therapeutically useful molecules.

Discovery of novel TRPV1 modulators through machine learning-based molecular docking and molecular similarity searching

The transient receptor potential vanilloid 1 (TRPV1) channel belongs to the transient receptor potential channel superfamily and participates in many physiological processes. TRPV1 modulators (both agonists and antagonists) can effectively inhibit pain caused by various factors and have curative effects in various diseases, such as itch, cancer, and cardiovascular diseases. Therefore, the development of TRPV1 channel modulators is of great importance. In this study, the structure-based virtual screening and ligand-based virtual screening methods were used to screen compound databases respectively. In the structure-based virtual screening route, a full-length human TRPV1 protein was first constructed, three molecular docking methods with different precisions were performed based on the hTRPV1 structure, and a machine learning-based rescoring model by the XGBoost algorithm was constructed to enrich active compounds. In the ligand-based virtual screening route, the ROCS program was used for 3D shape similarity searching and the EON program was used for electrostatic similarity searching. Final 77 compounds were selected from two routes for in vitro assays. The results showed that 8 of them were identified as active compounds, including three hits with IC50 values close to capsazepine. In addition, one hit is a partial agonist with both agonistic and antagonistic activity. The mechanisms of some active compounds were investigated by molecular dynamics simulation, which explained their agonism or antagonism.

Molecular Simulations of Conformational Transitions within the Insulin Receptor Kinase Reveal Consensus Features in a Multistep Activation Pathway

Modulating the transitions between active and inactive conformations of protein kinases is the primary means of regulating their catalytic activity, achieved by phosphorylation of the activation loop (A-loop). To elucidate the mechanism of this conformational activation, we applied the string method to determine the conformational transition path of insulin receptor kinase between the active and inactive conformations and the corresponding free-energy profiles with and without A-loop phosphorylation. The conformational change was found to proceed in three sequential steps: first, the flipping of the DFG motif of the active site; second, rotation of the A-loop; finally, the inward movement of the αC helix. The main energetic bottleneck corresponds to the conformational change in the A-loop, while changes in the DFG motif and αC helix occur before and after A-loop conformational change, respectively. In accordance with this, two intermediate states are identified, the first state just after the DFG flipping and the second state after the A-loop rotation. These intermediates exhibit structural features characteristic of the corresponding inactive and active conformations of other protein kinases. To understand the impact of A-loop phosphorylation on kinase conformation, the free energies of A-loop phosphorylation were determined at several states along the conformational transition path using the free-energy perturbation simulations. The calculated free energies reveal that while the unphosphorylated kinase interconverts between the inactive and active conformations, A-loop phosphorylation restricts access to the inactive conformation, thereby increasing the active conformation population. Overall, this study suggests a consensus mechanism of conformational activation between different protein kinases.

Deep Boosted Molecular Dynamics (DBMD): Accelerating molecular simulations with Gaussian boost potentials generated using probabilistic Bayesian deep neural network

We have developed a new Deep Boosted Molecular Dynamics (DBMD) method. Probabilistic Bayesian neural network models were implemented to construct boost potentials that exhibit Gaussian distribution with minimized anharmonicity, thereby allowing for accurate energetic reweighting and enhanced sampling of molecular simulations. DBMD was demonstrated on model systems of alanine dipeptide and the fast-folding protein and RNA structures. For alanine dipeptide, 30ns DBMD simulations captured up to 83-125 times more backbone dihedral transitions than 1μs conventional molecular dynamics (cMD) simulations and were able to accurately reproduce the original free energy profiles. Moreover, DBMD sampled multiple folding and unfolding events within 300ns simulations of the chignolin model protein and identified low-energy conformational states comparable to previous simulation findings. Finally, DBMD captured a general folding pathway of three hairpin RNAs with the GCAA, GAAA, and UUCG tetraloops. Based on Deep Learning neural network, DBMD provides a powerful and generally applicable approach to boosting biomolecular simulations. DBMD is available with open source in OpenMM at https://github.com/MiaoLab20/DBMD/.

Synthesis, in vitro biological screening and docking study of benzo[d]oxazole bis Schiff base derivatives as a potent anti-Alzheimer agent

We have synthesized benzo[d]oxazole derivatives (1-21) through a multistep reaction. Alteration in the structure of derivatives was brought in the last step via using various substituted aromatic aldehydes. In search of an anti-Alzheimer agent, all derivatives were evaluated against acetylcholinesterase and butyrylcholinesterase enzyme under positive control of standard drug donepezil (IC₅₀ = 0.016 ± 0.12 and 4.5 ± 0.11 µM) respectively. In case of acetylcholinesterase enzyme inhibition, derivatives 8, 9 and 18 (IC₅₀ = 0.50 ± 0.01, 0.90 ± 0.05 and 0.3 ± 0.05 µM) showed very promising inhibitory potentials. While in case of butyrylcholinesterase enzyme inhibition, most of the derivatives like 6, 8, 9, 13, 15, 18 and 19 (IC₅₀ = 2.70 ± 0.10, 2.60 ± 0.10, 2.20 ± 0.10, 4.25 ± 0.10, 3.30 ± 0.10, 0.96 ± 0.05 and 3.20 ± 0.10 µM) displayed better inhibitory potential than donepezil. Moreover, derivative 18 is the most potent one among the series in both inhibitions. The binding interaction of derivatives with the active gorge of the enzyme was confirmed via a docking study. Furthermore, the binding interaction between derivatives and the active site of enzymes was correlated through the SAR study. Structures of all derivatives were confirmed through spectroscopic techniques such as ¹H-NMR, ¹³C-NMR and HREI-MS, respectively.Communicated by Ramaswamy H. Sarma.

An effort to find new α-amylase inhibitors as potent antidiabetics compounds based on indole-based-thiadiazole analogs

Inhibition of α-amylase enzyme is of key significance for the therapy of diabetes mellitus (DM). Numerous indole-based compounds have earlier been described for broad range of bioactivities. From our previous study, we knew that indole and thiadiazole are potent inhibitors of diabetics II. We design the hybrid molecules of them and synthesized 18 derivatives of indole-based-thiadiazole (1-18). All synthesized compounds were characterized using different spectroscopic methods and evaluated for their α-amylase inhibitory activities. All synthetic compounds, except 4, 13, 15 and 16, were found to be strongly active (IC₅₀ values in the range of 0.80 ± 0.05 - 9.30 ± 0.20 µM) than the standard drug, acarbose (IC₅₀ = 11.70 ± 0.10 µM). Nevertheless, compound 18 was found to be inactive. The modes of binding interactions of five most active compounds 2, 3, 5, 10 and 17 were also studies through molecular docking study. In brief, current study identifies a novel class of α-amylase inhibitors which can be further studied for the treatment of hyperglycemia and obesity.Communicated by Ramaswamy H. Sarma.

Markov State Modeling Analysis Captures Changes in the Temperature-Sensitive N-Terminal and β-Turn Regions of the p53 DNA-Binding Domain

The transcription factor p53 is one of the most widely studied cancer proteins. Its temperature-sensitive nature suggests reduction in functionality at physiological temperatures. Temperature-induced conformational variations and their impact on its functional ability still remain unexplored. A total of 20.8 μs molecular dynamics simulations of wildtype p53 in the apo and the DNA-bound states have been performed at 300 K and 310 K. Further, Markov State Modeling (MSM) analyses were performed, considering C_α-C_α distances as reaction coordinates. Filtering of these distances based on correlation with the time-independent components (tICs) resulted in 16 and 32 distances for apo and DNA-bound systems, respectively. Individual MSM analyses using these filtered distances were performed for both p53 systems. These C_α-C_α residue pairs belonged to the N-terminal, S6/7 β-turn, loop L2, loop L3, and hydrophobic core residues. At physiological temperatures, apo-p53 exhibits exposure of its hydrophobic core, where the temperature-sensitive hotspot residues were also located. This exposure was the result of the S6/7 β-turn and N-terminal moving apart. In the DNA-bound p53 system, loop L1 attains an open conformation at physiological temperatures, which weakens the DNA binding. It is already known that p53 mutants that lack DNA binding also tend to show similar conformational variations. The S6/7 β-turn along with the already known functionally important loop L2 may pose as regions to be targeted to overcome the loss in binding of temperature-sensitive wildtype p53. Rescue strategies directed toward these temperature-sensitive regions may be useful to recuperate its strong binding at physiological temperatures.

Gaussian Accelerated Molecular Dynamics in OpenMM

Gaussian accelerated molecular dynamics (GaMD) is a computational technique that provides both unconstrained enhanced sampling and free energy calculations of biomolecules. Here, we present the implementation of GaMD in the OpenMM simulation package and validate it on model systems of alanine dipeptide and RNA folding. For alanine dipeptide, 30 ns GaMD production simulations reproduced free energy profiles of 1000 ns conventional molecular dynamics (cMD) simulations. In addition, GaMD simulations captured the folding pathways of three hyperstable RNA tetraloops (UUCG, GCAA, and CUUG) and binding of the rbt203 ligand to the HIV-1 Tar RNA, both of which involved critical electrostatic interactions such as hydrogen bonding and base stacking. Together with previous implementations, GaMD in OpenMM will allow for wider applications in simulations of proteins, RNA, and other biomolecules.

Application of molecular dynamics simulation in biomedicine

Molecular dynamics (MD) simulation has been widely used in the field of biomedicine to study the conformational transition of proteins caused by mutation or ligand binding/unbinding. It provides some perspectives those are difficult to find in traditional biochemical or pathological experiments, for example, detailed effects of mutations on protein structure and protein-protein/ligand interaction at the atomic level. In this review, a broad overview on conformation changes and drug discovery by MD simulation is given. We first discuss the preparation of protein structure for MD simulation, which is a key step that determines the accuracy of the simulation. Then, we summarize the applications of commonly used force fields and MD simulations in scientific research. Finally, enhanced sampling methods and common applications of these methods are introduced. In brief, MD simulation is a powerful tool and it can be used to guide experimental study. The combination of MD simulation and experimental techniques is an a priori means to solve the biomedical problems and give a deep understanding on the relationship between protein structure and function.

GLOW: A Workflow Integrating Gaussian-Accelerated Molecular Dynamics and Deep Learning for Free Energy Profiling

We introduce a Gaussian-accelerated molecular dynamics (GaMD), deep learning (DL), and free energy profiling workflow (GLOW) to predict molecular determinants and map free energy landscapes of biomolecules. All-atom GaMD-enhanced sampling simulations are first performed on biomolecules of interest. Structural contact maps are then calculated from GaMD simulation frames and transformed into images for building DL models using a convolutional neural network. Important structural contacts are further determined from DL models of attention maps of the structural contact gradients, which allow us to identify the system reaction coordinates. Finally, free energy profiles are calculated for the selected reaction coordinates through energetic reweighting of the GaMD simulations. We have also successfully demonstrated GLOW for the characterization of activation and allosteric modulation of a G protein-coupled receptor, using the adenosine A₁ receptor (A₁AR) as a model system. GLOW findings are highly consistent with previous experimental and computational studies of the A₁AR, while also providing further mechanistic insights into the receptor function. In summary, GLOW provides a systematic approach to mapping free energy landscapes of biomolecules. The GLOW workflow and its user manual can be downloaded at http://miaolab.org/GLOW.

Adaptive Synthesis of Functional Amphiphilic Dendrons as a Novel Approach to Artificial Supramolecular Objects

Supramolecular structures, such as micelles, liposomes, polymerosomes or dendrimerosomes, are widely studied and used as drug delivery systems. The behavior of amphiphilic building blocks strongly depends on their spatial distribution and shape of polar and nonpolar component. This report is focused on the development of new versatile synthetic protocols for amphiphilic carbosilane dendrons (amp-CS-DDNs) capable of self-assembly to regular micelles and other supramolecular objects. The presented strategy enables the fine modification of amphiphilic structure in several ways and also enables the facile connection of a desired functionality. DLS experiments demonstrated correlations between structural parameters of amp-CS-DDNs and the size of formed nanoparticles. For detailed information about the organization and spatial distribution of amp-CS-DDNs assemblies, computer simulation models were studied by using molecular dynamics in explicit water.

Carbosilane Glycodendrimers for Anticancer Drug Delivery: Synthetic Route, Characterization, and Biological Effect of Glycodendrimer-Doxorubicin Complexes

The complexity of drug delivery mechanisms calls for the development of new transport system designs. Here, we report a robust synthetic procedure toward stable glycodendrimer (glyco-DDM) series bearing glucose, galactose, and oligo(ethylene glycol)-modified galactose peripheral units. In vitro cytotoxicity assays showed exceptional biocompatibility of the glyco-DDMs. To demonstrate applicability in drug delivery, the anticancer agent doxorubicin (DOX) was encapsulated in the glyco-DDM structure. The anticancer activity of the resulting glyco-DDM/DOX complexes was evaluated on the noncancerous (BJ) and cancerous (MCF-7 and A2780) cell lines, revealing their promising generation- and concentration-dependent effect. The glyco-DDM/DOX complexes show gradual and pH-dependent DOX release profiles. Fluorescence spectra elucidated the encapsulation process. Confocal fluorescence microscopy demonstrated preferential cancer cell internalization of the glyco-DDM/DOX complexes. The conclusions were supported by computer modeling. Overall, our results are consistent with the assumption that novel glyco-DDMs and their drug complexes are very promising in drug delivery and related applications.

Binding Analysis Using Accelerated Molecular Dynamics Simulations and Future Perspectives

Biomolecular recognition such as binding of small molecules, nucleic acids, peptides and proteins to their target receptors plays key roles in cellular function and has been targeted for therapeutic drug design. Molecular dynamics (MD) is a computational approach to analyze these binding processes at an atomistic level, which provides valuable understandings of the mechanisms of biomolecular recognition. However, the rather slow biomolecular binding events often present challenges for conventional MD (cMD), due to limited simulation timescales (typically over hundreds of nanoseconds to tens of microseconds). In this regard, enhanced sampling methods, particularly accelerated MD (aMD), have proven useful to bridge the gap and enable all-atom simulations of biomolecular binding events. Here, we will review the recent method developments of Gaussian aMD (GaMD), ligand GaMD (LiGaMD) and peptide GaMD (Pep-GaMD), which have greatly expanded our capabilities to simulate biomolecular binding processes. Spontaneous binding of various biomolecules to their receptors has been successfully simulated by GaMD. Microsecond LiGaMD and Pep-GaMD simulations have captured repetitive binding and dissociation of small-molecule ligands and highly flexible peptides, and thus enabled ligand/peptide binding thermodynamics and kinetics calculations. We will also present relevant application studies in simulations of important drug targets and future perspectives for rational computer-aided drug design.

Synthesis of indole derivatives as diabetics II inhibitors and enzymatic kinetics study of α-glucosidase and α-amylase along with their in-silico study

In this study, we have investigated a series of indole-based compounds for their inhibitory study against pancreatic α-amylase and intestinal α-glucosidase activity. Inhibitors of carbohydrate degrading enzymes appear to have an essential role as antidiabetic drugs. All analogous exhibited good to moderate α-amylase (IC₅₀ = 3.80 to 47.50 μM), and α-glucosidase inhibitory interactions (IC₅₀ = 3.10-52.20 μM) in comparison with standard acarbose (IC₅₀ = 12.28 μM and 11.29 μM). The analogues 4, 11, 12, 15, 14 and 17 had good activity potential both for enzymes inhibitory interactions. Structure activity relationships were deliberated to propose the influence of substituents on the inhibitory potential of analogues. Docking studies revealed the interaction of more potential analogues and enzyme active site. Further, we studied their kinetic study of most active compounds showed that compounds 15, 14, 12, 17 and 11 are competitive for α-amylase and non- competitive for α-glucosidase.

Gaussian accelerated molecular dynamics (GaMD): principles and applications

Gaussian accelerated molecular dynamics (GaMD) is a robust computational method for simultaneous unconstrained enhanced sampling and free energy calculations of biomolecules. It works by adding a harmonic boost potential to smooth biomolecular potential energy surface and reduce energy barriers. GaMD greatly accelerates biomolecular simulations by orders of magnitude. Without the need to set predefined reaction coordinates or collective variables, GaMD provides unconstrained enhanced sampling and is advantageous for simulating complex biological processes. The GaMD boost potential exhibits a Gaussian distribution, thereby allowing for energetic reweighting via cumulant expansion to the second order (i.e., "Gaussian approximation"). This leads to accurate reconstruction of free energy landscapes of biomolecules. Hybrid schemes with other enhanced sampling methods, such as the replica exchange GaMD (rex-GaMD) and replica exchange umbrella sampling GaMD (GaREUS), have also been introduced, further improving sampling and free energy calculations. Recently, new "selective GaMD" algorithms including the ligand GaMD (LiGaMD) and peptide GaMD (Pep-GaMD) enabled microsecond simulations to capture repetitive dissociation and binding of small-molecule ligands and highly flexible peptides. The simulations then allowed highly efficient quantitative characterization of the ligand/peptide binding thermodynamics and kinetics. Taken together, GaMD and its innovative variants are applicable to simulate a wide variety of biomolecular dynamics, including protein folding, conformational changes and allostery, ligand binding, peptide binding, protein-protein/nucleic acid/carbohydrate interactions, and carbohydrate/nucleic acid interactions. In this review, we present principles of the GaMD algorithms and recent applications in biomolecular simulations and drug design.

Mechanism of zinc ejection by disulfiram in nonstructural protein 5A

Hepatitis C virus (HCV) is a notorious member of the Flaviviridae family of enveloped, positive-strand RNA viruses. Non-structural protein 5A (NS5A) plays a key role in HCV replication and assembly. NS5A is a multi-domain protein which includes an N-terminal amphipathic membrane anchoring alpha helix, a highly structured domain-1, and two intrinsically disordered domains 2-3. The highly structured domain-1 contains a zinc finger (Zf)-site, and binding of zinc stabilizes the overall structure, while ejection of this zinc from the Zf-site destabilizes the overall structure. Therefore, NS5A is an attractive target for anti-HCV therapy by disulfiram, through ejection of zinc from the Zf-site. However, the zinc ejection mechanism is poorly understood. To disclose this mechanism based on three different states, A-state (NS5A protein), B-state (NS5A + Zn), and C-state (NS5A + Zn + disulfiram), we have performed molecular dynamics (MD) simulation in tandem with DFT calculations in the current study. The MD results indicate that disulfiram triggers Zn ejection from the Zf-site predominantly through altering the overall conformation ensemble. On the other hand, the DFT assessment demonstrates that the Zn adopts a tetrahedral configuration at the Zf-site with four Cys residues, which indicates a stable protein structure morphology. Disulfiram binding induces major conformational changes at the Zf-site, introduces new interactions of Cys39 with disulfiram, and further weakens the interaction of this residue with Zn, causing ejection of zinc from the Zf-site. The proposed mechanism elucidates the therapeutic potential of disulfiram and offers theoretical guidance for the advancement of drug candidates.

Designing a multi-epitope vaccine against Mycobacteroides abscessus by pangenome-reverse vaccinology

Mycobacteroides abscessus (Previously Mycobacterium abscessus) is an emerging microorganism of the newly defined genera Mycobacteroides that causes mainly skin and tissue diseases in humans. The recent availability of total 34 fully sequenced genomes of different strains belonging to this species has provided an opportunity to utilize this genomics data to gain novel insights and guide the development of specific antimicrobial therapies. In the present study, we collected collectively 34 complete genome sequences of M. abscessus from the NCBI GenBank database. Pangenome analysis was conducted on these genomes to understand the genetic diversity and to obtain proteins associated with its core genome. These core proteins were then subjected to various subtractive filters to identify potential antigenic targets that were subjected to multi-epitope vaccine design. Our analysis projected the open pangenome of M. abscessus containing 3443 core genes. After applying various stepwise filtration steps on the core proteins, a total of four potential antigenic targets were identified. Utilizing their constituent CD4 and CD8 T-cell epitopes, a multi-epitope based subunit vaccine was computationally designed. Sequence-based analysis as well as structural characterization revealed the immunological effectiveness of this designed vaccine. Further molecular docking, molecular dynamics simulation and binding free energy estimation with Toll-like receptor 2 indicated strong structural associations of the vaccine with the immune receptor. The promising results are encouraging and need to be validated by additional wet laboratory studies for confirmation.

Distal lysine (de)coordination in the algal hemoglobin THB1: A combined computer simulation and experimental study

THB1 is a monomeric truncated hemoglobin from the green alga Chlamydomonas reinhardtii. In the absence of exogenous ligands and at neutral pH, the heme group of THB1 is coordinated by two protein residues, Lys53 and His77. THB1 is thought to function as a nitric oxide dioxygenase, and the distal binding of O₂ requires the cleavage of the Fe-Lys53 bond accompanied by protonation and expulsion of the lysine from the heme cavity into the solvent. Nuclear magnetic resonance spectroscopy and crystallographic data have provided dynamic and structural insights of the process, but the details of the mechanism have not been fully elucidated. We applied a combination of computer simulations and site-directed mutagenesis experiments to shed light on this issue. Molecular dynamics simulations and hybrid quantum mechanics/molecular mechanics restrained optimizations were performed to explore the nature of the transition between the decoordinated and lysine-bound states of the ferrous heme in THB1. Lys49 and Arg52, which form ionic interactions with the heme propionates in the X-ray structure of lysine-bound THB1, were observed to assist in maintaining Lys53 inside the protein cavity and play a key role in the transition. Lys49Ala, Arg52Ala and Lys49Ala/Arg52Ala THB1 variants were prepared, and the consequences of the replacements on the Lys (de)coordination equilibrium were characterized experimentally for comparison with computational prediction. The results reinforced the dynamic role of protein-propionate interactions and strongly suggested that cleavage of the Fe-Lys53 bond and ensuing conformational rearrangement is facilitated by protonation of the amino group inside the distal cavity.

Control of distal lysine coordination in a monomeric hemoglobin: A role for heme peripheral interactions

THB1 is a monomeric truncated hemoglobin (TrHb) found in the cytoplasm of the green alga Chlamydomonas reinhardtii. The canonical heme coordination scheme in hemoglobins is a proximal histidine ligand and an open distal site. In THB1, the latter site is occupied by Lys53, which is likely to facilitate Fe(II)/Fe(III) redox cycling but hinders dioxygen binding, two features inherent to the NO dioxygenase activity of the protein. TrHb surveys show that a lysine at a position aligning with Lys53 is an insufficient determinant of coordination, and in this study, we sought to identify factors controlling lysine affinity for the heme iron. We solved the "Lys-off" X-ray structure of THB1, represented by the cyanide adduct of the Fe(III) protein, and hypothesized that interactions that differ between the known "Lys-on" structure and the Lys-off structure participate in the control of Lys53 affinity for the heme iron. We applied an experimental approach (site-directed mutagenesis, heme modification, pH titrations in the Fe(III) and Fe(II) states) and a computational approach (MD simulations in the Fe(II) state) to assess the role of heme propionate-protein interactions, distal helix capping, and the composition of the distal pocket. All THB1 modifications resulted in a weakening of lysine affinity and affected the coupling between Lys53 proton binding and heme redox potential. The results supported the importance of specific heme peripheral interactions for the pH stability of iron coordination and the ability of the protein to undergo redox reactions.

Dihydroquinazolin-4(1H)-one derivatives as novel and potential leads for diabetic management

A variety of dihydroquinazolin-4(1H)-one derivatives (1-37) were synthesized via "one-pot" three-component reaction scheme by treating aniline and different aromatic aldehydes with isatoic anhydride in the presence of acetic acid. Chemical structures of compounds were deduced by different spectroscopic techniques including EI-MS, HREI-MS, ¹H-, and ¹³C-NMR. Compounds were subjected to α-amylase and α-glucosidase inhibitory activities. A number of derivatives exhibited significant to moderate inhibition potential against α-amylase (IC₅₀ = 23.33 ± 0.02-88.65 ± 0.23 μM) and α-glucosidase (IC₅₀ = 25.01 ± 0.12-89.99 ± 0.09 μM) enzymes, respectively. Results were compared with the standard acarbose (IC₅₀ = 17.08 ± 0.07 μM for α-amylase and IC₅₀ = 17.67 ± 0.09 μM for α-glucosidase). Structure-activity relationship (SAR) was rationalized by analyzing the substituents effects on inhibitory potential. Kinetic studies were implemented to find the mode of inhibition by compounds which revealed competitive inhibition for α-amylase and non-competitive inhibition for α-glucosidase. However, in silico study identified several important binding interactions of ligands (synthetic analogues) with the active site of both enzymes.

Post-translational modifications in collagen type I of bone in a mouse model of aging

The fracture resistance of cortical bone and matrix hydration are known to decline with advanced aging. However, the underlying mechanisms remain poorly understood, and so we investigated levels of matrix proteins and post-translational modifications (PTM) of collagen I in extracts from the tibia of 6-mo. and 20-mo. old BALB/c mice (female and male analysis done separately). Liquid chromatography-tandem mass spectrometry (LC-MS/MS) analysis revealed that the levels of collagen I deamidation at specific asparagine (Asn) and glutamine (Gln) residues significantly increased with age. Other non-enzymatic PTMs such as carboxymethylation of lysine (CML) were detected as well, but the relative abundance did not vary with age. No significant age-related differences in the abundance of hydroxylysine glycosylation sites were found, but hydroxylation levels at a few of the numerous lysine and proline hydroxylation sites significantly changed by a small amount with age. We performed molecular modeling and dynamics (MD) simulations for three triple helical fragments representing collagen I regions with prominent age-dependent increases in deamidation as identified by LC-MS/MS of male extracts. These 3 fragments included deamidated Asn and Gln residues as follows: 1) an Asn⁴²⁸ site of the α2(I) chain in which deamidation levels increased from 4.4% at 6-mo. to 8.1% at 20-mo., 2) an Asn⁹⁸³ site of the α2(I) chain with a deamidation increase from 18.3% to 36.8% with age and an Asn¹⁰⁵² site of the α1(I) chain with consistent deamidation levels of ~60% across the age groups, and 3) a Gln⁴¹⁰ site of the α1(I) chain that went from no detectable deamidation at 6-mo. to 2.7% at 20-mo. and a neighboring Asn⁴²¹ site of the same chain with an age-related deamidation increase from 3.6% to 16.3%. The deamidation levels at these sites inversely correlated with an estimate of toughness determined from three-point bending tests of the femur mid-diaphysis. MD revealed that the sidechains become more negatively charged at deamidated sites and that deamidation alters hydrogen bonding with water along the collagen backbone while increasing water interactions with the aspartic and glutamic acid sidechains. Our findings suggest a new mechanism of the age-dependent reduction in the fracture resistance of cortical bone whereby deamidation of Asn and Glu residues redistributes bound water within collagen I triple helix.

A novel multiscale scheme to accelerate atomistic simulations of bio-macromolecules by adaptively driving coarse-grained coordinates

All-atom molecular dynamics (MD) simulations of bio-macromolecules can yield relatively accurate results while suffering from the limitation of insufficient conformational sampling. On the other hand, the coarse-grained (CG) MD simulations efficiently accelerate conformational changes in biomolecules but lose atomistic details and accuracy. Here, we propose a novel multiscale simulation method called the adaptively driving multiscale simulation (ADMS)-it efficiently accelerates biomolecular dynamics by adaptively driving virtual CG atoms on the fly while maintaining the atomistic details and focusing on important conformations of the original system with irrelevant conformations rarely sampled. Herein, the "adaptive driving" is based on the short-time-averaging response of the system (i.e., an approximate free energy surface of the original system), without requiring the construction of the CG force field. We apply the ADMS to two peptides (deca-alanine and Ace-GGPGGG-Nme) and one small protein (HP35) as illustrations. The simulations show that the ADMS not only efficiently captures important conformational states of biomolecules and drives fast interstate transitions but also yields, although it might be in part, reliable protein folding pathways. Remarkably, a ∼100-ns explicit-solvent ADMS trajectory of HP35 with three CG atoms realizes folding and unfolding repeatedly and captures the important states comparable to those from a 398-µs standard all-atom MD simulation.

Reformulation of the self-guided molecular simulation method

Self-guided molecular/Langevin dynamics (SGMD/SGLD) simulation methods were developed to enhance conformational sampling through promoting low frequency motion of molecular systems and have been successfully applied in many simulation studies. Quantitative understanding of conformational distribution in SGLD has been achieved by separating microscopic properties according to frequency. However, a missing link between the guiding factors and conformational distributions makes it highly empirical and system dependent when choosing the values of the guiding parameters. Based on the understanding that molecular interactions are the source of energy barriers and diffusion friction, this work reformulates the equation of the low frequency motion to resemble Langevin dynamics. This reformulation leads to new forms of guiding forces and establishes a relation between the guiding factors and conformational distributions. We call simulations with these new guiding forces the generalized self-guided molecular/Langevin dynamics (SGMDg/SGLDg). In addition, we present a new way to calculate low frequency properties and an efficient algorithm to implement SGMDg/SGLDg that minimizes memory usage and inter-processor communication. Through example simulations with a skewed double well system, an argon fluid, and a cryo-EM map flexible fitting case, we demonstrate the guiding effects on conformational distributions and conformational searching.

In vivo analgesic, anti-inflammatory, and sedative activity and a molecular docking study of dinaphthodiospyrol G isolated from Diospyros lotus

BACKGROUND

Analgesic, anti-inflammatory, and sedative drugs are available with potential side effects such as peptic ulcer and addiction among other things. In this regard, research is underway to find safe, effective, and economical drugs free of these side effects. In this study, an isolated natural product from Diospyros lotus, was tested for the aforementioned bioactivities.

OBJECTIVES

To evaluate analgesic, anti-inflammatory, and sedative potential of D. lotus extracts in animal paradigms using BALB/c mice as experimental model.

METHODS

Analgesic, anti-inflammatory and sedative activities of dinaphthodiospyrol G (1) isolated from the chloroform fraction of D. lotus were evaluated using different experimental procedures. Anti-inflammatory effect was evaluated using the carrageenan and histamine-induced paw edema, whereas the antinociceptive effect was quantified by means of the hot plate analgesiometer. On the other hand, the sedative effect was determined using animal assay for screening the locomotors effects of compound 1. Compound 1 was also subjected to molecular modeling studies against cyclooxygenase enzymes.

RESULTS

Results from this investigation showed that the extract is devoid of anti-inflammatory and antinociceptive potentials but has a significant sedative effect, whereas the tested compound exhibited 55.23 and 78.34% attenuation in paw edema by carrageenan and histamine assays, respectively. A significant (p < 0.001) and dose-dependent antinociceptive and sedative effects were demonstrated by the isolated compound. Molecular docking and dynamics simulation studies of the isolated compound against cyclooxygenase enzyme indicated that compound 1 forms specific interactions with key residues in the active site of the target receptor, which validates the potential use of the isolated compound as cyclooxygenase inhibitor.

CONCLUSIONS

Compound 1 exhibited remarkable analgesic, anti-inflammatory, and sedative activities. These findings strongly justify the traditional use of D. lotus in the treatment of inflammation, pain, and insomnia.

Enhanced sampling and free energy calculations for protein simulations

Molecular dynamics simulation is a powerful computational technique to study biomolecular systems, which complements experiments by providing insights into the structural dynamics relevant to biological functions at atomic scale. It can also be used to calculate the free energy landscapes of the conformational transitions to better understand the functions of the biomolecules. However, the sampling of biomolecular configurations is limited by the free energy barriers that need to be overcome, leading to considerable gaps between the timescales reached by MD simulation and those governing biological processes. To address this issue, many enhanced sampling methodologies have been developed to increase the sampling efficiency of molecular dynamics simulations and free energy calculations. Usually, enhanced sampling algorithms can be classified into methods based on collective variables (CV-based) and approaches which do not require predefined CVs (CV-free). In this chapter, the theoretical basis of free energy estimation is briefly reviewed first, followed by the reviews of the most common CV-based and CV-free methods including the presentation of some examples and recent developments. Finally, the combination of different enhanced sampling methods is discussed.

Novel N-Acyl-1H-imidazole-1-carbothioamides: Design, Synthesis, Biological and Computational Studies

The present study reports the convenient synthesis, spectroscopic characterization, bio-assays and computational evaluation of a novel series of N-acyl-1H-imidazole-1-carbothioamides. The screened derivatives displayed excellent antioxidant activity, moderate antibacterial and antifungal potential. The screened derivatives were found to be highly biocompatible against hRBCs. Molecular docking ascertained the mechanism and mode of action towards the molecular target delineating that ligands and complexes were stabilized at the active site by electrostatic and hydrophobic forces in accordance to the corresponding experimental results. Docking simulation provided additional information about the possibilities of inhibitory potential of the compounds against RNA. Computational evaluation predicted that N-acyl-1H-imidazole-1-carbothioamides 5c and 5g can serve as potential surrogates for hit to lead generation and design of novel antioxidant and antibacterial agents.

DFMD: Fast and Effective DelPhiForce Steered Molecular Dynamics Approach to Model Ligand Approach Toward a Receptor: Application to Spermine Synthase Enzyme

Here we report a novel approach, the DelPhiForce Molecular Dynamics (DFMD) method, for steered molecular dynamics simulations to model receptor-ligand association involving charged species. The main purpose of developing DFMD is to simulate ligand's trajectory toward the receptor and thus to predict the "entrance" of the binding pocket and conformational changes associated with the binding. We demonstrate that the DFMD is superior compared with molecular dynamics simulations applying standard cut-offs, provides correct binding forces, allows for modeling the ligand approach at long distances and thus guides the ligand toward the correct binding spot, and it is very fast (frequently the binding is completed in <1 ns). The DFMD is applied to model the binding of two ligands to a receptor (spermine synthase) and it is demonstrated that it guides the ligands toward the corresponding pockets despite of the initial ligand's position with respect to the receptor. Predicted conformational changes and the order of ligand binding are experimentally verified.

Thermal Stability of Globins: Implications of Flexibility and Heme Coordination Studied by Molecular Dynamics Simulations

Proteins are sensitive to temperature, and abrupt changes in the normal temperature conditions can have a profound impact on both structure and function, leading to protein unfolding. However, the adaptation of certain organisms to extreme conditions raises questions about the structural features that permit the structure and function of proteins to be preserved under these adverse conditions. To gain insight into the molecular basis of protein thermostability in the globin family, we have examined three representative examples: human neuroglobin, horse heart myoglobin, and Drosophila hemoglobin, which differ in their melting temperatures and coordination states of the heme iron in the absence of external ligands. In order to elucidate the possible mechanisms that govern the thermostability of these proteins, microsecond-scale classical molecular dynamics simulations were performed at different temperatures. Structural fluctuations and essential dynamics were analyzed, indicating that the flexibility of the CD region, which includes the two short C and D helixes and the connecting CD loop, is directly related to the thermostability. We observed that a larger inherent flexibility of the protein produces higher thermostability, probably concentrating the thermal fluctuations observed at high temperature in flexible regions, preventing unfolding. Globally, the results of this work improve our understanding of thermostability in the globin family.

Comparison of the umbrella sampling and the double decoupling method in binding free energy predictions for SAMPL6 octa-acid host-guest challenges

We calculate the absolute binding free energies of tetra-methylated octa-acids host-guest systems as a part of the SAMPL6 blind challenge (receipt ID vq30p). We employed two different free energy simulation methods, i.e., the umbrella sampling (US) and double decoupling method (DDM). The US method was used with the weighted histogram analysis method (WHAM) (US-WHAM scheme). In the DDM scheme, Hamiltonian replica-exchange method (HREM) was combined with the Bennett acceptance ratio (BAR) (HREM-BAR scheme). We obtained initial binding poses via molecular docking using GalaxyDock-HG program, which is developed for the SAMPL challenge. The root mean square deviation (RMSD) and the mean absolute deviations (MAD) using US-WHAM scheme were 1.33 and 1.02 kcal/mol, respectively. The MAD was the top among all submissions, however the correlation with respect to experiment was unexceptional. While the RMSD and MAD via HREM-BAR scheme were greater than US-WHAM scheme, (i.e., 2.09 and 1.76 kcal/mol), their correlations were slightly better than US-WHAM. The correlation between the two methods was high. Further discussion on the DDM method can be found in a companion paper by Han et al. (receipt ID 3z83m) in the same issue.

Conformational Selection of a Polyproline Peptide by Ebola Virus VP30

An adaptive temperature-based replica-exchange simulation of a peptide extracted from the Ebola virus nucleoprotein containing a polyproline sequence motif is reported. The simulation results of applying the CHARMM36m force field with a generalized Born solvent model is presented. Conformational heterogeneity is described by potentials of mean force (PMFs) for a set of reaction coordinates that define the topological fold space. Starting from an extended backbone conformation of the peptide observed in an X-ray crystallographic assembly with the Ebola virus protein VP30, the PMFs report a conformational landscape populated by chain excursions to collapsed states with limited transitions to either an extended fold or a canonical polyproline type II helix. Clustering of the conformations and applying an elastic network interpolation model yield a multistep pathway of conformational selection that minimizes the net transition-state cost from the population hub to the bound state. Related difference between the pathway endpoints taken from the PMFs reveal a significant free-energy penalty in reaching a population shift. To evaluate sequence fitness of the Ebola virus peptide in generating probability distributions, two human sequence variants are modeled and are found to produce profiles that show extensive deviations, thus suggesting either dissimilar binding mechanisms or the lack of recognition by VP30.

Hydronium Ions Accompanying Buried Acidic Residues Lead to High Apparent Dielectric Constants in the Interior of Proteins

Internal ionizable groups are known to play important roles in protein functions. A mystery that has attracted decades of extensive experimental and theoretical studies is the apparent dielectric constants experienced by buried ionizable groups, which are much higher than values expected for protein interiors. Many interpretations have been proposed, such as water penetration, conformational relaxation, local unfolding, protein intrinsic backbone fluctuations, etc. However, these interpretations conflict with many experimental observations. The virtual mixture of multiple states (VMMS) simulation method developed in our lab provides a direct approach for studying the equilibrium of multiple chemical states and can monitor p K_a values along simulation trajectories. Through VMMS simulations of staphylococcal nuclease (SNase) variants with internal Asp or Glu residues, we discovered that cations were attracted to buried deprotonated acidic groups and the presence of the nearby cations were essential to reproduce experimentally measured p K_a values. This finding, combined with structural analysis and validation simulations, suggests that the proton released from a deprotonation process stays near the deprotonated group inside proteins, possibly in the form of a hydronium ion. The existence of a proton near a buried charge has many implications in our understanding of protein functions.

Atomic Resolution Cryo-EM Structure of β-Galactosidase

The advent of direct electron detectors has enabled the routine use of single-particle cryo-electron microscopy (EM) approaches to determine structures of a variety of protein complexes at near-atomic resolution. Here, we report the development of methods to account for local variations in defocus and beam-induced drift, and the implementation of a data-driven dose compensation scheme that significantly improves the extraction of high-resolution information recorded during exposure of the specimen to the electron beam. These advances enable determination of a cryo-EM density map for β-galactosidase bound to the inhibitor phenylethyl β-D-thiogalactopyranoside where the ordered regions are resolved at a level of detail seen in X-ray maps at ∼ 1.5 Å resolution. Using this density map in conjunction with constrained molecular dynamics simulations provides a measure of the local flexibility of the non-covalently bound inhibitor and offers further opportunities for structure-guided inhibitor design.

Molecular dynamics simulations of lipid nanodiscs

A lipid nanodisc is a discoidal lipid bilayer stabilized by proteins, peptides, or polymers on its edge. Nanodiscs have two important connections to structural biology. The first is associated with high-density lipoprotein (HDL), a particle with a variety of functionalities including lipid transport. Nascent HDL (nHDL) is a nanodisc stabilized by Apolipoprotein A-I (APOA1). Determining the structure of APOA1 and its mimetic peptides in nanodiscs is crucial to understanding pathologies related to HDL maturation and designing effective therapies. Secondly, nanodiscs offer non-detergent membrane-mimicking environments and greatly facilitate structural studies of membrane proteins. Although seemingly similar, natural and synthetic nanodiscs are different in that nHDL is heterogeneous in size, due to APOA1 elasticity, and gradually matures to become spherical. Synthetic nanodiscs, in contrast, should be homogenous, stable, and size-tunable. This report reviews previous molecular dynamics (MD) simulation studies of nanodiscs and illustrates convergence and accuracy issues using results from new multi-microsecond atomistic MD simulations. These new simulations reveal that APOA1 helices take 10-20 μs to rearrange on the nanodisc, while peptides take 2 μs to migrate from the disc surfaces to the edge. These systems can also become kinetically trapped depending on the initial conditions. For example, APOA1 was trapped in a biologically irrelevant conformation for the duration of a 10 μs trajectory; the peptides were similarly trapped for 5 μs. It therefore remains essential to validate MD simulations of these systems with experiments due to convergence and accuracy issues. This article is part of a Special Issue entitled: Emergence of Complex Behavior in Biomembranes edited by Marjorie Longo.

Ensemble Docking in Drug Discovery

Ensemble docking corresponds to the generation of an "ensemble" of drug target conformations in computational structure-based drug discovery, often obtained by using molecular dynamics simulation, that is used in docking candidate ligands. This approach is now well established in the field of early-stage drug discovery. This review gives a historical account of the development of ensemble docking and discusses some pertinent methodological advances in conformational sampling.

Hamiltonian Monte Carlo with Constrained Molecular Dynamics as Gibbs Sampling

Compared to fully flexible molecular dynamics, simulations of constrained systems can use larger time steps and focus kinetic energy on soft degrees of freedom. Achieving ergodic sampling from the Boltzmann distribution, however, has proven challenging. Using recent generalizations of the equipartition principle and Fixman potential, here we implement Hamiltonian Monte Carlo based on constrained molecular dynamics as a Gibbs sampling move. By mixing Hamiltonian Monte Carlo based on fully flexible and torsional dynamics, we are able to reproduce free energy landscapes of simple model systems and enhance sampling of macrocycles.

On the Helix Propensity in Generalized Born Solvent Descriptions of Modeling the Dark Proteome

Intrinsically disordered proteins that populate the so-called “Dark Proteome” offer challenging benchmarks of atomistic simulation methods to accurately model conformational transitions on a multidimensional energy landscape. This work explores the application of parallel tempering with implicit solvent models as a computational framework to capture the conformational ensemble of an intrinsically disordered peptide derived from the Ebola virus protein VP35. A recent X-ray crystallographic study reported a protein-peptide interface where the VP35 peptide underwent a folding transition from a disordered form to a helix-β-turn-helix topological fold upon molecular association with the Ebola protein NP. An assessment is provided of the accuracy of two generalized Born solvent models (GBMV2 and GBSW2) using the CHARMM force field and applied with temperature-based replica exchange dynamics to calculate the disorder propensity of the peptide and its probability density of states in a continuum solvent. A further comparison is presented of applying an explicit/implicit solvent hybrid replica exchange simulation of the peptide to determine the effect of modeling water interactions at the all-atom resolution.

Understanding ligand-receptor non-covalent binding kinetics using molecular modeling

Kinetic properties may serve as critical differentiators and predictors of drug efficacy and safety, in addition to the traditionally focused binding affinity. However the quantitative structure-kinetics relationship (QSKR) for modeling and ligand design is still poorly understood. This review provides an introduction to the kinetics of drug binding from a fundamental chemistry perspective. We focus on recent developments of computational tools and their applications to non-covalent binding kinetics.

Origin of pKa Shifts of Internal Lysine Residues in SNase Studied Via Equal-Molar VMMS Simulations in Explicit Water

Protein internal ionizable groups can exhibit large shifts in pK_a values. Although the environment and interaction changes have been extensively studied both experimentally and computationally, direct calculation of pK_a values of these internal ionizable groups in explicit water is challenging due to energy barriers in solvent interaction and in conformational transition. The virtual mixture of multiple states (VMMS) method is a new approach designed to study chemical state equilibrium. This method constructs a virtual mixture of multiple chemical states in order to sample the conformational space of all states simultaneously and to avoid crossing energy barriers related to state transition. By applying VMMS to 25 variants of staphylococcal nuclease with lysine residues at internal positions, we obtained the pK_a values of these lysine residues and investigated the physics underlining the pK_a shifts. Our calculation results agree reasonably well with experimental measurements, validating the VMMS method for pK_a calculation and providing molecular details of the protonation equilibrium for protein internal ionizable groups. Based on our analyses of protein conformation relaxation, lysine side chain flexibility, water penetration, and the microenvironment, we conclude that the hydrophobicity of the microenvironment around the lysine side chain (which affects water penetration differently for different protonation states) plays an important role in the pK_a shifts.

Unconstrained Enhanced Sampling for Free Energy Calculations of Biomolecules: A Review

Free energy calculations are central to understanding the structure, dynamics and function of biomolecules. Yet insufficient sampling of biomolecular configurations is often regarded as one of the main sources of error. Many enhanced sampling techniques have been developed to address this issue. Notably, enhanced sampling methods based on biasing collective variables (CVs), including the widely used umbrella sampling, adaptive biasing force and metadynamics, have been discussed in a recent excellent review (Abrams and Bussi, Entropy, 2014). Here, we aim to review enhanced sampling methods that do not require predefined system-dependent CVs for biomolecular simulations and as such do not suffer from the hidden energy barrier problem as encountered in the CV-biasing methods. These methods include, but are not limited to, replica exchange/parallel tempering, self-guided molecular/Langevin dynamics, essential energy space random walk and accelerated molecular dynamics. While it is overwhelming to describe all details of each method, we provide a summary of the methods along with the applications and offer our perspectives. We conclude with challenges and prospects of the unconstrained enhanced sampling methods for accurate biomolecular free energy calculations.

Investigation of Structural Dynamics of Enzymes and Protonation States of Substrates Using Computational Tools

This review discusses the use of molecular modeling tools, together with existing experimental findings, to provide a complete atomic-level description of enzyme dynamics and function. We focus on functionally relevant conformational dynamics of enzymes and the protonation states of substrates. The conformational fluctuations of enzymes usually play a crucial role in substrate recognition and catalysis. Protein dynamics can be altered by a tiny change in a molecular system such as different protonation states of various intermediates or by a significant perturbation such as a ligand association. Here we review recent advances in applying atomistic molecular dynamics (MD) simulations to investigate allosteric and network regulation of tryptophan synthase (TRPS) and protonation states of its intermediates and catalysis. In addition, we review studies using quantum mechanics/molecular mechanics (QM/MM) methods to investigate the protonation states of catalytic residues of β-Ketoacyl ACP synthase I (KasA). We also discuss modeling of large-scale protein motions for HIV-1 protease with coarse-grained Brownian dynamics (BD) simulations.

A Virtual Mixture Approach to the Study of Multistate Equilibrium: Application to Constant pH Simulation in Explicit Water

Chemical and thermodynamic equilibrium of multiple states is a fundamental phenomenon in biology systems and has been the focus of many experimental and computational studies. This work presents a simulation method to directly study the equilibrium of multiple states. This method constructs a virtual mixture of multiple states (VMMS) to sample the conformational space of all chemical states simultaneously. The VMMS system consists of multiple subsystems, one for each state. The subsystem contains a solute and a solvent environment. The solute molecules in all subsystems share the same conformation but have their own solvent environments. Transition between states is implicated by the change of their molar fractions. Simulation of a VMMS system allows efficient calculation of relative free energies of all states, which in turn determine their equilibrium molar fractions. For systems with a large number of state transition sites, an implicit site approximation is introduced to minimize the cost of simulation. A direct application of the VMMS method is for constant pH simulation to study protonation equilibrium. Applying the VMMS method to a heptapeptide of 3 ionizable residues, we calculated the pKas of those residues both with all explicit states and with implicit sites and obtained consistent results. For mouse epidermal growth factor of 9 ionizable groups, our VMMS simulations with implicit sites produced pKas of all 9 ionizable groups and the results agree qualitatively with NMR measurement. This example demonstrates the VMMS method can be applied to systems of a large number of ionizable groups and the computational cost scales linearly with the number of ionizable groups. For one of the most challenging systems in constant pH calculation, SNase Δ+PHS/V66K, our VMMS simulation shows that it is the state-dependent water penetration that causes the large deviation in lysine66’s pKa.

Computer simulation plays an important role to understand molecular systems and has been applied to problems of increasing complexity. Multistate equilibrium is a fundamental concept behind the structure and function of biological systems. Due to the limit in computing resources and lack of good alternative methods, computer simulation has been conducted for systems in a single state, sampling from one state to another to infer equilibrium properties. This sequential approach has been successful in many cases such as protonation equilibrium with implicit solvation model. However, state transition is difficult when explicit solvent is used for more accurate solvation description. Many efforts have been dedicated to overcome this difficulty. Analogous to real multistate systems, we proposed a virtual mixture of multiple states (VMMS) to directly simulate the equilibrium. State transitions are replaced by changes in state molar fractions. Mimicking a test tube environment, all states are simulated in parallel to equilibrate with each other. Application to constant pH simulation in explicit water demonstrates the capability of this method. It is expected that the VMMS method will find more applications in biological problems related to the equilibrium of competing states.