Introduction

Based on a set of ordinary differential equations (ODE) describing the kinetics of the cells differentiation, we designed a model to find the best differentiation rate.
First we developed a deterministic algorithm based on the ordinary differential equations solutions.
Then we find out that a stochastic algorithm could be an other solution to solve our problem.

For system involving large cell counts, the ordinary differential equations model give an accurate representation of the behavior. But with small cell counts, the stochastic and discrete method has a significant influence on the observed behaviour.
These reasons led us to write both a deterministic program based on the mass action law and a stochastic program based on the Gillespie’s stochastic simulation algorithm (SSA). With these two programs we obtain an accurate analysis of the vitamin production.

Deterministic model

Parameters

We conceived a simple model with the minimum number of parameters.
We find seven significant parameters for our model.

• \(t\) : fermentation period.
• \(k_{1}\) : rate constant of the mother cell differentiation.
• \(k_{2}\) : rate constant of the mother cell doubling time.
• \(k_{3}\) : rate constant of the differentiated cell doubling time.
• \(k_{4}\) : rate constant of the differentiated cell vitamin production.
• \(MC_{0}\) : initial number of mother cells in the medium.
• \(DC_{0}\) : initial number of differentiated cells in the medium.

Kinetic equations

We wrote four simple kinetic equations.
\[ \begin{align} MC \xrightarrow[]{k_{1}} DC \\ MC \xrightarrow[]{k_{2}} 2.MC \\ DC \xrightarrow[]{k_{3}} 2.DC \\ DC \xrightarrow[]{k_{4}} DC + Vitamin \end{align} \]

Formal mathematical solution

Translation in ordinary differential equations

