Pattern selection in reaction diffusion systems

Turing's theory of pattern formation has been used to describe self-organisation in many biological, chemical and physical systems. However, while conditions sufficient for the existence of patterns are known, the general nonlinear mechanisms responsible for pattern selection are not. Here, we show that the movement and final position of peaks within a Turing pattern are determined by the flow of mass through the system and the minimisation of the total mass of the fast species. The apparent success, for particular parameters choice, of linear stability analysis in predicting the final pattern is found to be a coincidence. Furthermore, mass minimisation also underlies a competition instability that results in peak-coarsening of quasi-steady state patterns and correctly predicts the number of peaks at the true steady-state.


INTRODUCTION
Pattern formation occurs in a huge variety of natural and living systems [1], from chemical reactions [2,3] to living cells [4][5][6] to environmental patterns [7]. In systems described by reaction-diffusion equations, the formation of spatially periodic patterns can be explained by the Turing instability, in which patterns emerge due to the presence of two or more interacting components that diffuse (or are transported) at different rates [8][9][10][11]. The resulting patterns are multi-stable in that several different stable patterns can be obtained from the same set of parameters, albeit, for incompletely understood reasons, with different frequencies.
Sufficient conditions for pattern formation can be determined in the so-called Turing or linear regime, in which a spatially uniform stable steady state becomes linearly unstable to spatial perturbations in the presence of diffusion [9]. Consider the following one-dimensional system over the spatial domain [0, L] with reflexive boundary conditions. The evolution of any small perturbation from a spatially uniform steady state is given by its decomposition into its Fourier modes e λnt cos( nπ L x), n ∈ N, where Re(λ n ) is the growth rate. Then, the uniform steady state is laterally unstable if any mode n has a positive growth rate Re(λ n ) > 0 (see Fig. 1a and Supplementary Information). This linear analysis does not take into account higher-order terms that become important as the pattern evolves. So although we can determine which modes are initially excited, we can not predict, in general, which pattern is finally obtained. Furthermore, many studies have examined the effect of external constraints such as fixed boundary conditions, parameter ramps, template patterns, geometry and deformable boundaries [12][13][14][15], especially near onset (see [10,16] for reviews). Here, however, we study unforced systems away from onset with reflexive or periodic boundary conditions. Then, the naive expectation is that the fastest growing mode will dominate the final pattern and therefore set the wavelength. For certain models and parameters choices this is indeed the case, but not generically. Furthermore, there are several indications, described below, that even when the prediction from linear stability holds, there are other mechanisms acting in parallel to specify the set of obtainable patterns.
First, final patterns have a well-defined wavelength with regularly positioned peaks. By that we mean, that the peaks of the pattern are found at the same locations as for some fundamental mode n. Furthermore, only this mode and its harmonics n, 2n, 3n, . . . contribute to the pattern (necessarily due to periodicity) (Fig. 1b). This is despite neighbouring modes n ± 1 generally having similar growth rates, which would be expected to lead to aperiodic patterns. What underlies this 'exclusion principle' [12] is not known. The regular positioning of peaks is even found in the singular limit D v → 0 in which the pattern consists of well-separated point-like spikes, which necessarily contain contributions from a large number of fundamental modes (see [17] for a review).
Secondly, regular positioning is typically maintained even during domain growth [11,18]. While some models exhibit peak insertion and others peak splitting, in either case peaks move towards a regular positioned configuration. Thirdly, as already stated, the linear 'prediction' does not always hold. A mode that is more slowly growing initially can nonetheless dominate the final pattern. This can occur through a coarsening or competition effect in which a transient pattern evolves to a stable pattern with fewer peaks (Fig. 1c). Like domain growth, coarsening occurs outside the linear regime since a pattern has already been formed. However, the final pattern obtained nonetheless consists of periodic well-positioned peaks, as described above. A notable special case occurs in mass-conserved two-variable systems in which the initial pattern coarsens eventually to a single peak or halfpeak [19,20].
Together, these observations indicate that the observed regular positioning of peaks is not due to the predicted growth of some fundamental cosine modes, which in any case is only valid close to the spatially uniform solution.
Rather, some non-linear mechanism must be responsible, with any agreement with the prediction of linear stability analysis being likely coincidental. In the following, we show that the flow of mass through the system is responsible for the regular positioning of peaks. By flow, we mean something more than simply the flux through the system. In a diffusive non-mass-conserving system, 'molecules' enter the system, diffuse and either leave the system or are converted to another species. This combination of diffusion and turnover results, as we shall see, in peak positioning due to the concept of flux-balance [21,22]. This also explains why the peaks in massconserving two-variable reaction-diffusion systems do not move: there is no flow to drive the movement. Importantly, we show that flux-balance is superseded by a more fundamental and generalisable principle -mass minimisation. In this viewpoint, the peaks of a Turing pattern are positioned so as to minimise the total mass of the fast diffusing species. Furthermore, we show that mass minimisation correctly predicts the number of peaks of the final stable pattern after any coarsening of quasi-stable patterns has taken place. The principle of mass minimisation is therefore an incredibly simple yet powerful concept for the understanding of the non-linear behaviour of pattern-forming systems.

