Tutorial 2 - cellmodeller/CellModeller GitHub Wiki

In this tutorial, we will cover intracellular ODE solving - allowing for the simulation of chemical and gene dynamics inside cells.

Tutorial 2a: simple gene expression

In this tutorial, let's look at simple expression of a single gene at a constant rate. Let's open Examples/Tutorial_2/Tutorial_2a.py, and let's run this model in the GUI.

We have now imported a simple integrator into this model: from CellModeller.Integration.CLEulerIntegrator import CLEulerIntegrator, the CLEulerIntegrator. This module deals with ODE solving inside cells.

Inside setup, we now have to initialize this module, and link it to the simulator:

integ = CLEulerIntegrator(sim, 1, max_cells)
...
sim.init(biophys, regul, None, integ)

The integrator must be initialized with the 3 arguments: firstly the simulator object, sim in order to connect the module to the simulator, then the number of chemical species being simulated, in this case 1, and the maximum number of cells. Note that the max_cells is also passed into the biophysics module too.

Let's say that this is a protein, and that each cell produces the protein at the same, constant rate. We specify the initial concentration of our 1 protein in the init function. For visualization, the red colour of the cells is linked to the species in the update function with cell.color = [numpy.clip(cell.species[0]/6.0,0.0,1.0), 1.0, 0.1], so that green cells have little protein, and yellow have lots.

The actual ODEs are specified in the specRateCL() function. Note that this code is passed into an OpenCL file which is compiled at runtime, and thus is actually written in C and not in python (hence you must use semicolons ; on each line). This means that the code can be solved on the GPU. Let's look at the equations:

const float k1 = 2.f;
float x0 = species[0];
rates[0] = k1;

So, we have a very simple system here. We have one species, x0, produced inside every cell at an arbitrary rate of 2. The integrator also automatically dilutes the species based on the actual volume growth of each cell. Without degradation, dilution balances this production rate to establish an equilibrium based on the effective growth rate for each cell.

This equilibrium between growth rate and production leads to emergent spatial patterns within the colony: Tutorial 2a

As the cells in the center of the colony are prevented from growing quickly by physical forces, they do dilute the protein and build up the concentration.

As well as species concentrations, various variables can be passed into the species rate calculation and can be used with the following names:
Cell position pos[n] (where n=0,1,2 for x,y,z)
Cell area area
Cell volume volume
Cell type cellType
Effective cell growth effGrowth (this is the actual amount grown averaged over the cells lifetime)

Using this we can change these emergent spatial patterns. For example, if the multiply the production rate by the effective growth, the protein concentration becomes much more homogenous across a colony.

Tutorial 2b: simple oscillations

In this tutorial, we will make a simple oscillator inside cells to demonstrate one possible increase in complexity in the intracellular network. Let's open Examples/Tutorial_2/Tutorial_2b.py and run it in the GUI.

This model is set up in a similar way to the previous tutorial, except now we have 2 species initialized in the integrator.

integ = CLEulerIntegrator(sim, 2, max_cells)

In this case, the cells contain a simple 2 component oscillator, taken from Guantes and Poyatos PLoS Comput Biol. 2006 Mar; 2(3). (A full derivation and explanation of the equations can be found in the reference).

Tutorial 2b topology

The circuit topology is shown in the figure above. This network is made up of two genes, one producing an activator, A, and an inhibitor, B. Protein A upregulates both components, and B which downregulates the activator. The equations solved in the model are also shown in the figure above. Note that the protein dilution from growth adds extra issues in analyzing how these equations behave in bacteria, which we can explicitly simulate here. The ‘specRate()’ function specifies the equation in the model, and is in the model file like this:

def specRateCL():
    return '''
    const float k1 = 2.f;
    const float k2 = 2.f;
    const float delta = 1.f;
    const float rho = 1.f;
    const float sigma = 1.f;
    
    float x = species[0];
    float y = species[1];
    rates[0] = delta*(k1*(1+rho*x*x)/(1+x*x+sigma*y*y)-x);
    rates[1] = delta*k2*(1+rho*x*x)/(1+x*x)-y;
    '''

Here we see the equations in the figure above written out explicitly, with chemicals x and y as species[0] and [1], and their rates, with constants set to arbitrary values.

If you run the model, you will see the cells oscillate. As the simulation progresses, the oscillations become increasingly out of sync due to stochastic effects changing effective growth rates and decoupling the oscillations. Tutorial 2b

By linking oscillating proteins to growth rates (eg. by setting cell.growthRate = cell.species[0] in the update function) it is possible to generate interesting feedbacks between dilution and the gene network that can lead to interesting effects.

Now, hopefully you are comfortable with simulating intracellular chemical dynamics in CellModeller. Let's move on to Tutorial 3, where we look at intercellular chemical dynamics.

⚠️ **GitHub.com Fallback** ⚠️