Difference between revisions of "Team:Paris Bettencourt/Modeling"

 
(45 intermediate revisions by 4 users not shown)
Line 1: Line 1:
 
{{Paris_Bettencourt/header}}
 
{{Paris_Bettencourt/header}}
 
{{Paris_Bettencourt/menu}}
 
{{Paris_Bettencourt/menu}}
−
{{Paris_Bettencourt/notebookBanner}}
+
{{Paris_Bettencourt/banner|page_id=modeling|page_name=Modeling}}
  
 
<html>
 
<html>
−
Work in progress :)
+
 
 +
<h1>Modeling</h1>
 +
 
 +
<h2>Introduction</h2>
 +
 
 +
<div class="column-left">
 +
Based on a set of ordinary differential equations (ODE) describing the kinetics of cell differentiation, we designed a model to find the best differentiation rate.
 +
<br />
 +
First we developed a deterministic algorithm based on the ordinary differential equation.
 +
<br />
 +
Then we found out that a stochastic algorithm could be another solution to solve our problem.
 +
<br />
 +
<br />
 +
</div>
 +
 
 +
<div class="column-right">
 +
For a system involving large cell counts, the ordinary differential equations model gives an accurate representation of the behavior. But with small cell counts, the stochastic and discrete method has a significant influence on the observed behaviour.
 +
<br />
 +
This led us to write both a deterministic program based on the mass action law and a stochastic program based on Gillespie’s stochastic simulation algorithm (SSA). With these two programs we obtain an accurate analysis of vitamin production.
 +
<br />
 +
</div>
 +
 
 +
<h2>Deterministic model</h2>
 +
<h3>Parameters</h3>
 +
We conceived a simple model with the minimum number of parameters.
 +
<br />
 +
We found seven significant parameters for our model.
 +
<br />
 +
<ul>
 +
<li>\(t\) : fermentation period.</li>
 +
<li>\(k_{1}\) : rate constant of the mother cell differentiation.</li>
 +
<li>\(k_{2}\) : rate constant of the mother cell doubling time.</li>
 +
<li>\(k_{3}\) : rate constant of the differentiated cell doubling time.</li>
 +
<li>\(k_{4}\) : rate constant of the differentiated cell vitamin production.</li>
 +
<li>\(MC_{0}\) : initial number of mother cells in the medium.</li>
 +
<li>\(DC_{0}\) : initial number of differentiated cells in the medium.</li>
 +
</ul>
 +
 
 +
<h3>Kinetic equations</h3>
 +
We wrote four simple kinetic equations.
 +
<br />
 +
\[
 +
\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}
 +
\]
 +
 
 +
<h3>Formal mathematical solution</h3>
 +
<h4>Translation in ordinary differential equations</h4>
 +
We used the law of mass action to write the ordinary differential equations based on kinetic equations \((1)\), \((2)\), \((3)\) and \((4)\).
 +
<br />
 +
We found the following equations.
 +
<br />
 +
\[
 +
\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}
 +
\]
 +
 
 +
<h4>Mathematical resolution</h4>
 +
Simple ordinary differential equations' resolution methods are used to find the solutions of the previous equations \((5)\), \((6)\), and \((7)\).
 +
<br />
 +
We expressed the time evolution of the mother cell, differentiated cell and vitamin number.
 +
<br />
 +
\[
 +
\begin{align}
 +
MC(t) = MC_{0}.e^{(k_{2} - k_{1}).t}
 +
\end{align}
 +
\]
 +
 
 +
<br />
 +
 
 +
\[
 +
\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}
 +
\]
 +
 
 +
<br />
 +
 
 +
\[
 +
\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}
 +
\]
 +
 
 +
</br>
 +
<h3>Vitamin optimization</h3>
 +
Our goal is to optimize the vitamin production. We can only change three parameters : \(k_{1}\), \(MC_0\) and \(DC_0\).
 +
<br />
 +
\(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 vitamins are produced.
 +
<br />
 +
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}\).
 +
<br />
 +
We try to find the best \(k_{1}\). To this end we maximize the vitamin function numerically with MATLAB.
 +
<br />
 +
As an example, let's consider the following parameters.
 +
<br />
 +
<ul>
 +
<li>\(t = 10 \phantom{t} (hours)\)</li>
 +
<li>\(k_{2} = 0.66 \phantom{t} (hours^{-1})\)</li>
 +
<li>\(k_{3} = 0.33 \phantom{t} (hours^{-1})\)</li>
 +
<li>\(k_{4} = 1 \phantom{t} (hours^{-1})\)</li>
 +
<li>\(MC_0 = 5 \phantom{t} (cells)\)</li>
 +
<li>\(DC_0 = 0 \phantom{t} (cells)\)</li>
 +
</ul>
 +
We obtain the following graph with a simple MATLAB program available on the software section.
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/4/45/OptimizeK1.png" title="Deterministic vitamin optimization" alt="Deterministic vitamin optimization" style="align:center;">
 +
<p style="text-align:center"><b>Figure 1 :</b> Deterministic vitamin optimization.</p>
 +
<br />
 +
As you can see, we are able to find the best \(k_{1}\) to optimize the vitamin production.
 +
<br />
 +
We noticed that two different \(k_{1}\) exist to optimize the maximum number of differentiated cell \(DC\) or the maximum number of \(Vitamin\).
 +
<br />
 +
In our case, we chose of course the \(k_{1}\) that optimizes the vitamin production <i>i.e.</i> \(k_{1} = 0.207 \phantom{t} (hours^{-1})\).
 +
<h3>Deterministic evolution of mother cell, differentiated cell and vitamin numbers</h3>
 +
We wrote a deterministic algorithm with MATLAB using the previous solutions \((8)\), \((9)\) and \((10)\). For those interested, the source code is available on the software section.
 +
<br />
 +
In order to optimize the vitamin production, we use the same parameters as previously and set \(k_{1}\) with the previous resulting value <i>i.e.</i> \(k_{1} =
 +
0.207 \phantom{t} (hours^{-1})\).
 +
<br />
 +
We obtain the following graph.
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/f/f7/DeterministicEvolution.png" title="Deterministic evolution of the system" alt="Deterministic evolution of the
 +
system" style="align:center;">
 +
<p style="text-align:center"><b>Figure 2 :</b> Deterministic evolution of the system.</p>
 +
<br />
 +
As you can see, the mother cell, the differentiated cell and the vitamin numbers follow an exponential law of time.
 +
<br />
 +
This result seems relevant. The model does not take into account the cells' death and the nutrients present in the medium.
 +
</d>
 +
