The proneural wave in the Drosophila optic lobe is driven by an excitable reactiondiffusion mechanism
Abstract
In living organisms, selforganised waves of signalling activity propagate spatiotemporal information within tissues. During the development of the largest component of the visual processing centre of the Drosophila brain, a travelling wave of proneural gene expression initiates neurogenesis in the larval optic lobe primordium and drives the sequential transition of neuroepithelial cells into neuroblasts. Here, we propose that this ‘proneural wave’ is driven by an excitable reactiondiffusion system involving epidermal growth factor receptor (EGFR) signalling interacting with the proneural gene l’sc. Within this framework, a propagating transition zone emerges from molecular feedback and diffusion. Ectopic activation of EGFR signalling in clones within the neuroepithelium demonstrates that a transition wave can be excited anywhere in the tissue by inducing signalling activity, consistent with a key prediction of the model. Our model illuminates the physical and molecular underpinnings of proneural wave progression and suggests a generic mechanism for regulating the sequential differentiation of tissues.
https://doi.org/10.7554/eLife.40919.001Introduction
The development of multicellular organisms relies on a multitude of transient coordination processes that provide the spatiotemporal cues for cell fate decisionmaking and thereby ensure that tissues are specified with the correct size, pattern and composition (Perrimon et al., 2012; Oates et al., 2012; Sato et al., 2013). In one strategy, largescale patterning is engineered by selforganised concentration waves of biomolecular fate determinants that travel across tissues through intercellular exchange and the regulation of gene expression. Such travelling waves, which are viable carriers of spatiotemporal information, are a ubiquitous feature of developmental pattern formation, where they arise through different underlying mechanisms, from coordinated intracellular oscillations (Oates et al., 2012; Jörg et al., 2015; Hubaud et al., 2017; Verd et al., 2018) to selforganised reactiondiffusion processes (Lubensky et al., 2011; FormosaJordan et al., 2012; Fried et al., 2016; Gavish et al., 2016; Corson et al., 2017).
During the development of the fruitfly Drosophila melanogaster, a propagating wave of gene expression orchestrates the patterning of the largest component of the visual processing centre: neuroepithelial cells in the optic lobe of the larval brain divide symmetrically, expanding the progenitor pool, and then undergo a sequential transition into asymmetrically dividing neuroblasts, which generate the neurons of the medulla (Figure 1a,b) (Egger et al., 2007; Yasugi et al., 2008; Egger et al., 2010; Yasugi et al., 2010). The transition from neuroepithelial cells to neuroblasts occurs at a ‘transition zone’ that sweeps from one side of the optic lobe to the other and is marked by the expression of the proneural gene lethal of scute (l’sc) (Yasugi et al., 2008). This ‘proneural wave’ is controlled by the coordinated action of different signalling pathways: epidermal growth factor receptor (EGFR) signalling and DeltaNotch signalling (Figure 1c–f) (Yasugi et al., 2010). The sequential nature of the transition is crucial to generate populations of cells of different developmental ages that give rise to a diverse array of terminally differentiated medulla neurons (Li et al., 2013; Sato et al., 2013; Suzuki et al., 2013; Erclik et al., 2017). The transition zone exhibits localised EGFR signalling as well as expression of l’sc (Figure 1b,c). Absence of EGFR signalling leads to loss of the differentiation wave, indicating that EGFR signalling is a key component for proneural wave progression (Yasugi et al., 2010). The neuroepithelium exhibits low levels of Notch signalling activity (Egger et al., 2011; Weng et al., 2012). However, Notch activity peaks directly before the transition from neuroepithelial cell to neuroblast, drops during the transition and then is restored upon neuroblast transformation (Figure 1c) (Contreras et al., 2018). In addition, coordinating roles are played by the JAK/STAT and FatHippo pathways, which are broadly expressed in the neuroepithelium and prevent premature and ectopic transition of the neuroepithelium (Yasugi et al., 2008; Yasugi et al., 2010; Wang et al., 2011a; Reddy et al., 2010; Kawamori et al., 2011; Weng et al., 2012; Tanaka et al., 2018).
The question of how the specific functional feedbacks of EGFR signalling and proneural gene expression generate a localised propagating transition zone requires a mechanistic explanation of wave progression based on molecular feedbacks and signalling cascades. Such a description should explain (i) the dynamic nature of the wave, (ii) the emergence of a localised transition zone with spatially confined expression of the proneural gene l’sc and (iii) the specific profiles of gene expression and signalling activity around the transition zone. Moreover, the nature and function of the interaction of these components with DeltaNotch signalling, more commonly associated with lateral inhibition of neighbouring cells, is poorly understood, see Appendix 3. While a recent effort of a phenomenological description of the proneural wave (Sato et al., 2016) has started to model the coarsegrained aspects of proneural wave progression, the emergence of some major characteristics of the wave (such as spatially confined proneural gene expression in a localised transition zone) has not been addressed. Here we propose a model of signalling activity and proneural gene expression that describes the emergence of the proneural wave. Within this framework, the neuroepithelium behaves as an excitable medium in which changes in gene expression at the tissue boundary initiate a spontaneous wave of signalling activity that effects the transition from neuroepithelium to neuroblasts.
Results
Travelling front model of EGFR signalling activity
To develop the model, we first considered interactions between L’sc expression and associated signalling pathways within the transition zone. Previously, it was proposed that sequential induction of EGFR signalling is responsible for the progression of the proneural wave (Yasugi et al., 2008). EGFR signalling activates the expression of L’sc (Yasugi et al., 2010), which is sufficient to drive the neuroepithelium to neuroblast (NE to NB) transition (Yasugi et al., 2008; Contreras et al., 2018). The EGFR is activated by binding the secreted form of its ligand, Spitz. Secreted Spitz is generated by cleavage of a membranebound precursor by the transmembrane protease, Rhomboid (Klämbt, 2000). EGFR, Rhomboid and secreted Spitz together form an autocrine positive feedback loop (Figure 1d) (Wiley et al., 2003; Sato et al., 2013). In a first step, we noted that the dynamics of EGFR signalling alone has features that are sufficient to enable such a sequential induction and produce a travelling front of EGFR signalling activity; a feature notably absent in recent attempts to model the proneural wave, which also require further components to stabilise the propagating EGFR signalling front (Sato et al., 2016). In a minimal model based on the EGFR/Rhomboid/Spitz positive feedback loop, EGFR signalling activity is represented by a single component ‘E’ (e.g., the local cellular concentration of the active form of Spitz) that is diffusible between cells and involves the aforementioned positive feedback (Figure 1d and Figure 2a; Appendix 1),
Here ${\varphi}^{\mathrm{E}}$ denotes the local strength of E activity, $\eta $ is the effective diffusion constant, $\mu $ is the gain rate in signalling activity due to positive feedback, $k$ is the decay rate, and $h(\varphi )={\varphi}^{n}/(1+{\varphi}^{n})$ is a Hill function parameterising the nonlinear positive feedback. Generically, reactiondiffusion systems involving diffusion and selfactivation are known to support travelling bistable fronts that leave behind an elevated signalling state (Figure 2a; Video 1; Appendix 1) (Muratov and Shvartsman, 2004; Keener and Sneyd, 2009; Graham et al., 2010; Tayar et al., 2015). EGFR signalling is a natural candidate for being a key driver of the proneural wave and experimental evidence has shown that EGFR signalling is both necessary and sufficient for wave progression (Yasugi et al., 2010).
Travelling pulse model of EGFR signalling and proneural gene expression
However, notably, the EGFR/Rhomboid/Spitz feedback loop does not remain active in the wake of the travelling wave, but remains spatially confined as Rhomboid is expressed only transiently in the travelling transition zone (Yasugi et al., 2010; Sato et al., 2013). Therefore, in a second step, we considered the influence of the proneural gene l’sc, represented by a second component ‘L’ in our model. Elevated EGFR signalling activates L’sc expression (Yasugi et al., 2010), which is sufficient to drive the NE to NB transition (Yasugi et al., 2008; Contreras et al., 2018). In this minimal model, L’sc downregulates EGFR signalling as a consequence of the transition, leading to an indirect negative feedback (Figure 1f). The corresponding reactiondiffusion system for the local strength of E and L activity, ${\varphi}^{\mathrm{E}}$ and ${\varphi}^{\mathrm{L}}$, is given by
where ${\mu}_{i}$ (with $i=\mathrm{E},\mathrm{L}$) indicate production rates and ${k}_{i}$ denote decay rates. Simulations of Equation 2 demonstrate that this type of feedback is sufficient to describe a travelling localised pulse of signalling activity and L’sc expression through the tissue (Figure 2b,c; Video 2; Appendix 2). Notably, the dynamics of L’sc (L) alter the bistable signalling behaviour of EGFR signalling (E) into an excitable one: once sufficiently perturbed by diffusion from an adjacent cell with an elevated signalling state, the intracellular reaction dynamics produces a transient expression pulse that downregulates itself as a result of the NE to NB cell fate transition (Appendix 2).
Integrated model of the proneural wave to include EGFRL’scNotch interactions
We next aimed to develop a more refined model that could be challenged by experiment and compared with previously published data. Such a model necessarily includes DeltaNotch signalling, which has been shown to influence how long cells remain in the L’sc expressing state (Yasugi et al., 2010; Wang et al., 2011b; Weng et al., 2012). As a mediator of lateral inhibition, DeltaNotch signalling is often associated with the emergence of ‘saltandpepper’like patterns of cell fate (Bray, 2006; Shaya and Sprinzak, 2011). However, this pattern is not seen during proneural wave progression and the reasons for this are not clear (Egger et al., 2010; PérezGómez et al., 2013; Sato et al., 2016). To address this question, we extended our minimal model to include canonical DeltaNotch interactions (Figure 2d) (Collier et al., 1996; Bray, 2006; Simakov and Pismen, 2013): (i) transactivation of Notch by Delta, (ii) downregulation of Delta by Notch within the same cell and (iii) cisinhibition (downregulation of Notch by Delta in the same cell). The model incorporates interactions between the DeltaNotch signalling pathway, EGFR signalling and L’sc expression, namely, upregulation of Delta through EGFR signalling (Yasugi et al., 2010), upregulation of EGFR signalling through Notch signalling (Yasugi et al., 2010), downregulation of L’sc through Notch signalling (Reddy et al., 2010) and downregulation of Notch expression through L’sc (Egger et al., 2010). Despite these complex interactions, the functional ‘module’ comprising EGFR signalling and L’sc expression still remains the driver of the wave (green box in Figure 2d), while DeltaNotch signalling acts to provide further timing cues for the transition and to prevent premature differentiation (pink box in Figure 2d) (Egger et al., 2010; Reddy et al., 2010; Yasugi et al., 2010). The integrated model also includes an explicit representation of the cell state dynamics during the transition between neuroepithelial cell to neuroblast. The NE to NB transition is promoted by L’sc (Yasugi et al., 2008) and downregulation of Notch in the presence of EGFR signalling (Yasugi et al., 2010; Weng et al., 2012). The mathematical details of the refined model are given in Appendix 3.
Congruence with experimental data
In addition to the emergence of a propagating transition zone, the integrated model also yielded predictions on the spatial profiles of signalling activity and gene expression (Figure 2e; Figure 4a; Video 3; Video 4; Appendix 4), which were in striking agreement with features observed in prior experiments. First, EGFR signalling, as well as L’sc and Delta expression, was found to be elevated only in the transition zone (Figure 1c) (Egger et al., 2010; Yasugi et al., 2010). Second, a peak of Notch activity is observed slightly in advance of the transition zone (pink arrowheads in Figure 2e) followed by a sharp drop in Notch activity (Figure 1c) (Egger et al., 2011; OriharaOno et al., 2011; Weng et al., 2012; Contreras et al., 2018). According to the model, lateral inhibition promotes highDelta/lowNotch and lowDelta/highNotch states in adjacent cells, and leads to a travelling ‘laterally inhibited’ cell state as the wave progresses. By contrast, the drop in Notch levels in transitioning cells arises due to cisinhibition in our model as a consequence of Delta binding to Notch within the same cell, as has been shown experimentally (Reddy et al., 2010; Weng et al., 2012; Contreras et al., 2018).
To challenge the model further, we checked whether the documented effects of misregulating EGFR signalling, Notch signalling or L’sc expression (Yasugi et al., 2010) could be reproduced. For example, clones in which EGFR signalling has been constitutively activated tend to advance the transition zone within the clone, while the absence of EGFR signalling leads to loss of the proneural wave (Figure 3a,b) (Yasugi et al., 2010). To this end, we simulated the model dynamics on a twodimensional hexagonal lattice that mimics the topology of the neuroepithelium. Consistent with experiment, the model captured the acceleration, delay or loss of the proneural wavefront within a clone depending upon its genetic makeup (Figure 3a–f; Appendix 6).
Low levels of Notch signalling activity are observed both in the neuroepithelium and in neuroblasts but not at the transition zone (Figure 1c). It has been shown experimentally that Notch activity is required to maintain neuroepithelial cell fate and the loss of Notch results in premature transformation into neuroblasts (Egger et al., 2010; Ngo et al., 2010; Reddy et al., 2010; Yasugi et al., 2010; OriharaOno et al., 2011; Wang et al., 2011b; PérezGómez et al., 2013). Intriguingly, despite the observation of active Notch signalling, there is no evidence of lateral inhibition in the neuroepithelium. Lateral inhibition causes neighbouring cells to acquire complementary cell fates and so results in the emergence of a ‘saltandpepper pattern’ of Notch signalling activity (Bray, 2006; Shaya and Sprinzak, 2011). The reason for the absence of saltandpepper Notch signalling in the neuroepithelium is not clear. Notably, our model predicts that the basal level of Notch activity observed in the the neuroepithelium could be the reason for the suppression of lateral inhibition patterns outside the transition zone. In our model, basal levels of Notch activity in the neuroepithelium lead to a spatially homogeneous ‘oversaturation’ that prevents the Delta levels from rising before being activated by EGFR signalling. This is the case even in the presence of biochemical fluctuations (Figure 4a). However, if basal Notch levels are lowered to small values compared to the threshold levels for activation and inhibition in our model, we indeed recapitulate the saltandpepper patterns that are a consequence of lateral inhibition (Figure 4b). An analytical argument for the suppression of lateral inhibitions through basal Notch activity is given in Appendix 5.
We tested this prediction of the model by lowering Notch levels in the neuroepithelium but we did not observe ‘saltandpepper’ patterns of Delta/Notch expression within clones expressing an RNAi against Notch (Figure 4c). However, the absence of the emergence of lateral inhibition is likely due to the complete loss of detectable Notch in cells expressing the RNAi (Figure 4d), while the reduction of Notch levels in the model prediction is more subtle (Figure 4b). Referring to the ‘phase diagram’ in Figure 4e, it can be seen that both basal and Deltaregulated Notch activity need to be in the appropriate range for lateral inhibition patterns to occur, which is difficult to achieve experimentally. Furthermore, our model entails that Notch downregulation is a necessary (but not generally sufficient) condition for inducing saltandpepper patterns.
Dependence of wave speed and transition zone width on kinetic rate parameters
Next, we asked which aspects of the signalling and gene expression changes in our model have the largest effect on two important features of the system: the speed of the proneural wave and the width of the transition zone. To this end, we performed a sensitivity analysis on the kinetic rate parameters, as detailed in Appendix 7.
This analysis, based on the socalled ‘Morris method’ (Morris, 1991; Campolongo et al., 2007; Wu et al., 2013), entails a resampling of the parameter space of the model and yields three indices for each probed parameter. These indices indicate the impact of each parameter on the assessed output. The Morris indices $m$ and ${m}^{*}$ describe the impact of the respective parameter on the output, with $m$ including positive and negative effects (which may cancel each other as the parameter is varied) and ${m}^{*}$ the overall absolute effect (Campolongo et al., 2007). The third index $\sigma $ measures the nonlinearity of the parameter/output relation and/or interactions with other parameters (Wu et al., 2013). Detailed definitions of the respective indices are given in Appendix 7.
Here we probed the effects of the kinetic rate parameters of EGFR signalling, DeltaNotch signalling and L’sc expression on the propagation speed of the proneural wave and the width of the transition zone. As expected, this analysis showed that the diffusion, gain and decay rate of EGFR signalling ($\eta $, ${\mu}_{\mathrm{E}}$ and ${k}_{\mathrm{E}}$) are the key regulators of wave speed, since EGFR signalling constitutes the driver of the wave in our model (Figure 5a). In contrast, L’sc gene expression (parametrised by ${\mu}_{\mathrm{L}}$ and ${k}_{\mathrm{L}}$) had almost no effect on wave speed since its dynamics is ‘pulled’ by the EGFR signalling front. Interestingly, the basal gain rate of Notch signalling ($\beta $) as well as its decay rate (${k}_{\mathrm{N}}$) play another prominent role in setting the wave speed. This is consistent with experimental data showing that Notch signalling promotes EGFR signalling at the transition zone, leading to a reinforcement of activation as the wave arrives at undifferentiated cells (Yasugi et al., 2010). Considering the width of the transition zone (Figure 5b), we found that while all parameters had some effect, L’sc expression clearly had a tightening effect as it promoted differentiation and therefore leads to a faster termination of differentiation. In contrast, basal Notch signalling in the epithelium (described by $\beta $ and ${k}_{\mathrm{N}}$) tended to enlarge the width of the transition zone in our sensitivity analysis. Indeed, a recent study showed that overexpression of Notch at the transition zone extends the width of the L’sc stripe and delays the transformation into neuroblasts (Contreras et al., 2018), providing support for this prediction of the model. Subsequently, it follows that a complementary prediction of the model would be that a reduction of Notch at the transition zone would decrease the width of the L’sc stripe. We tested this prediction by knocking down Notch at the transition zone in clones (Figure 6). Within clones expressing Notch RNAi, the proneural wave was not only accelerated (Figure 3), as observed previously in Notch mutant clones (Egger et al., 2010; Reddy et al., 2010; Yasugi et al., 2010), but the width of the transition zone also appeared smaller (yellow arrowheads in Figure 6). In summary, this sensitivity analysis suggests that EGFR and Notch signalling are the key regulators of wave speed and width of the transition zone while L’sc expression provides an additional acceleration of the transition of the wave and therefore has a negative influence on the width of the transition zone.
Ectopic excitation of the transition in vivo
A signature feature of the model dynamics is that, through interactions between L’sc and EGFR signalling, the neuroepithelium functions as an excitable medium. As such, the model predicts that local induction of EGFR signalling would initiate a circular (targetlike) transition wave (Figure 7a; Video 5). To test whether a transition wave in the neuroepithelium could be excited at a position remote from the proneural wave, we induced clones expressing a downstream effector of the EGFR signalling pathway, Pointed P1 (PntP1; Figure 7b,c; Materials and methods). We found upregulation of L’sc at clonal boundaries and expression of Dpn within the clone, suggesting the ectopic generation of neuroblasts within the epithelium, that is, a NE to NB transition (Figure 7b,c). Our results agree with previous experiments showing the induction of neuroblasts within the neuroepithelium in response to ectopic EGFR signalling (Yasugi et al., 2010) and are in striking agreement with model simulations based on the same perturbation (Figure 7a; Video 5; Supplementary Text).
Discussion
Our findings suggest that the proneural wave involves the activation of an excitatory pulse of signalling activity and gene expression, giving rise to a tightlyregulated propagating transition zone. In contrast with Turingbased activatorinhibitor mechanisms (Turing, 1952), which typically comprise fast diffusible inhibitors, the reactiondiffusion system described here is based on strictly local inhibition (Figure 2b,c). The role of sequential patterning by the proneural wave is to ensure the correct timing and composition of the neuroblast population (Bertet et al., 2017). A similar process of sequential patterning occurs during the progression of the morphogenetic furrow in the Drosophila eye (Roignant and Treisman, 2009; Lubensky et al., 2011; FormosaJordan et al., 2012; Wartlick et al., 2014; Fried et al., 2016; Gavish et al., 2016). However, the progression of the morphogenetic furrow also entails transient growth as well as subsequent photoreceptor patterning and differentiation to generate ommatidia (Sato et al., 2013).
In our model, which focuses on the driving mechanism behind the proneural wave, we have refrained from considering additional signalling pathways that neither play key roles in driving the proneural wave nor exhibit strong signatures of bidirectional feedbacks (in contrast to EGFR and DeltaNotch signalling). These include the JAK/STAT and Hippo pathways, which serve important roles in modulating proneural wave progression but are not actively involved in propagating the transition zone through a reactiondiffusionlike mechanism (Yasugi et al., 2008; Reddy et al., 2010; Yasugi et al., 2010; Kawamori et al., 2011; Wang et al., 2011a; Weng et al., 2012).
On a mechanistic level, the excitable propagation behaviour illuminated here provides a mechanism to capture the transient and localised activity of the proneural gene l’sc and EGFR signalling, as well as robustness against fluctuating signalling activity and gene expression. In contrast to a recent model (Sato et al., 2016), our model also implies that neither differentiation nor proneural gene expression is required for the transient stabilisation of EGFR signalling activity that is required to advance the wave front. In the context of vertebrate somitogenesis, intracellular excitability has recently been suggested to underlie the emergence of genetic oscillations (Hubaud et al., 2017). The appearance here of excitability in the context of a propagating front of gene expression suggests that such a mechanism may serve more widely as a generic and robust strategy to achieve sequential transition waves in developing tissues.
Materials and methods
Fly strains
Request a detailed protocolFlies were raised on standard cornmeal medium at 25°C. Strains used were:
yw, hsFLP; FRT40A, tubGAL80/CyO, ActGFP; tubPGAL4,UASmCD8GFP/TM6B
w; FRT40A; UASPntP1/TM6B
w; FRT40A/CyO; N RNAi/TM6B (N RNAi lines used were BL33611 and BL33616)
HLHmgammaGFP (Almeida and Bray, 2005) was used to report active Notch signalling.
UASPntP1 clones and N RNAi clones
Request a detailed protocolTo generate clones in the developing optic lobe, larvae were collected 48–50 hr after egg laying and were heat shocked for 20 min at 37°C. Larvae were then dissected 50 or 60 hr after clone induction.
Immunohistochemistry
Request a detailed protocolLarval brains were dissected in PBS and fixed for 20 min at room temperature in 4% formaldehyde and fixation buffer (PBS, 5 mM MgCl2, 0.5 mM EGTA). After fixation, brains were rinsed and washed in 0.3% PBS Triton X100 (PBT). Samples were blocked with 10% normal goat serum (NGS) in 0.3% PBT at room temperature and incubated with the primary antibodies overnight at 4°C. Brains were then washed in 0.3% PBT and incubated with the secondary antibodies overnight at 4°C. Brains were washed in 0.3% PBT and mounted in Vectashield (Vector Laboratories, Burlingame, CA, USA). The following primary antibodies and dilutions were used: guinea pig antiDpn (1:10,000) and rat antiL’sc (1:5,000) (Caygill and Brand, 2017), chicken antiGFP (1:2,000) from Abcam, mouse antiDelta (1:50, C594.9B) from DSHB and mouse antiNotch (intracellular domain, ICD) (1:50, C17.9C6) from DSHB. Fluorescently conjugated secondary antibodies Alexa405, Alexa488, Alexa546 and Alexa633 (all 1:200) from Life Technologies.
Images were acquired with a Leica TCS SP8 confocal microscope (Leica Microsystems, Wetzlar, Germany) and analysed with Fiji (Schindelin et al., 2012). Figures and illustrations were assembled using Adobe Photoshop CS3 and Adobe Illustrator CS3 (Adobe Systems, San Jose, CA, USA).
Appendix 1
Bistable EGFR signalling fronts
We describe the details of our theoretical model by building up an integrated model of the proneural wave starting from a simple system of EGFR signalling activity. Reactiondiffusion models with different types of genetic interactions and correspondingly different phenomenologies have been used to describe the dynamics of EGFR signalling in other contexts such as Drosophila oogenesis (Shvartsman et al., 2002; Zartman et al., 2009; Simakov et al., 2012; Fauré et al., 2014) and the Drosophila retina (Pennington and Lubensky, 2010).
Model formulation
To illustrate the mechanism of wave propagation, we first show how EGFR signalling activity alone may give rise to a propagating signalling front. As outlined in the main text, we consider the dynamics of one signalling component E, which encapsulates the collective dynamics of the EGFR/Rhomboid/Spitz network (see below) and, hence, involves positive feedback and effective diffusion on the tissue level (see Figure 1c and Figure 2a). We describe the component E by a continuous field that defines the signalling activity, $\varphi =\varphi (\mathbf{\mathbf{x}},t)$, where $\mathbf{\mathbf{x}}$ denotes the position coordinate within the tissue and $t$ the time. (Throughout this text, we consider onedimensional systems for illustration purposes and twodimensional systems to describe the proneural wave in vivo. Hence, the position variable $\mathbf{\mathbf{x}}$ and the nabla operator $\nabla $ refer to $\mathbf{\mathbf{x}}=({x}_{1},{x}_{2})$ and $\nabla =(\partial /\partial {x}_{1},\partial /\partial {x}_{2})$ in the case of two dimensions and $\mathbf{\mathbf{x}}=x$ and $\nabla =\partial /\partial x$ in the case of one dimension.) The corresponding reactiondiffusion system is given by
where $\eta $ is the effective diffusion constant and $r(\varphi )$ is the reaction term describing the intracellular feedback,
Here, $\mu $ denotes the gain rate in signalling activity due to positive feedback, $k$ the degradation rate and $h$ is a monotonically increasing function describing the nonlinear positive feedback with $\mathrm{\Phi}$ being the threshold activity. Here, we choose a function of the Hill type (Novák and Tyson, 2008),
where $n$ is the ‘Hill exponent’ characterising the nonlinearity of the feedback. Systems of the type given by Equations 3–5 are known to exhibit solutions in the form of a propagating front for appropriate parameters (Keener and Sneyd, 2009). For illustrative purposes, we initially consider the dynamics of signalling activity in one spatial dimension on a domain of length $\mathrm{\ell}$ with ‘noflux’ boundary conditions, ${(\partial \varphi /\partial x)}_{x=0}=0={(\partial \varphi /\partial x)}_{x=\mathrm{\ell}}$. A numerical example of the system specified by Equations 3–5 is shown in Figure 2a. A localised concentration beyond a certain threshold concentration at the boundary $x=0$ initiates a travelling front of E, which propagates at constant velocity, leaving behind elevated levels of sustained signalling activity. (For the example shown in Figure 2a, we choose an initial condition of the form $\varphi (x,0)={\mathrm{e}}^{{x}^{2}/{x}_{0}^{2}}$ with ${x}_{0}^{2}=5\eta /k$.)
This behaviour is familiar in reactiondiffusion systems with reaction terms of the form given by Equation 4 and can be understood by studying the local reaction dynamics in a phase portrait: Appendix 1—figure 1a shows the reaction term $r$, which describes the local net loss/gain rate of signalling activity E. For sufficiently large gain rates, the function has three equilibrium points ${\varphi}_{i}^{*}$ for which the net loss/gain rate vanishes (dots in Appendix 1—figure 1a). Therefore, once the system has reached such an equilibrium, it remains there until perturbed. The two equilibrium points ${\varphi}_{0}^{*}$ and ${\varphi}_{2}^{*}$ are stable (black dots), that is small perturbations of the signalling activity will drive the system back to these equilibrium points—higher levels lead to a net loss whereas smaller levels lead to a net gain (black arrows in Appendix 1—figure 1a). The two equilibria are separated by an unstable equilibrium ${\varphi}_{1}^{*}$ (white dot). If the system is only slightly perturbed around $\varphi ={\varphi}_{0}^{*}$ (black dot), the reaction dynamics will drive the system back towards this equilibrium, because degradation is stronger than selfactivation in this region, $r<0$. However, if the signalling activity rises above ${\varphi}_{1}^{*}$ (e.g., due to diffusion from the neighboring cell), selfactivation will outcompete degradation, $r>0$, and the signalling activity will be driven towards the other stable equilibrium ${\varphi}_{2}^{*}$. This reactiondiffusion mechanism is known as a bistable front (Keener and Sneyd, 2009).
Further insight into the dynamics of the front propagation can be gained using analytical techniques (Keener and Sneyd, 2009). Nondimensionalising the model by rescaling $\varphi \to \varphi /\mathrm{\Phi}$, $t\to kt$, and $x\to \sqrt{k/\eta}x$, the corresponding reactiondiffusion system in one dimension is given by
Here, $\lambda =\mu /(k\mathrm{\Phi})$ is the only (dimensionless) parameter and $\rho $ is the dimensionless reaction term.
Anticipating that the system gives rise to a travelling front with velocity $c$ and a stationary front profile $\varphi $, we make the ansatz $\varphi (x,t)=\varphi (xct)$ in Equation 6, which yields the ordinary differential equation
where the prime denotes the derivative with respect to the argument $y=xct$ of $\varphi $.
Multiplying this equation by ${\varphi}^{\prime}$ and integrating from $\mathrm{\infty}$ to $\mathrm{\infty}$ yields an implicit expression for the velocity of the front,
where $\mathrm{\Delta}U={\mathrm{lim}}_{y\to \mathrm{\infty}}[U(\varphi (y))U(\varphi (y))]$ with the potential $U(\varphi )$ defined by $\rho =\partial U/\partial \varphi $ (see Appendix 1—figure 1b). Analytical solutions to Equation 8 only exist for special functional forms of the reaction term $\rho $. To obtain analytical insights into the wave speed, we consider the limiting case of large Hill exponents, $n\to \mathrm{\infty}$ in Equation 5. This corresponds to a switchlike response of the gain rate once the activation threshold is reached. In this case, the reaction term acquires a piecewise linear functional form,
where $\mathrm{\Theta}$ is the Heaviside step function with the convention $\mathrm{\Theta}(0)=1/2$ (see Appendix 1—figure 1c). This function has three equilibria for $\lambda >1$. For this reaction term, the front profile and front velocity can be calculated analytically via Equation 7 and Equation 8. The front profile is given by
where we have fixed the arbitrary position of the comoving reference frame by imposing the condition $\varphi (0)=1$. Appendix 1—figure 1e shows an example of the front solution Equation 10. The velocity of the front is given by
Expanding in the limit of large $\lambda $, dropping all orders higher than ${\lambda}^{1/2}$ and restoring the original parameter dependence by multiplying with the velocity scale $\sqrt{\eta k}$, we find
Equation 12 implies that the front velocity increases with increasing diffusion constant $\eta $ and synthesis rate $\mu $ as well as decreasing degradation rate $k$ and decreasing activation threshold $\mathrm{\Phi}$. Appendix 1—figure 1d shows that Equation 11 is a good approximation for the velocity of the front even when compared to numerical simulations with a finite Hill exponent $n$.
A onecomponent model reproduces the key features of more detailed models
Previously, we sought to capture the spatiotemporal dynamics of EGFR signalling by a single component, described by the signalling activity $\varphi $, even though EGFR signalling comprises three components: the EGF receptor, the transmembrane factor Rhomboid and the EGFR ligand Spitz (see Figure 1d). To show that the system above entails the essential features of more detailed descriptions, we now consider the kinetics of these three components and show that a corresponding model gives rise to the same qualitative dynamics. We consider a dimensionally reduced reactiondiffusion system where ${\psi}^{\mathrm{E}}$ denotes the concentration of bound EGFRs, ${\psi}^{\mathrm{R}}$ represents the concentration of Rhomboid, and ${\psi}^{\mathrm{S}}$ is the concentration of the secreted active form of Spitz (sSpi),
Here, ${\lambda}_{\mathrm{E}}$ is the binding rate of EGFRs, ${\lambda}_{\mathrm{R}}$ is the synthesis rate of Rhomboid, and ${\lambda}_{\mathrm{S}}$ is the secretion rate of Spitz; the functions ${H}_{\mathrm{E}}$ and ${H}_{\mathrm{R}}$ describe saturation of EGFR binding and saturation of the Rhomboid synthesis rate, respectively, and are assumed to be of the same qualitative form as the function $h$, given by Equation 5. Only the equation for the component S contains a diffusion term as only Spitz is secreted. The negative terms account for receptor unbinding and degradation of gene products; for simplicity, we have chosen identical unit decay rates and concentration thresholds. Numerical examples of Equation 13 show that the positive feedback of E, R, and S can generate a bistable front, see Appendix 1—figure 2.
Since each component is required to activate the next one in the loop, all three components show the same behaviour in terms of their localisation. Moreover, since the activation of one component is sufficient to activate the entire loop, diffusibility of S confers a diffusionlike effect to the entire E–R–S system. Therefore, the parameter $\mu $ of the onecomponent model Equation 3 and Equation 4 is to be interpreted as a composite parameter that characterises the strength of the positive feedback of the entire cycle. Likewise, the parameter $k$ represents the combined effects of degradation and receptor unbinding, and $\eta $ describes the effective diffusibility, mediated by diffusion of component S. Formally, a connection to the onecomponent model Equation 6 can be established when the dynamics of E and R is sufficiently fast such that they are always close to their stationary values determined by $\partial {\psi}^{\mathrm{E}}/\partial t=0$ and $\partial {\psi}^{\mathrm{R}}/\partial t=0$ while the front is travelling. In this case, we can solve for ${\psi}^{\mathrm{E}}$ and ${\psi}^{\mathrm{R}}$ and obtain a closed equation for $\varphi \equiv {\psi}^{\mathrm{S}}$,
where $\lambda \equiv {\lambda}_{\mathrm{S}}$ and $h(\varphi )={\lambda}_{\mathrm{R}}{H}_{\mathrm{R}}({\lambda}_{\mathrm{E}}{H}_{\mathrm{E}}(\varphi ))$. Note that Equation 14 is formally equivalent to Equation 6; if ${H}_{\mathrm{E}}$ and ${H}_{\mathrm{R}}$ are sigmoidal functions with the qualitive shape of Equation 5, then so is $h$.
Appendix 2
Excitable dynamics from EGFR–proneural interactions
While the minimal model Equations 3–5 leads to the emergence of a front that leaves behind an elevated signalling state, the transition zone of the proneural wave is characterised by localised EGFR signalling and proneural gene expression. We now show how such a travelling localised pulse arises if we take into account interactions between EGFR signalling and the proneural gene l’sc. The interactions between EGFR signalling and L’sc are schematically represented in Figure 1f and Figure 2b: the component E activates the expression of L while the component L effectively downregulates E as a consequence of the transition that is induced. The corresponding reactiondiffusion equations for the two fields ${\varphi}^{\mathrm{E}}$ and ${\varphi}^{\mathrm{L}}$ are given by
with the reaction terms
where we have introduced the Hill function $\overline{h}(\varphi )\equiv 1h(\varphi )$ which describes an inhibitory effect. For simplicity, we have also considered identical unit concentration thresholds for activation and inhibition. Here, ${\lambda}_{\mathrm{E}}$ and ${\lambda}_{\mathrm{L}}$ are the gain rates of the components E and L, respectively. Note that in the absence of L (${\varphi}^{\mathrm{L}}=0$), the reaction term ${\rho}_{\mathrm{E}}$ reduces to the one of the onecomponent model given in Equation 6, ${\rho}_{\mathrm{E}}({\varphi}^{\mathrm{E}},0)=\rho ({\varphi}^{\mathrm{E}})$ (Appendix 2—figure 1). Again, we consider a finite onedimensional domain of length $\mathrm{\ell}$ and noflux boundary conditions, ${(\partial {\varphi}^{\mathrm{E}}/\partial x)}_{x=0}=0={(\partial {\varphi}^{\mathrm{E}}/\partial x)}_{x=\mathrm{\ell}}$ and ${(\partial {\varphi}^{\mathrm{L}}/\partial x)}_{x=0}=0={(\partial {\varphi}^{\mathrm{L}}/\partial x)}_{x=\mathrm{\ell}}$.
Figure 2b in the main text displays a numerical example of the system Equation 15 and Equation 16. Again, starting from a localised concentration of the component E at $x=0$, a localised pulse of E and L travels through the system at a constant speed. (For the example shown in Figure 2b, we choose an initial condition of the form $\varphi (x,0)={\mathrm{e}}^{{x}^{2}/{x}_{0}^{2}}$ with ${x}_{0}^{2}=10\eta /k$.) Again, this behaviour can be elucidated by studying the local reaction dynamics, which is now given by the vector field $\mathbf{\mathbf{F}}=({\rho}_{\mathrm{E}},{\rho}_{\mathrm{L}})$, defined by Equation 16) (Appendix 2—figure 1). Note that now only the point $({\varphi}^{\mathrm{E}},{\varphi}^{\mathrm{L}})=(0,0)$ is a stable equilibrium. The black and white dots indicate points with ${\varphi}^{\mathrm{L}}=0$ that satisfy ${\rho}_{\mathrm{E}}=0$ as in the onecomponent system (compare to Appendix 1—figure 1a). However, note that only the point $(0,0)$ also satisfies ${\rho}_{\mathrm{L}}=0$. Again, diffusion from a neighboring cell will increase the levels of E. When the concentration level marked by the white point is exceeded, the reaction dynamics will elevate the levels of E and L until a turning point is reached when downregulation of E is sufficiently strong to suppress its positive selffeedback. Thus, the system finally returns to the fixed point $({\varphi}^{\mathrm{E}},{\varphi}^{\mathrm{L}})=(0,0)$ which marks the end of the pulse. Hence, the model Equation 15 and Equation 16 leads to propagation of a pulse of E and L at a defined speed.
Appendix 3
Integrated model of the proneural wave
The system described in Appendix 2 invokes a core mechanism giving rise to a propagating transition zone. As motivated in the main text, we now extend our model to include Delta–Notch interactions, based on the classical description advanced by (Collier et al., 1996) and also include cisinhibition (Sprinzak et al., 2010; Shaya and Sprinzak, 2011).
An earlier attempt to describe the proneural wave advanced by (Sato et al., 2016) has focused on a phenomenological description of the proneural wave, which, already in its general approach, differs from our model. The corresponding model is, in major parts, built around observed patterns of gene expression and patterning phenomena, rather than starting from the molecular underpinnings: there, proneural gene expression is not an independent dynamic ingredient; rather, the cell state (NE or NB) and the state of proneural gene expression is combined in a single abstract variable, making it impossible to address the effects of alterations in proneural gene expression independently of the NE to NB transition. Hence, in (Sato et al., 2016), the general driving mechanism of the wave is fundamentally different from our model: it relies on EGFR signalling inducing a change in cell state and, conversely, a change in cell state transiently driving EGFR signalling activity through a phenomenological prescription. In contrast to our model, this renders autocrine EGFR signalling without intrinsic bistability and therefore without the ability to autonomously drive the wavefront. Moreover, the resulting model is unstable with respect to additive fluctuations in gene expression and signalling activity; slight misexpression or perturbations in signalling activity will result in an immediate premature differentiation, as discussed in (Sato et al., 2016). DeltaNotch interactions are incorporated as a subsystem without intrinsic multistability usually found necessary in attempts to describe the emergence of lateral inhibition phenomena (Collier et al., 1996).
Since Delta–Notch interactions can give rise to lateral inhibition, that is stable lowDelta/highNotch and highDelta/lowNotch states in adjacent cells, instead of a continuous description of the tissue, we now consider a lattice where each lattice site represents a cell. We consider the signalling and gene activities ${\varphi}_{\mathbf{\mathbf{x}}}^{i}(t)$ with $i=\mathrm{E},\mathrm{L},\mathrm{D},\mathrm{N}$ describing EGFR signalling activity, L’sc expression, and Delta and Notch, respectively. The index $\mathbf{x}$ indicates the lattice site. Furthermore, we introduce the cell state ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}(t)$ which takes values from 0 to 1, where ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}=0$ indicates that cell $\mathbf{x}$ is a neuroepithelial cell and ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}=1$ indicates that cell $\mathbf{x}$ is a neuroblast. Figure 2d shows the regulatory network of the model. In contrast to the effective inhibition of signalling by the proneural gene l’sc that we considered before, we now include the more realistic shutdown of signalling as a consequence of differentiation. Moreover, motivated by the presence of low levels of Notch signalling in the neuroepithelium and the neuroblasts (Egger et al., 2010; OriharaOno et al., 2011), we include a basal source of Notch that is independent of transactivation by Delta in adjacent cells. As in the previous sections, we here consider identical reference concentrations for activation and inhibition and rescale all concentrations by this reference concentration, so that the fields ${\varphi}_{\mathbf{\mathbf{x}}}^{i}$ with $i=\mathrm{E},\mathrm{L},\mathrm{D},\mathrm{N}$ are dimensionless. Moreover, we consider identical degradation constants for all four components and rescale time by the degradation time. The corresponding dynamical equations are given by
Here, the parameters ${\mu}_{i}$ denote gain rates, ${k}_{i}$ denote decay rates, $\beta $ denotes the basal gain rate of the component N, and the ${\mathrm{\Phi}}_{i}$ denote threshold levels for positive and negative feedbacks. The operators $\widehat{\mathrm{\Delta}}$ and $\widehat{\mathrm{\Sigma}}$ are the lattice Laplacian and the sum over concentrations of neighboring lattice sites, respectively, and defined by ${[\widehat{\mathrm{\Delta}}\varphi ]}_{\mathbf{\mathbf{x}}}={\sum}_{y\in {U}_{\mathbf{\mathbf{x}}}}({\varphi}_{\mathbf{\mathbf{y}}}{\varphi}_{\mathbf{\mathbf{x}}})$ and ${[\widehat{\mathrm{\Sigma}}\varphi ]}_{\mathbf{\mathbf{x}}}={\sum}_{y\in {U}_{\mathbf{\mathbf{x}}}}{\varphi}_{\mathbf{\mathbf{y}}}$, where ${U}_{\mathbf{\mathbf{x}}}$ is the set of neighbours of site $x$. The dynamics of the cell state ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}$ is given by
where the function $V$ is a ‘potential’ for the cell state which ensures that in the absence of signalling and proneural gene expression, ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}$ has two stable equilibria: ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}=0$ (neuroepithelium) and ${\mathrm{\Omega}}_{\mathbf{\mathbf{x}}}=1$ (neuroblast) (see Appendix 3—figure 1). The qualitative features of our model do not depend on the exact choice of $V$. The term ${f}_{\mathrm{int}}$ acts as a ‘force’ that triggers the transition from one state to the other depending on the local signalling activity and proneural gene expression. The functional form of ${f}_{\mathrm{int}}$ is based on the observations that (i) proneural gene expression seems to be sufficient but not necessary for the transition while (ii) Notch expression seems to keep cells in their proneural state (Yasugi et al., 2010). Therefore, we here choose
The cell state potential $V$ has wells at $\mathrm{\Omega}=0$ (neuroepithelial state) and $\mathrm{\Omega}=1$ (neuroblast state) (see Appendix 3—figure 1), so that both states are stable and a transition occurs when the hill in between both states can be overcome. This choice of ${f}_{\mathrm{int}}$ leads to a transition of the cell state from neuroepithelial ($\mathrm{\Omega}=0$) to neuroblast ($\mathrm{\Omega}=1$) when (i) L exceeds the threshold level ${\mathrm{\Phi}}_{\mathrm{int}}$ and/or (ii) E exceeds and N drops below the threshold level.
Appendix 4
Robustness of the model against fluctuations and disorder
Robustness against biochemical noise
Gene expression and biochemical reactions typically suffer from fluctuations due to small numbers of molecules involved (Tsimring, 2014). To achieve reliable morphological results, any biochemical mechanism governing morphogenetic processes must be robust against such types of noise. From an analytical point of view, stability of the zero signalling fixed point in our model (see Appendix 2—figure 1) ensures that a cell does not differentiate prematurely due to a certain degree of noise in proneural gene expression or fluctuating signalling activity. To numerically demonstrate that our model is robust against fluctuations in molecule concentrations, we performed simulations of the system in the presence of biochemical fluctuations in all four components. The noisy system is given by Equations 17–19 with each dynamical equation replaced according to
where $\gamma $ denotes the noise strength and ${\xi}_{i}$ denotes Gaussian white noise characterised by the expectation values $\u27e8{\xi}_{i}(t)\u27e9=0$ and $\u27e8{\xi}_{i}(t){\xi}_{j}({t}^{\prime})\u27e9={\delta}_{ij}\delta (t{t}^{\prime})$. Furthermore, ${\delta}_{ij}$ denotes the Kronecker delta and $\delta (t)$ the Dirac delta distribution.
Appendix 4—figure 1 shows numerical examples of the system for different noise strengths $\gamma $. In these examples, the wave robustly travels from the left to the right for small and intermediate noise levels (as compared to the gain rate ${\mu}_{\mathrm{E}}$) (Appendix 4—figure 1a,b), whereas premature differentiation is only observed for large noise levels that introduce fluctuations comparable to physiological concentrations (Appendix 4—figure 1c). Finally, random differentiation throughout the tissue only occurs if the system is dominated by fluctuations (Appendix 4—figure 1d).
Robustness against lattice defects
To test whether the proposed mechanism is robust with respect to disordered lattice structures, we considered a sitediluted version of our model, in which randomly chosen sites in the hexagonal lattice were ‘deactivated’, that is removed from the dynamics and the coupling topology. Simulations showed that while the presence of disorder leads to a local distortion of the propagating front profile, the overall mechanism remains intact and leads to robust and sequential differentiation of the neuroepithelial tissue (Appendix 4—figure 2a). Even regions which are partially ‘shielded’ by defects eventually become differentiated as the diffusionmediated propagation has no intrinsic directionality and is able to reach such regions when the wave has surrounded the corresponding region (Appendix 4—figure 2). Moreover, the overall phenomenology of mutant and transgenic clones remains intact on such imperfect lattices, as illustrated using the ‘Notch upregulation’ clone (cf. Figure 3f and Appendix 4—figure 2b).
Appendix 5
Suppression of lateral inhibition patterns
An analytical argument demonstrating how basal Notch levels suppress lateral inhibition can be made in a simple picture involving Delta–Notch interactions in two cells $\mathbf{x}=1,2$,
where $\overline{\mathbf{x}}$ refers to the respective other cell and $\overline{h}(\varphi )\equiv 1h(\varphi )$ with $h$ being the Hill function Equation 5, as before. Here, $\lambda $ is the gain rate for Notch and $\beta $ indicates the basal production. For simplicity, we consider a linear positive feedback in the dynamics of Notch and the Hill exponent $n=2$ for $\overline{h}$. In steady state, where $\mathrm{d}{\varphi}_{\mathbf{\mathbf{x}}}^{\mathrm{D}}/\mathrm{d}t=0=\mathrm{d}{\varphi}_{\mathbf{\mathbf{x}}}^{\mathrm{N}}/\mathrm{d}t$, we can eliminate ${\varphi}_{1}^{\mathrm{N}}$ and ${\varphi}_{2}^{\mathrm{N}}$ to obtain
where $\mathrm{\Gamma}(\varphi )=\overline{h}(\beta +\lambda \varphi )$. From this, it follows that both ${\varphi}_{1}^{\mathrm{D}}$ and ${\varphi}_{2}^{\mathrm{D}}$ satisfy ${\varphi}_{\mathbf{\mathbf{x}}}^{\mathrm{D}}=\mathrm{\Gamma}(\mathrm{\Gamma}({\varphi}_{\mathbf{\mathbf{x}}}^{\mathrm{D}}))$. Among other solutions, this equation has two solutions of the form ${\varphi}_{\pm}=p\pm \sqrt{q}$ with
which, in the case that they are real, correspond to the lowDelta/highNotch and highDelta/lowNotch states in adjacent cells as they satisfy ${\varphi}_{\pm}=\mathrm{\Gamma}({\varphi}_{\mp})$.
However, bistability only exists if both solutions are real, which is the case for $q>0$. From Equation 23, we find that this corresponds to
Figure 4e shows the corresponding phase diagram for the occurrence of lateral inhibition in the twocell system. Therefore, a basal production term, if strong enough, can prevent lateral inhibition patterns.
Appendix 6
Description of mutant and transgenic clones
Simulation of clones
To capture experimental scenarios in which mutant or transgenic clones were induced in the neuroepithelium, we modified the model Equation 17 accordingly. In our model, a clone is defined by a cell or a group of cells within the simulated tissue which has altered kinetic rate parameters or a different initial condition, depending on the type of experimental perturbation (Figure 3). For the case of downregulation of proneural gene or signalling factors, the synthesis or binding rate of the respective gene or receptor in the clone cells is decreased or set to zero, as indicated in the caption of Figure 3. For the case of upregulation, which in all considered cases corresponds to a constitutively active gene, we added a source term to the clone that leads to constant synthesis and furthermore set the initial condition of the clone to the elevated steadystate concentration of the respective gene or signalling activity.
Simulation of Figure 7a
To simulate the effects of clones in which EGFR signalling is constitutively activated within the neuroepithelium (Figure 7), we simulated the model on a hexagonal lattice with circular boundaries (Video 4). The initial condition was set to ${{\varphi}^{\mathrm{E}}}_{t=0}={\mu}_{\mathrm{E}}/2$ in a onedimensional array of cells in the outermost cell layer that have angles between $\pi /3$ and $5\pi /3$ as measured from the center of the circular lattice. Moreover, we arbitrarily selected four lattice sites and endowed them with a constant production of E by adding the term ${\mu}_{\mathrm{E}}\overline{h}(2{\mathrm{\Omega}}_{\mathbf{\mathbf{x}}})$ to the reaction dynamics of the component E in Equation 17; this mimics the constitutively active EGFR signalling. Figure 7a shows a snapshot at time $t=14$. All parameter values are given in Appendix 3—table 1.
Appendix 7
Sensitivity analysis of the model
Morris method for global sensitivity analyses
To test the sensitivity of key observables on model parameters, we here employ the socalled Morris method, a widely used method for global sensitivity analyses (Morris, 1991; Campolongo et al., 2007; Wu et al., 2013). To briefly summarise, for a fixed model observable $\mathcal{O}$ and a given set of parameters ${\theta}_{1},\mathrm{\dots},{\theta}_{n}$, the Morris method consists in repeatedly sampling a discretised parameter space (or subspace of interest) of the model and, for each parameter $i$, calculating the socalled ‘elementary effects’
That is the finitedifference quotient of the output with respect to the parameter, given a finite step size $\mathrm{\Delta}$ that is chosen adequately; for standard choices and further details on the method, see, for example (Wu et al., 2013). This sampling procedure yields a distribution ${P}_{i}(e)$ of elementary effects for each parameter $i$, from which the following sensitivity indices are computed,
where ${\u27e8\cdot \u27e9}_{i}$ denotes the expectation value under the distribution ${P}_{i}$. The interpretation of these indices is given in the main text.
Probed observables and parameters
As output observables we here choose the linear propagation velocity $v$ of the proneural wave and the width $w$ of the transition zone. We formally define these quantities as follows for the latticebased full proneural wave model Equations 17–19. By $x$, we denote the extension of the system in the direction of the travelling wave. We define ${\overline{\mathrm{\Omega}}}_{x}={\mathrm{\ell}}_{\u27c2}^{1}{\sum}_{y}{\mathrm{\Omega}}_{xy}$ as the average of the cell state $\mathrm{\Omega}$ in the direction perpendicular to the wave, where ${\mathrm{\ell}}_{\u27c2}$ is the extension of the lattice in the perpendicular direction. A wave with constant velocity leads to a proportionally linear increase in number of neuroblasts, so that the wave velocity $v$ (in lattice sites per unit time) is given by
Practically, we determine $v$ as the slope obtained from a linear fit of ${\sum}_{x}{\overline{\mathrm{\Omega}}}_{x}$ in the linear regime.
We define the width $w$ of the transition zone via the spatial spread of transitioning cells, that is those with $0<\mathrm{\Omega}<1$. Formally, this width can be defined as
The discrete derivative, ${\overline{\mathrm{\Omega}}}_{x}{\overline{\mathrm{\Omega}}}_{x1}$, which measures the steepness of the profile, is nonzero only in the transition zone. For example, for a Gaussian profile of the discrete derivative, ${\overline{\mathrm{\Omega}}}_{x}{\overline{\mathrm{\Omega}}}_{x1}\propto {\mathrm{e}}^{x/2{\sigma}^{2}}$ with a variance ${\sigma}^{2}\gg 1$, Equation 28 yields $w=\sigma $. To avoid confounding effects by initial and boundary conditions in model simulations, we use the temporal median of $w$ as a proxy for the width of the transition zone.
We compute the Morris indices Equation 26 for the dependence of $v$ and $w$ on the kinetic and diffusion properties of the integrated model, that is on the parameters $\eta $, $\beta $ and ${\mu}_{i}$, ${k}_{i}$ for $i=\mathrm{E},\mathrm{L},\mathrm{D},\mathrm{N}$ and allow them to vary between the 0.2fold and 5fold reference value given in Appendix 3—table 1, while keeping all other parameters fixed to their values given in Appendix 3—table 1. 55000 parameter samples were used to compute expectation values.
Data availability
All data generated or analysed during this study are included in the manuscript.
References

