## Abstract

Herein we present a framework to characterize different sources of protein expression variability in Turing patterned tissues. In this context, we introduce the concept of granular noise to account for the unavoidable fluctuations due to finite cell-size effects and show that the nearest-neighbours autocorrelation function provides the means to measure it. To test our findings, we perform *in silico* experiments of growing tissues driven by a generic activator–inhibitor dynamics. Our results show that the relative importance of different sources of noise depends on the ratio between the characteristic size of cells and that of the pattern domains and on the ratio between the pattern amplitude and the effective intensity of the biochemical fluctuations. Importantly, our framework provides the tools to measure and distinguish different stochastic contributions during patterning: granularity versus biochemical noise. In addition, our analysis identifies the protein species that buffer the stochasticity the best and, consequently, it can help to determine key instructive signals in systems driven by a Turing instability. Altogether, we expect our study to be relevant in developmental processes leading to the formation of periodic patterns in tissues.

## 1. Introduction

Tissue patterning sets the developmental roadmap that provides positional information to cells and confers their unique identities [1–4]. Thus, as a response to distinct protein expression levels in a primordium, cells commit to different fates, may undergo apoptosis, increase/decrease the proliferation rate, or change their division mode to ultimately put into action the developmental plan that shapes the organism. As for the mechanisms of pattern formation, short-range signalling elicits local responses that may propagate in the tissue, e.g. Notch–Delta interactions leading to lateral inhibition/induction [5,6]. On the other hand, long-range signalling driven by diffusive molecules, morphogens, enables positional information to cells at a larger, tissue-level scale and induces pattern formation by following either the French flag model or other mechanisms [4,7]. In particular, Alan Turing proposed a long-range patterning mechanism where interacting species (activators and inhibitors) with different diffusive properties lead to periodic protein expression profiles [8,9]. Turing's proposal meant a breakthrough in mathematical biology; however, besides its success in explaining animal coating [9–11], its relevance in the field of development remained elusive for many years. Notably, during the last decade a number of developmental structures have been found to be shaped by this mechanism [12–14].

The characteristic size of the domains patterned by diffusive processes is controlled by the physical constraints set by the size of morphogen molecules, their production and degradation rates, and the time scales associated with exocytosis and endocytosis among other factors [15,16]. As a result, the reported ratio between the size of patterned domains and that of cells is typically of the order of [12,14]. Moreover, cell growth leads to size doubling that leads to perturbations of this ratio by around 20%. These facts raise the intriguing question of how cells deal with aliasing-like effects. By ‘aliasing-like effects’, we mean the consequences in terms of protein variability derived from the sampling process since, as illustrated in figure 1, cells must perform an analogue-to-digital conversion of the, otherwise idealized, continuous periodic pattern into a discretized version due to finite cell-size effects. Here we coin the term *granular noise* to describe such variability. We point out that the granular effect relies on the hypothesis that proteins are well mixed at the single-cell level and, hence, cells cannot ‘sample’ protein numbers with a subcellular resolution. This hypothesis is sustained by the fact that patterning depends ultimately on morphogens that are driven by a homogenizing transport mechanism: diffusion.

During the last decade, much progress has been attained in the field of noise in gene expression at the single-cell level [17–26]. On the other hand, a number of studies have pointed out the ordering role of fluctuations in spatially extended systems [27–34]. In the context of developmental processes, the role played by cell-to-cell variability, modelled as biochemical noise, in morphogen patterned systems to produce robust and precise positional information has been explored [35]. Also, the effects of embryo size variability, the interactions that help to attenuate the associated noise and the scaling properties of positional information in growing tissues have been recently addressed [36,37].

However, the identification, the characterization and the function of different sources of variability at the collective tissue level within a biological context remain barely explored [38–41]. Notably, to the best of our knowledge, granularity effects have been neglected altogether. Moreover, so far there are no tools that allow us to distinguish among different stochastic contributions, i.e. to measure independently the variability arising from biochemical noise and that from the, unavoidable, granularity. Major difficulties in performing such analysis include the competing effects of various sources of variability at different spatio-temporal scales, the dynamical character of growing tissues and the role played by cell communication that effectively propagates noise in the local cellular environment. Here, we address these questions and present a framework to study gene expression variability, and differentiate among various contributions, in Turing-like patterned tissues (i.e. patterns showing a periodic motif). To that end, we introduce a formalism based on the autocorrelation function that can be easily implemented when processing experimental data. We illustrate our findings by means of a generic activator–inhibitor system and perform realistic numerical simulations of growing tissues using a vertex-model approach that includes mechanical interactions between cells, cell-cycle variability, binomial partition of molecules between daughter cells following division and biochemical noise. Our results confirm that the autocorrelation function is a robust method to compute the granular noise and to estimate also the levels of biochemical fluctuations. The applicability of our methodology to specific developmental processes patterned by the Turing instability is straightforward and we provide simple guidelines to assess and quantify the importance of granularity versus other protein variability sources.

