Computing Committor Functions for the Study of Rare Events Using Deep Learning

TL;DR

This paper presents a deep learning-based method to efficiently compute committor functions, enabling the study of rare transition events in high-dimensional complex systems with rough energy landscapes.

Contribution

It introduces a novel computational approach combining deep learning, data sampling, and feature engineering to overcome challenges in calculating committor functions for realistic systems.

Findings

Achieves good performance on complex benchmark problems

Overcomes curse of dimensionality in rare event computation

Provides an alternative practical method for high-dimensional systems

Abstract

The committor function is a central object of study in understanding transitions between metastable states in complex systems. However, computing the committor function for realistic systems at low temperatures is a challenging task, due to the curse of dimensionality and the scarcity of transition data. In this paper, we introduce a computational approach that overcomes these issues and achieves good performance on complex benchmark problems with rough energy landscapes. The new approach combines deep learning, data sampling and feature engineering techniques. This establishes an alternative practical method for studying rare transition events between metastable states in complex, high dimensional systems.

Tables

Table 1. Table 1 : Comparison of the neural network approximation q θ subscript 𝑞 𝜃 q_{\theta} and the FEM solution q 𝑞 q for the extended Mueller potential. The errors are computed on 100 100 100 transition states sampled from the dynamics ( 25 ). The statistics of the error (mean ± plus-or-minus \pm deviation) is based on 10 10 10 independent runs. The network denoted by 10 ​ - ​ 20 ​ - ​ 20 ​ - ​ 1 10 - 20 - 20 - 1 10\text{-}20\text{-}20\text{-}1 contains 2 2 2 hidden layers with 20 20 20 nodes on each hidden layer; similarly for other networks.

sampling	data size	network	RMSE	MAE
method
artificial temperature	${10}^5$	$10\text{-}20\text{-}1$	$0.0607\pm 0.0201$	$0.0488\pm 0.0176$
		$10\text{-}20\text{-}20\text{-}1$	$0.0636\pm 0.0288$	$0.0519\pm 0.0240$
		$10\text{-}20\text{-}20\text{-}20\text{-}1$	$0.0520\pm 0.0157$	$0.0425\pm 0.0135$
	$2\times {10}^5$	$10\text{-}20\text{-}1$	$0.0456\pm 0.0070$	$0.0368\pm 0.0058$
		$10\text{-}20\text{-}20\text{-}1$	$0.0427\pm 0.0079$	$0.0345\pm 0.0065$
		$10\text{-}20\text{-}20\text{-}20\text{-}1$	$0.0416\pm 0.0111$	$0.0338\pm 0.0099$
	$4\times {10}^5$	$10\text{-}20\text{-}1$	$0.0261\pm 0.0062$	$0.0213\pm 0.0054$
		$10\text{-}20\text{-}20\text{-}1$	$0.0353\pm 0.0076$	$0.0297\pm 0.0074$
		$10\text{-}20\text{-}20\text{-}20\text{-}1$	$0.0309\pm 0.0100$	$0.0256\pm 0.0095$
meta- dynamics	${10}^4$	$10\text{-}20\text{-}1$	$0.0737\pm 0.0169$	$0.0572\pm 0.0135$
	$2\times {10}^4$	$10\text{-}20\text{-}1$	$0.0532\pm 0.0094$	$0.0441\pm 0.0081$
	$4\times {10}^4$	$10\text{-}20\text{-}1$	$0.0345\pm 0.0087$	$0.0285\pm 0.0081$

Equations

$$p(x)=\frac{1}{\mathcal{Z}}{e}^{-\beta V(x)},$$

$$\dot{x}(t)=-\nabla V(x(t))+\sqrt{2{\beta }^{-1}}\eta (t),$$

$${\tau }_A(x)$$

$${\tau }_B(x)$$

$$q(x)=\text{Prob}\{{\tau }_B(x)<{\tau }_A(x)\},$$

