Turing Patterns

I’d like here to share the reoprt I prepared for my final project for the course of Numerical Analysis of PDE here at KAUST. All the code to reproduce the image here presented can be found on this dedicated repo. We are interested in solving numerically the PDE system proposed by Alan Turing in 1952 to model the pattern formation on the tissue of some animals such as zebras, ladybugs etc. In particular we will focus our attention on the following diffusion-reaction coupled PDE system:

(1)   \begin{equation*}\begin{cases}\partial_t u = \delta_1(\partial^2_x u + \partial^2_y u)+\alpha u(1-\tau_1v^2)+v(1-\tau_2u)\\\partial_t v = \delta_2(\partial^2_x v + \partial^2_y v)+\beta v(1-\frac{\alpha\tau_1}{\beta}uv)+u(\gamma+\tau_2v)\end{cases}\end{equation*}

and we can separate the diffusion terms from the reactions term, obtaining,

(2)   \begin{equation*}\begin{cases}\partial_t u = \nabla^2\cdot u + f(u,v) \\\partial_t v = \nabla^2\cdot v + g(u,v)\end{cases}\end{equation*}

where the functions f and g are clearly defined as:

(3)   \begin{gather*}f(u,v) = \alpha u(1-\tau_1v^2)+v(1-\tau_2u),\g(u,v) = \beta v(1-\frac{\alpha\tau_1}{\beta}uv)+u(\gamma+\tau_2v).\end{gather*}

To begin we will focus our attention on solving the above coupled system of autonomous ordinary differential equations representing the reaction part of the Turing model. In particular we will solve these implementing the 4th order Runge-Kutta method. First we investigate numerically the stability of the method, to do this we kept fixed the initial conditions u_0 = 1.0, v=0.5, and vary the number of points in the temporal discretization for the values of the parameters \alpha,\beta,\gamma,\delta_1,\delta_2,\tau_1 and \tau_2 in the first line of Parameters Table. In particular from Figure 1 one observe that the method appear to be stable if we chose more then 100 points in the temporal discretization. Now we investigate the long term behaviour of the quantitates u and v for different values in Parameters Table. As we can observe from Figure 2, for different values of initial data and of the parameters all the ODEs in the long term engage in an alternating oscillatory behaviour, similar to the one observed in simpler non linear ODE model such as Lotka-Volterra equations. Studying the linearization of the systems in the equilibrium points, one will found that the Jacobian has only negative eigenvalues. In particular numerically we notice that all the equilibrium points of the system are unstable.

\mathbf{\delta_1}\mathbf{\delta_2}\mathbf{\tau_1}\mathbf{\tau_2}\mathbf{\alpha}\mathbf{\beta}\mathbf{\gamma}
0.002250.00450.020.20.899-0.91-0.899
0.0010.00450.020.20.899-0.91-0.899
0.002250.00450.020.21.9-0.91-1.9
0.002250.00452.020.02.0-2.0-2.0
0.002250.00213.50.00.899-0.91-0.899
0.002250.00450.020.21.9-0.85-1.9
0.002250.00452.020.02.0-0.91-2.0
Parameters Table The table contains the parameters for which we simulate the Turing patterns model.

The next step will be to solve the diffusion equation, in particular we implemented a second order five points finite difference scheme to approximate the second derivative. In particular the matrix:

(4)   \begin{equation*}A = \frac{1}{h^2} \begin{bmatrix}-2 & 1 & & \;\;1\\1 & \ddots & \ddots & \\& \ddots & \ddots & \;\;1\\1 & & 1 & -2\end{bmatrix}\end{equation*}

approximate the one dimensional second derivative with periodic boundary condition, then the approximation of the Laplacian is obtain by the Kronecker product K=\delta_i (A\otimes A), since the same mesh is taken in the horizontal and vertical direction. Now we would like to compute theoretically how small to chose the time marching step in order to obtain a stable integration in time. In order to answer this question we notice that A is a circulant matrix and therefore its eigenvalues can be easily computed,

