Cellular compartmentalisation and receptor promiscuity as a strategy for accurate and robust inference of position during morphogenesis

  1. Krishnan S Iyer
  2. Chaitra Prabhakara
  3. Satyajit Mayor  Is a corresponding author
  4. Madan Rao  Is a corresponding author
  1. Simons Center for the Study of Living Machines, National Center for Biological Sciences - TIFR, India
  2. National Center for Biological Sciences - TIFR, India

Abstract

Precise spatial patterning of cell fate during morphogenesis requires accurate inference of cellular position. In making such inferences from morphogen profiles, cells must contend with inherent stochasticity in morphogen production, transport, sensing and signalling. Motivated by the multitude of signalling mechanisms in various developmental contexts, we show how cells may utilise multiple tiers of processing (compartmentalisation) and parallel branches (multiple receptor types), together with feedback control, to bring about fidelity in morphogenetic decoding of their positions within a developing tissue. By simultaneously deploying specific and nonspecific receptors, cells achieve a more accurate and robust inference. We explore these ideas in the patterning of Drosophila melanogaster wing imaginal disc by Wingless morphogen signalling, where multiple endocytic pathways participate in decoding the morphogen gradient. The geometry of the inference landscape in the high dimensional space of parameters provides a measure for robustness and delineates stiff and sloppy directions. This distributed information processing at the scale of the cell highlights how local cell autonomous control facilitates global tissue scale design.

Editor's evaluation

The manuscript introduces a compelling theoretical framework to investigate architectures of signal processing. The predictions of the computational model have been convincingly validated with data from fly wing precursor tissues. The work is important and will be highly valuable to biological physicists and developmental biologists interested in morphogenesis and pattern formation.

https://doi.org/10.7554/eLife.79257.sa0

Introduction

Precise positioning of cell fates and cell fate boundaries in a developing tissue is of vital importance in ensuring a correct developmental path (reviewed in Tkačik and Gregor, 2021; Wolpert, 2016). The required positional information is often conveyed by concentration gradients of secreted signalling molecules, or morphogens (reviewed in Tabata and Takei, 2004; Briscoe and Small, 2015). Typically, a spatially varying input morphogen profile is translated into developmentally meaningful transcriptional outputs. Morphogen profile measurements, across several signalling contexts, show that the gradients are inherently noisy Houchmandzadeh et al., 2002; Gregor et al., 2007a; Kicheva et al., 2007; Bollenbach et al., 2008; Zagorski et al., 2017. However, precision of the signalling output should be robust to inherent genetic or environmental fluctuations in the concentrations of the ligands and receptors engaged in translating the positional information. For example, the noisy profile of the morphogen Bicoid (Bcd) that activates hunchback (hb) in the early Drosophila embryo Gregor et al., 2007a; Gregor et al., 2007b , and the expression of gap genes that activate pair-rule genes Dubuis et al., 2013; Petkova et al., 2019 result in cell fate boundaries that are positioned to a remarkable accuracy of about one cell’s width. This points to a local, cell autonomous morphogenetic decoding that is precise and robust to various sources of noise Kerszberg and Wolpert, 2007; Kerszberg, 2004; Jaeger et al., 2004.

Cell autonomous decoding of noisy morphogen profiles includes reading of morphogen concentration, followed by cellular processing, finally leading to inference in the form of transcriptional readout. Several strategies have been proposed to ensure precision in output (reviewed in Barkai and Shilo, 2009; Lander et al., 2009): feedbacks such as self-enhanced morphogen degradation Eldar et al., 2002; Eldar et al., 2003, spatial and temporal averaging Gregor et al., 2007a, use of two opposing gradients McHale et al., 2006, pre-steady state patterning Bergmann et al., 2007 and serial transcytosis Bollenbach et al., 2007.

Most cell signalling systems have regulatory mechanisms that fine-tune signalling by controlling ligand-specific receptor interactions Rogers and Schier, 2011. Ligands such as TGF β/BMP Mueller and Nickel, 2012, Jiang and Cong, 2016, D’Souza et al., 2008, show promiscuous interactions with different receptors. Chen and Schier, 2002; Sick et al., 2006 or sequestering components within the extracellular matrix Marjoram and Wright, 2011 or interactions with binding receptors such as heparin sulphate proteoglycans (HSPGs) Baeg et al., 2001; Baeg et al., 2004; Yan and Lin, 2009 can control availability of the ligand. Additionally, the multiple endocytic pathways that operate at the plasma membrane can control the extent of signalling Bökel and Brand, 2014; Di Fiore and von Zastrow, 2014. These examples argue for distributed information processing within the cell.

In this paper, we show how cellular compartmentalisation, a defining feature of multicellularity, provides a compelling realisation of such distributed cellular inference. We show that compartmentalisation together with multiple receptors, receptor promiscuity and feedback control, ensure precision and robustness in positional inference from noisy morphogen profiles during development. Compartments associated with specific chemical (e.g. lipids, proteins/enzymes) and physical (e.g. pH) environments, have been invoked as regulators of biochemical reactions during cellular signalling and development Ellisdon and Halls, 2016; Omerovic et al., 2007; Omerovic and Prior, 2009; Shilo and Schejter, 2011; Bökel and Brand, 2014. Deploying promiscuous receptors against a morphogen, in addition to its specific receptor, is a strategy to buffer variations in morphogen levels. These observations provide the motivation for a general conceptual framework for morphogenetic decoding based on a multi-tiered, multi-branched information channel. While our framework has broader applicability, we will, for clarity, use the terminology of Wingless signalling in Drosophila wing imaginal disc Hemalatha et al., 2016.

Conceptual framework and quantitative models

We pose the task of morphogenetic decoding as a problem in local, cell autonomous inference of position from a morphogen input (Figure 1), where each cell acts as an information/inference channel with the following information flow:

  1. ‘reading’ of the morphogen input by receptors on the cell surface,

  2. ‘processing’ by various cellular mechanisms such as receptor trafficking, secondary messengers, feedback control, and

  3. ‘inference’ of the cell’s position in the form of a transcriptional readout.

Schematic of information processing in the developing tissue.

(a) A morphogen is produced by a specific set of cells (blue), and secreted into the lumen surrounding the tissue. Due to stochasticity of the production and transport processes, the morphogen concentration received by the rest of the cells is contaminated by extrinsic noise, which defines a distribution of morphogen concentration along the y-direction at any position x. (b) The route from morphogens to a developmental outcome requires each cell to read, process and infer its position. This task is further complicated by the stochasticity of the reading and processing steps themselves, that lead to intrinsic noise. (c) The problem of robust inference of position can be considered in a channel framework. The positional information is noisily encoded in the local morphogen (ligand) concentrations, p(L|x). The cells receive this as input and process it into a less noisy output to ensure robustness in inferred positions.

At a phenomenological level, reading of the morphogen input is associated with the binding of the morphogen ligand to various receptors with varying degree of specificity, leading to the notion that the information channel describing positional inference must possess multiple branches. Furthermore, the multiple processing steps associated with compartmentalisation of cellular biochemistry and/or signal transduction modules, for example phosphorylation states, provide the motivation for invoking multiple tiers in the channel architecture. At an abstract level, one may think of the branch-tier architecture of the cellular processing as a bipartite Markovian network/graph Hartich et al., 2014, with a fast direction (involving multiple branches) consisting of ligand-bound and unbound states along with chemical state changes, and a slower direction (involving multiple tiers) consisting of intracellular transport, fission and fusion, characterised by energy-utilising processes or a flux imbalance. A general developmental context with multiple morphogens may involve several such bipartite Markov networks/graphs with different receptors (or branches) in parallel. Some of these receptors could be shared between different morphogens. We refer to signalling receptors as those which transduce a signal upon binding to their specific morphogen ligand and non-signalling receptors as those that participate in the signalling pathway without directly eliciting a signalling response. At the end of processing, each individual cell may pool information from the various branches for the final inference of position, i.e. a transcriptional readout (Figure 2).

Schematic for the branch-tier channel architecture.

Branches correspond to different receptor types and tiers denote the layers of compartmentalisation used in cellular processing. Cellular processing associated with each receptor type (here, branches 1 and 2) is depicted by a generic Markov network. The gray and brown planes depict the tiers in the two branches respectively (here, tiers 1, 2, and 3 in each branch). The bi-directional in-plane purple arrows correspond to faster transitions between receptor states, e.g. bound/unbound, and the green bi-directional arrows depict slower transitions involving intracellular transport driven by flux-imbalanced processes. There may exist several feedback control loops (red ━┥ arrows) in the network. Ligand concentration L drives one or several reaction rates in such Markov networks as in Harvey et al., 2020. The output θ is a collection f of several signalling states (purple nodes) from one or many branches. The statistics of the output θ then enables inference of position.

The task of achieving a precise inference is complicated by the noise in morphogen input arising from both production and transport processes, and by the stochasticity of the reading and processing steps; thus the inference must be robust to the extrinsic and intrinsic sources of noise. The use of feedback control mechanisms is a common strategy to bring about robustness in the context of morphogen gradient formation and sensing Averbukh et al., 2017. Motivated by this, in Section ‘Quantitative models for cellular reading and processing’ we consider different feedback controls in conjunction with the tiers and branches. With these three elements to the channel architecture, the task of morphogenetic decoding can be summarised in the following objective.

Objective:
Given a noisy ligand input distribution at position x, i.e. p(L|x), what are the requirements on the reading (number of receptor types and receptor concentrations) and processing steps (number of tiers and feedback type) such that the positional inference is precise and robust to extrinsic and intrinsic noise?

Mathematical framework

Figure 1 describes information processing during development across a two dimensional tissue of nx, ny cells in x and y directions, respectively. The direction of morphogen gradient is taken to be along x, with the morphogen source to the left of x=0. Each cell is endowed with a chemical reaction network (CRN) with the same multi-tiered, multi-branched architecture with feedbacks described previously, that reads a noisy input L(x,y) (morphogen concentration) and produces an ‘output’ (biochemical ‘signal’) θ(x,y) that is also noisy. Here, we choose to construct the noisy morphogen profile in the following manner: for a given position x[0,1], cells along the y-direction see different amounts of ligand coming from the same input distributionp(L|x),

(1) p(L|x)=22πσL2(x)Exp[(LμL(x))22σL2(x)](1+Erf[μL(x)2σL(x)])1.

characterised parametrically by a mean μL(x) and standard deviation σL(x). Experimental data can be fit to this distribution Equation 1 (or another distribution suitable for the specific experimental system) to obtain the parameters. Here, we consider an exponentially decaying mean μL and standard deviation σL.

(2) μL(x)=Aex/λ
(3) σL(x)=μL(x)

Alternatively, one could choose a different parametrisation consistent with experimental observations for a morphogen profile with a monotonically decaying mean. The values of A,λ chosen for our analysis are listed in Table 1. The corresponding output distribution p(θ|x) can be used to infer the cell’s position. Since we do not know the precise functional relationship between the output and inferred position, we invoke Bayes rule MacKay and Mac Kay, 2003, as in previous work Tkačik et al., 2015, to infer the cell’s position,

(4) p(x|θ)=p(θ|x)p(x)p(θ)

where p(θ)=01dxp(θ|x)p(x) and p(x) is the prior distribution which we take to be uniform over a tissue of unit length, p(x)=1. We quantify precision in the inference by the local inference error, σX(x). For each position x, the inferred position x of cells along the y-direction is taken to be the maximum a posteriori estimate,

(5) x(x,y)=argmaxx~p(x~|θ(x,y))

where we use x~ to differentiate from the true position x. From this, the local and average inference error can be computed.

(6) σX2(x)=(xx)2y
(7) σ¯X=01σX(x)p(x)dx

where the average in Equation 6 is over cells in the y-direction. The logic behind this definition of the inference error is that development of the tissue relies on the precision in the inference of cells’ positions throughout the tissue. However, there may be tissue developmental contexts, where only the positions of certain regions or cell fate boundaries need to be specified with any precision (as in the case of short-range morphogen gradients like Nodal Liu et al., 2022). The definition of inference error may be readily extended to incorporate such specifications (see Section ‘Choice of objective function’).

Table 1
Parameters associated with rates, feedback and receptor profiles along with their range of values.

The chemical rate values used in numerical analysis are scaled by the unbinding rate ru,κu taken to be 1. The corresponding experimental values have been taken from Lauffenburger and Linderman, 1996 where available.

ParameterSymbolNumerical valuesExperimental values
Chemical rates
Signalling branch
Unbinding rateru10.34 min-1
Binding raterb0.1–1 nM-10.072 nM-1min-1
Degradation raterd0.001–0.010.0022 min-1
Internalisation raterI0.1–10.03–0.3 min-1
Recycling raterR0.1–10.058 min-1
Non-signalling branch
Unbinding rateκu1-
Binding rateκb0.1–1 nM-1-
Degradation rateκd0.001–0.01-
Internalisation rateκI0.1–1-
Recycling rateκR0.1–1-
Conjugation rateκC0.1–1 nM-1-
Splitting rateκS0.1–1-
Feedback control
Amplificationα0.1-10-
Feedback Sensitivityγ0-1
nM-1
-
Feedback strengthn0-5-
Receptor control
Signalling receptors
Hill coefficienta0-5-
Minimum concentrationA050-250 nM-
Maximum concentrationA0+A150-500 nM-
Position of half-maximumA20.01-1-
Non-signalling receptors
Hill coefficientb0-5
Minimum concentrationB050-250 nM-
Maximum concentrationB0+B150-500 nM-
Position of half-maximumB20.01-1-
Ligand input
Maximum concentrationA30 nM-
Decay lengthλ0.2-0.5-
Number of cells along x-directionnx101-
Number of cells along y-directionny101-

We have been motivated to use the maximum a posteriori (MAP) estimate in Equation 5 by its successful use in previous studies in Drosophila embryo Dubuis et al., 2013; Petkova et al., 2019; Tkačik et al., 2015 and, more importantly, that it is a local estimate not requiring the computation of p(θ) (which is independent of x). We have checked that a different definition of the inference error, which does not use the MAP estimate and takes into account the entire distribution p(x|x),

σx2(x)=01dx(xx)2p(x|x)

leads to the same qualitative results.

Quantitative models for cellular reading and processing

In order to calculate the probability of the inferred position given the output p(x|θ) and hence the inference error σ¯X, one needs to know the prior p(x) and the input-output relation giving rise to the output distribution p(θ|x) in Equation 4. While a uniform prior may be justified by a homogeneous distribution of cells in the developing tissue at the stage considered, the input-output relation needs to be developed using a specific model based on the general channel design principles described previously. Thus, we will take each cell to be equipped with a chemical reaction network (CRN) that has up to two receptor types both of which bind the ligand on the cell surface but only one is signalling competent Hemalatha et al., 2016; Tabata and Takei, 2004. This latter aspect breaks the symmetry between the receptor types and hence the branches, a point that we will revisit in Section ‘Asymmetry in branched architecture: promiscuity of non-signalling receptors’. In multi-tier architectures, the bound states of both the receptors are internalised and shuttled through several compartments. The last compartment allows for a conjugation reaction between the two receptors (as in the case of Wingless and Dpp Hemalatha et al., 2016; Zhu et al., 2020). The signalling states, defined by all the bound states of the signalling receptor, contribute to the output. Within this schema, we consider control mechanisms on the surface receptor concentrations and in the chemical reactions downstream to binding on the surface (i.e. on internalisation, shuttling, conjugation, etc). We formulate the control on processing steps as a feedback/feedforward regulation from one of the signalling species in the CRN. On the other hand, the control of surface receptors is considered in the form of an open-loop control by allowing receptor profiles to vary within certain bounds, as described below. The key parameters are chemical rate parameters describing the rates of various reactions in the CRN, receptor parameters describing the receptor concentration profiles, feedback topology in the CRN that is a combination of actuator and rate under regulation, control parameters describing the strength and sensitivity of the feedback/feedforward. With these parameters specified, an input-output relation, calculated as a tier-wise weighted sum of all signalling states, can then be used to infer the cell’s position by Equation 4.

