Spots, Stripes and the Two Chemicals That Might Explain Both

Somebody Guessed This Before Anybody Could Look

Turing published the idea in 1952 [1]. No assay then in existence could measure a morphogen, and the machines that might have integrated his equations were still being soldered together, so he had nothing to test the claim with except the claim itself; what he did have was a question about how a ball of identical cells becomes an animal with a front and a back, and a hunch that the answer might be arithmetic rather than instruction.

His claim was strange enough to deserve stating slowly. Take a sheet of tissue in which two substances are spread perfectly evenly. Let them react with each other. Let them also diffuse, one faster than the other. Turing's result is that an arrangement this bland can be stable against every uniform disturbance while remaining unstable against a patterned one, so that the sheet, left entirely alone, has no option but to break its own symmetry. Spots appear. Their spacing is inherited from nothing at all, and the chemistry and the two diffusion rates fix it between them.

Diffusion is supposed to smooth things out. Smoothing is the entire reputation of diffusion, and a drop of ink in water has never once organised itself into a grid, however long anybody has sat and watched it. Turing's point is that a second diffusing species changes what smoothing means, so that under the right conditions the smoothing of one substance becomes, physically, the sharpening of the other.

He was dead two years after the paper appeared. He never saw the chemistry confirmed. Confirmation took until 1990, when the Bordeaux group produced a sustained Turing pattern in a gel reactor running the chlorite-iodide-malonic acid reaction [5], and until 1991, when Ouyang and Swinney watched the same system cross from a uniform state into hexagons and stripes [6]. Thirty-eight years is a long while for a prediction to sit unverified. He was not there for any of it.

We can do one thing he could not, which is run the equations cheaply, thousands of times over, and often enough that every number below arrives with an error bar attached to it. That much is the whole of this study. We are not testing whether leopards work this way; we are testing whether the mechanism does what he said it does, and then asking what separates a field of spots from a field of stripes.

What The Grid Actually Holds

We used the Gray-Scott system, a two-variable reaction-diffusion model that grew out of work on autocatalysis in stirred tank reactors [2] and became a standard numerical test bed after Pearson mapped its behaviour in 1993 [3]. A substrate \(u\) is fed in at rate \(F\) and consumed by an autocatalyst \(v\), which is removed at rate \(F + k\):

The domain is a 256 by 256 square with periodic edges and unit grid spacing, the Laplacian is the five-point stencil, and time advances by forward Euler at \(\Delta t = 1\). Diffusion coefficients are \(D_u = 0.16\) and, in the standard setting, \(D_v = 0.08\), so the substrate spreads twice as fast as the autocatalyst, which is the whole of the asymmetry the mechanism needs. Float64 throughout, except in the parameter sweep, where float32 is accurate enough and four times quicker.

Here is what the model does not contain. Two chemical species stand in for a signalling network that in any real tissue runs to dozens of them. Nothing inside it moves except by diffusion, and that matters, because zebrafish stripes are made by pigment cells that chase and kill each other rather than by morphogens seeping through tissue [12]. The square is flat and fixed, where an embryo grows and curves, and growth alone changes which pattern you get [14][16]. The reaction terms are exact mass action with a cubic autocatalytic step, chosen because a person can analyse them by hand, not because any cell anywhere does this. No noise term appears anywhere in it, and randomness enters exactly once, in the initial condition, never again.

Nothing in this article was observed, in an organism or a gel reactor or anywhere else that has a smell. Every number below is printed output from one program, run once, from one seed, and the raw text file is linked at the top of the page. When we say "measured" we mean measured from the model's own output, the way you measure the frequency of heads by flipping a coin two hundred times.

The Condition, Written Out

Turing's condition is not vague. Four inequalities carry the whole of it, and for this system every one of them can be written down in closed form, which is exactly why we chose Gray-Scott over something more biological.

Before the algebra, the picture. Suppose a patch of the sheet happens to run slightly rich in the autocatalyst, which makes more of itself and eats the substrate around it as it goes; the substrate, being the faster mover of the two, is drained out of a wide neighbourhood, while the autocatalyst stays more or less where it was made and piles up locally. The rich patch grows and starves a ring around itself at the same time, and once the sheet is dotted with patches each sitting inside its own exhausted ring, the spacing between them is whatever the two diffusion lengths can agree on. Turing's four inequalities are that story, written as algebra.

A non-trivial uniform steady state satisfies \(uv = F+k\), so \(u\) solves \(Fu^2 - Fu + (F+k)^2 = 0\) and

