Turing Numerics

Module 4: Space and Pattern Formation

This is the prediction step of the laboratory. You turn the theory of Turing Conditions and Dispersion Relation into three small functions: equilibrium, Jacobian, and Turing test (Turing 1952; murray2003mathematical?). The fourth, the mode scan, follows on Unstable Modes in 1D. You will reuse all four in 2D and in the assignment.

The workflow is always the same:

Find the homogeneous equilibrium.
Compute the Jacobian at that equilibrium.
Compute the Laplacian eigenvalues from the domain and the boundary conditions.
Build the mode matrix $A_n = \gamma J - k_n^2 D$ for each mode.
Compute the eigenvalues of each $2 \times 2$ matrix and keep the modes whose largest real part is positive.

Compute the Equilibrium

The steady state from Linearization is one line of code each.

def gierer_meinhardt_equilibrium(a=0.40, b=1.00):
    u_star = (a + 1) / b
    v_star = u_star**2
    return u_star, v_star

Write the Jacobian at the Fixed Point

Exercise

Use the general Jacobian from Linearization and evaluate it at $(u_*, v_*)$.

def gierer_meinhardt_jacobian(a=0.40, b=1.00):
    u_star, v_star = gierer_meinhardt_equilibrium(a, b)
    # TODO: return fu, fv, gu, gv

Solution

def gierer_meinhardt_jacobian(a=0.40, b=1.00):
    u_star, v_star = gierer_meinhardt_equilibrium(a, b)
    fu = -b + 2 * u_star / v_star
    fv = -(u_star**2) / (v_star**2)
    gu = 2 * u_star
    gv = -1.0
    return fu, fv, gu, gv

Simplified, this is $f_u = b(1-a)/(1+a)$, $f_v = -b^2/(1+a)^2$, $g_u = 2(1+a)/b$, $g_v = -1$. You can also write the entries directly in $a$ and $b$:

def gierer_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

Implement the Turing Test

Exercise

Translate the four inequalities of (murray2003mathematical?) into Python. Use & instead of and so the function works on arrays.

import numpy as np

def is_turing_instability(a, b, d):
    fu, fv, gu, gv = gierer_meinhardt_jacobian(a, b)
    det = fu * gv - fv * gu

    # TODO: write the four conditions
    cond1 =
    cond2 =
    cond3 =
    cond4 =
    return cond1 & cond2 & cond3 & cond4

Solution

def is_turing_instability(a, b, d):
    fu, fv, gu, gv = gierer_meinhardt_jacobian(a, b)
    det = fu * gv - fv * gu
    cond1 = (fu + gv) < 0
    cond2 = det > 0
    cond3 = (d * fu + gv) > 0
    cond4 = (d * fu + gv) ** 2 > 4 * d * det
    return cond1 & cond2 & cond3 & cond4

The course package writes the last two conditions as the single inequality $d f_u + g_v > 2\sqrt{d \det J}$, which is equivalent.

If your implementation is vectorized, it works with scalars and with arrays. Validate it carefully: you reuse it in 2D and in the assignment.

Plot the Turing Space

Exercise

Complete the code to recreate the diagram on Turing Space. Mark the points $(a, d) = (0.4, 20)$ and $(0.4, 30)$. Which one lies inside the instability region?

import matplotlib.pyplot as plt
import numpy as np

arr_a = np.linspace(0, 1, 500)
arr_d = np.linspace(0, 100, 500)
mesh_a, mesh_d = np.meshgrid(arr_a, arr_d)
mask_turing =  # TODO: evaluate the Turing test on the mesh, with b = 1

fig, ax = plt.subplots(figsize=(7, 4))
ax.contourf(mesh_a, mesh_d, mask_turing)
# TODO: mark the two test points
ax.set_xlabel("a")
ax.set_ylabel("d")
ax.set_title("Turing Space")
plt.show()

Solution

Pass the whole mesh to the vectorized test. The point $(0.4, 30)$ lies inside the region and $(0.4, 20)$ lies outside, because the boundary at $a = 0.4$ sits at $d_c \approx 26.24$.

mask_turing = is_turing_instability(mesh_a, 1.0, mesh_d)

fig, ax = plt.subplots(figsize=(7, 4))
ax.contourf(mesh_a, mesh_d, mask_turing)
ax.plot([0.4], [30], "o", color="tab:red", label="Case A: d = 30")
ax.plot([0.4], [20], "s", color="tab:gray", label="Case B: d = 20")
ax.set_xlabel("a")
ax.set_ylabel("d")
ax.set_title("Turing Space")
ax.legend()
plt.show()

The same diagram appears in the interactive 1D Gierer-Meinhardt simulator, so the algebra and the simulation refer to the same parameter plane.

References

Turing, Alan M. 1952. “The Chemical Basis of Morphogenesis.” Philosophical Transactions of the Royal Society B 237 (641): 37–72. https://doi.org/10.1098/rstb.1952.0012.