Cellular Reading via surface receptors

In the framework described previously, we consider the morphogen ligand as an external input to the receiving cells, outside the cellular information processing channel. The signal and noise of this external input are captured by the distribution Equation 1. This implicitly assumes that there is no feedback control from the output to the ligand input, that is no ‘sculpting’ of the morphogen ligand profile. We revisit this point in the Discussion. Given a distribution of the morphogen input, we address the local, cell autonomous morphogenetic decoding that allows the cells to tune their reading dynamically.

We subject the local, cellular reading to an open-loop control on total (ligand bound plus unbound) surface availability of the signalling ψ and non-signalling ϕ receptors. This implies that for each evaluation of inference error within the optimisation routine (see Section ‘Performance of the Channel Architectures’), the local surface receptor levels are held constant in time through a chemostat (see Appendix 1). In our analysis, we consider a family of monotonic (increasing or decreasing in x and independent of y) receptor profiles, which for convenience we take to be of the Hill form (Figure 3), that is either

(8) Monotonically increasing in x:fA(x)=A0+A1xaA2a+xaor
(9) Monotonically decreasing in x:fB(x)=B0+B11+(x/B2)b

The range of values for these parameters considered in the numerical analysis are listed in Table 1. Therefore, when considering ψ(x) to be monotonically increasing in x, we parametrise it with fA. It follows that in a one-branch channel, there are two possibilities: ψ{fA,fB} while in a two-branch channel, there are a total of four possibilities: (ψ,ϕ){fA,fB}×{fA,fB}. This allows us to simulate the ‘reading’ step performed by the cells (see Figure 1b).

Family of receptor profiles fA (monotonically increasing in x) and fB (monotonically decreasing in x) with an interpretation of function parameters (Equations 8; 9).

The total surface concentrations of both signalling and non-signalling receptors are taken from these families of receptor profiles.

Note that we are not fixing a receptor profile but taking it from a class of monotonic profiles (including a uniform profile), over which we vary to determine the optimal inference (see Section ‘Performance of the Channel Architectures’ below). Further, in the optimisation scheme (Section ‘Performance of the Channel Architectures’), we allow the receptor concentrations to vary over the space of all monotonically increasing, decreasing or flat profiles, and do not encode the positional information in the receptor profiles. Monotonicity implicitly assumes a spatial correlation in the receptor concentrations across cells – we return to this point in the Discussion.

Dynamics of processing in a single-tier channel

In a single tier channel, all processing is restricted to the cell surface. We represent the bound state of the signalling receptor as R(1) and that of the non-signalling receptor as S(1). The conjugated state is represented by Q(1). The CRN for such a system with one and two branches is shown in Figure 4a. Rates associated with these reactions are listed in Table 1. The differential equations that describe the binding, unbinding, conjugation, splitting and degradation reactions of the receptors are given by

(10) tR(1)=rbL(ψ(x)R(1))(ru+rd)R(1)κCR(1)S(1)+κSQ(1)
(11) tS(1)=κbL(ϕ(x)S(1))(κu+κd)S(1)κCR(1)S(1)+κSQ(1)
(12) tQ(1)=κCR(1)S(1)κSQ(1)

The steady-state output θ, defined as the sum of all the ligand-bound signalling states, is given by θ=R(1)+Q(1). Note that to describe the 1-branch system, we simply set all rates κ to zero.

Examples of channel architectures with single and multiple tiers, and upto two branches.

Signalling receptors in the bound state (colour purple) from each of the tiers contribute to the cellular output. The interpretation of the arrows is shown in the legend.

Dynamics of processing in a multi-tier channel

In a multi-tiered channel, the receptors go through additional steps of processing before generating an output. We represent the bound state of a receptor in k-th tier of the first branch as R(k), that of the second branch as S(k), and the conjugate species that forms in the last nT-th tier as Q(nT). The CRN for such a system with nT tiers is shown in Figure 4b. Rates associated with these reactions are listed in Table 1. The differential equations that describe the binding, unbinding, trafficking, recycling, conjugation, splitting and degradation reactions of the receptors are given by

(13) tR(1)=rbL(ψ(x)R(1))(ru+rd+rI)R(1)+rRR(2)
(14) tS(1)=κbL(ϕ(x)S(1))(κu+κd+κI)S(1)+κRS(2)
(15) tR(nT)=rIR(nT1)(rR+rd)R(nT)κCR(nT)S(nT)+κSQ(nT)
(16) tS(nT)=κIS(nT1)(κR+κd)S(nT)κCR(nT)S(nT)+κSQ(nT)
(17) tQ(nT)=κCR(nT)S(nT)κSQ(nT)

The output, realised from all the ligand-bound signalling states, now becomes θ=wnTQ(nT)+k=1nTwkR(k) at steady state with wk, such that kwk=1, representing the weight allotted to the tier (according to the mean residence time in the tier, for instance). For details regarding the setup of Equations 10–17 refer to Appendix 1. These differential equations for single-tiered and multi-tiered systems are to be augmented by stochastic contributions from both extrinsic and intrinsic sources. Extrinsic noise is a consequence of stochasticity of the ligand concentration presented to the cell, Lp(L|x), and enters the equations as a source term. On the other hand, intrinsic noise is a consequence of copy-number fluctuations in the CRNs that characterise the channel, and are treated using chemical master equations (CMEs) Sengupta, 2008.

Feedback Control

We consider all rates in the CRN, except the ligand binding and unbinding rates, as potentially under feedback regulation. Any chemical rate r{rI,κI,κC,.} that is under feedback control actuated by the node R{R(1),S(1),.} is modelled as.

(18) r+=r0(1+αRnγn+Rn)if under positive feedback
(19) r=r01+(γR)nif under negative feedback

with r0 as the reference value of the chemical rate in the absence of feedback. The range of values for amplification α, feedback sensitivity γ and feedback strength n are listed in Table 1. Figure 5 shows the different categories of possible feedback controls. We discuss the heuristics underlying the feedback controls in Appendix 2.

Schematic of feedback types.

(A) In a one-branch channel, feedbacks are considered on internalisation rates or degradation rates. (B) A second branch in the channel opens up the possibilities of (a) intra-branch and (b) inter-branch, (i) intra-tier and (ii) inter-tier feedbacks.

Performance of the channel architectures

With the model in place, we address the Objective discussed previously, by studying the performance of different channel architectures, i.e. number of tiers and branches, and feedback topology. We define a vector v belonging to a parameter space V of the channel parameters related to chemical rates, receptor profiles and feedback (see Table 1). While the chemical rates and feedback parameters are the same in all cells, the receptor profile parameters help define the receptor concentrations at each cell position x,y. For a given morphogen input distribution p(L|x) and a channel architecture under consideration, the optimisation can be stated as

(20) σ¯Xopt=minvVσ¯X(v;p(L|x))

and implemented by the following algorithm, the details of which are presented in Appendix 3.

Optimisation scheme
  1. Fix a morphogen input distribution for each position, p(L|x) using Equation 1.

  2. Define the channel architecture hierarchically, i.e. first declare the number of tiers and branches in the channel, and then choose a feedback topology (as in Figure 5).

  3. Optimise the average inference error Equation 20 w.r.t. to the channel parameters vV within the bounds provided in Table 1. We use a gradient independent method viz. Pattern Search algorithm for this step (implemented in MATLAB). For every poll (iteration) of the Pattern Search, we evaluate the average inference error σ¯X using the steady-state outputs of the equations corresponding to the CRN under optimisation that is Equations 10–17. The steady state solution is obtained analytically when possible or solved using ODE15s (MATLAB) algorithm.

  4. Repeat Step 3 until all feedback topologies under consideration are exhausted.

  5. Repeat Steps 2 and 3 until all channel architectures are scanned.

Results

As discussed previously, cells of a developing tissue face both extrinsic as well as intrinsic sources of noise. We first look at the issue of extrinsic noise in the morphogen input (described by Equation 1). The output then is a deterministic function of the morphogen input and parameters of the channel i.e. receptor concentrations, feedback topology, chemical rates and feedback parameters. The range of values considered for these parameters is listed in Table 1, consistent with the timescale separation between the rates of chemical reactions and transport as discussed in Section ‘Conceptual framework and quantitative models’. We apply the numerical analysis and the optimisation algorithm outlined in Section ‘Performance of the Channel Architectures’ to determine the design characteristics of ‘reading’ (receptor profiles) and ‘processing’ (tiers and feedback control) steps. Later, we check how channels, optimised in the reading and processing steps to deal with extrinsic noise, respond to intrinsic noise and what roles the elements of channel architecture play there. All the essential results are presented in this section and the reader may look up the appendices for further details.

Branched architecture with multiple receptors provides accuracy and robustness to extrinsic noise

We begin with architectures comprising single-tiered channels with one and two branches. Such architectures are similar in design to the classic picture of ligand-receptor kinetics Lauffenburger and Linderman, 1996; Alberts et al., 2017, but also to the self-enhanced degradation models for robustness of morphogen gradients Eldar et al., 2003. Before we proceed, it helps to recall a simple heuristic regarding signal discrimination. Appendix 4—figure 1 illustrates that precision in positional inference requires both that the output variance at a given position be small and that the mean output at two neighbouring positions be sufficiently different.

Let us first consider a minimal architecture of a one-tier one-branch channel without feedback control on any of the reaction rates. The output of this channel, here R(1), is a monotonic, saturating function of the input, with the surface receptor concentration setting the asymptote. As in Appendix 5—figure 1a, if the receptor concentrations decrease with mean ligand input, i.e. increases with distance from source (fA in Figure 3a; ), the outputs for different input ranges overlap significantly. On the other hand, if the receptor concentrations increase with mean input (fB in Figure 3b), the outputs overlap to a lesser degree (see Appendix 5—figure 1b). Thus within this minimal architecture, the inference error is optimised when the receptor concentrations increase with the mean input.

Introducing a feedback in this one-tier one-branch architecture, either on receptor levels or degradation rate, only partially reduces the inference errors (Figure 6a, c). As seen in Figure 6d, this is because the surface receptor concentration ψ sets both the asymptote and the steepness of the input-output functions, resulting in significant overlaps between outputs at neighbouring positions. The receptor control introduces a competition between robustness of the output to input noise and sensitivity to systematic changes in the mean input (see Appendix 4).

Characteristics of an optimised (a) one-tier one-branch channel with only the signalling receptor and feedback.

The optimised channel shows a moderately strong positive feedback on the degradation rate. (b) The optimal output is obtained when (b, inset) the total (bound plus unbound) signalling receptor concentration profile decreases away from the source. (c) Local inference errors in this optimised channel show a reduction compared to the expected inference errors from ligand with no cellular processing (i.e. reading directly from the free ligand). The minimum average inference error in this channel is σ¯X8%, which corresponds to 8 cells’ width. The dashed line denotes a local inference error of one cell’s width 1/nx. (d) The input-output relations in this channel are monotonically increasing sigmoid functions saturating at only large values of input. The solid lines correspond to the input-output relations at selected positions x=0.25,0.5,0.75, shaded with the same colour as the position-markers in (b inset, coloured rectangles). The signalling ψ(x) receptor concentration is mentioned in the legend. For a fixed distribution of ligand input (Equation 1), the range of input values recorded by the receptors at the selected positions gives rise to a range of outputs (circles). It is clear that neighbouring positions have significant overlaps in their outputs. The optimised parameter values for the plots in (b–d) can be found in Table 2 under the column corresponding to nT=1,nB=1,r+=rd(1).

Including a non-signalling receptor ϕ via an additional branch in the channel architecture opens up several new possibilities of feedback controls, in addition to providing an extra tuning variable. Now, as opposed to the one-tier one-branch case, an inter-branch feedback control (Figure 7a) results in an input-output relation with a sharp rise followed by a saturation (Figure 7d). By appropriately placing the receptors at spatial locations that receive different input, as shown by black arrow in Figure 7d, one can cleanly separate out the cellular outputs in neighbouring positions. For a detailed description see Appendix 6. This mitigates the above-mentioned tension between robustness to input noise and sensitivity to systematic changes in the mean input to a considerable extent (see Appendix 4—figure 2).

Results of optimisation of (a) one-tier two-branch channel.

(b) The output profile (with standard error in shaded region) corresponding to the (inset) optimised signalling (blue) and non-signalling (red) receptor profiles. The optimal signalling receptor now increases away from the source as opposed to the situation in the optimal one-tier one-branch channel (Figure 6). On the other hand, the optimal non-signalling receptor decreases away from the source. (c) The local inference error σX(x) is reduced throughout the tissue, when compared to the expected inference errors from ligand with no processing. (d) The input-output relations at selected positions x=0.25,0.5,0.75 (in the direction of the black arrow) are shown as solid lines, shaded with the same colour as the position-markers in (b inset, coloured rectangles). The signalling ψ(x) and non-signalling ϕ receptor concentrations are mentioned in the legend. For a fixed distribution of ligand input (Equation 1), the range of input values recorded by the receptors at the selected positions gives rise to a range of outputs (circles). Tuning of input-output relations through receptor concentrations reduces output variance and minimises overlaps in the outputs of neighbouring cell cohorts. The optimised parameter values for the plots in (b–d) can be found in Table 2 under the column corresponding to nT=1,nB=2,r-=κC.

As seen in Figure 7c, the two-branch architecture with inter-branch feedback leads to a dramatic reduction in the inference errors, to reach one cell’s width precision at most spatial locations in the tissue.

We would like to highlight two unexpected features of the optimised two-branch architecture. (i) The signalling and non-signalling receptors present opposing optimal profiles – a consequence of the negative inter-branch feedback. (ii) The optimal non-signalling receptor decreases away from the source, indicating that the non-signalling receptor ‘reads’ the ligand input, while the signalling receptor increases away from the source, buffering the noise in the output (Figure 7). A heuristic understanding of the opposing optimal receptor profiles is provided in Appendix 7. In contrast, in the one-branch architectures, it is the signalling receptor that does the reading and buffering.

Table 2
Values of rates, feedback and receptor control parameters obtained after optimising the different channel architectures with nT tiers and nB branches.

The optimised values of the chemical rates quoted below are scaled by the unbinding rate ru, κu taken to be 1. The symbols r- and r+ denote positive and negative feedbacks, respectively, on the rates following the equals sign; {} implies absence of feedback.

