The Relation between Approximation in Distribution and Shadowing in Molecular Dynamics

Molecular dynamics refers to the computer simulation of a material at the atomic level. An open problem in numerical analysis is to explain the apparent reliability of molecular dynamics simulations. The difficulty is that individual trajectories computed in molecular dynamics are accurate for only short time intervals, whereas apparently reliable information can be extracted from very long-time simulations. It has been conjectured that long molecular dynamics trajectories have low-dimensional statistical features that accurately approximate those of the original system. Another conjecture is that numerical trajectories satisfy the shadowing property: that they are close over long time intervals to exact trajectories but with different initial conditions. We prove that these two views are actually equivalent to each other, after we suitably modify the concept of shadowing. A key ingredient of our result is a general theorem that allows us to take random elements of a metric space that are close in distribution and embed them in the same probability space so that they are close in a strong sense. This result is similar to the Strassen-Dudley Theorem except that a mapping is provided between the two random elements. Our results on shadowing are motivated by molecular dynamics but apply to the approximation of any dynamical system when initial conditions are selected according to a probability measure.Comment: 21 pages, final version accepted in SIAM Dyn Sy

Long Time Scale Ensemble Methods in Molecular Dynamics: Ligand–Protein Interactions and Allostery in SARS-CoV-2 Targets

We subject a series of five protein-ligand systems which contain important SARS-CoV-2 targets, 3-chymotrypsin-like protease (3CLPro), papain-like protease, and adenosine ribose phosphatase, to long time scale and adaptive sampling molecular dynamics simulations. By performing ensembles of ten or twelve 10 μs simulations for each system, we accurately and reproducibly determine ligand binding sites, both crystallographically resolved and otherwise, thereby discovering binding sites that can be exploited for drug discovery. We also report robust, ensemble-based observation of conformational changes that occur at the main binding site of 3CLPro due to the presence of another ligand at an allosteric binding site explaining the underlying cascade of events responsible for its inhibitory effect. Using our simulations, we have discovered a novel allosteric mechanism of inhibition for a ligand known to bind only at the substrate binding site. Due to the chaotic nature of molecular dynamics trajectories, regardless of their temporal duration individual trajectories do not allow for accurate or reproducible elucidation of macroscopic expectation values. Unprecedentedly at this time scale, we compare the statistical distribution of protein-ligand contact frequencies for these ten/twelve 10 μs trajectories and find that over 90% of trajectories have significantly different contact frequency distributions. Furthermore, using a direct binding free energy calculation protocol, we determine the ligand binding free energies for each of the identified sites using long time scale simulations. The free energies differ by 0.77 to 7.26 kcal/mol across individual trajectories depending on the binding site and the system. We show that, although this is the standard way such quantities are currently reported at long time scale, individual simulations do not yield reliable free energies. Ensembles of independent trajectories are necessary to overcome the aleatoric uncertainty in order to obtain statistically meaningful and reproducible results. Finally, we compare the application of different free energy methods to these systems and discuss their advantages and disadvantages. Our findings here are generally applicable to all molecular dynamics based applications and not confined to the free energy methods used in this study

Long Time Scale Ensemble Methods in Molecular Dynamics: Ligand-Protein Interactions and Allostery in SARS-CoV-2 Targets

We subject a series of five protein–ligand systems which contain important SARS-CoV-2 targets, 3-chymotrypsin-like protease (3CLPro), papain-like protease, and adenosine ribose phosphatase, to long time scale and adaptive sampling molecular dynamics simulations. By performing ensembles of ten or twelve 10 μs simulations for each system, we accurately and reproducibly determine ligand binding sites, both crystallographically resolved and otherwise, thereby discovering binding sites that can be exploited for drug discovery. We also report robust, ensemble-based observation of conformational changes that occur at the main binding site of 3CLPro due to the presence of another ligand at an allosteric binding site explaining the underlying cascade of events responsible for its inhibitory effect. Using our simulations, we have discovered a novel allosteric mechanism of inhibition for a ligand known to bind only at the substrate binding site. Due to the chaotic nature of molecular dynamics trajectories, regardless of their temporal duration individual trajectories do not allow for accurate or reproducible elucidation of macroscopic expectation values. Unprecedentedly at this time scale, we compare the statistical distribution of protein–ligand contact frequencies for these ten/twelve 10 μs trajectories and find that over 90% of trajectories have significantly different contact frequency distributions. Furthermore, using a direct binding free energy calculation protocol, we determine the ligand binding free energies for each of the identified sites using long time scale simulations. The free energies differ by 0.77 to 7.26 kcal/mol across individual trajectories depending on the binding site and the system. We show that, although this is the standard way such quantities are currently reported at long time scale, individual simulations do not yield reliable free energies. Ensembles of independent trajectories are necessary to overcome the aleatoric uncertainty in order to obtain statistically meaningful and reproducible results. Finally, we compare the application of different free energy methods to these systems and discuss their advantages and disadvantages. Our findings here are generally applicable to all molecular dynamics based applications and not confined to the free energy methods used in this study

