From PDE to tissue

In this post, I’ll show how to convert a simple 1D PDE model into a 2D multi-scale tissue model in a few steps.

Excitable media

Excitable media are a type of spatial nonlinear dynamical systems that have the ability to propagate waves and are distinguished by a refractory period, during which a new wave cannot pass. Think of a forest where a forest fire can spread (the wave), but a new fire of fire cannot pass until the tree have grown back (the refractory period).

A famous example of excitability in biological cells is the aggregation of social amoeba (Dictyostelium discoideum). In tissue biology, excitability is observed in the generation of a spike of transmembrane potential (action potential) by a neuron or cardiac cell, induced by a depolarizing perturbation of a resting state. In fact, a number of pathologies in the brain (epilepsia) and the heart (ventricular fibrillations, tachycardia) can be modelled as excitable media.

In this post, I will show how to construct a simple multi-scale model of an 3D excitable tissue, starting from a 1D PDE model.

1D PDE Barkley model

Derived from the seminal Hodgkin-Huxley and FitzHugh-Nagumo models, the Barkley model is probably the simplest continuous model for excitable media. It replaces the cubic term in the FitzHugh-Nagumo model with a piece-wise linear term as a simplification that enables fast simulation.

The equations for the Barkley model are: \begin{aligned} \frac{\partial u}{\partial t}&=\frac{1}{\epsilon}u(1-u) \cdot \left( u - \frac{v+b}{a} \right) + \nabla^2 u \\ \frac{\partial v}{\partial t}&=u-v \\ \end{aligned}

where the variable $u$ represents the fast dynamics and controls wave propagation and $v$ represents the slow dynamics controlling refractoriness. The parameter $\epsilon$ determines the timescale separation between the two dynamics, larger $a$ gives wider waves, increasing the ratio $b/a$ gives larger an excitability threshold.

To simulate this model with Morpheus, we first define two Fields for u and v in the Global section. These define scalar field, spatially resolved continuous variables over the whole lattice domain. Second, we define a System of differential equations (DiffEqn) for $\partial u/{\partial}t$ and $\partial v/{\partial}t$ and enter the reaction term of the equations above. Moreover, we specify the constants e, a and b. What is left is the diffusion rate for u, which is specified in Field/Diffusion.

See tutorial_excitable/1_1D_continuous.xml to try yourself.

Once the core of the model is set up, you’ll have defined a linear lattice with periodic boundary in Space and define the simulation duration in Time. To visualize the results, you can add a Logger and a Logger/Plot that records and visualized the values of u and v over space and time.

As initial condition, we let u be non-zero in a patch in the middle of the lattice, and v be non-zero in the left part. This forces the wave to propagate to the right, as shown below.

Discrete diffusion

Now, we want to reproduce this continuous spatial model with a discrete spatial model where each lattice site represents a biological cell. Later, we will expand this to 2D and add cellular dynamics using CPM.

To do this, we take our System as defined in Global before, but move it a CellType in the CellTypes section. You can just copy/paste this. Instead of a scalar field, we now want every cell to compute the reaction part for itself. That is, we convert the PDE to a ODE system.

To spatially couple the ODEs, we use a NeighborhoodReporter that looks at the value of u in the neighboring cells and return the arithmetic mean of these values, which is assigned to the new Property u_n.

In the reaction term of $du/dt$, we now add Du*(0.5*u_n - 0.5_u) which defines the discritized version of diffusion, which can be interpreted as an undirected cell-to-cell transport term. The property Du is theset to the same value as the diffusion rate.

See tutorial_excitable/2_1D_discrete.xml to try yourself.

The initial condition of the system consist of the spatial distribution of cells and their initial values in CellPopulations. You can use the InitCellLattice to fill the entire lattice with cells: one cell per lattice site. To specify the initial values, you use the InitProperty plugin and specify the values according to the cell center position cell.center and the lattice size s like this for u:

