Automated sampling assessment for molecular simulations using the effective sample size

Platinum (II) complex of isopentyl glycine ligand: DNA binding, molecular dynamic, and anticancer activity against breast cancer

In this paper, we performed thorough experimental and theoretical calculations to examine the interaction between Pt derivative, as an anticancer, and ct-DNA. The mode of DNA binding with [Pt(NH3)2(Isopentylgly)]NO3, where Isopentylgly is Isopentyl glycine, was evaluated by various spectroscopic methods, docking, and molecular dynamics simulation studies. UV-Vis and fluorescence spectroscopic titration results and CD spectra of DNA-drug showed this interaction is via groove binding. Also, thermal stability studies or DNA melting temperature changes (ΔTm), as well as the quenching emissions monitoring proved it. Also, the thermodynamic parameter and binding constant displayed that complex-DNA formation is a spontaneous process, and H-binding and also groove binding were found to be the main forces. Theoretical studies stated [Pt(NH3)2(Isopentylgly)]NO3-DNA formation occurs on C-G center on DNA, along with rising DNA-compound stability. IC50 value against the human breast cell line probably is due to the Isopentyl glycine ligand in the structure of the Pt compound, and it was obtained more than cisplatin and less than carboplatin against the MCF7 cell.Communicated by Ramaswamy H. Sarma.

Spectroscopy and molecular dynamic study of the interaction of calf thymus DNA by anticancer Pt complex with butyl glycine ligand

Despite the past few decades since the discovery of anticancer drugs, there is still no definitive treatment for its treatment. Cisplatin is a chemotherapy medication used to treat some cancers. In this research, the DNA binding affinity of Pt complex with butyl glycine ligand was studied by various spectroscopy methods and simulation studies. Fluorescence and UV-Vis spectroscopic data showed groove binding in ct-DNA-[Pt(NH₃)₂(butylgly)]NO₃ complex formation by the spontaneous process. The results were also confirmed by small changes in CD spectra and thermal study (T_m), as well as the quenching emission of [Pt(NH₃)₂(butylgly)]NO₃ complex on DNA. Finally, thermodynamic and binding parameters displayed that hydrophobic forces are the main forces. Based on docking simulation, [Pt(NH₃)₂(butylgly)]NO₃ could bind to DNA and via minor groove binding on C-G center on DNA, formed a stable DNA complex.

Molecular and Biological Investigation of Isolated Marine Fungal Metabolites as Anticancer Agents: A Multi-Target Approach

Cancer is the leading cause of death globally, with an increasing number of cases being annually reported. Nature-derived metabolites have been widely studied for their potential programmed necrosis, cytotoxicity, and anti-proliferation leading to enrichment for the modern medicine, particularly within the last couple of decades. At a more rapid pace, the concept of multi-target agents has evolved from being an innovative approach into a regular drug development procedure for hampering the multi-fashioned pathophysiology and high-resistance nature of cancer cells. With the advent of the Red Sea Penicillium chrysogenum strain S003-isolated indole-based alkaloids, we thoroughly investigated the molecular aspects for three major metabolites: meleagrin (MEL), roquefortine C (ROC), and isoroquefortine C (ISO) against three cancer-associated biological targets Cdc-25A, PTP-1B, and c-Met kinase. The study presented, for the first time, the detailed molecular insights and near-physiological affinity for these marine indole alkaloids against the assign targets through molecular docking-coupled all-atom dynamic simulation analysis. Findings highlighted the superiority of MEL's binding affinity/stability being quite in concordance with the in vitro anticancer activity profile conducted via sulforhodamine B bioassay on different cancerous cell lines reaching down to low micromolar or even nanomolar potencies. The advent of lengthy structural topologies via the metabolites' extended tetracyclic cores and aromatic imidazole arm permitted multi-pocket accommodation addressing the selectivity concerns. Additionally, the presence decorating polar functionalities on the core hydrophobic tetracyclic ring contributed compound's pharmacodynamic preferentiality. Introducing ionizable functionality with more lipophilic characters was highlighted to improve binding affinities which was also in concordance with the conducted drug-likeness/pharmacokinetic profiling for obtaining a balanced pharmacokinetic/dynamic profile. Our study adds to the knowledge regarding drug development and optimization of marine-isolated indole-based alkaloids for future iterative synthesis and pre-clinical investigations as multi-target anticancer agents.

