Turing Instability

Pattern Formation in Reaction-Diffusion

Turing instability explains when a uniform steady state becomes unstable to spatial perturbations, leading to patterns. We will implement the checks and visualize the instability region.

Figure 1: Turing space for the Gierer-Meinhardt model (a vs d).

Conditions

Let the Jacobian at the fixed point be

\[J=\begin{pmatrix} f_u & f_v \\ g_u & g_v \end{pmatrix}. \tag{1}\]

A common set of conditions for Turing instability is:

\[ \begin{aligned} &f_u + g_v < 0, \\ &f_u g_v - f_v g_u > 0, \\ &g_v + d f_u > 2\sqrt{d\,(f_u g_v - f_v g_u)}. \end{aligned} \tag{2}\]

Implement a Turing test

Create a function that returns True if all conditions are satisfied. In the reference script the Jacobian entries are hard-coded at the fixed point.

def giere_meinhardt_jacobian(a=0.40, b=1.00):
    # TODO: compute fu, fv, gu, gv at the fixed point
    return fu, fv, gu, gv


def is_turing_instability(a, b, d):
    fu, fv, gu, gv = giere_meinhardt_jacobian(a, b)
    # TODO: implement the three conditions
    return cond1 & cond2 & cond3

Hint: Turing instability functions (click to expand)

def giere_meinhardt_jacobian(a=0.40, b=1.00):
    fu = 2 * b / (a + 1) - b
    fv = -((b / (a + 1)) ** 2)
    gu = 2 * (a + 1) / b
    gv = -1.0
    return fu, fv, gu, gv


def is_turing_instability(a, b, d):
    fu, fv, gu, gv = giere_meinhardt_jacobian(a, b)
    nabla = fu * gv - fv * gu
    cond1 = (fu + gv) < 0
    cond2 = nabla > 0
    cond3 = (gv + d * fu) > (2 * np.sqrt(d) * np.sqrt(nabla))
    return cond1 & cond2 & cond3

Plot the Turing space

Plot the region of \((a, d)\) where the instability holds.

import numpy as np
import matplotlib.pyplot as plt

arr_a = np.linspace(0, 1, 1000)
arr_d = np.linspace(0, 100, 1000)
mesh_a, mesh_d = np.meshgrid(arr_a, arr_d)
mask_turing = is_turing_instability(mesh_a, b, mesh_d)

fig, ax = plt.subplots(figsize=(6, 4))
ax.contourf(mesh_a, mesh_d, mask_turing)
ax.set_xlabel("a")
ax.set_ylabel("d")
ax.set_title("Turing Space")

Ensure your function works with arrays.

Unstable spatial modes

The temporal eigenvalues for spatial mode \(\lambda_n\) are the eigenvalues of

\[A_n = J - \lambda_n \;\mathrm{diag}(1, d). \tag{3}\]

You can scan modes to find which ones are unstable:

import numpy as np

def find_unstable_spatial_modes(a=0.40, b=1.00, d=30.0, length=40.0, num_modes=10):
    # TODO: build Jacobian and scan spatial modes
    return unstable_modes

Hint: Spatial mode scan (click to expand)

def find_unstable_spatial_modes(a=0.40, b=1.00, d=30.0, length=40.0, num_modes=10):
    fu, fv, gu, gv = giere_meinhardt_jacobian(a, b)
    jac = np.array([[fu, fv], [gu, gv]])
    n_values = np.arange(num_modes)
    max_eigs = np.zeros(num_modes)

    for n in n_values:
        lambda_n = (n * np.pi / length) ** 2
        a_n = jac - lambda_n * np.diag([1, d])
        sigma1, sigma2 = np.linalg.eigvals(a_n)
        max_eigs[n] = max(sigma1.real, sigma2.real)

    sorted_indices = np.argsort(max_eigs)[::-1]
    unstable_modes = sorted_indices[max_eigs[sorted_indices] > 0]
    return unstable_modes

Interaction extension (optional)

Make the plot interactive: clicking on the Turing diagram updates \((a, d)\) and restarts the animation of \(v(x)\).

from matplotlib.backend_bases import MouseEvent

def mouse_click(event: MouseEvent):
    if event.inaxes == ax_turing:
        a = event.xdata
        d = event.ydata
        # TODO: reset u,v and restart the animation

fig.canvas.mpl_connect("button_press_event", mouse_click)

Tip

You can find a complete interactive implementation in amlab/pdes_1d.