Computational Modeling of (De)-Solvation Effects and Protein Flexibility in Protein-Ligand Binding using Molecular Dynamics Simulations

Water is a crucial participant in virtually all cellular functions. Evidently, water molecules in the binding site contribute significantly to the strength of intermolecular interactions in the aqueous phase by mediating protein-ligand interactions, solvating and de-solvating both ligand and protein upon protein-ligand dissociation and association. Recently many published studies use water distributions in the binding site to retrospectively explain and rationalize unexpected trends in structure-activity relationships (SAR). However, traditional approaches cannot quantitatively predict the thermodynamic properties of water molecules in the binding sites and its associated contribution to the binding free energy of a ligand. We have developed and validated a computational method named WATsite to exploit high-resolution solvation maps and thermodynamic profiles to elucidate the water molecules’ potential contribution to protein-ligand and protein-protein binding. We have also demonstrated the utility of the computational method WATsite to help direct medicinal chemistry efforts by using explicit water de-solvation. In addition, protein conformational change is typically involved in the ligand-binding process which may completely change the position and thermodynamic properties of the water molecules in the binding site before or upon ligand binding. We have shown the interplay between protein flexibility and solvent reorganization, and we provide a quantitative estimation of the influence of protein flexibility on desolvation free energy and, therefore, protein-ligand binding. Different ligands binding to the same target protein can induce different conformational adaptations. In order to apply WATsite to an ensemble of different protein conformations, a more efficient implementation of WATsite based on GPU-acceleration and system truncation has been developed. Lastly, by extending the simulation protocol from pure water to mixed water-organic probes simulations, accurate modeling of halogen atom-protein interactions has been achieved

An investigation of molecular dynamics for simple liquids

This thesis contains work expanding the theoretical understanding of molecular dynamics used to aid the study of simple liquids. It does so by focusing on investigating forces, which govern the dynamics of manybody systems. We loosely address three questions: How can we categorise force distributions? What can we gauge from force data? When do forces obey Newton’s third law? The first of these questions is addressed using statistical mechanics to derive standardised moments of the force distribution for a simple LennardJones liquid in both 1d and 3d with the aid of molecular dynamics. To answer the second question, we introduce the notions of force spaces and configurations spaces, and look at equivalence of these. We begin the investigation using the harmonic potential, and develop homotopy continuation methods for non-linear forces like Lennard-Jones. Convergent behaviour and limitations are explored for many-body systems, and a general two-body direct inversion is developed and implemented. The final question is entrenched in classical potential theory, and approached through work focusing on understanding the functional dependence of the interatomic potential. We develop theorems and provide corresponding constructive proofs concluding that potentials which obey certain symmetries can be described by distances, as opposed to positions. This enables us to understand when forces display reciprocity

Optimisation of Non-Canonical Hamiltonian Monte Carlo

This thesis contributes novel developments in the optimisation and theoretical basis of non-canonical Hamiltonian Monte Carlo (NCHMC). NCHMC is a generalisation of Hamiltonian Monte Carlo (HMC) which has gained some interest in recent years for its potential for improved mixing over HMC. Despite the large number of adjustable parameters in the skew-symmetric matrix that govern the dynamics of NCHMC, prior work has not focused on the optimisation of the sampler over these parameters. Additionally, prior work in optimising the numerical integration required by the sampler, by replacing implicit integration with explicit integration, has introduced a source of bias. This work provides methods to address these issues. It also provides a novel derivation for the form the dynamics must adhere to using the global balance condition and demonstrates how NCHMC can be equivalently described by the irreversible jump (I-Jump) sampler. A method to optimise the skew-symmetric matrix via gradient descent is proposed, which is straightforward to implement and capable of producing significant reductions in the chain autocorrelations. For ill-conditioned unimodal densities in particular, the improvement in effective sample size (ESS) is very large. However, this is not necessarily the case for multimodal distributions, with the small improvements in ESS being outweighed by the additional computation time required. A method for explicit, volume-preserving numerical integration of non-canonical Hamiltonian dynamics is presented, on which a computationally fast procedure for generating unbiased samples with NCHMC is built. Unfortunately, experiments demonstrate this sampling method generally performs worse than when using the implicit generalised Störmer-Verlet scheme, as the additional constraints of the sampler force either a smaller step size in the discretised dynamics, or a larger number of rejected samples in the Metropolis-Hastings steps. Whilst not outperforming implicit integration, the results from the method still provide valuable insight into the cost of using an explicit integrator with bias correction