Biomolecular Simulations Research Articles

All-Atom Biomolecular Simulation in the Exascale Era.

Exascale supercomputers have opened the door to dynamic simulations, facilitated by AI/ML techniques, that model biomolecular motions over unprecedented length and time scales. This new capability holds the potential to revolutionize our understanding of fundamental biological processes. Here we report on some of the major advances that were discussed at a recent CECAM workshop in Pisa, Italy, on the topic with a primary focus on atomic-level simulations. First, we highlight examples of current large-scale biomolecular simulations and the future possibilities enabled by crossing the exascale threshold. Next, we discuss challenges to be overcome in optimizing the usage of these powerful resources. Finally, we close by listing several grand challenge problems that could be investigated with this new computer architecture.

Log-periodic oscillations as real-time signatures of hierarchical dynamics in proteins.

The time-dependent relaxation of a dynamical system may exhibit a power-law behavior that is superimposed by log-periodic oscillations. D. Sornette [Phys. Rep. 297, 239 (1998)] showed that this behavior can be explained by a discrete scale invariance of the system, which is associated with discrete and equidistant timescales on a logarithmic scale. Examples include such diverse fields as financial crashes, random diffusion, and quantum topological materials. Recent time-resolved experiments and molecular dynamics simulations suggest that discrete scale invariance may also apply to hierarchical dynamics in proteins, where several fast local conformational changes are a prerequisite for a slow global transition to occur. Employing entropy-based timescale analysis and Markov state modeling to a simple one-dimensional hierarchical model and biomolecular simulation data, it is found that hierarchical systems quite generally give rise to logarithmically spaced discrete timescales. By introducing a one-dimensional reaction coordinate that collectively accounts for the hierarchically coupled degrees of freedom, the free energy landscape exhibits a characteristic staircase shape with two metastable end states, which causes the log-periodic time evolution of the system. The period of the log-oscillations reflects the effective roughness of the energy landscape and can, in simple cases, be interpreted in terms of the barriers of the staircase landscape.

Widespread alteration of protein autoinhibition in human cancers

Autoinhibition is a prevalent allosteric regulatory mechanism in signaling proteins. Reduced autoinhibition underlies the tumorigenic effect of some known cancer drivers, but whether autoinhibition is altered generally in cancer remains elusive. Here, we demonstrate that cancer-associated missense mutations, in-frame insertions/deletions, and fusion breakpoints are enriched within inhibitory allosteric switches (IASs) across all cancer types. Selection for IASs that are recurrently mutated in cancers identifies established and unknown cancer drivers. Recurrent missense mutations in IASs of these drivers are associated with distinct, cancer-specific changes in molecular signaling. For the specific case of PPP3CA, the catalytic subunit of calcineurin, we provide insights into the molecular mechanisms of altered autoinhibition by cancer mutations using biomolecular simulations, and demonstrate that such mutations are associated with transcriptome changes consistent with increased calcineurin signaling. Our integrative study shows that autoinhibition-modulating genetic alterations are positively selected for by cancer cells.

A Machine Learning Force Field for Bio-Macromolecular Modeling Based on Quantum Chemistry-Calculated Interaction Energy Datasets.

Accurate energy data from noncovalent interactions are essential for constructing force fields for molecular dynamics simulations of bio-macromolecular systems. There are two important practical issues in the construction of a reliable force field with the hope of balancing the desired chemical accuracy and working efficiency. One is to determine a suitable quantum chemistry level of theory for calculating interaction energies. The other is to use a suitable continuous energy function to model the quantum chemical energy data. For the first issue, we have recently calculated the intermolecular interaction energies using the SAPT0 level of theory, and we have systematically organized these energies into the ab initio SOFG-31 (homodimer) and SOFG-31-heterodimer datasets. In this work, we re-calculate these interaction energies by using the more advanced SAPT2 level of theory with a wider series of basis sets. Our purpose is to determine the SAPT level of theory proper for interaction energies with respect to the CCSD(T)/CBS benchmark chemical accuracy. Next, to utilize these energy datasets, we employ one of the well-developed machine learning techniques, called the CLIFF scheme, to construct a general-purpose force field for biomolecular dynamics simulations. Here we use the SOFG-31 dataset and the SOFG-31-heterodimer dataset as the training and test sets, respectively. Our results demonstrate that using the CLIFF scheme can reproduce a diverse range of dimeric interaction energy patterns with only a small training set. The overall errors for each SAPT energy component, as well as the SAPT total energy, are all well below the desired chemical accuracy of ~1 kcal/mol.