Such a state exists only when \(F \ge 4(F+k)^2\). Evaluating the reaction Jacobian there and simplifying with \(uv = F+k\) gives a matrix with a pleasant amount of cancellation in it:

Without diffusion the state is stable exactly when \(\mathrm{tr}\,J < 0\) and \(\det J > 0\), which together reduce to the single requirement \(v^{*2} > \max(F, k)\). Add diffusion and a perturbation of wavenumber \(q\) grows when

So the four conditions are: \(\mathrm{tr}\,J < 0\), \(\det J > 0\), \(B > 0\), and \(B^2 > 4 D_u D_v \det J\). The first two say the uniform state survives being nudged uniformly; the last two say the same state does not survive being nudged unevenly, one patch of the sheet pushed up while the patch beside it goes down.

The equal-diffusion case falls out in two lines and this is the part worth reading twice. Put \(D_u = D_v = D\). Then \(B = D\,(f_u + g_v) = D\,\mathrm{tr}\,J\). Condition (i) says \(\mathrm{tr}\,J < 0\), so \(B < 0\), so condition (iii) fails, for every positive \(D\) and every parameter pair whose uniform state was stable to begin with. No Turing instability exists when the two chemicals spread at the same speed. Nothing here is a numerical finding we could have got wrong. The algebra is an identity. The equal-diffusion runs reported below check our code rather than Turing.

The Arithmetic

Work the validation point by hand. \(F = 0.0740\), \(k = 0.0620\), so \(F + k = 0.1360\) and \(4(F+k)^2 = 4 \times 0.018496 = 0.073984\), which is just under \(F = 0.074\). A steady state exists, barely. The discriminant \(1 - 0.073984/0.074 = 1 - 0.999784 = 0.000216\), and its square root is \(0.014704\). Two roots follow. Take the lower, because its steady state is the one that holds still without diffusion:

Square the second one: \(v^{*2} = 0.076209\). Now everything is subtraction.

\(\mathrm{tr}\,J = k - v^{*2} = 0.062000 - 0.076209 = -0.014209\). Negative, so condition (i) holds. \(\det J = (F+k)(v^{*2} - F) = 0.136 \times (0.076209 - 0.074000) = 0.136 \times 0.002209 = 0.000300\). Positive, so condition (ii) holds. The uniform state is stable against a uniform kick, and it is stable by a slim margin, because \(v^{*2}\) beats \(F\) by only two parts in a thousand.

Then the diffusion term. \(B = D_u(F+k) - D_v(v^{*2} + F) = 0.16 \times 0.136 - 0.08 \times 0.150209 = 0.021760 - 0.012017 = 0.009743\). Positive, so condition (iii) holds. And \(B^2 = 9.4932 \times 10^{-5}\) against \(4 D_u D_v \det J = 0.0512 \times 0.000300385 = 1.5380 \times 10^{-5}\), a factor of 6.17, so condition (iv) holds too.

The critical wavenumber is where \(h\) bottoms out, at \(q_c^2 = B / (2 D_u D_v) = 0.009743303 / 0.0256 = 0.380598\), so \(q_c = 0.616926\) and the critical wavelength is \(2\pi / 0.616926 = 10.18\) cells. The fastest-growing mode is a different question, because it maximises the eigenvalue rather than minimising \(h\), and it comes out at \(q_{\max} = 0.489998\), a wavelength of \(2\pi/0.489998 = 12.82\) cells, growing at \(0.014986\) per unit time. That growth rate is an e-folding every 66.7 time units, which on a grid stepping at \(\Delta t = 1\) means the pattern is visible after a few hundred steps and not before.

One correction before we compare anything, because the grid cannot hold \(q_{\max}\) exactly and never could. The nearest available mode on a 256-cell periodic square is \((9, 18)\), with radius \(\sqrt{81 + 324} = \sqrt{405} = 20.1246\) and wavelength \(256 / 20.1246 = 12.7207\) cells. The continuum answer sits 0.1021 cells above the discrete one, 0.80% high, and every comparison from here on uses the discrete number, because that is the only one the simulation could conceivably produce. Four numbers, one of them a wavelength, and not one of them adjustable afterwards.

Mode By Mode

A phase diagram is a weak test of a theory, because a picture that looks like the right picture can be produced by half a dozen mechanisms that have nothing to do with each other. Linear analysis predicts an exact growth rate for every single Fourier mode the grid can hold, so we went and measured every mode we could reach. A single mode with a single predicted growth rate is a claim the theory cannot wriggle out of, which is what makes it worth the twelve separate runs.

