In this tutorial, we will look at a simulation of a reaction-diffusion system described by a discretized partial differential equation used as a simple model animal coat pattern formation. The use of reactionâ€“diffusion systems to explain pattern formation in biology may be traced back to Alan Turing's seminal 1952 paper (The Chemical Basis of Morphogenesis).

We will use HoloViews to investigate how the simulated reaction-diffusion system evolves over time for one initial state and then again over varying values of one of the simulation parameters. In the process, it will become clear how easy it is to leverage HoloViews to explore the behaviour of existing analysis or simulation code.

The simulation code for this tutorial has been adapted directly
from the freely available
recipe in the
**IPython Interactive Computing and Visualization Cookbook** by
Cyrille Rossant:

As in the original tutorial, we first need to define a discretization of the Laplacian operator:

In [1]:

```
import numpy as np
def laplacian(Z, dx):
"""
Function to computes the discrete Laplace operator of
a 2D variable on the grid (using a five-point stencil
finite difference method.)
"""
Ztop = Z[0:-2,1:-1]
Zleft = Z[1:-1,0:-2]
Zbottom = Z[2:,1:-1]
Zright = Z[1:-1,2:]
Zcenter = Z[1:-1,1:-1]
return (Ztop + Zleft + Zbottom + Zright - 4 * Zcenter) / dx**2
```

*(time, np.array)* tuples:

In [2]:

```
def reaction_diffusion(a=2.8e-4, b=5e-3, tau=0.1, k=-0.005, samples=10):
"""
We simulate the PDE with the finite difference method.
The samples value is the number of equally spaced samples
to collect over the total simulation time T.
"""
size = 100 # size of the 2D grid
dx = 2./size # space step
T = 10.0 # total time
dt = 4.5 * dx**2 # simulation time step
n = int(T/dt)
result = []
U = np.random.rand(size, size)
V = np.random.rand(size, size)
sample_times = [int(el) for el in np.linspace(0, n, samples)]
for i in range(n):
# We compute the Laplacian of u and v.
deltaU = laplacian(U, dx=dx)
deltaV = laplacian(V, dx=dx)
# We take the values of u and v inside the grid.
Uc = U[1:-1,1:-1]
Vc = V[1:-1,1:-1]
# We update the variables.
U[1:-1,1:-1], V[1:-1,1:-1] = \
Uc + dt * (a * deltaU + Uc - Uc**3 - Vc + k), \
Vc + dt * (b * deltaV + Uc - Vc) / tau
# Neumann conditions: derivatives at the edges
# are null.
for Z in (U, V):
Z[0,:] = Z[1,:]
Z[-1,:] = Z[-2,:]
Z[:,0] = Z[:,1]
Z[:,-1] = Z[:,-2]
if i in sample_times:
result.append((i * dt,U.copy()))
return result
```

The array *U* represents the concentration of some compound
involved in pigmentation, whereas the *V* array represents some
other substance that reacts with the first compound to impede
pigmentation. The partial differential equation for the evolution
of these two compounds over time may be described by these two
coupled equations:

The code above is a fairly direct translation of the original
recipe except we
now collect the simulation results over time as numpy arrays
instead of simply mutating the *U* array until the final state is
reached.

To visualize the results of the simulation using the default
parameter, we first load the ipython extension and import the
`HoloMap`

and
`Image`

classes:

In [3]:

```
import holoviews as hv
hv.notebook_extension()
```

Lets run the first simulation which should complete after a few seconds:

In [4]:

```
sim1 = reaction_diffusion()
```

In [5]:

```
%%opts Image (cmap='copper')
hv.HoloMap({time: hv.Image(array) for (time, array) in sim1}, kdims=['Time'])
```

Out[5]:

*'a'* parameter in addition to the
3-dimensions we explored above:

In [6]:

```
a_values = np.linspace(2.8e-4, 6e-4, 3)
```

*'a'* parameter with each
numpy array:

In [7]:

```
sim2 = [((time, a_value), array) for a_value in a_values
for (time, array) in reaction_diffusion(a=a_value)]
```

In [8]:

```
%%opts Image (cmap='copper')
hv.HoloMap({key: hv.Image(array) for (key,array) in sim2}, kdims=['Time', 'A'])
```

Out[8]: