# Long-range morphogen gradient formation by cell-to-cell signal propagation

Johanna Dickmann, Jochen Rink, Frank Jülicher

► To cite this version:

Johanna Dickmann, Jochen Rink, Frank Jülicher. Long-range morphogen gradient formation by cell-to-cell signal propagation. Physical Biology, 2022, 19 (6), pp.066001. 10.1088/1478-3975/ac86b4. hal-03838337

HAL Id: hal-03838337

https://hal.inrae.fr/hal-03838337

Submitted on 3 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Long-range morphogen gradient formation by cell-to-cell signal propagation

To cite this article: Johanna E M Dickmann et al 2022 *Phys. Biol.* **19** 066001

View the [article online](#) for updates and enhancements.

You may also like

- Impulsive signaling model of cytoneme-based morphogen gradient formation
  Hyunjoong Kim and Paul C Bressloff

- Partial differential equations for self-organization in cellular and developmental biology
  R E Baker, E A Gaffney and P K Maini

- Roadmap for the multiscale coupling of biochemical and mechanical signals during development
  Pierre-François Lenne, Edwin Munro, Idse Heemskerk et al.

<figure><img src="image_2.png" /></figure>
OPEN ACCESS

<figure><img src="image_3.png" /></figure>

CrossMark

RECEIVED

13 April 2022

REVISED

12 July 2022

ACCEPTED FOR PUBLICATION

3 August 2022

PUBLISHED

7 September 2022

Original content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.
Any further distribution
of this work must
maintain attribution to
the author(s) and the
title of the work, journal
citation and DOI.

<figure><img src="image_4.png" /></figure>

PAPER

# Long-range morphogen gradient formation by cell-to-cell signal propagation

Johanna E M Dickmann¹,⁴ ID, Jochen C Rink²,⁵ ID and Frank Jülicher¹,³,* ID

¹ Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Straße 38, 01187 Dresden, Germany

² Max Planck Institute for Multidisciplinary Sciences, Am Faßberg 11, 37077 Göttingen, Germany

³ Cluster of Excellence, Physics of Life, TU Dresden, 01307 Dresden, Germany

* Author to whom any correspondence should be addressed.

⁴ Present address: Laboratoire Reproduction et Développement des Plantes, ENS de Lyon, INRAE, CNRS, UCB Lyon 1, F-69364, Lyon, France.

⁵ Former affiliation: Max Planck Institute of Molecular Cell Biology and Genetics, Pfotenhauerstraße 108, 01307 Dresden, Germany.

E-mail: johanna.dickmann@ens-lyon.fr, jochen.rink@mpinat.mpg.de and julicher@pks.mpg.de

**Keywords:** morphogen gradients, tissue patterning, long-range patterning, cellular signaling, Wnt signaling, development, planarians

## Abstract

Morphogen gradients are a central concept in developmental biology. Their formation often involves the secretion of morphogens from a local source, that spread by diffusion in the cell field, where molecules eventually get degraded. This implies limits to both the time and length scales over which morphogen gradients can form which are set by diffusion coefficients and degradation rates. Towards the goal of identifying plausible mechanisms capable of extending the gradient range, we here use theory to explore properties of a cell-to-cell signaling relay. Inspired by the millimeter-scale *wnt*-expression and signaling gradients in flatworms, we consider morphogen-mediated morphogen production in the cell field. We show that such a relay can generate stable morphogen and signaling gradients that are oriented by a local, morphogen-independent source of morphogen at a boundary. This gradient formation can be related to an effective diffusion and an effective degradation that result from morphogen production due to signaling relay. If the secretion of morphogen produced in response to the relay is polarized, it further gives rise to an effective drift. We find that signaling relay can generate long-range gradients in relevant times without relying on extreme choices of diffusion coefficients or degradation rates, thus exceeding the limits set by physiological diffusion coefficients and degradation rates. A signaling relay is hence an attractive principle to conceptualize long-range gradient formation by slowly diffusing morphogens that are relevant for patterning in adult contexts such as regeneration and tissue turn-over.

## 1. Introduction

Morphogens and morphogen gradients are key concepts in developmental biology. Morphogens are defined as secreted signaling molecules that spread through the tissue and specify cell fate in a concentration-dependent manner [1]. The decrease in morphogen concentrations as a function of distance from the site of secretion (the source) is referred to as morphogen gradient. The existence of morphogens has been predicted [2–4] and studied theoretically [3, 5, 6], decades before they were first observed experimentally [7, 8]. The current concept combines the ideas of ‘form producers’ introduced by Turing [3], and the idea of positional information

introduced by Wolpert [4] with a large body of experimental observations [7–18]. The emergence of the morphogen concept highlights the importance of theory in the elucidation and understanding of morphogenesis.

Morphogen gradients can form when morphogens are secreted in a local source and spread by diffusion through a target tissue where they are also degraded [5, 12, 14]. This model can account for observations of morphogen gradients in different model organisms including the fruit fly embryo [12–14, 19] and the zebrafish embryo [15, 16]. In these systems, morphogen gradients typically reach length scales in the range of tens to a few hundreds of micrometers [7, 12–16, 20, 21]. Salient biophysical
parameters, i.e. diffusion coefficient and morphogen degradation rate, have been estimated in quantitative experiments [12–16, 21]. These measured parameter values can explain the formation of morphogen gradients spanning the observed tissue dimensions within the relevant developmental time interval. Tissue patterning by morphogen gradients does not only occur during development on the micrometer dimensions of early embryos, but can also occur at larger spatial dimensions of adult tissues. Examples for patterning of adult tissues by soluble signaling molecules include the organization of mammalian liver lobules [22, 23], axolotl regeneration [24], as well as regeneration and steady-state turn-over of planarian flatworms [25–31]. These examples show that signaling gradients are relevant to pattern tissues on millimeter length scales and possibly even up to centimeter length scales and thus require mechanisms capable of gradient formation over post-embryonic length scales.

In principle, morphogen gradients formed by diffusion and degradation can achieve arbitrary decay lengths, either by increasing the diffusion coefficient or by increasing the morphogen life time. However, constraints on parameter values (e.g., the diffusion coefficient of proteins in aqueous solution [21]) or the typical time scales of biological processes (e.g., observed growth or regeneration times) impose limits on the system's dimensions in which gradient formation by diffusion and degradation is feasible. For instance, it has been suggested that extracellular diffusion of Fgf8 ($D = 53 \pm 8 \ \mu m^2 s^{-1}$ [15]) together with the estimated degradation rate ($k = D/\lambda^2 = (1.3 \pm 0.3) \times 10^{-3} \ s^{-1}$ [15]) gives rise to an Fgf8 gradient with a decay length of $\lambda = 197 \pm 7 \ \mu m$ [15] during embryonic development of the zebrafish on a time scale of around 13 min, thus matching the dimension of the embryonic structure to be patterned and the time scale relevant for development. However, many morphogens are lipid-modified [32–34] and therefore poorly soluble in the aqueous tissue environment, which results in slow diffusion [35–37]. This is in particular the case for Wnt [34] and it has been debated how far Wnt can spread by diffusion [35, 36, 38–42]. Already when assuming a diffusion coefficient one order of magnitude smaller than the above cited value for Fgf8, molecular life times of the order of months would be required for the generation of the millimeter-scale Wnt signaling gradient that organizes the main body axis in planarians [25–27, 30, 31]. It has recently been suggested that morphogens could be lost by moving out of the tissue, challenging long-range gradient formation [20]. This has been shown in an engineered system to be reduced by receptor binding [20]. Thus, both the slow diffusivity of morphogens and their potential loss out of the tissue during diffusive spreading pose the fundamental question of what other means of

gradient formation could be suited to generate such long-rage gradients.

Directed molecular transport is one concept to increase the decay length of gradients that has been examined [43]. A directed motion of molecules through the system (drift) can result in faster and more long-ranged gradient formation as compared to diffusive spreading, see appendix A. Drift could be caused by active intracellular transport processes, or by extracellular fluid flows, generated for example by coordinated cilia beating or by coordinated contractions [43]. An alternative mechanism for signal spreading is represented by cytonemes [44]. These thin, actin-dependent membranous structures extend from the plasma membrane of the signaling cell and can contact a signal receiving-cell more than 40 µm away [41, 44]. Cytonemes have therefore been suggested to contribute to signaling gradient formation on the micrometer range [41]. An alternative mechanism for gradient formation has been introduced in the ‘bucket brigade’ model by Kerszberg and Wolpert in 1998 [45]. In this model, morphogens move between transmembrane receptors along and between cells. Receptors recently bound by signaling molecules become refractory, preventing ‘backward’ spreading of the morphogen. Further, the morphogen can be handed over from receptors on one cell to those on another [45]. This leads to the formation of gradients of receptor-bound morphogens. Kerszberg and Wolpert briefly discuss whether the hand-over of morphogens between cells could be replaced by the idea of positive feedback in which new morphogens are generated in response to a morphogen binding to a receptor [45]. They conclude that this would prevent the morphogen concentration from decaying away from the local source and would thus not be a feasible concept of gradient formation [45].

In this paper, we revisit the idea of positive feedback as a means of gradient formation. Inspired by the evidence for Wnt-mediated *wnt* expression and the ensuing *wnt* expression gradients in planarians [27, 29, 31, 46] we explore how tissue-scale morphogen gradients can form by signaling relays. We develop a model in which the extracellular morphogen gradient and resulting intracellular gradient of signaling activity are connected by a positive feedback loop (i.e. Wnt-mediated *wnt* expression). This describes morphogen-mediated morphogen production and introduces a signaling relay by which morphogens can propagate in the tissue. In fact, the whole tissue effectively becomes a signaling-dependent source due to the positive feedback. A signaling-independent source of morphogen at the boundary initiates the signaling relay (inspired by the putative role of a tail-tip expressed Wnt ligand in planarians [31, 46]). In contrast to the work by Kerszberg and Wolpert [45], our results show that a
signaling relay can induce stable morphogen and signaling gradients. Moreover, we show that spreading of the signal by a cell-to-cell relay leads to effective diffusion of the morphogen through the system. That is, individual molecules do not have to travel long distances but are instead produced anew further down the system in a feedback-dependent manner. Furthermore, the signaling relay leads to an effective degradation rate smaller than the molecular degradation rate. This way, a cell-to-cell signaling relay can efficiently generate long-range morphogen and signaling gradients for physiologically relevant molecular diffusion coefficients and degradation rates on relevant time scales. Moreover, when this feedback-dependent production is combined with secretion polarity in the plane of the cell field, the cell-to-cell relay additionally leads to the emergence of an effective drift resembling directed transport of morphogens through the system in the absence of molecular transport. This can further increase the decay length of the signaling and morphogen gradients. Finally, the range of the gradients can be set by regulating the feedback strength, allowing the system to adjust gradient range without the need to change diffusivity or molecular life time of the morphogen. Altogether, this makes cell-to-cell signaling relays an attractive concept to explore in systems with millimeter-scale patterns as observed during regeneration such as in planarians, hydra, zebrafish, and salamanders.

## 2. Gradient formation by cell-to-cell signaling relay

In this paper, we focus on morphogen gradient formation by a cell-to-cell relay where individual cells transmit the signaling information to their neighbors. In particular, our approach takes the coupling between intracellular levels of signaling activity and the extracellular morphogen levels into account. This differs from scenarios where intracellular signaling gradients are strictly downstream of extracellular morphogen levels. To capture the coupling, we account for individual cells that react to the morphogen concentration in the extracellular space by increasing their intracellular signaling activity, captured as the concentration of an intracellular signaling molecule downstream of the morphogen. This increased signaling activity eventually leads to the production and secretion of more morphogen, increasing the extracellular morphogen concentration. This way, the information is relayed from cell to cell as the whole cell field becomes a signaling-dependent source. The level of intracellular signaling activity and concentration of extracellular morphogen are thus interdependent. The gradient is positioned and oriented by a localized, signaling-independent source of morphogen at the proximal boundary. From this signaling-independent source,

the morphogen profile propagates via the cell-to-cell signaling relay towards the distal end. Note that the feedback is present along the whole length of the system. That is, the cells at the proximal boundary of the system serve a dual function: they produce morphogen firstly in a signaling-independent manner by virtue of being the signaling independent source and secondly in a signaling-dependent manner depending on their intracellular signaling level like all the other cells in the cell field.

The schematic of the cell-based model is shown in figure 1. We denote by $b_n$ the signaling activity in cell $n$, where $n = 0, \dots, N-1$, and $N$ is the number of cells along the length of the system. Thus, $b_n$ denotes the concentration of an intracellular signaling molecule that is downstream of the morphogen. Similarly, we introduce the morphogen concentration in the extracellular space $a_n$, where $n$ refers to the extracellular space between the cells $n-1$ and $n$. For simplicity we use a discrete description of the extracellular space with concentrations $a_n$ averaged within a cell length. The dynamic equations for these concentrations read

$$ \partial_t a_n = D \frac{a_{n-1} - 2a_n + a_{n+1}}{\delta^2} - k_A a_n + s_A \theta_{w-n} + f(b_{n-1}, b_n), \quad (1) $$

$$ \partial_t b_n = s_B - k_B g(a_n, a_{n+1}) b_n. \quad (2) $$

Here, $D$ denotes the diffusion coefficient, $\delta$ the width of a cell, and $k_A$ the rate of morphogen loss which could be due to degradation, internalization or leakage. For simplicity, we refer to $k_A$ as the degradation rate of the morphogen. The constant rate of production in the signaling-independent source region (shaded in red in figure 1), is denoted $s_A$. The cells that are part of this signaling-independent source region at the proximal end ($n=0$) of the system are specified by the function $\theta_{w-n}$, where $\theta_{w-n} = 1$ inside the source for $n < w$, and $\theta_{w-n} = 0$ outside the source for $n > w$. At the source boundary, $n = w$, we use $\theta_{w-n} = p$ where $p$ describes a cellular asymmetry of molecular secretion, see below. The distal end is at $n = N$. In equation (1), the positive feedback gives rise to an additional source $f(b_{n-1}, b_n)$ in the whole cell field that depends on the signaling activities in the two adjacent cells, $b_{n-1}$ and $b_n$. We refer to this as the signaling-dependent source. Equation (2) describes the dynamics of the signaling activity. Here, $s_B$ denotes the rate of production of the intracellular signaling molecule, $k_B$ is a degradation rate and $g(a_n, a_{n+1})$ describes the positive feedback regulating the degradation of the intracellular signaling molecule by the morphogen concentrations in the adjacent extracellular spaces $a_n$ and $a_{n+1}$. The boundary conditions are discussed in appendix B. The signaling-dependent source $f$, as well as the regulation of the degradation of the intracellular signaling molecule $g$ are given by
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Schematic representation of the model geometry in one dimension. The extracellular morphogen $a_n$ (depicted as red dots) is constantly produced in the signaling-independent source region shaded in red at the proximal end ($n = 0$). Moreover, it is produced throughout the system, which becomes a signaling-dependent source due to the signaling relay. It spreads across the extracellular spaces, where it is also degraded, towards the distal end ($n = N$). The intracellular signaling activity $b_n$ (depicted as blue dots) is confined to the cell interior and the intracellular signaling molecules do not move between cells of size $\delta$. The dynamics of the system is given by equations (1) and (2).</figcaption></figure>

$$f(b_{n-1}, b_n) = \alpha p \frac{b_{n-1}}{c_B + b_{n-1}} + \alpha(1-p) \frac{b_n}{c_B + b_n}, \quad (3)$$

$$g(a_n, a_{n+1}) = \frac{c_A}{c_A + a_n + a_{n+1}}, \quad (4)$$

where we use the Hill activator function to model the positive feedback of the intracellular signaling molecule $b_n$ on morphogen production and $c_B$ is a Hill activation threshold, i.e. the concentration of $b_n$ at which half-maximal activation is reached. Similarly, $g$ is the Hill inhibitor function and $c_A$ is the Hill inhibitor coefficient at which half-maximal inhibition is reached. Here, we use Hill exponents of 1 because we consider cells that respond linearly to weak stimuli before the response saturates at high stimuli. The maximal signaling-induced secretion rate is denoted $\alpha$. Equation (3) describes polarized morphogen secretion. Such polarized secretion might for instance arise from distinct activities of the two lateral sides of the cell to secrete morphogen. This secretion polarity $p$ is associated with cell polarity in the tissue akin to planar cell polarity. In the case of non-polar cells $p = \frac{1}{2}$. Extreme secretion asymmetries correspond to $p = 0$ and $p = 1$, respectively.

Equations (1)–(4) describe a positive feedback between morphogen concentrations and concentrations of the intracellular signaling molecule, where the parameter $\alpha$ sets the overall feedback strength. This feedback consists of two elements: Equation (3) contains a positive feedback element stimulating morphogen secretion with increasing intracellular signaling activity, the signaling-dependent source. Equation (4) also is a positive feedback element, increasing signaling activity for increased morphogen levels. This feedback element acts on the degradation of the intracellular signaling molecule and thus is positive due to inhibition of inhibition, a principle commonly observed in signaling pathways.

Our model is inspired by Wnt signaling. It is known that high extracellular Wnt concentrations prevent intracellular $\beta$-catenin degradation [47] as captured by equation (2), corresponding to inhibition

of inhibition. Equation (1) is motivated by observations in planarian flatworms suggesting that high intracellular $\beta$-catenin levels lead to increased expression of several Wnt-signaling pathway components, including several Wnt ligands [31].

The cell-to-cell relay involving positive feedback is conceptually different from gradient formation by diffusion and degradation. The relay turns the whole system into a signaling-dependent source. Thus, individual signaling molecules do not have to move a distance spanning the entire gradient. They can form a long-range gradient even if they only move between neighboring cells. In the following, we analyze the

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 2.** Steady-state concentration profiles of the morphogen (a) and the intracellular signaling activity (b). The signaling-independent source region is shaded in red. Circles denote the numerically obtained steady-state solution to the non-linear equations, lines indicate the analytical approximation to the steady state. The insets show the same profiles on a log-y scale. Note that we subtracted the plateau value reached at the distal end of the system in order to visualize the exponential decay towards this value. Parameters: table H1, set A in the appendix.</figcaption></figure>
impact of the positive feedback on the shape of the gradient and the dynamics of gradient formation.

### 3. Regulation of the gradient range by signaling feedback

#### 3.1. Morphogen and signaling gradients at steady state

We first consider steady-state concentration profiles in order to discuss the range of the morphogen and signaling gradients. We quantify the range of the gradients as the decay length $\lambda$ over which they decay. Note that here we measure all lengths relative to the cell size $\delta$ and $\lambda$ is therefore the number of cells corresponding to the gradient range. We can then analyze the impact of the positive feedback on the gradient range. The steady-state solutions to equations (1) and (2) are denoted $a_n^*$ and $b_n^*$ and obey $\partial_t a_n^* = 0$,

$\partial_t b_n^* = 0$. We can obtain these solutions numerically, see figure 2 and appendix C. However, to discuss the impact of the feedback, we use an approximation which we can obtain analytically. Note that at steady-state, equation (2) implies that the profile of the intracellular signaling activity is linearly dependent on the morphogen concentration profile:

$$b_n^* = \frac{s_B}{k_B c_A} (c_A + a_n^* + a_{n+1}^*). \quad (5)$$

Thus, at steady state, the signaling activity-gradient can be obtained from the morphogen gradient using equation (5).

We focus our discussion on the morphogen gradient $a_n^*$ in the following. We obtain an approximation for the steady state by matching simplified solutions inside the source region with $n = 0, \dots, w$, and outside of the source region with $n = w + 1, \dots, n$.

$$a_n^* = \begin{cases} a_c^{\text{in}} + C_1^{\text{in}} \exp(n/\lambda_1^{\text{in}}) + C_2^{\text{in}} \exp(-n/\lambda_2^{\text{in}}) & \text{for } n = 0, \dots, w, \\ a_c^{\text{out}} + C_1^{\text{out}} \exp(n/\lambda_1^{\text{out}}) + C_2^{\text{out}} \exp(-n/\lambda_2^{\text{out}}) & \text{for } n = (w+1), \dots, N. \end{cases} \quad (6)$$

Here, the gradient amplitudes $C_1^{\text{in}}$, $C_2^{\text{in}}$, $C_1^{\text{out}}$, and $C_2^{\text{out}}$ are determined by the boundary conditions, as well as the matching condition at the source boundary, see appendix D. The decay lengths $\lambda_{1,2}^{\text{in, out}}$ can be expressed explicitly, see appendix D. We focus on the decay of the gradient outside the source, which occurs over a range $\lambda_2^{\text{out}}$. This term dominates for large system sizes outside of the source (see insets in figure 2) and therefore sets the gradient range

$$\lambda := \lambda_2^{\text{out}}. \quad (7)$$

The decay length $\lambda$ can be related to effective transport coefficients which we discuss in the next section. Based on these effective transport coefficients, we then discuss the regulation of the gradient range by the signaling-feedback strength in section 3.3. We show in particular, that the decay length is increased compared to gradient formation by diffusion and degradation, see figure 3.

#### 3.2. Effective coefficients of molecular transport

The model for gradient formation by cell-to-cell signaling relay introduced in the last sections describes the kinetics of molecular transport and signaling at the scale of cells. In this model, a gradient with a long range emerges revealing that the system generates concentration patterns at larger scales. In the spirit of a hydrodynamic theory, molecular concentration profiles can be captured by an effective continuum description at scales much larger than

cells. In this limit, it has the form of a convection–diffusion–degradation equation:

$$\partial_t C = D_{\text{eff}} \partial_x^2 C - \partial_x(v_{\text{eff}} C) - k_{\text{eff}} C + s\theta(S-x). \quad (8)$$

Here, $D_{\text{eff}}$ is an effective diffusion coefficient, $k_{\text{eff}}$ is an effective degradation rate, and the velocity $v_{\text{eff}}$ describes an effective drift. Production with rate s is limited to a local source of width S, which corresponds to the signaling-independent source in the cell-based model. The spatial coordinate in the direction along which the gradient forms is denoted x, and $\theta$ denotes the Heaviside function, where $\theta(y) = 0$ for $y < 0$ and $\theta(y) = 1$ for $y \ge 0$. In equation (8), $0 < x \le L$ with a tissue of size L and the source at the proximal boundary. Positive values of v indicate drift in positive x-direction. We can estimate effective transport coefficients as well as the effective degradation rate by comparing the large scale dynamics of our cell-based model to the corresponding behavior of the continuum model given by equation (8), see appendix E. We obtain

$$D_{\text{eff}} \approx \phi \left( \frac{D}{\delta^2} + \frac{1}{2} \alpha z^{\text{out}} \right), \quad (9)$$

$$k_{\text{eff}} \approx \phi(k_{\text{A}} - 2\alpha z^{\text{out}}), \quad (10)$$

$$v_{\text{eff}} \approx \phi(2p - 1)\alpha z^{\text{out}}, \quad (11)$$

where
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 3.** The decay length of the morphogen gradients reaches a maximum close to the critical feedback strength $\alpha^{\text{crit}}$ in the relay model. (a) Decay length $\lambda$ in cells as a function of the relative feedback strength $\alpha/\alpha^{\text{crit}}$. (b) Decay length $\lambda$ in cells as a function of the relative distance from the critical feedback strength $\delta\alpha/\alpha^{\text{crit}}$, where $\delta\alpha = \alpha^{\text{crit}} - \alpha$. The log–log plot exposes the scaling relationship stated in equation (15). Solid blue lines: decay length according to the analytical approximation to the steady-state solution (equation (13)), blue dots: decay length based on the numerical steady-state solution (appendix F), dashed blue line: decay length in the limit of a linear feedback response, i.e. in the absence of feedback saturation, purple line: decay length in the limit of $\alpha = 0$, i.e. a diffusion–degradation model, vertical black line: $\alpha = \alpha^{\text{crit}}$. Parameters: table H1, set B in the appendix.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 4.** Polarized secretion towards the distal end ($v_{\text{eff}} > 0$) increases the decay length of the morphogen gradients at steady state. Decay length $\lambda$ in cells as a function of the relative distance from the critical feedback strength $\delta\alpha/\alpha^{\text{crit}}$, where $\delta\alpha = \alpha^{\text{crit}} - \alpha$. Solid lines: decay length according to the analytical approximation to the steady-state solution (equation (13)), dots: decay length based on the numerical steady-state solution (appendix F), dashed lines: decay length in the limit of a linear feedback response, i.e. in the absence of feedback saturation. Blue: absence of secretion polarity, $p = 0.5$, green: preferred secretion towards the distal end, $p = 0.99$, red: preferred secretion towards the proximal end, $p = 0.01$. Parameters: table H1, set C in the appendix.</figcaption></figure>

$$z^{\text{out}} = \frac{c_A c_B k_B}{s_B} \left( c_A + \frac{c_A c_B k_B}{s_B} + 2a_c^{\text{out}} \right)^{-2} . \quad (12)$$

Here $a_c^{\text{out}}$ is given in appendix D, see equation (D.3), and the coefficient $\phi$ is given in appendix E, see equation (E.24). Note that the effective diffusion coefficient and degradation rate contain the molecular diffusion coefficient $D/(\delta^2)$ and the molecular degradation rate $k_A$, respectively. Both are modulated by the additive contribution from the relay $\alpha z^{\text{out}}$, where $\alpha$ is the feedback strength and $z^{\text{out}}$, with units of inverse concentration, characterizes the linear response of the morphogen system at the homogeneous reference state, see appendix D. They are further modulated by the multiplicative coefficient $\phi$. In contrast, the effective drift is purely due to polarized secretion of molecules with polarity $p$ in the complete absence of a molecular drift. This implies that molecules appear to be transported through the system, when actually they are produced and secreted with secretion polarity $p$. Accordingly, the effective drift vanishes in the absence of secretion polarity at $p = \frac{1}{2}$. The polarized secretion $2p - 1$ is again modulated by $\alpha z^{\text{out}}$ and $\phi$.

Using these effective quantities, the decay length of the morphogen gradient at steady-state is given as:

$$\lambda \approx \frac{2D_{\text{eff}}}{(v_{\text{eff}}^2 + 4D_{\text{eff}}k_{\text{eff}})^{1/2} - v_{\text{eff}}} . \quad (13)$$

Based on this expression for the decay length, we can discuss the impact of the signaling relay on the gradient range, based on the effective quantities we identified above. We begin this discussion by focusing on the influence of the feedback strength $\alpha$.

### 3.3. Regulation of gradient range by signaling-feedback strength

The positive feedback strength $\alpha$ can regulate the gradient range by modulating all effective quantities introduced above. In the absence of positive feedback, $\alpha = 0$, the model becomes a diffusion–degradation system and thus the decay length simplifies to the well-known expression $\lambda = \delta^{-1}(D/k)^{1/2}$ [14], see appendix A, equation (A.5). In the presence of positive feedback, $\alpha > 0$, the range of the gradient is increased as compared to the case $\alpha = 0$, see figure 3. In particular, the decay length exhibits a maximum near a feedback strength $\alpha^{\text{crit}}$. Here, $\alpha^{\text{crit}}$ corresponds to a critical point that occurs for a linear feedback response when $b_n/(c_B + b_n)$ is replaced by $b_n/c_B$, which is a linear function of $b_n$. Such a linear regime occurs in the limit when the Hill activation threshold $c_B$ is much larger than the signaling activity $b_n$. In this case, the concentration grows without bounds for $\alpha > \alpha^{\text{crit}}$ with

$$\alpha^{\text{crit}} = \frac{c_A c_B k_A k_B}{2s_B}, \quad (14)$$

and no steady state exists. In this linear response regime, $\lambda$ diverges at the critical feedback strength. For $p = \frac{1}{2}$, $\lambda$ increases near the critical point as

$$\lambda \propto (\delta\alpha)^{-1/2}, \quad (15)$$
where $\delta\alpha = \alpha - \alpha^{\text{crit}}$, see dashed line in figure 3(b). Because the Hill function $b_n/(c_B + b_n)$ is non-linear, $\lambda$ remains finite. The maximal value reached near $\alpha^{\text{crit}}$ scales with $c_B$ as

$$ \lambda \approx (c_B)^{1/4} \left[ \frac{4D/\delta^2 + k_A}{8k_A(s_B/k_B)^{1/2}} \right]^{1/2} . \quad (16) $$

In the presence of secretion polarity oriented towards the distal end, $p > \frac{1}{2}$, $\lambda$ diverges near the critical point as

$$ \lambda \propto (\delta\alpha)^{-1}. \quad (17) $$

It is capped by the non-linearities at the maximal value given by

$$ \lambda \approx (c_B)^{1/2} \frac{2p - 1}{4(s_B/k_B)^{1/2}}. \quad (18) $$

The maximal $\lambda$ corresponds to a minimal value of $k_{\text{eff}}$ which never reaches 0. Correspondingly, as a function of $\alpha$, $\lambda$ passes through a maximum near $\alpha^{\text{crit}}$, see figure 3(a), equations (16) and (18). Note that the increase of gradient range due to the signaling relay is particularly strong close to $\alpha^{\text{crit}}$, see figure 3(b).

Taken together, the positive feedback in the relay increases the decay length of the steady-state signaling gradient compared to diffusion and degradation, see figure 3. Moreover, the range of the signaling gradient can be regulated solely by changes in the feedback strength $\alpha$. This is different from gradient formation by diffusion and degradation alone, where the decay length of the gradient is regulated by changes of the degradation rate or the diffusion coefficient. We next focus on the impact of secretion polarity on gradient range.

### 3.4. The influence of the emergent morphogen drift on gradient range

In the signaling relay, an emergent drift with effective velocity $v_{\text{eff}}$ arises as a consequence of polarized secretion of morphogen, with $p \neq \frac{1}{2}$, see equation (11). Thus, the emergent morphogen drift $v_{\text{eff}}$ occurs in the absence of an actual directed transport of molecules through the system. In particular, the emergent drift has positive values for secretion polarity oriented towards the distal end i.e. for $p > \frac{1}{2}$. We observe that in this case the emergent drift increases the decay length of the steady-state signaling gradient reached at the critical feedback strength $\alpha^{\text{crit}}$, see figure 4, compare equations (16) and (18). This dependence of decay length on effective drift is also captured by the continuum description, see equation (8).

The emergent morphogen drift also influences the sensitivity of the system to changes in the feedback strength as the feedback approaches its critical strength, compare equations (15) and (17). Thus, for preferred secretion towards the distal end, the decay length increases more rapidly as the feedback strength approaches its critical value as compared to the case without effective drift $p = \frac{1}{2}$, see figure 4.

<figure><img src="image_11.png" /><figcaption>Figure 5. The relaxation time of the signaling-relay mechanism reaches a finite value at the critical feedback strength. Relaxation time $\tau_{\text{max}}$ as a function of the relative distance from the critical feedback strength $\delta\alpha/\alpha^{\text{crit}}$, where $\delta\alpha = \alpha^{\text{crit}} - \alpha$. Dots: relaxation time obtained by linearizing the dynamics around the numerically obtained steady-state solution, solid lines: relaxation time obtained by linearizing the dynamics around the analytical approximation to the steady-state solution, dashed lines: the relaxation time in the limit of a linear feedback response, i.e. in the absence of feedback saturation. Blue: absence of secretion polarity, $p = 0.5$, green: preferred secretion towards the distal end, $p = 0.99$. $\delta\alpha = \alpha - \alpha^{\text{crit}}$. Parameters: table H1, set C in the appendix.</figcaption></figure>

Taken together, polarized secretion can increase the decay length of the steady-state gradients and can make the decay length more responsive to changes in the feedback strength. We next investigate how the feedback strength influences the dynamics of the system, in particular the time needed to build the gradient.

## 4. Dynamics of gradient formation with cell-to-cell signaling relay

In order to exert a patterning function, a morphogen gradient has to form. Therefore, apart from the decay length of the steady-state profile, the time it takes to form a profile is an important property of gradient formation. We use the slowest relaxation time of the system to reach steady state as a measure of the time it takes to form a gradient. We define the slowest relaxation time as the slowest exponential relaxation of concentrations close to steady state

$$ a_n(t) = a_n^* + \delta a_n(t), \quad (19) $$

$$ b_n(t) = b_n^* + \delta b_n(t), \quad (20) $$

where $a_n(t)$ and $b_n(t)$ denote the concentration profiles of the morphogen and the intracellular signaling molecule at time $t$, $a_n^*$ and $b_n^*$ denote the respective steady-state concentration profiles and $\delta a_n(t)$ and $\delta b_n(t)$ the respective time-dependent deviations from it. Using the numerically determined steady-state profiles, we can then obtain the relaxation modes numerically, for details see appendix G, see figure 5 circles. We observe that the relaxation time increases as $\alpha$ increases and that the dependence on $\alpha$ is steepest close to $\alpha^{\text{crit}}$, see figure 6. Compared to the case $p = \frac{1}{2}$, the relaxation time close to $\alpha^{\text{crit}}$ becomes shorter when secretion polarity is oriented towards
<figure><img src="image_12.png" /></figure>

**Figure 6.** The relaxation time of the relay mechanism reaches a local maximum close to the critical feedback strength. (a) The relaxation time as a function of the relative feedback strength. (b) Blow-up of the boxed region in (a). Dots: relaxation time obtained by linearizing the dynamics around the numerically obtained steady-state solution, solid lines: relaxation time obtained by linearizing the dynamics around the analytical approximation to the steady-state solution, dotted lines: relaxation time in the limit of a linear feedback response, i.e. in the absence of feedback saturation, vertical black lines: critical feedback strength. Blue: absence of secretion polarity, *p* = 0.5, green: preferred secretion towards the distal end, *p* = 0.99. Parameters: table H1, set C in the appendix.

<figure><img src="image_13.png" /></figure>

**Figure 7.** The trade-off between the decay length and the relaxation time of the relay model for changing feedback strength. Dots: relaxation time based on dynamics linearized around the numerical steady-state solution and the decay length based on the numerical steady-state solution. Solid lines: relaxation time based on the dynamics linearized around the analytical approximation to the steady-state solution and the decay length according to the analytical approximation to the steady-state solution (equation (13)). Dashed lines: relaxation time and the decay length in the limit of a linear feedback response, i.e. in the absence of feedback saturation. Blue: absence of secretion polarity, *p* = 0.5, green: preferred secretion towards the distal end, *p* = 0.99. Black lines above and below the curves: scaling behavior of τmax ∝ λ² and τmax ∝ λ, respectively. τmax ∝ λ² is the scaling behavior of gradient formation by diffusion and degradation alone. The region shaded in red corresponds to gradient decay lengths and time scales relevant for planarian regeneration, i.e. 0.125 mm to 0.5 cm and 1 h to 14 days, respectively. Parameters: table H1, set C in the appendix.

the distal end, $p > \frac{1}{2}$, see figures 5 and 6. We compare these results to the relaxation modes obtained based on the approximation of the steady-state profile discussed in section 3.1, equation (6). This approximation works well away from the critical feedback strength but underestimates the relaxation time near $\alpha^{\text{crit}}$, figure 5 solid lines. Using a linear feedback, $b_n/c_B$, provides a good approximation for the relaxation time away from the critical feedback strength, but leads to a divergence of the relaxation time at $\alpha^{\text{crit}}$, see figures 5 and 6 dashed lines.

Taken together, the relaxation time slows down with increasing feedback strength $\alpha$ and reaches a local maximum near $\alpha^{\text{crit}}$ where it stays finite due to the non-linearities in the feedback, see figure 5. However, different from the decay length, the relaxation time is very asymmetrical around $\alpha^{\text{crit}}$, reaching higher levels for $\alpha > \alpha^{\text{crit}}$ as compared to $\alpha < \alpha^{\text{crit}}$. Secretion polarity towards the distal end can speed up the relaxation times close to $\alpha^{\text{crit}}$, see figures 5 and 6.

## 5. Scaling of decay length and relaxation time in gradient formation

Patterning of tissues requires a range of gradient decay lengths from tens of micrometers up to the order of millimeters and potentially even centimeters. In order to exert a patterning function, those gradients

have to be formed on time scales compatible with the respective patterning process, for instance days to weeks in case of planarian regeneration. In this section, we discuss the scaling of relaxation time with decay length as the decay length is increased by changing the feedback strength.

As discussed above, the decay length of the gradients at steady state is largest for feedback strengths close to the critical one at which it reaches a maximum, see figure 3(a). In contrast, the relaxation time reaches a local maximum close to $\alpha^{\text{crit}}$ and stays high at values $\alpha > \alpha^{\text{crit}}$, see figure 6. Therefore, generating gradients is fastest at feedback strengths below and up to the critical feedback strength, see figure 7. Figure 7 shows the scaling of relaxation time $\tau_{\text{max}}$ with decay length $\lambda$. We observe that for $p = \frac{1}{2}$, i.e. vanishing effective drift $v_{\text{eff}}$, the relaxation time scales as $\tau_{\text{max}} \propto \lambda^\nu$ for $\alpha < \alpha^{\text{crit}}$, with $\nu \approx 2$, see figure 7. In other words, doubling the decay length increases the relaxation time by a factor of four. This is the same scaling relationship observed for gradient formation by diffusion and degradation, see appendix A, equation (A.8). Note that for gradient formation by diffusion and degradation, the decay length and the relaxation time are regulated by changes in the degradation rate, whereas in the signaling relay discussed here, they are regulated by changes in the feedback strength $\alpha$.

For secretion polarity oriented distally, $p > \frac{1}{2}$, the scaling relationship is altered, with $\nu < 2$, see figure 7. As $\lambda$ reaches a maximum, deviations from these
scaling relations occur. These are most pronounced when estimated using our approximations, see solid lines in figure 7. Taken together, the signaling relay combined with secretion polarity outperforms the scaling of relaxation time with the decay length that is reached by diffusion and degradation alone. Thus, the decay length of the steady-state profiles can be increased up to the millimeter range while staying within a biologically relevant relaxation times of days, see figure 7.

# 6. Discussion

## 6.1. Formation of long-range morphogen gradients by a signaling relay

In this paper, we explore the consequences of incorporating morphogen-mediated morphogen production into morphogen gradient formation. Previous models of morphogen gradient formation assume the release of morphogen from a spatially restricted source region and morphogen gradient formation by diffusion and degradation within the adjacent cell field [5, 12–16]. Motivated by recent observations in planarian flatworms (discussed in more detail below), our model incorporates a signaling relay in form of positive feedback of the morphogen on its own production across the entire cell field into the gradient formation process. This positive feedback loop consists of two elements that jointly mediate morphogen-mediated morphogen production. Morphogen production is positively regulated by intracellular morphogen-induced signaling activity, which in turn depends on the extracellular morphogen concentration, see figure 1 and equations (1) and (2). Within the cell field, morphogen-mediated morphogen production facilitates signal propagation between neighboring cells. Gradient initiation and orientation are set by a local morphogen source at the proximal boundary of the cell field. The source cells constitutively release morphogen in a signaling-independent manner, analogous to the role of the source cells in previous models [5, 12–16]. However, in our model the source cells additionally produce morphogen in a morphogen signaling-dependent manner, similar to all the cells in the cell field. While the signaling-independent morphogen production primes gradient formation analogous to the role of the source in previous gradient formation models, morphogen-mediated morphogen production in our model effectively extends the source region across the entire cell field. A further important distinction to gradient formation by diffusion and degradation is that in our model, the morphogen does not necessarily need to travel the entire length of the gradient, as the cell-to-cell signaling relay can provide an additional signal propagation route even in the absence of diffusion. As further elaborated below, signaling relay

mechanisms may thus extend the lengths and time scales over which patterns can be formed under physiological diffusion coefficients, while also providing potentially new mechanisms for pattern scaling.

## 6.2. The role of the feedback strength

The strength $\alpha$ of the morphogen-dependent morphogen production feedback loop is a key parameter in our model. It is assumed to be spatially homogeneous across the entire cell field. Note that despite this spatially homogeneous feedback strength, morphogen production may still be spatially modulated due to spatial differences in extracellular morphogen concentrations, see equations (1) and (2). Morphogen production controlled by a positive feedback loop thus forms a self-organized system that can generate signaling gradients and morphogen gradients.

The key element to achieve steady-state concentration gradients in the presence of a positive feedback loop is feedback saturation. In the absence of feedback saturation, there exists a critical point in the feedback strength at which the production due to the positive feedback would exceed the degradation and no steady-state would be reached. In the presence of feedback saturation, the system reaches a steady state at this critical feedback strength $\alpha^{crit}$ with several interesting consequences on the system's properties. Firstly, close to the critical feedback strength, the decay length of the gradient profile becomes maximal, see figure 3. Secondly, for feedback strengths above the critical one, gradient formation is inefficient, as the relaxation time stays high for $\alpha > \alpha^{crit}$, see figure 6, while the decay length decreases, see figures 3 and 7. Note that the effect of the signaling relay on both decay length and relaxation time is particularly pronounced close to the critical feedback strength, see figures 3 and 6. Thirdly, above the critical feedback strength, steady-state concentration profiles become rather flat with a high base level of morphogen due to the positive feedback and a shallow spatial profile, see figures D1–D4 in the appendix. We speculate that in the regime $\alpha > \alpha^{crit}$, the signaling relay may lead to traveling fronts, similar to those reported in reference [48].

Given that the focus of this study is on the efficient formation of long-range morphogen and signaling gradients, we focus our discussion on feedback strengths $\alpha \le \alpha^{crit}$. In this regime, the signaling relay affords long decay lengths, see figures 3 and 4 and gradient formation on biologically relevant time scales of hours to days, see figures 6 and 7.

In summary, our results demonstrate that the formation of steady-state gradients is indeed possible in systems incorporating morphogen-mediated morphogen production. The strength of the feedback has important consequences on the decay length of the morphogen gradients formed, as well as on the time scale of their formation.
### 6.3. Signaling relay increases the decay length of morphogen and signaling gradients

An interesting feature of our signaling relay model is that it affords the extension of gradient decay lengths beyond the decay length set by diffusion and degradation alone. We make a conservative choice of the molecular diffusion coefficient of the morphogen of $1 \mu m^2 s^{-1}$ to take into account the slow diffusivity of lipid-modified molecules, and use a degradation rate of the morphogen of $10^{-3} s^{-1}$, like the one measured for Wg [14] or Fgf8 [15]. With these numbers, gradients formed by diffusion and degradation alone would reach a decay length of $31 \mu m$. In contrast, we show that with these same numbers, the signaling relay can generate gradients with decay lengths of up to about 950 $\mu m$, see figure 4. Note that, in order to achieve such decay lengths by diffusion and degradation alone, molecular life times of days would be required. In contrast, the signaling relay described here can build such long-range gradients using short-lived molecules with life times of minutes to hours.

This remarkable decay length-extension effect is due to the fact that morphogen-mediated morphogen production throughout the tissue uncouples the decay length of the gradient from the diffusion distance of individual morphogen molecules. The morphogen profile still follows a transport equation (equation (8)) which is characterized by effective transport coefficients (equations (9) and (11)) and an effective degradation rate (equation (10)). We can understand the increased gradient decay length in terms of these effective transport coefficients and degradation rate: the effective diffusion of the morphogen is increased compared to its molecular diffusion coefficient (equation (9)), and the effective degradation rate of the morphogen (equation (10)) is smaller than the molecular degradation rate. Both effects are due to the production of morphogen across the whole cell field in response to the signaling relay and lead to an increase in the decay length of the gradients. Importantly, the resulting signaling pattern can extend beyond the limits set by the molecular diffusion coefficient and degradation rate of the morphogen. Taken to its extreme, this means that with the signaling relay, morphogen gradients can form even in the limit of no molecular diffusion (e.g., the morphogen remaining associated with the surface of the producing cell), exclusively due to the positive feedback and the cell-to-cell relay.

In addition to effective diffusion and effective degradation, the production of morphogen across the whole cell field can lead to an effective drift of morphogen through the system (equation (11)) if the morphogen is preferentially secreted towards the distal end of the system, $p > \frac{1}{2}$. This drift emerges due to the secretion polarity even though individual molecules do not exhibit drift. Biological mechanisms

that could achieve such polarized secretion include planar cell polarity [49]. The effective drift further increases the decay length of the gradients formed, see figure 4. For a diffusion coefficient of $1 \mu m^2 s^{-1}$ and a degradation rate of the morphogen of $10^{-3} s^{-1}$, the signaling relay reaches decay lengths of up to 1770 $\mu m$ in the presence of secretion polarity, see figure 4. Thus, secretion polarity generates an additional lever to increase the range of gradients formed by a signaling relay.

Taken together, a signaling relay allows to reach decay lengths of the order of millimeters even for slowly diffusing molecules. This is one to two orders of magnitude larger than what is achieved by diffusion and degradation alone, see figures 3 and 4.

### 6.4. Efficient formation of long-range morphogen gradients

Apart from the decay length of a morphogen gradient, the time it takes to form the profile is an important characteristic of any patterning mechanism, since the profile has to form on time scales relevant for the patterning process in question. The relaxation time of a morphogen gradient, defined as the slowest exponential relaxation time close to the steady state, is a good measure of the dynamics of the gradient. In general, increasing the decay length of gradients requires longer times compared to the formation of more short-range profiles. Accordingly, we find that for a cell-to-cell signaling relay the relaxation time slows down with increasing feedback strength $\alpha$, see figures 5 and 6 as the decay length increases, for $\alpha < \alpha^{crit}$, see figures 3 and 4. To quantify this trade-off between the decay length and the relaxation time, we analyzed the scaling between the decay length and the relaxation time, $\tau_{max} \propto \lambda^{\nu}$ for $\alpha < \alpha^{crit}$. In the signaling relay-model, the relaxation time scales with the decay length with an exponent of $\nu \simeq 2$ in the absence of secretion polarity, see figure 7. This is the same scaling relationship observed for gradient formation by diffusion and degradation. However, in the signaling relay-model in the presence of secretion polarity towards the distal end, $p > \frac{1}{2}$ this scaling exponent is smaller than 2, see figure 7. We note that in order to generate long-range gradients on efficient time scales, the degradation rate of the intracellular signaling molecule $k_B$ must be larger than the degradation rate of the morphogen $k_A$.

Taken together, positive feedback combined with polarized secretion outperforms the scaling relationship of the relaxation time with the decay length reached by diffusion and degradation alone. This way, the decay lengths of the steady-state morphogen and signaling gradients can be increased into the millimeter range while staying within a biologically relevant relaxation time of hours to days, see figure 7, for a physiologically relevant choice of parameters.
## 6.5. Scaling and robustness of gradients formed by a signaling relay

The signaling relay does not only lead to an increased decay length compared to a gradient formed by diffusion and degradation alone, it additionally provides a new way of regulating the decay length of the gradient: the signaling relay modulates the effective transport coefficients and the effective degradation rate that together determine the decay length of the gradient. Thus, the decay length can be tuned by exclusively changing the feedback strength.

The presence of a critical point in the signaling relay provides a sensitive response to the feedback strength: small relative changes in the feedback strength $\alpha$ in the vicinity of the critical point can generate changes in the decay length of more than one order of magnitude in our signaling relay-model, see figures 3 and 4. In contrast, in order to scale the decay length of a profile generated by diffusion and degradation across one order of magnitude, either the diffusion coefficient or the degradation rate would have to change by two orders of magnitude. The sensitive feedback in our model could be relevant for scaling, i.e. the increase of the decay length proportional to system size in a growing tissue [19, 50]. In particular, when scaling is achieved through a control of the feedback strength linked to the size of the system [50–52] a sensitive feedback could facilitate robust and reliable scaling. Indeed, a sensitive feedback implies an efficient response in adjusting the gradient to small changes in system size.

However, a sensitive feedback may also raise questions of the robustness of the system to small perturbations. Those questions are based on the fact that biological processes are inherently noisy. It has been suggested that robustness then implies an insensitivity to changes in parameters [53, 54]. Interestingly, careful experimental observations of the bicoid gradient in *Drosophila* embryos show that gradient formation is very precise and in fact quantitatively reproducible between embryos, despite being sensitive to the mRNA levels of the morphogen, i.e. to the production rate [55]. This suggests that gradient formation is precisely controlled to ensure robustness at every step of gradient formation [55]. In the context of our model of gradient formation by a signaling relay, robustness could be achieved through an efficient control of a sensitive feedback.

## 6.6. Conclusion and outlook

The concept of patterning by morphogen gradients has mostly been established in few molecularly tractable developmental model systems, e.g. *Drosophila* [7, 8, 12, 14, 19, 55] or zebrafish embryos [15, 16]. In these systems, as in other embryonic systems, patterning occurs on small length scales and is well explained by diffusion and degradation [14, 15]. However, patterning processes are not restricted to embryonic development and the associated small length

scales; patterning processes are also required in adult animals, with regeneration as one notable example. Limb regeneration in adult salamanders [56] or whole body regeneration inclusive of the brain in planarian flatworms [25–31, 57, 58] not only entail spectacular examples of complex tissue patterning processes, but also involve length scales up to the centimeter range. As we show in figure 3 the decay length based on diffusion and degradation is much shorter than in the presence of a signaling relay, suggesting that the biological utility of morphogen gradients formed by diffusion and degradation is strongly diminished at large scales. Our signaling relay mechanism via morphogen-mediated morphogen expression thus provides an attractive concept for bridging scales, see figures 4 and 7.

Our model is inspired by the posterior-to-anterior Wnt signaling gradient in planarian flatworms that patterns the planarian anteroposterior axis and that is established via Wnt-ligand expression gradients [25–31, 46]. Moreover, the expression gradient-forming Wnt ligands [27–29, 31, 46] are themselves expressed in a Wnt-signaling dependent manner [28, 29, 31] within the sheet-like body wall musculature [59–61]. Gradient formation is likely initiated by Wnt1, which is expressed at the very tail tip [27, 46] in an at least partially Wnt-independent manner [28, 31]. Jointly, these data motivate signal propagation across cell fields (sheet-like body wall musculature) from a partially signaling-independent source at the boundary (*wnt1* expression at the tail tip) via morphogen-dependent-morphogen expression (*Wnt*-dependent *wnt* expression of multiple Wnt-ligands). Though many mechanistic aspects of the signaling relay remain to be experimentally demonstrated, the uncoupling of pattern decay length from the diffusion distance of individual molecules provides an attractive proposition for generating the up to centimeter-scale *wnt* expression and signaling gradients that are observed in planarians.

Interestingly, morphogen expression gradients, that are a general hallmark of our signaling relay model, are also observed in other systems. Wnt expression gradients of different lengths emanate from the hypostome in adult Hydra [62]. They are further present during bud formation and reform in a sequential manner during head regeneration [62]. The hypostome-restricted expression of other *wnt* genes might be related to the signaling-independent source [62, 63]. Moreover, *wnt* expression gradients are observed in overlapping domains during embryonic development of the sea anemone *Nematostella vectensis* [64]. Furthermore, a signaling relay might also be involved in TGF-β signaling-dependent mesoderm induction during early embryonic development of Xenopus, as evidenced by mRNA injection experiments: the signal is only propagated if the clone of the signal-injected blastomer touches the one of the ligand-injected blastomer [65]. We predict that other
examples of signaling relays may be discovered in biological contexts that require patterning over large length scales, such as adult body part regeneration or in continuously growing organisms. Moreover, it will be interesting to revisit cases in which morphogen gradients formed by diffusion and degradation are believed to organize patterning. For instance, it could be investigated whether the shape of the Wg (the main *Drosophila* Wnt) gradient in the *Drosophila* wing imaginal disk changes in the absence of the *Drosophila* β-catenin homologue Amardillo. wg actually has a TCF binding site in its regulatory region [66] and could thus be regulated by intracellular signaling (Amardillo) levels. In fact, evidence for positive feedback of Wg on its own expression exists in intersegmental patterning during *Drosophila* development [67]. Beyond long-range pattern formation, pattern formation via poorly diffusible morphogens provides a further conceptual utility of signaling relays/morphogen-mediated morphogen expression in biological systems. Interestingly, both Wnt [34] and Hh [32, 33] as two of the most prominent morphogens are modified by lipid modifications that render them poorly diffusible [35–37]. In our model of a signaling relay, pattern formation can occur even in the limit of non-diffusible signals, i.e., the signal remaining associated with the producing cell, see equation (9). One observation that might be particularly interesting to revisit in this context is the finding that membrane-tethered Wg can still mediate patterning of the developing *Drosophila* wing imaginal disk [38]. This raises the question whether a Wg-induced Wg signaling relay could operate at low concentration levels.

Finally, the facile scalability of signaling relay patterns via modulation of the feedback strength $\alpha$, see figure 3, constitutes an intriguing property of our model. Pattern scaling remains an unsolved problem in biology, as on the one hand side the concomitance of patterning with growth necessitates scalable patterns, yet it is still debated how morphogen diffusion and degradation parameters could be modified dynamically, particularly in a self-organized manner [21, 50, 68]. Here, the adjustment of pattern length scales via the tuning of intracellular feedback loops as common element of signal transduction pathways constitutes a new concept that warrants experimental testing.

## Acknowledgments

We thank Jonathan Bauermann, Charlie Duclut, Steffen Werner, Daniel Aguilar-Hidalgo, Alexander Mietke, Felix Meigel, Hanh Vu, Mario Ivanković and James Cleland for insightful discussions.

## Data availability statement

The data that support the findings of this study are available upon reasonable request from the authors.

## Appendix A. Continuum theory of signaling gradient formation

In this section, we recapitulate the continuum theory of gradient formation including drift. The diffusion-degradation model for gradient formation that has successfully been applied to explain the Dpp gradient in the wing imaginal disc [14] is contained in this model as the special case of no-drift. In particular, the dynamics of morphogen concentration $C$ are given by secretion in a local source of width $S$ with rate $s$, diffusion with diffusion coefficient $D$, degradation with rate $k$, and drift or advective transport of molecule concentration with drift speed $v$,

$$ \partial_t C = D \partial_x^2 C - \partial_x(vC) - kC + s\theta(S-x), \quad (A.1) $$

where $x$ denotes the spatial coordinate in the dimension along which the gradient is formed, and $\theta$ denotes the Heaviside function, where $\theta(y) = 0$ for $y < 0$ and $\theta(y) = 1$ for $y \ge 0$. In equation (A.1), $0 < x \le L$ with a tissue of size $L$ and the source at the proximal boundary. Positive values of $v$ indicate transport in positive $x$-direction. Thus, for a profile decreasing in positive $x$-direction, $\partial_x C < 0$, the drift term leads to drift down the concentration gradient.

This dynamics leads to steady-state concentration profiles of the form

$$ C^* = \begin{cases} C_1^{\text{in}} \exp(x/\lambda_1) + C_2^{\text{in}} \exp(-x/\lambda_2) & 0 \le x \le S \\ C_1^{\text{out}} \exp(x/\lambda_1) + C_2^{\text{out}} \exp(-x/\lambda_2) & S < x \le L \end{cases}, \quad (A.2) $$

where the decay lengths are given by

$$ \lambda_1 = \frac{2D}{(v^2 + 4Dk)^{1/2} + v}, \quad (A.3) $$

$$ \lambda_2 = \frac{2D}{(v^2 + 4Dk)^{1/2} - v}, \quad (A.4) $$

Note that in the absence of drift, $v = 0$, this expression simplifies to the well-known decay length of a diffusion-degradation system that we denote by the index $v = 0$:

$$ \lambda_{v=0} = \sqrt{\frac{D}{k}}. \quad (A.5) $$

For no-flux boundary conditions, $-D\partial_x C(0) + vC(0) = -D\partial_x C(L) + vC(L) = 0$, and differentiability, i.e. matching value and flux at point $S$, where the source and the non-source regions meet, $C_1^{\text{out}}$ is very small and goes to zero as the system size $L$ increases. Therefore, the profile outside of the source
<figure><img src="image_14.png" /><figcaption>Figure A1. Comparison of the continuous diffusion-degradation model with and without drift. Comparison of steady-state profile decay lengths of the diffusion degradation model with drift (λ) to the one without drift (λv=0) (a), of the respective relaxation times (b), and the trade-off between decay length and relaxation time (λ/τ and λv=0/τv=0, respectively) (c). All plots as a function of a dimensionless parameter relating the influence of drift with speed v that is varied to the one of diffusion with diffusion coefficient D and degradation with degradation rate k, that are set to values measured in the fly wing disc (D = 0.1 µm², k = 10⁻³ s⁻¹ [14]). The dashed lines indicate the absence of drift, v = 0.</figcaption></figure>

region is given by an exponentially decreasing profile with the decay length λ := λ₂.

Comparing the decay lengths in the presence and absence of drift, we note that the decay length is increased in a system with drift for positive drift speeds, see figure A1(a). The relaxation time of gradient formation, defined as the relaxation time of the most slowly relaxing eigenmode, that is a good measure of how long the system takes to reach steady state, is given by

$$ \tau = \left( \frac{v^2}{4D} + D \left( \frac{\pi n}{L} \right)^2 + k \right)^{-1} . \quad (A.6) $$

Again, the time-scale of the diffusion-degradation model is contained as the case v = 0:

$$ \tau_{v=0} = \left( D \left( \frac{\pi n}{L} \right)^2 + k \right)^{-1} . \quad (A.7) $$

We can thus appreciate, that the system relaxes faster in the presence of drift, see figure A1(b) for a comparison in the limit of an infinitely large system (L → ∞).

There is a trade-off between the decay length of the steady-state profile and the relaxation time of gradient formation. That is, the formation of long-range steady state profiles takes increasingly long. In particular, in the absence of drift, the relaxation time scales as

$$ \tau_{v=0} \propto \lambda_{v=0}^2 . \quad (A.8) $$

As we discussed above, in the presence of drift, the steady-state decay length in a system with positive drift speeds is increased compared to a system without drift (figure A1(a)), and the relaxation time of gradient formation is decreased (figure A1(b)). Thus, the trade-off between the decay length of the steady state profile and the relaxation time, that can be quantified as the ratio of the decay length and the relaxation time λ/τ is markedly improved, see figure A1(c). This trade-off is a measure of how long it takes to form a gradient of a given decay length. Thus, a gradient with an equal decay length is reached faster in the presence of drift or a gradient with a longer decay length is reached in the same time.

While drift of morphogen concentration through the system may very well explain the formation of

signaling gradients on decay lengths in the millimeter range on biologically relevant time scales of hours to days, it leaves open the question of the origin of drift. Drift may arise due to concerted cilia beating or muscle contraction [43], but neither can explain the formation of long-range gradients in the context of regeneration in which the respective structures would have to be built first.

## Appendix B. Boundary conditions

We consider the case in which the system ends with an extracellular space 0 at the proximal end and an extracellular space N at the distal end, see figure 1. We thus need to specify the boundary conditions ∂_t a₀ and ∂_t a_N. We assume no diffusive flux at these boundaries. Moreover, we choose the volumes of the extracellular spaces at the system boundaries to be smaller than those in the bulk, in particular

$$ V_0 = (1 - p)V, \quad (B.1) $$

$$ V_N = pV, \quad (B.2) $$

where V denotes the volume of the extracellular spaces in the bulk, i.e. V_n = V for n = 1, ..., N − 1. This leads to the following dynamics at the system boundaries:

$$ \partial_t a_0 = \frac{D}{(1-p)} \frac{-a_0 + a_1}{\delta^2} - k_A a_0 + s_A + \alpha \frac{b_0}{c_B + b_0}, \quad (B.3) $$

$$ \partial_t a_N = \frac{D}{p} \frac{a_{N-1} - a_N}{\delta^2} - k_A a_N + \alpha \frac{b_{N-1}}{c_B + b_{N-1}}. \quad (B.4) $$

Note that with this choice of volumes at the boundaries, the secretion polarity in the production terms is balanced at the boundaries. To see this, consider that the cells produce individual molecules rather than concentrations and that thus the production rates measured in concentration per time depend on the volume of the extracellular spaces into which the produced molecules are secreted.
Note further that this choice of boundary conditions gives rise to flat profiles at the boundaries of the system, i.e. vanishing spatial derivative at the boundaries, see figure 2. Finally, note that in the presence of secretion polarity, $p \neq \frac{1}{2}$, this choice of boundary conditions is not the same as no-flux boundary conditions.

## Appendix C. Numerical steady-state solution

We can obtain the steady-state solution to the signaling relay model given by equations (1) and (2) numerically by starting from an initial condition and computing the time evolution for long times. We choose the average of the first ($a_0$ and $b_0$, respectively) and the last value ($a_N$ and $b_{N-1}$, respectively) of the analytical approximation to the steady-state solution as initial conditions. In order to obtain the concentration profiles of the morphogen and the intracellular signaling molecule at time $t$, $a^t$, and $b^t$, respectively:

$$ \mathbf{a}^t = \begin{bmatrix} a_0^t \\ \vdots \\ a_N^t \end{bmatrix}, \qquad \mathbf{b}^t = \begin{bmatrix} b_0^t \\ \vdots \\ b_{N-1}^t \end{bmatrix}, \qquad (C.1) $$

where the upper index refers to the time step and bold italic fond indicates vector notation, we need their concentration profiles at time $t-1$, $a^{t-1}$, and $b^{t-1}$. We use a forward (explicit) Euler method to compute the (non-linear) dynamics of the concentration profile of the intracellular signaling molecule. Specifically, we add the change in concentration during one time-step of length $\Delta t$ to the current concentration profile $b^{t-1}$:

$$ \mathbf{b}^t = \mathbf{b}^{t-1} + \Delta t \mathbf{G}(\mathbf{a}^{t-1}, \mathbf{b}^{t-1}), \qquad (C.2) $$

where the differential equation for the signaling activity $\mathbf{G}(\mathbf{a}^{t-1}, \mathbf{b}^{t-1})$ is defined according to equation (2):

$$ [\mathbf{G}(\mathbf{a}^{t-1}, \mathbf{b}^{t-1})]_n = s_B - k_B \frac{c_A}{c_A + a_n^{t-1} + a_{n+1}^{t-1}} b_n^{t-1}. \qquad (C.3) $$

To compute the dynamics of the morphogen, we use an implicit-explicit Euler method [69]. In particular, we evaluate the linear parts (diffusion and degradation, described by the matrix $N$) implicitly and the non-linear parts (i.e. the local, signaling-independent source, as well as the feedback, i.e. the signaling-dependent source, $\mathbf{F}(\mathbf{b}^{t-1})$) explicitly (see reference [69] for details on the method). This results in the following dynamics:

$$ \mathbf{a}^t = (1 - \Delta t \mathbf{N})^{-1} \cdot [\mathbf{a}^{t-1} + \Delta t \mathbf{F}(\mathbf{b}^{t-1})]. \qquad (C.4) $$

The linear part is given by diffusion and degradation:

$$ N_n = \left( 0 \dots 0 \frac{D}{\delta^2} - 2 \frac{D}{\delta^2} - k_A \frac{D}{\delta^2} 0 \dots 0 \right). \qquad (C.5) $$

The non-linear part of the dynamics of the signaling molecule, $\mathbf{F}(\mathbf{b}^{t-1})$, is given by the constant production in the local, signaling-independent source with rate $s_A$ and the positive feedback according to equation (3). This results in:

$$ [\mathbf{F}(\mathbf{b}^{t-1})]_{n=1,...,w-1} = s_A + p \alpha \frac{b_{n-1}^{t-1}}{c_B + b_{n-1}^{t-1}} + (1-p) \alpha \frac{b_n^{t-1}}{c_B + b_n^{t-1}}, \qquad (C.6) $$

$$ [\mathbf{F}(\mathbf{b}^{t-1})]_w = p s_A + p \alpha \frac{b_{w-1}^{t-1}}{c_B + b_{w-1}^{t-1}} + (1-p) \alpha \frac{b_w^{t-1}}{c_B + b_w^{t-1}}, \qquad (C.7) $$

$$ [\mathbf{F}(\mathbf{b}^{t-1})]_{n=w+1,...,N-1} = p \alpha \frac{b_{n-1}^{t-1}}{c_B + b_{n-1}^{t-1}} + (1-p) \alpha \frac{b_n^{t-1}}{c_B + b_n^{t-1}}. \qquad (C.8) $$

The boundary conditions are specified in the same way according to equations (B.3) and (B.4). This way, we can compute the dynamic solution to our model and in particular, for long times, the steady-state solution. We verify that the steady-state is reached by observing that the profiles no longer change with time.

## Appendix D. Analytical approximation to the steady state

In short, we find an analytical approximation to the solution of the steady-state equation of the morphogen (equation (1)) using equation (5)) by linearizing it around the piece-wise constant steady-state solution given by $a_c^{\text{in}}$ and $a_c^{\text{out}}$ inside ($n=0, \dots, w$) and outside ($n=w+1, \dots, N$) of the source region, respectively. This piece-wise constant solution is defined by the constants solving the steady-state equations of the bulk of the source and the non-source region, respectively. We describe this procedure in more detail in this section:

