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

  1. Yonghyun Song
  2. Changbong Hyeon  Is a corresponding author
  1. Korea Institute for Advanced Study, Republic of Korea

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.

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 xb, σx2(xb), is defined as the product between the variance of the morphogen concentration, σρ2(xb), and the squared inverse slope of the morphogen profile, (xρ(x))x=xb-2.

We begin by defining a target boundary, x=xb, where the morphogen concentration is at its critical threshold value ρbρ(xb). 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 ρ^b, the nucleus can estimate the location of target boundary, such that x^b=x^b(ρ^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=xb with respect to its true value ρb is given by σρ2(xb)[(ρ^b-ρb)2=(1/N)i=1N(ρ^b,i-ρb)2]. In other words, σρ2(xb) represents the inherent variability of the morphogen profile, which can be determined experimentally by performing multiple measurements, ρ^b,i (i=1,2,N). For the nucleus to make a single measurement, ρ^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, ρ(x^b)ρ(xb)+(ρ(x^b)/x^b)|x^b=xb(x^b-xb)+, allows one to relate the error in measured concentration with the positional error in estimating the target boundary σx2(xb)[(x^b-xb)2=(1/N)i=1N(x^b,i-xb)2] as follows:

(1) σx2(xb)(ρ(x^b)x^b)x^b=xb-2σρ2(xb),

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 0xL, 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 (vcell) and length (lcell) at the interval between x and x+lcell by ρ(x,t)[conc#/vcell]. 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 jin[conc×lcell/time]. At all positions, the morphogen molecules are depleted with rate kd [time-1], while spreading across the cells with the diffusivity D[lcell2/time] (Figure 2a). For L/lcell1 (i.e. at the continuum limit), the spatiotemporal dynamics of the morphogen is described using the reaction-diffusion equation

(2) tρ(x,t)=Dx2ρ(x,t)-kdρ(x,t),

with boundary conditions -Dxρ(x,t)|x=0=jin and Dxρ(x,t)|x=L=0. Then, the concentration profile at steady state is obtained as

(3) ρss(x)=jinDkdcosh(Lxλ)sinh(Lλ)jinDkdex/λ,

where λ=(D/kd) is the characteristic decay length determined by the diffusivity of morphogen (D) and the depletion rate (kd). 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
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(xb), 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 xb/L. The color scale indicates the trade-off product πo(xb) computed for each pair of λ/L and xb/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 αoL/lcell. 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)(ρ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,

(4) τ(x)=0R(x,t)𝑑t=kd-1f(x;λ),

where f(x;λ)12[1+Lλcoth(Lλ)-L-xλtanh(L-xλ)] is a dimensionless quantity determined by x and λ. Then, vcelljin/lcell×τ(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

(5) 𝒞(x)=αo(vcelljinlcellkd)f(x;λ)αo(vcelljinlcellkd)12(1+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 αo494×42×103 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 × 108 Bcd molecules are produced (Drocco et al., 2012). Thus, the cost of generating the Bcd profile is on the order of 1012 ATPs, which is a small fraction of the total energy budget of Drosophila embryogenesis (∼7 × 1016 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):

(6) ϵ2(x)σx2(x)x2=(logρss(x)logx)2σρ2(x)ρ^(x)2=(λ3x2lcell)(vcelljinlcellkd)1sinh(Lλ)cosh(Lxλ)[sinh(Lxλ)]2(λ3x2lcell)(vcelljinlcellkd)1ex/λ,

where ρ^(x) represents the morphogen concentration estimated at x from a measurement. For large system size (Lλ), ϵ2(x)(λ3/x2)ex/λ. Note that the multiple measurements of ρ^(x) yield the mean concentration profile ρ^(x)=(1/N)i=1Nρ^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. vcell2σρ2(x)=vcellρ^(x)=vcellρ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(xb)7×10-4 for Bcd to ϵ2(xb)1.5×10-2 for Dpp, where xb(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 nm is the number of morphogen molecules produced for a certain time duration then 𝒞(x)nm and σρ2(x)/ρ^(x)21/nm, the latter of which simply arises from the central limit theorem, or can be rationalized based on the Berg-Purcell result (δρ^/ρ^1/Daρ^τ) (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, nm corresponds to (vcelljin/lcellkd) (the morphogen molecules produced for the time duration kd-1). The nm-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,

(7) πo(x;λ)𝒞(x;λ)×ϵ2(x;λ)π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=xb, 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 xb, 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(xb;λ) is close to its lower bound πo,min(xb;λ), the system is cost-effective at generating precise morphogen profiles at xb. Generally, πo,min(xb;λ) increases monotonically with xb, 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 xb, 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=xb, that is, πo,min(xb), such that

(8) λmin(xb)=argminλπo(xb;λ).

The optimal decay length, λmin(xb), also increases monotonically with the target boundary xb; for xb<L/2, λmin is well approximated by the linear relationship λmin(xb)0.43xb (See Equation 40 for the derivation).

It is of great interest to compare xb and λ of real morphogen profiles against the optimal λmin(xb) that one can predict from Equation 8. For the morphogen profiles of Bcd, Wg, Hh, and Dpp, we first estimated the possible range of xb 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 λBcd=100μm, λWg=6μm, λHh=8μm, and λDpp=20μm (Houchmandzadeh et al., 2002; Kicheva et al., 2007; Wartlick et al., 2011). Figure 2b shows the pairs of xb 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(xb). 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(xb;λ), is minimized to its lower bound, πo,min(xb;λ) (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 xb0.4L, πo(xb;λBcd)αo(L/lcell)0.6 (Figure 2b) where L/lcell50, and the experimentally reported precision is ϵ2(xb)7×10-4 (Gregor et al., 2007). Thus, the cost associated with the Bcd profile is C(xb;λBcd)=πo(xb;λBcd)/ϵ2(xb;λBcd)4×104αo, or equivalently, the thermodynamic cost of synthesizing and degrading 4 × 104 Bcd molecules. This is orders of magnitude smaller than the earlier estimate of ∼5 × 108 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(xb) 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

(9) mi(x)=1T0T𝑑tx-ax+aρ^i(y,t)𝑑y,

where ρ^i(y,t) represents an estimate of molecular count at position y from an i-th independent measurement. Repeated measurements (N1) lead to the mean value of molecular counts

(10) m(x)=1Ni=1Nmi(x)=2jinkde-x/λsinh(a/λ).

In this case, jin and kd are in the units of #/time and time-1, respectively. With a sufficiently long measurement time (Tkd-1), the variance of m(x) can be approximated to Fancher and Mugler, 2020

(11) σm2(x)1Ni=1N(mi(x)-m(x))2=4jine-x/λsinh(a/λ)Tkd2(1-(λ/a)cosh(a/λ)sinh(a/λ)+12ea/λ(λ/a)sinh(a/λ)).

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

(12) ϵT2(x)(logm(x)logx)-2σm2(x)m(x)2=(λ/x)2ex/λjinT(e-a/λ+3e3a/λ-4(1+a/λ)ea/λ2(e2a/λ-1)2).

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

(13) 𝒞T=αojinT.

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, jinT cancels off, which yields the following expression of the cost-precision trade-off:

(14) πT(x;λ)𝒞T×ϵT2(x;λ)πT,min(x;λ).

With a given target boundary x=xb and the sensor size a, the value of πT(xb;λ) is solely determined by λ, analogously to πo(xb;λ). Thus, the lower bound of πT(xb;λ), i.e., πT,min(xb;λ), can be determined simply by tuning λ to a value that minimizes πT(xb;λ), that is, λmin(xb). Both πT,min(xb;λ) and λmin(xb) increase monotonically with xb. In particular, λmin(xb) is found in the range of (xb-a)/2λmin(xb)xb/2 (see Equations 42 and 43 for derivation). The boundary of the inaccessible region in Figure 3b constrains πT(xb;λ), 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(xb)0.1) at xb=60a, the characteristic length must be λ30a, which demands that the system synthesize at least ∼182 morphogen molecules (Figure 3b).

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

(a) Schematic of the model. (b) The optimal trade-off product πT(xb) obtained numerically with respect to the location of the target boundary position normalized by the sensor size (xb/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 xb/a. The color scale indicates the trade-off product πT(xb) computed for each pair of λ/a and xb/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 (xb), 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. Tkd-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 (λ=D/kd) of all four morphogen profiles are close to their respective λmin(xb) 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)/(2a) 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 (jin/Dkd)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, ϵT2(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 ϵT2(x) measured by a single receptor.

The cost associated with the morphogen profile in the case of point measurement, vcelljinlcell-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, jinT, integrated over a time interval Tkd-1, is identical for all positions. The latter assumption is required in the derivation of the expression of ϵT2(x) (see Appendix 1 of Fancher and Mugler, 2020). We additionally demand Tτ(x) in order for jinT to represent the total cost of morphogen profile formation and maintenance. However, T may not necessarily be larger than either kd-1 or τ(x). For instance, the Bcd profile degrades and stabilizes at time scales of 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: λmin0.43xb for the point measurement, and λmin0.5xb 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 (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 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 (ϵrev2(x)). The trade-off among the thermodynamic cost, speed of formation (τrev-1(x)), and precision is quantified by the product of the three quantities,

(15) πo,rev(x)S˙totτrev(x)ϵrev2(x)πo,rev,min(x).

The lower bound, πo,rev,min(xb), increases monotonically with xb (Appendix 1—figure 1).

With S˙totτrev(xb) 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 ρ(xb) at the target boundary x=xb through the combination of morphogen synthesis rate (jin), diffusivity (D) and depletion rate (kd) (inset in Figure 4). A steeper morphogen profile with larger jin and smaller λ(=D/kd) (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 jin 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(0.43-0.50)xb, which suggests that the target boundary of the exponentially decaying morphogen profile, c(x)=c0e-x/λ, should be formed at c(xb)0.1c0. 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(xb)0.1c0 (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.

Optimal concentration profile of morphogens.

(a) Three possible morphogen profiles (i) (red), (ii) (cyan), (iii) (brown) with different λ’s (λ(i)<λ(ii)<λ(iii)), generated with different values of morphogen influxes jin (jin(i)>jin(ii)>jin(iii)). The morphogen concentration of the three profiles coincide at xb, giving rise to the same threshold value ρ(xb) but different positional errors (ϵ2(xb;λ)). (b) The precision of three possible morphogen profiles (i) (red), (ii) (cyan), (iii) (brown) with different λ’s (λ(i)<λ(ii)<λ(iii)), generated with different values of morphogen influxes jin (jin(i)>jin(iii)>jin(ii)). The λ values are identical to those with matching colors in (a), but the red and brown curves are generated with different jin values from those in (a). The cost associated with each morphogen profile are shown in units of αoL/lcell. (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=xb. The trade-off product is minimized to πo=πo(ii)αo(L/lcell)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 (S˙tot). 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
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(xb), obtained numerically by minimizing πo,rev(xb) with respect to λ and γ, with βout. The trade-off product is normalized by L/lcell. (c,d) The optimal characteristic decay length λmin (b) and γmin (c) corresponding to πo,rev,min(xb).

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 (vcell) and size (lcell) at the interval between x and x+lcell by the concentration, ρ(x,t)[conc#/vcell]. At the left boundary (see Appendix 1—figure 1a), the morphogen is injected with rate jin[conc×lcell/time], and taken out by the corresponding reverse reaction with rate vout[lcell/time]. At all positions, the morphogen molecules are depleted with rate kd[time-1] and synthesized by the corresponding reverse reaction with rate rs[conc/time], and spread across space with diffusivity D[lcell2/time]. At the continuum limit (Llcell), the dynamics of the morphogen concentration is described by the partial differential equation

(16) tρ(x,t)=Dx2ρ(x,t)-kdρ(x,t)+rs,

with two boundary conditions at x=0 and x=L, Dxρ(x,t)|x=0=jinvoutρ(0,t) and Dxρ(x,t)|x=L=0. λD/kd is the characteristic length of the dynamics associated with Equation 16. The steady state concentration is expressed as

(17) ρss(x)=rskd(γ-1)cosh[L-xλ]βout-1sinh[Lλ]+cosh[Lλ]=ρss,0(x)+rskd,

where ρss,0(x) denotes the position-dependent portion of the concentration profile. We introduce three dimensionless parameters: βinjin/(λrs), βoutvout/(λkd), and their ratio γβ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)limtρ(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:

(18) ρ^(x,s)=1s[(jin(s+kd)-rsvout)cosh[(L-x)(s+kd)/D](s+kd)sinh[L(s+kd)/D](D(s+kd)/D+voutcoth[L(s+kd)/D])+rs(s+kd)].

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

(19) lims0(sρ^(x,s)-ρ(x,0))=lims00𝑑te-stdρ(x,t)dt=0𝑑tdρ(x,t)dt=ρ(x,)-ρ(x,0)=ρss(x).

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

(20) R(x,t)ρss(x)-ρ(x,t)ρss(x)-ρ(x,0)=1-ρ(x,t)ρss(x),

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

(21) τrev(x)=0R(x,t)𝑑t.

Since the Laplace transform of the relaxation function is R^(x,s)=0e-stR(x,t)𝑑t=s-1-ρ^(x,s)/ρss(x), the local accumulation time at position x, τrev(x)=lims0R^(x,s), is obtained as

(22) τrev(x)=𝒩/

with

𝒩2βout(coth[L/λ])2+csch[Lλ]{[(γ3)+βout(γ1)Lλ]cosh[Lxλ]+2βout1sinh[Lλ](γ1)Lxλsinh[Lxλ]}+coth[Lλ]csch[Lλ]([γLλ(2βout+Lλ)]cosh[Lxλ]+4sinh[Lλ]βout(γ1)Lxλsinh[Lxλ])

and

2kd(1+βoutcoth[Lλ])(βout-1+coth[Lλ]+(γ-1)cosh[L-xλ]csch[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 rs and vout are negligibly small. At the limits of rs0 and vout0, with the substitutions βout=vout/(λkd) and γ=jinkd/(voutrs), Equation 22 is simplified to

(23) τ(x)=limrs0vout0τrev(x)=12kd(1+Lλcoth[Lλ]Lxλtanh[Lxλ]).

Precision

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

(24) ϵrev2(x)σx2(x)x2=(logρss(x)logx)2σρ2(x)ρss(x)2=λ2kdvcellx2rs(βout1+coth[Lλ])(βout1+coth[Lλ]+(γ1)cosh[Lxλ]csch[Lλ])(csch[Lxλ]sinh[Lλ])2(γ1)2.

The expression for the squared relative error in the irreversible process of morphogen dynamics can be obtained at the limits of rs0 and vout0, with the substitutions βout=vout/(λkd) and γ=jinkd/(voutrs):

(25) ϵ2(x)=limrs0vout0ϵrev2(x)=λ3lcellx2(vcelljinlcellkd)1sinh[Lλ]cosh[Lxλ](sinh[Lxλ])2.

Thermodynamic cost

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

(26) S˙tot=1lcell0Ldx[S˙in(x)+S˙dep(x)+S˙diff(x)],

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

(27) S˙in(x)2vcell(jin-voutρss(0))δ(x)log(jinvoutρss(0)),
(28) S˙dep(x)vcell(kdρss(x)-rs)log(kdρss(x)rs),
(29) S˙diff(x)vcellD[xρss(x)][xlnρss(x)]=vcellD[xρss(x)]2ρss(x),

with the Boltzmann constant set to unity (kB=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 limlcell0Dlcell-2(ρss(x+lcell)-ρss(x))ln[ρss(x+lcell)/ρss(x)], in which Dρss(x+lcell)/lcell2 and Dρss(x)/lcell2 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 (S˙tot=0), the cells cannot infer any spatial information from it.

To obtain a simplified expression of 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

(30) 0L2δ(x)[jin-voutρss(0)]𝑑x=0L[kdρss(x)-rs]𝑑x.

Second, xρss(x)=xρss,0(x) and x2ρss(x)=x2ρss,0(x)=λ-2ρss,0(x), with ρss,0 defined in Equation 17, lead to the following equality:

(31) 0L[ρss,0(x)lnρss(x)]𝑑x=λ2([xρss,0(x)]lnρss(x)|0L-0L[[xρss(x)]2ρss(x)]𝑑x).

Then, S˙tot can be expressed as follows:

(32) S˙tot=1lcell0L(S˙in(x)+S˙dep(x)+S˙diff(x))dx=vcellkdlcell0L(ρss,0(x)lnjinkdvoutrs+ρss,0(x)lnρss(x)ρss(0)+λ2[xρss(x)]2ρss(x))dx=vcellkdlcell(0Lρss,0(x)lnγdx0Lρss,0(x)lnρss(0)dx+0L[ρss,0(x)lnρss(x)+λ2[xρss(x)]2ρss(x)]dx)=vcellkdlcell(0Lρss,0(x)lnγdxλ2[xρss,0(x)]lnρss(0)|0L+λ2[xρss,0(x)]lnρss(x)|0L)=vcellkdlcelllnγ0Lρss,0(x)dx=vcellλrslcell(γ1)lnγβout1+coth[Lλ].

We used 30 and 31 to obtain the second and fourth lines of Equation 32, respectively. The boundary condition of the system, 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 S˙tot (Equation 32) is normalized by lnγ, and rs and vout are negligibly small, we obtain

(33) limrs0vout0S˙totlnγ=vcelljinlcell.

Comparison of Equation 33 with the expression in Equation 5 enables us to relate αo with the free energy cost, Tlnγ. Thus, the thermodynamic cost associated with the morphogen profile is expressed as C=(αovcelljin/lcell)τ (Equation 5) where τ is the time required to generate the profile, and αovcelljin/lcell approximates the rate of dissipation, TS˙tot.

Trade-off product

The product of the local accumulation time (τrev(x), Equation 22), the squared relative positional error (ϵrev2(x), Equation 24), and the entropy production rate (S˙tot, Equation 32) leads to the position-dependent trade-off product πo,rev(x),

(34) πo,rev(x;λ,βout,γ)=S˙totτrev(x)ϵrev2(x)=12Llcellλ3Lx2csch[(Lx)/λ]sinh[L/λ]lnγ(γ1)(βoutcoth[L/λ]+1)×{(1γ)Lxλ+(γ3+βout(γ1)Lλ)coth[Lxλ]+4cosh[Lλ]csch[Lxλ]+coth[Lλ](βout(1γ)Lxλ+((γ1)Lλ2βout)coth[Lxλ]+2βoutcosh[Lλ]csch[Lxλ])+2βout1csch[Lxλ]sinh[Lλ]}.

With fixed L, lcell, and x=xb, the lower bound of πo,rev(xb;λ,βout,γ), πrev,min(xb), can be obtained by minimizing πo,rev(xb;λ,βout,γ) with respect to γ, βout, and λ. For a fixed value of γ>1, the partial derivative of πo,rev(xb;λ,βout,γ) with respect to βout is always smaller than 0,

(35) βoutπo,rev(xb;λ,βout,γ)=Llcellλ3L(xb)2lnγ(csch[Lxbλ])2×((sinh[L/λ])2βout2(γ1)+cosh[(Lxb)/λ]sinh[L/λ](2L/λ+sinh[2L/λ])4(βoutcosh[L/λ]+sinh[L/λ])2)<0.

Thus, πo,rev(xb) is minimized at the limit of large βout (or, equivalently, large βin=γβout), with the following expression:

(36) limβoutπo,rev(xb;λ,βout,γ)=14Llcellλ3L(xb)2(csch[(Lxb)/λ])2tanh[L/λ]lnγ(γ1)×{2+2cosh[2Lλ]2cosh[2Lxbλ]2cosh[xbλ]+(γ1)(xbλsinh[2Lxbλ]+(2Lxbλ)sinh[xbλ])}.

With the characteristic decay length set to λ=xb, it can be shown that limxb0βoutπo,rev(xb;xb,βout,γ)=0. Since πo,rev,min(xb)πo,rev(xb;xb,βout,γ) by the definition of πo,rev,min(xb), πo,rev,min(0)=0.

For xb>0, the optimal trade-off product, πrev,min(xb), and the corresponding values of λmin(xb) and γmin(xb), were numerically obtained (Appendix 1—figure 1). πo,rev,min(xb) 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(xb;λ,βout,γ) close to πo,rev,min(xb), the system can rapidly generate precise morphogen profiles at xb with low thermodynamic cost. In general, πo,rev,min(xb), λmin(xb), and γmin(xb) all increase monotonically with xb (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 vout,rs (see Equation 35), while the irreversible model is equivalent to the reversible version at the limits of vout,rs0.

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 vcelljin/lcell (Equation 33),

(37) 𝒞(x)=αo(vcelljinlcellkd)12(1+Lλcoth[Lλ]-L-xλtanh[L-xλ]),

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

(38) ϵ2(x)=λ3lcellx2(vcelljinlcellkd)-1sinh[Lλ]cosh[L-xλ](sinh[L-xλ])2.

Then, the trade-off product reads

(39) πo(x;λ)=𝒞(x)ϵ2(x)=αo2λ3lcellx2 (1+Lλcoth[Lλ]Lxλtanh[Lxλ])sinh[Lλ]cosh[Lxλ](sinh[Lxλ])2.

For fixed x=xb, πo(xb;λ) 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

(40) πo(x;λ)αo2λlcell[(λx)+(λx)2]ex/λ.

For a given position x=xb, λπo(xb;λ)|λ=λmin(xb)=0 determines λmin(xb)(13-1)xb/60.43xb (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

(41) πT(x;λ)=𝒞TϵT2(x)=αojinT(logm(x)logx)2σm2(x)m(x)2=αojinTλ2x2[2jinkdex/λsinh(a/λ)]2=[xxm(x)]2×4jinex/λTkd2sinh(a/λ)(1λcosh(a/λ)sinh(a/λ)+a2ea/λλsinh(a/λ))=σm2(x)=αo(λx)2ex/λsinh(a/λ)(1λcosh(a/λ)sinh(a/λ)+a2ea/λλsinh(a/λ))=αo(λx)2ex/λea/λ+3e3a/λ4(aλ1+1)ea/λ2(e2a/λ1)2,

The πT(x;λ) minimizing λ at x=xb, namely λmin(xb), is found at the interval (xb-a)/2λmin(xb)xb/2, since

(42) λπT(xb;λ)|λ=xb/2=αoae2+2a/xb(1+3e4a/xb)xb2(e4a/xb-1)3(sinh[4axb]-4axb)0,

and

(43) λπT(xb;λ)|λ=(xb-a)/2=αo8a2e2(xb+a)/(xb-a)xb2(xb-a)(e4a/(xb-a)-1)2([2axb-a]-1-coth[2axb-a])0.

For Figure 3 in the main text, the exact values of πT,min(xb) and λmin(xb) were obtained by numerically minimizing Equation 41 with respect to λ.

Appendix 2

Distributed morphogen source

Appendix 2—figure 1
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.05L, and (d) λp=0.1L. The color scale indicates the trade-off product πo normalized by αoL/lcell. The characteristic length scale λd is defined as λd0L𝑑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

(44) tρ(x,t)=Dx2ρ(x,t)-kdρ(x,t)+Jin(x),

with the boundary conditions xρ(x,t)|x=0,L=0. The production term, Jin(x) is assumed to be exponentially distributed with a characteristic length λp,

(45) Jin(x)=jinλpcosh[L-xλp]sinh[Lλp],

so that 0LJin(x)𝑑x=jin. The solution to Equation 44 in the Laplace domain is given by

(46) ρ^(x,s)=jinsD(kd+s)DD-(kd+s)λp2(cosh[L-xD/(kd+s)]sinh[LD/(kd+s)]-λpkd+sDcosh[L-xλp]sinh[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):

(47) ρss(x)=jinDkdλ2λ2-λp2(cosh[L-xλ]sinh[Lλ]-λpλcosh[L-xλp]sinh[Lλp]),
(48) ϵ2(x)=ρss(x)vcell(xxρss(x))2=kdλvcelljin(λ2-λp2)(cosh[L-xλ]sinh[Lλ]-λpλcosh[L-xλp]sinh[Lλp])x2[(sinh[L-xλ]sinh[Lλ]-sinh[L-xλp]sinh[Lλp])]2,

with λD/kd. Next, the characteristic timescale to approach the steady state profile (τ(x), Equation 22) can be written as

(49) τ(x)=𝒩/,

 with

𝒩(13λp2λ2+Lλ(1λp2λ2)coth[Lλ])cosh[Lxλ]sinh[Lλ]+2λp3λ3cosh[Lxλp]sinh[Lλp]Lxλ(1λp2λ2)sinh[Lxλ]sinh[Lλ],

and

2kd(1+λp2λ2)(λpλcosh[xλp]coth[Lλp]λpλsinh[xλp]cosh[Lxλ]sinh[Lλ]).

Then, the trade-off product is expressed as

(50) π(x;λ,λp)=αo(vcelljinlcell)τ(x)ϵ2(x),

with the proportionality constant αo. For fixed x=xb, L, and λp, the trade-off product can be minimized by tuning λ to the optimal value, λmin(xb). 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 λd0L𝑑x(ρss(x)-ρss(L))/(ρss(0)-ρss(L)).

Dynamics on a sphere

Appendix 2—figure 2
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 R2/λ2=10. The solid black line was computed from evaluating ρss(z)=lims0sρ^(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/lcell). The characteristic length scale λs is defined as λs01𝑑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
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. ϕρ=0, rρ=0, and r=R) the morphogen concentration satisfies the following evolution equation:

(51) ρ(θ,t)dt=D[1R2sinθθ(sinθρ(θ,t)θ)]-kdρ(θ,t)+jinδ(cosθ),

with boundary conditions at θρ|θ=0,π=0. The morphogen dynamics is fully specified by the system size R[lcell], the depletion rate kd[1/time], the diffusion coefficient D[lcell2/time], and the injection rate jin[conc/time]. 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 ρ^(z,s)=0𝑑te-stρ(z,t), can be obtained by solving

(52) z((1-z2)ρ^(z,s)z)-R2(kd+s)Dρ^(z,s)+R2jinDsδ(z)=0.

The solution to the differential equation of this type can be expressed as ρ^(z,s)=l=0cl(s)Pl(z), where Pl(z) (l=0,1,) is the Legendre polynomial, satisfying Pl(x)ddx[(1-x2)dPl(x)dx]=-l(l+1)Pl(x), with the orthogonality condition -11Pl(x)Pm(x)𝑑x=22l+1δlm. To find the coefficient cl(s), we rearrange Equation 52 into

-R2jinDsδ(z)=-l=0(l(l+1)+R2(kd+s)D)cl(s)Pl(z).

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

(53) ρ^(z,s)=R2jinDsl=0(2l+1)2(R2(kd+s)D+l(l+1))Pl(0)Pl(z).

The steady state concentration can be expressed as ρss(z)=lims0sρ^(z,s).

The average relaxation time can be calculated using τ(z)t=0t(-dR(t)dt)𝑑t=0R(t)𝑑t=lims0R^(z,s), which yields

(54) τ(z)=kd-1𝒯(z,λ/R),

with λD/kd and

𝒯(z,λ/R)=R2λ2l=0(2l+1)2(R2/λ2+l(l+1))2Pl(0)Pl(z)l=0(2l+1)2(R2/λ2+l(l+1))Pl(0)Pl(z).

The number of morphogen molecules produced from the cells at the equator per unit time is vcell(2πR/lcell)jin. Thus, the thermodynamic cost to approach the steady state concentration at z is

(55) C(z)=αo2πRvcelljinlcellkd𝒯(z,λ/R),

where αo is the proportionality constant.

Next, using ρss(z)=lims0sρ^(z,s) (Equation 53), one obtains the local precision of the morphogen profile at position z,

(56) ϵ2(z)=ρss(z)vcell(zzρss(z))2=kdvcelljin(z,λ/R),

where

(z,λ/R)λ2R2l=0(2l+1)2(R2/λ2+l(l+1))Pl(0)Pl(z)z2(l=0(2l+1)2(R2/λ2+l(l+1))Pl(0)zPl(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

(57) π(z;λ/R)C(z)×ϵ2(z)=αo2πRlcell𝒯(z,λ/R)(z,λ/R).

The length scale λs used in the plot is defined as λsz=0z=1𝑑z(ρss(z)-ρss(1))/(ρss(0)-ρss(1)). At the target position z=zb, 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 λBcd=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 xb=(140188)μm.

  • Wg. The characteristic decay length of Wg was set to λWg=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 xb=(1121)μm.

  • Hh. The characteristic decay length of Hh was set to λHh=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 xb=(1530)μm.

  • Dpp. The characteristic decay length of Dpp was set to λDpp=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 xb=(3654)μ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 (kd) of Cyclops and Squint have been measured through the ectopic expressions of the fluorescently tagged constructs. The resulting characteristic decay lengths (λD/kd) are λcyclops=76μm and λsquint=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, zb 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, λfgf8=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 zb 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.

Time scales of Bcd, Wg, Hh, and Dpp

The expression for the variance of the space-time-averaged signal, σm2(x) (Equation 11), is derived with the assumption of a long measurement time, Tkd-1 (Appendix 1 of Fancher and Mugler, 2020). Additionally, in order for Tjin 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τ(xb)). For Bcd, Wg, Hh, and Dpp, Lλ and xb/λ2, leading to the approximate expression for the characteristic time of τ(xb)1.5kd-1 (Equation 4). We assume that the maximum measurement time, Tmax, is set by the cell doubling time, after which the fates of the daughter cells must be re-established.

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

Article and author information

Author details

  1. 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
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8713-1294
  2. 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
    ORCID icon "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.

Version history

  1. Preprint posted: April 14, 2021 (view preprint)
  2. Received: May 4, 2021
  3. Accepted: August 13, 2021
  4. Accepted Manuscript published: August 17, 2021 (version 1)
  5. Version of Record published: September 22, 2021 (version 2)

Copyright

© 2021, Song and Hyeon

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

Metrics

  • 1,177
    Page views
  • 203
    Downloads
  • 4
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Yonghyun Song
  2. Changbong Hyeon
(2021)
Cost-precision trade-off relation determines the optimal morphogen gradient for accurate biological pattern formation
eLife 10:e70034.
https://doi.org/10.7554/eLife.70034

Share this article

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

Further reading

    1. Developmental Biology
    2. Neuroscience
    Tariq Zaman, Daniel Vogt ... Michael R Williams
    Research Article

    The cell-type-specific expression of ligand/receptor and cell-adhesion molecules is a fundamental mechanism through which neurons regulate connectivity. Here, we determine a functional relevance of the long-established mutually exclusive expression of the receptor tyrosine kinase Kit and the trans-membrane protein Kit Ligand by discrete populations of neurons in the mammalian brain. Kit is enriched in molecular layer interneurons (MLIs) of the cerebellar cortex (i.e., stellate and basket cells), while cerebellar Kit Ligand is selectively expressed by a target of their inhibition, Purkinje cells (PCs). By in vivo genetic manipulation spanning embryonic development through adulthood, we demonstrate that PC Kit Ligand and MLI Kit are required for, and capable of driving changes in, the inhibition of PCs. Collectively, these works in mice demonstrate that the Kit Ligand/Kit receptor dyad sustains mammalian central synapse function and suggest a rationale for the affiliation of Kit mutation with neurodevelopmental disorders.

    1. Developmental Biology
    2. Neuroscience
    Smrithi Prem, Bharati Dev ... Emanuel DiCicco-Bloom
    Research Article

    Autism spectrum disorder (ASD) is defined by common behavioral characteristics, raising the possibility of shared pathogenic mechanisms. Yet, vast clinical and etiological heterogeneity suggests personalized phenotypes. Surprisingly, our iPSC studies find that six individuals from two distinct ASD-subtypes, idiopathic and 16p11.2 deletion, have common reductions in neural precursor cell (NPC) neurite outgrowth and migration even though whole genome sequencing demonstrates no genetic overlap between the datasets. To identify signaling differences that may contribute to these developmental defects, an unbiased phospho-(p)-proteome screen was performed. Surprisingly despite the genetic heterogeneity, hundreds of shared p-peptides were identified between autism subtypes including the mTOR pathway. mTOR signaling alterations were confirmed in all NPCs across both ASD-subtypes, and mTOR modulation rescued ASD phenotypes and reproduced autism NPC associated phenotypes in control NPCs. Thus, our studies demonstrate that genetically distinct ASD subtypes have common defects in neurite outgrowth and migration which are driven by the shared pathogenic mechanism of mTOR signaling dysregulation.