<h2>Stochastic model</h2>
 +
 
 +
<h3>Parameters</h3>
 +
We used more parameters than the deterministic model. Indeed, because of the nature of the algorithm, we need a constant rate probability distribution.
 +
<br />
 +
In order to do a simulation with laboratory data, we modelize the time probability distribution with a normalized Gaussian distribution using the <a href="http://fr.mathworks.com/help/stats/normrnd.html">normrnd</a> MATLAB function.
 +
<br />
 +
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
 +
<ul>
 +
<li>\(\mathcal P_{i}(t)\) : time \(t\) probability.</li>
 +
<li>\(\mu_{i}\) : time mean.</li>
 +
<li>\(\sigma_{i}\) : time standard deviation.</li>
 +
</ul>
 +
 
 +
<br />
 +
 
 +
Rate constants for the time distribution \(i\) are calculated using the following 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.
 +
 
 +
<br />
 +
The vitamin number is computed by numerical integration using the rate constant \(k_{4}\).
 +
<br />
 +
We have ten parameters for this model.
 +
<br />
 +
<ul>
 +
<li>\(t\) : fermentation period.</li>
 +
<li>\(\mu_{1}\) : mother cell differentiation time mean.</li>
 +
<li>\(\sigma_{1}\) : mother cell differentiation time standard deviation.</li>
 +
<li>\(\mu_{2}\) : mother cell doubling time mean.</li>
 +
<li>\(\sigma_{2}\) : mother cell doubling time standard deviation.</li>
 +
<li>\(\mu_{3}\) : differentiated cell doubling time mean.</li>
 +
<li>\(\sigma_{3}\) : differentiated cell doubling time standard deviation.</li>
 +
<li>\(k_{4}\) : rate constant of the differentiated cell vitamin production.</li>
 +
<li>\(MC_0\) : initial number of mother cells in the medium.</li>
 +
<li>\(DC_0\) : initial number of differentiated cells in the medium.</li>
 +
</ul>
 +
 
 +
<h3>Algorithm</h3>
 +
 
 +
<h4>Introduction</h4>
 +
 
 +
In the program, each cell is a structure containing three elements.
 +
<ul>
 +
<li>The cell type.
 +
<ul>
 +
<li>1 : mother cell.</li>
 +
<li>2 : differentiated cell.</li>
 +
</ul>
 +
</li>
 +
<li>The next event type.
 +
<ul>
 +
<li>1 : a mother cell becomes a differentiated cell.</li>
 +
<li>2 : a mother cell divides into two mother cells.</li>
 +
<li>3 : a differentiated cell divides into two differentiated cells.</li>
 +
 
 +
</ul>
 +
</li>
 +
<li>The time before the next event.</li>
 +
 
 +
</ul>
 +
 
 +
<br />
 +
The time \(T_{i}\) for an event is chosen thanks to the equation \((11)\).
 +
<br />
 +
There is a competition between the first two events concerning the mother cell.
 +
<br />
 +
Either the mother cell differentiates or it divides into two differentiated cells.
 +
<br />
 +
To do the choice we compare \(T_{1}\) and \(T_{2}\). The smallest time is the next event time.
 +
<br />
 +
<img src="https://static.igem.org/mediawiki/2015/7/72/GeneralStochasticAlgorithm.png"  height="420" style="float:right;" alt="Algorithm flowchart" title="Algorithm flowchart" >
 +
 
 +
<h4>Step 1 : Initialize the cells.</h4>
 +
First of all we have to generate the first cells in accordance with \(MC_0\) and \(DC_0\).
 +
<br />
 +
For each cell we set the cell type, the next event and the time before the next event.
 +
<br />
 +
<h4>Step 2 : Choose the next event among all cells.</h4>
 +
We have an array containing all the cells. Each cell contains the time before its next event (divide or differentiate).
 +
<br />
 +
The next event time processed by the algorithm is the smallest cell next event time.
 +
<br />
 +
<h4>Step 3 : Do the next event.</h4>
 +
Now that we know the next event type and the concerned cell (only one cell is involved here), we compute this event.
 +
<br />
 +
As explained above, we have three choices.
 +
<br />
 +
If the cell is a mother cell, the cell can divide or differentiate.
 +
<br />
 +
If the cell is a differentiated cell, the cell can only divide.
 +
<br />
 +
<br />
 +
We define a new event and time event for the cell and for the new cell in case of division.
 +
<br  />
 +
We also performed the number of vitamin calculus by numerically integrating the number of \(DC\).
 +
<br />
 +
We use the following formula.
 +
 
 +
\[
 +
\begin{align}
 +
Vitamin(t) = Vitamin(t-T_{i}) + k_{4}.DC(t).T_{i}
 +
\end{align}
 +
\]
 +
 
 +
where
 +
<ul>
 +
<li>\(Vitamin(t)\) : number of vitamin at the instant \(t\).</li>
 +
<li>\(DC(t)\) : number of differentiated cells in the medium at the instant \(t\).</li>
 +
<li>\(T_{i}\) : time of the next event being performed.</li>
 +
<li>\(k_{4}\) : rate constant of the differentiated cell vitamin production.</li>
 +
</ul>
 +
 
 +
<h4>Step 4 : Update the simulation time.</h4>
 +
For each cell except the concerned cell, define the cell next event time with the following formula.
 +
 
 +
\[
 +
\begin{align}
 +
T_{event \phantom{t} cell} = T_{event \phantom{t} cell} - T_{i}
 +
\end{align}
 +
\]
 +
 
 +
where
 +
<ul>
 +
<li>\(T_{i}\) : time of the next event being performed.</li>
 +
<li>\(T_{event \phantom{t} cell}\) : cell next time event.</li>
 +
</ul>
 +
 
 +
<h4>Step 5 : Go to step 2 until the simulation time is superior to the fermentation period.</h4>
 +
We do this process until the simulation time is superior to the fermentation period.
 +
The number of cells increases with an exponential law. The computational power increases with the same law.
 +
<br />
 +
Do not try to process a long simulation time if your computer is not powerful.
 +
<h4>Step 6 : Show the results.</h4>
 +
Here is a graph generated by the stochastic algorithm.
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/1/14/StochasticEvolution.png" title="Stochastic evolution of the system" alt="Stochastic evolution of the system" style="align:center;">
 +
<p style="text-align:center"><b>Figure 3 :</b> Stochastic evolution of the system.</p>
 +
 
 +
With the following distribution of event time probability.
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/9/97/StochasticEventTimeProbabilityDistribution.png" title="Distribution of event time probability" alt="Distribution of event time probability" style="align:center;">
 +