Regulation of postembryonic neuroblasts by Drosophila grainyheadMechanisms of Development 122:1282–1293.https://doi.org/10.1016/j.mod.2005.08.004

BookThe Developmental Origin of Cell Type Diversity in the Drosophila Visual SystemCham, Switzerland: Springer International Publishing.

Notch signalling: a simple pathway becomes complexNature Reviews Molecular Cell Biology 7:678–689.https://doi.org/10.1038/nrm2009

An effective screening design for sensitivity analysis of large modelsEnvironmental Modelling & Software 22:1509–1518.https://doi.org/10.1016/j.envsoft.2006.10.004

Pattern formation by lateral inhibition with feedback: a mathematical model of deltanotch intercellular signallingJournal of Theoretical Biology 183:429–446.https://doi.org/10.1006/jtbi.1996.0233

A discrete model of Drosophila eggshell patterning reveals cellautonomous and juxtacrine effectsPLOS Computational Biology 10:e1003527.https://doi.org/10.1371/journal.pcbi.1003527

Regulation of neuronal differentiation at the neurogenic wavefrontDevelopment 139:2321–2329.https://doi.org/10.1242/dev.076406

A model of the Spatiotemporal dynamics of Drosophila eye disc developmentPLOS Computational Biology 12:e1005052.https://doi.org/10.1371/journal.pcbi.1005052

