Software

Demo

Demo 1

We created a demo simulation to present some of the CeCe simulator abilities. This simulation contains two different paths for input cells with different size. As you can see, bottom flow velocity is lower than in the upper path so yeast cells have time to produce enough signal and react to it. The connection between tanks shows how the simulator handle difference of pressure.

Demo 2

The second demo shows how yeast colony cells reacts to contaminous cells which produce specific signal. Contaminous cells are added randomly in specific area and starts producing signal (not visible). When yeast cells detect some specific amount of signal in surrounding environment they start to produce RFP molecules.

This simulation was requested by Chalmers iGEM team in order to confirm their Study in the Scarlett detection design.

Download

Simulate yourself, both simulation files can be downloaded here.

Architecture

The simulator is designed to be highly modular. The simulator core contains almost no functionality and missing functionality is provided by plugins. This design allows to extend simulator functionality by adding plugins written by someone else. Plugins are loaded on demand by simulation file so unnecessary functionality is not used.

Simulation is computed by stepping over independent iterations. Each iteration is defined by time step and all modules and objects (provided by plugins) are updated. Module is unit that is updated in each iteration and provides some global functionality (e.g. diffusion). Only one of same type is allowed to be in simulation. Objects are simulation entities that provides some local functionality (e.g. Yeast cells).

