Cellular compartmentalisation and receptor promiscuity as a strategy for accurate and robust inference of position during morphogenesis
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.sa0Introduction
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 pairrule 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 selfenhanced 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, presteady state patterning Bergmann et al., 2007 and serial transcytosis Bollenbach et al., 2007.
Most cell signalling systems have regulatory mechanisms that finetune signalling by controlling ligandspecific receptor interactions Rogers and Schier, 2011. Ligands such as TGF $\beta $/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 multitiered, multibranched 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:
‘reading’ of the morphogen input by receptors on the cell surface,
‘processing’ by various cellular mechanisms such as receptor trafficking, secondary messengers, feedback control, and
‘inference’ of the cell’s position in the form of a transcriptional readout.
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 branchtier architecture of the cellular processing as a bipartite Markovian network/graph Hartich et al., 2014, with a fast direction (involving multiple branches) consisting of ligandbound and unbound states along with chemical state changes, and a slower direction (involving multiple tiers) consisting of intracellular transport, fission and fusion, characterised by energyutilising 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 nonsignalling 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).
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(Lx)$, 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 n_{x}, n_{y} 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 multitiered, multibranched architecture with feedbacks described previously, that reads a noisy input $L(x,y)$ (morphogen concentration) and produces an ‘output’ (biochemical ‘signal’) $\theta (x,y)$ that is also noisy. Here, we choose to construct the noisy morphogen profile in the following manner: for a given position $x\in [0,1]$, cells along the $y$direction see different amounts of ligand coming from the same input distribution$p(Lx)$,
characterised parametrically by a mean ${\mu}_{L}(x)$ and standard deviation ${\sigma}_{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 ${\mu}_{L}$ and standard deviation ${\sigma}_{L}$.
Alternatively, one could choose a different parametrisation consistent with experimental observations for a morphogen profile with a monotonically decaying mean. The values of $A,\lambda $ chosen for our analysis are listed in Table 1. The corresponding output distribution $p(\theta 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,
where $p(\theta )={\int}_{0}^{1}dxp(\theta 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, ${\sigma}_{X}(x)$. For each position $x$, the inferred position ${x}^{\ast}$ of cells along the $y$direction is taken to be the maximum a posteriori estimate,
where we use $\stackrel{~}{x}$ to differentiate from the true position $x$. From this, the local and average inference error can be computed.
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 shortrange 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’).
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(\theta )$ (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}^{\ast}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}^{\ast}\theta )$ and hence the inference error ${\overline{\sigma}}_{X}$, one needs to know the prior $p(x)$ and the inputoutput relation giving rise to the output distribution $p(\theta 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 inputoutput 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 nonsignalling receptors’. In multitier 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 openloop 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 inputoutput relation, calculated as a tierwise 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 openloop control on total (ligand bound plus unbound) surface availability of the signalling $\psi $ and nonsignalling $\varphi $ 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
The range of values for these parameters considered in the numerical analysis are listed in Table 1. Therefore, when considering $\psi (x)$ to be monotonically increasing in $x$, we parametrise it with ${f}_{A}$. It follows that in a onebranch channel, there are two possibilities: $\psi \in \{{f}_{A},{f}_{B}\}$ while in a twobranch channel, there are a total of four possibilities: $(\psi ,\varphi )\in \{{f}_{A},{f}_{B}\}\times \{{f}_{A},{f}_{B}\}$. This allows us to simulate the ‘reading’ step performed by the cells (see Figure 1b).
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 singletier 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 nonsignalling 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
The steadystate output $\theta $, defined as the sum of all the ligandbound signalling states, is given by $\theta ={R}^{(1)}+{Q}^{(1)}$. Note that to describe the 1branch system, we simply set all rates $\kappa $ to zero.
Dynamics of processing in a multitier channel
In a multitiered 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 ${n}_{T}$th tier as ${Q}^{({n}_{T})}$. The CRN for such a system with ${n}_{T}$ 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
The output, realised from all the ligandbound signalling states, now becomes $\theta ={w}_{{n}_{T}}{Q}^{({n}_{T})}+{\sum}_{k=1}^{{n}_{T}}{w}_{k}{R}^{(k)}$ at steady state with w_{k}, such that ${\sum}_{k}{w}_{k}=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 singletiered and multitiered 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, $L\sim p(Lx)$, and enters the equations as a source term. On the other hand, intrinsic noise is a consequence of copynumber 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\in \{{r}_{I},{\kappa}_{I},{\kappa}_{C},\mathrm{\dots}.\}$ that is under feedback control actuated by the node $R\in \{{R}^{(1)},{S}^{(1)},\mathrm{\dots}.\}$ is modelled as.
with r_{0} as the reference value of the chemical rate in the absence of feedback. The range of values for amplification $\alpha $, feedback sensitivity $\gamma $ 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.
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 $\overrightarrow{v}$ belonging to a parameter space $\mathcal{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(Lx)$ and a channel architecture under consideration, the optimisation can be stated as
and implemented by the following algorithm, the details of which are presented in Appendix 3.
Optimisation scheme 


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 singletiered channels with one and two branches. Such architectures are similar in design to the classic picture of ligandreceptor kinetics Lauffenburger and Linderman, 1996; Alberts et al., 2017, but also to the selfenhanced 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 onetier onebranch 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 (${f}_{A}$ in Figure 3a; ), the outputs for different input ranges overlap significantly. On the other hand, if the receptor concentrations increase with mean input (${f}_{B}$ 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 onetier onebranch 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 $\psi $ sets both the asymptote and the steepness of the inputoutput 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).
Including a nonsignalling receptor $\varphi $ 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 onetier onebranch case, an interbranch feedback control (Figure 7a) results in an inputoutput 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 abovementioned tension between robustness to input noise and sensitivity to systematic changes in the mean input to a considerable extent (see Appendix 4—figure 2).
As seen in Figure 7c, the twobranch architecture with interbranch 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 twobranch architecture. (i) The signalling and nonsignalling receptors present opposing optimal profiles – a consequence of the negative interbranch feedback. (ii) The optimal nonsignalling receptor decreases away from the source, indicating that the nonsignalling 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 onebranch architectures, it is the signalling receptor that does the reading and buffering.
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 twotier twobranch architectures, one with interbranch feedback on the internalisation rate of the nonsignalling receptors ${\kappa}_{I}$ and the other on the conjugation rate ${\kappa}_{C}$, that have comparable inference errors (Figure 8b, c). Both the receptor profiles and the inputoutput relations of these two optimised twotier twobranch channels are qualitatively similar (Appendix 8—figure 1).
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 signaltonoise ratio (SNR) Stoeger et al., 2016. Further, by making the output $\theta $ a multivariate 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 onetier twobranch and twotier twobranch 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 onetier onebranch channel in the presence of intrinsic noise, can be traced to the instabilities of steadystate trajectories of the two signalling species ${R}^{(1)}$ and ${Q}^{(1)}$ driven by the nonlinear 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 twotier twobranch architecture (Figure 9g–i), the fluctuations in the signalling species are more tempered, the interbranch 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.
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 interbranch 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 (nonsignalling) 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 nonsignalling 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 nonsignalling receptors
Before comparing the theoretical results with experiments, we comment on the implications for the cellular control of the signalling $\psi $ and nonsignalling $\varphi $ receptors. In the twobranch architecture, the symmetry between the signalling and nonsignalling receptors is broken by the interbranch feedback and the definition of output $\theta $, 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 ${\overline{\sigma}}_{X}$ in the $\psi \varphi $ plane around the optimal point. We compute the eigenvalues of the local curvature of ${\overline{\sigma}}_{X}(\mathrm{\Delta}\psi ,\mathrm{\Delta}\varphi )$ around the optimal point ($\mathrm{\Delta}\psi =\mathrm{\Delta}\varphi =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 $\psi $ and $\varphi $ axes, respectively. This implies that while the signalling receptor is under tight cellular control, the control on the nonsignalling receptor is allowed to be sloppy. A similar feature is observed in the contour plots for the robustness measure $\chi $ (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 $\psi $, while the control on $\varphi $ can be lax.
This sloppiness in the levels of nonsignalling 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 nonsignalling 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 nonsignalling receptor, being promiscuous, is nonspecific. This, as we see below, is the case with the Heparan sulfate proteoglycans (HSPGs) such as Dally and Dallylike protein (Dlp) that participate in the Wingless (Wg) and Decapentaplegic (Dpp) signalling networks Lin and Perrimon, 2000; RomanovaMichaelides 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 onetier twobranch channel (shown in Figure 7a with optimum channel parameters listed in Table 2, Table 3) in presence of extrinsic noise.
To address the geometry of the local fidelity landscape around the optimum, we compute (i) percent changes in inference error ${\overline{\sigma}}_{X}$ due to perturbations in channel parameters (Figure 10a), and (ii) the eigenspectrum of the Fisher information metric (FIM, Figure 10b). The FIM ${g}_{\mu \nu}$ is evaluated in the logparameter space as Transtrum et al., 2015.
where, $\overrightarrow{v}\in \mathcal{V}$ is the channel parameter vector, and ${x}_{i},{y}_{j}$ 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 ${\overline{\sigma}}_{X}\le 2.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 $\gamma ,n$. This implies that variation of the feedback parameters from the optimum would result in significant changes in the inferred positions. Perturbing conjugation ${\kappa}_{C}$ and splitting ${\kappa}_{S}$ rates simultaneously (see eigenvector 16) does not produce any notable change to the inferred positions (eigenvalue $\sim {10}^{13}$). Further, perturbations to channel parameters other than the feedback parameters (eigenvectors 7–16) produce marginal changes in inferred positions.
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 onetier twobranch channel architecture with 2^{16} spacefilling initial points in the 16dimensional parameter space of this architecture. We then define the low inference error states as those channel parameters ${\overrightarrow{v}}^{\text{opt}}$ that yield ${\overline{\sigma}}_{X}\le 2\%$. This cutoff, which equals $\lceil \frac{{\overline{\sigma}}_{X}}{0.01}\rceil $, corresponds to declaring as equivalent all the inference errors ${\overline{\sigma}}_{X}$ that lie between one and two cells’ widths. Consistent with the local analyses, we find that the frequency distribution of optimal feedback parameters $\gamma ,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 nonsignalling 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 nonsignalling receptors as described in Section ‘Asymmetry in branched architecture: promiscuity of nonsignalling 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 ${\overrightarrow{v}}^{\text{opt}}$ in the parameter space $\mathcal{V}$ along the eigenvectors of the Hessian of ${\overline{\sigma}}_{X}$, defined as
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 $\overrightarrow{v}}^{\phantom{\rule{thinmathspace}{0ex}}\text{opt}}\in \mathcal{V$ 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.
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 ${N}_{p}$ segmented cell identities as.
where $\delta $ and $\mathrm{\Theta}$ denote the Kroneckerdelta and Heavisidetheta functions respectively, ${\xi}_{i}$ is the position of the boundary between the ith and (i+1)th segments, and $g$ is a function that maps position (actual or inferred) to a segment, that is $g:[0,1]\to \{1,2,\mathrm{\dots},{N}_{p}\}$ with ${N}_{p}$ as the total number of segments.
We optimise onetier onebranch, onetier twobranch and twotier twobranch channel architectures for the inference error as defined in Equation 23 with ${N}_{p}=4$ and equally spaced boundaries located at positions ${\xi}_{1}=0.25,{\xi}_{2}=0.5,{\xi}_{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 twobranch 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.
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 onetoone mapping to the twotier twobranch channel defined above, thus providing an ideal experimental system for a realisation of the ideas presented here (Figure 14a, b).
Wingless (Wg) is secreted by a line of cells (1–3 cells) at the dorsoventral 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, Frizzled2 (DFz2), initiates signal transduction pathway and nuclear translocation of $\beta $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; FranchMarro 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 Wgbound 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 twotier twobranch 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 coreceptors 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 celltocell 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 twotier twobranch 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 nonsignalling 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 twobranch 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, CAAXGFP (expressed using ubiquitin promoter), and observed that the CV of CAAXGFP 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 pointwise manner:
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 multibranch multitier channel architecture with feedback control. This ensures a separation between the reading of the morphogen input and buffering against noise.
Having a combination of signalling and nonsignalling receptors in the channel can significantly improve the performance of cells in their positional decoding.
For a monotonically decaying morphogen input, the signalling and nonsignalling receptors exhibit spatially varying profiles with the signalling receptor increasing away from the source and the nonsignalling receptor decreasing away from the source. This implies that the nonsignalling receptor ‘reads’ the morphogen input, while the signalling receptor buffers against noise.
The performance of the multitier multibranch channels is enhanced by having a feedback from the signalling branch to the nonsignalling branch. Along with control on the levels of signalling receptor, this interbranch feedback provides buffering against extrinsic noise.
Having a multitier architecture (cellular compartmentalisation) tempers the effects of intrinsic noise in the channel by stabilising fluctuations of the output at steadystate.
The optimisation shows that the characteristics of the signalling receptor are tightly controlled, whereas those of the nonsignalling receptor are flexible. This implies that the signalling receptor is specific whereas the nonsignalling receptor is promiscuous.
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’).
The efficacy of interbranch feedback control is enabled by having a conjugated state corresponding to a confluence of the signalling and nonsignalling branches in a common compartment.
Our analysis demonstrates how local, cell autonomous control can facilitate the optimisation of a tissuelevel 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 twotier twobranch channel. In the experiments, we use signaltonoise ratio (SNR) of the output as a proxy for robustness of inference. Perturbation of the architecture, i.e. removal of the nonsignalling 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 highdimensional optimisation Transtrum et al., 2015; Yadav et al., 2022.
Such a geometrical approach also provides insight on the differences between the signalling and the nonsignalling 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 nonsignalling receptor do not affect the inference errors significantly. This gives rise to the notion of stiff and sloppy directions of control  with nonsignalling 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 nonspecific receptor can potentially facilitate crosstalks between them. A sloppy control on nonspecific 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 tradeoffs 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 intercellular 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 RomanovaMichaelides 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 closedloop 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 twobranch channel architecture, described in Section "Quantitative models for cellular reading and processing" and Figure 4, with any number of tiers ${n}_{T}$ are: production of receptors, bindingunbinding 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 firstorder ordinary differential equations (ODEs):
${R}_{j}^{(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, ${S}_{j}^{(k)}$ denotes the same for the nonsignalling receptor. ${Q}^{({n}_{T})}$ 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 $\kappa $. 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 tierindex 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 cellautonomous control on the total cell surface receptor concentrations, in a manner akin to the openloop control Stengel, 1994, such that the sum of free and bound signalling receptors is $\psi (x)$ and that of the nonsignalling receptor is $\varphi (x)$. This implies that the openloop 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 ${R}_{0}^{(1)}$ and ${R}_{1}^{(1)}$ in Eq. A, with now an openloop control realised through ${r}_{p}(t)$ on ${R}_{0}^{(1)}$.
At any position $x$, if the total receptor concentration ${R}_{0}^{(1)}(x,t)+{R}_{1}^{(1)}(x,t)=\psi (x)$ are to stay constant in time (through a chemostat), ${r}_{p}(t)$ must satisfy
We emphasize that here we do not consider a closedloop 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 openlooplike control on total surface receptor concentrations (Section ‘Quantitative models for cellular reading and processing’), ${R}_{0}^{(1)}$ and ${S}_{0}^{(1)}$ are replaced by $\psi {R}_{1}^{(1)}$ and $\varphi {S}_{1}^{(1)}$, respectively. We then drop the subscripts from ${R}_{1}^{(k)}$ and ${S}_{1}^{(k)}$ to simplify notation.
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 steadystate solutions of the above equations, without feedback controls on chemical rates, are given by,
Onetier channel:
Twotier channel:
For channels with ${n}_{T}>2$, the solutions are written in the following recursive form Channel with ${n}_{T}$ tiers:
where ${\stackrel{~}{r}}_{d}^{(k)}$ are in turn evaluated recursively as
The general form of Equations 26–28 holds upon introduction of interbranch 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 Hillform 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 rebind 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
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,$\frac{{r}_{d}{r}_{I}}{{r}_{R}+{r}_{d}}$ and $\frac{{\kappa}_{d}{\kappa}_{I}}{{\kappa}_{R}+{\kappa}_{d}}$ remain small compared to r_{u}, and $\psi \sim \frac{{r}_{p}}{{r}_{d}^{(1)}},\varphi \sim \frac{{\kappa}_{p}}{{\kappa}_{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’.
Considering the choices of chemical rate under control $r\in \{{r}_{I},{r}_{R},{r}_{d},{\kappa}_{I},{\kappa}_{R},{\kappa}_{d}\}$ and the actuating nodes $R\in \{{R}^{(1)},{R}^{(2)},\mathrm{\dots}\}$, 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 steadystate output of a channel without the feedback controls on chemical rates is instructive in this categorisation. For instance, consider the output $\theta $ for a twotier twobranch channel with ${n}_{T}=2$.
This form of the output suggests that feedbacks from ${R}^{(1)}$ (term outside the bracket) on any of the rates in $\frac{{r}_{I}}{{r}_{R}+{r}_{d}}$ and $\frac{{\kappa}_{C}}{{\kappa}_{S}}\frac{{\kappa}_{I}}{{\kappa}_{R}+{\kappa}_{d}}$ could crosscorrelate ${R}^{(1)}$ with the other signalling species ${R}^{(2)}$ and ${Q}^{(2)}$. Thus, a negative feedback on ${r}_{I}$ (${\kappa}_{I}$) would be qualitatively equivalent to a positive feedback on ${r}_{R}$ (${\kappa}_{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 ${r}_{R}$ from ${R}^{(1)}$ that may have nonunique solutions. Similarly, feedbacks from different but related actuators could be interchanged. For example, the feedback on ${\kappa}_{I}$ from ${R}^{(2)}$, i.e. ${\kappa}_{I}\to \frac{{\kappa}_{I}}{1+{({\gamma}^{\prime}{R}^{(2)})}^{n}}$, could be written in terms of the feedback on ${\kappa}_{I}$ from ${R}^{(1)}$, i.e. ${\kappa}_{I}\to \frac{{\kappa}_{I}}{1+{(\gamma {R}^{(1)})}^{n}}$, by a simple transformation on feedback sensitivity ${\gamma}^{\prime}=\frac{{r}_{I}}{{r}_{R}+{r}_{d}}\gamma $. 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 $\overrightarrow{v}\in \mathcal{V}$, belonging to a parameter space $\mathcal{V}$, is highly nonconvex (see Section ‘Performance of the Channel Architectures’). Moreover, the derivative of ${\overline{\sigma}}_{X}$ with respect to $\overrightarrow{v}$ is illdefined due to finite sampling of input distributions. With these considerations, we chose a derivativefree (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 ${n}_{init}=32$ in the parameter space $\mathcal{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 ${n}_{init}=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 steadystate 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 explicitimplicit tau leaping algorithm Sandmann, 2009 to determine the steadystate 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 tauleaping algorithms (${n}_{a},{n}_{d}$), (ii) decide the number of reactions when Gillespie algorithm is selected (n_{b}), and (iii) various threshold (${n}_{c},\u03f5$). The values chosen for these are ${n}_{a}=10$, ${n}_{b}=10$ and ${n}_{b}=100$ if implicit tauleaping algorithm was used in the previous step, ${n}_{c}=10$, ${n}_{d}=100$, $\delta =0.05$ and $\u03f5=0.1$.
Appendix 4
Robustness, sensitivity, and tradeoffs
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).
With this heuristic understanding of an ideal output, we define two local measures: (i) $\chi $ 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 $\theta $ to the coefficient of variation in the input $L$ for the same cohort of cells.
Thus, robustness of the output increases with decreasing values of $\chi $; (ii) $\xi $ 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 $\u27e8\theta \u27e9$ to the relative change in the mean input ${\mu}_{L}$ for the same cell cohort.
The angular brackets in the above equation denote averaging over cells belonging to the same cohort. Thus, higher sensitivity implies higher values of $\xi $. Precision in positional inference, in terms of these two local measures, implies simultaneous minimisation of $\chi $ and ${\xi }^{1}$.
We calculate $\chi $ and $\xi $ 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).
Furthermore, a onetier twobranch channel with an interbranch 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 $\chi $. To maintain lower values of $\chi $ (higher robustness), the signalling receptors $\psi $ must decrease with increasing mean of the input ${\mu}_{L}$ (Appendix 4—figure 3). Increasing the nonsignalling receptors $\varphi $ as a function of mean input ${\mu}_{L}$ helps separate the mean outputs $\u27e8\theta \u27e9$ at neighbouring positions and thus aids in increasing the sensitivity (Appendix 4—figure 4). This is consistent with the receptor profiles in the optimised onetier twobranch channel (Figure 7b of the main text).
Appendix 5
Inputoutput 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 6
Robustness due to interbranch feedback
Here we describe the logic behind the working of an interbranch feedback control. Consider the case of a onetier twobranch channel with the feedback on conjugation rate ${\kappa}_{C}$ (Figure 7a of the main text) actuated by the ligand bound state of the signalling receptor in the first tier ${R}^{(1)}$. The steadystate solution to Equations 10–12 corresponding to this channel with ${\kappa}_{C}\to \frac{{\kappa}_{C}}{1+{(\gamma {R}^{(1)})}^{n}}$ is given by.
Therefore, the output can be expressed in terms of $L,\psi $ and $\varphi $ as follows.
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 onetier twobranch 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 $\psi $. Meanwhile, the conjugate species ${Q}^{(1)}$ has a nonmonotonic 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 anticorrelated behaviour of ${R}^{(1)}$ and ${Q}^{(1)}$ due to the feedback results in the output $\theta \equiv {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 $\psi $ and nonsignalling $\varphi $ 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 inputoutput relations in Figure 7d of the main text).
Appendix 7
Heuristic for the influence of feedback action on receptor profiles
The optimised onetier twobranch channel showed a monotonically increasing profile of the signalling receptor $\psi $, and a monotonically decreasing profile of the nonsignalling receptor $\varphi $ (Figure 7b). Here, we provide an understanding for this counterintuitive result using the forms of the channel output $\theta $ and variation in the output ${{\partial}_{L}\theta }_{\psi ,\varphi}$.
Consider the output $\theta $ for a onetier twobranch channel with interbranch feedback on the conjugation rate ${\kappa}_{C}$. Equation 37 gives the explicit form for the output of this channel in terms of input ligand concentration $L$, signalling $\psi $ and nonsignalling $\varphi $ receptor concentrations, feedback strength $n$, feedback sensitivity $\gamma $ and binding rates ${r}_{b},{\kappa}_{b}$. The last three are small parameters in the optimised channel i.e. ${r}_{b}\simeq 0.1,{\kappa}_{b}\simeq 0.05,\gamma \simeq 0.2$. Note that the conjugation and splitting rates are absorbed into $\varphi $, i.e. $\frac{{\kappa}_{C}}{{\kappa}_{S}}\varphi \to \varphi $.
We now look at the output in the limit of high ligand input $L\simeq 2{r}_{b}^{1}\simeq {\kappa}_{b}^{1}$, and low ligand input $L\simeq 0.1{r}_{b}^{1}$. The output forms in these limits are given by Equation 38 and Equation 39.
To analyse the variation of channel output $\theta $ with ligand input $L$, we compute ${\partial}_{L}\theta $ at fixed receptor concentrations $\psi ,\varphi $.
In the limits of high and low ligand input, the form simplifies to
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 $\psi \sim 5{r}_{b}^{1}$ and higher levels of the nonsignalling receptor $\varphi \sim 50{r}_{b}^{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 $\psi \sim 15{r}_{b}^{1}$ and lower levels of the nonsignalling receptor $\varphi \sim 15{r}_{b}^{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 onetier twobranch channel with an interbranch feedback, the optimum signalling receptor $\psi $ has a monotonically increasing spatial profile while the optimum nonsignalling receptor $\varphi $ has a monotonically decreasing spatial profile. This qualitatively explains the result obtained for the optimised onetier twobranch channel in Section ‘Branched architecture with multiple receptors provides accuracy and robustness to extrinsic noise’, Figure 7.
high $L$  low $L$  

output θ  variation ${\mathrm{\partial}}_{L}\theta$  output θ  variation ${\mathrm{\partial}}_{L}\theta$  
low ψ, low Φ  293.95  0.59  23.75  165.93 
high ψ, low Φ  156.53  1.67  22.52  1.12 
low ψ, high Φ  902.07  0.69  67.5  543.46 
high ψ, high Φ  288.43  1.69  40.06  –25.19 
Appendix 8
Results of optimisation of twotier twobranch channels
As discussed in Section ‘Branched architecture with multiple receptors provides accuracy and robustness to extrinsic noise’, additional tiers in a twobranch channel only marginally reduced the inference errors due to extrinsic noise when compared to the optimised onetier twobranch channel (Figure 8). Here, we show that both the receptor profiles and the inputoutput relations of these two optimised twotier twobranch channels are qualitatively similar (Appendix 8—figure 1).
Appendix 9
Additional tiers dampen the effects of feedback to provide stability
Stochastic simulations of the chemical reaction network (CRN) corresponding to the optimised onetier twobranch channel (Figure 9a of the main text) show large fluctuations in the time trajectories of signalling species ${Q}^{(1)}$ about its steadystate 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 nonsignalling receptor $\varphi $. Such amplified fluctuations are absent in the twotier twobranch channel (Figure 9g–i of the main text). Here, we provide heuristics of the differing behaviours in the optimised onetier and twotier channels by analysing a simpler set of CRNs (Appendix 9—figure 1) with the essential elements of these channels.
The dynamics for the first, twospecies CRN is given by.
and that of the second, 3species CRN by
with constant $c$ playing the role of $\psi $, thus providing an upper bound to species 1. In both these sets of equations, we consider a feedback k_{3} such that ${k}_{3}\to \frac{{k}_{0}}{1+{(\gamma {s}_{1})}^{n}}$ where k_{0} represents the reference value of k_{3} in absence of feedback. Appendix 9—figure 2 shows the phase portrait for the dynamics of the 2species CRN when ${k}_{0}\gg {k}_{1}$ with a moderately strong feedback $n=2$. This is representative of ${\kappa}_{C}{S}^{(1)}\gg {r}_{b}L$ in the optimised onetier twobranch channel (Figure 9a of the main text). The nullcline ${\dot{s}}_{1}=0$ (dashed, green curve) acts as a separatrix for the behaviour of this system: if due to fluctuations in species s_{1} the system (${s}_{1},{s}_{2}$) crosses the nullcline from its steadystate (pink point), it sets out on a large trajectory (black line) such that s_{1} remains close to zero while s_{2} 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 nonlinearity of the separatrix is due to the feedback from s_{1} on the rate k_{3} that couples to s_{1} in Equation 43. Higher production rates k_{1} (akin to higher values of ligand $L$ in the optmised channel) bring the steadystate 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 threespecies CRN removes the nonlinearity in ${\dot{s}}_{1}=0$ and provides buffering against this effect (Equation 44).
Appendix 10
Twotier twobranch channel with no feedback
Here we show that the twotier twobranch channel without any interbranch feedback control has fundamentally different optimisation characteristics and a poorer positional inference. Crucially, the optimised profiles of both signalling $\psi $ and nonsignalling $\varphi $ receptors are monotonically decreasing away from the source (Appendix 10—figure 1b, inset).
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 nonsignalling 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 interbranch 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 nonsignalling receptors (Equations 8; 9). Note that this implicitly assumes spatial correlations in the surface concentrations of receptors across the cells.
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 onetier twobranch channel architecture is optimised for input distributions with different decay lengths of mean ligand input $\lambda $ (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 13
Response in inference error due to perturbations in receptor concentrations
The definition of cellular output $\theta $ in a branched architecture involves making a distinction between the signalling and nonsignalling receptors. In conjunction, the direction of the interbranch feedback is from the signalling to the nonsignalling 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 ${\overline{\sigma}}_{X}$ due to the receptor profiles in Appendix 13—figure 1b, c resulting from perturbations in the receptor control parameters A_{2} and B_{2} (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 $\mathrm{\Delta}\psi ,\mathrm{\Delta}\varphi $ from the optimal receptor profile (blue curves in Appendix 13—figure 1b, c), which is computed as follows.
where ${A}_{2}^{\ast}$ and ${B}_{2}^{\ast}$ are the optimum values, and $\delta {A}_{2},\delta {B}_{2}$ 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 ${\overline{\sigma}}_{X}(\mathrm{\Delta}\psi ,\mathrm{\Delta}\varphi )$ around the optimal point ($\mathrm{\Delta}\psi =\mathrm{\Delta}\varphi =0$, red point in Appendix 13—figure 1a) in the $\psi \varphi $ plane has eigenvalues ${\lambda}_{1}\simeq 0.016,{\lambda}_{2}\simeq 6.1$ corresponding to the eigenvectors that are nearly parallel to $\mathrm{\Delta}\varphi $ and $\mathrm{\Delta}\psi $ 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 $\psi $ than to changes in the nonsignalling receptor $\varphi $, implying a stiff direction of control along the former and a sloppy direction of control along the latter.
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 ${\overline{\sigma}}_{X}\le 2\%$. Of the 16 channel parameters, we showed six parameters corresponding to ligand binding rates to the signalling and nonsignalling 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 ${\overline{\sigma}}_{X}\le 2\%$. 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 15
Experimental methods
Fly stocks, endocytic assays, and imaging
Fly stocks used in this study are wildtype w1118, c5GAL4, CAAXGFP (Kyoto DGRC  109823), Port et al., 2014 and UASmyrGarzE740K. 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 (WgAF568) 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 laserscanning confocal microscope using a 60 X/1.42 NA oil objective (with acquisition XY pixel dimensions of $0.138\mu m$ and Z stacks of size $0.5\mu m$).
Analysis
The slicebyslice 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 apicobasal height of cells across the domeshaped 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 databased 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 databased 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.
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/morphogeneticdecoding; copy archived at swh:1:rev:59c184aaa46c5e769a95ea39b48e911d0fc4fd5b).
References

Dealing with noise: the challenge of buffering biological variabilityCurrent Opinion in Systems Biology 1:69–74.https://doi.org/10.1016/j.coisb.2016.12.011

Robust generation and decoding of morphogen gradientsCold Spring Harbor Perspectives in Biology 1:a001990.https://doi.org/10.1101/cshperspect.a001990

Endocytosis and signaling during developmentCold Spring Harbor Perspectives in Biology 6:a017020.https://doi.org/10.1101/cshperspect.a017020

Morphogen transport in epitheliaPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 75:011901.https://doi.org/10.1103/PhysRevE.75.011901

Endocytosis, signaling, and beyondCold Spring Harbor Perspectives in Biology 6:a016865.https://doi.org/10.1101/cshperspect.a016865

Compartmentalization of GPCR signalling controls unique cellular responsesBiochemical Society Transactions 44:562–567.https://doi.org/10.1042/BST20150236

Ras/Mapk signaling from endomembranesMolecular Oncology 3:297–307.https://doi.org/10.1016/j.molonc.2009.06.004

CellCell junctions organize structural and signaling networksCold Spring Harbor Perspectives in Biology 10:a029181.https://doi.org/10.1101/cshperspect.a029181

Stochastic thermodynamics of bipartite systems: transfer entropy inequalities and a maxwell’s demon interpretationJournal of Statistical Mechanics 1:E2016.https://doi.org/10.1088/17425468/2014/02/P02016

`` direct search’’ solution of numerical and statistical problemsJournal of the ACM 8:212–229.https://doi.org/10.1145/321062.321069

Novel regulation of Wnt signaling at the proximal membrane levelTrends in Biochemical Sciences 41:773–783.https://doi.org/10.1016/j.tibs.2016.06.003

Noise, delays, robustness, canalization and all thatCurrent Opinion in Genetics & Development 14:440–445.https://doi.org/10.1016/j.gde.2004.06.001

The measure of success: constraints, objectives, and tradeoffs in morphogenmediated patterningCold Spring Harbor Perspectives in Biology 1:a002022.https://doi.org/10.1101/cshperspect.a002022

BookReceptors: Models for Binding, Trafficking, and SignalingOxford University Press on Demand.

BookInformation Theory, Inference and Learning AlgorithmsCambridge university press.

Promiscuity and specificity in BMP receptor activationFEBS Letters 586:1846–1859.https://doi.org/10.1016/j.febslet.2012.02.043

Ras proteins: paradigms for compartmentalised and isoformspecific signallingCellular and Molecular Life Sciences 64:2575–2589.https://doi.org/10.1007/s0001800771338

Compartmentalized signalling: Ras proteins and signalling nanoclustersThe FEBS Journal 276:1817–1825.https://doi.org/10.1111/j.17424658.2009.06928.x

Morphogen gradients: from generation to interpretationAnnual Review of Cell and Developmental Biology 27:377–407.https://doi.org/10.1146/annurevcellbio092910154148

BookModeling Biomolecular Networks: An Introduction to Systems BiologyOxford University Press.

Morphogens, their identification and regulationDevelopment 131:703–712.https://doi.org/10.1242/dev.01043

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

Positional information and pattern formationCurrent Topics in Developmental Biology 117:597–608.https://doi.org/10.1016/bs.ctdb.2015.11.008

Shaping morphogen gradients by proteoglycansCold Spring Harbor Perspectives in Biology 1:a002493.https://doi.org/10.1101/cshperspect.a002493
Decision letter

Thomas GregorReviewing Editor; Princeton University, United States

Naama BarkaiSenior Editor; Weizmann Institute of Science, Israel

Marcin ZagorskiReviewer; 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 suboptimal 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 tradeoffs, 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 spatiallyfluctuating 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, nonspecific 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 networklike schematic to indicate reaction and flux imbalances could improve presentation.
– In lines 8587, 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 wavelike 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 quiteconsuming 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 fourtiered systems? What sort of computational cost does adding one extra tier there?
https://doi.org/10.7554/eLife.79257.sa1Author 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 suboptimal 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 tradeoffs, 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 nonsignalling 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 onetier twobranch channel for this analysis which yields a minimum inference error of $\overline{\sigma}x\sim 1.9\mathrm{\%}$ i.e. approximately two cell widths. We see that the inference error does not change significantly with most parameters, i.e. it remains within $\overline{\sigma}x\le 2.2\mathrm{\%}$ 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. $\overline{\sigma}x\le 2\mathrm{\%}$. We find that the frequency distribution of optimal feedback parameters whose inference error lies within the range $1.9\mathrm{\%}\le \overline{\sigma}x\le 2\mathrm{\%}$ 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 spatiallyfluctuating 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 steadystate 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 steadystate 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, nonspecific 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 nonsignalling 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 twotier twobranch 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 networklike 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 8587, 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 wavelike 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 quiteconsuming 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 steadystate 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 colocalized 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 colocalization of downstream adaptor proteins to endosomes containing WgDFz2, 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 fourtiered 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 threetier 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.sa2Article and author information
Author details
Funding
Department of Atomic Energy, Government of India (RTI4006)
 Krishnan S Iyer
 Chaitra Prabhakara
 Satyajit Mayor
 Madan Rao
Simons Foundation (287975)
 Madan Rao
DBTWellcome 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 DBTWellcome Trust India alliance (IA/M/15/1/502018).
Senior Editor
 Naama Barkai, Weizmann Institute of Science, Israel
Reviewing Editor
 Thomas Gregor, Princeton University, United States
Reviewer
 Marcin Zagorski, Jagiellonian University, Poland
Version history
 Preprint posted: March 30, 2022 (view preprint)
 Received: April 5, 2022
 Accepted: January 14, 2023
 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
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Developmental Biology
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 centrioleduplication pathways. Collectively, our data suggest a differential requirement for mouse SAS‑6 in centriole formation or integrity depending on PLK4 and centrosome composition.

 Developmental Biology
 Neuroscience
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 reexamined 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. Highresolution 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 circuitspecific recognition. Nerve cellspecific 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.