Improving the Efficiency of Variationally Enhanced Sampling with Wavelet-Based Bias Potentials

Collective variable-based enhanced sampling methods are routinely used on systems with metastable states, where high free energy barriers impede the proper sampling of the free energy landscapes when using conventional molecular dynamics simulations. One such method is variationally enhanced sampling (VES), which is based on a variational principle where a bias potential in the space of some chosen slow degrees of freedom, or collective variables, is constructed by minimizing a convex functional. In practice, the bias potential is taken as a linear expansion in some basis function set. So far, primarily basis functions delocalized in the collective variable space, like plane waves, Chebyshev, or Legendre polynomials, have been used. However, there has not been an extensive study of how the convergence behavior is affected by the choice of the basis functions. In particular, it remains an open question if localized basis functions might perform better. In this work, we implement, tune, and validate Daubechies wavelets as basis functions for VES. The wavelets construct orthogonal and localized bases that exhibit an attractive multiresolution property. We evaluate the performance of wavelet and other basis functions on various systems, going from model potentials to the calcium carbonate association process in water. We observe that wavelets exhibit excellent performance and much more robust convergence behavior than all other basis functions, as well as better performance than metadynamics. In particular, using wavelet bases yields far smaller fluctuations of the bias potential within individual runs and smaller differences between independent runs. Based on our overall results, we can recommend wavelets as basis functions for VES.

What Markov State Models Can and Cannot Do: Correlation versus Path-Based Observables in Protein-Folding Models

Markov state models (MSMs) have been widely applied to study the kinetics and pathways of protein conformational dynamics based on statistical analysis of molecular dynamics (MD) simulations. These MSMs coarse-grain both configuration space and time in ways that limit what kinds of observables they can reproduce with high fidelity over different spatial and temporal resolutions. Despite their popularity, there is still limited understanding of which biophysical observables can be computed from these MSMs in a robust and unbiased manner, and which suffer from the space-time coarse-graining intrinsic in the MSM model. Most theoretical arguments and practical validity tests for MSMs rely on long-time equilibrium kinetics, such as the slowest relaxation time scales and experimentally observable time-correlation functions. Here, we perform an extensive assessment of the ability of well-validated protein folding MSMs to accurately reproduce path-based observable such as mean first-passage times (MFPTs) and transition path mechanisms compared to a direct trajectory analysis. We also assess a recently proposed class of history-augmented MSMs (haMSMs) that exploit additional information not accounted for in standard MSMs. We conclude with some practical guidance on the use of MSMs to study various problems in conformational dynamics of biomolecules. In brief, MSMs can accurately reproduce correlation functions slower than the lag time, but path-based observables can only be reliably reproduced if the lifetimes of states exceed the lag time, which is a much stricter requirement. Even in the presence of short-lived states, we find that haMSMs reproduce path-based observables more reliably.

Massive-Scale Binding Free Energy Simulations of HIV Integrase Complexes Using Asynchronous Replica Exchange Framework Implemented on the IBM WCG Distributed Network

To perform massive-scale replica exchange molecular dynamics (REMD) simulations for calculating binding free energies of protein-ligand complexes, we implemented the asynchronous replica exchange (AsyncRE) framework of the binding energy distribution analysis method (BEDAM) in implicit solvent on the IBM World Community Grid (WCG) and optimized the simulation parameters to reduce the overhead and improve the prediction power of the WCG AsyncRE simulations. We also performed the first massive-scale binding free energy calculations using the WCG distributed computing grid and 301 ligands from the SAMPL4 challenge for large-scale binding free energy predictions of HIV-1 integrase complexes. In total there are ∼10000 simulated complexes, ∼1 million replicas, and ∼2000 μs of aggregated MD simulations. Running AsyncRE MD simulations on the WCG requires accepting a trade-off between the number of replicas that can be run (breadth) and the number of full RE cycles that can be completed per replica (depth). As compared with synchronous Replica Exchange (SyncRE) running on tightly coupled clusters like XSEDE, on the WCG many more replicas can be launched simultaneously on heterogeneous distributed hardware, but each full RE cycle requires more overhead. We compared the WCG results with that from AutoDock and more advanced RE simulations including the use of flattening potentials to accelerate sampling of selected degrees of freedom of ligands and/or receptors related to slow dynamics due to high energy barriers. We propose a suitable strategy of RE simulations to refine high throughput docking results which can be matched to corresponding computing resources: from HPC clusters, to small or medium-size distributed campus grids, and finally to massive-scale computing networks including millions of CPUs like the resources available on the WCG.