The protocol: start from the uniform steady state, add a single cosine of wavenumber \((m, 0)\) with amplitude \(10^{-7}\), run 400 steps, and measure how fast that one mode's amplitude grows between step 200 and step 400. The first 200 steps go in the bin, so that the decaying eigenvector has died away and what remains to be measured is the growing one, alone and uncontaminated. The analytic comparison is the logarithm of the larger eigenvalue of \(I + \Delta t\,(J + \mu D)\), where \(\mu\) is the exact symbol of the five-point Laplacian rather than the continuum \(-q^2\).

Twelve modes, \(m = 12\) down to \(m = 28\), wavelengths from 21.3 cells to 9.1. The worst relative disagreement across all twelve was \(1.08 \times 10^{-5}\). Linear theory is not approximately right here. Right to better than one part in ninety thousand, the analysis describes the grid's behaviour to a precision the grid itself can barely represent.

Figure 1 carries the whole Turing condition in a single picture. The curve is positive over a band of wavenumbers and negative everywhere else, which is what "stable to uniform perturbation, unstable to patterned perturbation" looks like once somebody plots it. The band runs from 7.359 cells to 35.020 cells, and outside it every mode shrinks.

The dashed line running across the same figure is the equal-diffusion case. Set \(D_u = D_v\) and no mode anywhere has a positive growth rate; the best any mode can manage is \(-0.007104\) per unit time, which is decay. The line is flat in the figure because we plotted the ceiling rather than the curve, and the point of it is that the ceiling sits below zero.

Take the difference away and the pattern dies

Theory says the pattern cannot exist when the two chemicals diffuse at the same rate, so a simulation that produced one anyway would mean our code was wrong rather than our reasoning. This control tests us, not Turing.

Three values of the shared coefficient, \(D = 0.16\), \(0.12\) and \(0.08\), three replicates each, every run starting from the steady state plus noise of amplitude \(10^{-3}\). The reaction Jacobian has eigenvalues \(-0.0071044 \pm 0.0158086 i\), a complex pair, so the decay is an oscillation inside a shrinking envelope with a period of 397.5 time units.

The right comparison is subtler than it looks. We got it wrong first. The slowest-decaying mode overall is the uniform one at \(q = 0\), decaying at \(-0.0071044\), but amplitude here means the standard deviation of \(v\), which subtracts the mean and therefore cannot see the uniform mode at all. The mode to compare against is instead the slowest non-uniform one the grid can hold, \((1, 0)\), decaying at \(-0.0071000\) for \(D = 0.16\), and the nine runs had to match that one. Nine measured decay rates, worst disagreement with the analytic value: 2.78%.

After 3000 steps the largest surviving amplitude anywhere in the nine runs was \(6.469 \times 10^{-14}\), which is the float64 round-off floor of the scheme and nothing else. The unequal-diffusion runs saturate at 0.1553, putting the ratio between the two at \(4.2 \times 10^{-13}\). The pattern does not weaken when you equalise the diffusion coefficients. It ceases to exist.

Working Notes From The Club Table

A Phase Diagram With A Hole Where The Theory Is

Thirteen values of \(F\) crossed with thirteen values of \(k\), 169 points, each integrated 6000 steps from Pearson's protocol: the trivial state \(u = 1\), \(v = 0\) disturbed by finite-amplitude patches [3]. Every finished field was classified automatically by four numbers, with the thresholds written down and fixed before the sweep ran, so that nobody could tune a boundary to taste afterwards.

A field with \(\mathrm{sd}(v) < 0.010\) is uniform. A field still drifting by more than \(2.5 \times 10^{-4}\) per step over its last 500 steps is chaos, because it never settles down into anything. Of what remains, global alignment above 0.35 counts as aligned stripes and local anisotropy below 0.55 counts as spots, with everything else falling into maze. Local anisotropy is the gradient structure tensor of \(v\) smoothed over half a wavelength, reported as \(|\lambda_1 - \lambda_2| / (\lambda_1 + \lambda_2)\), which runs from 0 for a round blob to 1 for a perfect stripe; global alignment is the same quantity computed once for the whole field, so a labyrinth pointing everywhere at once scores near zero.

The census: 123 uniform, 14 spots, 18 maze, 14 chaotic, 0 stripes. Two entries there need explaining. One of them is the point of this section.

