Costprecision tradeoff relation determines the optimal morphogen gradient for accurate biological pattern formation
Abstract
Spatial boundaries formed during animal development originate from the prepatterning 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 diffusiondepletion 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 costeffectiveness 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 prepatterning of tissues with biomolecules, called morphogens, that instruct cells to acquire distinct cell fates in a concentrationdependent 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 cellfate 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 tradeoff between the cost and precision of the morphogen profile formation in the framework of the reactiondiffusion 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 anteriorposterior axis of the Drosophila embryo through the Bicoid (Bcd) gradient. In the 2hr postfertilization 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 anteriorend of the embryo by maternally deposited bcd mRNAs. Its subsequent diffusion and degradation engender an exponentially decaying profile of Bcd concentration (Driever and NüssleinVolhard, 1988a; Driever and NüssleinVolhard, 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 dorsalventral (DV) and anteriorposterior (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 (JafarNejad 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 spaltmajor (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 singlecell width emerge from the coordinated actions of multiple patterning events (Abouchar et al., 2014).
The positional information, encoded in the concentration profile $\rho (x)$ (red line in Figure 1), is decoded to yield the target gene expression profile, $g(\rho (x))$ (blue line in Figure 1). The morphogen profile, $\rho (x)$, and the corresponding gene expression level, $g(\rho (x))$, together specify the cell fate in a morphogen concentrationdependent 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.
We begin by defining a target boundary, $x={x}_{\mathrm{b}}$, where the morphogen concentration is at its critical threshold value ${\rho}_{\mathrm{b}}\equiv \rho ({x}_{\mathrm{b}})$. The nuclei exposed to morphogen concentrations higher than ${\rho}_{\mathrm{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 ${\widehat{\rho}}_{\mathrm{b}}$, the nucleus can estimate the location of target boundary, such that ${\widehat{x}}_{\mathrm{b}}={\widehat{x}}_{\mathrm{b}}({\widehat{\rho}}_{\mathrm{b}})$. For instance, the Bcd concentration at each nucleus at position $x$ ($\rho (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 positiondependent expression profile of hb ($g(\rho (x))$) (Gregor et al., 2007). The error (variance) associated with the measured concentration at $x={x}_{\mathrm{b}}$ with respect to its true value ${\rho}_{\mathrm{b}}$ is given by ${\sigma}_{\rho}^{2}({x}_{\mathrm{b}})[\equiv \u27e8{({\widehat{\rho}}_{\mathrm{b}}{\rho}_{\mathrm{b}})}^{2}\u27e9=(1/N){\sum}_{i=1}^{N}{({\widehat{\rho}}_{\mathrm{b},i}{\rho}_{\mathrm{b}})}^{2}]$. In other words, ${\sigma}_{\rho}^{2}({x}_{\mathrm{b}})$ represents the inherent variability of the morphogen profile, which can be determined experimentally by performing multiple measurements, ${\widehat{\rho}}_{\mathrm{b},i}$ ($i=1,2,\mathrm{\cdots}N$). For the nucleus to make a single measurement, ${\widehat{\rho}}_{\mathrm{b},i}$, it would integrate the morphogeninduced 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, $\rho ({\widehat{x}}_{\mathrm{b}})\approx \rho ({x}_{\mathrm{b}})+{\left(\partial \rho ({\widehat{x}}_{\mathrm{b}})/\partial {\widehat{x}}_{\mathrm{b}}\right)}_{{\widehat{x}}_{\mathrm{b}}={x}_{\mathrm{b}}}({\widehat{x}}_{\mathrm{b}}{x}_{\mathrm{b}})+\mathrm{\cdots}$, allows one to relate the error in measured concentration with the positional error in estimating the target boundary ${\sigma}_{x}^{2}({x}_{\mathrm{b}})[\equiv \u27e8{({\widehat{x}}_{\mathrm{b}}{x}_{\mathrm{b}})}^{2}\u27e9=(1/N){\sum}_{i=1}^{N}{({\widehat{x}}_{\mathrm{b},i}{x}_{\mathrm{b}})}^{2}]$ as follows:
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 ‘costeffectiveness’ 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 costminimizing λ (Emberly, 2008). Here, we expand this argument with an indepth 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 morphogensensing 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) Spacetimeaveraged 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 spacetimeaveraged 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 productionindependent quantity by taking the product of total cost and the positional error to quantify the tradeoff between the cost and precision of the morphogen profile, and show that the tradeoff product can be minimized when the morphogen gradient’s characteristic length, λ, is properly selected. We evaluate the costeffectiveness 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 onedimensional (1D) array of cells in the domain $0\le x\le 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}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}$) and length ($l}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}$) at the interval between $x$ and $x+{l}_{\mathrm{cell}}$ by $\rho (x,t)[\mathrm{conc}\equiv \mathrm{\#}/{v}_{\mathrm{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}$[\mathrm{conc}\times {l}_{\mathrm{cell}}/\mathrm{time}]$. At all positions, the morphogen molecules are depleted with rate ${k}_{\mathrm{d}}$ [${\mathrm{time}}^{1}$], while spreading across the cells with the diffusivity $D[{l}_{\mathrm{cell}}^{2}/\mathrm{time}]$ (Figure 2a). For $L/{l}_{\mathrm{cell}}\gg 1$ (i.e. at the continuum limit), the spatiotemporal dynamics of the morphogen is described using the reactiondiffusion equation
with boundary conditions ${D{\partial}_{x}\rho (x,t)}_{x=0}={j}_{\mathrm{in}}$ and ${D{\partial}_{x}\rho (x,t)}_{x=L}=0$. Then, the concentration profile at steady state is obtained as
where $\lambda =(\sqrt{D/{k}_{\mathrm{d}}})$ is the characteristic decay length determined by the diffusivity of morphogen ($D$) and the depletion rate (${k}_{\mathrm{d}}$). For large system size ($L\gg \lambda $), ${\rho}_{\mathrm{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 reactiondiffusion model.
(i) With a function $R(x,t)\equiv ({\rho}_{\mathrm{ss}}(x)\rho (x,t))/({\rho}_{\mathrm{ss}}(x)\rho (x,0))$, which captures the evolution of the morphogen profile from $\rho (x,t=0)=0$ to the steady state value ${\rho}_{\mathrm{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,
where $f(x;\lambda )\equiv \frac{1}{2}\left[1+\frac{L}{\lambda}\mathrm{coth}\left(\frac{L}{\lambda}\right)\frac{Lx}{\lambda}\mathrm{tanh}\left(\frac{Lx}{\lambda}\right)\right]$ is a dimensionless quantity determined by $x$ and λ. Then, ${v}_{\mathrm{cell}}{j}_{\mathrm{in}}/{l}_{\mathrm{cell}}\times \tau (x)$ quantifies the total number of morphogens produced over the time scale of $\tau (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
where $C(x)$ increases linearly with $x$ for large system size ($L\gg \lambda $) (Berezhkovskii et al., 2010).
The proportionality constant, ${\alpha}_{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 ${\alpha}_{o}\approx 494\times 4\approx 2\times {10}^{3}$ ATPs, where we assume that each peptidebond formation requires 4 ATPs (Lynch and Marinov, 2015). In the first ∼2 hr of Drosophila embryogenesis when the anteriorposterior axis patterning takes place, ∼5 × 10^{8} Bcd molecules are produced (Drocco et al., 2012). Thus, the cost of generating the Bcd profile is on the order of 10^{12} ATPs, which is a small fraction of the total energy budget of Drosophila embryogenesis (∼7 × 10^{16} 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):
where $\widehat{\rho}(x)$ represents the morphogen concentration estimated at $x$ from a measurement. For large system size ($L\gg \lambda $), ${\u03f5}^{2}(x)\propto ({\lambda}^{3}/{x}^{2}){e}^{x/\lambda}$. Note that the multiple measurements of $\widehat{\rho}(x)$ yield the mean concentration profile $\u27e8\widehat{\rho}(x)\u27e9=(1/N){\sum}_{i=1}^{N}{\widehat{\rho}}_{i}(x)$, which is equivalent to ${\rho}_{\mathrm{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}_{\mathrm{cell}}^{2}{\sigma}_{\rho}^{2}(x)={v}_{\mathrm{cell}}\u27e8\widehat{\rho}(x)\u27e9={v}_{\mathrm{cell}}{\rho}_{\mathrm{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, ${\u03f5}^{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 ${\u03f5}^{2}({x}_{\mathrm{b}})\approx 7\times {10}^{4}$ for Bcd to ${\u03f5}^{2}({x}_{\mathrm{b}})\approx 1.5\times {10}^{2}$ for Dpp, where ${x}_{\mathrm{b}}\approx (0.40.5)L$ for both systems (Gregor et al., 2007; Bollenbach et al., 2008).
There is a tradeoff between $\mathcal{C}(x)$ and ${\u03f5}^{2}(x)$. If ${n}_{\mathrm{m}}$ is the number of morphogen molecules produced for a certain time duration then $\mathcal{C}(x)\propto {n}_{\mathrm{m}}$ and ${\sigma}_{\rho}^{2}(x)/{\u27e8\widehat{\rho}(x)\u27e9}^{2}\propto 1/{n}_{\mathrm{m}}$, the latter of which simply arises from the central limit theorem, or can be rationalized based on the BergPurcell result ($\delta \widehat{\rho}/\widehat{\rho}\sim 1/\sqrt{Da\widehat{\rho}\tau}$) (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}_{\mathrm{m}}$ corresponds to $\left({v}_{\mathrm{cell}}{j}_{\mathrm{in}}/{l}_{\mathrm{cell}}{k}_{\mathrm{d}}\right)$ (the morphogen molecules produced for the time duration ${k}_{\mathrm{d}}^{1}$). The n_{m}independent tradeoff 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,
where the λ dependence of the tradeoff product is made explicit for the discussion that follows.
We are mainly concerned with the precision at the boundary position, $x={x}_{\mathrm{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}_{\mathrm{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 tradeoff product ${\pi}_{o}({x}_{\mathrm{b}};\lambda )$ is close to its lower bound ${\pi}_{o,\mathrm{min}}({x}_{\mathrm{b}};\lambda )$, the system is costeffective at generating precise morphogen profiles at ${x}_{\mathrm{b}}$. Generally, ${\pi}_{o,\mathrm{min}}({x}_{\mathrm{b}};\lambda )$ increases monotonically with ${x}_{\mathrm{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}_{\mathrm{b}}$, the value of the tradeoff product is determined solely by the decay length, λ. By tuning λ, one can find the value of ${\lambda}_{\mathrm{min}}$ that minimizes the tradeoff product at $x={x}_{\mathrm{b}}$, that is, ${\pi}_{o,\mathrm{min}}({x}_{\mathrm{b}})$, such that
The optimal decay length, ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})$, also increases monotonically with the target boundary ${x}_{\mathrm{b}}$; for ${x}_{\mathrm{b}}<L/2$, ${\lambda}_{\mathrm{min}}$ is well approximated by the linear relationship ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})\approx 0.43{x}_{\mathrm{b}}$ (See Equation 40 for the derivation).
It is of great interest to compare ${x}_{\mathrm{b}}$ and λ of real morphogen profiles against the optimal ${\lambda}_{\mathrm{min}}({x}_{\mathrm{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}_{\mathrm{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 $\lambda}_{\mathrm{B}\mathrm{c}\mathrm{d}}=100\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$, $\lambda}_{\mathrm{W}\mathrm{g}}=6\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$, $\lambda}_{\mathrm{H}\mathrm{h}}=8\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$, and $\lambda}_{\mathrm{D}\mathrm{p}\mathrm{p}}=20\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$ (Houchmandzadeh et al., 2002; Kicheva et al., 2007; Wartlick et al., 2011). Figure 2b shows the pairs of ${x}_{\mathrm{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, ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})$. Taken together, our findings indicate that the concentration profiles of the four morphogens are effectively formed under the condition that the costprecision tradeoff, ${\pi}_{o}({x}_{\mathrm{b}};\lambda )$, is minimized to its lower bound, ${\pi}_{o,\mathrm{min}}({x}_{\mathrm{b}};\lambda )$ (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}_{\mathrm{b}}\approx 0.4L$, ${\pi}_{o}({x}_{\mathrm{b}};{\lambda}_{\mathrm{Bcd}})\approx {\alpha}_{o}(L/{l}_{\mathrm{cell}})0.6$ (Figure 2b) where $L/{l}_{\mathrm{cell}}\approx 50$, and the experimentally reported precision is ${\u03f5}^{2}({x}_{\mathrm{b}})\approx 7\times {10}^{4}$ (Gregor et al., 2007). Thus, the cost associated with the Bcd profile is $C({x}_{\mathrm{b}};{\lambda}_{\mathrm{Bcd}})={\pi}_{o}({x}_{\mathrm{b}};{\lambda}_{\mathrm{Bcd}})/{\u03f5}^{2}({x}_{\mathrm{b}};{\lambda}_{\mathrm{Bcd}})\approx 4\times {10}^{4}{\alpha}_{o}$, or equivalently, the thermodynamic cost of synthesizing and degrading 4 × 10^{4} Bcd molecules. This is orders of magnitude smaller than the earlier estimate of ∼5 × 10^{8} 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 ${\alpha}_{o}$ so that $C({x}_{\mathrm{b}})$ represents the overall cost to produce the concentration profile of nuclear Bcd.
Spacetimeaveraged 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 ${\u03f5}^{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 spacetimeaveraging.
To incorporate the scenario of spacetimeaveraging, we assume that the cell at position $x$ uses a sensor with size $a$ to detect the morphogen concentration over the space interval $(xa,x+a)$ for time $T$ (Figure 3a). Then the spacetimeaveraged molecular count is written as
where ${\widehat{\rho}}_{i}(y,t)$ represents an estimate of molecular count at position $y$ from an $i$th independent measurement. Repeated measurements ($N\gg 1$) lead to the mean value of molecular counts
In this case, $j}_{\mathrm{i}\mathrm{n}$ and ${k}_{\mathrm{d}}$ are in the units of $\mathrm{\#}/\mathrm{time}$ and ${\mathrm{time}}^{1}$, respectively. With a sufficiently long measurement time ($T\gg {k}_{\mathrm{d}}^{1}$), the variance of $m(x)$ can be approximated to Fancher and Mugler, 2020
Analogously to Equation 6, the squared relative error (precision) at position $x$ can be quantified as follows,
The cost of maintaining the morphogen profile at steady states is proportional to ${j}_{\mathrm{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
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}_{\mathrm{in}}T$ cancels off, which yields the following expression of the costprecision tradeoff:
With a given target boundary $x={x}_{\mathrm{b}}$ and the sensor size $a$, the value of ${\pi}_{T}({x}_{\mathrm{b}};\lambda )$ is solely determined by λ, analogously to ${\pi}_{o}({x}_{\mathrm{b}};\lambda )$. Thus, the lower bound of ${\pi}_{T}({x}_{\mathrm{b}};\lambda )$, i.e., ${\pi}_{T,\mathrm{min}}({x}_{\mathrm{b}};\lambda )$, can be determined simply by tuning λ to a value that minimizes ${\pi}_{T}({x}_{\mathrm{b}};\lambda )$, that is, ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})$. Both ${\pi}_{T,\mathrm{min}}({x}_{\mathrm{b}};\lambda )$ and ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})$ increase monotonically with ${x}_{\mathrm{b}}$. In particular, ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})$ is found in the range of $({x}_{\mathrm{b}}a)/2\le {\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})\le {x}_{\mathrm{b}}/2$ (see Equations 42 and 43 for derivation). The boundary of the inaccessible region in Figure 3b constrains ${\pi}_{T}({x}_{\mathrm{b}};\lambda )$, 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% (${\u03f5}_{T}({x}_{\mathrm{b}})\approx 0.1$) at ${x}_{\mathrm{b}}=60a$, the characteristic length must be $\lambda \approx 30a$, which demands that the system synthesize at least ∼182 morphogen molecules (Figure 3b).
As we have done for ${\pi}_{o}$ in the point measurement, we can evaluate ${\pi}_{T}$ of the four morphogen profiles at their target boundaries (${x}_{\mathrm{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\gg {k}_{\mathrm{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 ($\lambda =\sqrt{D/{k}_{\mathrm{d}}}$) of all four morphogen profiles are close to their respective ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})$ values. Thus, from the spacetimeaveraged measurement of morphogen profiles, the four morphogen profiles are also formed under the conditions of nearly optimal costprecision tradeoff (Figure 3b).
Discussion
Comparison of the tradeoff products, ${\pi}_{o}$ and ${\pi}_{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 spacetimeaveraging, defined as $m(x)/(2a)$ at the limit of small sensor size ($a/\lambda \ll 1$) (Equation 10) is identical to the one in the point measurement in the limit of $L\gg \lambda $ (Equation 3), both yielding $({j}_{\mathrm{in}}/\sqrt{D{k}_{\mathrm{d}}}){e}^{x/\lambda}$. On the one hand, ${\u03f5}^{2}(x)$, defined through ${\rho}_{\mathrm{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, ${\u03f5}_{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, ${\u03f5}^{2}(x)$ can be interpreted as the short time limit of ${\u03f5}_{T}^{2}(x)$ measured by a single receptor.
The cost associated with the morphogen profile in the case of point measurement, ${v}_{\mathrm{cell}}{j}_{\mathrm{in}}{l}_{\mathrm{cell}}^{1}\tau (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 spacetimeaveraging, ${j}_{\mathrm{in}}T$, integrated over a time interval $T\gg {k}_{\mathrm{d}}^{1}$, is identical for all positions. The latter assumption is required in the derivation of the expression of ${\u03f5}_{T}^{2}(x)$ (see Appendix 1 of Fancher and Mugler, 2020). We additionally demand $T\gg \tau (x)$ in order for ${j}_{\mathrm{in}}T$ to represent the total cost of morphogen profile formation and maintenance. However, $T$ may not necessarily be larger than either ${k}_{\mathrm{d}}^{1}$ or $\tau (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 costprecision tradeoff involving the pattern formation is presumably at work in between the two scenarios. Thus, of great significance is the finding that ${\pi}_{o}(x;\lambda )$ and ${\pi}_{T}(x;\lambda )$ for the two limiting models can be effectively minimized by similar λ values at a given position: ${\lambda}_{\mathrm{min}}\approx 0.43{x}_{\mathrm{b}}$ for the point measurement, and ${\lambda}_{\mathrm{min}}\approx 0.5{x}_{\mathrm{b}}$ for the spacetimeaveraging.
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 reactiondiffusion model of morphogen dynamics, we extend our discussion of costprecision tradeoff on the models that include unidirectional irreversible steps in the synthesis and depletion by considering a more general reaction model with bidirectional 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}}_{\mathrm{tot}}$) or the thermodynamic cost for the formation of morphogen gradient and clarifies the physical meaning of ${\alpha}_{o}$ (for the relation between $C$ (Equation 5) and ${\dot{S}}_{\mathrm{tot}}$, see Appendix 1 and Equations 32, 33). Similar to the point measurement, we derive expressions for the relaxation time to the steady state (${\tau}_{\mathrm{rev}}(x)$) and the positional error (${\u03f5}_{\mathrm{rev}}^{2}(x)$). The tradeoff among the thermodynamic cost, speed of formation ($\sim {\tau}_{\mathrm{rev}}^{1}(x)$), and precision is quantified by the product of the three quantities,
The lower bound, ${\pi}_{o,\mathrm{rev},\mathrm{min}}({x}_{\mathrm{b}})$, increases monotonically with ${x}_{\mathrm{b}}$ (Appendix 1—figure 1).
With ${\dot{S}}_{\mathrm{tot}}{\tau}_{\mathrm{rev}}({x}_{\mathrm{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 tradeoff relation between the entropy production and the precision of a currentlike 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 modelspecific and positiondependent. The costprecision tradeoff 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 currentlike observable with an oddparity with timereversal; 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 morphogeninduced currents through signaling pathways, such as transcription, as the output observable to which TUR can be directly applied.
Costeffectiveness of precise morphogen profile formation
There are multiple possibilities of achieving the same critical threshold concentration $\rho ({x}_{\mathrm{b}})$ at the target boundary $x={x}_{\mathrm{b}}$ through the combination of morphogen synthesis rate ($j}_{\mathrm{i}\mathrm{n}$), diffusivity ($D$) and depletion rate (${k}_{\mathrm{d}}$) (inset in Figure 4). A steeper morphogen profile with larger $j}_{\mathrm{i}\mathrm{n}$ and smaller $\lambda \phantom{\rule{veryverythickmathspace}{0ex}}(=\sqrt{D/{k}_{\mathrm{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}_{\mathrm{i}\mathrm{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 costprecision tradeoff product, which leads to the formation of costeffective morphogen profiles. In other words, morphogen profiles with optimal characteristic length ${\lambda}_{\mathrm{min}}$ achieve a desired precision when the cost is minimal (Figure 4b). Specifically, ${\lambda}_{\mathrm{min}}\approx (0.430.50){x}_{\mathrm{b}}$, which suggests that the target boundary of the exponentially decaying morphogen profile, $c(x)={c}_{0}{e}^{x/\lambda}$, should be formed at $c({x}_{\mathrm{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}_{\mathrm{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.
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 costprecision tradeoff 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 tradeoff 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 tradeoff bound (Appendix 2—figure 2b). The large tradeoff 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; AlmuedoCastillo 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 costprecision tradeoff 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 costprecision tradeoff relation for the localized synthesisdiffusiondepletion model of the morphogen dynamics in a 1D array of cells. Appendix 2 derives modified tradeoff 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 reactiondiffusion model of morphogen dynamics
Here, we provide detailed derivations of the costprecision tradeoff relation for the localized synthesisdiffusiondepletion model of the morphogen dynamics in a 1D array of cells. We begin with the reversible reactiondiffusion 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}}_{\mathrm{t}\mathrm{o}\mathrm{t}$). Through the discussion of reversible process, we provide mathematical details on the derivation of the costprecision tradeoff relation for the irreversible reactiondiffusion process of morphogen synthesis and depletion presented in the main text. Next, we derive the conditions that lead to the minimum tradeoff product for the irreversible process with point ($\pi}_{o$) and spacetimeaveraged ($\pi}_{T$) measurements.
We explore the tradeoffs 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}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}$) and size ($l}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}$) at the interval between $x$ and $x+{l}_{\mathrm{cell}}$ by the concentration, $\rho (x,t)[\mathrm{conc}\equiv \mathrm{\#}/{v}_{\mathrm{cell}}]$. At the left boundary (see Appendix 1—figure 1a), the morphogen is injected with rate $j}_{\mathrm{i}\mathrm{n}$$[\mathrm{conc}\times {l}_{\mathrm{cell}}/\mathrm{time}]$, and taken out by the corresponding reverse reaction with rate $v}_{\mathrm{o}\mathrm{u}\mathrm{t}$$[{l}_{\mathrm{cell}}/\mathrm{time}]$. At all positions, the morphogen molecules are depleted with rate ${k}_{\mathrm{d}}[{\mathrm{time}}^{1}]$ and synthesized by the corresponding reverse reaction with rate ${r}_{\mathrm{s}}[\mathrm{conc}/\mathrm{time}]$, and spread across space with diffusivity $D[{l}_{\mathrm{cell}}^{2}/\mathrm{time}]$. At the continuum limit ($L\gg {l}_{\mathrm{cell}}$), the dynamics of the morphogen concentration is described by the partial differential equation
with two boundary conditions at $x=0$ and $x=L$, $D{\mathrm{\partial}}_{x}\rho (x,t){}_{x=0}={j}_{\mathrm{i}\mathrm{n}}{v}_{\mathrm{o}\mathrm{u}\mathrm{t}}\rho (0,t)$ and $D{\mathrm{\partial}}_{x}\rho (x,t){}_{x=L}=0$. $\lambda \equiv \sqrt{D/{k}_{\mathrm{d}}}$ is the characteristic length of the dynamics associated with Equation 16. The steady state concentration is expressed as
where ${\rho}_{\mathrm{ss},0}(x)$ denotes the positiondependent portion of the concentration profile. We introduce three dimensionless parameters: ${\beta}_{\mathrm{in}}\equiv {j}_{\mathrm{in}}/(\lambda {r}_{\mathrm{s}})$, ${\beta}_{\mathrm{out}}\equiv {v}_{\mathrm{out}}/(\lambda {k}_{\mathrm{d}})$, and their ratio $\gamma \equiv {\beta}_{\mathrm{in}}/{\beta}_{\mathrm{out}}$ which is related with the thermodynamic cost (see Equation 32). We only consider the regime with $\gamma >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, ${\rho}_{\mathrm{ss}}(x)\equiv {lim}_{t\to \mathrm{\infty}}\rho (x,t)$, from the initial profile $\rho (x,0)=0$. The time evolution of $\rho (x,t)$ can be analyzed by solving the Equation 16 in Laplace domain:
The expression of the steady state profile (Equation 17) is obtained using
To characterize the relaxation dynamics to the steady state profile, we define the relaxation function
which monotonically decays from 1 to 0 at all positions as $t\to \mathrm{\infty}$. The associated mean relaxation time at position $x$, ${\tau}_{\mathrm{rev}}(x)$, is given by
Since the Laplace transform of the relaxation function is $\widehat{R}(x,s)={\int}_{0}^{\mathrm{\infty}}{e}^{st}R(x,t)\mathit{d}t={s}^{1}\widehat{\rho}(x,s)/{\rho}_{\mathrm{ss}}(x)$, the local accumulation time at position $x$, ${\tau}_{\mathrm{rev}}(x)={lim}_{s\to 0}\widehat{R}(x,s)$, is obtained as
with
and
${\tau}_{\mathrm{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}_{\mathrm{s}}$ and $v}_{\mathrm{o}\mathrm{u}\mathrm{t}$ are negligibly small. At the limits of ${r}_{\mathrm{s}}\to 0$ and ${v}_{\mathrm{out}}\to 0$, with the substitutions ${\beta}_{\mathrm{out}}={v}_{\mathrm{out}}/(\lambda {k}_{\mathrm{d}})$ and $\gamma ={j}_{\mathrm{in}}{k}_{\mathrm{d}}/({v}_{\mathrm{out}}{r}_{\mathrm{s}})$, Equation 22 is simplified to
Precision
The local measure of precision at $x$ is defined as the squared relative positional error (Equation 1),
The expression for the squared relative error in the irreversible process of morphogen dynamics can be obtained at the limits of ${r}_{\mathrm{s}}\to 0$ and ${v}_{\mathrm{out}}\to 0$, with the substitutions ${\beta}_{\mathrm{out}}={v}_{\mathrm{out}}/(\lambda {k}_{\mathrm{d}})$ and $\gamma ={j}_{\mathrm{in}}{k}_{\mathrm{d}}/({v}_{\mathrm{out}}{r}_{\mathrm{s}})$:
Thermodynamic cost
The total entropy production rate required to maintain the steady state morphogen profile is given by
where the three components of the positiondependent entropy production rate are
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 $\underset{{l}_{\mathrm{cell}}\to 0}{lim}D{l}_{\mathrm{cell}}^{2}\left({\rho}_{\mathrm{ss}}(x+{l}_{\mathrm{cell}}){\rho}_{\mathrm{ss}}(x)\right)\mathrm{ln}\left[{\rho}_{\mathrm{ss}}(x+{l}_{\mathrm{cell}})/{\rho}_{\mathrm{ss}}(x)\right]$, in which $D{\rho}_{\mathrm{ss}}(x+{l}_{\mathrm{cell}})/{l}_{\mathrm{cell}}^{2}$ and $D{\rho}_{\mathrm{ss}}(x)/{l}_{\mathrm{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}}_{\mathrm{tot}}=0$), the cells cannot infer any spatial information from it.
To obtain a simplified expression of ${\dot{S}}_{\mathrm{tot}}$, we use the following two equalities. First, at steady state, the net injection of the morphogen at the leftboundary must be equal to the net depletion of the morphogen occurring at all positions, so that
Second, ${\partial}_{x}{\rho}_{\mathrm{ss}}(x)={\partial}_{x}{\rho}_{\mathrm{ss},0}(x)$ and ${\partial}_{x}^{2}{\rho}_{\mathrm{ss}}(x)={\partial}_{x}^{2}{\rho}_{\mathrm{ss},0}(x)={\lambda}^{2}{\rho}_{\mathrm{ss},0}(x)$, with ${\rho}_{\mathrm{ss},0}$ defined in Equation 17, lead to the following equality:
Then, ${\dot{S}}_{\mathrm{tot}}$ can be expressed as follows:
We used 30 and 31 to obtain the second and fourth lines of Equation 32, respectively. The boundary condition of the system, ${{\partial}_{x}\rho (x)}_{x=L}=0$, and the integral of ${\rho}_{\mathrm{ss},0}(x)$ leads to the last line.
The reversible model motivates a physical interpretation of the proportionality constant ${\alpha}_{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}}_{\mathrm{tot}}$ (Equation 32) is normalized by $\mathrm{ln}\gamma $, and ${r}_{\mathrm{s}}$ and $v}_{\mathrm{o}\mathrm{u}\mathrm{t}$ are negligibly small, we obtain
Comparison of Equation 33 with the expression in Equation 5 enables us to relate ${\alpha}_{o}$ with the free energy cost, $T\mathrm{ln}\gamma $. Thus, the thermodynamic cost associated with the morphogen profile is expressed as $C=\left({\alpha}_{o}{v}_{\mathrm{cell}}{j}_{\mathrm{in}}/{l}_{\mathrm{cell}}\right)\tau $ (Equation 5) where τ is the time required to generate the profile, and ${\alpha}_{o}{v}_{\mathrm{cell}}{j}_{\mathrm{in}}/{l}_{\mathrm{cell}}$ approximates the rate of dissipation, $T{\dot{S}}_{\mathrm{tot}}$.
Tradeoff product
The product of the local accumulation time (${\tau}_{\mathrm{rev}}(x)$, Equation 22), the squared relative positional error (${\u03f5}_{\mathrm{rev}}^{2}(x)$, Equation 24), and the entropy production rate (${\dot{S}}_{\mathrm{tot}}$, Equation 32) leads to the positiondependent tradeoff product ${\pi}_{o,\mathrm{rev}}(x)$,
With fixed $L$, $l}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}$, and $x={x}_{\mathrm{b}}$, the lower bound of ${\pi}_{o,\mathrm{rev}}({x}_{\mathrm{b}};\lambda ,{\beta}_{\mathrm{out}},\gamma )$, ${\pi}_{\mathrm{rev},\mathrm{min}}({x}_{\mathrm{b}})$, can be obtained by minimizing ${\pi}_{o,\mathrm{rev}}({x}_{\mathrm{b}};\lambda ,{\beta}_{\mathrm{out}},\gamma )$ with respect to γ, ${\beta}_{\mathrm{out}}$, and λ. For a fixed value of $\gamma >1$, the partial derivative of ${\pi}_{o,\mathrm{rev}}({x}_{\mathrm{b}};\lambda ,{\beta}_{\mathrm{out}},\gamma )$ with respect to ${\beta}_{\mathrm{out}}$ is always smaller than 0,
Thus, ${\pi}_{o,\mathrm{rev}}({x}_{\mathrm{b}})$ is minimized at the limit of large ${\beta}_{\mathrm{out}}$ (or, equivalently, large ${\beta}_{\mathrm{in}}=\gamma {\beta}_{\mathrm{out}}$), with the following expression:
With the characteristic decay length set to $\lambda ={x}_{\mathrm{b}}$, it can be shown that $\underset{\genfrac{}{}{0ex}{}{{x}_{\mathrm{b}}\to 0}{{\beta}_{\mathrm{o}\mathrm{u}\mathrm{t}}\to \mathrm{\infty}}}{lim}{\pi}_{o,\mathrm{r}\mathrm{e}\mathrm{v}}({x}_{\mathrm{b}};{x}_{\mathrm{b}},{\beta}_{\mathrm{o}\mathrm{u}\mathrm{t}},\gamma )=0$. Since ${\pi}_{o,\mathrm{rev},\mathrm{min}}({x}_{\mathrm{b}})\le {\pi}_{o,\mathrm{rev}}({x}_{\mathrm{b}};{x}_{\mathrm{b}},{\beta}_{\mathrm{out}},\gamma )$ by the definition of ${\pi}_{o,\mathrm{rev},\mathrm{min}}({x}_{\mathrm{b}})$, ${\pi}_{o,\mathrm{rev},\mathrm{min}}(0)=0$.
For ${x}_{\mathrm{b}}>0$, the optimal tradeoff product, ${\pi}_{\mathrm{rev},\mathrm{min}}({x}_{\mathrm{b}})$, and the corresponding values of ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})$ and ${\gamma}_{\mathrm{min}}({x}_{\mathrm{b}})$, were numerically obtained (Appendix 1—figure 1). ${\pi}_{o,\mathrm{rev},\mathrm{min}}({x}_{\mathrm{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 ${\pi}_{o,\mathrm{rev}}({x}_{\mathrm{b}};\lambda ,{\beta}_{\mathrm{out}},\gamma )$ close to ${\pi}_{o,\mathrm{rev},\mathrm{min}}({x}_{\mathrm{b}})$, the system can rapidly generate precise morphogen profiles at ${x}_{\mathrm{b}}$ with low thermodynamic cost. In general, ${\pi}_{o,\mathrm{rev},\mathrm{min}}({x}_{\mathrm{b}})$, ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})$, and ${\gamma}_{\mathrm{min}}({x}_{\mathrm{b}})$ all increase monotonically with ${x}_{\mathrm{b}}$ (Appendix 1—figure 1b,c,d).
A cautionary remark is in place. A direct onetoone comparison between the tradeoff products of the reversible and irreversible models should be avoided, since ${\pi}_{o,\mathrm{rev}}$ is minimized at the limits of ${v}_{\mathrm{out}},{r}_{\mathrm{s}}\to \mathrm{\infty}$ (see Equation 35), while the irreversible model is equivalent to the reversible version at the limits of ${v}_{\mathrm{out}},{r}_{\mathrm{s}}\to 0$.
Minimum tradeoff product of the point measurement for the irreversible reactiondiffusion process
We provide the mathematical details involving the tradeoff product of the point measurement (${\pi}_{o}$, Equation 7 and Figure 2 in the main text). With the expressions derived in the Appendix 1, Reversible reactiondiffusion model of morphogen dynamics, the thermodynamic cost of generating the steady state morphogen concentration at position $x$ (Equation 5) amounts to the product among ${\alpha}_{o}$, $\tau (x)$ (Equation 23), and ${v}_{\mathrm{cell}}{j}_{\mathrm{in}}/{l}_{\mathrm{cell}}$ (Equation 33),
and the squared relative positional error, which is identical to Equation 25, is
Then, the tradeoff product reads
For fixed $x={x}_{\mathrm{b}}$, ${\pi}_{o}({x}_{\mathrm{b}};\lambda )$ is numerically minimized with respect to λ to give rise to the solid black line in Figure 2c in the main text. When $L/\lambda \gg 1$, the tradeoff product (Equation 39) is approximated to
For a given position $x={x}_{\mathrm{b}}$, ${{\partial}_{\lambda}{\pi}_{o}({x}_{\mathrm{b}};\lambda )}_{\lambda ={\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})}=0$ determines ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})\approx (\sqrt{13}1){x}_{\mathrm{b}}/6\approx 0.43{x}_{\mathrm{b}}$ (gray dotted line in Figure 2c of the main text).
Minimum tradeoff product of the spacetimeaveraged measurement for the irreversible reactiondiffusion process
As in the main text, the tradeoff product from the spacetimeaveraged measurement is given by
The ${\pi}_{T}(x;\lambda )$ minimizing λ at $x={x}_{\mathrm{b}}$, namely ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})$, is found at the interval $({x}_{\mathrm{b}}a)/2\le {\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})\le {x}_{\mathrm{b}}/2$, since
and
For Figure 3 in the main text, the exact values of ${\pi}_{T,\mathrm{min}}({x}_{\mathrm{b}})$ and ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})$ were obtained by numerically minimizing Equation 41 with respect to λ.
Appendix 2
Distributed morphogen source
We consider a case where the source of morphogen production is distributed over the space in the irreversible reactiondiffusion model of morphogen dynamics. The spacetime evolution of the morphogen concentration is described by
with the boundary conditions ${{\partial}_{x}\rho (x,t)}_{x=0,L}=0$. The production term, ${J}_{\mathrm{in}}(x)$ is assumed to be exponentially distributed with a characteristic length ${\lambda}_{\mathrm{p}}$,
so that ${\int}_{0}^{L}{J}_{\mathrm{in}}(x)\mathit{d}x={j}_{\mathrm{in}}$. The solution to Equation 44 in the Laplace domain is given by
By following a procedure similar to that shown in the section Reversible reactiondiffusion model of morphogen dynamics, we can obtain the following expressions for the steady state profile (${\rho}_{\mathrm{ss}}(x)$, Equation 19), and the precision (${\u03f5}^{2}(x)$, Equation 24):
with $\lambda \equiv \sqrt{D/{k}_{\mathrm{d}}}$. Next, the characteristic timescale to approach the steady state profile ($\tau (x)$, Equation 22) can be written as
with
and
Then, the tradeoff product is expressed as
with the proportionality constant ${\alpha}_{o}$. For fixed $x={x}_{\mathrm{b}}$, $L$, and ${\lambda}_{\mathrm{p}}$, the tradeoff product can be minimized by tuning λ to the optimal value, ${\lambda}_{\mathrm{min}}({x}_{\mathrm{b}})$. As ${\lambda}_{\mathrm{p}}$ increases, the optimal tradeoff product and the associated optimal characteristic length ${\lambda}_{\mathrm{min}}(x)$ decrease (Appendix 2—figure 1). In Appendix 2—figure 1, the characteristic length scale ${\lambda}_{\mathrm{d}}$ is defined as ${\lambda}_{\mathrm{d}}\equiv {\int}_{0}^{L}\mathit{d}x\left({\rho}_{\mathrm{ss}}(x){\rho}_{\mathrm{ss}}(L)\right)/\left({\rho}_{\mathrm{ss}}(0){\rho}_{\mathrm{ss}}(L)\right)$.
Dynamics on a sphere
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,\theta ,\varphi )$ is denoted by $\rho (r,\theta ,\varphi )$, 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}_{\varphi}\rho =0$, ${\partial}_{r}\rho =0$, and $r=R$) the morphogen concentration satisfies the following evolution equation:
with boundary conditions at ${{\partial}_{\theta}\rho }_{\theta =0,\pi}=0$. The morphogen dynamics is fully specified by the system size $R[{l}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}]$, the depletion rate ${k}_{\mathrm{d}}\phantom{\rule{thinmathspace}{0ex}}[1/\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}]$, the diffusion coefficient $D[{l}_{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}^{2}/\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}]$, and the injection rate ${j}_{\mathrm{i}\mathrm{n}}[\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{c}/\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{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=\mathrm{cos}\theta $ and the initial condition $\rho (z,t)=0$, the solution of Equation 51 in the Laplace domain with $\widehat{\rho}(z,s)={\int}_{0}^{\mathrm{\infty}}\mathit{d}t{e}^{st}\rho (z,t)$, can be obtained by solving
The solution to the differential equation of this type can be expressed as $\widehat{\rho}(z,s)={\sum}_{l=0}^{\mathrm{\infty}}{c}_{l}(s){P}_{l}(z)$, where ${P}_{l}(z)$ ($l=0,1,\mathrm{\dots}$) is the Legendre polynomial, satisfying $\mathcal{L}{P}_{l}(x)\equiv \frac{d}{dx}\left[(1{x}^{2})\frac{d{P}_{l}(x)}{dx}\right]=l(l+1){P}_{l}(x)$, with the orthogonality condition ${\int}_{1}^{1}{P}_{l}(x){P}_{m}(x)\mathit{d}x=\frac{2}{2l+1}{\delta}_{lm}$. To find the coefficient ${c}_{l}(s)$, we rearrange Equation 52 into
By using the orthogonality condition of Legendre polynomials, we obtain the solution of Equation 52,
The steady state concentration can be expressed as ${\rho}_{\mathrm{ss}}(z)={lim}_{s\to 0}s\widehat{\rho}(z,s)$.
The average relaxation time can be calculated using $\tau (z)\equiv \u27e8t\u27e9={\int}_{0}^{\mathrm{\infty}}t\left(\frac{dR(t)}{dt}\right)\mathit{d}t={\int}_{0}^{\mathrm{\infty}}R(t)\mathit{d}t={lim}_{s\to 0}\widehat{R}(z,s)$, which yields
with $\lambda \equiv \sqrt{D/{k}_{\mathrm{d}}}$ and
The number of morphogen molecules produced from the cells at the equator per unit time is ${v}_{\mathrm{cell}}(2\pi R/{l}_{\mathrm{cell}}){j}_{\mathrm{in}}$. Thus, the thermodynamic cost to approach the steady state concentration at $z$ is
where ${\alpha}_{o}$ is the proportionality constant.
Next, using ${\rho}_{\mathrm{ss}}(z)={lim}_{s\to 0}s\widehat{\rho}(z,s)$ (Equation 53), one obtains the local precision of the morphogen profile at position $z$,
where
Finally, the tradeoff product, which is plotted in Appendix 2—figure 2b. Eb as a function of $z$ and ${\lambda}_{s}$, is defined as
The length scale ${\lambda}_{\mathrm{s}}$ used in the plot is defined as ${\lambda}_{\mathrm{s}}\equiv {\int}_{z=0}^{z=1}\mathit{d}z\left({\rho}_{\mathrm{ss}}(z){\rho}_{\mathrm{ss}}(1)\right)/\left({\rho}_{\mathrm{ss}}(0){\rho}_{\mathrm{ss}}(1)\right)$. At the target position $z={z}_{\mathrm{b}}$, the tradeoff product can be minimized by tuning λ to the optimal value ${\lambda}_{\mathrm{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 $\lambda}_{\mathrm{B}\mathrm{c}\mathrm{d}}=100\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$ (Houchmandzadeh et al., 2002). The length of the embryo across the anteriorposterior 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\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$. The sensor radius $a=2.6\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{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}_{\mathrm{b}}=(140188)\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$.
Wg. The characteristic decay length of Wg was set to $\lambda}_{\mathrm{W}\mathrm{g}}=6\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$ (Kicheva et al., 2007). The overall size of the DV axis patterned by Wg was set to $L=70\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{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}_{\mathrm{b}}=(1121)\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$.
Hh. The characteristic decay length of Hh was set to $\lambda}_{\mathrm{H}\mathrm{h}}=8\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$ (Figure S2B of Wartlick et al., 2011). The overall size of the domain patterned by Hh was set to $L=100\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{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}_{\mathrm{b}}=(1530)\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$.
Dpp. The characteristic decay length of Dpp was set to $\lambda}_{\mathrm{D}\mathrm{p}\mathrm{p}}=20\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$ (Kicheva et al., 2007). The overall size of the domain patterned by Dpp was set to $L=80\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ (Figure 2H of Bakker et al., 2020). The boundary of Dpp target salm expression was set to $x}_{\mathrm{b}}=(3654)\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$ (Figure 2H of Bakker et al., 2020). The radius of the sensor for Dpp and Hh mediated patterning was set to $a=2\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{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}_{\mathrm{d}}$) of Cyclops and Squint have been measured through the ectopic expressions of the fluorescently tagged constructs. The resulting characteristic decay lengths ($\lambda \equiv \sqrt{D/{k}_{\mathrm{d}}}$) are $\lambda}_{\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{s}}=76\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m$ and $\lambda}_{\mathrm{s}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{n}\mathrm{t}}=178\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{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\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$ (Figure 1 of Nowak et al., 2011). Overall, for Nodal morphogens Cyclops and Squint, ${z}_{\mathrm{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, $\lambda}_{\mathrm{f}\mathrm{g}\mathrm{f}8}=197\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{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}_{\mathrm{b}}$ in the range of $[(70\mu \mathrm{m}70\mu \mathrm{m})/200\mu \mathrm{m},(120\mu \mathrm{m}70\mu \mathrm{m})/200\mu \mathrm{m}]=[0.0,0.25]$.
Time scales of Bcd, Wg, Hh, and Dpp
The expression for the variance of the spacetimeaveraged signal, ${\sigma}_{m}^{2}(x)$ (Equation 11), is derived with the assumption of a long measurement time, $T\gg {k}_{\mathrm{d}}^{1}$ (Appendix 1 of Fancher and Mugler, 2020). Additionally, in order for $T{j}_{\mathrm{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\gg \tau ({x}_{\mathrm{b}})$). For Bcd, Wg, Hh, and Dpp, $L\gg \lambda $ and ${x}_{\mathrm{b}}/\lambda \approx 2$, leading to the approximate expression for the characteristic time of $\tau ({x}_{\mathrm{b}})\approx 1.5{k}_{\mathrm{d}}^{1}$ (Equation 4). We assume that the maximum measurement time, ${T}_{\mathrm{max}}$, is set by the cell doubling time, after which the fates of the daughter cells must be reestablished.
Bcd: $T}_{\mathrm{m}\mathrm{a}\mathrm{x}}[\approx 10\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{i}\mathrm{n$(Foe and Alberts, 1983)] $<{k}_{\mathrm{d}}^{1}[\approx 1\phantom{\rule{thinmathspace}{0ex}}\mathrm{h}\mathrm{r}$(Drocco et al., 2011)].
Wg: $T}_{\mathrm{m}\mathrm{a}\mathrm{x}}[\approx 5\phantom{\rule{thinmathspace}{0ex}}\mathrm{h}\mathrm{r$(Martín et al., 2009)] $\gg {k}_{\mathrm{d}}^{1}[\approx 15\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{i}\mathrm{n}$(Kicheva et al., 2007)].
Hh: $T}_{\mathrm{m}\mathrm{a}\mathrm{x}}[\approx 5\phantom{\rule{thinmathspace}{0ex}}\mathrm{h}\mathrm{r$(Martín et al., 2009)] $\gg {k}_{\mathrm{d}}^{1}[\approx 5\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{i}\mathrm{n}$(Nahmad and Stathopoulos, 2009)].
Dpp: $T}_{\mathrm{m}\mathrm{a}\mathrm{x}}[\approx 5\phantom{\rule{thinmathspace}{0ex}}\mathrm{h}\mathrm{r$(Martín et al., 2009)] $>{k}_{\mathrm{d}}^{1}[\approx 1.4\phantom{\rule{thinmathspace}{0ex}}\mathrm{h}\mathrm{r}$(Kicheva et al., 2007)].
For Dpp, Hh, and Wg, the spacetimeaveraged measurements may be performed for sufficiently long durations to represent the costprecision tradeoff by ${\pi}_{T}$ (Equation 14). For Bcd, ${\pi}_{o}$, which is formulated for the point measurement, may better represent the pertinent costprecision tradeoff.
Data availability
All data analyzed in this study are from the figures of previously published works, shown in Figure 2S1 and Appendix 2 Figure 3. The reference associated with each panel is provided in the respective figure legend.
References

Fly wing vein patterns have spatial reproducibility of a single cellJournal of the Royal Society Interface 11:20140443.https://doi.org/10.1098/rsif.2014.0443

Scaleinvariant patterning by sizedependent inhibition of nodal signallingNature Cell Biology 20:1032–1042.https://doi.org/10.1038/s4155601801557

Thermodynamic uncertainty relation for biomolecular processesPhysical Review Letters 114:1158101.https://doi.org/10.1103/PhysRevLett.114.158101

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

Formation of morphogen gradients: local accumulation timePhysical Review E 83:051906.https://doi.org/10.1103/PhysRevE.83.051906

Physics of chemoreceptionBiophysical Journal 20:193–219.https://doi.org/10.1016/S00063495(77)855446

Measurement and perturbation of morphogen lifetime: effects on gradient shapeBiophysical Journal 101:1807–1815.https://doi.org/10.1016/j.bpj.2011.07.025

Optimizing the readout of morphogen gradientsPhysical Review E 77:041903.https://doi.org/10.1103/PhysRevE.77.041903

Role of spatial averaging in the precision of gene expression patternsPhysical Review Letters 103:258101.https://doi.org/10.1103/PhysRevLett.103.258101

Information thermodynamics of turing patternsPhysical Review Letters 121:108301.https://doi.org/10.1103/PhysRevLett.121.108301

Dissipation bounds all SteadyState current fluctuationsPhysical Review Letters 116:120601.https://doi.org/10.1103/PhysRevLett.116.120601

Shaping a morphogen gradient for positional precisionBiophysical Journal 99:697–707.https://doi.org/10.1016/j.bpj.2010.04.073

Grand canonical markov model: a stochastic theory for open nonequilibrium biochemical networksThe Journal of Chemical Physics 124:044110.https://doi.org/10.1063/1.2165193

Energetic costs, precision, and transport efficiency of molecular motorsThe Journal of Physical Chemistry Letters 9:513–520.https://doi.org/10.1021/acs.jpclett.7b03197

Dynamics of the dorsal morphogen gradientPNAS 106:21707–21712.https://doi.org/10.1073/pnas.0912395106

Formation of a morphogen gradient: acceleration by degradationThe Journal of Physical Chemistry Letters 2:1502–1505.https://doi.org/10.1021/jz2004914

Robust and precise morphogenmediated patterning: tradeoffs, constraints and mechanismsJournal of the Royal Society, Interface 12:20141041.https://doi.org/10.1098/rsif.2014.1041

3 minutes to precisely measure morphogen concentrationPLOS Genetics 14:e1007676.https://doi.org/10.1371/journal.pgen.1007676

Interpretation of the FGF8 morphogen gradient is regulated by endocytic traffickingNature Cell Biology 13:153–158.https://doi.org/10.1038/ncb2155

Precision of hunchback expression in the Drosophila embryoCurrent Biology 22:2247–2252.https://doi.org/10.1016/j.cub.2012.09.051

Contribution of increasing plasma membrane to the energetic cost of early zebrafish embryogenesisMolecular Biology of the Cell 31:520–526.https://doi.org/10.1091/mbc.E19090529

Modeling of Wntmediated tissue patterning in vertebrate embryogenesisPLOS Computational Biology 16:e1007417.https://doi.org/10.1371/journal.pcbi.1007417

Morphogen profiles can be optimized to buffer against noisePhysical Review E 80:041902.https://doi.org/10.1103/PhysRevE.80.041902

Mathematical models of morphogen gradients and their effects on gene expressionWiley Interdisciplinary Reviews: Developmental Biology 1:715–730.https://doi.org/10.1002/wdev.55

Energy budget of Drosophila embryogenesisCurrent Biology 29:R566–R567.https://doi.org/10.1016/j.cub.2019.05.025

Thermodynamic cost, speed, fluctuations, and error reduction of biological copy machinesThe Journal of Physical Chemistry Letters 11:3136–3143.https://doi.org/10.1021/acs.jpclett.0c00545

Thermodynamic uncertainty relation to assess biological processesThe Journal of Chemical Physics 154:130901.https://doi.org/10.1063/5.0043671

Development of morphogen gradient: the role of dimension and discretenessThe Journal of Chemical Physics 140:085102.https://doi.org/10.1063/1.4866453

New model for understanding mechanisms of biological signaling: direct transport via cytonemesThe Journal of Physical Chemistry Letters 7:180–185.https://doi.org/10.1021/acs.jpclett.5b02703

The many bits of positional informationDevelopment 148:dev176065.https://doi.org/10.1242/dev.176065

Fundamental limits to position determination by concentration gradientsPLOS Computational Biology 3:e78.https://doi.org/10.1371/journal.pcbi.0030078

The chemical basis of morphogenesis. 1953Bulletin of Mathematical Biology 52:153–197.https://doi.org/10.1016/S00928240(05)800084

Positional information and the spatial pattern of cellular differentiationJournal of Theoretical Biology 25:1–47.https://doi.org/10.1016/S00225193(69)800160
Decision letter

Gordon J BermanReviewing Editor; Emory University, United States

Naama BarkaiSenior Editor; Weizmann Institute of Science, Israel

Mariela PetkovaReviewer; Harvard University, United States

Anatoly KolomeiskyReviewer; Rice University, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
This manuscript presents a theoretical investigation on costprecision relations in the formation of biological signaling profiles (also known as morphogen gradients) that are responsible for the proper development of biological organisms. It is argued that the thermodynamic cost is an important factor in the morphogen gradient formation, and several signaling molecule gradients in Drosophila appear to be optimized according to the presented theoretical estimates. This is one of the first quantitative studies of the optimization of the development morphogen gradients, and its main arguments are strong and convincing. One of the main conclusions of this work is that the complex task of biological signaling processes seems to be optimized, and the paper should be interesting to a wide range of scientists working on both experimental and theoretical aspects of understanding the fundamental origin of optimization in biological processes.
Decision letter after peer review:
Thank you for submitting your article "Costprecision tradeoff relation determines the optimal morphogen gradient for accurate biological pattern formation" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Mariela Petkova (Reviewer #2); Anatoly Kolomeisky (Reviewer #3).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential Revisions:
1) The main assumption in the paper is that the cost of creating a morphogen gradient is something that the biological system cares about. However, it is not shown that the cost of creating the morphogen gradient is a significant portion of the total energy budget of the biological system, and therefore potentially necessary to be optimized. The choice of tradeoff function was presented without supporting arguments that either reflect biological expectations or demonstrate that the results are robust to other tradeoff constructs. Importantly, the paper argues that the agreement between the optimal gradient shape predicted by the tradeoff framework and the measurements in real biological systems indicates that thermodynamic cost is a key constraint determining the morphogen shape. However, perturbed systems such as increasing and decreasing the dosage of Bicoid still produce viable embryos (Liu et al., PNAS, 2013). The latter observation suggests that both producing different amount of morphogen molecules and employing regulatory mechanisms to correct the effects of perturbations are not energetically costlimiting to an embryo. Justifying that this cost is significant with respect to the cells' energy budget would greatly aid the manuscript.
2) Similarly, the authors choose a product of thermodynamic cost and the squared relative error in the position as a tradeoff function that they analyze. My question is how robust are the obtained results? The choice of the tradeoff function is not unique and it does not have any fundamental foundation, to the best of my understanding. For instance, one defines X as a thermodynamic cost and Y as relative error in position, then not only X*Y, but also (X+Y), X*Y*Y, X*X*Y and many other combinations also might be viewed as tradeoff functions. But they would lead to a different prediction on the optimized minimal length and other related issues. Then how to choose and why?
3) The work constraints the patterning problem to a gradient which produces a single target boundary, but validates the results with measurements in systems where gradients are known to regulate many downstream patterns. The result of the tradeoff analysis is that the location of the target boundary determines the shape of the gradient. How is this reconciled with the understanding that, for example, the Bicoid gradient is read out by multiple downstream genes along the entire length of the fly embryo? (Hannon et al., eLife, 2017).
4) The description of the error above Equation 1 is somewhat confusing in two different ways: i) definitions are given in terms of mathematical variables like N and the subscript i that are not explicitly explained. The reader can figure out the intended meaning via later discussions, but it would be helpful to make these terms clear from the start. ii) More broadly, the notion of "measurement" used in the paper is a little abstract, and is likely to confuse biologicallyoriented readers. Before delving into the mathematical definition of measurement errors, could the authors provide a summary of how measurements could hypothetically be carried out via concrete biological mechanisms? For example, they mention that spacetimeaveraged measurements might be carried out by cell surface receptors. But what mechanisms could lead to instantaneous point measurements?
5) In line 181 the authors state that "the two limiting models can be effectively minimized by similar λ values at a given position is of great significance." However comparing Figures 2(c) and 3(c) it is a hard to directly see this, because in one case the yaxis is normalized as λ/L and in the other case as λ/a. Could the authors somehow make this similarity in predicted λ values clearer? Maybe adding a twin y axis scale to one of the figures, showing λ in terms of the other normalization?
6) The discussion of the reversible model in the appendix is quite nice, particularly how you can derive all the irreversible results directly from the reversible model by taking an appropriate limit. Particularly interesting is that the cost defined in terms of entropy production rate times characteristic time, if normalized by ln(γ) [using Equation 33] would seem to give you the morphogen production cost of Equation 37, up to a factor of α_{0}, assuming the irreversible limit where ln(γ) > infinity. So it seems that the entropybased cost is closely related to the other definition of cost, and the authors could consider highlighting this fact in the main text. Is there any qualitative way of understanding why ln(γ) is the appropriate normalization factor?
7) All derivations are mode for the finite length L. However, as argued in the specific calculations L>>\λ, so why not to derive the results for L>>1. In this case, all the formulas are much simpler and there is a clear physical interpretation for all the quantities. Why the bulkier and less transparent expressions are used in the text? If authors wanted to explain things better to a broader readership and to emphasize the physics they better use a simple analytical framework. Or they need to explain why they used a finite L.
8) It would be better to describe fully how the characteristic lengths were estimated. This is because the reported numbers do not fully agree with what is reported in experiments. For example, for BCD Berezhkovskii and Shvartsman report λ=75 nm, not 100nm as given in this paper. Why are such different results are reported?
Reviewer #1:
There has been a recent surge of interest in trying to understand the tradeoff between the precision achieved by nonequilibrium processes and the underlying thermodynamic costs. One of the most interesting arenas for this exploration has been biology, and in the current paper, the authors have discovered a novel form of this tradeoff in the context of pattern formation during animal development.
Their main focus is a 1D mathematical model describing how morphogens are produced, diffuse, and depleted within the cell, forming a gradient that divides the cell into two spatial regions with different subsequent developmental trajectories. Establishing a sharp boundary between these regions comes at a certain cost in terms of morphogen production, and the authors provide an elegant theoretical argument that the boundary precision and cost bounded by an "uncertainty principle" type relationship, where their product is always greater than a certain systemspecific quantity. The result is reminiscent of, but qualitatively different from, the socalled thermodynamic uncertainty principle that has generated a lot of excitement in the nonequilibrium statistical physics community in the last few years.
Finding such a tradeoff is interesting in itself, but is it biologically relevant? In this respect the paper makes two striking claims: (i) making the system optimally cost effective in achieving a certain precision predicts a certain characteristic decay length for the morphogen concentration profile; (ii) this prediction is directly testable using data from earlier Drosophila development experiments. The authors' analysis shows how four morphogens (Hh, Wg, Bcd, Dpp) to a good approximation fulfill the criteria for optimality in this sense. This conclusion seems robust to different mechanisms by which cells can "measure" the morphogen profile to translate it into cell fate decisions, whether these measurements are done instantaneously at precise locations or averaged over time and space.
Overall the results of the paper are well supported by the theoretical analysis and data. The optimality of the four Drosophila morphogens suggests that costeffectiveness has played an evolutionary role in shaping morphogenesis, at least in this organism. One complicating factor in this conclusion (which to their credit, is explored by the authors) is that optimality may not be achieved in systems with more complex geometries, like zebrafish. How and when thermodynamic costs are decisive in biological contexts is a fascinating question for future study. This paper is a significant contribution to that discussion and a wonderful example of physicsinspired theory helping us interpret wellstudied biological systems in a completely new way.
Reviewer #2:
This paper assumes that the thermodynamic cost of producing patterning molecules constraints the events which occur during development. It takes a minimalistic view of a general class of problems in developmental systems: how a morphogen gradient specifies a boundary location in a downstream gene expression pattern. The authors evaluate the cost of generating a morphogen gradient against the precision with which the gradient can specify the target boundary. Making more morphogen molecules costs more, but reduces variability in the gradient and improves the precision with which a target boundary can be set up. The authors develop a tradeoff framework for two gradient readout limits: instantaneous and spatiotemporally averaged. In both cases, the optimal solution for the shape of the morphogen gradient is consistent with experimental measurements in multiple biological systems. The key claim of the work is that this agreement between prediction and measurement is evidence that developing systems employ such tradeoffs.
Strengths:
The tradeoff framework developed in this paper has three major strengths. First, it is simple and generalizable to multiple systems and geometries: developing embryos in multiple species as well as gradientmediated tissue patterning. Second, it depends only on geometric features such as length, boundary location and gradient shape. These can be readily measured in real biological systems and thus test the framework's predictions. Lastly, it captures how the dynamics of reading out the morphogen gradient affects the tradeoff between morphogen cost and target boundary precision. The key result of the predictions is that the same shape of the morphogenic gradient optimizes two readout strategies which are at the two extremes of biological function – instantaneous readout, and spatiotemporal averaging.
Weaknesses:
The weakness in this paper lies in reconciling the work's assumptions and results with their interpretations in the context of complex biological systems.
The main assumption in the paper is that the cost of creating a morphogen gradient is something that the biological system cares about. However, it is not shown that the cost of creating the morphogen gradient is a significant portion of the total energy budget of the biological system, and therefore potentially necessary to be optimized.
Next, the work constraints the patterning problem to a gradient which produces a single target boundary, but validates the results with measurements in systems where gradients are known to regulate many downstream patterns. The result of the tradeoff analysis is that the location of the target boundary determines the shape of the gradient. How is this reconciled with the understanding that, for example, the Bicoid gradient is read out by multiple downstream genes along the entire length of the fly embryo? (Hannon et al.,2017)
Finally, the authors present only one possibility for the tradeoff function and evaluate it against measurements in unperturbed systems. The choice of tradeoff function was presented without supporting arguments that either reflect biological expectations or demonstrate that the results are robust to other tradeoff constructs. Importantly, the paper argues that the agreement between the optimal gradient shape predicted by the tradeoff framework and the measurements in real biological systems indicates that thermodynamic cost is a key constraint determining the morphogen shape. However, perturbed systems such as increasing and decreasing the dosage of Bicoid still produce viable embryos (Liu et al., 2013). The latter observation suggests that both producing different amount of morphogen molecules and employing regulatory mechanisms to correct the effects of perturbations are not energetically costlimiting to an embryo.
References:
Hannon, C.E., Blythe, S.A. and Wieschaus, E.F., 2017. Concentration dependent chromatin states induced by the bicoid morphogen gradient. eLife, 6, p.e28275.
Liu, F., Morrison, A.H. and Gregor, T., 2013. Dynamic interpretation of maternal inputs by the Drosophila segmentation gene network. Proceedings of the National Academy of Sciences, 110(17), pp.67246729
Reviewer #3:
This is a very elegant theoretical investigation on the microscopic features of the formation of biological signaling profiles. A comprehensive theoretical model that estimates the thermodynamic cost and the accuracy is presented. The theoretical predictions are utilized for analyzing several morphogens in Drosophila. The strengths of the paper are the following: (i) it addresses a fundamental question of optimization of biological systems from the fundamental point of view, (ii) the approach is fully quantitative, (iii) the theoretical framework is explored for estimating the optimization of morphogen gradients in the real biological system, and (iv) several extensions for more realistic situations of the morphogen gradients formation are outlined. All these features make me believe that the conclusions of the paper are justified. However, despite the fact that the paper is well written and the overall results seem convincing, there are several issues that need to be clarified. The weaknesses (not very significant) include: (i) the estimates of lengths do not fully agree with some of the reported lengths in the literature, (ii) the robustness of the costprecision tradeoff relation is not critically evaluated, and (iii) some references to the existing literature are missing.
https://doi.org/10.7554/eLife.70034.sa1Author response
Essential Revisions:
1) The main assumption in the paper is that the cost of creating a morphogen gradient is something that the biological system cares about. However, it is not shown that the cost of creating the morphogen gradient is a significant portion of the total energy budget of the biological system, and therefore potentially necessary to be optimized. The choice of tradeoff function was presented without supporting arguments that either reflect biological expectations or demonstrate that the results are robust to other tradeoff constructs. Importantly, the paper argues that the agreement between the optimal gradient shape predicted by the tradeoff framework and the measurements in real biological systems indicates that thermodynamic cost is a key constraint determining the morphogen shape. However, perturbed systems such as increasing and decreasing the dosage of Bicoid still produce viable embryos (Liu et al., PNAS, 2013). The latter observation suggests that both producing different amount of morphogen molecules and employing regulatory mechanisms to correct the effects of perturbations are not energetically costlimiting to an embryo. Justifying that this cost is significant with respect to the cells' energy budget would greatly aid the manuscript.
Thank you for raising these important points. As demonstrated in Liu et al., it is difficult to predict the effect of perturbing the Bcd profile on the organism’s fitness, because the embryonic anterior patterning is governed by a large network of coregulating signaling molecules. Additionally, as detailed in the revised manuscript, the thermodynamic cost of generating the Bcd morphogen profile is only a small fraction of the total energy budget of embryogenesis. 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., Dev. Cell 2019; Rodenfels el al., MBoC 2020; Song et al., Curr. Biol. 2019). The generation of each morphogen profile is an indispensable developmental process that contributes to the total energy budget, if only by a small amount. In this light, our theory proposes a quantitative framework to evaluate the costprecision tradeoff 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. It would be fruitful in future work to assess the costprecision tradeoff of more comprehensive models by incorporating the thermodynamic costs of other components of the signaling network underlying biological pattern formation. These discussion points are added throughout the revised manuscript.
– Rodenfels, J., Neugebauer, K. M., and Howard, J. (2019). Heat Oscillations Driven by the Embryonic Cell Cycle Reveal the Energetic Costs of Signaling. Dev. Cell, 48(5):646–658.
– Rodenfels, J., Sartori, P., Golfier, S., Nagendra, K., Neugebauer, K., and Howard, J. (2020). Contribution of increasing plasma membrane to the energetic cost of early zebrafish embryogenesis. MBoC, 31(7):520–526.
– Song, Y., Park, J. O., Tanner, L., Nagano, Y., Rabinowitz, J. D., and Shvartsman, S. Y. (2019). Energy budget of Drosophila embryogenesis. Curr. Biol., 29(12):R566–R567
Our response to the choice of tradeoff function and its robustness is given in 2. (please see below).
As mentioned above, the revised manuscript contains a more detailed discussion on the cost of generating the morphogen profile. The proportionality constant, ${\alpha}_{o}$, in the expression of the cost $C\equiv {\alpha}_{o}\left({v}_{\mathrm{\text{cell}}}{j}_{\mathrm{\text{in}}}\mathrm{\text{/}}{l}_{\mathrm{\text{cell}}}\right)\tau $ (Equation 5) represents the thermodynamic cost (free energy) involving 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 for the synthesis and degradation of a single Bcd protein composed of 494 amino acids is estimated to be on the order of ${\alpha}_{o}\approx 494\times 4\approx 2\times {10}^{3}$ ATPs, where we assume that each peptidebond formation requires 4 ATPs (Lynch et al., 2015). During the first ~2 hours of Drosophila embryogenesis when the anteriorposterior axis patterning takes place, about ~$5\times {10}^{8}$ Bcd molecules are produced (Drocco et al., 2012). Thus, the cost of generating the Bcd profile is on the order of 10^{12} ATPs, which is only a small fraction of the total energy budget of Drosophila embryogenesis ~$7\times {10}^{16}\phantom{\rule{0.222em}{0ex}}$ ATPs (Song et al., 2019).
However, the energetic cost of numerous individual molecular processes comprise the overall energy budget of an organism. Thus, it is of great interest to quantify the energetic efficiency of each molecular process. For instance, molecular motors (e.g. kinesin) and biomassproducing enzymes (e.g. ribosome) balance multiple functional features, such as the thermodynamic cost, reaction speed, and the fluctuations of the reactions. In light of the thermodynamic uncertainty relation, the product $Q\equiv \Delta S\left(t\right)\frac{\u27e8\delta X{\left(t\right)}^{2}\u27e9}{{\u27e8X\left(t\right)\u27e9}^{2}}$ quantifies the tradeoffs among the thermodynamic cost ($\Delta S\left(t\right)$), speed (${\u27e8X\left(t\right)\u27e9}^{2}$), and the fluctuations ($\u27e8\delta X{\left(t\right)}^{2}\u27e9$) of the timedependent observable $X\left(t\right)$, which can be defined as the steps taken by the molecular motor or the number of peptidebonds synthesized by the ribosome. Previous works by Hyeon and colleagues have shown that $Q$ of certain molecular processes are close to the physical lower bound of $2{k}_{B}$ ($Q\approx 715{k}_{B}$ for the molecular motors, and $Q\approx 4550{k}_{B}$ for E. coli ribosome), which suggests that those molecular processes have evolved to operate costeffectively (Hwang and Hyeon, JPCL 2018; Song and Hyeon, JPCL 2020; Song and Hyeon, JCP 2021). Similarly, in our study, it is remarkable to find that the fruit fly embryos are nearly optimized with respect to the costprecision tradeoff product ($\pi $), which suggests that both the cost and the precision are key physical features that determine the shape of morphogen profiles at the time scale of evolution.
While it is possible to assign explicit values to the cost of generating the morphogen profile, $C$, the interpretation of our theoretical result in the context of the biological system must be done with care. For the Bcd profile with the boundary position ${x}_{b}\approx 0.4L$, ${\pi}_{o}({x}_{b},{\lambda}_{B\mathrm{\text{cd}}})\approx {\alpha}_{o}\left(L\mathrm{\text{/}}{l}_{\mathrm{\text{cell}}}\right)0.6$ (Figure 2b) where $L\mathrm{\text{/}}{l}_{\mathrm{\text{cell}}}\approx 50,$ and the experimentally reported precision is ${\u03f5}^{2}\left({x}_{b}\right)\approx 7\times {10}^{4}$ (Gregor et al., 2007). Thus, the cost associated with the Bcd profile is $C({x}_{b},{\lambda}_{B\mathrm{\text{cd}}})={\pi}_{o}({x}_{b},{\lambda}_{B\mathrm{\text{cd}}})\mathrm{\text{/}}{\u03f5}^{2}({x}_{b},{\lambda}_{B\mathrm{\text{cd}}})\approx \phantom{\rule{0.222em}{0ex}}40000{\alpha}_{o}$, or equivalently, the thermodynamic cost of synthesizing and degrading 40000 Bcd molecules. This is orders of magnitude smaller than the estimate of $5\times {10}^{8}$ based on the direct measurement of the 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 the 3D volume including the cytoplasm and the nuclei. Additionally, the duration to generate the steady state morphogen profile is longer than the characteristic time, $\tau $. However, assuming that the cytoplasmic and nuclear concentrations are proportional to each other, and that the time for the morphogen profile to reach steady states is proportional to $\tau $, it is possible to adjust the proportionality constant ${\alpha}_{o}$ so that $C\left({x}_{b}\right)$ represents the overall cost to produce the concentration profile of nuclear Bcd.
– Lynch, M. and Marinov, G. K. (2015). The bioenergetic costs of a gene. Proc. Natl. Acad. Sci. USA, 112(51):15690– 15695.
– Gregor, T., Tank, D. W., Wieschaus, E. F., and Bialek, W. (2007). Probing the Limits to Positional Information. Cell, 130(1):153–164.
– Drocco, J. A., Wieschaus, E. F., and Tank, D. W. (2012). The synthesisdiffusiondegradation model explains Bicoid gradient formation in unfertilized eggs. Phys. Biol., 9(5):055004.
– Hwang, W. and Hyeon, C. (2018). Energetic Costs, Precision, and Transport Efficiency of Molecular Motors. J. Phys. Chem. Lett., 9(3):513–520.
– Song, Y. and Hyeon, C. (2020). Thermodynamic cost, speed, fluctuations, and error reduction of biological copy machines. J. Phys. Chem. Lett., 11:3136–3143.
– Song, Y. and Hyeon, C. (2021). Thermodynamic uncertainty relation to assess the efficiency of biological processes. J. Chem. Phys., 154:130901.
2) Similarly, the authors choose a product of thermodynamic cost and the squared relative error in the position as a tradeoff function that they analyze. My question is how robust are the obtained results? The choice of the tradeoff function is not unique and it does not have any fundamental foundation, to the best of my understanding. For instance, one defines X as a thermodynamic cost and Y as relative error in position, then not only X*Y, but also (X+Y), X*Y*Y, X*X*Y and many other combinations also might be viewed as tradeoff functions. But they would lead to a different prediction on the optimized minimal length and other related issues. Then how to choose and why?
The main reason for defining the tradeoff by $\pi \equiv C\times {\u03f5}^{2}$ was the following:
There are two quantities of interest in the problem of costprecision tradeoff of morphogen profile formation. First, researchers interested in the precision of morphogen concentration at the target boundary have been measuring the quantity which amounts to “the relative error of morphogen concentration” defined to be proportional to the ratio of concentration fluctuation to its mean value, namely, ${\sigma}_{\rho}/\u27e8\rho \u27e9$. In fact, this quantity, often used in statistics, is proportional to $1/{n}_{m}^{1/2}$, which is nothing but the outcome of central limit theorem. Second, the cost of producing morphogens ($C$), which is the main concern in our study, is expected to be proportional to the number of morphogens being produced (${n}_{\mathrm{\text{m}}}).$ – Let us clarify here that the cost ($C)$ we refer to in our study is specific to the morphogen, not the whole cellular budget. – Then, when the cost of morphogen production is large, the precision of morphogen profile at the target boundary is high (or the relative error of morphogen concentration at the target boundary is small). Conversely, when the cost is small, the precision is low (or the error is large). This means that there is a tradeoff between the cost and precision. The simplest and generic construction for discussing the tradeoff relationship between the two quantities, independently from ${n}_{\mathrm{\text{m}}}$ (i.e. the system size) is to consider the product between the cost ($C$) and the squared value of ${\sigma}_{\rho}/\u27e8\rho \u27e9$.
Let us emphasize again that $\pi =C\times {\u03f5}^{2}$ , which we called the tradeoff product, is a dimensionless quantity, independent of the number of morphogens being produced (${n}_{\mathrm{\text{m}}}$).
Other constructions suggested by the referee(s), such as X+Y, X*Y*Y, and X*X*Y do not fulfill this requirement, and the minimization of these constructions leads to physically uninterpretable conclusions.
In fact, this form of the costprecision tradeoff was motivated by the thermodynamic uncertainty relation (TUR), which defines the product between the squared relative error of the time integrated currentlike observable $X\left(t\right)$ and the total entropy production for $t$, $\mathrm{\text{\Delta S}}\left(t\right)$, namely, $Q\equiv \Delta S\left(t\right)\frac{\u27e8\delta X{\left(t\right)}^{2}\u27e9}{{\u27e8X\left(t\right)\u27e9}^{2}}$. In the case of TUR, the two quantities of interest are a function of time t, such that $\Delta S\left(t\right)\sim t$, $\frac{\u27e8\delta X{\left(t\right)}^{2}\u27e9}{{\u27e8X\left(t\right)\u27e9}^{2}}\sim 1/t$, and the dependence of time t is cancelled off when the product of the two is taken. What is remarkable in TUR is that there exists a model independent universal bound in the product, $Q\ge 2{k}_{B}$, where ${k}_{B}$ is the Boltzmann constant, and the costeffectiveness of a dynamical process generated in nonequilibrium is dictated by how close the parameter Q is to its minimal bound.
Regarding the robustness of our result, we would like to emphasize this: For the scenarios of the point measurement and the spacetime averaging, we addressed the costprecision tradeoff of the morphogen gradient formation by employing slightly different definitions of the cost and the precision. We obtained similar conclusions as to the optimal characteristic lengths, ${\lambda}_{min}\approx 0.43{x}_{\mathrm{\text{b}}}\phantom{\rule{0.222em}{0ex}}$ for the point measurement and ${\lambda}_{min}\approx 0.5{x}_{\mathrm{\text{b}}}$ for the spacetimeaveraged measurement, both of which are the outcomes of minimizing the respective tradeoff product.
3) The work constraints the patterning problem to a gradient which produces a single target boundary, but validates the results with measurements in systems where gradients are known to regulate many downstream patterns. The result of the tradeoff analysis is that the location of the target boundary determines the shape of the gradient. How is this reconciled with the understanding that, for example, the Bicoid gradient is read out by multiple downstream genes along the entire length of the fly embryo? (Hannon et al., eLife, 2017).
Our theory implies that the precision of the morphogen induced spatial boundary can be optimized costeffectively at only one spatial position. For the morphogens Bcd, Dpp, and Wg, our analysis focused on the target gene expressions whose quantitative characterization revealed a relatively sharp boundary (Figure 2S1). In the case of Hh, we included the full range of spatial boundaries formed by the three target genes, dpp, col, and en (Figure 2S1). As the referees point out, however, it is conceivable that the precision of the morphogen profile is essential at multiple boundary positions. 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 required to form the profile. This discussion is added to the revised manuscript.
4) The description of the error above Equation 1 is somewhat confusing in two different ways: i) definitions are given in terms of mathematical variables like N and the subscript i that are not explicitly explained. The reader can figure out the intended meaning via later discussions, but it would be helpful to make these terms clear from the start. ii) More broadly, the notion of "measurement" used in the paper is a little abstract, and is likely to confuse biologicallyoriented readers. Before delving into the mathematical definition of measurement errors, could the authors provide a summary of how measurements could hypothetically be carried out via concrete biological mechanisms? For example, they mention that spacetimeaveraged measurements might be carried out by cell surface receptors. But what mechanisms could lead to instantaneous point measurements?
Thank you for this suggestion. In the revised manuscript, we clarified that Bcd concentration at each nucleus 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 its positiondependent expression profile. Additionally, cells with many receptors positioned across the cell surface may perform spacetimeaveraging by integrating the signal from morphogenbound receptors for an extended time interval. Theoretically, cells may control the duration of the measurement by modulating the availability of morphogensensing receptors, or by tuning the overall transcription rate of the target gene. Then, the point measurement is effectively interpreted as the short time limit of the spacetime averaged counterpart with a single receptor. Alternatively, one can define the point measurement of the morphogen profile as an image of the morphogen profile obtained from a fixed sample.
5) In line 181 the authors state that "the two limiting models can be effectively minimized by similar λ values at a given position is of great significance." However comparing Figures 2(c) and 3(c) it is a hard to directly see this, because in one case the yaxis is normalized as λ/L and in the other case as λ/a. Could the authors somehow make this similarity in predicted λ values clearer? Maybe adding a twin y axis scale to one of the figures, showing λ in terms of the other normalization?
Our main finding can be summarized by the expression of the optimal characteristic length $\lambda $ that minimizes the tradeoff product at a given target boundary. In Figure 3c of the revised manuscript, we have overlaid the linear approximation of the optimal characteristic length ${\lambda}_{\mathrm{\text{min}}}$ from the point measurement model with ${\lambda}_{\mathrm{\text{min}}}$ from the model with spacetimeaveraging, so that the readership can compare the two results sidebyside. Additionally, in the discussion, we clarified that ${\lambda}_{\mathrm{\text{min}}}$ of the two models were similar, in that ${\lambda}_{\mathrm{\text{min}}}\approx 0.43L$ in the former and ${\lambda}_{\mathrm{\text{min}}}\approx 0.5L$ in the latter.
6) The discussion of the reversible model in the appendix is quite nice, particularly how you can derive all the irreversible results directly from the reversible model by taking an appropriate limit. Particularly interesting is that the cost defined in terms of entropy production rate times characteristic time, if normalized by ln(γ) [using Equation 33] would seem to give you the morphogen production cost of Equation 37, up to a factor of α_{0}, assuming the irreversible limit where ln(γ) > infinity. So it seems that the entropybased cost is closely related to the other definition of cost, and the authors could consider highlighting this fact in the main text. Is there any qualitative way of understanding why ln(γ) is the appropriate normalization factor?
We appreciate your suggestion. The reversible model motivates a physical interpretation of the proportionality constant, ${\alpha}_{o}$, which we defined in the irreversible model of morphogen dynamics as the number of ATPs hydrolyzed for the synthesis and degradation of a morphogen molecule. When the entropy production rate (${\stackrel{\u0307}{S}}_{\mathrm{\text{tot}}}$, Equation 32) is normalized by $ln\left(\gamma \right)$, and ${v}_{\mathrm{\text{out}}}$ and ${r}_{\mathrm{\text{s}}}$ are negligibly small, we obtain the equality ${lim}_{{r}_{s},{v}_{\mathrm{\text{out}}}\to 0}{\stackrel{\u0307}{S}}_{\mathrm{\text{tot}}}/ln\gamma ={v}_{\mathrm{\text{cell}}}{j}_{\mathrm{\text{in}}}/{l}_{\mathrm{\text{cell}}}$ (Equation 33). Comparison of this with the expression of $C$ (Equation 5) enables us to relate ${\alpha}_{o}$ with the free energy cost, $Tln\left(\gamma \right)$. Thus, the thermodynamic cost associated with the morphogen profile is expressed as $C\equiv \left({\alpha}_{o}{v}_{\mathrm{\text{cell}}}{j}_{\mathrm{\text{in}}}\mathrm{\text{/}}{l}_{\mathrm{\text{cell}}}\right)\tau $ (Equation 5) where $\tau $ is the time required to generate the profile and ${\alpha}_{o}{v}_{\mathrm{\text{cell}}}{j}_{\mathrm{\text{in}}}\mathrm{\text{/}}{l}_{\mathrm{\text{cell}}}$ approximates the rate of total dissipation, $T{\stackrel{\u0307}{S}}_{\mathrm{\text{tot}}}$. In the revised manuscript, we direct the readership to the Appendix, where this discussion has been added.
7) All derivations are mode for the finite length L. However, as argued in the specific calculations L>>λ, so why not to derive the results for L>>1. In this case, all the formulas are much simpler and there is a clear physical interpretation for all the quantities. Why the bulkier and less transparent expressions are used in the text? If authors wanted to explain things better to a broader readership and to emphasize the physics they better use a simple analytical framework. Or they need to explain why they used a finite L.
The revised manuscript includes the more transparent expressions at the limit of large $L$.
8) It would be better to describe fully how the characteristic lengths were estimated. This is because the reported numbers do not fully agree with what is reported in experiments. For example, for BCD Berezhkovskii and Shvartsman report λ=75 nm, not 100nm as given in this paper. Why are such different results are reported?
In the plots, we used the characteristic lengths reported from each experimental study. In the case of the Bcd profile, there have been slight differences depending on the reference that reports the characteristic length. However, relatively small discrepancies of the characteristic lengths do not change the qualitative conclusions of our work. For more clarity, we added guiding lines in Figure 2S1 with the experimental measurements of the morphogen profiles, which show that the approximations of exponentially decaying profiles with the reported characteristic lengths fits the raw data faithfully.
https://doi.org/10.7554/eLife.70034.sa2Article and author information
Author details
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.
Senior Editor
 Naama Barkai, Weizmann Institute of Science, Israel
