In this example, we will use Multicell to simulate the self-organization of a geometrical Turing pattern (Turing 1952; Note about other proposals and ways to produce spatial patterns), based on equations developed by Gierer and Meinhardt (Gierer and Meinhardt 1972). These equations describe a simple molecular mechanism that involves two chemical species, an activator and a repressor. The activator activates itself, as well as the repressor. The repressor represses the activator. Both species are diffusible, the activator within a short-range, and the repressor within a longer range.
Despite its simplicity, this mechanism has been successfully used to explain the formation of many different molecular spatial patterns in tissue structures (Meinhardt and Gierer 1974; Meinhardt book 1982; Kondo 2010;). In this section, we will implement the Gierer-Meinhardt equations (Gierer and Meinhardt 1972).
In [1]:
%matplotlib notebook
In [2]:
import multicell
import numpy as np
In [3]:
sim = multicell.simulation_builder.generate_cell_grid_sim(20, 20, 1, 1e-3)
In [4]:
sim.register_cell_variable("a")
sim.register_cell_variable("h")
The concentrations of a will be computed automatically for all cells. However, we will be going to use their squares multiple times per time step. To avoid raising the vector c_a to the square multiple times, we define a computed variable c_a2 that will be computed once per time step. The equation of c_a2 is defined using a Python function, which is then registered using the method register_computed_variable of the Simulation object.
In [5]:
def c_a2(c_a, **kwargs):
return c_a**2
sim.register_computed_variable("c_a2", c_a2)
In [6]:
sim.set_constants({"mu_a": 1e-1, "mu_h": 2e-1, "rho_a": 1., "rho_h": 1., "q": 1., "H": 0.35, "A": 0., "D_h": 5., "D_a": 0.025})
In [7]:
def da_dt(simulation, a, c_a, c_a2, c_h, D_a, mu_a, rho_a, A, q, adjacency_matrix, **kwargs):
return simulation.diffusion(D_a, c_a, adjacency_matrix) + rho_a * c_a2 / c_h / (1 + q**2 * c_a2) - mu_a * a + A
The formula of $\dfrac{dh}{dt}$ is similarly built, except for the fact that the variable synthesis term is only a-dependent.
In [8]:
def dh_dt(simulation, h, c_a2, c_h, D_h, mu_h, rho_h, H, adjacency_matrix, **kwargs):
return simulation.diffusion(D_h, c_h, adjacency_matrix) + rho_h * c_a2 - mu_h * h + H
sim.set_ODE("a", da_dt)
sim.set_ODE("h", dh_dt)
In [9]:
sim.initialize_cell_variables()
We initialize initial quantities of matter to values that would be close to their steady state if there was no diffusion (i.e. if cells were independent). As cell volumes are all slightly different (the grid is noisy), concentrations will also all be slightly different and there is no need to randomize the initial quantities of matter.
In [10]:
a0 = np.full(sim.n_cells, 0.789)
h0 = np.full(sim.n_cells, 4.863)
sim.set_cell_variable("a", a0)
sim.set_cell_variable("h", h0)
In [11]:
sim.set_duration(3200)
sim.set_time_steps(10, "linear")
In [12]:
sim.register_renderer(multicell.rendering.MatplotlibRenderer, "c_a", {"max_cmap": 1.3, "view": (90, -90), "axes": False})
In [13]:
sim.renderer.display("c_a")
In [14]:
sim.simulate()
In [15]:
sim.set_duration(1e7)
sim.set_time_steps(1)
sim.simulate()
Small differences in concentrations amplified and propagated into a pattern.