We used the law of mass action to write the ordinary differential equations based on kinetic equations \((1)\), \((2)\), \((3)\) and \((4)\).
We found the following equations.
\[ \begin{align} \frac{dMC}{dt}(t) = (k_{2} - k_{1}).MC(t) \end{align} \] \[ \begin{align} \frac{dDC}{dt}(t) = k_{1}.MC(t) + k_{3}.DC(t) \end{align} \] \[ \begin{align} \frac{dVitamin}{dt}(t) = k_{4}.DC(t) \Rightarrow Vitamin(t) = k_{4}.\int_{0}^{t}{DC(t').dt'} \end{align} \]

Mathematical resolution

Simple ordinary differential equations resolution methods are used to find the solutions of the previous equations \((5)\), \((6)\), and \((7)\).
We expressed the time evolution of the mother cell, differentiated cell and vitamin number.
\[ \begin{align} MC(t) = MC_{0}.e^{(k_{2} - k_{1}).t} \end{align} \]
\[ \begin{align} DC(t) = \begin{cases} DC_{0}.e^{k_{3}.t} + MC_{0}.\frac{k_{1}}{k_{2}-k_{1}-k_{3}}.(e^{(k_{2} - k_{1}).t} - e^{k_{3}.t}) & \mbox{if } k_{1} \ne k_{3} - k_{2} \\ \\ (MC_{0}.(k_{2} - k_{3}).t + DC_{0}).e^{k_{3}.t} & \mbox{if } k_{1} = k_{3} - k_{2} \end{cases} \end{align} \]
\[ \begin{align} Vitamin(t) = \begin{cases} DC_{0}.\frac{k_{4}}{k_{3}}.(e^{k_{3}.t} - 1) + MC_{0}.\frac{k_{4}.k_{1}}{(k_{2}-k_{1}-k_{3}).(k_{2} - k_{1})}.(e^{(k_{2} - k_{1}).t} - 1) + MC_{0}.\frac{k_{4}.k_{1}}{(k_{2}-k_{1}-k_{3}).k_{3}}.(1 - e^{k_{3}.t}) & \mbox{if } k_{1} \ne k_{3} - k_{2} \\ \\ DC_{0}.\frac{k_{4}}{k_{3}}.(e^{k_{3}.t} - 1) + MC_{0}.\frac{k_{4}}{k_{3}}.(k_{3}.t - 1).(\frac{k_{2}}{k_{3}} - 1).(e^{k_{3}.t} - 1) + MC_{0}.\frac{k_{4}}{k_{3}}.(\frac{k_{2}}{k_{3}} - 1) & \mbox{if } k_{1} = k_{3} - k_{2} \end{cases} \end{align} \]

Vitamin optimization

Our goal is to optimize the vitamin production. We can only change three parameters : \(k_{1}\), \(MC_0\) and \(DC_0\).
\(t\), \(k_{2}\), \(k_{3}\) and \(k_{4}\) are constants. \(MC_0\) and \(DC_0\) are not relevant parameters. It seems logical that the more cells are in the medium, the more vitamin are producted.
However, it could be interesting to compare different \(\frac{MC_0}{DC_0}\) ratios but we prefer to focus on the rate constant \(k_{1}\).
We try to find the best \(k_{1}\). To this end we maximize the vitamin function numerically with MATLAB.
As an example, lets consider the following parameters.

• \(t = 10 \phantom{t} (hours)\)
• \(k_{2} = 0.66 \phantom{t} (hours^{-1})\)
• \(k_{3} = 0.33 \phantom{t} (hours^{-1})\)
• \(k_{4} = 1 \phantom{t} (hours^{-1})\)
• \(MC_0 = 5 \phantom{t} (cells)\)
• \(DC_0 = 0 \phantom{t} (cells)\)

We obtain the following graph with a simple MATLAB program available here.


As you can see, we are able to find the best \(k_{1}\) to optimize the vitamin production.
We noticed that two differents \(k_{1}\) exists to optimize the maximum number of differentiated cell \(DC\) or the maximum number of \(Vitamin\).
In our case, we chose of course the \(k_{1}\) that optimize the vitamin production ie \(k_{1} = 0.207\).

Deterministic evolution of mother cell, differentiated cell and vitamin numbers

We wrote a deterministic algorithm with MATLAB using the previous solutions \((8)\), \((9)\) and \((10)\). For those interested, the source code is available here.
In order to optimize the vitamin production, we use the same parameters as previously and set \(k_{1}\) with the previous resulting value ie \(k_{1} = 0.207\).
We obtained the following graph.


As you can see, the mother cell, differentiated cell and vitamin numbers follow an exponential law of time.
This result seems relevant. The model does not take into account the cells death and the nutrients present in the medium.

Stochastic model

Parameters

We used more parameters than the deterministic model. Indeed, because of the nature of the algorithm, we need a constant rate probability distribution.
In order to do a simulation with laboratory data, we modelize the time probability distribution with a normalized gaussian distribution using the normrnd MATLAB function.
For each time distribution \(i\), we use the following time probability distribution : \[ \begin{align} \mathcal P_{i}(t) = \frac{1}{\sigma_{i}.\sqrt{2\pi}}.e^{\frac{-(t - \mu_{i})^2}{2.\sigma_{i}^2}} \end{align} \] where

• \(\mathcal P_{i}(t)\) : time \(t\) probability.
• \(\mu_{i}\) : time mean.
• \(\sigma_{i}\) : time standard deviation.

Rate constants for the time distribution \(i\) is calculated using this formula : \[ \begin{align} k_{i} = \frac{ln(2)}{t} \end{align} \] The time \(t\) is a random variable with a probability \(\mathcal P_{i}(t)\) as explained above.
The vitamin number is computed by numerical integration using the rate constant \(k_{4}\).
We have ten parameters for this model.

• \(t\) : fermentation period.
• \(\mu_{1}\) : mother cell differentiation time mean.
• \(\sigma_{1}\) : mother cell differentiation time standard deviation.
• \(\mu_{2}\) : mother cell doubling time mean.
• \(\sigma_{2}\) : mother cell doubling time standard deviation.
• \(\mu_{3}\) : differentiated cell doubling time mean.
• \(\sigma_{3}\) : differentiated cell doubling time standard deviation.
• \(k_{4}\) : rate constant of the differentiated cell vitamin production.
• \(MC_0\) : initial number of mother cells in the medium.
• \(DC_0\) : initial number of differentiated cells in the medium.

Algorithm

Introduction

In the program, each cell is a structure containing three elements :

• the cell type :
  • 1 : mother cell.
  • 2 : differentiated cell.
• the next event type :
  • 1 : a mother cell become a differentiated cell.
  • 2 : a mother cell divide into two mother cells.
  • 3 : a differentiated cell divide into two differentiated cells.
• the time before the next event.

The time \(T_{i}\) for an event is chosen thanks to the equation \((11)\).
There is a competition between the first two events concerning the mother cell.
Either the mother cell differentiate, or it divide into two differentiated cells.
To do the choice, we compare \(T_{1}\) and \(T_{2}\). The smallest time is the next event time.

Step 1 : Initialize the cells

First of all we have to generate the first cells in accordance with \(MC_0\) and \(DC_0\).
For each cell we set the cell type, the next event and the time before the next event.

Step 2 : Choose the next event among all cells

We have an array containing all the cells. Each cell contain the time before his next event (divise or differentiate).
The next event time processed by the algorithm is the smallest cell next event time.

Step 3 : Do the next event

Now that we know the next event type and the concerned cell (only one cell is involved here), we compute this event.
As explained above, we have three choices.
If the cell is a mother cell, the cell can divise or differentiate.
If the cell is a differentiate cell, the cell can only divise.

We define a new event and time event for the cell and for the new cell in case of division.
We also performed the number of vitamin calculus by numerically integrating the number of \(DC\).
We use the following formula. \[ \begin{align} Vitamin(t) = Vitamin(t-T_{i}) + k_{4}.DC(t).T_{i} \end{align} \] where

• \(Vitamin(t)\) : number of vitamin at the instant \(t\)
• \(DC(t)\) : number of differentiated cells in the medium at the instant \(t\).
• \(T_{i}\) : time of the next event being performed.
• \(k_{4}\) : rate constant of the differentiated cell vitamin production.

Step 4 : Update the simulation time

For each cell except the concerned cell, define the cell next event time with the following formula : \[ \begin{align} T_{event cell} = T_{event cell} - T_{i} \end{align} \] where

• \(T_{i}\) : time of the next event being performed.
• \(T_{event cell}\) : cell next time event.

Step 5 : Go to step 2 until the simulation time is superior to the fermentation period.

We do this processus until the simulation time is superior to the fermentation period. Of course, the number of cells increase with an exponential law. The computational power increase too with the same law.
Do not try to processed a long simulation time if your computer is not powerful.

Step 6 : Show the results.

Here are a graph generated by the stochastic algorithm.
With the following distribution of event time probability.

Vitamin optimization

Like the deterministic model, our goal is to maximize the vitamin production.
In order to get relevant results, we average by doing a large number of simulations specified by the \(Averaging\) \(number\) parameter.
We can only change three parameters : \(mu_{1}\), \(MC_0\) and \(DC_0\).
\(t\), \(\sigma_{1}\), \(mu_{2}\), \(\sigma_{2}\), \(mu_{3}\), \(\sigma_{3}\) and \(k_{4}\) are constants.
We try to find the best \(mu_{1}\) to maximize the vitamin production.
As an example we can consider the following parameters.

• \(t = 5 \phantom{t} (hours)\)
• \(Averaging\) \(number = 500\)
• \(\sigma_{1} = 0.1 \phantom{t} (hours)\)
• \(\mu_{2} = 1.0502 \phantom{t} (hours)\)
• \(\sigma_{2} = 0.1 \phantom{t} (hours)\)
• \(\mu_{3} = 2.1004 \phantom{t} (hours)\)
• \(\sigma_{3} = 0.1 \phantom{t} (hours)\)
• \(k_{4} = 1 \phantom{t} (hours^{-1})\)
• \(MC_0 = 5 \phantom{t} (cells)\)
• \(DC_0 = 0 \phantom{t} (cells)\)

We obtained the following graph.
Here is a deterministic vitamin optimization.

We can compare the two results by superimposing the two graphs.

We performed cross-correlation and correlation coefficient between the two graphs with the xcorr and corr2 MATLAB functions.
As described in the official MATLAB website, the \(DC\) and \(Vitamin\) correlation coefficients are computed with the following formula.

\[ \begin{align} DC \phantom{t} correlation \phantom{t} coefficient = \frac{\sum\limits^{}_{k_{1}} (DC_{d}(k_{1}) - \overline{DC_{d}}).(DC_{s}(k_{1}) - \overline{DC_{s}})} {\sqrt{\left(\sum\limits^{}_{k_{1}} (DC_{d}(k_{1}) - \overline{DC_{d}})^2 \right).\left(\sum\limits^{}_{k_{1}} (DC_{s}(k_{1}) - \overline{DC_{s}})^2 \right)}} \end{align} \]
\[ \begin{align} Vitamin \phantom{t} correlation \phantom{t} coefficient = \frac{\sum\limits^{}_{k_{1}} (Vitamin_{d}(k_{1}) - \overline{Vitamin_{d}}).(Vitamin_{s}(k_{1}) - \overline{Vitamin_{s}})} {\sqrt{\left(\sum\limits^{}_{k_{1}} (Vitamin_{d}(k_{1}) - \overline{Vitamin_{d}})^2 \right).\left(\sum\limits^{}_{k_{1}} (Vitamin_{s}(k_{1}) - \overline{Vitamin_{s}})^2 \right)}} \end{align} \]
The \(DC\) and \(Vitamin\) cross-correlation are computed with these equations.

\[ \begin{align} DC \phantom{t} cross\text{-}correlation(k_{1}') = \sum\limits^{}_{k_{1}} DC_{d}(k_{1}).DC_{s}(k_{1} + k_{1}') \end{align} \] \[ \begin{align} Vitamin \phantom{t} cross\text{-}correlation(k_{1}') = \sum\limits^{}_{k_{1}} Vitamin_{d}(k_{1}).Vitamin_{s}(k_{1} + k_{1}') \end{align} \] where

• \(DC_{d}(k_{1})\) : deterministic \(DC\) function of \(k_{1}\).
• \(DC_{s}(k_{1})\) : stochastic \(DC\) function of \(k_{1}\).
• \(Vitamin_{d}(k_{1})\) : deterministic \(Vitamin\) function of \(k_{1}\).
• \(Vitamin_{s}(k_{1})\) : stochastic \(Vitamin\) function of \(k_{1}\).

Conclusion

As we can see, the stochastic and the deterministic models prodvide a priori differents results. It seems logical because of the statistical nature of the problem.
However, the comparison of the two graphs give us useful and accurate informations about the vitamin optimization.
Because of small cell counts, the stochastic and discrete method has a significant influence on the observed behaviour.
In vivo, the cells population likely follow the stochastic model. We are waiting for laboratory results to confirm or disprove our model.
We are able to maximize the vitamin production by determining the best rate constant \(k_{1}\) and therefore the best doubling time \(\mu_{1}\) to set the differentiation rate.

MATLAB algorithm

In order to make this model and code accessible, understandable and editable by everyone, we have ceated a special wiki page.
Feel free to download the source code, modified it and make it accessible for everyone.

Bibliography

• Gillespie, D.T. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of Computational Physics 22, 403-434 (1976).
• Gillespie, D.T. Exact Stochastic Simulation Of Coupled Chemical-Reactions. Journal of Physical Chemistry 81, 2340-2361 (1977).
• Gillespie, D.T. A rigorous derivation of the chemical master equation. Physica A: Statistical Mechanics and its Applications 188, 404-425 (1992).
• Cao, Y., Gillespie, D.T. & Petzold, L.R. The slow-scale stochastic simulation algorithm. The Journal of Chemical Physics 122, 14116 (2005).
• Gillespie, D.T., Hellander, A. & Petzold, L.R. Perspective: Stochastic algorithms for chemical kinetics. The Journal of Chemical Physics 138, (2013).