## 2. Results

As for the various contributions to cell-to-cell variability, finite cell-size effects in protein variability are relevant at length scales of the order of *Ω* (cell's size) and set protein concentrations *quanta*: a characteristic protein concentration difference between neighbouring cells that results in a discontinuous protein concentration profile (figure 1). In addition, at the single-cell level, the effects of the biochemical noise further contribute to increasing the variability by a factor *σ*_{z} (the effective intensity of the biochemical fluctuations). Finally, patterning leads to protein concentration differences of order *δ*_{z} (pattern's amplitude) over length scales of order *λ* (pattern's wavelength).

To filter out the influence of the global patterning on the concentration variability from our analyses, we make use of the nearest-neighbours autocorrelation function
2.1where *z*_{n} stands for the protein concentration levels of species *Z* at cell *n* and the average, 〈•〉, is performed over cells (see below for considerations about time averaging). By estimating the concentration differences between neighbouring cells, *G*_{z} effectively measures the mutual information in the local cellular environment. Note that *G*_{z} also provides an (indirect) estimation of the value of the protein quanta, Δ*z* = [〈(*z*_{n+1} − *z*_{n})^{2}〉]^{1/2}/〈*z*_{n}〉.

Note that if there is no biochemical noise, either in the case of constant functions (no patterning) or in the case that there is patterning but it is a continuous solution, , then *z*_{n+1}∼ *z*_{n} and consequently *G*_{z} ∼ 0 (maximal similarity in the local cellular environment). On the other hand, when biochemical fluctuations are dominant, *z*_{n} behaves as a spatial white noise, 〈*z*_{n}*z*_{m}〉 ∼ *δ*_{nm}/*Ω* (*δ*_{nm} being the Kronecker delta) and *G*_{z} ∼ 1 (minimal similarity among neighbouring cells).

If further sampling over configurations, i.e. noise realizations, is needed (see below) then 〈•〉 accounts for a time average too. Ideally, the sampling frequency for time averages must be smaller than the typical protein turnover rate (to ensure that the pattern has reached stationary conditions) but larger than the characteristic frequency for pattern remodelling (to render enough statistics). Alternatively, if the focus of the study is to characterize the non-equilibrium fluctuations during the process of pattern formation, the analysis should be performed using a sampling frequency larger than the typical protein turnover rate. As for pattern remodelling, as the size of the tissue increases due to cellular growth more of the pattern's domains, i.e. wavelengths, can fit into the size-increasing primordium, thus leading to rearrangements (see the electronic supplementary material, movies S1–S4). Assuming exponential growth conditions, remodelling occurs at a rate ∼ (*τ*log(*λ*/*Ω*))^{−1} : *τ* being the average cell-cycle duration. Notably, pattern remodelling relies on the ability of cells to change their expression profile depending on their local environment as experimentally reported [14].