Biophysical Modelling for Insight into Oxygen Diffusion, Distribution, and Measurement.

Building on an extensive history of physiological and systems-oriented modelling, my group and others have recently used molecular simulation studies to understand oxygen (O2) transport and localisation. Molecular simulations enable biophysical insight into processes difficult to study with experiments alone and are sometimes described as a "computational microscope." Our work has emphasised lipid membrane contributions to oxygen diffusion and uptake, suggesting that lipid-based pathways along membranes and lipid deposits are likely to accelerate diffusive transport through cells and tissues. Moreover, the lipid and fluid fractions of the tissue are expected to be primary determinants of the local oxygen partial pressure (pO2) as well as the oxygen permeability. Measurements using molecular probes can be influenced by the local molecular environment, due to differential solubility of both the probe and the oxygen molecules in various components of the cell's complex solvent system. The biomolecular simulation work complements experimental studies, which enable evaluation of the models' accuracy and their applicability to real biological systems. Further work is needed to assess fully the possible influence of nanoscale crowders and obstacles (especially protein molecules) on tissue-level diffusive transport of oxygen. Likewise, water-rich carbohydrate layers, such as the glycocalyx, should be evaluated as potential barriers to oxygen transport. Insights gained through biophysical modelling studies could be broadly relevant to clinical phenomena affected by tissue oxygenation, such as tumour radiotherapy, ischaemia, neuropathy, and wound healing.

Machine-learned molecular mechanics force fields from large-scale quantum chemical data.

The development of reliable and extensible molecular mechanics (MM) force fields-fast, empirical models characterizing the potential energy surface of molecular systems-is indispensable for biomolecular simulation and computer-aided drug design. Here, we introduce a generalized and extensible machine-learned MM force field, espaloma-0.3, and an end-to-end differentiable framework using graph neural networks to overcome the limitations of traditional rule-based methods. Trained in a single GPU-day to fit a large and diverse quantum chemical dataset of over 1.1 M energy and force calculations, espaloma-0.3 reproduces quantum chemical energetic properties of chemical domains highly relevant to drug discovery, including small molecules, peptides, and nucleic acids. Moreover, this force field maintains the quantum chemical energy-minimized geometries of small molecules and preserves the condensed phase properties of peptides and folded proteins, self-consistently parametrizing proteins and ligands to produce stable simulations leading to highly accurate predictions of binding free energies. This methodology demonstrates significant promise as a path forward for systematically building more accurate force fields that are easily extensible to new chemical domains of interest.

Development and Comprehensive Benchmark of a High-Quality AMBER-Consistent Small Molecule Force Field with Broad Chemical Space Coverage for Molecular Modeling and Free Energy Calculation.

Biomolecular simulations have become an essential tool in contemporary drug discovery, and molecular mechanics force fields (FFs) constitute its cornerstone. Developing a high quality and broad coverage general FF is a significant undertaking that requires substantial expert knowledge and computing resources, which is beyond the scope of general practitioners. Existing FFs originate from only a limited number of groups and organizations, and they either suffer from limited numbers of training sets, lower than desired quality because of oversimplified representations, or are costly for the molecular modeling community to access. To address these issues, in this work, we developed an AMBER-consistent small molecule FF with extensive chemical space coverage, and we provide Open Access parameters for the entire modeling community. To validate our FF, we carried out benchmarks of quantum mechanics (QM)/molecular mechanics conformer comparison and free energy perturbation calculations on several benchmark data sets. Our FF achieves a higher level of performance at reproducing QM energies and geometries than two popular open-source FFs, OpenFF2 and GAFF2. In relative binding free energy calculations for 31 protein-ligand data sets, comprising 1079 pairs of ligands, the new FF achieves an overall root-mean-square error of 1.19 kcal/mol for ΔΔG and 0.92 kcal/mol for ΔG on a subset of 463 ligands without bespoke fitting to the data sets. The results are on par with those of the leading commercial series of OPLS FFs.

