Paris/Stochastic model

From 2007.igem.org

• Main Page
• Introduction
• Design Process
• Modeling
• Construction Process
• Results
• Perspectives
• Acknowledgments

In this last part of the models section, we are developing a stochastic simulation of the microscopic model. The major contribution is to handle in a stochastic context a dynamic and heterogeneous population of bacteria.

Introduction

In 1977, Gillespie developed an exact Simulation Stochastic 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 approach 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) that 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 [http://dx.doi.org/10.1016/j.biosystems.2006.12.009 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 compartment. 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:

 DAPAp + DAP ->K_DAPApI DAPAp_i
 DAPAp_i     ->K_DAPApA DAPAp + DAP
 DAPAp       ->K_Cre    DAPAp + Cre

DAPAp promoter activity for independent differentiation:

 DAPAp       ->K_Cre    DAPAp + Cre

Differentiation:

 LOXP_Box    + Cre ->K_Diff LOXP_BoxCre
 LOXP_BoxCre + Cre ->K_Diff DAP_Box

DAP production:

 DAP_Box ->K_DAPiP DAP_Box + DAP

DAP import and export:

 [DAP, ...]Bact  ->K_DAPEx DAP + [...]Bact
 DAP + [...]Bact ->K_DAPIm [DAP, ...]Bact

Bacteria macroscopic evolutions (division and death):

 [nDAP, ...]BactG ->K_Div   [(n/2)DAP, ...]sub>BactG</sub> + [(n/2)DAP, ...]sub>BactG</sub>
 [nDAP, ...]Bact  ->K_Death nDAP

MGS Implementation

As presented [http://dx.doi.org/10.1016/j.biosystems.2006.12.009 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 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 ne 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 [http://dx.doi.org/10.1016/j.biosystems.2006.12.009 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 simulation 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