The impact of cell size on morphogen gradient precision

A recently developed numerical framework estimates how much variability in and between morphogen gradients can be accounted for by cell-to-cell variability reported for morphogen production, decay and diffusion (Vetter and Iber, 2022). In this article, we extend the model to take a different perspective on the precision of gradient-based patterning in cellular tissues. We analyse the impact of various length scales present in the epithelium, such as the cell diameter and source size, as well as spatial averaging, on morphogen gradient variability. The findings suggest that positional accuracy is higher, the narrower the cells and the larger the morphogen source.

As a new source of noise, we introduced cell size variability. As the cell area distributions in the Drosophila larval and prepupal wing discs, and in the mouse neural tube resemble log-normal distributions (Sánchez-Gutiérrez et al., 2016; Guerrero et al., 2019), we drew individual cell areas A_i independently of a log-normal distribution with a prescribed mean μ_A and a coefficient of variation CV_A. This allowed us to evaluate the impact of cell-to-cell variability in the production, degradation and diffusion rates p_i, d_i and D_i, as well as in the cell cross-sectional areas A_i, on gradient variability (Fig. 1B).

Cell-to-cell variability in the cross-sectional cell area A does not affect the gradient variability as long as CV_A<1 (Fig. 2C,D). Only for extreme cell area variability exceeding 1 does the variability in λ grow (Fig. 2C). However, we are not aware of any reported CV_A>1 (Guerrero et al., 2019; Kokic et al., 2019 preprint; Gómez et al., 2021; Bocanegra-Moreno et al., 2023). Consequently, cell size has a considerable impact on gradient variability, while physiological levels of variability in the cell area do not contribute to gradient imprecision. A larger source or gradient length reduces only the amplitude variability, but does not affect the decay length variability (Fig. 2E-H). Amplitude and gradient decay length variability is reduced in a source that is composed of many cells with a small mean diameter (see supplementary Materials and Methods for further details, Fig. S5). The parameter values in all reported simulations correspond to those reported for the mouse neural tube (μ_λ=20 μm, μ_δ=5 μm, L_s=5μ_δ and L_p=50μ_δ), unless stated otherwise. At these values, source sizes above 25 μm and gradient decay lengths above 20 μm barely reduce amplitude variability. Sonic hedgehog (SHH) in the neural tube is secreted from both the notochord and the floor plate, while bone morphogenetic protein (BMP) is secreted from both the ectoderm and the roof plate. Intriguingly, while the SHH-secreting notochord shrinks over time, it still measures about 30 μm in width by the 5-somite stage (Imuta et al., 2014), and the SHH-secreting floor plate then emerges in the ventral part of the neural tube and widens over time (Kicheva et al., 2014). The gradient length remains constant at about μ_λ=20 μm (Cohen et al., 2015; Zagorski et al., 2017), the largest value for which the positional error remains small at a large distance (12μ_λ=240 μm) from the source. The source size thus assumes the smallest value and the gradient decay length the largest value for which morphogen gradient variability remains small.

In the case of rectangular rather than circular cell areas, cells are confined to the interval [x₀−r, x₀+r]. The theoretically predicted shift is then approximately 0.052 μm in the mouse neural tube (see supplementary Materials and Methods for further details, Fig. S2) or 1% of the cell diameter. This agrees with the shift we measured in our simulations, Δx=0.0523±0.0001 μm (mean±s.e.m.), confirming that spatial averaging of an exponential gradient results in a higher average concentration than centroid readout. Kinetic and area variability both increase Δx (Fig. 3B), but it remains small enough (small fractions of a cell diameter) to be neglected in the analysis of tissue patterning under biological conditions where r/λ≪1. Linear gradients (Wolpert, 1969) would not result in any shift at all.