We analyze the steady-state of the morphogen concentration in order to determine the piece-wise constant steady-state solution. Consider that at steady state, $b_n^*$ is linearly dependent on $a_n^*$ (equation (5)). Inserting this relationship (equation (5)) into the steady state equation for $a_n^*$ (equation (1)), we obtain

$$ 0 = D \frac{a_{n-1} - 2a_n + a_{n+1}}{\delta^2} - k_A a_n + s_A \theta_{w-n} + \alpha p \frac{s_B}{k_B} \left( 1 + \frac{a_{n-1}^* + a_n^*}{c_A} \right) \times \left[ c_B + \frac{s_B}{k_B} \left( 1 + \frac{a_{n-1}^* + a_n^*}{c_A} \right) \right]^{-1} $$
$$ + \alpha(1-p)\frac{s_B}{k_B}\left(1+\frac{a_n^*+a_{n+1}^*}{c_A}\right) \times \left[c_B+\frac{s_B}{k_B}\left(1+\frac{a_n^*+a_{n+1}^*}{c_A}\right)\right]^{-1}. \quad (D.1) $$

Based on this equation, we define $a_c^{\text{in}}$ and $a_c^{\text{out}}$ as the constants solving equation (D.1) for $n = 1, \dots, w-1$ and $n = w+1, \dots, N-1$, respectively. We obtain