THE MODEL
We introduce the following exploratory onedimensional system, inspired by our recent model of bacterial condensin [23,24], written in terms of the variables u = u(x, t) and v = v(x, t), defined over the spatial domain [−L/2, L/2], with reflexive boundary conditions, all parameters non-negative and D v < D u . While superficially similar to the some of the classic Turing models such as the Brusselator [25] and Schnakenberg [26] models, this model has some notable properties that make some analyses easier. Firstly, in the absence of diffusion, it has a single fixed point that is stable for all parameter values. This means that the stability diagram of the system is particularly simple. There are only two regions, specified by a single inequality: one in which the spatially uniform solution is stable and another in which it is Turing unstable (Fig. 1d). There are no oscillatory instabilities. Secondly, the model has the property that every steady state solution has the same  (2) with default parameters and L = 4. Note that the growth rate at n = 0 is negative -the system is not generically mass conserving. b. The Fourier decomposition of the obtained two-peak pattern (inset). c. An example kymogrpah showing pattern development starting for random perturbation of the uniform state. While mode n = 7 dominates initially, the pattern coarsens down to two peaks, dominated by mode n = 4. d. The Turing space of the model (blue shaded region).
total concentration Thirdly, the model has the form of a mass-conserving Turing system with additional terms: a global source term, cδ, and two depletion terms, δu and δv. By writing the source term as cδ, we can change δ, the turnover rate, while leaving the total steady state concentration c fixed. We obtain a mass-conserved Turing model when δ = 0 and the limit δ → 0 is well defined as long as we constrain the total initial mass to be the same as the steady-state mass, i.e. C(0) = c. The condition for a Turing instability is most easily stated by non-dimensionalising the system and introducing the dimensionless parameters a = βc 2 γ , b = δ γ , Γ = γL 2 Dv , d = Du Dv (see supplemental text for details). As can be seen in Figure 1d for typically choices of the diffusivity ratio d, we require b 1 for patterning, i.e. the timescale of mass flow through system, 1/δ, must be much longer than the timescale underlying the Turing instability 1/γ. Note that Γ does not enter the condition for a Turing instability but does affect which modes have positive growth rate. Peak movement and regular positioning depend on flux through the system. a. The system is initialised with a peak away from middomain. The peak subsequently moves to mid-domain. b.
The centroid of the peak (blue line) plotted as function of simulation time. The orange dashed line is an exponential fit. Inset: (Top right) The rate of movement obtained from fitting the centroid to an exponential as in b shows a linear dependence on the turnover rate δ. (Bottom right) Peak velocity is linear in peak position. c. A single peak in the mass conserved limit δ = 0 can be positioned anywhere on the domain. No peak movement is observed d. The mass conserved system exhibits complete coarsening. Irrespective of how many peaks there are initially, the pattern eventually coarsens to a single peak, the position of which depends on which peak of the initial pattern has not coarsened. In d Γ = 19200.
Numerically solving the system, we found that it indeed produces regularly positioned peaks. We also observed that, like the model it is based on [23], it exhibits a competition or coarsening instability in that the final dominant mode has a shorter wavelength than predicted by linear stability analysis. For our default parameter set with L = 4 (Γ = 4800), linear stability predicts (Fig. 1a) that the pattern consists of four peaks (or valleys) (mode n = 8) whereas the obtained steady-state pattern most frequently consists of two peaks (mode n = 4) (Fig. 1b). While multiple peaks often form initially, consistent with the linear prediction, coarsening rapidly occurs, leaving mis-positioned peaks that then move slowly towards opposite quarter positions, while maintaining their shape (Fig. 1c).
To examine the movement of peaks in more detail, we focused on the case of a single peak (n = 2), typically obtained for L = 2 (Γ = 1200). Examining the movement of the peak (Fig. 2a), we found that it moves to mid-domain exponentially in time (Fig. 2b), indicating the peak velocity is linearly proportional to its displacement from mid-domain. This was the case whether the system was initialised with a random perturbation of the uniform state or with a peak preformed somewhere on the domain. Furthermore, the rate of movement was found to be directly proportional to the turnover rate δ (or equivalently cδ the flux through the system per unit length). In the mass-conserved limit, δ = 0, no peak movement was observed (Fig. 2c). Furthermore, as in other massconserving models [20], a coarsening process occurs such that the final steady-state pattern always consists of a single peak (Fig.2d). Importantly, since peaks do not move, the position of this final peak is determined by whichever peak of the transient state remains after coarsening. If the system is initialised with a single preformed peak, then that peak does not move and constitutes a stable pattern (Fig. 2c). Thus, the mass-conserved model has a continuum of single-peak stable states and thus regular positioning is not an intrinsic property of the system. Overall, these results demonstrate a connection between peak movement towards the regular positioned configuration and the flow of mass through the system.