Parameter (Symbol)Value obtained in the optimised channel (nT,nB)
(1,1)(1,2)(2,2)(2,2)(2,2)
r+=rd(1)r-=κCr-=κIr-=κCr-={}
Chemical rates
Signalling branch
Binding rate (rb,nM-1)0.08980.09490.09320.08930.0787
Degradation rate in tier 1 (rd(1))0.00130.00810.00860.00980.0038
Degradation rate in tier 2 (rd(2))--0.00660.00870.0016
Internalisation rate (rI)--0.05310.07840.0363
Recycling rate (rR)--0.06810.03590.0758
Non-signalling branch
Binding rate (κb,nM-1)-0.05900.09540.08350.0288
Degradation rate in tier 1 (κd(1))-0.00860.0010.00430.0068
Degradation rate in tier 2 (κd(2))--0.00370.00310.0033
Internalisation rate (κI)--0.07410.08460.0559
Recycling rate (κR)--0.01230.01340.0998
Conjugation rate (κC,nM-1)-0.99260.98230.97220.6019
Splitting rate (κS)-0.12850.15450.13500.7512
Feedback control
Amplification (α)3.2085----
Feedback Sensitivity (γ)0.24910.18310.55350.8259-
Feedback strength (n)2.68252.36832.09532.1880-
Tier-wise weights
weight of tier 1 (w1)110.00180.12320.9259
weight of tier 2 (w2)--0.99820.87680.0741
Receptor control
Signalling receptors
Hill coefficient (a)4.92311.99743.83633.52513.3835
Minimum concentration (A0,nM)51.813051.096069.694051.977051.2
Maximum concentration (A0+A1,nM)298.283290.356304.114134301
Position of half-maximum (A2)0.47520.78180.94050.83440.4091
Non-signalling receptors
Hill coefficient (b)-4.89511.08021.74723.1821
Minimum concentration (B0,nM)-192.32248.69192.494.1850
Maximum concentration (B0+B1,nM)-442489.77441.67305
Position of half-maximum (B2)-0.74280.51770.31960.0902

Tiered architecture with compartmentalisation adds robustness to intrinsic noise

We next investigate the effects of addition of tiers (compartments) on the inference errors. Our optimisation shows there are two distinct optimised two-tier two-branch architectures, one with inter-branch feedback on the internalisation rate of the non-signalling receptors κI and the other on the conjugation rate κC, that have comparable inference errors (Figure 8b, c). Both the receptor profiles and the input-output relations of these two optimised two-tier two-branch channels are qualitatively similar (Appendix 8—figure 1).

Performance of the optimised two-branch channels with increasing numbers of tiers.

(a) Minimum average inference error σ¯X in two-branch architectures with increasing number of tiers nT. The dashed line corresponds to a local inference error of one cell’s width 1/nx. (b,c) Results of optimisation of two-tier two-branch channels with inter-branch feedback. These two architectures perform equally well: local inference errors in both the channels (blue dots) are low throughout the tissue (with average inference errors 1.6% and 1.7%) as compared to a case with no processing of ligand prior to inference (black dots). Note that the local inference errors in the optimised channels increase towards the end of the tissue due to lower ligand concentrations. The dashed line corresponds to a local inference error of one cell’s width 1/nx. The optimised parameter values for the plots in (b–c) can be found in Table 2 under the column corresponding to nT=2,nB=2,r-=κI and nT=2,nB=2,r-=κC, respectively.

It would seem that addition of further tiers, that is more than two, would lead to further improvement in the inference. However, in both these optimised architectures, addition of tiers leads only to a marginal reduction of inference errors (Figure 8a) while invoking a cellular cost. Of course, extensions of our model that involve modification of the desired output could favour the addition of more tiers. For instance, additional tiers could facilitate signal amplification or improvement in robustness to input noise through an increase in signal-to-noise ratio (SNR) Stoeger et al., 2016. Further, by making the output θ a multi-variate function of the tier index (compartment identity) one can multitask the various cellular outcomes (as in Ras/MAPK signalling Fehrenbacher et al., 2009 or with GPCR compartmentalisation Ellisdon and Halls, 2016).

So far, we have only considered noise due to fluctuations in the morphogen profile, that is extrinsic noise. Given that we are considering a distributed channel, intrinsic noise due to low copy numbers of the reacting species in the CRN will have a significant influence on the inference. As discussed in Section ‘Conceptual framework and quantitative models’ and Appendix 3, we solve the stochastic chemical master equations (CMEs) to compute the output distributions and the positional inference. It is here that we find that the addition of tiers contribute significantly to reducing inference errors. A comparison of the one-tier two-branch and two-tier two-branch channel architectures (Figure 9a and b) optimised for extrinsic noise, shows that in the presence of intrinsic noise, additional tiers lead to significantly lower inference errors (Figure 9c). The large inference errors seen in the one-tier one-branch channel in the presence of intrinsic noise, can be traced to the instabilities of steady-state trajectories of the two signalling species R(1) and Q(1) driven by the non-linear feedback (Figure 9d–f). This effect is more prominent for larger values of ligand concentrations, that is closer to the source at x=0. On the other hand, we find that in the two-tier two-branch architecture (Figure 9g–i), the fluctuations in the signalling species are more tempered, the inter-branch feedback leads to a mutual damping of the fluctuations of the signalling species from the two branches. Details of this heuristic argument appear in Appendix 9.

Robustness to intrinsic noise in (a) one-tier two-branch (1T2B) channel and (b) two-tier two-branch (2T2B) channel architectures, previously optimised for extrinsic noise alone.

(c) A comparison of local inference errors due to intrinsic noise shows consistently better performance in the case of a two-tier two-branch channel (red dots). (d-f) Sample steady-state trajectories of the signalling species R(1) (blue) and Q(1) (red) of a one-tier two-branch channel (purple nodes in (a)) at positions x=0.1,0.4,0.7, respectively. (g–i) Sample steady-state trajectories of the signalling species R(2) (blue) and Q(2) (red) of a two-tier two-branch channel (purple nodes in (b)) at positions x=0.1,0.4,0.7, respectively. The optimised parameter values for the plots in (c,d–f,g–i) can be found in Table 2 under the column corresponding to nT=1,nB=2,r-=κC and nT=2,nB=2,r-=κC, respectively.

In summary, we find that the nature of the channel architectures play a significant role in robustness of morphogenetic decoding to both extrinsic and intrinsic sources of noise. Of the three elements to the channel architecture - branches, tiers, and feedback control, we find that a branched architecture can significantly reduce inference errors by employing an inter-branch feedback and a control on its local receptor concentrations. For this, the receptor concentration profiles required to minimise inference errors are such that the concentration of signalling (non-signalling) receptor should decrease (increase) with mean morphogen input. Crucially, in the absence of feedback, performance of the channel diminishes and the optimised receptor profiles both decrease away from the source (Appendix 10—figure 1). Further, we show in Appendix 11—figure 1 that having uniform profiles for the signalling and non-signalling receptors, with or without uncorrelated noise, fares poorly in terms of inference capability. This provides a posteriori justification for the monotonicity in receptor profiles. Addition of tiers can help in further bringing down inference errors due to extrinsic noise, but with diminishing returns. An additional tier, however, does provide a buffering role for feedback when dealing with intrinsic noise. We note that these qualitative conclusions remain unaltered for different morphogen input characteristics, that is input noise and morphogen decay lengths (see Appendix 12).

Asymmetry in branched architecture: promiscuity of non-signalling receptors

Before comparing the theoretical results with experiments, we comment on the implications for the cellular control of the signalling ψ and non-signalling ϕ receptors. In the two-branch architecture, the symmetry between the signalling and non-signalling receptors is broken by the inter-branch feedback and the definition of output θ, the latter taken to be a function only of the signalling states R(k) and Q(k) (Section ‘Conceptual framework and quantitative models’, purple nodes in Figure 7a and Figure 8b and c). What are the phenotypic implications of this asymmetry? In Appendix 13—figure 1, we plot the contours of average inference errors σ¯X in the ψ-ϕ plane around the optimal point. We compute the eigenvalues of the local curvature of σ¯X(Δψ,Δϕ) around the optimal point (Δψ=Δϕ=0). The difference in the magnitudes of these eigenvalues, as discussed in Appendix 13, immediately describes stiff and sloppy directions Transtrum et al., 2015 along the ψ and ϕ axes, respectively. This implies that while the signalling receptor is under tight cellular control, the control on the non-signalling receptor is allowed to be sloppy. A similar feature is observed in the contour plots for the robustness measure χ (defined as the ratio of coefficients of variation in the output to that in the input). Appendix 4—figure 3 shows that for any given input distribution, reduction in output variance requires a stricter control on ψ, while the control on ϕ can be lax.

This sloppiness in the levels of non-signalling receptor would manifest at a phenotypic level in the context of multiple morphogen inputs as in the case of Drosophila imaginal disc Tabata and Takei, 2004. Participation of the same non-signalling receptor in the different signalling networks would imply its promiscuous interactions with all ligands. The signalling receptors, therefore, are specific for the various ligands while the non-signalling receptor, being promiscuous, is non-specific. This, as we see below, is the case with the Heparan sulfate proteoglycans (HSPGs) such as Dally and Dally-like protein (Dlp) that participate in the Wingless (Wg) and Decapentaplegic (Dpp) signalling networks Lin and Perrimon, 2000; Romanova-Michaelides et al., 2021.

Geometry of fidelity landscape

The above section and Appendix 13 motivate us to study the changes in the inference error upon perturbations of all the channel parameters. We therefore discuss the nature of optima in terms of the local geometry of the fidelity landscape around the optimum, and the geometry of the low inference error states. We work with the case of the optimised one-tier two-branch channel (shown in Figure 7a with optimum channel parameters listed in Table 2, Table 3) in presence of extrinsic noise.

Table 3
Values of chemical rates and feedback parameters obtained after optimising the two-tier two-branch channel with inter-branch feedback on the internalisation rate κI of the non-signalling branch, keeping the receptor profiles spatially uniform, with and without uncorrelated noise.

The optimised values of the chemical rates quoted below are scaled by the unbinding rate ru, κu taken to be 1.

Parameter (Symbol)Optimised value
uniform receptor profilesuniform receptor profiles with uncorrelated noise
Chemical rates
Signalling branch
Binding rate (rb,nM-1)0.09220.0782
Degradation rate in tier 1 (rd(1))0.00890.0041
Degradation rate in tier 2 (rd(2))0.00920.0095
Internalisation rate (rI)0.02250.0611
Recycling rate (rR)0.04030.0971
Non-signalling branch
Binding rate (κb,nM-1)0.04640.0265
Degradation rate in tier 1 (κd(1))0.00350.0045
Degradation rate in tier 2 (κd(2))0.00710.0068
Internalisation rate (κI)0.020.0513
Recycling rate (κR)0.09890.0770
Conjugation rate (κC,nM-1)0.76050.7579
Splitting rate (κS)0.70380.3036
Feedback control
Amplification (α)--
Feedback Sensitivity (γ)0.09390.1946
Feedback strength (n)4.63100.6202
Tier-wise weights
weight of tier 1 (w1)0.00460.2875
weight of tier 2 (w2)0.99540.7125

To address the geometry of the local fidelity landscape around the optimum, we compute (i) percent changes in inference error σ¯X due to perturbations in channel parameters (Figure 10a), and (ii) the eigenspectrum of the Fisher information metric (FIM, Figure 10b). The FIM gμν is evaluated in the log-parameter space as Transtrum et al., 2015.

(21) gμν=xiyjx(M(xi,yj),v)lnvμx(M(xi,yj),v)lnvν

where, vV is the channel parameter vector, and xi,yj are the indices of cells that run along the x- and y-directions. As shown in Figure 10a, we see that the inference error does not change significantly (up to 20% change with most parameters), that is it remains within σ¯X2.2%. Varying the feedback strength n, however, drives a much stronger deviation from the minimum. Similarly, as seen from the heat map (Figure 10b), eigenvectors with the larger eigenvalues (index 1–6) have an appreciable component of the feedback parameters γ,n. This implies that variation of the feedback parameters from the optimum would result in significant changes in the inferred positions. Perturbing conjugation κC and splitting κS rates simultaneously (see eigenvector 16) does not produce any notable change to the inferred positions (eigenvalue 10-13). Further, perturbations to channel parameters other than the feedback parameters (eigenvectors 7–16) produce marginal changes in inferred positions.

Geometry of the fidelity landscape around the optimum.

(a) Percent changes in the inference error upon perturbations in the channel parameters (as described in Table 1) around the optimum for one-tier two-branch channel (optimised σ¯X=1.9%). For most perturbations, the inference error deviates by up to 20% of the optimum i.e. the inference error σ¯X remains below 2.2%. (b) Left: eigen spectrum of the Fisher information metric (FIM, see Equation 21) around the global minimum of σ¯X, Right: weight of the different channel parameters in the eigenvectors of FIM, obtained from projecting each eigenvector along the channel parameter axes. The index 1 corresponds to the eigenvector with the largest eigenvalue and the index 16 corresponds to the eigenvector with the smallest eigenvalue.

Moving now from a local to global analysis of the fidelity landscape, we run the optimisation algorithm (Section ‘Performance of the Channel Architectures’) on the one-tier two-branch channel architecture with 216 space-filling initial points in the 16-dimensional parameter space of this architecture. We then define the low inference error states as those channel parameters vopt that yield σ¯X2%. This cutoff, which equals σ¯X0.01, corresponds to declaring as equivalent all the inference errors σ¯X that lie between one and two cells’ widths. Consistent with the local analyses, we find that the frequency distribution of optimal feedback parameters γ,n is narrowly distributed about the global optimum (Figure 11a). As shown in Figure 11a, the parameters corresponding to forward and backward rates are skewed towards the upper and lower bounds of the allowed parameter range, respectively. We see that the optimal binding rates in the non-signalling branch (Figure 11a) are more broadly distributed across the permissible range than the optimal binding rates in the signalling branch, which are concentrated towards the upper bound of the permissible range. This again reflects the promiscuity of the non-signalling receptors as described in Section ‘Asymmetry in branched architecture: promiscuity of non-signalling receptors’. All other optimal parameters corresponding to degradation rates, minimum and maximum receptor values and steepness of the receptor profiles, show a very broad spread over this range (Appendix 14—figure 1). To explore the topography of the low inference error landscape, we evaluate the components of the ‘position vectors’ of these minima vopt in the parameter space V along the eigenvectors of the Hessian of σ¯X, defined as

(22) hμν=2σ¯X(M,v)vμvν

where M stands for the entire morphogen profile and we have assumed a Euclidean metric. As shown in Figure 11b and c, components of the ‘position vector’ of the minima voptV lie predominantly along the sloppy directions of the Hessian that is along the eigenvectors with small eigenvalues. This suggests that geometry of the low inference error landscape resembles a deep valley, which is shallow along the several sloppy directions and steep along the few stiff directions.

Geometry of the low inference error landscape defined by channels within a band σ¯X2% about the global minimum.

(a) Frequency distributions of optimised channel parameters in the low inference error landscape. Here we show the ligand binding rates of the signalling and non-signalling receptors, conjugation and splitting rates, and feedback sensitivity and feedback strength parameters. The distributions of the other optimised channel parameters are shown in Appendix 14. (b) Eigenvalues of the Hessian hμν (see Equation 22) of σ¯X around the global minimum. (c) Components of the normalised ‘position vectors’ of the minima voptV along the eigenvectors of the Hessian hμν, obtained from projecting each position vector along the eigenvector of the Hessian. Here, position vectors in the parameter space V are defined by the usual Euclidean metric.

Choice of objective function

The objective function as defined in Equation 7 gave equal weight to inference errors at all positions x along the tissue, driving the inference error to reduce at all positions simultaneously. In certain developmental contexts, the objective could be to partition the tissue into cell identity segments (reviewed in Briscoe and Small, 2015). In such a case, the partition boundaries would need to be sharp Gregor et al., 2007a that is only the errors at the segment boundaries would need to be minimised. We show that even with this choice of objective function, the qualitative results for the optimal channel architectures remain unaltered. We define the inference error for a tissue with Np segmented cell identities as.