To separate the effects arising from different sources, we analyse first the case when the biochemical noise is negligible and show that we can provide an accurate theoretical estimation of the background granularity and its confidence bounds. If a periodic pattern develops, then *z*_{n} can be approximated as *z*_{n} = *z*_{0} + *δ*_{z}cos (*q***nΩ*) along, at least, one of the symmetry axes, where *z*_{0} is the average protein concentration and *q** = 2*π*/*λ* is the most unstable Fourier mode. Under these conditions, *G*_{z} can be estimated (see Material and methods)
2.2where *Ω* is the characteristic mean cell size and *σ*_{Ω} is its standard deviation. The error band of *G*_{z} sets the confidence bounds of the theoretical prediction and helps to (i) spot pattern misalignment with sampling boxes and (ii) assess if granularity is the dominant source of protein variability. Expression (2.2) determines the, unavoidable, background variability among cells in a local environment due to finite-size effects: the so-called *granular* noise. Note that this expression is independent of the expression levels of protein species and just depends on geometrical constraints: the ratio between the pattern wavelength and the cell size. Importantly, these quantities, together with *σ*_{Ω}, can be characterized and, consequently, one can compute the basal value of the autocorrelation in the absence of biochemical noise, *G*_{z}|_{σ2z = 0}. As shown below, the relative difference of the estimated autocorrelation, *G*_{z}|_{σ2z = 0}, from its actual, measured, value in the presence of noise, *G*_{z}|_{σ2z≠0}, allows us to distinguish between stochastic sources of protein variability.

In regards of the protein variability due to the biochemical noise, in reaction–diffusion systems the *effective* magnitude of the fluctuations depends on the regulatory interactions between species, which can either enhance or damp out the noise, and on the diffusive process, which propagates the noise among cells in their local environment and also contributes to average out the fluctuations. We consider that the overall effect of this source of noise can be described as a random additive contribution, *η*^{z}_{n}, to the protein concentration: *z*_{n} = *z*_{0} + *δ*_{z}cos(*q***nΩ*) + *η*^{z}_{n}. Thus, we expect the fluctuations to be more relevant at locations where cos(*q***nΩ*) ≃ 0 and the ratio between the noise and the pattern solution is maximal: the domain boundaries of the pattern. If there is enough statistics then either spatial or time averages, assuming ergodicity, sample the pattern configuration space (noise realizations): . In that case (see Material and methods)
2.3Consequently, the *noise-to-signal* ratio can be written as
2.4where *G*_{z}|_{σ2z = 0} is the background granular noise, which can be estimated by measuring *λ* and *Ω* by using equation (2.2), and *G*_{z}|_{σ2z≠0} is the actual nearest-neighbour autocorrelation measured in the presence of biochemical noise using equation (2.1). Table 1 summarizes the most relevant parameters defined in our analysis to estimate the granular noise and measure the effective value of the biochemical fluctuations.

To test and illustrate our findings, we perform computer simulations of growing tissues driven by a generic activator–inhibitor protein dynamics (see Material and methods). Figure 2 shows the regulatory interactions between species and the region where a Turing instability develops and patterns the growing tissue. We stress that the formalism introduced herein does not depend on the specific details of the model (e.g. the number of species or the functional form of the equations) but simply on the existence of a patterning solution with a well-defined periodicity regardless of the mechanism. We simulate the regulatory interactions, the cell's biomechanics and growth/division effects (including a binomial redistribution of proteins between daughter cells following a division) by means of a vertex model approach (see Material and methods). We consider parameter sets with (i) distinct aspect ratios between the pattern domain size and the cell size and (ii) different pattern amplitudes (figure 2*b*).

In order to assess whether our theoretical estimation of *G*_{z}|_{σ2z = 0} captures correctly the background granularity, we first performed numerical simulations of the growing tissue in the absence of biochemical fluctuations. Figure 3 shows snapshots of the growing, patterned tissue for different *Ω*/*λ* ratios (see the electronic supplementary material, movies S1 and S3). To quantify *G*_{z} using equation (2.1), we sampled the tissue at regular time intervals (approx. 40 frames per cell cycle). We show the importance of sampling the periodicity of the patterned solution correctly by sampling cells from the tissue along perpendicular, fixed stripes with a large aspect ratio: approximately one cell diameter width and whole-tissue length (figure 3*b*). Note that in the case of stripe-like patterns it is possible to find a large degree of alignment of the pattern with a sampling box (thus masking the pattern periodicity). In our analysis, the cell's characteristic diameter is estimated as the square root of the apical area of the sampled cells.