POINT SINKS
To explore this connection in more detail, we turn to a toy model involving diffusion and point sinks. We consider the steady-state diffusion equation for a variable A = A(x) over a one-dimensional domain of length L in the presence of global source and decay terms as well as n localised point sinks at positions x = (x 1 , · · · , x n ) (each with rate µ): We take the domain to be [−L/2, L/2] and impose zeroflux boundary conditions. As before, we write the global source term in terms of the decay rate δ and a concentration c, which is the steady-state concentration in the absence of the point sinks. A simpler system without the decay term and with perfect points sinks (i.e. µ → ∞) was used by Ietswaart et al. to model the positioning of plasmids within rod-shaped bacterial cells [22]. They found that the gradient differential across each sink vanishes if and only if the sinks are regularly positioned and, therefore, if sinks were to move up the concentration gradient, they would be regularly positioned. We will extend this result to the more complicated case of equation (3).
Note that the presence of the decay term introduces an additional length scale D/δ into the system, namely the distance that a molecule of A would diffuse (in the absence of any point sinks) before it decays. We refer to this as the length-scale of diffusion. It is small when either diffusion is slow or the decay rate (turnover) δ is fast. We can write the solution to (3) as where G(x; x i ) is the modified Green's function defined by where the dimensionless parameter κ = L δ D is the ratio of the length of domain to the length-scale of diffusion. The coefficients µ i = µ i (x) are determined by the linear algebraic conditions where we have defined a second dimensionless parameter λ = µ δ , the ratio of the sink and background decay rates. The quantities µ i have a simple interpretation. They are directly related to J i , the flux leaving the system through each sink where J i = |D dA dx | and the − and + subscripts refer to the diffusive flux from the left and right respectively. We also define the flux differential across each sink as Note the total mass (concentration) of A in the system is readily given by where the second term solely describes the effect of the point sinks.
We can now investigate what would in happen in this system if sinks were to move up the gradient of A. As in Ietwaart et al., we can determine the configurations for which the flux differentials are all zero. In the supplementary text, we prove that this occurs uniquely for regularly positioned sinks, Interestingly, we also show that the regular positioned configurationx i is the unique stationary point of the mass M , i.e.
Thus if sinks move up the concentration gradient (in the direction of greatest flux), they will be regularly positioned as this is the configuration for which the fluxes into each sink from either side balance. Furthermore this configuration minimizes the total mass of the system. In other words the sinks are positioned so as to 'consume' mass at the greatest rate. This connection between regular positioning and mass minimisation appears to be generalisable and we have observed numerically that it holds for spatial sinks i.e. if the delta function in equation (3) is replaced by a peak-shaped spatial function such as a Gaussian function or sech 2 (x), then the total mass is minimised when the sink is centred at mid-domain.
Let us consider the case of a single sink, n = 1, in more detail. We focus on the regime κ 1 in which the diffusive length-scale is much longer than the domain size. We expand in κ to find first and then Hence, if κ 1, then the flux differential across the sink depends linearly on its relative displacement from the mid-domain. For strong sinks (λ 1), the proportionality factor is linear in δ, just as we observed for the Turing system (Fig. 2b). As κ increases, the flux-differential becomes inflected about x 1 = 0 (Fig. 3a). We can think of this heuristically as follows. If the diffusive length-scale is much shorter than the domain size (κ 1), then only particles initially created near the sink will fall into it. Therefore the flux-differential is only significantly nonzero close to the boundaries (or another sink). In essence, the geometry sensing of the system breaks down. On the other hand, when the diffusive length-scale is much longer than the domain size (κ 1), particles can explore the entire domain before decaying and so the flux-differential across the sink reflects its position on the domain, with the fluxes into the sink from either side balancing at middomain. The relevance of this dependence on κ to Turing systems will be made clear later.
We can make sink movement explicit by specifying the sink velocities. We make the simplest choice that their velocity is directly proportional to the flux-differential ∆J i , From the results above, the steady-state solution then consists of regularly positioned sinks and this is also the unique configuration that minimises the total mass M . This holds for any dependence of the velocity on ∆J i such that dxi dt = 0 if and only if ∆J i = 0. If we assume that sinks move on a much slower timescale than that of diffusion, then we can use the steady state solution for A given in equation (4) to solve the dynamic system (10). We find that increasing δ, which shortens the diffusive length-scale while also increasing the flux through the system, leads to faster sink movement (Fig.  3b,c) reminiscent of the Turing system (Fig. 2b). On the other hand, if we decrease D, which decreases the diffusive length-scale without affecting the flux through the system, the sink moves more slowly towards mid-domain (Fig. S1).