(23) σX2(x)=(1δg(x),g(x))(xx)2ywhereg(x)=1+i=1NpΘ(xξi)

where δ and Θ denote the Kronecker-delta and Heaviside-theta functions respectively, ξi is the position of the boundary between the i-th and (i+1)-th segments, and g is a function that maps position (actual or inferred) to a segment, that is g:[0,1]{1,2,,Np} with Np as the total number of segments.

We optimise one-tier one-branch, one-tier two-branch and two-tier two-branch channel architectures for the inference error as defined in Equation 23 with Np=4 and equally spaced boundaries located at positions ξ1=0.25,ξ2=0.5,ξ3=0.75 along the x-axis. As before, this optimisation suggests that an additional branch aids in reducing the inference errors due to extrinsic noise (compare Figure 12b and d), with similar opposing receptor profiles as in Section ‘Branched architecture with multiple receptors provides accuracy and robustness to extrinsic noise’. Tiers play only a moderate role in reducing the inference errors further in a two-branch channel (compare Figure 12d and f). However, just as with the previous objective function, an additional tier provides substantial robustness to intrinsic noise as shown in Figure 13c.

Robustness to extrinsic noise with a different choice of objective function (Equation 23) for one-tier one-branch channel (a–b), one-tier two-branch channel (c–d) and two-tier two-branch channel (e–f).

(a) Profile of the signalling receptor for (a, inset) the optimised one-tier one-branch channel. (b) Corresponding inference errors due to extrinsic noise in the optimised one-tier one-branch channel. (c) Profiles of the signalling (blue) and non-signalling (red) receptor for (c, inset) the optimised one-tier two-branch channel. (d) Corresponding inference errors due to extrinsic noise in the optimised one-tier two-branch channel. Errors are predominantly located around the segment boundaries at x=0.25,0.5,0.75 and still increase in the direction of reducing morphogen concentrations. (e) Profiles of the signalling (blue) and non-signalling (red) receptor for (e, inset) the optimised two-tier two-branch channel. (f) Corresponding inference errors due to extrinsic noise in the optimised two-tier two-branch channel. Note that the errors here are predominantly around the segment boundaries (x=0.25,0.5,0.75) and diminished compared to the one-tier two-branch channel in (d).

Robustness to intrinsic noise, with a different choice of objective function (Equation 23), in (a) one-tier two-branch (1T2B) channel and (b) two-tier two-branch (2T2B) channel architectures, previously optimised for extrinsic noise alone.

(c) A comparison of local inference errors of the two optimised channels in (a,b) in presence of intrinsic noise. Even for this choice of objective function, the two-tier channel shows consistently better performance.

Experimental verification in the Drosophila Wg signalling system

The phenomenology of the morphogen reading and processing of Wg in the wing imaginal disc of Drosophila melanogaster Hemalatha et al., 2016 suggests a one-to-one mapping to the two-tier two-branch channel defined above, thus providing an ideal experimental system for a realisation of the ideas presented here (Figure 14a, b).

Figure 14 with 1 supplement see all
Comparison of theoretical results with experimental observations on Wg signalling system in Drosophila wing imaginal disc.

(a) Schematic of the cellular processes involved in Wg signalling, showing the two endocytic routes for the receptors (see text for further description). (b) Two-tier two-branch channel architecture corresponding to the Wg signalling system. (c) Schematic describing the XY view of wing disc. The vertical brown stripe marks the Wg producing cells. Horizontal green stripes mark the regions in wing disc used for analysis. See Experimental Methods (Appendix 15) for more information. (d) Coefficient of variation (CV) of CAAX-GFP intensity profiles, expressed in wing discs, as a function of (normalized) distance from producing cells (n=4). (e) Coefficient of variation in the output of the optimised two-tier two-branch channel (blue), and upon perturbation (orange) via removal of the non-signalling branch, implemented by setting all rates in the non-signalling branch κ to zero. The optimised parameter values for the plot can be found in Table 2 under the column corresponding to nT=2,nB=2,r-=κC. (f) CV of intensity profiles of endocytosed Wg in control wing discs (C5GAL4Xw1118; blue; n=4) and discs where CLIC/GEEC endocytic pathway is removed using UAS-myr-garz-DN (C5GAL4XUAS-myr-garz-DN; orange; n=5).

Figure 14—source data 1

CV of intensity measurements of CAAX GFP and endocytosed Wg, in control and myr-Garz-DN, as a function of distance from producing cells in individual samples of wing imaginal discs.

https://cdn.elifesciences.org/articles/79257/elife-79257-fig14-data1-v1.xlsx

Wingless (Wg) is secreted by a line of cells (1–3 cells) at the dorso-ventral boundary and forms a concentration gradient across the receiving cells Neumann and Cohen, 1997. Receiving cells closer to the production domain show higher Wg signalling while those farther away have lower Wg signalling Neumann and Cohen, 1997. Several cell autonomous factors influence reading and processing of the morphogen Wg in the receiving cells. Binding of Wg to its signalling receptor, Frizzled-2 (DFz2), initiates signal transduction pathway and nuclear translocation of β-catenin which further results in activation of Wg target genes (reviewed in Clevers and Nusse, 2012). In addition to the signalling receptor, binding receptors such as Heparin Sulphate Proteoglycans (HSPGs) – Dally and Dlp also contribute to Wg signalling Baeg et al., 2001; Franch-Marro et al., 2005. Further, the two receptors follow distinct endocytic pathways Hemalatha et al., 2016: while, DFz2 enters cells via the Clathrin Mediated Endocytic pathway (CME), Wg also enters cells independent of DFz2, possibly by binding to HSPGs, through CLIC/GEEC (CG) endocytic pathway. The two types of vesicles, containing Wg bound to different receptors, merge in common early endosomes Hemalatha et al., 2016. However, only DFz2 receptors in their Wg-bound state, both at the cell surface and early endosomes, are capable of generating a downstream signal leading to positional inference through a transcriptional readout Tsuda et al., 1999. This phenomenology is faithfully recapitulated in our two-tier two-branch channel architecture (Figure 14b) in which DFz2 and HSPG receptors play the role of the two branches. The conjugated state ‘Q’ represents a combination of the readings from the two branches, possibly realised by the co-receptors HSPGs that bind Kirkpatrick et al., 2006; Capurro et al., 2008 and present Hemalatha et al., 2016 diffusible ligands to signalling receptors (either on the cell surface or within endosomes).

Since an experimental measurement of positional inference error poses difficulties, we measure the cell-to-cell variation in the signalling output for a given position x as a proxy for inference error (Appendix 4—figure 3). Larger the variation, higher is the inference error. This is calculated as coefficient of variation (CV, Appendix 15) in the output across cells in the y-direction (Figure 14c).

Let us first discuss the results from the theoretical analysis. The optimised two-tier two-branch channel (Figure 14b) shows that the magnitude and the fluctuations in the coefficients of variation are small, with a slight increase with position (blue, Figure 14f). This is consistent with the low inference error associated with the optimised channel (Figure 8b). Upon perturbing this channel via removal of the non-signalling branch, the magnitude and fluctuations in the signalling output variation increases significantly (orange, Figure 14e). This qualitative feature of the coefficient of variation in the optimised two -tier two-branch channel is replicated in the Wg measurements of wild type cells.

In the experiments, we first established the method by determining the CV of a uniformly distributed signal, CAAX-GFP (expressed using ubiquitin promoter), and observed that the CV of CAAX-GFP is relatively uniform in x, the distance from Wg producing cells (Figure 14d). In order to study the steady state distribution of Wg within a cell and within the endosomes, we performed a long endocytic pulse (1 hr) with fluorescently labelled antibody against Wg Hemalatha et al., 2016; Prabhakara et al., 2022. Following this, we estimated the CV of the Wg endocytic profile as a function of x (Figure 14f, and Figure 14—figure supplement 1).

We assessed the CV of endocytosed Wg under two conditions: one, where the endocytic pulse of Wg is captured by the two branches and two tiers (control condition), and another, where we disengage one of the tiers by inhibiting the second endocytic pathway using a genetically expressed dominant negative mutant of Garz, a key player in the CG endocytic pathway Gupta et al., 2009. This perturbation has little or no effect on the functioning of the CME or the levels of the surface receptors that are responsible for Wg endocytosis (Hemalatha et al., 2016; Prabhakara et al., 2022). As predicted by the theory (Figure 14e), CV in the control shows a slight increase with position (Figure 14f) with fluctuations about the mean profile being small. In the perturbed condition, with the CG endocytic pathway disengaged, we find the CV shows a steeper increase with x and has larger fluctuations about the mean profile.

In principle, the coefficient of variation of the output is affected by all the microscopic stochastic processes that intersect with Wg signalling network in the wing imaginal disc and in the ligand input. Therefore, one has to be careful about interpreting the changes in the coefficient of variation of the output, based on the such perturbation experiments. Notwithstanding, this qualitative agreement between theory and experiment is encouraging.

Discussion

In this paper, we have posed the problem of spatial patterning of cell fates in a developing tissue as a local, cell autonomous morphogenetic decoding that ensures precise inference of position, that is robust to extrinsic and intrinsic noise. We treat the cells as inference channels capable of reading and processing the morphogen input. We describe the architecture of the inference channels in terms of three elements: branches (number of receptor types), tiers (number of compartments) and feedbacks. We ask for properties of the inference channel architectures that allow for precision and robustness in the task of morphogenetic decoding of cellular position.

Key results

Taking an information theoretic and systems biology approach, we have addressed the issue of accurate and robust morphogenetic decoding of position. For convenience, we summarise our key results in a point-wise manner:

  1. The main result is that given a noisy morphogen input, cells in a developing tissue can achieve low inference error of their positions by deploying a more elaborate multi-branch multi-tier channel architecture with feedback control. This ensures a separation between the reading of the morphogen input and buffering against noise.

  2. Having a combination of signalling and non-signalling receptors in the channel can significantly improve the performance of cells in their positional decoding.

  3. For a monotonically decaying morphogen input, the signalling and non-signalling receptors exhibit spatially varying profiles with the signalling receptor increasing away from the source and the non-signalling receptor decreasing away from the source. This implies that the non-signalling receptor ‘reads’ the morphogen input, while the signalling receptor buffers against noise.

  4. The performance of the multi-tier multi-branch channels is enhanced by having a feedback from the signalling branch to the non-signalling branch. Along with control on the levels of signalling receptor, this inter-branch feedback provides buffering against extrinsic noise.

  5. Having a multi-tier architecture (cellular compartmentalisation) tempers the effects of intrinsic noise in the channel by stabilising fluctuations of the output at steady-state.

  6. The optimisation shows that the characteristics of the signalling receptor are tightly controlled, whereas those of the non-signalling receptor are flexible. This implies that the signalling receptor is specific whereas the non-signalling receptor is promiscuous.

  7. Analysis of the geometry of the fidelity landscape reveals that the channel parameters corresponding to feedback, binding rates and profile of the signalling receptor are stiff, while the rest of the channel parameters are sloppy (elaborated in Section ‘Geometry of the inference error landscape: implications for control’).

  8. The efficacy of inter-branch feedback control is enabled by having a conjugated state corresponding to a confluence of the signalling and non-signalling branches in a common compartment.

  9. Our analysis demonstrates how local, cell autonomous control can facilitate the optimisation of a tissue-level task, here morphogenetic decoding of cellular position.

Our theoretical predictions are compared with experimental observations from Wg morphogen system of Drosophila wing imaginal disc. We first show that Wg signalling in the experimental system is equivalent to a two-tier two-branch channel. In the experiments, we use signal-to-noise ratio (SNR) of the output as a proxy for robustness of inference. Perturbation of the architecture, i.e. removal of the non-signalling branch, results in reduction of SNR. In a forthcoming manuscript, we will provide a detailed verification of the predicted opposing receptor profiles.

Geometry of the inference error landscape: implications for control

We have explored the local geometry of the fidelity landscape around the optimum, and the global geometry of the low inference error states, by perturbing channel parameters and concentration profiles of the receptors.

The local geometry of the fidelity landscape is studied using the Fisher information metric. This shows that steepest variation in the inference error comes from moving along the feedback parameters while perturbations to other channel parameters produces only marginal changes. Further, we explore the global geometry using the spectrum of the Hessian of the inference error. We find that the topography of the low inference error landscape resembles a ravine or a deep valley, which is shallow along the several sloppy directions and steep along the few stiff directions, the latter being predominantly along the feedback parameters. This dimensional reduction appears to be a recurring feature of such high-dimensional optimisation Transtrum et al., 2015; Yadav et al., 2022.

Such a geometrical approach also provides insight on the differences between the signalling and the non-signalling receptors, which shows up in the extent to which they influence inference errors in the neighbourhood of the optimum. Slight changes in the signalling receptor away from the optimum lead to a sharp increase in inference error while similar changes in the non-signalling receptor do not affect the inference errors significantly. This gives rise to the notion of stiff and sloppy directions of control - with non-signalling receptor placed under sloppy control. In a context with multiple morphogen ligands setting up the different coordinate axes (e.g. Wg, Dpp and Hh in imaginal discs Lin, 2004; Lin and Perrimon, 2000), the non-specific receptor can potentially facilitate cross-talks between them. A sloppy control on non-specific receptor would allow for accommodation of robustness in the outcomes of the different morphogens. This could potentially be tested in experiments.

Future directions

We end our discussion with a list of tasks that we would like to take up in the future. First, the information processing framework established here is very general. Obvious extensions of our models, such as adding more branches, tiers and chemical states, will not lead to qualitatively new features. However, one may alter the objective function – for instance, in the case of short range morphogens like Nodal Liu et al., 2022, only the positions of certain regions (closer to the morphogen source) or cell fate boundaries need to be specified with any precision. To this end, we have analysed another objective function which partitions the tissue into cell identity segments. The qualitative features of the optimised channel architectures remain unaltered. Depending on the developmental context, one might explore other objective functions. This would be a task for a future investigation.

Next, our optimisation study ignores cellular costs due to compartmentalisation, additional receptors and implementation of feedback controls, and thus possible trade-offs between cellular economy and precision in inference. Nevertheless, the observation that addition of extra tiers beyond two provides only marginal improvements to inference, already suggests a balance between precision and cellular costs.

Third, our theoretical result that the optimised surface receptor profiles are either monotonically increasing or decreasing from the morphogen source, suggests that the surface receptor concentrations are spatially correlated across cells. Such correlations could have a mechanochemical basis, either via cell surface tension that could in turn affect internalisation rates Thottacherry et al., 2018 or inter-cellular communication through cell junction proteins Garcia et al., 2018 or from adaptive feedback mechanisms between the output and receptor concentrations Barkai and Leibler, 1997. We emphasize that in the current optimisation scheme, we have allowed the receptor concentrations to vary over the space of all monotonically increasing, decreasing or flat profiles, and have not encoded the positional information in the receptor profiles.

Finally, we have considered the morphogen ligand as an external input to the receiving cells, outside the cellular information processing channel. There is no feedback from the output to the receptors and thus no ‘sculpting’ of the morphogen ligand profile. Morphogen ligand profiles (e.g. Dpp Romanova-Michaelides et al., 2021) are set by the dynamics of morphogen production at the source, diffusion via transcytosis and luminal transport, and degradation via internalisation. These cellular processes are common to both the reading and processing modules in our channel architecture. This would suggest a dynamical coupling and feedback between reading and ligand internalisation, which naturally introduces closed-loop controls on the surface receptors and a concomitant sculpting of the morphogen profile.