Best Practices for Quantification of Uncertainty and Sampling Quality in Molecular Simulations [Article v1.0]

The quantitative assessment of uncertainty and sampling quality is essential in molecular simulation. Many systems of interest are highly complex, often at the edge of current computational capabilities. Modelers must therefore analyze and communicate statistical uncertainties so that "consumers" of simulated data understand its significance and limitations. This article covers key analyses appropriate for trajectory data generated by conventional simulation methods such as molecular dynamics and (single Markov chain) Monte Carlo. It also provides guidance for analyzing some 'enhanced' sampling approaches. We do not discuss systematic errors arising, e.g., from inaccuracy in the chosen model or force field.

Tunable Coarse Graining for Monte Carlo Simulations of Proteins via Smoothed Energy Tables: Direct and Exchange Simulations

Many commonly used coarse-grained models for proteins are based on simplified interaction sites and consequently may suffer from significant limitations, such as the inability to properly model protein secondary structure without the addition of restraints. Recent work on a benzene fluid (Lettieri S.; Zuckerman D. M.J. Comput. Chem.2012, 33, 268-275) suggested an alternative strategy of tabulating and smoothing fully atomistic orientation-dependent interactions among rigid molecules or fragments. Here we report our initial efforts to apply this approach to the polar and covalent interactions intrinsic to polypeptides. We divide proteins into nearly rigid fragments, construct distance and orientation-dependent tables of the atomistic interaction energies between those fragments, and apply potential energy smoothing techniques to those tables. The amount of smoothing can be adjusted to give coarse-grained models that range from the underlying atomistic force field all the way to a bead-like coarse-grained model. For a moderate amount of smoothing, the method is able to preserve about 70-90% of the α-helical structure while providing a factor of 3-10 improvement in sampling per unit computation time (depending on how sampling is measured). For a greater amount of smoothing, multiple folding-unfolding transitions of the peptide were observed, along with a factor of 10-100 improvement in sampling per unit computation time, although the time spent in the unfolded state was increased compared with less smoothed simulations. For a β hairpin, secondary structure is also preserved, albeit for a narrower range of the smoothing parameter and, consequently, for a more modest improvement in sampling. We have also applied the new method in a "resolution exchange" setting, in which each replica runs a Monte Carlo simulation with a different degree of smoothing. We obtain exchange rates that compare favorably to our previous efforts at resolution exchange (Lyman E.; Zuckerman D. M.J. Chem. Theory Comput.2006, 2, 656-666).

Sampling Assessment for Molecular Simulations Using Conformational Entropy Calculations

The extent and significance of conformational sampling is a major factor determining the reliability of long-scale molecular simulations of large and flexible biomolecules. Although several methods have been proposed to quantify the effective sample size of molecular simulations by transforming root mean squared distances between pairs of configurations into statistical/probabilistic quantities, there is still no standard technique for measuring the size of sampling. In this work, we study conformational entropy (Sconform) as a purely informational and probabilistic measure of sampling that does not require the adoption of any clustering protocol or distance metric between configurations. In addition Sconform, which is calculated from the probability mass functions associated with discretized dihedral angles, offers other potential advantages for sampling assessment (e.g., global character, thermodynamic significance, relationship with internal degrees of freedom, etc.). The utility of Sconform in sampling assessment is illustrated by carrying out test calculations on configurations produced by two extended molecular dynamics simulations, namely, a 2.0 μs trajectory of a highly flexible 17-residue peptide and the trajectory data set of the 1.0 ms bovine pancreatic trypsin inhibitor simulation provided by the D. E. Shaw research group.

