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

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

.

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.

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.

# Add CPM cells

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.

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.

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

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

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