Continuum theory of gene expression waves during vertebrate segmentationNew Journal of Physics 17:093042.https://doi.org/10.1088/13672630/17/9/093042

Fat/Hippo pathway regulates the progress of neural differentiation signaling in the Drosophila optic lobeDevelopment, Growth & Differentiation 53:653–667.https://doi.org/10.1111/j.1440169X.2011.01279.x

BookMathematical Physiology (2nd ed)Springer Science+Business Media.https://doi.org/10.1007/9780387758473

EGF receptor signalling: the importance of presentationCurrent Biology 10:R388–R391.https://doi.org/10.1016/S09609822(00)004851

Signal propagation and failure in discrete autocrine relaysPhysical Review Letters 93:500–504.https://doi.org/10.1103/PhysRevLett.93.118101

Design principles of biochemical oscillatorsNature Reviews Molecular Cell Biology 9:981–991.https://doi.org/10.1038/nrm2530

Switch and template pattern formation in a discrete reactiondiffusion system inspired by the Drosophila eyeThe European Physical Journal E 33:129–148.https://doi.org/10.1140/epje/i2010106476

Signaling mechanisms controlling cell fate and embryonic patterningCold Spring Harbor Perspectives in Biology 4:a005975.https://doi.org/10.1101/cshperspect.a005975