<p style="text-align:center"><b>Figure 4 :</b> Distribution of event time probability.</p>
 +
 
 +
<h3>Vitamin optimization</h3>
 +
Like the deterministic model our goal is to maximize the vitamin production.
 +
<br />
 +
In order to get relevant results, we average by doing a large number of simulations specified by the \(Averaging\) \(number\) parameter.
 +
<br />
 +
We can only change three parameters : \(\mu_{1}\), \(MC_0\) and \(DC_0\).
 +
<br />
 +
\(t\), \(\sigma_{1}\), \(\mu_{2}\), \(\sigma_{2}\),  \(\mu_{3}\), \(\sigma_{3}\) and \(k_{4}\) are constants.
 +
<br />
 +
We try to find the best \(\mu_{1}\) to maximize the vitamin production.
 +
<br />
 +
As an example we can consider the following parameters.
 +
<br />
 +
<ul>
 +
<li>\(t = 5 \phantom{t} (hours)\)</li>
 +
<li>\(Averaging\) \(number = 500\)</li>
 +
<li>\(\sigma_{1} = 0.1 \phantom{t} (hours)\)</li>
 +
<li>\(\mu_{2} = 1.0502 \phantom{t} (hours)\)</li>
 +
<li>\(\sigma_{2} = 0.1 \phantom{t} (hours)\)</li>
 +
<li>\(\mu_{3} = 2.1004 \phantom{t} (hours)\)</li>
 +
<li>\(\sigma_{3} = 0.1 \phantom{t} (hours)\)</li>
 +
<li>\(k_{4} = 1 \phantom{t} (hours^{-1})\)</li>
 +
<li>\(MC_0 = 5 \phantom{t} (cells)\)</li>
 +
<li>\(DC_0 = 0 \phantom{t} (cells)\)</li>
 +
</ul>
 +
<br />
 +
We obtain the following graph.
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/0/0b/OptimizeK1Stochastic.png" title="Stochastic vitamin optimization" alt="Stochastic vitamin optimization" style="align:center;">
 +
<p style="text-align:center"><b>Figure 5 :</b> Stochastic vitamin optimization.</p>
 +
<br />
 +
Here is a deterministic vitamin optimization.
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/c/ca/OptimizeK1StochasticDeterministic.png" title="Deterministic vitamin optimization" alt="Deterministic vitamin optimization" style="align:center;">
 +
<p style="text-align:center"><b>Figure 6 :</b> Deterministic vitamin optimization.</p>
 +
<br />
 +
We can compare the two results by superimposing the two graphs.
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/3/3c/CompareStochasticAndDeterministic.png" title="Comparison of deterministic and stochastic vitamin optimization" alt="Comparison of deterministic and stochastic vitamin optimization" style="align:center;">
 +
<p style="text-align:center"><b>Figure 7 :</b> Comparison of deterministic and stochastic vitamin optimization.</p>
 +
<br />
 +
We performed a cross-correlation and found correlation coefficients between the two graphs with the <a href="http://fr.mathworks.com/help/signal/ref/xcorr.html">xcorr</a> and <a href="http://fr.mathworks.com/help/images/ref/corr2.html">corr2</a> MATLAB functions.
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/7/76/CrossCorrelation.png" title="Deterministic and stochastic cross-correlation and correlation coefficients" alt="Deterministic and stochastic cross-correlation and correlation coefficients" style="align:center;">
 +
<p style="text-align:center"><b>Figure 8 :</b> Deterministic and stochastic cross-correlation and correlation coefficients.</p>
 +
<br />
 +
As described in the <a href="http://fr.mathworks.com/help/images/ref/corr2.html">official MATLAB website</a>, the \(DC\) and \(Vitamin\) correlation coefficients are computed with the following formula.
 +
<br />
 +
<br />
 +
 
 +
\[
 +
\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}
 +
\]
 +
 
 +
<br />
 +
 
 +
\[
 +
\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}
 +
\]
 +
 
 +
<br />
 +
The \(DC\) and \(Vitamin\) cross-correlation are computed with these equations.
 +
<br />
 +
<br />
 +
 
 +