The easy one is the zero, and the reason for it is geometric rather than chemical. Nothing in a periodic square started from isotropic noise tells a stripe which way to point, so what forms is a labyrinth: locally striped, globally directionless. To check that aligned stripes are a real solution rather than something the model cannot do, we broke the symmetry by hand and seeded a one-dimensional cosine at two Turing points. Both reached a global alignment of exactly 1.000 after 4000 steps, at amplitudes 0.1597 and 0.1671. Aligned stripes are perfectly stable once they exist. They are simply never chosen. A zebra is not a periodic square, and that difference is the whole reason real animals carry oriented stripes while our model produces mazes; work on how gradients and tissue anisotropy orient Turing patterns addresses exactly this gap [18].

The hard one is the dashed ring in Figure 2. Across the entire 169-point window, exactly one point satisfies the strict Turing conditions. We then scanned a much finer grid, 10,556 points covering \(F\) from 0.005 to 0.12 and \(k\) from 0.030 to 0.075, and found 41 Turing-unstable points, 0.4% of them, confined to \(F\) above 0.049 and \(k\) between 0.054 and 0.062. The spots and mazes that made Gray-Scott famous mostly live somewhere else. They are nucleated by finite-amplitude patches from a different steady state, which is the bistable behaviour mapped in the 1990s [10], worked out as a bifurcation structure more recently [9], and the reason the same system also produces self-replicating spots [8]. Two different routes to a pattern, sharing a figure.

Restricting to the 25 points where the minority phase covers at least 10% of the area, so that the pattern genuinely fills the domain rather than surviving as three lonely spots, the wavelength is 12.817 ± 2.878 cells. Over that set \(F\) varies by a factor of 8.2 and \(k\) by a factor of 1.53, and the wavelength moves by 22.5% either side of its mean. A diffusion-set length scale should behave exactly like that. The reaction rates change the shape of the pattern. The diffusion coefficients change its size.

The Wavelength We Got Is Not The Wavelength We Predicted

We would rather not have had to write this section, which is why we wrote it more carefully than any other part of the article, and why it carries more arithmetic than argument.

Twenty-four independent replicates at the validation point, each starting from the uniform steady state plus white noise of amplitude \(10^{-4}\), each with its own random stream spawned from the master seed. We read the dominant wavelength twice: once at \(t = 261\), while the pattern is still growing and the dynamics are still linear, and once at \(t = 2500\), when it has stopped changing altogether. Each reading takes the peak of the mode-count-normalised radial power spectrum of \(v\), with parabolic interpolation between neighbouring bins to beat the bin width down a little.

In the linear phase: 13.6371 ± 0.1372 cells against a predicted 13.5612. The difference is \(+0.0759\) cells, \(+0.56\%\), and z = +0.55. Linear theory stands or falls on that one number, and it stands.

At saturation: 13.8870 ± 0.0877 cells against the fastest-growing mode's 12.7207. The difference is \(+1.1662\) cells, \(+9.17\%\), and z = +13.30. Nothing about that is a rounding error, and we are not going to bury it somewhere near the references and hope nobody reads that far.

Here is what we think is going on. An assumption failed, not Turing's analysis. The growth rates in "Mode By Mode" match to one part in ninety thousand, so the linear theory is not wrong; what is wrong is the unstated step from "this mode grows fastest" to "this is the mode you end up with." The unstable band runs from 7.36 to 35.02 cells, and across the middle of it the growth rate is so nearly flat that at the realised wavelength of 14.22 cells the rate is 0.014551, only 2.91% below the peak of 0.014986. Nonlinear saturation reselects inside that almost level band, and it reselects longer. Weakly nonlinear analysis of exactly that reselection is the reason hexagons and stripes hold different stability ranges in chemical Turing systems [6], and it is not a small correction.

We checked three explanations that turned out to be wrong, and each one deserves listing, because each is the sort of thing that quietly rescues a bad result if you let it.

Was it the Fourier binning? On a 256 grid the bins near the peak are 0.6737 cells apart, 5.26% of the wavelength, which is uncomfortably close to the size of the effect. So we re-ran the whole protocol on a 512 by 512 grid, where the bins are 0.3282 cells apart, 2.56%, and the saturated wavelength came out at 13.8186 ± 0.1150, which is 0.0683 cells from the 256 answer, 0.47 standard errors. The finer mesh did not move it.

Was it the peak finder? We recomputed the same fields with the spectral first moment instead of the interpolated peak, a completely different statistic, and it gives 13.7260 ± 0.0542, still 7.90% above the fastest mode. Same conclusion.