$$ a_c^{\text{in}} = \frac{1}{2} \left[ \frac{s_A}{k_A} + \frac{\alpha}{k_A} - \frac{c_A}{2} \left( \frac{c_B k_B}{s_B} + 1 \right) \right] + \sqrt{ \frac{1}{4} \left[ \frac{s_A}{k_A} + \frac{\alpha}{k_A} - \frac{c_A}{2} \left( \frac{c_B k_B}{s_B} + 1 \right) \right]^2 + \frac{c_A}{2k_A} \left[ s_A \left( \frac{c_B k_B}{s_B} + 1 \right) + \alpha \right] }, \quad (D.2) $$

$$ a_c^{\text{out}} = \frac{1}{2} \left[ \frac{\alpha}{k_A} - \frac{c_A}{2} \left( \frac{c_B k_B}{s_B} + 1 \right) \right] + \sqrt{ \frac{1}{4} \left[ \frac{\alpha}{k_A} - \frac{c_A}{2} \left( \frac{c_B k_B}{s_B} + 1 \right) \right]^2 + \frac{\alpha c_A}{2k_A} } \quad (D.3) $$

We then obtain the approximation to the steady-state solution by linearizing equation (D.1) around $a_c^{\text{in}}$ inside the source region for $n = 0, \dots, w$, and around $a_c^{\text{out}}$ outside of the source region for $n = w+1, \dots, N$, respectively. To this end, we express $a_n$ as a deviation from these constants:

