How to learn from inconsistencies: Integrating molecular simulations with experimental data
Simone Orioli, Andreas Haahr Larsen, Sandro Bottaro, Kresten Lindorff-Larsen
HHow to learn from inconsistencies: Integratingmolecular simulations with experimental data
Simone Orioli , *, Andreas Haahr Larsen , , *, Sandro Bottaro , and Kresten Lindorff-Larsen Structural Biology and NMR Laboratory & Linderstrøm-Lang Centre for ProteinScience, Department of Biology, University of Copenhagen, Copenhagen, Denmark. Structural Biophysics, Niels Bohr Institute, Faculty of Science, University of Copen-hagen, Copenhagen, Denmark. Present address: Structural Bioinformatics andComputational Biochemistry Unit, Department of Biochemistry, University of Ox-ford, Oxford, United Kingdom. Atomistic Simulations Laboratory, Istituto Italianodi Tecnologia, Genova, Italy. ∗ These authors contributed equally to this work
Keywords : Molecular simulations; Integration with experiments; Force fields; Maxi-mum Entropy; Bayesian Methods; Time-dependent; Time-resolved
Abstract
Molecular simulations and biophysical experiments can be used to pro-vide independent and complementary insights into the molecular origin ofbiological processes. A particularly useful strategy is to use molecular sim-ulations as a modelling tool to interpret experimental measurements, andto use experimental data to refine our biophysical models. Thus, explicitintegration and synergy between molecular simulations and experimentsis fundamental for furthering our understanding of biological processes.This is especially true in the case where discrepancies between measuredand simulated observables emerge. In this chapter, we provide an overviewof some of the core ideas behind methods that were developed to improvethe consistency between experimental information and numerical predic-tions. We distinguish between situations where experiments are used torefine our understanding and models of specific systems, and situationswhere experiments are used more generally to refine transferable models.We discuss different philosophies and attempt to unify them in a singleframework. Until now, such integration between experiments and simula-tions have mostly been applied to equilibrium data, and we discuss morerecent developments aimed to analyse time-dependent or time-resolveddata. a r X i v : . [ phy s i c s . c h e m - ph ] D ec ontents Introduction
Molecular mechanics simulations of biological systems have matured over thelast decades, to a state where reliable predictions and interpretations of biolog-ical phenomena can be achieved [1, 2]. Simulations can readily be used eitherbefore experiments to suggest hypotheses and design experiments, or after ex-periments to analyse and interpret the data. With simulations it is possible toprobe timescales and spatial details that are, yet, impossible to access exper-imentally, and thus they provide a unique tool to study, e.g. conformationaltransitions [3–7], signalling events [8], ion transport [8, 9], and many other phe-nomena. Since most experimental observables are averaged over time and a largenumber of molecules, simulations can help resolve the underlying structural anddynamical distribution with atomistic spatial accuracy and femtosecond timeresolution.As simulations are in essence a theoretical framework, experimental verifica-tion is pivotal and there are still, in many cases, significant differences betweenexperimental data and the corresponding observables calculated from state-of-the art simulations [10]. These discrepancies are usually due to imperfect forcefields [11–14], insufficient sampling [15] or inaccurate forward models [16]. Evenwith accurate forward models and robust sampling, however, the quality ofthe simulation results are bound by the accuracy of the underlying force field[13]. Indeed, molecular mechanics is intrinsically limited by the fact that it em-ploys classical approximations for phenomena that are known to happen on aquantum scale. Force fields have been developed to approximate these quantuminteractions within a classical framework, in such a way to balance between com-putational simplicity and accuracy: thus, an exact match is not to be expectedand synergy between simulations and experiments remains of fundamental im-portance. Substantial efforts of the biophysical community have been thereforedirected towards the development of methods that correct the results of simu-lations by accounting for experimental information, either in a system specificor generalisable fashion. In this chapter we provide a self-contained overview ofthe principles and strategies that have been applied to infer conformational en-sembles, i.e. collections of structures reflecting the dynamism and plasticity ofa molecule, as well as to improve force fields using experimental data. In otherwords, we will make a clear cut distinction between system-specific and generalforce field corrections: in the former case, simulations are performed accordingto experimental conditions, and if inconsistencies emerge, the resulting ensembleis corrected according to the available experimental information; in the latter3ase, instead, discrepancies between simulations and experimental observablesare used to actively improve the physical potential employed to describe thesystem. We note that the goal of the chapter is not to be comprehensive of allthe literature that has been previously published: rather, it focuses on somecornerstone principles and methods that summarise the efforts of the commu-nity. For this reason we also focus mostly on the methods and theory, and onlyprovide few examples of applications.This chapter is meant for readers experienced in molecular simulations andcomfortable with the fundamentals of statistics. It is particularly suited forscientists that approach the combination of experiments with simulations forthe first time or readers interested in a non-technical overview of the field.The chapter is comprised of seven sections, each one concerned with a dif-ferent aspect of the literature on the subject. Section 2 is dedicated to thedescription of some of the main strategies to obtain consistency between sim-ulation and data by manipulating the ensemble after simulations have beenperformed. Differently, in section 3 we will discuss how consistency can beachieved by introducing a system-specific empirical energy term in the forcefield. In this case, the refinement step occurs before the actual simulation, asthe experimental bias will guide the sampling towards only the relevant regionsof conformational space. In section 4 we discuss how to employ experimentaldata to refine the physical description of macromolecules, i.e. the force field,instead of the conformational ensemble of a particular system. Like the biasedforce field approaches described in section 3, these methods also adjust the forcefield before the production simulation. In section 5 we discuss how to combinetime-dependent and time-resolved data with simulations. This is particularlyrelevant, as computational power [17–20] and enhanced sampling techniques[3, 21] have pushed the reachable timescales to the ones resolved by some ex-perimental techniques. In section 6 we discuss some of the most importantoverall challenges arising from combining experimental data and simulations,and, finally, in section 7 we conclude by summarising similarities and differ-ences between the various approaches. We thus aim to give an overview of thepresent state of the field as well as pointing out key aspects where further collec-tive effort is needed to improve the predictive power of combined computationaland experimental methods for structural studies of biological systems.4
Reweighting strategies
Consistency between simulations and experimental data can in many cases beachieved by reweighting a trajectory (or a set of trajectories) obtained with agiven force field (Fig. 1). Let us assume to perform a molecular dynamics (MD)or Monte Carlo (MC) simulation and to save N snapshots, X i , . . . , X N , fromthe trajectory for analysis. Each frame in the simulation is given an initialweight that reflects the population of that structure as predicted by the forcefield. We denote the initial set of weights, w = w , w , ..., w N , the referencedistribution and we assume them to be non-negative and normalised as prob-abilities, (cid:80) i w i = 1. When performing a standard MD/MC simulation andunder the assumption that enough sampling has been performed, the initialweights are constant, w i = 1 /N , as such simulations generate conformationsdistributed according to the Boltzmann distribution. When using enhancedsampling techniques that do not directly sample the Boltzmann distribution or,in the case of shorter off-equilibrium simulations, the initial weights are in gen-eral non-uniform [22]. Given the N snapshots and the reference distribution, astatic observable can, in the simplest case, be calculated as a weighted ensembleaverage: (cid:104) O calc (cid:105) ( w ) = N (cid:88) i =1 w i O calc ( X i ) . (1)Some experimental quantities (e.g. NMR relaxation rates and other inherentlytime-dependent observables) cannot generally be expressed as linear ensembleaverages over individual structures, i.e. employing Eq. (1). Therefore, through-out the manuscript the validity of Eq. (1) will be taken as a working assumptionuntil section 5, where we will discuss the case of time-dependent observables.In a reweighting process, the weights are modified until the calculated ob-servables become consistent with corresponding experimental observables, O exp j ,where j = 1 , . . . , M . We deonote the reweighted set of weights, w = w , w , ..., w N ,the optimal distribution [23]. Unfortunately, many possible sets of weights canlead to the same averages, and thus to consistency with data. I.e. it is an ill-posed problem, and strategies are therefore needed to regularise the problem,so the most optimal of these distributions can be selected. There exist severalstrategies to regularise and solve the resulting optimisation problem [24], andwe here focus on three of them: the principle of maximum entropy (MaxEnt,section 2.1), the principle of maximum parsimony (MaxPars, section 2.2) andBayesian inference (MaxPrior, section 2.3). Each method encodes a different5 xperimental information R e w e i g h t i n g Bare FF Probability density MaxEnt/MaxPriorMaxParsSampling Probability densityExperiment-Biased FFSampling
Analysis AnalysisSimulation Simulation
Figure 1: Overview of methods to achieve consistency between simulation andexperimental data for the distribution of a generic static observable. All methodsstart with an initial force field E ff ( x ) (‘bare’ FF, blue line) which is eitheremployed as is or modified in accordance to some experimentally driven bias(experiment-biased FF, red dashed line). In the upper case, the biased force fieldis sampled to directly generate an observable distribution which is coherent withexperimental data. In the lower case, the observable distribution is subsequentlyreweighted against experimental data using either the MaxEnt/MaxPrior orMaxPars principles to yield some reweighted distributions (red dashed line).6hilosophy for what it means for a set of weights to be optimal , which translatesinto different optimal distributions after reweighting. The principle of maximum entropy (MaxEnt) [25] states that the optimal dis-tribution is the one adding the least amount of information given some imposed constraints , e.g. (cid:104) O calc (cid:105) = O exp . In other words, the solution with highest en-tropy and still consistent with the data is considered optimal. In this context,entropy is quantified via the relative entropy which we define as the negative ofthe Kullback-Leibler divergence [26]: S ( w ) = − N (cid:88) i =1 w i log (cid:18) w i w i (cid:19) , (2)where w i and w i are the weights of frame X i before and after reweighting.Other expressions for the entropy are possible [27, 28], but Eq. (2) is usuallypreferred [29–32]. As the distributions of initial and refined weights diverges, S ( w ) becomes smaller, and vice versa. In the simplest case, the set of weightsthat maximises the entropy is chosen by the only constraint that the weightsare normalised as probabilities, i.e. we look for weights for which the function g ( w ) = 1 − (cid:80) i w i is zero. This problem can be solved by constructing aLagrangian: L ( w , λ ) = S ( w ) − λ g ( w ) , (3)where λ is a Lagrange multiplier. The solution of the set of equations ∇L = 0 (4)gives the highest entropy weights which fulfil the normalisation condition, whichare simply w = w . It is clear then, that the problem is interesting onlywhen more constraints are introduced. For example, exact consistency witha given set of experimental data can be imposed via functions of the form g j ( w ) = O exp j − (cid:104) O calc j (cid:105) ( w ). The Lagrangian in Eq. (3), in the case of M experimental constraints and the normalisation constraint, becomes: L ( w , λ ) = S ( w ) − M (cid:88) j =0 λ j g j ( w ) , (5)where a Lagrange multiplier λ j is introduced for each constraint. The solution ofEq. (4) for the Lagrangian in Eq. (5) provides weights that ensure consistency7etween the reweighted trajectory and the experimental data. The optimalweights are uniquely given in terms of the M Lagrange multipliers [33–37]: w i ( λ ) = w i exp (cid:104) − (cid:80) Mj =0 λ j O calc j ( X i ) (cid:105)(cid:80) Nk =0 w k exp (cid:104) − (cid:80) Mj =0 λ j O calc j ( X k ) (cid:105) , (6)where O calc j ( X i ) is the value of the j -th observable, as calculated from the i -th frame in the simulation. Note that the denominator just ensures propernormalisation of the weights.Fitting data tightly to gain full consistency by imposed constraints might,however, be unrealistic and lead to overfitting because of errors in the data andin the model [34]. Indeed, experimental observables are only known with somelimited certainty and, likewise, forward models used to calculate observablesfrom simulations have an associated uncertainty [16]. Systematic errors in thedata lead to further disagreement with the model [23, 38, 39], and finally, thefinite number of structures N used to compute the theoretical observable (seeEq. (1)) introduces additional error to the estimated averages [39]. All theseeffects add up to an effective uncertainty for each observable, σ j . Therefore it isgenerally preferable to apply restraints rather than constraints, i.e. to imposeconsistency with the data only within the estimated uncertainty (restraining)rather than exactly (constraining). This translates into relaxing the equality (cid:104) O calc (cid:105) = O exp to a similarity (cid:104) O calc (cid:105) ∼ O exp within a given threshold.There are two main strategies to introduce restraints in the MaxEnt methodusing Lagrange multipliers. The first one consists in modifying the constraintsin such a way that the calculated and experimental observables are allowed todiffer by some small quantity ε , i.e. g σj ( w ) = ( O exp j + ε ( λ j , σ j )) − (cid:104) O calc j (cid:105) ( w ). ε depends on the effective uncertainty, σ j and the Lagrange multiplier associatedto the constraint, λ j [36, 40, 41]. The second way to include experimental errorswas introduced by Gull and Daniel [42] and, in the assumption of normallydistributed errors, describes the discrepancy between data and model via the χ distribution: χ ( w ) = M (cid:88) j =1 (cid:32) O exp j − (cid:104) O calc j (cid:105) ( w ) σ j (cid:33) , (7)The expectation value of χ is equal to the number of degrees of freedom ν , soa reduced χ can be defined such that its expectation value is unity, χ r = 1 ν χ . (8)8xperimental restraints can thus be included by imposing the constraint func-tion g χ r ( w ) = χ r ( w ) − θ − ,and whose role we discuss in further detail below. This treatment is correct inthe assumption of normally distributed errors and for datasets that guarantee χ r = 1. The latter assumption is not safe, as point estimates can deviate fromthe χ r expectation value: more details on this are provided in section 6. Toovercome this problem, the optimal distribution can be obtained by minimisingthe equation [23, 29, 32] T MaxEnt ( w ) = χ ( w ) − θS ( w ) , (9)where θ is treated as an adjustable parameter rather than a Lagrange multiplier,and controls how tightly the data should be fitted. The solution is then saidto be regularised by the entropy, i.e. at a fixed value of θ the optimal solutionis the one with the highest entropy S ( w ) among the ones that minimise thediscrepancy with the data, χ ( w ).After reweighting, some of the weights may be close to zero and the cor-responding structures will thus effectively be ignored in the calculation of thenew averages. MaxEnt ensures that the optimised ensemble preserves as manystructures from the reference one as possible. In this sense, the entropy term S can be interpreted via a more intuitive quantity, φ eff , which represents theeffective fraction of frames used in the reweighted ensemble compared to theinitial ensemble [32]: φ eff ( w ) = exp[ S ( w )] . (10)If the reference distribution is unaltered, φ eff ( w = w ) = 1. In the oppositeextreme all weight is given to a few frames of the simulation and φ eff ( w ) (cid:28) K = (cid:16)(cid:80) Ni =1 w i (cid:17) (cid:80) Ni =1 w i . (11)When a single frame is dominant over all the others, Eq. (11) returns K = 1,while in a situation where no particular frame is preferred one has K = N . The principle of Maximum Parsimony (MaxPars), also referred to as
Occam’srazor , is another strategy that can help choosing among several models consis-9ent with the data. As we shall discuss later, many different interpretations ofOccam’s razor exist. In this section, the principle is interpreted as follows: theoptimal distribution coherent with a given set of data is the one providing thesmallest possible ensemble while still fitting the data. The basic idea is that if amodel (ensemble) with few parameters (structures) can explain the data, thereis no reason to further complicate the model by introducing more parameters.The above principle can for example be quantified by the Akaike informationcriterion [45]: AIC = 2 n − L, (12)where the optimal solution is found by simultaneously maximising the likeli-hood, L , that quantifies agreement with experimental data, and minimising thenumber of model parameters, n . While the number n is effectively equal tothe number of frames in the ensemble, it is customary to associate it to thenumber of non-zero weights, so n (cid:28) N . For normally distributed experimentalrestraints one has log L = − χ ( w ) /
2, so: T MaxPars ( w , n ) = AIC = χ ( w ) + θn, (13)where θ = 2 if Eq. (12) is directly applied. θ was introduced here to emphasisesimilarity with the MaxEnt methods, i.e. that the difference is in the form ofthe regularisation term.The Akaike information criterion was directly applied by Bowerman et al.[46] to find a minimal ensemble of tri-ubiquitin, using MD simulations andSAXS data. Bouma et al. [47], on the other hand, used θ as a free adjustableparameter to tune the strength of the MaxPars regularisation term. In mostmethods however, the optimal n is found by increasing it by incremental stepsand, for each step, finding the ensemble with lowest χ ( w ). The optimal n is then found when χ ( w ) has converged, either judged by manual assessment[48, 49] or by an automatic convergence criterion [50–52], thus avoiding havingto set an explicit value for θ . An alternative method to MaxEnt and MaxPars is provided by Bayesian statis-tics [23, 37, 53, 54]. While MaxEnt regularises the negative log-likelihood bythe entropy, and MaxPars by parsimony, Bayesian inference employs the priorfor the same goal. For this reason we will refer to this approach as
MaxPrior .10n Bayesian inference all available information is expressed by means of prob-abilities, as captured in Bayes’ theorem: P ( w | O exp , σ ) = P ( O exp | w , σ ) P ( O exp , σ ) P ( w ) . (14) P ( w ) is called the prior and it quantifies the information known about the sys-tem prior to the introduction of experimental data. In case of reweighting, itrepresents the prior probability associated to the weights, and typically comesfrom the distribution encoded in an energy function. The term P ( O exp | w , σ ),known as the likelihood L ( w ), represents the probability of measuring the ex-perimental observables, given the set of weights and uncertainties. Assumingnormally distributed errors, L ( w ) is given as: P ( O exp | w , σ ) = L ( w ) ∝ exp (cid:18) − χ ( w ) (cid:19) . (15)Note that this is not the only possible expression for the likelihood. More com-plicated expressions, even non-analytical ones, can occur if different types ofdata are combined or when sources of error cannot be assumed to be Gaussian[55, 56]. The term on the left-hand side of Eq. (14) is known as the posteriorand represents the probability of the weights after experimental data have beenconsidered. Finally, the term in the denominator is treated as a normalisa-tion constant in which case Bayes’ theorem takes the form of a proportionalityrelation: P ( w | O exp , σ ) ∝ P ( O exp | w , σ ) P ( w ) . (16)Given Eq. (16), Bayesian inference defines the optimal distribution by followinga two-step procedure: first, the information from data is included, i.e. a like-lihood function is defined; second, all other forms of data are included, i.e. aprior on the weights is provided. If we assume the employed force field to rep-resent the best estimate of our prior knowledge on the system, then the priorprobability P ( w ) must decrease as w deviates from w . This observation cane.g. be quantified through the relative entropy S ( w ) from Eq. 2 [32, 57, 58], inwhich case the prior takes the form: P ( w ) ∝ exp (cid:18) θ S ( w ) (cid:19) . (17)The factor 1 / T MaxPrior ( w ), that must be minimised to find the optimaldistribution, w : − P ( w | O exp , σ ) ∝ T MaxPrior ( w ) = χ ( w ) − θS ( w ) . (18)It is important to note that, differently from the case of MaxEnt where θ wasempirically used to tune the strength of the regularisation term, in a Bayesianframework the interpretation of θ is clear. Indeed, it accounts fot the uncertaintyon the weights, σ w i , and thus effectively on the force field, i.e. θ = ( σ w i ) − .This means that if the uncertainty on the weights, data and model were knownaccurately, θ could be exactly determined. Unfortunately, the uncertainty onthe force field parameters, and thus on the weights, is generally unknown and θ must therefore be treated as a hyperparameter of the model. In practice, θ plays a double-role: it expresses our trust in the force field and at the same timeit can be used to compensate for over- or underestimated experimental errors,as well as errors in the forward model [23, 64]. See more on the determinationof θ in section 6.The usefulness of Bayes theorem in refining structural ensembles from sim-ulations and data is clear from the amount of studies on the subject, all includ-ing the word Bayesian in their title [23, 31, 38, 39, 50, 57, 59–67] (see also theoverview in [68]). The methodology has become so widespread in the field, thatone could argue that researchers should start stating in the title when Bayesianmethods are not used, rather than the opposite. We note also that differentmethods may differ substantially in how and the extent to which they apply theBayesian formalism including whether priors are defined over all parameters andwhether these are integrated out during the procedures.Several methods combine MaxPars, MaxEnt and MaxPrior approaches (seean overview in Table 1 in the review by Bonomi et al. [68]). Bayesian priorscan, e.g., be constructed to prefer minimal ensembles [63, 69], thus combiningMaxPrior with MaxPars. Additional information can also be used together withprior information from the simulation, e.g. secondary structure restraints, likein the case where MD simulated data are fitted into low-resolution cryo-electronmicroscopy density maps [70]. Therefore, the Bayesian approach should be seenmore as a toolbox than an actual principle for selection among ensembles.The Bayesian framework also provides an interpretation of the principle of12arsimony, which is different from the Akaike Information Criterion, where themodel with fewer parameters is considered the simplest one. In a Bayesiansetting, the principle of maximum parsimony is a matter of surprise : the fur-ther a posterior result is from the prior, the more surprising it is. This leadsto a Bayesian Occam’s term: a measure of the distance between the referencedistribution (prior) and the optimal distribution (posterior) [71, 72]. The prin-ciple is implicitly built into Bayes theorem. In that sense, the S ( w ) term in theBayesian inference and MaxEnt methods can be interpreted as an Occam’s term,and both these methods could claim to fulfil the principle of maximum parsi-mony. However, to keep notation clear, we will use the Akaike interpretation ofOccam’s razor. Different interpretations of the concept of optimal distribution can lead to dra-matically different conformational ensembles. Indeed, MaxEnt includes as manyconfigurations of the initial ensemble as possible, while MaxPars tries to includeonly the ones that are strictly necessary to model the data. The behaviour ofMaxPrior, on the other hand, depends on the prior employed. In most ap-proaches however, a term similar to MaxEnt is used, preferring solutions thatare consistent with the reference distribution [68], so that MaxPrior will lead tosolutions with as many frames as possible. Therefore, we will only distinguishbetween MaxEnt and MaxPars in the following, when discussing which methodis more appropriate under what circumstances. To our knowledge, no system-atic and direct comparison of the different reweighting methods have been made.However, the methods have some clear and important differences that will bediscussed below.
A relevant point of comparison between the methods is how easily the resultsare interpreted. The results of the MaxPars method are easier to interpret andvisualise, as the optimal ensemble usually contains only 2-5 representative struc-tures [65, 73], whereas the MaxEnt method may result in an optimal ensemblewith thousands of structures, whereof many are similar.13 .4.2 General applicability
As argued by Bonomi et al. [68] and Ravera et al. [24], both methods can be usedwhen the system’s free energy landscape shows a few distinct and well-definedminima. In that case, the ensemble is represented well by a few structures, onefor each free energy minimum and weighted by its corresponding depth (seeFig. 1). On the other hand, if the energy landscape representing the ensembleis flat or high-dimensional, i.e. in case of high-entropy systems, a few structuresprovide a poor description of the ensemble, even if they may fit the data. Soin the case of high-entropy systems, the MaxEnt method should be preferredover MaxPars. A simplistic example of a high entropy system is a regulardie[24, 33, 35]. For an unbiased die, all six faces, which represent the states ofthe system, are equally likely. Let us consider the experimental observation thatthe average result of throwing the die is 3.5. The application of MaxPars to thisproblem would lead to the conclusion that the die is well described by two statesonly, e.g states and with weights w = w = 0 .
5, or states and withweights w = 0 .
25 and w = 0 .
75, or any other pair of states consistent withthe observation. This result is ambiguous, as many combinations are equally optimal , and it is also a poor description of the regular die. MaxEnt, on theother hand, would lead to a uniform distribution for all outcomes. More relevantexamples of high-entropy systems are multi-domain proteins with flexible linkers[74], intrinsically disordered proteins [75], and unfolded proteins [76], which areall described poorly by a few distinct states. Thus, the MaxEnt method ismore versatile, in the sense that it gives reliable results also for high-entropysystems [24, 68]. The dimensionality of the system also plays an importantrole. For reasonably small systems, a description in terms of low-dimensionalfree-energy landscapes might be possible and relevant. I.e. a few structurescan adequately represent the system. However, for high-dimensional biologicalsystems the assumption that a few coordinates, or states, can properly describethe collective dynamics of the molecule might be ventured. Recent studies havehighlighted how also the dynamics of some fast-folding proteins might be toocomplex to be represented in fewer than ∼
10 dimensions [77]. Therefore, wesuggest that MaxEnt, or related methods, should be considered the method ofchoice when no safe assumptions can be made about the free energy landscapeof the ensemble. 14 .4.3 Imperfect force fields
The choice of the reweighting method also depends on the quality of the sim-ulations used to generate the reference distribution w . Indeed, the results ofMaxEnt and MaxPrior (with e.g. Gaussian [59–62] and Dirichlet priors [63])reweighting will strongly depend on the quality of the simulation, because theoptimal distribution is expected not to deviate too much from w . On the otherhand, MaxPars reweighting depends less on the quality of the initial distribu-tion, as it freely picks the simulated structures that best represent the data.This decoupling from the initial simulation is both the strength and the weak-ness of the MaxPars method. It implies that, in principle, a realistic referencedistribution is not necessarily preferred over a poor one, as long as some goodrepresentative structures have been sampled (Fig. 2). Consequently, approxi-mate but efficient simulation methods, e.g. Monte Carlo sampling with implicitwater, can be used in conjunction with the MaxPars method (e.g. Rosetta [78]),allowing for a fast exploration of the conformational space and better overallsampling without the need for enhanced sampling methods. Such approximatesimulation methods are evidently less suitable for the MaxEnt and MaxPriorreweighting and therefore, more costly explicit solvent MD simulations shouldbe employed. On the other hand, MaxPars does not benefit much from thequalified prior information of a good force field. We want to stress here thatone of the main advantages of using accurate simulations as a prior is whenthe provided experimental dataset is sparse. In that case, the refined ensemblesobtained from MaxEnt/MaxPrior are generally more reliable than MaxPars, asthey are regularised by the simulated model. Reweighting methods, as we presented them here, are essentially optimisationproblems. As such, they are cursed by the problem of dimensionality [79] or,in other words, the higher the dimension of the parameter space the more com-putationally challenging it becomes to find the true minimum of the functionalof interest. Several strategies have been employed to reduce the computationalburden of reweighting.In MaxEnt, for example, the principal numerical challenge is to find the setof weights that minimise Eq. (5). This can become highly nontrivial when N islarge, e.g. for tens of thousands of structures. The computational complexitycan nonetheless be alleviated by optimising the values of the M Lagrange mul-tipliers (Eq. (6)) rather than the N weights, as usually M (cid:28) N [23, 29, 32]. We15 axEnt Experimental information
Sampling
Probability density
SimulationSimulation A n a l y s i s Simulation Analysis R e w e i g h t i n g Inexact FF
Experiment-Biased FF
Unreliable FF
SamplingPoor sampling
Probability density MaxPars
Figure 2: Lower two rows: The two cases of poor sampling of an inexact forcefield and the one of full sampling of an unreliable force field can both lead to a re-sult where the distribution of an observable of interest is poorly estimated. Theperformance of subsequent MaxEnt/MaxPrior and MaxPars reweighting differ,as a poor prior inevitably leads to a poor reweighted distribution in the case ofMaxEnt/MaxPrior, while MaxPars can still select some relevant conformationsfrom the initial unreliable pool. Top row: In both cases, experiment-biasedsimulations can help alleviate the problem.16ote however that, depending on the implementation, the optimisation of the N weights can be faster than the one of the Lagrange multipliers [37]. Moreover,the functional in Eq. 18 is strictly convex for finite θ [80], so the correspondingoptimisation can be carried out in a straightforward fashion by means of highlyefficient routines such as limited memory Broyden-Fletcher-Goldfarb-Shanno(L-BFGS) algorithm [81].In MaxPars, the numerically challenging part is provided by the test of the N ! / ( n !( N − n )!) combinations of n -sized ensembles out of N structures. Anaive approach where all the possible combinations are tested is intractable forrealistic situations as, for example, in the case of n = 5 and N = 1000 thenumber of combinations becomes C ( n, N ) ∼ . Therefore, the best fittingcombinations of n structures has to be determined with some other strategies,e.g. via Monte Carlo minimisation [82] or more complex algorithms such as theone implemented in ASTEROID [83]. Other strategies limit the search to thesubset of structures that best match the experimental data individually [52],which however violates the principle that the goodness of fit should only beassessed only against the ensemble average (Eq. (1)), and not against the fit tosingle structures.As a final remark, we note that both MaxPars and MaxEnt methods canbenefit from clustering, where structures are grouped according to their struc-tural resemblance, to reduce N before reweighting [29, 60, 61]. In section 2 we focused on three possible principles, MaxEnt, MaxPars andMaxPrior, to reweight the results of simulations a posteriori . However, as weargued in section 1, this is not the only possible approach. Instead of changingthe simulated ensemble after the simulation has been performed, the simulationcan be guided by adding a bias to the underlying force field, E ff ( X ), employedin the simulation to ensure consistency with the experimental data. In thissection we provide an overview of the methods that have been applied withthis goal, which are either based on the MaxEnt principle (section 3.1), on theaddition of empirical energy terms (section 3.2), or on the MaxPrior principle(section 3.3). Finally, in section 3.4 we compare and discuss the differencesbetween experiment-biased simulations and a posteriori reweighting strategies,discussing their possible advantages and shortcomings.17 .1 Maximum Entropy In section 2.1, we discussed how Maximum Entropy can be used to obtainconsistency between data and simulations by optimising a set of Lagrange mul-tipliers (see Eq. (6) and Ref. [36]). The mathematical formulation of theprinciple makes it immediate to extend it to the problem of experiment-biasedsimulations. Starting from an unbiased force field E ff ( X ) (Fig. 1), the idea isto perturb it as little as possible under the constraint/restraint of consistencywith some given experimental data. Maximum Entropy principle can be ex-pressed in terms of the probability of a conformation X given the new forcefield P ( X ) ∝ exp( − βE ( X )), where β = ( k B T ) − with k B being the Boltzmannconstant and T the temperature of the system, and the probability of the sameconfiguration obtained from the unbiased force field P ff ( X ) ∝ exp( − βE ff ( X ))[34, 41]: S [ P ( X )] = − (cid:90) dXP ( X ) log (cid:18) P ( X ) P ff ( X ) (cid:19) . (19)Practical examples of this approach are provided by experiment directed sim-ulations (EDS) [84] and experiment directed metadynamics (EDM) [85]. Theeffective error from experiment and forward model can also be accounted for inthe Lagrangian framework, thus only restraining the solution within the givenuncertainty [40, 41]. Rather than approaching with MaxEnt, experimental restraints can be directlyincluded as empirical penalty terms [86]: E ( X ) = E ff ( X ) + E exp ( (cid:104) O calc (cid:105) , O exp , σ )= E ff ( X ) + M (cid:88) j =1 θ j h j ( (cid:104) O calc j (cid:105) , O exp j , σ j ) , (20)where the sum is carried out over M restraints. θ j are the force constantsassociated to each restraint, while the specific choice for the functional formof each h j depends on the distribution of the corresponding observable. Theconstraints are usually expressed as squared residuals [39, 86] or, when errorsare taken into account, as χ terms [68] (Eq. (7)).For concreteness and simplicity, we assume here normally distributed ex-perimental observables. We also assume that the same force constant can beused for all observables. This is a good assumption, e.g. when the observables18ome from SAXS data [38], but cannot generally be assumed when more thanone experimental techniques are combined. With this assumption, the energyfunction reduces to: E ( X ) = E ff ( X ) + θχ ( (cid:104) O calc (cid:105) , O exp , σ ) . (21)Just as in reweighting approaches, parameter θ in Eq. (21) takes into accountunknown uncertainties of the force field parameters as well as unknown or im-perfectly determined errors in the data and the forward model. Therefore, θ is generally not known and determining it is one of the key challenges of themethod [23, 37].Another issue in the practical implementation of this approach is that theaverage (cid:104) O calc (cid:105) can only be determined after the simulation is over. Thus therestraints need to be applied iteratively [35]. As an alternative approach to ob-tain an ensemble averages at each point in the simulation, replica methods havebeen introduced [23, 87–89], where N independent replicas of the same systemare simulated and observables are calculated as ensemble averages from these.Interestingly, replica experiment-biased methods mathematically converge to aMaximum Entropy constrained solution as N → ∞ and θ → ∞ , as discussed byPitera & Chodera [90], Roux & Weare [34] and Cavalli et al. [91], and reviewedby Boomsma et al. [35]. Very recently, K¨ofinger et al. [37] combined reweightingand experiment-biased methods to simultaneously ensure a large ensemble (asin MaxEnt) is considered and that all relevant states are visited by adding abiasing energy term to the force field. Just as in the case of MaxEnt, Bayesian inference, or MaxPrior principle (seesection 2.3), can be also employed to bias a priori molecular simulations ratherthan just reweighting them a posteriori . The method known as Metainference[39] is an implementation of this principle: it employs a Bayesian approachto quantify how much the prior is modified by the introduction of noisy andheterogeneous sources of data. In its essence, the strategy works by running N replicas of the system and guiding the sampling by means of a log-posteriorscoring function: s ( X, σ ) = − N (cid:88) i =1 log P ( X i , σ i ) + ∆ ( X ) N (cid:88) i =1 σ i , (22)where σ takes into account all the sources of errors, P ( X i , σ i ) is the prior andthe factor ∆ ( X ) estimates the deviation between the predicted observables19nd the experimental ones. It can be shown that in the single-replica limit, N = 1, Eq. (22) reduces to Eq. (18). Therefore, the log-prior plays the role ofan effective entropy term, while the second component of s ( X, σ ) is a χ termcomputed over the replicas. It is interesting to notice here that the equivalencebetween Eq. (22) and Eq. (18) is valid only in the θ = 1 case: this comesfrom the fact that the method proposes a model to take into account all thesources of error [39], which means that θ is not expected to be a free parameteranymore. We notice that Eq. (22) is only valid for Gaussian error sources, andmore complicated and complete expressions can be obtained in the general case(see the discussion in section Materials and Methods in Ref. [39] and Ref. [92]).Metainference has been effectively combined with Metadynamics [93–95] inits parallel bias formulation [92, 96]: this synergy enables one to explore the con-figuration space in an efficient way while simultaneously sampling conformationsthat are coherent with experimental data. Reweighting methods can be used with many different types of simulations andforce fields, as the reweighting process is independent from the simulation andsampling (Fig. 1). This makes it a rather adaptable module-like tool, with theinput being the trajectory and the experimental data only [32, 36, 44]. Also, newexperimental data can easily be incorporated in the reweighting process withouthaving to re-run the conformational sampling. In the experiment-biased sim-ulations, the implementation is more specific, as the empirical energy term isan integrated part of the simulation [38]. Decisions about types of experimentsand force field constants have therefore to be taken before the simulation. It isstill possible, however, to reweight experiment-biased simulations a posteriori to remove or add sets of experimental data, though the procedure is techni-cally somewhat more challenging than reweighting MD trajectories; indeed, theapplied experimental bias has to be estimated and subtracted from the simula-tion before further reweighting [23, 44]. For the same reason, experiment-biasedsimulations are often not carried out together with enhanced sampling tech-niques. A notable exception is provided by metadynamics [97], which has beencombined with metainference [39, 92] to increase its efficiency, and in principlemetadynamics (and other enhanced sampling methods) could be used in othermethods for experiment-biased simulations.20 .4.2 Forward models
When a trajectory is reweighted a posteriori , the forward model is only evalu-ated on the frames that are to be reweighted, which are typically only a smallfraction of the frames generated during the simulation. Also, as long as theobservables are calculated according to Eq. (1), the calculations are done inde-pendently of one another and may thus be easily parallelised. For these reasons,the forward model can be of high complexity (e.g. quantum calculations canbe carried out on the ensemble structures [98]). This is not the case, however,of experiment-biased simulations: when implemented within a molecular simu-lation, forward models are evaluated and differentiated at (almost) every step,so they have to be sufficiently simple to assure computational efficiency [44]and their gradients have to be known analytically. Therefore, for complex for-ward models reweighting might be more applicable than a priori experimentalbias. In some cases, such as for NMR chemical shifts, it is possible to employ afast forward model [99] that is almost as accurate as more refined and complexmodels [100], which would be too complex to be computed at each step in asimulation. In other cases, it might not be possible to derive sufficiently compu-tationally efficient and accurate forward models. We suggest that a possibilitywould be to use simpler and less accurate models to bias the simulations andthen reweight a posteriori the simulation with the more realistic models.
Reweighting methods and experiment-biased simulation methods may in prac-tice perform differently in cases where the force field provides a relatively poordescription of the system’s conformational space. With a poor force field, somerelevant states may be rarely or never visited when running an unbiased simu-lation (Fig. 2). Consequently, reweighting methods may fail in predicting thecorrect average of an observable. An indication that this is happening is usu-ally provided by the fraction of effective frames φ eff becoming close to 0 and bythe reweighted distributions of relevant observables being skewed towards theexperimental average (Fig. 2). In principle, this can be overcome by sufficientsampling, but it might be computationally very expensive and practically un-realistic. Experiment-biased force fields, on the other hand, can better providereasonable results even when a rather poor unperturbed force field is used as ba-sis, as the empirical energy term will alter the energy landscape such that evenrelevant states with high energies in the unbiased force field become reachable(Fig. 2). The effect of the added empirical energy term can be monitored by21omparison with a control simulation without the additional energy term ( θ = 0in Eq. (21)). In sections 2 and 3 we described techniques to refine simulations in a system-specific manner by including experimental data. From a different perspective,substantial progress has been made when using experimental data to improvethe force field as a general and transferable predictive tool. In this section wefocus on this point and review some of the fundamental advancements in thesubject.
As widely discussed in literature and assessed by empirical knowledge, the qual-ity of the physical description provided by force fields is fundamental for theaccuracy of biophysical simulations. The reliability of these simulations indeeddepends critically on the ability of the underlying physical description to effec-tively model all the relevant inter-atomic interactions. After decades of forcefields development [101–104], MD simulations have reached a high level of re-liability and the ever-growing amount of experimental observations calls for asystematic and detailed comparison of theoretical predictions, coming from MDsimulations, with the available data.Unfortunately, force fields do not always provide results that are in per-fect agreement with experimental findings. To understand the reasons and thesources of these emerging discrepancies, we shall first recall how a force field isusually designed. For a comprehensive introduction on the subject we refer thereader to Chapter 1 of this book, and here we instead focus on how experimentaldata may be used in force field parameterization.The typical force field is composed by two fundamental elements: (i) A func-tional form E ( X ). This element embeds our physical understanding of molecu-lar processes by providing a classical parametrisation of the Born-Oppenheimerenergy surface [105]. The typical functional form of the force fields used forbiomolecular simulations (with minor re-adjustments among the different inter-22retations) is given by [106]: E ( X ) = E bonded ( X ) + E non-bonded ( X )= (cid:88) bonds k b ( r − r ) + (cid:88) angles k θ ( θ − θ ) + (cid:88) dihedrals k φ [1 + cos( nφ − φ )]+ (cid:88) i,j ∈ non-bonded (cid:34) A ij r ij − B ij r ij + q i q j r ij (cid:35) , (23)where X is a configuration of the system and r , θ and φ are functions of theatomic coordinates and q are atomic charges. k b , k θ and k φ denote the dif-ferent strengths of the interactions and are usually tensors, as they depend onthe specific group of atoms involved in the interaction; (ii) A set of parameters ξ = ( r , θ , φ , k b , k θ , k φ , . . . ). The functional expression of the force field de-pends on the choice of these parameters which set, for example, the strengthof interactions, equilibrium distances and angles. To determine their values, itis necessary to fit them against known experimental and quantum mechanical(i.e. ab initio) properties. The specific choice of these properties depends onthe philosophy underlying the force field development.Historically, force field development has always been a daunting task be-cause of its technical complexity and the amount of time, experimental dataand simulations required [107–110]. To give an example of this, let us focuson the AMBER class of force fields. In its first version [104], bonded angleswere fit to the vibrational frequencies of single amino acids or small molecules,in order to reproduce experimental frequencies; fixed charges were fit to repro-duce the results of quantum calculations [111–113]. Lennard-Jones parameters,instead, were set in such a way to reproduces enthalpies of vaporisation anddensities in organic liquids [101]. Finally, dihedral and torsional angles were fitto reproduce quantum calculations of single amino acids or experimental barrierheights of small molecules. With the increasing quality of experimental and abinitio data, modern and widely used versions of the AMBER force fields [114]have reached a high level of complexity. Nonetheless, the underlying generalphilosophy remained the same: fit the force field parameters to ab initio andexperimental data of single amino acids or small molecules and compare theresults against data available for larger systems. Force field ff19SB, obtainedby only employing ab initio simulations, constitutes a notable exception [115].Despite its historical significance and great success, this procedure shows someimportant limitations: (i) It is difficult, within this approach, to improve theforce field parameters by making use of discrepancies with respect to experimen-tal data on large biomolecules; (ii) errors induced by the forward models used to23alculate experimental observables are not consistently taken into account, anderrors on the force field parameters are not estimated and/or not used. In thelast decade, many new strategies for force field parametrisation have been de-veloped, that less strictly follow the traditional approaches and instead embracea more Bayesian approach [116–124]. These efforts, combined with the growinginterest in automated force field parametrisation (as exemplified by the develop-ment of the ForceBalance framework [121, 122, 124]), partially solve the issuesreported in point (ii). A deeper discussion on point (ii) will be left for section 6.In this section we focus on point (i) or, more precisely, describe and summarisethe strategies available in the literature to optimise force field parameters usingdiscrepancies between MD simulations and experimental observations on longerpeptides or even entire proteins. As we will see, these methods are tightly con-nected to common reweighting algorithms (see section 2) but, when applied toa wide set of molecules, provide a an answer of more general purpose than thesole ensemble refinement. NMR data have been widely and successfully used in combination with ensemblereweighting strategies of proteins [30, 32, 48, 57, 58, 61, 62, 64, 83, 125–135];other approaches employed various sources of experimental data, e.g. for smallmolecules. This fact opens to the appealing possibility to use NMR data onfull proteins to directly optimise force fields. Such an approach is more generalthan ensemble reweighting, in the sense that experimental NMR data can becollected for a larger number of proteins and used to increase the quality of theforce field rather than optimising the description of a single system of interest.To understand how this could be possible in practice, we need first to realise thatthe general scheme of ensemble reweighting methods basically comprises of threesteps: (i) Generation of an ensemble of structures with a given force field; (ii)Calculation of experimental observables from the ensemble; (iii) Determinationof weights, associated to each structure, that better describe the experimentalobservables.If this procedure is repeated for a large set of proteins, the informationcarried by the optimised weights can be in principle used to increase the qualityof the underlying physical model, i.e. the values of the force field fit parameters,rather than just the single molecule ensemble (Fig. 3). We note, however, thatthis approach is less flexible than, e.g., MaxEnt where the weights of all thesimulation frames can be fine tuned. When optimising a force field, one is24nevitably constrained by its functional form and only limited adjustments canbe made.In the next sections we review the different strategies that have been de-veloped to compute and use the weights in force field refinement for proteins(section 4.2.1) and RNA (section 4.2.2).
Let us assume we can access a given set of N snapshots X i obtained from asampling strategy of choice (MD or Monte Carlo) and generated using a givenforce field E ff , defined by a set of parameters ξ . These snapshots can be usedto calculate the ensemble average of some observables of interest, (cid:104) O calc j (cid:105) , forwhich we possess experimental information O exp j . The goal is to determine aset of force field parameters ξ that decreases the discrepancy between experi-mental and simulated averages. This discrepancy can be estimated via the χ (Eq. (7)), whose minimisation with respect to the force field parameters wouldmaximise the compatibility between simulations and the experiments of choice.The search for the χ minimum can in principle be done in a brute force fashion,by infinitesimally perturbing the set of force field parameters by a quantity δξ i and re-simulating the system of interest. However, this strategy is particularlyinefficient, as it corresponds to a search in a high-dimensional parameter space,where the evaluation of each new trial force field requires a full re-sampling ofthe system’s conformational ensemble. To avoid this, it is possible to use someideas coming from statistical mechanics. In principle, if we knew the exact forcefield E exact describing our system, we could reweight each configuration usingthe Boltzmann relationship: w exact i = w i e − β ( E exact ( X i ) − E ff ( X i )) . (24)Eq. (24) is the same idea behind the method of free energy perturbation [136]devised by Zwanzing to determine the changes in the free energy of a systemwhen a perturbing potential is introduced. In this case, force field E ff playsthe role of the unperturbed potential of the system, while the exact force field E exact is the contribution to the total energy provided by the perturbation.Even though E exact can never be known, Eq. (24) can be used as the basis of anefficient optimisation strategy [116, 137–139]. In this method, a force field E old is iteratively optimised to obtain a new one, E new , that is used to recomputethe weights of the parent MD simulation via: w new i = w old i e − β ( E new ( X i ) − E old ( X i )) . (25)25 xperimental informationHigh overlap with old FF RefinementCalculation of new weightsLow overlap with old FFForward modelInitial FF Convergence New FFWeightsSamplingSimulationRefined FF i-th FFObservable Figure 3: Schematic representation of the force field optimisation process. First,an initial force field E ff ( X ) is employed to sample an ensemble of configurationsof some systems of interest. A forward model is then applied to the computedframes, to reproduce the known experimental average of a set of relevant ob-servables. If the computed and the experimental observables disagree, furtherexperimental information is used to refine the current force field parameters anddefining a new force field E i ( X ). A weight is associated to each frame obtainedfrom the previous sampling using Eq. (25). If the overlap between the old andthe new force field is high (i.e. most of the weights do not change substantially)the new weights are used to refine the estimation of the observable’s average,which in turn is used to further refine the force field parameters and to estimatenew weights. If the overlap is instead small, a further round of conformationalsampling is carried out with the new force field E i ( X ). The process is repeateduntil some convergence criterion is satisfied.26he new weights w new i are used to refine the estimation of the simulated averageby means of Eq. (1). E new force field is obtained by minimising the χ , wherethe specific strategy for minimisation can change (e.g. Levenberg-Marquardtprocedure [116, 140], simplex minimisation [137] or simulated annealing followedby simplex minimisation [138]). The strength of the reweighting strategy in Eq.(25) is that it does not require one to sample new conformations every timethe force field parameters are updated. Nonetheless, it only works under theassumption that minimal perturbations of the force field are introduced by the χ minimisation and therefore that the snapshots X i could have been reasonablysampled also by means of the E new force field. For this reason, a sensibleoverlap between E old and E new force fields is expected to exist, and particularcare is dedicated to the assessment of this overlap [138, 141], which can bequantified, e.g. through φ eff (Eq. (10)). A re-sampling of the full conformationalensemble needs to be repeated every time the overlap between the two forcefields becomes too small: when φ eff < ε , with ε (cid:28) α , C β and C (cid:48) carbon atoms of 4 trial proteins and subsequently benchmarking it against a testset of 18 proteins of different topologies. The obtained results show an average27mprovement of the comparison with the experimental data for the proteins inthe test set, and was later refined further [138].More recently Chen et al. [139] extended the method from Norgaard etal. [116] to use Markov State Models (refer to section 5.1 for more details onMarkov State Models) as an intermediate step to calculate observables from sim-ulations and parameterize the conformational landscape. Despite the differentframework, the general philosophy of the method, named ODEM (Observable-driven Design of Effective Molecular models), remains unaltered. ODEM hasbeen applied to design a C α − C β coarse-grained model of protein FIP35 whichis able to reproduce relevant pair-distance distributions measured by FRET. Molecular simulations may also be used to study the conformational landscapesand dynamics of RNA molecules. Unfortunately, the precision of RNA forcefields is still limited and not comparable to the one reached by force fieldsdesigned for proteins [150–153]. The approach of integrative structural biologytherefore becomes a powerful tool to enhance the comprehension of fundamentalprocessed governing RNA dynamics and indeed the effectiveness of a posteriori reweighting on RNA simulations has been established by several works [32, 154–157]. The lack of a common strategy to increase the reliability of RNA forcefields from first principles [158–160], however, makes it particularly appealingto resort to approaches that exploit discrepancies between MD simulations andexperimental data to refine force fields. Here we review a recent study [134]that tackles this problem by proposing a likelihood minimisation scheme. Letus assume the system of interest is described by a force field E ff ( X ) and acorresponding Boltzmann probability distribution P ( X ) ∝ e − βE ff ( X ) . For thesame system, M experimental data have also been collected. In this work, theauthors seek for an optimised probability density P ( X, µ ) ∝ P ( X ) e − β (cid:80) Ni =0 µ i f i ( X ) (26)where now the force field is expressed by an expansion on a basis f i ( X ), whereeach function is associated to a weight µ i . We stress that weights µ have notto be interpreted as the previously introduced weights w , as the former arearbitrarily normalised and not interpreted as probability densities. Each of the N terms helps in reproducing the M experimental constraints, but, differentlyfrom MaxEnt methods, f i ( X ) do not represent the forward model connectingconfiguration X to an experimental measurements. Rather, f i ( X ) can be genericfunctions and N is sought in such a way that N (cid:28) M . While the analytic form28f the basis functions is enforced at the beginning of the optimisation process,the value of the weights λ i needs to be determined through the minimisation of afunction describing the discrepancy between the simulated and the experimentalobservables. We note that in this strategy the functional form of the force fieldis actively modified by introducing basis functions f i ( X ) that are not explicitlyincluded in E ff ( X ). Alternatively, this can be seen as a way to associate non-zeroweights to terms in the force field that have an effective null coupling constant.If functions f i ( X ) are already part of the potential, instead, the correspondingweights µ i amount for a refinement of the interaction strength. Supposing tobe interested in M observables O calc j , the target for the minimisation takes theform of a regularised error function T ( µ ) = T ( (cid:104) O calc1 (cid:105) ( µ ) , . . . , (cid:104) O calc M (cid:105) ( µ )) + θ | µ | θ ≥ (cid:104) O j (cid:105) ( µ ) are computed in the refined ensemble (cid:104) O j (cid:105) ( µ ) = 1 Z (cid:90) dX O calc j ( X ) P ( X, µ ) (28)The error function is designed to enforce both equalities and inequalities, i.e. (cid:104) O calc j (cid:105) ( µ ) = O exp j and (cid:104) O calc j (cid:105) ( µ ) < O exp j and the optimal set of weights µ isobtained as the one minimising T . The strength θ of the regularisation term iskey to the strategy: for θ → ∞ the error function does not feel the contribu-tion of the experimental constraints and thus the potential resulting from theoptimisation is just the original one E ff ( X ). Instead, for θ → θ has tobe determined or alternative statistical tools have to be applied to integrateout the variable. In this work the authors employ cross-validation to determinethe optimal θ , but this is a general issue common to many ensemble optimisa-tion methods. We refer therefore to section 6 and references therein for furtherdetails on the choice of the optimal prior strength.The RNA systems chosen for this application were four tetranucleotidesand two tetraloops, and NOEs (Nuclear Overhouser Effect) as well as scalarcoupling NMR data were used in the force field refinement procedure. Thebasis functions were chosen to be sines and cosines acting on torsional angles,while the employed guess force field was Amber ff99bsc0 + χ OL3 with the OPCwater model [104, 161–163]. 29
Matching time-dependent and time-resolveddata
When dealing with systems at equilibrium, the average of an observable (cid:104) O (cid:105) isconsistent with many possible distributions of the same observable p ( O ). Aswe discussed in section 2, MaxEnt, MaxPars and MaxPrior principles can helpdiscriminate, among the several possibilities, which one is the optimal distri-bution. Suppose, however, we can access a good description of some relevantconformations of a molecule (by experimental evidences or previous sampling)and we are now interested in understanding the dynamics of interconversionbetween them. From an MD perspective, one can hope to see interconversionshappening by running several trajectories starting from the collected relevantconfigurations. If extensive sampling is achieved and enough state transitionsare captured, one can for example use the trajectories to build a Markov StateModel for the system, and obtain the desired information about the inter-statedynamics [164–166]. However, full sampling of the conformational dynamicsmight be challenging to achieve, especially when the desired motions happen inthe µ s timescale or above. Moreover, as discussed in detail in section 4, limi-tations in the force field employed in the simulations might lead to conclusionswhich disagree with key experimental data, even in the case of excellent sam-pling. This may be the case, for example, for NMR spin relaxation experiments:spin relaxation rates are sensitive to both short- and long-range dynamics of aprotein, occurring from the picosecond to the nanosecond timescale [167]. Fordynamical observables like NMR spin relaxation rates, methods developed toincrease the consistency between computed and experimental values of staticobservables might be of only limited help and different strategies are needed.In this section we describe some recent efforts in this direction and how it ispossible to use simulations to match experimental knowledge on time-dependent observables (i.e. observables that depend on time because of fast processes hap-pening in the system and cannot therefore be expressed by means of Eq. (1), e.g.NOEs with spin diffusion [135]) and time-resolved ones (i.e. observables that areat equilibrium locally in time, but have been monitored for long timescales andcan therefore be expressed as in Eq. (1) by including a time label, e.g. time-resolved SAXS). The following sections will focus on applications concerningMaximum Entropy and Likelihood estimation applied to Markov State Models[168, 169] (section 5.1), the principle of Maximum Caliber [170] (section 5.2)and average block selection [171] (section 5.3) (Fig. 4).30 ugmented Markov Models
20 1 w p ik
20 1 w p ik Experimental information(a)
Maximum Caliber Principle Average Block Selection (b) (c)Experimental information w w w w Figure 4: Graphical summary of methods implemented to recover time-dependent and time-resolved data. (a) Augmented Markov Models employ ex-perimental information to refine a guess Markov State Model and consequentlybetter estimate kinetic observables of interests; (b) Maximum Caliber, in itsexperiment-biased formulation, is employed to bias the sampling to match sometime-resolved experimental quantities; (c) Average block selection is used toassign weights to sub-trajectories and better estimate time-dependent experi-mental data. 31 .1 Maximum Entropy and Likelihood in dynamical sys-tems
In a kinetic model, the equilibrium distribution w can be determined from atransition probability matrix T ( τ ), i.e. a matrix whose elements p ik providethe probabilities for the system to be found in state k at time t + τ , given itwas in the state i at time t . T ( τ ), together with the structures of the statesamong which transitions happen, is usually known as a Markov State Model (MSM) [165, 168]. In this section, we will discuss a method to refine a MSM byadding static (not time-dependent or time-resolved) experimental information.This framework, that balances data coming from simulations and experimentalaverages, is called Augmented Markov Model (AMM) [172].The construction of an AMM starts from the definition of an initial MSM, T ( τ ), with transition probabilities p ik . The MSM is typically built by followinga multi-step procedure [165, 168]: (i) relevant features (e.g. dihedral angles,contact maps etc.) of the system are selected. The dynamics of these featuresneeds to play a key role in the conformational transitions under consideration,so their identification might not be trivial. This task can become easier withthe help of automatic selection tools [173, 174]; (ii) the features are used togenerate a lower-dimensional representation of the system dynamics [175, 176];(iii) the lower dimensionality [177] allows one to cluster similar configurationsin the dimensionally reduced space. The N clusters will define the states of themodel; (iv) finally, statistical tools [168] are employed to estimate the transitionmatrix T ( τ ). Among other possibilities, the transition probability matrix canbe obtained by maximisation of the likelihood [165, 168]: L ( T ) ∝ (cid:89) i,k ( p ik ) c ik , (29)where c ik is the number of transitions occurring between state i and k , which areobtained by explicitly counting the transitions between states in the simulation[178, 179]. The resulting transition matrix has to satisfy further criteria, e.g.it has to be row- or column-stochastic, and detailed balance can be enforcedexplicitly. Maximisation of Eq. (29) is equivalent to minimising the negativelog-likelihood: S ( T ) = − N (cid:88) i,k c ik log p ik , (30)We note that while Eq. 30 resembles an entropy term, but it cannot represent aproper one, as the counts c ik enter Eq. 29 un-normalised. However, S ( T ) can32e used in a similar way as a regularising prior when including experimental in-formation, as we describe further below. The equilibrium distribution of the sys-tem, usually referred to as the stationary probability distribution of the MSM, isobtained from the transition probability matrix by solving the eigenvalue equa-tion T ( τ ) w = w . For more information about Markov State Models, werefer the reader to some excellent reviews on the subject [164–166, 168, 180].Given the initial Markov State Model T ( τ ), the framework of AugmentedMarkov Models employs the MaxEnt principle and likelihood maximisation tobuild an optimised MSM, T ( τ ). Let us assume we can access both the expec-tation values of M experimental observables and an MSM built on simulationdata, with the reference stationary distribution w . To obtain the optimal dis-tribution w , one can apply MaxEnt principle, as reported in Eq. (6), whichbridges the model distribution with the experimental one by means of the pro-vided experimental averages. Rather than estimating the weights correspondingto the optimised distribution w , the MaxEnt formulation allows one to optimisethe set of Lagrange multipliers λ (Eq. (6)). As usual, the Lagrange multipliersare obtained by enforcing constraints on the experimental averages. In orderto account for the statistical errors in both sampling and experiments, one canintroduce the so-called Augmented Markov Model likelihood , assuming Gaussianerrors: L ( w , T ) ∝ (cid:89) i,k ( p ik ) c ik e − χ ( w ) . (31)The term in the parentheses is the MSM likelihood, Eq. (29), while the secondterm incorporates the Gaussian error model proposed for the observables. Thecorresponding negative log-likelihood is thus proportional to: T AMM ( w , T ) = χ ( w ) − θS ( T ) , (32)where S ( T ) is the entropy-like term introduced in Eq. (30). Note that w is uniquely determined by T by solving the eigenvalue equation T ( τ ) w = w ,so Eq. (32) is effectively a minimisation problem for the matrix elements p ik alone. We stress that we have introduced Eq. (32) only to show the similaritywith the reweighting methods for equilibrium ensembles, Eq. (9) and (13), andthat it does not directly represent the way AMMs are practically implemented.Rather, one obtains an AMM that optimally balances between the informationfrom simulation and the experimental data (Fig. 4a) by employing fixed-pointiteration algorithm to maximise Eq. 31 with respect to the unknowns p ik and theweights w (see the supplementary information of Ref. [172] for further details).Olsson et al. [172] applied this method to two 1 ms simulations of ubiquitin.33hey built AMMs using NMR scalar and residual dipolar couplings, i.e. staticdata, and their results were compared against NMR relaxation dispersion data.The results of such experiments can either be calculated directly from a longmolecular dynamics simulation [18, 181] or indirectly from a MSM [182]. Theauthors found that the AMM that had been optimized against the experimentaldata was in overall better agreement this indpendent data compared to the rawMSM. This is a strong indication of a non-trivial fact: reweighting an MSM withrespect to experimental equilibrium observables helps increase the reliability inthe prediction of time-dependent data. The principle of Maximum Caliber (MaxCal) can be regarded as a generalisationof MaxEnt, and enables one to compute the probabilities associated to dynam-ical pathways, rather than probabilities (weights) associated with equilibriumstates. To infer pathway probabilities, MaxCal principle maximises a path en-tropy [183] defined over all the possible pathways, constrained to reproduce agiven dynamical observable. The path entropy is defined as: S [ p ( γ )] = − (cid:88) γ p ( γ ) log p ( γ ) , (33)where p ( γ ) is the probability that the system follows the structural path γ .Each path is considered as a collection of configurations, γ = { X γ , . . . , X γT } ,labelled by a discrete time index t = 0 , . . . , T . Maximisation of Eq. (33),together with M constraints of the form: g j [ p ( γ )] = 0 j = 1 , . . . , M (34)yields the optimal distribution p ( γ ) over pathways. Just like in the MaxEntcase, constraints are usually enforced by Lagrange multipliers, and typicallythey concern the calculated dynamical observables O calc ( X γt ) and correspondingexperimental values at the same time point t , O exp t : (cid:88) γ p ( γ ) O calc ( X γt ) − O exp t = 0 (35)The normalisation of the probability distribution of pathways is also enforced: (cid:88) γ p ( γ ) − . (36)If it is possible to introduce an estimate for the time-dependent error σ t , thenthe notion of χ can be extended to a time-dependent form as well, χ = χ t .34n this way, by analogy with MaxEnt, MaxCal principle can also be expressedin a regularised fashion as T MaxCal [ p ( γ )] = χ t [ p ( γ )] − θ t S [ p ( γ )] (37)where the parameter θ of Eq. (9) acquires a dependence on time. All theterms in Eq. (37) are functionals of the probability density in the space ofpathways: as such, the minimisation of T MaxCal [ p ( γ )] would require a searchin path space, which is knowingly a hard task [184–187]. To our knowledge,indeed, MaxCal principle has never been used in the acceptation of Eq. (37).Instead, the principle of Maximum Caliber can be employed to run restrainedMD simulations, where restraints are provided by time-resolved experimentaldata (Fig. 4b).Note, when comparing simulations with time-resolved experiments, a singlesimulation cannot be used to generate an ensemble as in the MaxEnt method(section 2), due to the non-equilibrium conditions. Therefore, replicas areneeded to describe ensemble development over time (see also section 3). Capelliet al. [188] introduced such replica-averaging implementation of MaxCal, soensembles can be determined at each time t . In this formulation, Eq. (35) and(36) are employed to constrain the minimisation of the path entropy, togetherwith the path probability density normalisation and two more equations. Thefirst one constrains the system’s diffusion constant D :12∆ t (cid:88) γ p ( γ )[ X γt +1 − X γt ] − D = 0 , (38)while the other constrains the standard deviation σ t of the observable of interestat time t , obtained by averaging it over all the N replicas (cid:88) γ i p ( γ i ) (cid:16) (cid:104) O calc ( X γ i t ) (cid:105) − O exp t (cid:17) − σ t = (cid:88) γ i p ( γ i ) ξ t − σ t = 0 , (39)where (cid:80) γ i is used to specify that the sum is carried out over all the possiblepathways of all replicas and the average (cid:104)·(cid:105) is computed over replicas. The ex-pression of the path entropy, together with all the constraints (in this particularcase, Eq. (35), (36), (38) and (39)) imposed via Lagrange multipliers is called caliber . The result of caliber maximisation provides the optimal distribution ofpathways: p ( γ ) = 1 Z exp − (cid:88) t,i ( ν it (cid:2) X γt +1 − X γt (cid:3) + λ it O ( X i,γt ) + µ it ξ t ) , (40)35here Z is a normalisation factor, ν it , λ it and µ it are replica- and time-dependentLagrange multipliers. Despite the complicated expression of the optimal pathprobability distribution, it is possible to prove that, in the assumption wherethe system’s dynamics is Brownian and by analogy with MaxEnt [188–190], theMaxCal distribution in Eq. (40) can be sampled by adding a time-dependent,harmonic bias potential to the underlying force field: E exp ( γ i , t ) = N k (cid:0) (cid:104) O calc ( X γt ) (cid:105) − O exp t (cid:1) , (41)where the constant k is used to tune the strength of the interaction. This ap-proach was applied to the case of the second hairpin of protein G B1 domainby employing synthetically generated time-resolved SAXS data [188].Rather than using MaxCal principle as an experiment-biasing technique,it is also possible to employ it for reweighting, in particular for the case ofMSMs [191–193]. Suppose we built a MSM from a set of trajectories and wewant to predict how a given non-equilibrium observable O ik changes in thetransition between state i and state k . A possible example would be the changein a spectroscopic signal (e.g. F¨orster resonance energy transfer or circulardichroism) when passing from a partially unfolded state i to a helical state k ofa protein. The average, computed over many transitions, is provided by: (cid:104) O calc (cid:105) ( w ) = (cid:88) i,k w i p ik O calc ik , (42)where w is the stationary probability distribution of the MSM. Let us supposewe also possess experimental information of such observable, O exp and that thecomputed average does not match the expected result. In this case, we can usethe MaxCal principle to select among the models with the correct average [192].To do so, first the entropy is built as: S ( T, w ) = − (cid:88) i,k w i p ik log (cid:18) p ik p ik (cid:19) , (43)where w is the optimal stationary distribution and p ik are the refined transitionprobabilities. Note that Eq. (43) is a version of Eq. (33) in the case of discretepathways. The entropy has to be maximised together with the constraint thatthe function g ( w ) = (cid:104) O calc (cid:105) ( w ) − O exp (44)is zero, where (cid:104) O calc (cid:105) ( w ) is the average over transitions introduced in Eq. (42).Moreover, w i and p ik are interdependent because of probability conservation,36o three more constraints naturally emerge: (cid:88) i w i p ik − w k = 0 (cid:88) k w i p ik − w i = 0 (cid:88) i,k w i p ik − . (45)The minimisation of the entropy in Eq. (43), with the constraints in Eq. (44)and (45), can be carried out with the method of Lagrange multipliers. It yields: p ik = ηφ k φ − i M ik ( λ, p ik , O calc ik ) (46)where M ik ( λ, p ik , O calc ik ) is a non-Hermitian matrix, φ is its only right-eigenvectorhaving only positive elements (whose existence and uniqueness is guaranteed bythe Perron-Frobenius theorem), η the corresponding eigenvalue and λ is the La-grange multiplier controlling the constraint in Eq. (44). The positivity of theelements of vector φ is necessary to guarantee the positivity of the transitionprobabilities. Eq. (46) provides a new set of transition probabilities definingan optimal transition probability matrix T ( τ ). Effectively, similarly to whathappened in the case of AMMs, the introduction of experimental informationactively modifies the transition probabilities between states, and interconver-sions can become more or less favourable depending on the cases. Finally, theoptimal stationary distribution is obtained as: w i = ψ i φ i , (47)where ψ is the left-eigenvector of matrix M ik ( λ, p ik , O calc ik ) corresponding to theeigenvalue η . Matrix M ik ( λ, p ik , O calc ik ) only depends on the Lagrange multi-plier λ , which is used as a parameter, and on quantities that can be readilycomputed from the initial MSM, i.e. the transition probabilities and the ob-servables of interest. Therefore, the optimal transition probabilities and equi-librium probability distribution can be obtained by single value decompositionof M ik ( λ, p ik , O calc ik ) for different values of λ : the optimal Lagrange multiplierwill be the one for which the constraint in Eq. (44), estimated a posteriori fromthe updated MSM, holds true.In the original work, the method was tested on a toy model for a growthfactor activation pathway, showing how the incorporation of experimental infor-mation can help correcting the transition probabilities in a guess MSM. Despitethe method, as illustrated here, does not incorporate experimental errors, it ispossible to take them into account in a Bayesian fashion [192]. Time-dependent data can also be applied directly in the ensemble refinement bymeans of time-dependent average block selection. Suppose we have N MD tra-37ectories, started from different conformations, and a set of M time-dependentexperimental observables O exp j . As anticipated at the beginning of this section,by time-dependent we here mean observables that cannot be calculated usingEq. (1) because the values depend on the underlying dynamics as well as on thestationary distribution. To increase the accuracy of the MD predictions, giventhe experimental data, and preserve information about dynamical processes,Salvi et al. [171] proposed to divide each trajectory into B blocks representingsubsequent time-windows. The total number of blocks is then N × B , each oneassociated to a weight w b . The block-weights are then optimised by minimisingthe residual sum: R ( w ) = M (cid:88) j =1 (cid:32) O exp j − N · B (cid:88) b =1 w b O calc bj (cid:33) (48)where O calc bj is the j -th observable computed from the b -th block. In later studies,the authors included errors in the expression by the usual χ ( w ) expression [194].Division in blocks allows to determine how much each block contributes to thetime-dependent experimental signal. For example, blocks with null weights areexcluded and can be interpreted as non-physical artefacts. Despite the similaritywith other approaches used for equilibrium observables, this method shows twoimportant differences: on the one hand, weights are not associated to a singlestructure of an ensemble, but rather to sub-trajectories (Fig. 4c). It is evident,then, that the method can be effective only in the case where the experimentaltimescales of interest can be sampled within a single trajectory. On the otherhand the method as it stands is unregularised, i.e. it assumes that the optimalsolution is obtained by minimising the sum of residuals without adding anyrestraint. This is a crucial difference with respect to MaxEnt, MaxPars andMaxPrior strategies: one would expect nonetheless the method to benefit fromregularisation terms. The addition of such terms would be advisable in furtherapplications, but it is not clear what kind of prior would be needed in thiscase. Because of the fact that Eq. (48) associates weights to pathways and notconformations, it might be possible to explore a connection with a posteriori reweighting using MaxCal principle.In its original formulation, the method (called ABSURD — Average BlockSelection Using Relaxation Data) [171] was designed to deal in particular withNMR spin relaxation rates, but its application can be extended to any experi-mental source sharing similar timescales in a straightforward way. ABSURD wasapplied to the C-terminal domain of the nucleoprotein of Sandai virus. A cumu-lative time of approximately 6.5 µ s of MD was sampled, employing two different38ater models. While the reproduction of spin relaxation rates by MD simula-tions alone was imperfect, the ABSURD-optimised trajectories were found inmuch better agreement with relaxation data on a wide spectrum of timescales,even in the simplest case where a single experimental rate was employed in theoptimisation. In sections 2–5 we focused on major efforts that have been done to merge in-formation coming from experiments and simulations to refine conformationalensembles or improve existing force fields. Despite many important steps for-ward in the last decade, some major challenges remain. In this section we willpoint out some of these and try to draw a possible paths for solutions. In section6.1 we will focus on the problem of setting the parameter θ in a robust way. Insection 6.2 we will discuss the option of combining force field corrections andreweighting for a given system. In section 6.3 we discuss possibilities and obsta-cles to employ kinetic data to reweight equilibrium ensembles. Finally, section6.4 is dedicated to a discussion on the next generation of force fields. In the MaxEnt and MaxPars reweighting methods, and equivalent methods bias-ing on the fly, the prior knowledge coming from simulations is balanced againstexperimental data by tuning the parameter θ (Eqs. (9), (13) and (21)). In prin-ciple, the balance could be exactly known if the effective uncertainties on theweights, force field, sampling, and on the both the calculated and experimentalobservables were known (as also discussed in section 3.3). However, the uncer-tainties on the weights are usually unknown and the experimental uncertaintystems from different sources, only some of which are typically known or easy toestimate. We can subdivide the sources to the experimental uncertainty in thefollowing categories: (i) Statistical errors on experimental data; (ii) statisticalerrors in the calculation of average observables from a limited amount of struc-tures (Eq. (1)); (iii) systematic errors on experimental data; (iv) inaccuracyof the forward model. For normally distributed errors, these add up to a totalvariance [23]. In the following we will shortly discuss each of them.The statistical error on experimental data is usually estimated by repetitivemeasurements and counting statistics, and it is included in most approaches,e.g. as standard deviation in the χ . Over- or underestimated errors may be39etected by visual good fits having χ r much lower or higher than unity.The statistical error on the mean of the observables calculated from the en-semble can easily be taken into account [68]. This is important when the ensem-ble is small, which is typically the case for MaxPars methods and experiment-biased approaches with replicas.Systematic errors in experimental data are generally difficult to take intoaccount, as their magnitude and nature are usually unknown and may be sys-tem specific (e.g. incorrect buffer subtraction in SAXS [38]). Shevchuk andHub [195] treated the systematic errors with Bayesian statistics as a nuisance parameter, i.e. an unknown (and uninteresting) parameter that should be de-termined together with the model parameters. Bonomi et al. [39] discussedthe more general case of outliers, and showed that an approach that combinesreweighting and experiment-biased force fields is more robust against outliersthan other related methods. Similar methods are also discussed by K¨ofinger etal. [37].Inaccuracies in the forward model are likewise highly non-trivial to estimate.While negligible in some cases, they are the dominant error source in others,e.g. NMR chemical shifts, where they can be orders of magnitude greater thanthe statistical error. SAXS forward models share similar problems [16]. Thesource of the model inaccuracy might also come from neglecting an intrinsictime-dependency of the observable (e.g. spin diffusion or dynamic effects in theestimation of NOEs [135]).In summary, as long as errors in the force field, systematic errors in theexperimental data, and errors on the forward model continue to be challenging,if not impossible, to estimate accurately, the parameter θ remains necessary asan effective scaling of the total error σ . A key challenge is therefore to determineit. A simple way is to tune θ until χ r reaches unity (Eq. (8)). This is, however,not generally a valid approach [196]. First, the number of degrees of freedom ν is ill-defined when reweighting is concerned. Conventionally, ν is estimatedas M − k , where M is the number of data points and k is the number of fittedparameters. In MaxEnt, for example, k is given by the number of Lagrangemultipliers, k = M , therefore ν is effectively zero and χ r has consequently nomeaning. Second, although the expectation value of χ r is one, point estimatesare typically different from unity and χ r -distributed. Therefore, this methodmay provide incorrect balance between data and simulation. Another strategyto determine θ is to plot S ( w ) against χ ( w ) and look for an elbow in the curve[32, 37]. Indeed, if plotted with double-logarithmic axis, χ ( S ) often becomes anL-shaped curve, and the optimal θ can be estimated by finding the kink of the40urve [197]. A third strategy is cross-validation, i.e. fitting the optimal ensembleto a separate set of data that has not been used in the analysis [64, 86, 134,198, 199]. The reweighting or experiment-biased simulation may then be donefor several values of θ to monitor when the goodness of fit to the unused datastarts to decrease. While intuitively appealing, the practical implementationmay be difficult. First, the data needs to be divided into independent subsets,which may sometimes be difficult for highly correlated or interdependent data.Second, as different sources of data may report on very different aspects of thesystem, they may in practice not be useful for cross validation [199]. A fourthstrategy for determining θ is strictly Bayesian, and θ is in that context treatedas a nuisance parameter. The posterior probabilities at each value of θ arecalculated and used to find the most probable value of θ [72], or to integrateout θ completely [23]. In section 3.4 we mentioned that reweighting might fail in reproducing theexperimental averages when the prior (force field) used is too inaccurate, becauserelevant states are either poorly sampled or not sampled at all. A pragmatic,and potentially transferable, solution to this problem would be to identify theparameters in the force field responsible for the incorrect behaviour and modifythem slightly in order to move closer to the expected averages. Such rescalingapproaches have proved succesful for specific systems and force fields, e.g. foradjusting the protein-protein interaction strength in the coarse-grained Martiniforce field [200, 201] or protein-water interactions for simulations of disorderedproteins [202]. A subsequent reweighting of the obtained trajectories mightcorrect for inconsistencies that are not related to the identified and re-scaledparameter, and could therefore lead to better consistency with the experimentaldata.Applying such approaches, however, raises some relevant questions. Forexample, to what extent should one correct the force field? Considerable modi-fications to a force field would typically require the need to iteratively re-samplethe system under consideration, an exercise that might become computation-ally expensive. Moreover, after any substantial reparametrisation the force fieldshould be benchmarked against other data and ab initio calculations, as donein the original parametrisation of the force field. Therefore, should one perhapsmodify the force field just enough to allow for reweighting? And, more gener-41lly, should force field corrections and reweighting be employed together at all?Standard force field reparameterization effectively corresponds to reweightinga conformational ensemble, but the extent is limited by the functional formof the energy function, and the simultaneous consideration of data on othersystems. A concurrent application of reweighting and force field optimizationcould violate the Bayesian principle of having well-defined prior and likelihoodin the reweighting process, as the re-scaled prior used for reweigting has alreadybeen adjusted against experimental data. The data is, so to say, used twicein such protocol. We do not have an answer to the aforementioned questions,nonetheless we believe they provide important points of discussion for futureapplications of reweighting, and that Bayesian methods for force field parame-terization would make it easier to merge these two different strategies.
As discussed in section 5.1, Augmented Markov Models are a framework basedon Maximum Entropy and Maximum Likelihood to include experimental infor-mation in a Markov State Model. Olsson et al. [172] used AMMs to predictNMR relaxation dispersion data, by employing only static experimental data toreweight the model. This fact leads to intriguing questions which, to our knowl-edge, have been largely unanswered in literature. As a matter of principle,it should be possible to directly employ kinetic data for ensemble reweight-ing, just like it is done with the ABSURD method [171]. What would happenthen to equilibrium observables? Would the amount of provided information beenough to increase the accuracy of simulated averages with respect to experi-mental equilibrium quantities? What would be the optimal framework to testthis hypothesis? While the first question has no clear answer yet, concerningthe latter we believe that a suitable framework would be the one of MarkovState Models, because of their intrinsic ability to encode kinetic information[172, 191]. We note however that MSMs, by construction, ignore the fastestdynamical timescales of the system, which can however be important for sometypes of experimental measurements. Therefore, a careful choice of kinetic datais recommended.
Force field parametrisation has been historically guided by a combination ofchemical and biological intuition, ab initio quantum calculation of small moleculesand trial and error approaches [108]. This approach has been proven to be ex-42remely successful [203] but, as the ever growing amount of experimental in-formation calls for force fields which are easily improvable once new data arecollected, the complexity and the amount of expertise required in this procedurehave made it somewhat impractical. We discussed some advancements in sec-tion 4.2, but none of the applications we reviewed provided an ultimate solution.Indeed, suggested improvements are usually small adjustments (e.g. modifica-tions in the backbone dihedral angles terms [137, 138]), though in some casesmore extensive changes have been introduced using such fitting to experiments[13, 121, 122, 124]. Therefore, in this section we want to discuss some principleson which a new generation of force fields could be build.One of the main challenges in force field development lies in how to in-clude new data, which is potentially conflicting with old information, and useit to improve the force field after it has already been parametrised. Auto-matic reparametrisation would make it extremely easy to modify and updateforce field once such new experimental data become available. The develop-ment of a Bayesian formalism [116] and later the more systematic ForceBalance[121, 122, 124] framework, for example, goes exactly in this direction and showsthat it is possible in principle to approach to the problem in an automatedfashion. It would be also important to assign some level of trust to force fieldparameters, such that it can be assessed to what degree new data should al-ter their values. At the same time, a fully Bayesian approach would involvedistribution of force fields, rather than point estimates, and thus parallel simu-lations could be used to integrate out force field uncertainty. We also stress thatsuch developments would ideally be carried concurrently with the constructionof worldwide accessible and curated databases of experimental and simulationdata. In order to include new data in a more automatic fashion, there have to bedata quality checks and consensus on experimental and forward model errors,as this all affects how much the new data should be able to alter the existingparametrisation (in case of inconsistency). At the same time, it is also impor-tant to keep in mind that such models should ideally capture well-understoodphysical effects, and that lack of agreement with experiments might indicateimportant effects that are missing from the functional form or parameter com-bining rules [110].The molecular mechanics force fields were developed as a classical parametri-sation of the Born-Oppenheimer energy surface that balances computationalspeed and precision and the choices of functional forms of the force field termshas always been guided by chemical intuition. This has led to families of forcefields with differences in the parametrisation [106] and to the proliferation of43ig sets of parameters needed to accommodate empirical choices (see, for ex-ample, the discussion on atom types in Refs. [204, 205]). We expect the nextgeneration of force fields to be more flexible with respect to specific choices ofparametrisation: more specifically, tools are required to define not only the forcefield parameters, but also the functional shapes of each term in a data-drivenfashion. These developments go hand in hand with a robust estimation of bothstatistical and systematic errors: the two sources of errors need to be fully de-coupled and it should be clear when poor estimates are due to poor trainingdatasets.Some of these advancements have already been applied to the SMIRNOFF99Frosstforce field and the Open Force Field Toolkit [204], developed by the Open ForceField Initiative [206]. For now, SMIRNOFF99Frosst has been tested only on awide set of pharmaceutically relevant small molecules, but it represents nonethe-less an important step towards a new generation of force fields.
Much research in the field of structural biology and molecular biophysics re-quires one to integrate several heterogeneous sources of data, coming both fromdiverse experimental techniques and simulations. Experiments are employed tocharacterise thermodynamic and kinetic quantities of the system under studyand simulations aid the interpretation of or complement these results thanks totheir high spatial and time resolutions. However, technological [17–20] and theo-retical [3, 21] advancements in the field of molecular mechanics have highlightedthe existence of discrepancies between simulations and experiments [11–16]. Inthis review, we have approached this issue through a specific perspective: in-consistencies between computational and experimental results carry informationthat can be systematically extracted and exploited to improve our understand-ing of biochemical entities and more general biophysical models. To pursue thisperspective, we focused on the cornerstone ideas that have guided the optimi-sation of system-specific conformational ensembles (sections 2 and 3) or generalpurpose force fields (section 4) against experimental information. We provideda summary of alternative strategies and discussed their strengths and short-comings depending on the level of trust one places in the force field used tocarry out the simulations. For each method, we critically assessed how mucha realistic scenario deviates from the ideal version of the approach, examiningthe challenges that arise when the systems grow in size and complexity. Westressed that many sources of data can be employed and different strategies are44eeded depending on whether observables can or cannot be represented as alinear combination of a given forward model applied to the single structures inthe ensemble (see Eq. (1) and section 5). Despite the differences in the techni-cal implementation, however, we argued that most of the approaches presentedhere share a common framework, rooted in the minimisation of the functional T = χ − θR, (49)where χ is the negative log-likelihood in the case of normally distributed errors, θ is a hyperparameter that balances between simulations and experimental dataand R is a regularisation term. Depending on the method, R assumes differentfunctional shapes: in MaxEnt, it becomes a cross-entropy term, Eq. (9); inMaxPars, it reduces a parsimony term, i.e. the negative number of conforma-tions with non-zero associated weight, Eq. (13); in MaxPrior, the regularisationterm is provided by the log-prior, Eq. (18); in AMMs R represents the logarithmof the MSM likelihood, Eq. (32); in MaxCal the regularisation term takes theform of a path-entropy, Eq. (37); finally, average block selection is unregularisedand so R = 0, Eq. (48). Therefore, the different philosophies that distinguishamong the several approaches are encoded in the choice of the regularisationterm.As a final remark, we devoted section 6 to the challenges that the communityis facing and we anticipate some open questions that we believe will captureexperts’ attention in the incoming years. Acknowledgements
We acknowledge support by a grant from the Lundbeck Foundation to theBRAINSTRUC structural biology initiative, the NordForsk Nordic Neutron Sci-ence Programme, the Carlsberg Foundation, a grant from the Velux Founda-tions and a Hallas-Møller Stipend from the Novo Nordisk Foundation. We wouldlike to thank members of the Linderstrøm-Lang Centre for Protein Science fornumerous discussions on these topics, and thank Giovanni Bussi, Ramon Cre-huet, Clemens Kauffmann and Simon Olsson for insightful comments on themanuscript.
References [1] P. E. M. Lopes, O. Guvench, and A. D. MacKerell,
Current Status ofProtein Force Fields for Molecular Dynamics Simulations , pp. 47–71. New45ork, NY: Springer New York, 2015.[2] S. Bottaro and K. Lindorff-Larsen, “Biophysical experiments andbiomolecular simulations: A perfect match?,”
Science , vol. 361, no. 6400,pp. 355–360, 2018.[3] T. Maximova, R. Moffatt, B. Ma, R. Nussinov, and A. Shehu, “Principlesand overview of sampling methods for modeling macromolecular structureand dynamics,”
PLoS Computational Biology , vol. 12, no. 4, p. e1004619,2016.[4] S. A. Adcock and J. A. McCammon, “Molecular dynamics: surveyof methods for simulating the activity of proteins,”
Chemical Reviews ,vol. 106, no. 5, pp. 1589–1615, 2006.[5] J. L. Klepeis, K. Lindorff-Larsen, R. O. Dror, and D. E. Shaw, “Long-timescale molecular dynamics simulations of protein structure and func-tion,”
Current Opinion in Structural Biology , vol. 19, no. 2, pp. 120–127,2009.[6] C. Abrams and G. Bussi, “Enhanced sampling in molecular dynamics us-ing metadynamics, replica-exchange, and temperature-acceleration,”
En-tropy , vol. 16, no. 1, pp. 163–199, 2013.[7] R. Elber, “Perspective: Computer simulations of long time dynamics,”
The Journal of Chemical Physics , vol. 144, no. 6, p. 060901, 2016.[8] S. A. Hollingsworth and R. O. Dror, “Molecular dynamics simulation forall,”
Neuron , vol. 99, no. 6, pp. 1129–1143, 2018.[9] C. Maffeo, S. Bhattacharya, J. Yoo, D. Wells, and A. Aksimentiev, “Mod-eling and simulation of ion channels,”
Chemical Reviews , vol. 112, no. 12,pp. 6250–6284, 2012.[10] W. F. van Gunsteren, X. Daura, N. Hansen, A. E. Mark, C. Oosten-brink, S. Riniker, and L. J. Smith, “Validation of molecular simulation:an overview of issues,”
Angewandte Chemie International Edition , vol. 57,no. 4, pp. 884–902, 2018.[11] S. Rauscher, V. Gapsys, M. J. Gajda, M. Zweckstetter, B. L. de Groot, andH. Grubm¨uller, “Structural ensembles of intrinsically disordered proteinsdepend strongly on force field: a comparison to experiment,”
Journal ofChemical Theory and Computation , vol. 11, no. 11, pp. 5513–5524, 2015.4612] J. Henriques, C. Cragnell, and M. Skep¨o, “Molecular dynamics simulationsof intrinsically disordered proteins: force field evaluation and comparisonwith experiment,”
Journal of Chemical Theory and Computation , vol. 11,no. 7, pp. 3420–3431, 2015.[13] P. Robustelli, S. Piana, and D. E. Shaw, “Developing a molecular dynam-ics force field for both folded and disordered protein states,”
Proceedingsof the National Academy of Sciences , vol. 115, no. 21, pp. E4758–E4766,2018.[14] P. S. Nerenberg and T. Head-Gordon, “New developments in force fieldsfor biomolecular simulations,”
Current Opinion in Structural Biology ,vol. 49, pp. 129–138, 2018.[15] R. C. Bernardi, M. C. Melo, and K. Schulten, “Enhanced sampling tech-niques in molecular dynamics simulations of biological systems,”
Biochim-ica et Biophysica Acta (BBA)-General Subjects , vol. 1850, no. 5, pp. 872–877, 2015.[16] T. N. Cordeiro, P. C. Chen, A. De Biasio, N. Sibille, F. J. Blanco, J. S.Hub, R. Crehuet, and P. Bernad´o, “Disentangling polydispersity in thePCNA-p15PAF complex, a disordered, transient and multivalent macro-molecular assembly,”
Nucleic Acids Research , vol. 45, no. 3, pp. 1501–1515, 2017.[17] S. Piana, K. Lindorff-Larsen, and D. E. Shaw, “Atomic-level descriptionof ubiquitin folding,”
Proceedings of the National Academy of Sciences ,vol. 110, no. 15, pp. 5915–5920, 2013.[18] K. Lindorff-Larsen, P. Maragakis, S. Piana, and D. E. Shaw, “Picosecondto millisecond structural dynamics in human ubiquitin,”
The Journal ofPhysical Chemistry B , vol. 120, no. 33, pp. 8313–8320, 2016.[19] V. A. Voelz, M. J¨ager, S. Yao, Y. Chen, L. Zhu, S. A. Waldauer, G. R.Bowman, M. Friedrichs, O. Bakajin, L. J. Lapidus, et al. , “Slow unfolded-state structuring in acyl-coa binding protein folding revealed by simulationand experiment,”
Journal of the American Chemical Society , vol. 134,no. 30, pp. 12565–12577, 2012.[20] G. R. Bowman, V. A. Voelz, and V. S. Pande, “Atomistic folding simu-lations of the five-helix bundle protein λ
6- 85,”
Journal of the AmericanChemical Society , vol. 133, no. 4, pp. 664–667, 2010.4721] C. Camilloni and F. Pietrucci, “Advanced simulation techniques for thethermodynamic and kinetic characterization of biological systems,”
Ad-vances in Physics: X , vol. 3, no. 1, p. 1477531, 2018.[22] H. Wu, F. N¨uske, F. Paul, S. Klus, P. Koltai, and F. No´e, “Varia-tional koopman models: Slow collective variables and molecular kineticsfrom short off-equilibrium simulations,”
The Journal of Chemical Physics ,vol. 146, no. 15, p. 154104, 2017.[23] G. Hummer and J. K¨ofinger, “Bayesian ensemble refinement by replicasimulations and reweighting,”
The Journal of Chemical Physics , vol. 143,p. 243150, 2015.[24] E. Ravera, L. Sgheri, and C. Luchinat, “A critical assessment of methodsto recover information from averaged data,”
Physical Chemistry ChemicalPhysics , vol. 18, pp. 5686–5701, 2016.[25] E. T. Jaynes, “Information Theory and Statistical Mechanics,”
The Phys-ical Review , vol. 106, no. 4, pp. 620–630, 1957.[26] S. Kullback and R. Leibler, “On Information and Sufficiency,”
Annals ofMathematical Statistics , vol. 22, no. 1, pp. 79–86, 1951.[27] S. Hansen and J. Pedersen, “A comparison of three different methods foranalysing small-angle scattering data,”
Journal of Applied Crystallogra-phy , vol. 24, pp. 541–548, 1991.[28] J. Skilling,
Maximum Entropy and Bayesian Methods . Springer Nether-lands, 1989.[29] B. R´ozycki, Y. C. Kim, and G. Hummer, “SAXS ensemble refinement ofESCRT-III CHMP3 conformational transitions,”
Structure , vol. 19, no. 1,pp. 109–116, 2011.[30] H. T. A. Leung, O. Bignucolo, R. Aregger, S. A. Dames, A. Mazur,S. Bern`eche, and S. Grzesiek, “A Rigorous and Efficient Method toReweight Very Large Conformational Ensembles Using Average Exper-imental Data and to Determine Their Relative Information Content,”
Journal of Chemical Theory and Computation , vol. 12, no. 1, pp. 383–394, 2016.[31] K. Reichel, L. S. Stelzl, J. K¨ofinger, and G. Hummer, “Precision DEERDistances from Spin-Label Ensemble Refinement,”
Journal of PhysicalChemistry Letters , vol. 9, no. 19, pp. 5748–5752, 2018.4832] S. Bottaro, G. Bussi, S. D. Kennedy, D. H. Turner, and K. Lindorff-Larsen,“Conformational ensembles of rna oligonucleotides from integrating nmrand molecular simulations,”
Science Advances , vol. 4, no. 5, p. eaar8521,2018.[33] E. T. Jaynes,
Where Do We Stand on Maximum Entropy? In: Papers onProbability, Statistics and Statistical Physics . Springer, Dordrecht, 1978.[34] B. Roux and J. Weare, “On the statistical equivalence of restrained-ensemble simulations with the maximum entropy method,”
The Journalof Chemical Physics , vol. 138, no. 8, p. 084107, 2013.[35] W. Boomsma, J. Ferkinghoff-Borg, and K. Lindorff-Larsen, “CombiningExperiments and Simulations Using the Maximum Entropy Principle,”
PLoS Computational Biology , vol. 10, no. 2, p. e1003406, 2014.[36] A. Cesari, S. Reißer, and G. Bussi, “Using the Maximum Entropy Prin-ciple to Combine Simulations and Solution Experiments,”
Computation ,vol. 6, no. 1, p. 15, 2018.[37] J. K¨ofinger, L. S. Stelzl, K. Reuter, C. Allande, K. Reichel, and G. Hum-mer, “Efficient ensemble refinement by reweighting,”
Journal of ChemicalTheory and Computation , vol. 15, no. 5, pp. 3390–3401, 2019.[38] R. Shevchuk and J. S. Hub, “Bayesian refinement of protein structures andensembles against SAXS data using molecular dynamics,”
PLoS Compu-tational Biology , vol. 13, no. 10, pp. 1–27, 2017.[39] A. Cavalli, M. Bonomi, C. Camilloni, and M. Vendruscolo, “Metainfer-ence: A Bayesian inference method for heterogeneous systems,”
ScienceAdvances , vol. 2, no. 1, pp. e1501177–e1501177, 2016.[40] A. Cesari, A. Gil-Ley, and G. Bussi, “Combining simulations and solu-tion experiments as a paradigm for rna force field refinement,”
Journal ofChemical Theory and Computation , vol. 12, no. 12, pp. 6192–6200, 2016.[41] D. B. Amirkulova and A. D. White, “Recent advances in maximum en-tropy biasing techniques for molecular dynamics,”
ArXiv , p. 1902.02252v1,2019.[42] S. F. Gull and G. J. Daniell, “Image reconstruction from incomplete andnoisy data,”
Nature , vol. 272, pp. 686–690, 1978.4943] H. Wiegand, “Kish, l.: Survey sampling. john wiley & sons, inc., newyork, london 1965, ix + 643 s., 31 abb., 56 tab., preis 83 s.,”
BiometrischeZeitschrift , vol. 10, no. 1, pp. 88–89, 1968.[44] R. Rangan, M. Bonomi, G. T. Heller, A. Cesari, G. Bussi, and M. Ven-druscolo, “Determination of structural ensembles of proteins: restrainingvs reweighting,”
Journal of Chemical Theory and Computation , vol. 14,no. 12, pp. 6632–6641, 2018.[45] H. A. I. Akaike, “A New Look at the Statistical Model Identification,”
IEEE Transaction on Automatic Control , vol. 19, no. 6, pp. 716–723,1974.[46] S. Bowerman, A. S. J. B. Rana, A. Rice, G. H. Pham, E. R. Strieter, andJ. Wereszczynski, “Determining Atomistic SAXS Models of Tri-UbiquitinChains from Bayesian Analysis of Accelerated Molecular Dynamics Simu-lations Samuel,”
Journal of Chemical The , vol. 13, no. 6, pp. 2418–2429,2017.[47] E. Boura, B. R´o˙zycki, D. Z. Herrick, H. S. Chung, J. Vecer, W. A.Eaton, D. S. Cafiso, G. Hummer, and J. H. Hurley, “Solution struc-ture of the escrt-i complex by small-angle x-ray scattering, epr, and fretspectroscopy,”
Proceedings of the National Academy of Sciences , vol. 108,no. 23, pp. 9437–9442, 2011.[48] Y. Chen, S. L. Campbell, and N. V. Dokholyan, “Deciphering proteindynamics from NMR data using explicit structure sampling and selection,”
Biophysical Journal , vol. 93, no. 7, pp. 2300–2306, 2007.[49] D. M. Francis, B. R¨a, D. Koveal, G. Hummer, R. Page, and W. Peti,“Structural basis of p38 regulation by hematopoietic tyrosine phos-phatase,”
Nature Chemical Biology , vol. 7, no. 12, pp. 916–924, 2011.[50] P. Cossio and G. Hummer, “Bayesian analysis of individual electronmicroscopy images: Towards structures of dynamic and heterogeneousbiomolecular assemblies,”
Journal of Structural Biology , vol. 184, no. 3,pp. 427–437, 2013.[51] K. Berlin, C. A. Casta˜neda, D. Schneidman-Duhovny, A. Sali, A. Nava-Tudela, and D. Fushman, “Recovering a representative conformationalensemble from underdetermined macromolecular structural data,”
Journalof the American Chemical Society , vol. 135, no. 44, pp. 16595–16609, 2013.5052] D. Schneidman-Duhovny, M. Hammel, J. A. Tainer, and A. Sali, “FoXS ,FoXSDock and MultiFoXS : Single-state and multi-state structural mod-eling of proteins and their complexes based on SAXS profiles,”
NucleicAcids Research , vol. 44, no. Web server issue, pp. 424–429, 2016.[53] W. Rieping, “Inferential Structure Determination,”
Science , vol. 309,no. 5732, pp. 303–306, 2005.[54] S. Olsson, J. Frellsen, W. Boomsma, K. V. Mardia, and T. Hamelryck,“Inference of structure ensembles of flexible biomolecules from sparse, av-eraged data,”
PLOS ONE , vol. 8, pp. 1–7, 11 2013.[55] R. Dutta, Z. F. Brotzakis, and A. Mira, “Bayesian calibration of force-fields from experimental data: Tip4p water,”
The Journal of ChemicalPhysics , vol. 149, no. 15, p. 154110, 2018.[56] P. Pernot and F. Cailliez, “A critical review of statistical calibra-tion/prediction models handling data inconsistency and model inade-quacy,”
AIChE Journal , vol. 63, no. 10, pp. 4642–4665, 2017.[57] K. A. Beauchamp, V. S. Pande, and R. Das, “Bayesian energy landscapetilting: Towards concordant models of molecular ensembles,”
BiophysicalJournal , vol. 106, no. 6, pp. 1381–1390, 2014.[58] D. H. Brookes and T. Head-Gordon, “Experimental inferential structuredetermination of ensembles for intrinsically disordered proteins,”
Journalof the American Chemical Society , vol. 138, no. 13, pp. 4530–4538, 2016.[59] C. K. Fisher, A. Huang, and C. M. Stultz, “Modeling intrinsically disor-dered proteins with Bayesian statistics,”
Journal of the American Chem-ical Society , vol. 132, no. 42, pp. 14919–14927, 2010.[60] A. Sethi, D. Anunciado, J. Tian, D. M. Vu, and S. Gnanakaran, “Deduc-ing conformational variability of intrinsically disordered proteins from in-frared spectroscopy with Bayesian statistics,”
Chemical Physics , vol. 422,pp. 143–155, 2013.[61] X. Xiao, N. Kallenbach, and Y. Zhang, “Peptide conformation analysisusing an integrated bayesian approach,”
Journal of Chemical Theory andComputation , vol. 10, no. 9, pp. 4152–4159, 2014.[62] Y. Ge and V. A. Voelz, “Model Selection Using BICePs: A BayesianApproach for Force Field Validation and Parameterization,”
Journal ofPhysical Chemistry B , vol. 122, no. 21, pp. 5610–5622, 2018.5163] W. Potrzebowski, J. Trewhella, and I. Andre, “Bayesian inference of pro-tein conformational ensembles from limited structural data,”
PLoS Com-putational Biology , vol. 14, no. 12, p. e1006641, 2018.[64] S. Bottaro, T. Bengtsen, and K. Lindor, “Integrating Molecular Simula-tion and Experimental Data : A Bayesian / Maximum Entropy reweight-ing approach,” bioRxiv , 2018.[65] K. S. Molnar, M. Bonomi, R. Pellarin, G. D. Clinthorne, G. Gonzalez,S. D. Goldberg, M. Goulian, A. Sali, and W. F. Degrado, “Cys-Scanningdisulfide crosslinking and bayesian modeling probe the transmembranesignaling mechanism of the histidine kinase, PhoQ,”
Structure , vol. 22,no. 9, pp. 1239–1251, 2014.[66] M. Mechelke and M. Habeck, “Bayesian weighting of statistical potentialsin nmr structure calculation,”
PloS one , vol. 9, no. 6, p. e100197, 2014.[67] L. D. Antonov, S. Olsson, W. Boomsma, and T. Hamelryck, “Bayesian in-ference of protein ensembles from SAXS data,”
Physical Chemistry Chem-ical Physics , vol. 18, no. 8, pp. 5832–5838, 2016.[68] M. Bonomi, G. T. Heller, C. Camilloni, and M. Vendruscolo, “Princi-ples of protein structural ensemble determination,”
Current Opinion inStructural Biology , vol. 42, pp. 106–116, 2017.[69] C. K. Fisher, O. Ullman, and C. M. Stultz, “Efficient construction of dis-ordered protein ensembles in a Bayesian framework with optimal selectionof conformations,”
Pac Symp Biocomput , pp. 82–93, 2012.[70] S. Kirmizialtin, J. Loerke, E. Behrmann, C. M. Spahn, and K. Y. San-bonmatsu,
Using molecular simulation to model high-resolution cryo-EMreconstructions , vol. 558. Elsevier Inc., 1 ed., 2015.[71] D. J. C. MacKay, “Bayesian model comparison and backprop nets,” in
Advances in Neural Information Processing Systems 4 (J. E. Moody, S. J.Hanson, and R. P. Lippmann, eds.), pp. 839–846, Morgan-Kaufmann,1992.[72] A. H. Larsen, L. Arleth, and S. Hansen, “Analysis of small-angle scatteringdata using model fitting and Bayesian regularization,”
Journal of AppliedCrystallography , vol. 51, no. 4, pp. 1151–1161, 2018.5273] M. Pelikan, G. L. Hura, and M. Hammel, “Structure and flexibility withinproteins as identified through small angle x-ray scattering,”
General Phys-iology and Biophysics , vol. 28, no. 2, p. 174, 2009.[74] C. Vogel, M. Bashton, N. D. Kerrison, C. Chothia, and S. A. Teichmann,“Structure, function and evolution of multidomain proteins,”
CurrentOpinion in Structural Biology , vol. 14, no. 2, pp. 208–216, 2004.[75] H. J. Dyson and P. E. Wright, “Intrinsically unstructured proteins andtheir functions,”
Nature Reviews Molecular Cell Biology , vol. 6, no. 3,pp. 197–208, 2005.[76] K. Lindorff-Larsen, N. Trbovic, P. Maragakis, S. Piana, and D. E. Shaw,“Structure and dynamics of an unfolded protein examined by moleculardynamics simulation,”
Journal of the American Chemical Society , vol. 134,no. 8, pp. 3787–3791, 2012.[77] A. Rodriguez, M. d’Errico, E. Facco, and A. Laio, “Computing the free en-ergy without collective variables,”
Journal of Chemical Theory and Com-putation , vol. 14, no. 3, pp. 1206–1215, 2018.[78] R. F. Alford, A. Leaver-Fay, J. R. Jeliazkov, M. J. O’Meara, F. P. DiMaio,H. Park, M. V. Shapovalov, P. D. Renfrew, V. K. Mulligan, K. Kappel,J. W. Labonte, M. S. Pacella, R. Bonneau, P. Bradley, R. L. Dunbrack,R. Das, D. Baker, B. Kuhlman, T. Kortemme, and J. J. Gray, “TheRosetta All-Atom Energy Function for Macromolecular Modeling and De-sign,”
Journal of Chemical Theory and Computation , vol. 13, pp. 3031–3048, 2017.[79] R. E. Bellman,
Adaptive control processes: a guided tour , vol. 2045. Prince-ton university press, 2015.[80] N. Agmon, Y. Alhassid, and R. D. Levine, “An algorithm for findingthe distribution of maximal entropy,”
Journal of Computational Physics ,vol. 30, no. 2, pp. 250–258, 1979.[81] R. Malouf, “A comparison of algorithms for maximum entropy parame-ter estimation,” in
Proceedings of the 6th conference on Natural languagelearning-Volume 20 , pp. 1–7, Association for Computational Linguistics,2002.[82] P. Bernad´o, E. Mylonas, M. V. Petoukhov, M. Blackledge, and D. I. Sver-gun, “Structural characterization of flexible proteins using small-angle53-ray scattering,”
Journal of the American Chemical Society , vol. 129,no. 17, pp. 5656–5664, 2007.[83] G. Nodet, L. Salmon, V. Ozenne, S. Meier, M. R. Jensen, and M. Black-ledge, “Quantitative description of backbone conformational sampling ofunfolded proteins at amino acid resolution from NMR residual dipolarcouplings,”
Journal of the American Chemical Society , vol. 131, no. 49,pp. 17908–17918, 2009.[84] A. D. White and G. A. Voth, “Efficient and minimal method to bias molec-ular simulations with experimental data,”
Journal of Chemical Theory andComputation , vol. 10, no. 8, pp. 3023–3030, 2014.[85] A. D. White, J. F. Dama, and G. A. Voth, “Designing Free Energy Sur-faces That Match Experimental Data with Metadynamics,”
Journal ofChemical Theory and Computation , vol. 11, no. 6, pp. 2451–2460, 2015.[86] K. Lindorff-Larsen, R. B. Best, M. A. DePristo, C. M. Dobson, andM. Vendruscolo, “Simultaneous determination of protein structure anddynamics,”
Nature , vol. 433, no. 7022, pp. 128–132, 2005.[87] M. Vendruscolo, “Determination of conformationally heterogeneous statesof proteins,”
Current Opinion in Structural Biology , vol. 17, no. 1, pp. 15–20, 2007.[88] B. T. Burnley, P. V. Afonine, P. D. Adams, and P. Gros, “Modellingdynamics in protein crystal structures by ensemble refinement,” eLife ,vol. 1, p. e00311, 2012.[89] E. J. Levin, D. A. Kondrashov, G. E. Wesenberg, and G. N. J. Phillips,“Ensemble refinement of protein crystal structures: validation and appli-cation,”
Structure , vol. 15, no. 9, pp. 1040–1052, 2007.[90] J. W. Pitera and J. D. Chodera, “On the use of experimental observationsto bias simulated ensembles,”
Journal of Chemical Theory and Computa-tion , vol. 8, no. 10, pp. 3445–3451, 2012.[91] A. Cavalli, C. Camilloni, and M. Vendruscolo, “Molecular dynamics simu-lations with replica-averaged structural restraints generate structural en-sembles according to the maximum entropy principle,”
The Journal ofChemical Physics , vol. 138, no. 9, p. 094112, 2013.5492] T. L¨ohr, C. Camilloni, M. Bonomi, and M. Vendruscolo, “A practicalguide to the simultaneous determination of protein structure and dynam-ics using metainference,” arXiv preprint arXiv:1901.08030 , 2019.[93] A. Laio and F. L. Gervasio, “Metadynamics: a method to simulate rareevents and reconstruct the free energy in biophysics, chemistry and ma-terial science,”
Reports on Progress in Physics , vol. 71, no. 12, p. 126601,2008.[94] A. Barducci, M. Bonomi, and M. Parrinello, “Metadynamics,”
Wiley In-terdisciplinary Reviews: Computational Molecular Science , vol. 1, no. 5,pp. 826–843, 2011.[95] L. Sutto, S. Marsili, and F. L. Gervasio, “New advances in metadynamics,”
Wiley Interdisciplinary Reviews: Computational Molecular Science , vol. 2,no. 5, pp. 771–779, 2012.[96] J. Pfaendtner and M. Bonomi, “Efficient sampling of high-dimensionalfree-energy landscapes with parallel bias metadynamics,”
Journal ofChemical Theory and Computation , vol. 11, no. 11, pp. 5062–5067, 2015.[97] M. Bonomi, C. Camilloni, and M. Vendruscolo, “Metadynamic metain-ference: enhanced sampling of the metainference ensemble using metady-namics,”
Scientific Reports , vol. 6, p. 31232, 2016.[98] A. Ianeselli, S. Orioli, G. Spagnolli, P. Faccioli, L. Cupellini, S. Jurinovich,and B. Mennucci, “Atomic detail of protein folding revealed by an abinitio reappraisal of circular dichroism,”
Journal of the American ChemicalSociety , vol. 140, no. 10, pp. 3674–3682, 2018.[99] K. J. Kohlhoff, P. Robustelli, A. Cavalli, X. Salvatella, and M. Vendrus-colo, “Fast and accurate predictions of protein nmr chemical shifts from in-teratomic distances,”
Journal of the American Chemical Society , vol. 131,no. 39, pp. 13894–13895, 2009.[100] Y. Shen and A. Bax, “Sparta+: a modest improvement in empirical nmrchemical shift prediction by means of an artificial neural network,”
Journalof Biomolecular NMR , vol. 48, no. 1, pp. 13–22, 2010.[101] W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, “Development andtesting of the opls all-atom force field on conformational energetics andproperties of organic liquids,”
Journal of the American Chemical Society ,vol. 118, no. 45, pp. 11225–11236, 1996.55102] C. Oostenbrink, A. Villa, A. E. Mark, and W. F. Van Gunsteren, “Abiomolecular force field based on the free enthalpy of hydration and sol-vation: the gromos force-field parameter sets 53a5 and 53a6,”
Journal ofComputational Chemistry , vol. 25, no. 13, pp. 1656–1676, 2004.[103] A. D. MacKerell Jr, D. Bashford, M. Bellott, R. L. Dunbrack Jr, J. D.Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, et al. , “All-atom empirical potential for molecular modeling and dynamics studies ofproteins,”
The Journal of Physical Chemistry B , vol. 102, no. 18, pp. 3586–3616, 1998.[104] W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M.Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, and P. A. Kollman,“A second generation force field for the simulation of proteins, nucleicacids, and organic molecules,”
Journal of the American Chemical Society ,vol. 117, no. 19, pp. 5179–5197, 1995.[105] M. Tuckerman,
Statistical mechanics: theory and molecular simulation .Oxford University Press, 2010.[106] O. Guvench and A. D. MacKerell,
Comparison of Protein Force Fields forMolecular Dynamics Simulations , pp. 63–88. Totowa, NJ: Humana Press,2008.[107] J. W. Ponder and D. A. Case, “Force fields for protein simulations,” in
Advances in Protein Chemistry , vol. 66, pp. 27–85, Elsevier, 2003.[108] X. Zhu, P. E. Lopes, and A. D. MacKerell Jr, “Recent developments andapplications of the charmm force fields,”
Wiley Interdisciplinary Reviews:Computational Molecular Science , vol. 2, no. 1, pp. 167–185, 2012.[109] P. Dauber-Osguthorpe and A. T. Hagler, “Biomolecular force fields: wherehave we been, where are we now, where do we need to go and how do weget there?,”
Journal of Computer-aided Molecular Design , vol. 33, no. 2,pp. 133–203, 2019.[110] A. T. Hagler, “Force field development phase ii: Relaxation of physics-based criteria. . . or inclusion of more rigorous physics into the represen-tation of molecular energetics,”
Journal of Computer-aided Molecular De-sign , vol. 33, no. 2, pp. 205–264, 2019.[111] C. I. Bayly, P. Cieplak, W. Cornell, and P. A. Kollman, “A well-behavedelectrostatic potential based method using charge restraints for deriving56tomic charges: the resp model,”
The Journal of Physical Chemistry ,vol. 97, no. 40, pp. 10269–10280, 1993.[112] A. Jakalian, B. L. Bush, D. B. Jack, and C. I. Bayly, “Fast, efficientgeneration of high-quality atomic charges. am1-bcc model: I. method,”
Journal of Computational Chemistry , vol. 21, no. 2, pp. 132–146, 2000.[113] A. Jakalian, D. B. Jack, and C. I. Bayly, “Fast, efficient generation of high-quality atomic charges. am1-bcc model: Ii. parameterization and valida-tion,”
Journal of Computational Chemistry , vol. 23, no. 16, pp. 1623–1641,2002.[114] K. Lindorff-Larsen, S. Piana, K. Palmo, P. Maragakis, J. L. Klepeis, R. O.Dror, and D. E. Shaw, “Improved side-chain torsion potentials for the am-ber ff99sb protein force field,”
Proteins: Structure, Function, and Bioin-formatics , vol. 78, no. 8, pp. 1950–1958, 2010.[115] C. Tian, K. Kasavajhala, K. Belfon, L. Raguette, H. Huang, A. Migues,J. Bickel, Y. Wang, J. Pincay, Q. Wu, et al. , “ff19sb: Amino-acid specificprotein backbone parameters trained against quantum mechanics energysurfaces in solution,”
ChemRxiv , 2019.[116] A. B. Norgaard, J. Ferkinghoff-Borg, and K. Lindorff-Larsen, “Experimen-tal parameterization of an energy function for the simulation of unfoldedproteins,”
Biophysical Journal , vol. 94, no. 1, pp. 182–192, 2008.[117] R. B. Best and G. Hummer, “Optimized molecular dynamics force fieldsapplied to the helix- coil transition of polypeptides,”
The Journal of Phys-ical Chemistry B , vol. 113, no. 26, pp. 9004–9015, 2009.[118] F. Cailliez and P. Pernot, “Statistical approaches to forcefield calibra-tion and prediction uncertainty in molecular simulation,”
The Journal ofChemical Physics , vol. 134, no. 5, p. 054124, 2011.[119] F. Rizzi, H. N. Najm, B. J. Debusschere, K. Sargsyan, M. Salloum,H. Adalsteinsson, and O. M. Knio, “Uncertainty quantification in md sim-ulations. part ii: Bayesian inference of force-field parameters,”
MultiscaleModeling & Simulation , vol. 10, no. 4, pp. 1460–1492, 2012.[120] P. Angelikopoulos, C. Papadimitriou, and P. Koumoutsakos, “Bayesianuncertainty quantification and propagation in molecular dynamics simu-lations: a high performance computing framework,”
The Journal of Chem-ical Physics , vol. 137, no. 14, p. 144103, 2012.57121] L.-P. Wang, J. Chen, and T. Van Voorhis, “Systematic parametrization ofpolarizable force fields from quantum chemistry data,”
Journal of Chem-ical Theory and Computation , vol. 9, no. 1, pp. 452–460, 2012.[122] L.-P. Wang, T. J. Martinez, and V. S. Pande, “Building force fields: Anautomatic, systematic, and reproducible approach,”
The Journal of Phys-ical Chemistry Letters , vol. 5, no. 11, pp. 1885–1891, 2014.[123] S. Wu, P. Angelikopoulos, C. Papadimitriou, R. Moser, and P. Koumout-sakos, “A hierarchical bayesian framework for force field selection inmolecular dynamics simulations,”
Philosophical Transactions of the RoyalSociety A: Mathematical, Physical and Engineering Sciences , vol. 374,no. 2060, p. 20150032, 2016.[124] L.-P. Wang, K. A. McKiernan, J. Gomes, K. A. Beauchamp, T. Head-Gordon, J. E. Rice, W. C. Swope, T. J. Mart´ınez, and V. S. Pande, “Build-ing a more predictive protein force field: a systematic and reproducibleroute to amber-fb15,”
The Journal of Physical Chemistry B , vol. 121,no. 16, pp. 4023–4039, 2017.[125] A. M. Bonvin, R. Boelens, and R. Kaptein, “Time-and ensemble-averageddirect noe restraints,”
Journal of Biomolecular NMR , vol. 4, no. 1,pp. 143–149, 1994.[126] J.-r. Huang and S. Grzesiek, “Ensemble calculations of unstructured pro-teins constrained by rdc and pre data: a case study of urea-denaturedubiquitin,”
Journal of the American Chemical Society , vol. 132, no. 2,pp. 694–705, 2009.[127] O. F. Lange, N.-A. Lakomek, C. Far`es, G. F. Schr¨oder, K. F. Walter,S. Becker, J. Meiler, H. Grubm¨uller, C. Griesinger, and B. L. De Groot,“Recognition dynamics up to microseconds revealed from an rdc-derivedubiquitin ensemble in solution,”
Science , vol. 320, no. 5882, pp. 1471–1475,2008.[128] S. Olsson, B. R. V¨ogeli, A. Cavalli, W. Boomsma, J. Ferkinghoff-Borg,K. Lindorff-Larsen, and T. Hamelryck, “Probabilistic determination ofnative state ensembles of proteins,”
Journal of Chemical Theory and Com-putation , vol. 10, no. 8, pp. 3484–3491, 2014.[129] S. Esteban-Mart´ın, R. B. Fenwick, and X. Salvatella, “Refinement ofensembles describing unstructured proteins using nmr residual dipolar58ouplings,”
Journal of the American Chemical Society , vol. 132, no. 13,pp. 4626–4632, 2010.[130] R. Scheek, A. Torda, J. Kemmink, and W. Van Gunsteren, “Structuredetermination by nmr: The modeling of nmr parameters as ensemble av-erages,” in
Computational aspects of the study of biological macromoleculesby nuclear magnetic resonance spectroscopy , pp. 209–217, Springer, 1991.[131] A. B. Mantsyzov, A. S. Maltsev, J. Ying, Y. Shen, G. Hummer, andA. Bax, “A maximum entropy approach to the study of residue-specificbackbone angle distributions in α -synuclein, an intrinsically disorderedprotein,” Protein Science , vol. 23, no. 9, pp. 1275–1290, 2014.[132] A. B. Mantsyzov, Y. Shen, J. H. Lee, G. Hummer, and A. Bax, “Mera: awebserver for evaluating backbone torsion angle distributions in dynamicand disordered proteins from nmr data,”
Journal of Biomolecular NMR ,vol. 63, no. 1, pp. 85–95, 2015.[133] J. Graf, P. H. Nguyen, G. Stock, and H. Schwalbe, “Structure and dy-namics of the homologous series of alanine peptides: a joint molecular dy-namics/nmr study,”
Journal of the American Chemical Society , vol. 129,no. 5, pp. 1179–1189, 2007.[134] A. Cesari, S. Bottaro, K. Lindorff-Larsen, P. Ban´aˇs, J. Sponer, andG. Bussi, “Fitting corrections to an rna force field using experimentaldata,”
Journal of Chemical Theory and Computation , 2019.[135] F. Vasile and G. Tiana, “Determination of structural ensembles of flexiblemolecules in solution from nmr data undergoing spin diffusion,”
Journalof Chemical Information and Modeling , 2019.[136] R. W. Zwanzig, “High-temperature equation of state by a perturbationmethod. i. nonpolar gases,”
The Journal of Chemical Physics , vol. 22,no. 8, pp. 1420–1426, 1954.[137] D.-W. Li and R. Br¨uschweiler, “Nmr-based protein potentials,”
Ange-wandte Chemie International Edition , vol. 49, no. 38, pp. 6778–6780, 2010.[138] D.-W. Li and R. Br¨uschweiler, “Iterative optimization of molecular me-chanics force fields from nmr data of full-length proteins,”
Journal ofChemical Theory and Computation , vol. 7, no. 6, pp. 1773–1782, 2011.59139] J. Chen, J. Chen, G. Pinamonti, and C. Clementi, “Learning effectivemolecular models from experimental observables,”
Journal of ChemicalTheory and Computation , vol. 14, no. 7, pp. 3849–3858, 2018.[140] S. A. Teukolsky, B. P. Flannery, W. Press, and W. Vetterling, “Numericalrecipes in c,”
SMR , vol. 693, no. 1, pp. 59–70, 1992.[141] R. Br¨uschweiler, “Collective protein dynamics and nuclear spin relax-ation,”
The Journal of Chemical Physics , vol. 102, no. 8, pp. 3396–3403,1995.[142] C. H. Bennett, “Efficient estimation of free energy differences from montecarlo data,”
Journal of Computational Physics , vol. 22, no. 2, pp. 245–268,1976.[143] M. R. Shirts and J. D. Chodera, “Statistically optimal analysis of sam-ples from multiple equilibrium states,”
The Journal of Chemical Physics ,vol. 129, no. 12, p. 124105, 2008.[144] J. R. Gillespie and D. Shortle, “Characterization of long-range structurein the denatured state of staphylococcal nuclease. i. paramagnetic relax-ation enhancement by nitroxide spin labels,”
Journal of Molecular Biology ,vol. 268, no. 1, pp. 158–169, 1997.[145] M. A. Lietzow, M. Jamin, H. J. Dyson, and P. E. Wright, “Mappinglong-range contacts in a highly unfolded protein,”
Journal of MolecularBiology , vol. 322, no. 4, pp. 655–662, 2002.[146] Q. Yi, M. L. Scalley-Kim, E. J. Alm, and D. Baker, “Nmr characteriza-tion of residual structure in the denatured state of protein l,”
Journal ofMolecular Biology , vol. 299, no. 5, pp. 1341–1351, 2000.[147] K. Teilum, B. B. Kragelund, and F. M. Poulsen, “Transient structureformation in unfolded acyl-coenzyme a-binding protein observed by site-directed spin labelling,”
Journal of Molecular Biology , vol. 324, no. 2,pp. 349–357, 2002.[148] J. R. Gillespie and D. Shortle, “Characterization of long-range structure inthe denatured state of staphylococcal nuclease. ii. distance restraints fromparamagnetic relaxation and calculation of an ensemble of structures,”
Journal of Molecular Biology , vol. 268, no. 1, pp. 170–184, 1997.60149] V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg, and C. Sim-merling, “Comparison of multiple amber force fields and development ofimproved protein backbone parameters,”
Proteins: Structure, Function,and Bioinformatics , vol. 65, no. 3, pp. 712–725, 2006.[150] D. E. Condon, S. D. Kennedy, B. C. Mort, R. Kierzek, I. Yildirim, andD. H. Turner, “Stacking in rna: Nmr of four tetramers benchmark molec-ular dynamics,”
Journal of Chemical Theory and Computation , vol. 11,no. 6, pp. 2729–2742, 2015.[151] C. Bergonzo, N. M. Henriksen, D. R. Roe, and T. E. Cheatham, “Highlysampled tetranucleotide and tetraloop motifs enable evaluation of commonrna force fields,”
RNA , vol. 21, no. 9, pp. 1578–1590, 2015.[152] P. Kuhrova, R. B. Best, S. Bottaro, G. Bussi, J. Sponer, M. Otyepka, andP. Banas, “Computer folding of rna tetraloops: identification of key forcefield deficiencies,”
Journal of Chemical Theory and Computation , vol. 12,no. 9, pp. 4534–4548, 2016.[153] S. Bottaro, P. Ban´aˇs, J. Sponer, and G. Bussi, “Free energy landscape ofgaga and uucg rna tetraloops,”
The Journal of Physical Chemistry Letters ,vol. 7, no. 20, pp. 4032–4038, 2016.[154] A. N. Borkar, M. F. Bardaro, C. Camilloni, F. A. Aprile, G. Varani,and M. Vendruscolo, “Structure of a low-population binding intermedi-ate in protein-rna recognition,”
Proceedings of the National Academy ofSciences , vol. 113, no. 26, pp. 7171–7176, 2016.[155] M. Krepl, M. Blatter, A. Cl´ery, F. F. Damberger, F. H. Allain, andJ. Sponer, “Structural study of the fox-1 rrm protein hydration revealsa role for key water molecules in rrm-rna recognition,”
Nucleic Acids Re-search , vol. 45, no. 13, pp. 8046–8063, 2017.[156] P. Podbevˇsek, F. Fasolo, C. Bon, L. Cimatti, S. Reißer, P. Carninci,G. Bussi, S. Zucchelli, J. Plavec, and S. Gustincich, “Structural deter-minants of the sine b2 element embedded in the long non-coding rnaactivator of translation as uchl1,”
Scientific Reports , vol. 8, no. 1, p. 3189,2018.[157] H. Kooshapur, N. R. Choudhury, B. Simon, M. M¨uhlbauer, A. Jussupow,N. Fernandez, A. N. Jones, A. Dallmann, F. Gabel, C. Camilloni, et al. ,61Structural basis for terminal loop recognition and stimulation of pri-mirna-18a processing by hnrnp a1,”
Nature Communications , vol. 9, no. 1,p. 2479, 2018.[158] A. H. Aytenfisu, A. Spasic, A. Grossfield, H. A. Stern, and D. H. Math-ews, “Revised rna dihedral parameters for the amber force field improverna molecular dynamics,”
Journal of Chemical Theory and Computation ,vol. 13, no. 2, pp. 900–915, 2017.[159] D. Tan, S. Piana, R. M. Dirks, and D. E. Shaw, “Rna force field withaccuracy comparable to state-of-the-art protein force fields,”
Proceedingsof the National Academy of Sciences , vol. 115, no. 7, pp. E1346–E1355,2018.[160] P. Kuhrova, V. Mlynsky, M. Zgarbova, M. Krepl, G. Bussi, R. B. Best,M. Otyepka, J. Sponer, and P. Banas, “Improving the performance ofthe amber rna force field by tuning the hydrogen-bonding interactions,”
Journal of Chemical Theory and Computation , vol. 15, no. 5, pp. 3288–3305, 2019.[161] A. P´erez, I. March´an, D. Svozil, J. Sponer, T. E. Cheatham III, C. A.Laughton, and M. Orozco, “Refinement of the amber force field for nucleicacids: improving the description of α / γ conformers,” Biophysical Journal ,vol. 92, no. 11, pp. 3817–3829, 2007.[162] M. Zgarbov´a, M. Otyepka, J. ˇSponer, A. Ml´adek, P. Ban´aˇs, T. E.Cheatham III, and P. Jurecka, “Refinement of the cornell et al. nucleicacids force field based on reference quantum chemical calculations of gly-cosidic torsion profiles,”
Journal of Chemical Theory and Computation ,vol. 7, no. 9, pp. 2886–2902, 2011.[163] S. Izadi, R. Anandakrishnan, and A. V. Onufriev, “Building water models:a different approach,”
The Journal of Physical Chemistry Letters , vol. 5,no. 21, pp. 3863–3871, 2014.[164] V. S. Pande, K. Beauchamp, and G. R. Bowman, “Everything you wantedto know about markov state models but were afraid to ask,”
Methods ,vol. 52, no. 1, pp. 99–105, 2010.[165] G. R. Bowman, V. S. Pande, and F. No´e,
An introduction to Markovstate models and their application to long timescale molecular simulation ,vol. 797. Springer Science & Business Media, 2013.62166] B. E. Husic and V. S. Pande, “Markov state models: From an art toa science,”
Journal of the American Chemical Society , vol. 140, no. 7,pp. 2386–2396, 2018.[167] A. G. Palmer III, “Nmr characterization of the dynamics of biomacro-molecules,”
Chemical Reviews , vol. 104, no. 8, pp. 3623–3640, 2004.[168] J.-H. Prinz, H. Wu, M. Sarich, B. Keller, M. Senne, M. Held, J. D.Chodera, C. Sch¨utte, and F. No´e, “Markov models of molecular kinetics:Generation and validation,”
The Journal of Chemical Physics , vol. 134,no. 17, p. 174105, 2011.[169] C. R. Schwantes, R. T. McGibbon, and V. S. Pande, “Perspective: Markovmodels for long-timescale biomolecular dynamics,”
The Journal of Chem-ical Physics , vol. 141, no. 9, p. 09B201 1, 2014.[170] E. T. Jaynes, “The minimum entropy production principle,”
Annual Re-view of Physical Chemistry , vol. 31, no. 1, pp. 579–601, 1980.[171] N. Salvi, A. Abyzov, and M. Blackledge, “Multi-timescale dynamics in in-trinsically disordered proteins from nmr relaxation and molecular simula-tion,”
The Journal of Physical Chemistry Letters , vol. 7, no. 13, pp. 2483–2489, 2016.[172] S. Olsson, H. Wu, F. Paul, C. Clementi, and F. No´e, “Combining experi-mental and simulation data of molecular processes via augmented markovmodels,”
Proceedings of the National Academy of Sciences , vol. 114, no. 31,pp. 8265–8270, 2017.[173] G. R. Bowman, K. A. Beauchamp, G. Boxer, and V. S. Pande, “Progressand challenges in the automated construction of markov state models forfull protein systems,”
The Journal of Chemical Physics , vol. 131, no. 12,p. 124101, 2009.[174] M. M. Sultan, G. Kiss, D. Shukla, and V. S. Pande, “Automatic selectionof order parameters in the analysis of large scale molecular dynamics sim-ulations,”
Journal of Chemical Theory and Computation , vol. 10, no. 12,pp. 5217–5223, 2014.[175] G. P´erez-Hern´andez, F. Paul, T. Giorgino, G. De Fabritiis, and F. No´e,“Identification of slow molecular order parameters for markov model con-struction,”
The Journal of Chemical Physics , vol. 139, no. 1, p. 07B604 1,2013. 63176] C. R. Schwantes and V. S. Pande, “Improvements in markov state modelconstruction reveal many non-native interactions in the folding of ntl9,”
Journal of Chemical Theory and Computation , vol. 9, no. 4, pp. 2000–2009, 2013.[177] F. No´e and C. Clementi, “Kinetic distance and kinetic maps from molecu-lar dynamics simulation,”
Journal of Chemical Theory and Computation ,vol. 11, no. 10, pp. 5002–5011, 2015.[178] M. K. Scherer, B. Trendelkamp-Schroer, F. Paul, G. P´erez-Hern´andez,M. Hoffmann, N. Plattner, C. Wehmeyer, J.-H. Prinz, and F. No´e,“Pyemma 2: A software package for estimation, validation, and analy-sis of markov models,”
Journal of Chemical Theory and Computation ,vol. 11, no. 11, pp. 5525–5542, 2015.[179] M. P. Harrigan, M. M. Sultan, C. X. Hern´andez, B. E. Husic, P. Eastman,C. R. Schwantes, K. A. Beauchamp, R. T. McGibbon, and V. S. Pande,“Msmbuilder: statistical models for biomolecular dynamics,”
BiophysicalJournal , vol. 112, no. 1, pp. 10–15, 2017.[180] F. No´e and E. Rosta, “Markov models of molecular kinetics,”
The Journalof Chemical Physics , vol. MMMK, p. 190401, 2019.[181] Y. Xue, J. M. Ward, T. Yuwen, I. S. Podkorytov, and N. R. Skryn-nikov, “Microsecond time-scale conformational exchange in proteins: us-ing long molecular dynamics trajectory to simulate nmr relaxation dis-persion data,”
Journal of the American Chemical Society , vol. 134, no. 5,pp. 2555–2562, 2012.[182] S. Olsson and F. No’e, “Mechanistic models of chemical exchange inducedrelaxation in protein nmr,”
Journal of the American Chemical Society ,vol. 139, no. 1, pp. 200–210, 2016.[183] S. Press´e, K. Ghosh, J. Lee, and K. A. Dill, “Principles of maximumentropy and maximum caliber in statistical physics,”
Reviews of ModernPhysics , vol. 85, no. 3, p. 1115, 2013.[184] P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, “Transitionpath sampling: Throwing ropes over rough mountain passes, in the dark,”
Annual Review of Physical Chemistry , vol. 53, no. 1, pp. 291–318, 2002.64185] P. Eastman, N. Grønbech-Jensen, and S. Doniach, “Simulation of proteinfolding by reaction path annealing,”
The Journal of Chemical Physics ,vol. 114, no. 8, pp. 3823–3841, 2001.[186] A. C. Pan, D. Sezer, and B. Roux, “Finding transition pathways usingthe string method with swarms of trajectories,”
The Journal of PhysicalChemistry B , vol. 112, no. 11, pp. 3432–3440, 2008.[187] E. Weinan, W. Ren, and E. Vanden-Eijnden, “String method for the studyof rare events,”
Physical Review B , vol. 66, no. 5, p. 052301, 2002.[188] R. Capelli, G. Tiana, and C. Camilloni, “An implementation of themaximum-caliber principle by replica-averaged time-resolved restrainedsimulations,”
The Journal of Chemical Physics , vol. 148, no. 18, p. 184114,2018.[189] B. Roux and J. Weare, “On the statistical equivalence of restrained-ensemble simulations with the maximum entropy method,”
The Journalof Chemical Physics , vol. 138, no. 8, p. 02B616, 2013.[190] A. Cavalli, C. Camilloni, and M. Vendruscolo, “Molecular dynamics simu-lations with replica-averaged structural restraints generate structural en-sembles according to the maximum entropy principle,”
The Journal ofChemical Physics , vol. 138, no. 9, p. 03B603, 2013.[191] P. D. Dixit and K. A. Dill, “Caliber corrected markov modeling (c2m2):Correcting equilibrium markov models,”
Journal of Chemical Theory andComputation , vol. 14, no. 2, pp. 1111–1119, 2018.[192] P. D. Dixit, “Communication: Introducing prescribed biases in out-of-equilibrium markov models,”
The Journal of Chemical Physics , vol. 148,no. 9, p. 091101, 2018.[193] M. Bause, T. Wittenstein, K. Kremer, and T. Bereau, “Microscopicreweighting for non-equilibrium steady states dynamics,” arXiv preprintarXiv:1907.08480 , 2019.[194] N. Salvi, A. Abyzov, and M. Blackledge, “Solvent-dependent segmentaldynamics in intrinsically disordered proteins,”
Science Advances , vol. 5,no. 6, p. eaax2348, 2019.[195] R. Shevchuk and J. S. Hub, “Bayesian refinement of protein structuresand ensembles against saxs data using molecular dynamics,”
PLoS Com-putational Biology , vol. 13, no. 10, p. e1005800, 2017.65196] R. Andrae, T. Schulze-Hartung, and P. Melchior, “Dos and don’ts of re-duced chi-squared,” arXiv , p. 1012.374v1, 2010.[197] P. C. Hansen, “The l-curve and its use in the numerical treatment of in-verse problems,” in
Computational Inverse Problems in Electrocardiology ,WIT Press, 2000.[198] P.-c. Chen, R. Shevchuk, F. M. Strnad, C. Lorenz, L. Karge, R. Gilles,A. M. Stadler, J. Hennig, and J. S. Hub, “Combined small-angle x-ray andneutron scattering restraints in molecular dynamics simulations,”
Journalof Chemical Theory and Computation , vol. 0, p. 0, 0.[199] R. Crehuet, P. J. B. Jorro, K. Lindorff-Larsen, and X. Salvatella,“Bayesian-maximum-entropy reweighting of idps ensembles based on nmrchemical shifts,” bioRxiv , p. 689083, 2019.[200] A. C. Stark, C. T. Andrews, and A. H. Elcock, “Toward optimized po-tential functions for protein-protein interactions in aqueous solutions: Os-motic second virial coefficient calculations using the MARTINI coarse-grained force field,”
Journal of Chemical Theory and Computation , vol. 9,no. 9, pp. 4176–4185, 2013.[201] M. Javanainen, H. Martinez-Seara, and I. Vattulainen, “Excessive aggre-gation of membrane proteins in the martini model,”
PLOS ONE , vol. 12,pp. 1–20, 11 2017.[202] R. B. Best, W. Zheng, and J. Mittal, “Balanced protein–water interac-tions improve properties of disordered proteins and non-specific proteinassociation,”
Journal of chemical theory and computation , vol. 10, no. 11,pp. 5113–5124, 2014.[203] A. Hospital, J. R. Go˜ni, M. Orozco, and J. L. Gelp´ı, “Molecular dynamicssimulations: advances and applications,”
Advances and Applications inBioinformatics and Chemistry: AABC , vol. 8, p. 37, 2015.[204] D. L. Mobley, C. C. Bannan, A. Rizzi, C. I. Bayly, J. D. Chodera, V. T.Lim, N. M. Lim, K. A. Beauchamp, D. R. Slochower, M. R. Shirts, et al. ,“Escaping atom types in force fields using direct chemical perception,”
Journal of Chemical Theory and Computation , vol. 14, no. 11, pp. 6076–6092, 2018.[205] C. Zanette, C. C. Bannan, C. I. Bayly, J. Fass, M. K. Gilson, M. R. Shirts,J. D. Chodera, and D. L. Mobley, “Toward learned chemical perception66f force field typing rules,”