Cost-precision trade-off relation determines the optimal morphogen gradient for accurate biological pattern formation

Abstract
Introduction
Results
Discussion
Materials and methods
Appendix 1
Appendix 2
Appendix 3
Data availability
References
Article and author information
Metrics

Abstract

Spatial boundaries formed during animal development originate from the pre-patterning of tissues by signaling molecules, called morphogens. The accuracy of boundary location is limited by the fluctuations of morphogen concentration that thresholds the expression level of target gene. Producing more morphogen molecules, which gives rise to smaller relative fluctuations, would better serve to shape more precise target boundaries; however, it incurs more thermodynamic cost. In the classical diffusion-depletion model of morphogen profile formation, the morphogen molecules synthesized from a local source display an exponentially decaying concentration profile with a characteristic length λ. Our theory suggests that in order to attain a precise profile with the minimal cost, λ should be roughly half the distance to the target boundary position from the source. Remarkably, we find that the profiles of morphogens that pattern the Drosophila embryo and wing imaginal disk are formed with nearly optimal λ. Our finding underscores the cost-effectiveness of precise morphogen profile formation in Drosophila development.

Introduction

Complex spatial structures are shaped across the entire body from an ensemble of initially identical cells during animal development. The emergence of spatial structures is generally linked to the pre-patterning of tissues with biomolecules, called morphogens, that instruct cells to acquire distinct cell fates in a concentration-dependent manner (Turing, 1990; Wolpert, 1969; Crick, 1970). At the molecular level, the positional information encoded in the local morphogen concentration is translated by the cells into the expression of specific genes associated with cell-fate determination. In order to generate reproducible spatial organization, the concentration profile of morphogens must be stable against the noisy background signals inherent to the cellular environment.

Along with quantitative measurements of morphogen gradients (Gregor et al., 2007; Kicheva et al., 2007; Bollenbach et al., 2008; Kanodia et al., 2009; He et al., 2010; Zagorski et al., 2017; Petkova et al., 2019), a number of theoretical studies have been devoted to understanding the precision and speed by which morphogen gradients are formed and interpreted (Tostevin et al., 2007; Emberly, 2008; Erdmann et al., 2009; Berezhkovskii et al., 2010; Saunders and Howard, 2009; Morishita and Iwasa, 2009; Kolomeisky, 2011; Tkačik et al., 2015; Lo et al., 2015; Desponds et al., 2020; Fancher and Mugler, 2020); however, generation of morphogen gradient, which breaks the spatial symmetry, incurs thermodynamic cost, and the relation of this cost with the precision and speed of morphogen gradient formation has rarely been addressed except for a few cases (Emberly, 2008). Creating and maintaining morphogen gradients require an influx of energy (Falasco et al., 2018), which is a limited resource for biological systems (Ilker and Hinczewski, 2019), particularly at the stage of embryonic development (Rodenfels et al., 2019; Rodenfels et al., 2020; Song et al., 2019). In the present work, we study the trade-off between the cost and precision of the morphogen profile formation in the framework of the reaction-diffusion model of localized synthesis, diffusion, and depletion (SDD model) (Tostevin et al., 2007; Gregor et al., 2007; Kicheva et al., 2007; Bollenbach et al., 2008; Shvartsman and Baker, 2012; Teimouri and Kolomeisky, 2014).

The precision associated with morphogen gradients are perhaps best exemplified by the patterning of the anterior-posterior axis of the Drosophila embryo through the Bicoid (Bcd) gradient. In the 2-hr post-fertilization of fruit fly embryogenesis, a uniform field of ∼6000 nuclei characterizes the periphery of the shared cytoplasm (Foe and Alberts, 1983). Bcd, a transcription factor, is produced from the anterior-end of the embryo by maternally deposited bcd mRNAs. Its subsequent diffusion and degradation engender an exponentially decaying profile of Bcd concentration (Driever and Nüsslein-Volhard, 1988a; Driever and Nüsslein-Volhard, 1988b), which is translated into the anterior expression of the target gene, hunchback (hb) (Struhl et al., 1989). Nuclei at around the middle section of the embryo can infer their relative spatial positions by detecting the concentration of Bcd with an exquisite precision of ∼ a single nucleus width (Gregor et al., 2007).

Other quantitatively characterized morphogens include Wingless (Wg), Hedgehog (Hh), and Decapentaplegic (Dpp), which pattern the dorsal-ventral (DV) and anterior-posterior (AP) axes of the wing imaginal disk in the fly larvae. Wg, a member of the Wnt signaling pathway, spreads through the wing disk from a narrow band of cells at the DV boundary (Zecca et al., 1996; Neumann and Cohen, 1997). The Wg concentration profile leads to the differential activation of multiple target genes, and, in particular, induces the sharp expression boundary of senseless (sens) a few cells away from the DV boundary (Jafar-Nejad et al., 2006; Bakker et al., 2020). Similarly, Hh patterns the AP axis by spreading from the posterior to the anterior side of the wing disk, affecting the expression of multiple genes. Its downstream transcriptional regulation gives rise to a strip of 8–10 dpp expressing cells that form the AP boundary (Tabata and Kornberg, 1994; Strigini and Cohen, 1997). Dpp, which spreads out from the localized production at the AP boundary, further patterns the AP axis. Major changes in the expression level of Dpp target genes, such as spalt-major (salm), occur at ∼50% of the total length of the domain patterned by Dpp (Nellen et al., 1996; Entchev et al., 2000). Overall, Drosophila wing structures with the spatial precision of a single-cell width emerge from the coordinated actions of multiple patterning events (Abouchar et al., 2014).

The positional information, encoded in the concentration profile $ρ (x)$ (red line in Figure 1), is decoded to yield the target gene expression profile, $g (ρ (x))$ (blue line in Figure 1). The morphogen profile, $ρ (x)$ , and the corresponding gene expression level, $g (ρ (x))$ , together specify the cell fate in a morphogen concentration-dependent manner. In what follows, we will motivate a quantitative expression for the positional error associated with the ‘boundary’ positions, where the target gene expression profile displays a sharp change.

Figure 1

Download asset Open asset

Positional information transfer by the morphogen gradient.

(Top) The specification of the anterior region of the fruit fly embryo. The uniformly distributed nuclei (gray circles) are subjected to different concentrations of the morphogen (red dots) in the local environment, which leads to the anterior expression of the target gene (blue shade). (Bottom) The red and blue lines respectively depict the morphogen profile, $ρ (x)$ , and the target gene expression, $g (x)$ , which together specify cell fate. The squared positional error at the boundary $x_{b}$ , $σ_{x}^{2} (x_{b})$ , is defined as the product between the variance of the morphogen concentration, $σ_{ρ}^{2} (x_{b})$ , and the squared inverse slope of the morphogen profile, ${(\partial_{x} ρ (x))}_{x = x_{b}}^{- 2}$ .

We begin by defining a target boundary, $x = x_{b}$ , where the morphogen concentration is at its critical threshold value $ρ_{b} \equiv ρ (x_{b})$ . The nuclei exposed to morphogen concentrations higher than $ρ_{b}$ would adopt an anterior cell fate. Prior to measuring the local morphogen concentration, each nucleus has no information regarding its own position (Figure 1). After measuring ${\hat{ρ}}_{b}$ , the nucleus can estimate the location of target boundary, such that ${\hat{x}}_{b} = {\hat{x}}_{b} ({\hat{ρ}}_{b})$ . For instance, the Bcd concentration at each nucleus at position $x$ ( $ρ (x)$ ) can be detected in terms of the frequency of Bcd binding to regulatory sequences in DNA. Higher local concentration of Bcd leads to a more frequent expression of the target gene, resulting in the position-dependent expression profile of hb ( $g (ρ (x))$ ) (Gregor et al., 2007). The error (variance) associated with the measured concentration at $x = x_{b}$ with respect to its true value $ρ_{b}$ is given by $σ_{ρ}^{2} (x_{b}) [\equiv ⟨ {({\hat{ρ}}_{b} - ρ_{b})}^{2} ⟩ = (1 / N) \sum_{i = 1}^{N} {({\hat{ρ}}_{b, i} - ρ_{b})}^{2}]$ . In other words, $σ_{ρ}^{2} (x_{b})$ represents the inherent variability of the morphogen profile, which can be determined experimentally by performing multiple measurements, ${\hat{ρ}}_{b, i}$ ( $i = 1, 2, \dots N$ ). For the nucleus to make a single measurement, ${\hat{ρ}}_{b, i}$ , it would integrate the morphogen-induced signal for a short time interval during which the local morphogen concentration remains effectively constant. Then, the Taylor expansion of the concentration at the target boundary, $ρ ({\hat{x}}_{b}) \approx ρ (x_{b}) + {(\partial ρ ({\hat{x}}_{b}) / \partial {\hat{x}}_{b}) |}_{{\hat{x}}_{b} = x_{b}} ({\hat{x}}_{b} - x_{b}) + \dots$ , allows one to relate the error in measured concentration with the positional error in estimating the target boundary $σ_{x}^{2} (x_{b}) [\equiv ⟨ {({\hat{x}}_{b} - x_{b})}^{2} ⟩ = (1 / N) \sum_{i = 1}^{N} {({\hat{x}}_{b, i} - x_{b})}^{2}]$ as follows:

σ_{x}^{2} (x_{b}) ≃ {(\frac{\partial ρ ({\hat{x}}_{b})}{\partial {\hat{x}}_{b}})}_{{\hat{x}}_{b} = x_{b}}^{- 2} σ_{ρ}^{2} (x_{b}),

such that the morphogen profiles exhibiting a large gradient and small fluctuations give rise to small positional errors.

Of fundamental importance is to address how naturally occurring morphogen profiles, tasked with the transfer of positional information, are formed under limited amount of resources. In an earlier work, Emberly considered the total morphogen content at steady state as a proxy for the cost, and evaluated the ‘cost-effectiveness’ of the exponentially decaying Bcd profile with a characteristic decay length λ (Emberly, 2008). Remarkably, it was shown that for a given positional error, the Bcd profile is shaped with a nearly cost-minimizing λ (Emberly, 2008). Here, we expand this argument with an in-depth treatment of the dynamics of morphogen profile formation and the thermodynamic cost involving the formation and maintenance of precise morphogen profiles.

Results

The cost of transferring the positional information should include the cellular resources used to generate the steady state profile prior to the measurement, as well as the resources to maintain the profile during the measurement. Theoretically, cells may control the act of measurement by modulating the availability of morphogen-sensing receptors, or by tuning the overall transcription rate of the target gene. In this context, we consider two limiting scenarios. (i) Point measurement, in which the morphogen profile is measured instantaneously; (ii) Space-time-averaged measurement, in which the morphogen profile is measured over a finite space and for a time interval $T$ . In the former, the associated cost is defined as the amount of morphogen molecules produced while the profile approaches to the steady state. In the latter, we assume a space-time-averaged measurement carried out for a long time duration. Then, the morphogen produced over the measurement can be approximated as the total cost required to create and maintain the morphogen gradient. For both of the limiting scenarios, we will show that the cost of generating morphogen profiles and the precision of the profiles are counterbalanced. We consider a morphogen production-independent quantity by taking the product of total cost and the positional error to quantify the trade-off between the cost and precision of the morphogen profile, and show that the trade-off product can be minimized when the morphogen gradient’s characteristic length, λ, is properly selected. We evaluate the cost-effectiveness of morphogen profiles patterning the fruit fly embryo and wing disk, which have been quantitatively characterized in a number of studies (Gregor et al., 2007; Kicheva et al., 2007; Bollenbach et al., 2008; Nahmad and Stathopoulos, 2009; Drocco et al., 2011; Wartlick et al., 2011; Perry et al., 2012; Zhou et al., 2012; Chaudhary et al., 2019; Petkova et al., 2019; Bakker et al., 2020).

Point measurement

