Empirical optimization of molecular simulation force fields by Bayesian inference

The demands on the accuracy of force fields for classical molecular dynamics simulations are steadily growing as larger and more complex systems are studied over longer times. One way to meet these growing demands is to hand over the learning of force fields and their parameters to machines in a systematic (semi)automatic manner. Doing so, we can take full advantage of exascale computing, the increasing availability of experimental data, and advances in quantum mechanical computations and the calculation of experimental observables from molecular ensembles. Here, we discuss and illustrate the challenges one faces in this endeavor and explore a way forward by adapting the Bayesian inference of ensembles (BioEn) method [Hummer and Köfinger, J. Chem. Phys. (2015)] for force field parameterization. In the Bayesian inference of force fields (BioFF) method developed here, the optimization problem is regularized by a simplified prior on the force field parameters and an entropic prior acting on the ensemble. The latter compensates for the unavoidable over simplifications in the parameter prior. We determine optimal force field parameters using an iterative predictor–corrector approach, in which we run simulations, determine the reference ensemble using the weighted histogram analysis method (WHAM), and update the force field according to the BioFF posterior. We illustrate this approach for a simple polymer model, using the distance between two labeled sites as the experimental observable. By systematically resolving force field issues, instead of just reweighting a structural ensemble, the BioFF corrections extend to observables not included in ensemble reweighting. We envision future force field optimization as a formalized, systematic, and (semi)automatic machine-learning effort that incorporates a wide range of data from experiment and high-level quantum chemical calculations, and takes advantage of exascale computing resources.


Introduction
The predictive power of molecular dynamics (MD) simulations relies on the accuracy and the computational efficiency of the applied force fields describing the molecular interactions. Widely used semi-empirical force fields aim to strike a balance between accuracy and efficiency [1,2]. The latter restricts their computational complexity. Force fields are parameterized in part by quantum mechanical calculations and in part by fitting calculated observables to experimental data. For decades, molecular force fields have been adapted and re-parameterized to fit an ever-growing resource of theoretical and experimental information [3]. These cumulative efforts have dramatically increased the predictive power of MD simulations [4].
The advent of exascale computing makes it possible -and ultimately necessary -to hand over the task of force field optimization to machines in a formala e-mail: juergen.koefinger@biophys.mpg.de (corresponding author) b e-mail: gerhard.hummer@biophys.mpg.de ized, systematic, and (semi)automatic manner [5][6][7][8][9][10][11]. However, the complexity of molecular force fields, their high-dimensional parameter spaces, and the interdependence of the parameters pose difficult challenges for automated force field optimization. Whereas this problem is well recognized for classical mechanical energy functions, a quite similar problem arises in quantum mechanical modeling. The latter also involves multiple model choices (say, which level of description or which basis set to use) or extensive parameterizations (e.g., in density functional theory). Intrinsically, force field calibration is a problem of statistical inference, where we update our current understanding of molecular interactions, e.g., from quantum mechanical calculations for molecular fragments, with uncertain information from various sources [9,12]. Quantum mechanical calculations, experimental measurements, MD simulations, and the calculations of experimental observables from molecular models are all subject to systematic or sampling errors. In a Bayesian formulation, a prior distribution on the force field parameters encodes our current state of knowledge which is then updated with new experimental data through the likelihood. Unfortunately, it is difficult to define such priors. The complex and often strong interdependence of force field parameters is hard to quantify. Consequently, without good priors, even small changes in the parameters can lead to highly non-physical conformations in Boltzmann sampling.
We illustrate these challenges and explore the use of maximum entropy [13][14][15][16] and Bayesian [17][18][19] ensemble refinement methods to compensate for the lack of good priors on the force field parameters. In ensemble refinement, we use experimental data to refine the statistical weights of an ensemble of structures. Our starting point is the Bayesian inference of ensembles (BioEn) method [18]. In BioEn, an entropic prior acts on the statistical weights of the structures in an ensemble such that refinement changes the weights only minimally to better fit the experimental data. Properties of the system that are independent of the data remain untouched.
Here, we extend BioEn by expressing the statistical weights as functions on selected force field parameters and by adding a prior on these parameters. Unavoidable deficiencies in the parameter prior are compensated by the entropic prior on the ensemble of structures. In this Bayesian inference of force fields (BioFF) method, we iteratively optimize the resulting BioFF posterior by incorporating the results of MD simulations with original and newly optimized force field parameters to ensure convergence. To illustrate the BioFF method, we infer the value of the single force field parameter of a simple polymer model using a single label distance as experimental data.
The article is organized as follows. In Sect. 2 we briefly review ensemble refinement in the context of force field parameterization and introduce the BioFF method. We present method details in Sect. 3 and a simple proof of principle for a two-dimensional polymer in Sect. 4. We discuss the results in a broader context in Sect. 5 and present a tentative outlook on the future of force field optimization and conclusions in Sect. 6.