Reviewing Editor
 Gordon J Berman, Emory University, United States
Reviewers
 Mariela Petkova, Harvard University, United States
 Anatoly Kolomeisky, Rice University, United States
Publication history
 Preprint posted: April 14, 2021 (view preprint)
 Received: May 4, 2021
 Accepted: August 13, 2021
 Accepted Manuscript published: August 17, 2021 (version 1)
 Version of Record published: September 22, 2021 (version 2)
Copyright
© 2021, Song and Hyeon
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.
Metrics

 1,060
 Page views

 168
 Downloads

 3
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Developmental Biology
 Neuroscience
The formation of neural circuits requires extensive interactions of cellsurface proteins to guide axons to their correct target neurons. Transcellular interactions of the adhesion G proteincoupled receptor latrophilin2 (Lphn2) with its partner teneurin3 instruct the precise assembly of hippocampal networks by reciprocal repulsion. Lphn2 acts as a repulsive receptor in distal CA1 neurons to direct their axons to proximal subiculum, and as a repulsive ligand in proximal subiculum to direct proximal CA1 axons to distal subiculum. It remains unclear if Lphn2mediated intracellular signaling is required for its role in either context. Here, we show that Lphn2 couples to Gα_{12/13 }in heterologous cells; this coupling is increased by constitutive exposure of the tethered agonist. Specific mutations of Lphn2's tethered agonist region disrupt its G protein coupling and autoproteolytic cleavage, whereas mutating the autoproteolytic cleavage site alone prevents cleavage but preserves a functional tethered agonist. Using an in vivo misexpression assay, we demonstrate that wildtype Lphn2 misdirects proximal CA1 axons to proximal subiculum and that Lphn2 tethered agonist activity is required for its role as a repulsive receptor in axons. By contrast, neither tethered agonist activity nor autoproteolysis was necessary for Lphn2's role as a repulsive ligand in the subiculum target neurons. Thus, tethered agonist activity is required for Lphn2mediated neural circuit assembly in a contextdependent manner.

 Developmental Biology
 Evolutionary Biology
The gill skeleton of cartilaginous fishes (sharks, skates, rays, and holocephalans) exhibits a striking anterior–posterior polarity, with a series of fine appendages called branchial rays projecting from the posterior margin of the gill arch cartilages. We previously demonstrated in the skate (Leucoraja erinacea) that branchial rays derive from a posterior domain of pharyngeal arch mesenchyme that is responsive to Sonic hedgehog (Shh) signaling from a distal gill arch epithelial ridge (GAER) signaling centre. However, how branchial ray progenitors are specified exclusively within posterior gill arch mesenchyme is not known. Here, we show that genes encoding several Wnt ligands are expressed in the ectoderm immediately adjacent to the skate GAER, and that these Wnt signals are transduced largely in the anterior arch environment. Using pharmacological manipulation, we show that inhibition of Wnt signalling results in an anterior expansion of Shh signal transduction in developing skate gill arches, and in the formation of ectopic anterior branchial ray cartilages. Our findings demonstrate that ectodermal Wnt signalling contributes to gill arch skeletal polarity in skate by restricting Shh signal transduction and chondrogenesis to the posterior arch environment and highlights the importance of signalling interactions at embryonic tissue boundaries for cell fate determination in vertebrate pharyngeal arches.