As shown in figure 3, the quantification of the background granular noise is in agreement with the theoretical estimation, equation (2.2). As expected, finite cell-size effects are more pronounced as the ratio *Ω*/*λ* increases, that is, as the ratio of the pattern wavelength to the cell size decreases. Note also that the granular noise has an *extrinsic* character [19] since, as predicted, it affects all protein species equally and independently of their expression levels, i.e. *G*_{u} ≃ *G*_{v}. Thus, in a periodic patterning situation, if the pattern wavelength and the characteristic cell size are determined, then equation (2.2) accounts for the variability in protein concentration due to finite cell-size effects and allows us to distinguish among different noisy sources as shown below.

In order to explore the competing effects between granular noise and intrinsic/extrinsic fluctuations, we implement a stochastic dynamics for proteins *U* and *V* and study different noise-to-signal ratios (see Material and methods). Figure 4 shows the quantification of *G*_{z}|_{σ2z}_{≠0} as measured by equation (2.1), the theoretical estimations of *G*_{z}|_{σ2z = 0}, the noise-to-signal ratio, and the patterning domains of *U* and *V* (see the electronic supplementary material, movies S2 and S4). As the levels of *G*_{z} indicate, the protein concentration profiles of the activator (slow-diffusive) species, *U*, are noisier than those of *V* and these stochastic effects are specially relevant at the domain boundaries. In addition, the effects due to fluctuations are more significant as the background granular noise decreases. These results reveal that the noise-to-signal ratio is dominant in (i) slow versus fast diffusing protein species, because noise averages out more easily in the second case, (ii) low versus high background granular noise situations, because otherwise the latter prevails, and (iii) activator versus inhibitor species, because the positive auto-regulation of the former contributes to amplify the noise. Note also that *G*_{z} deviates from the extrinsic behaviour as predicted by equation (2.3).

## 3. Discussion

We have shown that by using the nearest-neighbours autocorrelation function the variability due to global patterning effects can be filtered out. Such variability of order *δ*_{z} can be estimated by Fourier analysis (amplitude of the most unstable mode). Aliasing effects and biochemical noisy contributions (intrinsic and/or extrinsic fluctuations) can be separated and determined by means of *G*_{z} when combining theoretical expressions and measurements. As for the applicability of our approach, experimental results have revealed that the ratio *λ*/*Ω*, pattern wavelength versus cell size, can be small and patterning domains may comprise around 8–10 cells [12–14]. Thus, we predict that the granular noise contributes significantly to protein variability. In the particular case of the limb bud, recent results have shown that this primordium is patterned by a Turing instability resulting from the interactions between two morphogens, *Bmp* and *Wnt*, and the skeletal marker *Sox*9 [14]. A quick assessment indicates that *λ* ∼ 10*Ω* and consequently the background granularity is of the order of around 20%. Evidence supporting the buffer effect of diffusion over the fluctuations can also be found in the limb bud. Both morphogens, *Bmp* and *Wnt*, lack a self-regulatory motif that in turn helps to damp the biochemical fluctuations, as shown in our simulations (see below for further considerations about the instructive role of inhibitors in Turing systems). However, the values of the diffusion coefficients satisfy that *D*_{Bmp} ≫ *D*_{Wnt} and one would expect, according to our analysis, that *Wnt* has a noiser profile than *Bmp*. The experimentally reported profiles of these species during limb bud patterning [14] seem to be in agreement with this prediction.

Importantly, we have presented our results using a methodology that is feasible in experimental situations. In particular, since the dimensionality of the pattern is not relevant for estimating these effects, by determining the pattern wavelength and by sampling the size and protein expression levels for (i) a few cells, , along a direction that captures the best the pattern periodicity and (ii) a small number of time points in the tissue dynamics, , then the background granularity and the biochemical noise can be specified. Further sampling can provide better, more accurate results but the fact that the method works with a limited statistics reveals its robustness. We acknowledge that the experimental quantification of protein concentrations at the single-cell level requires an accurate calibration of the fluorescence intensity [42]. However, we stress that the background granularity levels, as prescribed by the theoretical expression equation (2.2), are independent of the protein levels. Thus, as long as the fluorescence levels of protein species are normalized, our method provides a way to compute the relative noise-to-signal ratio and weight the importance of biochemical noise with respect to granularity.

