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

Decision letter

  1. Gordon J Berman
    Reviewing Editor; Emory University, United States
  2. Naama Barkai
    Senior Editor; Weizmann Institute of Science, Israel
  3. Mariela Petkova
    Reviewer; Harvard University, United States
  4. 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 α0, 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αo(vcelljin/lcell)τ (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 αo494×42×103 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×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 only a small fraction of the total energy budget of Drosophila embryogenesis ~7×1016 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ΔS(t)δX(t)2X(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 2kB (Q715kB for the molecular motors, and Q4550kB 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 xb0.4L, πo(xb,λBcd)αo(L/lcell)0.6 (Figure 2b) where L/lcell50, and the experimentally reported precision is ϵ2(xb)7×104 (Gregor et al., 2007). Thus, the cost associated with the Bcd profile is C(xb,λBcd)=πo(xb,λBcd)/ϵ2(xb,λBcd)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×108 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(xb) 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 πC×ϵ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/nm1/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 (nm). – 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 nm (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×ϵ2 , which we called the trade-off product, is a dimensionless quantity, independent of the number of morphogens being produced (nm).

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ΔS(t)δX(t)2X(t)2. In the case of TUR, the two quantities of interest are a function of time t, such that ΔS(t)t, δX(t)2X(t)21/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, Q2kB, where kB 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, λmin0.43xb for the point measurement and λmin0.5xb 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 λmin0.43L in the former and λmin0.5L in the latter.

6) The discussion of the reversible model in the appendix is quite nice, particularly how you can derive all the irreversible results directly from the reversible model by taking an appropriate limit. Particularly interesting is that the cost defined in terms of entropy production rate times characteristic time, if normalized by ln(γ) [using Equation 33] would seem to give you the morphogen production cost of Equation 37, up to a factor of α0, assuming the irreversible limit where ln(γ) -> infinity. So it seems that the 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 (Ṡtot, Equation 32) is normalized by ln(γ), and vout and rs are negligibly small, we obtain the equality limrs,vout0Ṡtot/lnγ=vcelljin/lcell (Equation 33). Comparison of this with the expression of C (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 total dissipation, TṠ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

  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.

Senior Editor

  1. Naama Barkai, Weizmann Institute of Science, Israel

Reviewing Editor

  1. Gordon J Berman, Emory University, United States

Reviewers

  1. Mariela Petkova, Harvard University, United States
  2. Anatoly Kolomeisky, Rice University, United States

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,145
    Page views
  • 200
    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

Further reading

    1. Chromosomes and Gene Expression
    2. Developmental Biology
    Virginia L Pimmett, Mounia Lagha
    Insight

    Imaging experiments reveal the complex and dynamic nature of the transcriptional hubs associated with Notch signaling.

    1. Cell Biology
    2. Developmental Biology
    Simon Schneider, Andjela Kovacevic ... Hubert Schorle
    Research Article

    Cylicins are testis-specific proteins, which are exclusively expressed during spermiogenesis. In mice and humans, two Cylicins, the gonosomal X-linked Cylicin 1 (Cylc1/CYLC1) and the autosomal Cylicin 2 (Cylc2/CYLC2) genes, have been identified. Cylicins are cytoskeletal proteins with an overall positive charge due to lysine-rich repeats. While Cylicins have been localized in the acrosomal region of round spermatids, they resemble a major component of the calyx within the perinuclear theca at the posterior part of mature sperm nuclei. However, the role of Cylicins during spermiogenesis has not yet been investigated. Here, we applied CRISPR/Cas9-mediated gene editing in zygotes to establish Cylc1- and Cylc2-deficient mouse lines as a model to study the function of these proteins. Cylc1 deficiency resulted in male subfertility, whereas Cylc2-/-, Cylc1-/yCylc2+/-, and Cylc1-/yCylc2-/- males were infertile. Phenotypical characterization revealed that loss of Cylicins prevents proper calyx assembly during spermiogenesis. This results in decreased epididymal sperm counts, impaired shedding of excess cytoplasm, and severe structural malformations, ultimately resulting in impaired sperm motility. Furthermore, exome sequencing identified an infertile man with a hemizygous variant in CYLC1 and a heterozygous variant in CYLC2, displaying morphological abnormalities of the sperm including the absence of the acrosome. Thus, our study highlights the relevance and importance of Cylicins for spermiogenic remodeling and male fertility in human and mouse, and provides the basis for further studies on unraveling the complex molecular interactions between perinuclear theca proteins required during spermiogenesis.