Parametric Pattern Selection in a Reaction-Diffusion Model

We compare spot patterns generated by Turing mechanisms with those generated by replication cascades, in a model one-dimensional reaction-diffusion system. We determine the stability region of spot solutions in parameter space as a function of a natural control parameter (feed-rate) where degenerate patterns with different numbers of spots coexist for a fixed feed-rate. While it is possible to generate identical patterns via both mechanisms, we show that replication cascades lead to a wider choice of pattern profiles that can be selected through a tuning of the feed-rate, exploiting hysteresis and directionality effects of the different pattern pathways.


S1 Connection with the Gray-Scott model
If we consider Eq. (1) without diffusion, we have By rescaling a → (k 3 /F )a, b → (k 3 /F )b, t → (k 1 F 2 /k 2 3 )t, we obtain a rescaled Gray-Scott model of the form introduced by Pearson [1]. where and GS refers to Gray-Scott. It is apparent that k GS takes on unphysical negative values for k 3 > k 2 , and thus a one-to-one correspondence between the two models is strictly valid for k 2 ≥ k 3 .

S2 Hopf analysis
To determine the stability of the steady state solutions of Eq. (1) we perform a linear stability analysis and examine the eigenvalues of the Jacobian matrix which are: where  a sign of a Hopf bifurcation which occurs at the curve, Above this curve state 3 is stable whereas below it may display oscillations. For k 3 = k 2 /2, F SN = F H and the saddle-node curve collides with the Hopf curve at a Bogdanov-Takens (BT) point beyond which the Hopf curve ceases to exist. Thus for k 3 > k 2 /2 state 3 is always stable and does not display any oscillations. For k 3 < k 2 /2, below the Hopf curve there may be stable oscillations depending on whether the bifurcation is subcritical or supercritical. In the subcritical case state 3 undergoes a bifurcation from a stable focus surrounded by an unstable limit cycle to an unstable focus (and only then state 1 remains as an attractor). In the supercritical case, a stable focus converts into an unstable focus surrounded by a stable limit cycle such that the concentrations may display stable oscillations. The point separating the super-and subcritical bifurcations is referred to as a Generalized Hopf (GH) point, located at k 3 = k 2 /4. The Hopf curve along with the BT and GH points are shown in Fig. S1.

S3 Turing analysis
To determine the conditions for the Turing instability, we employ the ansatz for i = 1, 2 (corresponding to a, b) where λ is the eigenvalue and q a wavenumber. Inserting Eq. (S7) into Eq. (1) we obtain an equation identical to Eq. (S4), however, now with For the equilibrium solutions to be unstable, we require that the real part of the eigenvalues change sign from negative to positive. It is straightforward to show that α > 0 for any non-zero positive value for q, therefore we look for a non-zero positive value for q that makes θ = 0 such that λ = 0 (cf. Eq. (S4)). State 1 is always stable for any q value so we focus on state 3. We notice that θ is a parabola in the quantity q 2 that opens upwards, being positive for q 2 = 0 (for temporally stable solutions) and positive for large q 2 . A condition for linear instability in the presence of diffusion is then obtained by determining the point at which the minimum of the parabola first becomes negative. Setting the derivative of θ with respect to q 2 to zero, we learn that the minimum occurs at the wavenumber q c given by We also require that this wavenumber be real (guaranteed by F > F SN ) and positive, leading to the condition where σ = D b /D a . Substituting q c into θ and setting it to zero then gives the Turing curve: The positivity of q c together with (S10), enforces the inequality (S12) At this point F T = F SN and the Turing curve collides with the saddle-node curve, beyond which it ceases to exist. For σ = 1 this happens at the BT point, while σ can be adjusted, such that there is overlap (or lack thereof) with the Hopf regime. F T is shown as a blue curve in Fig. S1. For the parameter set k 2 = 1.3, k 3 = 1.5, D a = 1, D b = 50, we find F T = 4.5116 and q c = 0.6897. Therefore, a Turing pattern with wavenumber q c corresponds to a n-spot pattern with n = Lq c /2π spots in a system of size L, in our case (L = 200) to a spot number of n = 21.9544.
While at threshold, there is an unique unstable mode with wavenumber q c , above threshold, there is whole band of unstable wavenumbers, i.e., modes with a positive growth rate (real part of the eigenvalue λ of the linear stability analysis). The most unstable mode q max can be numerically determined through the condition d dq Re(λ) = 0. (S13) The resulting q max can be transformed into a number of spots by n T = Lq max /2π, where L is the size of the medium. It is important to point out that not every unstable wavenumber q (or n) corresponds to a stable pattern. A weakly-nonlinear analysis can determine the set of wavenumbers that give rise to stable n-spot patterns.

