Turing patterns in RD systems

Turing patterns are what is called by mathematicians a class of finite-wavelength, stationary solutions which can develop from homogeneous initial conditions following local systems of parabolic reaction-diffusion equations. Barbaric terms aside, this is what the common of mortals call "spots" and "stripes" that can be observed on the fur coat of zebras, giraffe, leopards, cats and other animal species (see figure [1]). In this rather short blog post, I intend to numerically implement the classical system of two partial differential equations derived by Alan Turing in his original paper "The Chemical basis of morphogenesis" published in 1952 using Python. The structure is three-folds: In the first part I derive mathematical equations of finite differences schemes used in the code. A second part is dedicated to the code implementation in Python using Numpy and Matplotlib. Finally, the last part is dedicated to what can one obtain by playing with parameters of the model.

Deriving Equations of the scheme

For the sake of simplicity, I choose the easy way by using an Euler-Explicit method for the time discretization while using 1D Finite-Differences for each equation. However, if one wishes, they could use a more creative method such as RK-4 for time as an example. We recall that the system proposed by Turing is the following one $$\begin{aligned} \partial_t U = d_1 \partial_{xx}^2 U + f(U, V) & \quad \text{in } \Omega \\[0.7em] \partial_t V = d_2 \partial_{xx}^2 V + g(U, V) & \quad \text{in } \Omega\\[0.7em] \partial_x U = 0, \partial_x V = 0 & \quad \text{in }\partial\Omega, \end{aligned}$$ Where $U$ and $V$ denote two chemical quantities interacting with each other and $\Omega = [0, 1]$ is the domain in which the problem lies. Conventionally, $U$ is considered to be a chemical activator, which is mathematically translated by $\partial_U f(U, V) > 0$ and $V$ an inhibitor, i.e. $\partial_V g(U, V) < 0$. The two functions $f$ and $g$ are what defines the kinetics of the reaction between $U$ and $V$. They are assumed to be continuous and sufficiently smooth: $f, g \in \mathcal{C}^k(\Omega, \mathbb{R})$ with $k \in \mathbb{N}$. We keep them as general as we can for now, but we will later replace them by actual expressions fitting our modelling needs.

Space-discretization

To be able to run simulations of such model, one has to move to a discrete world, since Python does not understand the notion of continuity. Let us start with space and consider $\Omega_D$ to be a discretized mesh of $\Omega$ defined by the series of points $\{x_0, ..., x_{N-1}\}$ with the convention $x_i = i/(N-1) = ih$, i.e. $(h = 1/(N-1))$. Let $u \in \mathcal{C^2(\Omega, \mathbb{R})}$, we perform classical Taylor expansions to the left and right of an arbitrary point $x_i \in \Omega_D$. $$\begin{aligned} u(x_i + h) = u(x_i) + h u'(x_i) + \frac{h^2}{2} u''(x_i) + O(h^3), \\[1em] u(x_i - h) = u(x_i) - h u'(x_i) + \frac{h^2}{2} u''(x_i) + O(h^3). \end{aligned}$$ Summing both lines makes odd order terms vanish on the right-hand side, leaving us with $$u(x_i + h) + u(x_i - h) = 2u(x_i) + h^2 u''(x_i) + O(h^3).$$ or equivalently $$\frac{u(x_i - h) - 2u(x_i) + u(x_i + h)}{h^2} = u''(x_i) + O(h).$$ We can then drop the "$O(h)$" notation for now and consider that the quantity on the left-hand side is a nice approximation of the second-order derivative of $u$ around $x_i$. In a quest of lighter notations, I will now write $u(x_i) = u_i$ so that we can rewrite both equations of the system at each point of the mesh $\Omega_D$: $$ \begin{aligned} \frac{d}{dt} U_i(t) &= d_1 \frac{U_{i-1}(t) - 2U_i(t) + U_{i+1}(t)}{h^2} + f(U_i(t), V_i(t)) \\[1em] \frac{d}{dt} V_i(t) &= d_2 \frac{V_{i-1}(t) - 2V_i(t) + V_{i+1}(t)}{h^2} + g(U_i(t), V_i(t)). \end{aligned} $$ Still, we have to deal with boundary conditions, but luckily the one-dimensional case is easy to deal with. An approximation of the first order derivative on both endpoints is enough to get discrete boundary conditions $$\partial_x u(x_i) \approx \frac{u(x_i) - u(x_{i+h})}{h}.$$ Plugging in $i=0$ and $i=N-1$ (since $\partial\Omega = \{0, 1\} = \{ x_0, x_{N-1}\}$) yields $$0 = \partial_x U(0, t) = \frac{U_0(t) - U_1(t)}{h}.$$ In other words, $U_0(t) = U_1(t)$ as well as $U_{N-1}(t) = U_{N-2}(t)$. With that, we're all done regarding the space discretization, it is left to deal with time discretization.

Time discretization

hello