$$\{\begin{array}{l}\nabla V\cdot \nabla q-{\beta }^{-1}\Delta q=0,x\in \Omega \setminus (A\cup B),\\ q(x)=0,x\in \partial A;q(x)=1,x\in \partial B,\end{array}$$

$$\underset{q}{\min }\frac{1}{Z}{\int }_{\Omega \setminus (A\cup B)}|\nabla q(x){|}^2{e}^{-\beta V(x)}dx,$$

$$q(x)=0,x\in \partial A;q(x)=1,x\in \partial B.$$

$$q(x)=\left(1-{\chi }_A(x)\right)\left[(1-{\chi }_B(x))\tilde{q}(x)+{\chi }_B(x)\right],x\in \Omega \setminus (A\cup B),$$

$$\chi_{A}(x)=\left\{\begin{array}[]{ll}1,&x\in\partial A,\\ 0,&x\in\Omega\setminus A^{\epsilon},\end{array}\right.\quad\chi_{B}(x)=\left\{\begin{array}[]{ll}1,&x\in\partial B,\\ 0,&x\in\Omega\setminus B^{\epsilon},\end{array}\right.$$

$$\begin{array}{r}{A}^{ϵ}=\left\{x\in \Omega :\underset{y\in A}{\inf }|x-y|\le ϵ\right\},B^{ϵ}=\left\{x\in \Omega :\underset{y\in B}{\inf }|x-y|\le ϵ\right\}.\end{array}$$

$$q(x)=\tilde{q}(x),\text{for }x\in \Omega \setminus (A^{ϵ}\cup B^{ϵ}).$$

$$\underset{\tilde{q}}{\min }\frac{1}{Z}{\int }_{\Omega \setminus (A\cup B)}|\nabla q(x){|}^2{e}^{-\beta V(x)}dx,$$

$$q_{\theta }(x)=\left(1-{\chi }_A(x)\right)\left[(1-{\chi }_B(x)){\tilde{q}}_{\theta }(x)+{\chi }_B(x)\right],x\in \Omega \setminus (A\cup B),$$

$$\tilde{q}(x)=F_{L+1}\circ F_L\circ \cdots \circ F_2\circ F_1(x),$$

$$q(x)\approx f(z_1(x),...,z_m(x)).$$

$$\underset{\theta }{\mathrm{argmin}}\frac{1}{Z}{\int }_{\Omega \setminus (A\cup B)}|{\nabla }_x{q}_{\theta }(x){|}^2{e}^{-\beta V(x)}dx,$$

$$\frac{1}{Z}{\int }_{\Omega \setminus (A\cup B)}|{\nabla }_x{q}_{\theta }(x){|}^2{e}^{-\beta V(x)}dx={\text{E}}_{X\sim p_{\beta }}\left[|{\nabla }_x{q}_{\theta }(X){|}^2\right].$$

$${\text{E}}_{X\sim p_{\beta }}\left[|{\nabla }_x{q}_{\theta }(X){|}^2\right]\approx \frac{1}{N}\sum _{k=1}^N|{\nabla }_x{q}_{\theta }(X^{(k)}){|}^2,$$

$$\underset{\theta }{\mathrm{argmin}}{\int }_{\Omega \setminus (A\cup B)}|{\nabla }_x{q}_{\theta }(x){|}^2{e}^{-(\beta -{\beta }^{\prime })V(x)}{p}_{{\beta }^{\prime }}(x)dx,$$

$$\underset{\theta }{\mathrm{argmin}}{\int }_{\Omega \setminus (A\cup B)}|{\nabla }_x{q}_{\theta }(x){|}^2{e}^{\beta V_G(x)}{p}_G(x)dx,$$

$$V_{G,t}(x)=\sum _{t^{\prime }=0,\tau ,2\tau ,...}^{t^{\prime }<t}w\exp \left(-\sum _{i=1}^d\frac{(S_i(x)-S_i(x_{t^{\prime }}){)}^2}{2{\sigma }_i^2}\right),$$

$$V(x)=V_m(x_1,x_2)+\frac{1}{2{\sigma }^2}\sum _{i=3}^{10}{x}_i^2,x\in {\mathbb{R}}^{10}$$

$$V_m(x_1,x_2)$$

$${\chi }_A(x)=\frac{1}{2}-\frac{1}{2}\tanh [1000(|(x_1,x_2)-a{|}^2-(r+0.02{)}^2)],x\in {\mathbb{R}}^{10}$$

$$\dot{x}=-\nabla \left(V(x)+V_q(x)\right)+\sqrt{2{\beta }^{-1}}\eta ,$$

$$=\sqrt{\frac{1}{N_s}\sum _{k=1}^{N_s}{\left(q_{\theta }(x^{(k)})-q(x^{(k)})\right)}^2},$$

$$=\frac{1}{N_s}\sum _{k=1}^{N_s}\left|q_{\theta }(x^{(k)})-q(x^{(k)})\right|,$$

$$A$$

$$B$$

Full text

I Introduction

Understanding transition events between metastable states is of great importance in the applied sciences. Well-known examples of the transition events include nucleation events during phase transitions, conformational changes of bio-molecules, dislocation dynamics in crystalline solids, etc. The long time scale associated with these events is a consequence of the disparity between the effective thermal energy and typical energy barrier of the systems. The dynamics proceeds by long waiting periods around metastable states followed by sudden jumps from one state to another. For this reason, the transition event is called rare event. The main objective in the study of rare events is to understand the transition mechanism, such as the transition pathway and transition states. Some numerical methods have been proposed for this purpose, among which the well-known ones include the nudged elastic band method, Jónsson, Mills, and Jacobsen (1998) the string method, E, Ren, and Vanden-Eijnden (2002, 2007); Ren et al. (2005) the action-based method, Olender and Elber (1996) and the transition path sampling technique, Bolhuis et al. (2002); Dellago, Bolhuis, and Geissler (2002) accelerated molecular dynamics, Voter (1997) etc.

One object that plays an important role in understanding the transition event is the committor function. This is a function which is defined in the configuration or phase space and describes the progress of the transition. Most of the interesting information regarding the transition can be extracted from the committor function. E, Ren, and Vanden-Eijnden (2005); Ren et al. (2005); E and Vanden-Eijnden (2010) For example, the transition states lie in the region where the committor value is around $1/2$; the transition path ensemble and the transition rate can be computed as well, based on the committor function and the equilibrium probability distribution of the system. Thus, developing efficient numerical methods to compute the committor function is an important problem for understanding rare events.

The committor function has a very simple mathematical description - it satisfies the backward Kolmogorov equation. However, it is very difficult to compute in practice, due to the curse of dimensionality. For example, the equation needs to be solved in the configuration space of dimension $3N$ for a molecule consisting of $N$ atoms. Traditional numerical approaches, such as finite difference or finite element methods, are computationally prohibited even when $N=2$ or 3, and obviously out of the question for realistic physical systems. Alternative methods have been proposed. In the transition path sampling technique, Bolhuis et al. (2002); Dellago, Bolhuis, and Geissler (2002) the committor function is computed using Monte Carlo methods. In the finite-temperature string method, Ren et al. (2005) the problem is reformulated into an equivalent variational form, which is minimized in a particular function space. In the work of Ref. Lai and Lu, 2018, the Kolmogorov equation for the committor function is solved using a point cloud discretization.

More recently, it is proposed to compute the committor function using neural networks. Khoo, Lu, and Ying (2019) Satisfactory numerical results were obtained for a model problem in $10$ dimensions. The success of this approach hinges on the availability of data, especially in the transition state region where the committor function changes sharply from 0 to 1. As the term “rare event” suggests, such data is rarely available. In Ref. Khoo, Lu, and Ying, 2019, the data was generated by solving the underlying Langevin dynamics at the physical temperature. While this works well when the physical temperature is high in which case transitions are easily observed, it becomes less efficient when the temperature is low, i.e. the case of rare events.

In the current work, we propose to use importance sampling techniques to overcome this difficulty. We consider two methods. In the first one, we generate the data using the Langevin dynamics at an artificial temperature. This temperature is high enough so that the thermal energy becomes close or even comparable to the energy barrier. In this case, the transition event between the metastable states becomes less rare or even frequent. This enables us to collect sufficient amount of data in the transition state region. The difference between the physical temperature and the artificial temperature is accounted for by the likelihood ratio in the objective function to be minimized. In the second method, we sample the data using metadynamics. The relevant potential wells on the energy surface are first filled using localized Gaussian functions. This effectively lowers the energy barrier between the metastable states. The data can then be sampled efficiently following the Langevin dynamics on the modified potential.

In addition to these sampling methods, in this work we also introduce new methods to impose the boundary conditions for the committor function, and employ collective variables as the input features for the neural network.

The paper is organized as follows. The background and problem formulation are presented in Section II. In Section III, we introduce the key ingredients of our method, including the method of imposing the boundary conditions, data sampling techniques, and the neural network architecture. Numerical results for two examples are presented in Section IV. Finally, we draw conclusions in Section V.

II Background

The typical starting point of transition modeling is the specification of a potential energy function $V:\Omega\subset\mathbb{R}^{n}\rightarrow\mathbb{R}$, which takes as inputs the microscopic configuration $x$ of the system (e.g. the positions of the constituent atoms of a bio-molecule). The subset $\Omega$ is the configuration space under consideration. Through standard statistical mechanics arguments, the probability density of the system’s configuration is determined by the potential energy function $V$ via the Boltzmann-Gibbs distribution

[TABLE]

where $\beta=1/k_{B}T$ is the inverse temperature and $\mathcal{Z}=\int_{\Omega}e^{-\beta V(x)}dx$ is the normalization factor, or partition function.

To study dynamical properties, one may consider the noisy gradient flow induced by $V$ in the form of the over-damped Langevin equation:

[TABLE]

where $\eta$ is a white noise. One can check that Eq. (2) has the invariant distribution in Eq. (1). Through this dynamical model, it is clear how one can then introduce the notion of stability of the system: local minima, or a collection of local minima, of $V$ correspond to metastable configurations; at low temperatures such that the thermal energy is much lower than the energy barrier, the system will remain in these configurations for exponentially long times before making transition to another such configuration. Our principal goal is to study the transition dynamics between two metastable configurations.

For two given distinct metastable regions $A,B\subset\Omega$, $A\cap B=\varnothing$, consider the first hitting times

[TABLE]

The committor function $q:\Omega\rightarrow[0,1]$ is defined by

[TABLE]

i.e. $q(x)$ is the probability that the system initiated at $x$ first reaches $B$ rather than $A$. It can be shown that much of the vital information regarding transition pathways and rates can be extracted from the committor function. E, Ren, and Vanden-Eijnden (2005); E and Vanden-Eijnden (2010)

Obviously, $q$ takes the value 0 in $A$ and 1 in $B$. In the domain $\Omega\setminus(A\cup B)$, $q$ satisfies the backward Kolmogorov equation with Dirichlet boundary conditions, Gardiner (1997)

[TABLE]

where $\nabla=(\frac{\partial}{\partial x_{1}},...,\frac{\partial}{\partial x_{n}})$ is the gradient operator, $\Delta=\sum_{i=1}^{n}\frac{\partial^{2}}{\partial x_{i}^{2}}$ is the Laplace operator, $\partial A$ and $\partial B$ denote the boundary of $A$ and $B$, respectively. In addition, we impose the boundary condition $\nabla q\cdot\mathbf{n}=0$ on $\partial\Omega$. Eq. (5) has an equivalent variational formation. Specifically, the committor function is also the solution to the variational problem

[TABLE]

where $Z=\int_{\Omega\setminus(A\cup B)}e^{-\beta V(x)}dx$, subject to the boundary conditions

[TABLE]

The main challenge in using Eq. (5) to compute $q$ is the curse of dimensionality: as the dimension $n$ of the configuration $x$ increases, the computational complexity of classical finite difference and finite element methods increase exponentially, thereby prohibiting the use of these methods for realistic models. Although the variational formulation somewhat ameliorate this issue, at low temperatures (large $\beta$), which is the regime of interest for rare transition events, solving the minimization problem in (6) becomes more challenging as the measure $Z^{-1}e^{-\beta V}$ becomes more “singular”, i.e. the density is more concentrated near the minima of the potential energy.

In the next section, we introduce our proposed method of combining deep learning and data sampling to efficiently compute $q$ for high-dimensional systems and at low temperatures.

III Methods

The numerical method is based on the variational formulation (6)-(7). We first reduce this constrained minimization problem into an unconstrained problem by choosing a particular form for $q(x)$. Then we use a deep neural network to parameterize the function to be minimized, thereby convert the minimization problem into an unsupervised learning problem. This is followed by the discussion of the key issue of the method - the sampling of training data for the neural network.

III.1 Imposing the Boundary Conditions

In the earlier work, Khoo, Lu, and Ying (2019) the boundary conditions were imposed by minimizing a modified functional with an additional penalty term. While this is straightforward to implement, the introduction of an additional penalty parameter may require a careful tuning. Moreover, it requires sampling additional data on the boundaries. We proceed differently by introducing a particular form of $q$ that naturally satisfies the boundary conditions (7). Specifically, we consider the composite form

[TABLE]

where $\chi_{A}$ and $\chi_{B}$ are smooth functions such that $\chi_{A}(x)|_{\partial A}=1$, $\chi_{A}(x)|_{\partial B}=0$ and $\chi_{B}(x)|_{\partial B}=1$. One can easily verify that $q(x)$ in Eq. (8) satisfies the boundary conditions $q(x)|_{\partial A}=0$ and $q(x)|_{\partial B}=1$.

A convenient choice for $\chi_{A}$ and $\chi_{B}$ is to use the mollified indicator functions

[TABLE]

where $A^{\epsilon}$ and $B^{\epsilon}$ are two sets expanded from $A$ and $B$, respectively

[TABLE]

Away from $\partial A$ and $\partial B$, $\chi_{A}$ and $\chi_{B}$ changes smoothly from 1 to 0 in a region of width $\epsilon$, respectively. With this choice, $q$ agrees with $\tilde{q}$ outside $A^{\epsilon}$ and $B^{\epsilon}$,

[TABLE]

The constrained minimization problem (6)-(7) now reduces to an unconstrained problem

[TABLE]

where $q(x)$ takes the form (8).

III.2 Parameterization of the Committor Function

We use a deep neural network to parameterize the function $\tilde{q}$, hence approximating the committor function $q(x)$ as

[TABLE]

where $\{\theta\}$ are trainable variables of the neural network. Specifically, given input vector $x\in\Omega\setminus(A\cup B)$, the neural network with $L$ hidden layers is defined as

[TABLE]

where “$\circ$” denotes the composition of functions, $F_{i}(\xi)=\Phi_{i}(W_{i}\xi+b_{i})$ for $i=1,...,L+1$ and $\Phi_{i}$’s, $W_{i}$’s and $b_{i}$’s refer to the activation functions, weight matrices and bias vectors, respectively. Since the range of the committor function is [0, 1], we use a sigmoid activation for the output layer, i.e. $\Phi_{L+1}(\xi)=1/(1+\exp(-\xi))$.

Very often, only a few coarse-grained variables, called the collective variables and denoted by $(z_{1}(x),z_{2}(x),...,z_{m}(x))$ with $m\ll n$, play a major role in the transition event. For example, in conformational changes of bio-molecules, it is often that only a few torsion or bond angles are sufficient to characterize the transitions between metastable states. In this situation, the committor function can be well approximated by a function of these collective variables

[TABLE]

This is a form of dimensional reduction or feature engineering, and the choice of $\left\{z_{i}\right\}$ often allows one to build some degree of physical knowledge into the parameterization of $q$. To that end, we shall make use of collective variables in the design of the first input transformation layer of the neural network. The architecture of the network is summarized in Fig. 1.

With the neural network parameterization, the minimization problem (12) can be posed as an unsupervised learning problem

[TABLE]

where $q_{\theta}(x)$ takes the form (13).

III.3 Sampling Training Data for the Neural Network

The objective function in Eq. (16) is in the form of an expectation

[TABLE]

A conventional approach to solve this type of learning problem is to use the stochastic gradient descent (SGD) algorithm, in which the expectation is approximated by sample average

[TABLE]

where $X^{(k)}$’s are independent samples from the distribution $p_{\beta}(x)=Z^{-1}e^{-\beta V(x)}$.

The sampling of $p_{\beta}$ at low temperature $T$ (large $\beta$) poses a serious challenge. Naive sampling using the dynamics (2) at low temperatures yields very few transition events and hence the distribution $p_{\beta}$ is not explored well enough to solve (16). Specifically, with very few observed transition events, the majority of the sampled data will be clustered near the metastable states $A$ and $B$, with very few distributed in the region of our interest, i.e. the transition state region which lies in between $A$ and $B$. This will lead to a poor estimate of the committor function.

The difficulty is caused by the disparity between the thermal energy $k_{B}T$ and the energy barrier between the metastable states. Next, we propose two sampling methods to overcome this difficulty. The first one is to use an artificially increased temperature, and the other one is based on metadynamics which modifies the potential to lower the energy barrier.

III.3.1 Sampling Data at Artificial Temperature

We consider another form of the problem (16)

[TABLE]

where $\beta^{\prime}=1/k_{B}T^{\prime}$, with $T^{\prime}>T$ and $p_{\beta^{\prime}}(x)={Z^{\prime}}^{-1}e^{-\beta^{\prime}V(x)}$ in which $Z^{\prime}$ is the normalization factor.

Solving the new problem (19) requires sampling $p_{\beta^{\prime}}$, for example, by simulating the Langevin dynamics (2) at the high temperature $T^{\prime}$. This becomes much more efficient when $T^{\prime}$ is sufficiently large, as energy barriers are much easier to cross at higher temperatures. In essence, this is a form of importance sampling, and the function $e^{-(\beta-\beta^{\prime})V(x)}$ is the likelihood ratio associated with the change of measure. Note that the choice of importance sampling measure $p_{\beta^{\prime}}$ is not arbitrary. Its purpose is to produce enough sample points in the a priori unknown transition region for us to estimate the committor function at low temperature $T$. Generic sampling measures (e.g. uniform) will not achieve this goal in moderately high dimensions.

III.3.2 Sampling Data Using Metadynamics

In this approach, an external potential $V_{G}$ is added to the potential $V$ to lower the energy barrier thereby facilitate transitions between the metastable states. Correspondingly, we consider the following equivalent form of the problem (16)

[TABLE]

where $p_{G}(x)=Z_{G}^{-1}e^{-\beta\left[V(x)+V_{G}(x)\right]}$ is the equilibrium distribution for the modified potential and $Z_{G}$ is the normalization factor. We use metadynamics to construct the external potential $V_{G}$.

Metadynamics was developed to enhance the sampling in molecular dynamics simulations and construct the free energy landscapes. Laio and Parrinello (2002); Barducci, Bonomi, and Parrinello (2011) The idea is to fill the potential wells containing the metastable states by depositing localized Gaussian functions following the dynamics. We denote the frequency of deposition by $\tau^{-1}$ and use Gaussian functions in collective variables. Then at time $t$, the external potential that was added to $V$ is given by

[TABLE]

where $w$ and $\sigma_{i}$’s are parameters that control the height and width of the Gaussian respectively, $x_{t}$ denotes the trajectory on the time-dependent potential $V+V_{G,t}$, and $\left\{S_{i}(x),i=1,\dots,d\right\}$ are some coarse-grained variables. Since the system evolves towards local minima of the modified potential, the basins containing the metastable states are filled after sufficiently long time; as a result, the transitions take place more frequently on the modified potential compared to the original system. Thus metadynamics provides an efficient way to sample transition events even at low temperatures.

In our simulation, we first fill the potential wells of $V$ using Eq. (21) up to a certain time $t^{*}$ when the system can easily hop over the potential barrier. This gives the external potential $V_{G}$ in Eq. (20): $V_{G}(x)=V_{G,t^{*}}(x)$. Afterwards, we sample data for the distribution $p_{G}(x)$ using the Langevin dynamics with the modified potential $V(x)+V_{G}(x)$. These data are used to solve the problem (20).

IV Numerical Examples

We demonstrate the effectiveness of the proposed numerical method using two benchmark problems: one is the Mueller potential extended to high dimensions, the other is the isomerization of alanine dipeptide.

IV.1 Extended Mueller Potential

We first consider the Mueller potential embedded in the $10$-dimensional ($10$D) space, Khoo, Lu, and Ying (2019)

[TABLE]

where $V_{m}(x_{1},x_{2})$ is the rugged Mueller potential in two dimensions (2D),

[TABLE]

and a harmonic potential is used in each of the other eight dimensions. The parameters $\gamma$ and $k$ control the roughness of the energy landscape, and $\sigma$ controls the scale of the quadratic terms. In this example, we use $\gamma=9,\ k=5,\ \sigma=0.05$. The other parameters are taken from Ref. Lai and Lu, 2018.

We first use the finite element method (FEM) to solve the backward Kolmogorov equation (5) for the 2D rugged Mueller potential $V_{m}$ on $\tilde{\Omega}=[-1.5,1]\times[-0.5,2]$ by the solver FreeFem++. Hecht (2012) The numerical solution is denoted by $q_{m}$. The discretization error in $q_{m}$, which is on the order of $10^{-5}$, is much smaller than the numerical error in the neural network approximation. Therefore, we neglect the discretization error and treat $q_{m}$ as the “exact” solution. Then the “exact” solution for the committor function in $10$D is given by $q(x)=q_{m}(x_{1},x_{2})$, $x\in\Omega=\{x:(x_{1},x_{2})\in\tilde{\Omega},x\in\mathbb{R}^{10}\}$. Fig. 2 (top panel) shows the contour plots of the 2D rugged Mueller potential and the committor function $q_{m}$ at $k_{B}T=10$ (the energy barrier from $A$ to $B$ is about $100$).

Next, we compute the neural network approximation to the committor function in the $10$D space using the method proposed in section III. The Mueller potential (23) has two local minima around $a=(-0.558,1.441)$ and $b=(0.623,0.028)$, respectively. We take the two metastable sets $A$ and $B$ as the cylinders centered at $(x_{1},x_{2})=a$ and $(x_{1},x_{2})=b$ respectively, each with radius $r=0.1$. The function $\chi_{A}(x)$ is constructed as

[TABLE]

and similarly for $\chi_{B}(x)$. These two functions satisfy (9) approximately.

The data are sampled using the two different methods discussed in section III.3. In the first one, we generate the data at the artificial temperature $k_{B}T^{\prime}=20$ by solving the Langevin equation (2) using the Euler-Maruyama scheme with the time step $\Delta t=10^{-5}$. In the second method, we use metadynamics to generate the data. We first bias the coordinates $x_{1}$ and $x_{2}$ (i.e. $S_{1}(x)=x_{1}$, $S_{2}(x)=x_{2}$ in Eq. (21)) by adding $2000$ Gaussian functions with height $w=5$ and width $\sigma_{1}=\sigma_{2}=0.05$ into the potential, one for every 500 time steps. Then a set of data are sampled by simulating the Langevin dynamics on the modified potential with the time step $\Delta t=10^{-5}$.

In both methods, we take one sample for every $100$ time steps, and only keep those data points with coordinates $x\in\Omega\setminus(A\cup B)$. Of these data, $70\%$ serves as the training dataset and the other $30\%$ serves as the validation dataset. The neural network used in this example is fully connected, and the hyperbolic tangent (tanh) function is used as the activation function in the hidden layers. We use the package TensorFlow Abadi et al. (2015) with Adam optimizer Kingma and Ba (2015) to train the network by minimizing (19) and (20) respectively at the physical temperature $k_{B}T=10$.

To quantitatively evaluate our results, we compare neural network approximation $q_{\theta}$ with the FEM result $q$. For rare events considered here, we are interested in the iso-surface of the committor function where $q\approx 1/2$, in particular, the region near the minima of $V$ on this surface. This region represents the bottleneck to the transition and is known as the transition state region. The committor function in this region is also the most difficult to compute, as data here is typically scarce. Therefore, we shall focus our evaluation of the numerical results in this region. To this end, we carry out the constrained sampling in the transition state region by simulating the dynamics

[TABLE]

where the additional potential $V_{q}(x)=\frac{1}{2}\kappa\left(q_{\theta}(x)-\frac{1}{2}\right)^{2}$ with $\kappa=3\times 10^{4}$ is to constrain the system on the $1/2$-isocommittor surface $\Gamma_{1/2}=\{x\in\Omega:q_{\theta}(x)=1/2\}$. After equilibration, we sample $100$ points: $\left\{x^{(k)}:k=1,\dots,100\right\}$. These points, projected on the $(x_{1},x_{2})$ plane, are shown in Fig. 2 (middle and lower panels). The middle panels shows the result obtained using artificial temperature, and the lower panel shows the result obtained using data sampled from metadynamics. The region where these points are clustered is the transition state region. Also shown in the two panels are comparisons of the $1/2$-isosurface of $q_{\theta}$ with that of $q$ obtained from the FEM calculation. Evidently these two agree well in the transition state region. Certain discrepancies occur away from the transition state region, due to the lack of training data there; nevertheless, those regions are irrelevant to the transition events.

We also compute the root-mean-square error (RMSE) and the mean-absolute error (MAE) between $q_{\theta}$ and the FEM solution $q$ at the sampled transition states $x^{(k)},k=1,...,100$,

[TABLE]

where $N_{s}=100$. The results for different dataset and difference choice of network structure are reported in Table 1. Each error shown in the table is computed from $10$ independent runs; each run includes sampling the data, training the neural network, and sampling the transition states. We observe that the numerical results are insensitive to the number of hidden layers, but the accuracy improves for larger set of training data. We also observe that the sampling method based on metadynamics is more efficient than the method of artificial temperature. To achieve roughly the same accuracy, the latter requires one order of more data as compared to the method based on metadynamics. This is due to the fact that a larger portion of training data from metadynamics are concentrated in the transition state region.

IV.2 Alanine Dipeptide

In this example, we study the isomerization process of the alanine dipeptide in vacuum at $T=$300\text{,}\mathrm{K}$$. The isomerization of alanine dipeptide has been the subject of several theoretical and computational studies, Apostolakis, Ferrara, and Caflisch (1999); Bolhuis, Dellago, and Chandler (2000); Ma and Dinner (2005); Ren et al. (2005) therefore it serves as a good benchmark problem for the proposed method.

The molecule consists of 22 atoms and has a simple chemical structure, yet it exhibits some of the important features common to biomolecules. Figure 3 shows the stick and ball representation of the molecule (upper panel) and its adiabatic energy landscape on the plane of the two torsion angles $\phi$ and $\psi$ (lower panel). The molecule has two metastable conformers $C_{7eq}$ and $C_{ax}$ located around (-$$,$$) and ($$,-$$), respectively. Accordingly, the metastable sets $A$ and $B$ are chosen as

[TABLE]

Our goal is to compute the committor function for transitions between $A$ to $B$ and sample the transition states at the temperature $T=$300\text{,}\mathrm{K}$$.

We use the two sampling methods discussed in section III.3 to generate the training data. In the method of artificial temperature, we raise the temperature to $T^{\prime}=$800\text{,}\mathrm{K} and use the package NAMD Phillips *et al.* ([2005](#bib.bib22)) to simulate the Langevin dynamics with the time step $\Delta t=$0.5\text{\,}\mathrm{fs}. In the method of metadynamics, we bias the two torsion angles $(\phi,\psi)$ (i.e. $S_{1}(x)=\phi$, $S_{2}(x)=\psi$ in Eq. (21)). We add one Gaussian function with height $w=$0.1\text{,}\mathrm{kcal}\text{,}{\mathrm{mol}}^{-1}$$ and width \sigma_{1}=\sigma_{2}=$$$ for every 500 time steps. In total, 10^{4}$Gaussian functions are added to fill the potential wells containing the metastable sets$A$and$B. To save the computational cost, Fiorin, Klein, and Hénin ([2013](#bib.bib23)) the biased potential and its gradient are accumulated and computed on a uniform grid with mesh size [math] for \phi$and$\psi$.

In Fig. 4, we show the data sampled at the temperature $T^{\prime}=$800\text{,}\mathrm{K} (middle panel) and from metadynamics on the modified potential at the physical temperature $T=$300\text{\,}\mathrm{K} (lower panel), respectively. These data are sampled from $10^{7}$ steps of Langevin dynamics, one for every 100 steps. We observe that the data sampled from metadynamics has better quality in the sense that a larger portion of the data are in the transition state region. As a comparison, we also run the dynamics for the same number of time steps on the original potential at the temperature $T=$300\text{,}\mathrm{K}, and plot the sampled data in the figure (upper panel). It is seen that at the temperature $T=$300\text{\,}\mathrm{K}, no transition from $C_{7eq}$ to $C_{ax}$ is observed within the simulation time; consequently, no training data is collected in the transition state region, the region of our interest. This demonstrates the benefit of our proposed sampling methods.

In the design of the network structure, we use feature engineering via collective variables. As conformational changes of bio-molecules can often be well described by a few collective variables such as torsion angles, we design the first layer of our network to extract a collection of such torsion angles from the molecule, whose sines and cosines are then fed into the subsequent hidden layers as extracted features. Note that these include the four torsion angles shown in Fig. 3 that were identified as being adequate in describing the transition in earlier work (e.g. Ref. Maragliano et al., 2006). But we do not a priori assume that we know this precise information, and a redundant description is supplied.

We train the network at $T=$300\text{,}\mathrm{K} using $2\times 10^{6}$ data points sampled at $T^{\prime}=$800\text{\,}\mathrm{K} and $10^{5}$ data points sampled from metadynamics, respectively. All the data are sampled outside the sets $A$ and $B$. Of these data, $70\%$ serves as the training data and $30\%$ are used for validation. With these data, we use the package TensorFlow Abadi et al. (2015) with Adam optimizer Kingma and Ba (2015) to minimize the objective functions in (19) and (20), respectively. The computation is terminated when the validation error no longer decreases.

For this high-dimensional problem, we cannot afford to compute the committor function using the FEM method as we did in the first example. In order to check the accuracy of the numerical results, in particular, whether the $1/2$-isocommittor surface really locates the transition states, we carry out the constrained Langevin dynamics simulation on the $1/2$-isocommittor surface at $T=$300\text{,}\mathrm{K}$$. Following the dynamics, we collect $100$ states. The committor values of these states are then computed directly using Langevin dynamics. Specifically, we generate $200$ trajectories initiating from each of these states with random initial velocities and estimate the probability (i.e. the committor value) of the system first reaching $B$ rather than $A$ (c.f. (4)). If adequate accuracy is achieved, we would expect the computed probabilities to cluster around $1/2$. We carry out two independent runs, including sampling data, training the network and sampling the transition states, for each of the two sampling methods. Fig. 5 shows the distribution of the committor values at the sampled transition states. Results in (a)-(d) and (e)-(f) are obtained using 9 and 41 torsion angles in the feature layer of the network, respectively. In all the results, the committor values cluster around $q=1/2$, indicating the transition states are correctly identified. We also tested neural networks with fewer nodes in the hidden layer and obtained similar results.

V Conclusion

In this paper, we introduced a method for computing the committor function at low temperatures. The committor function characterizes rare transition events between metastable states. The main idea of the method is to combine deep learning with efficient sampling methods in order to overcome the curse of dimensionality associated with realistic systems and the scarcity of transition data at low temperatures. We also incorporated collective variables into the network structure as a form of crude feature selection to improve the efficiency of learning.

We considered two sampling methods: the method of artificial temperature and the method based on metadynamics. In the numerical examples considered in this work, the method of metadynamics outperformed the sampling method at artificial temperature. The former produced training data of better quality in the sense that a larger portion of the data lie in the transition state region. As a result, less amount of data were needed to achieve the same accuracy as in the sampling method using higher temperatures. Nevertheless, metadynamics requires a suitable choice of collective variables to bias, while sampling at higher temperature does not need this information.

Our method is demonstrated to be effective on a relatively simple example involving the rugged Mueller potential, as well as a more complex benchmark example of the alanine dipeptide molecule. This provides an alternative approach to study complex systems with rough energy landscapes. We intend to apply the method to more complex systems in the future work.

Acknowledgements.

The work of Ren was supported in part by Singapore MOE AcRF grants R-146-000-232-112 (Tier 2) and R-146-000-267-114 (Tier 1), and the NSFC grant (No. 11871365).

reference