Calculate which ripples will grow

Choose reactions we can follow completely

Let $u,v$ be dimensionless concentrations, with reference units also chosen for time and length. Use $f=a-u+u^2v$ and $g=b-u^2v$ , with constant supplies $a,b>0$ . The conversion term $u^2v$ is autocatalytic: more of the first substance accelerates its local production.

A uniform equilibrium has no diffusion flux. Adding the two reaction equations cancels conversion and gives $u_*=a+b$ . Substitution into the second equation gives $v_*=b/(a+b)^2$ . Supply and consumption continue even though the concentrations remain constant.

f(u_*,v_*)=g(u_*,v_*)=0,\qquad (u_*,v_*)=\left(a+b,\frac{b}{(a+b)^2}\right)

A star denotes the uniform equilibrium. This is an open reaction system; supply and consumption change the sum of concentrations.

Reading: J. Schnakenberg

Organize the response to small changes

Write $u=u_*+\delta u$ and $v=v_*+\delta v$ . Keeping only terms linear in the disturbance replaces each reaction curve by its local tangent. The four partial derivatives form the Jacobian. Each column asks how both reaction rates respond to a small increase of one concentration.

For a disturbance with fixed component ratio and time dependence $e^{\lambda t}$ , substitution gives an eigenvalue problem. Negative $\operatorname{Re}\lambda$ means decay; positive values mean growth. An imaginary part describes oscillation. The trace is the sum of the diagonal entries; the determinant is their product minus the product of the other two entries. Negative trace and positive determinant make both eigenvalue real parts negative.

J=\begin{pmatrix}-1+2b/(a+b)&(a+b)^2\\-2b/(a+b)&-(a+b)^2\end{pmatrix}

Differentiate each reaction with respect to each concentration, then evaluate at equilibrium.

\lambda^2-\operatorname{tr}(J)\lambda+\det J=0,\qquad \operatorname{tr}(J)<0,\quad\det J>0

The inequalities ensure stability of the uniform reaction system. Equality requires further analysis.

Diffusion changes the same matrix

Consider a spatial component $e^{i\mathbf q\cdot\mathbf x}$ . Its wave number is $q=|\mathbf q|$ and its wavelength is $2\pi/q$ . The Laplacian multiplies this component by $-q^2$ , so diffusion subtracts $D_uq^2$ and $D_vq^2$ from the two diagonal entries.

Both changes are negative. Yet the reacting variables are coupled, so their diagonal terms alone do not determine stability. Unequal diffusion changes the relative response of the two components and can give their coupled system a positive eigenvalue.

Write $J=(J_{ij})$ . If the reactions are stable, diffusion makes the trace more negative. Instability can then occur when the determinant becomes negative, giving one positive and one negative real eigenvalue. Expanding the determinant produces a quadratic in $q^2$ .

A(q)=J-q^2\begin{pmatrix}D_u&0\\0&D_v\end{pmatrix},\qquad \lambda_+(q)=\max\operatorname{Re}\operatorname{eig}A(q)

The growth spectrum associates each wave number with its fastest linear growth rate. Boundaries determine which wave numbers are available.

\det A(q)=\det J-(D_vJ_{11}+D_uJ_{22})q^2+D_uD_vq^4

A sufficiently large middle coefficient lets the determinant become negative over a finite interval.

\operatorname{Re}\lambda_+(q)

Analytic reference · dimensionless

Shaded interval: positive growth

q

Original calculation for the worked-example parameters. This is a continuum linear prediction, separate from the experiment’s recorded growth points.

Collect the Turing conditions

Define $B=D_vJ_{11}+D_uJ_{22}$ . To cross zero at positive $q^2$ , the quadratic needs $B>0$ and a positive discriminant. Combining these with local stability gives the Turing conditions for two species with constant diagonal diffusion.

If $D_u=D_v=D$ , every reaction eigenvalue simply loses $Dq^2$ . Equal diffusion cannot destabilize a locally stable equilibrium. This provides a useful independent limit.

A continuous unstable interval does not guarantee that a small finite box contains an allowed mode in it. The next chapter introduces boundaries and grids, then estimates growth from actual simulation records.

\operatorname{tr}J<0,\quad\det J>0,\quad B>0,\quad B^2>4D_uD_v\det J

These conditions predict an unstable nonzero wave number in continuous space. The full nonlinear dynamics determine the final pattern.

Reading: A. M. Turing

Worked example

Parameters that select an intermediate scale

Choose $a=0.1,b=0.9,D_u=0.05,D_v=1$ . The equilibrium is $(1,0.9)$ , with Jacobian $\left(\begin{smallmatrix}0.8&1\\-1.8&-1\end{smallmatrix}\right)$ , trace $-0.2$ and determinant $1$ . The reaction system is stable.