\[
 +
\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
 +
<br />
 +
<ul>
 +
<li>\(DC_{d}(k_{1})\) : deterministic \(DC\) function of \(k_{1}\).</li>
 +
<li>\(DC_{s}(k_{1})\) : stochastic \(DC\) function of \(k_{1}\).</li>
 +
<li>\(Vitamin_{d}(k_{1})\) : deterministic \(Vitamin\) function of \(k_{1}\).</li>
 +
<li>\(Vitamin_{s}(k_{1})\) : stochastic \(Vitamin\) function of \(k_{1}\).</li>
 +
</ul>
 +
 
 +
<h2>Conclusion</h2>
 +
 
 +
As we can see, the stochastic and the deterministic models provide <i>a priori</i> different results.
 +
It seems logical because of the statistical nature of the problem.
 +
<br />
 +
However, the comparison of the two graphs gives us useful and accurate information about the vitamin optimization.
 +
<br />
 +
Because of small cell counts, the stochastic and discrete method has a significant influence on the observed behaviour.
 +
<br />
 +
<i>In vivo</i>, the cell population likely follows the stochastic model.
 +
We are waiting for laboratory results to confirm or disprove our model.
 +
<br />
 +
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.
 +
<br />
 +
<br />
 +
<br />
 +
 
 +
<h1>Software</h1>
 +
 
 +
<h2>Introduction</h2>
 +
As explained in the modeling section, this program models the cell population evolution with mother cell differentiation and cell division.
 +
<br />
 +
In order to make this model and code accessible, understandable and editable by everyone, we have created this software section.
 +
<br />
 +
We focus here on the stochastic program.
 +
 
 +
<h2>Files</h2>
 +
The MATLAB code is available on <a href="https://github.com/iGEMParisBettencourt2015/modeling" title="Source code in GitHub" >GitHub</a> under the GNU General Public License version 3.0.
 +
<br />
 +
Feel free to download the source code, modify it and make it accessible for everyone.
 +
 
 +
<br />
 +
Here is a quick explanation of the different files <i>i.e.</i> the different functions.
 +
<h3>Stochastic algorithm</h3>
 +
<ul>
 +
<li>\(timeEvolutionStochastic\) : calculate the time evolution.</li>
 +
<li>\(checkInput\) : check the input and show an error message if the input is not good.</li>
 +
<li>\(elmPixel\) : find the coordinates of a pixel around another pixel.</li>
 +
<li>\(errorProgramm\) : display an error and stop the program.</li>
 +
<li>\(freePosition\) : create a list of free pixels around one pixel.</li>
 +
<li>\(getComputerName\) : get the computer name.</li>
 +
<li>\(getSizeMat\) : find the matrix size.</li>
 +
<li>\(imageToMatrix\) : convert an image into a matrix.</li>
 +
<li>\(initCells\) : initialize the cells.</li>
 +
<li>\(initialization\) : initialize the program.</li>
 +
<li>\(initTimeEvent\) : calculate the different time.</li>
 +
<li>\(initVariables\) : initialize the variables.</li>
 +
<li>\(isBinary\) : check if a number is 0 or 1.</li>
 +
<li>\(isPositiveInteger\) : check if a number is a positive integer (1) or not (0).</li>
 +
<li>\(makeCells\) : make the cells.</li>
 +
<li>\(randEvent\) : generate the next random event.</li>
 +
<li>\(randPosition\) : find a free position for a cell.</li>
 +
<li>\(saveInitialParameters\) : save the initial parameters in a text file.</li>
 +
<li>\(saveResult\) : save and plot the results.</li>
 +
<li>\(showAndSaveEventTimeDistribution\) : show and save the event time probability distribution.</li>
 +
<li>\(showPopulation\) : show the cells matrix.</li>
 +
<li>\(writeParameters\) : write a text file with the parameters.</li>
 +
<li>\(launchStochasticProgramm\) : launch the stochastic time evolution simulation.</li>
 +
<li>\(k1OptimizationStochastic\) : calculate the number of differentiated cells and vitamin with different \(k_{1}\) in a stochastic way.</li>
 +
<li>\(k1OptimizationStochasticAndDeterministic\) : superimpose the deterministic and stochastic models.</li>
 +
</ul>
 +
<h3>Deterministic algorithm</h3>
 +
<ul>
 +
<li>\(timeEvolutionDeterministic\) : show the time evolution by numerical calculus and equations.</li>
 +
<li>\(k1OptimizationDeterministic\) : calculate the number of differentiated cells with different \(k_{1}\) in a deterministic way.</li>
 +
</ul>
 +
<h2>Input</h2>
 +
In order to launch a simulation, you must define five variables.
 +
<ul>
 +
<li>Initial cells.</li>
 +
<li>Rate constants.</li>
 +
<li>Fermentation period.</li>
 +
<li>Action.</li>
 +
<li>Folder name.</li>
 +
</ul>
 +
<h3>Cells input</h3>
 +
You can choose to generate a random or a predetermined position for the initial cells.
 +
<br />
 +
You have two cell types : mother cells and differentiated cells.
 +
<h4>Random initial position</h4>
 +
The input array must contain four integer values in this order.
 +
<ul>
 +
<li>\(MC_0\) : initial number of mother cells in the medium.</li>
 +
<li>\(DC_0\) : initial number of differentiated cells in the medium.</li>
 +
<li>\(sizeMatCell\) : cells box size (in pixels).</li>
 +
<li>\(sizeMat\) : real box size (in pixels).</li>
 +
</ul>
 +
<br />
 +
Warning : \(sizeMatCell\) must be larger than \(sizeMat\).
 +
<br />
 +
<br />
 +
 
 +
In the program, the input variable has the following pattern.
 +
<br />
 +
\(input = [MC_0, DC_0, sizeMatCell, sizeMat]\)
 +
<br />
 +
<br />
 +
Here is an example.
 +
<br />
 +
<ul>
 +
<li>\(MC_0 = 5 \phantom{t} (cells)\)</li>
 +
<li>\(DC_0 = 0 \phantom{t} (cells)\)</li>
 +
<li>\(sizeMatCell = 100 \phantom{t} (pixels)\)</li>
 +
<li>\(sizeMat = 250 \phantom{t} (pixels)\)</li>
 +
</ul>
 +
 
 +
<h4>Predetermined initial position</h4>
 +
In order to make an easy predetermined initial position input, we wrote a code that understands images.
 +
<br />
 +
The input must contain the image path.
 +
<br />
 +
It is very easy to design a custom image <i>i.e.</i> a custom spatial cell distribution. You must follow these four rules (we work only with the image red component)
 +
<ul>
 +
<li>the image must be square and in jpeg format.</li>
 +
<li>255 means no cell.</li>
 +
<li>128 means one differentiated cell.</li>
 +
<li>0 means one mother cell.</li>
 +
</ul>
 +
0, 128 and 255 are red components of any colours that you want because the MATLAB code take into account only the image red component.
 +
<br />
 +
If you want, you can use the following hexadecimal colors codes.
 +
<ul>
 +
<li>#FFFFFF (white) means no cell.</li>
 +
<li>#800000 (red) means one differentiated cell.</li>
 +
<li>#000000 (black) means one mother cell.</li>
 +
</ul>
 +
<br />
 +
Here is an example only with mother cells.
 +
<br />
 +
<br />
 +
<img src="https://static.igem.org/mediawiki/2015/9/9a/IgemImg.png" style="width:40%;height:40%;" alt="Example of predetermined cells initial position with mother cells" title="Example of predetermined cells initial position with mother cells" >
 +
<p><b>Figure 9 :</b> Example of predetermined cells initial position with mother cells.</p>
 +
<br />
 +
Here is another example with mother cells (black) and differentiated cells (red).
 +
<br />
 +
<br />
 +
<img src="https://static.igem.org/mediawiki/2015/a/a4/IgemParisBettencourtImg.png" style="width:40%;height:40%;" alt="Example of predetermined cells initial position with mother cells (black) and differentiated cells (red)" title="Example of predetermined cells initial position with mother cells (black) and differentiated cells (red)" >
 +
<p><b>Figure 10 :</b> Example of predetermined cells initial position with mother cells (black) and differentiated cells (red).</p>
 +
<br />
 +
 
 +
<h3>Rate constants</h3>
 +
You have to define seven variables.
 +
<ul>
 +
<li>\(\mu_{1}\) : mother cell differentiation time mean.</li>
 +
<li>\(\sigma_{1}\) : mother cell differentiation time standard deviation.</li>
 +
<li>\(\mu_{2}\) : mother cell doubling time mean.</li>
 +
<li>\(\sigma_{2}\) : mother cell doubling time standard deviation.</li>
 +
<li>\(\mu_{3}\) : differentiated cell doubling time mean.</li>
 +
<li>\(\sigma_{3}\) : differentiated cell doubling time standard deviation.</li>
 +
<li>\(k_{4}\) : rate constant of the differentiated cell vitamin production.</li>
 +
</ul>
 +
<br />
 +
In the program, the rate constants variable has the following pattern.
 +
<br />
 +
\(constantRate = [\mu_{1}, \sigma_{1}, \mu_{2}, \sigma_{2}, \mu_{3}, \sigma_{3}, k_{4}]\)
 +
<br />
 +
<br />
 +
Here is an example.
 +
<br />
 +
<ul>
 +
<li>\(\mu_{1} = 1 \phantom{t} (hours)\)</li>
 +
<li>\(\sigma_{1} = 0.1 \phantom{t} (hours)\)</li>
 +
<li>\(\mu_{2} = 1.0502 \phantom{t} (hours)\)</li>
 +
<li>\(\sigma_{2} = 0.1 \phantom{t} (hours)\)</li>
 +
<li>\(\mu_{3} = 2.1004 \phantom{t} (hours)\)</li>
 +
<li>\(\sigma_{3} = 0.1 \phantom{t} (hours)\)</li>
 +
<li>\(k_{4} = 1 \phantom{t} (hours^{-1})\)</li>
 +
</ul>
 +
 
 +
<h3>Fermentation period</h3>
 +
The fermentation period \(t\) is a simple scalar variable.
 +
<br />
 +
You can use \(t = 5 \phantom{t} (hours)\) for example.
 +
<h3>Action</h3>
 +
In order to choose what you want to do, you must define three Boolean variables.
 +
<ul>
 +
<li>\(createFolder\) : create a folder in the disk to save information.</li>
 +
<li>\(showAnimation\) : show the cell population animation. Warning : requires a lot of computing power.</li>
 +
<li>\(plotGraph\) : plot the results in graphs.</li>
 +
</ul>
 +
 
 +
<br />
 +
In the program, the action variable has the following pattern.
 +
<br />
 +
\(action = [createFolder, showAnimation, plotGraph\)]
 +
<br />
 +
<br />
 +
In the following example we want to create a folder and plot the graphs.
 +
<br />
 +
<ul>
 +
<li>\(createFolder = 1\)</li>
 +
<li>\(showAnimation = 0\)</li>
 +
<li>\(plotGraph = 1\)</li>
 +
</ul>
 +
<br />
 +
In the following example we want to show the cell animation.
 +
<br />
 +
<ul>
 +
<li>\(createFolder = 0\)</li>
 +
<li>\(showAnimation = 1\)</li>
 +
<li>\(plotGraph = 0\)</li>
 +
</ul>
 +
<br />
 +
In the following example we want to show the cell animation and save the images in a folder.
 +
<br />
 +
With the images you can create gif animations.
 +
<br />
 +
<ul>
 +
<li>\(createFolder = 1\)</li>
 +
<li>\(showAnimation = 1\)</li>
 +
<li>\(plotGraph = 0\)</li>
 +
</ul>
 +
 
 +
<h3>Folder name</h3>
 +
The \(folderName\) variable is a string containing the folder name where results folder are stored in the disk.
 +
<br />
 +
You can use \(folderName =\) '\(Results\)' for example.
 +
 
 +
<h2>Output</h2>
 +
Four types of output are generated.
 +
<ul>
 +
<li>raw data : all the useful information.</li>
 +
<li>folder creation : if chosen in the \(action\) variable a folder is created.</li>
 +
<li>graphs : if chosen in the \(action\) variable graphs are generated.</li>
 +
<li>animation : if chosen in the \(action\) variable cell animation is generated.</li>
 +
</ul>
 +
 
 +
<h3>Raw data</h3>
 +
<ul>
 +
<li>\(mainTic\) : start of the program.</li>
 +
<li>\(MC_0\) : initial number of mother cells in the medium.</li>
 +
<li>\(DC_0\) : initial number of differentiated cells in the medium.</li>
 +
<li>\(pasDistributionTimeEvent\) : resolution for the graph time event distribution plot (see the \(initVariables\) function).</li>
 +
<li>\(frameRate\) : animation frame rate (see the \(initVariables\) function).</li>
 +
<li>\(timeShowPopulation\) : simulation time between two animations (see the \(initVariables\) function).</li>
 +
<li>\(simulationTime\) : fermentation period \(t\).</li>
 +
<li>\(resultFolderName\) : folder where the results folder are created.</li>
 +
<li>\(\mu_{1}\) : mother cell differentiation time mean.</li>
 +
<li>\(\sigma_{1}\) : mother cell differentiation time standard deviation.</li>
 +
<li>\(\mu_{2}\) : mother cell doubling time mean.</li>
 +
<li>\(\sigma_{2}\) : mother cell doubling time standard deviation.</li>
 +
<li>\(\mu_{3}\) : differentiated cell doubling time mean.</li>
 +
<li>\(\sigma_{3}\) : differentiated cell doubling time standard deviation.</li>
 +
<li>\(k_{4}\) : rate constant of the differentiated cell vitamin production.</li>
 +
<li>\(cVitamin\) : vitamin counter.</li>
 +
<li>\(time\) : main timer containing the simulation time.</li>
 +
<li>\(createFolder\) : action boolean.</li>
 +
<li>\(showPopulationBool\) : action boolean.</li>
 +
<li>\(plotGraph\) : action boolean.</li>
 +
<li>\(input\) : input array.</li>
 +
<li>\(cTime1\) : time distribution counter 1 (see the \(randEvent\) function).</li>
 +
<li>\(cTime2\) : time distribution counter 2 (see the \(randEvent\) function).</li>
 +
<li>\(cTime3\) : time distribution counter 3 (see the \(randEvent\) function).</li>
 +
<li>\(counterEvent\) : event counter.</li>
 +
<li>\(timerShow\) : time animation counter.</li>
 +
<li>\(sizeMat\) : real box size (in pixels).</li>
 +
<li>\(nextEvent\) : next event type.</li>
 +
<li>\(time1Array\) : time distribution array 1 (see the \(initTimeEvent\) function).</li>
 +
<li>\(time2Array\) : time distribution array 2 (see the \(initTimeEvent\) function).</li>
 +
<li>\(time3Array\) : time distribution array 3 (see the \(initTimeEvent\) function).</li>
 +
<li>\(cell\) : structure containing all the cells' information.</li>
 +
<li>\(sizeMatCell\) : cells box size (in pixels).</li>
 +
<li>\(timeNextEvent\) : time before the next event.</li>
 +
<li>\(saveTime\) : array containing the simulation time evolution.</li>
 +
<li>\(saveNbrMC\) : array containing the mother cells time evolution.</li>
 +
<li>\(saveNbrDC\) : array containing the differentiated cells time evolution.</li>
 +
<li>\(saveNbrVitamin\) : array containing the vitamin time evolution.</li>
 +
<li>\(folderName\) : folder where the results are stored. The pattern is the following : \(resultFolderName/DD-MM-YY\_hh\_mm\_ss\_mss\).</li>
 +
<li>\(imgType\) : image type (see the \(initialization\) function).</li>
 +
<li>\(parametersFileName\) : parameters file name where the parameters are saved (see the \(initialization\) function).</li>
 +
<li>\(parametersDoc\) : file object where the parameters are saved.</li>
 +
<li>\(maxMC\) : maximum number of mother cells.</li>
 +
<li>\(maxDC\) : maximum number of differentiated cells.</li>
 +
<li>\(maxVitamin\) : maximum number of vitamins.</li>
 +
<li>\(timeSpent\) : computation time (in seconds).</li>
 +
 
 +
</ul>
 +
<h3>Folder creation</h3>
 +
If the \(action\) variable \(createFolder\) is set, the program will save important results in this folder.
 +
<br />
 +
<ul>
 +
<li>\(parameters.txt\) : text file containing different simulation parameters and results.</li>
 +
<li>\(data.mat\) : MATLAB file containing the raw data structure. To load the file, use this formula : \(load('data.mat')\).</li>
 +
</ul>
 +
 
 +
<h3>Graphs</h3>
 +
If the \(action\) variable \(plotGraph\) is set, the program will show the graphs.
 +
<br />
 +
Moreover, if the \(action\) variable \(createFolder\) is set, the program will save the graphs in the folder.
 +
 
 +
<h3>Animation</h3>
 +
If the \(action\) variable \(showAnimation\) is set, the program will show the cell animation.
 +
<br />
 +
Moreover, if the \(action\) variable \(createFolder\) is set, the program will save the images in the folder.
 +
 
 +
<h2>Examples</h2>
 +
We present here some results obtained with the program. For more details about the model sees the modeling section.
 +
 
 +
<h3>Time evolution of mother cells, differentiated cells and vitamin</h3>
 +
<h4>Graphs</h4>
 +
As explained in the modeling section the program give these results.
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/1/14/StochasticEvolution.png" title="Stochastic evolution of the system" alt="Stochastic evolution of the system" style="align:center;">
 +
<p style="text-align:center"><b>Figure 11 :</b> Stochastic evolution of the system.</p>
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/0/0b/OptimizeK1Stochastic.png" title="Stochastic vitamin optimization" alt="Stochastic vitamin optimization" style="align:center;">
 +
<p style="text-align:center"><b>Figure 12 :</b> Stochastic vitamin optimization.</p>
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/9/97/StochasticEventTimeProbabilityDistribution.png" title="Distribution of event time probability" alt="Distribution of event time probability" style="align:center;">
 +
<p style="text-align:center"><b>Figure 13 :</b> Distribution of event time probability.</p>
 +
 
 +
<h4>Animation</h4>
 +
Here is an example with mother cells (black) and differentiated cells (red).
 +
<br />
 +
This image is the cell input.
 +
<br />
 +
<br />
 +
<img style="width:40%;height:40%;" src="https://static.igem.org/mediawiki/2015/a/a4/IgemParisBettencourtImg.png" alt="Example of predetermined cells initial position with mother cells (black) and differentiated cells (red)" title="Example of predetermined cells initial position with mother cells (black) and differentiated cells (red)" >
 +
<p><b>Figure 14 :</b> Example of predetermined cells initial position with mother cells (black) and differentiated cells (red).</p>
 +
<br />
 +
The result is the following animation after a conversion in a gif file.
 +
<br />
 +
The mother cells are in orange and the differentiated cells are in yellow.
 +
<br />
 +
<br />
 +
<img style="width:40%;height:40%;" src="https://static.igem.org/mediawiki/2015/7/73/IgemParisBettencourtGif.gif" alt="Animation example with predetermined cells initial position with mother cells and differentiated cells" title="Animation example with predetermined cells initial position with mother cells and differentiated cells" >
 +
<p><b>Figure 15 :</b> Animation example with predetermined cells initial position with mother cells and differentiated cells.</p>
 +
<br />
 +
 
 +
<h3>Vitamin optimization</h3>
 +
This program is designed to maximize the vitamin production. Here are some results.
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/0/0b/OptimizeK1Stochastic.png" title="Stochastic vitamin optimization" alt="Stochastic vitamin optimization" style="align:center;">
 +
<p style="text-align:center"><b>Figure 16 :</b> Stochastic vitamin optimization.</p>
 +
 
 +
<h3>Deterministic and stochastic vitamin optimization comparison</h3>
 +
Here is a comparison between the two models.
 +
<br />
 +
<br />
 +
<img width="100%" src="https://static.igem.org/mediawiki/2015/3/3c/CompareStochasticAndDeterministic.png" title="Comparison of deterministic and stochastic vitamin optimization" alt="Comparison of deterministic and stochastic vitamin optimization" style="align:center;">
 +
<p style="text-align:center"><b>Figure 17 :</b> Comparison of deterministic and stochastic vitamin optimization.</p>
 +
 
 +
<h2>Conclusion</h2>
 +
Feel free to use and modify this program. The source code is available on <a href="https://github.com/iGEMParisBettencourt2015/modeling" title="Source code in GitHub"> GitHub : iGEM Paris Bettencourt 2015</a> under the <a href="http://www.gnu.org/licenses/gpl.txt" alt="GNU General Public License" >GNU General Public License version 3.0</a>.
 +
For more information concerning the model see the modeling section.
 +
 
 +
 
 +
 
 
</html>
 
</html>
  
 
{{Paris_Bettencourt/footer}}
 
{{Paris_Bettencourt/footer}}

Latest revision as of 18:44, 28 October 2015

Modeling

Introduction

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

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

Deterministic model

Parameters

We conceived a simple model with the minimum number of parameters.
We found 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 vitamins are produced.
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, let's 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 on the software section.

Deterministic vitamin optimization

Figure 1 : Deterministic vitamin optimization.


As you can see, we are able to find the best \(k_{1}\) to optimize the vitamin production.
We noticed that two different \(k_{1}\) exist 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 optimizes the vitamin production i.e. \(k_{1} = 0.207 \phantom{t} (hours^{-1})\).

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 on the software section.
In order to optimize the vitamin production, we use the same parameters as previously and set \(k_{1}\) with the previous resulting value i.e. \(k_{1} = 0.207 \phantom{t} (hours^{-1})\).
We obtain the following graph.

Deterministic evolution of the system

Figure 2 : Deterministic evolution of the system.


As you can see, the mother cell, the differentiated cell and the 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\) are calculated using the following 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 becomes a differentiated cell.
    • 2 : a mother cell divides into two mother cells.
    • 3 : a differentiated cell divides 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 differentiates or it divides into two differentiated cells.