Appendix 1

Setting up the dynamical equations

The cellular processes involved in a general two-branch channel architecture, described in Section "Quantitative models for cellular reading and processing" and Figure 4, with any number of tiers nT are: production of receptors, binding-unbinding of ligand with the receptors, intracellular transport between cellular compartments, conjugation and splitting in the final tier, and degradation of the chemical species involved. The mass action kinetics for these cellular processes are written as the following set of first-order ordinary differential equations (ODEs):

(24) tR0(1)=rp(rbL+rd(1))R0(1)+ruR1(1)tR1(1)=rbLR0(1)(ru+rd(1)+rI(1))R1(1)+rR(1)R1(2)tS0(1)=κp(κbL+κd(1))S0(1)+κuS1(1)tS1(1)=κbLS0(1)(κu+κd(1)+κI(1))S1(1)+κR(1)S1(2)tR1(2)=rI(1)R1(1)(rR(1)+rd(2)+rI(2))R1(2)+rR(2)R(3)tS1(2)=κI(1)S1(1)(κR(1)+κd(2)+κI(2))S1(2)κR(2)S(3)tR1(nT)=rI(nT1)R1(nT1)(rR(nT1)+rd(nT))R1(nT)κCR1(nT)S1(nT)+κSQ(nT)tS1(nT)=κI(nT1)S1(nT1)(κR(nT1)+κd(nT))S1(nT)κCR1(nT)S1(nT)+κSQ(nT)tQ(nT)=κCR1(nT)S1(nT)κSQ(nT)

Rj(k) denotes species in the first branch, associated with the signalling receptor, in tier k and in bound (j=1) or unbound (j=0) states. Likewise, Sj(k) denotes the same for the non-signalling receptor. Q(nT) denotes the conjugated species in the final tier. Chemical rates related to the first branch are denoted by r while those in the second branch are denoted by κ. The processes of production, binding, unbinding, internalisation, recycling, conjugation, splitting and degradation are represented by the subscripts p,b,u,I,R,C,S,d, respectively (Table 1). The superscript on degradation rate indicates the index of the tier in which the degradation takes place. The superscript over internalisation rate stands for the tier-index of the species internalised. Superscript over the recycling rate follows that of the corresponding internalisation rate.

As described in Section ‘Quantitative models for cellular reading and processing’, we consider local cell-autonomous control on the total cell surface receptor concentrations, in a manner akin to the open-loop control Stengel, 1994, such that the sum of free and bound signalling receptors is ψ(x) and that of the non-signalling receptor is ϕ(x). This implies that the open-loop control actuated on the production (or secretion) of the cell surface receptors balances the net flux of receptors away from the surface. Consider the dynamics of R0(1) and R1(1) in Eq. A, with now an open-loop control realised through rp(t) on R0(1).

tR0(1)=rp(t)(rbL+rd(1))R0(1)+ruR1(1)tR1(1)=rbLR0(1)(ru+rd(1)+rI(1))R1(1)+rR(1)R1(2)

At any position x, if the total receptor concentration R0(1)(x,t)+R1(1)(x,t)=ψ(x) are to stay constant in time (through a chemostat), rp(t) must satisfy

rp(t)=rd(1)(R0(1)+R1(1))+rI(1)R1(1)rR(1)R1(2)

We emphasize that here we do not consider a closed-loop version of the receptor control. However, it can be realised through an integral control Goodwin et al., 2001, also known as perfect adaptation Barkai and Leibler, 1997. With open-loop-like control on total surface receptor concentrations (Section ‘Quantitative models for cellular reading and processing’), R0(1) and S0(1) are replaced by ψ-R1(1) and ϕ-S1(1), respectively. We then drop the subscripts from R1(k) and S1(k) to simplify notation.

(25) tR(1)=rbL(ψ(x)R(1))(ru+rd(1)+rI(1))R(1)+rR(1)R(2)tS(1)=κbL(ϕ(x)S(1))(κu+κd(1)+κI(1))S(1)+κR(1)S(2)tR(nT)=rI(nT1)R(nT1)(rR(nT1)+rd(nT))R(nT)κCR(nT)S(nT)+κSQ(nT)tS(nT)=κI(nT1)S(nT1)(κR(nT1)+κd(nT))S(nT)κCR(nT)S(nT)+κSQ(nT)tQ(nT)=κCR(nT)S(nT)κSQ(nT)

To keep the notation simple, we avoid explicitly indicating superscripts over chemical rate symbols when presenting the dynamical equations in the main text (Equations 10–17). In the numerical analysis, however, the degradation, internalisation and recycling rates of different tiers are treated separately. The steady-state solutions of the above equations, without feedback controls on chemical rates, are given by,

One-tier channel:

(26) R(1)=rbLrbL+ru+rd(1)ψS(1)=κbLκbL+κu+κd(1)ϕQ(1)=κCκSR(1)S(1)

Two-tier channel:

(27) R(1)=rbLrbL+ru+rd(1)+rd(2)rI(1)rR(1)+rd(2)ψS(1)=κbLκbL+κu+κd(1)+κd(2)κI(1)κR(1)+κd(2)ϕR(2)=rI(1)rR(1)+rd(2)R(1)S(2)=κI(1)κR(1)+κd(2)S(1)Q(2)=κCκSR(2)S(2)

For channels with nT>2, the solutions are written in the following recursive form Channel with nT tiers:

(28) R(1)=rbLrbL+ru+r~d(nT)ψR(k)=rI(k1)rR(k1)+r~d(nTk+1)R(k1)fork=2,3,...,nTS(1)=κbLκbL+κu+κ~d(nT)ϕS(k)=κI(k1)κR(k1)+κ~d(nTk+1)S(k1)fork=2,3,...,nTQ(nT)=κCκSR(nT)S(nT)

where r~d(k) are in turn evaluated recursively as

r~d(1)=rd(n)r~d(nTk+1)=rd(k)(1+rI(k)rR(k)+r~d(k+1)r~d(k+1)rd(k))fork=nT1,nT2,...,1κ~d(1)=κd(n)κ~d(nTk+1)=κd(k)(1+κI(k)κR(k)+κ~d(k+1)κ~d(k+1)κd(k))fork=nT1,nT2,...,1

The general form of Equations 26–28 holds upon introduction of inter-branch feedbacks (refer Figure 5) and therefore the set of ODEs (Equation 25) can be solved analytically by appropriately replacing the rates under feedback with the Hill-form functions discussed in Section ‘Quantitative models for cellular reading and processing’. Cases with other types of feedbacks need to be solved numerically (Appendix 3).

In writing these dynamical equations, we have made two assumptions, simply as a matter of convenience. One may note that only the ligand bound states of receptors are transported between tiers. This is done with the consideration that residence times of unbound receptors within a compartment, other than the cell surface (first tier), is very small that is receptors in enclosed compartments re-bind the ligand quickly after unbinding due to small volume of the compartment. Further, the open loop control on receptors need not be strictly on surface receptors (tier 1). In principle, the control could be on the production rate of receptors. This does not change the solution in any substantial manner. For instance, in the case of two tiers, the solution would be

(29) R1(1)=rbLrbL(1+rI(1)rR(1)+rd(2))+ru+rd(1)+rd(2)rI(1)rR(1)+rd(2)rprd(1)
(30) S1(1)=κbLκbL(1+κI(1)κR(1)+κd(2))+κu+κd(1)+κd(2)κI(1)κR(1)+κd(2)κpκd(1)

and the rest of the species would follow hierarchically as in the case with surface receptor control (see Equation 28). As long as the deviations in the denominator,rdrIrR+rd and κdκIκR+κd remain small compared to ru, and ψrprd(1),ϕκpκd(1), the steady solutions for the two cases (control on surface receptors versus total receptors) will not differ significantly.

Appendix 2

Heuristics and choice of feedback controls

The figure below illustrates the effect of feedback parameters in Equations 18; 19 discussed under Section ‘Quantitative models for cellular reading and processing’.

Appendix 2—figure 1
Hill functions representing positive (left) and negative (right) feedback control on chemical rates r actuated by chemical species R.

The effect of feedback amplification α, sensitivity γ and strength n are indicated. r0 denotes the reference value of the chemical rate r in absence of feedback.

Considering the choices of chemical rate under control r{rI,rR,rd,κI,κR,κd} and the actuating nodes R{R(1),R(2),}, there are a large number of possible feedbacks for any given channel architecture. However, categorising the different feedbacks according to their effect can reduce the redundancies and aid numerical analysis. The steady-state output of a channel without the feedback controls on chemical rates is instructive in this categorisation. For instance, consider the output θ for a two-tier two-branch channel with nT=2.

(31) θ=wnTQ(nT)+k=1nTwkR(k)=wR(1)(1+1wwrIrR+rd(1+κCκSκIκR+κdS(1)))

This form of the output suggests that feedbacks from R(1) (term outside the bracket) on any of the rates in rIrR+rd and κCκSκIκR+κd could cross-correlate R(1) with the other signalling species R(2) and Q(2). Thus, a negative feedback on rI (κI) would be qualitatively equivalent to a positive feedback on rR (κR). Beside helping to reduce the number of CRNs that need to be considered, this observation also allows one to forego the numerical analysis of a CRN with feedback on rR from R(1) that may have non-unique solutions. Similarly, feedbacks from different but related actuators could be interchanged. For example, the feedback on κI from R(2), i.e. κIκI1+(γR(2))n, could be written in terms of the feedback on κI from R(1), i.e. κIκI1+(γR(1))n, by a simple transformation on feedback sensitivity γ=rIrR+rdγ. Note that Equation 31 also indicates the preferred direction of feedbacks, that is with the signalling receptors as actuators.

Appendix 3

A note on numerical methods

Optimisation over the channel parameter vector vV, belonging to a parameter space V, is highly non-convex (see Section ‘Performance of the Channel Architectures’). Moreover, the derivative of σ¯X with respect to v is ill-defined due to finite sampling of input distributions. With these considerations, we chose a derivative-free (gradient independent) optimisation algorithm viz. Pattern Search algorithm Hooke and Jeeves, 1961. The implementation of this algorithm in MATLAB can be found as patternsearch function under the Global Optimization Toolbox.

With this pattern search algorithm, we perform the optimisation routine as described in Section ‘Performance of the Channel Architectures’ in two rounds. First, with a certain number of initial points ninit=32 in the parameter space V, we determine the advantageous feedback topologies (CRNs) that is, that give lower inference errors. In the next round, we go through the same optimisation routine for the advantageous feedback topologies determined through the the first round, but now with ninit=320 to converge to the true global minimum. Using parallel computing clusters, the cost of this computation in terms of real running time can be brought down to 12–48 hr.

Ideally, this optimisation ought to be done with extrinsic and intrinsic noises considered together. However, we restrict optimisation to the case of extrinsic noise as the computational cost associated with steady-state solutions of chemical master equations (CMEs), in the case of intrinsic noise, is rather high. Therefore, we evaluate the inference error due to intrinsic noise only of those CRNs that were optimised in the context of extrinsic noise previously. For this, we solve the CMEs using adaptive explicit-implicit tau leaping algorithm Sandmann, 2009 to determine the steady-state outputs and thus inference error. Implementation of this algorithm requires the definition of some numerical parameters that (i) help switch between Gillespie and implicit/explicit tau-leaping algorithms (na,nd), (ii) decide the number of reactions when Gillespie algorithm is selected (nb), and (iii) various threshold (nc,ϵ). The values chosen for these are na=10, nb=10 and nb=100 if implicit tau-leaping algorithm was used in the previous step, nc=10, nd=100, δ=0.05 and ϵ=0.1.

Appendix 4

Robustness, sensitivity, and trade-offs

Precision in positional inference requires both that the output variance at a given position be small (i.e. the output is robust to the noise in input) and that the mean output at two neighbouring positions be sufficiently different (the output is sensitive to systematic changes in mean input).

Appendix 4—figure 1
Schematic of an ideal output profile.

The monotonically decreasing curve is the mean output profile. Due to various sources of noise, the output at any point will have a distribution around the local mean. If neighbouring cohorts of cells are to accurately distinguish their inferred positions from their outputs (in the Bayesian sense), the output distributions must have a small overlap (shaded region under the distribution curves).

With this heuristic understanding of an ideal output, we define two local measures: (i) χ measures the noise in the output of a cell cohort as compared to the noise in the received input. We define it as the ratio of coefficient of variation in the output θ to the coefficient of variation in the input L for the same cohort of cells.

(32) χCVθCVL.

Thus, robustness of the output increases with decreasing values of χ; (ii) ξ measures the sensitivity of the output to a systematic change in the input. We define it as the ratio of relative change in the mean output θ to the relative change in the mean input μL for the same cell cohort.

(33) ξΔθθ/ΔμLμL

The angular brackets in the above equation denote averaging over cells belonging to the same cohort. Thus, higher sensitivity implies higher values of |ξ|. Precision in positional inference, in terms of these two local measures, implies simultaneous minimisation of χ and |ξ|-1.

We calculate χ and ξ at all positions x in the various optimised channels. We find that addition of a branch allows for better robustness and alleviates the tension between sensitivity to systematic change in mean input and robustness to input noise (Appendix 4—figure 2).

Appendix 4—figure 2
Robustness-Sensitivity plots for nx-1=100 cohorts in the optimised channel architectures.

The robustness-sensitivity objective is to reach the origin along both the axes. Optimised single-branch architectures (labelled 1T1B and 2T1B) show the two measures as conflicting objectives, such that improvement in robustness is achieved at the expense of sensitivity beyond a certain point. This conflict is absent in the optimised two-branch architectures (labelled 1T2B and 2T2B). An additional tier in two-branch architectures can further improve the two local measures. Note: the choice of coordinate axes reflects the requirement of simultaneous minimisation of χ (Equation 32) and ξ-1 (Equation 33) and the allowed tolerance in ξ-1 when χ is below a certain level.

Furthermore, a one-tier two-branch channel with an inter-branch feedback control (Figure 7a of the main text) shows a preference towards certain concentrations of the two receptors, which provide the cellular output robustness to input noise, i.e. minimise χ. To maintain lower values of χ (higher robustness), the signalling receptors ψ must decrease with increasing mean of the input μL (Appendix 4—figure 3). Increasing the non-signalling receptors ϕ as a function of mean input μL helps separate the mean outputs θ at neighbouring positions and thus aids in increasing the sensitivity (Appendix 4—figure 4). This is consistent with the receptor profiles in the optimised one-tier two-branch channel (Figure 7b of the main text).

Appendix 4—figure 3
Contour plots of χ (see Eq.Equation 32) in the optimised one-tier two-branch channel with inter-branch feedback control (Figure 7a of the main text) showing the preferred receptor combinations (deep blue) for different values of mean input μL.

Green dots denote the receptor concentrations in the optimised channel (Figure 7b,inset of the main text) at positions corresponding to the values of mean input μL indicated above the contour plots. The optimised parameter values for the plots can be found in Table 2 under the column corresponding to nT=1,nB=2,r-=κC.

Appendix 4—figure 4
Contour plots of mean output θ in the optimised one-tier two-branch channel with inter-branch feedback control (Figure 7a of the main text).