Spatial and temporal averaging can reduce the positional error of morphogen gradients (Berg and Purcell, 1977). Previously, these mechanisms have been mainly analysed at the level of the morphogen readouts – typically transcription factors (TFs) – which are averaged by diffusion between nuclei (Houchmandzadeh et al., 2002; Bialek and Setayeshgar, 2005; Gregor et al., 2007; Erdmann et al., 2009; Sokolowski and Tkačik, 2015; Ellison et al., 2016; Mugler et al., 2016). This is easily possible in a syncytium, as present in the early Drosophila embryo, but the role of TF diffusion in increasing patterning precision has remained controversial (Jaeger and Verd, 2020). In an epithelium, nuclei are separated by cell membranes such that the averaging of morphogen-induced factors would require transport between cells, a complex and slow process with many additional sources of molecular noise (Entchev and González-Gaitán, 2002; Lander et al., 2002). However, epithelial cells potentially can reap the benefits of spatial averaging by averaging the morphogen signal over their surface (Fig. 4A, green). Receptors may either be dispersed on the apical cell surface or along the baso-lateral surface, or, in the case of hormones, be limited to nuclei (Saitoh et al., 2013; Zhang et al., 2019). In the last case, morphogen receptors would be limited to a small patch, which could either be randomly positioned (Fig. 4A, blue) or located at the centroid of the cell (Fig. 4A, red). In the mouse neural tube, the SHH receptor PTCH1 is restricted to a cilium located on the apical surface (Saade et al., 2013). The range of spatial averaging then depends on the cilium length and flexibility, rather than the cross-sectional cell area (Fig. 4A, purple). We sought to analyse how the different spatial averaging strategies without crosstalk between neighbouring cells affect the variability of gradients and thus the positional error.

Although the mean cell diameter μ_δ greatly affects the concentration variability CV_C, the readout strategy has only a moderate impact (Fig. 4B). The difference is most pronounced for large cells (μ_δ=μ_λ), where the sensed morphogen variability is largest if the cellular readout point is randomly placed (Fig. 4B, blue). Readout at the centroid or averaged over the entire cell yield similar sensed concentration variabilities. This is understandable because the theoretical considerations above predict only a small shift. In addition, a cilium that averages the gradient concentration over larger regions than a single cell area barely reduces the sensed variability (Fig. 4C).

In summary, larger cross-sectional cell diameters increase the variability of the morphogen concentration profiles, while spatial averaging over the cell surface barely reduces the gradient variability. Spatial averaging may, however, counteract detection noise at low morphogen concentrations far away from the source. It is currently unknown over what distance morphogen gradients operate. At a distance 12λ from the source, for example, exponential concentrations will have declined by e¹²≈160-thousand-fold. At such low levels, detection noise may dominate readout variability unless removed by spatial averaging.

The Decapentaplegic (Dpp) morphogen gradient in the Drosophila wing imaginal disc defines the position of several veins in the adult wing (Fig. 5A). Thus, the anterior-most limits of the Dpp source and the Dpp target gene spalt (sal) define the positions of the third (L3) and second (L2) longitudinal veins in the anterior compartment, respectively (Sturtevant et al., 1997; Bollenbach et al., 2008; Restrepo et al., 2014; Tripathi and Irvine, 2022), while the fifth longitudinal (L5) wing vein forms at the border between the expression domains of optomotor-blind (omb) and brinker (brk) in the posterior compartment (Cook et al., 2004). The Dpp readout positions scale with the total length of the uniformly expanding patterning domain, such that the anterior position of the Sal-domain boundary remains roughly at 40-45% of the anterior domain length L_a, while the posterior Omb domain boundary remains at approximately 50% of the posterior domain length L_p (Bollenbach et al., 2008; Wartlick et al., 2011; Hamaratoglu et al., 2011; Restrepo et al., 2014). The gradient readout positions scale with the length of the patterning domain, because both the gradient length, λ, and the gradient amplitude, C₀, increase dynamically with the expanding tissue (Wartlick et al., 2011; Hamaratoglu et al., 2011; Fried and Iber, 2014, 2015) (Fig. 5B). On their own, the increases in μ_λ and in μ_x would lower the precision of the readout substantially over time (Eq. 6). However, the Dpp source widens in parallel, keeping the μ_λ/L_s ratio at about 0.69 (Fig. 5B). Moreover, the apical cell diameter μ_δ shrinks threefold close to the source from 4.5 to 1.5 μm (Aegerter-Wilmsen et al., 2012; Corrigall et al., 2007; Escudero et al., 2011; Legoff et al., 2013; Kokic et al., 2019 preprint), which somewhat balances the increase in μ_x over time. Plugging these dynamics into our model, the simulations showed that the positional error at μ_x=0.4L_a increases from 2.9 μm to 4.3 μm over developmental time (Fig. 5C, orange diamonds). If no compensation were taking place, the positional error would increase to about 6.5 μm in the same time period (Fig. 5C, blue circles).