Was it the prediction? Partly, and this is the one that moved, which took us a humbling while to accept. Our first comparison for the linear phase used \(\exp(2\gamma t)\) as the power in each mode, which sat 6.98% away and cost us two evenings of suspicion aimed at the simulation. Running the scheme's own one-step propagator forward for white-noise initial data, and pushing the result through our own Fourier binning, gives 13.5612 and the 0.56% agreement quoted above. The simulation had been right the whole time. The prediction was the sloppy half.

Turn The Diffusion Ratio And Watch The Size Change

This part was fun. Hold \(D_u\) at 0.16, bring \(D_v\) down step by step, and the ratio \(r = D_u/D_v\) climbs from just above its critical value all the way to 5, which is the closest thing to a knob this system has. Theory says the pattern should get finer, and that for large \(r\) the critical wavenumber approaches \(q_c^2 \to r\,g_v / (2 D_u)\), a wavelength falling as \(r^{-1/2}\). The physical reason is unmysterious: the substrate's reach decides how wide a neighbourhood a single patch can starve, so shortening the autocatalyst's reach shortens the territory each patch can hold, and the spots crowd in against one another. Three replicates at each of eleven ratios.

One complication: the growth rate varies by a factor of thirty across this sweep, from 0.001492 at \(r = 1.45\) to 0.046583 at \(r = 5.00\), so reading every run at the same clock time would compare a pattern that has barely started against one that has long since saturated. The reading time is therefore chosen separately at each ratio to put the amplitude at about 0.005 when we look, which is why the \(t\) column below runs from 2623 steps down to 84.

Across all eleven ratios the mean absolute difference between club and theory is 0.3526 cells, which is 3.36% of the predicted wavelength, with a mean signed difference of only \(+0.0206\) cells. The largest single gap is 0.8225 cells at \(r = 2.30\). The z column is not the right thing to read row by row, and we have printed it only so that nobody accuses us of hiding it: with three replicates the standard error is itself badly determined, so a row can post a z of 25.86 purely because its three runs happened to land on top of each other. The mean absolute difference is the honest summary.

Fitting a power law to the linear-phase wavelengths gives an exponent of \(p = -0.6778 \pm 0.0359\) for the club measurement against \(-0.6873\) for the linear theory evaluated the same way. The gap is \(+0.0095\), which is 0.26 standard errors. Both sit well away from the asymptotic \(-1/2\), and they should: the exact relation is \(q_c^2 = (r g_v + f_u)/(2 D_u)\), and with \(f_u = -0.150\) the offset still matters at a ratio of a few. Club and theory bend by the same amount, which is the thing being tested.

Then the sweep produced something we had not planned for. Below \(r \approx 1.8\) the pattern grows to a visible amplitude and the entire field then collapses onto the other stable state, \(u = 1\), \(v = 0\), leaving a blank square. The two open circles in Figure 5 are runs whose linear phase matched theory beautifully and whose final state was nothing whatsoever, a blank square where a pattern had been. The Turing bifurcation here is subcritical. Nothing small and patterned exists for the growing mode to settle onto, so the overshoot escapes to the far attractor and stays there. Linear stability analysis says nothing about this and cannot, because it is a statement about infinitesimal perturbations and this is a statement about where they go when they stop being infinitesimal.

The Strongest Objection

Put the best version of the case against us. It goes like this.

You have shown that a particular pair of partial differential equations, integrated on a square, produces blobs whose spacing follows a formula you derived from the same equations. Arithmetic agreeing with itself. It has no bearing on animals. Leopards are not squares and morphogens are not two, so you have demonstrated nothing whatever about biology except that a certain equation happens to be solvable. Worse, pattern-matching is a notoriously weak form of evidence: the number of mechanisms that can make evenly spaced spots is large, and chemotaxis and cell-cell contact rules both do it without a single diffusing morphogen anywhere in the tissue, which means a match between your field and a photograph of an animal carries almost no information about the animal.

Most of that is correct and we accept it. A matching pattern is not proof that an animal uses this mechanism, and anybody who tells you the leopard's spots are explained because a simulation made similar-looking spots is selling you a picture. Murray's 1981 pre-pattern argument [4] produced strikingly animal-like coats. The similarity settled nothing.