Biomolecular simulation to elucidate small-molecule modulation of mechanosensor protein.

Biomolecular simulation to elucidate small-molecule modulation of mechanosensor protein.

Identification of novel STAT3 inhibitors for liver fibrosis, using pharmacophore-based virtual screening, molecular docking, and biomolecular dynamics simulations

The signal transducer and activator of transcription 3 (STAT3) plays a fundamental role in the growth and regulation of cellular life. Activation and over-expression of STAT3 have been implicated in many cancers including solid blood tumors and other diseases such as liver fibrosis and rheumatoid arthritis. Therefore, STAT3 inhibitors are be coming a growing and interesting area of pharmacological research. Consequently, the aim of this study is to design novel inhibitors of STAT3-SH3 computationally for the reduction of liver fibrosis. Herein, we performed Pharmacophore-based virtual screening of databases including more than 19,481 commercially available compounds and in-house compounds. The hits obtained from virtual screening were further docked with the STAT3 receptor. The hits were further ranked on the basis of docking score and binding interaction with the active site of STAT3. ADMET properties of the screened compounds were calculated and filtered based on drug-likeness criteria. Finally, the top five drug-like hit compounds were selected and subjected to molecular dynamic simulation. The stability of each drug-like hit in complex with STAT3 was determined by computing their RMSD, RMSF, Rg, and DCCM analyses. Among all the compounds Sa32 revealed a good docking score, interactions, and stability during the entire simulation procedure. As compared to the Reference compound, the drug-like hit compound Sa32 showed good docking scores, interaction, stability, and binding energy. Therefore, we identified Sa32 as the best small molecule potent inhibitor for STAT3 that will be helpful in the future for the treatment of liver fibrosis.

Pipelines for automating compliance-based elimination and extension (PACE2): a systematic framework for high-throughput biomolecular materials simulation workflows

The formation of biomolecular materials via dynamical interfacial processes, such as self-assembly and fusion, for diverse compositions and external conditions can be efficiently probed using ensemble Molecular Dynamics (MD). However, this approach requires many simulations when investigating a large composition phase space. In addition, there is difficulty in predicting whether each simulation will yield biomolecular materials with the desired properties or outcomes and how long each simulation will run. These difficulties can be overcome by rules-based management systems, including intermittent inspection, variable sampling, and premature termination or extension of the individual MD simulations. Automating such a management system can significantly improve runtime efficiency and reduce the burden of organizing large ensembles of MD simulations. To this end, a computational framework, the Pipelines for Automating Compliance-based Elimination and Extension (PACE2), is proposed for high-throughput ensemble biomolecular materials simulations. The PACE2 framework encompasses Candidate pipelines, where each pipeline includes temporally separated simulation and analysis tasks. When a MD simulation is completed, an analysis task is triggered, which evaluates the MD trajectory for compliance. Compliant simulations are extended to the next MD phase with a suitable sample rate to allow additional, detailed analysis. Non-compliant simulations are eliminated, and their computational resources are reallocated or released. The framework is designed to run on local desktop computers and high-performance computing resources. Preliminary scientific results enabled by the use of PACE2 framework are presented, which demonstrate its potential and validates its function. In the future, the framework will be extended to address generalized workflows and investigate composition-structure-property relations for other classes of materials.

NNP/MM: Accelerating Molecular Dynamics Simulations with Machine Learning Potentials and Molecular Mechanics.

Machine learning potentials have emerged as a means to enhance the accuracy of biomolecular simulations. However, their application is constrained by the significant computational cost arising from the vast number of parameters compared with traditional molecular mechanics. To tackle this issue, we introduce an optimized implementation of the hybrid method (NNP/MM), which combines a neural network potential (NNP) and molecular mechanics (MM). This approach models a portion of the system, such as a small molecule, using NNP while employing MM for the remaining system to boost efficiency. By conducting molecular dynamics (MD) simulations on various protein-ligand complexes and metadynamics (MTD) simulations on a ligand, we showcase the capabilities of our implementation of NNP/MM. It has enabled us to increase the simulation speed by ∼5 times and achieve a combined sampling of 1 μs for each complex, marking the longest simulations ever reported for this class of simulations.