S4 Spot analysis
In our spot analysis, we follow the work done by Muratov and Osipov [2,3], who-starting from the Gray-Scott model-have extensively studied pattern formation in reaction-diffusion systems of the form: where a is the activator and b is the inhibitor/substrate. The functions q, Q are the interaction terms, θ is a constant parameter (such as the feed-rate F ), τ a 1,2 are their respective time-scales and l a 1,2 are their characteristic length scales of variation. They show that in general, for localized pattern to emerge, the limit l a l b must be satisfied. This is interpreted as the length scale of the activator a varying on a much smaller scale than the inhibitor/substrate b. The formation of spikes/spots occurs in the background of the uniform state 1 and when it does, the size of the spot is of O(l a ). Furthermore the stability of this pattern to perturbations is determined by the relative time scale of variation τ a /τ b . Therefore, the spot is characterized by a sharp peak of the activator, its variation in space must be much smaller than the inhibitor. When its decay rate is smaller than that of the inhibitor, it will be stable. However when it decays faster than the inhibitor (or indeed on the same order) then various kinds of instabilities may manifest themselves, among them spot replication.
To compare our system with the general form (S14), we first rescale Eq. (1) via the following transformations: Note that the steady-state value of b in state 1 is F/k 3 , thus the rescaling of the equations is with S-5 respect to deviations from state 1. After rescaling we get The characteristic time scale for each species are: We also notice a characteristic length scale for the species a 1,2 defined as which can now be interpreted as a diffusion length. Focusing our attention on a, b and comparing (S14) with (S17) we see that: We find it convenient to define the following dimensionless parameters: where σ = D b /D a . In the limitl 1, we can re-write the equations for a, b as, (S19)

S4.1 Example solution
In general, different types of spike solutions (in terms of its profile) will exist in the different regimes of the limitl 1. As an example, we provide a simple case which is accurate in the limit l θ 2 1. In one dimension and in the steady-state Eq. (S19) reduces tõ Since b varies on the order of unity, and a varies on the order ofl 1, one can separate scales inside and outside the spike. Assuming that within the spike b = b s (s denotes spike) is roughly a constant, substituting it into the first equation in (S21) and solving for a we get, where a(0) is the amplitude of the spike. Away from the spike, a = 0 (since this is in the background of state 1). Substituting this into the second equation in (S21) and solving for b we have Matching this solution with the condition b(0) = b s we get the following expressions: where θ c = 12l. We see thus that if θ < θ c , the spike solution does not exist, whereas for θ ≥ θ c there are two solutions. The one corresponding to the positive sign has larger amplitude and is always stable. Whereas the lower amplitude solution is always unstable. Along similar lines for the critical feed rates F SN,H,T (corresponding to the saddle-node, Hopf and Turing bifurcations in the previous section) we can define a critical feed rate F sp for the formation of spots. Using the condition θ ≥ θ c we find, such that spots exist for F ≥ F sp . In Fig. S1 we plot F sp as a red curve for qualitative purposes, keeping in mind that this is an approximate solution and breaks down with increasing k 3 .

S5 Expanding fronts
An alternative mechanism leading to the formation of Turing patterns is provided by expanding fronts. In particular, a perturbation of state 1 within the Turing regime, can lead to an expanding front that in its wake leaves a stationary periodic spot pattern (similar to what has been observed in the Gray-Scott model [4]). In Fig. S2, we show three examples of this for a fixed F = 2.96.
In the left panel, we see the result of perturbing state 1, where a slowly expanding front leaves a spot array in its wake. In the central panel a small localized perturbation to state 3 leads to a much faster expanding front than in the previous case. This is because state 3 is unstable in the Turing regime, while state 1 is stable. In the right panel we show the result of a slightly different perturbation to state 3 leading to an array with different n than in the middle panel.
The figures thus illustrate that the asymptotic state of patterns generated by this mechanism is quite sensitive to initial conditions. An averaging over an exhaustive state of initial conditions is thus required to determine the statistical properties of the dynamics. However, this is beyond the scope of this paper.