$$ a_n = a_c^{\text{in}} + \delta a_n \quad \text{for: } n = 0, \dots, w, \quad (D.4) $$

$$ a_n = a_c^{\text{out}} + \delta a_n \quad \text{for: } n = w+1, \dots, N, \quad (D.5) $$

and linearize for small $\delta a_n$. We thus obtain a set of linear equations:

$$ \begin{aligned} \partial_t a_0^* = 0 &= \frac{D}{(1-p)} \frac{-\delta a_0^* + \delta a_1^*}{\delta^2} - k_A \delta a_0^* \\ &\quad + \alpha z^{\text{in}}(\delta a_0^* + \delta a_1^*), \quad (D.6) \\ \text{for } n \in [1, \dots, (w-1)]: \end{aligned} $$

$$ \begin{aligned} \partial_t a_n^* = 0 &= D \frac{\delta a_{n-1}^* - 2\delta a_n^* + \delta a_{n+1}^*}{\delta^2} - k_A \delta a_n^* \\ &\quad + p \alpha z^{\text{in}}(\delta a_{n-1}^* + \delta a_n^*) \\ &\quad + (1-p) \alpha z^{\text{in}}(\delta a_n^* + \delta a_{n+1}^*), \quad (D.7) \end{aligned} $$

$$ \begin{aligned} \partial_t a_w^* = 0 &= D \frac{-a_c^{\text{in}} + a_c^{\text{out}}}{\delta^2} - k_A a_c^{\text{in}} + p s_A \\ &\quad + p \alpha (c_A + 2a_c^{\text{in}}) \\ &\quad \times \left[ c_A + \frac{c_A c_B k_B}{s_B} + 2a_c^{\text{in}} \right]^{-1} \\ &\quad + (1-p) \alpha (c_A + a_c^{\text{in}} + a_c^{\text{out}}) \\ &\quad \times \left[ c_A + \frac{c_A c_B k_B}{s_B} + a_c^{\text{in}} + a_c^{\text{out}} \right]^{-1} \\ &\quad + D \frac{\delta a_{w-1}^* - 2\delta a_w^* + \delta a_{w+1}^*}{\delta^2} - k_A \delta a_w^* \\ &\quad + p \alpha z^{\text{in}}(\delta a_{w-1}^* + \delta a_w^*) \\ &\quad + (1-p) \alpha (\delta a_w^* + \delta a_{w+1}^*) \frac{c_A c_B k_B}{s_B} \\ &\quad \times \left[ c_A + \frac{c_A c_B k_B}{s_B} + a_c^{\text{in}} + a_c^{\text{out}} \right]^{-2}, \quad (D.8) \end{aligned} $$

