Statistical Learning of Reduced Kinetic Monte Carlo
Models of Complex Chemistry from MolecularDynamics
Qian Yang, Carlos A. Sing-Long, Enze Chen, Muralikrishna Raju,
Evan J. Reed
Stanford University
number of molecules

CH4250MD
train 12 ps
train 25 ps
train 50 ps
train 100 ps

200
150
100
50

0200

400

600

time (ps)February 26, 2018
AI at SLAC Seminar

Molecular dynamics simulations (existing or new)contain a wealth of data.
I

1000s of atoms

I

1000s of reactions

I

100s of moleculesI

1000s of timesteps perpicosecond simulation

Can we do more with this complex, expensive DATA?
3

We propose using statistical learning to learn a chemical
reaction network from atomistic simulation data.

IWe discover that it’s possible to build kinetic Monte Carlo (KMC)
models from a single or few MD simulations that can extrapolate
the dynamics of the chemical system more than 10x in time and inchemical space.

IWe develop a new data-driven method that reduces thousands ofreactions to fewer than a hundred in a matter of minutes.

4Outline

Building KMC models from MD simulations
Extrapolation in Time and Chemical Space

Reduction from thousands to dozens of reactions

Summary5We project molecular dynamics data onto a model
consisting of elementary reactions and reaction rates.
CH4200150
10050H2100
mean

number of molecules

number of molecules

250

0200

400

mean

80
60
402060000200time (ps)400

600time (ps)C2H6+number of molecules

30
mean

25
20
1510

Elementary Reactions and Reaction Rates
500200400

600

time (ps)Fig. 2 Six independent molecular dynamics simulations of the same
system under the same thermodynamic conditions, resulting insomewhat different molecular concentration trajectories due to different
initial velocities. Each colored line corresponds to a projected molecularconcentration trajectory derived from a single molecular dynamics
simulation. The dotted black line corresponds to the mean of these
trajectories. We see that although the general trends are the sameacross simulations, the number of molecules at any given time can6fluctuate between simulations. The average difference (in the root meansquare sense) between each molecular concentration trajectory and theBond duration is used to control bias-variance tradeo↵
in our model
RMSE (molecules)

15CH4H2
C2H6

1050
10 -2

10 -1

10 0

10 1(ps)3500
original+signal+from+bond+length+criteria

all reactions
unimolecularbimoleculartrimolecularcomplexincreasing++!

bondednot+bonded
bondednot+bonded

number of reactions

3000
2500
2000
1500
1000
500
0time

00.1

0.2

0.3

0.4

0.5(ps)7

We model chemical reaction networks as a set ofelementary reactions and reaction rates governed by
the chemical master equation.
Reactions are random events. We associate with every reaction apropensity function, aj (x), such that aj (x)dt gives the probability ofthat reaction occurring in the time interval [t, t + dt).
III

Unimolecular reactions X1 ! products: aj (X) = kj X1

Bimolecular reactions 2X1 ! products: aj (X) = kj X1 (X11)

Bimolecular reactions X1 + X2 ! products: aj (X) = kj X1 X2

The chemical master equation is a system of ODEs that gives theprobability P(x, t) of being in a particular state X(t) = x in molecularconcentration space at time t.
8

We use Maximum Likelihood Estimation to estimatereaction rates
In a given MD simulation with M unique species and R unique reactions,
we observe at every timestep t:
1. vector of molecule counts X(t) = (x1 (t), ..., xM (t))
2. vector of reaction counts R(t) = (r1 (t), ..., rR (t))
For t 1, we can consider X(t) as a function of X(0) and all of theprevious R(t0 ) for t0 < t, so that our set of observations for timesteps 1
to T are(X(0), R(0), ..., R(T 1))
For N independent MD simulations, we estimate the vector of reaction
rate constants k = (k1 , ..., kR ) by maximizing the likelihood of theobservations:k⇤ = max L(k; (X(0), R(0), ..., R(T 1))1 , ..., (X(0), R(0), ..., R(T 1))N )k9

We simulate the chemical master equation using kinetic
Monte CarloGiven a reaction network with thousands of reactions and reaction rates,Gillespie Stochastic Simulation is equivalent to kinetic Monte Carlo andcan be used to simulate the chemical master equation exactly.
At each step:
I Choose the next reaction based on the propensity functions foreach reaction.
I Choose the time until the next reaction.

number of molecules

250
CH4H2C2H6MD

200150

Our KMC simulations take minutes
to model the MD data that tookweeks to generate.

1005000200400time (ps)

600

11

Outline

Building KMC models from MD simulations
Extrapolation in Time and Chemical Space

Reduction from thousands to dozens of reactions

Summary16

We show that the chemical reaction network can besignificantly reduced when describing target molecules.
CH4
number of molecules

250
molecular dynamics
kinetic Monte Carlo, 2000 reactions
kinetic Monte Carlo, 75 reactions

200150100
500200400600time (ps)Model reduction of complex systems is a formidable technical challenge:
I Sensitivity Analysis
I Mixed Integer Programming17

We represent chemical reaction networks as a linear
system in the set of reactions
Consider the toy reaction network:
A+B !C(1)

2C ! A(3)

A ! 2C(2)

Then according to the reaction rate equations, we haved[A]
= k1 [A][B] k2 [A] + k3 [C] ([C] 1)
dt
d[B]= k1 [A][B]dt
d[C]= k1 [A][B] + 2k2 [A] 2k3 [C] ([C] 1)
dt
We would like to choose reaction rate coefficients to match the first andsecond moments of the stochastic model at the sampled X(t). Themoments are nonlinear in molecular concentrations, but linear inreaction rate coefficients.

18

We represent chemical reaction networks as a linear
system in the set of reactions
23
AAt
1 4 t+⌧
E Bt+⌧ Bt 5 = µt+⌧
⌧Ct+⌧ Ct
|{z}
2yt

32
3
2
3
1
1
1
= k1 4 1 5 At Bt + k2 4 0 5 At + k3 4 0 5 Ct (Ct 1)
12223232
31
1 1
At B t 00k1
5 4 k2 50 54 0
At0=4 1 012200 Ct (Ct 1)k3|{z}|{z} | {z }RDtkA similar expression linear in the vector of ki ’s can be obtained for1cov (yt ) = HDt k⌧

19

We represent model reduction of chemical reaction
networks as a linear system subject to L1 regularization

23 2
3µ2 | X1
RD16 ⌃2 |X1 7 6 HDt 77 67166
7 = 6 ... 7 k...
6767
⌧ 4
µn+1 |Xn 5 4 RDn 5⌃n+1 |XnHDnsample nonlinear data

Then after scaling the system by k from the MLE, we have thenonnegative LASSO (least absolute selection and shrinkage operator):
minkAccbk2

s.t.kck1  ,c020