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

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: September 22, 2021 (This version)
Accepted Manuscript published: August 17, 2021 (Go to version)
Accepted: August 13, 2021
Received: May 4, 2021
Preprint posted: April 14, 2021 (Go to version)

1. Of interest
Dynamic compartmentalization of the pro-invasive transcription factor NHR-67 reveals a role for Groucho in regulating a proliferative-invasive cellular switch in C. elegans

Taylor N Medwig-Kinney, Brian A Kinney ... David Q Matus

Research Article Dec 1, 2023
Further reading

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

Abstract

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

Introduction

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

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

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

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

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

Figure 1

Download asset Open asset

Positional information transfer by the morphogen gradient.

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

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

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

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

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

Results

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

Point measurement

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

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

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

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

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

Figure 2 with 1 supplement see all

Download asset Open asset

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Space-time-averaged measurement

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Figure 3

Download asset Open asset

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

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

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

Discussion

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

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

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

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

The entropic cost of forming precise morphogen profiles

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

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

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

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

Cost-effectiveness of precise morphogen profile formation

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

Figure 4

Download asset Open asset

Optimal concentration profile of morphogens.

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

Concluding remarks

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

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

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

Materials and methods

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

Appendix 1

Reversible reaction-diffusion model of morphogen dynamics

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

Appendix 1—figure 1

Download asset Open asset

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

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

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

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

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

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

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

Speed

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

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

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

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

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

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

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

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

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

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

with

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

and

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

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

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

Precision

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

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

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

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

Thermodynamic cost

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Trade-off product

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Then, the trade-off product reads

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

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

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

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

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

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

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

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

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

and

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

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

Appendix 2

Distributed morphogen source

Appendix 2—figure 1

Download asset Open asset

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

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

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

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

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

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

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

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

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

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

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

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

τ (x) = 𝒩 / ℳ,

with

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

and

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

Then, the trade-off product is expressed as

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

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

Dynamics on a sphere

Appendix 2—figure 2

Download asset Open asset

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

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

Appendix 2—figure 3

Download asset Open asset

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

where

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

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

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

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

Appendix 3

Length scales of Bcd, Wg, Hh, and Dpp

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

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

Length scales of Cyclops, Squint, and Fgf8

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

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

Appendix 3—table 1

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

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

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

Time scales of Bcd, Wg, Hh, and Dpp

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

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

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

Data availability

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

References