$$ \begin{aligned} \partial_t a_{w+1}^* = 0 &= D \frac{a_c^{\text{in}} - a_c^{\text{out}}}{\delta^2} - k_A a_c^{\text{out}} \\ &\quad + p \alpha (c_A + a_c^{\text{in}} + a_c^{\text{out}}) \\ &\quad \times \left[ c_A + \frac{c_A c_B k_B}{s_B} + a_c^{\text{in}} + a_c^{\text{out}} \right]^{-1} \\ &\quad + (1-p) \alpha (c_A + 2a_c^{\text{out}}) \\ &\quad \times \left[ c_A + \frac{c_A c_B k_B}{s_B} + 2a_c^{\text{out}} \right]^{-1} \\ &\quad + D \frac{\delta a_w^* - 2\delta a_{w+1}^* + \delta a_{w+2}^*}{\delta^2} \\ &\quad - k_A \delta a_{w+1}^* + p \alpha (\delta a_w^* + \delta a_{w+1}^*) \\ &\quad \times \frac{c_A c_B k_B}{s_B} \\ &\quad \times \left[ c_A + \frac{c_A c_B k_B}{s_B} + a_c^{\text{in}} + a_c^{\text{out}} \right]^{-2} \\ &\quad + (1-p) \alpha z^{\text{out}} (\delta a_{w+1}^* + \delta a_{w+2}^*), \quad (D.9) \end{aligned} $$