To do the choice we compare \(T_{1}\) and \(T_{2}\). The smallest time is the next event time.
Algorithm flowchart

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 contains the time before its next event (divide 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 divide or differentiate.
If the cell is a differentiated cell, the cell can only divide.

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 \phantom{t} cell} = T_{event \phantom{t} cell} - T_{i} \end{align} \] where
  • \(T_{i}\) : time of the next event being performed.
  • \(T_{event \phantom{t} cell}\) : cell next time event.

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

We do this process until the simulation time is superior to the fermentation period. The number of cells increases with an exponential law. The computational power increases with the same law.
Do not try to process a long simulation time if your computer is not powerful.

Step 6 : Show the results.

Here is a graph generated by the stochastic algorithm.

Stochastic evolution of the system

Figure 3 : Stochastic evolution of the system.

With the following distribution of event time probability.

Distribution of event time probability

Figure 4 : 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 obtain the following graph.

Stochastic vitamin optimization

Figure 5 : Stochastic vitamin optimization.


Here is a deterministic vitamin optimization.

Deterministic vitamin optimization

Figure 6 : Deterministic vitamin optimization.


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

Comparison of deterministic and stochastic vitamin optimization

Figure 7 : Comparison of deterministic and stochastic vitamin optimization.


We performed a cross-correlation and found correlation coefficients between the two graphs with the xcorr and corr2 MATLAB functions.