Finding semirigid domains in biomolecules by clustering pair-distance variations

Dynamic variations in the distances between pairs of atoms are used for clustering subdomains of biomolecules. We draw on a well-known target function for clustering and first show mathematically that the assignment of atoms to clusters has to be crisp, not fuzzy, as hitherto assumed. This reduces the computational load of clustering drastically, and we demonstrate results for several biomolecules relevant in immunoinformatics. Results are evaluated regarding the number of clusters, cluster size, cluster stability, and the evolution of clusters over time. Crisp clustering lends itself as an efficient tool to locate semirigid domains in the simulation of biomolecules. Such domains seem crucial for an optimum performance of subsequent statistical analyses, aiming at detecting minute motional patterns related to antigen recognition and signal transduction.

Simultaneous Computation of Dynamical and Equilibrium Information Using a Weighted Ensemble of Trajectories

Equilibrium formally can be represented as an ensemble of uncoupled systems undergoing unbiased dynamics in which detailed balance is maintained. Many nonequilibrium processes can be described by suitable subsets of the equilibrium ensemble. Here, we employ the “weighted ensemble” (WE) simulation protocol [Huber and Kim, Biophys. J.1996, 70, 97–110] to generate equilibrium trajectory ensembles and extract nonequilibrium subsets for computing kinetic quantities. States do not need to be chosen in advance. The procedure formally allows estimation of kinetic rates between arbitrary states chosen after the simulation, along with their equilibrium populations. We also describe a related history-dependent matrix procedure for estimating equilibrium and nonequilibrium observables when phase space has been divided into arbitrary non-Markovian regions, whether in WE or ordinary simulation. In this proof-of-principle study, these methods are successfully applied and validated on two molecular systems: explicitly solvated methane association and the implicitly solvated Ala4 peptide. We comment on challenges remaining in WE calculations.

Molecular dynamics simulations on the Tre1 G protein-coupled receptor: exploring the role of the arginine of the NRY motif in Tre1 structure

BACKGROUND

The arginine of the D/E/NRY motif in Rhodopsin family G protein-coupled receptors (GPCRs) is conserved in 96% of these proteins. In some GPCRs, this arginine in transmembrane 3 can form a salt bridge with an aspartic acid or glutamic acid in transmembrane 6. The Drosophila melanogaster GPCR Trapped in endoderm-1 (Tre1) is required for normal primordial germ cell migration. In a mutant form of the protein, Tre1sctt, eight amino acids RYILIACH are missing, resulting in a severe disruption of primordial germ cell development. The impact of the loss of these amino acids on Tre1 structure is unknown. Since the missing amino acids in Tre1sctt include the arginine that is part of the D/E/NRY motif in Tre1, molecular dynamics simulations were performed to explore the hypothesis that these amino acids are involved in salt bridge formation and help maintain Tre1 structure.

RESULTS

Structural predictions of wild type Tre1 (Tre1+) and Tre1sctt were subjected to over 250 ns of molecular dynamics simulations. The ability of the model systems to form a salt bridge between the arginine of the D/E/NRY motif and an aspartic acid residue in transmembrane 6 was analyzed. The results indicate that a stable salt bridge can form in the Tre1+ systems and a weak salt bridge or no salt bridge, using an alternative arginine, is likely in the Tre1sctt systems.

CONCLUSIONS

The weak salt bridge or lack of a salt bridge in the Tre1sctt systems could be one possible explanation for the disrupted function of Tre1sctt in primordial germ cell migration. These results provide a framework for studying the importance of the arginine of the D/E/NRY motif in the structure and function of other GPCRs that are involved in cell migration, such as CXCR4 in the mouse, zebrafish, and chicken.

Modeling complex biological systems: From solution chemistry to membranes and channels

