# Paris/Stochastic model

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.

## 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 `Rmu` 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 `cmu` 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 `R1` ... `RN`, this probability allows to compute the probability law p(`mu`,`tau`) of the next reaction to occur is `Rmu` 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 `R0` ... `RN` can occur at a given time `t`, the SSA can be expressed as follows:

``` ```

``` SSA (S,t,R0,...,RN) = { 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 ad-hoc 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,R0,...,RN) = { (S',t') := SSA (S,t,R0,...,RN) foreach compartment M of S do let (M',tM) = ExtSSA (M,t,R0,...,RN) in if tM < t' then S' := replace M by M' in S t' := tM 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 LOX-Cre 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_BoxCre` the intermediary state ``` ```

``` LOXP_Box + Cre ->K_Diff LOXP_BoxCre LOXP_BoxCre + 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 so-called bag, a particular data structure corresponding to multi-sets (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 well-mixed 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 in `inBact` 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 top-level is a prototyped and generic system that clearly not provides the efficiency of an ad-hoc 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 life-size project.

On the other hand, the simulation is not volume-dependent. 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 spatial model to simulate the spatial organization of cells
• the Gillespie based simulation to simulate the chemical reactions occurring in cells