The proposed method relies on several assumptions that are worth commenting on. First, we have assumed that a pattern with a well-defined wavelength develops. In any particular patterning situation, a Fourier analysis on the protein profile, as obtained, for example, by fluorescence microscopy, will reveal the validity of this assumption. We point out that, even if several Fourier modes are relevant, the methodology remains valid. However, the analytical expressions to estimate the background granularity, equation (2.2), and the noise-to-signal ratio, equation (2.4), would need revision. In any case, if the pattern solution is a collection of waves then the calculations can be carried out easily (Material and methods). Second, we have assumed that, besides the pattern's rearrangement dynamics due to cell growth, stationary conditions apply. In this regard, we notice that if the amplitude of the pattern is modulated in time then our results still hold. However, we stress that time averaging must be done carefully to ensure that the pattern configuration space is properly sampled. Third, we have considered that the overall effect of the biochemical fluctuations renders a Gaussian-like, white behaviour that is independent of the value of the pattern amplitude. This approximation is motivated by the central limit theorem [43] and makes it possible to obtain analytical expressions for the noise-to-signal ratio. However, it may raise questions about its fidelity since other statistical properties for the noise are certainly possible [21,41]. Yet, the background granular noise is independent of the statistical properties of the noise since it merely depends on the geometrical constraints of the problem (cell size and pattern wavelength). Consequently, the difference between the nearest-neighbours autocorrelation and the estimated theoretical background granular noise still provides an estimation of the levels of biochemical fluctuations.

As for additional implications of our study, our results indicate that variability due to the biochemical noise is larger for cells at domain boundaries, especially for activator species. According to this, fate commitment, as a result of patterning, would be more difficult for those cells. In addition, boundary lines among cell populations would be wiggly, thus suggesting that the formation of developmental structures, e.g. the fingering pattern in the case of the limb bud [14], could be challenged. However, the robustness displayed by biological systems argues against these ideas. This raises the intriguing question about the mechanisms that are able to buffer noisy effects in tissue patterning. In this regard, we note that the level of variability is reduced in inhibitor species (figure 4). The latter suggests a possible instructive role for this species in fate decision-making (less noisy). In addition, it may explain the recurrent feedback found in developmental boundaries between patterning and cell mechanical properties to refine compartmentalization [44–49].

We expect our study to be relevant in developmental processes leading to the formation of periodic patterns in tissues when cell-to-cell variability needs to be characterized to, for example, better understand fate decision-making. The extension of the proposed method to other patterning situations, e.g. morphogen gradient profiles, or the implementation of feedback between signalling and cell mechanics to buffer the fluctuation effects are promising avenues of research that can shed further light on the role of noise in developmental patterning systems. Work in these directions is in progress.

## 4. Material and methods

### 4.1. Tissue simulations

The tissue dynamics is implemented in our simulations using a vertex model approach [50]. Each cell in the tissue is represented by a discrete set of points: the vertexes that define its shape. The energy associated with a vertex *i* reads
where the sums indexed by *α* and 〈*ij*〉 run, respectively, over the cells *α* and vertices *j* sharing vertex *i*. is the cell apical area, *K*_{α} is proportional to the Young modulus, *Λ*_{ij} is a line tension that weights the cell adhesion effects (*l*_{ij} being the length of the edge connecting neighbouring vertices *i* and *j*), *Γ*_{α} is a term that accounts for the contraction effect of the actomyosin cortical ring and *L*_{α} is the cell perimeter (see [44] for details). The parameters used in our simulations are (dimensionless): , *K* = 1, *Γ* = 2 × 10^{−2} and *Λ* = 10^{−2}. In addition we consider that in the tissue periphery *Λ* = 5 × 10^{−2}. The cell-cycle duration, *τ*, is stochastic in our simulations and follows the rule *τ* = *ɛt*_{det} + (1 − *ɛ*) *t*_{sto}, where *t*_{det} is a deterministic time scale that accounts for a mean cell-cycle duration in the absence of mechanical stress and *t*_{sto} is a random variable that accounts for the variability of cell-cycle duration, . The parameter *ɛ* ∈ [0, 1] weights the stochasticity of the cell-cycle duration. In our simulations, *ɛ* = 0.8.

Cellular growth is implemented by prescribing the following dynamics for the *preferred* apical cell area (dimensionless): cells start to grow their apical area linearly (towards doubling) at the moment they reach the middle of the cell-cycle progression [44]. With respect to the cleavage orientation, the code evaluates the inertia tensor of the cell with respect to its centre of mass assuming that a proper representation of the former is a polygonal set of rods, i.e. the cell edges. The principal inertia axes indicate the symmetry axes of the cell: the longest axis of the cell is orthogonal to the largest principal inertia axis. Cells divide following the Hertwig rule: the cleavage plane is perpendicular to the longest axis [51]. Cleavage is assumed to happen instantaneously. In our *in silico* experiments, tissues are typically grown from approximately 2.5 × 10^{3} to approximately 1.2 × 10^{4} cells (approx. two cell cycles). Transients effects in the tissue dynamics due to the initial conditions are eliminated from our analysis by discarding data from the first cell cycle.