(2014) Fly wing vein patterns have spatial reproducibility of a single cell
Journal of the Royal Society Interface 11:20140443.
https://doi.org/10.1098/rsif.2014.0443
- PubMed
- Google Scholar
(2018) Scale-invariant patterning by size-dependent inhibition of nodal signalling
Nature Cell Biology 20:1032–1042.
https://doi.org/10.1038/s41556-018-0155-7
- PubMed
- Google Scholar
(2020) The WG and DPP morphogens regulate gene expression by modulating the frequency of transcriptional bursts
eLife 9:e56076.
https://doi.org/10.7554/eLife.56076
- PubMed
- Google Scholar
1. Barato AC
2. Seifert U
(2015) Thermodynamic uncertainty relation for biomolecular processes
Physical Review Letters 114:1158101.
https://doi.org/10.1103/PhysRevLett.114.158101
- Google Scholar
(2010) How long does it take to establish a morphogen gradient?
Biophysical Journal 99:L59–L61.
https://doi.org/10.1016/j.bpj.2010.07.045
- PubMed
- Google Scholar
(2011) Formation of morphogen gradients: local accumulation time
Physical Review E 83:051906.
https://doi.org/10.1103/PhysRevE.83.051906
- Google Scholar
1. Berg HC
2. Purcell EM
(1977) Physics of chemoreception
Biophysical Journal 20:193–219.
https://doi.org/10.1016/S0006-3495(77)85544-6
- PubMed
- Google Scholar
(2008) Precision of the Dpp gradient
Development 135:1137–1146.
https://doi.org/10.1242/dev.012062
- PubMed
- Google Scholar
1. Briscoe J
2. Small S
(2015) Morphogen rules: design principles of gradient-mediated embryo patterning
Development 142:3996–4009.
https://doi.org/10.1242/dev.129452
- PubMed
- Google Scholar
1. Chaudhary V
2. Hingole S
3. Frei J
4. Port F
5. Strutt D
6. Boutros M
(2019) Robust wnt signaling is maintained by a wg protein gradient and Fz2 receptor activity in the developing Drosophila wing
Development 146:dev174789.
https://doi.org/10.1242/dev.174789
- PubMed
- Google Scholar
1. Chen W
2. Huang H
3. Hatori R
4. Kornberg TB
(2017) Essential basal cytonemes take up hedgehog in the Drosophila wing imaginal disc
Development 144:3134–3144.
https://doi.org/10.1242/dev.149856
- PubMed
- Google Scholar
1. Crick F
(1970) Diffusion in embryogenesis
Nature 225:420–422.
https://doi.org/10.1038/225420a0
- Google Scholar
1. Dechant A
2. Sasa SI
(2020) Fluctuation-response inequality out of equilibrium
PNAS 117:6430–6436.
https://doi.org/10.1073/pnas.1918386117
- PubMed
- Google Scholar
(2020) A mechanism for hunchback promoters to readout morphogenetic positional information in less than a minute
eLife 9:e49758.
https://doi.org/10.7554/eLife.49758
- PubMed
- Google Scholar
1. Driever W
2. Nüsslein-Volhard C
(1988a) A gradient of bicoid protein in Drosophila embryos
Cell 54:83–93.
https://doi.org/10.1016/0092-8674(88)90182-1
- PubMed
- Google Scholar
1. Driever W
2. Nüsslein-Volhard C
(1988b) The bicoid protein determines position in the Drosophila embryo in a concentration-dependent manner
Cell 54:95–104.
https://doi.org/10.1016/0092-8674(88)90183-3
- PubMed
- Google Scholar
1. Drocco JA
2. Grimm O
3. Tank DW
4. Wieschaus E
(2011) Measurement and perturbation of morphogen lifetime: effects on gradient shape
Biophysical Journal 101:1807–1815.
https://doi.org/10.1016/j.bpj.2011.07.025
- PubMed
- Google Scholar
(2012) The synthesis-diffusion-degradation model explains bicoid gradient formation in unfertilized eggs
Physical Biology 9:055004.
https://doi.org/10.1088/1478-3975/9/5/055004
- PubMed
- Google Scholar
(2015) Response to nodal morphogen gradient is determined by the kinetics of target gene induction
eLife 4:e05042.
https://doi.org/10.7554/eLife.05042
- Google Scholar
1. Emberly E
(2008) Optimizing the readout of morphogen gradients
Physical Review E 77:041903.
https://doi.org/10.1103/PhysRevE.77.041903
- Google Scholar
(2000) Gradient formation of the TGF-beta homolog Dpp
Cell 103:981–992.
https://doi.org/10.1016/S0092-8674(00)00200-2
- PubMed
- Google Scholar
(2009) Role of spatial averaging in the precision of gene expression patterns
Physical Review Letters 103:258101.
https://doi.org/10.1103/PhysRevLett.103.258101
- PubMed
- Google Scholar
(2018) Information thermodynamics of turing patterns
Physical Review Letters 121:108301.
https://doi.org/10.1103/PhysRevLett.121.108301
- PubMed
- Google Scholar
1. Fancher S
2. Mugler A
(2020) Diffusion vs. direct transport in the precision of morphogen readout
eLife 9:e58981.
https://doi.org/10.7554/eLife.58981
- PubMed
- Google Scholar
1. Foe VE
2. Alberts BM
(1983) Studies of nuclear and cytoplasmic behaviour during the five mitotic cycles that precede gastrulation in Drosophila embryogenesis
Journal of Cell Science 61:31–70.
https://doi.org/10.1242/jcs.61.1.31
- PubMed
- Google Scholar
(2016) Dissipation bounds all Steady-State current fluctuations
Physical Review Letters 116:120601.
https://doi.org/10.1103/PhysRevLett.116.120601
- PubMed
- Google Scholar
(2007) Probing the limits to positional information
Cell 130:153–164.
https://doi.org/10.1016/j.cell.2007.05.025
- PubMed
- Google Scholar
(2017) Concentration dependent chromatin states induced by the bicoid morphogen gradient
eLife 6:e28275.
https://doi.org/10.7554/eLife.28275
- PubMed
- Google Scholar
1. He F
2. Saunders TE
3. Wen Y
4. Cheung D
5. Jiao R
6. ten Wolde PR
7. Howard M
8. Ma J
(2010) Shaping a morphogen gradient for positional precision
Biophysical Journal 99:697–707.
https://doi.org/10.1016/j.bpj.2010.04.073
- PubMed
- Google Scholar
1. Heuett WJ
2. Qian H
(2006) Grand canonical markov model: a stochastic theory for open nonequilibrium biochemical networks
The Journal of Chemical Physics 124:044110.
https://doi.org/10.1063/1.2165193
- PubMed
- Google Scholar
1. Horowitz JM
2. Gingrich TR
(2017) Proof of the finite-time thermodynamic uncertainty relation for steady-state currents
Physical Review E 96:020103.
https://doi.org/10.1103/PhysRevE.96.020103
- PubMed
- Google Scholar
1. Horowitz JM
2. Gingrich TR
(2020) Thermodynamic uncertainty relations constrain non-equilibrium fluctuations
Nature Physics 16:15–20.
https://doi.org/10.1038/s41567-019-0702-6
- Google Scholar
(2002) Establishment of developmental precision and proportions in the early Drosophila embryo
Nature 415:798–802.
https://doi.org/10.1038/415798a
- PubMed
- Google Scholar
1. Huang H
2. Kornberg TB
(2015) Myoblast cytonemes mediate wg signaling from the wing imaginal disc and Delta-Notch signaling to the air sac primordium
eLife 4:e06114.
https://doi.org/10.7554/eLife.06114
- PubMed
- Google Scholar
1. Hwang W
2. Hyeon C
(2018) Energetic costs, precision, and transport efficiency of molecular motors
The Journal of Physical Chemistry Letters 9:513–520.
https://doi.org/10.1021/acs.jpclett.7b03197
- PubMed
- Google Scholar
1. Hyeon C
2. Hwang W
(2017) Physical insight into the thermodynamic uncertainty relation using brownian motion in tilted periodic potentials
Physical Review E 96:012156.
https://doi.org/10.1103/PhysRevE.96.012156
- PubMed
- Google Scholar
1. Ilker E
2. Hinczewski M
(2019) Modeling the growth of organisms validates a general relation between metabolic costs and natural selection
Physical Review Letters 122:2381011.
https://doi.org/10.1103/PhysRevLett.122.238101
- Google Scholar
(2006) Senseless and daughterless confer neuronal identity to epithelial cells in the Drosophila wing margin
Development 133:1683–1692.
https://doi.org/10.1242/dev.02338
- PubMed
- Google Scholar
(2009) Dynamics of the dorsal morphogen gradient
PNAS 106:21707–21712.
https://doi.org/10.1073/pnas.0912395106
- PubMed
- Google Scholar
(2007) Kinetics of morphogen gradient formation
Science 315:521–525.
https://doi.org/10.1126/science.1135774
- PubMed
- Google Scholar
1. Kolomeisky AB
(2011) Formation of a morphogen gradient: acceleration by degradation
The Journal of Physical Chemistry Letters 2:1502–1505.
https://doi.org/10.1021/jz2004914
- Google Scholar
1. Li P
2. Markson JS
3. Wang S
4. Chen S
5. Vachharajani V
6. Elowitz MB
(2018) Morphogen gradient reconstitution reveals hedgehog pathway design principles
Science 360:543–548.
https://doi.org/10.1126/science.aao0645
- PubMed
- Google Scholar
(2013) Dynamic interpretation of maternal inputs by the Drosophila segmentation gene network
PNAS 110:6724–6729.
https://doi.org/10.1073/pnas.1220912110
- PubMed
- Google Scholar
1. Lo WC
2. Zhou S
3. Wan FY
4. Lander AD
5. Nie Q
(2015) Robust and precise morphogen-mediated patterning: trade-offs, constraints and mechanisms
Journal of the Royal Society, Interface 12:20141041.
https://doi.org/10.1098/rsif.2014.1041
- PubMed
- Google Scholar
1. Lucas T
2. Tran H
3. Perez Romero CA
4. Guillou A
5. Fradin C
6. Coppey M
7. Walczak AM
8. Dostatni N
(2018) 3 minutes to precisely measure morphogen concentration
PLOS Genetics 14:e1007676.
https://doi.org/10.1371/journal.pgen.1007676
- PubMed
- Google Scholar
1. Lynch M
2. Marinov GK
(2015) The bioenergetic costs of a gene
PNAS 112:15690–15695.
https://doi.org/10.1073/pnas.1514974112
- PubMed
- Google Scholar
(2009) Cell competition, growth and size control in the Drosophila wing imaginal disc
Development 136:3747–3756.
https://doi.org/10.1242/dev.038406
- PubMed
- Google Scholar
1. Morishita Y
2. Iwasa Y
(2009) Accuracy of positional information provided by multiple morphogen gradients with correlated noise
Physical Review E 79:061905.
https://doi.org/10.1103/PhysRevE.79.061905
- Google Scholar
1. Müller P
2. Rogers KW
3. Jordan BM
4. Lee JS
5. Robson D
6. Ramanathan S
7. Schier AF
(2012) Differential diffusivity of nodal and lefty underlies a reaction-diffusion patterning system
Science 336:721–724.
https://doi.org/10.1126/science.1221920
- PubMed
- Google Scholar
1. Müller P
2. Rogers KW
3. Yu SR
4. Brand M
5. Schier AF
(2013) Morphogen transport
Development 140:1621–1638.
https://doi.org/10.1242/dev.083519
- PubMed
- Google Scholar
1. Nahmad M
2. Stathopoulos A
(2009) Dynamic interpretation of hedgehog signaling in the Drosophila wing disc
PLOS Biology 7:e1000202.
https://doi.org/10.1371/journal.pbio.1000202
- PubMed
- Google Scholar
1. Nellen D
2. Burke R
3. Struhl G
4. Basler K
(1996) Direct and long-range action of a DPP morphogen gradient
Cell 85:357–368.
https://doi.org/10.1016/S0092-8674(00)81114-9
- PubMed
- Google Scholar
1. Neumann CJ
2. Cohen SM
(1997) Long-range action of wingless organizes the dorsal-ventral Axis of the Drosophila wing
Development 124:871–880.
https://doi.org/10.1242/dev.124.4.871
- PubMed
- Google Scholar
1. Nowak M
2. Machate A
3. Yu SR
4. Gupta M
5. Brand M
(2011) Interpretation of the FGF8 morphogen gradient is regulated by endocytic trafficking
Nature Cell Biology 13:153–158.
https://doi.org/10.1038/ncb2155
- PubMed
- Google Scholar
1. Oginuma M
2. Moncuquet P
3. Xiong F
4. Karoly E
5. Chal J
6. Guevorkian K
7. Pourquié O
(2017) A gradient of glycolytic activity coordinates FGF and wnt signaling during elongation of the body Axis in amniote embryos
Developmental Cell 40:342–353.
https://doi.org/10.1016/j.devcel.2017.02.001
- PubMed
- Google Scholar
1. Oginuma M
2. Harima Y
3. Tarazona OA
4. Diaz-Cuadros M
5. Michaut A
6. Ishitani T
7. Xiong F
8. Pourquié O
(2020) Intracellular pH controls WNT downstream of glycolysis in amniote embryos
Nature 584:98–101.
https://doi.org/10.1038/s41586-020-2428-0
- PubMed
- Google Scholar
1. Perry MW
2. Bothma JP
3. Luu RD
4. Levine M
(2012) Precision of hunchback expression in the Drosophila embryo
Current Biology 22:2247–2252.
https://doi.org/10.1016/j.cub.2012.09.051
- PubMed
- Google Scholar
(2019) Optimal decoding of cellular identities in a genetic network
Cell 176:844–855.
https://doi.org/10.1016/j.cell.2019.01.007
- PubMed
- Google Scholar
(2017) Finite-time generalization of the thermodynamic uncertainty relation
Physical Review E 96:012101.
https://doi.org/10.1103/PhysRevE.96.012101
- PubMed
- Google Scholar
(2019) Heat oscillations driven by the embryonic cell cycle reveal the energetic costs of signaling
Developmental Cell 48:646–658.
https://doi.org/10.1016/j.devcel.2018.12.024
- PubMed
- Google Scholar
(2020) Contribution of increasing plasma membrane to the energetic cost of early zebrafish embryogenesis
Molecular Biology of the Cell 31:520–526.
https://doi.org/10.1091/mbc.E19-09-0529
- PubMed
- Google Scholar
(2020) Optogenetic investigation of BMP target gene expression diversity
eLife 9:e58641.
https://doi.org/10.7554/eLife.58641
- PubMed
- Google Scholar
1. Rosenbauer J
2. Zhang C
3. Mattes B
4. Reinartz I
5. Wedgwood K
6. Schindler S
7. Sinner C
8. Scholpp S
9. Schug A
(2020) Modeling of Wnt-mediated tissue patterning in vertebrate embryogenesis
PLOS Computational Biology 16:e1007417.
https://doi.org/10.1371/journal.pcbi.1007417
- PubMed
- Google Scholar
1. Saunders TE
2. Howard M
(2009) Morphogen profiles can be optimized to buffer against noise
Physical Review E 80:041902.
https://doi.org/10.1103/PhysRevE.80.041902
- PubMed
- Google Scholar
1. Shvartsman SY
2. Baker RE
(2012) Mathematical models of morphogen gradients and their effects on gene expression
Wiley Interdisciplinary Reviews: Developmental Biology 1:715–730.
https://doi.org/10.1002/wdev.55
- PubMed
- Google Scholar
1. Song Y
2. Park JO
3. Tanner L
4. Nagano Y
5. Rabinowitz JD
6. Shvartsman SY
(2019) Energy budget of Drosophila embryogenesis
Current Biology 29:R566–R567.
https://doi.org/10.1016/j.cub.2019.05.025
- PubMed
- Google Scholar
1. Song Y
2. Hyeon C
(2020) Thermodynamic cost, speed, fluctuations, and error reduction of biological copy machines
The Journal of Physical Chemistry Letters 11:3136–3143.
https://doi.org/10.1021/acs.jpclett.0c00545
- PubMed
- Google Scholar
1. Song Y
2. Hyeon C
(2021) Thermodynamic uncertainty relation to assess biological processes
The Journal of Chemical Physics 154:130901.
https://doi.org/10.1063/5.0043671
- PubMed
- Google Scholar
1. Strigini M
2. Cohen SM
(1997) A hedgehog activity gradient contributes to AP axial patterning of the Drosophila wing
Development 124:4697–4705.
https://doi.org/10.1242/dev.124.22.4697
- PubMed
- Google Scholar
(1989) The gradient morphogen bicoid is a concentration-dependent transcriptional activator
Cell 57:1259–1273.
https://doi.org/10.1016/0092-8674(89)90062-7
- PubMed
- Google Scholar
1. Tabata T
2. Kornberg TB
(1994) Hedgehog is a signaling protein with a key role in patterning Drosophila imaginal discs
Cell 76:89–102.
https://doi.org/10.1016/0092-8674(94)90175-9
- PubMed
- Google Scholar
1. Teimouri H
2. Kolomeisky AB
(2014) Development of morphogen gradient: the role of dimension and discreteness
The Journal of Chemical Physics 140:085102.
https://doi.org/10.1063/1.4866453
- Google Scholar
1. Teimouri H
2. Kolomeisky AB
(2016) New model for understanding mechanisms of biological signaling: direct transport via cytonemes
The Journal of Physical Chemistry Letters 7:180–185.
https://doi.org/10.1021/acs.jpclett.5b02703
- PubMed
- Google Scholar
(2015) Positional information, positional error, and readout precision in morphogenesis: a mathematical framework
Genetics 199:39–59.
https://doi.org/10.1534/genetics.114.171850
- PubMed
- Google Scholar
1. Tkačik G
2. Gregor T
(2021) The many bits of positional information
Development 148:dev176065.
https://doi.org/10.1242/dev.176065
- PubMed
- Google Scholar
(2004) Patched controls the hedgehog gradient by endocytosis in a dynamin-dependent manner, but this internalization does not play a major role in signal transduction
Development 131:2395–2408.
https://doi.org/10.1242/dev.01102
- PubMed
- Google Scholar
(2007) Fundamental limits to position determination by concentration gradients
PLOS Computational Biology 3:e78.
https://doi.org/10.1371/journal.pcbi.0030078
- PubMed
- Google Scholar
1. Turing AM
(1990) The chemical basis of morphogenesis. 1953
Bulletin of Mathematical Biology 52:153–197.
https://doi.org/10.1016/S0092-8240(05)80008-4
- PubMed
- Google Scholar
(2011) Dynamics of Dpp signaling and proliferation control
Science 331:1154–1159.
https://doi.org/10.1126/science.1200037
- PubMed
- Google Scholar
1. Wolpert L
(1969) Positional information and the spatial pattern of cellular differentiation
Journal of Theoretical Biology 25:1–47.
https://doi.org/10.1016/S0022-5193(69)80016-0
- Google Scholar
1. Yu SR
2. Burkhardt M
3. Nowak M
4. Ries J
5. Petrásek Z
6. Scholpp S
7. Schwille P
8. Brand M
(2009) Fgf8 morphogen gradient forms by a source-sink mechanism with freely diffusing molecules
Nature 461:533–536.
https://doi.org/10.1038/nature08391
- PubMed
- Google Scholar
1. Zagorski M
2. Tabata Y
3. Brandenberg N
4. Lutolf MP
5. Tkačik G
6. Bollenbach T
7. Briscoe J
8. Kicheva A
(2017) Decoding of position in the developing neural tube from antiparallel morphogen gradients
Science 356:1379–1383.
https://doi.org/10.1126/science.aam5887
- PubMed
- Google Scholar
(1996) Direct and long-range action of a wingless morphogen gradient
Cell 87:833–844.
https://doi.org/10.1016/S0092-8674(00)81991-1
- PubMed
- Google Scholar
1. Zhou S
2. Lo WC
3. Suhalim JL
4. Digman MA
5. Gratton E
6. Nie Q
7. Lander AD
(2012) Free extracellular diffusion creates the Dpp morphogen gradient of the Drosophila wing disc
Current Biology 22:668–675.
https://doi.org/10.1016/j.cub.2012.02.065
- PubMed
- Google Scholar