$$ \begin{aligned} \partial_t a_n^* = 0 &= D \frac{\delta a_{n-1}^* - 2\delta a_n^* + \delta a_{n+1}^*}{\delta^2} - k_A \delta a_n^* \\ &\quad + p \alpha z^{\text{out}}(\delta a_{n-1}^* + \delta a_n^*) + (1-p) \\ &\quad \times \alpha z^{\text{out}}(\delta a_n^* + \delta a_{n+1}^*), \quad (D.10) \end{aligned} $$

$$ \begin{aligned} \partial_t a_N^* = 0 &= \frac{D}{p} \frac{\delta a_{N-1}^* - \delta a_N^*}{\delta^2} - k_A \delta a_N^* + \alpha z^{\text{out}} \\ &\quad \times (\delta a_{N-1}^* + \delta a_N^*) \quad (D.11) \end{aligned} $$

where

$$ z^{\text{in}} = \frac{c_A c_B k_B}{s_B} \left( c_A + \frac{c_A c_B k_B}{s_B} + 2a_c^{\text{in}} \right)^{-2}, \quad (D.12) $$

$$ z^{\text{out}} = \frac{c_A c_B k_B}{s_B} \left( c_A + \frac{c_A c_B k_B}{s_B} + 2a_c^{\text{out}} \right)^{-2}. \quad (D.13) $$

Note that this linearization (equation (D.10)) defines $z^{\text{out}}$ (equation (12)), repeated here (equation (D.13))
<figure><img src="image_15.png" /><figcaption>**Figure D1.** Steady-state concentration profiles of the morphogen (a) and the signaling activity (b) for indicated values of the feedback strength $\alpha$. The signaling-independent source region is shaded in red. Circles denote the numerically obtained steady-state solution to the non-linear equations, lines indicate the analytical approximation to the steady state. Parameters: table H1, set B in the appendix.</figcaption></figure>
for convenience. Note further that the term for the signaling relay-dependent morphogen production by cell w at steady-state is given by $b_w^*/(c_B + b_w^*)$. When re-writing this using equation (5), it contains $a_w^*$ and $a_{w+1}^*$. Therefore, it needs to be linearized around the sum of $a_c^{\text{in}}$ and $a_c^{\text{out}}$. Thus, for $a_w^*$ and $a_{w+1}^*$, the constant term does not vanish (see equations (D.8) and (D.9)) as it contains both $a_c^{\text{in}}$ and $a_c^{\text{out}}$.

We can solve these linearized steady-state equations using an exponential ansatz and obtain the solution given in equation (6). The constants in equation (6) are defined by the boundary conditions at $n=0$ and $n=N$, as well as by the source/non-source interface at $n=w$ and $n=w+1$, equations (D.6), (D.8), (D.9), and (D.11). Plugging the exponential ansatz (6) into equations (D.6), (D.8), (D.9), and (D.11), we obtain:

$$ 
\begin{gathered}
\begin{pmatrix} \partial_t a_0^* \\ \partial_t a_w^* \\ \partial_t a_{w+1}^* \\ \partial_t a_N^* \end{pmatrix} = 0 = \begin{pmatrix} M_{11} & M_{12} & 0 & 0 \\ M_{21} & M_{22} & M_{23} & M_{24} \\ M_{31} & M_{32} & M_{33} & M_{34} \\ 0 & 0 & M_{43} & M_{44} \end{pmatrix} \\ \cdot \begin{pmatrix} C_1^{\text{in}} \\ C_1^{\text{out}} \\ C_2^{\text{out}} \\ C_2^{\text{in}} \end{pmatrix} + \begin{pmatrix} s_1 \\ s_2 \\ s_3 \\ s_4 \end{pmatrix} \qquad (D.14) \\ = M \cdot c + s. \qquad (D.15)
\end{gathered}
$$

Thus, the constants in c are given as:

$$ \Leftrightarrow c = M^{-1} \cdot (-s). \qquad (D.16) $$

We obtain the decay lengths by expanding the exponents in $1/\lambda$. This expansion is valid for decay lengths much larger than one cell, i.e. $\lambda \gg 1$, and thus the long decay lengths we are interested in. We find

$$ \lambda_1^{\text{in}} \approx \frac{\frac{2D}{\delta^2} + \alpha z^{\text{in}}}{\sqrt{[(2p-1)\alpha z^{\text{in}}]^2 + 2\left[\frac{2D}{\delta^2} + \alpha z^{\text{in}}\right][k_A - 2\alpha z^{\text{in}}] + (2p-1)\alpha z^{\text{in}}}}, \qquad (D.17) $$

$$ \lambda_2^{\text{in}} \approx \frac{\frac{2D}{\delta^2} + \alpha z^{\text{in}}}{\sqrt{[(2p-1)\alpha z^{\text{in}}]^2 + 2\left[\frac{2D}{\delta^2} + \alpha z^{\text{in}}\right][k_A - 2\alpha z^{\text{in}}] - (2p-1)\alpha z^{\text{in}}}}, \qquad (D.18) $$

$$ \lambda_1^{\text{out}} \approx \frac{\frac{2D}{\delta^2} + \alpha z^{\text{out}}}{\sqrt{[(2p-1)\alpha z^{\text{out}}]^2 + 2\left[\frac{2D}{\delta^2} + \alpha z^{\text{out}}\right][k_A - 2\alpha z^{\text{out}}] + (2p-1)\alpha z^{\text{out}}}}, \qquad (D.19) $$

$$ \lambda_2^{\text{out}} \approx \frac{\frac{2D}{\delta^2} + \alpha z^{\text{out}}}{\sqrt{[(2p-1)\alpha z^{\text{out}}]^2 + 2\left[\frac{2D}{\delta^2} + \alpha z^{\text{out}}\right][k_A - 2\alpha z^{\text{out}}] - (2p-1)\alpha z^{\text{out}}}}. \qquad (D.20) $$

We can compare this approximation to the steady-state profiles obtained numerically, see figures D1–D4. Note that the analytical approximation to the steady state is particularly good for feedback levels away from the critical feedback strength. Note further that close to the critical feedback strength, the approximation is better for $b_n/c_B \ll 1$, compare figure D1 ($c_B = 1.66 \times 10^8$ nM)

and figure D3 ($c_B = 1.66 \times 10^3$ nM). Finally, note that for $\alpha \gg \alpha^{\text{crit}}$, the positive feedback on the production $f(b_{n-1}, b_n)$ (equation (3)) essentially reduces to an additional constant production with strength $\alpha$ since in this case, $b_n \gg c_B$. Therefore, the concentrations are elevated throughout the system and we need to look at a baseline-corrected profile in order to see the decay length of the profile, see
<figure><img src="image_16.png" /><figcaption>**Figure D2.** Baseline-corrected steady-state concentration profiles of the morphogen (a) and the signaling activity (b) for indicated values of the feedback strength $\alpha$. In order to better visualize the gradient shape of the profiles depicted in figure D1, we subtracted the contractions $a_N$ (a) and $b_{N-1}$ (b), respectively from the entire profile. The signaling-independent source region is shaded in red. Circles denote the numerically obtained steady-state solution to the non-linear equations, lines indicate the analytical approximation to the steady state. Parameters: table H1, set B in the appendix.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure D3.** Steady-state concentration profiles of the morphogen (a) and the signaling activity (b) for indicated values of the feedback strength $\alpha$. The signaling-independent source region is shaded in red. Circles denote the numerically obtained steady-state solution to the non-linear equations, lines indicate the analytical approximation to the steady state. Parameters: table H1, set C in the appendix, $p = 0.50$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure D4.** Baseline-corrected steady-state concentration profiles of the morphogen (a) and the signaling activity (b) for indicated values of the feedback strength $\alpha$. In order to better visualize the gradient shape of the profiles depicted in figure D3, we subtracted the contractions $a_N$ (a) and $b_{N-1}$ (b), respectively from the entire profile. The signaling-independent source region is shaded in red. Circles denote the numerically obtained steady-state solution to the non-linear equations, lines indicate the analytical approximation to the steady state. Parameters: table H1, set C in the appendix, $p = 0.50$.</figcaption></figure>
figures D2 and D4. In this case, the decay length of the profile is dependent on the diffusion and the degradation, whereas the feedback contributes to overall high levels of morphogen and signaling molecule.

## Appendix E. Linearized dynamic equations

In order to calculate the effective transport coefficients and loss rate, defined by equation (8), we first determine eigenmodes of the dynamics linearized around the homogeneous steady-state solution given by $a_c^{in}$, $a_c^{out}$, $b_c^{in}$, $b_c^{out}$, see appendix D. The $b_c^{in/out}$ are obtained from the $a_c^{in/out}$ using equation (5). This results in:

$$b_c^{in} = \frac{s_B}{c_A k_B} (c_A + 2a_c^{in}) \quad (E.1)$$

$$b_c^{out} = \frac{s_B}{c_A k_B} (c_A + 2a_c^{out}). \quad (E.2)$$

We discuss the dynamics outside of the source region. The dynamics inside of the source region can be obtained in the same way using the respective constants $a_c^{in}$ and $b_c^{in}$. Outside of the source region, the dynamics linearized around the homogeneous steady state solution is given by

$$\partial_t \delta a_n \approx D \frac{\delta a_{n-1} - 2\delta a_n + \delta a_{n+1}}{\delta^2} - k_A \delta a_n + p\alpha \times \frac{c_B \delta b_{n-1}}{(c_B + b_c^{out})^2} + (1-p)\alpha \frac{c_B \delta b_n}{(c_B + b_c^{out})^2}, \quad (E.3)$$

$$\begin{aligned} \partial_t \delta b_n &\approx s_B - k_B \frac{c_A}{c_A + 2a_c^{\text{out}}} b_c^{\text{out}} - k_B \frac{c_A}{c_A + 2a_c^{\text{out}}} \delta b_n \\ &\quad + k_B \frac{c_A}{(c_A + 2a_c^{\text{out}})^2} b_c^{\text{out}} (\delta a_n + \delta a_{n+1}) \\ &= \frac{1}{c_A + 2a_c^{\text{out}}} [s_B(\delta a_n + \delta a_{n+1}) - c_A k_B \delta b_n], \end{aligned} \quad (E.4)$$

where we used equation (E.2) in the last step. We use the ansatz

$$\delta a_n(t) = A \exp(iq\delta n) \exp(-i\omega t), \quad (E.5)$$

$$\delta b_n(t) = B \exp(iq\delta n) \exp(-i\omega t) \quad (E.6)$$

for the eigenmodes of the system, where $q$ specifies an inverse wavelength, $\delta$ denotes the cell width, $\omega$ denotes a relaxation rate, $t$ specifies the time, $n$ the spatial index of the cell or extracellular space, and $A, B \in C$ denote the amplitudes of the eigenmodes that can be $q$-dependent. With this ansatz, the dynamics simplify to the eigenvalue problem

$$-i\omega \begin{pmatrix} A \\ B \end{pmatrix} = L \begin{pmatrix} A \\ B \end{pmatrix}, \quad (E.7)$$

where