The relative patterning precision, as quantified by the coefficient of variation CV_x=σ_x/μ_x, has even been reported to increase during development, as the CV of the distance between the L2 and L3 veins in the adult fly is only half (CV_x=0.08) that of the anterior-most Sal domain boundary (CV_x=0.16) (Bollenbach et al., 2008). How this increase in precision is achieved has remained elusive. In light of Eqn 6, (Fig. 5D), such that the decreasing CV_x in adult stages could at least partly be a consequence of the increase in μ_x=0.4L_a between the stage when the precision of the Sal domain boundary was measured and the termination of Dpp-dependent patterning. The asymptotic relationship may thus provide an explanation of how the relative precision of patterning increases during Drosophila wing disc development.

Our theoretical considerations and simulations above are based on statistical independence between adjacent cells. To examine the effect of spatial correlations, we performed additional simulations in which this assumption was relaxed. We introduced a maximal degree of spatial correlation between neighbouring cells, given a certain degree of intercellular variability CV_k, by sorting the kinetic parameters p_i, d_i and D_i in ascending or descending order along the patterning axis after they had been drawn from their respective probability distributions, and then solved the reaction-diffusion problem (Eqn 2). The square-root increase of the positional error with the mean cell diameter remains intact in the presence of such spatial correlations between cells (see supplementary Materials and Methods for further details, Fig. S3), with a slightly smaller prefactor. Because any physiological level of cell-to-cell correlation that preserves CV_k will lie somewhere between the uncorrelated and the maximally correlated extremes, the impact of such a form of spatial correlation on patterning precision can be expected to be minimal, and our findings also remain valid in presence of spatial correlations.

An additional form of inter-cellular correlation may occur if nearby cells stem from the same lineage and, as such, may have correlated kinetic properties. In its most extreme form, neighbouring cells may share all their molecular parameters, p, d and D, effectively becoming one wider joint cell in our model. We can use our results for cell-autonomous noise to predict the dependency of patterning precision on the number of adjacent cells sharing their kinetic properties, N. As the effective cell diameter simply becomes Nμ_δ, the positional error will scale as ⁠. In this sense, the mean cell diameter μ_δ in our formulas may be interpreted as an effective spatial distance over which morphogen kinetics are shared, proportional to a spatial correlation length in the tissue, if any.