Objects within simulation are driven by physical engine ([http://box2d.org Box2D]) so interactions between objects looks almost realistic. Physical engine use different time step (engine is designed for games and mostly expect simulation step around 33ms) than is used by simulation and physical engine uses coefficient to convert units between those two systems. The reason why are those time step separated is simulations can require larger time step (seconds, minutes, ...) - cell growth is slow process that takes hours and have same time step doesn't make sense in meaning of simulation.

Visualization

The CeCe simulator has ability to visualize given simulation (OpenGL). Final scene visualization is done by combination of modules and objects visualization (if they provide it). Visualization is pretty fast so powerful graphics card is not required even integrated GPU is able to render visualization in real time.

Modeling

Default simulator package contains several plugins that offer some additional functionality which is described below. More detailed overview how the simulator has been used in our project can be found on modeling page.

Diffusion

The diffusion plugin is used to simulate diffusion where regions with high concentration are spread to regions with lower concentration. Our model do this by distributing grid cell concentration between surrounding grid cells. Following formula are used to compute signal distribution matrix (in our case it's 3x3):

$$M_{i,j} = \frac{1}{4 \pi D dt} \cdot e^{\frac{d_{i,j}^2}{4 D dt}}\;i,\,j = 0, 1, 2$$

where \(M\) is mapping matrix, \(D\) is diffusion rate, \(dt\) is time step and \(d_{i,j}\) is cell distance from source cell (\(M_{1,1}\)). For each diffusion grid cell the \(M\) matrix is used to distribute cell signal between surrounding grid cells.

$$G_{x,y+i,j-1,1}^{n+1} = G_{x,y-1,1}^{n} (1 - R_d dt) \bar{M_{i,j}},\;i,\,j = 0, 1, 2;\;x = 0, ... N;\;y = 0, ... M$$

where \(G\) is diffusion grid, \(R_d\) is degradation rate, \(\bar{M_{i,j}}\) is normalized \(M_{i,j}\) and \(N\), \(M\) are size of diffusion grid.

Streamlines

For simulating complex fluid systems the Lattice Boltzman method became de facto standard last years. This relative new method works in two steps:

1. Collision step
2. Streaming step

Collision step compute all distribution functions of node's mass and streaming step move that specific mass to other node (in defined direction).

The simulator implements square lattice LBGK model (D2Q9) which is expressed as (collision step):

$$f_i(\bar{x} + \bar{e_i}, t + 1) - f_i(\bar{x}, t) = - \frac{1}{\tau}[f_i(\bar{x}, t) - f^{(eq)}_i(\bar{x}, t)], i = 0, 1, ..., 8$$

where \(f_i\) is distribution function in direction \(i\), \(\bar{x}\) is current position, \(\bar{e_i}\) is position change in direction \(i\), \(\tau\) is relaxation time and \(f^{(eq)}_i\) is equilibrium distribution function in direction \(i\). The direction \(i\) represents 2D direction of momentum. Node density \(\rho\) and macroscopic flow velocity \(\bar{u}\) can be calculated:

$$\sum_{i=0}^{8}f_i = \rho$$ $$\sum_{i=1}^{8}f_i e_i = \rho \bar{u}$$

The equilibrium distribution functions \(f^{(eq)}_i\) depends only on local density and velocity. In our case we chose to calculate it in following way:

$$f^{(eq)}_i = \rho \cdot w_i [1 - 3(\bar{e_i} \cdot \bar{u}) + \frac{9}{2}(\bar{e_i} \cdot \bar{u})^2 - \frac{3}{2} \bar{u} \cdot \bar{u}]$$

where \(w_i\) is direction weight:

$$w_i = \begin{cases} \frac{4}{9} & i = 0\\ \frac{1}{9} & i = 1,2,3,4\\ \frac{1}{36} & i = 5,6,7,8 \end{cases}$$

After collision step is performed streaming step takes place.

$$f_i(\bar{x} + \bar{e_i}, t + 1) = f_i(\bar{x}, t + 1)$$

This means the distribution function value is streamed to next node in that specific direction.

Boundary conditions

Boundary conditions are a quite complex problem. The problem arises from the fact that there no physical intuition on the behaviour of the velocities on boundaries. Choosing the right boundary conditions play an important role in simulation stability. In our simulator we implement boundary conditions defined by Zou and He^[ZouHe1997].

Inlet

At inlet boundary condition we have no information about behaviour before inlet. For each inlet node we compute density and then use this density and inlet velocity \(\bar{u}\) to compute distribution functions.

$$\rho = \frac{1}{1 - \bar{u}_i} (f_{0,2,4} + 2 f_{3,6,7})$$

Outlet

At outlet boundary condition there is no information about behaviour after outlet. We compute outlet velocity and then use this velocity with \(\rho = 1\) to compute distribution functions.

$$u = -1 + \frac{1}{\rho} (f_{0,2,4} + 2 f_{1,5,8})$$

Stochastic reactions

In stochastic reaction each reaction is defined by chemical components and a rate. To simulate a reaction, you have to know how much of each reactant you have plus how probable it is for the reaction to happen. The chemical components can be divided into two classes:

• Non-diffusive
• Diffusive

Non diffusive components are proteins, markers and genes which do not leave the cell. They may be treated as being present only inside the cell’s center, i.e., programmatically speaking, each biological cell remembers how many molecules of each component there are.

Diffusive components are a bit tricky, because they are not linked to the cell interior but rather to the environment itself. During each step a routine is run to determine how many molecules of a diffusive component are currently present in close cell's surrounings.

In some simulations you may need some reactions inside the cell, but not diffusive components. In order to save your resources diffusive reactions are an optional extension, which may or may not included in your simulation.

Propensity

Each reaction is described by the reactants, the rate, and when talking stochastically, the probability that the reaction happens in the following moments. Usually, a reaction is described by a probability measure called propensity, which is the probability that the reaction fires in the infinitesimal time dt. To simulate a system of chemical reactions, we need to know how to compute the propensity for simple reactions

Birth

\(\varnothing \Rightarrow k_A \Rightarrow A\), propensity \(k_A\)

Death

\(A \Rightarrow d_A \Rightarrow \varnothing\), propensity \(d_A A\)

Two component

\(A + B \Rightarrow k_{AB} \Rightarrow C\), propensity \(k_{AB} A B\)

Probability of birth reaction is constant, it depends only on its rate, but probability that something dies increase with the amount of things that can die (The more A's, the higher probability of any of them dying) Similarly, the probability that something collides increases with amounts of both things that can collide.

It is often useful to write down reactions in terms of stoichiometric matrix. This matrix’s rows correspond with reactions and columns with species. Each row says how much of a change in species numbers there is if the linked reaction happens. For example consider the following reaction system:

• \(\varnothing \Rightarrow k_A \Rightarrow A\)
• \(B \Rightarrow d_B \Rightarrow \varnothing\)
• \(A + B \Rightarrow k_{AB} \Rightarrow C\)

The underlying stoichiometric matrix is $$\begin{bmatrix}+1 & 0 & 0\\0 & -1 & 0\\ -1 & -1 & +1\end{bmatrix}$$ The number at position i, j says how much of a change in number of j-th species molecules there is when the i-th reaction fires.

Gillespie First reaction method

It can be proved, that under well-mixed solution and thermal equilibrium assumptions the reaction times have an exponential distribution in time. The simplest way to simulate a stochastic system is the Gillespie first reaction method. This methods consists of the following step in each iteration

1. Compute propensities of all reactions in the system
2. Generate a uniform random variable (0,1)
3. Calculate the random time of next reaction from the following formula $$ \tau = -\dfrac{1}{\sum\nolimits Propenisty} log[rand[0,1]]$$
4. Determine which reaction happened – the index of reaction is the index of the first entry in the cumulative sum of propensities that is greater than $$rand[0,1] \sum\nolimits Propensity$$
5. Update the system’s state according to the stoichiometry matrix

Tau-leaping and Dependency graph

Tau leaping is an enhancement of the Gillespie algorithm, where propensities are updated only once per fixed timestep instead of after each reaction. This method makes Gillespie algorithm more efficient in simulation of larger systems. On the other hand, this approximation doesn't fully satisfy the idea of Gillespie algorithm, because some computations are based on outdated propensities.

To have an efficient simulation together with up-to-date propensities we decided to include the Dependency graph. This improvement makes recomputation of propensities more effective, because propensities are no longer updated all by once, but only those which got changed when particular reaction has occured. For example, when reaction with expression of A has happened, only reactions which propensity is function of A are updated. Dependency graph removes disadvantages of tau-leaping and keeps propensities accurate.

Download

Simulator binaries are available to download for following platforms:

• Windows 64bit
• Mac OS X
• Ubuntu/Mint x64

Package contains simulator binary, examples, manual and tutorial (HTML).

Source code

The CeCe simulator source code is available on [http://github.com/GeorgievLab/CeCe GitHub].

References

1. ↑ Qisu Zou, Xiaoyi He (1997) On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Phys. Fluids, 9, 1592-1598