Complex biological systems are intimately linked to their environment, a very crowded and equally complex solution compartmentalized by fluid membranes. Modeling such systems remains challenging and requires a suitable representation of these solutions and their interfaces. Here, we focus on particle-based modeling at an atomistic level using molecular dynamics (MD) simulations. As an example, we discuss important steps in modeling the solution chemistry of an ion channel of the ligand-gated ion channel receptor family, a major target of many drugs including anesthetics and addiction treatments. The bacterial pentameric ligand-gated ion channel (pLGIC) called GLIC provides clues about the functional importance of solvation, in particular for mechanisms such as permeation and gating. We present some current challenges along with promising novel modeling approaches.

Relaxation estimation of RMSD in molecular dynamics immunosimulations

Molecular dynamics simulations have to be sufficiently long to draw reliable conclusions. However, no method exists to prove that a simulation has converged. We suggest the method of "lagged RMSD-analysis" as a tool to judge if an MD simulation has not yet run long enough. The analysis is based on RMSD values between pairs of configurations separated by variable time intervals Δt. Unless RMSD(Δt) has reached a stationary shape, the simulation has not yet converged.

Tunable, mixed-resolution modeling using library-based Monte Carlo and graphics processing units

Building on our recently introduced library-based Monte Carlo (LBMC) approach, we describe a flexible protocol for mixed coarse-grained (CG)/all-atom (AA) simulation of proteins and ligands. In the present implementation of LBMC, protein side chain configurations are pre-calculated and stored in libraries, while bonded interactions along the backbone are treated explicitly. Because the AA side chain coordinates are maintained at minimal run-time cost, arbitrary sites and interaction terms can be turned on to create mixed-resolution models. For example, an AA region of interest such as a binding site can be coupled to a CG model for the rest of the protein. We have additionally developed a hybrid implementation of the generalized Born/surface area (GBSA) implicit solvent model suitable for mixed-resolution models, which in turn was ported to a graphics processing unit (GPU) for faster calculation. The new software was applied to study two systems: (i) the behavior of spin labels on the B1 domain of protein G (GB1) and (ii) docking of randomly initialized estradiol configurations to the ligand binding domain of the estrogen receptor (ERα). The performance of the GPU version of the code was also benchmarked in a number of additional systems.

Will molecular dynamics simulations of proteins ever reach equilibrium?

We show that conformational entropies calculated for five proteins and protein-ligand complexes with dihedral-distribution histogramming, the von Mises approach, or quasi-harmonic analysis do not converge to any useful precision even if molecular dynamics (MD) simulations of 380-500 ns length are employed (the uncertainty is 12-89 kJ mol(-1)). To explain this, we suggest a simple protein model involving dihedrals with effective barriers forming a uniform distribution and show that for such a model, the entropy increases logarithmically with time until all significantly populated dihedral states have been sampled, in agreement with the simulations (during the simulations, 52-70% of the available dihedral phase space has been visited). This is also confirmed by the analysis of the trajectories of a 1 ms simulation of bovine pancreatic trypsin inhibitor (31 kJ mol(-1) difference in the entropy between the first and second part of the simulation). Strictly speaking, this means that it is practically impossible to equilibrate MD simulations of proteins. We discuss the implications of such a lack of strict equilibration of protein MD simulations and show that ligand-binding free energies estimated with the MM/GBSA method (molecular mechanics with generalised Born and surface-area solvation) vary by 3-15 kJ mol(-1) during a 500 ns simulation (the higher estimate is caused by rare conformational changes), although they involve a questionable but well-converged normal-mode entropy estimate, whereas free energies estimated by free-energy perturbation vary by less than 0.6 kJ mol(-1) for the same simulation.

Equilibrium sampling in biomolecular simulations

Equilibrium sampling of biomolecules remains an unmet challenge after more than 30 years of atomistic simulation. Efforts to enhance sampling capability, which are reviewed here, range from the development of new algorithms to parallelization to novel uses of hardware. Special focus is placed on classifying algorithms--most of which are underpinned by a few key ideas--in order to understand their fundamental strengths and limitations. Although algorithms have proliferated, progress resulting from novel hardware use appears to be more clear-cut than from algorithms alone, due partly to the lack of widely used sampling measures.