Pattern formation in the Drosophila eye discThe International Journal of Developmental Biology 53:795–804.https://doi.org/10.1387/ijdb.072483jr

Waves of differentiation in the fly visual systemDevelopmental Biology 380:1–11.https://doi.org/10.1016/j.ydbio.2013.04.007

Fiji: an opensource platform for biologicalimage analysisNature Methods 9:676–682.https://doi.org/10.1038/nmeth.2019

From Notch signaling to finegrained patterning: Modeling meets experimentsCurrent Opinion in Genetics & Development 21:732–739.https://doi.org/10.1016/j.gde.2011.07.007

Modeling and computational analysis of EGF receptormediated cell communication in Drosophila oogenesisDevelopment 129:2577–2589.

Noise in biologyReports on Progress in Physics 77:026601–026630.https://doi.org/10.1088/00344885/77/2/026601

The chemical basis of morphogenesisPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 237:37–72.https://doi.org/10.1098/rstb.1952.0012

Role of JAK/STAT signaling in neuroepithelial stem cell maintenance and proliferation in the Drosophila optic lobeBiochemical and Biophysical Research Communications 410:714–720.https://doi.org/10.1016/j.bbrc.2011.05.119

Changes in Notch signaling coordinates maintenance and differentiation of the Drosophila larval optic lobe neuroepitheliaDevelopmental Neurobiology 72:1376–1390.https://doi.org/10.1002/dneu.20995

