Biological Patterning

When you look at patterning for animals we typically think of it with respect to stripes or spots. But why are some things striped and other things spotted? We will look into some models that help describe these regular and self-repairing patterns that we see in the animals around us (and at the zoo). Hopefully we can build an intuition about what causes these patterns and make some cool pictures along the way.

Reaction Diffusion Models

The reaction-diffusion model is a theoretical model used to explain self-regulated pattern formation in biology. It is one of the best-known theoretical models to explain self-regulated pattern formation in the developing animal embryo as well as expressed patterns on the fur. These models can create many spatial patterns that are seen in the natural world. The reaction-diffusion model describes a simple system of 2 substances interacting with each other to generate emergent patterns.

The Activator substance increases its own activity as well as the activity of the inhibitor. The inhibitor substance decreases the activity of the activator. Because these substances diffuse at different speeds, patterns can be created. The behavior of these systems depends on the diffusion properties of the individual morphogens. The activator typically diffuses slower than the inhibitor. This system can replicate most biological patterns seen in nature and is a standard theory for biological pattern formation.

In simulation studies, classic gradient boundary conditions are used to add realism to the generated pattern. With this model, an initially homogenous concentration can form repeating patterns. It is estimated that 2-component Reaction-Diffusion models are insufficient to explain complex phenotype expressions in biology. By adding a second substance, the diffusivity of the system is no longer a problem, and bistable patterns can be created. With these models, the cells don’t have positional information about their place, but are still able to form predictable patterns. In reaction-diffusion models, ligands can diffuse via active diffusion to control the morphogens range.

Reaction-diffusion models have been shown to explain the signaling of morphogens across a large range of distances. These distances range from 5um during spindle assembly to 200um for embryo development in zebrafish. They also show non-equilibrium behaviors with a stable source of energy, such as the BZ reaction. Diffusion can couple isolated parts for oscillators with non-trivial phase differences.

Turing Patterns

Turing patterns are named after mathematician Alan Turing, who first proposed the idea that these patterns could arise spontaneously from the interaction of chemical substances in the 1950s. In the years since, Turing patterns have been observed in a wide range of biological systems, including the development of embryos. Turing patterns are a special case of the reaction-diffusion Turing model that create stationary patterns. It is a nonlinear wave that is maintained by a dynamic equlibrium of a system. The wavelength of Turing patterns are determined by the diffusion rates and interactions of the particles without the need for positional information.

One interesting property is that Turing patterns can regenerate even after disturbances are introduced. A limitless variety of spatial patterns can be produced by tuning the morphogen parameters and boundary conditions. A short-range activator coupled with a long-range inhibitor creates stable patterns. Turing patterns have been observed in Zebrafish stripes, The Feather buds of chicks, the hair follicles of mice, and the spots on flower petals.

The two main principles needed for a Turing pattern are local-auto autoactivation and long-range lateral inhibition but more complex patterning can behave like a 2-species model. Notice that with the Activator/Inhibitor model the concentrations are in-phase with each other and for the substrate-depletion model the concentrations are out-f-phase with eachother. This will be important to determine some of the observed dynamics of our system later on.

1D Reaction Diffusion Model

Let’s explore this a little bit more by creating our own model and tweaking some of its parameters to start to build an intuition about these systems. The reaction-diffusion model that we ill use uses 2 morphogens, one which is an inhibitor and one which is an activator.

Our physical domain is a discreteized 1-dimensional field that can be represented by a series of 2-dimensional vectors. We will use a length of 100 and a discretization distance of 1.

The next thing that we need to define is the update equations for our simulation. The first equation is the approximation of the 2nd derivative and its inputs are the 3 spatial values that we need to approximate the 2nd derivative.

The next equation is the equation that describes how the activator is affected by the inhibitor and itself. \(R_a\) is the reaction rate, \(A\) is the activator concentration and \(I\) is the inhibitor concentration.

Next we have the equation for the inhibitor

These can then be used in a typical activator-inhibitor function which includes the diffusion coefficients \(D_A\) and \(D_I\)

Now we can look at the integration methods that our system will use.

In the early experiments we can use Euler’s method but for more advanced simulations we should use a more accurate integration method.

Next we set the parameter values for the simulation, including number of iterations.

Then we populate the initial random fields.

Then we simply step through the simulation time and update our concentrations.

Finally we run the simulation and view the result.

…. And that result is super boring. The activator overpowers the inhibitor almost entirely. If we roll up our previous equations into a single function we can sweep through the parameters to see if we can find anything useful.

After running through different parameters a few times, we can see that most of our combinations don’t work. They reach a sate where either the activator or the inhibitor entirely dominate the space as seen in two examples below.

As I looked through some of the generated results, a few of them had features that jumped out. Our goal was to find parameters for a reaction-diffusion model that mimics patterns that we see in nature. As stated above, the relative concentrations of the morphogens can be used to affect the expression of patterns in tissue such as color or function. Therefore We should be looking for results where the concentration of activators is higher than the inhibitors in some areas, and lower in others.

Below are some examples of some results and my thoughts about what might cause them

Random Results

Don’t know what’s happening here, but I think the parameters might be too small and cause the simulation to be unstable.

This next one is showing the behavior of the activator and inhibitor being out of phase with eachother. This is what we would expect of a substrate-depletion method but not an activator inhibitor method. The only problem is that the distance between the two concentrations is too large, so there would be no expression visible.

These next ones have fewer amplitude changes in the concentrations, but they are closer together. This might be a good canditate to expand the parameters to see if we can get them to interact in a more desirable way.

Now wer’re getting somewhere. This example shows the concentration inversion that we will need to generate patterns. It has large amplitudes and the gradients are close together.

For these next ones we start to see some problems arise. The first and foremost is that we can see an oscillatory behavior in the gradients. This artifact shows the limitations of our simple integration method. There is something about these parameter sets that makes our simulation unstable, therefore we either need to decrease the timestep or use a different integration method to get better results.

All of the first-pass results that were interesting are saved on github Here if you would like to look at them yourself.

Improvements

This project was a good jump into the shallow end of simulating reaction-diffusion models. It was not extensive, and there are many improvements that we may review at a later date.

1) Integration Methods:
Euler’s method is relatively crude and inaccurate. Using other methods such as Runge-Kutta, or Adams-Bashforth, would make for a more accurate simulation. Some adaptive integration methods might allow for larger step sizes, thereby reducing the computation time required. I have a previous example of the classic Runge-Kutta method in this previous article.

2) Boundary Conditions:
I implemented static boundary conditions by not changing the values at the edge of the simulation domain. This is not a very good method for this type of simulation, and an improvement would be to use a more suitable method.

3) Timesteps and Simulation Time:
Another thing that could maybe help is to experiment with the simulation time and the size of the timestep. Without this it is difficult to tell if the results we got were representative of a steady state or if they are just transient snapshots. Adding an additional condition that stops the simulations if the residuals converge to a steady state will be more beneficial

4) Change the Model Equations
From the results that we’re getting it looks like I implemented a substrate/depletion model and not an activator-inhibitor model as the activator and inhibitor are out of phase with eachother. This makes me think that I might need to review my model update equations, but it was useful in illustrating the very basics of implementing these types of systems.

2-Dimensional Models

I extended the model to 2-dimensions to see what types of patterns that wet got and here are some of the results!

What about Spots?

Just as a bonus, we can also generate cheetah-like spots with these types of models. The code for this can be found here.

Sources

Keep In Touch

Twitter Linkedin-in