The contours move downward along the axis of the non-signalling receptor ϕ. Therefore, as the preferred values of signalling receptor ψ decrease with mean ligand input μL (Appendix 4—figure 3), non-signalling receptor concentration ϕ needs to increase with μL to ensure that the mean output θ is a monotonically increasing function of.μL Green dots denote the receptor concentrations in the optimised channel (Figure 7b,inset of the main text) at positions corresponding to the values of mean input μL indicated above the contour plots. The optimised parameter values for the plots can be found in Table 2 under the column corresponding to.nT=1,nB=2,r-=κC

Appendix 5

Input-output relations in a minimal channel

Here, we try to understand why having an appropriate spatial gradient of receptors helps in separating the outputs in nearby cells, thus facilitating accurate positional inference. This is best done in the simplest case of a minimal channel that is with one tier and one branch. This channel ‘reads’ the ligand input through binding of the ligand on only one, signalling receptor.

Appendix 5—figure 1a shows that setting receptor conentrations to decrease with mean ligand input creates an unfavourable scenario. Low receptor availability for higher ligand concentrations ensures a saturation of the output at lower values (blue). On the other hand, higher receptor availability for ligand concentrations is impractical as the output remains low despite increase in potential saturation point (purple). This causes large overlaps of the outputs at neighbouring positions. On the other hand, Appendix 5—figure 1b shows that the output are better separated with receptor profiles that increase with the mean ligand input and thus would enable a better positional inference (Appendix 4).

Appendix 5—figure 1
Input-Output function of a minimal channel shows the importance of choosing the correct receptor profiles.

Appendix 6

Robustness due to inter-branch feedback

Here we describe the logic behind the working of an inter-branch feedback control. Consider the case of a one-tier two-branch channel with the feedback on conjugation rate κC (Figure 7a of the main text) actuated by the ligand bound state of the signalling receptor in the first tier R(1). The steady-state solution to Equations 10–12 corresponding to this channel with κCκC1+(γR(1))n is given by.

(34) R(1)=rbLrbL+ru+rdψS(1)=κbLκbL+κu+κdϕ
(35) Q(1)=κCκSrbLψrbL+ru+rdκbLϕκbL+κu+κd11+(γrbLrbL+ru+rdψ)n

Therefore, the output can be expressed in terms of L,ψ and ϕ as follows.

(36) θ=rbLψrbL+ru+rd(1+κCκSκbLϕκbL+κu+κd11+(γrbLrbL+ru+rdψ)n)

Appendix 6—figure 1 describes the behaviour of the signalling species R(1) and Q(1) as given by Equations 34; 35 where the parameters are taken from the optimised one-tier two-branch channel. R(1) increases monotonically with the ligand input L (blue curve, Appendix 6—figure 1) and saturates at a value set by the signalling receptor ψ. Meanwhile, the conjugate species Q(1) has a non-monotonic behaviour (yellow curve, Appendix 6—figure 1): for very low values of input, Q(1) rises sharply due to absence of the feedback effect from the small values of R(1), but then decreases with further increase in input as the value of the feedback actuator R(1) rises. This anti-correlated behaviour of R(1) and Q(1) due to the feedback results in the output θR(1)+Q(1) being a more stable function of the input for an intermediate range of the input (region around the cusp in the green curve, Appendix 6—figure 1). Modulating the signalling ψ and non-signalling ϕ receptors allows for placement of the stability region (cusp) in accordance with the range of ligand input received at any position x, thus tempering the noise in the output (see input-output relations in Figure 7d of the main text).

Appendix 6—figure 1
Concentrations of the signalling species R(1),Q(1) and total cellular output θ in the optimised one-tier two-branch channel with ψ=164nM,ϕ=318nM.

The optimised chemical rates and feedback parameters for the above plot can be found in Table 2 under the column corresponding to nT=1,nB=2,r-=κC.

Appendix 7

Heuristic for the influence of feedback action on receptor profiles

The optimised one-tier two-branch channel showed a monotonically increasing profile of the signalling receptor ψ, and a monotonically decreasing profile of the non-signalling receptor ϕ (Figure 7b). Here, we provide an understanding for this counter-intuitive result using the forms of the channel output θ and variation in the output Lθ|ψ,ϕ.

Consider the output θ for a one-tier two-branch channel with inter-branch feedback on the conjugation rate κC. Equation 37 gives the explicit form for the output of this channel in terms of input ligand concentration L, signalling ψ and non-signalling ϕ receptor concentrations, feedback strength n, feedback sensitivity γ and binding rates rb,κb. The last three are small parameters in the optimised channel i.e. rb0.1,κb0.05,γ0.2. Note that the conjugation and splitting rates are absorbed into ϕ, i.e. κCκSϕϕ.

(37) θ=LψL+rb1[1+LϕL+κb111+(γLψL+rb1)n]

We now look at the output in the limit of high ligand input L2rb-1κb-1, and low ligand input L0.1rb-1. The output forms in these limits are given by Equation 38 and Equation 39.

(38) θhigh L=23ψ+13ϕψ1+(23γψ)n
(39) θlow L=rbψ+rbψκbϕ1+(rbγψ)n

To analyse the variation of channel output θ with ligand input L, we compute Lθ at fixed receptor concentrations ψ,ϕ.

(40) Lθ|(ψ,ϕ)=rb1ψ(L+rb1)2[1+LϕL+κb111+(γLψL+rb1)n]+ϕL+κb1LψL+rb111+(γLψL+rb1)n[κb1L+κb1+nrb1L+rb111+(γLψL+rb1)n]

In the limits of high and low ligand input, the form simplifies to

(41) Lθ|(ψ,ϕ)high L=rbψ9+118ψϕ1+(23γψ)n(rb+3κbrbn1+(23γψ)n)
(42) Lθ|(ψ,ϕ)low L=rbψ+rbψκbϕ1+(rbγψ)n(2n1+(rbγψ)n)

We provide numerical estimates for the output and the variation in the output, in the high and low ligand input limits, with all possible receptor profile combinations (Figure 3). As shown in the table below (highlighted in bold), it is optimal to have lower levels of the signalling receptor ψ5rb-1 and higher levels of the non-signalling receptor ϕ50rb-1 in the high L region (close to source). On the other hand, in the low L region (far from source), it is optimal to have higher levels of the signalling receptor ψ15rb-1 and lower levels of the non-signalling receptor ϕ15rb-1. These receptor combinations in the two limits on the input (highlighted rows in the table below) help maximally separate the output at the two ends of the tissue while keeping the variation in the output low at both ends. Taking the receptor profiles as monotonic functions of position, this would imply that for a one-tier two-branch channel with an inter-branch feedback, the optimum signalling receptor ψ has a monotonically increasing spatial profile while the optimum non-signalling receptor ϕ has a monotonically decreasing spatial profile. This qualitatively explains the result obtained for the optimised one-tier two-branch channel in Section ‘Branched architecture with multiple receptors provides accuracy and robustness to extrinsic noise’, Figure 7.

high L
low L
output θvariation Lθoutput θvariation Lθ
low ψ, low Φ293.950.5923.75165.93
high ψ, low Φ156.531.6722.521.12
low ψ, high Φ902.070.6967.5543.46
high ψ, high Φ288.431.6940.06–25.19

Appendix 8

Results of optimisation of two-tier two-branch channels

As discussed in Section ‘Branched architecture with multiple receptors provides accuracy and robustness to extrinsic noise’, additional tiers in a two-branch channel only marginally reduced the inference errors due to extrinsic noise when compared to the optimised one-tier two-branch channel (Figure 8). Here, we show that both the receptor profiles and the input-output relations of these two optimised two-tier two-branch channels are qualitatively similar (Appendix 8—figure 1).

Appendix 8—figure 1
Results of optimisation of (a,d) two-tier two-branch channels.

(b,e) The output profiles (with standard error in shaded region) and (insets) the corresponding optimised signalling (blue) and non-signalling (red) receptor profiles. (c,f) The input-output relations at selected positions x=0.25,0.5,0.75 are shown as solid lines, shaded with the same colour as the position-markers (coloured rectangles in b,e insets). The signalling ψ(x) and non-signalling ϕ receptor concentrations are mentioned in the legend. For a fixed distribution of ligand input (Equation 1), the range of input values recorded by the receptors at the selected positions gives rise to a range of outputs (circles). Tuning of input-output relations through receptor concentrations reduces output variance and minimises overlaps in the outputs of neighbouring cell cohorts. The optimised parameter values for the plots in (b–c,e–f) can be found in Table 2 under the column corresponding to nT=2,nB=2,r-=κI and nT=2,nB=2,r-=κC, respectively.

Appendix 9

Additional tiers dampen the effects of feedback to provide stability

Stochastic simulations of the chemical reaction network (CRN) corresponding to the optimised one-tier two-branch channel (Figure 9a of the main text) show large fluctuations in the time trajectories of signalling species Q(1) about its steady-state mean (Figure 9d–f of the main text). Q(1) increases due to absence of feedback from small values of R(1). However, beyond some amount of increase in Q(1), the trajectories veer back towards their mean values. The state of small R(1) and increasing Q(1) can be maintained by replenishment of R(1) from binding reaction followed by immediate conjugation with large pools of S(1) due to high availability of the non-signalling receptor ϕ. Such amplified fluctuations are absent in the two-tier two-branch channel (Figure 9g–i of the main text). Here, we provide heuristics of the differing behaviours in the optimised one-tier and two-tier channels by analysing a simpler set of CRNs (Appendix 9—figure 1) with the essential elements of these channels.

Appendix 9—figure 1
Two- and three-species CRNs with production k1 and degradation k2 rates of species 1 (s1) and inter-conversion k3,k4 rates between the “signalling" (output generating) species (in purple box).

These rates mimic the binding rb, unbinding ru and conjugation-splitting κC,κS rates respectively in the optimised one-tier two-branch and two-tier two-branch channels (Figure 9a and b of the main text). Consistent with this mapping, the feedback is from species 1 on k3. The three-species CRN has additional rates k3,k4 mimicking internalisation rI and recycling rR rates, respectively, in the optimised two-tier two-branch channel (Figure 9b of the main text). In both cases, output is the sum of the last two nodes in the purple box.

The dynamics for the first, two-species CRN is given by.

(43) s1˙=k1(cs1)(k2+k3(s1))s1+k4s2s2˙=k3(s1)s1k4s2

and that of the second, 3-species CRN by

(44) s1˙=k1(cs1)(k2+k3)s1+k4s2s2˙=k3s1(k4+k3(s1))s2+k4s3s3˙=k3(s1)s2k4s3

with constant c playing the role of ψ, thus providing an upper bound to species 1. In both these sets of equations, we consider a feedback k3 such that k3k01+(γs1)n where k0 represents the reference value of k3 in absence of feedback. Appendix 9—figure 2 shows the phase portrait for the dynamics of the 2-species CRN when k0k1 with a moderately strong feedback n=2. This is representative of κCS(1)rbL in the optimised one-tier two-branch channel (Figure 9a of the main text). The nullcline s˙1=0 (dashed, green curve) acts as a separatrix for the behaviour of this system: if due to fluctuations in species s1 the system (s1,s2) crosses the nullcline from its steady-state (pink point), it sets out on a large trajectory (black line) such that s1 remains close to zero while s2 grows fast until the system turns back towards the steady state. This is similar to the trajectories of R(1) and Q(1) discussed earlier. The non-linearity of the separatrix is due to the feedback from s1 on the rate k3 that couples to s1 in Equation 43. Higher production rates k1 (akin to higher values of ligand L in the optmised channel) bring the steady-state closer to the separatrix, making crossing of the separatrix due to fluctuations more likely. Having an additional node in between the actuator species and the controlled rate in the three-species CRN removes the non-linearity in s˙1=0 and provides buffering against this effect (Equation 44).

Appendix 9—figure 2
Phase portrait for Equation 43 with k0k1.

The red and green dashed lines are the nullclines s2˙=0 and s1˙=0, respectively. The pink dot denotes the steady-state solution (stable fixed point) of this system. The solid black line is a trajectory with initial point at the origin. Note that at steady-state, s1s2. If s1 drops beyond the green nullcline due to a fluctuation, the system goes back to the steady state through a long trajectory with s1 essentially remaining close to zero for a long period of time while s2 increases dramatically. Parameter values for the plot: k1=k2=k4=1,k3=100,c=50,γ=0.5,n=2.

Appendix 10

Two-tier two-branch channel with no feedback

Here we show that the two-tier two-branch channel without any inter-branch feedback control has fundamentally different optimisation characteristics and a poorer positional inference. Crucially, the optimised profiles of both signalling ψ and non-signalling ϕ receptors are monotonically decreasing away from the source (Appendix 10—figure 1b, inset).

Appendix 10—figure 1
Results of optimisation of (a) two-tier two-branch channel with no feedback on rates.

(b) The output profile (with standard error in shaded region) corresponding to the (inset) optimised signalling (blue) and non-signalling (red) receptor profiles. (c) The local inference error σX(x) is only marginally reduced throughout the tissue, when compared to the expected inference errors from ligand with no processing. The dashed line corresponds to a local inference error of one cell’s width 1/nx. (d) The input-output relations in this channel are monotonically increasing sigmoid functions saturating at only large values of input. The solid lines correspond to the input-output relations at selected positions, x=0.25,0.5,0.75 shaded with the same colour as the position-markers in (b inset, coloured rectangles). The signalling ψ(x) and non-signalling ϕ(x) receptor concentrations are mentioned in the legend. For a fixed distribution of ligand input (Equation 1), the range of input values recorded by the receptors at the selected positions gives rise to a range of outputs (circles). It is clear that neighbouring positions have significant overlaps in their outputs. The optimised parameter values for the plots in (b–d) can be found in Table 2 under the column corresponding to.nT=2,nB=2,r-={}.

Appendix 11

Uniform receptor profiles with uncorrelated noise

Note that in arriving at the optimised channel characteristics for a given morphogen profile, we go through all possible monotonic receptor profiles, including flat profiles. The optimised receptor profiles show a spatial gradient. Here, we ask why can’t a flat receptor profile (possibly modified by noise) infer positions accurately from a noisy morphogen gradient? Following the arguments in Appendix 4, we reason that if the morphogen gradient was not corrupted by noise, then flat receptor profiles would have sufficed to infer positions accurately. It is because one wants to discriminate between morphogen concentrations in neighbouring cells in a noisy background that there is a need for a spatial variation in the receptor profiles.

To demonstrate this, we consider uniform spatial profiles, with or without uncorrelated noise, for both the signalling and non-signalling receptors (Appendix 11—figure 1b, c), and optimise the rates and feedback parameters anew (Table 3) to show that this leads to a higher inference error compared to the optimal (Appendix 11—figure 1a). In fact, the inference error in these cases, even with an inter-branch feedback, is only marginally smaller than a channel with no processing of the ligand (black dots in Appendix 11—figure 1a). The inference from flat receptor profiles reflects the noise in morphogen gradient itself. This provides the motivation for choosing monotonically increasing or decreasing profiles for both the signalling and non-signalling receptors (Equations 8; 9). Note that this implicitly assumes spatial correlations in the surface concentrations of receptors across the cells.

Appendix 11—figure 1
Behaviour of inference error on varying receptor profiles.