As for the protein dynamics, we assume each cell to be a well-stirred system, where spatial effects can be disregarded. Protein concentration values in each cell are obtained by integrating numerically the discretized model equations using an Euler algorithm. We take into account cell growth (dilution effects) to determine protein concentrations. Also, proteins are distributed binomially between daughter cells as a consequence of a division event. The morphogen diffusion process is mimicked by an out-of-lattice, discretized Laplacian operator that conserves the number of molecules [52],
where *Z*_{i} stands for the number of proteins at cell *i*, is its area, the sum runs over its nearest neighbours *j*, *L*_{ij} is the length of the membrane shared among cells and *r*_{ij} is the distance between cell centres. This definition captures rigorously the discretization of the Laplacian operator as long as the well-mixed hypothesis holds (see Introduction).

Biochemical noise effects in protein species *Z* at cell *n* and time *t*, *ξ*^{z}_{n}(*t*), are implemented by means of additive uncorrelated Gaussian fluctuations: 〈*ξ*^{z}_{n}(*t*)〉 = 0 and 〈*ξ*^{z}_{n}(*t*)*ξ*^{z′}_{m}(*t*′)〉 = (*σ*^{2}_{ξ}/*Ω*)*δ*_{nm}*δ*_{zz′}*δ*(*t* − *t*′), where *σ*^{2}_{ξ} is the intensity of the noise. Our simulations explore the conditions *σ*_{ξ} < *δ*_{z} (black circle in figure 1*b*) and *σ*_{ξ} ∼ *δ*_{z} (white circle in figure 1*b*). We do not consider the case *σ*_{ξ} > *δ*_{z} since it leads to unrealistic fluctuation-controlled expression profiles in which patterning does not play a key role. In all cases *σ*_{ξ} = 10^{−1}. We point out that the value of *σ*_{ξ} does not represent the *effective* intensity of the fluctuations *σ*_{z} since the latter depends on the regulatory interactions and on the morphogen diffusion process. That is, *ξ*_{n} are the *input* fluctuations with intensity *σ*_{ξ} that we use to represent the biochemical noise at the single-cell level in the simulations, and *η*_{n} are the *output* fluctuations with intensity *σ*_{z} that, as a result of the interactions and diffusion, the proteins effectively experience and we register using the autocorrelation method. The code for carrying out the numerical simulations is provided as the electronic supplementary material.

### 4.2. Analytical calculations

Given the expression of the autocorrelation function, equation (2.2), and taking into account that 〈*z*_{n+1}〉 = 〈*z*_{n}〉, then *G*_{z} can be rewritten as
4.1Thus, to characterize *G*_{z} the moments 〈*z*_{n}〉, 〈*z*^{2}_{n}〉 and 〈*z*_{n+1}*z*_{n}〉 need to be estimated. Given the functional form of the Turing patterning solution, *z*_{n} = *z*_{0} + *δ*_{z} cos(*q***nΩ*), and the following identities (valid if *N* = *L*/*Ω* ≫ 1, *L* being the size of the tissue), 〈cos(*q***nΩ*)〉 ≃ 0, and 〈cos(*q***nΩ*)cos(*q**(*n* + 1)*Ω*)〉 ≃ cos(*q***Ω*)/2, those moments read
4.2
4.3
4.4Substituting the above expressions into equation (4.1), we obtain *G*_{z} = 1 − cos(*q***Ω*). The confidence bounds for *G*_{z} are determined by the variation of the parameters *q** and *Ω*. Yet, Turing patterning tissues show a precise wavelength and the main source of variability is the cell size, which may show variations up to 100% due to size doubling during cell-cycle progression. Thus, given the standard deviation of the cell size in the tissue, *σ*_{Ω} , the error propagates to *G*_{z} as ± |∂*G*_{z}/∂*Ω*|*σ*_{Ω} = ± *σ*_{Ω}*q**|sin(*q***Ω*)|.