Replica Exchange with Hybrid Tempering Efficiently Samples PGLa Peptide Binding to Anionic Bilayer.

We evaluated the utility of a variant of the replica exchange method, a replica exchange with hybrid tempering (REHT), for all-atom explicit water biomolecular simulations and compared it with a more traditional replica exchange with the solute tempering (REST) algorithm. As a test system, we selected a 21-mer antimicrobial peptide PGLa binding to an anionic DMPC/DMPG lipid bilayer. Application of REHT revealed the following binding mechanism. Due to the strong hydrophobic moment, the bound PGLa adopts an extensive helical structure. The binding free energy landscape identifies two major bound states, a metastable surface bound state and a dominant inserted state. In both states, positively charged PGLa amino acids maintain electrostatic interactions with anionic phosphate groups by rotating the PGLa helix around its axis. PGLa binding causes an influx of anionic DMPG and an efflux of zwitterionic DMPC lipids from the peptide proximity. PGLa thins the bilayer and disorders the adjacent fatty acid tails. Deep invasion of water wires into the bilayer hydrophobic core is detected in the inserted peptide state. The analysis of charge density distributions indicated that peptide positive charges are nearly compensated for by lipid negative charges and water dipole ordering, whereas ions play no role in peptide binding. Thus, electrostatic interactions are the key energetic factor in binding cationic PGLa to an anionic DMPC/DMPG bilayer. Comparison of REHT and REST shows that due to exclusion of lipids from tempered partition, REST lags behind REHT in peptide equilibration, particularly, with respect to peptide insertion and helix acquisition. As a result, REST struggles to provide accurate details of PGLa binding, although it still qualitatively maps the bimodal binding mechanism. Importantly, REHT not only equilibrates PGLa in the bilayer faster than REST, but also with less computational effort. We conclude that REHT is a preferable choice for studying interfacial biomolecular systems.

Binding kinetics of ten small-molecule drug candidates on SARS-CoV-2 3CLpro revealed by biomolecular simulations

3CL protease (3CLpro) is the main protease (Mpro) found in severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which cuts the coronavirus polyprotein at eleven conserved sites and is essential for the virus replication. Therefore, 3CLpro has been widely used as a promising drug target. Many small-molecule drug candidates targeting 3CLpro have been proposed to inhibit the virus replication. In this work, we aim to reveal detailed interactions between ten small candidate molecules with extensive attention and 3CLpro using molecular docking and molecular dynamics simulations. First, we identified the possible binding sites of these candidate molecules on 3CLpro via molecular docking. Then, a series of 100 ns all-atom molecular dynamics simulations of strongest binding modes were performed to further evaluate the dynamical interactions between the molecules and 3CLpro in detail. Last, the binding free energy of these molecules on 3CLpro was calculated using MM/PBSA calculation, where the contribution of key amino acids was highlighted. The binding kinetics revealed in this work may provide useful insights into the action mechanism and applicability of these small-molecule drug candidates.

Toward a Benchmark for Markov State Models: The Folding of HP35.

Adopting a 300 μs long MD trajectory of the folding of villin headpiece (HP35) by D. E. Shaw Research, we recently constructed a Markov state model (MSM) based on inter-residue contacts. The model reproduces the folding time and predicts that the native basin and unfolded region consist of metastable substates that are structurally well-characterized. Recognizing the need to establish well-defined benchmark problems, we study to what extent and in what sense this MSM can be employed as a reference model. Hence, we test the robustness of the MSM by comparing it to models that use alternative combinations of features, dimensionality reduction methods, and clustering schemes. The study suggests some main characteristics of the folding of HP35 that should be reproduced by other competitive models. Moreover, the discussion reveals which parts of the MSM workflow matter most for the considered problem and illustrates the promises and pitfalls of state-based models for the interpretation of biomolecular simulations.

MDSubSampler: a posteriori sampling of important protein conformations from biomolecular simulations.