--- title: "Turing Numerics" subtitle: "Module 4: Space and Pattern Formation" format: html --- This is the prediction step of the laboratory. You turn the theory of [Turing Conditions](turing-conditions.qmd) and [Dispersion Relation](dispersion-relation.qmd) into three small functions: equilibrium, Jacobian, and Turing test [@turing1952chemical; @murray2003mathematical, Section 2.3]. The fourth, the mode scan, follows on [Unstable Modes in 1D](unstable-modes-1d.qmd). You will reuse all four in 2D and in the assignment. The workflow is always the same: - Find the homogeneous equilibrium. - Compute the Jacobian at that equilibrium. - Compute the Laplacian eigenvalues from the domain and the boundary conditions. - Build the mode matrix $A_n = \gamma J - k_n^2 D$ for each mode. - Compute the eigenvalues of each $2 \times 2$ matrix and keep the modes whose largest real part is positive. ## Compute the Equilibrium The steady state from [Linearization](linearization.qmd) is one line of code each. ```python def gierer_meinhardt_equilibrium(a=0.40, b=1.00): u_star = (a + 1) / b v_star = u_star**2 return u_star, v_star ``` ## Write the Jacobian at the Fixed Point ::: {.callout-note} ## Exercise Use the general Jacobian from [Linearization](linearization.qmd) and evaluate it at $(u_*, v_*)$. ```python def gierer_meinhardt_jacobian(a=0.40, b=1.00): u_star, v_star = gierer_meinhardt_equilibrium(a, b) # TODO: return fu, fv, gu, gv ``` ::: ::: {.callout-tip collapse="true"} ## Solution ```python def gierer_meinhardt_jacobian(a=0.40, b=1.00): u_star, v_star = gierer_meinhardt_equilibrium(a, b) fu = -b + 2 * u_star / v_star fv = -(u_star**2) / (v_star**2) gu = 2 * u_star gv = -1.0 return fu, fv, gu, gv ``` Simplified, this is $f_u = b(1-a)/(1+a)$, $f_v = -b^2/(1+a)^2$, $g_u = 2(1+a)/b$, $g_v = -1$. You can also write the entries directly in $a$ and $b$: ```python def gierer_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 ``` ::: ## Implement the Turing Test ::: {.callout-note} ## Exercise Translate the four inequalities of @murray2003mathematical [Section 2.3] into Python. Use `&` instead of `and` so the function works on arrays. ```python import numpy as np def is_turing_instability(a, b, d): fu, fv, gu, gv = gierer_meinhardt_jacobian(a, b) det = fu * gv - fv * gu # TODO: write the four conditions cond1 = cond2 = cond3 = cond4 = return cond1 & cond2 & cond3 & cond4 ``` ::: ::: {.callout-tip collapse="true"} ## Solution ```python def is_turing_instability(a, b, d): fu, fv, gu, gv = gierer_meinhardt_jacobian(a, b) det = fu * gv - fv * gu cond1 = (fu + gv) < 0 cond2 = det > 0 cond3 = (d * fu + gv) > 0 cond4 = (d * fu + gv) ** 2 > 4 * d * det return cond1 & cond2 & cond3 & cond4 ``` The course package writes the last two conditions as the single inequality $d f_u + g_v > 2\sqrt{d \det J}$, which is equivalent. ::: If your implementation is vectorized, it works with scalars and with arrays. Validate it carefully: you reuse it in 2D and in the assignment. ## Plot the Turing Space ::: {.callout-note} ## Exercise Complete the code to recreate the diagram on [Turing Space](turing-space.qmd). Mark the points $(a, d) = (0.4, 20)$ and $(0.4, 30)$. Which one lies inside the instability region? ```python import matplotlib.pyplot as plt import numpy as np arr_a = np.linspace(0, 1, 500) arr_d = np.linspace(0, 100, 500) mesh_a, mesh_d = np.meshgrid(arr_a, arr_d) mask_turing = # TODO: evaluate the Turing test on the mesh, with b = 1 fig, ax = plt.subplots(figsize=(7, 4)) ax.contourf(mesh_a, mesh_d, mask_turing) # TODO: mark the two test points ax.set_xlabel("a") ax.set_ylabel("d") ax.set_title("Turing Space") plt.show() ``` ::: ::: {.callout-tip collapse="true"} ## Solution Pass the whole mesh to the vectorized test. The point $(0.4, 30)$ lies inside the region and $(0.4, 20)$ lies outside, because the boundary at $a = 0.4$ sits at $d_c \approx 26.24$. ```python mask_turing = is_turing_instability(mesh_a, 1.0, mesh_d) fig, ax = plt.subplots(figsize=(7, 4)) ax.contourf(mesh_a, mesh_d, mask_turing) ax.plot([0.4], [30], "o", color="tab:red", label="Case A: d = 30") ax.plot([0.4], [20], "s", color="tab:gray", label="Case B: d = 20") ax.set_xlabel("a") ax.set_ylabel("d") ax.set_title("Turing Space") ax.legend() plt.show() ``` ::: The same diagram appears in the [interactive 1D Gierer-Meinhardt simulator](https://bam-applied-math-lab.streamlit.app/Gierer_Meinhardt_1D), so the algebra and the simulation refer to the same parameter plane.