A kinetic Monte Carlo algorithm for Markovian dynamics of multi-state systems

Clear[k12, k13, k21, k23, k31, k32]

The general model we are trying to study consists of N discrete states, j =1,2, ..., N. These could be energy levels of a molecule, distinct conformations of a polymer chain, or adsorption sites on a surface, on which an atom is performig its random walk. The kinetics is described by 1-st order rate equations, with a rate constant (Subscript[k, ij]) for the transition i -> j. The rate constants are contained in the matrix  (Subscript[k, ij])  . Typically N would be very large but for now assume N=3.

n = 3

kMatrix = {{0, k12, k13}, {k21, 0, k23}, {k31, k32, 0}}

kMatrix //MatrixForm

If p = (p1, p2, ... pN)  is a vector describing the probabilities of being in states 1, ..., N, then it obeys the differential equation (standard 1st order kinetics):  d(Subscript[p, i])/dt=-(∑_(j=1)^N)(Subscript[k, ij])(Subscript[p, i])+(∑_(j=1)^N)(Subscript[k, ji])(Subscript[p, j]), or, in matrix form:  dp/dt = Kp, where K is a matrix or the form:   (Subscript[K, ii])=-(∑_(j=1)^N)(Subscript[k, ij])  and (Subscript[K, ij])=(Subscript[k, ji])

K = Transpose[kMatrix]; For[i = 1, i ≤ n, i++, K[[i, i]] = -Sum[kMatrix[[i, j]], {j, 1, n}]];

K // MatrixForm

Notice that the sum of the elements in each column is equal to zero.

The solution of the equation for the probability vector as a function of time can be written in terms of the matrix exponent: p(t) = exp(Kt)  p(0)

Clear[p]; p[t_, p0_] := MatrixExp[t K] . p0

For example, suppose we start in state 1 at t=0, then p0 = {1,0,0}

p0 = {1,0,0}

Let' also select a specific set of rates:

k12 = k13 = k21 = k31 = k23=k32=1

K // MatrixForm

Here is the time dependence of the solution:

p[t, p0]

Plot[{p[t, p0][[1]], p[t, p0][[2]], p[t, p0][[3]]}, {t, 0, 3}, PlotRange -> All]

A physically reasonable matrix K must satisfy 3 properties:

Property 1. It has one zero eigenvalue, λ=0

Eigenvalues[K]

This means that there is a vector (Subscript[p, 0])such that K(Subscript[p, 0])=0. Since dp/dt = Kp, the vector (Subscript[p, 0])is the steady-state or equilibrium solution for the kinetics

Property 2. Eigenvalues of K are real numbers. If the eigenvalues were complex, we would have oscillatory solutions for p(t). That would be very entertaining but it never happens (except in nonlinear systems that are not considered here). The origin of this property can be traced back to the fact that the rate matrix is "nearly symmetric". Specifically, the rates satisfy the detailed balance condition (Subscript[k, ij])(Subscript[p, i])=(Subscript[k, ji])(Subscript[p, j]), where (Subscript[p, i])~exp(-(Subscript[βG, i])) is the equilibrium probability for being in state i.

Property 3. K has no positive eigenvalues. If there's a positive eigenvalue, the solution would blow up at t  ⟶∞

So far there was no kinetic Monte Carlo involved. Now consider the same problem from a single-molecule perspective. We have one molecule that jumps among its possible states j=1, ..., N. If, at certain time, the molecule happens to be in a state i, during the next dt seconds it can make a transition to any other state j with a probability equal to (Subscript[k, ij]) dt.  As a result we have a random walk in the j-space, which is controlled by the K matrix. Kinetic Monte Carlo is a method for generating such a random walk on a computer.

suppose at certain time t my random walk has arrived at state i.  In KMC, one first generates the time tNext = t +Δt  when the next jump occurs. While I'm in state i, the overall probability to make a transition to ANY other state is determined by the total escape rate, which is equal to -(Subscript[K, ii]). More precisely, the probability of surviving in the state i during the next Δt seconds is given by exp((Subscript[K, ii])Δt). It turns out that to generate Δt that has such survival probaility all one needs is to solve the equation  exp((Subscript[K, ii])Δt) = r, where r is a uniform random number from 0 to 1.

Clear[survivalTime]; survivalTime[i_] := Log[Random[Real, {0, 1}]]/K[[i,i]]

survivalTime[2]

Now that we have determined the time to the next jump we need to decide what state we are going to jump to. The probability to jump to state j is proportional to the rate constant (Subscript[k, ij]) so all we need to select j according to these probabilities. This is achieved by the "roulette wheel" procedure. Here is a function that does it:

Clear[selectOne];
selectOne[probs_]:=Module[{lo=1,hi=Length[probs],mid,r},cumProb = FoldList[Plus, 0, probs] // Rest;  r=Random[] Last[cumProb]; While[lo≠hi-1,mid=Round[(lo+hi)/2];
If[r<cumProb[[mid]],hi=mid,lo=mid]];
If[r≤cumProb[[lo]],lo,hi]]

selectOne[{1, 1}]

Table[selectOne[{1, 1}], {40}]

given an unnormalized list of probabilities  probs = {p1, p2...pN} it will choose a number 1... N with a probability proportional to pN. For example, selectOne[{1, 1}] will choose1 or 2 with equal probabilities.  selectOne[{1, 1, 4}] will produce 1, 2, or 3 such that 3 will apear 4 times more likely than 1 or 2

Table[selectOne[{1, 1, 4}], {40}]

Now let's combine all this into an algorithm for a single KMC step that, given initial state i, tells us at after what time deltaT the next step will happen and what will be the new state j

Clear[kmcStep]; kmcStep[i_] := {survivalTime[i], selectOne[kMatrix[[i]]]}

kmcStep[2]

Finally, by repeatedly invoking this procedure we generate a KMC "trajectory"

initialState = 1

iter = 0; time = 0; state = initialState;trajectoryList = {};  While[iter < 20, blah = kmcStep[state]; time = time + blah[[1]]; nextState = blah[[2]];trajectoryList = Append[trajectoryList, {time, state}];  trajectoryList = Append[trajectoryList, {time, nextState}]; state = nextState; iter += 1; ]

trajectoryList;

Plot trajecory  (state as a function of time)

ListPlot[trajectoryList,Joined->True]