Decision letter

Gordon J Berman

Reviewing Editor; Emory University, United States
Naama Barkai

Senior Editor; Weizmann Institute of Science, Israel
Mariela Petkova

Reviewer; Harvard University, United States
Anatoly Kolomeisky

Reviewer; 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 cost-precision 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 "Cost-precision trade-off 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 trade-off function was presented without supporting arguments that either reflect biological expectations or demonstrate that the results are robust to other trade-off constructs. Importantly, the paper argues that the agreement between the optimal gradient shape predicted by the trade-off 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 cost-limiting 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 trade-off function that they analyze. My question is how robust are the obtained results? The choice of the trade-off 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 trade-off 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 trade-off 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 biologically-oriented 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 space-time-averaged 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 y-axis 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 α₀, assuming the irreversible limit where ln(γ) -> infinity. So it seems that the entropy-based 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 non-equilibrium 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 system-specific quantity. The result is reminiscent of, but qualitatively different from, the so-called thermodynamic uncertainty principle that has generated a lot of excitement in the non-equilibrium 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 cost-effectiveness 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 physics-inspired theory helping us interpret well-studied 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 trade-off 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 trade-offs.

Strengths:

The trade-off 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 gradient-mediated 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 trade-off 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 trade-off 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 trade-off function and evaluate it against measurements in unperturbed systems. The choice of trade-off function was presented without supporting arguments that either reflect biological expectations or demonstrate that the results are robust to other trade-off constructs. Importantly, the paper argues that the agreement between the optimal gradient shape predicted by the trade-off 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 cost-limiting 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.6724-6729

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 cost-precision trade-off relation is not critically evaluated, and (iii) some references to the existing literature are missing.