In the case where we assume a random perturbation, *η*, of the patterning solution due to the biochemical noise, the functional form of the protein concentration profile now reads
4.5By invoking the central limit theorem [43], here we implement a Gaussian approximation and assume that *η*^{z}_{n} satisfies the distribution
where is a normalization constant and *σ*^{2}_{z}/*Ω* is the *effective* intensity of the biochemical noise of species *Z*. In addition, we consider that the fluctuations are uncorrelated in space and time: 〈*η*_{n}*η*_{m}〉 = *σ*^{2}_{z}(*δ*_{nm}/2*Ω*)*δ*(*t* − *t*′), where *δ*_{ij} stands for the Kronecker delta and *δ*(*s*) for the Dirac delta. Note that the cell size, *Ω*, is required in the definition of the noise autocorrelation such that in the continuous limit the autocorrelation satisfies a white noise character in space: . Formally, one may choose a definition of the noise intensity leaving out the cell size *Ω* and integrating the latter in the correlation properties of the fluctuations (i.e. in the delta function). Yet, herein we preferred to explicitly include this dependence of the cell size effectively in the definition of the biochemical fluctuations. The Gaussian approximation is valid as long as the intensity of the noise does not overcome the patterning solution, that is, *σ*_{z} < *z*_{0} − *δ*_{z}. Otherwise the noise could lead to unphysical values (negative) of the protein concentration. Given *ρ*(*η*^{z}_{n}) the probability distribution for *z*_{n} reads
4.6In this case, the moments 〈*z*_{n}〉, 〈*z*^{2}_{n}〉 and 〈*z*_{n+1}*z*_{n}〉, by averaging over noise realizations, read
4.7
4.8
4.9Consequently,
4.10and
4.11

We notice that in the case where the assessment of the patterning solution reveals more than one unstable Fourier mode, e.g. *z*_{n} = *z*_{0} + *δ*^{(1)}_{z}cos(*q**_{(1)}*nΩ*) + *δ*^{(2)}_{z}cos(*q**_{(2)}*nΩ*), the above calculations are straightforward. The generalization of our framework to other patterning situations, e.g. morphogen gradients, can also be easily implemented. By taking into account the morphogen decay length, *λ*, the protein profile reads [15], *z*_{n} = *z*_{0} e^{−nΩ/λ}. Thus, the moments 〈*z*_{n}〉, 〈*z*^{2}_{n}〉, and 〈*z*_{n+1}*z*_{n}〉 can be estimated and the background granularity reads
As in the case of Turing patterning, the background granularity goes to zero as the ratio *Ω*/*λ* does.

### 4.3. Activator–inhibitor Turing patterning systems

For the sake of simplicity, we restrict our analysis to a system with two coupled reaction–diffusion equations in one spatial dimension. Yet, the analysis presented below can be easily generalized. Thus, the following dimensionless system of equations describes the reaction and diffusion terms of two protein species *U* and *V*:
4.12where *u* = *u*(*x*, *t*) and *v* = *v*(*x*, *t*) represent protein concentrations. The reaction terms *f* and *g* are supposed to have a single equilibrium point, , such that *f*(*u*^{0}, *v*^{0}) = *g*(*u*^{0}, *v*^{0}) = 0 and the point is assumed to be stable in the absence of diffusion. That is, by defining and , where *z* stands for either the field *u* or the field *v*, the following conditions hold:
4.13A Fourier analysis reveals that a mode, *q**≠0, becomes destabilized and a pattern develops if [9,33]. It is easy to prove that in terms of the signs of the entries of the Jacobian matrix
only four out of the possible eight options can lead to a Turing pattern instability:
Here we choose the first option where species *U* stands for an activator and *V* for an inhibitor: {*f*_{u}, *f*_{v}} > 0, {*g*_{u}, *g*_{v}} < 0. If *f*_{u} < |*g*_{v}|, and *f*_{u}|*g*_{v}| < *f*_{v}|*g*_{u}|, then conditions (4.11) are satisfied. In addition, if *f*_{u} = *f*_{v} = *a*, *g*_{u} = − 2 and *g*_{v} = − 1 a pattern develops if
The modelling equations require nonlinear terms to saturate to a finite value the amplitude of the pattern, which otherwise diverges. By choosing a cubic non-linearity and , we obtain the following modelling equations describing an activator–inhibitor system:
4.14

## Data accessibility

The simulation code and compilation instructions for generating all data shown in the manuscript are available at https://doi.org/10.6084/m9.figshare.5146564.

## Competing interests

I declare I have no competing interests.

## Funding

No funding has been received for this article.

## Footnotes

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.3852112.

- Received May 2, 2017.
- Accepted August 2, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.