{'fraction_non_alphanumeric': 0.02121085326209405, 'fraction_numerical': 0.004763590717435386, 'mean_word_length': 5.114938004329856, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Towards empirical force fields that match experimental observables

Biomolecular force fields have been traditionally derived based on a mixture of reference quantum chemistry data and experimental information obtained on small fragments. However, the possibility to run extensive molecular dynamics simulations on larger systems achieving ergodic sampling is paving the way to directly using such simulations along with solution experiments obtained on macromolecular systems. Recently, a number of methods have been introduced to automatize this approach. Here we review these methods, highlight their relationship with machine learning methods, and discuss the open challenges in the field.


I. INTRODUCTION
Classical molecular dynamics (MD) simulations at the atomistic scale offer a unique opportunity to model the conformational dynamics of biomolecular systems. Being able to reveal mechanisms at spatial and temporal scales that are difficult to observe experimentally, MD simulations are often seen as a computational microscope. 1 In the past years, they have been applied to study problems ranging from protein folding 2 and aggregation 3 to RNA-protein interactions, 4,5 transmembrane proteins dynamics, 6 and full viruses, 7 bacteria, 8 or organelles. 9 The capability of MD simulations to reproduce and predict experimental results is limited by the statistical errors arising from the finite length of simulations and by the systematic errors resulting from the inaccuracies of the underlying models. Interactions are often modeled using empirically parametrized force fields that allow timescales of the order of the microsecond to be routinely simulated. Importantly, the two sources of error mentioned above are deeply intertwined, because only systematic errors that are larger than statistical errors can be detected by comparison with reference experimental results. Indeed, in the past 20 years, the use of special purpose hardware, 10 optimized software, 11,12 and enhanced sampling methods, 13,14 has significantly reduced the statistical errors, thereby allowing force fields inaccuracies to be detected and largely alleviated. In spite of this, empirical force fields are still far from perfect and in some cases are poorly predictive. For instance, it is not trivial to have force fields capable of simultaneously describing correctly folded, disordered, single-chain proteins or protein complexes, 15,16 to correctly predict RNA structure from sequence-only information across a wide range of structural motifs, 17 or to reproduce experimental kinetics in ligand-receptor systems. 18 Solution experiments are optimally suited for validation of force fields, since they provide information about transiently populated structures as well, and they have traditionally been used in this sense. Nevertheless, several approaches have enabled solution experiments to be used directly during forcefield fitting, on a par with available quantum chemistry data. The aim of this perspective is to review these approaches, highlight their relationship with machine learning methods, and discuss the open challenges in the field. We will use here as paradigmatic examples some of the force fields that are most used for simulating biomolecular systems, namely AMBER, 19 CHARMM, 20 OPLS, 21 and GROMOS. 22 All the mentioned force fields share a common functional form, including bond stretching, angle potentials, torsional potentials, Lennard-Jones, and electrostatic interactions: The parameters (k b ; r 0 ; k a ; a 0 ;V n ; σ ; ε; q) are derived from small fragments in advance and depend on the atom type and its chemical environment. Polarizable force fields (such as AMOEBA 23 and a variant of CHARMM 24 ), reactive force fields (such as ReaxFF 25 ), and semi-empirical methods (such as DFTB 26 ) have different functional forms but similar considerations can be applied. The parameters in Eq. 1 are derived with a variety of different procedures that depend on the specific force field and are summarized in Table I. In particular, some of the parameters are derived from quantum chemistry calculations performed at a varying level of accuracy, in a bottom-up spirit. Other parameters are instead derived from experimental data, either using spectroscopy experiments, databases of crystallographic structures, or other gas-phase or solution-phase experiments, in a top-down spirit. The reliability of a force field largely depends on the accuracy of the employed reference data. For instance, a force field fitted purely on quantum chemistry data cannot provide results that are more accurate than the reference method. However, this limit can be surpassed if multiple sources of data are combined. As an additional and perhaps even more important source of error, one should take into account that reference data used in force-field fitting, either computational or experimental ones, are obtained studying systems that are necessarily not identical to those that one wants to simulate later. For instance, torsional parameters and partial charges in Traditional force-field parametetrization and risks associated to extrapolations. The left panel shows the traditional procedure used for force-field parametetrization. Parameters are obtained from calculations or experiments on small molecules or fragments. Simulations are then validated for their capability to maintain the native structure of a macromolecule or against solution experiments. Since fitting and validation are done on different types of systems, there is a large risk associated to extrapolation. The right panels show typical errors observed when fitting a function. The horizontal axis represents a configurational coordinate (e.g., a dihedral angle) and the vertical axis an observable that is used for fitting (e.g., the energy of the system). The true function is shown as a solid line, and the available reference data are shown as grey points. Lines fitted on the reference data using a simple (linear) and a complex (quadratic) form are shown as dotted and dashed, respectively. Data are collected on a narrow (upper panel) or wide (lower panel) range of configurations. The simple model represents a force field with too few parameters or with an incorrect functional form. When fitted on a narrow range of configurations (upper panel) it reproduces well the true function. However, it fails the extrapolation to the right part of the graph. When fitted on a wide range of configurations (lower panel) the intrinsic limited transferability of the model emerges from the error observed on the fitted points. The complex model represents a force field with more parameters and a more physical functional form. When fitted on a narrow range of configurations (upper panel) it can lead to significant overfitting. Conversely, when fitted on a wide range of configurations (lower panel) it reproduces well the true function on the entire range of configurations.
the AMBER force field are traditionally obtained using quantum chemistry calculations in small fragments of up to a few dozen atoms, typically including a couple of aminoacids, but are later used to simulate oligopeptides or full protein domains (see Fig. 1, left panel). Similarly, Lennard-Jones parameters in the OPLS force field are obtained from vaporization calorimetry of pure organic liquids such as tetrahydro-furan, pyridine or benzene, but then applied to cases where the analyzed compounds are only portions of a sugar, nucleobase, or aminoacid respectively. The reliability of a force field when used in a context different from the one in which it was parametrized depends on the transferability of the functional form in Eq. 1 (see Fig. 1, right panels). Given the very large gap between the size and complexity of the systems used for parameter fitting and the systems to which force fields are applied, it appears almost a miracle that current force fields are, for instance, capable of correctly identifying the folded state of a protein. 2 It is interesting to look at a few anecdotal examples to better understand how this is possible. The traditional AM-BER force field for nucleic acids has been used for several years before it was realized that sufficiently long simulations could lead to a transition to experimentally unobserved rotamers in the α and γ torsions of DNA backbone. 27,28 Following this empirical observation, a joint effort of several groups lead to the parmbsc0 reparameterization of DNA backbone, 28 where the parameters corresponding to these two torsional angles were fitted against quantum chemistry calculations. A similar episode occurred later with the χ OL3 corrections, derived to counteract the occurrence of ladder-like structures in RNA. 29 For protein systems, one of the most important additions after the initial development of the CHARMM force field has been the introduction of empirical corrections maps (CMAP) 30 , that deviate from the functional form of Eq. 1 by the presence of coupling terms between consecutive torsional angles. These corrections were fitted on quantum chemistry data, but required also a heuristic adjustment to fix the typical values of torsional angles in α-helical and β -sheet regions. As a further example, empirical adjustments of the AMBER and CHARMM force fields were performed respectively in Refs. 31 and 32, where solution data on short oligopeptides were used to optimize backbone dihedrals so as to reproduce helix-coil transitions.
A general trend that can be seen is that experimental data on macromolecular systems (e.g., nucleic acids duplexes or protein domains) are typically used for validation, whereas the parameters are fitted on either theoretical or experimental information available for much smaller systems. Nonetheless, the observation of failures in macromolecular systems is the only way to detect which precise parameters should be corrected. The last two mentioned works, 31,32 instead, report direct fitting of parameters on simulations of short oligomers.

III. RECENT APPROACHES FOR FITTING FORCE FIELD PARAMETERS ON EXPERIMENTAL DATA
A number of approaches have been introduced to allow fitting force fields directly on experimental data taken on macromolecular systems rather than on small fragments, all of them following a flowchart similar to the one illustrated in Fig. 2.
Since solution experiments often report results that are averaged over an ensemble of copies of the same molecule, these methods are typically designed to enforce ensemble averages rather than instantaneous values. Norgaard et al. 33 introduced  2. Schematic representation of a force-field fitting procedure. Initial parameters are tuned based on both quantum chemistry data and experimental data on small systems (e.g., individual residues). Molecular dynamics simulations are then performed on macromolecular systems. Reweighting is used to optimize force field parameters in order to maximize the agreement with a set of available data including experiments on macromolecular systems. In principle, this second stage might include also quantum chemistry data and experimental data on small systems. Even when not explicitly used at this stage, the initial set of quantum chemistry and experimental data is still playing a role for all the parameters that are not further adjusted. Even on the adjusted parameters, the information about the initial force field remains present if regularization terms are included.
an approach where a force field is iteratively refined until agreement with experiment is obtained. At each iteration, a simulation is performed and the force field parameters are optimized by assigning new weights to the visited conformations. Thus, through such reweighting procedure, one can predict what result would be obtained using these slightly modified parameters. At some point, when the refined force field and the initial one become too different, it is necessary to iterate the procedure performing a new simulation. The method was applied to the refinement of a coarse-grained model of a protein and fitted against paramagnetic relaxation enhancement experiments. Li et al. 34 showed how to refine an allatom protein force field using chemical shifts and full-length protein simulations. A common trait of all these methods is that even small changes in force field parameters can make the resulting ensemble very different from the original one making the reweighting procedure less accurate. In Ref. 34, a local reweighting procedure was introduced to alleviate this issue. This procedure is based on the heuristic observation that the ensemble of conformations accessible to a residue is maximally affected by the parameters used for that residue and, to a lesser extent, by the parameters used for the other (possibly identical) residues. Since this is an approximation, a subsequent simulation performed with the corrected force field was necessary to validate the modification. Refs. 35 and 36 used a similar automatic procedure to optimize water models. Interestingly, they realized that a straightforward fitting procedure might lead to overfitting and showed how a regularization term can be included in order to alleviate this issue. Chen et al. 37 used a force-field fitting procedure to develop a coarse-grained model for proteins based on synthetic experimental data generated using a long trajectory obtained with an atomistic model. 38 Finally, Cesari et al. 39 introduced a procedure to refine atomistic force fields where heterogenous systems and types of experimental data are used to refine the AMBER RNA force field. Enhanced sampling techniques are employed to ergodically sample the conformational space for a number of RNA tetramers and hairpin loops and a regularization term is used in the fitting scheme to maintain the refined force field close to the initial one. The weight of the regularization term is chosen with a cross-validation procedure aimed at maximizing the transferability of the parameters.
It is important to recognize the difference between the mentioned approaches, that are meant to generate transferable force-field parameters, and methods meant to improve the agreement with experiment for a specific system for which data are available. 40 This second class includes a variety of approaches such as Bayesian schemes 41,42 and methods based on the maximum entropy principle. 43,44 In the maximum entropy formalism the number of free parameters is equal to the number of experimental datapoints. For instance, in a homogenous polymer, each of the monomers will feel a different correction that makes its structure as compatible as possible with experiments. Since the number of parameters is very high, regularization methods can be used and tuned with a cross-validation procedure (see, e.g., Ref. 45 dition, if a polymer of a different length needs to be simulated, new experimental data should be obtained. In forcefield fitting procedures, the chemical structure of the investigated molecule is a priori used to reduce the number of parameters. For instance, in a homogenous polymer, each of the monomers will feel the same correction (although perhaps terminal monomers might be treated differently 46 ). On the one hand, this allows to encode a large amount of information in the specific choice of the functional form employed. This type of information is similar to the one that is included when atoms are classified in types in order to obtain their parameters. 47 On the other hand, it significantly reduces the number of parameters potentially making the resulting force field transferable. Ref. 48 used an hybrid approach were maximum entropy restraints were used but kept by construction constant across chemically equivalent parts of the system. For a recent comparison of approaches taken from both classes, see Ref. 49.
Besides the discussed systematic approaches, that report methodological improvements aimed at optimizing parameters based on experimental data, a number of recently developed force fields include terms that were chosen based on the result of MD simulations on systems of different complexity and their capability to reproduce experimental data. so as to provide correct salt-brigde interactions. 50 Finally, two recent variants of the AMBER RNA force field contain corrections on hydrogen bonds obtained scanning a series of parameters and minimizing the discrepancy with solution experiment for RNA oligomers. 46,51

IV. THE MACHINE LEARNING LESSON: HOW TO AVOID OVERFITTING
Overfitting is a ubiquitous problem when fitting procedures are done in a blind manner. The prototypical cases are machine learning algorithms where functions of arbitrary complexity, supported by no or little physical understanding, are used to fit empirical data. The machine learning community has thus developed a number of tools that can be used to avoid or at least alleviate this issue.
Many different machine learning techniques exist and are typically based on a common framework. 52 The basic ingredient is a dataset made up of a matrix X of independent variables (data) and a set Y of dependent variables (labels). Next, a set of models is proposed to map X into Y with best accuracy. A model is defined by a set of parameters plus a set of hyperparameters. This splitting is guided by computational convenience such that inference can be approached in a multi-level fashion: typically, model parameters are found by solving an optimization problem at fixed hyperparameters, that on the other hand are preferably scanned over a discrete scale. This double approach is more easily understood when another basic ingredient of machine learning is introduced, that is the cost function. The cost function is used to estimate the performance of a model, and while it is usually a continuous function of the model parameters, it can have a non-trivial dependence on the hyperparameters. For example, the set of hyperparameters can include the architecture of the model, the optimization algorithm used to find the optimal model parameters, the functional form of the cost function itself, etc. Since the sets of parameters and hyperparameters defining models are fitted against a finite set of examples {X, Y}, overfitting can easily occur. In the limit of fitting on an infinite amount of data, the only limitation of a model would be determined by its complexity. In this limit a too much simple model would underfit the data, leading to a bias in the result. This bias can be decreased by increasing the model complexity. But since in general we deal with datasets of finite size, increasing the complexity of the model would result in a large contribution to the error (variance) due to the sampling. A too much complex model would overfit the data, thus having a seriously low performance on new independent data. The search for the model with the optimal tradeoff between bias and variance (i.e. between under-and over-fitting) follows two directions. One is to split the dataset into a training and a cross-validation set, prior to analysis. Model parameters are fitted against data in the training set, and afterwards the optimized model is validated against the validation set data not included in the training procedure. This procedure is usually referred to as crossvalidation (Fig. 4). The other is to reduce the risk of overfitting by means of regularization techniques, the most common consisting in adding terms to the cost function, that prevent the model parameters from reaching values extremely adapted to the dataset. This comes at the cost of increasing the number of hyperparameters (e.g., the relative size of training and cross-validation sets, their composition, coefficients of regularization terms, etc.) that continue to be affected by risk of overfitting. Even if a close solution to this problem is not established yet, overfitting should be taken into account for each level of inference (for both parameters and hyperparameters). The most straightforward way to deal with this multi-level risk of overfitting is to a priori split the dataset into three subsets: in addition to the standard training and cross-validation subsets, an independent test set is introduced. The training set is used to fit the optimal values of parameters at fixed hyperparameters; optimal hyperparameters are then fitted against the cross-validation set. Eventually, the performance of the model defined by the optimal parameters and hyperparameters is evaluated on the test set. A more robust approach consists in nested cross-validation 53 , in which parameters and hyperparameters are optimized on a single dataset, but the criterion used to optimize model parameters (training) is different from the optimization criterion used for hyperparameters (model selection). Also in this case, validation of the selected opti-FIG. 5. Free-energy landscapes when using reweighting and iterative simulations. The native state, compatible with experimental information, is in the middle (N). Two metastable non-native states that, based on experimental information, are supposed to have a low population, are also shown (U1 and U2). In the original force field (panel a), the N state is sampled, but the most stable state is U1. During the reweighting procedure, the force field learns how to improve the agreement with experiments by disfavoring U1. However, since U2 was never observed in this simulation, there is nothing that prevents it to be stable when using the refined force field. Once a simulation is performed with the refined force field (panel b), the state U2 appears with large population leading to disagreement with experiment. In principle, if a reweighting is performed using only the second simulation, state U1 might appear again with an incorrect population. Only a reweighting where both simulations are combined, and thus all the possible states can be observed, is capable of generating a force field that correctly sets N as the global free-energy minimum and U1 and U2 as metastable states with low population (panel c). mized model against new data that, importantly, has not been used to adjust neither parameters nor hyperparameters, is best practice.

V. OVERFITTING IN FORCE FIELD DEVELOPMENT
Force-field fitting procedures can be interpreted as machine learning methods where the parameters are the optimized coefficients and data and labels are a mixture of information obtained from both quantum chemistry calculations and various experimental techniques. One should thus pay attention to overfitting. Whenever overfitting occurs, transferability of the Force-field fitting can be based on a combination of quantum chemistry data, experimental data on small systems and experimental data on macromolecular systems. All data can be used in parameter fitting, and leave-one-out cross-validation (see Fig. 4) can be used to help the choice of hyperparameters. To this end, it is necessary that a separate set of data, either theoretical or experimental, is left out until the very end of the procedure to validate the transferability of the model. This separate data should not be used to take any decision during the fitting procedure, or the information leak might make the final validation not truly independent.
force field to a different case might be compromised. As already discussed, if parameters are only fitted on small systems, their transferability to larger systems might be limited. The other phenomenon that can be observed in reweighting methods is the subtle overfitting on the analyzed trajectory.
In particular, if parameters are derived to match experimental data by reweighting a trajectory that is not sufficiently long, they might not work correctly on another trajectory obtained using the same force field but with different initial conditions. In addition, since reweighting schemes can only modulate the weight of states that have been explored but cannot predict the population of states that have not been observed (see, e.g., Ref. 54 for a comparison of restraining and reweighting when used to implement the maximum entropy principle), the only way to detect these problems is to keep the target force field as close as possible to the original one with some form of regularization and to then perform a new simulation once parameters have been optimized (see Fig. 5). Every other decision taken in the path should be included in the list of hyperparameters. Coefficients controlling regularization terms used in the optimization, that control the relative weight of initial force field and of experimental data, are naturally considered hyperparameters. The functional form of the force field itself is a hyperparameter. The analogous of these hyperparameters in the training of neural network are regularization terms or early stopping criteria and network architecture, respectively. 55 In force-field fitting, a number of ad-ditional hyperparameters might be used whose control might be more or less explicit. For instance, the so-called forward models used to calculate experimental observables from MD trajectories contain a number of parameters. If the training is done to reproduce the energetics of quantum chemistry calculations, the set of structures used for fitting and their relative weights are to be considered as hyperparameters. Even the precise quantum methods used to compute the total energy might contain a number of hidden parameters (e.g., the possibility to use either implicit or explicit solvent or the method used to solve the many-body Schrödinger equation).
If hyperparameters are chosen a priori based on some independent intuition or information, for instance the fact that a given quantum chemistry method is more accurate than another one, then this extra information will be encoded in the final result, improving the quality of the resulting model. However, if hyperparameters are optimized by monitoring the performance of the force field on a specific system, then this system will implicitly become part of the training set. Thus, the resulting model should be validated against a separate system (Fig. 6). A practical example would be if different variants of a force field are derived using three different quantum chemistry methods, then the best method is chosen evaluating the performance of the stability of the native structure of a specific system. Unless there are other independent evidences that the selected quantum chemistry method is better than the other ones, this choice should be considered as fitted on the specific system and should be then validated on an independent one. Therefore, as a final remark, all the decisions taken in the process should be critically evaluated in this respect.

VI. CRITICAL ISSUES AND OPEN CHALLENGES
The recent works done in adjusting force field parameters including experimental data suggests that this is a promising field that will lead to important improvements in the future. There are however a number of critical issues that one should carefully consider.
First, we suggest that all input data should be considered at the same level, irrespectively of being obtained from experiment or from quantum chemistry calculations. All types of input data can indeed be equivalently used for training or for validation, taking into account their relative errors and the different information content. Particularly valuable are data obtained on systems as close as possible to those that one is interested in simulating. Less weight instead should be given to data obtained in very different conditions (e.g., without solvent) or on systems that are too simple to be considered as representative (e.g., individual aminoacids or nucleotides). As an exception to this general rule, one should consider that different types of data typically give access to the energetics in different portions of the conformational space. For instance, solution experiments on macromolecular systems are valuable in providing the relative stability of structures that can be distinguished using some probe. Quantum chemistry calculations are instead valuable when states are difficult to be distinguished in the experiment, or when probing rarely visited states (such as transition states).
Reference data should be obtained in conditions as realistic as possible. For what concerns experimental data, one should carefully consider the specific conditions in which experiments are carried out, and prefer experiments performed in conditions that can be reproduced in MD simulations. Ideally, specific experiments might be designed and performed in order to facilitate force field development. When instead basing the fits on quantum chemistry calculations, one should consider the importance of the solvent. Additionally, errors in the experimental data should be taken into account, as well as errors in the forward models used to connect structures with experiments. This is also true for errors in the quantum chemistry calculations.
Taking inspiration from the machine learning community, it is fundamental to learn how to avoid overfitting. In particular, overfitting on specific systems should be avoided and this can be achieved by including as heterogenous as possible systems in the dataset. Similarly, overfitting should be avoided on specific trajectories. To this end, separate validation simulations can be run or robust estimates of the statistical errors can be pursued. Regularization terms can be used to tune model complexity thus reducing the impact of overfitting. Validation should be made on data that are obtained in an as independent as possible manner.
Finally, the current functional form (Eq. 1) might be too limited to be usable on a wide range of cases. Increasing the complexity of the model might help in this respect. Complexity can be introduced by physical insight (e.g., polarizable force fields) or by blind learning of non-linear models (e.g., neural network potentials). Nonetheless, one should keep in mind that, whenever complexity is increased, overfitting has more chance to appear. In this respect, for a fixed number of parameters, the more physical the functional form is, the less it will tend to overfit. Interestingly, neural networks are now routinely used to fit bottom up potentials where the training data can be generated by computational methods and can then be easily made very abundant. [56][57][58][59] These approaches are however typically designed to be trained on very small systems or chemical groups, and their applicability to macromolecular systems has not been showed yet. It is thus still to be seen if neural network potentials can be used fruitfully when force fields are directly fitted on experimental data. ACKNOWLEDGMENTS Sandro Bottaro and Lillian Chong are acknowledged for carefully reading the manuscript and providing useful suggestions. Max Bonomi is greatly acknowledged for carefully reading the manuscript and reporting enlightening feedbacks. Data sharing not applicable -no new data generated