The impact of cell size on morphogen gradient precision

ABSTRACT Tissue patterning during embryonic development is remarkably precise. Here, we numerically determine the impact of the cell diameter, gradient length and the morphogen source on the variability of morphogen gradients. We show that the positional error increases with the gradient length relative to the size of the morphogen source, and with the square root of the cell diameter and the readout position. We provide theoretical explanations for these relationships, and show that they enable high patterning precision over developmental time for readouts that scale with expanding tissue domains, as observed in the Drosophila wing disc. Our analysis suggests that epithelial tissues generally achieve higher patterning precision with small cross-sectional cell areas. An extensive survey of measured apical cell areas shows that they are indeed small in developing tissues that are patterned by morphogen gradients. Enhanced precision may thus have led to the emergence of pseudostratification in epithelia, a phenomenon for which the evolutionary benefit had so far remained elusive.

In this supplementary document, we theoretically show that averaging morphogen concentrations over a spatial region (such as cell areas) can shift the effective readout position compared to point-like readout, and we derive the corresponding shift ∆x analytically for isotropic and for rectangular cell shapes. We focus on exponential morphogen gradients here as they arise in systems with diffusion-driven morphogen transport and uniform linear degradation, but note that the developed formalism can be applied directly also to other gradient shapes. Moreover, the impact of spatial correlation of the kinetic cell parameters on the positional error, the choice of the kinetic parameter distribution and the effect of cell number in the source domain are discussed.

Readout in a continuous domain
Consider an exponential morphogen concentration gradient with concentration C0 at the source at x = 0, and decay length λ. Assuming a continuous readout based on a threshold concentration C θ = C(x θ ), a positional identity boundary forms at position This mechanism allows for gradient-based tissue patterning, where individual patterning domains are delineated by different boundary positions x θ resulting from different readout thresholds C θ .

Readout in a tissue of isotropic cells
For morphogen readout in a cellular tissue, we consider several different cases in a unified description. Cells can either sense the morphogen concentration at a singular point, averaged over a spatial region with radius r about that point (which may or may not be smaller than a cell), or as an average concentration over the entire cell area. We denote this readout region by Ω (Fig. S1). The average concentration in Ω is Assuming that the averaging domain is circular (i.e., the cell areas have no orientational bias) in a two-dimensional tissue cross section or surface, we can approximate Ω as a disk with radius r about a center point (x0, 0): . Averaging an exponential morphogen concentration (blue) over a local region such as the cell area leads to a larger readout concentration (green) than taking the concentration at the middle of the region (red). To compensate for this effect, the readout position shifts downhill (away from the source) by a distance ∆x from x θ to x0.
In the case where the concentration is averaged over the entire cell area, r is the effective cell radius. The average concentration thus becomes is the modified Bessel function of the first kind for integer n.
The series converges very quickly if r ≪ λ, such that higher order terms in r/λ can be dropped. Substitution and expansion of the Bessel function yields Thus, the mean concentration ⟨C⟩ is larger than the one in the middle of the readout domain, C(x0), and this deviation increases with larger readout regions and shorter gradient decay lengths. If threshold-based readout operates on the averaged concentration, we must have C θ = ⟨C⟩. Therefore, The location of domain boundaries is shifted down the concentration gradient by the distance as shown in Fig. S1. Notably, the shift is independent of both the gradient amplitude C0 and the concentration threshold C θ for an exponential gradient. Therefore, it is the same for all readout positions in the pattern if the averaging radius r and the decay length λ are spatially invariant, such that all domain boundaries are shifted equally by this averaging effect. Eq. S2 is plotted in Fig. S2. Using the power series expansion of the natural logarithm, the boundary shift can be expanded to For a mean cell radius of r = 2.5 µm and a gradient decay length of λ = 20 µm, the shift is ∆x ≈ 0.039 µm.
By combining Eqs. S1 and S2, we find the mean domain boundary position at

Readout in a tissue of rectangular cells
We now derive the downhill shift ∆x also for rectangular cell areas, effectively rendering the problem one-dimensional. This scenario corresponds to a tissue composed of cuboidal cells in which the morphogen gradient forms in a direction perpendicular to one of the cells' axes. In this case, Averaging over the cell area thus gives Requiring again that the readout threshold be the average concentration, C θ = ⟨C⟩, yields The shift in the readout position then follows as Eq. S4 is plotted in Fig. S2. For a mean cell radius of r = 2.5 µm (which in this case corresponds to the half-width of the rectangular cells) and a gradient decay length of λ = 20 µm, the shift is ∆x ≈ 0.052 µm.
In analogy to Eq. S3, the mean domain boundary position is found at in tissues composed of rectangular cells.

