Motivation

Until now, intercellular interactions has been difficult to model. And because we also use microfluidics as a tool to confirm our designs, we made a decision to create an instrument which allows easy and completely boundless simulation of cell communication in any microfluidic environment. We decided to name the instrument CeCe, as abbreviation of Cell-Cell interaction.

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) and use coefficient to convert units between those two systems. The reason why are those time step separated is the simulator 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.

Demo

Modeling

Default simulator package contains several plugins that offer some additional functionality. 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 current 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 - Lattice-Boltzman method

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 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