The proneural wave in the Drosophila optic lobe is driven by an excitable reaction-diffusion mechanism

In living organisms, self-organised 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 reaction-diffusion 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.


Introduction
The development of multicellular organisms relies on a multitude of transient coordination processes that provide the spatiotemporal cues for cell fate decision-making 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, large-scale patterning is engineered by self-organised 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 self-organised reaction-diffusion processes (Lubensky et al., 2011;Formosa-Jordan 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 Delta-Notch 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 increases just before the transition zone (marked by L'sc, cyan), drops during the transition and then increases again in neuroblasts. Bottom: The proneural wave, characterised by expression of L'sc as well as EGF receptor (EGFR) and Notch signalling activity, sequentially converts the neuroepithelium into neuroblasts. (d) EGFR signalling in the transition zone activates expression of the transmembrane protein Rhomboid, which in turn cleaves the membrane-tethered form of the EGFR ligand Spitz (mSpi) to generate its active secreted form (sSpi). (The shaded region depicts an individual cell in the neuroepithelium.) sSpi can bind to the EGFR on the same cell and neighbouring cells. (e) Delta-Notch signalling is a contactdependent signalling pathway active in both the neuroepithelium and the neuroblasts. The Delta ligand binds to Notch receptors on adjacent cells upon which their intracellular domain (NICD) is cleaved. The NICD regulates target genes, which, in turn, affects expression of Delta. (f) Active EGFR signalling promotes the expression of L'sc within the same cell, which is sufficient for the neuroepithelium to neuroblast transition and which in turn downregulates EGFR signalling. DOI: https://doi.org/10.7554/eLife.40919.002 by the JAK/STAT and Fat-Hippo 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 Delta-Notch 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 coarse-grained 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.

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 membrane-bound 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 f E denotes the local strength of E activity, h is the effective diffusion constant, is the gain rate in signalling activity due to positive feedback, k is the decay rate, and hðfÞ ¼ f n =ð1 þ f n Þ is a Hill function parameterising the nonlinear positive feedback. Generically, reaction-diffusion systems involving diffusion and self-activation 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 reaction-diffusion system for the local strength of E and L activity, f E and f L , is given by (2) where i (with i ¼ E; 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 EGFR-L'sc-Notch 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 Delta-Notch 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, Delta-Notch signalling is often associated with the emergence of 'salt-and-pepper'-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érez-Gó mez et al., 2013;Sato et al., 2016). To address this question, we extended our minimal model to include canonical Delta-Notch interactions ( Figure 2d) (Collier et al., 1996;Bray, 2006;Simakov and Pismen, 2013): (i) trans-activation of Notch by Delta, (ii) downregulation of Delta by Notch within the same cell and (iii) cis-inhibition (downregulation of Notch by Delta in the same cell). The model incorporates interactions between the Delta-Notch 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 Delta-Notch 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;Orihara-Ono et al., 2011;Weng et al., 2012;Contreras et al., 2018). According to the model, lateral inhibition promotes high-Delta/low-Notch and low-Delta/high-Notch 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 cis-inhibition 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;Orihara-Ono et al., 2011;Wang et al., 2011b;Pérez-Gó 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 'salt-and-pepper pattern' of Notch signalling activity (Bray, 2006;Shaya and Sprinzak, 2011). The reason for the absence of salt-and-pepper 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 salt-and-pepper 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 'salt-and-pepper' 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 Delta-regulated 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 salt-and-pepper 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 so-called '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 s measures the non-linearity 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, Delta-Notch 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 (h, E and k 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 L and k 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 (b) as well as its decay rate (k 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 b and k 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 (target-like) 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 tightly-regulated propagating transition zone. In contrast with Turing-based activator-inhibitor mechanisms (Turing, 1952), which typically comprise fast diffusible inhibitors, the reaction-diffusion 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 . 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;Formosa-Jordan 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 . 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 Delta-Notch 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 reaction-diffusion-like 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).

L'sc Dpn
Delta L'sc a b L'sc Dpn L'sc N(ICD) Figure 6. As in Notch mutant clones (Egger et al., 2010;Reddy et al., 2010;Yasugi et al., 2010), the proneural wave is accelerated in clones expressing Notch RNAi, see also  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.

UAS-Pnt-P1 clones and N RNAi clones
To 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.
Video 5. Travelling proneural wave with ectopic activation of EGFR signalling within clones. The movie shows a simulation of the proneural wave model Equations 17-19 on the same lattice as in Video 3 but with four single-cell clones with constitutively active EGFR signalling. The movie corresponds to the snapshots shown in Figure 7a. DOI: https://doi.org/10.7554/eLife.40919.012 gain rate of signalling activity E. For sufficiently large gain rates, the function has three equilibrium points f Ã 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 f Ã 0 and f Ã 2 are stable (black dots), that is small perturbations of the signalling activity will drive the system back to these equilibrium pointshigher 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 f Ã 1 (white dot). If the system is only slightly perturbed around f ¼ f Ã 0 (black dot), the reaction dynamics will drive the system back towards this equilibrium, because degradation is stronger than self-activation in this region, r<0. However, if the signalling activity rises above f Ã 1 (e.g., due to diffusion from the neighboring cell), self-activation will outcompete degradation, r>0, and the signalling activity will be driven towards the other stable equilibrium f Ã 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). Non-dimensionalising the model by rescaling f ! f=F, t ! kt, and x ! ffiffiffiffiffiffiffiffi k=h p x, the corresponding reaction-diffusion system in one dimension is given by Here, l ¼ =ðkFÞ is the only (dimensionless) parameter and is the dimensionless reaction term.
Anticipating that the system gives rise to a travelling front with velocity c and a stationary front profile f, we make the ansatz fðx; tÞ ¼ fðx À ctÞ in Equation 6, which yields the ordinary differential equation where the prime denotes the derivative with respect to the argument y ¼ x À ct of f.
Multiplying this equation by f 0 and integrating from À¥ to ¥ yields an implicit expression for the velocity of the front, where DU ¼ lim y!¥ ½UðfðyÞÞ À UðfðÀyÞÞ with the potential UðfÞ defined by ¼ ÀqU=qf (see Appendix 1-figure 1b). Analytical solutions to Equation 8 only exist for special functional forms of the reaction term . To obtain analytical insights into the wave speed, we consider the limiting case of large Hill exponents, n ! ¥ 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 Q is the Heaviside step function with the convention Qð0Þ ¼ 1=2 (see Appendix 1- figure 1c). This function has three equilibria for l>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 co-moving reference frame by imposing the condition fð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 l, dropping all orders higher than l À1=2 and restoring the original parameter dependence by multiplying with the velocity scale ffiffiffiffiffi ffi hk p , we find Equation 12 implies that the front velocity increases with increasing diffusion constant h and synthesis rate as well as decreasing degradation rate k and decreasing activation threshold F. 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 one-component model reproduces the key features of more detailed models Previously, we sought to capture the spatio-temporal dynamics of EGFR signalling by a single component, described by the signalling activity f, 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 reaction-diffusion system where E denotes the concentration of bound EGFRs, R represents the concentration of Rhomboid, and S is the concentration of the secreted active form of Spitz (sSpi), Here, l E is the binding rate of EGFRs, l R is the synthesis rate of Rhomboid, and l S is the secretion rate of Spitz; the functions H E and H 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 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 diffusion-like effect to the entire E-R-S system. Therefore, the parameter 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 h describes the effective diffusibility, mediated by diffusion of component S. Formally, a connection to the one-component 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 q E =qt ¼ 0 and q R =qt ¼ 0 while the front is travelling. In this case, we can solve for E and R and obtain a closed equation for f S , where l l S and hðfÞ ¼

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 reaction-diffusion equations for the two fields f E and f L are given by (15) with the reaction terms where we have introduced the Hill function hðfÞ 1 À hðfÞ which describes an inhibitory effect. For simplicity, we have also considered identical unit concentration thresholds for activation and inhibition. Here, l E and l L are the gain rates of the components E and L, respectively. Note that in the absence of L (f L ¼ 0), the reaction term E reduces to the one of the one-component model given in Equation 6, E ðf E ; 0Þ ¼ ðf E Þ (Appendix 2- figure  1). Again, we consider a finite one-dimensional domain of length ' and no-flux boundary  Figure 2b, we choose an initial condition of the form fðx; 0Þ ¼ e Àx 2 =x 2 0 with x 2 0 ¼ 10h=k.) Again, this behaviour can be elucidated by studying the local reaction dynamics, which is now given by the vector field F ¼ ð E ; L Þ, defined by Equation 16) (Appendix 2- figure 1). Note that now only the point ðf E ; f L Þ ¼ ð0; 0Þ is a stable equilibrium. The black and white dots indicate points with f L ¼ 0 that satisfy E ¼ 0 as in the one-component system (compare to Appendix 1- figure 1a). However, note that only the point ð0; 0Þ also satisfies 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 self-feedback. Thus, the system finally returns to the fixed point ðf E ; f 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.
Here, the parameters i denote gain rates, k i denote decay rates, b denotes the basal gain rate of the component N, and the F i denote threshold levels for positive and negative feedbacks. The operatorsD andŜ are the lattice Laplacian and the sum over concentrations of neighboring lattice sites, respectively, and defined by ½Df x ¼ P y2Ux ðf y À f x Þ and ½Ŝf x ¼ P y2Ux f y , where U x is the set of neighbours of site x. The dynamics of the cell state 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, x has two stable equilibria: x ¼ 0 (neuroepithelium) and 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 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 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 ¼ 0 (neuroepithelial state) and ¼ 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 int leads to a transition of the cell state from neuroepithelial ( ¼ 0) to neuroblast ( ¼ 1) when (i) L exceeds the threshold level F int and/or (ii) E exceeds and N drops below the threshold level. All other parameters are given in Appendix 3-table 1. The system given by Equations 17-20 was simulated on a hexagonal lattice with circular geometry with a radius of 15 lattice sites. Initial conditions were localised elevated levels of E in those outer boundary cells that have angles between p=3 and 5p=3 as measured from the center of the circular lattice. The respective simulation panels show snapshots of the activity of EGFR signalling (green), L'sc expression (blue), Delta activity (magenta) and Notch activity (pink), as well as the cell state , for which black indicates neuroepithelium ( ¼ 0) and red indicates neuroblasts ( ¼ 1). Colour intensity indicates the local gene expression levels or signalling activities, respectively. The snapshots show the time t ¼ 12:5.

Robustness against lattice defects
To test whether the proposed mechanism is robust with respect to disordered lattice structures, we considered a site-diluted 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 diffusion-mediated 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