Measure a growth rate from a simulation record

The box first selects the wavelengths

In a square of side $L$ , opposite edges are identified: periodic boundaries. A wave must return to its starting phase after crossing the box. Each direction can therefore contain only an integer number of periods.

Label modes by integers $n_x,n_y$ . The zero mode changes every position together. Other modes vary in space. In the default box $L=16$ , transverse mode $5$ has wave number about $1.96$ , inside the unstable continuous interval from the preceding chapter.

A no-flux wall imposes zero normal concentration gradient and admits a different family of standing waves. Boundaries are part of the model: changing them can change the selected mode with identical reaction kinetics.

q_x=\frac{2\pi n_x}{L},\qquad q_y=\frac{2\pi n_y}{L},\qquad q^2=q_x^2+q_y^2

Allowed wave numbers in a periodic square, using the dimensionless length scale of this model.

Reading: A. M. Turing

Extract one component from the whole image

Multiply every grid value by the complex conjugate of a chosen wave and average over the domain. Matching phases add; other complete periods cancel. This gives the Fourier coefficient $\widehat{\delta u}_{n_x,n_y}$ . Its modulus separates amplitude from the phase change caused by translating a pattern.

Growth points come from saved tiny-disturbance simulations, with theoretical eigenvalues drawn separately as reference lines. The colorful main field evolves until its spots are visible. These separate stages test early growth and display the later pattern.

Preparing the two concentration disturbances in the ratio of a growing eigenvector reduces transients. Perturbing only one species generally excites both growing and decaying components, so the earliest interval need not follow a single exponential.

\widehat{\delta u}_{n_x,n_y}(t)=\frac{1}{N^2}\sum_{j,k=0}^{N-1}\delta u_{jk}(t)\exp\!\left[-\frac{2\pi i}{N}(n_xj+n_yk)\right]

Project the disturbance relative to the uniform equilibrium. Each direction has the same number of grid points.

Fit the interval that remains linear

For $A(t)=A_0e^{\lambda t}$ , the natural logarithm is a straight line. Fit $\log A(t)$ against time to estimate its slope. Several samples reveal curvature or transients more clearly than a ratio of two frames.

A very early window may mix several components; a late window includes nonlinear reaction changes. Move both ends of the window and check whether the slope remains stable. Inspect the amplitudes as well: a tiny signal near rounding error can give a misleading line.

In a deterministic simulation without observation noise, fitting residuals reflect transients, nonlinearity and numerical error. They are not uncertainty across biological specimens. Studying noise requires a specified random process, seeds and independent repeated records.

\widehat\lambda=\frac{\sum_r(t_r-\bar t)(\log A_r-\overline{\log A})}{\sum_r(t_r-\bar t)^2}

The ordinary least-squares slope. Amplitudes must be positive and belong to a resolvable linear interval.

A grid has its own spectrum

The five-point stencil approximates the Laplacian using neighboring values. With spacing $h=L/N$ , a wave receives the discrete eigenvalue below instead of exactly $-q^2$ . Short waves show the largest difference.

First compare fitted points with theory for the same grid, which tests time integration and record processing. Then refine the grid and compare with continuous theory. Reducing the time step cannot remove spatial error from an inadequate grid.

Keep the physical box, reaction parameters and mode numbers fixed during refinement. Record the grid, step size, boundaries and fitting window so that another reader can repeat the same comparison.

q_h^2=\frac{4}{h^2}\left[\sin^2\!\left(\frac{q_xh}{2}\right)+\sin^2\!\left(\frac{q_yh}{2}\right)\right],\qquad A_h=J-q_h^2\operatorname{diag}(D_u,D_v)

Continuous and grid references use the eigenvalues of their respective spatial operators. They approach one another as the grid is refined.

Worked example

A first estimate from two frames

Suppose saved amplitudes at $t=0.5$ and $t=1.5$ are $0.00120$ and $0.00155$ . Their separation is one time unit, giving $\log(0.00155/0.00120)\approx0.256$ . These are synthetic values illustrating the calculation, separate from the experiment’s saved record.

Agreement with a predicted value is a starting point. Inspect intermediate samples and repeat on a finer grid. If the later amplitude begins to flatten, shorten the window: fitting through saturation underestimates the early linear rate.

Try the reasoning

Trace a plotted point back to its record

Choose a growing mode and calculate its allowed wave number from the box length.
Compare the fitted point with grid and continuum references. Identify which discrepancy requires a smaller time step and which requires a finer grid.
Use the stable preset as a control. The same logarithmic analysis should give a negative slope for a decaying mode.

Exercises and explanations

What happens to wave number if the box doubles in size but the mode number stays fixed?
Read the explanation
Wave number halves and wavelength doubles. This moves the mode along the growth spectrum. To test convergence at the original wave number, keep the box and mode fixed while refining the grid.
Why is a clearly developed spot pattern unsuitable for testing an early linear growth rate?
Read the explanation
Its concentrations have left the neighborhood of equilibrium. Reaction slopes change and modes interact. The full nonlinear dynamics govern saturation.
The discrete reference agrees with the record but the continuum reference does not. What should you change first?
Read the explanation
Refine the spatial grid while holding the box, parameters and target mode fixed. The solver may be accurately solving a coarse discrete problem.

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.