(5)   \begin{equation*}\lambda_p = \frac{e^{-2\pi i p h}-2+e^{2\pi i p h}}{h^2}=\frac{-2+\cos(2\pi i ph)}{h^2}.\end{equation*}

Furthermore since circulant matrix commute one has that the eigenvalues of the Kronecker product K are given by:

(6)   \begin{equation*}\lambda_{p,q}=\delta_i\frac{-2+2\cos(2\pi i ph)+2\cos(2\pi i qh)}{h^2}\end{equation*}

and this tell us that the eigenvalue of B lie on the real line and attains as maximum value 0 and as minimum value -\delta_i\frac{8}{h^2}. Now for the numerical method to be stable we wont k, the marching step in time, to be such that -\delta_i\frac{8}{h^2}k lies with in the stability region for the 4-th order Runge-Kutta method. Now from Figure 3 we notice that this imply that the following inequality has to be verified,

(7)   \begin{gather*}-\delta_i \frac{8}{h^2}k \geq -z\k \leq \frac{z}{8}\frac{1}{\delta_i} h^2\end{gather*}

where -z is the where the negative real line meets the border of the stability region, marked in red in Figure 3. If we take -z\approx 2.5 we get that for the first value in the Table ?? that we should take a time marching step k\approx 0.027225. The simulation of the diffusive phenomena obtain with this numerical scheme can be observed in Figure 4.

Now is time to simulate both the reaction and diffusion term together, the result of this simulation can be seen in Figure 5. Using the 4th order Runge-Kutta method for t=100 the running time is approximatively 3 minutes.

To over come the CFL condition and have a faster method we resort to a splitting technique, i.e. we solve first the diffusion part of the PDE using an implicit Euler method. The implicit Euler scheme is absolutely stable therefore we can use much larger time steps. Then we use the solution obtain as initial data to the reaction part of the system and this is solved using the 4-th order Runge-Kutta method previously implemented. This is advantageous because as we have seen in the beginning the time step required to solve in a stable manner \eqref{eq:ReactionODE} is much larger then the one needed to solve in a stable manner the diffusive part of the model. The result of this simulation at t=150, obtained using a time marching step k=0.3 are shown in Figure 6. From Figure 6 we observe that even starting from initial data without any sort of pattern such as the one shown in Figure 4 one obtains solutions with patterns such as the one presented in different animals tissue.

Here will follow the code to obtain all figures presented in this report.

import numpy as np
from math import sqrt
import scipy.integrate
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from scipy import sparse
from scipy.sparse.linalg import eigs
from scikits.umfpack import spsolve #UMFPACK If you work in serial.
from tqdm import trange, tqdm
from nodepy.runge_kutta_method import *

class ODESol:
    def __init__(self,timesteps,timestep,U):
        self.t = timesteps;
        self.h = timestep;
        self.y = U;
def RK4(F,T,U0,arg,N):
    tt = np.linspace(T[0],T[1],N);
    h = tt[1]-tt[0];
    U = np.zeros([len(U0),N]);
    U[:,0] = U0;
    for i in trange(0,N-1):
        Y1 = U[:,i];
        Y2 = U[:,i] + 0.5*h*F(tt[i],Y1,arg);
        Y3 = U[:,i] + 0.5*h*F(tt[i]+0.5*h,Y2,arg);
        Y4 = U[:,i] + h*F(tt[i]+0.5*h,Y3,arg);
        U[:,i+1] = U[:,i]+(h/6)*(F(tt[i],Y1,arg)+2*F(tt[i]+ 0.5*h,Y2,arg)+2*F(tt[i]+ 0.5*h,Y3,arg)+F(tt[i]+h,Y4,arg))
    sol = ODESol(tt,h,U);
    return sol;