Theory
We briefly introduce ensemble refinement methods and review their use for force field optimization before presenting the BioFF method. We do so step-by-step: We first present the BioFF posterior in its general form as a probability density function of the force field parameters, which we then adapt for discrete ensembles. For optimization, we reweight the discrete weights as function of the force field parameters. The key to the accuracy and efficiency of BioFF by reweighting is to pool all MD simulations generated for different force fields and calculate new reference weights using the binless weighted histogram analysis method (WHAM) [20][21][22]. Equivalently, one can use the multistate Bennett acceptance ratio (MBAR) method [23,24]. Binless WHAM has to be executed only once for each new simulation added to the pool. Note that pooling all simulations for different force field parameters using binless WHAM also distinguishes BioFF from previous methods pursuing similar goals [25][26][27]. Pooling and priors stabilize the optimization procedure and minimize the risk of sampling non-physical conformations. We sketch the BioFF-by-reweighting procedure at the end of the section.

Background and preliminaries
In ensemble refinement, we adapt the statistical weights of the structures in an ensemble to obtain better agreement with experimentally measured ensemble averages [15,18,28]. By properly accounting for uncertainties and balancing the information encoded in the structural ensemble and in the experimental data, the resulting ensemble is a better representation of the true ensemble underlying the experiment. The structural ensemble is usually generated in free or restrained MD simulations. The refinement of such simulation ensembles by reweighting ensures convergence with ensemble size [18].
Many of the current approaches to ensemble refinement that appear to be different on the surface are actually quite similar and will in principle provide nearly the same optimal ensembles if we encode the same prior knowledge [13,[17][18][19][29][30][31][32]. In the following, we use the formalism of the BioEn method [18] to introduce ensemble refinement. We then adapt BioEn for force field parameterization.
Consider a potential energy function U (x|c) for conformations x that depends on a vector of m force field parameters c i , c = (c 1 , . . . , c m ). Structures are distributed according to the normalized Boltzmann distribution (1) β = 1/k B T is the inverse temperature and k B is Boltzmann's constant. For uncorrelated Gaussian errors, the BioEn posterior is given by the functional [18] where the Kullback-Leibler divergence [33] is given by and where is the chi-squared defined as the sum of the squared errors in units of the estimated standard errors. p(x|c 0 ) is the reference ensemble or c 0 -ensemble. Y k are the M experimental data points, y k (x) are the corresponding calculated observable values for structure x, and σ k is the combined theoretical and experimental error. θ expresses our confidence in the reference ensemble [13,18]. By maximizing the BioEn posterior given by Eq. (2) with respect to p(x), we obtain the BioEn optimal ensemble determined self-consistently by the M generalized forces Angular brackets indicate the average over all structures in the optimal ensemble, i.e., y k = dx p(x) y k (x) for a normalized p(x). Thus, the mean force to restore the reference distribution is exactly balanced by the mean force to fit the data [18]. The generalized forces correspond to Lagrange multipliers in the corresponding expressions for Eq. (5) in the maximum entropy formulations of ensemble refinement [14,28,31,34,35]. The generalized forces parameterizing the optimal ensemble can be efficiently determined by gradient-based optimization methods [36,37]. The optimal ensemble corresponds to a Boltzmann ensemble ∝ exp[−βU (x|c 0 ) − βΔU (x)], where the energy is given by the sum of the potential energy of the reference ensemble βU (x|c 0 ) and a correction term which is linear in the calculated observables y α k = y k (x α ) [14,17,18].
For well sampled ensembles, the bias energy given by Eq. (7) can be used to update force field parameters. However, in general, it is not guaranteed that we can adapt the force field parameters to fit βΔU (x). This correction may not be adequately representable given the functional form of the force field. In the special case of 3 J couplings, Cesari et al. exploited the functional similarity of the Karplus equation to dihedral potentials to calibrate RNA force fields [31]. More recently, Cesari et al. generalized this approach and replaced the observables y k (x) in the expression for the optimal ensemble, Eq. (5), by a potential energy correction function and treated the generalized forces as fitting parameters [38]. A regularized error function keeps the calculated observables close to the measured values and the fitting parameters small. Cesari et al. tune the latter using a parameter to balance experiment and MD simulations similar to the BioEn confidence parameter.
In any systematic approach to force field refinement [5,27,31,38,39], we have to regularize the problem and balance the information from experiment and MD simulations for the given uncertainties. BioEn has both features already built in. In the following, we adapt BioEn for systematic force field refinement.

The BioFF posterior
For force field refinement, we restrict the function space of p(x) in the BioEn posterior, Eq. (2), to probability distributions parameterized by c, i.e., to p(x|c) ∝ exp[−βU (x|c)]. Note that c denotes the parameters we want to optimize and that p(x|c) will usually depend on additional force field parameters. Doing so, the BioEn posterior becomes a function of c. Introducing an additional prior p 0 (c) on the force field parameters, we obtain the BioFF negative log-posterior as The BioFF optimal solution satisfies where prime ( ) denotes a gradient with respect to c. We show next how to optimize the BioFF posterior for a given ensemble of structures, which is at the core of the BioFF-by-reweighting procedure presented below.

BioFF optimization by reweighting
Let us assume we performed a reference simulation using the force field parameters c 0 . For unbiased MD simulations, each of the N structures in this reference ensemble has a statistical weight w  and thus satisfy Note that w α (c) and w The gradient of these weights, which is needed to evaluate the gradient of the negative log-posterior, Eq. (11), is given by The calculation of the Hessian is straightforward. The Hessian is needed for optimization using the (L-)BFGS algorithm, for example. The evaluation of Eq. (13) for the weights w α (c) is computationally cheap such that numerical optimization can be performed efficiently.

Data for multiple different systems
We note that Eqs. (10) and (11) can readily be generalized to incorporate M s data from multiple different systems s = 1, . . . , S sampled with N s conformations each, as will be the case in most force field optimizations. The negative log-likelihood then becomes For each system s, the weights satisfy Eqs. (13) and (14) with w α , x α , and N replaced by w αs , x αs , and N s , respectively. If systems differ also in the thermodynamic state, e.g., temperature, then Eq. (13) and following have to be adapted accordingly.

Data for individual conformations from quantum mechanical calculations
In many practical cases, one also has data Y j (x j ) (j = M + 1, . . . , M + M * ) for individual conformations x j , such as forces or potential energies from quantum mechanical calculations. Note that x j can refer to the same conformation for different values of j. We can include such data in BioFF by adding a term χ 2 * /2 to the negative log-posterior, where y j (x j |c) is the corresponding calculation for the force field with parameters c, and σ j 2 is the combined squared uncertainty of data and calculation. χ 2 * /2 can be added to L(c) in Eq. (10). Its gradient with respect to c is This expression can be added to the gradient L (c) in Eq. (11) to account for data reporting on the properties of individual conformations.

The BioFF-by-reweighting method
In the following, we sketch an iteration procedure to determine the BioFF optimal parameters by reweighting ( Fig. 1).
At iteration 0, we start with a single simulation using force field parameters c i=0 and calculate all M observables y k (x α ) for all N 0 structures in this initial estimate for the reference ensemble. The reference weights are given by w (0) α = 1/N 0 . We minimize Eq. (10) numerically to obtain new parameters c i=1 .
1. We run an additional simulation for the newly optimized force field parameters c i and generate N i structures in the c i -ensemble.
α from all MD simulations at parameters c 0 , . . . , c i using binless WHAM [21,22,24]. We calculate all observables y k (x α ) for the new simulation ensemble. 3. We minimize Eq. (10) numerically to obtain new parameters c i+1 . 4. We check for convergence of the force field parameters. If not converged, we continue with step 1 and the new force field parameters c i+1 . If converged, Fig. 1 Flowchart of the BioFF-by-reweighting optimization procedure. The input (green) are experimental or theoretical data and a structural ensemble generated in a molecular simulation using the reference force field. The calculation of the reference weights and observables, optimization of the BioFF posterior, and running of an additional simulation for the newly optimized parameters are iterated until convergence (blue). The result are optimal force field parameter values (red) c i ≈ c i+1 are the BioFF optimal parameters and we are finished.
In the presented procedure, we start with an unbiased simulation using the reference force field. One can, however, also use any unbiased or biased simulations and different force fields than the reference force field and apply binless WHAM to calculate the reference weights w (0) α in the c 0 -ensemble.

BioFF implementation
The ultimate goal of full automatization of force field optimization requires careful consideration of steps that would otherwise be handled by experienced practitioners. In particular, when reweighting using Eq. (13) we have to ensure that the reference ensemble covers the reweighted ensemble sufficiently. If the population of structures in the overlap region between reference ensemble and reweighted ensembles is too low, then the reweighed ensemble will suffer from artifacts. As a measure of the coverage between the reference and reweighted ensemble we use Kish's effective sample size [40,41] given by In contrast to Rangan et al. [41], we use the effective sample size N itself and not the relative size N /N as a measure. Our iterative optimization procedure ensures good coverage of the optimal ensemble even far from the reference ensemble, which grows with each iteration. Consequently, an absolute measure for the quality of sampling is needed because, relatively, the coverage may get smaller with each iteration.
To be able to use common optimization libraries, we include a threshold on the effective sample size N given in Eq. (18) in the objective function. To account for the reweighting limit in the objective function used for optimization, we introduce with a suitable chosen N * > 0 and k > 0. This function approaches one for N N * . For N approaching N * , it diverges exponentially. Defining the objective function for optimization as we thus ensure that we do not extrapolate further than supported by the data. Iterating optimization and simulation at the new optimal parameter sets, we make sure that O(c) and L(c) converge to the same optimum. As a check we could rerun the final optimization for O(c) = L(c).

Polymer model
We use a simple polymer model to illustrate the principles of the BioFF method. The two-dimensional polymer consists of a string of n beads, with positions r i and r i+1 of neighboring beads separated by a unit of length. That is, |v i | = 1, where the n−1 bond vectors are given Δφ j and ' ' indicates the transpose of the column vector v i .
We introduce a potential energy acting on Δφ i such that the Boltzmann distribution corresponds to the von Mises distribution given by where I 0 (κ i ) is the modified Bessel function of order zero. The parameters κ i > 0 determine the stiffness of the polymer at positions i. μ i introduce directional biases.
The probability of a polymer with n beads is given by where we introduced Δφ = (Δφ 1 , . . . , Δφ n−1 ), κ = (κ 1 , . . . , κ n−1 ), and μ = (μ 1 , . . . , μ n−1 ). Without distance-dependent bead interactions, we can efficiently generate uncorrelated polymer conformations by drawing n − 1 angles from the von Mises distributions given by Eq. (22). We can include distance-dependent bead interactions and sample conformations using Markov chain Monte Carlo simulations and trial moves in the angles Δφ i . Despite its simplicity, the 2(n − 1) parameters of the polymer model facilitate rich structural diversity. Stiff, ordered parts have large κ i -values. Disordered regions, linkers, and short, flexible hinges, have low κ i -values. Turns can be introduced as short stretches with large κ i -values and non-zero μ i -values of the same sign.
The computational efficiency of generating polymer conformations and the simplicity of their visualization in two dimensions make the polymer model an attractive system for method development. We use this polymer model here to sketch the BioFF optimization procedure.

Calculation details
As a simple example, we use the polymer model introduced above with n = 100 beads and κ i = κ and μ i = 0 for i = 1, . . . , n−1. We use the distance between bead 1 and n/2 = 50 as the experimental observable. For validation, we additionally monitor the end-to-end distance given by the distance between bead 1 and n.
We use the same model to generate synthetic experimental data for a chosen value κ exp = 10 of the single force field parameter by drawing N = 10000 random conformations according to Eq. (23). We calculate the experimental value of the label distance as an average over this ensemble.
For the initial reference simulation, we set the value of our force field parameter, corresponding to c 0 , to κ 0 = 20. We draw N 0 = 4000 conformations, which form the initial reference ensemble and for which we calculate the label distances. In contrast to the sampling in MD simulations, the structures in the ensemble are strictly uncorrelated. Correlations would have to be taken into account in the WHAM evaluation.
To detect the limit of the reweighting, we set N * = 300 and k = 10 in Eq. (19). As we will see below, this choice is rather conservative in the example con-sidered here. Lower values of N * or higher values of k will decrease the number of iterations needed for convergence.
To detect convergence, we demand that for at least two consecutive iteration steps the relative difference between the previous and the actual value of κ is smaller than 0.05, i.e., |κ i − κ i−1 |/κ i < 0.05.
For selected values of the confidence parameter θ = 10, 1, and 0.1, we optimize Eq. (20) using the Nelder-Mead method [42] for the reference simulation and w (0) α = 1/N , which returns a new value κ 1 . For this parameter value, we run a simulation, i.e., we draw N 1 = 4000 new structures, and calculate the label distances for all conformations. We then perform binless WHAM [21,22,24] to obtain new values of the reference weights w (0) α for the combined N = N 0 + N 1 structures. With this new set of simulations, we again optimize the objective function given by Eq. (20). We iterate this procedure until we detect convergence. Note that we do not use a prior p 0 (κ) in this example.

Results
We obtain convergence for all values of θ (Fig. 2). For the smallest value of the confidence parameter θ = 0.1, the BioFF optimal value (κ ≈ 10.3) comes close to the value underlying the synthetic experimental data (κ exp = 10).
For decreasing values of the confidence parameter θ, BioFF optimal distance distributions converge to the experimental distribution (Figs. 3 and 4). Reflecting the fact that BioEn is less restrained in its optimization than BioFF, we find that the calculated average label distance of the BioEn optimal solution is consistently closer to the experimental value than the BioFF solution for the same θ value (Fig. 3). However, by correcting for the underlying force field error, BioFF also improves the statistics of the end-to-end distance, which we do not use in the refinement (Fig. 4). By contrast, BioEn gives only small improvements for the end-to-end distance even at large values of θ. Indeed, the BioEnoptimal average end-to-end distance is approximately the same for all θ-values examined here. By contrast, the BioFF end-to-end distance distribution and their average converge towards the reference with decreasing θ even though they are not used in the refinement.
Note that the difference between the BioEn optimal weights and the BioFF optimal weights increases with decreasing values of θ. This difference is quantified by the Kullback-Leibler divergence S KL between optimal weights and the reference weights (see figure legends in Fig. 3). For θ = 0.1, S KL ≈ 0.4 for the BioEn optimal weights and S KL ≈ 14.7 for the BioFF optimal weights. The reason for this large difference between the S KL values is that the label distance distributions of the reference ensemble and the experimental ensemble have a large overlap and thus little BioEn reweighting is needed to obtain agreement. However, the potential energy distributions overlap only very little (Fig. 2), such that the BioFF optimal weights are quite different from the reference weights. By combining the ensembles obtained for different parameters c i with binless WHAM, issues with poor overlap are effectively avoided in the BioFF optimization.

Discussion
BioFF combines a simplified prior on the force field parameters and an entropic prior on the statistical weights of the sampled conformations. For good param-eter priors, the entropic prior can be neglected. A good prior properly accounts for the interdependence of the parameters and their uncertainties. For poor parameter priors, the entropic prior takes care that we do not move too far from the reference ensemble. Indeed, the entropic prior prefers parameter values that keep the Kullback-Leibler divergence to the reference ensemble small.
The BioFF method can be viewed as an umbrella of methods distinguished by the prior on the force field parameters and the value of the confidence parameter θ. If we set θ = 0 then the prior on the force field parameters is the only a-priori regularization of the parame- ter optimization problem. In this case, the prior has to have an appropriate functional form with accurately known uncertainties to properly account for the complex interdependence of parameters. If we use an uninformative prior on the force field parameters and θ > 0 then the a-priori regularization enters through the functional dependence of the weights on these parameters, Eq. (13), and the Kullback-Leibler divergence. If additionally the potential energy contains terms of the same functional form as the observables, Eq. (7), then the BioEn optimal ensemble and BioFF optimal ensemble are the same [31].
The BioEn optimal ensemble poses a limit for the BioFF optimal ensemble. For the same value of θ, BioEn Fig. 4 Distribution of the end-to-end distance and mean values (vertical lines) for the initial reference simulation (κ0 = 20, green), the synthetic experimental data (κexp = 10, black dashed lines), BioEn (blue), and BioFF(orange). The end-to-end distance was not used for refinement. The confidence in the reference ensemble decreases top to bottom (θ = 10, 1, 0.1) will essentially always provide a better fit to the experimental data used in refinement while staying closer to the reference ensemble than BioFF. In BioEn, the Kullback-Leibler divergence restrains the optimal statistical weights. In BioFF, the functional dependence of the statistical weights on the force field parameters and the force field parameter prior provide additional restraints. As a consequence, the less restrained BioEn posterior in Eq. (2) evaluated for the BioFF optimal weights is lower than or, at best, equal to that for the BioEn optimal weights.
Typically, one will select certain force field parameters for optimization based on prior knowledge. Experience with the force field and physical insight on which structural features are probed by the experimental observables guide this process. For example, if an atomistic force field provides an ensemble of intrinsically disordered proteins that is too compact then likely candidates are dihedral angles and the water-protein interaction [43][44][45].
In a more systematic approach, force field parameters can be identified as candidates for optimization by a sensitivity analysis. Given a set of ensembles corresponding to U (x|c), we can use binless WHAM to estimate the sensitivity of the estimated observable averages y k c as a function of c to the different force field parameters c j by Here, Δc j is the uncertainty in coefficient c j as encoded, e.g., in a Gaussian prior. To linear order, the dimensionless |J kj | reports on the improvement achievable for observable k by varying parameter j. Indeed, we have Δc j ∂χ 2 /∂c j = 2 k J kj , i.e., the sum of the J kj is twice the change in χ 2 associated with a small change Δc j . In general, we select parameters that sensitively impact χ 2 . It remains to be seen to what extent the complexity of the force fields, the large number of their parameters, their poorly quantified couplings, and the importance of solvent-mediated interactions and the required Boltzmann averaging will require fine tuning of the optimization approach and human intervention. BioFF optimization enables us to learn about the interdependence of parameters. For the BioFF optimal solutions, we can perform a sensitivity analysis by evaluating the Hessian of the log-posterior at the optimum. Performing this analysis for different sets of parameters, we are able to quantify how informative the experimental data are with respect to individual parameters. Importantly, the off-diagonal elements provide information about the interdependence of the parameters. This information is conditioned on the molecular systems we are simulating and the observables that have been measured in experiments and used for parameter refinement.
The reweighting approach in BioFF is more efficient for soft degrees of freedom but it can also handle stiffer degrees of freedom. An example for a soft degree is the dispersion energy parameter of the Lennard-Jones interaction. Adapting this parameter, mainly the depth of the minimum changes. In general, conformations will be part of both the reference and the reweighted ensemble with significant albeit different weights. By contrast, the diameter of a Lennard-Jones particle is a stiff degree of freedom due to the diverging repulsive part of the interaction. A small decrease of the diameter tends to introduce new conformations not previously sampled in the reference simulation. Consequently, reweighting is less efficient in predicting ensembles for decreased diameters. In BioFF, the iteration of parameter optimization for a given ensemble and running new MD simulations for optimized parameters will converge eventually. Con-vergence can be sped up by softening stiff degrees of freedom, e.g., by replacing the diverging repulsive interaction by a soft-core potential of finite height, and using binless WHAM to obtain the proper reference weights.

Conclusions and outlook
Optimized force fields from diverse data and systems. Force field parameterization is a complex process with broad implications for the MD simulation community. In light of the many possible pitfalls, the optimization has to be performed with great care. Seemingly small corrections of certain force field parameters driven by data for one system may cause unexpected, negative consequences in MD simulations of other systems. In practice, it is therefore important to include a wide range of systems and data in the optimization process and to test the results on systems and data not used in the optimization process.
Covering diverse systems and data requires careful accounting for the respective errors. However, uncertainty quantification is a recognized challenge [9,46,47], as we are dealing with errors in experiments, quantum mechanical calculations, statistical sampling, and calculation of observables. The respective errors and error models will have to be carefully and regularly re-assessed to minimize the risk that data with poorly modeled errors skew the optimization process.
In this context, it is also worth noting that coarsegraining may be considered a force field learning and optimization problem [48,49]. In this vein, Chen et al. used experimental inputs to learn coarse-grained models by constructing a Markov state model and expressing the likelihood as a function of the model parameters [48].
Force field optimization as a community effort. The requirement to cover a wide range of systems and data makes force field improvement a community effort. Data and computing resources of individual researchers are necessarily limited. To include as wide a range of inputs and tests as possible, data should be pooled, and optimization and testing tasks should be distributed and ultimately collected. Indeed, ultimate validation happens in the field by scientists applying these force field and checking for consistency.
Practitioners choose force fields depending on the molecular system at hand, their research question, their resources, and their prior experience with molecular force fields. They can compensate for known deficiencies by applying biases during MD simulations or by adapting their analysis, for example. However, for structurally and chemically heterogeneous systems, large system sizes, and long simulation times, it becomes more difficult to find a force field that satisfies all demands. Therefore, the development of computationally efficient force fields that accurately account for a wide range of data for diverse systems becomes critical.
Systematic and automated force field optimization. We can only meet the steadily growing demands on the accuracy of molecular force fields by applying systematic and comprehensive approaches to force field parameterization. We envision that in a systematic and (semi)automatic process, force field parameters are validated and updated and their uncertainties are quantified. Experimental data and their uncertainties would be collected in an open data base for wide ranges of methods, systems, and experimental conditions. Importantly, the results of high-level quantum chemical calculations can be readily included in the process, both at the level of priors (e.g., by defining tolerable variations of certain parameters) and of data (by entering into χ 2 ; see Sect. 2.5). For each of the systems, molecular models are built. Uncertainties in the modeling are accounted for by building multiple models for the same system and weighing them according to prior knowledge in the subsequent refinement. One then runs (un)biased MD simulations for these systems, combines them with previous runs, and quantifies the sampling errors. From these ensembles, one calculates the experimental observables and quantifies the uncertainties in these calculations. By performing Bayesian inference of the parameters, e.g., using BioFF or similar methods, the different sources of errors can be taken into account. Validation of the new parameters involves systems, observables, and data not included in the optimization.
Following this rough outline, resources could be set aside at a super-computing center to regularly perform force field parameterization, validation, and uncertainty quantification. Users can submit experimental data, molecular models, quantum mechanical calculations, and trajectories. With this steadily growing wealth in data, models, and MD simulations one could refine force fields with different emphasis. Specialized force fields can be parameterized for protein folding, liquid-liquid phase separations, protein-ligand binding, and so on. For structurally and chemically heterogeneous systems, as we encounter in (sub)cellular structures, we need to parameterize general-purpose force fields. The comparison of general-purpose force fields with specialized force fields also enables us to identify the limitations stemming from the chosen functional form of the force fields, i.e., model adequacy [9].
BioFF produces experimentally refined simulation models. The BioFF approach also addresses a fundamental problem in MD simulations, namely whether to accept or reject the simulation model as a faithful description of the system of interest. Quantitative comparisons of MD simulations to experiment often end in a frustrating experience, in which for a given system some of the observables agree well with experiment, yet others not so well or not at all. One is then left with the decision whether the simulation model should be discarded as a whole (on the basis of partial but significant disagreements) or accepted for further analysis. This decision is a serious one since one of the goals of molecular simulation studies is to gain insight beyond what can be deduced readily from experiment alone. By providing a means to incorporate the experimental information into the model, frameworks like BioFF allow one to refine the original simulation model "on the fly" to account for the measurements. By providing an experimentally refined model, BioFF, like BioEn, offers a possible route out of the dilemma of accepting or rejecting the original model.
The outcome of BioFF is not only an empirically optimized force field but also an ensemble of molecular conformations that can then be analyzed further. As such, BioFF can be viewed as an enhanced sampling method. Seen from this perspective, the procedure addresses the dual challenges of obtaining better simulation models and of gaining insight into molecular systems on the basis of MD simulations. These simulations capture "all" available experimental information and remain solidly based on physics, which enters, e.g., through quantum mechanical calculations used in the parameterization.
Improving the functional form of force fields. BioFF, like many other approaches [5,27,31,38,39], tries to work around one main challenge in force field refinement: the lack of good force field priors. Ideally, priors not only act on the parameters but also on the functional form of the force fields. Going forward with any of these methods, we should be able to better understand force fields, the interdependence of their parameters and their uncertainties, and the limitations due to their functional forms. This effort should also lead to better priors and methods to develop force fields along the way. Python 3 source code using Numpy [50], Scipy [51], Matplotlib [52], and Numba [53] and Jupyter notebooks [54] to perform BioFF for the simple example case presented here can be found at https://github.com/ bio-phys/BioFF. This code uses the open-source BioEn optimization software (https://github.com/bio-phys/ BioEn) and an open-source binless WHAM implementation (https://github.com/bio-phys/binless-wham).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.