(a) Inference error profiles due to extrinsic noise in the optimised two-tier two-branch channels with optimal receptor profiles (blue), uniform receptor profiles (red) and uniform receptor profiles with uncorrelated noise (orange). Having uniform receptor profiles simply reflects the noise in the ligand input (black). (b) Signalling receptor profiles corresponding to the cases in (a). (c) Non-signalling receptor profiles corresponding to the cases in (a). In (b,c), the mean uniform receptor concentration is set to the mid point of the optimised receptor profile while the strength of the uncorrelated noise is 0.1 times the mean. The chemical rates, receptor profile parameters and feedback parameters for the optimised two-tier two-branch channel can be found in Table 2 under the column corresponding to nT=2,nB=2,r-=κI. The optimised chemical rates and feedback parameters for the two-tier two-branch channels with uniform receptors can be found in Table 3.

Appendix 12

Dependence of inference errors on input characteristics

The general qualitative features of the optimised channels remain invariant to changes in input characteristics. We find the same feedback topology and qualitative results when the one-tier two-branch channel architecture is optimised for input distributions with different decay lengths of mean ligand input λ (Appendix 12—figure 1). Additionally, lowering noise in the ligand input reduces the inference error of the optimised channel and extends the region of robustness in the tissue (Appendix 12—figure 2).

Appendix 12—figure 1
Optimisation of one-tier two-branch channels for extrinsic noise with varying mean input decay lengths λ (Equation 2).

(a) The channel architecture with inter-branch feedback shows the lowest inference error for all values of λ considered. (b) The minimum local and average inference errors decrease with λ. (c) Optimised profiles of the signalling receptors are increasing functions of x for the different values of λ. (d) Optimised profiles of the non-signalling receptors are decreasing functions of x for the different values of λ.

Appendix 12—figure 2
The one-tier two-branch channel optimised for ligand distribution with standard deviation equal to the square root of mean, σL=μL, and decay length λ=0.3 shows smaller inference error for lower levels of input noise.

The dashed line corresponds to a local inference error of one cell’s width 1/nx. Note that the point at which local inference error departs away from 1% (one cell width error) extends further away from the source.

Appendix 13

Response in inference error due to perturbations in receptor concentrations

The definition of cellular output θ in a branched architecture involves making a distinction between the signalling and non-signalling receptors. In conjunction, the direction of the inter-branch feedback is from the signalling to the non-signalling branch (Figure 8b of the main text). This gives rise to the possibility of an asymmetric response of the optimised channel to perturbations in the two receptors around the optimal point. Appendix 13—figure 1a shows the average inference error σ¯X due to the receptor profiles in Appendix 13—figure 1b, c resulting from perturbations in the receptor control parameters A2 and B2 (Table 1, Equations 8; 9) around their optimal values. Each perturbed receptor profile (black curves in Appendix 13—figure 1b, c) leads to a net deviation Δψ,Δϕ from the optimal receptor profile (blue curves in Appendix 13—figure 1b, c), which is computed as follows.

(45) Δψ=dx(ψ(x;A2+δA2)-ψ(x;A2))
(46) Δϕ=dx(ϕ(x;B2+δB2)-ϕ(x;B2))

where A2 and B2 are the optimum values, and δA2,δB2 are the perturbations in the receptor control parameters. This is simply the signed area between the optimised and perturbed receptor profiles. Note that the perturbed receptor profiles are such that they maintain the nature of monotnonicity. The local curvature of σ¯X(Δψ,Δϕ) around the optimal point (Δψ=Δϕ=0, red point in Appendix 13—figure 1a) in the ψ-ϕ plane has eigenvalues λ10.016,λ26.1 corresponding to the eigenvectors that are nearly parallel to Δϕ and Δψ axes, respectively (white arrows in Appendix 13—figure 1a). This indicates that the inference error is much more sensitive to changes in the signalling receptor ψ than to changes in the non-signalling receptor ϕ, implying a stiff direction of control along the former and a sloppy direction of control along the latter.

Appendix 13—figure 1
Response of average inference error in the optimised two-tier two-branch channel (Figure 8b of the main text) to changes in receptor profiles shows the stiff-sloppy directions of control on receptors.

(a) Contours of average inference error σ¯X as functions of the net deviation from the optimal receptor profiles, as defined in Equations 45; 46. The white arrows indicate the directions of eigenvectors of the local Hessian (curvature) around the optimum (red point). The shorter arrow corresponds to the smaller eigenvalue (sloppy direction) while the longer arrow corresponds to the larger eigenvalue (stiff direction). (b,c) The allowed perturbations in receptor profiles, ψ(x) and ϕ(x) (black) around the optimal receptor profiles (blue), maintaining the nature of monotonicity. The optimised parameter values for the plots can be found in Table 2 under the column corresponding to nT=2,nB=2,r-=κI.

Appendix 14

Distribution of optimum channel parameters

In Section "Geometry of fidelity landscape", we commented on the nature of low inference error landscape as defined by optimum channel parameters that yield an average inference error σ¯X2%. Of the 16 channel parameters, we showed six parameters corresponding to ligand binding rates to the signalling and non-signalling receptors, conjugation and splitting rates, and feedback sensitivity and feedback strength. These parameters were stiff, that is small changes in these parameters led to strong variantions in the inference error. For completeness, here we present the frequency distributions of all the optimum channel parameters that yield an inference error of σ¯X2%. While some of these parameters are narrowly distributed about the upper or lower bounds of their permissible ranges, others are more broadly distributed across the range.

Appendix 14—figure 1
Frequency distributions of optimum channel parameters yielding σ¯X2%.

Symbols below each panel represent channel parameters listed in Table 1.

Appendix 15

Experimental methods

Fly stocks, endocytic assays, and imaging

Fly stocks used in this study are wildtype w1118, c5-GAL4, CAAX-GFP (Kyoto DGRC - 109823), Port et al., 2014 and UAS-myr-Garz-E740K. Flies were reared on corn flour medium containing sugar, yeast and agar along with antibacterial and antifungal agents. Flies were grown in 25°C incubators with 12 hr light/dark cycles for experiments and otherwise maintained at 18°C or 22°C. Third instar larval wing discs were dissected in Grace’s live imaging media Dye et al., 2017. Dissected discs were incubated with labelled Wg antibodies (Wg-AF-568) on ice for 45 min. Discs were transferred to room temperature for indicated time of pulse, washed with ice cold 1XPBS buffer and fixed using 4% PFA (5 min on ice +15 min at room temperature). Discs were then mounted in imaging chambers and imaged on a FV3000 laser-scanning confocal microscope using a 60 X/1.42 NA oil objective (with acquisition XY pixel dimensions of 0.138μm and Z stacks of size 0.5μm).

Analysis

The slice-by-slice images of the dome shaped wing disc from the confocal microscope were transformed into images from the outer most surface of the wing disc. This allowed us to compare the intensities of Wg from similar apico-basal height of cells across the dome-shaped epithelial. Wg production plane (Plane Q) and a perpendicular plane (Plane R) were defined. Intensity of different probes in curved tissues is affected by sample geometry and imaging depth. A data-based correction matrix was constructed using a uniform marker – CAAX GFP expressed uniformly under a ubiquitin promoter. Intensity for each disc, for each probe, was corrected using this data-based correction matrix. Detailed experimental and analysis methodology for extracting gradients from a curved tissue is described in Prabhakara et al., 2022. For computing the coefficient of variation, 18 bins parallel to Plane R were defined and intensities at different distances from the production plane was computed (schematic in Figure 14c of the main text). Mean and standard deviation (SD) for each disc across multiple parallel bins at different distances from the producing cells x were used to estimate the coefficient of variation.

(47) CV(x)=SD(x)Mean(x)

Computation of CV was done using intensity collected from the apical region of cells (20% of the entire length of cells). Normalized distance from the production plane is represented in all plots. Here, we have considered dorsal and ventral gradients to be equivalent.

Data availability

Figure 14 and Figure 14—figure supplement 1 contain source data in xlsx format. Modelling code and numerical data is available on Gitlab (https://gitlab.com/eoskrish/morphogenetic-decoding; copy archived at swh:1:rev:59c184aaa46c5e769a95ea39b48e911d0fc4fd5b).

References

  1. Book
    1. Goodwin GC
    2. Graebe SF
    3. Salgado ME
    (2001)
    Control System Design
    Prentice Hall Upper Saddle River.
  2. Book
    1. Lauffenburger DA
    2. Linderman JJ
    (1996)
    Receptors: Models for Binding, Trafficking, and Signaling
    Oxford University Press on Demand.
  3. Book
    1. MacKay DJ
    2. Mac Kay DJ
    (2003)
    Information Theory, Inference and Learning Algorithms
    Cambridge university press.
  4. Conference
    1. Sandmann W
    (2009) IEEE
    2009 Winter Simulation Conference.
    https://doi.org/10.5555/1995456
  5. Book
    1. Sengupta AM
    (2008)
    Modeling Biomolecular Networks: An Introduction to Systems Biology
    Oxford University Press.
  6. Book
    1. Stengel RF
    (1994)
    Optimal Control and Estimation
    Courier Corporation.

Decision letter

  1. Thomas Gregor
    Reviewing Editor; Princeton University, United States
  2. Naama Barkai
    Senior Editor; Weizmann Institute of Science, Israel
  3. Marcin Zagorski
    Reviewer; Jagiellonian University, Poland

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Cellular compartmentalisation and receptor promiscuity as a strategy for accurate and robust inference of position during morphogenesis" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor. The following individual involved in the review of your submission has agreed to reveal their identity: Marcin Zagorski (Reviewer #2).

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:

Both reviewers and reviewing editor agreed that no additional data is needed to support the major conclusions. Therefore a revision focused on improving clarity and presentation will be sufficient. The following two major concerns should be addressed; additional detailed comments and suggestions in the reviews below should be considered for clarity and overall readability.

1) One important point to clarify is whether the presented solutions (parameters identified) are representative of a larger class of sub-optimal or optimal solutions. This can be checked by perturbing the solutions. Presumably, the authors have investigated this point already as in the SM there are parts about, robustness, sensitivity and trade-offs, so improving presentation in the main part should solve that. This will also shed more light on the generality of the obtained results.

2) The important question to address is whether the choice of Eq. 5, which describes how the cell evaluates its position, influences other results in the manuscript.

Reviewer #1 (Recommendations for the authors):

– Extrinsic noise coming from the stochastic ligand production at the source is included in the model by letting cells along the direction parallel to the source boundary experience a spatially-fluctuating level of ligand, with statistics given in (1). From what we understand, this noise is quenched in space. Is it clear to what extent this is equivalent to the more realistic annealed disorder originating from Brownian motion and stochastic source production?

– In some cases, non-specific receptors (e.g. Dlp, Hanandal, Development, 2005) have been shown to increase in expression levels away from the morphogen source. Can the authors comment on this observation in light of their model?

– In appendix L, the local inference error for a (optimised) graded receptor expression is compared to that of a uniform receptor expression pattern. We find that the way this result is presented slightly misleading since the expression level in the uniform case is not optimised; the two setups should be compared after optimisation.

Reviewer #2 (Recommendations for the authors):

– In my opinion, Figure 2 is misleading and is disrupting the flow of the manuscript. Almost the same information is conveyed in Figure 4 and Figure 5. Further Figure 2 suggests some very regular arrangement of the nodes (regular topology of signalling architecture), which is not the case. Presenting a model with tiers and branches, or some different network-like schematic to indicate reaction and flux imbalances could improve presentation.

– In lines 85-87, different timescales of signalling processing are mentioned corresponding to branches (fast) and tiers (slow), but this aspect of regulation seems to be not discussed in the later parts of the paper. It might be worth drawing this analogy again when discussing how noise is integrated by different architectures. Possibly, there is the separation of the timescales or signalling integration takes place on the same timescales.

– Line 98, Figure 5 is mentioned just after Figure 2 (line 91). Please amend the ordering or rephrase so figures are referenced in the text as they appear in the paper.

– Line 114, I would avoid statements "of course any choice consistent with experimental observations would do". For instance, would the model work for stationary wave-like patterns of morphogens that could emerge via Turing mechanism (Green and Sharpe, Development 2015)? Basically, I would rephrase providing a description of acceptable ligand profiles (e.g. monotonically decreasing).

– The optimization procedure employed is computationally quite-consuming as the intrinsic noise is calculated by solving a chemical master equation. Can the model be solved without directly solving CME? What are the differences? If the differences with and without CME are small this might help to have faster optimization and hence more explicitly explore the space of available signalling architectures that result in optimal or close to the optimal solution.

– The intrinsic and extrinsic noises are varied separately as far as I understand. Are there any arguments/heuristics that would indicate the resulting global solutions would be the same if two types of noise would be varied simultaneously?

– In appendix C, the screening is described for the model parameters, but it is still slightly obscure what process was carried out.

– In table 2, in the middle column, we see that the weight of tier 2 is orders of magnitude greater than that for tier 1. In this case, what role does the first tier really play? And does such a strong weighting make sense biologically?

– I don't really understand the need for appendix E, but I might be missing something.

– In Appendix G, Figure 8, we see that adding a third tier does not reduce the minimum average inference error very much at all. Were any simulations done for four-tiered systems? What sort of computational cost does adding one extra tier there?

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

Author response

Essential revisions:

Both reviewers and reviewing editor agreed that no additional data is needed to support the major conclusions. Therefore a revision focused on improving clarity and presentation will be sufficient. The following two major concerns should be addressed; additional detailed comments and suggestions in the reviews below should be considered for clarity and overall readability.

1) One important point to clarify is whether the presented solutions (parameters identified) are representative of a larger class of sub-optimal or optimal solutions. This can be checked by perturbing the solutions. Presumably, the authors have investigated this point already as in the SM there are parts about, robustness, sensitivity and trade-offs, so improving presentation in the main part should solve that. This will also shed more light on the generality of the obtained results.

We thank the reviewers/editor for suggesting this. Indeed the original manuscript has a detailed discussion on sensitivity and robustness to variations in receptor concentrations (Section “Asymmetry in branched architecture: promiscuity of non-signalling receptors”, Appendix 13 and paragraph 6 of Discussion in the revised manuscript). We have now included an additional Subsection “Geometry of fidelity landscape” of Results, where we have described the geometry of the fidelity landscape (inference error) around the optimum and the low inference error states. This has involved additional computation and we have presented our detailed results on this as Figure 10 and Figure 11.

Briefly, we measure (i) percent change in inference error upon perturbation of the channel parameters, and (ii) the Fisher information metric (FIM) around the optimum. We use the optimised one-tier two-branch channel for this analysis which yields a minimum inference error of σ¯x1.9% i.e. approximately two cell widths. We see that the inference error does not change significantly with most parameters, i.e. it remains within σ¯x2.2% margin, except for the feedback strength n. The eigenspectrum of the FIM shows only a few stiff directions (eigenvectors with high eigenvalues), which have a significant component along the feedback parameters.

We now explicitly connect this detailed perturbative analysis in channel parameters to the discussion on perturbations of receptor profiles in Appendix 13.

To address the geometry of low inference error states, we look at the frequency distributions of channel parameters in a small cutoff region around the minimum identified, i.e. σ¯x2%. We find that the frequency distribution of optimal feedback parameters whose inference error lies within the range 1.9%σ¯x2% is narrowly distributed about the global optimum. Parameters corresponding to forward and backward rates are skewed towards the upper and lower bounds of the allowed parameter range respectively. All other optimal parameters corresponding to degradation rates, minimum and maximum receptor values and steepness of the receptor profiles, show a very broad spread over this range. Further, position vector of the optimum parameters lie predominantly along the sloppy directions of the Hessian evaluated at the global minimum. This suggests that geometry of the low inference error landscape resembles a narrow valley, which is shallow along several sloppy directions and steep along the few stiff directions.