https://doi.org/10.7554/eLife.70034.sa1

Author 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 trade-off function was presented without supporting arguments that either reflect biological expectations or demonstrate that the results are robust to other trade-off constructs. Importantly, the paper argues that the agreement between the optimal gradient shape predicted by the trade-off 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 cost-limiting 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 co-regulating 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 cost-precision trade-off of individual morphogen profiles, which allows us to show that the morphogen profiles of fruit fly development are nearly optimal, as opposed to those of the zebrafish embryo. It would be fruitful in future work to assess the cost-precision trade-off 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 trade-off 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, $α_{o}$ , in the expression of the cost $C \equiv α_{o} (v_{cell} j_{in} / l_{cell}) τ$ (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 $α_{o} \approx 494 \times 4 \approx 2 \times 10^{3}$ ATPs, where we assume that each peptide-bond formation requires 4 ATPs (Lynch et al., 2015). During the first ~2 hours of Drosophila embryogenesis when the anterior-posterior 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¹² ATPs, which is only a small fraction of the total energy budget of Drosophila embryogenesis ~ $7 \times 10^{16}$ 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 biomass-producing 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 Δ S (t) \frac{⟨ δ X {(t)}^{2} ⟩}{{⟨ X (t) ⟩}^{2}}$ quantifies the trade-offs among the thermodynamic cost ( $Δ S (t)$ ), speed ( ${⟨ X (t) ⟩}^{2}$ ), and the fluctuations ( $⟨ δ X {(t)}^{2} ⟩$ ) of the time-dependent observable $X (t)$ , which can be defined as the steps taken by the molecular motor or the number of peptide-bonds 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 7 - 15 k_{B}$ for the molecular motors, and $Q \approx 45 - 50 k_{B}$ for E. coli ribosome), which suggests that those molecular processes have evolved to operate cost-effectively (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 cost-precision trade-off product ( $π$ ), 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.4 L$ , $π_{o} (x_{b}, λ_{B cd}) \approx α_{o} (L / l_{cell}) 0.6$ (Figure 2b) where $L / l_{cell} \approx 50,$ and the experimentally reported precision is $ϵ^{2} (x_{b}) \approx 7 \times 10^{- 4}$ (Gregor et al., 2007). Thus, the cost associated with the Bcd profile is $C (x_{b}, λ_{B cd}) = π_{o} (x_{b}, λ_{B cd}) / ϵ^{2} (x_{b}, λ_{B cd}) \approx 40000 α_{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, $τ$ . 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 $τ$ , it is possible to adjust the proportionality constant $α_{o}$ so that $C (x_{b})$ represents the overall cost to produce the concentration profile of nuclear Bcd.

– 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 synthesis-diffusion-degradation 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 trade-off function that they analyze. My question is how robust are the obtained results? The choice of the trade-off 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 trade-off 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 trade-off by $π \equiv C \times ϵ^{2}$ was the following:

There are two quantities of interest in the problem of cost-precision trade-off 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, $σ_{ρ} / ⟨ ρ ⟩$ . 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_{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 trade-off between the cost and precision. The simplest and generic construction for discussing the trade-off relationship between the two quantities, independently from $n_{m}$ (i.e. the system size) is to consider the product between the cost ( $C$ ) and the squared value of $σ_{ρ} / ⟨ ρ ⟩$ .

Let us emphasize again that $π = C \times ϵ^{2}$ , which we called the trade-off product, is a dimensionless quantity, independent of the number of morphogens being produced ( $n_{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 cost-precision trade-off was motivated by the thermodynamic uncertainty relation (TUR), which defines the product between the squared relative error of the time integrated current-like observable $X (t)$ and the total entropy production for $t$ , $ΔS (t)$ , namely, $Q \equiv Δ S (t) \frac{⟨ δ X {(t)}^{2} ⟩}{{⟨ X (t) ⟩}^{2}}$ . In the case of TUR, the two quantities of interest are a function of time t, such that $Δ S (t) \sim t$ , $\frac{⟨ δ X {(t)}^{2} ⟩}{{⟨ X (t) ⟩}^{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 \geq 2 k_{B}$ , where $k_{B}$ is the Boltzmann constant, and the cost-effectiveness 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 space-time averaging, we addressed the cost-precision trade-off 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, $λ_{min} \approx 0.43 x_{b}$ for the point measurement and $λ_{min} \approx 0.5 x_{b}$ for the space-time-averaged measurement, both of which are the outcomes of minimizing the respective trade-off 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 trade-off 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 cost-effectively 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 2-S1). In the case of Hh, we included the full range of spatial boundaries formed by the three target genes, dpp, col, and en (Figure 2-S1). 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 biologically-oriented 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 space-time-averaged 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 position-dependent expression profile. Additionally, cells with many receptors positioned across the cell surface may perform space-time-averaging by integrating the signal from morphogen-bound receptors for an extended time interval. Theoretically, cells may control the duration of the measurement by modulating the availability of morphogen-sensing 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 space-time 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 y-axis 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 $λ$ that minimizes the trade-off product at a given target boundary. In Figure 3c of the revised manuscript, we have overlaid the linear approximation of the optimal characteristic length $λ_{min}$ from the point measurement model with $λ_{min}$ from the model with space-time-averaging, so that the readership can compare the two results side-by-side. Additionally, in the discussion, we clarified that $λ_{min}$ of the two models were similar, in that $λ_{min} \approx 0.43 L$ in the former and $λ_{min} \approx 0.5 L$ 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 α₀, assuming the irreversible limit where ln(γ) -> infinity. So it seems that the entropy-based 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, $α_{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 ( ${\dot{S}}_{tot}$ , Equation 32) is normalized by $ln (γ)$ , and $v_{out}$ and $r_{s}$ are negligibly small, we obtain the equality ${lim}_{r_{s}, v_{out} \to 0} {\dot{S}}_{tot} / ln γ = v_{cell} j_{in} / l_{cell}$ (Equation 33). Comparison of this with the expression of $C$ (Equation 5) enables us to relate $α_{o}$ with the free energy cost, $T ln (γ)$ . Thus, the thermodynamic cost associated with the morphogen profile is expressed as $C \equiv (α_{o} v_{cell} j_{in} / l_{cell}) τ$ (Equation 5) where $τ$ is the time required to generate the profile and $α_{o} v_{cell} j_{in} / l_{cell}$ approximates the rate of total dissipation, $T {\dot{S}}_{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 2-S1 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.sa2

Article and author information

Author details

Yonghyun Song

Korea Institute for Advanced Study, Seoul, Republic of Korea

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

Competing interests
No competing interests declared

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

Korea Institute for Advanced Study, Seoul, Republic of Korea

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

For correspondence
hyeoncb@kias.re.kr

Competing interests
No competing interests declared

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

Funding

Korea Institute for Advanced Study (CG067102)

Yonghyun Song

Korea Institute for Advanced Study (CG035003)

Changbong Hyeon

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

Acknowledgements

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

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

Version 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

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.