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; Murray 2003, vol. 18, sec. 2.3). 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 Murray (2003, vol. 18, sec. 2.3) 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

Murray, James D. 2003. Mathematical Biology II: Spatial Models and Biomedical Applications. 3rd ed. Vol. 18. Interdisciplinary Applied Mathematics. New York: Springer. https://doi.org/10.1007/b98869.

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.