Molecular dynamics (MD) simulations have become routine tools for the study of protein dynamics and function. Thanks to faster GPU-based algorithms, atomistic and coarse-grained simulations are being used to explore biological functions over the microsecond timescale, yielding terabytes of data spanning multiple trajectories, thereby extracting relevant protein conformations without losing important information is often challenging. We present MDSubSampler, a Python library and toolkit for a posteriori subsampling of data from multiple trajectories. This toolkit provides access to uniform, random, stratified, weighted sampling, and bootstrapping sampling methods. Sampling can be performed under the constraint of preserving the original distribution of relevant geometrical properties. Possible applications include simulations post-processing, noise reduction, and structures selection for ensemble docking. MDSubSampler is freely available at https://github.com/alepandini/MDSubSampler, along with guidance on installation and tutorials on how it can be used.

Optimized Quantum Drude Oscillators for Atomic and Molecular Response Properties.

The quantum Drude oscillator (QDO) is an efficient yet accurate coarse-grained approach that has been widely used to model electronic and optical response properties of atoms and molecules as well as polarization and dispersion interactions between them. Three effective parameters (frequency, mass, and charge) fully characterize the QDO Hamiltonian and are adjusted to reproduce response properties. However, the soaring success of coupled QDOs for many-atom systems remains fundamentally unexplained, and the optimal mapping between atoms/molecules and oscillators has not been established. Here we present an optimized parametrization (OQDO) where the parameters are fixed by using only dipolar properties. For the periodic table of elements as well as small molecules, our model accurately reproduces atomic (spatial) polarization potentials and multipolar dispersion coefficients, elucidating the high promise of the presented model in the development of next-generation quantum-mechanical force fields for (bio)molecular simulations.

Free Energy Calculations and Conformational Analysis of Dibenzo-30-crown-10 with Sm2+, Eu2+, and Three Halide Salts in THF Using the AMOEBA Force Field.

Crown ether complexes have been tailored for use in industrial separations of lanthanides (Ln) as a part of rare earth mining and refining. Dibenzo-30-crown-10 (DB30C10) is one of the most efficient complexants for the separation of rare earth mixtures based on the cation size. To understand the origin of this complexation, molecular dynamics (MD) simulations of DB30C10 have been performed using different combinations of divalent Sm and Eu and three halide salts Cl-, Br-, and I- in tetrahydrofuran (THF) solvent. DB30C10 was parameterized here for the polarizable atomic multipole optimized energetics for biomolecular simulation (AMOEBA) force field, and the existing parameters of THF, Sm2+, and Eu2+ were employed from our previous efforts. The large conformational fluctuations present in the DB30C10 systems were found to be dependent both on the identity of the lanthanide and halide complexes. For Cl- and Br- systems, there were no observed conformational changes at 200 ns, while in I- systems, there were two conformational changes with Sm2+ and one with Eu2+ within that same timeframe. In SmI2-DB30C10, there were three stages of conformational changes. In the first stage, the molecule is unfolded, in the second stage, the molecule is partly folded, and finally, in the third stage, the molecule is completely folded. Lastly, the Gibbs binding free energies of DB30C10 with SmBr2 and EuBr2 have been computed, which resulted in nearly identical ΔGcomp values for each lanthanide with Sm2+ being slightly more favorable. Considering the folding mechanism of the SmI2 system with DB30C10, the Gibbs binding free energies of DB30C10 and dicyclohexano-18-crown-6 (DCH18C6) with SmI2 were calculated separately and compared to probe their complexation affinities, in which the former was found to be more favorable.

Simulating water adsorption in metal–organic frameworks with open metal sites using the 12-6-4 Lennard–Jones potential

ABSTRACT Metal–organic frameworks (MOFs) with coordinatively unsaturated open metal sites (OMSs) are promising sorbent materials in atmospheric water harvesting (AWH) systems at low relative humidity (RH) due to their strong interactions with water molecules. However, accurate computational modelling of water adsorption in those materials are challenging, as standard force fields (FFs) based on the 12-6 Lennard–Jones (L–J) potential cannot properly describe the water–OMS interactions. In biomolecular simulations, the 12-6-4 L–J potential has been successfully used to model metal ion–water interactions in which the extra 1/r 4-term is for charge-induced dipole electrostatics. In this work, we adopted this strategy and used the 12-6-4 L–J potential to model water–OMS interactions in MOFs for low RH AWH applications. Without any modifications, the parameters used for aqueous metal ions were able to greatly improve the accuracy in predicting water adsorption isotherms in MOF-74, which highlights the simplicity and transferability of this method.

