Paris/Stochastic model
From 2007.igem.org
In this last part of the models section, we are developing a stochastic simulation of the microscopic model. The major contribution is to handle a dynamic and heterogeneous population of bacteria in a stochastic context.
Contents 
Introduction
In 1977, Gillespie developed an exact Stochastic Simulation Algorithm (SSA) dedicated to the simulation of homogeneous chemical systems. This method was recently used in many applications for the simulation of biological systems. A good point of this approach is that it allows to handle biochemical systems where numbers of molecules are low and that cannot be well characterized by classical approaches using differential equations and chemical concentrations. Nevertheless this method requires strong hypotheses about the spatial homogeneity of molecules distribution. Extensions of Gillespie's SSA have been proposed to deal with compartments.
As our system is composed of a growing and heterogeneous population of bacteria, we propose to use this extension to simulate it. In the following paragraphs, we first detail the extended SSA we use and then we present some samples generated by our implementation using the set of parameters found in the numerical analysis of the model. Note that the main contribution here is in the development of the simulation algorithm.
Extended SSA
Gillespie's SSA
From a computational point of view, the Gillespie SSA relies on a discrete events simulation of chemical reactions between individual molecules. A reaction R_{mu}
like A + B > C
occurs when reactants A
and B
meet with enough energy to produce a molecule C
. The probability that this reaction occurs during an infinitesimal time, is proportional to the number of molecules A
and B
in the system (dependence on the concentration) and a coefficient c_{mu}
called the stochastic constant (corresponding to the reaction kinetic) when reactants are uniformly distributed in the system (intuitively it means that each couple of molecules has the same probability to meet). Assuming a system composing of molecules that can interact with respect to reactions R_{1}
... R_{N}
, this probability allows to compute the probability law p(mu
,tau
) of the next reaction to occur is R_{mu}
after tau
units of time.
The algorithm developed by Gillespie consists in iterating the drawing of the next reaction mu
to occur together with its reaction time tau
using p(mu
,tau
), and in making the system evolves. Assuming that S
is chemical solution where reactions R_{0}
... R_{N}
can occur at a given time t
, the SSA can be expressed as follows:
SSA (S,t,R_{0},...,R_{N}) = {
compute mu and tau using p(mu,tau)
compute S' by applying mu on S
return (S',t+tau)
}
Handling Membranes
In the SSA, molecules have to be uniformly distributed in space. It cannot be straightly used to deal with systems that exhibit a complex spatial organization, as the system we are interested in. Taking space into account means modifying the probability law p(mu,tau) in order to take care of the localization of the reaction. Obviously, it cannot be computed in practice for any kind of organization, but it is possible to develop adhoc algorithms for specific organizations.
We propose to consider a population of bacteria as a nested membranes system, called compartments. The environment where bacteria live is represented by a compartment that contains molecules and bacteria, and a bacterium is then represented by a compartment that contains molecules. In each compartment, the contained elements are assumed to be uniformly distributed in space.
Extended Algorithm
The extension relies on the fact that compartments elements are uniformly organized, meaning the SSA can be used in each compartment. Intuitively, at each iteration of the new algorithm, the standard SSA is evaluated on each compartment of the system ; the compartment with the smallest reaction time is updated (with respect to the corresponding reaction) and the other compartments remain unchanged. This procedure is naturally recursive and can be expressed as follows:
ExtSSA (S,t,R_{0},...,R_{N}) = {
(S',t') := SSA (S,t,R_{0},...,R_{N})
foreach compartment M of S do
let (M',t_{M}) = ExtSSA (M,t,R_{0},...,R_{N}) in
if t_{M} < t'
then
S' := replace M by M' in S
t' := t_{M}
done
return (S',t')
}
This algorithm is taken from here.
SMB in Terms of Compartments and Chemical Reactions
We start from the same representation of the system as presented in the ODE based simulation. Obviously, bacteria are represented by compartments that contain molecules as shown on the picture of the previous section. Let now consider the figure next to this text in order to sum up the different molecules we have to deal with:
 DAP: it is of course one of the major components of the system. It can be located in the bacteria or out of them in the environment.
 Cre: it is the second main molecule to be considered. It remains in cells and cannot be exported or imported.
 DNA: we focus only on the two parts of the DNA we are interested in:
 Both types of cell are able to produce Cre. This encoding region is under the control of a constitutive promoter. We then consider a molecule DAPAp that is an abstract representation of this part of the DNA. Thus every cell contains a DAPAp molecule. Following the two possible behaviors we consider, DAPAp can be inhibited or not by DAP (here is controlled the DAP (in)dependence of the differentiation). If the case occurs, the inhibited form of DAPAp is named DAPAp_i.
 The germ cells differentiate using a Cre controlled recombination of a DNA part. This part have two different states as the considered cell is either differentiated (somatic cell) or not (germ cells). See here for a biological description of this mechanism. We can call these states LOX_Box for germ cells (the DNA part uses a LOXCre mechanism to recombine) and DAP_Box for somatic cells (after recombination the DNA is able to produce DAP). Once again each cell has to contain one molecule of DAP_Box or one molecule of LOX_Box.
The other components are of course bacteria compartments. The two kinds of compartment differ only by the presence of either LOX_Box for germ cells or DNA_Box for soma cells. Using these molecules, we are able to model the dynamics by chemical reactions.
DAP and Cre degradation:
Cre >_{K_CreD} .
DAP >_{K_DAPD} .
DAPAp promoter activity for DAP dependence differentiation: these rules prevent the production of Cre when DAP inhibits DAPAp
DAPAp + DAP >_{K_DAPApI} DAPAp_i
DAPAp_i >_{K_DAPApA} DAPAp + DAP
DAPAp >_{K_Cre } DAPAp + Cre
DAPAp promoter activity for independent differentiation: DAPAp is never inhibited
DAPAp >_{K_Cre } DAPAp + Cre
Differentiation: note that the differentiation is irreversible. Molecules of Cre have to bind the both LOX sites of the LOX_BOX. We name LOXP_Box_{Cre}
the intermediary state
LOXP_Box_{ } + Cre >_{K_Diff} LOXP_Box_{Cre}
LOXP_Box_{Cre} + Cre >_{K_Diff} DAP_Box
DAP production: this only happens in somatic cells, germ cell does not contain DAP_Box. Therefore, only somatic cells are to feed germ cells
DAP_Box >_{K_DAPiP} DAP_Box + DAP
DAP import and export: these rules are applied in the environment. The square brackets represent compartments (for example [...]_{BactS}
is a somatic cell)
[DAP, ...]_{Bact} >_{K_DAPEx} DAP + [...]_{Bact}
DAP + [...]_{Bact} >_{K_DAPIm} [DAP, ...]_{Bact}
Bacteria macroscopic evolutions (division and death): when a bacterium dies, its internal DAP molecules are released in the environment. As DAP is essential for germ cells division, we can envisage that the value of the kinetic constant K_Div
is a function of the number n
of DAP molecules in the cell
[n.DAP, ...]_{BactG} >_{K_Div } [(n/2).DAP, ...]_{BactG} + [(n/2).DAP, ...]_{BactG}
[n.DAP, ...]_{Bact } >_{K_Death} n.DAP
MGS Implementation
As presented here, the MGS language provides a specific rules application strategy implementing the SSA. We propose to implement the algorithm ExtSSA
using MGS. The program (available here) is briefly described with the representation of the system, the encoding of chemical reactions and a proposition of optimization for ExtSSA
.
State and Structure
First of all, our system handle different types of molecules. They are represented in MGS by symbols (back quoted strings):

`DAP
: DAP molecules

`Cre
: Cre molecules

`DAPAp
: DAP sensitive promoter (inhibited form:`DAPAp_i
)

`LOXP_Box
: DNA configuration of germ cells (Cre bounding form:`LOXP_Box_Cre
)