if( cell.center.x>=s.x/2-5 and cell.center.x<=s.x/2+5, 1, 0 )


And similarly for v:

if( cell.center.x<=s.x/2, 1, 0 )


To visualize the behavior of this system, you can use a Gnuplotter and colorcode the cells using the u and v properties and see something like this:

Convert to 2D

Converting the model above to a 2D model is trivial in Morpheus. In the Space section, you only need to change size and structure of the Lattice to e.g. 20,20,0 and square.

The most complicated part is to set up the initial conditions in order to get spiral waves. Here, we define u to be non-zero in a square in the middle of the lattice and v to be non-zero in the top-right quarter of the lattice. This suffices to obtain spirals. You can do this with an expression in InitProperty similarly to the one above.

Without further changes, you will see a simulation with popagating spiral waves, in a discrete lattice of square cells.

See tutorial_excitable/3_2D_discrete.xml to try yourself.

To move away from these regular and static cells and add motility and cell mechanics, we will now include cellular Potts sampling to our model.

This is actually surprisingly easy. In our existing CellTypes section, we only need to add VolumeConstraint and a SurfaceConstraint that together make up the energy fuction for our cellular Potts model (CPM) to maintain the certain ratio between the size and perimeters of the cells.

In addition, we need to specify some CPM-specific parameters in the CPM section of the model. In particular, we need to add a MonteCarloSampler and specify the temperature of the MetropolisKinetics which defines the amount of stochasticity in membrane fluctuations.

Another important parameter here is the MCSDuration. Typically, and by default, set to 1.0 to define one simulation time step as one Monte Carlo step (where on average every lattice points has been updated). Here, however, we want the simulation time to be determined by the dynamics of the Barkley model and scale the CPM model relative to this timescale. Here, we set the MCSDuration to 0.1 such that we simulate 10 CPM Monte Carlo steps for each simulation time step.

See tutorial_excitable/4_2D_CPM.xml to try yourself.

Since our cells will now become larger than a single lattice site, we obviously need to increase the size of the Lattice and change the initial configuration. We’ll use a InitRectangle to distribute our cell seeds over the lattice and InitVoronoi to render a Voronoi tesselation from these seeds.

Once set up, you will see a simulation like the one below. Note that the velocity and shape of the wave are not deterministic anymore because wave propagation is mediated by cell-cell contacts and these are governed by the stochastic cellular Potts model.

Couple signals to biomechanics

In the simulation above, the CPM dynamics affect the wave propagation. However, the wave itself does not affect CPM. Now, we’ll add some juice and realize a mutual coupling between the intracellular dynamics and tissue mechanics.

Here, we will assume two types of coupling. First, we assume cells contract based on the concentration of v. To implement this, we can use a mathematical expression instead of a value for the cell’s target area (VoumeConstraint/target):

Second, we assume that cells become more adherent to each other as a function of the concentration of u. This can be modelled using the HomophilicAdhesion plugin under CPM/Interaction/Contact. This modulates the cell-contact energy based on binding between two identical adhesives on adjacent cells, here represented by u.

Finally, to make the effect visually more dramatic, we do not simulate a confluent tissue, but only a fragment. I have done this with the InitCellObjects initializer that let’s you initialize cells are geometric objects such as boxes or circles.

As you can see, during the propagation of the wave, cell contract and become smaller, due to the effect of the intracellular model variables on the cell mechanics. On the level of the tissue, we see cells moving back and forth in wavy patterns, which is most clearly visible at the tissue boundary.

See tutorial_excitable/5_2D_coupled.xml to try yourself.

Conclusion

In this post, I have shown how to convert a 1D PDE model into a 2D multi-scale tissue model in Morpheus using 3 steps:

1. Convert a PDE model to a cell-discrete diffusion model,
2. Add motility and cell mechanics using cellular Potts sampling,
3. Couple intracellular dynamics and tissue mechanics.

I hope each step was clear enough to inspire you to experiment around!