Deterministic and stochastic cross-correlation and correlation coefficients

Figure 8 : Deterministic and stochastic cross-correlation and correlation coefficients.


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 provide a priori different results. It seems logical because of the statistical nature of the problem.
However, the comparison of the two graphs gives us useful and accurate information 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 cell population likely follows 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.


Software

Introduction

As explained in the modeling section, this program models the cell population evolution with mother cell differentiation and cell division.
In order to make this model and code accessible, understandable and editable by everyone, we have created this software section.
We focus here on the stochastic program.

Files

The MATLAB code is available on GitHub under the GNU General Public License version 3.0.
Feel free to download the source code, modify it and make it accessible for everyone.
Here is a quick explanation of the different files i.e. the different functions.

Stochastic algorithm

  • \(timeEvolutionStochastic\) : calculate the time evolution.
  • \(checkInput\) : check the input and show an error message if the input is not good.
  • \(elmPixel\) : find the coordinates of a pixel around another pixel.
  • \(errorProgramm\) : display an error and stop the program.
  • \(freePosition\) : create a list of free pixels around one pixel.
  • \(getComputerName\) : get the computer name.
  • \(getSizeMat\) : find the matrix size.
  • \(imageToMatrix\) : convert an image into a matrix.
  • \(initCells\) : initialize the cells.
  • \(initialization\) : initialize the program.
  • \(initTimeEvent\) : calculate the different time.
  • \(initVariables\) : initialize the variables.
  • \(isBinary\) : check if a number is 0 or 1.
  • \(isPositiveInteger\) : check if a number is a positive integer (1) or not (0).
  • \(makeCells\) : make the cells.
  • \(randEvent\) : generate the next random event.
  • \(randPosition\) : find a free position for a cell.
  • \(saveInitialParameters\) : save the initial parameters in a text file.
  • \(saveResult\) : save and plot the results.
  • \(showAndSaveEventTimeDistribution\) : show and save the event time probability distribution.
  • \(showPopulation\) : show the cells matrix.
  • \(writeParameters\) : write a text file with the parameters.
  • \(launchStochasticProgramm\) : launch the stochastic time evolution simulation.
  • \(k1OptimizationStochastic\) : calculate the number of differentiated cells and vitamin with different \(k_{1}\) in a stochastic way.
  • \(k1OptimizationStochasticAndDeterministic\) : superimpose the deterministic and stochastic models.