$$L = \begin{pmatrix} \frac{D}{\delta^2} [\exp(-iq\delta) - 2 + \exp(iq\delta)] - k_A & \frac{\alpha c_B}{(c_B + b_c^{\text{out}})^2} [p \exp(-iq\delta) + 1 - p] \\ \frac{s_B}{c_A + 2a_c^{\text{out}}} [1 + \exp(iq\delta)] & - \frac{c_A k_B}{c_A + 2a_c^{\text{out}}} \end{pmatrix}. \quad (E.8)$$

From the resulting characteristic polynomial, we obtain two relaxation rates $\omega_1$ and $\omega_2$ which depend on the wave number, i.e. the inverse wavelength $q$:

$$\omega_{1,2} = \frac{1}{2} i \gamma \pm \frac{1}{2} \sqrt{i^2 \gamma^2 + \eta}, \quad (E.9)$$

$$\gamma = \frac{D}{\delta^2} [\exp(-iq\delta) - 2 + \exp(iq\delta)] - k_A - \frac{c_A k_B}{c_A + 2a_c^{\text{out}}}, \quad (E.10)$$

$$\eta = -4C \left\{ \exp(iq\delta) \left[ \frac{D}{\delta^2} + \alpha z^{\text{out}}(1-p) \right] + \exp(-iq\delta) \left[ \frac{D}{\delta^2} + \alpha z^{\text{out}} p \right] - \frac{2D}{\delta^2} - k_A + \alpha z^{\text{out}} \right\}, \quad (E.11)$$

where $z^{\text{out}}$ is defined by equation (12) and $C = (c_A k_B)/(c_A + 2a_c^{\text{out}})$. We aim to analyze the dynamics in the vicinity of the steady state. The steady state of
the system is contained in these dynamics as the limit of $\omega = 0$. We see that $\omega_2 = \frac{1}{2}i\gamma - \frac{1}{2}\sqrt{i^2\gamma^2 + \eta} = 0$ if $\eta = 0$. Thus, the steady-state and therefore the slow dynamics are contained in the relaxation rate $\omega_2$. We hence focus our further analysis on this relaxation rate. We identify the two roots of $\omega_2$ determined by $\eta = 0$. We start by expanding the exponential

functions for small values of $iq\delta$, i.e. long wavelengths, using

$$ \exp(iq\delta) \approx 1 + iq\delta + \frac{(iq\delta)^2}{2}. \qquad (E.12) $$

We obtain the approximations

$$ \gamma \approx \frac{D}{\delta^2}(iq\delta)^2 - k_A - C, \qquad (E.13) $$

$$ \eta \approx -4C \left[ (iq\delta)^2 \left( \frac{D}{\delta^2} + \frac{\alpha z^{\text{out}}}{2} \right) + (iq\delta) \left( (1-2p)\alpha z^{\text{out}} \right) - k_A + 2\alpha z^{\text{out}} \right]. \qquad (E.14) $$

Based on these, we can then calculate the inverse steady-state decay length $iq_{1,2}\delta$ defined by the steady-state condition $\eta = 0$. We obtain

$$ iq_{1,2}\delta \approx \frac{(2p-1)\alpha z^{\text{out}} \pm \sqrt{[(2p-1)\alpha z^{\text{out}}]^2 - 4\left(\frac{D}{\delta^2} + \frac{\alpha z^{\text{out}}}{2}\right)\left(-k_A + 2\alpha z^{\text{out}}\right)}}{2\left(\frac{D}{\delta^2} + \frac{\alpha z^{\text{out}}}{2}\right)}. \qquad (E.15) $$

Note that $q_{1,2}$ are equivalent to the steady-state decay lengths obtained in appendix D, equations (D.19) and (D.20):

$$ \lambda_{\text{out}}^{1} = \frac{1}{iq_{1}\delta}, \qquad (E.16) $$

$$ \lambda_{\text{out}}^{2} = -\frac{1}{iq_{2}\delta}, \qquad (E.17) $$

where the difference in sign is due to the definition of signs in equation (6). We can then expand $\omega_2$ around its roots, i.e. the inverse steady-state decay lengths $q_{1,2}$. To this end, we define

$$ x_1 = iq\delta - iq_1\delta, \qquad (E.18) $$

$$ x_2 = iq\delta - iq_2\delta. \qquad (E.19) $$

We can then express $\eta$ as

$$ \eta \approx -4C \left[ \frac{D}{\delta^2} + \frac{\alpha z^{\text{out}}}{2} \right] x_1 x_2, \qquad (E.20) $$

and $\gamma$ as

$$ \gamma \approx \frac{D}{\delta^2}x_1x_2 + \frac{D}{\delta^2}[(iq_1\delta + iq_2\delta)](iq\delta) - D(iq_1\delta)(iq_2\delta) - k_A - C. \qquad (E.21) $$

Note that the steady state now corresponds to $x_1 = 0$ or $x_2 = 0$. We aim to expand $w$ for wavenumbers in the vicinity of the steady-state, i.e. small $x_1x_2$. The hydrodynamic limit can be further approximated considering $x_2$ small such that $iq\delta \approx iq_2\delta$. With these simplifications, $\gamma$ can be approximated as

$$ \gamma \approx \frac{D}{\delta^2}x_1x_2 + \frac{D}{\delta^2}[(iq_1\delta + iq_2\delta)](iq_2\delta) - D(iq_1\delta)(iq_2\delta) - k_A - C. \qquad (E.22) $$

Using equations (E.18)–(E.22), we can then expand $\omega_2$ to lowest order in $x_1x_2$ as

$$ w_2 \approx i\phi \left( \frac{D}{\delta^2} + \frac{\alpha z^{\text{out}}}{2} \right) x_1 x_2, \qquad (E.23) $$

where

$$ \phi = \frac{C \left( \frac{D}{\delta^2} + \frac{\alpha z^{\text{out}}}{2} \right)^2}{\frac{D}{\delta^2} \alpha z^{\text{out}} (2p-1) \psi + \left( \frac{D}{\delta^2} + \frac{\alpha z^{\text{out}}}{2} \right) \left[ \frac{D}{\delta^2} (4\alpha z^{\text{out}} + C) + \alpha z^{\text{out}} (k_A + C) \right]}, \qquad (E.24) $$

where

$$ \psi = -\alpha z^{\text{out}}(2p-1) + \sqrt{[(2p-1)\alpha z^{\text{out}}]^2 - 4\left(\frac{D}{\delta^2} + \frac{\alpha z^{\text{out}}}{2}\right)(-k_A + 2\alpha z^{\text{out}})}. \qquad (E.25) $$
With

$$x_1x_2 = (iq\delta)^2 - (iq\delta_1 + iq_2\delta)(iq\delta) + (iq_1\delta)(iq_2\delta), \quad (\text{E.26})$$

$$i\omega_2 \approx \phi \left\{ \left[ -\left( \frac{D}{\delta^2} + \frac{\alpha z^{\text{out}}}{2} \right) \right] (iq\delta)^2 + [(2p-1)\alpha z^{\text{out}}] (iq\delta) + (k_A - 2\alpha z^{\text{out}}) \right\} \quad (\text{E.27})$$

equation (E.27) can be compared to the dispersion-relation of a convection-diffusion-degradation equation (equation (8)) which reads

$$i\omega = D_{\text{eff}}q^2 + iv_{\text{eff}}q + k_{\text{eff}}. \quad (\text{E.28})$$

By comparing equations (E.27) and (E.28), we can identify the effective transport coefficients and loss rate as the coefficients of powers of q. They are given in equations (9)-(11).

## Appendix F. Decay length based on numerical steady-state solution

In order to compute the decay length based on the numerical steady-state solution, we assume that the profile outside of the source region follows an exponential decay added on to a constant plateau or baseline. To reveal its decay length, we numerically compute the second derivative of the numerical steady-state profile. We obtain a baseline-corrected numerical steady-state profile by subtracting the baseline or plateau value (given by the value at the distal end of the numerical steady-state profile) from all points of the profile. Subsequently, we point-wise divide this baseline-corrected profile by the second derivative of the profile. We take the square root of these values and use their median of in the non-source region as the decay length based on the numerical steady-state profile.

## Appendix G. Relaxation times

In order to obtain the relaxation modes of the signaling relay, we analyze the dynamics of the system close to the steady state, see equations (19) and (20). We approximate the dynamics of the system close to the steady state by linearizing the dynamic equations of the morphogen (equation (1)) and the intracellular signaling molecule concentration (equation (2)) around the steady-state solution of the system given by $a_n^*$ and $b_n^*$. We obtain:

$$\partial_t \delta a_n \approx D \frac{\delta a_{n-1} - 2\delta a_n + \delta a_{n+1}}{\delta^2} - k_A \delta a_n + p \alpha c_B \frac{\delta b_{n-1}}{(c_B + b_{n-1}^*)^2} + (1-p) \alpha c_B \frac{\delta b_n}{(c_B + b_n^*)^2}, \quad (\text{G.1})$$

$$\partial_t \delta b_n \approx \frac{1}{c_A + a_n^* + a_{n+1}^*} \times [s_B(\delta a_n + \delta a_{n+1}) - c_A k_B \delta b_n] \quad (\text{G.2})$$

in the bulk of the system. The linearized dynamic equations at the boundaries can be obtained in the same way.

We can obtain the steady-state profiles $a_n^*$ and $b_n^*$ numerically, see appendix C. Alternatively, we can use the approximation to the steady-state profiles discussed in section 3.1, equation (6), using equations (5) and (6) as $b_n^*$ and $a_n^*$, respectively. The approach for analyzing the dynamics stays the same, the only difference is which profiles are used to linearize the dynamics.

Either way, the resulting system of linear equations can be expressed in matrix form as:

$$\partial_t \delta \mathbf{x} = \mathbf{M} \cdot \delta \mathbf{x}, \quad \text{where} \quad \delta \mathbf{x} = \begin{bmatrix} \delta a_0 \\ \delta b_0 \\ \vdots \\ \delta a_{N-1} \\ \delta b_{N-1} \\ \delta a_N \end{bmatrix}, \quad (\text{G.3})$$

and M describes the dynamics of the system according to equations (G.1) and (G.2) and the respective linearized equations at the boundaries. The dynamic solution to this system is given by the sum of all eigenmodes:

$$\delta \mathbf{x}(t) = \sum_{i=0}^{2N+1} K_i e^{-t/\tau_i} \mathbf{v}_i^r, \quad (\text{G.4})$$