We begin by defining the dynamics of morphogen profile formation in a system composed of a one-dimensional (1D) array of cells in the domain $0 \leq x \leq L$ , where $L$ is the system size (see Figures 1 and 2(a)). The term ‘cell’ refers to the spatial grid in which to define the local positional error and the cost. We denote the amount of morphogen in the cell of the volume ( $v_{c e l l}$ ) and length ( $l_{c e l l}$ ) at the interval between $x$ and $x + l_{cell}$ by $ρ (x, t) [conc \equiv # / v_{cell}]$ . Unless otherwise specified, ρ refers to the ensemble averaged morphogen concentration. At the left boundary (see Figure 2a), the morphogen is injected at a constant flux of j_in $[conc \times l_{cell} / time]$ . At all positions, the morphogen molecules are depleted with rate $k_{d}$ [ ${time}^{- 1}$ ], while spreading across the cells with the diffusivity $D [l_{cell}^{2} / time]$ (Figure 2a). For $L / l_{cell} ≫ 1$ (i.e. at the continuum limit), the spatiotemporal dynamics of the morphogen is described using the reaction-diffusion equation

\partial_{t} ρ (x, t) = D \partial_{x}^{2} ρ (x, t) - k_{d} ρ (x, t),

with boundary conditions $- {D \partial_{x} ρ (x, t) |}_{x = 0} = j_{in}$ and ${D \partial_{x} ρ (x, t) |}_{x = L} = 0$ . Then, the concentration profile at steady state is obtained as

ρ_{s s} (x) = \frac{j_{i n}}{\sqrt{D k_{d}}} \frac{c o s h (\frac{L - x}{λ})}{s i n h (\frac{L}{λ})} \approx \frac{j_{i n}}{\sqrt{D k_{d}}} e^{- x / λ},

where $λ = (\sqrt{D / k_{d}})$ is the characteristic decay length determined by the diffusivity of morphogen ( $D$ ) and the depletion rate ( $k_{d}$ ). For large system size ( $L ≫ λ$ ), $ρ_{ss} (x)$ is simply an exponentially decaying profile with characteristic length λ. In what follows, we will define (i) the cost, and (ii) the precision associated with this reaction-diffusion model.

Figure 2 with 1 supplement see all

Download asset Open asset

Cost-precision trade-off associated with the point measurement.

(a) Schematic of the model. (b) The position-dependent lower bound of the trade-off product $π_{o, \min} (x_{b})$ , obtained numerically. The gray hashed area represent the inaccessible regions. The trade-off product of the morphogen profiles of Bcd, Wg, Hh, and Dpp are shown in the respective insets. (c) The black line denotes the optimal characteristic decay length ( $λ_{\min}$ ) with respect to the position $x_{b} / L$ . The color scale indicates the trade-off product $π_{o} (x_{b})$ computed for each pair of $λ / L$ and $x_{b} / L$ values. The grea dotted line depicts the linear approximation of $λ_{\min}$ at large $L$ . In (b) and (c), the depicted trade-off product is normalized by $α_{o} L / l_{cell}$ . The parameters for the naturally occurring morphogen profiles are further described in Appendix 3, Length scales of Bcd, Wg, Hh, and Dpp, and Appendix 3—table 1.

(i) With a function $R (x, t) \equiv (ρ_{ss} (x) - ρ (x, t)) / (ρ_{ss} (x) - ρ (x, 0))$ , which captures the evolution of the morphogen profile from $ρ (x, t = 0) = 0$ to the steady state value $ρ_{ss} (x)$ (Equation 3), the mean local accumulation time at $x$ that characterizes the average time scale for establishing the steady state profile is calculated as Berezhkovskii et al., 2010,

τ (x) = \int_{0}^{\infty} R (x, t) 𝑑 t = k_{d}^{- 1} f (x; λ),

where $f (x; λ) \equiv \frac{1}{2} [1 + \frac{L}{λ} \coth (\frac{L}{λ}) - \frac{L - x}{λ} \tanh (\frac{L - x}{λ})]$ is a dimensionless quantity determined by $x$ and λ. Then, $v_{cell} j_{in} / l_{cell} \times τ (x)$ quantifies the total number of morphogens produced over the time scale of $τ (x)$ . The thermodynamic cost of producing the morphogens, which is effectively the driving force of pattern formation, is expected to be proportional to this number (Lynch and Marinov, 2015), such that

𝒞 (x) = α_{o} (\frac{v_{cell} j_{in}}{l_{cell} k_{d}}) f (x; λ) \approx α_{o} (\frac{v_{cell} j_{in}}{l_{cell} k_{d}}) \frac{1}{2} (1 + \frac{x}{λ}),

where $C (x)$ increases linearly with $x$ for large system size ( $L ≫ λ$ ) (Berezhkovskii et al., 2010).

The proportionality constant, $α_{o}$ , amounts to the thermodynamic cost of synthesizing and degrading a single morphogen molecule, which can be quantified by the number of ATPs hydrolyzed in the process. For instance, the thermodynamic cost required to translate and degrade a single Bcd protein composed of 494 amino acids is on the order of $α_{o} \approx 494 \times 4 \approx 2 \times 10^{3}$ ATPs, where we assume that each peptide-bond formation requires 4 ATPs (Lynch and Marinov, 2015). In the first ∼2 hr of Drosophila embryogenesis when the anterior-posterior axis patterning takes place, ∼5 × 10⁸ Bcd molecules are produced (Drocco et al., 2012). Thus, the cost of generating the Bcd profile is on the order of 10¹² ATPs, which is a small fraction of the total energy budget of Drosophila embryogenesis (∼7 × 10¹⁶ ATPs) (Song et al., 2019). However, the energy budget of developing systems remains largely uncharacterized, and there is certainly a lack of understanding on how much of the cost can be attributed to 'housekeeping processes', as opposed to the cost of generating new spatial structures (Rodenfels et al., 2019; Song et al., 2019). The generation of each morphogen profile is an indispensable developmental process, although its thermodynamic cost is relatively small compared to the total energy budget.

(ii) The local measure of precision at position $x$ can be quantified by the squared relative error in the positional measurement of morphogen profile (Equation 1):

\begin{array}{ll} ϵ^{2} (x) & \equiv \frac{σ_{x}^{2} (x)}{x^{2}} = {(\frac{\partial \log ρ_{s s} (x)}{\partial \log x})}^{- 2} \frac{σ_{ρ}^{2} (x)}{⟨ \hat{ρ} (x) ⟩^{2}} \\ = (\frac{λ^{3}}{x^{2} l_{c e l l}}) {(\frac{v_{c e l l} j_{i n}}{l_{c e l l} k_{d}})}^{- 1} \frac{s i n h (\frac{L}{λ}) c o s h (\frac{L - x}{λ})}{{[s i n h (\frac{L - x}{λ})]}^{2}} \\ \approx (\frac{λ^{3}}{x^{2} l_{c e l l}}) {(\frac{v_{c e l l} j_{i n}}{l_{c e l l} k_{d}})}^{- 1} e^{x / λ}, \end{array}

where $\hat{ρ} (x)$ represents the morphogen concentration estimated at $x$ from a measurement. For large system size ( $L ≫ λ$ ), $ϵ^{2} (x) \propto (λ^{3} / x^{2}) e^{x / λ}$ . Note that the multiple measurements of $\hat{ρ} (x)$ yield the mean concentration profile $⟨ \hat{ρ} (x) ⟩ = (1 / N) \sum_{i = 1}^{N} {\hat{ρ}}_{i} (x)$ , which is equivalent to $ρ_{ss} (x)$ at steady states given in Equation 3. We use the fact that the probability distribution of the steady state concentration profile obeys the Poisson statistics (Heuett and Qian, 2006) (i.e. $v_{cell}^{2} σ_{ρ}^{2} (x) = v_{cell} ⟨ \hat{ρ} (x) ⟩ = v_{cell} ρ_{ss} (x)$ ), since all the reactions in the system are of zeroth or first order. Although the present study mainly focuses on the mathematical form of the precision, $ϵ^{2} (x)$ is also bounded by physical constraints such as the size of the morphogen sensing machineries that bind and measure the local concentration (Tostevin et al., 2007). Experimentally determined values of the precision range from $ϵ^{2} (x_{b}) \approx 7 \times 10^{- 4}$ for Bcd to $ϵ^{2} (x_{b}) \approx 1.5 \times 10^{- 2}$ for Dpp, where $x_{b} \approx (0.4 - 0.5) L$ for both systems (Gregor et al., 2007; Bollenbach et al., 2008).

There is a trade-off between $𝒞 (x)$ and $ϵ^{2} (x)$ . If $n_{m}$ is the number of morphogen molecules produced for a certain time duration then $𝒞 (x) \propto n_{m}$ and $σ_{ρ}^{2} (x) / {⟨ \hat{ρ} (x) ⟩}^{2} \propto 1 / n_{m}$ , the latter of which simply arises from the central limit theorem, or can be rationalized based on the Berg-Purcell result ( $δ \hat{ρ} / \hat{ρ} \sim 1 / \sqrt{D a \hat{ρ} τ}$ ) (Berg and Purcell, 1977) where $a$ is the radius of the volume in which a receptor detects the morphogen and τ is the detection time. Increasing the overall morphogen content reduces the morphogen profile’s positional error at the expense of a larger thermodynamic cost. In Equations 5 and 6, $n_{m}$ corresponds to $(v_{cell} j_{in} / l_{cell} k_{d})$ (the morphogen molecules produced for the time duration $k_{d}^{- 1}$ ). The n_m-independent trade-off between the cost of generating the morphogen profile and the squared relative error in the position of morphogen profile can be quantified by taking the product of the two quantities,

π_{o} (x; λ) \equiv 𝒞 (x; λ) \times ϵ^{2} (x; λ) \geq π_{o, \min} (x; λ),

where the λ dependence of the trade-off product is made explicit for the discussion that follows.

We are mainly concerned with the precision at the boundary position, $x = x_{b}$ , where a sharp change in the downstream gene expression, $g (x)$ , is observed. The inequality in Equation 7 constrains the properties of the morphogen profile at $x_{b}$ , specifying either the minimal cost of generating the morphogen profile for a given positional error or the minimal error in the morphogen profile for a given thermodynamic cost (Figure 2). In other words, when the trade-off product $π_{o} (x_{b}; λ)$ is close to its lower bound $π_{o, \min} (x_{b}; λ)$ , the system is cost-effective at generating precise morphogen profiles at $x_{b}$ . Generally, $π_{o, \min} (x_{b}; λ)$ increases monotonically with $x_{b}$ , which signifies that boundaries farther away from the origin require the synthesis of more morphogen molecules to achieve comparable positional error (Figure 2b). For a fixed system size $L$ and position $x_{b}$ , the value of the trade-off product is determined solely by the decay length, λ. By tuning λ, one can find the value of $λ_{\min}$ that minimizes the trade-off product at $x = x_{b}$ , that is, $π_{o, \min} (x_{b})$ , such that

λ_{\min} (x_{b}) = \arg \min_{λ} π_{o} (x_{b}; λ) .

The optimal decay length, $λ_{\min} (x_{b})$ , also increases monotonically with the target boundary $x_{b}$ ; for $x_{b} < L / 2$ , $λ_{\min}$ is well approximated by the linear relationship $λ_{\min} (x_{b}) \approx 0.43 x_{b}$ (See Equation 40 for the derivation).

It is of great interest to compare $x_{b}$ and λ of real morphogen profiles against the optimal $λ_{\min} (x_{b})$ that one can predict from Equation 8. For the morphogen profiles of Bcd, Wg, Hh, and Dpp, we first estimated the possible range of $x_{b}$ from the experimentally measured expression profiles of the respective target genes, hb, sens, dpp and salm (Torroja et al., 2004; Perry et al., 2012; Bakker et al., 2020). The range of the inferred boundary positions where the target gene expressions undergo relatively sharp changes are shown on the right column of Figure 2—figure supplement 1 and listed in Appendix 3—table 1. The morphogen concentration profiles obtained from experiments are shown on the left column of Figure 2—figure supplement 1, where we overlaid the exponential profiles with characteristic lengths $λ_{B c d} = 100 μ m$ , $λ_{W g} = 6 μ m$ , $λ_{H h} = 8 μ m$ , and $λ_{D p p} = 20 μ m$ (Houchmandzadeh et al., 2002; Kicheva et al., 2007; Wartlick et al., 2011). Figure 2b shows the pairs of $x_{b}$ and λ normalized by the system size, $L$ , which reveals that the λ’s of naturally occurring morphogen profiles are close to their respective optimal values, $λ_{\min} (x_{b})$ . Taken together, our findings indicate that the concentration profiles of the four morphogens are effectively formed under the condition that the cost-precision trade-off, $π_{o} (x_{b}; λ)$ , is minimized to its lower bound, $π_{o, \min} (x_{b}; λ)$ (Figure 2b) (see the Appendix 3, Length scales of Bcd, Wg, Hh, and Dpp, for more details on the relevant parameters of the morphogen profiles).

Our theoretical result must be interpreted with care. For the Bcd profile with the target boundary $x_{b} \approx 0.4 L$ , $π_{o} (x_{b}; λ_{Bcd}) \approx α_{o} (L / l_{cell}) 0.6$ (Figure 2b) where $L / l_{cell} \approx 50$ , and the experimentally reported precision is $ϵ^{2} (x_{b}) \approx 7 \times 10^{- 4}$ (Gregor et al., 2007). Thus, the cost associated with the Bcd profile is $C (x_{b}; λ_{Bcd}) = π_{o} (x_{b}; λ_{Bcd}) / ϵ^{2} (x_{b}; λ_{Bcd}) \approx 4 \times 10^{4} α_{o}$ , or equivalently, the thermodynamic cost of synthesizing and degrading 4 × 10⁴ Bcd molecules. This is orders of magnitude smaller than the earlier estimate of ∼5 × 10⁸ based on the experimental measurement of total Bcd synthesis in the embryo (Drocco et al., 2012). The discrepancy between the two numbers most likely arises because our theory simplifies the dynamics of the nuclear concentration of Bcd in one dimension, whereas the direct measurement of Bcd synthesis accounts for all the molecules in 3D volume including the cytoplasm and the nuclei (Gregor et al., 2007; Drocco et al., 2012). Additionally, the duration to generate the steady state morphogen profile is longer than the characteristic time, τ. However, assuming that the cytoplasmic and nuclear concentrations are proportional to each other, and that the time for the morphogen profile to reach steady state is proportional to τ, it is possible to adjust the proportionality constant $α_{o}$ so that $C (x_{b})$ represents the overall cost to produce the concentration profile of nuclear Bcd.

Space-time-averaged measurement

The positional information available from the morphogen profile varies depending on the method by which the cell senses and interprets the local morphogen concentration. Our previous definition of $ϵ^{2} (x)$ (Equation 6) represents the positional error from a single independent measurement of the morphogen concentration at a specific location in space. However, organisms may further improve the positional information of the morphogen gradient through space-time-averaging.

To incorporate the scenario of space-time-averaging, we assume that the cell at position $x$ uses a sensor with size $a$ to detect the morphogen concentration over the space interval $(x - a, x + a)$ for time $T$ (Figure 3a). Then the space-time-averaged molecular count is written as

m_{i} (x) = \frac{1}{T} \int_{0}^{T} 𝑑 t \int_{x - a}^{x + a} {\hat{ρ}}_{i} (y, t) 𝑑 y,

where ${\hat{ρ}}_{i} (y, t)$ represents an estimate of molecular count at position $y$ from an $i$ -th independent measurement. Repeated measurements ( $N ≫ 1$ ) lead to the mean value of molecular counts

⟨ m (x) ⟩ = \frac{1}{N} \sum_{i = 1}^{N} m_{i} (x) = \frac{2 j_{in}}{k_{d}} e^{- x / λ} \sinh (a / λ) .

In this case, $j_{i n}$ and $k_{d}$ are in the units of $# / time$ and ${time}^{- 1}$ , respectively. With a sufficiently long measurement time ( $T ≫ k_{d}^{- 1}$ ), the variance of $m (x)$ can be approximated to Fancher and Mugler, 2020

σ_{m}^{2} (x) \equiv \frac{1}{N} \sum_{i = 1}^{N} {(m_{i} (x) - ⟨ m (x) ⟩)}^{2} = \frac{4 j_{in} e^{- x / λ} \sinh (a / λ)}{T k_{d}^{2}} (1 - \frac{(λ / a) \cosh (a / λ) \sinh (a / λ) + 1}{2 e^{a / λ} (λ / a) \sinh (a / λ)}) .

Analogously to Equation 6, the squared relative error (precision) at position $x$ can be quantified as follows,

ϵ_{T}^{2} (x) \equiv {(\frac{\partial \log ⟨ m (x) ⟩}{\partial \log x})}^{- 2} \frac{σ_{m}^{2} (x)}{{⟨ m (x) ⟩}^{2}} = \frac{{(λ / x)}^{2} e^{x / λ}}{j_{in} T} (\frac{e^{- a / λ} + 3 e^{3 a / λ} - 4 (1 + a / λ) e^{a / λ}}{2 {(e^{2 a / λ} - 1)}^{2}}) .

The cost of maintaining the morphogen profile at steady states is proportional to $j_{in} T$ , which simply represents the total number of morphogens produced for $T$ , such that the total thermodynamic cost for producing morphogens for time $T$ is

𝒞_{T} = α_{o} j_{in} T .

In the product between the net amount of morphogen molecules synthesized and the positional error from the averaged signal, the number of morphogens synthesized for time $T$ , $j_{in} T$ cancels off, which yields the following expression of the cost-precision trade-off:

π_{T} (x; λ) \equiv 𝒞_{T} \times ϵ_{T}^{2} (x; λ) \geq π_{T, \min} (x; λ) .

With a given target boundary $x = x_{b}$ and the sensor size $a$ , the value of $π_{T} (x_{b}; λ)$ is solely determined by λ, analogously to $π_{o} (x_{b}; λ)$ . Thus, the lower bound of $π_{T} (x_{b}; λ)$ , i.e., $π_{T, \min} (x_{b}; λ)$ , can be determined simply by tuning λ to a value that minimizes $π_{T} (x_{b}; λ)$ , that is, $λ_{\min} (x_{b})$ . Both $π_{T, \min} (x_{b}; λ)$ and $λ_{\min} (x_{b})$ increase monotonically with $x_{b}$ . In particular, $λ_{\min} (x_{b})$ is found in the range of $(x_{b} - a) / 2 \leq λ_{\min} (x_{b}) \leq x_{b} / 2$ (see Equations 42 and 43 for derivation). The boundary of the inaccessible region in Figure 3b constrains $π_{T} (x_{b}; λ)$ , suggesting that there is a minimal morphogen production for a given positional error. For instance, in order to suppress the positional error down to 10% ( $ϵ_{T} (x_{b}) \approx 0.1$ ) at $x_{b} = 60 a$ , the characteristic length must be $λ \approx 30 a$ , which demands that the system synthesize at least ∼182 morphogen molecules (Figure 3b).

Figure 3

Download asset Open asset

Cost-precision trade-off associated with the space-time-averaged measurement.

(a) Schematic of the model. (b) The optimal trade-off product $π_{T} (x_{b})$ obtained numerically with respect to the location of the target boundary position normalized by the sensor size ( $x_{b} / a$ ). The gray hashed area represent the inaccessible regions. Shown are the trade-off products of the morphogens profiles of Bcd, Wg, Hh, and Dpp. (c) The black line denotes the optimal characteristic decay length ( $λ_{\min}$ ) with respect to the position $x_{b} / a$ . The color scale indicates the trade-off product $π_{T} (x_{b})$ computed for each pair of $λ / a$ and $x_{b} / a$ values. The gray dotted line depicts the linear approximation of $λ_{\min}$ . The white dotted line depicts the linear approximation of $λ_{\min}$ for the point measurement model. In (b) and (c), the depicted trade-off product is normalized by $α_{o}$ . The parameters for the naturally occurring morphogen profiles are further described in the Appendix 3, Length scales of Bcd, Wg, Hh, and Dpp, and Appendix 3—table 1.

As we have done for $π_{o}$ in the point measurement, we can evaluate $π_{T}$ of the four morphogen profiles at their target boundaries ( $x_{b}$ ), by using the cell size as a proxy for the sensor size ( $a$ ), and by assuming that the morphogen profile is measured for a sufficiently long time (i.e. $T ≫ k_{d}^{- 1}$ , the validity of which is further discussed in the Appendix 3, Time scales of Bcd, Wg, Hh, and Dpp). The boundary position and characteristic length associated with each morphogen, normalized by the respective cell size, are shown in Figure 3c (see Appendix 3, Length scales of Bcd, Wg, Hh, and Dpp and Appendix 3—table 1 for more detail). We find that the characteristic decay lengths ( $λ = \sqrt{D / k_{d}}$ ) of all four morphogen profiles are close to their respective $λ_{\min} (x_{b})$ values. Thus, from the space-time-averaged measurement of morphogen profiles, the four morphogen profiles are also formed under the conditions of nearly optimal cost-precision trade-off (Figure 3b).

Discussion

Comparison of the trade-off products, $π_{o}$ and $π_{T}$

Although the two models of biological pattern formation seem to differ in the definitions of both the precision and the associated cost, the measure of precision in the two models are in fact limiting expressions of each other. The concentration detected in the measurement with space-time-averaging, defined as $m (x) / (2 a)$ at the limit of small sensor size ( $a / λ ≪ 1$ ) (Equation 10) is identical to the one in the point measurement in the limit of $L ≫ λ$ (Equation 3), both yielding $(j_{in} / \sqrt{D k_{d}}) e^{- x / λ}$ . On the one hand, $ϵ^{2} (x)$ , defined through $ρ_{ss} (x)$ , is experimentally quantifiable through repeated measurements of the morphogen profile from the images of fixed samples (Kicheva et al., 2007; Gregor et al., 2007; Bollenbach et al., 2008). On the other hand, $ϵ_{T}^{2} (x)$ , defined through $m (x)$ , may represent the precision accessible to cells that integrate the signal over time from multiple receptors positioned across the cell surface. Effectively, $ϵ^{2} (x)$ can be interpreted as the short time limit of $ϵ_{T}^{2} (x)$ measured by a single receptor.

The cost associated with the morphogen profile in the case of point measurement, $v_{cell} j_{in} l_{cell}^{- 1} τ (x)$ , is a quantity that reflects the number of morphogen molecules required to reach steady state at position $x$ . In contrast, the cost of morphogen production with space-time-averaging, $j_{in} T$ , integrated over a time interval $T ≫ k_{d}^{- 1}$ , is identical for all positions. The latter assumption is required in the derivation of the expression of $ϵ_{T}^{2} (x)$ (see Appendix 1 of Fancher and Mugler, 2020). We additionally demand $T ≫ τ (x)$ in order for $j_{in} T$ to represent the total cost of morphogen profile formation and maintenance. However, $T$ may not necessarily be larger than either $k_{d}^{- 1}$ or $τ (x)$ . For instance, the Bcd profile degrades and stabilizes at time scales of $\sim 1$ hour (Berezhkovskii et al., 2010; Drocco et al., 2012), but must be measured within 10 min by the nuclei (Lucas et al., 2018) (see further discussion in the Appendix 3, Time scales of Bcd, Wg, Hh, and Dpp). The energetic cost of naturally occurring morphogen profiles is likely determined in between the two limiting cases.

For actual biological systems, the cost-precision trade-off involving the pattern formation is presumably at work in between the two scenarios. Thus, of great significance is the finding that $π_{o} (x; λ)$ and $π_{T} (x; λ)$ for the two limiting models can be effectively minimized by similar λ values at a given position: $λ_{\min} \approx 0.43 x_{b}$ for the point measurement, and $λ_{\min} \approx 0.5 x_{b}$ for the space-time-averaging.

The entropic cost of forming precise morphogen profiles

Biological pattern formation by morphogen gradients is a process operating out of equilibrium (Falasco et al., 2018), in which thermodynamic cost is incurred to generate a morphogen gradient with minimal error for the transfer of positional information against stochastic fluctuations. In the Appendix 1, Reversible reaction-diffusion model of morphogen dynamics, we extend our discussion of cost-precision trade-off on the models that include uni-directional irreversible steps in the synthesis and depletion by considering a more general reaction model with bi-directional reversible kinetics where the forward and reverse rates are well defined at every elementary process (Appendix 1—figure 1. A and Equation 16). Such a model allows us to quantify the entropy production rate ( ${\dot{S}}_{tot}$ ) or the thermodynamic cost for the formation of morphogen gradient and clarifies the physical meaning of $α_{o}$ (for the relation between $C$ (Equation 5) and ${\dot{S}}_{tot}$ , see Appendix 1 and Equations 32, 33). Similar to the point measurement, we derive expressions for the relaxation time to the steady state ( $τ_{rev} (x)$ ) and the positional error ( $ϵ_{rev}^{2} (x)$ ). The trade-off among the thermodynamic cost, speed of formation ( $\sim τ_{rev}^{- 1} (x)$ ), and precision is quantified by the product of the three quantities,

π_{o, rev} (x) \equiv {\dot{S}}_{tot} τ_{rev} (x) ϵ_{rev}^{2} (x) \geq π_{o, rev, \min} (x) .

The lower bound, $π_{o, rev, \min} (x_{b})$ , increases monotonically with $x_{b}$ (Appendix 1—figure 1).

With ${\dot{S}}_{tot} τ_{rev} (x_{b})$ representing the entropy production during the time scale associated with the morphogen profile formation, Equation 15 is reminiscent of the thermodynamic uncertainty relation (TUR), a fundamental trade-off relation between the entropy production and the precision of a current-like output observable for dynamical processes generated in nonequilibrium with a universal bound (Barato and Seifert, 2015; Gingrich et al., 2016; Hyeon and Hwang, 2017; Hwang and Hyeon, 2018; Horowitz and Gingrich, 2020; Song and Hyeon, 2020). However, unlike TUR which has a model independent lower bound, the lower bound of Equation 15 is model-specific and position-dependent. The cost-precision trade-off of pattern formation discussed in this study fundamentally differs from that of TUR (Pietzonka et al., 2017; Horowitz and Gingrich, 2017; Dechant and Sasa, 2020; Song and Hyeon, 2021), in that the concentration profile of morphogen is not a current-like observable with an odd-parity with time-reversal; however, still of great significance is the discovery of the underlying principle that the pattern formation is quantitatively bounded by the dissipation. For future work, it would be of great interest to consider employing the morphogen-induced currents through signaling pathways, such as transcription, as the output observable to which TUR can be directly applied.

Cost-effectiveness of precise morphogen profile formation

There are multiple possibilities of achieving the same critical threshold concentration $ρ (x_{b})$ at the target boundary $x = x_{b}$ through the combination of morphogen synthesis rate ( $j_{i n}$ ), diffusivity ( $D$ ) and depletion rate ( $k_{d}$ ) (inset in Figure 4). A steeper morphogen profile with larger $j_{i n}$ and smaller $λ (= \sqrt{D / k_{d}})$ (red in Figure 4a) leads to a more precise boundary but it incurs a higher thermodynamic cost. The opposite case with a morphogen profile with smaller $j_{i n}$ and larger λ gives rise to a less precise boundary with a lower cost (brown in Figure 4a). As the main result, we show that λ can be tuned to minimize the cost-precision trade-off product, which leads to the formation of cost-effective morphogen profiles. In other words, morphogen profiles with optimal characteristic length $λ_{\min}$ achieve a desired precision when the cost is minimal (Figure 4b). Specifically, $λ_{\min} \approx (0.43 - 0.50) x_{b}$ , which suggests that the target boundary of the exponentially decaying morphogen profile, $c (x) = c_{0} e^{- x / λ}$ , should be formed at $c (x_{b}) \approx 0.1 c_{0}$ . Remarkably, the λ’s of naturally occurring morphogen profiles in fruit fly development are close to their respective optimal values, and we confirm that the target boundaries of biological pattern from those profiles are identified at $c (x_{b}) \approx 0.1 c_{0}$ (see Figure 2—figure supplement 1). These findings lend support to the hypothesis that along with the reduction in the positional error, the thermodynamic costisone of the key physical constraints in designing the morphogen profiles.

Figure 4

Download asset Open asset

Optimal concentration profile of morphogens.

(a) Three possible morphogen profiles (i) (red), (ii) (cyan), (iii) (brown) with different λ’s ( $λ^{(i)} < λ^{(i i)} < λ^{(i i i)}$ ), generated with different values of morphogen influxes $j_{i n}$ ( $j_{in}^{(i)} > j_{in}^{(ii)} > j_{in}^{(iii)}$ ). The morphogen concentration of the three profiles coincide at $x_{b}$ , giving rise to the same threshold value $ρ (x_{b})$ but different positional errors ( $ϵ^{2} (x_{b}; λ)$ ). (b) The precision of three possible morphogen profiles (i) (red), (ii) (cyan), (iii) (brown) with different λ’s ( $λ^{(i)} < λ^{(i i)} < λ^{(i i i)}$ ), generated with different values of morphogen influxes $j_{i n}$ ( $j_{in}^{(i)} > j_{in}^{(iii)} > j_{in}^{(ii)}$ ). The λ values are identical to those with matching colors in (a), but the red and brown curves are generated with different $j_{i n}$ values from those in (a). The cost associated with each morphogen profile are shown in units of $α_{o} L / l_{cell}$ . (c) The diagram of the trade-off product associated with the point measurement, $π_{o} (x; λ)$ , plotted with respect to $x$ and λ. The black line indicates the optimal decay length, $λ_{\min}$ at position $x$ . Shown on the diagram are the trade-off product $π_{o}$ ’s for the three cases shown in (a) and (b). (d) The value of $π_{o}$ as a function of λ at $x = x_{b}$ . The trade-off product is minimized to $π_{o} = π_{o}^{(ii)} \approx α_{o} (L / l_{cell}) 0.4$ with $π_{o}^{(ii)} < π_{o}^{(iii)} < π_{o}^{(i)}$ .

Concluding remarks

The classical SDD model offers a simple and powerful framework to study basic properties of morphogen dynamics (Gregor et al., 2007; Kicheva et al., 2007; Bollenbach et al., 2008; Emberly, 2008; Berezhkovskii et al., 2010; Teimouri and Kolomeisky, 2014; Fancher and Mugler, 2020). The molecular mechanisms underlying the formation of morphogen profiles are, however, much more complex than those discussed here. Even the seemingly simple diffusive spreading of the morphogen can originate from many different mechanisms (Müller et al., 2013). In fact, among the four biological examples shown in Figures 2 and 3, it has been suggested that Wg and Hh spread over the space through active transport mechanisms rather than through passive diffusion (Huang and Kornberg, 2015; Teimouri and Kolomeisky, 2016; Chen et al., 2017; Fancher and Mugler, 2020; Rosenbauer et al., 2020). Furthermore, the morphogens considered in the present study induce the expression of multiple target gene expressions, potentially leading to multiple spatial boundaries (Torroja et al., 2004; Hannon et al., 2017; Bakker et al., 2020). While our theory suggests that the cost-precision trade-off can only be optimized at a single target boundary, one can consider an extension of our study in which a weighted average of the precision of many boundary positions is balanced against the total cost of generating the morphogen profile.

Modified trade-off relations in systems with more complex geometries, such as the 1D model with a distributed morphogen source, and the morphogen dynamics on a sphere, can be conceived as well (see Appendix 2). For the latter case, the λ values, for the morphogens inducing the endoderm and mesoderm of zebrafish embryos, are found far greater than those leading to the trade-off bound (Appendix 2—figure 2b). The large trade-off product values may either simply indicate that the thermodynamic cost is not necessarily a key physical constraint or reflect the presence of other molecular players in a more complex mechanism establishing the zebrafish germ layers. (Nowak et al., 2011; Müller et al., 2012; Dubrulle et al., 2015; Almuedo-Castillo et al., 2018).

Generally, multiple morphogen profiles relay combinatorial input signals that determine the expressions of an array of downstream target genes and their subsequent interactions (Briscoe and Small, 2015; Tkačik and Gregor, 2021). For instance, the anterior patterning of the Drosophila embryo depends on the dynamic interpretation of maternal input signals from bcd, nanos, and torso (Liu et al., 2013). Nevertheless, each morphogen profile is a fundamental component of biological pattern formation, and its thermodynamic cost must be taken into account in the energy budget of a developing organism. In this light, our theory proposes a quantitative framework to evaluate the cost-precision trade-off of individual morphogen profiles, which allows us to show that the morphogen profiles of fruit fly development are nearly optimal, as opposed to those of the zebrafish embryo. For future work, it would be of great interest to extend our approach to more complex tissue patterning mechanisms, for which key insights are being generated from ongoing studies on the fruit fly, zebrafish, chicken, and mouse (Dubrulle et al., 2015; Oginuma et al., 2017; Zagorski et al., 2017; Li et al., 2018; Petkova et al., 2019; Rogers et al., 2020; Oginuma et al., 2020).

Materials and methods

Methods are described in Appendices 1–3. Appendix 1 derives the expressions associated with the cost-precision trade-off relation for the localized synthesis-diffusion-depletion model of the morphogen dynamics in a 1D array of cells. Appendix 2 derives modified trade-off relations in a 1D model with a distributed morphogen source, and the morphogen dynamics occurring on a sphere. Appendix 3 details how the length and time scales of the naturally occurring morphogen profiles are obtained from literature.

Appendix 1

Reversible reaction-diffusion model of morphogen dynamics

Here, we provide detailed derivations of the cost-precision trade-off relation for the localized synthesis-diffusion-depletion model of the morphogen dynamics in a 1D array of cells. We begin with the reversible reaction-diffusion process of morphogen profile formation, where the forward and reverse rates are well defined at every elementary step, which enables the calculation of the total entropy production rate ( ${\dot{S}}_{t o t}$ ). Through the discussion of reversible process, we provide mathematical details on the derivation of the cost-precision trade-off relation for the irreversible reaction-diffusion process of morphogen synthesis and depletion presented in the main text. Next, we derive the conditions that lead to the minimum trade-off product for the irreversible process with point ( $π_{o}$ ) and space-time-averaged ( $π_{T}$ ) measurements.

Appendix 1—figure 1

Download asset Open asset

Speed-precision-cost trade-off for the reversible morphogen dynamics.

(a) Schematic of the reversible reaction-diffusion dynamics of morphogen profile formation. (b) The solid black line depicts the optimal trade-off product $π_{o, rev, \min} (x_{b})$ , obtained numerically by minimizing $π_{o, rev} (x_{b})$ with respect to λ and γ, with $β_{out} \to \infty$ . The trade-off product is normalized by $L / l_{cell}$ . (**c,d**) The optimal characteristic decay length $λ_{\min}$ (b) and $γ_{\min}$ (c) corresponding to $π_{o, rev, \min} (x_{b})$ .

We explore the trade-offs among the speed, precision, and entropy production in the reversible process of morphogen dynamics. Similarly to the irreversible version presented in the main text, the dynamics is defined on a 1D array of cells. We denote the amount of morphogen in the cell of volume ( $v_{c e l l}$ ) and size ( $l_{c e l l}$ ) at the interval between $x$ and $x + l_{cell}$ by the concentration, $ρ (x, t) [conc \equiv # / v_{cell}]$ . At the left boundary (see Appendix 1—figure 1a), the morphogen is injected with rate $j_{i n}$ $[conc \times l_{cell} / time]$ , and taken out by the corresponding reverse reaction with rate $v_{o u t}$ $[l_{cell} / time]$ . At all positions, the morphogen molecules are depleted with rate $k_{d} [{time}^{- 1}]$ and synthesized by the corresponding reverse reaction with rate $r_{s} [conc / time]$ , and spread across space with diffusivity $D [l_{cell}^{2} / time]$ . At the continuum limit ( $L ≫ l_{cell}$ ), the dynamics of the morphogen concentration is described by the partial differential equation

\partial_{t} ρ (x, t) = D \partial_{x}^{2} ρ (x, t) - k_{d} ρ (x, t) + r_{s},

with two boundary conditions at $x = 0$ and $x = L$ , $- D \partial_{x} ρ (x, t) |_{x = 0} = j_{i n} - v_{o u t} ρ (0, t)$ and $D \partial_{x} ρ (x, t) |_{x = L} = 0$ . $λ \equiv \sqrt{D / k_{d}}$ is the characteristic length of the dynamics associated with Equation 16. The steady state concentration is expressed as

ρ_{ss} (x) = \underset{= ρ_{ss, 0} (x)}{\underset{⏟}{\frac{r_{s}}{k_{d}} \frac{(γ - 1) \cosh [\frac{L - x}{λ}]}{β_{out}^{- 1} \sinh [\frac{L}{λ}] + \cosh [\frac{L}{λ}]}}} + \frac{r_{s}}{k_{d}},

where $ρ_{ss, 0} (x)$ denotes the position-dependent portion of the concentration profile. We introduce three dimensionless parameters: $β_{in} \equiv j_{in} / (λ r_{s})$ , $β_{out} \equiv v_{out} / (λ k_{d})$ , and their ratio $γ \equiv β_{in} / β_{out}$ which is related with the thermodynamic cost (see Equation 32). We only consider the regime with $γ > 1$ , where a net positive current of morphogen is supplied at the left boundary. In what follows, we will derive the quantitative expressions of the (i) speed, (ii) precision, and (iii) thermodynamic cost associated with the morphogen profile.

Speed

The concentration profile varies with time to reach its steady state, $ρ_{ss} (x) \equiv {lim}_{t \to \infty} ρ (x, t)$ , from the initial profile $ρ (x, 0) = 0$ . The time evolution of $ρ (x, t)$ can be analyzed by solving the Equation 16 in Laplace domain:

\hat{ρ} (x, s) = \frac{1}{s} [\frac{(j_{in} (s + k_{d}) - r_{s} v_{out}) \cosh [(L - x) \sqrt{(s + k_{d}) / D}]}{(s + k_{d}) \sinh [L \sqrt{(s + k_{d}) / D}] (D \sqrt{(s + k_{d}) / D} + v_{out} \coth [L \sqrt{(s + k_{d}) / D}])} + \frac{r_{s}}{(s + k_{d})}] .

The expression of the steady state profile (Equation 17) is obtained using

lim_{s \to 0} (s \hat{ρ} (x, s) - ρ (x, 0)) = lim_{s \to 0} \int_{0}^{\infty} 𝑑 t e^{- s t} \frac{d ρ (x, t)}{d t} = \int_{0}^{\infty} 𝑑 t \frac{d ρ (x, t)}{d t} = ρ (x, \infty) - ρ (x, 0) = ρ_{ss} (x) .

To characterize the relaxation dynamics to the steady state profile, we define the relaxation function

R (x, t) \equiv \frac{ρ_{ss} (x) - ρ (x, t)}{ρ_{ss} (x) - ρ (x, 0)} = 1 - \frac{ρ (x, t)}{ρ_{ss} (x)},

which monotonically decays from 1 to 0 at all positions as $t \to \infty$ . The associated mean relaxation time at position $x$ , $τ_{rev} (x)$ , is given by

τ_{rev} (x) = \int_{0}^{\infty} R (x, t) 𝑑 t .

Since the Laplace transform of the relaxation function is $\hat{R} (x, s) = \int_{0}^{\infty} e^{- s t} R (x, t) 𝑑 t = s^{- 1} - \hat{ρ} (x, s) / ρ_{ss} (x)$ , the local accumulation time at position $x$ , $τ_{rev} (x) = {lim}_{s \to 0} \hat{R} (x, s)$ , is obtained as

τ_{rev} (x) = 𝒩 / ℳ

with

\begin{array}{ll} 𝒩 & \equiv 2 β_{o u t} {(\coth [L / λ])}^{2} \\ + c s c h [\frac{L}{λ}] {[(γ - 3) + β_{o u t} (γ - 1) \frac{L}{λ}] c o s h [\frac{L - x}{λ}] + 2 β_{o u t}^{- 1} s i n h [\frac{L}{λ}] - (γ - 1) \frac{L - x}{λ} s i n h [\frac{L - x}{λ}]} \\ + c o t h [\frac{L}{λ}] c s c h [\frac{L}{λ}] ([γ \frac{L}{λ} - (2 β_{o u t} + \frac{L}{λ})] c o s h [\frac{L - x}{λ}] + 4 s i n h [\frac{L}{λ}] - β_{o u t} (γ - 1) \frac{L - x}{λ} s i n h [\frac{L - x}{λ}]) \end{array}

and

ℳ \equiv 2 k_{d} (1 + β_{out} \coth [\frac{L}{λ}]) (β_{out}^{- 1} + \coth [\frac{L}{λ}] + (γ - 1) \cosh [\frac{L - x}{λ}] csch [\frac{L}{λ}]) .

$τ_{rev}^{- 1} (x)$ quantifies the speed at which the morphogen profile is established at position $x$ (Berezhkovskii et al., 2010; Berezhkovskii et al., 2011). The expression for the characteristic time in the irreversible morphogen dynamics is obtained when $r_{s}$ and $v_{o u t}$ are negligibly small. At the limits of $r_{s} \to 0$ and $v_{out} \to 0$ , with the substitutions $β_{out} = v_{out} / (λ k_{d})$ and $γ = j_{in} k_{d} / (v_{out} r_{s})$ , Equation 22 is simplified to

\begin{array}{ll} τ (x) = lim_{\binom{r_{s} \to 0}{v_{o u t} \to 0}} τ_{r e v} (x) = \frac{1}{2 k_{d}} (1 + \frac{L}{λ} \coth [\frac{L}{λ}] - \frac{L - x}{λ} \tanh [\frac{L - x}{λ}]) \end{array} .

Precision

The local measure of precision at $x$ is defined as the squared relative positional error (Equation 1),

\begin{array}{ll} ϵ_{r e v}^{2} (x) & \equiv \frac{σ_{x}^{2} (x)}{x^{2}} = {(\frac{\partial \log ρ_{s s} (x)}{\partial \log x})}^{- 2} \frac{σ_{ρ}^{2} (x)}{ρ_{s s} (x)^{2}} \\ = \frac{λ^{2} k_{d}}{v_{c e l l} x^{2} r_{s}} \frac{(β_{o u t}^{- 1} + c o t h [\frac{L}{λ}]) (β_{o u t}^{- 1} + c o t h [\frac{L}{λ}] + (γ - 1) c o s h [\frac{L - x}{λ}] c s c h [\frac{L}{λ}]) {(c s c h [\frac{L - x}{λ}] s i n h [\frac{L}{λ}])}^{2}}{{(γ - 1)}^{2}} . \end{array}

The expression for the squared relative error in the irreversible process of morphogen dynamics can be obtained at the limits of $r_{s} \to 0$ and $v_{out} \to 0$ , with the substitutions $β_{out} = v_{out} / (λ k_{d})$ and $γ = j_{in} k_{d} / (v_{out} r_{s})$ :

ϵ^{2} (x) = lim_{\binom{r_{s} \to 0}{v_{o u t} \to 0}} ϵ_{r e v}^{2} (x) = \frac{λ^{3}}{l_{c e l l} x^{2}} {(\frac{v_{c e l l} j_{i n}}{l_{c e l l} k_{d}})}^{- 1} s i n h [\frac{L}{λ}] c o s h [\frac{L - x}{λ}] {(s i n h [\frac{L - x}{λ}])}^{- 2} .

Thermodynamic cost

The total entropy production rate required to maintain the steady state morphogen profile is given by

{\dot{S}}_{t o t} = \frac{1}{l_{c e l l}} \int_{0}^{L} d x [{\dot{S}}_{i n} (x) + {\dot{S}}_{d e p} (x) + {\dot{S}}_{d i f f} (x)],

where the three components of the position-dependent entropy production rate are

{\dot{S}}_{in} (x) \equiv 2 v_{cell} (j_{in} - v_{out} ρ_{ss} (0)) δ (x) \log (\frac{j_{in}}{v_{out} ρ_{ss} (0)}),

{\dot{S}}_{dep} (x) \equiv v_{cell} (k_{d} ρ_{ss} (x) - r_{s}) \log (\frac{k_{d} ρ_{ss} (x)}{r_{s}}),

{\dot{S}}_{diff} (x) \equiv v_{cell} D [\partial_{x} ρ_{ss} (x)] [\partial_{x} \ln ρ_{ss} (x)] = v_{cell} D \frac{{[\partial_{x} ρ_{ss} (x)]}^{2}}{ρ_{ss} (x)},

with the Boltzmann constant set to unity ( $k_{B} = 1$ ). Equations 27 and 28 represent the entropy production rates from the net morphogen current flowing in and out of the system, respectively. Similarly, Equation 29 arrises from the diffusive flux (Falasco et al., 2018). Briefly, Equation 29 can be understood as $lim_{l_{cell} \to 0} D l_{cell}^{- 2} (ρ_{ss} (x + l_{cell}) - ρ_{ss} (x)) \ln [ρ_{ss} (x + l_{cell}) / ρ_{ss} (x)]$ , in which $D ρ_{ss} (x + l_{cell}) / l_{cell}^{2}$ and $D ρ_{ss} (x) / l_{cell}^{2}$ amount to the forward and reverse currents between adjacent spatial compartments in the discretized spatial representation. It is noteworthy that the asymmetry of concentrations generates a flow of morphogen current between the neighboring cells, and contributes to the entropy production. If the morphogen profile is uniform in space, which occurs at the detailed balance condition ( ${\dot{S}}_{tot} = 0$ ), the cells cannot infer any spatial information from it.

To obtain a simplified expression of ${\dot{S}}_{tot}$ , we use the following two equalities. First, at steady state, the net injection of the morphogen at the left-boundary must be equal to the net depletion of the morphogen occurring at all positions, so that

\int_{0}^{L} 2 δ (x) [j_{in} - v_{out} ρ_{ss} (0)] 𝑑 x = \int_{0}^{L} [k_{d} ρ_{ss} (x) - r_{s}] 𝑑 x .

Second, $\partial_{x} ρ_{ss} (x) = \partial_{x} ρ_{ss, 0} (x)$ and $\partial_{x}^{2} ρ_{ss} (x) = \partial_{x}^{2} ρ_{ss, 0} (x) = λ^{- 2} ρ_{ss, 0} (x)$ , with $ρ_{ss, 0}$ defined in Equation 17, lead to the following equality:

\int_{0}^{L} [ρ_{ss, 0} (x) \ln ρ_{ss} (x)] 𝑑 x = λ^{2} ({[\partial_{x} ρ_{ss, 0} (x)] \ln ρ_{ss} (x) |}_{0}^{L} - \int_{0}^{L} [\frac{{[\partial_{x} ρ_{ss} (x)]}^{2}}{ρ_{ss} (x)}] 𝑑 x) .

Then, ${\dot{S}}_{tot}$ can be expressed as follows:

\begin{array}{ll} {\dot{S}}_{t o t} & = \frac{1}{l_{c e l l}} \int_{0}^{L} ({\dot{S}}_{i n} (x) + {\dot{S}}_{d e p} (x) + {\dot{S}}_{d i f f} (x)) d x \\ = \frac{v_{c e l l} k_{d}}{l_{c e l l}} \int_{0}^{L} (ρ_{s s, 0} (x) \ln \frac{j_{i n} k_{d}}{v_{o u t} r_{s}} + ρ_{s s, 0} (x) \ln \frac{ρ_{s s} (x)}{ρ_{s s} (0)} + λ^{2} \frac{[\partial_{x} ρ_{s s} (x)]^{2}}{ρ_{s s} (x)}) d x \\ = \frac{v_{c e l l} k_{d}}{l_{c e l l}} (\int_{0}^{L} ρ_{s s, 0} (x) \ln γ d x - \int_{0}^{L} ρ_{s s, 0} (x) \ln ρ_{s s} (0) d x + \int_{0}^{L} [ρ_{s s, 0} (x) \ln ρ_{s s} (x) + λ^{2} \frac{[\partial_{x} ρ_{s s} (x)]^{2}}{ρ_{s s} (x)}] d x) \\ = \frac{v_{c e l l} k_{d}}{l_{c e l l}} (\int_{0}^{L} ρ_{s s, 0} (x) \ln γ d x - λ^{2} [\partial_{x} ρ_{s s, 0} (x)] \ln ρ_{s s} (0) |_{0}^{L} + λ^{2} [\partial_{x} ρ_{s s, 0} (x)] \ln ρ_{s s} (x) |_{0}^{L}) \\ = \frac{v_{c e l l} k_{d}}{l_{c e l l}} \ln γ \int_{0}^{L} ρ_{s s, 0} (x) d x \\ = \frac{v_{c e l l} λ r_{s}}{l_{c e l l}} \frac{(γ - 1) \ln γ}{β_{o u t}^{- 1} + c o t h [\frac{L}{λ}]} . \end{array}

We used 30 and 31 to obtain the second and fourth lines of Equation 32, respectively. The boundary condition of the system, ${\partial_{x} ρ (x) |}_{x = L} = 0$ , and the integral of $ρ_{ss, 0} (x)$ leads to the last line.

The reversible model motivates a physical interpretation of the proportionality constant $α_{o}$ , which we previously defined in the irreversible model of morphogen dynamics (Equation 2) as the number of ATPs hydrolyzed for the synthesis and depletion of a morphogen molecule (Equation 5). When ${\dot{S}}_{tot}$ (Equation 32) is normalized by $\ln γ$ , and $r_{s}$ and $v_{o u t}$ are negligibly small, we obtain

lim_{\binom{r_{s} \to 0}{v_{o u t} \to 0}} \frac{{\dot{S}}_{t o t}}{\ln γ} = \frac{v_{c e l l} j_{i n}}{l_{c e l l}} .

Comparison of Equation 33 with the expression in Equation 5 enables us to relate $α_{o}$ with the free energy cost, $T \ln γ$ . Thus, the thermodynamic cost associated with the morphogen profile is expressed as $C = (α_{o} v_{cell} j_{in} / l_{cell}) τ$ (Equation 5) where τ is the time required to generate the profile, and $α_{o} v_{cell} j_{in} / l_{cell}$ approximates the rate of dissipation, $T {\dot{S}}_{tot}$ .

Trade-off product

The product of the local accumulation time ( $τ_{rev} (x)$ , Equation 22), the squared relative positional error ( $ϵ_{rev}^{2} (x)$ , Equation 24), and the entropy production rate ( ${\dot{S}}_{tot}$ , Equation 32) leads to the position-dependent trade-off product $π_{o, rev} (x)$ ,

\begin{array}{ll} π_{o, r e v} (x; λ, β_{o u t}, γ) & = {\dot{S}}_{t o t} τ_{r e v} (x) ϵ_{r e v}^{2} (x) \\ = \frac{1}{2} \frac{L}{l_{c e l l}} \frac{λ^{3}}{L x^{2}} \frac{c s c h [(L - x) / λ] s i n h [L / λ] \ln γ}{(γ - 1) (β_{o u t} c o t h [L / λ] + 1)} \\ \times {(1 - γ) \frac{L - x}{λ} + (γ - 3 + β_{o u t} (γ - 1) \frac{L}{λ}) c o t h [\frac{L - x}{λ}] + 4 c o s h [\frac{L}{λ}] c s c h [\frac{L - x}{λ}] \\ + c o t h [\frac{L}{λ}] (β_{o u t} (1 - γ) \frac{L - x}{λ} + ((γ - 1) \frac{L}{λ} - 2 β_{o u t}) c o t h [\frac{L - x}{λ}] \\ + 2 β_{o u t} c o s h [\frac{L}{λ}] c s c h [\frac{L - x}{λ}]) \\ + 2 β_{o u t}^{- 1} c s c h [\frac{L - x}{λ}] s i n h [\frac{L}{λ}]} . \end{array}

With fixed $L$ , $l_{c e l l}$ , and $x = x_{b}$ , the lower bound of $π_{o, rev} (x_{b}; λ, β_{out}, γ)$ , $π_{rev, \min} (x_{b})$ , can be obtained by minimizing $π_{o, rev} (x_{b}; λ, β_{out}, γ)$ with respect to γ, $β_{out}$ , and λ. For a fixed value of $γ > 1$ , the partial derivative of $π_{o, rev} (x_{b}; λ, β_{out}, γ)$ with respect to $β_{out}$ is always smaller than 0,

\begin{array}{ll} \partial_{β_{o u t}} π_{o, r e v} (x_{b}; λ, β_{o u t}, γ) \\ = - \frac{L}{l_{c e l l}} \frac{λ^{3}}{L (x_{b})^{2}} \ln γ {(c s c h [\frac{L - x_{b}}{λ}])}^{2} \\ \times (\frac{{(s i n h [L / λ])}^{2}}{β_{o u t}^{2} (γ - 1)} + \frac{c o s h [(L - x_{b}) / λ] s i n h [L / λ] (2 L / λ + \sinh [2 L / λ])}{4 {(β_{o u t} c o s h [L / λ] + s i n h [L / λ])}^{2}}) < 0. \end{array}

Thus, $π_{o, rev} (x_{b})$ is minimized at the limit of large $β_{out}$ (or, equivalently, large $β_{in} = γ β_{out}$ ), with the following expression:

\begin{array}{cc} lim_{β_{o u t} \to \infty} π_{o, r e v} (x_{b}; λ, β_{o u t}, γ) = \frac{1}{4} \frac{L}{l_{c e l l}} \frac{λ^{3}}{L (x_{b})^{2}} \frac{{(c s c h [(L - x_{b}) / λ])}^{2} t a n h [L / λ] \ln γ}{(γ - 1)} \\ \times {2 + 2 c o s h [\frac{2 L}{λ}] - 2 c o s h [\frac{2 L - x_{b}}{λ}] - 2 c o s h [\frac{x_{b}}{λ}] \\ + (γ - 1) (\frac{x_{b}}{λ} s i n h [\frac{2 L - x_{b}}{λ}] + (\frac{2 L - x_{b}}{λ}) s i n h [\frac{x_{b}}{λ}])} . \end{array}

With the characteristic decay length set to $λ = x_{b}$ , it can be shown that $lim_{\binom{x_{b} \to 0}{β_{o u t} \to \infty}} π_{o, r e v} (x_{b}; x_{b}, β_{o u t}, γ) = 0$ . Since $π_{o, rev, \min} (x_{b}) \leq π_{o, rev} (x_{b}; x_{b}, β_{out}, γ)$ by the definition of $π_{o, rev, \min} (x_{b})$ , $π_{o, rev, \min} (0) = 0$ .

For $x_{b} > 0$ , the optimal trade-off product, $π_{rev, \min} (x_{b})$ , and the corresponding values of $λ_{\min} (x_{b})$ and $γ_{\min} (x_{b})$ , were numerically obtained (Appendix 1—figure 1). $π_{o, rev, \min} (x_{b})$ constrains dynamical properties of the morphogen profile, specifying, for instance, the minimum entropy production rate for a given speed and precision (Appendix 1—figure 1b). With a small $π_{o, rev} (x_{b}; λ, β_{out}, γ)$ close to $π_{o, rev, \min} (x_{b})$ , the system can rapidly generate precise morphogen profiles at $x_{b}$ with low thermodynamic cost. In general, $π_{o, rev, \min} (x_{b})$ , $λ_{\min} (x_{b})$ , and $γ_{\min} (x_{b})$ all increase monotonically with $x_{b}$ (Appendix 1—figure 1b,c,d).

A cautionary remark is in place. A direct one-to-one comparison between the trade-off products of the reversible and irreversible models should be avoided, since $π_{o, rev}$ is minimized at the limits of $v_{out}, r_{s} \to \infty$ (see Equation 35), while the irreversible model is equivalent to the reversible version at the limits of $v_{out}, r_{s} \to 0$ .

Minimum trade-off product of the point measurement for the irreversible reaction-diffusion process

We provide the mathematical details involving the trade-off product of the point measurement ( $π_{o}$ , Equation 7 and Figure 2 in the main text). With the expressions derived in the Appendix 1, Reversible reaction-diffusion model of morphogen dynamics, the thermodynamic cost of generating the steady state morphogen concentration at position $x$ (Equation 5) amounts to the product among $α_{o}$ , $τ (x)$ (Equation 23), and $v_{cell} j_{in} / l_{cell}$ (Equation 33),

𝒞 (x) = α_{o} (\frac{v_{cell} j_{in}}{l_{cell} k_{d}}) \frac{1}{2} (1 + \frac{L}{λ} \coth [\frac{L}{λ}] - \frac{L - x}{λ} \tanh [\frac{L - x}{λ}]),

and the squared relative positional error, which is identical to Equation 25, is

ϵ^{2} (x) = \frac{λ^{3}}{l_{cell} x^{2}} {(\frac{v_{cell} j_{in}}{l_{cell} k_{d}})}^{- 1} \frac{\sinh [\frac{L}{λ}] \cosh [\frac{L - x}{λ}]}{{(\sinh [\frac{L - x}{λ}])}^{2}} .

Then, the trade-off product reads

\begin{array}{ll} π_{o} (x; λ) & = 𝒞 (x) ϵ^{2} (x) \\ = \frac{α_{o}}{2} \frac{λ^{3}}{l_{c e l l} x^{2}} (1 + \frac{L}{λ} \coth [\frac{L}{λ}] - \frac{L - x}{λ} \tanh [\frac{L - x}{λ}]) \frac{\sinh [\frac{L}{λ}] \cosh [\frac{L - x}{λ}]}{{(s i n h [\frac{L - x}{λ}])}^{2}} . \end{array}

For fixed $x = x_{b}$ , $π_{o} (x_{b}; λ)$ is numerically minimized with respect to λ to give rise to the solid black line in Figure 2c in the main text. When $L / λ ≫ 1$ , the trade-off product (Equation 39) is approximated to

π_{o} (x; λ) \approx \frac{α_{o}}{2} \frac{λ}{l_{cell}} [(\frac{λ}{x}) + {(\frac{λ}{x})}^{2}] e^{x / λ} .

For a given position $x = x_{b}$ , ${\partial_{λ} π_{o} (x_{b}; λ) |}_{λ = λ_{\min} (x_{b})} = 0$ determines $λ_{\min} (x_{b}) \approx (\sqrt{13} - 1) x_{b} / 6 \approx 0.43 x_{b}$ (gray dotted line in Figure 2c of the main text).

Minimum trade-off product of the space-time-averaged measurement for the irreversible reaction-diffusion process

As in the main text, the trade-off product from the space-time-averaged measurement is given by

\begin{array}{ll} π_{T} (x; λ) & = 𝒞_{T} ϵ_{T}^{2} (x) \\ = α_{o} j_{i n} T {(\frac{\partial \log ⟨ m (x) ⟩}{\partial \log x})}^{- 2} \frac{σ_{m}^{2} (x)}{⟨ m (x) ⟩^{2}} \\ = α_{o} j_{i n} T \underset{= [x \partial_{x} ⟨ m (x) ⟩]^{- 2}}{\underset{⏟}{\frac{λ^{2}}{x^{2} {[\frac{2 j_{i n}}{k_{d}} e^{- x / λ} \sinh (a / λ)]}^{2}}}} \\ \times \underset{= σ_{m}^{2} (x)}{\underset{⏟}{\frac{4 j_{i n} e^{- x / λ}}{T k_{d}^{2}} \sinh (a / λ) (1 - \frac{λ \cosh (a / λ) \sinh (a / λ) + a}{2 e^{a / λ} λ \sinh (a / λ)})}} \\ = α_{o} {(\frac{λ}{x})}^{2} \frac{e^{x / λ}}{\sinh (a / λ)} (1 - \frac{λ \cosh (a / λ) \sinh (a / λ) + a}{2 e^{a / λ} λ \sinh (a / λ)}) \\ = α_{o} {(\frac{λ}{x})}^{2} e^{x / λ} \frac{e^{- a / λ} + 3 e^{3 a / λ} - 4 (a λ^{- 1} + 1) e^{a / λ}}{2 {(e^{2 a / λ} - 1)}^{2}}, \end{array}

The $π_{T} (x; λ)$ minimizing λ at $x = x_{b}$ , namely $λ_{\min} (x_{b})$ , is found at the interval $(x_{b} - a) / 2 \leq λ_{\min} (x_{b}) \leq x_{b} / 2$ , since

{\partial_{λ} π_{T} (x_{b}; λ) |}_{λ = x_{b} / 2} = α_{o} \frac{a e^{2 + 2 a / x_{b}} (1 + 3 e^{4 a / x_{b}})}{x_{b}^{2} {(e^{4 a / x_{b}} - 1)}^{3}} (\sinh [\frac{4 a}{x_{b}}] - \frac{4 a}{x_{b}}) \geq 0,

and

{\partial_{λ} π_{T} (x_{b}; λ) |}_{λ = (x_{b} - a) / 2} = α_{o} \frac{8 a^{2} e^{2 (x_{b} + a) / (x_{b} - a)}}{x_{b}^{2} (x_{b} - a) {(e^{4 a / (x_{b} - a)} - 1)}^{2}} ({[\frac{2 a}{x_{b} - a}]}^{- 1} - \coth [\frac{2 a}{x_{b} - a}]) \leq 0 .

For Figure 3 in the main text, the exact values of $π_{T, \min} (x_{b})$ and $λ_{\min} (x_{b})$ were obtained by numerically minimizing Equation 41 with respect to λ.

Appendix 2

Distributed morphogen source

Appendix 2—figure 1

Download asset Open asset

Cost-precision trade-off in the morphogen dynamics with distributed synthesis.

(a) Schematic of the distributed synthesis of the morphogen. The black solid line represents the morphogen synthesis rate which is roughly exponentially distributed with characteristic decay length $λ_{p}$ . (**b–d**) The trade-off product ( $π_{o} (x, λ_{d})$ ) for (b) $λ_{p} = 0$ , (c) $λ_{p} = 0.05 L$ , and (d) $λ_{p} = 0.1 L$ . The color scale indicates the trade-off product $π_{o}$ normalized by $α_{o} L / l_{cell}$ . The characteristic length scale $λ_{d}$ is defined as $λ_{d} \equiv \int_{0}^{L} 𝑑 x (ρ_{ss} (x) - ρ_{ss} (L)) / (ρ_{ss} (0) - ρ_{ss} (L))$ .

We consider a case where the source of morphogen production is distributed over the space in the irreversible reaction-diffusion model of morphogen dynamics. The space-time evolution of the morphogen concentration is described by

\partial_{t} ρ (x, t) = D \partial_{x}^{2} ρ (x, t) - k_{d} ρ (x, t) + J_{in} (x),

with the boundary conditions ${\partial_{x} ρ (x, t) |}_{x = 0, L} = 0$ . The production term, $J_{in} (x)$ is assumed to be exponentially distributed with a characteristic length $λ_{p}$ ,

J_{in} (x) = \frac{j_{in}}{λ_{p}} \frac{\cosh [\frac{L - x}{λ_{p}}]}{\sinh [\frac{L}{λ_{p}}]},

so that $\int_{0}^{L} J_{in} (x) 𝑑 x = j_{in}$ . The solution to Equation 44 in the Laplace domain is given by

\hat{ρ} (x, s) = \frac{j_{in}}{s \sqrt{D (k_{d} + s)}} \frac{D}{D - (k_{d} + s) λ_{p}^{2}} (\frac{\cosh [\frac{L - x}{\sqrt{D / (k_{d} + s)}}]}{\sinh [\frac{L}{\sqrt{D / (k_{d} + s)}}]} - \frac{λ_{p} \sqrt{k_{d} + s}}{\sqrt{D}} \frac{\cosh [\frac{L - x}{λ_{p}}]}{\sinh [\frac{L}{λ_{p}}]}) .

By following a procedure similar to that shown in the section Reversible reaction-diffusion model of morphogen dynamics, we can obtain the following expressions for the steady state profile ( $ρ_{ss} (x)$ , Equation 19), and the precision ( $ϵ^{2} (x)$ , Equation 24):

ρ_{ss} (x) = \frac{j_{in}}{\sqrt{D k_{d}}} \frac{λ^{2}}{λ^{2} - λ_{p}^{2}} (\frac{\cosh [\frac{L - x}{λ}]}{\sinh [\frac{L}{λ}]} - \frac{λ_{p}}{λ} \frac{\cosh [\frac{L - x}{λ_{p}}]}{\sinh [\frac{L}{λ_{p}}]}),

ϵ^{2} (x) = \frac{ρ_{ss} (x)}{v_{cell} {(x \partial_{x} ρ_{ss} (x))}^{2}} = \frac{k_{d} λ}{v_{cell} j_{in}} \frac{(λ^{2} - λ_{p}^{2}) (\frac{\cosh [\frac{L - x}{λ}]}{\sinh [\frac{L}{λ}]} - \frac{λ_{p}}{λ} \frac{\cosh [\frac{L - x}{λ_{p}}]}{\sinh [\frac{L}{λ_{p}}]})}{x^{2} {[(\frac{\sinh [\frac{L - x}{λ}]}{\sinh [\frac{L}{λ}]} - \frac{\sinh [\frac{L - x}{λ_{p}}]}{\sinh [\frac{L}{λ_{p}}]})]}^{2}},

with $λ \equiv \sqrt{D / k_{d}}$ . Next, the characteristic timescale to approach the steady state profile ( $τ (x)$ , Equation 22) can be written as

τ (x) = 𝒩 / ℳ,

with

𝒩 \equiv (1 - 3 \frac{λ_{p}^{2}}{λ^{2}} + \frac{L}{λ} (1 - \frac{λ_{p}^{2}}{λ^{2}}) c o t h [\frac{L}{λ}]) \frac{c o s h [\frac{L - x}{λ}]}{s i n h [\frac{L}{λ}]} + 2 \frac{λ_{p}^{3}}{λ^{3}} \frac{c o s h [\frac{L - x}{λ_{p}}]}{s i n h [\frac{L}{λ_{p}}]} - \frac{L - x}{λ} (1 - \frac{λ_{p}^{2}}{λ^{2}}) \frac{s i n h [\frac{L - x}{λ}]}{s i n h [\frac{L}{λ}]},

and

ℳ \equiv 2 k_{d} (- 1 + \frac{λ_{p}^{2}}{λ^{2}}) (\frac{λ_{p}}{λ} c o s h [\frac{x}{λ_{p}}] c o t h [\frac{L}{λ_{p}}] - \frac{λ_{p}}{λ} s i n h [\frac{x}{λ_{p}}] - \frac{c o s h [\frac{L - x}{λ}]}{s i n h [\frac{L}{λ}]}) .

Then, the trade-off product is expressed as

π (x; λ, λ_{p}) = α_{o} (\frac{v_{cell} j_{in}}{l_{cell}}) τ (x) ϵ^{2} (x),

with the proportionality constant $α_{o}$ . For fixed $x = x_{b}$ , $L$ , and $λ_{p}$ , the trade-off product can be minimized by tuning λ to the optimal value, $λ_{\min} (x_{b})$ . As $λ_{p}$ increases, the optimal trade-off product and the associated optimal characteristic length $λ_{\min} (x)$ decrease (Appendix 2—figure 1). In Appendix 2—figure 1, the characteristic length scale $λ_{d}$ is defined as $λ_{d} \equiv \int_{0}^{L} 𝑑 x (ρ_{ss} (x) - ρ_{ss} (L)) / (ρ_{ss} (0) - ρ_{ss} (L))$ .

Dynamics on a sphere

Appendix 2—figure 2

Download asset Open asset

Cost-precision trade-off of the morphogen dynamics on a sphere.

(a) (Top) Schematic of the morphogen synthesized from the equator of the spherical shell (red region). (Bottom) The relative steady state concentration profile $ρ_{ss} (z)$ when $R^{2} / λ^{2} = 10$ . The solid black line was computed from evaluating $ρ_{ss} (z) = {lim}_{s \to 0} s \hat{ρ} (z, s)$ up to $l = 84$ (Equation 53). The dotted red line was obtained by numerically solving Equation 51. (b) The trade-off product ( $π_{o} (z, λ_{s})$ ) for pairs of target position $z$ and characteristic length scale $λ_{s}$ values, normalized by $α_{o} (2 π R / l_{cell})$ . The characteristic length scale $λ_{s}$ is defined as $λ_{s} \equiv \int_{0}^{1} 𝑑 z (ρ_{ss} (z) - ρ_{ss} (1)) / (ρ_{ss} (0) - ρ_{ss} (1))$ . The parameters for the naturally occurring morphogen profiles are further described in Appendix 3, Length scales of Cyclops, Squint, and Fgf8, and Appendix 3—table 1.

Appendix 2—figure 3

Download asset Open asset

Length scales of naturally occurring morphogen profiles in the zebrafish embryo.

(a) (Left) Schematic of the Nodal ligands Cyclops and Squint patterning the vegetal-animal axis of the zebrafish embryo. The Nodal signals activate Smad2, which affects the expression of short-range and long-range target genes. (Right) The expression profile of Nodal target genes *gsc* and *ntl* obtained by in situ hybridization. The profiles are from the zebrafish embryos 5 hr post fertilization, when the embryo covers roughly 50% of the yolk. The red lines denote the width of the target gene expression domain away from the vegetal margin. (b) (Left) Expression of *fgf8* in 60% epiboly whole-mount embryos obtained by in situ hybridization. Animal side of the embryo is to the top, and dorsal is to the right. (Right) The expression profiles of the Fgf8 target genes *spry4*, *spry2*, and *pea3*. The red lines denote the width of the target gene expression domain away from the vegetal margin. Scale bars = 100 $μ m$ . (A) is reproduced from Figure 2 Dubrulle et al., 2015. (B) is reproduced from Figure 1 Nowak et al., 2011. Reprinted by permission from Springer Nature: Springer Nature, Nature Cell Biology, Interpretation of the FGF8 morphogen gradient is regulated by endocytic trafficking, Nowak, M., MacHate, A., Yu, S. R., Gupta, M., and Brand, M. 2011, Nat. Cell Biol., 13(2):153–158 Copyright 2011. It is not covered by the CC-BY 4.0 licence and further reproduction of this panel would need permission from the copyright holder.

© 2015, Dubrulle et al. Appendix 2—figure 3A is reproduced from Figure 2 Dubrulle et al. 2015, with permission from The Company of Biologists Creative Commons Attribution 4.0 International Public License.

© 2011, Springer Nature. Appendix 2—figure 3B is reproduced from Figure 1 Nowak et al., 2011. Reprinted by permission from Springer Nature: Nowak et al., 2011 Copyright 2011. It is not covered by the CC-BY 4.0 licence and further reproduction of this panel would need permission from the copyright holder.Creative Commons Attribution 4.0 International Public License.

We consider the morphogen profile formation on the surface of a sphere. The profile formations of Cyclops, Squint, and Fgf8 in the developing zebrafish embryo at ∼50% epiboly pertain to this situation. When the concentration of the morphogen in a cell with spherical coordinates $(r, θ, ϕ)$ is denoted by $ρ (r, θ, ϕ)$ , with the symmetry along the azimuthal angle and the condition of the morphogen dynamics being confined to the spherical shell of radius $r = R$ (i.e. $\partial_{ϕ} ρ = 0$ , $\partial_{r} ρ = 0$ , and $r = R$ ) the morphogen concentration satisfies the following evolution equation:

\frac{\partial ρ (θ, t)}{d t} = D [\frac{1}{R^{2} \sin θ} \frac{\partial}{\partial θ} (\sin θ \frac{\partial ρ (θ, t)}{\partial θ})] - k_{d} ρ (θ, t) + j_{in} δ (\cos θ),

with boundary conditions at ${\partial_{θ} ρ |}_{θ = 0, π} = 0$ . The morphogen dynamics is fully specified by the system size $R [l_{c e l l}]$ , the depletion rate $k_{d} [1 / t i m e]$ , the diffusion coefficient $D [l_{c e l l}^{2} / t i m e]$ , and the injection rate $j_{i n} [c o n c / t i m e]$ . The last term represents the influx of morphogens due to their synthesis by the ‘band’ of cells at the equator (Appendix 2—figure 2a). With the substitution $z = \cos θ$ and the initial condition $ρ (z, t) = 0$ , the solution of Equation 51 in the Laplace domain with $\hat{ρ} (z, s) = \int_{0}^{\infty} 𝑑 t e^{- s t} ρ (z, t)$ , can be obtained by solving

\frac{\partial}{\partial z} ((1 - z^{2}) \frac{\partial \hat{ρ} (z, s)}{\partial z}) - \frac{R^{2} (k_{d} + s)}{D} \hat{ρ} (z, s) + \frac{R^{2} j_{in}}{D s} δ (z) = 0 .

The solution to the differential equation of this type can be expressed as $\hat{ρ} (z, s) = \sum_{l = 0}^{\infty} c_{l} (s) P_{l} (z)$ , where $P_{l} (z)$ ( $l = 0, 1, \dots$ ) is the Legendre polynomial, satisfying $ℒ P_{l} (x) \equiv \frac{d}{d x} [(1 - x^{2}) \frac{d P_{l} (x)}{d x}] = - l (l + 1) P_{l} (x)$ , with the orthogonality condition $\int_{- 1}^{1} P_{l} (x) P_{m} (x) 𝑑 x = \frac{2}{2 l + 1} δ_{l m}$ . To find the coefficient $c_{l} (s)$ , we rearrange Equation 52 into

- \frac{R^{2} j_{in}}{D s} δ (z) = - \sum_{l = 0}^{\infty} (l (l + 1) + \frac{R^{2} (k_{d} + s)}{D}) c_{l} (s) P_{l} (z) .

By using the orthogonality condition of Legendre polynomials, we obtain the solution of Equation 52,

\hat{ρ} (z, s) = \frac{R^{2} j_{in}}{D s} \sum_{l = 0}^{\infty} \frac{(2 l + 1)}{2 (\frac{R^{2} (k_{d} + s)}{D} + l (l + 1))} P_{l} (0) P_{l} (z) .

The steady state concentration can be expressed as $ρ_{ss} (z) = {lim}_{s \to 0} s \hat{ρ} (z, s)$ .

The average relaxation time can be calculated using $τ (z) \equiv ⟨ t ⟩ = \int_{0}^{\infty} t (- \frac{d R (t)}{d t}) 𝑑 t = \int_{0}^{\infty} R (t) 𝑑 t = {lim}_{s \to 0} \hat{R} (z, s)$ , which yields

τ (z) = k_{d}^{- 1} 𝒯 (z, λ / R),

with $λ \equiv \sqrt{D / k_{d}}$ and

𝒯 (z, λ / R) = \frac{R^{2}}{λ^{2}} \frac{\sum_{l = 0}^{\infty} \frac{(2 l + 1)}{2 {(R^{2} / λ^{2} + l (l + 1))}^{2}} P_{l} (0) P_{l} (z)}{\sum_{l = 0}^{\infty} \frac{(2 l + 1)}{2 (R^{2} / λ^{2} + l (l + 1))} P_{l} (0) P_{l} (z)} .

The number of morphogen molecules produced from the cells at the equator per unit time is $v_{cell} (2 π R / l_{cell}) j_{in}$ . Thus, the thermodynamic cost to approach the steady state concentration at $z$ is

C (z) = α_{o} 2 π R \frac{v_{cell} j_{in}}{l_{cell} k_{d}} 𝒯 (z, λ / R),

where $α_{o}$ is the proportionality constant.

Next, using $ρ_{ss} (z) = {lim}_{s \to 0} s \hat{ρ} (z, s)$ (Equation 53), one obtains the local precision of the morphogen profile at position $z$ ,

ϵ^{2} (z) = \frac{ρ_{ss} (z)}{v_{cell} {(z \partial_{z} ρ_{ss} (z))}^{2}} = \frac{k_{d}}{v_{cell} j_{in}} ℰ (z, λ / R),

where

ℰ (z, λ / R) \equiv \frac{\frac{λ^{2}}{R^{2}} \sum_{l = 0}^{\infty} \frac{(2 l + 1)}{2 (R^{2} / λ^{2} + l (l + 1))} P_{l} (0) P_{l} (z)}{z^{2} {(\sum_{l = 0}^{\infty} \frac{(2 l + 1)}{2 (R^{2} / λ^{2} + l (l + 1))} P_{l} (0) \partial_{z} P_{l} (z))}^{2}} .

Finally, the trade-off product, which is plotted in Appendix 2—figure 2b. Eb as a function of $z$ and $λ_{s}$ , is defined as

π (z; λ / R) \equiv C (z) \times ϵ^{2} (z) = α_{o} \frac{2 π R}{l_{cell}} 𝒯 (z, λ / R) ℰ (z, λ / R) .

The length scale $λ_{s}$ used in the plot is defined as $λ_{s} \equiv \int_{z = 0}^{z = 1} 𝑑 z (ρ_{ss} (z) - ρ_{ss} (1)) / (ρ_{ss} (0) - ρ_{ss} (1))$ . At the target position $z = z_{b}$ , the trade-off product can be minimized by tuning λ to the optimal value $λ_{\min}$ , which is plotted in black in Appendix 2—figure 2b.

Appendix 3

Length scales of Bcd, Wg, Hh, and Dpp

Here we explain the length scales associated with the naturally occurring morphogen profiles of Drosophila (Figure 2—figure supplement 1, Appendix 3—table 1).

Bcd. The characteristic decay length of Bcd was set to $λ_{B c d} = 100 μ m$ (Houchmandzadeh et al., 2002). The length of the embryo across the anterior-posterior axis is 475 µm, and the Bcd profile decays exponentially starting at 50 µm from the anterior end (Houchmandzadeh et al., 2002). Thus, we set the system length to $L = 425 μ m$ . The sensor radius $a = 2.6 μ m$ was obtained from the image of nuclei in Figure 3A of Gregor et al., 2007. The Bcd target Hb expression displays a sharp decline in the range of (190–238) µm (Figure 3C of Perry et al., 2012). Since the Bcd profile decays exponentially starting at 50 µm from the anterior end, we set $x_{b} = (140 - 188) μ m$ .
Wg. The characteristic decay length of Wg was set to $λ_{W g} = 6 μ m$ (Kicheva et al., 2007). The overall size of the DV axis patterned by Wg was set to $L = 70 μ m$ , based on the Hoeschst staining image in Figure 1E of Chaudhary et al., 2019. Although the Wg dynamics reported in Kicheva et al., 2007 is not from the endogenous concentration profile, we assume that the wild type profile forms with similar length and time scales. While Wg induces the expression of multiple genes, including distalless and vestigial, we picked sens for our analysis because its expression profile was quantitatively shown to undergo a sharp transition. sens expression decreases rapidly at (15–25) µm away from the DV boundary (Figure 1G of Bakker et al., 2020). Assuming that Wg is secreted by the two stripes of cells with radius 2 µm at the DV boundary (Figure 2A from Chaudhary et al., 2019), we set the target boundary to $x_{b} = (11 - 21) μ m$ .
Hh. The characteristic decay length of Hh was set to $λ_{H h} = 8 μ m$ (Figure S2B of Wartlick et al., 2011). The overall size of the domain patterned by Hh was set to $L = 100 μ m$ (Figure 6C of Torroja et al., 2004). The Hh target genes include en, col, and dpp, whose expression profiles were obtained from Figure 5 of Torroja et al., 2004. The expression domains of en, col, and dpp, respectively, span 15 µm, 25 µm, and 30 µm away from the Hh producing cells. Accordingly, the range of boundary positions associated with Hh was set to $x_{b} = (15 - 30) μ m$ .
Dpp. The characteristic decay length of Dpp was set to $λ_{D p p} = 20 μ m$ (Kicheva et al., 2007). The overall size of the domain patterned by Dpp was set to $L = 80 μ m$ (Figure 2H of Bakker et al., 2020). The boundary of Dpp target salm expression was set to $x_{b} = (36 - 54) μ m$ (Figure 2H of Bakker et al., 2020). The radius of the sensor for Dpp and Hh mediated patterning was set to $a = 2 μ m$ , same as for Wg, since they diffuse through the same wing imaginal disk. While we only focus on the boundary defined by salm for simplicity, other genes affected by Dpp, including omb, brk, and dad, also display changes in their spatial expression profiles at locations similar to that of salm (Bakker et al., 2020).

Length scales of Cyclops, Squint, and Fgf8

In zebrafish embryos, the morphogens Cyclops, Squint, and Fgf8 are responsible for the induction of the endoderm and the mesoderm. All three morphogens are synthesized at the embryonic margin, and diffuse toward the animal pole (Appendix 2—figure 3a). For simplicity, we will only consider zebrafish embryos around the ∼50% epiboly stage, when the embryo encompasses the top half of the spherical yolk. The length scales associated with Cyclops, Squint, and Ffg8 are also given in Appendix 3—table 1.

Cyclops and Squint. The diffusion coefficient ( $D$ ) and depletion rate ( $k_{d}$ ) of Cyclops and Squint have been measured through the ectopic expressions of the fluorescently tagged constructs. The resulting characteristic decay lengths ( $λ \equiv \sqrt{D / k_{d}}$ ) are $λ_{c y c l o p s} = 76 μ m$ and $λ_{s q u i n t} = 178 μ m$ (Müller et al., 2012). The target boundary of Nodal signaling was set by the expression domain width of the target genes gsc (50 µm) and ntl (125 µm) (Figure 1f of Dubrulle et al., 2015). The width of the domain of Nodal morphogen production was assumed to be 25 µm, following the kinetic model describing Nodal dynamics in Dubrulle et al., 2015. The radius of the zebrafish embryo was set to $R = 200 μ m$ (Figure 1 of Nowak et al., 2011). Overall, for Nodal morphogens Cyclops and Squint, $z_{b}$ ranges in [(50 µm – 25 µm)/200 µm, (125 µm – 25 µm)/200 µm] = [0.125, 0.5].
Fgf8. The characteristic decay length of Fgf8, $λ_{f g f 8} = 197 μ m$ , was measured from the ectopic expression of the fluorescently tagged constructs (Yu et al., 2009). Fgf8 affects the expression of target genes such as spry4, spry2, and pea3. The target boundary of Fgf8 was set by the expression domain width of the target genes spry4 (70 µm), spry2 (85 µm), and pea3 (120 µm) (Figure 1 of Nowak et al., 2011). Subtracting the width of the Fgf8 expression domain (70 µm) and dividing by the radius of the embryo (200 µm) amount to $z_{b}$ in the range of $[(70 μ m - 70 μ m) / 200 μ m, (120 μ m - 70 μ m) / 200 μ m] = [0.0, 0.25]$ .

Appendix 3—table 1

Various length scales associated with the naturally occurring morphogen profiles found in the fruit fly (*) and zebrafish (†) embryos.

For each plot of the morphogen profile in Figure 2—figure supplement 1, the solid red line denotes the exponentially decaying profile with the associated characteristic length $λ$ . Further details on how each value was obtained from the associated reference are provided in the text of Appendix 3.

	$L^{*}$ or $R^{†}$ (µm)	$λ [\equiv \sqrt{D / k_{d}}]$ (µm)	Range of $x_{b}^{*}$ or $z_{b}^{†}$ (µm)	$a$ (µm)
Bcd*	425 Houchmandzadeh et al., 2002	100 Houchmandzadeh et al., 2002	140–188 Perry et al., 2012	2.6 Gregor et al., 2007
Wg*	70 Chaudhary et al., 2019	6 Kicheva et al., 2007	11–21 Bakker et al., 2020	2 Chaudhary et al., 2019
Hh*	100 Torroja et al., 2004	8 Wartlick et al., 2011	15–30 Torroja et al., 2004	2 Chaudhary et al., 2019
Dpp*	80 Bakker et al., 2020	20 Kicheva et al., 2007	36–54 Bakker et al., 2020	2 Chaudhary et al., 2019
Cyclops†	200 Nowak et al., 2011	76 Dubrulle et al., 2015	25–100 Dubrulle et al., 2015	10 Yu et al., 2009
Squint†	200 Nowak et al., 2011	178 Dubrulle et al., 2015	25–100 Dubrulle et al., 2015	10 Yu et al., 2009
Fgf8†	200 Nowak et al., 2011	197 Yu et al., 2009	0–50 Nowak et al., 2011	10 Yu et al., 2009

Time scales of Bcd, Wg, Hh, and Dpp

The expression for the variance of the space-time-averaged signal, $σ_{m}^{2} (x)$ (Equation 11), is derived with the assumption of a long measurement time, $T ≫ k_{d}^{- 1}$ (Appendix 1 of Fancher and Mugler, 2020). Additionally, in order for $T j_{in}$ to represent the total cost of establishing and maintaining the morphogen profile, we require the measurement time to be much longer than the time scale of reaching the steady state (i.e. $T ≫ τ (x_{b})$ ). For Bcd, Wg, Hh, and Dpp, $L ≫ λ$ and $x_{b} / λ \approx 2$ , leading to the approximate expression for the characteristic time of $τ (x_{b}) \approx 1.5 k_{d}^{- 1}$ (Equation 4). We assume that the maximum measurement time, $T_{\max}$ , is set by the cell doubling time, after which the fates of the daughter cells must be re-established.

Bcd: $T_{m a x} [\approx 10 m i n$ (Foe and Alberts, 1983)] $< k_{d}^{- 1} [\approx 1 h r$ (Drocco et al., 2011)].
Wg: $T_{m a x} [\approx 5 h r$ (Martín et al., 2009)] $≫ k_{d}^{- 1} [\approx 15 m i n$ (Kicheva et al., 2007)].
Hh: $T_{m a x} [\approx 5 h r$ (Martín et al., 2009)] $≫ k_{d}^{- 1} [\approx 5 m i n$ (Nahmad and Stathopoulos, 2009)].
Dpp: $T_{m a x} [\approx 5 h r$ (Martín et al., 2009)] $> k_{d}^{- 1} [\approx 1.4 h r$ (Kicheva et al., 2007)].

For Dpp, Hh, and Wg, the space-time-averaged measurements may be performed for sufficiently long durations to represent the cost-precision trade-off by $π_{T}$ (Equation 14). For Bcd, $π_{o}$ , which is formulated for the point measurement, may better represent the pertinent cost-precision trade-off.

Data availability

All data analyzed in this study are from the figures of previously published works, shown in Figure 2-S1 and Appendix 2 Figure 3. The reference associated with each panel is provided in the respective figure legend.

References

(2014) Fly wing vein patterns have spatial reproducibility of a single cell
Journal of the Royal Society Interface 11:20140443.

https://doi.org/10.1098/rsif.2014.0443
- PubMed
- Google Scholar
(2018) Scale-invariant patterning by size-dependent inhibition of nodal signalling
Nature Cell Biology 20:1032–1042.

https://doi.org/10.1038/s41556-018-0155-7
- PubMed
- Google Scholar
(2020) The WG and DPP morphogens regulate gene expression by modulating the frequency of transcriptional bursts
eLife 9:e56076.

https://doi.org/10.7554/eLife.56076
- PubMed
- Google Scholar
1. Barato AC
2. Seifert U
(2015) Thermodynamic uncertainty relation for biomolecular processes
Physical Review Letters 114:1158101.

https://doi.org/10.1103/PhysRevLett.114.158101
- Google Scholar
(2010) How long does it take to establish a morphogen gradient?
Biophysical Journal 99:L59–L61.

https://doi.org/10.1016/j.bpj.2010.07.045
- PubMed
- Google Scholar
(2011) Formation of morphogen gradients: local accumulation time
Physical Review E 83:051906.

https://doi.org/10.1103/PhysRevE.83.051906
- Google Scholar
1. Berg HC
2. Purcell EM
(1977) Physics of chemoreception
Biophysical Journal 20:193–219.

https://doi.org/10.1016/S0006-3495(77)85544-6
- PubMed
- Google Scholar
(2008) Precision of the Dpp gradient
Development 135:1137–1146.

https://doi.org/10.1242/dev.012062
- PubMed
- Google Scholar
1. Briscoe J
2. Small S
(2015) Morphogen rules: design principles of gradient-mediated embryo patterning
Development 142:3996–4009.

https://doi.org/10.1242/dev.129452
- PubMed
- Google Scholar
1. Chaudhary V
2. Hingole S
3. Frei J
4. Port F
5. Strutt D
6. Boutros M
(2019) Robust wnt signaling is maintained by a wg protein gradient and Fz2 receptor activity in the developing Drosophila wing
Development 146:dev174789.

https://doi.org/10.1242/dev.174789
- PubMed
- Google Scholar
1. Chen W
2. Huang H
3. Hatori R
4. Kornberg TB
(2017) Essential basal cytonemes take up hedgehog in the Drosophila wing imaginal disc
Development 144:3134–3144.

https://doi.org/10.1242/dev.149856
- PubMed
- Google Scholar
1. Crick F
(1970) Diffusion in embryogenesis
Nature 225:420–422.

https://doi.org/10.1038/225420a0
- Google Scholar
1. Dechant A
2. Sasa SI
(2020) Fluctuation-response inequality out of equilibrium
PNAS 117:6430–6436.

https://doi.org/10.1073/pnas.1918386117
- PubMed
- Google Scholar
(2020) A mechanism for hunchback promoters to readout morphogenetic positional information in less than a minute
eLife 9:e49758.

https://doi.org/10.7554/eLife.49758
- PubMed
- Google Scholar
1. Driever W
2. Nüsslein-Volhard C
(1988a) A gradient of bicoid protein in Drosophila embryos
Cell 54:83–93.

https://doi.org/10.1016/0092-8674(88)90182-1
- PubMed
- Google Scholar
1. Driever W
2. Nüsslein-Volhard C
(1988b) The bicoid protein determines position in the Drosophila embryo in a concentration-dependent manner
Cell 54:95–104.

https://doi.org/10.1016/0092-8674(88)90183-3
- PubMed
- Google Scholar
1. Drocco JA
2. Grimm O
3. Tank DW
4. Wieschaus E
(2011) Measurement and perturbation of morphogen lifetime: effects on gradient shape
Biophysical Journal 101:1807–1815.

https://doi.org/10.1016/j.bpj.2011.07.025
- PubMed
- Google Scholar
(2012) The synthesis-diffusion-degradation model explains bicoid gradient formation in unfertilized eggs
Physical Biology 9:055004.

https://doi.org/10.1088/1478-3975/9/5/055004
- PubMed
- Google Scholar
(2015) Response to nodal morphogen gradient is determined by the kinetics of target gene induction
eLife 4:e05042.

https://doi.org/10.7554/eLife.05042
- Google Scholar
1. Emberly E
(2008) Optimizing the readout of morphogen gradients
Physical Review E 77:041903.

https://doi.org/10.1103/PhysRevE.77.041903
- Google Scholar
(2000) Gradient formation of the TGF-beta homolog Dpp
Cell 103:981–992.

https://doi.org/10.1016/S0092-8674(00)00200-2
- PubMed
- Google Scholar
(2009) Role of spatial averaging in the precision of gene expression patterns
Physical Review Letters 103:258101.

https://doi.org/10.1103/PhysRevLett.103.258101
- PubMed
- Google Scholar
(2018) Information thermodynamics of turing patterns
Physical Review Letters 121:108301.

https://doi.org/10.1103/PhysRevLett.121.108301
- PubMed
- Google Scholar
1. Fancher S
2. Mugler A
(2020) Diffusion vs. direct transport in the precision of morphogen readout
eLife 9:e58981.

https://doi.org/10.7554/eLife.58981
- PubMed
- Google Scholar
1. Foe VE
2. Alberts BM
(1983) Studies of nuclear and cytoplasmic behaviour during the five mitotic cycles that precede gastrulation in Drosophila embryogenesis
Journal of Cell Science 61:31–70.

https://doi.org/10.1242/jcs.61.1.31
- PubMed
- Google Scholar
(2016) Dissipation bounds all Steady-State current fluctuations
Physical Review Letters 116:120601.

https://doi.org/10.1103/PhysRevLett.116.120601
- PubMed
- Google Scholar
(2007) Probing the limits to positional information
Cell 130:153–164.

https://doi.org/10.1016/j.cell.2007.05.025
- PubMed
- Google Scholar
(2017) Concentration dependent chromatin states induced by the bicoid morphogen gradient
eLife 6:e28275.

https://doi.org/10.7554/eLife.28275
- PubMed
- Google Scholar
1. He F
2. Saunders TE
3. Wen Y
4. Cheung D
5. Jiao R
6. ten Wolde PR
7. Howard M
8. Ma J
(2010) Shaping a morphogen gradient for positional precision
Biophysical Journal 99:697–707.

https://doi.org/10.1016/j.bpj.2010.04.073
- PubMed
- Google Scholar
1. Heuett WJ
2. Qian H
(2006) Grand canonical markov model: a stochastic theory for open nonequilibrium biochemical networks
The Journal of Chemical Physics 124:044110.

https://doi.org/10.1063/1.2165193
- PubMed
- Google Scholar
1. Horowitz JM
2. Gingrich TR
(2017) Proof of the finite-time thermodynamic uncertainty relation for steady-state currents
Physical Review E 96:020103.

https://doi.org/10.1103/PhysRevE.96.020103
- PubMed
- Google Scholar
1. Horowitz JM
2. Gingrich TR
(2020) Thermodynamic uncertainty relations constrain non-equilibrium fluctuations
Nature Physics 16:15–20.

https://doi.org/10.1038/s41567-019-0702-6
- Google Scholar
(2002) Establishment of developmental precision and proportions in the early Drosophila embryo
Nature 415:798–802.

https://doi.org/10.1038/415798a
- PubMed
- Google Scholar
1. Huang H
2. Kornberg TB
(2015) Myoblast cytonemes mediate wg signaling from the wing imaginal disc and Delta-Notch signaling to the air sac primordium
eLife 4:e06114.

https://doi.org/10.7554/eLife.06114
- PubMed
- Google Scholar
1. Hwang W
2. Hyeon C
(2018) Energetic costs, precision, and transport efficiency of molecular motors
The Journal of Physical Chemistry Letters 9:513–520.

https://doi.org/10.1021/acs.jpclett.7b03197
- PubMed
- Google Scholar
1. Hyeon C
2. Hwang W
(2017) Physical insight into the thermodynamic uncertainty relation using brownian motion in tilted periodic potentials
Physical Review E 96:012156.

https://doi.org/10.1103/PhysRevE.96.012156
- PubMed
- Google Scholar
1. Ilker E
2. Hinczewski M
(2019) Modeling the growth of organisms validates a general relation between metabolic costs and natural selection
Physical Review Letters 122:2381011.

https://doi.org/10.1103/PhysRevLett.122.238101
- Google Scholar
(2006) Senseless and daughterless confer neuronal identity to epithelial cells in the Drosophila wing margin
Development 133:1683–1692.

https://doi.org/10.1242/dev.02338
- PubMed
- Google Scholar
(2009) Dynamics of the dorsal morphogen gradient
PNAS 106:21707–21712.

https://doi.org/10.1073/pnas.0912395106
- PubMed
- Google Scholar
(2007) Kinetics of morphogen gradient formation
Science 315:521–525.

https://doi.org/10.1126/science.1135774
- PubMed
- Google Scholar
1. Kolomeisky AB
(2011) Formation of a morphogen gradient: acceleration by degradation
The Journal of Physical Chemistry Letters 2:1502–1505.

https://doi.org/10.1021/jz2004914
- Google Scholar
1. Li P
2. Markson JS
3. Wang S
4. Chen S
5. Vachharajani V
6. Elowitz MB
(2018) Morphogen gradient reconstitution reveals hedgehog pathway design principles
Science 360:543–548.

https://doi.org/10.1126/science.aao0645
- PubMed
- Google Scholar
(2013) Dynamic interpretation of maternal inputs by the Drosophila segmentation gene network
PNAS 110:6724–6729.

https://doi.org/10.1073/pnas.1220912110
- PubMed
- Google Scholar
1. Lo WC
2. Zhou S
3. Wan FY
4. Lander AD
5. Nie Q
(2015) Robust and precise morphogen-mediated patterning: trade-offs, constraints and mechanisms
Journal of the Royal Society, Interface 12:20141041.

https://doi.org/10.1098/rsif.2014.1041
- PubMed
- Google Scholar
1. Lucas T
2. Tran H
3. Perez Romero CA
4. Guillou A
5. Fradin C
6. Coppey M
7. Walczak AM
8. Dostatni N
(2018) 3 minutes to precisely measure morphogen concentration
PLOS Genetics 14:e1007676.

https://doi.org/10.1371/journal.pgen.1007676
- PubMed
- Google Scholar
1. Lynch M
2. Marinov GK
(2015) The bioenergetic costs of a gene
PNAS 112:15690–15695.

https://doi.org/10.1073/pnas.1514974112
- PubMed
- Google Scholar
(2009) Cell competition, growth and size control in the Drosophila wing imaginal disc
Development 136:3747–3756.

https://doi.org/10.1242/dev.038406
- PubMed
- Google Scholar
1. Morishita Y
2. Iwasa Y
(2009) Accuracy of positional information provided by multiple morphogen gradients with correlated noise
Physical Review E 79:061905.

https://doi.org/10.1103/PhysRevE.79.061905
- Google Scholar
1. Müller P
2. Rogers KW
3. Jordan BM
4. Lee JS
5. Robson D
6. Ramanathan S
7. Schier AF
(2012) Differential diffusivity of nodal and lefty underlies a reaction-diffusion patterning system
Science 336:721–724.

https://doi.org/10.1126/science.1221920
- PubMed
- Google Scholar
1. Müller P
2. Rogers KW
3. Yu SR
4. Brand M
5. Schier AF
(2013) Morphogen transport
Development 140:1621–1638.

https://doi.org/10.1242/dev.083519
- PubMed
- Google Scholar
1. Nahmad M
2. Stathopoulos A
(2009) Dynamic interpretation of hedgehog signaling in the Drosophila wing disc
PLOS Biology 7:e1000202.

https://doi.org/10.1371/journal.pbio.1000202
- PubMed
- Google Scholar
1. Nellen D
2. Burke R
3. Struhl G
4. Basler K
(1996) Direct and long-range action of a DPP morphogen gradient
Cell 85:357–368.

https://doi.org/10.1016/S0092-8674(00)81114-9
- PubMed
- Google Scholar
1. Neumann CJ
2. Cohen SM
(1997) Long-range action of wingless organizes the dorsal-ventral Axis of the Drosophila wing
Development 124:871–880.

https://doi.org/10.1242/dev.124.4.871
- PubMed
- Google Scholar
1. Nowak M
2. Machate A
3. Yu SR
4. Gupta M
5. Brand M
(2011) Interpretation of the FGF8 morphogen gradient is regulated by endocytic trafficking
Nature Cell Biology 13:153–158.

https://doi.org/10.1038/ncb2155
- PubMed
- Google Scholar
1. Oginuma M
2. Moncuquet P
3. Xiong F
4. Karoly E
5. Chal J
6. Guevorkian K
7. Pourquié O
(2017) A gradient of glycolytic activity coordinates FGF and wnt signaling during elongation of the body Axis in amniote embryos
Developmental Cell 40:342–353.

https://doi.org/10.1016/j.devcel.2017.02.001
- PubMed
- Google Scholar
1. Oginuma M
2. Harima Y
3. Tarazona OA
4. Diaz-Cuadros M
5. Michaut A
6. Ishitani T
7. Xiong F
8. Pourquié O
(2020) Intracellular pH controls WNT downstream of glycolysis in amniote embryos
Nature 584:98–101.

https://doi.org/10.1038/s41586-020-2428-0
- PubMed
- Google Scholar
1. Perry MW
2. Bothma JP
3. Luu RD
4. Levine M
(2012) Precision of hunchback expression in the Drosophila embryo
Current Biology 22:2247–2252.

https://doi.org/10.1016/j.cub.2012.09.051
- PubMed
- Google Scholar
(2019) Optimal decoding of cellular identities in a genetic network
Cell 176:844–855.

https://doi.org/10.1016/j.cell.2019.01.007
- PubMed
- Google Scholar
(2017) Finite-time generalization of the thermodynamic uncertainty relation
Physical Review E 96:012101.

https://doi.org/10.1103/PhysRevE.96.012101
- PubMed
- Google Scholar
(2019) Heat oscillations driven by the embryonic cell cycle reveal the energetic costs of signaling
Developmental Cell 48:646–658.

https://doi.org/10.1016/j.devcel.2018.12.024
- PubMed
- Google Scholar
(2020) Contribution of increasing plasma membrane to the energetic cost of early zebrafish embryogenesis
Molecular Biology of the Cell 31:520–526.

https://doi.org/10.1091/mbc.E19-09-0529
- PubMed
- Google Scholar
(2020) Optogenetic investigation of BMP target gene expression diversity
eLife 9:e58641.

https://doi.org/10.7554/eLife.58641
- PubMed
- Google Scholar
1. Rosenbauer J
2. Zhang C
3. Mattes B
4. Reinartz I
5. Wedgwood K
6. Schindler S
7. Sinner C
8. Scholpp S
9. Schug A
(2020) Modeling of Wnt-mediated tissue patterning in vertebrate embryogenesis
PLOS Computational Biology 16:e1007417.

https://doi.org/10.1371/journal.pcbi.1007417
- PubMed
- Google Scholar
1. Saunders TE
2. Howard M
(2009) Morphogen profiles can be optimized to buffer against noise
Physical Review E 80:041902.

https://doi.org/10.1103/PhysRevE.80.041902
- PubMed
- Google Scholar
1. Shvartsman SY
2. Baker RE
(2012) Mathematical models of morphogen gradients and their effects on gene expression
Wiley Interdisciplinary Reviews: Developmental Biology 1:715–730.

https://doi.org/10.1002/wdev.55
- PubMed
- Google Scholar
1. Song Y
2. Park JO
3. Tanner L
4. Nagano Y
5. Rabinowitz JD
6. Shvartsman SY
(2019) Energy budget of Drosophila embryogenesis
Current Biology 29:R566–R567.

https://doi.org/10.1016/j.cub.2019.05.025
- PubMed
- Google Scholar
1. Song Y
2. Hyeon C
(2020) Thermodynamic cost, speed, fluctuations, and error reduction of biological copy machines
The Journal of Physical Chemistry Letters 11:3136–3143.

https://doi.org/10.1021/acs.jpclett.0c00545
- PubMed
- Google Scholar
1. Song Y
2. Hyeon C
(2021) Thermodynamic uncertainty relation to assess biological processes
The Journal of Chemical Physics 154:130901.

https://doi.org/10.1063/5.0043671
- PubMed
- Google Scholar
1. Strigini M
2. Cohen SM
(1997) A hedgehog activity gradient contributes to AP axial patterning of the Drosophila wing
Development 124:4697–4705.

https://doi.org/10.1242/dev.124.22.4697
- PubMed
- Google Scholar
(1989) The gradient morphogen bicoid is a concentration-dependent transcriptional activator
Cell 57:1259–1273.

https://doi.org/10.1016/0092-8674(89)90062-7
- PubMed
- Google Scholar
1. Tabata T
2. Kornberg TB
(1994) Hedgehog is a signaling protein with a key role in patterning Drosophila imaginal discs
Cell 76:89–102.

https://doi.org/10.1016/0092-8674(94)90175-9
- PubMed
- Google Scholar
1. Teimouri H
2. Kolomeisky AB
(2014) Development of morphogen gradient: the role of dimension and discreteness
The Journal of Chemical Physics 140:085102.

https://doi.org/10.1063/1.4866453
- Google Scholar
1. Teimouri H
2. Kolomeisky AB
(2016) New model for understanding mechanisms of biological signaling: direct transport via cytonemes
The Journal of Physical Chemistry Letters 7:180–185.

https://doi.org/10.1021/acs.jpclett.5b02703
- PubMed
- Google Scholar
(2015) Positional information, positional error, and readout precision in morphogenesis: a mathematical framework
Genetics 199:39–59.

https://doi.org/10.1534/genetics.114.171850
- PubMed
- Google Scholar
1. Tkačik G
2. Gregor T
(2021) The many bits of positional information
Development 148:dev176065.

https://doi.org/10.1242/dev.176065
- PubMed
- Google Scholar
(2004) Patched controls the hedgehog gradient by endocytosis in a dynamin-dependent manner, but this internalization does not play a major role in signal transduction
Development 131:2395–2408.

https://doi.org/10.1242/dev.01102
- PubMed
- Google Scholar
(2007) Fundamental limits to position determination by concentration gradients
PLOS Computational Biology 3:e78.

https://doi.org/10.1371/journal.pcbi.0030078
- PubMed
- Google Scholar
1. Turing AM
(1990) The chemical basis of morphogenesis. 1953
Bulletin of Mathematical Biology 52:153–197.

https://doi.org/10.1016/S0092-8240(05)80008-4
- PubMed
- Google Scholar
(2011) Dynamics of Dpp signaling and proliferation control
Science 331:1154–1159.

https://doi.org/10.1126/science.1200037
- PubMed
- Google Scholar
1. Wolpert L
(1969) Positional information and the spatial pattern of cellular differentiation
Journal of Theoretical Biology 25:1–47.

https://doi.org/10.1016/S0022-5193(69)80016-0
- Google Scholar
1. Yu SR
2. Burkhardt M
3. Nowak M
4. Ries J
5. Petrásek Z
6. Scholpp S
7. Schwille P
8. Brand M
(2009) Fgf8 morphogen gradient forms by a source-sink mechanism with freely diffusing molecules
Nature 461:533–536.

https://doi.org/10.1038/nature08391
- PubMed
- Google Scholar
1. Zagorski M
2. Tabata Y
3. Brandenberg N
4. Lutolf MP
5. Tkačik G
6. Bollenbach T
7. Briscoe J
8. Kicheva A
(2017) Decoding of position in the developing neural tube from antiparallel morphogen gradients
Science 356:1379–1383.

https://doi.org/10.1126/science.aam5887
- PubMed
- Google Scholar
(1996) Direct and long-range action of a wingless morphogen gradient
Cell 87:833–844.

https://doi.org/10.1016/S0092-8674(00)81991-1
- PubMed
- Google Scholar
1. Zhou S
2. Lo WC
3. Suhalim JL
4. Digman MA
5. Gratton E
6. Nie Q
7. Lander AD
(2012) Free extracellular diffusion creates the Dpp morphogen gradient of the Drosophila wing disc
Current Biology 22:668–675.

https://doi.org/10.1016/j.cub.2012.02.065
- PubMed
- Google Scholar

Article and author information

Author details

Yonghyun Song

Korea Institute for Advanced Study, Seoul, Republic of Korea

Contribution
Conceptualization, Formal analysis, Validation, Investigation, Visualization, Writing - original draft, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-8713-1294
Changbong Hyeon

Korea Institute for Advanced Study, Seoul, Republic of Korea

Contribution
Conceptualization, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Writing - original draft, Project administration, Writing - review and editing

For correspondence
hyeoncb@kias.re.kr

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4844-7237

Funding

Korea Institute for Advanced Study (CG067102)

Yonghyun Song

Korea Institute for Advanced Study (CG035003)

Changbong Hyeon

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported by the KIAS individual Grants CG067102 (Y.S.) and CG035003 (C.H.) at Korea Institute for Advanced Study. We thank the Center for Advanced Computation in KIAS for providing computing resources.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.