COMPARISON WITH THE TURING SYSTEM
The similarity between moving point-sinks (Fig. 3) and the movement of peaks in a Turing pattern (Fig. 2) is striking. In the regime κ 1, a point sink that moves with a velocity proportional to the gradient, moves exponential to mid-domain since the gradient across the sink is a linear function of its displacement from middomain. This suggests that the exponential movement of a Turing peak may also be due to a dependence of the peak velocity on the flux-differential (in this case of the fast species across a peak of the slow species). Note that this comparison does not give us insight into Turing patterns containing a peak at the domain boundary and we therefore restrict ourselves in the following to patterns without peaks at the boundary. We introduce the following definition of the flux-differential into the peak of a single-peak Turing pattern: We initialised the system with a single peak and monitored ∆J s as a function of the peak position and velocity. We found that, like for point sinks, the flux-differential is, away from the domain boundaries, directly proportional to the displacement from mid-domain (Fig. 4a). Thus, peaks move with a velocity proportional to the flux-differential. It is not clear how to extend the above definition of the flux-differential to patterns with multiple peaks as well as to higher dimensions in which Turing patterns can consist of complex structures such as stripes, spirals and hexagons. However, the concept of mass minimisation is easy to generalise. When we examined the total mass (concentration) of u in the system, M = u(x, t) dx, we found that it decreases monotonically as the peak moves to mid-domain (Fig. 4a). Further, when we initialised the system with two peaks positioned at various location, we found similar behaviour (Fig.4c) suggesting that, for a given number of peaks, the regularly positioned configuration minimises the total mass of u, just as we have proven for point sinks in the previous section. Indeed, the trajectories show a remarkable similarity to those of moving point sinks (Fig.  3d).
To explore this analytically, we considered the singular limit D v D u in which the peaks in v take the form of narrow spikes or pulses of width order = D v /γ. Away from the spike v is approximately constant with a value v out that is much smaller than u. This limit allows the use of non-linear analysis methods to study the existence, stability and dynamics of Turing patterns [17]. Here, our goal is simply to derive an approximation for u in this limit by treating the spikes of v as Dirac delta functions a. Flux differential measured numerically using equation (11) for a single peak (orange) is a linear function of the peak position. The mass of the fast species M (blue) is plotted is minimised at mid-domain. b. The same quantities as in a but for the analytical expressions from the spike approximation (eq. (16) and M = c − ρ + ). c. Mass minimisation for two peaks. Trajectories of two peaks as they move towards opposite quarter positions (while lines). The contours and colour bar represent the mass M interpolated from trajectories. The mass is minimised for regular positioning d. Same as c but trajectories obtained from the approximation of peaks as spikes using equation (17). Parameters: Dv = 0.0012. L = 2 in a and b, L = 4 in c and d, otherwise default.
as described above.
We look for steady-state solutions consisting of n spikes at positions x 1 , · · · , x n . We assume that u changes slowly within each spike and so can be approximated by a constant u i and within each spike u i v. First we introduce the inner coordinate, y i = (x − x i )/ , within each spike. We then have the following system for the inner variable v i (y) In the outer region, each spike is approximated by a weighted Dirac delta function and we therefore replace the v and uv 2 terms by Dirac delta functions with weights w 1 and w 2 given by respectively. The equation for u in the outer region then becomes and where we have used that, according to the spike approximation, (2b)) to simplify the contribution away from the spike. Note the inverse dependence on u in the point sink term (which we call an inverted sink). This form is also obtained for other Turing systems with a uv 2 non-linearity, such as the Schnakenberg and Brusselator models [27,28].
Following the approach of the previous section, the solution to equation (12) is given by where the Green's function is defined as for point sinks but in terms of the corresponding dimensionless parameter κ = L δ Du , the ratio of the length of the domain to the diffusive length scale of u (henceforth κ replaces b in the set of dimensionless parameters of the system). The coefficients ρ i = ρ i (x) are now determined by the non-linear algebraic system where σ = ρ c 2 δ = 6 √ b+1 a √ Γ . The inverse dependence on u(x i ) makes solving this algebraic system challenging. For a general choice of sink positions x i , there are n coupled quadratic equations in ρ i , and therefore up to 2 n real solutions.
However, not all spike configurations are stable. Based on our numerical observations, stable solutions consist only of regularly positioned spikes of the same height (ρ i = ρ ), also referred to as symmetric spike solutions, similar to other models. The flux-differentials (defined analogously to equation (7)) of these solutions vanish, ∆J i (x) = 0, just like for the points sinks of the previous section.
Solving this system for a single arbitrarily positioned spike, we find two solutions corresponding to different spike amplitudes and corresponding masses M ± = c − ρ ± . Since we only ever observe spikes within large amplitudes, i.e. patterns in which almost all the mass of the system is contained with spikes, we assume that the low amplitude solution is unstable and an artefact of the approximation. We find that the flux-differential across the spike depends linearly on the spike position in the regime κ 1 just as we found for the non-inverted sink case in equation (9). Note also the linear dependence on δ when σ 1 (as it is for our parameter choice), consistent with our numerical observations (Fig. 2b). We also find that M , the total mass of u, is minimised at mid-domain. When we evaluated these quantities using the same parameters as for our numerical results, we found very similar qualitative profiles (Fig. 4b). However, the agreement was not quantitative, likely due to the nature of the approximation and/or because our solution is not sufficiently spike-like. We will see in the next section that our main result is unaffected.
Finally, we considered the case of two moving spikes. As before, spike movement is introduced by specifying that the spike velocities are proportional to the fluxdifferentials where ν is a new parameter. Together with (14), this defines a differential-algebraic system. For two spikes, it has up to four real solutions for each configuration (x 1 , x 2 ). As above, we initialise the system on the solution branch with the smallest mass M . In Fig. 4d, we show the mass M overlaid with several sample trajectories. The similarity to the numerical observation (Fig. 4c) is clear and in both cases, the steady-state solution, consisting of quarter position peaks/spikes, minimises the mass of u.
These results indicate that the movement of a peak in a developing Turing pattern is due to the flux-differential of the fast species into the peak. Furthermore, regular positioning is a result of the flow (creation, diffusion, decay) of mass through the system. It is the configuration for which all the flux-differentials balance or equivalently, the configuration that minimises the mass of the fast species.
We have also seen that the inverted (12) and noninverted (3) sink models have very similar properties modulo the multiplicity of solutions of the inverted model. In the next section we will see where they differ and the important role of the inverted term.