where $\mathbf{v}_i^r$ denotes a right eigenvector of M with eigenvalue $-1/\tau_i$, and the $\tau_i$'s denote the relaxation times of the system. We need to distinguish between the left and the right eigenvectors, since M is not symmetric and hence the right eigenvectors are not necessarily orthogonal to their transpose. Note that the system is comprised of both the cells and the extracellular spaces, gathered in $\delta \mathbf{x}$, and thus has one set of $2N+1$ relaxation times and corresponding eigenvectors. We obtain the eigenvectors and eigenvalues by numerically diagonalizing M. The coefficients $K_i$ are determined by the initial conditions, i.e. the initial concentration of morphogen ($a(0)$) and the initial concentration of the intracellular signaling molecule
<table><caption>**Table H1.** Parameters used to generate plots. Parameter set A is used in figure 2. Parameter set B is used in figures 3, D1, and D2. Parameter set C is used in figures 4–7, D3, and D4.</caption><thead><tr><th>Parameters</th><th>Parameter set A</th><th>Parameter set B</th><th>Parameter set C</th></tr></thead><tbody><tr><td>α (# molecules V<sup>-1</sup> s<sup>-1</sup>)</td><td>0.998 33 × α<sup>crit</sup></td><td>Varying</td><td>Varying</td></tr><tr><td>c<sub>A</sub> (nM)</td><td>1.66 × 10<sup>2</sup></td><td>1.66 × 10<sup>2</sup></td><td>1.66 × 10<sup>2</sup></td></tr><tr><td>c<sub>B</sub> (nM)</td><td>1.66 × 10<sup>8</sup></td><td>1.66 × 10<sup>8</sup></td><td>1.66 × 10<sup>3</sup></td></tr><tr><td>D (μm<sup>2</sup> s<sup>-1</sup>)</td><td>1.0</td><td>1.0</td><td>1.0</td></tr><tr><td>δ (μm)</td><td>10.0</td><td>10.0</td><td>10.0</td></tr><tr><td>k<sub>A</sub> (s<sup>-1</sup>)</td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>k<sub>B</sub> (s<sup>-1</sup>)</td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td><td>100</td></tr><tr><td>N</td><td>600</td><td>600</td><td>600</td></tr><tr><td>s<sub>A</sub> (# molecules V<sup>-1</sup> s<sup>-1</sup>)</td><td>20</td><td>20</td><td>1000</td></tr><tr><td>s<sub>B</sub> (# molecules W<sup>-1</sup> s<sup>-1</sup>)</td><td>4</td><td>4</td><td>4</td></tr><tr><td>p</td><td>0.5</td><td>0.5</td><td>See legend</td></tr><tr><td>V (μm<sup>3</sup>)</td><td>10<sup>2</sup></td><td>10<sup>2</sup></td><td>10<sup>2</sup></td></tr><tr><td>W (μm<sup>3</sup>)</td><td>10<sup>3</sup></td><td>10<sup>3</sup></td><td>10<sup>3</sup></td></tr><tr><td>w</td><td>60</td><td>60</td><td>60</td></tr></tbody></table>

(b(0)). Gathering these initial conditions, as well as the steady-state solution in respective vectors similar to what we did for δx:

$$ \mathbf{x}(0) = \begin{bmatrix} a_0(0) \\ b_0(0) \\ \vdots \\ a_{N-1}(0) \\ b_{N-1}(0) \\ a_N(0) \end{bmatrix}, \quad \mathbf{x}^* = \begin{bmatrix} a_0^* \\ b_0^* \\ \vdots \\ a_{N-1}^* \\ b_{N-1}^* \\ a_N^* \end{bmatrix}, \quad (G.5) $$

we obtain the $K_i$ according to:

$$ K_i = \frac{1}{r} [\mathbf{v}_i^1 \cdot (\mathbf{x}(0) - \mathbf{x}^*)]. $$

That is, we compute the $K_i$ by projecting the initial condition onto the corresponding elements of the left eigenvector $\mathbf{v}_i^1$ and normalizing to $\mathbf{v}_i^1 \cdot \mathbf{v}_i^r = r$. Note that we use the *left* eigenvectors to compute the $K_i$ as we used the *right* ones to obtain the eigenmodes.

Based on this approximation to the dynamic solution of the system close to the steady state, we can define the characteristic time scale of the relay mechanism. For a linear system like the one presented in equation (G.3), the dynamics of the system are governed by the slowest relaxation time $\tau_{max}$ for long times. In particular, this slowest relaxation time is a good estimate of how long it takes the linear system to reach its steady state. We thus define the slowest relaxation time

$$ \tau_{max} := \max\{\tau_i\} \quad (G.6) $$

as the time scale of the relay mechanism.

## Appendix H. Parameter sets

Table H1 summarizes the three parameter sets used to generate the plots in this paper.

## ORCID iDs

Johanna E M Dickmann https://orcid.org/0000-0002-0861-4440

Jochen C Rink https://orcid.org/0000-0001-6381-6742

Frank Jülicher https://orcid.org/0000-0003-4731-9185

## References

[1] Gilbert S F 2014 *Developmental Biology* 10th edn (Sunderland, Mass: Sinauer Associates, Inc.)

[2] Morgan T H 1905 'Polarity' considered as a phenomenon of gradation of materials *J. Exp. Zool.* **2** 495–506

[3] Turing A M 1952 The chemical basis of morphogenesis *Phil. Trans. R. Soc.* B **237** 37–72

[4] Wolpert L 1969 Positional information and the spatial pattern of cellular differentiation *J. Theor. Biol.* **25** 1–47

[5] Crick F 1970 *Diffusion in embryogenesis* *Nature* **225** 671

[6] Gierer A, Berking S, Bode H, David C N, Flick K, Hansmann G, Schaller H and Trenkner E 1972 Regeneration of Hydra from reaggregated cells *Nat. New Biol.* **239** 98–101

[7] Driever W and Nüsslein-Volhard C 1988 A gradient of bicoid protein in Drosophila embryos *Cell* **54** 83–93

[8] Driever W and Nüsslein-Volhard C 1988 The bicoid protein determines position in the Drosophila embryo in a concentration-dependent manner *Cell* **54** 95–104

[9] Green J B A and Smith J C 1990 Graded changes in dose of a Xenopus activin a homologue elicit stepwise transitions in embryonic cell fate *Nature* **347** 391–4

[10] Gurdon J B, Harger P, Mitchell A and Lemaire P 1994 Activin signalling and response to a morphogen gradient *Nature* **371** 487–92

[11] Entchev E V, Schwabedissen A and González-Gaitán M 2000 Gradient formation of the TGF-β homolog Dpp *Cell* **103** 981–92

[12] Gregor T, Bialek W, Van Steveninck R R d R, Tank D W and Wieschaus E F 2005 Diffusion and scaling during early embryonic pattern formation *Proc. Natl. Acad. Sci. USA* **102** 18403–7

[13] Gregor T, Wieschaus E F, McGregor A P, Bialek W and Tank D W 2007 Stability and nuclear dynamics of the bicoid morphogen gradient *Cell* **130** 141–52

[14] Kicheva A, Pantazis P, Bollenbach T, Kalaidzidis Y, Bittig T, Jülicher F and González-Gaitán M 2007 Kinetics of morphogen gradient formation *Science* **315** 521–5

[15] Yu S R, Burkhardt M, Nowak M, Ries J, Petrásek Z, Scholpp S, Schwille P and Brand M 2009 Fgf8 morphogen gradient
forms by a source-sink mechanism with freely diffusing
molecules *Nature* **461** 533–6

[16] Müller P, Rogers K W, Jordan B M, Lee J S, Robson D, Ramanathan S and Schier A F 2012 Differential diffusivity of nodal and lefty underlies a reaction-diffusion patterning system *Science* **336** 721–4

[17] Sagner A and Briscoe J 2019 Establishing neuronal diversity in the spinal cord: a time and a place *Development* **146** dev182154

[18] Delgado I et al 2020 Proximo-distal positional information encoded by an Fgf-regulated gradient of homeodomain transcription factors in the vertebrate limb *Sci. Adv.* **6** 1–10

[19] Wartlick O, Mumcu P, Kicheva A, Bittig T, Seum C, Jülicher F and González-Gaitán M 2011 Dynamics of Dpp signaling and proliferation control *Science* **331** 1154–9

[20] Stapornwongkul K S, de Gennes M, Cocconi L, Salbreux G and Vincent J-P 2020 Patterning and growth control *in vivo* by an engineered GFP gradient *Science* **370** 321–7

[21] Romanova-Michaelides M, Hadjivasiliou Z, Aguilar-Hidalgo D, Basagiannis D, Seum C, Dubois M, Jülicher F and Gonzalez-Gaitan M 2022 Morphogen gradient scaling by recycling of intracellular Dpp *Nature* **602** 287–93

[22] Halpern K B et al 2017 Single-cell spatial reconstruction reveals global division of labour in the mammalian liver *Nature* **542** 352–6

[23] Halpern K B et al 2018 Paired-cell sequencing enables spatial gene expression mapping of liver endothelial cells *Nat. Biotechnol.* **36** 962

[24] Nacu E, Gromberg E, Oliveira C R, Drechsel D and Tanaka E M 2016 Fgf8 and SHH substitute for anterior–posterior tissue interactions to induce limb regeneration *Nature* **533** 407–10

[25] Gurley K A, Rink J C and Alvarado A S 2008 *β*-catenin defines head versus tail identity during planarian regeneration and homeostasis *Science* **319** 323–7

[26] Iglesias M, Gomez-Skarmeta J L, Saló E and Adell T 2008 Silencing of Smed-β catenin1 generates radial-like hypercephalized planarians *Development* **135** 1215–21

[27] Petersen C P and Reddien P W 2008 Smed-β catenin-1 is required for anteroposterior blastema polarity in planarian regeneration *Science* **319** 327–30

[28] Sureda-Gómez M, Pascual-Carreras E and Adell T 2015 Posterior Wnts have distinct roles in specification and patterning of the planarian posterior region *Int. J. Mol. Sci.* **16** 26543–54

[29] Reuter H et al 2015 β-catenin-dependent control of positional information along the AP body axis in planarians involves a teashirt family member *Cell Rep.* **10** 253–65

[30] Sureda-Gómez M, Martin-Durán J M and Adell T 2016 Localization of planarian *β*-catenin-1 reveals multiple roles during anterior–posterior regeneration and organogenesis *Development* **143** 4149–60

[31] Stückemann T, Cleland J P, Werner S, Vu H T-K, Bayersdorf R, Liu S-Y, Friedrich B, Jülicher F and Rink J C 2017 Antagonistic self-organizing patterning systems control maintenance and regeneration of the anteroposterior axis in planarians *Dev. Cell* **40** 248–63

[32] Porter J A, Young K E and Beachy P A 1996 Cholesterol modification of hedgehog signaling proteins in animal development *Science* **274** 255–9

[33] Pepinsky R B et al 1998 Identification of a palmitic acid-modified form of human sonic hedgehog *J. Biol. Chem.* **273** 14037–45

[34] Willert K, Brown J D, Danenberg E, Duncan A W, Weissman I L, Reya T, Yates J R and Nusse R 2003 Wnt proteins are lipid-modified and can act as stem cell growth factors *Nature* **423** 448–52

[35] Greco V, Hannus M and Eaton S 2001 Argosomes: a potential vehicle for the spread of morphogens through epithelia *Cell* **106** 633–45

[36] Panáková D, Sprong H, Marois E, Thiele C and Eaton S 2005 Lipoprotein particles are required for hedgehog and wingless signalling *Nature* **435** 58–65

[37] Bartscherer K and Boutros M 2008 Regulation of Wnt protein secretion and its role in gradient formation *EMBO Rep.* **9** 977–82

[38] Alexandre C, Baena-Lopez A and Vincent J-P 2014 Patterning and growth control by membrane-tethered wingless *Nature* **505** 180–5

[39] Mulligan K A, Fuerer C, Ching W, Fish M, Willert K and Nusse R 2012 Secreted wingless-interacting molecule (swim) promotes long-range signaling by maintaining wingless solubility *Proc. Natl Acad. Sci. USA* **109** 370–7

[40] Mii Y and Taira M 2009 Secreted frizzled-related proteins enhance the diffusion of Wnt ligands and expand their signalling range *Development* **136** 4083–8

[41] Stanganello E, Hagemann A I H, Mattes B, Sinner C, Meyen D, Weber S, Schug A, Raz E and Scholpp S 2015 Filopodia-based Wnt transport during vertebrate tissue patterning *Nat. Commun.* **6** 1–14

[42] Pani A M and Goldstein B 2018 Direct visualization of a native Wnt *in vivo* reveals that a long-range Wnt gradient forms by extracellular dispersal *eLife* **7** 1–22

[43] Howard J, Grill S W and Bois J S 2011 Turing’s next steps: the mechanochemical basis of morphogenesis *Nat. Rev. Mol. Cell Biol.* **12** 400–6

[44] Ramírez-Weber F-A and Kornberg T B 1999 Cytonemes: cellular processes that project to the principal signaling center in Drosophila imaginal discs *Cell* **97** 599–607

[45] Kerszberg M and Wolpert L 1998 Mechanisms for positional signalling by morphogen transport: a theoretical study *J. Theor. Biol.* **191** 103–14

[46] Gurley K A, Elliott S A, Simakov O, Schmidt H A, Holstein T W and Alvarado A S 2010 Expression of secreted Wnt pathway components reveals unexpected complexity of the planarian amputation response *Dev. Biol.* **347** 24–39

[47] Nusse R and Clevers H 2017 Wnt/β-catenin signaling, disease, and emerging therapeutic modalities *Cell* **169** 985–99

[48] Dieterle P B, Min J, Irimia D and Amir A 2020 Dynamics of diffusive cell signaling relays *eLife* **9** 1–48

[49] Yang Y and Mlodzik M 2015 Wnt-frizzled/planar cell polarity signaling: cellular orientation by facing the wind (Wnt) *Annu. Rev. Cell Dev. Biol.* **31** 623–46

[50] Barkai N and Ben-Zvi D 2009 'Big frog, small frog'—maintaining proportions in embryonic development *FEBS J.* **276** 1196–207

[51] Ben-Zvi D and Barkai N 2010 Scaling of morphogen gradients by an expansion-repression integral feedback control *Proc. Natl Acad. Sci. USA* **107** 6924–9

[52] Dickmann J 2020 Formation of long-ranged morphogen gradients by cell-to-cell relay *PhD Thesis* TU Dresden, Dresden, Germany

[53] Eldar A, Rosin D, Shilo B-Z and Barkai N 2003 Self-enhanced ligand degradation underlies robustness of morphogen gradients *Dev. Cell* **5** 635–46

[54] England J L and Cardy J 2005 Morphogen gradient from a noisy source *Phys. Rev. Lett.* **94** 07810–1–4

[55] Gregor T, Tank D W, Wieschaus E F and Bialek W 2007 Probing the limits to positional information *Cell* **130** 153–64

[56] Otsuki L and Tanaka E M 2022 Positional memory in vertebrate regeneration: a century's insights from the salamander limb *Cold Spring Harb. Perspect. Biol.* **14** a040899

[57] Morgan T H 1898 Experimental studies of the regeneration of planaria maculata *Arch. Entwickelungsmechanik Org.* **7** 364–97

[58] Petersen C P and Reddien P W 2009 A wound-induced Wnt expression program controls planarian regeneration polarity *Proc. Natl Acad. Sci. USA* **106** 17061–6
[59] Witchley J N, Mayer M, Wagner D E, Owen J H and Reddien P W 2013 Muscle cells provide instructions for planarian regeneration *Cell Rep.* **4** 633–41

[60] Lucila Scimone M, Cote L E, Rogers T and Reddien P W 2016 Two FGFRL-Wnt circuits organize the planarian anteroposterior axis *eLife* **5** 1–19

[61] Scimone M L, Cote L E and Reddien P W 2017 Orthogonal muscle fibres have different instructive roles in planarian regeneration *Nature* **551** 623–8

[62] Lengfeld T et al 2009 Multiple Wnts are involved in hydra organizer formation and regeneration *Dev. Biol.* **330** 186–99

[63] Hobmayer B, Rentzsch F, Kuhn K, Happel C M, Von Laue C, Snyder P, Rothbächer U and Holstein T W 2000 Wnt signalling molecules act in axis formation in the diploblastic metazoan hydra *Nature* **407** 186–9

[64] Kusserow A et al 2005 Unexpected complexity of the Wnt gene family in a sea anemone *Nature* **433** 156–60

[65] Reilly K M and Melton D A 1996 Short-range signaling by candidate morphogens of the TGF-β family and evidence for a relay mechanism of induction *Cell* **86** 743–54

[66] Lessing D and Nusse R 1998 Expression of wingless in the Drosophila embryo: a conserved cis-acting element lacking conserved Ci-binding sites is required for patched-mediated repression *Development* **125** 1469–76

[67] Yoffe K B, Manoukian A S, Wilder E L, Brand A H and Perrimon N 1995 Evidence for engrailed-independent wingless autoregulation in Drosophila *Dev. Biol.* **170** 636–50

[68] Ben-Zvi D, Pyrowolakis G, Barkai N and Shilo B-Z 2011 Expansion-repression mechanism for scaling the Dpp activation gradient in *Drosophila* wing imaginal discs *Curr. Biol.* **21** 1391–6

[69] Ascher U M, Ruuth S J and Spiteri R J 1997 Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations *Appl. Numer. Math.* **25** 151–67