Quantum healthcare computing using precision based granular approach

Previously, doctors interpreted diseases and their outcomes according to their experience in diagnosis. However, with the rapid increase in technology and population, the task of examining the patient becomes cumbersome and sometimes human efforts produce inconsistent results. Several research is being done for healthcare in terms of improving visualization and accuracy by using machine learning models. The current research targets to explore quantum computing as a different way of processing information compared to classical computer systems such as the use of quantum bits (qubits) along with superposition and entanglement for extending the computation capabilities at an unprecedented level of thinking in the healthcare domain. Quantum computing systems provide exponential benefits in terms of high-speed processing, faster and easier diagnostic assistance, unimaginable reduction in processing throughput, and many more. An extensive comparative analysis of existing approaches has been made which benchmarks the need for quantum healthcare computing. The objective of this work is to interpret whether Quantum computers prove to be more trusted when it comes to patient diagnosis, and faster analysis leading to cost optimization. In order to accelerate patient diagnosis, different approaches have been presented. The authors have proposed a precision-based granular approach for patient diagnosis that incorporates diagnosing the disease with enhanced precision and granularity. It involves reporting symptoms by the patient, encountering by healthcare expert on multiple factors, precise examination, granular health status (understanding past and present medical history), followed by a precise intervention by understanding biomolecular simulations. The algorithm has been presented to describe the flow process for patient diagnosis modeling using quantum computing. It involves qubits initialization, pairing the values, assigning probabilistic values, cross-validation, and quantum circuit formation. Precision-based granular approach has been implemented for a scenario (consisting of medical parameters such as oxygen and heart rate level, with the functionality of diagnosing oxygen level and heart range which lies as either normal or not normal (high/low)). Precision-based granular approach deals specifically with the individual ‘biomolecular simulation by understanding variations in the individual body whereas the umbrella-based approach does not deal with specifically to individual mechanisms. Granular level of encounter is not possible in umbrella-based treatment. Python Jupyter notebook and IBM Composer tool is used for the implementation of results. Bloch sphere and computational state graph are obtained as an output for better visualization and understanding. Falcon r5.11H processor is used with the version of 1.0.24 of IBM Composer to simulate the experiment. The methodology using precision based granular approach provides timely encounter of disease along with umbrella diagnosis and precise treatment. The time is taken and frequency of qubits have been presented with promising results. The diagnosis process and optimizing cost efficiency can aid in an early detection of the disease.

Comparison of force fields to study the zinc-finger containing protein NPL4, a target for disulfiram in cancer therapy

Molecular dynamics (MD) simulations are a powerful approach to studying the structure and dynamics of proteins related to health and disease. Advances in the MD field allow modeling proteins with high accuracy. However, modeling metal ions and their interactions with proteins is still challenging. NPL4 is a zinc-binding protein and works as a cofactor for p97 to regulate protein homeostasis. NPL4 is of biomedical importance and has been proposed as the target of disulfiram, a drug recently repurposed for cancer treatment. Experimental studies proposed that the disulfiram metabolites, bis-(diethyldithiocarbamate)‑copper and cupric ions, induce NPL4 misfolding and aggregation. However, the molecular details of their interactions with NPL4 and consequent structural effects are still elusive. Here, biomolecular simulations can help to shed light on the related structural details. To apply MD simulations to NPL4 and its interaction with copper the first important step is identifying a suitable force field to describe the protein in its zinc-bound states. We examined different sets of non-bonded parameters because we want to study the misfolding mechanism and cannot rule out that the zinc may detach from the protein during the process and copper replaces it. We investigated the force-field ability to model the coordination geometry of the metal ions by comparing the results from MD simulations with optimized geometries from quantum mechanics (QM) calculations using model systems of NPL4. Furthermore, we investigated the performance of a force field including bonded parameters to treat copper ions in NPL4 that we obtained based on QM calculations.