Sensitivity analysis of infectious disease models: methods, advances and their applicationJournal of The Royal Society Interface 10:20121018.https://doi.org/10.1098/rsif.2012.1018
Article and author information
Author details
Funding
Wellcome (092096)
 David J Jörg
 Elizabeth E Caygill
 Anna E Hakes
Cancer Research UK (C6946/A14492)
 David J Jörg
 Elizabeth E Caygill
 Anna E Hakes
 Esteban G Contreras
 Andrea H Brand
 Benjamin D Simons
Biotechnology and Biological Sciences Research Council (Project Grant BB/L007800/1)
 Elizabeth E Caygill
 Andrea H Brand
Wellcome (Senior Investigator Award 103792)
 Anna E Hakes
 Esteban G Contreras
 Andrea H Brand
Royal Society (Darwin Trust Research Professorship)
 Andrea H Brand
Wellcome (Senior Investigator Award 098357)
 Benjamin D Simons
Royal Society (EP Abraham Research Professorship RP\R1\180165)
 Benjamin D Simons
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We are grateful to Pau FormosaJordan, Claude Desplan, François Schweisguth and members of the Simons and Brand labs for useful discussions.
Version history
 Received: August 9, 2018
 Accepted: February 8, 2019
 Version of Record published: February 22, 2019 (version 1)