Impact of spatial correlation on the positional error
In the main article, we assumed uncorrelated morphogen kinetics. Here, we demonstrate how spatial correlation affects the positional error. First, we consider total correlation, where all three kinetic parameters (p, d, D) are the same for all cells, but are still varied between different simulations (different tissues). In this limiting case, morphogen gradient variability occurs only between tissues, not within them. The positional error is significantly greater than with independent cells, and the square-root scaling is lost (Fig. S3, green triangles), because to the morphogen gradient, the tissue effectively appears like a homogeneous continuum with uniform properties. Next, we consider, as a second extreme case, a maximal degree of cell-to-cell correlation in the kinetic parameters, while preserving their probability distributions within the tissue. The kinetic cell parameters (pi, di, Di) are drawn individually and independently for each cell, but are then sorted along the patterning axis and assigned to the cells i, prior to solving the reaction-diffusion equation. Sorting does not affect the patterning precision appreciably, independent of the ordering (Fig. S3). In comparison to zero correlation, sorting slightly reduces the positional error-an effect that is most pronounced for larger cell diameters. But even with this maximal level of spatial cell-to-cell correlation, the square-root scaling of the positional error holds. Intermediate levels of spatial correlation can be expected to yield positional errors lying in between the curves for zero and maximal cell-to-cell correlation.  Fig. 4F). Blue (red) triangles correspond to the case with maximal spatial correlation at given cell-to-cell variability CV p,d,D , where the cell parameters were drawn from log-normal distributions and then sorted in descending (ascending) order. All simulations were repeated n = 10 3 times and the mean positional error ± SEM is plotted.

Choice of the kinetic parameter distribution
In the main article, we assumed log-normally distributed morphogen kinetics. In this section, we show that our results are largely independent of the probability distribution assumed for the kinetic parameters, provided that it meets certain physiological criteria: • The morphogen production rates, degradation rates and diffusivities must be strictly positive. This rules out a normal distribution.
• The probability density of near-zero kinetic parameters must vanish quickly, as otherwise no successful patterning can occur. For example, a tiny diffusion coefficient would not enable morphogen transport over biologically useful distances within useful time periods. This rules out a normal distribution truncated at zero, because very low diffusivities would occur rather frequently for such a distribution.
We repeated the simulations shown in Figs. 2A,B and 4F with a gamma distribution in place of the log-normal distribution. Among other distributions that are conceivable, a gamma distribution with appropriate shape parameter α and inverse scale parameter β fulfills the above criteria. In order to recover the mean and variance of the kinetic parameters, we set α k = 1/CV 2 k and β k = CV 2 k /µ k , where CV k is the coefficient of variation and µ k the mean value of a specific kinetic parameter k. As can be appreciated from Fig. S4, the results are not significantly altered by the specific choice of probability distribution, and our conclusions remain valid. The scaling exponents are consistent within statistical errors.

Effect of cell number in the source domain on gradient precision
In the main article, we showed that patterning precision increases with narrower cells and wider sources. These effects are coupledwider sources will be composed of more cells if the average cell diameter remains constant. In this section, we demonstrate that the positional error is mainly dominated by the cell diameter rather than the source size, and that the found scaling σx ∼ 1/Ls (Eq. 6) is largely due to higher cell numbers in wider sources.
Increasing the number of cells in a source of fixed length improves the precision of the morphogen gradient parameters  S4. Gradient variability and positional error under gamma-distributed morphogen kinetics. All simulations were repeated n = 10 3 times and the mean values ± SEM are plotted. A,B The same scaling laws for the gradient variability found for the gamma and log-normal distributions ( Fig. 2A,B) are consistent. C Different r eadout s trategies (identical to Fig. 4A). D Square-root scaling of the positional error with the cell diameter is found also with gamma-distributed morphogen kinetics. Symbol colours in D correspond to the different m orphogen s ensing s trategies i n C.
according to the asymptotic relationship where N cells is the number of cells in the source domain (Fig. S5A,B). They thus approximately follow the law of large numbers. The positional error decreases analogously with increased cell number in a source of fixed length (Fig. S5C). If, on the other hand, the number of source cells is fixed but the source size increases, the variability in the gradient parameters increases according to power laws (Fig. S5D,E), with exponents α = 0.510 ± 0.005 (Fig. S5D, blue curve) and β = 0.43±0.02 (Fig. S5E, blue curve), suggesting again CV λ,0 ∼ µ δs /Ls. A source composed of a fixed number of cells yields only a mildly greater positional error if its constituent cells have a larger average diameter, however (Fig. S5F). In these simulations, the mean cell diameter in the patterning domain was fixed. Thus, in order to achieve high spatial gradient precision, a morphogen source must have a large number of cells with small diameters, but the cell count is more decisive than the source length.
To study the competition of cell sizes between the source and patterning domain, we then changed the mean cell diameter separately in both subdomains, retaining the mean diameter in the other at a constant value. No further appreciable increase in gradient precision takes place once the mean cell diameter in the source subceeds the one in the patterning domain (µ δs < µ δp , Fig. S6). The mean cell diameter in the source has a limited impact on gradient precision (Fig. S6, pink symbols) compared to the mean diameter in the patterning domain (Fig. S6, yellow symbols). Overall, this suggests that a large number of narrow cells in both the source and patterning domain, but mainly in the latter, is advantageous for patterning precision.

Fit parameters
In Table S1, we list all functional relationships used to fit the data shown in the main article and this supplementary document, together with the fit parameters and their standard errors (SE).