MCB111: Mathematics in Biology (Fall 2024)


week 11:

From the Markov stochastic process to chemical concentration changes

In this notes, we deduce how the differential equations that control the change of concentrations of the different species in a set of chemical reactions as the averages of the master equation describing a stochastic Markov process underlying the reactions.

In other words, we are going to go from a Markov birth/death master equation describing a stochastic process, to a set of ordinary differental equations describing the behavior of the compounds concentrations changes with time.

An example

The chemical reaction

We have a system that converts two compounds A, B into two different ones C, and D. The reaction is given by

\[A+ B\overset{k}{\longrightarrow} C+D\]

In the w11 lecture we indicated that the differential equations describing the change of the concentrations of the four compounds over time is

\[\begin{aligned} \frac{d A}{d t} = \frac{d B}{d t} &= -k\ A B\\ \frac{d C}{d t} = \frac{d D}{d t} &= +k\ A B. \end{aligned}\]

Here we are going to deduce this results from the master equation describing the underlying stochastic Markov process.

The master equation

The reaction

\[A+ B\overset{k}{\longrightarrow} C+D\]

indicates that instantaneously one molecule of the A and B compounds disappear and a molecule of the C and D compounds appears with probability \(k\).

If the system has \((A, B, C, D)\) molecules of the four different compounds, in one small time interval \(\Delta t\), the number of molecules changes to \((A-1, B-1, C+1, D+1)\) with probability \(k A B \Delta t\) because any combination of one \(A\) and one \(B\) molecule can be used in the reaction that results in one \(C\) and one \(D\) more molecules.

In a time interval \(\Delta t\), the system reacts according to the following probability distribution,

\[\begin{aligned} (A, B, C, D) &\overset{k A B \Delta t}{\longrightarrow} (A-1, B-1, C+1, D+1) \\ (A, B, C, D) &\overset{1- k A B \Delta t}{\longrightarrow} (A, B, C, D). \end{aligned}\]

Then, the master equation following the description in the w10 lecture can be written as

\[\begin{aligned} P(A, B, C, D\mid t+\Delta t)&= \left. \begin{array}{rl} + k\, (A+1) (B+1) \Delta t & P(A+1, B+1, C-1, D-1\mid t)\\ + \left[1-k A B \Delta t\right] & P(A, B, C, D\mid t) \end{array} \right. \end{aligned}\]

And in the limit of \(\Delta t\rightarrow 0\) we obtain the differential equation

\[\begin{aligned} \frac{d P(A, B, C, D\mid t)}{dt}&= \left. \begin{array}{l} k\, (A+1)\, (B+1) P(A+1, B+1, C-1, D-1\mid t)\\ - k\, A\, B\, P(A, B, C, D\mid t) \end{array} \right. \end{aligned}\]

The deterministic solution

By taking averages over the number of \(A\) particles we obtain

\[\begin{aligned} \sum_{A, B, C, D} A \frac{d P(A, B, C, D\mid t)}{dt}&= \left. \begin{array}{l} + k\, \sum_{A, B, C, D} A\, (A+1)\, (B+1)\, P(A+1, B+1, C-1, D-1\mid t)\\ - k\, \sum_{A, B, C, D} A^2\, B\, P(A, B, C, D\mid t) \end{array} \right. \end{aligned}\]

which can be re-written as

\[\begin{aligned} \sum_{A, B, C, D} A \frac{d P(A, B, C, D\mid t)}{dt}&= \left. \begin{array}{l} + k\, \sum_{A, B, C, D} (A-1)\, A\, B\, P(A, B, C, D\mid t)\\ - k\, \sum_{A, B, C, D} A^2\, B\, P(A, B, C, D\mid t) \end{array} \right. \end{aligned}\]

or

\[\sum_{A, B, C, D} A \frac{d P(A, B, C, D\mid t)}{dt} = - k\, \sum_{A, B, C, D} A\, B\, P(A, B, C, D\mid t).\]

so that

\[\frac{d \langle A\rangle_t}{dt} = - k\, \langle A\, B\rangle_t.\]

Similarly, one can find the results for the change with time of the concentrations of the other species B, C and D.

The general case

A general multiple chemical reaction

In the w11 lecture, we indicated that for \(R\) reactions involving \(N\) substrates \(X_1\ldots X_N\), with rate constants \(k_1,\ldots k_R\),

\[\begin{aligned} n^1_1\ X_1 + \ldots + n^1_N\ X_N &\overset{k_1}{\longrightarrow} p^1_1 X_1 + \ldots + p^1_N X_N\\ \vdots\\ n^R_1\ X_1 + \ldots + n^R_N\ X_N &\overset{k_R}{\longrightarrow} p^R_1 X_1 + \ldots + p^R_N X_N \end{aligned}\]

the differential equation controlling the average number of particles per volume unit are given by

\[\begin{aligned} \frac{d X_i}{d t} &= \sum_{r=1}^R \eta^r_i\ v_r, \end{aligned}\]