`DAP_Box
: DNA configuration of somatic cells
As mentioned above, the system is represented by nested compartments. MGS provides the socalled bag, a particular data structure corresponding to multisets (sets where multiple occurrences of the same element are allowed) whose topology is a complete graph. In other words, each element of the bag is neighbor with all the others. From the dynamical systems point of view, the neighborhood is equivalent to a potential interaction: in bags, each element can react with all the others. This topology is then perfect to represent wellmixed chemical solutions and, of course, compartments. The system is composed of 3 types of compartments defined here:
 the environment of type
Env
where all bacteria are localized
 the germ cells of type
BactG
that must contain one molecule`LOXP_Box
 the somatic cells of type
BactS
that must contain one molecule`DAP_Box
The value (`DAP::`DAP_Box::BactS:())::(`LOXP_Box::BactG:())::`DAP::Env:()
is an example of a system state. It is composed of two bacteria, a germ cell and a somatic cell (the somatic contains one molecule of DAP), and an outside diffusing DAP molecule.
Dynamics
The chemical reactions are specified using transformation rules (see here). There are 14 different reactions embedding in two different transformations:
 transformation
inBact
: this transformation specifies the chemical reactions occurring in bacteria (independently if they are germ or somatic). The reactions are chemicals degradation, DAPAp promoter inhibition, differentiation and DAP production (for somatic cells only).
 transformation
inEnv
: this transformation specifies the chemical reactions occurring in the environment. Reactions are DAP degradation, DAP import and export (this rule actually appears ininBact
to improve the efficiency of the program and to noticeably reduce the execution time).
As an example, the following rule appears in inBact
:
`DAPAp, `DAP ={ C = sndOrder(K_DAPApI,1.0) }=> `DAPAp_i ;
It specifies the inhibition of the promoter DAPAp by a molecule of DAP: `DAPAp
and `DAP
reacts to produce `DAPAp_i
, the inhibited form of `DAPAp
. In the arrow, the parameter C
defines the stoachstic constant of this inhibition. It is computed using the function sndOrder
that takes as arguments