Table = []
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02,"tau2":0.2,"alpha":0.899, "beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.001,"delta2":0.0045,"tau1":0.02,"tau2":0.2,"alpha":0.899, "beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02,"tau2":0.2,"alpha":1.9, "beta":-0.91,"gamma":-1.9}]

Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":2.02,"tau2":0.0,"alpha":2.0, "beta":-0.91,"gamma":-2}]
Table = Table + [{"delta1":0.00105,"delta2":0.0021,"tau1":3.5,"tau2":0.0,"alpha":0.899, "beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02,"tau2":0.2,"alpha":1.9, "beta":-0.85,"gamma":-1.9}]
Table = Table + [{"delta1":0.00225,"delta2":0.0005,"tau1":2.02,"tau2":0.0,"alpha":2.0, "beta":-0.91,"gamma":-2}]

def Reaction(t,x,parameters):
    #The ODE is autonomus so we don't really need
    #the depende on time.
    u = x[0]; v = x[1]; #We grab the useful quantity.
    d1 = parameters["delta1"]; d2 = parameters["delta2"];
    t1 = parameters["tau1"]; t2 = parameters["tau2"];
    a = parameters["alpha"]; b = parameters["beta"]; g = parameters["gamma"];
    du = a*u*(1-t1*(v**2))+v*(1-t2*u);
    dv = b*v*(1+(a*t1/b)*(u*v))+u*(g+t2*v);
    return np.array([du,dv])

t0 = 0.                  # Initial time
u0 = np.array([0.5,0.5])# Initial values
tfinal = 100.              # Final time
dt_output=0.1# Interval between output for plotting
N=int(tfinal/dt_output)       # Number of output times
print(N)
tt=np.linspace(t0,tfinal,N)  # Output times
ODE = RK4(Reaction,[t0,tfinal],u0,Table[6],N);
uu=ODE.y
plt.plot(tt,uu[0,:],tt,uu[1,:])
plt.legend(["u","v"])
plt.show();
import numpy as np
from math import sqrt
import scipy.integrate
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from scipy import sparse
from scipy.sparse.linalg import eigs
from scikits.umfpack import spsolve #UMFPACK If you work in serial.
from tqdm import trange, tqdm
from nodepy.runge_kutta_method import *

class ODESol:
    def __init__(self,timesteps,timestep,U):
        self.t = timesteps;
        self.h = timestep;
        self.y = U;
def RK4(F,T,U0,arg,N):
    tt = np.linspace(T[0],T[1],N);
    h = tt[1]-tt[0];
    U = np.zeros([len(U0),N]);
    U[:,0] = U0;
    for i in trange(0,N-1):
        Y1 = U[:,i];
        Y2 = U[:,i] + 0.5*h*F(tt[i],Y1,arg);
        Y3 = U[:,i] + 0.5*h*F(tt[i]+0.5*h,Y2,arg);
        Y4 = U[:,i] + h*F(tt[i]+0.5*h,Y3,arg);
        U[:,i+1] = U[:,i]+(h/6)*(F(tt[i],Y1,arg)+2*F(tt[i]+ 0.5*h,Y2,arg)+2*F(tt[i]+ 0.5*h,Y3,arg)+F(tt[i]+h,Y4,arg))
    sol = ODESol(tt,h,U);
    return sol;

Table = []
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":0.899,"beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.001,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":0.899,"beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":1.9,"beta":-0.91,"gamma":-1.9}]

Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":2.02, "tau2":0.0,"alpha":2.0,"beta":-0.91,"gamma":-2}]
Table = Table + [{"delta1":0.00105,"delta2":0.0021,"tau1":3.5, "tau2":0.0,"alpha":0.899,"beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":1.9,"beta":-0.85,"gamma":-1.9}]
Table = Table + [{"delta1":0.00225,"delta2":0.0005,"tau1":2.02, "tau2":0.0,"alpha":2.0,"beta":-0.91,"gamma":-2}]

def laplacian_1D(m):
    em = np.ones(m)
    e1=np.ones(m-1)
    A = (sparse.diags(-2*em,0)+sparse.diags(e1,-1)+sparse.diags(e1,1))/((2/(m+1))**2);
    A[0,-1]=1/((2/(m+1))**2);
    A[-1,0]=1/((2/(m+1))**2);
    return A;
def laplacian_2D(m):
    I = np.eye(m)
    A = laplacian_1D(m)
    return sparse.kron(A,I) + sparse.kron(I,A)

def Diffusion(t,x,parameters):
    B = sparse.bmat([[parameters["delta1"]*A,None],[None,parameters["delta2"]*A]]);
    return B*x;
m=100
x=np.linspace(-1,1,m+2); x=x[1:-1]
y=np.linspace(-1,1,m+2); y=y[1:-1]
X,Y=np.meshgrid(x,y)
A=laplacian_2D(m)
#Generating the initial data
mu, sigma = 0, 0.5 # mean and standard deviation
u0 = np.random.normal(mu, sigma,m**2);
v0 = np.random.normal(mu, sigma,m**2)
#Plotting the initial data
plt.figure()
U0=u0.reshape([m,m])
plt.pcolor(X,Y,U0)
plt.colorbar();
plt.figure()
V0=v0.reshape([m,m])
plt.pcolor(X,Y,V0)
plt.colorbar();
TIndex = 0;
t0 = 0.0                  # Initial time
tfinal = 100            # Final time
dt_output=0.018          # Interval between output for plotting
N=int(tfinal/dt_output)       # Number of output times
print("CFL suggested time step",(2.5*(0.0198)**2)/(4*(Table[TIndex]["delta1"]+Table[TIndex]["delta2"])));
#ODE = scipy.integrate.solve_ivp(Diffusion,[t0,tfinal],np.append(u0,v0,axis=0),args=[Table[TIndex]],method='RK45',t_eval=tt,atol=1.e-10,rtol=1.e-10);
ODE = RK4(Diffusion,[t0,tfinal],np.append(u0,v0,axis=0),Table[TIndex],N);
uu=ODE.y
print("Step size of the Method",ODE.h);
print(ODE.h*eigs(sparse.bmat([[Table[TIndex]["delta1"]*A,None],[None,Table[TIndex]["delta2"]*A]]),k=10)[0])
RK44 = loadRKM('RK44')
RK44.plot_stability_region(bounds=[-5,1,-5,5])
plt.plot(ODE.h*eigs(sparse.bmat([[Table[TIndex]["delta1"]*A,None],[None,Table[TIndex]["delta2"]*A]]),k=10)[0].real,ODE.h*eigs(sparse.bmat([[Table[TIndex]["delta1"]*A,None],[None,Table[TIndex]["delta2"]*A]]),k=10)[0].imag,"b*")
ut = uu[0:m**2,-1];
vt = uu[m**2:,-1];
Ut=ut.reshape([m,m])
plt.figure()
plt.pcolor(X,Y,Ut)
plt.colorbar();
plt.figure()
Vt=vt.reshape([m,m])
plt.pcolor(X,Y,Vt)
plt.colorbar();
plt.show();
import numpy as np
from math import sqrt
import scipy.integrate
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from scipy import sparse
from scipy.sparse.linalg import eigs
from scikits.umfpack import spsolve #UMFPACK If you work in serial.
from tqdm import trange, tqdm
from nodepy.runge_kutta_method import *

class ODESol:
    def __init__(self,timesteps,timestep,U):
        self.t = timesteps;
        self.h = timestep;
        self.y = U;
def laplacian_1D(m):
    em = np.ones(m)
    e1=np.ones(m-1)
    A = (sparse.diags(-2*em,0)+sparse.diags(e1,-1)+sparse.diags(e1,1))/((2/(m+1))**2);
    A[0,-1]=1/((2/(m+1))**2);
    A[-1,0]=1/((2/(m+1))**2);
    return A;
def laplacian_2D(m):
    I = np.eye(m)
    A = laplacian_1D(m)
    return sparse.kron(A,I) + sparse.kron(I,A)
def Reaction(t,x,parameters):
    u = x[0:m**2]; v = x[m**2:]; #We grab the useful quantity.
    d1 = parameters["delta1"]; d2 = parameters["delta2"];
    t1 = parameters["tau1"]; t2 = parameters["tau2"];
    a = parameters["alpha"]; b = parameters["beta"]; g = parameters["gamma"];
    #Reaction
    du = a*u*(1-t1*(v**2))+v*(1-t2*u);
    dv = b*v*(1+(a*t1/b)*(u*v))+u*(g+t2*v);
    b = np.append(du,dv,axis=0);
    return b;
def SplitSolver(F,T,U0,arg,N):
    tt = np.linspace(T[0],T[1],N);
    h = tt[1]-tt[0];
    U = np.zeros([len(U0),N]);
    U[:,0] = U0;
    for i in trange(0,N-1):
        B = sparse.bmat([[arg["delta1"]*A,None],[None,arg["delta2"]*A]]);
        UStar = spsolve((sparse.identity(B.shape[0])-h*B),U[:,i])
        Y1 = UStar;
        Y2 = UStar + 0.5*h*F(tt[i],Y1,arg);
        Y3 = UStar + 0.5*h*F(tt[i]+0.5*h,Y2,arg);
        Y4 = UStar + h*F(tt[i]+0.5*h,Y3,arg);
        U[:,i+1] = UStar+(h/6)*(F(tt[i],Y1,arg)+2*F(tt[i]+ 0.5*h,Y2,arg)+2*F(tt[i]+ 0.5*h,Y3,arg)+F(tt[i]+h,Y4,arg))
    sol = ODESol(tt,h,U);
    return sol;
Table = []
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":0.899,"beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.001,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":0.899,"beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":1.9,"beta":-0.91,"gamma":-1.9}]

Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":2.02, "tau2":0.0,"alpha":2.0,"beta":-0.91,"gamma":-2}]
Table = Table + [{"delta1":0.00105,"delta2":0.0021,"tau1":3.5, "tau2":0.0,"alpha":0.899,"beta":-0.91,"gamma":-0.899}]
Table = Table + [{"delta1":0.00225,"delta2":0.0045,"tau1":0.02, "tau2":0.2,"alpha":1.9,"beta":-0.85,"gamma":-1.9}]
Table = Table + [{"delta1":0.00225,"delta2":0.0005,"tau1":2.02, "tau2":0.0,"alpha":2.0,"beta":-0.91,"gamma":-2}]

TIndex = int(input("What pattern you want to simulate ? "))
m=100
x=np.linspace(0,1,m+2); x=x[1:-1]
y=np.linspace(0,1,m+2); y=y[1:-1]
X,Y=np.meshgrid(x,y)
A=laplacian_2D(m)
#Generating the initial data
mu, sigma = 0, 0.5 # mean and standard deviation
u0 = np.random.normal(mu, sigma,m**2);
v0 = np.random.normal(mu, sigma,m**2)
#Plotting the initial data
plt.figure()
U0=u0.reshape([m,m])
plt.pcolor(X,Y,U0)
plt.colorbar();
plt.title("Initial Data u");
plt.figure()
V0=v0.reshape([m,m])
plt.pcolor(X,Y,V0)
plt.colorbar();
plt.title("Initial Data v");

t0 = 0.0                  # Initial time
tfinal = 150            # Final time
dt_output=0.3          # Interval between output for plotting
N=int(tfinal/dt_output)       # Number of output times
ODE = SplitSolver(Reaction,[t0,tfinal],np.append(u0,v0,axis=0),Table[TIndex],N);
uu=ODE.y
ut = uu[0:m**2,-1];
vt = uu[m**2:,-1];
Ut=ut.reshape([m,m])
plt.figure()
plt.pcolor(X,Y,Ut)
plt.colorbar();
plt.figure()
Vt=vt.reshape([m,m])
plt.pcolor(X,Y,Vt)
plt.colorbar();
plt.show();

Leave a Reply

Your email address will not be published. Required fields are marked *