Cell-specific morphogen production and decay rates, and local variability in morphogen transport rates have not yet been quantified in epithelial tissues. A spatial coupling of molecular noise in dividing cells would require a perfectly symmetric division of cell contents upon cell division and the absence of cell-intrinsic noise. Dpp-containing endosomes are indeed distributed equally upon cell division in the Drosophila wing disc (Bökel et al., 2006). However, no cellular system without intrinsic noise has so far been reported. Differences between genetically identical sister cells were first shown for bacterial cells (Elowitz et al., 2002), but have since also been demonstrated for mammalian cells, and pose a key challenge in synthetic biology (Raser and O'Shea, 2005; Zoller et al., 2015; Urban and Johnston, 2018). The coefficients of variation that we used are based on the reported variabilities of production and decay rates in single genetically identical cells in cell culture (Vetter and Iber, 2022).

There are further reasons why low spatial correlation of the kinetic parameters is to be expected. In pseudostratified epithelia, interkinetic nuclear migration (IKNM) introduces differences between cells as the cell cross-sectional areas change along the entire apical-basal axis over time (Gómez et al., 2021). As the tight junctions constitute a diffusion barrier between the apical and the basolateral domains, the apical receptor density between cells will change dynamically between cells if the apical receptor number is equal and fixed for all cells. To maintain the same receptor density, even though IKNM proceeds at different rates between neighbouring cells, as reflected in the different nuclear positions along the apical-basal axis (Gómez et al., 2021), the processes that balance receptor production and internalisation would need to be identical between neighbouring cells, although differences in cell and nuclear volumes may also need to be compensated for. The same holds for the glyocalyx and extracellular matrix, which define the speed of morphogen diffusion, or for fillipodia, in the case of cytoneme-based transport. In summary, the combination of an unequal distribution of cell components in cell division, differences in the relative surface area to cell and nuclear volume, and intrinsic noise in gene expression must be expected to lead to individual differences between neighbouring cells, even if they stem from the same lineage.

After finding that patterning precision is greater with narrower cells in our model, we collected mean apical cell areas for a wide range of tissues from the literature to check whether cell diameters are small in tissues that rely on gradient-based patterning (Fig. 6). In the chick (cNT) and mouse neural tube (mNT), where SHH, BMP and WNT gradients define the progenitor domain boundaries (Briscoe and Small, 2015), the mean apical cell areas are largely around 7 μm² and remain below 12 μm² (Escudero et al., 2011; Guerrero et al., 2019; Bocanegra-Moreno et al., 2023). The chick embryonic ectoderm (cEE) appears to be patterned by BMP gradients (Pera et al., 1999), with mean apical cell area just below 12 μm² (Escudero et al., 2011). In the Drosophila larval eye disc (dEYE), notum (dNP) and wing disc (dWL), Hedgehog (Hh), Decapentaplegic (Dpp) and Wg gradients pattern the epithelium (Tomoyasu et al., 2000; Cavodeassi et al., 2002; Briscoe and Small, 2015), with mean apical cell areas smaller than 7 μm² (Corrigall et al., 2007; Escudero et al., 2011; Aegerter-Wilmsen et al., 2012; Kokic et al., 2019 preprint). The mean apical cell areas of the wing disc increase through the pre-pupal stages (dWP and dPW), to approximately 18 μm² in the pupal stages (Escudero et al., 2011; Kokic et al., 2019 preprint); other measurements in the Drosophila wing disc (dWD) report mean apical cell areas from 0 to 16 μm² (Aegerter-Wilmsen et al., 2012). In the Drosophila eye antennal disc, no gradient-based patterning was described (dEA folded, mean apical cell areas of approximately 33 μm²; dEA non-folded, mean apical cell areas of approximately 39 μm²) (Ku and Sun, 2017). For the peripodal membrane (dPE10-24) of the Drosophila eye disc, no gradient-based patterning has been described and mean apical cell areas range from 85 μm² to more than 300 μm² (Kokic et al., 2019 preprint). In the Drosophila egg chamber (dEC), the mean apical cell areas decline from around 30 μm² at stage 2/3 to around 10 μm² by stage 6/7 (Finegan et al., 2019), consistent with reported gradient-based patterning at stage 6 (Osterfield et al., 2017); we did not find reports of earlier gradient-based patterning. Although gradients pattern the Drosophila blastoderm syncytium (Briscoe and Small, 2015), we are not aware of morphogen gradient readout during cellularisation. In the Drosophila embryo anterior pole (dEAP), the mean apical cell area is approximately 46 μm² and in the embryo trunk (dET) it is roughly 35 μm² (Rupprecht et al., 2017), much larger than in the neural tube or wing disc. Before cellularisation, the situation is different from that in an epithelium in that free diffusion in the inter-nuclear space of the syncytium likely counteracts any sharp transition in the kinetic parameters, as represented in our epithelial model, where cell membranes compartmentalise space. In the Drosophila L2 trachea (dL2 T), no gradients have been reported and the mean apical cell areas are greater than 200 μm² (Skouloudaki et al., 2019). In the mouse embryonic lung (mLUNG), no morphogen gradients have been reported, despite chemical patterning (Iber, 2021), and the mean apical cell area is approximately 19 μm² (Kadzik et al.,2014). Mean apical cell areas in the postnatal (P1-P21) cochlea are between 15 and 55 μm² (Etournay et al., 2010). In adult mouse retinal pigment epithelial (mRPE) cells, the mean apical cell areas exceed 200 μm² in young mice (P30) and increase to over 400 μm² in old mice (P720) (Kim et al., 2021). No gradient-based patterning has been reported in mouse outer hair cells (mOHC1-3; P1, P3, P5, P7.5); mean apical cell areas decrease from 35 μm² (P1) to 16 μm² (P7.5). No gradient-based patterning takes place in the inner hair cells (mIHC1; P1, P3, P6, P7.5); mean apical cell areas decrease from 54 μm² (P1) to 29 μm² (P7.5) (Etournay et al., 2010). No gradient-based patterning has been reported in the mouse ear epidermis (mEE), with mean apical cell areas of 1044 μm² (Yokouchi et al., 2016). The data thus confirm that apical cell areas are small in tissues that employ gradient-based patterning. Our theory makes no prediction about the apical areas in tissues that do not employ gradient-based patterning, but in all cases that we have checked, apical areas are larger and appear to further increase in later developmental stages and in adult animals.

We have shown that gradient precision decreases with increasing cross-sectional area of the patterned cells. Consistent with our prediction, apical surface areas are small in epithelia that employ gradient-based patterning. In curved domains, spatial precision will be higher on the inside, where the average cell diameter is smaller. In the mouse neural tube, the SHH-sensing cilium is indeed located on the inner, apical surface (Saade et al., 2013), while in the flat Drosophila imaginal discs, cells sense Hedgehog along the entire apical-basal axis (Gore et al., 2021). In the Drosophila wing disc, the apical cell diameters shrink in the centre of the domain, such that the apical areas are almost twofold smaller close to the source and increase roughly linearly (Corrigall et al., 2007; Widmann and Dahmann, 2009; Legoff et al., 2013; Bai et al., 2013). In the eye disc, the size gradient is even more pronounced, with tiny apical areas in the Dpp secreting morphogenetic furrow (Corrigall et al., 2007). The declining apical cell diameters have previously been attributed to a mechanical pressure feedback caused by growth (Hufnagel et al., 2007; Aegerter-Wilmsen et al., 2012). However, signalling by Dpp, the fly homolog of mammalian BMP2/4, has been shown to result in taller cells with smaller cross-sectional area in its patterning domain compared with other parts of the Drosophila wing and eye disc (Corrigall et al., 2007; Widmann and Dahmann, 2009; Legoff et al., 2013; Bai et al., 2013). Similarly, the morphogens SHH and WNT have been observed to increase cell height and reduce the cell cross-sectional area via their impact on actin polymerisation, myosin localisation and activity in the embryonic mouse neural tube and lung (Kadzik et al., 2014; Widmann and Dahmann, 2009; Gritli-Linde et al., 2002; Kondo and Hayashi, 2015; Chiang et al., 1996). Complementary to these observations, it would be interesting to test our hypothesis in experiments that alter cell shape using either a genetic or mechanical approach (Neufeld et al., 1998; Duda et al., 2019).

In light of our study, it is possible that the morphogen-dependent reduction in the cross-sectional cell area via positive modulation of cell height serves to enhance patterning precision. The precision advantage of small cell diameters may then have led to the emergence of pseudostratification in epithelial monolayers, a phenomenon that has so far remained unexplained. Our finding that wide cells and very large cell area variability are both detrimental to patterning precision indicate that there is potentially a window for epithelial pseudostratification in which patterning precision is optimal: High cell density benefits precision because cell diameters are small; however, with nuclei much wider than the average cell diameter (Gómez et al., 2021), precision would decline due to large area variability. It is remarkable that all the tissues we analysed seem to lie in the optimal range of this trade-off (Kokic et al., 2019 preprint).

We have revealed scaling relationships between the positional error, cell diameter, gradient decay length and source length (Eqn 6). In follow-up work, we found that these also hold for non-exponential gradients arising from non-linear morphogen degradation (Adelmann et al., 2023), as far as they were studied, and also in 2D tissue patterning (Long et al., 2023 preprint). These relationships predict that morphogen gradients remain highly accurate over very long distances, providing precise positional information even at a distance from the morphogen source. Our results are system-agnostic, and could thus apply widely in development. The compensation between cell diameter, gradient length, source size and readout location, which we have found here, allows a patterning system to tune its length scales to achieve a particular level of spatial precision. Our theoretical work suggests a potential evolutionary benefit for a developmental mechanism that regulates features such as the cell diameter or the ratio to maintain high patterning precision. A loss in precision due to a shift in readout position away from the morphogen source, for example, can be compensated for by narrower cells in the source or in the patterning domain. This allows developmental systems to maintain high patterning precision at readout positions that scale with a growing tissue domain.

Whether pre-steady-state gradients, as likely play a role in the patterning of the Drosophila wing disc (Fried and Iber, 2014), follow the same behaviour as discovered here for the steady state, remains an issue for future research. Assuming that they do, our results offer a potential explanation for the observed increase in relative patterning precision during wing disc development.

The patterning axis was constructed as follows: a random cell area A_i was drawn for cell i=1, and then converted to a diameter ⁠, which assumes that cell surfaces are roughly isotropic. This process was repeated for the next cells i=2, 3,... until their cumulated diameters matched the domain length L_s or L_p. To control the mean cell diameter μ_δ, cell areas were drawn with a mean value of for given μ_δ and CV_A, as follows from the transformation properties of log-normal random variables, such that ⁠. The patterning axis was then discretized into subintervals of length δ_i; the source and patterning domains were pasted together, such that x=0 marked the source boundary; random kinetic parameters p_i, d_i and D_i were drawn independently for each cell from log-normal distributions. The results reported in this work are largely independent of the specific choice of probability distribution, given that they do not allow for very small (or even negative) kinetic parameters, which would not be compatible with a successful morphogen transport and patterning process. A gamma distribution with the same mean and variance, for example, yields largely unchanged behaviour (see supplementary Materials and Methods for further details, Fig. S4).

We then solved Eqn 2 numerically on the discretized domain using Matlab's built-in fourth-order boundary value problem solver bvp4c (version R2020b). Continuity of the morphogen concentration and its flux was imposed at each cell boundary. Further technical details can be found in Vetter and Iber (2022). Each simulation was repeated n=10³ times with independent random parameters and cell areas.

We determined the amplitude C₀ and decay length λ for each numerically generated noisy morphogen gradient by fitting the deterministic solution to it. With no-flux boundaries, the gradient shapes are hyperbolic cosines that slightly deviate from a pure exponential at the far end (Vetter and Iber, 2022). We fitted these inside the patterning domain to obtain C₀ and λ after logarithmisation of the morphogen concentration, as detailed by Vetter and Iber (2022).

We thank Marco Meer for providing cell area data, and Fernando Casares and Nikolaos Doumpas for discussions.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.201702.reviewer-comments.pdf