K_DAPApI
, the kinetic constant of this reaction (whose value is defined here), and
 the volume of the cell (we assume here that the volume is always 1.0 ; volume dependence is not implemented yet).
Optimization
The previously described ExtSSA
is not straightly implemented. In fact, while stochastic simulations can deal with a small number of molecules, they often exhibit limitations as they are time greedy. On the first hand, events are triggered asynchronously: only one reaction is applied at each time step. As a consequence, a simulation requires a lot of time steps to be significant. On the other hand, choosing the next reaction requires to draw a lot of random numbers. Unfortunately, this operation is badly time consuming and increases the execution time.
We propose to decrease the number of random numbers generations by freezing cells that do not evolve. This improvement consists in using a scheduler that stocks cells with their reaction time in a buffer (the program is available at here). At each time step, the SSA is evaluated on the environment. If the returning reaction time is smaller than the smallest reaction time of the buffer, then we make the environment evolves and reaction times of cells touched by this update are recomputed and the buffer updated. On the other, if the environment evolution is slower, the fastest reaction of the buffer is applied and the corresponding cell is updated. This improved version of ExtSSA
is programmed here.
To give an idea of the efficiency of our improved algorithm, the evaluation of a 100000 steps evolution starting with a population of 20 cells to reach a population size of about 230, is about 290s on a standard computer. Moreover, note that the MGS toplevel is a prototyped and generic system that clearly not provides the efficiency of an adhoc programming of our improved ExtSSA
.
Some Results
The following figures was generated by the implementation of ExtSSA
we have just described. We start with 20 germ cells without external DAP production. On the left, the temporal evolution of the bacteria population size. The exponential growth of the population is once again stressed by the simulation. On the left, the graph presents the mean number of DAP molecules in a cell. As we can see, this never goes above 5 molecules. It justifies
the use of a stochastic process for realistic simulations.
Conclusion
We have developed a complex and generic framework for stochastic simulations of nested membranes systems following ideas presented here. Moreover, our implementation provides an optimization that allows to use this framework in a lifesize project.
On the other hand, the simulation is not volumedependent. As we do not take into account the evolution of the cells size, reactions rates are not yet well computed. This is a work in progress. Our simulation is also limited because of the homogeneous assumption of ExtSSA
in compartments: the somatic cells feed all the germ cells. But as previous simulations show, somatic cells only feed germ cells that surround them. The spatial distribution of cells is not here taken into account. A possible extension is to merge the three simulations to deal with all the aspects of the system:
 the cellular automaton to simulate DAP diffusion in the environment
 the spatial model to simulate the spatial organization of cells
 the Gillespie based simulation to simulate the chemical reactions occurring in cells