This has now been included as a part of Subsection “Geometry of fidelity landscape” in the Results section along with new Figures 10 and 11 in the revised manuscript.

2) The important question to address is whether the choice of Eq. 5, which describes how the cell evaluates its position, influences other results in the manuscript.

We thank the reviewers for this comment. We presume here that the question refers to Eq. 6 since this describes the choice of objective function and how the cell evaluates its position. On the other hand, Eq. 5 refers to a standard inference of position from an output distribution computed by the cell. We discuss the assumptions behind Eq. 5 in our response to Reviewer #1’s first question.

Indeed the original manuscript contains a discussion on the choice of objective function i.e. inference error (last paragraph of Section “Mathematical Framework” and subsection “Future directions” of the Discussion in the revised manuscript). To highlight this important point, we have now included Subsection “ Choice of objective function” of Results on the choice of objective function. With the alternate objective function which specifies positions of certain regions or cell fate boundaries, we observe that the general, qualitative results remain the same, i.e. an additional branch helps reduce the inference errors due to extrinsic noise and an additional tier helps with the same when faced with intrinsic noise. Along with these results, we see that the receptor profiles too remain qualitatively the same as before.

Reviewer #1 (Recommendations for the authors):

– Extrinsic noise coming from the stochastic ligand production at the source is included in the model by letting cells along the direction parallel to the source boundary experience a spatially-fluctuating level of ligand, with statistics given in (1). From what we understand, this noise is quenched in space. Is it clear to what extent this is equivalent to the more realistic annealed disorder originating from Brownian motion and stochastic source production?

In our problem setup, each cellular output is calculated from the steady-state distribution of the morphogen concentration in its bound and unbound states within the cell. This implicitly assumes that the relaxation timescale of morphogen concentrations to its steady-state is faster than the timescales associated with individual cellular processing. This is a fair assumption as the timescales associated with diffusive transport (0.1 – 0.01s) and extracellular degradation (seconds to minutes) of the morphogen are faster than cellular processes which are in the order of minutes to tens of minutes.

These estimates follow from FCS measurements on the diffusion of free Wg (manuscript in preparation by C Prabhakara et al) and for free Dpp [10] in Drosophila wing discs.

– In some cases, non-specific receptors (e.g. Dlp, Han et al, Development, 2005) have been shown to increase in expression levels away from the morphogen source. Can the authors comment on this observation in light of their model?

We thank the reviewer for pointing out this important paper by Han et al. These authors demonstrate that there are two non-signalling receptors, Dally and Dlp, involved in the Wg signalling, one of which Dlp interacts strongly with Wg and facilitates its sculpting. As mentioned in response to Reviewer #1 (comment #4), in the current manuscript, we have treated the morphogen profiles as independent of the receptor dynamics and hence cannot address the issue of morphogen sculpting. We will take up this important issue in a separate study.

– In appendix L, the local inference error for a (optimised) graded receptor expression is compared to that of a uniform receptor expression pattern. We find that the way this result is presented slightly misleading since the expression level in the uniform case is not optimised; the two setups should be compared after optimisation.

We have indeed optimised the chemical rates and feedback parameters in the two-tier two-branch channels with spatially uniform receptor profiles. We have modified the corresponding appendix (Appendix 11 in the revised manuscript) and included Table 3 in the revised manuscript to emphasize this.

In arriving at the optimised channel characteristics for a given morphogen profile, we go through all possible monotonic receptor profiles including flat profiles. Thus we know that any deviation from the optimum will lead to inferior positional inferences. However, the question posed in Appendix 11 of the revised manuscript is the following: “why can’t a flat receptor profile (modified by uncorrelated noise) infer positions accurately from a noisy morphogen gradient?”. If the morphogen gradient was not corrupted by noise, then surely flat receptor profiles would have sufficed. It is because one wants to discriminate between morphogen concentrations in neighbouring cells in a noisy background that the you need a dynamic and spatially varying profile for the receptors. As seen from Appendix 11—figure 1 (of the revised manuscript), the inference from flat receptor profiles reflect the noise in morphogen gradient itself.

We elaborate on this point in Appendix 11 in the revised manuscript.

Reviewer #2 (Recommendations for the authors):

– In my opinion, Figure 2 is misleading and is disrupting the flow of the manuscript. Almost the same information is conveyed in Figure 4 and Figure 5. Further Figure 2 suggests some very regular arrangement of the nodes (regular topology of signalling architecture), which is not the case. Presenting a model with tiers and branches, or some different network-like schematic to indicate reaction and flux imbalances could improve presentation.

We agree with the reviewer on this point. Thanks to the suggestions, we have now modified the Figure 2 in the revised manuscript. We hope updated figure is more appealing and conveys the point better.

– In lines 85-87, different timescales of signalling processing are mentioned corresponding to branches (fast) and tiers (slow), but this aspect of regulation seems to be not discussed in the later parts of the paper. It might be worth drawing this analogy again when discussing how noise is integrated by different architectures. Possibly, there is the separation of the timescales or signalling integration takes place on the same timescales.

The timescale separation in signal processing is between the direction of chemical reactions (receptor state transitions) of the same receptor type i.e. within the same branch (fast, horizontal) and the direction of transport corresponding to tiers (slow, vertical). This is reflected in our choice of parameter values (listed in Table 1, consistent with experimental values) over which all our optimisation is carried out. The unbinding rate, which corresponds to transition in receptor state, is an order of magnitude faster than the internalisation and recycling rates, which correspond to intracellular transport processes.

We have added a line to this effect in the first paragraph of the Results section in the revised manuscript.

– Line 98, Figure 5 is mentioned just after Figure 2 (line 91). Please amend the ordering or rephrase so figures are referenced in the text as they appear in the paper.

Thank you. We have attended to this point in the revised manuscript.

– Line 114, I would avoid statements "of course any choice consistent with experimental observations would do". For instance, would the model work for stationary wave-like patterns of morphogens that could emerge via Turing mechanism (Green and Sharpe, Development 2015)? Basically, I would rephrase providing a description of acceptable ligand profiles (e.g. monotonically decreasing).

Indeed, the morphogen profiles applicable here belong to the family of monotonically behaving functions. We have attended to this point in the revised manuscript.

– The optimization procedure employed is computationally quite-consuming as the intrinsic noise is calculated by solving a chemical master equation. Can the model be solved without directly solving CME? What are the differences? If the differences with and without CME are small this might help to have faster optimization and hence more explicitly explore the space of available signalling architectures that result in optimal or close to the optimal solution.

Solving the CME is the most accurate approach which is applicable in any parameter regime or any network topology. In some conditions such as low (additive) noise, other approaches e.g. solving an equivalent Langevin equation, may solve for steady-state faster with minor deviations from the CME solution. However, there is no guarantee that the intrinsic noise is small in any given network topology and across the parameter range.

– The intrinsic and extrinsic noises are varied separately as far as I understand. Are there any arguments/heuristics that would indicate the resulting global solutions would be the same if two types of noise would be varied simultaneously?

Yes, the extrinsic and intrinsic noises have been considered separately in our work. As the reviewer would appreciate, doing an optimisation considering intrinsic and extrinsic noises simultaneously is computationally exorbitant. Instead, for a fixed number of tiers and branches, we first optimise the channel architecture subject to extrinsic noise alone over all possible rates, receptor profiles and feedback topologies. For each of these optimised channels with fixed number of tiers and branches (e.g. as given by Table 2), we calculate the inference errors in face of intrinsic noise. By doing so, we compare robustness of each such optimised channel’s performance to intrinsic noise. This is appropriate considering the timescale of variation in the extracellular morphogen concentrations is taken to be faster than the cellular processing timescales, as discussed in our response to Reviewer #1’s comment #6.

This is the scope of our exercise. Were we to address the different question of what is the optimum channel characteristics subject to both extrinsic and intrinsic noises then we would have had to perform a simultaneous optimisation.

– In appendix C, the screening is described for the model parameters, but it is still slightly obscure what process was carried out.

Thank you for the suggestion. We have now made minor edits to Appendix 3 in the revised manuscript for clarity.

– In table 2, in the middle column, we see that the weight of tier 2 is orders of magnitude greater than that for tier 1. In this case, what role does the first tier really play? And does such a strong weighting make sense biologically?

As discussed in the original manuscript (Section “Quantitative models for cellular reading and processing” of the revised manuscript), the weights associated with a tier correspond to the residence time of the signal transducing molecules (or secondary messengers) in that tier. We optimise over these weights along with the other channel parameters. A large weight to the second tier in an optimised channel would imply that signal transduction takes place predominantly from the second tier. The first tier, nevertheless, is still important in generating the characteristics of the second tier e.g. through actuating the feedback (or feedforward) control.

The Wg signalling in Drosophila wing disc provides a concrete example. Hemalatha et al., 2016 [4], Seto et al., 2006 [11] and Rives et al., 2006 [12] have shown that endocytosis is important for signalling. Removal of different components of the endocytic network abrogates both short range as well as long range Wg signalling. Intracellular Wg endosomes often appear to be co-localized with the downstream signalling mediator Dishevelled (Dsh). However, a quantitative study comparing the extent of localization of Dsh at the plasma membrane versus endosomal membrane is currently missing. Given the relevance of the endocytic network in Wg signalling and co-localization of downstream adaptor proteins to endosomes containing Wg-DFz2, we believe that these endosomes (corresponding to the second tier) act as signalling hubs to transduce positional information.

– I don't really understand the need for appendix E, but I might be missing something.

In Appendix 5, we are trying to understand why having an appropriate spatial gradient of receptors helps in separating the outputs in nearby cells, thus facilitating accurate positional inference. This is best done in the simplest case of a minimal channel i.e. with one tier and one branch.

We have now made minor edits to Appendix 5 in the revised manuscript to reflect the above.

– In Appendix G, Figure 8, we see that adding a third tier does not reduce the minimum average inference error very much at all. Were any simulations done for four-tiered systems? What sort of computational cost does adding one extra tier there?

The number of network topologies increases as a quadratic function of the number of tiers. Further, the number of parameters to be optimised would also increase linearly with the number of tiers, resulting in a net cubic increase in computation time. Practically, even with parallelisation, the number of computing hours (in real time) would increase significantly with additional tiers. From our analysis of two- and three-tier channels, we would predict that further addition of tiers provides a marginal advantage in reduction of inference errors while incurring additional cellular costs.

References

https://doi.org/10.7554/eLife.79257.sa2

Article and author information

Author details

  1. Krishnan S Iyer

    Simons Center for the Study of Living Machines, National Center for Biological Sciences - TIFR, Bangalore, India
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0930-5164
  2. Chaitra Prabhakara

    National Center for Biological Sciences - TIFR, Bangalore, India
    Contribution
    Resources, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
  3. Satyajit Mayor

    National Center for Biological Sciences - TIFR, Bangalore, India
    Contribution
    Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    mayor@ncbs.res.in
    Competing interests
    Reviewing editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9842-6963
  4. Madan Rao

    Simons Center for the Study of Living Machines, National Center for Biological Sciences - TIFR, Bangalore, India
    Contribution
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    madan@ncbs.res.in
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6210-6386

Funding

Department of Atomic Energy, Government of India (RTI4006)

  • Krishnan S Iyer
  • Chaitra Prabhakara
  • Satyajit Mayor
  • Madan Rao

Simons Foundation (287975)

  • Madan Rao

DBT-Wellcome Trust India Alliance (IA/M/15/1/502018)

  • Satyajit Mayor

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

Acknowledgements

We thank Thomas Lecuit for insights during the course of investigation. We thank past and present members of the Simons Centre, especially Alkesh Yadav, Amit Kumar, Archishman Raju, Kabir Husain, Mukund Thattai and Sandeep Krishna, for critical inputs. In particular, we thank Archishman Raju for useful discussions on the geometric analysis of the fidelity landscape. We acknowledge support from the Department of Atomic Energy (India), under project no. RTI4006, and the Simons Foundation (Grant No. 287975). MR and SM acknowledge DST (India) for JC Bose Fellowships. SM acknowledges a Margadarshi Fellowship of DBT-Wellcome Trust India alliance (IA/M/15/1/502018).

Senior Editor

  1. Naama Barkai, Weizmann Institute of Science, Israel

Reviewing Editor

  1. Thomas Gregor, Princeton University, United States

Reviewer

  1. Marcin Zagorski, Jagiellonian University, Poland

Version history

  1. Preprint posted: March 30, 2022 (view preprint)
  2. Received: April 5, 2022
  3. Accepted: January 14, 2023
  4. Version of Record published: March 6, 2023 (version 1)

Copyright

© 2023, Iyer et al.

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

  • 842
    Page views
  • 97
    Downloads
  • 0
    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. Krishnan S Iyer
  2. Chaitra Prabhakara
  3. Satyajit Mayor
  4. Madan Rao
(2023)
Cellular compartmentalisation and receptor promiscuity as a strategy for accurate and robust inference of position during morphogenesis
eLife 12:e79257.
https://doi.org/10.7554/eLife.79257

Share this article

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

Further reading

    1. Developmental Biology
    Marta Grzonka, Hisham Bazzi
    Research Article

    SAS‑6 (SASS6) is essential for centriole formation in human cells and other organisms but its function in mouse is unclear. Here, we report that Sass6‑mutant mouse embryos lack centrioles, activate the mitotic surveillance cell death pathway and arrest at mid‑gestation. In contrast, SAS‑6 is not required for centriole formation in mouse embryonic stem cells (mESCs), but is essential to maintain centriole architecture. Of note, centrioles appeared after just one day of culture of Sass6‑mutant blastocysts, from which mESCs are derived. Conversely, the number of cells with centrosomes is drastically decreased upon the exit from a mESC pluripotent state. At the mechanistic level, the activity of the master kinase in centriole formation, PLK4, associated with increased centriolar and centrosomal protein levels, endow mESCs with the robustness in using SAS‑6‑independent centriole-duplication pathways. Collectively, our data suggest a differential requirement for mouse SAS‑6 in centriole formation or integrity depending on PLK4 and centrosome composition.

    1. Developmental Biology
    2. Neuroscience
    Athina Keramidioti, Sandra Schneid ... Charles N David
    Research Article

    The Hydra nervous system is the paradigm of a ‘simple nerve net’. Nerve cells in Hydra, as in many cnidarian polyps, are organized in a nerve net extending throughout the body column. This nerve net is required for control of spontaneous behavior: elimination of nerve cells leads to polyps that do not move and are incapable of capturing and ingesting prey (Campbell, 1976). We have re-examined the structure of the Hydra nerve net by immunostaining fixed polyps with a novel antibody that stains all nerve cells in Hydra. Confocal imaging shows that there are two distinct nerve nets, one in the ectoderm and one in the endoderm, with the unexpected absence of nerve cells in the endoderm of the tentacles. The nerve nets in the ectoderm and endoderm do not contact each other. High-resolution TEM (transmission electron microscopy) and serial block face SEM (scanning electron microscopy) show that the nerve nets consist of bundles of parallel overlapping neurites. Results from transgenic lines show that neurite bundles include different neural circuits and hence that neurites in bundles require circuit-specific recognition. Nerve cell-specific innexins indicate that gap junctions can provide this specificity. The occurrence of bundles of neurites supports a model for continuous growth and differentiation of the nerve net by lateral addition of new nerve cells to the existing net. This model was confirmed by tracking newly differentiated nerve cells.