COMPETITION AND PATTERN SELECTION
We have seen that in the mass-conserved limit δ → 0, the model exhibits a coarsening effect in which the only stable pattern is a singe peak (Fig. 2d). For non-zero δ we have also seen that the model exhibits incomplete coarsening. Linear stability predicts that mode n = 8 (four peaks) will dominate (Fig. 1a) but we most frequently obtain a steady state pattern dominated by mode n = 4 (Fig. 1b,c).
This motivated us to explore the connection between coarsening and the flow rate δ in more detail. We measured the distribution of steady state patterns obtained for different values of δ and compared against the prediction of linear instability (Fig 5a,b). We used periodic boundary conditions to avoid peaks on the boundary that are not described by the spike approximation. We found that for κ > ∼ 1 linearly stability analysis correctly predicts the dominant mode at steady state. However, for κ < ∼ 1, a coarsening process occurs and the steady-state pattern is dominated by a lower mode than that predicted. We explain this as follows. When the diffusive length-scale is longer than the domain size, all peaks compete for u molecules created across the domain. Whereas, when the length-scale is short, peaks only absorb molecules of u created within a distance given by the diffusive lengthscale and therefore compete less or not at all. Competition is exasperated by the fact that decreasing δ also decreases the total flux through the system (cδL).
We next applied the spike approximation developed in the previous section. We decreased D v from the default value so that the obtained pattern was reasonably spikelike (Fig. 5c) while at the same time not resulting in a very much enlarged Turing space (since we want to sweep over different values of δ). We considered only symmetric, regularly positioned spike solutions, which are the only observed steady-state solutions. For n spikes, we obtain two possible values of ρ , of which we take the larger, given rise to solution u(x) = c − ρ + i G(x;x i ) with mass M = c − nρ + . Note that for a real solution we must have a > 12κ b+1 Γ coth( κ 2n ). Therefore, for a given choice of parameters, there is an upper bound on the number of spikes that a solution can contain. However, in general, a solution exists for multiple values of n. Yet, numerically, we observe a very narrow distribution of the a. The number of peaks in the most frequent steady-state pattern is plotted as a function of a and κ. For each set of parameters, the most frequent pattern was obtained from 5 simulations each initialised with a different random perturbation from the uniform state. The simulations were ran for long enough to ensure the steady-state pattern was reached. b. The number of peaks in the mode with the greatest growth rate as predicted by linear stability analysis is plotted as a function of a and κ. Plots a and b are similar for κ 1. They disagree for κ ≤ 1, which indicates coarsening. c. Example of a steady state pattern in the spiky limit. d. Normalised total mass M/c = 1 − nρ + /c plotted as a function of n for different values of κ. There exists a critical n for which the mass is minimal. e. Normalised total mass M plotted as a function of κ for different numbers of spikes. As κ → 0, M is minimal for a single spike. f. The numerically obtained distribution of peak number at steady state for different values of κ (colour scale) overlaid with the prediction of the dominant pattern from linear stability (green triangles) and the prediction from mass minimisation (red circles). Mass minimisation correctly predicts the number of peaks at steady state. Data from 50 simulations for each parameter set. Parameters: Default values as in Figure 1 with L = 4 except c-f which use Dv = 0.006 (to make peaks narrower).
number of peaks. We hypothesised that mass minimisation might play a role. Indeed, when we examined the mass M of solutions consisting of different numbers of spikes, we found that the mass is minimal for a specific number of spikes (Fig. 5d). This could also be seen by plotting the mass as a function of κ for different values of n (Fig. 5e). The value of n at the minimum decreases with κ, with a single spike being minimal at κ → 0. The curves invert so that as κ is increased multiple spikes produce the lowest mass. To test if this critical number of spikes has any relevance to the obtained pattern, we compared the number of spikes predicted by mass minimisation against the distribution of patterns obtained numerically (starting from a small random perturbation around the uniform state). We found remarkable agreement (Fig. 5f). Mass minimisation correctly predicts the most frequent pattern obtained over the entire range of κ, with significant deviation only at the transition points and close to exiting the Turing regime at high κ. In comparison, the linear prediction only agrees for the highest values of κ. Remarkably, the prediction was also reasonably accurate for our default parameter set ( Fig. S3) even though the solution is less spike-like (Fig. 2). Thus, mass minimisation not only predicts where the peaks of a Turing pattern are positioned but also how many peaks there will be.