Now $B=0.75$ and the discriminant is $0.3625>0$ . The determinant is negative for approximately $1.216<q<3.677$ . At $q=2$ , the mode matrix has trace $-4.4$ and determinant $-1.2$ , giving a larger eigenvalue of approximately $0.258$ .

A sufficiently small pure growing mode therefore multiplies its amplitude by $e^{0.258}\approx1.29$ in one time unit. Once the disturbance changes the local reaction slopes appreciably, the full nonlinear equations must replace this linear prediction.

\det A(q)=1-0.75q^2+0.05q^4

The specified parameters reduce the stability calculation to this explicit polynomial.

Try the reasoning

Predict before inspecting the spectrum

Use the worked example to check uniform stability and predict the unstable wave-number interval.
Compare the theoretical spectrum with points estimated from the early simulation record. Distinguish continuous-space and discrete-grid references.
Compare the stable and oscillatory presets. Inspect the zero-wave-number behavior before explaining how each differs from the Turing case.

Exercises and explanations

Why does $J_{11}>0$ alone not prove uniform instability?
Read the explanation
The second concentration responds and feeds back onto the first. Eigenvalues describe their combined response. In the worked example the first diagonal entry is positive, yet both uniform eigenvalues have negative real parts.
What happens if both diffusion coefficients become $0.05$ ?
Read the explanation
Each reaction eigenvalue loses $0.05q^2$ . Its real part becomes more negative, so no diffusion-driven instability occurs.
Does the fastest linear mode uniquely determine whether the final pattern contains spots or stripes?
Read the explanation
No. Linear analysis predicts early growth and a preferred scale. Nonlinear mode interactions, boundaries and initial conditions affect the saturated pattern. The nonlinear equations must be evolved to study that state.

Course sources and calculation scope

The derivations, diagrams and numerical examples are authored for this course. Source locations below distinguish adopted equations, teaching extensions and the one exact published example recalculated here. Interactive fields and noisy observations are model-generated.

A. M. Turing, The Chemical Basis of Morphogenesis (1952), Philosophical Transactions B 237, 37–72
Read §4, pp. 42–43, for the two-cell numerical example and §6 for mode analysis. The reproduced result is early linear growth in §4. No source figure is copied; the Schnakenberg experiments use a different model.
J. Schnakenberg, Simple Chemical Reaction Systems with Limit Cycle Behaviour (1979), Journal of Theoretical Biology 81, 389–400
Original reaction-model reference. The dimensionless irreversible teaching form is specified here; its equilibrium, Jacobian and numerical examples are derived in the course. Parameters are not biological measurements.
D. Barkley, A Model for Fast Computer Simulation of Waves in Excitable Media (1991), Physica D 49, 61–70
Read the two-variable equations in §2 and the discussion of applicability. The experiment solves those continuous equations with its own discretization and teaching initial conditions; source images and rotation frequencies are not reproduced.
E. J. Crampin, E. A. Gaffney and P. K. Maini, Reaction and Diffusion on Growing Domains (1999), Bulletin of Mathematical Biology 61, 1093–1120
Read §3, equations 9–20, on conservation, transport and growth. The paper mainly studies one-dimensional growth. The two-dimensional uniformly expanding square is a separately derived teaching extension, not a reproduction of its mode-doubling figures.
M. Nishikawa, S. R. Naganathan, F. Jülicher and S. W. Grill, Controlling Contractile Instabilities in the Actomyosin Cortex (2017), eLife 6, e19595
Read Appendix equations 4–9 on density flux, active viscous stress, saturation and force balance. This teaching reduction uses one density, linear turnover and local stress inhibition. It does not fit embryo data or reproduce Figure 2. The source is CC BY.
A. L. Krause, M. A. Ellis and R. A. Van Gorder, Influence of Curvature, Growth, and Anisotropy on the Evolution of Turing Patterns on Growing Manifolds (online 2018; issue 2019)
Read equations 4–6 for moving-domain conservation, 14–15 for material coordinates and 44–45 for Schnakenberg reactions and Jacobian. The periodic square and parameters are teaching choices; ellipse geometries and Figures 1–4 are not reproduced.
L. Giomi, Geometry and Topology of Turbulence in Active Nematics, arXiv:1409.1555v2 (2015); Physical Review X 5, 031003
Frontier branch reading on two-dimensional incompressible active nematics, vortex statistics and defects. Those orientation and fluid equations belong to FR22; their results are outside the one-dimensional scalar cortex model.