Deterministic algorithm

  • \(timeEvolutionDeterministic\) : show the time evolution by numerical calculus and equations.
  • \(k1OptimizationDeterministic\) : calculate the number of differentiated cells with different \(k_{1}\) in a deterministic way.

Input

In order to launch a simulation, you must define five variables.
  • Initial cells.
  • Rate constants.
  • Fermentation period.
  • Action.
  • Folder name.

Cells input

You can choose to generate a random or a predetermined position for the initial cells.
You have two cell types : mother cells and differentiated cells.

Random initial position

The input array must contain four integer values in this order.
  • \(MC_0\) : initial number of mother cells in the medium.
  • \(DC_0\) : initial number of differentiated cells in the medium.
  • \(sizeMatCell\) : cells box size (in pixels).
  • \(sizeMat\) : real box size (in pixels).

Warning : \(sizeMatCell\) must be larger than \(sizeMat\).

In the program, the input variable has the following pattern.
\(input = [MC_0, DC_0, sizeMatCell, sizeMat]\)

Here is an example.
  • \(MC_0 = 5 \phantom{t} (cells)\)
  • \(DC_0 = 0 \phantom{t} (cells)\)
  • \(sizeMatCell = 100 \phantom{t} (pixels)\)
  • \(sizeMat = 250 \phantom{t} (pixels)\)

Predetermined initial position

In order to make an easy predetermined initial position input, we wrote a code that understands images.
The input must contain the image path.
It is very easy to design a custom image i.e. a custom spatial cell distribution. You must follow these four rules (we work only with the image red component)
  • the image must be square and in jpeg format.
  • 255 means no cell.
  • 128 means one differentiated cell.
  • 0 means one mother cell.
0, 128 and 255 are red components of any colours that you want because the MATLAB code take into account only the image red component.
If you want, you can use the following hexadecimal colors codes.
  • #FFFFFF (white) means no cell.
  • #800000 (red) means one differentiated cell.
  • #000000 (black) means one mother cell.

Here is an example only with mother cells.

Example of predetermined cells initial position with mother cells

Figure 9 : Example of predetermined cells initial position with mother cells.


Here is another example with mother cells (black) and differentiated cells (red).

Example of predetermined cells initial position with mother cells (black) and differentiated cells (red)

Figure 10 : Example of predetermined cells initial position with mother cells (black) and differentiated cells (red).


Rate constants

