key: cord-0486709-n0nwy9b9 authors: Wan, Shunzhou; Sinclair, Robert C.; Coveney, Peter V. title: Uncertainty Quantification in Classical Molecular Dynamics date: 2020-06-12 journal: nan DOI: nan sha: d6aff5450045f17e17ea9f58a03b83cb30a44c3b doc_id: 486709 cord_uid: n0nwy9b9 Molecular dynamics simulation is now a widespread approach for understanding complex systems on the atomistic scale. It finds applications from physics and chemistry to engineering, life and medical science. In the last decade, the approach has begun to advance from being a computer-based means of rationalising experimental observations, to producing apparently credible predictions for a number of real-world applications within industrial sectors such as advanced materials and drug discovery. However, key aspects concerning the reproducibility of the method have not kept pace with the speed of its uptake in the scientific community. Here, we present a discussion of uncertainty quantification for molecular dynamics simulation designed to endow the method with better error estimates that will enable the method to be used to report actionable results. The approach adopted is a standard one in the field of uncertainty quantification, namely using ensemble methods, in which a sufficiently large number of replicas are run concurrently, from which reliable statistics can be extracted. Indeed, because molecular dynamics is intrinsically chaotic, the need to use ensemble methods is fundamental and holds regardless of the duration of the simulations performed. We discuss the approach and illustrate it in a range of applications from materials science to ligand-protein binding free energy estimation. which is to be relied upon for taking actionable decisions and thus to become valuable in diverse applications, including inter alia industrial and clinical contexts. For that, we need uncertainty quantification (UQ), verification and validation (V&V), or VVUQ. But while careful control of uncertainty is the mainstay of weather forecasting, along with many branches of engineering and applied mathematics, it is rather rarely performed in disciplines such as physics and chemistry where much time is spent investigating matter at shorter length and time scales than the macroscopic ones of direct concern in many real world situations. A major use of MD simulation is to predict the binding affinity of a lead compound or drug with a protein target, of major relevance in drug discovery and personalised medicine. That target may be respectively either a generic protein or a sequence specific variant, reflecting the fact that individuals respond differently to a given drug based on their genetic make-up. The binding affinity, also known as the free energy of binding, is the single most important initial indicator of drug potency, and the most challenging to predict. Reproducibility is an intrinsic feature of the scientific method, whether experimental or computational. A method cannot be regarded as scientific if it does not yield the same result regardless of who performs it. what boundary conditions are to be used, and so on. As Michael Levitt has stated: "the art is to find an approximation simple enough to be computable, but not so simple that you lose the useful detail"5. In principle, a fine-grained model should produce more accurate predictions than a coarse-grained one, although in practice it is not always the case. In drug discovery approaches, for example, a ligand-protein model with explicit water molecules is usually better than one with implicit water. These choices all affect the outcome of a simulation, usually in a deterministic way. Biases in the interaction parameters chosen to represent the system can significantly influence the results; for example, different protein force fields favour different secondary structure types6, populating either helical or sheet-like structures within independent simulations. When the cause of such systematic errors can be identified, it can be reduced or even eliminated, as shown for example in recent simulations with state-of-the-art force fields7. The implementation of the model in an appropriate MD engine and the calibration of the engine can also influence the results. The thermodynamic conditions, such as constant volume or pressure in a closed system, must be specified. A few operational parameters need to be fine-tuned, including those for temperature and pressure couplings, for the calculation of long-range interactions, for the time step(s) used within the integration algorithms, and so on. Given the extreme sensitivity of Newtonian dynamics to initial conditions, two independent MD simulations will sample the microscopic states with different probabilities no matter how close the initial conditions used8. The difference thus produced in two simulations introduces a variation in results that can often be larger than the quantity of interest, making the results practically useless. Although the chaotic nature of molecular dynamics is mentioned in various textbooks9 it generally receives remarkably short shrift. Extensive studies we have performed in recent years show that molecular dynamics models indeed exhibit extreme sensitivity to initial conditions10-13. From our investigations, we observe that the properties one computes from molecular dynamics trajectories appear superficially to be described by a Gaussian random process (GRP) with a normal distribution denoted by N(,2), characterised by a  and standard deviation  (the square root of the variance in the data). Note, however, that a normal distribution cannot be assumed and in fact there are frequently significant deviations from this form in nonlinear systems14. Moreover, it is entirely possible that molecular dynamics may manifest a pathology which we recently discovered in the simulation of simple chaotic systems on digital computers15. This is caused by the limitations of the IEEE floating point numbers in describing the statistical behaviour of systems with such exquisite sensitivity: ensemble averages, designed to address random errors, also contribute substantial systematic errors to predicted properties15. In a recent study, the effect of box size on simulations of protein dynamics in water was reported16; the authors reported that calculated values of various properties, such as the stabilities of the unliganded and liganded states of human hemoglobin and the density or number of hydrogen bonds per water molecule, changed systematically with increase in box size; the authors maintained that a surprisingly large box of 150 Å was required to obtain meaningful results. Although at first sight this dependency on the box size appeared to be an example of systematic errors in the simulation, it was in fact caused by a lack of reproducibility in the study which became manifest when random errors were taken into account17. Indeed, the ensuing debate16-19 highlights the importance of setting up systems for simulation correctly and, more importantly, applying ensemble approaches to get statistically significant results. As we noted above, without first correctly handling the stochastic errors, it is not possible to assess correctly the nature/size of the systematic errors, and to interpret a finding convincingly. Although it was recognised more than two decades ago that one-off classical molecular dynamics simulations do not generate consistent protein conformations20,21, systematic investigation as to how to make these calculations reproducible has not been performed until recently. Considerable effort has been invested in the development of so-called "enhanced sampling protocols" in order to accelerate phase space sampling, their purpose being to make computed properties more reliable by demonstrating more rapid "convergence" of computed properties. However, in all practical cases it is quite impossible to calculate the (equilibrium or other) probability distribution function, against which expectation values would be calculated. Performing ensemble averaging of replicas of such enhanced sampling simulations shows that there is significant variance in such ensembles22. Ensemble-averaging is not just a practical consideration invoked in the repertoire of uncertainty quantification methods. When molecular dynamics is used, as it frequently is, to estimate thermodynamic properties, such as the free energy of a system, it should be recalled that the connection between microstates (generated by individual MD simulation trajectories) and thermodynamic properties is achieved using ensemble averages. This is true whether the system is in or out of equilibrium. The very common resort to perform socalled "long time averages" of a single microsystem appeals to the ergodic theorem, which is in fact only valid for long times at which this time average should converge to the ensemble average. In reality, that time interval would need to be of the order of a Poincaré recurrence time-a truly astronomical epoch-for the equality to hold. In practice, it is taken to be as long as authors deem to be reasonable. In the ergodic hierarchy of dynamical systems, those which approach and reach equilibrium must be at least mixing8. Mixing systems are ergodic, but the converse is not true. Mixing systems exhibit the tell-tale property of dynamical chaos: neighbouring trajectories, no matter how close, diverge exponentially, at a rate given by a Lyapounov exponent. Since we can never know the true initial conditions for a real system (which arise as a consequence of whatever it was doing before we started to observe it), we are obliged to formulate the approach to equilibrium in probabilistic terms. Indeed, even a single trajectory associated with a given initial condition becomes increasingly inaccurate as time passes, since the exquisite sensitivity of the dynamics means that round-off errors accruing during the time integration of the equations of motion inevitably put the system on orbits other than the one it began on. The lack of reproducibility thus stems primarily from the chaotic nature of classical molecular dynamics. We therefore focus on ensemble averaging, which is mandatory in statistical mechanics, for the convergence, reproducibility, reliability and uncertainty quantification of properties obtained from MD simulations. If we adjudge the system to be in a state of equilibrium, we can in addition perform time averaging, a procedure bereft of meaning out of equilibrium. Extensive studies we and others have performed in recent years10-13,22-33 confirm that the most effective and reliable computational route to reproducible binding free energies of ligands to proteins using MD simulation can be achieved using ensemble methods. The same conclusion has been drawn from MD simulations in other areas, including studies on materials applications such a graphene based systems34 on DNA nanopores using coarse-grained MD simulations35, on rate parameter estimation for binding kinetics36 and so on37,38. An ensemble approach employs a set of independent MD simulations, referred to as "replicas" both in statistical mechanics and within the uncertainty quantification domain, to obtain the required averages and associated. The key feature of such simulations is the use of ensembles and-for systems at equilibrium-time averaging. It is useful to recognise the stochastic nature of these simulations; it can be convenient to approximate the statistical properties of such ensembles as Gaussian random processes8. There is no theoretical means to establish the number of replicas required to produce low errors from ensemble simulation: the criterion for ensuring convergence of the ensemble average is to establish the number N of replicas required such that using N+1 of them makes no significant difference to the expectation values calculated. Of course, this can also be looked at another way, as amounting to a trade-off between the amount of computation one performs (which increases linearly with N) and the size of the error one is willing to tolerate (which reduces roughly as N-1/2). In general terms, the principles discussed above are applicable to all-atom MD and coarse-grained MD. The reduction in the number degrees of freedom achieved by grouping several atoms into single particles or pseudoatoms, as is done in coarse-grained MD, reduces the level of fluctuations in such systems. While this form of coarsening of the model's representation can typically lead to a decrease in accuracy, the benefits which accrue are an ability to study larger systems and for longer time periods; the reduced degrees of freedom also reduce the phase space that needs to be sampled. We find that smaller ensembles are required for coarse-grained MD than all-atom MD. This typically leads to an ensemble with around two-thirds the number replicas as compared to all-atom MD ensembles, as we have shown in recent work with graphene oxide dispersions39 and DNA nanopores in lipid bilyaers40. The requirement to simply run a large ensemble of replicas may sound trivial but it comes with significant overheads in terms of effort expended. Generally speaking, the numbers of replicas required is usually at least ten and frequently many more. So it will be clear that using ensemble methods greatly adds to the computational cost of a study and to the wall clock time unless one has access to modern high performance computers which are equipped with large number of nodes and cores. In such cases, in the time it takes to run one simulation, one can produce the output for all of them. There is then a lot more data to handle, and process. The key step from the overall technical perspective is to bring all these output trajectory data together and then perform the analysis on the aggregate of all of this. This can be done in a number of ways, but the one we generally use is called a bootstrap error41. Given N results from an ensemble of simulations, the bootstrap method involves calculating the distribution of means from resamples of size N from that original results. Many resamples are taken, typically greater than 10,000, are made with replacement. If the original sample is representative of the true distribution, this method can provide error bounds or confidence intervals on any calculated value. The bootstrap error behaves similarly to a standard error; indeed, it is meaningful for quantities that have non-Gaussian distributions. This is of practical value in cases where we do not know the distribution of the quantity of interest. Evidently, automation of some sort is necessary to manage the extra effort involved, and efficient sampling techniques are required to make these kinds of workflow possible. We have previously developed software to assist us in this task for the computation of binding free energies, the so-called "binding affinity calculator"42. This in turn led us to develop software for more general forms of uncertainty quantification, and to extend this to address verification and validation too. In particular, we have developed the VECMA Toolkit4, 43 , as an open source, open development project which is allowing us to apply these methods much more widely, to address uncertainty quantification in a set of diverse domains. For example, the toolkit includes a python library called EasyVVUQ, which permits users to instrument their own codes with capabilities to perform UQ using a wide range of methods, including quasi-Monte Carlo (the method described here), stochastic collocation, and polynomial chaos3,44. The use of ensembles in molecular dynamics simulation only started to be systematically and routinely viable since the advent of the petascale era (that is, in a little over the past ten years), as a result of the vast increase in the number of nodes and cores available on supercomputers. An instructive thing to do is to plot the frequency distribution of observables as it emerges from all the members of an ensemble, as this gives us an indication of the nature of the distributions we can expect. Figure 1 shows a set of examples of the data which typically emerges from these studies. It is clear from these plots that, while they are approximately Gaussian, they exhibit deviations from the standard bell curve expected on the basis that the variables are independent of one another, as one assumes in conventional statistics. Instead, we find that the distributions have a skewness associated with them, the asymmetry favouring the occurrence of values of the observable higher than the mean. The majority of the There are various approaches to estimate the magnitude of the binding free energy (a measure of how strong the interaction is between a ligand and its target protein), based on different theories and approximations45. The Ensemble end-point simulations To generate reliable, precise, and reproducible binding free energies from MMPBSA and MMGBSA approaches, we have proposed an ensemble based MD approach, named "enhanced sampling of molecular dynamics with approximation of continuum solvent (ESMACS)"12. This builds around the so-called MMPB(GB)SA method, including configurational entropy and free energy of association, but with important additional features to address reliability and reproducibility. Correctly accounting for entropic contributions is essential for reliably predicting binding free energies in cases where the ligands are diverse and/or flexible with many rotatable bonds. The contributions can be incorporated to the calculated free energies using normal mode (NMODE) approach or a variety of other options25. We have found a varying number of replicas may be required to achieve a desired level of precision; for many small molecule-protein systems, 25 replicas are typically required with 4 nanosecond production run for each replica, for ESMACS studies12,13. As already noted, the combination of the simulation length and the size of the ensemble provides a tradeoff between computational cost and precision. In collaboration with several pharmaceutical companies, we have used the ESMACS approach to investigate druglike small molecules bound to therapeutic targets23-25,33, and show that ESMACS is well suited for use in the initial hit-to-lead activities within drug discovery. The alchemical approach calculates relative binding free energies between two physical states which are linked by an "alchemical" path. A series of nonphysical steps are involved in the path. The two physical states can be a protein binding with two ligands, or a ligand binding with wild-type and mutant proteins. Along the alchemical path, some atoms change their chemical identities -appearing, disappearing or alchemically transforming from one to another. The alchemical methods thus have a more restricted domain of validity: they are applicable mainly to estimating small relative free energy changes for structures (drugs or proteins) which involve relatively minor (perturbative) variations. The alchemical method can also be used to calculate absolute binding free energies49. It is the equivalent of the relative free energy calculation when one of the ligands involved is replaced by nothing. Free energy calculations using such alchemical methods had rarely been used seriously in drug development projects until recently when Schrödinger Inc. released their "FEP+" simulation software for relative free energy calculations50. With the improved methodology, much of which is proprietary and thus not available for assessment, and the use of graphical processing units (GPUs), FEP+ has made a significant impact in the pharmaceutical industry within its domain of applicability51, although further evaluation is still needed on its accuracy and precision22,30,31. From the perspective of this study, however, it is interesting to observe that the methodology advocated is decidedly based on use of "one-off" simulations, so that any attempt to provide uncertainty quantification is entirely lacking here. As in many other approaches, an alchemical calculation certainly generates random errors, and very likely systematic errors too. As we have stated above, we need to correctly handle the stochastic errors before we can reliably estimate the possible systematic errors. Thermodynamic integration can produce significant differences from separated simulations, up to 1.58 kcal/mol and ~7 kcal/mol for relative and absolute binding free energies, respectively, for the cases tested using 5 independent simulations31. These simulations vary only in their initial As one example among many, the application of TIES to protein mutations provides insights underpinning the impact of the gatekeeper mutation of the FGFR-1 kinase on drug efficacy31. Using an ensemble based approach, we were able to quantify the uncertainties in the free energy calculations, and to compare the performance of different software and hardware for the calculation of the same free energy changes22. Ensemble approaches like TIES provide a reliable, rapid and inexpensive method for uncertainty quantification applied to both relative and absolute binding free energy calculations using alchemical methods22,31. Predicting the properties of modern advanced materials typically requires understanding the structure and the dynamical processes on an atomistic level. Many large-scale, macroscopic, engineering properties can be modelled using methods taken from continuum mechanics, such as the finite element method. However, diffusive processes, self-assembly, structural degradation, surface/interface characteristics, and many other quantities of interest to modern materials scientists and engineers, are all heavily influenced if not controlled by dynamical processes on the scale of atoms and molecules. This is particularly clear in the case of nanocomposites, in which one has to deal with a polymer matrix in which is embedded a nanomaterial such as graphene (and its oxide), carbon nanotubes, clays and such like. Graphene and other so-called two-dimensional materials have one dimension which is of the order of nanometers, and thus just one or a few atoms thick, yet they impart dramatically enhanced large-scale materials properties in the composites they produce. In such circumstances, it is clear that one must use MD as part of the range of techniques available for studying such complex systems. MD techniques are uniquely equipped to explore processes that occur on time and length scales of nanometers to microns, and nanoseconds to microseconds. For complex systems, especially those with anisotropic structures on the nanoscale, such as "soft matter systems", MD has proven useful for predicting the nanoscale structure and material properties. In the development of new materials, the chemical constituents are often among the first known aspects of the system, and the subsequent time necessary for developing useful applications is spent optimising the fabrication and processing for engineering tests. MD can often help to reduce if not remove many of these practical barriers and assess a material's suitability for a given purpose based only on its atomistic structure. Indeed, there is substantial interest in many areas of materials design in virtual testing using computer simulation to speed up the process from concept to real-world implementation, which currently takes of order twenty years, at a cost of many billions of dollars54. The challenge then becomes providing reliable computer based "in silico" predictions which reduce the need for expensive and timeconsuming experimental work. This puts a premium on providing tight error bars since these furnish a key measure of the confidence with which we can accept modelling results and use them to guide experimental work. By way of example, consider trying to predict a material's stiffness using molecular dynamics (see Figure 2) 55. This is one of the most fundamental applications of MD in materials science. A uniaxial stretch of such system is a fairly trivial simulation to perform; however, it poses several fundamental questions about the certainty we can have in a result of an MD simulation. As discussed above, several systematic uncertainties exist with this technique which are hard to quantify56, the most glaring example of which is the choice of force field used to represent the material57. The often quoted limitations of MD-such as finite size/time effects or structure generation-are also systematic errors. Using appropriate workflows we can quantify these effects to produce accurate results58. MD simulations in the condensed phase are typically performed by imposing periodic boundary conditions in all three spatial directions, which means that we only expect to simulate a comparatively small simulation cell to approximate the bulk properties. The size of this simulation cell has many implications for computational cost but, more importantly, the reliability of the scientific results it furnishes. Finite size effects and fluctuations can be expected to affect the outcome of a simulation. To measure the Young's modulus (YM) of a material system, the pressure exerted along one axis is sampled before and after (or during) imposition of a small strain. Since the instantaneous pressure of a molecular dynamics simulation can fluctuate by several GPa, it is necessary to average this value over a long sample period to measure the change in pressure due to an applied strain. In a recent study3, we considered an epoxy resin system, a thermosetting polymer. The investigation quantified the effect of specific MD parameters on the measured Young's modulus, including the system size, starting velocities, and polymer generation random seed. We found that the mean YM of an ensemble of simulations is independent of simulation size but below a box of size 4 nm there is a finite size effect which makes the system artificially stiffer (Figure 3 ). This effect only became evident after performing 300 replica simulations at each box size, with mean YM and standard deviation 3.4  1.9 GPa. The smallest box size gave a distribution of YMs with a significant skewness of -0.8; as the box size increased the skewness tended to zero. The distribution is effectively Gaussian because there are no long-range interactions present. The analysis was greatly facilitated by use of the EasyVVUQ software3. From this, it should be clear that one single measurement, at the low strains imposed here, is wholly inadequate to measure this property reliably. Another benefit of running ensemble simulations is that one can perform sensitivity analysis and thereby learn about the variance in properties arising due to different input variables and parameters. By applying the law of total variance, our study showed that the expected variance due to the polymer network generation was equal to that due to the starting velocities of the atoms in the simulation. In other words, the exact connectivity of monomer units was inconsequential, provided that the same crosslink density is achieved. However, one single simulation of such a polymer, is insufficient; the aforementioned ensemble of 300 replicas was required in order to be 95% confident of the size effect seen above. These approaches are broadly applicable. Running ensembles that give statistically significant results not only allows us to efficiently sample more phase space, therefore increasing the accuracy in the result, but will provides much more information that single simulations would not be able to tell us. It is often forgotten that the starting structure can itself be a major source of uncertainty in MD simulations; it is as true for materials as it is in biomedical simulation. The most straightforward way of achieving variation in an ensemble is to use different random numbers to seed the initial atomic velocities; however, generating different starting configurations should also be considered. The process of building the initial structure and coordinates of a system is not trivial and can often be the most time-consuming step in such research work (it may involve non-trivial polymer chain building for synthetic and biological macromolecular structures, the details about the microstructure of nanocomposites, and so on). The initial state of an MD simulation is inevitably artificial and therefore not itself a representation of the system of interest. Initial states but must be built such that, for example, if one wishes to study the equilibrium state, it will not take an inordinately long time to reach that state by MD simulation. The starting structure must be sampled sufficiently and carefully assessed to check that the unphysical starting conditions do not influence the production stage of a simulation. Polymer systems are manifestly difficult to build while diffusion in high molecular weight polymers can be extremely slow, so entanglements and anisotropic structures must be built carefully as good starting points. Numerous techniques exist to do this. Diffusive processes are so slow that 'sampling the phase space' with one long simulation is clearly impossible; instead several structure generations are essential to be confident in a result. Graphene oxide poses a different problem as its structure is that of an amorphous crystal61. Oxygen containing groups are present across the surface of graphene with some random distribution; in the presence of other materials, the precise distribution of these groups may influence their interaction. In the case of a graphene oxide dissolved in a polymer melt, it is not be sufficient to build one structure and generate an ensemble with different initial velocities; instead one must generate several graphene oxide structures to understand the system. The errors caused by inaccurate force fields can be completely catastrophic for the physics under investigation. A stable inaccuracy in a force field may give a binding energy within some error of the true value, but a more significant inaccuracy may result in divergent dynamical and/or structural regimes62. For example, when predicting the structure of crystals, a flawed force field may never produce certain crystal arrangements. Lennard-Jones forcefields are known to underestimate the friction between layered materials63. Sinclair et al. 34 found that the spherical symmetry of these potentials is too gross an approximation whilst simulating graphene bilayers so a new forcefield was developed. In experiment, graphene flakes are observed to show superlubric behaviour when propelled across a graphitic surface. Propelling flakes in this way is a chaotic process (in the technical sense that it is highly dependent on the initial conditions). In order to achieve an acceptable error which was comparable with experimental data on frictional properties, ensembles of 40 replicas were required. The distance travelled by a flake seems to follow a lognormal distribution, with a Fisher-Pearson coefficient of skewness of 1.4. Ensemble molecular dynamics simulation provides us with a powerful methodology that enables us to connect ergodic theory and uncertainty quantification, and to obtain reproducible results from simulations in a systematic and theoretically well-grounded manner. As evidenced in the case of ligand-protein binding free energy predictions, ensemble simulation-based approaches yield statistically robust, precise and reproducible, hence reliable results. Using ensemble methods, the errors in predictions can be systematically controlled, amenable to further reduction by increasing the number of replicas in an ensemble and in propitious circumstances too by extending the length of such simulations. Ensemble approaches are scalable: for example, they permit hundreds to thousands of binding affinities to be calculated per day, depending on the computing resources available. Computing capabilities are set to increase as the exascale era heaves in to view. In the near future, rapid, accurate and reliable predictions of materials properties may emerge that can be exploited in the aerospace and automotive industries; free energy prediction at high throughput will assist physicians in clinical decision making and medicinal chemists in directing compound synthesis in a routine manner. In sectors such as these, virtual certification and regulatory approval for the use of in silico methods will depend critically on the application of rigorous uncertainty quantification along the lines we have described. Phase Transition for a Hard Sphere System The Alan Turing Institute. Reliability and reproducibility in computational science: Implementing receptor sequences with selected lung cancer drugs Precise, and Reproducible Prediction of Peptide-MHC Binding Affinities from Molecular Dynamics That Correlate Well with Experiment Computing clinically relevant binding free energies of HIV-1 protease inhibitors Big data: the end of the scientific method A New Pathology in the Simulation of Chaotic Dynamical Systems on Digital Computers Valid molecular dynamics simulations of human Prediction of Bromodomain Inhibitors: A Computational Study Application of the ESMACS Binding Free Energy Protocol to a Multi-Binding Site Lactate Dehydogenase A Ligand Dataset Multidimensional Replica Exchange Molecular Dynamics Yields a Converged Ensemble of an RNA Tetranucleotide On the Structural Convergence of Biomolecular Simulations by Determination of the Effective Sample Size Robust Prediction of Resistance to Trimethoprim in Staphylococcus aureus Avoiding False Positive Conclusions in Molecular Simulation: The Importance of Replicas Ensemble-Based Replica Exchange Alchemical Free Energy Methods: The Effect of Protein Mutations on Inhibitor Binding Uncertainty Quantification in Alchemical Free Energy Methods Rapid, accurate, precise, and reliable relative free energy prediction using ensemble based thermodynamic integration Evaluation and Characterization of Trk Kinase Inhibitors for the Treatment of Pain: Reliable Binding Affinity Predictions from Theory and Computation Interactions: Friction, Superlubricity, and Exfoliation Ensemble-Based Steered Molecular Dynamics Predicts Relative Residence Time of A 2A Receptor Binders Toward High Fidelity Materials Property Prediction from Multiscale Modeling and Simulation Ensemble Docking in Drug Discovery Principles Governing Control of Aggregation and Dispersion of Graphene and Graphene Oxide in Polymer Melts Structure and dynamics of DNA nanopores An introduction to the bootstrap Automated molecular simulation based binding affinity calculator for ligand-bound HIV-1 proteases A Library for Verification, Validation and Uncertainty Quantification in High Performance Computing Advances in the calculation of binding free energies Simulations meet machine learning in structural biology Ligand Binding Affinities from MD Simulations How to obtain statistically converged MM/GBSA results Uncertainty Quantification in Alchemical Free Energy Methods Accurate and Reliable Prediction of Relative Ligand Binding Potency in Prospective Drug Discovery by Way of a Modern Free-Energy Calculation Protocol and Force Field Collaborating to improve the use of free-energy and other quantitative methods in drug discovery Replica Exchange with Solute Scaling: A More Efficient Version of Replica Exchange with Solute Tempering (REST2) Statistically optimal analysis of samples from multiple equilibrium states The economic benefits of publicly funded basic research: a critical review The Role of Graphene in Enhancing the Material Properties of Thermosetting Polymers Computation of elastic constants of solids using molecular simulation: comparison of constant volume and constant pressure ensemble methods Prediction of Mechanical Properties of Polymers with Various Force Fields Functional uncertainty quantification for isobaric molecular dynamics simulations and defect formation energies. Model Simul Mater Sci Eng The kinetics of homogeneous melting beyond the limit of superheating Assessing the inner core nucleation paradox with atomic-scale simulations Modeling Nanostructure in Graphene Oxide: Inhomogeneity and the Percolation Threshold Multi-objective optimization of interatomic potentials with application to MgO. Model Simul Mater Sci Eng A Gaussian treatment for the friction issue of Lennard-Jones potential in layered materials: Application to friction between graphene, MoS 2 , and black phosphorus PVC is grateful for many stimulating conversations with Dario Alfè, Bruce Boghosian, Alfons Hoekstra, Peter Sloot and Sauro Succi. All authors contributed to the concept and writing of the article. The authors have no competing interests. We are grateful for funding from the UK EPSRC for the UK High-End Computing Consortium (EP/R029598/1), from MRC for a Medical Bioinformatics grant (MR/L016311/1), the European Commission for the CompBioMed, CompBioMed2 and VECMA grants (numbers 675451, 823712 and 800925 respectively) and special funding from the UCL Provost.