In 1929 Paul Dirac (1, p. 714) stated,
The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble. It therefore becomes desirable that approximate practical methods of applying quantum mechanics should be developed, which can lead to an explanation of the main features of complex atomic systems without too much computation.
Ninety years later, this quote is still state of the art. However, in the past decade, new tools from the rapidly developing field of machine learning (ML) have started to make a significant impact on the development of approximate methods for complex atomic systems, bypassing the direct solution of "equations much too complicated to be soluble."
ML aims at extracting complex patterns and relationships from large data sets, to predict specific properties of the data. A classical application of ML is to the problem of image classification in which descriptive labels need to be associated with images that are presented in terms of sets of pixels. The machine is trained on a large number of examples and learns how to classify new images. The underlying idea is that there exists a complex relationship between the input (the pixels) and the output (the labels) that is unknown in its explicit form but can be inferred by a suitable algorithm. Clearly, such an operating principle can be very useful in the description of atomic and molecular systems as well. We know that atomistic configurations dictate the chemical properties, and the machine can learn to associate the latter with the former without solving first-principle equations, if presented with enough examples.
Although different ML tools are available and have been applied to molecular simulation [e.g., kernel methods (2)], here we mostly focus on the use of neural networks, now often synonymous with the term deep learning. We assume the reader has basic knowledge of ML, and we refer to the literature for an introduction to statistical learning theory (3,4) and deep learning (5,6).
One of the first applications of ML in chemistry has been to extract classical potential energy surfaces from quantum mechanical (QM) calculations, to efficiently perform molecular dynamics (MD) simulations that can incorporate quantum effects. The seminal work of Behler & Parrinello (7) in this direction has opened the way to a now rapidly advancing area of research (8)(9)(10)(11)(12)(13)(14)(15)(16)(17). In addition to atomistic force fields, effective molecular models at resolution coarser than atomistic can be designed by ML (18)(19)(20)(21). Analysis and simulation of MD trajectories have also been affected by ML-for instance, for the definition of optimal reaction coordinates (22)(23)(24)(25)(26)(27), the estimation of free energy surfaces (25,(28)(29)(30), the construction of Markov state models (24,26,31) and dynamic graphical models (32) of molecular kinetics, and the enhancement of MD sampling by learning bias potentials (33)(34)(35)(36)(37) or selecting starting configurations by active learning (38)(39)(40). Finally, ML can be used to generate samples from the equilibrium distribution of a molecular system without performing MD altogether, as proposed in the recently introduced Boltzmann generators (41). A selection of these topics are reviewed and discussed in this article.
All these different aspects of molecular simulation have evolved independently so far. For instance, ML-generated force fields have mostly been developed and applied on small molecules and ordered solids, while the analysis of MD trajectories is mostly relevant for the simulation of large, flexible molecules like proteins. In order to really revolutionize the field, these tools and methods need to evolve, become more scalable and transferable, and converge into a complete pipeline for the simulation and analysis of molecular systems. There are still some significant challenges toward this goal, as we discuss here, but considering the rapid progress of the past few years, we envision that in the near future ML will have transformed the way molecular systems are studied in silico.
This review focuses on physics-based ML approaches for molecular simulation. ML is also having a big impact in other areas of chemistry without involving a physics-based model-for example, by enabling direct attempts to predict physicochemical or pharmaceutical properties (42)(43)(44)(45) or to design materials and molecules with certain desirable properties by using generative learning (46)(47)(48)(49). We refer the interested reader to other recent reviews on the subject (50,51).
The review is organized as follows. We start by describing the most important ML problems and principles for molecular simulation (Section 2). A crucial aspect of the application of ML in molecular simulations is to incorporate physical constraints, and we discuss this for the most commonly used physical symmetries and invariances for molecular systems (Section 3). We then provide examples of specific ML methods and applications for molecular simulation tasks, focusing on deep learning and the neural network architectures involved (Section 4). We conclude by outlining open problems and pointing out possible approaches to their solution (Section 5).
In this section we discuss how several of the open challenges in the simulation of molecular systems can be formulated as ML problems and describe the recent efforts to address them.
MD and Markov chain Monte Carlo simulations employing classical force fields within the Born-Oppenheimer approximation constitute the cornerstone of contemporary atomistic modeling in chemistry, biology, and materials science. These methods perform importance sampling; that is, in the long run, they sample states x from the molecular system's equilibrium distribution, which has the general form μ(x) ∝ e -u(x) .
1.
The reduced potential u(x) contains terms that depend on the thermodynamic constraints (e.g., fixed temperature, pressure). In the canonical ensemble (fixed number of particles, volume, and temperature), u(x) = U (x)/k B T , where k B T is the thermal energy at temperature T. However, the predictive power of these simulations is only as good as the underlying potential energy surface (PES). Hence, predictive simulations of properties and functions of molecular systems require an accurate description of the global PES, U(x), where x indicates the nuclear Cartesian coordinates. All many-body interactions between electrons are encoded in the U (x) function. Although U (x) could be obtained on the fly using explicit ab initio calculations, more efficient approaches that can access long timescales are required to understand relevant phenomena in large molecular systems. A plethora of classical mechanistic approximations to U (x) exist, in which the parameters are typically fitted to a small set of ab initio calculations or experimental data (52)(53)(54)(55). Unfortunately, these classical approximations often suffer from a lack of transferability and can yield accurate results only close to the conditions (geometries) they have been fitted to.
Alternatively, sophisticated ML approaches that can accurately reproduce the global PES for elemental materials (7,9,14,(56)(57)(58) and small molecules (10,11,(58)(59)(60)(61)(62) have recently been developed (see Figure 1). Such methods learn a model of the PES, Û (x, θ ), where the parameters θ are optimized by energy matching and/or by force matching. In energy matching, an ML model is trained to minimize the loss function Constructing force field models with machine learning. (a) Reference geometries are sampled from a molecular dynamics trajectory that is sufficiently long to ensure optimal coverage of the configuration space. (b) For a small subset of geometries, energy and force labels are then computed at a high level of theory to form the training, validation, and test data sets. A globally consistent atom-atom assignment across the whole training set enables the identification and reconstructive exploitation of relevant spatial and temporal physical symmetries of the molecular dynamics. (c) The resulting machine learning model of the potential energy surface is finally used to speed up sampling-intensive path-integral molecular dynamics simulations at the accuracy of the reference electronic structure method (10). (d) An interatomic distance distribution, sampled with classical and path-integral molecular dynamics, is shown. Abbreviations: CCSD(T), coupled-cluster single double (triple); DFT, density functional theory.
where U i are energy values obtained by QM calculations at specific configurations (see Figure 1). For force matching, we compute the QM forces at specified configurations
and minimize the loss function:
The existing ML PES models are based either on nonlinear kernel learning (9,10,57,61) or on (deep) neural networks (14,56,62). Specific neural network architectures are discussed in Sections 4.1 and 4.2. Both approaches have advantages and limitations. The advantage of the kernel methods is that their convex nature yields a unique solution, whose behavior can be controlled outside of the training data. Neural networks are nonconvex ML models and are often harder to interpret and generalize outside of the training manifold.
Given the Cartesian coordinates of a molecule with N atoms, x ∈ R 3N , we define the collective coordinates by the encoding
where y ∈ R m and m is a small number. In ML terms, the coordinates y define a latent space.
In molecular simulations, the collective coordinates are often defined to describe the slowest processes of the system (63). In general, the mapping (or encoding) of E can be highly nonlinear. If the energy function U (x) associated with the atomistic representation is known, an important quantity to compute-for instance, to connect with experimental measurements-is the free energy of the system as a function of the collective variables. The free energy is defined as
where μ Y (y) is the marginal distribution on y of the equilibrium distribution μ(x) given by Equation 1:
The integral in Equation 7 is in practice impossible to compute analytically for high-dimensional systems, and several methods have been developed for its numerical estimation by enhancing the sampling of the equilibrium distribution in MD simulations of the system. The definition of the free energy can also be formulated as a learning problem: The aim is to optimize the parameters θ of a free energy function model F (y, θ ) such that Equations 6 and 7 are satisfied to a good approximation. Fulfilling these equations is usually referred to as enforcing thermodynamic consistency.
Using methods that estimate the free energy F λ (or its gradient F λ ) at a given set of points in the collective variables space y λ (λ = 1, . . . , M), one can use the free energy loss,
or the free energy gradient loss,
with these free energy estimates to reconstruct the entire surface F (y,θ ). Both kernel regression (28) and neural networks (30) have been used for this purpose. ML has also been used in combination with enhanced sampling methods to learn the free energy surface on the fly (25,29,(64)(65)(66)(67)(68)(69).
An alternative way to learn Equation 6 is by using force matching (70,71). It has been demonstrated (72) that, given the forces xU of the atomistic system collected along an MD trajectory, x t , t = 1, . . . , T, the gradient of a free energy model F (y,θ ) that best satisfies Equation 6 also minimizes the force-matching loss,
8 . where the term f y lmf is the local mean force,
that is, the projection of the atomistic force xU on the collective variable space through the mapping E.
In practice, the estimator in Equation 8 is very noisy: Because of the dimensionality reduction from x to y, multiple realizations of the projected force f y lmf can be associated with the same value of the collective coordinates y, and the minimum of the loss function in Equation 8 cannot go to zero. It can be shown by invoking statistical estimator theory that this loss can be broken down into bias, variance, and noise terms (20).
The use of coarse-grained models of complex molecular systems, such as proteins, presents an attractive alternative to the use of atomistic models that are very expensive to simulate (71,73).
The design of a coarse-grained model for a system with N atoms into a reduced representation with n effective beads starts with the definition of a mapping similar to Equation 5, where now y ∈ R 3n are the coordinates of the coarse-grained beads. In this case, the mapping is usually a linear function, as in general the beads are defined as a subset or a linear combination of sets of atoms in the original system, in a way that allows one to keep some information about the geometry of the molecule. For instance, in protein systems, a coarse-grained mapping can be defined by replacing all the atoms in a residue with a bead centered on the corresponding Cα atom. There is at present no general theory to define the optimal coarse-graining mapping for a specific system. A few methods have been proposed in this direction by optimizing the mapping E to preserve certain properties of the original system. Examples include the definition of a system's dynamical assembly units (74) or the use of an autoencoder to minimize the reconstruction error (75).
Once the coarse-graining mapping is given, several strategies exist to define the model energy function, either to match experimental observables (top-down) or to reproduce specific properties of the atomistic system (bottom-up) (71). If one wants to define a coarse-grained model that is thermodynamically consistent with the original model (Equations 6 and 7), then the definition of the energy function associated with a coarse-grained molecular model can be seen as a special case of the free energy learning discussed in the previous section (70). The effective energy of the coarse-grained model is the free energy defined by Equation 6, where y are now the coarse variables. Therefore, once the coarse-graining map is defined, the definition of the effective potential can also be seen as a learning problem. In particular, the force-matching loss function given by Equation 8 can be used to train the effective energy of the model from the atomistic forces. For a general linear mapping y = Ex where matrix E ∈ R 3n×3N clusters N atomic coordinates into n bead coordinates, the expression for the loss, Equation 8, reduces to
where the instantaneous coarse-grained force is given by f y (x) = (EE ) -1 Ef (x) of the atomistic forces f (x) (Equation 3). This loss function has been used to design coarse-grained force fields for different systems with both kernel methods (18) and deep neural networks (19,20).
Kinetics are the slow part of the dynamics. Due to the stochastic components in the MD integrator, for any trajectory emerging from a configuration x t at time t, there is a probability distribution of finding the molecule in configuration x t+τ at a later time t + τ :
We can express the transition density (Equation 10) as the action of the Markov propagator in continuous space and by its spectral decomposition (76,77):
The spectral decomposition can be read as follows: The evolution of the probability density p(x) can be approximated as the superposition of functions ψ. A second set of functions, φ, is required to compute the amplitudes of these functions.
In general, Equation 11 is a singular value decomposition with left and right singular functions φ k , ψ k and true singular values σ * k (77). The approximation is then a low-rank decomposition in which the small singular values are discarded. For the special case that dynamics are in thermal equilibrium, detailed balance (Equation 21) holds, and Equation 11 is an eigenvalue decomposition with the choices
Hence, Equation 11 simplifies: We need only one set of functions, the eigenfunctions ψ k . The true eigenvalues λ * k are real valued and decay exponentially in time τ with characteristic relaxation times t * i that are directly linked to kinetic experimental observables (78,79). The approximation in Equation 11 is due to truncating all terms with relaxation times shorter than t * n . A quite established approach is to learn molecular kinetics (Equation 11) from a given trajectory data set. To obtain a low-dimensional model of the molecular kinetics that is easy to interpret and analyze, this usually involves two steps: (a) finding a low-dimensional latent space representation of the collective variables, y = E(x), using the encoder E; and (b) learning a dynamical propagator P in that space:
A common approach in MD, but also in other fields such as dynamical systems and fluid mechanics, is to seek a pair (E, P) such that P is a small matrix that propagates state vectors in a Markovian (memory-less) fashion (77,78,(80)(81)(82)(83)(84)(85)(86)(87)(88). This is motivated by the spectral decomposition of dynamics (Equation 11): If τ is large enough to filter fast processes, a few functions are sufficient to describe the kinetics of the system, and if E maps to the space spanned by these functions, P can be a linear, Markovian model. If no specific constraints are imposed on P, then the minimum regression error of y t+τ = Py t , the variational approach of Markov processes (VAMP) (77,89), and maximum likelihood will all lead to the same estimator (90),
using the latent space covariance matrices
If E performs a one-hot encoding that indicates which state the system is in, then the pair (E, P) is called a Markov state model (78,(80)(81)(82)(83)(84), and P is a matrix of conditional probabilities to be in state j at time t + τ given that the system was in state i at time t.
Learning the embedding E is more difficult than learning P, as optimizing E by maximum likelihood or minimal regression error in latent space y leads to a collapse of E to trivial, uninteresting solutions (90). This problem can be avoided by following a variational approach to approximating the leading terms of the spectral decomposition (Equation 10) (77,89). The variational approach for conformation dynamics (89) states that for dynamics obeying detailed balance ( 21), the eigenvalues λ k of a propagator matrix P via any encoding y = E(x) are, in the statistical limit, lower bounds of the true λ * k . The VAMP variational principle is more general, as it does not require detailed balance (see Equation 21) and applies to the singular values σ k :
As VAMP is the more general principle, we can use it to define the loss function for estimating molecular kinetics models:
If dynamics obey detailed balance, we can use C τ τ = C 00 and plug in a symmetric estimate for C 0τ (26,90).
MD time steps are on the order of one femtosecond (10 -15 s), while configuration changes governing molecular function are often rare events that may take 10 -3 to 10 3 s. Even when the potential and the forces generated by single-protein folding and unfolding are evaluated quickly, simulating these processes by direct MD simulation may take years to centuries on a supercomputer. To avoid this sampling problem, ML methods can be employed to learn generation of equilibrium samples from μ(x) more efficiently, or even to generate statistically independent samples in one shot.
Learning to sample probability distributions is the realm of generative learning (91)(92)(93)(94). In the past few years, directed generative networks, such as variational autoencoders (92), generative adversarial networks (93), and flows (94,95), have received a particular surge of interest. Such networks are trained to transform samples from an easy-to-sample probability distribution, such as Gaussian noise, into realistic samples of the objects of interest. These methods have been used to draw photorealistic images (96,97), generate speech or music audio (98), and generate chemical structures to design molecules or materials with certain properties (46)(47)(48)(49).
If the aim is to learn a probability distribution from sampling data (density estimation) or to learn to efficiently sample from a given probability distribution [Boltzmann generation (41)], one typically faces the challenge of matching the probability distribution of the trained model with a reference.
Matching probability distributions can be achieved by minimizing probability distances. The most commonly used probability distance is the Kullback-Leibler (KL) divergence, also called relative entropy between two distributions q and p: KL(q p) = q(x) log q(x)log p(x) dx. 16.
In Boltzmann generation (41), we aim to efficiently sample the equilibrium distribution μ(x) of a many-body system defined by its energy function u(x) (Equation 1). We can learn the parameters θ of a neural network that generates samples from the distribution p X (x; θ) by making its generated distribution similar to the target equilibrium distribution. Choosing q ≡ p X (x) and p ≡ μ(x) in Equation 16, we can minimize the KL θ q X μ X with the energy loss (EL):
As the network samples from an energy surface u(x; θ) =log p X (x; θ) + const (Equation 1), this loss is performing energy matching, similar to what was discussed in Section 2.1. In order to evaluate the loss (Equation 17), we need not only to be able to generate samples x ∼ p X (x; θ) from the network but also to be able to compute p X (x; θ) for every sample x (see Section 4.5 for an example). The reverse case is density estimation. Suppose we have simulation data x, and we want to train the probability distribution p X (x; θ) to resemble the data distribution. We choose q(x) as the data distribution and p ≡ p X (x; θ) in Equation 16and exploit that E x∼data log q(x) is a constant that does not depend on the parameters θ. Then we can minimize the loss:
18.
This loss is the negative likelihood (NL) that the model generates the observed sample; hence, minimizing it corresponds to a maximum likelihood approach. Likelihood maximization is abundantly used in ML; in this review we discuss it in Section 4.5.
Compared to working with classical ML problems such as image classification, we have a decisive advantage when working with molecular problems: We know a lot of physical principles that restrict the possible predictions of our machine to the physically meaningful ones.
Let us start with a simple example to illustrate this. We are interested in predicting the potential energy and the atom-wise forces of the diatomic molecule O 2 with positions x 1 , x 2 in a vacuum (without external forces). Independent of the details of our learning algorithm, physics tells us that a few general rules must hold:
1. The energy is invariant when translating or rotating the molecule. We can thus arbitrarily choose the positions x 1 = (0, 0, 0); x 2 = (d, 0, 0), and the energy becomes a function of the interatomic distance only:
2. Energy is conserved. The energy U and the force of the molecule are related by F(x) = -∇U (x). Now we can compute the components of the force as
∂d , 0, 0). 3. Identical particles are indistinguishable. The energy is unchanged if we exchange the labels 1 and 2 of the identical oxygen atoms.
In ML, there are two principal approaches when dealing with such invariances or symmetries: (a) data augmentation and (b) building the invariances into the ML model.
Data augmentation means we learn invariances by heart by artificially generating more training data and applying the known invariances to it. For example, given a training data point for positions and energy/force labels, (x; U , f ), we can augment this data point with translation invariance of energy and force by adding more training data (x + x; U , f ) with random displacements x. Data augmentation makes the ML model more robust and helps us predict the invariances approximately. It is an important ML tool, as it is easy to do, although for certain invariances it is conceptually difficult or computationally expensive to hardwire those invariances into the ML model. However, data augmentation is statistically inefficient, as additional training data are needed, and inaccurate, because a network that does not have translation invariance hardwired into it will never predict that the energy is exactly constant when translating the molecule. This inaccuracy may lead to unphysical and potentially catastrophic predictions when such an energy model is inserted into an MD integrator.
The more accurate, statistically efficient, and elegant approach to incorporating physical constraints is to directly build them into the ML model. Doing so involves accounting for two related aspects that are both essential to making the learning problem efficient. First are equivariances:
The ML model should have the same invariances and symmetries as the modeled physics problem. Second is parameter sharing: Whenever there is an invariance or symmetry, this should be reflected by sharing model parameters in different parts of the network.
If we can hardwire the physical symmetries into the ML structure, we can reduce the dimensionality of the problem. In the example of the O 2 molecule above, the energy and force are oneand six-dimensional functions defined for each point x in six-dimensional configuration space. However, when we use the invariances described above, we have to learn only a one-dimensional energy function over a one-dimensional space, and we can compute the full force field from it. The learning problem has become much simpler because we now learn only in the manifold of physically meaningful solutions (Figure 2a).
An important concept related to invariance is equivariance. A function is equivariant if it transforms the same way as its argument. For example, the force is defined by the negative gradient of the energy, -∇U (x). If we rotate the molecule by applying the rotation R•, the energy is invariant, but the force is equivariant as it rotates in the same way as x (Figure 2b):
Equivariances are closely linked to convolutions in ML. For example, standard convolutions are translation invariant: Each convolution kernel is a feature detector that is applied to each pixel neighborhood (99). When that feature is present in the image, the convolved image will show a signal in the corresponding position (Figure 2c). When translating the input, the convolved image translates in the same way (Figure 2c). We briefly review below common invariances and equivariances useful for applications to molecular simulations and discuss strategies to incorporate them in an ML algorithm.
The force is equivariant to rotation but invariant to translation (Figure 2b; Section 3.4):
-∇U (Rx + T) = -R∇U (x).
Rototranslational invariance can be achieved by transforming the configuration x into rototranslationally invariant features, such as intramolecular distances or angles. Equivariance of the force can then be achieved by computing the force explicitly by a network layer that computes the derivatives with respect to x, as it is done, for example, in SchNet (100) and CGnet (coarse-graining neural network) (20) (Section 4.3). Note that periodic systems, such as crystals and explicit-solvent boxes, have translational but not rotational invariances and equivariances.
Physical quantities, such as QM energies, are invariant if we exchange the labels of identical atoms-for example, carbons. As the number of possible permutations increases exponentially with the number of identical particles, trying to learn permutation invariance by data augmentation is hopeless. Permutation invariance can be built into the ML model by choosing a functional form to compute the quantity of interest that is permutation invariant (see Reference 101 for the general conditions). Following the pioneering work in References 102-105, a specific choice that is common for networks that compute extensive quantities such as potential energy U (x) is to model them as a sum,
where U i is the contribution of the energy by the ith atom in its chemical environment (7,11,62,100). To account for the multibody character of QM, U i (x) must generally also be a multibody function. Equation 19 can be implemented by using separate networks for the individual contributions U i and adding up their results (Figure 3). The force resulting from Equation 19 is automatically permutation equivariant.
a Coordinates b Coordinates H2O
Total energy
Atom-specific neural network Oxygen network
Behler-Parrinello networks for learning quantum energies. (a) A Behler-Parrinello network computes the atomic energy of a single atom i. The system coordinates are mapped to rototranslationally invariant features describing the chemical environment of atom i. (b) A molecular system such as H 2 O is composed by employing a network copy for each atom. The parameters are shared between networks for same elements. The total energy is given by the sum of atomic energies, introducing permutation invariance.
Classical MD force fields define bonded interactions by assigning atoms to nodes in a bond graph, and thus, exchanging of individual atom pairs no longer preserves the energy. However, the energy is still invariant to the exchange of identical molecules, such as solvent, and it is important to take that into account for coarse graining (20) and generating samples from the equilibrium density (41).
A simple alternative to building permutation invariance into the ML function is to map all training and test data to a reference permutation. This can be done efficiently by so-called bipartite graph matching methods such as the Hungarian method (106) that are frequently used in recent learning models (10,41,107).
Closed physical systems are energy conserving, which implies that the force field is defined by the gradient of the energy (Equation 3). When learning potential energy, it can be a great advantage to use force information, because each N-atom configuration is associated with only one energy but 3N forces, which may result in superior data efficiency when force labels are used during learning (10,61). If we use supervised learning for coarse graining with thermodynamic consistency, we can only use forces, as no labels for the free energies are available (Section 2.3).
In these examples, we have labeled training data [x, f (x)] i . Using a network that directly predicts the forces would not guarantee that Equation 3 holds. Therefore, when Equation 3 holds, it is common to learn an energy U (x) and then compute the force in the network by using a gradient layer (20, 100) (Section 4.3). An alternative to ensure Equation 3 is gradient-domain ML (61).
In statistical mechanics (both equilibrium and nonequilibrium), we are interested in the probability of events. An important principle is thus probability conservation-that is, the sum over all events is 1. A common approach to encode the probability of classes or events with a neural network is to use a SoftMax output layer, where the activation of each neuron can be defined as
where j runs over all neurons in the layer. In this representation, u can be seen as a vector of energies giving rise to the Boltzmann probabilities y i (u); y i (u) are nonnegative, because of the exponential functions, and they sum up to 1.
In Markov state models of molecular kinetics, we would like to obtain a Markov transition matrix P that is stochastic; that is, p ij ≥ 0 for all elements, and j p ij = 1 for all i. If the encoder E(x) uses a SoftMax, then the estimator in Equation 14will result in a transition matrix P that conserves probability ( j p ij = 1 for all i). SoftMax can be exploited in VAMPnets to simultaneously learn an embedding E(x) from configurations to metastable states and a Markov transition matrix P (Section 4.4) (24).
Detailed balance connects thermodynamics and dynamics. In a dynamical system that evolves purely in thermal equilibrium-that is, without applying external forces-the equilibrium distribution μ(x) and the transition probability p(y | x) are related by the following symmetry:
21.
Thus, the unconditional probabilities of forward and backward trajectories are equal. Enforcing Equation 21 in kinetic ML models ensures the spectral decomposition (Equation 11) is real valued, which is useful for many analyses.
To learn a kinetic model that obeys detailed balance, estimators other than Equation 14 must be used; these typically enforce the unconditional transition probabilities p(x, y) = μ(x)p(y | x) to be symmetric (see Reference 90 for details).
A decisive advance in the performance of neural networks for classical computer vision problems such as text or digit recognition came with going from fully connected dense networks to convolutional neural networks (99). Convolution layers are equivariant, thus helping with the detection of an object independent of its location (Figure 2c), but the real efficiency of convolutional neural networks is due to parameter sharing.
In a classical dense network, all neurons of neighboring layers are connected and have independent parameters stored in a weight matrix W l . This leads to a problem with high-dimensional data. If we were to process relatively small images, say 100 × 100 = 10 4 pixels, and associate each pixel with a neuron, a single dense neural network layer l will have 10 4 × 10 4 = 10 8 parameters in W l . Not only would this be demanding in terms of memory and computing time but a network with so many independent parameters would also likely overfit and not be able to generalize to unknown data.
Convolutions massively reduce the number of independent parameters. A convolutional layer with a filter w acting on a one-dimensional signal x l-1 computes, before applying bias and nonlinearities,
In terms of an image, convolution applies the same filter w to every pixel neighborhood; in other words, all pixel transformations share the same parameters. Additionally, the filters w are usually much smaller than the signal x (often 3 × 3 for images), so each convolution has only a few parameters. For molecules, we extend this idea to continuous convolutions that use particle positions instead of pixels (Section 4.2).
Besides the sheer reduction of parameters, parameter sharing-for example, via convolution layers-is the keystone of transferability across chemical space. In convolutions of a grayscale image, we apply the same filters to every pixel, which implies translational equivariance but also means the same rules apply to all pixels. If we have multiple color channels, we have different channels in the filters as well, so pixels in the same color channels behave the same. In molecules we can apply the same idea to particle species. When one is learning energies from QM data, for example, every chemical element should be treated the same; that is, one should use the same convolution filters to sense its chemical environment. This treatment gives us a building-block principle that allows us to train on one set of molecules and make predictions for new molecules.
In this section, we present specific methods and neural network architectures that have been proposed to tackle the ML problems discussed in Section 2.
Behler-Parrinello networks are one of the first applications of ML in the molecular sciences (7). They aim at learning and predicting potential energy surfaces from QM data, and they combine all of the relevant physical symmetries and parameter sharing for this problem (Section 3).
First, the molecular coordinates are mapped to rototranslationally invariant features for each atom i (Figure 3a). In this step, the distances of neighboring atoms of a certain type, as well as the angles between two neighbors of certain types, are mapped to a fixed set of correlation functions that describe the chemical environment of atom i. These features are the input to a dense neural network that outputs one number, the energy of atom i in its environment. By design of the input feature functions, this energy is rototranslationally invariant. Parameters are shared between equivalent atoms-for example, all carbon atoms have the same network parameters for computing their atomic energies-but since the chemical environments will differ, their energies will differ. In a second step, the atomic energies are summed over all atoms of the molecule (Figure 3b). This second step, combined with parameter sharing, achieves permutation invariance, as explained in Section 3.4.2. Transferability is achieved due to parameter sharing but also because the summation principle allows one to grow or shrink the network to molecules of any size, including sizes that were never seen in the training data.
A related approach is Deep Potential net (62), where each atom is treated in a local coordinate frame that has the rotational and translational degrees of freedom removed.
Behler-Parrinello networks are traditionally trained by energy matching (Section 2.1) (7, 56), but they can be trained with force matching if a gradient layer is added to compute the conservative force (Section 2.1; Equation 3). The Behler-Parrinello method has been further developed in the ANI network, for example, by extending it to more advanced functions involving two neighbors (11). While Behler-Parrinello networks have mainly been used to make predictions of the same molecular system in order to run MD simulations unaffordable by direct ab initio QM MD (56), ANI has been trained on density functional theory and coupled-cluster data across a large chemical space (11,12,108).
One of the first deep learning architectures to learn to represent molecules or materials was the family of deep tensor neural networks (DTNNs) (14), with its recent addition SchNet (58,109). While in kernel-based learning methods (2,110) chemical compounds are compared in terms of prespecified kernel functions (8,(111)(112)(113), DTNN and its extension SchNet learn a multiscale representation of the properties of molecules or materials from large data sets.
DTNNs were inspired by the language processing approach word-to-vec (114), where the role of a word within its grammatical or semantic context is learned and encoded in a parameter vector. Likewise, DTNNs learn a representation vector for each atom within its chemical environment (Figure 4b, left). DTNN' s tensor construction algorithm then iteratively learns higher-order representations by first interacting with all pairwise neighbors, for example, extracting information implemented in the bond structure (Figure 4b, middle). By stacking such interaction layers deep, DTNNs can represent the structure and statistics of multibody interactions in deeper layers. As DTNNs are end-to-end trained to predict certain QM quantities, such as potential energies, they learn the representation that is relevant for the task of predicting these quantities (115).
SchNet (58) uses a deep convolutional neural network (116,117). Classically, convolutional neural networks were developed for computer vision using pixelated images, so they use discrete convolution filters. However, the particle positions of molecules cannot be discretized on a grid because QM properties such as the energy are highly sensitive to small position changes, such as the stretching of a covalent bond. For this reason, SchNet introduced continuous convolutions (58, 100), which are represented by filter-generating neural networks that map the rototranslationally invariant interatomic distances to the filter values used in the convolution (Figure 4b, right).
DTNNs and SchNet have both reached highly competitive prediction quality both across chemical compound space and across configuration space in order to simulate MD. In addition to their prediction quality, their scalability to large data sets (109) and their ability to extract novel chemical insights by means of their learned representation make the DTNNs increasingly popular research tools.
As mentioned in Section 2.3, ML has been used to define coarse-grained models for molecular systems. Both kernel methods (18) and deep neural networks (19,20) have been designed to minimize the force-matching loss, Equation 9, for given coarse-graining mappings for specific systems.
In both cases, it has been shown that the incorporation of physical constraints is crucial to the success of the model. The training data are obtained by means of atomistic molecular dynamic simulations, and regions of the configurational space that are physically forbidden, such as configurations with broken covalent bonds or overlapping atoms, are not sampled and not included in the training. Without additional constraints, the machine cannot make predictions far away from the training data and thus will not reliably predict that the energy should diverge when approaching physically forbidden regions.
Excluding high-energy states such as broken bonds or colliding atoms is different from enforcing physical symmetries as described in Section 3.4. Rather, it is about enforcing the correct asymptotic behavior of the energy when going toward an unphysical limit. CGnets proposed to achieve this by learning the difference from a simple prior energy that was defined to have the correct asymptotic behavior (20) (Figure 5a). The exact form of this prior energy is not essential for success, as the CGnet can correct the prior energy where training data are available. In Reference 20, the prior energy consisted of harmonic terms for bonds and angles of coarse-grained particles whose equilibrium values and force constants were obtained with Boltzmann inversion, as well as excluded volume terms in the form of (σ /r) c , where r is the interparticle distance and σ , c are hyperparameters.
As in Behler-Parrinello networks and SchNet, CGnet predicts a rototranslationally invariant energy as the first layer transforms the Cartesian coordinates into internal coordinates such as distances and angles (Figure 5a). Furthermore, CGnet predicts a conservative and rotationequivariant force field as the gradient of the total free energy with respect to input configuration x is computed self-consistently by the network (see Figure 5a). The network is trained by minimizing the force-matching loss of this prediction (Equation 8).
Figure 5b,c shows an application of CGnets to the coarse graining of the miniprotein chignolin, in which all solvent molecules are coarse grained away and the atoms of each residue are mapped to the corresponding Cα atom (Figure 5b). MD simulations performed with the force field predicted by the CGnet predict a free energy surface that is quantitatively similar to the free energy surface of the all-atom simulations, and they resolve the same metastable states (folded, unfolded, and misfolded). In contrast, a spline-based coarse-grained model where only two-body terms are included in the energy function cannot reproduce the all-atom free energy surface and does not even predict that folded and unfolded are separated metastable states. These results clearly illustrate the importance of multibody interactions in the coarse-grained energy-for example, surface or volume terms that can describe implicit solvation. While the spline model can be dramatically improved by adding suitable terms to the list of features, this is not necessary when using a deep neural network that automatically learns the required multibody terms. Similar conclusions can be obtained by using Gaussian approximation potentials as the ML model to capture multibody terms in coarse-grained energy functions (18).
VAMPnets (24) were introduced to replace the complicated and error-prone approach of constructing Markov state models by (a) searching for optimal features; (b) combining them into a low-dimensional representation y, for example, via time-lagged independent component analysis (118); (c) clustering y; (d) estimating the transition matrix P; and (e) coarse graining P. VAMPnets instead use a single end-to-end learning approach in which all of these steps are replaced by a deep neural network. This is possible because with the variational approach for conformational dynamics and VAMP principles (Section 2.4) (77,89), loss functions are available that are suitable to train the embedding E(x) and the propagator P simultaneously (see Section 2.4; Equation 13).
VAMPnets contain two network lobes representing the embedding E(x). These networks transform the molecular configurations found at a time delay τ along the simulation trajectories (Figure 6a). VAMPnets can be trained by minimizing the VAMP loss (Equation 15), which is meaningful for dynamics both with and without detailed balance (77). VAMPnets may, in general, use two distinct network lobes to encode the spectral representation of the left and right singular functions [which is important for nonstationary dynamics (119,120)]. Extended dynamic mode decomposition with dictionary learning (121) uses an architecture similar to that of VAMPnets but is optimized by minimizing the regression error in latent space. To avoid collapsing to trivial embeddings, such as constant functions (see Section 2.4), a suitable regularization must be employed (121).
While hyperparameter selection can be performed by minimizing the variational loss (Equation 15) on a validation set (122)(123)(124)(125), it is important to test the performance of a kinetic model on timescales beyond the training timescale τ . We can use the Chapman-Kolmogorov equation to test how well the learned model predicts longer times:
23.
A common way to implement this test is to compare the leading eigenvalues λ i (τ ) of the left-and right-hand sides (81,84,126). In Reference 24, parameters were shared between the VAMPnet nodes, and a unique embedding E(x) was thereby learned. When detailed balance is enforced while computing P(τ ) (Equation 21), the loss function automatically becomes a variational approach for conformational dynamics score (89). In this case, the embedding E(x) encodes the space of the dominant Markov operator eigenfunctions (24). This feature was extensively studied in state-free reversible VAMPnets (26).
To obtain a propagator that can be interpreted as a Markov state model, Mardt et al. VAMPnet and application to alanine dipeptide. (a) A VAMPnet (24). The network includes an encoder E that transforms each molecular configuration x t to a latent space of slow reaction coordinates y t and is trained on pairs (y t , y t+τ ) sampled from the molecular dynamics simulation using the VAMP score (77) (127,128). As a result, the propagator computed by Equation 14 conserves probability and is almost a transition matrix (Section 3.4.4), although it may have some negative elements with small absolute values.
The results described in Reference 24 (see, e.g., Figure 6) were competitive with and sometimes surpassed the state-of-the-art handcrafted Markov state model analysis pipeline. Given the rapid improvements in training efficiency and accuracy of deep neural networks seen in a broad range of disciplines, we expect end-to-end learning approaches such as VAMPnets to dominate the field eventually.
In Reference 31, a deep generative Markov state model was proposed that, in addition to the encoder E(x) and the propagator P, learns a generative part that samples the conditional distribution of configurations in the next time step. The model can be operated in a recursive fashion to generate trajectories to predict the system evolution from a defined starting state and to propose new configurations. The deep generative Markov state model was demonstrated to provide accurate estimates of the long-time kinetics and to generate valid distributions for small MD benchmark systems.
Boltzmann generators were introduced in Reference 41 to learn to sample equilibrium distributions (Equation 1). In contrast to standard generative learning, a Boltzmann generator does not attempt to learn the probability density from data but is trained to efficiently sample μ(x) ∝ e -u(x) using the dimensionless energy u(x) as an input. A Boltzmann generator consists of two parts:
1. a generative model that is trained to propose samples from a probability distribution p X (x) that is similar to μ(x) and that allows us to evaluate p X (x) (up to a constant) for every x, and 2. a reweighting procedure that takes proposals from p X (x) and generates unbiased samples from μ(x).
Boltzmann generators use a trainable generative network F zx that maps latent space samples z from a simple prior-for example, a Gaussian normal distribution-to samples x ∼ p X (x). Training is done by combining the energy-based training using the KL divergence (Equation 17) and maximum likelihood (Equation 18). For both training and reweighting, we need to be able to compute the probability p X (x) of generating a configuration x. This can be achieved by the change-of-variables equation if F zx is an invertible transformation (Figure 7a) (94,95). Such invertible networks are called flows due to the analogy of the transformed probability density with a fluid (94,129). In Reference 41, the nonvolume preserving transformations RealNVP were employed (129), but the development of more powerful invertible network architectures is an active field of research (97,130,131). By stacking multiple invertible blocks, a deep invertible neural network is obtained that can encode a complex transformation of variables.
Figure 7b-e illustrates the Boltzmann generator on a condensed matter model system that contains a bistable dimer in a box densely filled with repulsive solvent particles (Figure 7b). Opening or closing the dimer is a rare event, and it involves the collective rearrangement of the solvent particles due to the high particle density. Using short MD simulation in the open and closed states as an initialization, the Boltzmann generator can be trained to sample open, closed, and previously unseen transition states by generating Gaussian random variables in latent space (Figure 7c) and feeding them through the transformation F zx . Such samples have realistic structures and closeto-equilibrium energies (Figure 7d). Free energy differences can be computed by employing reweighting (Figure 7e). There is a direct relationship between the temperature of the canonical ensemble and the variance of the latent-space Gaussian of the Boltzmann generator (41). This allows us to learn to generate thermodynamics, such as the temperature-dependent free energy profiles, using a single Boltzmann generator (Figure 7e). Finally, as the latent space concentrates configurations of equilibrium probability around the origin, Boltzmann generators can be used to generate physically realistic reaction pathways by performing linear interpolations in latent space.
Despite rapid advances in the field of ML for molecular simulation, there are still significant open problems that need to be addressed, in all the areas discussed above.
To be practically useful, an ML model for both PES and atomic forces is needed that (a) can yield accuracy of 0.2-0.3 kcal/mol for the energy per functional group and about 1 kcal/(mol• Å) for the force per atom, (b) is not much more expensive to evaluate than classical force fields, (c) scales to large molecules such as proteins, and (d) is transferable to different covalent and noncovalent environments. Such a universal model does not exist yet.
Crucial steps toward items a and b have been recently taken by symmetrized gradient-domain machine learning (sGDML), a kernel-based approach to constructing molecular force fields (10,61,132,133). Currently, sGDML already enables MD simulations with electrons and nuclei treated at essentially exact QM level for molecules with up to 20-30 atoms.
Network-based approaches, such as SchNet and ANI, are better suited to items c and d, as they break down the energy in local interactions of atoms with their environment, enabling a buildingblock principle that is by design better scalable to molecules of different sizes and transferable across chemical space. However, these approaches do not reach the high accuracy in configuration space that sGDML does. Combining high accuracy in configuration and chemical space remains an active research topic.
The vast majority of approaches to make ML inference on molecular structures are based on local chemical information. Current neural networks for modeling molecular energies use the summation principle (e.g., Equation 19) to sum up local energies E i (x) of atom i with its neighbors. While multibody and long-ranged energies can be obtained by stacking multiple layers (14, 58, 100)-the working principle of deep convolution networks (116)-there are fundamental physical limits of this approach: Long-ranged interactions such as electrostatics cannot be cut off.
For classical point-charge models, long-ranged electrostatics methods have been developed, such as the Ewald summation method for periodic systems (134). One option is to combine shortranged ML models with such long-ranged electrostatics methods. To avoid double counting interactions, one must also predict atomic charges, which is an active field of research (135,136). An alternative option-currently unexplored territory-is to develop neural network structures for particle interactions that can compute long-ranged interactions by design.
In addition to electrostatics, van der Waals dispersion interactions can also have a substantial long-range character; that is, they can extend to separations of tens of nanometers or more in large molecular and nanoscale systems (137)(138)(139). Developing ML models that correctly treat the QM many-body nature of van der Waals interactions remains a difficult challenge to overcome (15).
With the availability of chemically transferable ML models that have quantum chemical accuracy, the next open problem is to sample metastable states and long-timescale kinetics. Although available ML models for predicting QM energies and forces are still significantly slower than force fields, the vast array of enhanced sampling methods and kinetic models (Sections 2.4 and 4.4) will likely allow us to explore the kinetics of quantum chemical systems on timescales of microseconds and beyond. A plethora of new physical insights that we cannot access with current MD force fields await us there. For example, what is the role of protonation dynamics in mediating protein folding or function?
An outstanding question in the design of coarse-grained models is that of transferability across chemical space. Bottom-up coarse-grained models are useful in practice if they can be parameterized on small molecules and then used to predict the dynamics of systems much larger than what is possible to simulate with atomistic resolution. It is not clear to what extent transferability of coarse-grained models can be achieved and how that depends on the coarse-graining mapping (140,141). Compared with the manual design of few-body free energy functionals, ML free energies can help with transferability, as they are able to learn the important multibody effects-for example, to model neglected solvent molecules implicitly (Section 4.3) (18)(19)(20).
It is natural to consider Behler-Parrinello networks or SchNet as a starting point for modeling transferable coarse-grained energies, but their application is nontrivial: It is a priori unclear what the interacting particles are in the coarse-grained model and how to define their types, as they are no longer given by the chemical element. Furthermore, these networks assume permutation invariance between identical particles, while classical MD force fields do not have permutation invariance of atoms within the same molecule. Therefore, particle network structures that can handle bonding graphs need to be developed.
While coarse-grained MD models may perform well in reproducing the thermodynamics of the atomistic system, they may fail in reproducing the kinetics. Existing approaches include adding fictitious particles (142) or training the coarse-grained model with spectral matching (143). There are indications that the kinetics can be approximated up to a global scaling factor in barrier-crossing problems when the barriers are well approximated (144), which could be achieved by identifying the slow reaction coordinates (63) and assigning more weight to the transition state in force matching or relative entropy minimization. This area of research is still underdeveloped.
Extensive properties such as potential energies can be well predicted across chemical space, as they can be conceptually broken down as a sum of parts that can be learned separately. This is not possible with intensive properties such as spectra or kinetics, and for this reason, the prediction of such properties is, as yet, far behind.
Generative networks, such as Boltzmann generators (Section 4.5), have been demonstrated to be able to generate physically realistic one-shot samples of model systems and proteins in implicit solvent (41). To scale to larger systems, we will need to build the invariances of the energy, such as the exchange of identical solvent particles, into the transformation and to include parameter sharing (Section 3.5), such that we can go beyond just sampling the probability density of one given system with energy u(x) and instead generalize from a data set of examples of one class of molecules, such as solvated proteins. To this end, equivariant networks with parameter sharing need to be developed for generative learning; these are, to date, not available.
Recently, the increasing popularity of explainable artificial intelligence methods (see, e.g., [145][146][147][148] has allowed us to gain insight into the inner workings of deep learning algorithms. In this manner, it has become possible to extract how a problem is solved by the deep model. This allows, for example, for the detection of so-called clever Hans solutions (147)-that is, nonsensical solutions relying on artifactual or nonphysical aspects in data. Combined with networks that learn a representation, such as DTNN/SchNet (14,100) and VAMPnets (24), these inspection methods may provide scientific insights into the mechanisms that give rise to the predicted physicochemical quantity and thereby fuel the development of new theories.
Annu. Rev. Phys. Chem. 2020.71:361-390. Downloaded from www.annualreviews.org Access provided by Rice University on 04/23/21. For personal use only.
www.annualreviews.org • Machine Learning for Molecular Simulation