What has moved the biological question forward is not resemblance but perturbation. Kondo and Asai's angelfish showed stripes rearranging as the fish grew, in the specific way a reaction-diffusion wave rearranges and a fixed pre-pattern cannot [7], which is a prediction coming true rather than a resemblance being noticed. Sick and colleagues predicted from a WNT/DKK model what should happen to hair follicle spacing if you over-express the inhibitor, then did it, and the spacing changed as predicted [11]. Sheth and colleagues knocked down Hox genes and got more, thinner digits, which is what a Turing mechanism with a shortened wavelength does and is difficult to get any other way [16]. Raspopovic and colleagues named the three interacting components in the limb bud [17]. Economou and colleagues cut a ruga out of a mouse palate. The neighbouring stripe bifurcated to fill the gap [14]. Zebrafish pigment cells turn out to implement the required short-range activation and long-range inhibition through cell contacts rather than diffusion, which is a Turing system in the mathematical sense with no morphogen in it at all [12][13].

So the honest position is narrow, and we would rather state it narrowly than let a reader walk away with something wider. We have confirmed that the mechanism works as advertised inside its own mathematics, that its central condition is exactly the condition Turing identified, and that removing the differential diffusion destroys it completely. Whether any particular animal uses it, our grid cannot say.

The second half of the objection deserves its own answer. Several modelling choices we made would have changed our numbers, and here is where.

Reading the wavelength only at saturation would have produced a headline of "linear theory is 9% wrong," which is a defensible thing to write and, we think, a misleading one. Reading it only during growth would have produced "linear theory is exact," which is defensible in the same way and misleading in the same way, and rather more flattering to us. We report both. Choosing Gray-Scott rather than an activator-inhibitor system in the Gierer-Meinhardt family [20] gave us closed-form conditions and a bistable background, and that bistability is precisely what produces the collapse at \(r < 1.8\) and the hole in the phase diagram; a supercritical system would have shown neither, and the phase-diagram section would have had nothing surprising in it. The periodic square guarantees zero aligned stripes in the census; a domain with an edge or a gradient across it would not, and we did not test that, so treat the zero as a fact about our boundary conditions rather than about the model. The five-point Laplacian is mildly anisotropic on a square lattice, which is why we stopped the ratio sweep at \(r = 5\) where the predicted wavelength is still 7.4 cells; a nine-point stencil would let the sweep run further, and we do not know what it would find there.

What We Would Tell Somebody Who Asked

Turing guessed right about a mechanism he had no instrument to check, and he died two years later without ever learning that he had guessed right.

Two chemicals, one reacting with the other, one diffusing faster, on a sheet with nothing written on it anywhere, and the uniform state holds against any uniform nudge while failing against an uneven one. The wavelength of that failure follows from four numbers and a pencil. Measured during growth, it came within 0.56% of the prediction, a z of 0.55 across 24 runs; set the two diffusion coefficients equal and the pattern falls to \(6.5 \times 10^{-14}\), the noise floor of double-precision arithmetic. Nothing here approximates the Turing condition or lands near it by coincidence. The grid obeys the Turing condition itself.

What decides spots against stripes has a longer answer than the question implies. In our sweep the feed and kill rates decide the area fraction: the spot fields put 11.8% of the domain into the high-\(v\) phase on average, the maze fields put 46.6%. Round objects when one phase is scarce. Connected ridges when the two phases sit near balance. But which of spots or stripes you actually observe also depends on things outside the reaction terms entirely, including whether anything in the domain ever tells a stripe which way to lie. Our domain does not, so it made 18 mazes and no stripes, and a zebra would have made something else entirely.

And the honest caveat, which should be the last thing anybody takes away. A model that makes a pattern like a leopard's is no evidence that the leopard makes its pattern that way; every biological result that actually moved this field did so by breaking something and predicting correctly what would break [11][14][16][17]. We broke our own model on purpose twice, once by equalising the diffusion coefficients and once by refining the mesh, and it behaved as the theory said it would. A computation can do no more than that.

Reproducing this

Python 3.12 and numpy. Nothing else. From the repository root:

The run that produced this article took 301.0 seconds of wall clock on a six-worker process pool. We re-ran it from scratch while drafting, on the same machine, and it took 285.9 seconds and produced a file byte-identical to the first apart from the six lines that print their own timings. Master seed 20250621, hard-coded at the top of the script; every random stream is spawned from it through numpy's SeedSequence, so the result does not depend on how many worker processes you give it. Recorded with numpy 2.4.2 on Python 3.12.3. Longest stage: the 169-point phase diagram, about 167 seconds. The 512 by 512 control comes next at about 66.

The interactive bench linked at the top of this page runs the same equations in a browser at lower resolution, with the feed and kill rates on sliders, and reports the predicted fastest-growing wavelength beside the one it measures from the running field. Its defaults are the parameters used in Figure 3 panel A.

References