where \(\eta^r_i\) is the stoichiometry of species \(X_i\) in reaction \(r\),

\[\eta^r_i = p^r_i-n^r_i.\]

and \(\nu_r\) is the mass action law reaction \(r\),

\[\nu_r(X_1\ldots X_N) = k_r\, X_1^{n_1^r}\ldots X_N^{n_N^r}.\]

Here we are going to deduce this results directly from the master equation.

The master equation

The master equation following the description in the w10 lecture can be written as

\[\begin{aligned} P(X_1\ldots\,X_N\mid t+\Delta t)&= \left. \begin{array}{rl} \sum_{r=1}^R w_r(X_1-\eta^r_1\ldots\,X_N-\eta^r_N) \Delta t & P(X_1-\eta^r_1\ldots\,X_N-\eta^r_N\mid t)\\ + \left[1- w_T(X_1\ldots\,X_N) \Delta t\right] & P(X_1\ldots\,X_N\mid t), \end{array} \right. \end{aligned}\]

where the propensity \(w_T(X_1\ldots\,X_N)\) is defined as the sum of the propensities for all reactions

\[w_T(X_1\ldots\,X_N) = \sum_{r=1}^R w_r(X_1\ldots\,X_N).\]

And in the limit of \(\Delta t\rightarrow 0\) we obtain the differential equation

\[\begin{aligned} \frac{d P(X_1\ldots\,X_N\mid t)}{dt}&= \left. \begin{array}{l} + \sum_{r=1}^R w_r(X_1-\eta^r_1\ldots\,X_N-\eta^r_N) P(X_1-\eta^r_1\ldots\,X_N-\eta^r_N\mid t)\\ - w_T(X_1\ldots\,X_N) \, P(X_1\ldots\,X_N\mid t) \end{array} \right. \end{aligned}\]

The deterministic solution

By taking averages over the number of \(X_i\) particles we obtain

\[\begin{aligned} \sum_{X_1\ldots\,X_N} X_i \frac{d P(X_1\ldots\,X_N\mid t)}{dt}&= \left. \begin{array}{l} + \sum_{X_1\ldots\,X_N} \sum_{r=1}^R X_i\, w_r(X_1-\eta^r_1\ldots\,X_N-\eta^r_N)\, P(X_1-\eta^r_1\ldots\,X_N-\eta^r_N\mid t)\\ - \sum_{X_1\ldots\,X_N} X_i w_T(X_1\ldots\,X_N)\, P(X_1\ldots\,X_N\mid t) \end{array} \right. \end{aligned}\]

which can be re-written as

\[\begin{aligned} \sum_{X_1\ldots\,X_N} X_i \frac{d P(X_1\ldots\,X_N\mid t)}{dt}&= \left. \begin{array}{l} + \sum_{X_1\ldots\,X_N} \sum_{r=1}^R (X_i+\eta_i^r)\, w_r(X_1\ldots\,X_N)\, P(X_1\ldots\,X_N\mid t)\\ - \sum_{X_1\ldots\,X_N} X_i w_T(X_1\ldots\,X_N)\, P(X_1\ldots\,X_N\mid t) \end{array} \right. \end{aligned}\]

that is

\[\begin{aligned} \sum_{X_1\ldots\,X_N} X_i \frac{d P(X_1\ldots\,X_N\mid t)}{dt}&= \left. \begin{array}{l} + \sum_{X_1\ldots\,X_N} X_i \sum_{r=1}^R w_r(X_1\ldots\,X_N)\, P(X_1\ldots\,X_N\mid t)\\ + \sum_{X_1\ldots\,X_N} \sum_{r=1}^R \eta_i^r w_r(X_1\ldots\,X_N)\, P(X_1\ldots\,X_N\mid t)\\ - \sum_{X_1\ldots\,X_N} X_i w_T(X_1\ldots\,X_N)\, P(X_1\ldots\,X_N\mid t) \end{array} \right. \end{aligned}\]

The first and third term cancel out, and we obtain,

\[\sum_{X_1\ldots\,X_N} X_i \frac{d P(X_1\ldots\,X_N\mid t)}{dt} = \sum_{r=1}^R \eta_i^r \sum_{X_1\ldots\,X_N} w_r(X_1\ldots\,X_N)\, P(X_1\ldots\,X_N\mid t) .\]

so that

\[\frac{d \langle X_i\rangle_t}{dt} = \sum_{r=1}^R \eta_i^r \langle w_r(X_1\ldots\,X_N)\rangle_t.\]

Introducing the mass action law for a reaction as the average of the propensities,

\[\nu _r(X_1\ldots\,X_N) = \langle w_r(X_1\ldots\,X_N)\rangle_t,\]

then we obtain the final result

\[\frac{d \langle X_i\rangle_t}{dt} = \sum_{r=1}^R \eta_i^r\, \nu _r(X_1\ldots\,X_N).\]