You have to define seven variables.
  • \(\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.

In the program, the rate constants variable has the following pattern.
\(constantRate = [\mu_{1}, \sigma_{1}, \mu_{2}, \sigma_{2}, \mu_{3}, \sigma_{3}, k_{4}]\)

Here is an example.
  • \(\mu_{1} = 1 \phantom{t} (hours)\)
  • \(\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})\)

Fermentation period

The fermentation period \(t\) is a simple scalar variable.
You can use \(t = 5 \phantom{t} (hours)\) for example.

Action

In order to choose what you want to do, you must define three Boolean variables.
  • \(createFolder\) : create a folder in the disk to save information.
  • \(showAnimation\) : show the cell population animation. Warning : requires a lot of computing power.
  • \(plotGraph\) : plot the results in graphs.

In the program, the action variable has the following pattern.
\(action = [createFolder, showAnimation, plotGraph\)]

In the following example we want to create a folder and plot the graphs.
  • \(createFolder = 1\)
  • \(showAnimation = 0\)
  • \(plotGraph = 1\)

In the following example we want to show the cell animation.
  • \(createFolder = 0\)
  • \(showAnimation = 1\)
  • \(plotGraph = 0\)

In the following example we want to show the cell animation and save the images in a folder.
With the images you can create gif animations.
  • \(createFolder = 1\)
  • \(showAnimation = 1\)
  • \(plotGraph = 0\)

Folder name

The \(folderName\) variable is a string containing the folder name where results folder are stored in the disk.
You can use \(folderName =\) '\(Results\)' for example.

Output

Four types of output are generated.
  • raw data : all the useful information.
  • folder creation : if chosen in the \(action\) variable a folder is created.
  • graphs : if chosen in the \(action\) variable graphs are generated.
  • animation : if chosen in the \(action\) variable cell animation is generated.

Raw data

  • \(mainTic\) : start of the program.
  • \(MC_0\) : initial number of mother cells in the medium.
  • \(DC_0\) : initial number of differentiated cells in the medium.
  • \(pasDistributionTimeEvent\) : resolution for the graph time event distribution plot (see the \(initVariables\) function).
  • \(frameRate\) : animation frame rate (see the \(initVariables\) function).
  • \(timeShowPopulation\) : simulation time between two animations (see the \(initVariables\) function).
  • \(simulationTime\) : fermentation period \(t\).
  • \(resultFolderName\) : folder where the results folder are created.
  • \(\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.
  • \(cVitamin\) : vitamin counter.
  • \(time\) : main timer containing the simulation time.
  • \(createFolder\) : action boolean.
  • \(showPopulationBool\) : action boolean.
  • \(plotGraph\) : action boolean.
  • \(input\) : input array.
  • \(cTime1\) : time distribution counter 1 (see the \(randEvent\) function).
  • \(cTime2\) : time distribution counter 2 (see the \(randEvent\) function).
  • \(cTime3\) : time distribution counter 3 (see the \(randEvent\) function).
  • \(counterEvent\) : event counter.
  • \(timerShow\) : time animation counter.
  • \(sizeMat\) : real box size (in pixels).
  • \(nextEvent\) : next event type.
  • \(time1Array\) : time distribution array 1 (see the \(initTimeEvent\) function).
  • \(time2Array\) : time distribution array 2 (see the \(initTimeEvent\) function).
  • \(time3Array\) : time distribution array 3 (see the \(initTimeEvent\) function).
  • \(cell\) : structure containing all the cells' information.
  • \(sizeMatCell\) : cells box size (in pixels).
  • \(timeNextEvent\) : time before the next event.
  • \(saveTime\) : array containing the simulation time evolution.
  • \(saveNbrMC\) : array containing the mother cells time evolution.
  • \(saveNbrDC\) : array containing the differentiated cells time evolution.
  • \(saveNbrVitamin\) : array containing the vitamin time evolution.
  • \(folderName\) : folder where the results are stored. The pattern is the following : \(resultFolderName/DD-MM-YY\_hh\_mm\_ss\_mss\).
  • \(imgType\) : image type (see the \(initialization\) function).
  • \(parametersFileName\) : parameters file name where the parameters are saved (see the \(initialization\) function).
  • \(parametersDoc\) : file object where the parameters are saved.
  • \(maxMC\) : maximum number of mother cells.
  • \(maxDC\) : maximum number of differentiated cells.
  • \(maxVitamin\) : maximum number of vitamins.
  • \(timeSpent\) : computation time (in seconds).

Folder creation

If the \(action\) variable \(createFolder\) is set, the program will save important results in this folder.
  • \(parameters.txt\) : text file containing different simulation parameters and results.
  • \(data.mat\) : MATLAB file containing the raw data structure. To load the file, use this formula : \(load('data.mat')\).

Graphs

If the \(action\) variable \(plotGraph\) is set, the program will show the graphs.
Moreover, if the \(action\) variable \(createFolder\) is set, the program will save the graphs in the folder.

Animation

If the \(action\) variable \(showAnimation\) is set, the program will show the cell animation.
Moreover, if the \(action\) variable \(createFolder\) is set, the program will save the images in the folder.

Examples

We present here some results obtained with the program. For more details about the model sees the modeling section.

Time evolution of mother cells, differentiated cells and vitamin

Graphs

As explained in the modeling section the program give these results.

Stochastic evolution of the system

Figure 11 : Stochastic evolution of the system.



Stochastic vitamin optimization

Figure 12 : Stochastic vitamin optimization.



Distribution of event time probability

Figure 13 : Distribution of event time probability.

Animation

Here is an example with mother cells (black) and differentiated cells (red).
This image is the cell input.

Example of predetermined cells initial position with mother cells (black) and differentiated cells (red)

Figure 14 : Example of predetermined cells initial position with mother cells (black) and differentiated cells (red).


The result is the following animation after a conversion in a gif file.
The mother cells are in orange and the differentiated cells are in yellow.

Animation example with predetermined cells initial position with mother cells and differentiated cells

Figure 15 : Animation example with predetermined cells initial position with mother cells and differentiated cells.


Vitamin optimization

This program is designed to maximize the vitamin production. Here are some results.

Stochastic vitamin optimization

Figure 16 : Stochastic vitamin optimization.

Deterministic and stochastic vitamin optimization comparison

Here is a comparison between the two models.

Comparison of deterministic and stochastic vitamin optimization

Figure 17 : Comparison of deterministic and stochastic vitamin optimization.

Conclusion

Feel free to use and modify this program. The source code is available on GitHub : iGEM Paris Bettencourt 2015 under the GNU General Public License version 3.0. For more information concerning the model see the modeling section.