DISCUSSION
One of the main challenges for the physics of pattern formation is the prediction of which pattern will be obtained, not only at onset, i.e. at entry into the parameter space giving patterns, but generically for any parameter values. While linear stability analysis can give a prediction for the dominant mode, non-linear effects mean that it is often inaccurate. Furthermore, linear analysis cannot explain the periodic nature of final patterns nor (with reflexive boundary conditions) the regular positioning of peaks, which occurs dynamically in several settings that are outside of the linear regime (e.g. domain growth, coarsening, initialised peaks).
Here we have developed a general principle of pattern selection, taking the one-dimensional reaction-diffusion (Turing) model in equation (2) as a model system. We have shown that the movement and positioning of peaks within the obtained Turing patterns is analogous to the behaviour of a diffusive system consisting of point sinks that move with a velocity proportional to the gradient. The flow of mass through such a point-sink system leads to sinks being positioned symmetrically and evenly spaced across the domain (regularly positioned), as this is the unique configuration for which the gradient across each sink vanishes. Importantly, we showed that this configuration also uniquely minimises the total mass of the diffuse species. In the Turing system, peaks (of the slowly diffusing species v) move toward the regularly positioned config-uration, with the same dynamics as point sinks and, in doing so, minimise the total mass of the fast species u.
In the singular limit D v D u , in which the peaks of the Turing pattern become narrow point-like spikes, an analytical approximation showed that u is indeed described by diffusion in the presence of point sinks but where the Dirac delta function terms have a 1/u pre-factor. This 'inverted' sink term leads to the system having a nontrivial dependence on the number of spikes and the rate of mass flow through the system. As a result, there is a critical number of spikes that minimises the mass of u, unlike in the non-inverted system in which more sinks necessarily results in lower mass. As the flow of mass through the system is decreased, the critical number of spikes decreases with a single spike being critical in the mass-conserved limit. Hence, mass minimisation explains the competition or coarsening effect observed in Turing systems. Indeed, we found that the prediction from mass minimisation matches the numerical observations almost perfectly. It also explains the complete coarsening down to a single peak observed in twocomponent mass-conserved systems [19,20]. Our previous three-component mass-conserved model [23], consisting of a two-component Turing system coupled linearly to a third species, did not display such complete coarsening. The role of mass flow provides the explanation. In the three-component mass-conserved system, there is still mass flow through the Turing subsystem and the rate of flow controls the strength of the competition, similar to the current model. Thus mass-minimisation appears to be a fundamental principle behind the behaviour of Turing systems and likely pattern forming systems in general.
Are our results applicable to other systems? If mass minimisation of the fast species is central to Turing patterning, then the mass of the fast species should be minimisable. In our model, the total mass of the system u + v dx at steady state is fixed. In the Gray-Scott model [29] a different positive linear combination of the two species is fixed, while in the Schnackenberg and Brusselator models it is the total mass of the slow species (see supplementary text). What happens then in models that have the mass of the fast species fixed? We consider the following class of systems: Note that at steady state the mass of u is fixed. It is straightforward to show (see supplementary information) that a system of this form cannot admit a Turing instability for any f , consistent with a general principle of mass minimisation. This result holds even if we replace u in the last term by any function g(u) with g (u) > 0, and a correspondingly different measure for the mass of u. More specifically, the outer equation for u obtained in the singular limit D v D u (equation (12)), has the same form as other systems of the substrate-depletion type such as the Schnakenberg [27] and Brusselator [28] and indeed both of these model exhibit the same peak movement towards regular positioning (see supplemental text). Substrate-inhibition models that have peaks of the two species overlapping, such as that of Gierer and Meinhardt [30], also exhibit peak movement towards regular positions. However, the outer equation of these models have a point source term rather than a point sink [31]. The effect of this on mass minimisation remains to be tested.
Our results indicate that both the position and number of peaks of a Turing pattern are positioned so as to minimise the mass of the fast species. There is therefore minimisation with respect to n continuous variables (the peak positions) and with respect to the discrete variable n itself. Might mass minimisation be able to predict which peaks coarsen and when? This remains to be seen but the answer may be connected to the changes in the existence and/or stability of the different solution branches of the differential-algebraic system (17). Finally, unlike flux-balance, mass minimisation extends naturally beyond point sinks and to higher dimensions. It may therefore be useful in the study of more complicated structures such as the stripes, hexagons and spots that appear in two dimensions.

Data availability statement
The data that support the plots within this paper and other findings of this study are available from the corresponding author upon request.

Code availability
All code used to generate the data in this study are available from the corresponding author upon request.