Copyright
© 2019, Jörg et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,999
 views

 300
 downloads

 14
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Developmental Biology
Ninein is a centrosome protein that has been implicated in microtubule anchorage and centrosome cohesion. Mutations in the human NINEIN gene have been linked to Seckel syndrome and to a rare form of skeletal dysplasia. However, the role of ninein in skeletal development remains unknown. Here, we describe a ninein knockout mouse with advanced endochondral ossification during embryonic development. Although the long bones maintain a regular size, the absence of ninein delays the formation of the bone marrow cavity in the prenatal tibia. Likewise, intramembranous ossification in the skull is more developed, leading to a premature closure of the interfrontal suture. We demonstrate that ninein is strongly expressed in osteoclasts of control mice, and that its absence reduces the fusion of precursor cells into syncytial osteoclasts, whereas the number of osteoblasts remains unaffected. As a consequence, nineindeficient osteoclasts have a reduced capacity to resorb bone. At the cellular level, the absence of ninein interferes with centrosomal microtubule organization, reduces centrosome cohesion, and provokes the loss of centrosome clustering in multinucleated mature osteoclasts. We propose that centrosomal ninein is important for osteoclast fusion, to enable a functional balance between boneforming osteoblasts and boneresorbing osteoclasts during skeletal development.

 Cell Biology
 Developmental Biology
The generation of distinct cell fates during development depends on asymmetric cell division of progenitor cells. In the central and peripheral nervous system of Drosophila, progenitor cells respectively called neuroblasts or sensory organ precursors use PAR polarity during mitosis to control cell fate determination in their daughter cells. How polarity and the cell cycle are coupled, and how the cell cycle machinery regulates PAR protein function and cell fate determination is poorly understood. Here, we generate an analog sensitive allele of CDK1 and reveal that its partial inhibition weakens but does not abolish apical polarity in embryonic and larval neuroblasts and leads to defects in polarisation of fate determinants. We describe a novel in vivo phosphorylation of Bazooka, the Drosophila homolog of PAR3, on Serine180, a consensus CDK phosphorylation site. In some tissular contexts, phosphorylation of Serine180 occurs in asymmetrically dividing cells but not in their symmetrically dividing neighbours. In neuroblasts, Serine180 phosphomutants disrupt the timing of basal polarisation. Serine180 phosphomutants also affect the specification and binary cell fate determination of sensory organ precursors as well as Baz localisation during their asymmetric cell divisions. Finally, we show that CDK1 phosphorylates SerineS180 and an equivalent Serine on human PAR3 in vitro.