I. LINEAR STABILITY ANALYSIS
We consider the following reaction-diffusion toy-model, It is easy to see from the above equations that total mass at steady state is the same for any set of initial conditions i.e., whereū andv are the steady state concentrations. We implement the following non-dimensionalisation: to obtain in terms of the dimensionless variables, Let us perform the linear stability analysis on the system. In the absence of diffusion there is a single fixed point The Jacobian at this point is given by where f and g are the reaction terms in (2a) and (2b) respectively. The trace and determinant of the Jacobian are easily found to be Since, T rJ < 0 and DetJ > 0, the homogeneous fixed point is always stable in the absence of diffusion for any set of parameters a, b, Γ. Following the standard approach, we then consider a spatial perturbation around the uniform state u = u 0 + δu, v = v 0 + δv. Expanding in terms of eigenfunctions of the Laplacian, we find that the perturbation w = δu δv evolves as w = k c k exp λ k t k where for each k, λ k are the eigenvalues of ΓJ − k 2 D. Hence any mode k (with k = nπ L for reflexive boundary conditions), that has an eigenvalue with positive real part will increase in amplitude i.e. the mode is unstable. In general there is a range of wave numbers k 2 1 < k 2 < k 2 2 having Re(λ(k 2 )) > 0. The condition for this to occur is which we evaluate to find This single inequality relating a, b and d determines the parameter values for which at least one mode k is unstable i.e. the Turing space.
Note that in limit b → 0, the system (2) becomes mass-conserved. We expand the Jacobian J around b = 0 to find The Jacobian (and hence the dispersion relation) becomes independent of b for b << 1. We can thus change b without significantly affecting the linear behaviour of the model.

II. POINT SINKS, FLUX BALANCE AND MASS MINIMISATION
In this section, we prove an important result described in the main text. Consider the introduced stationary equation We divide the equation by δ to obtain where κ = L δ D and µ i = λA(x i ) with λ = µ δ . Consider the Green's function equation defined by We then write solution of eq. (6) as with the same (reflexive) boundary conditions as A. Solving for W we find W (x) = c, so that the expression for A(x) is where the µ i = µ i (x) are determined by the algebraic equations A. The Green's function

The explicit form of the Green's function in equation is
The derivative of G(x; x i ) with respect to x is discontinuous at Note the following property Using this, the flux-differential ∆J i defined in the main text can be written as B. Properties of the Green's function at regular positioning Sinks are regular positioned when they are evenly spaced across the domain as defined bȳ wherex i is the position of i th sink and n is the total number of sinks.

Property I
Evaluating the Green's function at the sink position for regularly positioned sinks defines a symmetric matrix G ij := G(x i ;x j ). Consider the sum of the j column, where a = κ n and the last step follows from the identity, n j=1 cosh(a(j + m)) = csch( a 2 ) sinh( an 2 ) cosh( a 2 (2m + n + 1)) We can similarly define a matrix G + x by evaluating the derivative of the Green's function at regular positioning (G + x ) ij = G x (x + i ;x j ). Summing over the j th column we find using the similar identity, n j=1 sinh(a(j + m)) = csch( a 2 ) sinh( an 2 ) sinh( a 2 (2m + n + 1)).

Property II
Since the summation of G over any of its rows or columns is the same, the vector of 1s,ê, is an eigenvector of G. Evaluating the defining equations for the µ , equation (9), at regular positioning x =x, we obtain the matrix equation Sinceê is an eigenvector of λG + 1, we must have that i.e. all the µ i are identical at regular positioning or in other words, the profile of A is symmetric. We can sum over any row and use (14) to find .
C. Regular positioning and flux balance The flux differential across each sink is given by, We evaluate this expression at regularly positioning, x =x. First we know from equation (17) that all µ j are identical for regularly positioned sinks. Then from equation (15), it follows that immediately that the flux-differentials vanish at regular positioning To show that the regularly positioned configuration is the unique configuration for which the flux-differentials vanish, we perform a power series expansion of ∆J i in κ. It then suffices to show uniqueness for the κ 0 term. We first expand µ i and G (x i ; x j ) For the lowest order terms, we find first that G 0 (x i ; x j ) = 1. Inserting this into the defining equation for the µ i , we have Hence, all the flux-differentials ∆J i vanish uniquely for regularly positioned sinks x i =x i = L n i − L 2 ( 1 n + 1). If we add time-dependence to the system by specifying the sink velocities as being proportional to the their fluxdifferential Then regular positioning is the unique fixed point of the resultant dynamical system (specified by the differentialalgebraic system defined by equations (9), (12) and (21). Given our numerical observation and that the domain is bounded, we assume that this fixed point is stable.

D. Regular positioning and mass minimisation
The total mass (or rather concentration) of A is readily given by integrating equation (8) M We would like to show that the regularly positioned configuration is a stationary point of M i.e. we will show that Using equation (9), we can write Evaluating this expression at regular positioning, and writing C = j G(x i ;x j ) = κ 2 coth( κ 2n ) from equation (14), we obtain We have already seen in equation (17) that all the µ j are identical at regular positioning. Hence, we need only evaluate Inserting the definition of G(x i ; x j ) from equation (10) we have where the last line follows from noting that the summations are the same is in equation (15) but without the i = m term. We have therefore shown that regular positioning is a stationary configuration of the total mass To show that regular positioning is the unique stationary point, we proceed as in the previous section and perform a power series expansion of M , It then suffices to show uniqueness for the first non-trivial order in the expansion. For the Green's function we have We already saw that µ 0i = µ 0 = λ 1+nλ and hence M 0 is a constant. Inserting these into the equation for µ 2i , we obtain The derivative of M 2 is then proportional to which vanishes only for the regularly positioned configuration Hence we have shown that regular positioned sinks is the unique configuration for which the total mass, M , is stationary. Based on our numerical results, we assume that this configuration is generically a minimum. As in Figure 3b and c but where κ is changed via D rather than δ. In this case higher values of κ (b) leads to slower movement is expected by since the lower diffusive length-scale means that there is less flux exiting the system through through the sink (and the total flux through the system (cδL) is unchanged). A. Mass minimisation in the Turing system In the singular (spike) limit of the Turing system (1), we obtain the following equation for the u at steady-state  Figure 5a is reproduced with a changed colorbar for comparison. b. The number of peaks in the pattern minimising the total mass M is plotted as a function of a and κ. c. Contours of a measure; minima of peaks of the steady state patterns in a is plotted as a function a and κ. For, κ 1 the peaks have a high baseline. d. Same as c, but for the measure; peak width of the steady state patterns. For higher values of a, resulting peaks are very broad. e. A direct comparison of the entire Turing space in a and b is plotted. The spike approximation fails(black dots) for κ 1 and higher values a as peaks have a high baseline and are much broader, respectively( see inset). f. The numerically obtained distribution of peak number at steady state for different values of κ (colour scale) in the non-spiky limit overlaid with the prediction of the dominant pattern from linear stability (green triangles) and the prediction from mass minimisation (red circles). Mass minimisation correctly predicts the number of peaks at steady state, even though it is less accurate than the spike case in Figure 5f. Data from 50 simulations for each parameter set. Parameters: Default with L = 4. For f we sweep across κ (black line in a and b) for the default value of a = 3.75.
as discussed in the main text. Note the inverse dependence on the variable in the sink term. Following the same approach as in the previous section, the solution of this equation is given by with the Green's function as previously defined (with the equivalent dimensionless parameter κ = L δ D ) and where the ρ i are determined by the n algebraic equations where σ = ρ c 2 δ is a dimensionless parameter. Let us consider the case where the sinks are regularly positioned and of the same strength i.e. x =x and ρ i = ρ . We then have Solving for ρ , we find two solutions with corresponding total masses For a fixed set parameters, we can thus calculate the predicted total mass of u for a solution consisting of n regularly positioned spikes. In the main text, we show that, taking the ρ + solution, the value of n that minimizes M is an excellent predictor for the steady-state number of peaks obtained numerically (after any coarsening) even when the peaks are not very spike-like.

Single peak
Let us consider the general case (arbitrarily positioned) of a single spike. The solution is given by where, for a given spike position x 1 , ρ is determined from G(x 1 ; x 1 ) = κ 2 cosh( 2κx1 L ) sinh(κ) + coth(κ) The flux differential across the spike is given by Taylor expanding around κ = 0 we obtain Hence, we find a similar linear dependence of the flux-differential on the sink position as for the 'non-inverted' case discussed earlier.
B. M can not be fixed for all steady state solutions In the Turing system considered here, the total combined mass of the system at steady state is fixed Other systems such as the Brusselator and the Schnakenberg model have the integral ofv being fixed. In our model, the individual integrals ofū andv can vary and we have shown that the steady-sate reached by the system, after any coarsening, is determined by minimizing the integral ofū, the fast species. Since our results indicate that minimization of the mass (integral) of the fast species is the determinant of pattern selection, it suggests that a Turing patterns may not be possible in a model that has the mass of u fixed at steady state.
Let us consider a general class of Turing systems of this form with D v < D u . The functions f (u, v) and g(u) are arbitrary apart from the constraint that g (u) > 0. g(·)dx then acts a measure for u and all steady states have the same mass of u using this measure 1 L denote a stable fixed point of the homogeneous system by (u 0 , v 0 ). The corresponding Jacobian is given by, If the fixed point is stable we must have The latter relation implies that f v < 0. A necessary condition for a Turing instability is that However since D v < D u , this condition can never be satisfied. Hence, the general system (35) does not admit a Turing instability is consistent with our results mass minimisation of fast species underlies pattern selection.

C. Numerical Methods
The simulations were performed in a spatial lattice x ∈ [− L 2 , L 2 ] and time domain t ∈ [0, T ], where L is the length of the spatial domain and T , the total time. The MATLAB solver pdepe was used to solve the coupled partial differential equations. The simulations were performed with the following default parametric values (unless explicitly stated otherwise): The relative and absolute tolerances in the difference between two values of iteration were set to 10 −6 and 10 −12 respectively.
where r 1 , r 2 are vectors of random numbers. We now present an analysis of other reaction-diffusion models and show numerically that most of the features of the exploratory model in eq. (1) still holds. The number of peaks in the most frequent steady-state pattern is plotted as a function of a and κ. For each set of parameters, the most frequent pattern was obtained from 5 simulations each initialised with a different random perturbation from the uniform state. The simulations were ran for long enough to ensure the steady-state pattern was reached. b. The number of peaks in the mode with the greatest growth rate as predicted by linear stability analysis is plotted as a function of a and κ. Plots a and b are similar for κ 1. They disagree for κ ≤ 1, which indicates coarsening.
The general version of Brusselator is described by the following equations, The Brusselator model has a similiar form as the exploratory model in eq. (1). In the absence of diffusion it has a single fixed point. The model also has the form of a mass-conserving Turing system with additional linear terms. However, note that rather than total mass being fixed at steady state, it is the mass of v that is fixed