1. Computational and Systems Biology
  2. Developmental Biology
Download icon

A quantitative modelling approach to zebrafish pigment pattern formation

  1. Jennifer P Owen  Is a corresponding author
  2. Robert N Kelsh
  3. Christian A Yates  Is a corresponding author
  1. Department of Biology and Biochemistry and Department of Mathematical Sciences, University of Bath, Claverton Down, United Kingdom
Research Article
  • Cited 0
  • Views 576
  • Annotations
Cite this article as: eLife 2020;9:e52998 doi: 10.7554/eLife.52998

Abstract

Pattern formation is a key aspect of development. Adult zebrafish exhibit a striking striped pattern generated through the self-organisation of three different chromatophores. Numerous investigations have revealed a multitude of individual cell-cell interactions important for this self-organisation, but it has remained unclear whether these known biological rules were sufficient to explain pattern formation. To test this, we present an individual-based mathematical model incorporating all the important cell-types and known interactions. The model qualitatively and quantitatively reproduces wild type and mutant pigment pattern development. We use it to resolve a number of outstanding biological uncertainties, including the roles of domain growth and the initial iridophore stripe, and to generate hypotheses about the functions of leopard. We conclude that our rule-set is sufficient to recapitulate wild-type and mutant patterns. Our work now leads the way for further in silico exploration of the developmental and evolutionary implications of this pigment patterning system.

Introduction

Pattern formation - the process generating regular features from homogeneity - is a fascinating phenomenon that is as ubiquitous as it is diverse. It is a major aspect of developmental biology, with key exemplars including segmentation within the syncitial blastoderm of fruit flies (Clark and Peel, 2018), digit formation in the vertebrate limb (Tickle, 2006), and branching patterns in kidney and lung development (Davies, 2002).

Another key example, pigment pattern formation, the process generating functional and often beautiful distributions of pigment cells, represents a classic problem in both developmental and mathematical biology. Pigment patterns allow animals to distinguish between individuals within a group and identify those of different species and are an important characteristic for the survival of most animals in wild populations. Pigment patterns are striking. They form rapidly and, in many cases, autonomously, that is, the process relies on self-organisation and not internal body structures. Additionally, they often vary dramatically between even closely related species, therefore recognising similarities and differences in the development of these related species can allow us insight into the evolutionary change. Finally, pigment pattern formation is made experimentally tractable by the self-labelling nature of pigment cells.

The horizontal blue and gold stripes of zebrafish are now one of the best-studied examples of pigment pattern formation, especially at the level of underlying cellular mechanisms (Singh and Nüsslein-Volhard, 2015; Watanabe and Kondo, 2015; Patterson and Parichy, 2019). Zebrafish are amenable to observational studies, since all development takes place outside the mother and the skin is transparent. This, combined with the availability of multiple key mutants (affecting, for example, cell-type differentiation and patterning), and the development of innovative in vivo cell ablation and in vitro cell culture techniques, have provided a unique opportunity to investigate the cellular and molecular basis for pigment pattern formation experimentally (Eom et al., 2015; Budi et al., 2011; Ceinos et al., 2015; Yamanaka and Kondo, 2014; Eom and Parichy, 2017; Hamada et al., 2014; Fadeev et al., 2015; Inoue et al., 2014; Irion et al., 2014; Watanabe et al., 2006; Frohnhöfer et al., 2013; Parichy et al., 2009; Svetic et al., 2007; Mellgren and Johnson, 2006; Parichy et al., 2000b; Iwashita et al., 2006; Hirata et al., 2005; Kelsh et al., 1996; Lister et al., 1999; Parichy et al., 2000a; Maderspacher and Nüsslein-Volhard, 2003; Walderich et al., 2016; Patterson et al., 2014; McMenamin et al., 2014; Patterson and Parichy, 2013; Krauss et al., 2014; Parichy and Turner, 2003; Eom et al., 2012; Mahalwar et al., 2016; Mahalwar et al., 2014; Asai et al., 1999; Takahashi and Kondo, 2008).

The cellular composition of the stripes and how these become assembled has been well-described. Zebrafish generate, over a period of a few weeks and beginning around 3 weeks of age (Frohnhöfer et al., 2013; Parichy et al., 2009), a robust adult stripe pattern of alternating dark blue stripes and golden interstripes comprised of three different pigment-producing cell types: melanocytes, containing black melanin; xanthophores containing yellow and orange carotenoids and pteridines; and iridescent iridophores, containing guanine crystals within reflective platelets (Hirata et al., 2003Figure 1A).

WT stripe composition and development.

(A) An adult wild type (WT) fish. Stripes and interstripes are labelled according to their order of temporal appearance. X0 is the first interstripe to appear. 1D and 1V (D - Dorsal, V- Ventral) are the first two stripes to appear. X1D and X1V are the next two interstripes to appear and so on. Image reproduced from Frohnhöfer et al., 2013 and licensed under CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0). (B) Summary of pigment cell distribution in adult zebrafish. The cells in the xanthophore, S-iridophore, melanophore and L-iridophore layers consist of xanthophores and xanthoblasts, melanophores, S-iridophores and L-iridophores, respectively. Adapted from Hirata et al., 2003. (C) Schematic of WT patterns on the body of zebrafish. Stages PB, PR, SP, J+ correspond to developmental stages described in 3. Patterns form sequentially outward from the central interstripe, labelled X0, with additional dorsal stripes and interstripes labelled 1D, X1D, 2D, X2D from the centre (horizontal myoseptum) dorsally outward (similarly, ventral stripes and interstripes are labelled 1V, X1V, 2V, etc).

Of the iridophores, there are two types distinguished by their platelet distribution (Hirata et al., 2003); type L-iridophores, and type S-iridophores. Only type S-iridophores play a role in stripe formation (type L-iridophores appear later and are likely involved in pattern maintenance [Hirata et al., 2003]), and they appear in two different forms. In the light interstripes, S-iridophores appear in a ‘dense’ arrangement (dense S-iridophores), forming a continuous sheet, whilst in the dark stripes the cells are in a ‘loose’ arrangement and appear more widely spaced (loose S-iridophores) (Hirata et al., 2003; Fadeev et al., 2015). Pigment cells are found in the hypodermis below the dermis, and organised as layers of cells consistently stacked in the same order (Hirata et al., 2003). Starting from the deepest layer of the hypodermis just above the muscle and moving to the dermis, adult dark stripes consist of consecutive layers of L-iridophores, melanocytes, loose S-iridophores and xanthophores. Similarly, adult light interstripes are made up of layers of dense S-iridophores and xanthophores (Figure 1BHirata et al., 2003). The final striped pattern is generated by the self-organisation of xanthophores, melanocytes, loose and dense S-iridophores into the appropriate positions within the hypodermis.

Prior to the initiation of adult stripe formation zebrafish exhibit a larval pigment pattern, formed in the first 5 days of development. Embryonic pigment cells form a distinctive early larval pattern that is essentially complete by 5 days post-fertilisation (dpf) and remains unchanged until metamorphosis. This pattern consists of melanocytes in four stripes (dorsal to the central nervous system, within the horizontal myoseptum, dorsal to the gut and ventrally under the yolk; S-iridophores are found associated with three of these melanocyte stripes [Parichy et al., 2009; Frohnhöfer et al., 2013]). Xanthophores lie in a monolayer under the skin, filling the areas between the melanocytes above the CNS and extending ventrally to the level of the gut. Formation of the adult pattern involves replacement of melanocytes and S-iridophores with new cells derived from adult pigment stem cells. Early larval xanthophores dedifferentiate, forming unpigmented xanthoblasts that regain their proliferative ability and proceed to generate the adult xanthophores (McMenamin et al., 2014; Mahalwar et al., 2014; Parichy et al., 2000b); an unknown proportion of the latter may derive from de novo production from adult pigment stem cells (Kelsh et al., 2017). Xanthophore de-differentiation is complete by 21dpf, and early metamorphic melanocytes appear in a widely scattered distribution between 14 and 21dpf (Parichy et al., 2009), thus forming the initial metamorphic pattern (Figure 1C, stage PB). A key event in the initiation of adult pattern metamorphosis is the appearance of newly differentiated dense S-iridophores alongside the horizontal myoseptum. In response to the appearance of these S-iridophores, the first adult xanthophores are generated (Figure 1C, stage PR) by differentiation from xanthoblasts in this region, thus initiating the first interstripe, X0. Furthermore, metamorphic melanocytes begin to accumulate either side of this central interstripe, marking the first two stripes denoted 1D and 1V (Figure 1C, stage PR). Subsequently, S-iridophores proliferate rapidly and spread bidirectionally; at the edges of the interstripes they switch to a more scattered (less tightly-packed) form as they continue to spread dorsally and ventrally. Spreading loose S-iridophores transition back into dense S-iridophores at the locations of the future interstripes X1V and X1D (Figure 1C, stages PR to SP) (Mahalwar et al., 2014). Once S-iridophores aggregate at the next interstripe, the process starts again, that is, xanthophores differentiate in response to the dense S-iridophores and melanocytes accumulate either side of the new interstripe generating the subsequent stripe. This process of S-iridophore aggregation predetermining future interstripe locations and subsequent delamination in future stripe regions repeats until S-iridophores cover the domain and all stripes (between 4 and 5) and interstripes are fully formed (Figure 1C, stage J+).

In addition to the description of pattern development (Frohnhöfer et al., 2013), many studies have identified individual patterning mechanisms that contribute to stripe formation, although it is unclear whether these are sufficient to explain pattern formation. Stripe generation is complex and requires many interactions. During pattern metamorphosis, these interactions may determine cell birth (Mahalwar et al., 2014), cell death (Takahashi and Kondo, 2008), cell migration (Yamanaka and Kondo, 2014; Takahashi and Kondo, 2008; Patterson et al., 2014), long-distance communication, through stabilisation of elongated cellular projections (Eom and Parichy, 2017; Eom et al., 2015), as well as the shape transitions of S-iridophores (Fadeev et al., 2015). During this period, there is also simultaneous two-dimensional domain growth (Parichy et al., 2009). The pattern is formed by cell–cell interactions of all three pigment producing cell types: melanocytes, xanthophores and S-iridophores. Without any one of these cell types, pattern formation is disrupted (Frohnhöfer et al., 2013; Patterson and Parichy, 2013).

Mathematical modelling has been a complementary tool in assessing possible patterning mechanisms. Until the last few years, these studies have focused on melanocytes and xanthophores, neglecting S-iridophores. The most commonly used mathematical paradigm for stripe formation takes the form of a Turing reaction-diffusion model. In these representations, melanocytes and xanthophores diffuse and interact via a few long- and short-range ‘reactions’. This class of model typically rely on a small number of parameters which, upon being altered, can generate a diverse range of patterns. Minimal models such as these have the benefit that they are sometimes analytically tractable, allowing a deep understanding of the model. However, a potential limitation is that parameters do not always have a clear biological interpretation which, can sometimes make it difficult to link parameters to measurable data. In the context of zebrafish stripe formation, these models have not yet incorporated S-iridophores (Watanabe and Kondo, 2015; Kondo, 2017; Painter et al., 2015; Bloomfield et al., 2011; Binder and Simpson, 2013; Volkening and Sandstede, 2015; Kondo, 2017; Nakamasu et al., 2009; Moreira and Deutsch, 2005; Bullara and De Decker, 2015; Yamaguchi et al., 2007; Asai et al., 1999). They suggest that the role for iridophores is restricted to simply orienting stripes (Volkening and Sandstede, 2015; Nakamasu et al., 2009; Binder and Simpson, 2013). New biological observations demonstrate that S-iridophores play a fundamental role in body stripe formation (Singh and Nüsslein-Volhard, 2015; Frohnhöfer et al., 2013; Patterson and Parichy, 2013). In particular, it has been shown that without S-iridophores, spots of melanocyte aggregates form instead of stripes, which is contrary to what these Turing reaction-diffusion models predict. These findings have paved the way for more detailed modelling, such as that of Volkening and Sandstede, 2018, who demonstrated (using an off-lattice individual-based model) the need for understanding S-iridophore behaviour when representing all three cell-types. For these reasons, we consider an inclusive modelling approach, incorporating the crucial cell-type S-iridophores and the full range of interactions depicted above.

Here, in a bottom-up approach, we hypothesise that the current biological understanding is sufficient to explain the major aspects of pigment pattern development and construct a model to test this. In particular, we construct an agent-based model incorporating all three pigment cell-types and their documented cellular interactions. We use observations of a set of three mutants that each lack an individual cell-type, plus the three double mutant combinations lacking pairs of cell-types, to deduce the key rules likely underpinning S-iridophore dynamics. Combining these assumptions with experimentally verified biological mechanisms in the literature, we generate a working model of adult pattern formation. We then run simulations for wild type (WT) and these mutant fish. We show that in each case our model correctly predicts the patterns observed in vivo, and that pattern development displays multiple quantitative matches to that in vivo using a parameter sampling methodology to demonstrate the robustness of these patterns to parameter variation. In an independent test of the model, we simulate mutants with pigment pattern defects caused by changes other than to the presence of pigment cell-types, and show that these too are successfully matched in silico by our model.

Our work demonstrates that current biological understanding, alongside simple assumptions about S-iridophore behaviour, is sufficient to explain adult pigment pattern formation in WT and multiple mutants. Our work reinforces the growing realisation in the field that the previously neglected S-iridophores are crucial for stripe formation, suggests a minimum set of their rules, and reveals unexpected subtleties to the phenotypic impact of the well-studied leo mutant.

Materials and methods

Modelling overview

Request a detailed protocol

We build our model with direct reference to the known biology. We model five cell types as individual agents: melanocytes (M), xanthophores (X), xanthoblasts (Xb), the unpigmented precursor cell to xanthophores) and S-iridophores in either dense or loose form (Id, Il, respectively). These are the cells we deem from the literature to be crucial for successful pattern formation. We do not directly model L-iridophores, since these appear after the adult pattern is formed and are more likely involved in pattern maintenance (Frohnhöfer et al., 2013). Unlike previous models of stripe formation (Nakamasu et al., 2009; Bullara and De Decker, 2015; Volkening and Sandstede, 2015; Painter et al., 2015; Bloomfield et al., 2011; Volkening and Sandstede, 2018), we include xanthoblasts as an independent cell-type in our model. This is because the larval xanthoblasts appear principally by dedifferentiation of the embryonic xanthophores, and most metamorphic xanthophores arise from the larval xanthoblasts (Mahalwar et al., 2014; McMenamin et al., 2014; Budi et al., 2011; Singh et al., 2014; Dooley et al., 2013), whilst xanthoblasts that do not re-differentiate into xanthophores persist in the stripe regions where they play a role in consolidating melanocytes into stripes.

Zebrafish pattern formation generates distinct pigment cell layers in the hypodermis (Figure 1B) a melanocyte, xanthophore and S-iridophore layer (Hirata et al., 2003; Hirata et al., 2005). For consistency, we model each of the three layers as independent, two-dimensional lattice domains throughout pattern formation (Figure 2A).

Model setup and simulation.

(A) An example model setup. The domain is made up of three layers. Layer 𝑿 which contains yellow X (yellow circles) and unpigmented Xb (clear circles with black outline). Layer 𝑰 which contains silvery Id (white circles) and blue Il (blue circles). Layer 𝑴 which contains black M (black circles) only. Each lattice site on each of the respective layers contain at most one cell at any given time. The layers are stacked on top of each other as seen in real fish. We note that in our simulations the ordering of the layers does not play any role in determining pattern formation. (B) Schematics of model implementation. (1) The model is initialised as described in Appendix 1. (2) The model is checked for the requirements for a fixed event to occur (e.g. new cells differentiating at a given time). If there is a fixed event to be implemented, the algorithm moves to stage 3. Otherwise, the algorithm passes to stage 4. (3) The model is updated according to the fixed event and algorithm returns to stage 2. (4) The propensity functions - probabilities for all other events to occur α1(t),α15(t) are calculated. (5) Random numbers u0,u1U(0,1) are generated. (6) Numbers u0,u1 are used to determine the time τ until the next event and which event will be implemented. (7) The model configuration is updated. (8) The algorithm checks if the end criterion is satisfied, that is in the case of WT, that ΩSL=13.5mm. If so, the algorithm finishes. Otherwise the algorithm continues, returning to stage 2. (9) The simulation completes.

Agents representing X and Xb, M, Id and Il occupy lattice sites, within xanthophore, melanocyte and S-iridophore domains respectively (Figure 2A). To account for the different packing densities of the cell types, lattice sites within the xanthophore and S-iridophore model layers are half the width and length of melanocyte sites size. This packing density does not have an impact on pattern formation, but, is included for biological realism. (For more details, see Appendix 1). Within each layer, volume exclusion properties hold: no two agents can occupy the same site at any one time (i.e. cells do not overlap).

The system is initialised to represent a typical WT fish shortly after the start of adult pigment pattern development (≈25 dpf). We set the domain height to be 1 mm, since this is the approximate height of the fish at 25 dpf (Supplementary file 3 for details), we set the domain length to be 2 mm, representing approximately one-third of the full length, from the tip of the snout to the start of the tail, at 25 dpf, and thus equivalent to the trunk (Parichy et al., 2009). We populate the domain itself at t=0 as an approximation of the observed larval pattern at 25 dpf (Frohnhöfer et al., 2013). At this time, there is a central stripe of dense S-iridophores along the horizontal myoseptum, scattered melanocytes and de-differentiated xanthophores (xanthoblasts) scattered across the domain. We model this by populating the central three rows of the S-iridophore layer with dense S-iridophores, and by distributing melanocytes uniformly at random into sites within the melanocyte domain at density 0.04 and xanthoblasts uniformly at random into sites in the xanthophore domain at density 0.4.

The model is then updated according to the Gillespie algorithm (Gillespie, 1977). An overview of how the model is updated is given in Figure 2B and can be described as follows. At any given time t, the model is first assessed for meeting the criteria of a fixed event. Fixed events are all biologically determined events that occur once at a fixed time. For example at the start of pattern formation, the appearance of dense S-iridophores along the horizontal myoseptum is a fixed event. If the model meets the criteria, the fixed event occurs, is subsequently marked as complete and the simulation continues. If no fixed time event is to be implemented then one of 15 possible continuous time events is attempted. To do this, we treat all the potential actions, (for example cell birth or domain growth [as described in Section "Modelling assumptions"]), as individual ‘events’, each with an exponentially distributed waiting time which corresponds to their rate of occurrence (as specified in the literature Supplementary file 4). To update the model at any given time t=T, an exponentially distributed waiting time; τ is generated until the next possible ‘event’ occurs (based on the rates of all of the possible events). Next, a random number u1U(0,1) determines which event occurs based on the relative probability of each event occurring. Once an event is chosen, the domain is updated accordingly: if conditions required for that event to occur are met, the event is implemented, whereas if they are not then there is no change. Time is also now updated to t=T+τ. This process repeats until we reach the end of pattern metamorphosis, defined by the simulated field standard length reaching approximately 13.5 mm (Supplementary file 3). The stochastic nature of our algorithm means that in any given simulation, the final pattern and its individual development will be inherently different to any other simulation, just as in real fish. Events incorporated into our model include all processes involved in the self-organisation of pigment cells during pattern metamorphosis as well as uniform domain growth with rate 0.13 mm per day in horizontal axis and 0.033 mm per day in the vertical axis (Parichy et al., 2009). These events are described in more detail in Section "Modelling assumptions".

Cells interact in the fish skin at both short (neighbouring cells) and long (up to half a stripe width ≈0.25 mm) range, with interactions thought to use direct contact through cellular extensions (filopodia, dendrites, or longer airenemes). In our model, uniform disks, with radii on the order of the distance between 2 cells (≈0.04 mm) account for short-range interactions (Figure 3A–D), and an annulus with an outer radius of 0.24 mm (12 cells) and inner radius of 0.22 mm, (11 cells) represent long-range dynamics (Figure 3E–H). We allow cell interactions across different layers (as in real pattern formation). Cells that are chosen for movement can move into one of eight sites local to them. The probability of movement in one of the eight direction is biased according to how attracted or repelled the focal cell is to its local neighbours (Figure 3I–J). For more detail about how short- and long-range interactions are implemented see Appendix 1. See Supplementary files 4, 5 and 6 for a detailed justification of the rates, interaction types and parameter values, respectively.

Simulating short- and long-range interactions.

(AF) Comparing the number of cells in the short (AD) and long (EF) range distance (0.04 mm) from a central site, on different domain types. M are represented as black circles. X are represented as yellow circles. (AB) A visualisation of sites (marked in red) local to central melanocyte on the (A) melanocyte domain and (B) xanthophore domain.(CD) A visualisation of sites (marked in red) local to central xanthophore on the (C) xanthophore domain and (D) melanocyte domain. (EF) A visualisation of sites (marked in red) long-distance to central melanocyte on the (E) melanocyte domain and (F) xanthophore domain.(GH) A visualisation of sites (marked in red) long-distance to central xanthophore on the (G) xanthophore domain and (H) melanocyte domain. (I) Melanophore movement with respect to local neighbours. If a melanophore located in a central position (marked in dark grey) attempts to move, the melanophore will consider neighbours in all sites marked in red on the melanophore domain (left), xanthophore and S-iridophore domain (right). (J) Movement on xanthophore/iridophore domain with respect to local neighbours. If a xanthophore/iridophore located in a central position (marked in dark grey) attempts to move, the xanthophore/iridophore will consider neighbours in all sites marked in red on the xanthophore and S-iridophore domain (left), melanophore domain (right).

Modelling assumptions

In this section, we describe our modelling assumptions with regard to cell–cell interactions. These assumptions include all the known interactions between melanocytes, dense S-iridophores, loose S-iridophores, xanthophores and xanthoblasts, as well as some predictions about S-iridophore behaviour which have not been experimentally investigated in the literature. Apart from those involving S-iridophores, all the interactions and wherever possible their quantitative properties (strength, frequency etc) come directly from the literature, and are summarised in Figure 4G, 1-14, and described in Supplementary file 5. These include interactions influencing the movement, proliferation, differentiation and death of all cell types. These are represented explicitly in the model in as biologically realistic a manner as possible, at their determined rates.

‘Missing cell’ mutants and cell-cell interactions implemented.

(AF) Adult phenotype of the set of fish used to deduce S-iridophore interactions. (A) WT fish, (B) shady (shd), (C) nacre(nac), (D) pfeffer (pfe), (E) nac;pfe, (F) shd;pfe. (G) A representation of all the interactions implemented in the model. See Supplementary files 4 and 5 for a detailed justification of the rates and interaction types, respectively. See Supplementary file 6 for a detailed justification of all parameters. The direction of the arrow/inhibition sign indicates which cell is acting on which. For example in assumption 13, the blue arrow from melanophore to xanthophore indicates that melanophores attract (green background) xanthophores. The purple inhibition arrow indicates that the xanthophores repel melanophores. OR statements in assumptions 15–18 are logical OR statements. For example, conditons 15 and 16 indicate that a dense S-iridophore may become loose if there are melanophores in the short range OR, there are no xanthophores in the short range AND many xanthophores in the long range. (H) A visualisation of measurements HAA and SL. (I) Plots of SL versus days post-fertilisation (dpf) for 40 model realisations. Each coloured line is a single realisation. Black squares are SL versus days post-fertilisation extracted from Parichy et al., 2009 in real fish given in Supplementary file 3. Red diamonds are the mean SL versus days post-fertilisation from 40 simulations. (J) Plots of SL (mm) versus HAA (mm) for 40 model realisations. Each coloured line is a single realisation. Black squares are SL versus HAA extracted from Parichy et al., 2009 given in Supplementary file 3. Red diamonds are the mean HAA versus SL from 40 simulations. Error bars are one standard deviation. The model agrees well with the data in both cases. Images (AF) reproduced from Frohnhöfer et al., 2013 and licensed under CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0).

The interactions involving S-iridophores have not been well-characterised experimentally, so we have developed our own predictions based on the literature describing S-iridophore behaviour during pattern metamorphosis. It has been shown using clonal cell analysis that during pattern metamorphosis S-iridophores spread across the skin of the zebrafish bidirectionally by proliferation of existing cells (between once and twice per day) combined with quick migration (Mahalwar et al., 2014). We further predict that dense S-iridophores show a directional bias towards xanthophores in the short range. We propose that dense S-iridophores are attracted in the short range to xanthophores since they are highly associated with each other in each of the mutants and in WT (Frohnhöfer et al., 2013), and this mutual attraction may be important for interstripe consolidation. Furthermore, loose S-iridophores show a directional bias away from other loose S-iridophores in the short range. We propose that loose S-iridophores are repelled by other loose S-iridophores as this would facilitate the prompt spreading of loose S-iridophores across the stripe regions.

Interestingly, as the S-iridophores spread they switch between a loose and dense form, predetermining the positioning of stripes and interstripes consecutively. In the loose form S-iridophores are spread and stellate in appearance. In contrast, in the dense form, S-iridophores are compact. The transition between the two types appears to be interchangeable. When dense S-iridophores initially spread beyond the boundary of the first X0 interstripe, they can later change to loose type (Fadeev et al., 2015). Similarly when loose S-iridophores reach an interstripe region, they can aggregate and form dense S-iridophores. It is not clear exactly what causes these shape transitions physically, and this is not a question we address here. It has, however, been shown that loss of Tjp1a function in sbr mutants compromises the transition of S-iridophores from dense to loose state, suggesting that Tjp1a contributes to regulation of the molecular switch that regulates S-iridophore shape changes during their dispersal (Fadeev et al., 2015). We envisage iridoblasts as initially differentiating in a dense form along the horizontal myoseptum, proliferate, migrate and spread, later de-differentiating and then re-differentiating into the opposite form dependent on their location with respective to other cell types (melanocytes and xanthophores).

Here, we hypothesise how the cell types affect S-iridophore type (loose or dense). The cause of these transitions is largely unknown, however, it has been suggested to be dependent on signals from melanocytes and xanthophores transmitted by gap junctions (Irion et al., 2014; Fadeev et al., 2015). In order to investigate this, we consider a primary set of six mutants known to prevent the formation of one or more individual pigment cell-type. We use these to define the contribution and nature of S-iridophore interactions in pattern formation, by considering the outcomes in fish lacking each of the three cell types. Specifically we consider:

  • Mutants lacking S-iridophores: The gene shady (shd) encodes zebrafish leukocyte tyrosine kinase (Ltk) which plays a role in S-iridophore specification (Lopes et al., 2008). As a result, strong shd mutants lack all S-iridophore types. The resultant adult pattern consists of a widened X0 region of xanthophores, which are flanked dorsally and ventrally by melanocytes organised as spot-like clusters in a sea of xanthophores, forming broken stripes (Figure 4B).

  • Mutants lacking melanocytes: The gene nacre (nac) encodes the transcription factor Mitfa (Lister et al., 1999). nac mutants lack melanocytes throughout embryonic and larval development (Lister et al., 1999). As a result, stripes do not form properly and the adult phenotype consists of a prominent X0 interstripe of dense S-iridophores and xanthophores with irregular borders, accompanied by spots of dense S-iridophores and xanthophores ventrally. The rest of the flank is filled with loose form S-iridophores (Figure 4C).

  • Mutants lacking the xanthophore lineage: Gene pfeffer (pfe) (alternatively known as salz (sal)) encodes colony stimulating factor one receptor (csf1ra) that is expressed and required specifically in xanthophores (Frohnhöfer et al., 2013; Patterson and Parichy, 2013; Parichy et al., 2000b). In the adult fish of strong alleles, xanthophores are almost absent in embryos, and absent in adults. The resultant adult pattern consists of a spotted melanocyte stripe pigmentation of normal alignment which fades out into a ‘salt and pepper’-like pattern more posteriorly (i.e., in the tail) (Figure 4D). Melanocyte spots are associated with loose S-iridophores. In the regions lacking melanocyte aggregation (the ‘salt-and-pepper’ region), S-iridophores take a dense form, with melanocytes scattered at very low density, an arrangement never seen in WT patterns.

  • Double mutants of the aforementioned mutant types: nac;pfe, nac;shd and shd;pfe (Figure 4E,F depict the adult phenotypes of nac;pfe and shd;pfe respectively, there is no image available for shd;pfe). These mutants lack two of the aforementioned cell types. The resultant adult pattern is a uniform distribution of the remaining cell type (Frohnhöfer et al., 2013). These mutant phenotypes demonstrate that zebrafish stripe formation is not determined by an underlying pre-pattern, but instead is generated by cell-cell interaction.

Upon evaluating these mutants, we make the following deductions about S-iridophore shape transitions during pattern formation:

  • S-iridophores are initially dense and cannot change shape autonomously. This is based on observations of mutants nac;pfe which only contain S-iridophores and in which the adult phenotype consists of dense S-iridophores in a coherent sheet across the domain (Frohnhöfer et al., 2013). In contrast, pfe and nac both exhibit loose and dense S-iridophores (Frohnhöfer et al., 2013), suggesting that both melanocytes and xanthophores are capable of facilitating S-iridophore shape transitions.

  • Melanocytes in the short range promote the transition of dense to loose, conversely, a lack of melanocytes in the short range promotes the transition of S-iridophores from loose to dense. We propose these interactions for the following reasons. Firstly, in pfe and WT, dense S-iridophores are associated with lack of melanocytes, for example within the interstripes, whilst loose S-iridophores are associated with melanocytes, for example in the stripe region. Since we predict that melanocytes are required for dense S-iridophores transition the simplest assumption is that melanocytes promotes dense to loose transitions in the short range. Since loose S-iridophores can re-aggregate to dense form in pfe we assume that this signal is bidirectional and therefore a lack of melanocytes promotes the loose to dense transition.

  • Xanthophores in the long range and lack of xanthophores in the short range promotes the transition from dense to loose; conversely, a lack of xanthophores in the long range as well as many xanthophores in the short range, promotes the transition from loose to dense. We propose these interactions for the following reasons. In nac and WT, dense S-iridophores are associated with xanthophores, whilst loose S-iridophores are associated with a lack of xanthophores (Frohnhöfer et al., 2013). Since S-iridophores initially appear in dense form and become loose for example in nac, when there are xanthophores in a low density local to S-iridophores and high density in the far range, we predict it is this combination that promotes the transition of dense to loose in the long range. Since in nac, S-iridophores can transition back from loose to dense when the local xanthophore density is high and far xanthophore density is low, we assume the opposing interaction is also true (Frohnhöfer et al., 2013).

These descriptions are summarised in Figure 4G, 15-18. We note that in each of these cases variations of these interactions were already hypothesised by Frohnhöfer et al., 2013. However, since their predictions did not distinguish loose and dense S-iridophores and did not indicate transition mechanisms, their predictions though similar, are extended here to incorporate these differences. The predictions we describe are the simplest possible for generating the patterns observed in the aforementioned set of fish, upon removing any one of these interactions, the model fails to generate the robust patterns we will describe (Figure 11 and Section "Necessity of S-iridophore assumptions").

Comparing simulated fish with real data

Request a detailed protocol

In order to validate our model, we compare different aspects of our simulation (size, spatial distributions of cells, numbers of melanocytes, stripe and interstripe width) with real fish at different developmental stages. In real fish, developmental stages are categorised according to the standard length (SL) of the fish (Figure 4H; Supplementary file 3Parichy et al., 2009). For consistency, we calculate the ‘stage’ of our simulations using the length of our domain and a simple calculation to generate a simulated SL (see Appendix 1). This allows us to make a direct comparison between the range of sizes obtained in model simulations and the natural range in zebrafish HAA and SL. As a test of validity of this measure, Figure 4I and Figure 4J demonstrate 40 plots of simulated SL versus days post-fertilisation (dpf) and simulated HAA versus SL, respectively, compared with the averaged data (Parichy et al., 2009). These figures demonstrate that whilst growth rates are variable within simulations (as seen in real fish), the mean of our simulated rates approximately matches that in real fish.

Results

Modelling simulations

Having deduced this minimal rule-set from the literature and our further predictions from the phenotypes of the selected primary set of mutants, indicated in Section "Modelling overview" we use these as the basis for our model. We use our model to generate stochastic simulations of pigment pattern formation corresponding to the the period of adult metamorphic pigment pattern formation, during which the SL extends from 7.6 mm to 13.5 mm. We note that adult pattern metamorphosis and the appearance of metamorphic S-iridophores starts earlier, around 6–7 mm SL (Parichy et al., 2009). We initialise later at 7.6 mm as by this time the skin lying over the horizontal myoseptum is populated with an initial stripe of dense S-iridophores. We intialise our model accordingly to match this. Subsequently, we first assess the ability of the model to reproduce natural growth at a quantitative level, and then to generate the WT pigment pattern, both qualitatively and quantitatively. We then go on to simulate conditions corresponding to our primary set of mutants, considering the qualitative fit to the published patterns. To test robustness of the patterns, we provide a rigorous robustness analysis by carrying out one hundred repeats of the WT simulation and ‘missing cell’ mutants with perturbed parameter values chosen uniformly at random from the range 0.75–1.25 of their described value and show that in each case, the appropriate pattern is preserved. Finally, in a more rigorous test of the predictive power of our model, we explore three further mutant phenotypes that had not been considered in deriving the model’s rule-set.

Simulation of WT pattern

In this section, we compare qualitatively our simulations of WT fish. For WT simulations, the model rules are given in Section "Modelling overview". Figure 5A–D depicts WT development, while Figure 5A’–D’ (Video 1) shows a representative simulation using the model described by the rules in Section "Modelling overview". The simulations reproduce qualitatively most aspects of the biological pattern. The model is initialised at stage PB. At stage PR, we begin to see an accumulation of melanocytes either side of the initial interstripe and differentiation of new xanthophores. Furthermore, we see the development of 1D and 1V stripe regions and delamination of S-iridophores from dense to loose form at the edges of interstripe X0. At stage SP, we observe the spreading of loose S-iridophores across the two developing stripe regions. Finally at stage J+, we see three interstripes alternating with five dark stripes. The final pattern matches the stripes seen in the real WT fish and the cellular component of dark stripes (X, Il, M) and light interstripes (Id, X) matches the composition of pigment cells in real fish (Figure 1C). We emphasise that the simulations presented here (as well as in future sections) are a representative example of the model output.

Representative simulation of WT, shd, nac and pfe.

(AD), (EH), (IL), (MP) WT, shd, nac, pfe development respectively and (A’D’), (E’H’), (I’L’), (M’P’) corresponding model simulation. Red arrows: melanocytes in the skin layer, modelled in (A). White arrows indicate the embryonic pattern of melanocytes in four stripes, that are deeper than the skin level and are consequently not included in the model. See Videos 14 for representative examples of these simulations in movie format. Scale bar is 0.25 mm for all images. Experimental images (AD), (EH), (IL), (MP) reproduced from Frohnhöfer et al., 2013 and licensed under CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0).

Video 1
Simulated development of WT.

Robustness of the model

Due to the abundance of parameters and cell-cell interactions necessary to capture what is known biologically about zebrafish pigment pattern formation, it is not feasible to perform an exhaustive parameter sweep to demonstrate the robustness of the model. Instead, as a test of robustness, we perform a rigorous robustness analysis by carrying out one hundred repeats of the WT simulation with perturbed parameter values chosen uniformly at random from the range 0.75–1.25 of their described value. The precise value of each parameter is sampled uniformly from this region, independentally for each parameter and each repeat. Twenty of these randomly sampled repeats are given in Figure 6. We observe that for all one hundred repeats that small perturbations to the rates still generate consistent striping, demonstrating the robustness of the model.

Example WT simulations at stage J+ when the parameters governing the rate of proliferation, movement, differentiation and death are perturbed slightly.

Each square is an example WT simulation at stage J+ where each rate parameter is perturbed to 1+x times its normal value. The value x is chosen uniformly at random from the set [-0.25,0.25].

Simulation of ‘missing’ cell mutants

In the next four sections (3.1.2–3.1.2), we compare qualitatively our simulations of mutants lacking one or more cell types. For the case of generating these mutants, we simulate the same WT model except we remove the appropriate cell type from the initial conditions and turn off cell birth of that cell type to match the mutation. For example, in shd we remove S-iridophores from the initial conditions and switch off S-iridophore birth. No other changes are made. For more information about mutant implementation see Appendix 2. These mutants often display similar features to WT fish; however, some aspects of the stripe formation are incomplete. In order to describe these differences, with reference to pfe, shd and later in sbr, we define pseudo-stripes as the spots of melanophore aggregates that appear in a stripe-like orientation reminiscent of that in WT fish. We describe the pseudo-stripes in the order they appear as in WT fish. For example, we define pseudo-1D and pseudo-1V to be the first pseudo-stripes.

We demonstrate here the capability of our model simulations to reproduce qualitatively the pattern development of these mutants. For each mutant, we describe the initialisation of the model domain to match the fish at stage PB as well as all of the similarities observed between our model outputs and real fish at the following three developmental stages, PR, SP and J+. Finally, we describe the variation between many repeats of the model and how this correlates with real same-type siblings.

The shady mutant

At stage PB (Figure 5E,E’), we populate the domain with some randomly dispersed melanocytes at a lower density than that in WT (Frohnhöfer et al., 2013) and some randomly dispersed xanthoblasts at the same density as WT (Frohnhöfer et al., 2013). At stage PR (Figure 5F,F’), we observe some melanocytes beginning to differentiate in the usual 1D and 1V stripe regions. At stage SP (Figure 5G,G’), we observe the accumulation of melanocytes around the 1D and 1V stripe regions with a central stripe of xanthophores. Finally, at stage J+ (Figure 5H,H’), we observe two horizontal pseudo-stripes of melanocyte spots surrounded by xanthophores. We found that in 100 simulations, 100% of shd stage J+ mutants observed two pseudo-stripes (1D and 1V) just as in Figure 5H. See Video 2 for a movie of the simulation.

Video 2
Simulated development of shd.

Moreover, pseudo-stripes varied in how stripe-like they were as observed in real fish (Frohnhöfer et al., 2013). As a measure of this, we calculated the longest stretch of melanophores in a row without any significant breaks over 100 simulations. This gives an indication of the widest ‘spot’ or ‘pseudo-stripe’ of melanophores in a simulation. We found that on the average, the mean of widest spot width over one hundred simulations was 0.18 of the simulated length. The widest spot width in 100 simulations was 0.43 of the simulated length, demonstrating the variance in pseudo-stripe length without break.

The nacre mutant

At stage PB (Figure 5I,I’) we populate the domain such that there is an initial stripe of dense S-iridophores and randomly dispersed xanthoblasts at the same density as in WT. At stage PR (see Figure 5J,J’) we see the appearance of newly differentiated xanthophores associated with the dense S-iridophores in the initial X0 interstripe. At stage SP (Figure 5K,K’) we observe the switch of dense S-iridophores to loose form and the subsequent spreading of loose S-iridophores. Finally at stage J+ (Figure 5L,L’) we observe the jagged edges of the usually straight X0 and the formation of a second pseudo interstripe some distance below X0 just as in real nac fish (Figure 5L). See Video 3 for a movie of the simulation.

Video 3
Simulated development of nac.
The pfeffer mutant

At stage PB (Figure 5M,M’), we populate the domain with a central stripe of dense S-iridophores and randomly dispersed melanocytes at the same density as observed in WT (Frohnhöfer et al., 2013). At stage PR (Figure 5N,N’), we observe the arrival of melanocytes into the prospective 1D, 1V stripe regions that is less pronounced compared with WT simulations. At stage SP (Figure 5O,O’), we observe the accumulation of newly differentiated melanocytes into aggregates in prospective stripe regions 1D and 1V. Finally, at stage J+, (Figure 5P,P’) we observe the aggregation of melanocytes (associated with loose S-iridophores) into spots, surrounded by a sea of dense S-iridophores peppered with black melanocytes. In one hundred simulations, the median number of pseudo-stripes at stage J+ in these repeats was four, as WT. This is consistent with observations of pfe mutants, which typically show the same number of pseudo-stripes and -interstripes as WT fish (Frohnhöfer et al., 2013). We observe higher conservation of striping than in simulated shd mutants as observed in real fish. For example, in one hundred simulations, the average longest stretch of melanophores in any given simulation was 0.6 of the full length. See Video 4 for a movie of the simulation.

Video 4
Simulated development of pfe.

As a test of robustness, we perform a rigorous robustness analysis by carrying out one hundred repeats of the mutant simulations with perturbed parameter values chosen uniformly at random from the range 0.75–1.25 of their described value as in Section "Simulation of WT pattern". Ten of these randomly sampled repeats are given in Appendix 4—figure 1. We observe for all one hundred repeats and in all three mutants, that small perturbations to the rate parameters still generate consistent patterning, demonstrating the robustness of the model.

Double mutants; shd;pfe, shd;nac, nac;pfe

Lastly, we consider the double mutants. Figure 7A and Figure 7B depict adult patterns in shd;pfe and nac;pfe respectively. There is no image available for shd;nac adult or for the development of the aforementioned mutant phenotypes but it has been described in the literature that in all of the double mutants, the remaining cell type, by adulthood, fills the domain uniformly (Frohnhöfer et al., 2013). Figure 7A’–C’ show a representative simulation for the mutants shd;pfe, nac;pfe, shd;nac respectively. In all cases of our model simulations, we observe that by stage J+ the remaining cell type begins to fill the domain. For example, in nac;pfe S-iridophores in dense form cover most of the flank by stage J+ (Figure 7B’).

Other mutant fish simulations.

(AC) Real stage J+ shd;pfe, nac;pfe, shd;nac mutants, respectively. (A’C’) Model simulations at stage J+ for shd;pfe, nac;pfe, shd;nac mutants, respectively. (DF) Adult phenotype of selected (Drse, (E) sbr, (F) cho mutant. (GJ), (KN), (OQ) rse, sbr and cho development, respectively. (G’J’), (K’N’), (O’Q’) rse, sbr and cho model simulations, respectively. (R) seurat at stage J. (R’) seurat model simulation at stage J. Scale bar is 0.25 mm for all images. See Videos 57 for representative examples of these simulation in movie format. Experimental images (AC), (D), (F) (GJ) and (OQ) are reproduced from Frohnhöfer et al., 2013, (E), (KN) from Fadeev et al., 2015, (R) from Eom et al., 2012 and are all licensed under CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0).

Simulation of other mutants

In Section "Simulation of WT pattern", we demonstrated that our proposed model reproduces the WT, single and double mutant patterns and thus is sufficient to explain pattern formation in the skin. In this section, we perform a more stringent test of the model’s completeness, by asking whether it can successfully simulate the outcomes of a set of pigment pattern mutants which were not used to deduce the rules underpinning our model. Since we were particularly interested to test our predictions of the rules relating to S-iridophore interactions, our secondary set comprises mutants with S-iridophore-related phenotypes: rose (rse) homozygotes, which show a reduction of S-iridophore numbers; schachbrett (sbr) homozygotes, which show a delay in S-iridophore shape transitions from dense to loose and choker (cho) homozygotes, in which the absence of the horizontal myoseptum prevents the formation of the initial dense X0 band of dense S-iridophores (Figure 7D–F). In the next few sections (3.1.3–3.1.3), we demonstrate the capability of our model to reproduce quantitatively the patterns of these mutants. In each section, we first describe the nature of the mutation and the way in which we adapt our WT model to simulate the mutants. We describe the similarities of our model simulation with real fish at the different developmental stages considered.

The rse mutant

Rose (rse), encodes the Endothelin receptor B1a (Krauss et al., 2014) and has been shown to acts cell-autonomously in S-iridophores; homozygous mutants result in a reduction of S-iridophores to approximately 20% of that seen in WT (observed in stage PB and adult fish [Frohnhöfer et al., 2013]). Consequently, adult fish show two broken dark stripes (reduction from four) bordering a widened X0 interstripe region. (Figure 7D). To simulate the rse mutant, we changed the number of initial dense S-iridophores at stage PB to one fifth of its usual number as observed in real fish at stage PB (Figure 7G,G’Frohnhöfer et al., 2013). At stage PR (Figure 7H,H’), there is a strong reduction in melanocyte number compared to WT (Figure 7I,I’) and we observe that dense S-iridophores spread less far from the horizontal myoseptum. At stage SP, (Figure 7J,J’) the stripe boundaries at X0 are poorly defined, and dense S-iridophores are still largely associated with the X0 region, with only a few loose S-iridophores appearing at the dorsal and ventral margins. At stage J+ (Figure 7K,K’), the stripe boundaries at X0 are more distinct, but the dark stripes are thinner and partly fragmented. See Video 5 for a movie of the simulation.

Video 5
Simulated development of rse.
The sbr mutant

The sbr gene encodes Tight Junction Protein 1a (Tjp-1a), which is expressed cell autonomously in dense S-iridophores (but not in loose S-iridophores) and truncated in sbr mutants; in adult sbr mutants S-iridophore shape transition from dense to loose is delayed (Fadeev et al., 2015). As a result, adult fish exhibit interrupted dark stripes, generating a pattern reminiscent of a checkerboard (Figure 7E). Figure 7K–N depicts sbr early pattern development. During adult pigment pattern formation, differences from normal WT development are not seen until ≈10 mm SL (Fadeev et al., 2015), (SP stage) at which point instead of dense S-iridophores transitioning to loose S-iridophores at the edge of the interstripe, in sbr, dense S-iridophores remain dense and spread over melanocytes dorsally and ventrally bidirectionally . At later stages, some dense S-iridophores do switch to loose S-iridophores. See Video 6 for a movie of the simulation.

Video 6
Simulated development of sbr.

We interpret the sbr mutation as causing a delay in signaling driving the transition of S-iridophores from dense to loose S-iridophore. We model this by reducing the attempted rate of transitioning from dense to loose to a rate 40 × less than the rate of attempting loose to dense S-iridophore transition. Due to available data, we initialise the model for sbr at 7.5 mm SL to match that published regarding the real fish (Figure 7K,K’). At 9 mm (Figure 7L,L’) melanocytes begin to accumulate either side of the widened initial X0 interstripe. At 10.6 mm SL, we observe melanocytes that are associated with dense S-iridophores (white cells) and not just with loose S-iridophores (blue cells) as usually seen in WT at 10.6 mm ≈ stage SA (between stages SP and J+, Figure 5C–D). At 11.5 mm (Figure 7M,M’), melanocytes are organised into aggregates, approximately one stripe width in size, and only partially connected, thus forming a broken pseudo-stripe pattern.

The cho mutant

Homozygous cho mutant larvae lack the horizontal myoseptum (Svetic et al., 2007). As a result, dense S-iridophores are prevented from traveling via the horizontal myoseptum to generate the initial stripe of dense S-iridophores seen in WT at stage PB. Instead loose S-iridophores appear only later, at stage PR, uniformly across the domain. cho fish then proceed to develop a labyrinthine pigment pattern. Stripes and interstripes of normal width form in a parallel arrangement, but with orientation disrupted, with regions running vertically and horizontally and often strongly curved, sometimes branched and often interrupted (Figure 7F).

To model cho, we omitted the initial stripe of dense S-iridophores at the PB stage (Figure 7O,O’) and instead place 200 loose S-iridophores at random across the S-iridophore domain at stage PR (Figure 7P,P’). No other interactions were altered. At stage J+ (Figure 7Q,Q’), we see a pattern of normal width stripes and interstripes except with varying orientation, as seen in real cho fish. See Video 7 for a movie of the simulation.

Video 7
Simulated development of cho.
The seurat mutant

Homozygous seurat mutants develop fewer adult melanophores, thus forming irregular spots rather than stripes. This phenotype arises from lesions in the gene encoding Immunoglobulin superfamily member 11 (Igsf11) (Eom et al., 2012) which encodes a cell surface receptor (containing two immunoglobulin-like domains) that is expressed autonomously by the melanophore lineage. Igsfl1 promotes the migration and survival of these cells during adult stripe development as well as mediating adhesive interactions in vitro.

To model seurat, we reduced the rate at which melanocytes could differentiate to a twentieth of the usual rate. This was to reflect the inhibition of the migration of melanoblasts (precursors of melanophores) across the domain and increased the rate of attempted melanocyte death to one hundred times per day (usually once per day). No other interactions were altered. At stage J (Figure 7R,R’), we see a pattern of normal width stripes broken into spots with a reduced number of melanocytes. Melanocytes are associated with loose S-iridophores and xanthophores with dense S-iridophores, as seen in real seurat fish.

By modelling seurat and sbr, we can also make predictions about the phenotype of a double mutant seurat;sbr, shown in Appendix 5—figure 1. We predict that by stage J+, this mutant would be covered in dense S-iridophore and associated xanthophores, with a few melanocytes at the very dorsal and ventral region of the fish. We are not aware of a published description of the phenotype of this double mutant, so our prediction remains to be tested.

Regeneration experiments

In order to further validate our model, we test to see whether we observe similar behaviour when the cells are ablated and the pattern is left to regenerate. In 2007, Yamaguchi et al., 2007 ablated a rectangular window of pigment cells of adult zebrafish stripes and recorded the regeneration of pigment producing cells (Figure 8A–C). They found that after ablation, cells regenerated in a labyrinthine pattern. To model this ablation, we simulated WT development from stage PB to our latest simulation stage, J+ as seen in Figure 5D’. At stage J+, we simulate ablation by removing a square region of cells in the centre (horizontally) of three stripes and interstripes. We then observe the pattern regeneration 7, 14 and 21 days later in Figure 8C–E. At 14 days, we observe the production of irregular shaped spots of melanophores in the centre of the ablated region as seen in the ablated fish at day 7. At day 21, we observe a regeneration of the pattern where stripes are no longer oriented horizontally. In some regions, spots of melanophores are surrounded by xanthophores.

Stripe regeneration simulations.

(AC) Regeneration of new pigment producing cells 7, 14 and 21 days respectively after a small rectangular window of cells in the adult WT stripes are completely ablated (Patterson and Parichy, 2013). (DF) A representative simulation of the regeneration of an adult zebrafish 7, 14 and 21 days after a simulated ablation has occurred. Scale bar is 0.25 mm in all images. Experimental images (AC) are reproduced from Yamaguchi et al., 2007. Copyright (2007) National Academy of Sciences U.S.A.; these images not covered by the CC-BY 4.0 licence and further reproduction of this panel would need permission from the copyright holder.

© 2007 National Academy of Sciences USA. Experimental images (A–C) are reproduced from Yamaguchi et al., 2007, with permission of PNAS. These images not covered by the CC-BY 4.0 licence and further reproduction of this panel would need permission from the copyright holder.

In 2013, Patterson and Parichy, 2013 ablated a section of dense S-iridophores along the horizontal myoseptum using a S-iridophore-specific marker pnp4a:NTR+Mtz at the beginning of pattern metamorphosis (Figure 9A, stage PB). They then observed the subsequent pattern formation (Figure 9B). We simulate this by removing a section of dense S-iridophores from the horizontal myoseptum at stage PB (Figure 9C) and then simulating as normal. We observe the pattern at 10 days post ablation. Xanthophores are associated with the undamaged portion of the S-iridophore stripe and melanophores surround the damaged interstripe (Figure 9D). In both cases, our simulations closely approximate the published observations.

Stripe regeneration post S-iridophore ablation.

(A) S-iridophores are ablated using pnp4a:NTR+Mtz at stage PB (B) Stripe development in wild-type fish following S-iridophore ablation. Xanthophores are associated with S-iridophore interstripe. Melanophores aggregate around the new interstripe. (CD) A representative simulation of S-iridophore ablation and subsequent development. For clarity, dense S-iridophores are displayed in green in (C). Scale bar is 0.25 mm in all images. Experimental images (AB) are reproduced from Patterson and Parichy, 2013 and licensed under CC-BY 4.0 (http://creativecommons.org/licenses/by/4.0).

Quantitative analysis of simulations

In the next few sections (3.2.1)-(3.2.4), we test the consistency of our WT and mutant simulations, averaged over 100 simulations, with real published quantitative measures. We test four criteria using experimental data: (1) the number of melanophores in mutants at different stages with respect to WT at the same stage, (2) the average width of X0 interstripe for WT, rse, pfe, shd and nac at stage J+, (3) WT stripe straightness, (4) WT pigment cell density in stripes and interstripes.

Melanophore density across WT and mutant fish

Figure 10A is a comparison between the average number of melanocytes per ventral hemisegment for each mutant with the number of melanocytes at WT at stages PB, PR, SP, J and J+ using the data from Frohnhöfer et al., 2013. First, since we do not have simulated data for a whole hemi-segment we normalise our melanocyte numbers against WT numbers at each stage. We do this by, for each respective stage and each respective mutant, multiplying the number of simulated melanocytes at the given stage for the given mutant by the number of melanocytes at each stage in real WT fish and dividing by the number of melanocytes observed in our WT simulations at this stage. These comparisons are given in Figure 10A. We observe the same trends seen in real fish. Moreover, in all stages, except for stage SP, the simulated data falls within the error bars of the measured data in real fish. In particular, we observe that the number of melanocytes remains similar to WT in pfe until stage J+, similar to that in WT, whilst the number of melanocytes is significantly lower in shd and rse in comparison to WT just as in real fish.

Quantitative analyses.

(A) Number of melanocytes per ventral hemisegment. Quantification of melanocytes. A comparison of melanocyte numbers between mutants and WT seen in our simulations and data recreated from Frohnhöfer et al., 2013. For each stage, our simulated number of melanocytes was normalised by the WT number in data by Frohnhöfer et al., 2013. (The number of animals used for counting was at least five for each measurement point.) The number of simulations was 100. Error bars indicate the standard deviation. There is no error bar for WT, since WT data was used to normalise the data. (B) X0 interstripe width. A comparison of X0 stripe width in 100 simulations with X0 stripe width in real mutants (taken from representative images; Frohnhöfer et al., 2013). The darker bars depict experimental data, lighter bars depict our simulated data. Error bars indicate the standard deviation. Since there is only one image for each of the respective mutants, there is no error bar for actual (real experimental data). (CF) PCF for different cell types averaged over 100 simulations using the square uniform PCF (Gavagnin et al., 2018). The spatial pair correlation as a function of distance for (C) dense S-iridophores only (D) melanocytes only, (E) melanocytes and dense S-iridophores, (F) melanocytes and xanthophores. Experimental data in (A) is reproduced from Frohnhöfer et al., 2013.

WT stripe straightness

In real life, zebrafish stripes are quite straight, but are not necessarily perfectly so (see X0 in Figure 1A for example). To measure stripe straightness, we first generated a line representation (x) of the central interstripe (see Appendix 3). From this line x, we calculated the stripe straightness SS(x), measured by the ratio of the length of our line (L) to the straight line distance between the ends of it (C), i.e.

(1) SS(x)=CL.

The value of SS(x) lies between 0 and 1, since CL. For more information about how SS(x) is generated see Appendix 3. In 10 real WT stage J+ fish, the mean SS value was 0.98. In 100 simulations we observe a high albeit slightly lower mean SS value of 0.92 at stage J+, demonstrating good stripe straightness, that is close to that observed in real fish.

X0 interstripe width across WT and mutant fish

Finally, we compare the interstripe X0 width in our simulations with real data. We choose this interstripe as it is the only interstripe that the mutants we consider and WT have in common. We compare the width of this interstripe at J+ in our simulations (see Appendix 3 for detailed methods) with the observations of the corresponding J+ mutant in Frohnhöfer et al., 2013. We demonstrate in Figure 5B that, in all cases, the real J+ mutant X0 interstripe widths fall in the range of ±1 standard deviations of the simulated J+ X0 stripe width averaged over 100 simulations. This demonstrates consistency between our model and real data. Thus, our model demonstrates an excellent ability to quantitatively simulate the patterns of real fish.

Pigment cell density in WT

There are no published estimates of WT pigment cell density in each of the stripes at the J+ stage when our simulations end. However, our data are comparable to those of adult WT fish measured by Mahalwar et al., 2016, who observed that in the stripe regions there were approximately four times more xanthophores in the interstripe region than melanocytes in the stripe region. Furthermore, whilst the light interstripe were completely devoid of melanocytes, there was a low density of xanthophores in the stripe region. In our model simulations, we observe a mean of 4.01 times as many xanthophores in the interstripes than melanocytes in the stripes demonstrating good agreement. We also observe a low density of xanthophores in the stripe regions and negligible numbers of melanocytes in the interstripe regions.

Simulation reproducibility of pattern formation

To further test the accuracy of our model’s outputs, we compare the spatial correlation of different cell types at different distances. We use this measure as an objective test of whether the spatial distributions between cells we observe in our representative simulations, (i.e. the patterns generated), are consistent among different simulated outputs.

To measure spatial correlation, we use a pair correlation function (PCF). A PCF determines whether, given a spatial distribution of agents on a domain, the number of pairs of agents at a certain distance from each other are greater than or fewer than the number expected if the agents were distributed uniformly at random. For example, if the PCF value is unity for a certain distance, this indicates that there is no spatial correlation. If the PCF value is greater or less than unity for a certain distance then this indicates that there is spatial correlation or anticorrelation respectively at that distance. The PCF we employ is specific to on-lattice domains and is called the Square Uniform PCF (Gavagnin et al., 2018) adapted for multiple cell-types (see Appendix 3 and Dini et al., 2018). We describe the PCF as homotypic when we are measuring the spatial correlation of one cell-type and heterotypic when we consider the spatial correlation between two different cell-types. We choose this PCF as it uses a measure of distance that is complementary to the distance measurement in our simulations.

Figure 10C–F shows the square uniform PCF against distance for different mutants and cell types averaged over one hundred simulations. For each plot of a given PCF type (homotypic melanocytes, for example), we repeatedly simulate the relevant mutant to its final stage, compute the PCF of the resultant pattern and then average the PCFs over the number of repeats to give us an averaged PCF.

To interpret the data, consider the representative simulation for a WT fish at stage J+. In this example, we observe stripes (Figure 5D’). This is a consistent feature for all of the repeat simulations of WT at stage J+. To quantify average interstripe width (the distance vertically in mm from the top of any interstripe to the bottom) in our simulations we can consider the averaged homotypic dense S-iridophores PCF for WT in Figure 10C. We observe that this shows periodicity (sequential peaks and troughs) at different distances. These are a consequence of the striped pattern at J+ (Figure 5D’). Since dense S-iridophores occupy the interstripe regions in WT and not the stripe regions at J+, dense S-iridophores are spatially correlated at short distances, indicated by a positive value of the PCF at short distances. Conversely, they are anti-correlated at distances approximately one half, one and a half or two and a half stripe widths away, as these distances correspond to the relative positions of the dark stripes, which normally lack dense S-iridophores. We see troughs at these distances. The period of the PCF in this case thus quantifies an estimate for average interstripe width.

In the next few paragraphs, we test the reproducibility of different features that are observed in our representative simulations by considering a PCF of appropriate cell types averaged over one hundred repeats.

Real cho mutants and WT fish share similar interstripe width (Frohnhöfer et al., 2013)(also seen in our model; compare Figure 7Q’ and Figure 5D’). To test reproducibility we consider the homotypic PCF of dense S-iridophores for WT and cho in Figure 10C. For both cho and WT simulations, the averaged homotypic PCF for dense S-iridophores observes periodic behaviour with the same frequency, indicating maintenance of interstripe width between WT and cho in our model, consistent with observations in vivo.

In real shd at stage J+, there are two pseudo-stripes of melanocytes broken into spots, of a diameter that is approximately equal to the normal stripe width (Frohnhöfer et al., 2013) (also seen in our model; compare Figure 5H’ and Figure 5D’). To test whether this is consistent we consider the homotypic PCF of melanocytes for WT and shd in Figure 10D. The average stripe width of WT and the average aggregate size of shd can be approximated from the PCF as the distance related to the first trough, as this is the shortest distance at which melanocytes are most anticorrelated with other melanocytes. For both shd mutants and WT simulations, these are both approximately 0.3 mm.

In real pfe at stage J+, stripes and interstripes remain aligned and have the same width as in WT, except that stripes are broken into spots and some melanocytes lie ectopically in the usual interstripe region (Frohnhöfer et al., 2013) (also seen in our model; compare Figure 5D’ and Figure 5P’). To test reproducibility, we consider the heterotypic PCF of melanocytes and dense S-iridophores for WT and pfe in Figure 10E. For both the WT and pfe simulations, the averaged heterotypic PCF of melanocytes and dense S-iridophores displays periodic behaviour with the same period. However, in pfe the peaks and troughs are damped. We interpret this as follows. Firstly, this indicates that, in our model, stripe width is preserved between pfe and WT as the period of the PCF is the same. Moreover, as the peaks and troughs are damped in pfe, this indicates that, as seen in our representative simulation, some melanocytes tend to remain in the interstripe regions.

In real sbr at 11.5 mm SL, stripes are sometimes broken into spots of usual (vertical) width, but the overarching stripe pattern remains (Fadeev et al., 2015) (also seen in our model, compare Figure 7D’ and Figure 5N’). To test reproducibility, we consider the heterotypic PCF of melanocytes and xanthophores for WT and sbr in Figure 10F. The first peak of this PCF corresponds to the shortest distance at which melanocytes and xanthophores are most correlated, which is approximately the stripe width. For both sbr mutants and WT simulations these are both approximately 0.3 mm.

In these examples, using appropriate PCFs we have demonstrated the consistency of our simulations in generating patterns that match the qualitative differences we expect when we compare mutant fish with WT. We note that we have only provided the averaged PCF for the scenarios aforementioned for simplicity. For information about how the PCF is calculated, please see Appendix 3.

Iridophore assumptions and biological redundancy

Necessity of S-iridophore assumptions

For the less well-studied S-iridophore transitions, we analysed key mutant phenotypes to infer biologically realistic rules for these interactions, aiming to generate assumptions that were the simplest for pattern formation changes seen, but no simpler. These deductions are discussed in Section "Modelling overview". In Figure 11 B1–J3, we demonstrate the necessity of all of the assumptions for dense-to-loose, loose-to-dense S-iridophore transitions first outlined in Figure 4G, 15-18 for stripe formation by showing representative images (the model was run 50 times each) of simulations at stage J+ for fish lacking one or more S-iridophore transition mechanisms display patterns which diverge from those seen in real fish.

Representative simulations of the model with some of the S-iridophore interactions removed.

The first column displays a diagram of the S-iridophore interactions which remain (all other cell-cell interactions are unchanged). Columns 2–4 are representative simulations of WT, pfe and nac under these conditions. *This interaction is equivalent to removal of long range xanthophore inhibition in criterion 17. It is also equivalent to removal of short range xanthophore promotion in criterion 18.

First, we analyse simulations lacking one of the transition types loose-to-dense or dense-to-loose (Figure 11 B1–C3). Without a loose-to-dense transition (Figure 11 B1–B3), note how in all cases (WT, pfe, nac) only one pseudo-interstripe is preserved: the initial X0 interstripe. The X0 interstripe is surrounded by a sea of loose S-iridophores which have transitioned to loose either by promotion by xanthophores in the long range (nac Figure 11 B3), promotion by melanophores in the short range (pfe Figure 11 B2) or a combination of both (WT Figure 11 B1). This suggests that loose-to-dense interactions are important for generating subsequent interstripes. Interestingly, without loose-to-dense transitions WT fish demonstrate a striped pattern similar to Danio albolineatus (Parichy, 2006) suggesting a possible route of evolution between these fishes (also noted by Volkening and Sandstede, 2018). Without a dense-to-loose transition (Figure 11 C1–C3), the S-iridophores form a dense sheet over the entire domain with xanthophores and melanophores scattered across the domain. This demonstrates the necessity for S-iridophores to be able to transition between dense and loose form.

Next, we consider removing each of the criteria required for an S-iridophore transition one at a time (Figure 11 D1–J3). We first note that, in some cases, removal of an interaction in either nac or pfe results in loss of a transition type. These are not shown in Figure 11 for simplicity. In other scenarios, however, removal of an interaction leads to the uninhibited possibility of a transition in one of the single cell mutants. For example, consider Figure 11 E2. Removal of long-range xanthophore promotion in criterion 16, leads to the possibility of a transition from loose to dense provided that there are either melanophores in the short range or no xanthophores in the short range. Since in pfe there are are always no xanthophores in the short range, or indeed anywhere on the domain, S-iridophores are consistently promoted to dense type, with a non-zero rate, thus Figure 11 E2 is not distinguishable from Figure 11 C2 since in effect, in pfe the same interactions have been knocked out. We note that this is the case for one of either nac or pfe in Figure 11C,E–G,J.

In Figure 11 D1–D3, short-range xanthophore inhibition is removed from criteria 16. As this is exclusively a xanthophore-iridophore interaction this only effects WT and nac. Without the promotion of xanthophores in the short range, S-iridophores change from dense to loose in the interstripes, making the interstripes appear faded. In Figure 11 G1–G3, short-range melanophore inhibition is removed from criteria 18. Therefore S-iridophore can change from loose to dense when there are xanthophores in the short range. As a result interstripes become wider as they are unrestricted by local melanophore stripes. Finally, in Figure 11 I1–I3 criteria 18 is removed, so S-iridophores only change from loose to dense when there are no melanophores in the short range and simultaneously no xanthophores in the long range. In this case stripe integrity is lost in WT.

To summarise, all interactions are necessary for pattern formation in WT and single cell mutants, nac and pfe.

Biological redundancy

As part of this in-depth study, we have incorporated all of the interactions that we have identified from the literature. Consequently, there may be some in-built redundancy. However, we keep all interactions for the purposes of biological realism. In this section we explore the idea of biological redundancy by removing some interactions and observing the resultant simulated development.

First we consider movement. In real (and simulated) fish, melanocytes and xanthophores move 0.11 mm per day (Takahashi and Kondo, 2008) and 0.033 mm per day, respectively. Furthermore, in the short range their direction of movement is influenced by each other (Yamanaka and Kondo, 2014). Xanthophores chase melanocytes, which in turn, are repelled by xanthophores. In Figure 12C–D, we turn off the movement of xanthophores and melanocytes in WT and shd mutants and simulate to stage J+. We observe that in both cases, the pattern is conserved. This is not surprising since the cells do not move very quickly. However, there is a slight difference in shd wherein the interstripe width is slightly smaller than in shd without changes (Figure 12B). We suggest that this is due to the loss of chase-run dynamics observed between melanocytes and xanthophores. Next we consider the differentiation rate of melanocytes. We model this as being dependent on the number of dense S-iridophores and the number of xanthophores currently on the domain. This is because we assume that dense S-iridophores and xanthophores positively influence the rate of melanocyte birth in the long range. In Figure 12E, we change the differentiation rate so it is only dependent on the number of xanthophores currently on the domain and not the number of S-iridophores, effectively reducing the rate of differentiation of melanocytes in wild-type fish. The resultant pattern is still striped; however, it is less organised. In Figure 12F–G, we remove the mechanism allowing long-range communication between xanthoblasts and melanocytes. In shd our model predicts that without the consolidation of spots by xanthoblasts, melanocyte spots become more widely spaced. In Figure 12H, we remove xanthophore and xanthoblast proliferation. This limits the number of xanthophores to the number allocated at the start. Remarkably, our model predicts that stripe formation is largely preserved, however, interstripes are fainter due to the lack of xanthophores. Next we consider melanocyte death. In Figure 12I–J, we remove melanocyte death from WT and shd simulations. Whilst stripe formation is maintained in WT (Figure 12I), interstripes are littered with melanocytes. In shd (Figure 12J), melanocyte spots are more difficult to discern as melanocytes can be observed in the xanthophore regions. We predict that dense S-iridophores and xanthophores in the long-range and/or loose S-iridophores and the lack of xanthophores in the short range promote melanocyte differentiation. In Figure 12K we change the criteria for melanocyte differentiation, so that it does not depend on S-iridophores (only xanthophores). As a result the stripe pattern looses integrity. A similar phenotypic change happens when we change the criteria for melanocyte differentiation, so that it does not depend on xanthophores (only S-iridophores, Figure 12N). In Figure 12L,M, we change melanocyte and xanthophore movement so that they are no longer influenced by each other (no chase-run dynamics). As a result in WT (Figure 12L), stripe integrity is lost. For shd (Figure 12M) simulations, however, qualitatively there is not much difference (similarly to the case when movement is completely removed (Figure 12D) suggesting that movement does not play a significant role in generating spots. Rather, it is the death of melanocytes in xanthophore rich areas and xanthoblast pulling which are important (Figure 12G,J).

Representative simulation of the model with certain rules omitted.

(AB) WT and shd respectively with all rules included. (CD) No movement of melanocytes or xanthophores for WT and shd, respectively. (E) Differentiation of melanocytes only promoted by dense S-iridophores for WT fish that is the rate of melanocytes does not depend on the number of S-iridophores. (FG) No pull of melanocytes by xanthoblasts for WT and shd fish, respectively. (H) No proliferation of xanthophores for WT. (IJ) No melanocyte death for WT and shd fish. (K) Differentiation of melanocytes is not dependent on S-iridophores, (LM) No death of melanocytes for WT and shd (N) Differentiation of melanocytes is dependent on dense S-iridophores being present in the long distance only.

In summary, whilst removal of these interactions largely do not change the type of pattern generated (e.g. shd still generates spots, wild-type fish still generate stripes), and thus could be considered as biologically redundant for pattern generation, they appear to have large impacts on the integrity of the patterns formed. We thus, suggest that the retention of these interactions in vivo act as a buffer to protect the integrity of stripe formation in spite of stochastic variations in stripe patterning.

Model predictions

A major benefit of developing a fine-grained, cell-level model is the ability to perform in silico experiments that can be directly related to real-life equivalents. This not only allows us to explore parts of the system that may otherwise not be testable experimentally, giving us valuable insight into the biological system. It also gives us the ability to analyse dynamics of different hypothetical mechanisms before devoting expensive resources to experimental tests that could confirm theoretical findings.

In the next few sections, we focus on the ability of the model to make biologically testable predictions, demonstrating firstly, in Section "An in silico investigation into important mechanisms for controlling pattern formation", how we can use our model to explore important facets of successful pattern formation such as growth, domain size and initial conditions. Then in Section "An in silico investigation into the function of the leo gene", we give an example of how we can use our model to generate testable hypotheses about the leopard mutant.

An in silico investigation into important mechanisms for controlling pattern formation

Initial S-iridophore interstripe orientation alone does not determine the orientation of stripes and interstripes

Previously it has been hypothesised that the horizontal orientation of the initial S-iridophore interstripe (emerging from the horizontal myoseptum) that drives the organisation of subsequent stripes and interstripes horizontally. One way to test this hypothesis is to initialise the interstripe so it is oriented vertically instead of horizontally. If the initial S-iridophore interstripe does orient stripes and interstripes then we would expect to see the same pattern development we observe in WT fish, but rotated 90 degrees. That is, we would expect to see vertical bars across the domain at the time corresponding to stage J+. The position of the dense S-iridophores (a horizontal interstripe along the horizontal myoseptum in WT fish) at the start of pattern metamorphosis, is dictated by the fish’s anatomy and cannot be altered experimentally. However, we can simulate an altered iridophore initial distribution in silico by initialising the initial interstripe as a band of width three along the vertical axis instead of as a band of dense iridophores vertically (dorso-ventrally) down the centre of the domain of width three instead of as a band of width three, dense S-iridophores along the horizontal axis. We observe the subsequent pattern development at stages PR, SP and J+ in Figure 13A. Interestingly, instead of observing vertical bars, at stage J+ we observe a labyrinthine pattern. This demonstrates that, whilst the initial S-iridophore interstripe plays a role in orientating the pattern, it is not the only part of the initial condition that is important. Further observation of the model output reveals that, for a while, the pattern is oriented in a vertical pattern similar to the initial iridophore interstripe, but as growth continues, this pattern becomes reoriented into a horizontal form. This clearly reveals the impact of growth on pattern formation.

In silico investigation into important mechanisms for controlling pattern formation.

The simulated domains at stages PR, SP and J+ wherein the following are changed (A) the orientation of the initial S-iridophore interstripe, (B) the initial position of the S-iridophore interstripe, (C) the initial domain size (D) the initial domain so that it is populated with adult-width stripes and interstripes.

The position of the initial S-iridophore interstripe is important for successful pattern formation

Another way to better understand the role of the initial S-iridophore interstripe is to alter its position in the dorso-ventral axis so that it appears more ventrally on the initial domain. We might naively predict, given that it has been hypothesised that dense S-iridophores only orientate the stripes and interstripes, that in this case the final pattern would be the same as in WT fish. We simulate this in silico by initialising the initial interstripe to be one quarter of the way up the domain instead of half way. A typical pattern evolution for this initial condition is displayed in Figure 13B. We observe that subsequent stripes and interstripes still appear sequentially, either side of the initial interstripe, suggesting that the S-iridophores do play a role in the positioning of new stripes and interstripes. However, we do not observe usual WT patterning. In particular, stripes and interstripes exhibit more breaks compared to WT simulations. Moreover, developing stripes and interstripes become sequentially thinner as a result of the impact of domain growth. Once again, growth is the key factor: growth is centred at the middle of the domain and so when the initial stripe is not similarly centred, growth disrupts pattern formation in our model.

Initial domain size contributes to the number of stripes and interstripes

In order to test the role of domain size in pattern development, we initialise the domain so that it is three times as tall as in WT simulations. That is, we initialise the domain to be 2 mm × 3 mm instead of 2 mm × 1 mm. All other parameters remain unchanged, including the rate of growth in the horizontal and vertical axis. We present a typical example of subsequent pattern development in Figure 13C. At stage PR, pattern development is similar to the pattern seen in WT fish at the same stage. However, at stage SP, we observe more stripes than are observed even by stage J+ in WT fish. By stage J+, instead of three interstripes and four stripes as seen in WT, we observe six stripes and seven interstripes. This suggests that the initial domain height influences the number of stripes and interstripes that develop, provided that growth is uniform and centred.

Stripe insertion can occur on an initially striped domain

Kondo and Asai observed that, as the size of the marine angelfish Pomocanthus doubled, new stripes along the skin would develop between the old ones (Kondo and Asal, 1995). This phenomenon has not been observed in zebrafish, where new stripes and interstripes appear consecutively at the dorsal and ventral periphery. We hypothesise that this is likely related to either pattern maintenance mechanisms or the spatial localisation of growth. Here, we experiment with the model to see whether stripe insertion can occur when the domain is populated at stage PB with adult-width stripes and interstripes. The results of an example realisation with these initial conditions are given in Figure 13(D). We observe that, in this case, new interstripes do appear between pre-established stripes. This is because growth (which is centred in the middle of the domain) creates space within the middle of the already developed stripes and interstripes.

S-iridophores are more important to the generation of melanocytes than xanthophores

We also used the model to make some more subtle predictions. For example, in the case of melanocyte differentiation, which we model as being promoted in the long range by both xanthophores (from observation of ablation experiments; Kondo and Asal, 1995) and S-iridophores (from observations of pfe), there were no known parameter values for their relative strengths. We found using our model that by making the strength of S-iridophore promotion of melanocyte differentiation to be much greater than that of xanthophores, qualitatively and quantitatively the model simulated for WT, pfe and shd was greatly improved (see Figure 14A–B).

Substantial alterations in starting conditions generate major alterations in pattern formation.

The effects of different signalling strengths in pfe and shd. Representative images of J+ simulations of (A) shd and (B) pfe where the signalling strength of melanocyte differentiation in the long range is increased (from left to right). Increasing the signalling strength of the respective cell types decreases the width of the interstripe X0 and increases the number of pseudo-stripes from 2 to 4. For our simulations, we choose the weakest shown signalling strength for xanthophores, corresponding to the furthest left of (A) and the strongest shown signalling strength for S-iridophores, corresponding to the furthest right of (B). (C) Stripe straightness of our simulations correlates with the simulated HAA/SL ratio. Each circle is one simulation, green circles indicate simulations up to stage J, black circles indicate simulations of up to stage J+. There are one hundred simulations in total for each stage. The red line is the line of best fit.

Horizontal growth bias during development generates more tortuous stripes in WT fish

Interestingly, we also observed in our simulations that increased height-to-length ratio is correlated with stripes becoming more tortuous (R = −0.617, p<0.01, Figure 14C). This phenomenon is not something we can see as being consistent with real fish and thus suggests that some interactions may be missing regarding the maintenance of stripe and interstripe formation.

An in silico investigation into the function of the leo gene

The model provides testable hypotheses for cryptic functions of the leo gene

The gene leo encodes Connexin39.4 (Cx39.4) (Watanabe et al., 2006; Maderspacher and Nüsslein-Volhard, 2003; Irion et al., 2014). As a result, leo mutants display a leopard-like spotted pattern across the flank of the fish (Figure 15A–A’), instead of the usual striped pattern (Figure 1A). In this section, we aim to hypothesise key aspects of the leo mutations using our model alongside observations of relevant mutants.

Using the model to generate predictions regarding the underlying effects of the leo mutant.

(AD) Adult leo, nac, pfe and shd respectively. (A’D’) The flanks of leo, nac, pfe and shd respectively. (B’’D’’) The flanks of leo;nac, leo;pfe and leo;shd respectively. (E) Example simulations where combinations of hypotheses 1–3 in the text are implemented. (F) Example simulations where combinations of hypotheses 1–5 in the text are implemented. Note that simulations of 1–5 are similar to simulations for 2–5, indicating that hypotheses one is not necessary to predict the leo phenotype. (G) Adult leo;sbr (G’) Flank of leo;sbr. (HL) Adult mutants contrasting with (H’L’) simulations using hypotheses 1–5 from the main text. (HL) Flank of leo, leo;nac, leo;pfe, leo;shd and leo;sbr respectively. (H’L’) Simulated images of leo, leo;pfe, leo;nac, leo;shd and leo;sbr respectively at stage J+. (A), (A’), (G), (G’), (H) and (L) from Fadeev et al., 2015, (BD’’), (IK) from Irion et al., 2014 and are all licensed under CC-BY 4.0 (http://creativecommons.org/licenses/by/4.0).

Pattern formation is also altered in the double mutants leo;shd, leo;nac and leo;pfe when compared with shd, nac and pfe (Irion et al., 2014). For example, the flank of double mutant leo;nac is covered by xanthophores and dense S-iridophores (Figure 15B’’). This is in contrast to nac which also contains large patches of loose S-iridophores (Figure 15B–B’). Adult leo;pfe fish exhibit randomly distributed melanocytes instead of spots (Figure 15C–C’’). Finally, adult leo;shd exhibit an absence of melanocytes on the flank of the fish, instead, the flank of the fish is entirely covered with xanthophores. This is contrast to the melanocyte spots normally observed on shd (Figure 15D–D’’).

Connexins are involved in cell–cell communication and signalling. Since, Cx39.4 is required for normal function in melanocytes and xanthophores but not in S-iridophores (Watanabe et al., 2006; Maderspacher and Nüsslein-Volhard, 2003; Irion et al., 2014), this suggests that in leo, cell–cell communication between melanocyte and xanthophores may be disrupted. Moreover, from observation of the double mutants, it seems that leo presumably generate heteromeric gap junctions among and between melanophores and xanthophores, controlling S-iridophore shape transitions (Irion et al., 2014).

In order to investigate the influence of the leo gene, we first consider the individual cell–cell interactions that can be deduced from the literature and ask if these cell–cell interactions are sufficient for generating the pattern. So far, there has been one experimental study observing the individual behaviour of leo cells. Kondo and Watanabe, 2015 studied the movement in-vitro of leo melanocyte and xanthophore cells. They demonstrated that the leo melanocyte repulsive response to xanthophores was hardly observed in comparison to the marked repulsion in WT fish. This suggests that melanocyte repulsion from xanthophores is inhibited in leo (Kondo and Watanabe, 2015). We will refer to this as hypothesis one for the effects of the mutant leo.

  • Hypothesis 1:Melanocytes are not repelled by xanthophores.

We can simulate pattern development in this case by turning off melanocyte repulsion by xanthophores in our model. This is shown in Figure 15E in the column numbered 1. We observe that at J+ the pattern consists of thicker interstripes than in WT fish, but not spots. Simulating the shd phenotype (lack of S-iridophores) with hypothesis 1, also does not generate the pattern expected in leo;shd. Hypothesis one is also insufficient to explain the phenotype of leo;nac or leo;pfe. This is because there are either no xanthophores or no melanocytes respectively in these mutants and therefore no melanocyte-xanthophore interactions. From these observations, we conclude that hypothesis one alone is insufficient for leo pattern formation.

To deduce the other cell-cell interactions that a mutation in leo could affect, we look at the patterns of adult leo;nac (Figure 15B’’) and leo;shd (Figure 15D’’). Of the three double mutants, these are the easiest from which to deduce single cell-cell interaction changes that could generate the double mutant patterns.

The leo;nac mutant displays an absence of loose S-iridophores compared to nac, suggesting that signalling of xanthophores which promote S-iridophore transition from dense to loose is inhibited in leo. The leo;shd mutant displays an absence of melanocytes in the adult pattern compared to shd, suggesting that long-range survival signal sent by xanthophores to melanocytes is inhibited in leo. We propose two further potential hypotheses for the effects of the mutant leo.

  • Hypothesis 2: Xanthophores do not promote the survival of melanocytes in the long range.

  • Hypothesis 3: Xanthophores do not promote the change of S-iridophores from dense to loose in the long range.

In Figure 15E, we provide a table of results for all combinations of hypotheses 1–3. We comment that hypotheses 1–3 cannot be all-encompassing, as so far, none of our hypotheses effect melanocyte-iridophore interactions and thus cannot generate the phenotype leo;pfe. Meanwhile, we focus on being able to generate leo;nac and leo;shd.

In Figure 15E, we demonstrate, that none of the hypotheses alone can generate the leo pattern, nor can they generate more than one of the double mutant types. When hypotheses 1 and 2 are combined, striping is disrupted, there are more breaks than in WT and interstripes are wider than WT. When 1 and 3 are combined, stripes are the same as in WT, presumably due to compensation of other cell-cell interactions, however, the simulated shd pattern has less aggregation of melanocytes than typically observed in shd. When hypotheses 2 and 3 are combined we observe the expansion of the initial interstripe to the very edges of the domain dorso-laterally, where there are some melanocytes, unlike leo. However, in this case, our simulations of leo;nac and leo;shd match the real phenotype. We demonstrate that when all hypotheses are implemented we generate a unstable pattern of spots which, eventually disappear over time, leaving a domain consisting of dense S-iridophores and xanthophores. Also, we can replicate leo;nac and leo;shd.

Finally, we attempt to elucidate the melanocyte-iridophore cell-cell interactions that are affected in the leo mutant by considering the phenotype of adult leo;pfe (Figure 15C’–C’’). The leo;pfe mutant displays a random distribution of melanocytes below a field of S-iridophores, suggesting that leo affects directed differentiation of melanocytes by dense S-iridophores. This leads us to two extra hypotheses for the effects of the leo gene:

  • Hypothesis 4: Melanocytes lose death signals from local dense S-iridophores and, as a result, can differentiate in dense S-iridophore zones.

  • Hypothesis 5: Melanocytes lose directed signalling from S-iridophores and hence, in the absence of xanthophores differentiate randomly.

In Figure 15F, we provide a table of representative simulated patterns for all combinations of hypotheses 1–3 with hypotheses 4 and 5. Any of the combinations considered can generate the leo;shd and leo;nac mutant phenotypes. Hypotheses 1–4 and 1–3,5 cannot generate the leo spots. However, if we combine all five hypotheses then we can replicate the phenotypes of all considered mutants; leo, leo;pfe, leo;nac and leo;shd. These results suggest that hypotheses 1–5 are sufficient for explaining leo.

Necessity of the pre-existing hypothesis about leo

Next, we evaluate the necessity of hypothesis (1) Hypothesis 1 is the only assumption that comes directly from the literature (Yamanaka and Kondo, 2014). We show, using our model, in Figure 15F that there is no notable difference between the patterns generated with and without hypothesis 1 (i.e. hypotheses 1–5 versus hypotheses 2–5) other than that the spots may be slightly more irregular in when hypothesis 1 is included. This suggests that the cell–cell interaction mediating repulsion of melanocytes by xanthophores in leo, may not be necessary to generate the characteristic spots.

Hypotheses 1–5 can replicate the patterns of leo, leo;nac, leo;pfe and leo;shd

Figure 15H—K display the flanks of leo, leo;pfe, leo;nac and leo;shd respectively. Figure 15H’–K’ display the simulated patterns at stage J+ for in silico mutants leo, leo;pfe, leo;nac and leo;shd respectively where our assumptions for the function of leo are based on hypotheses 1–5 and, in the case of double mutants leo;pfe, leo;nac and leo;shd, our assumptions for pfe, nac and shd are as previously described in Section "Simulation of ‘missing’ cell mutants".

We observe, for our simulations of leo (Figure 15H’), that melanocyte spots are associated with loose S-iridophores in a sea of dense S-iridophores and xanthophores just as in real leo fish (Figure 15H). For our simulations of leo;nac (Figure 15I’) we observe a domain that is fully populated with xanthophores and dense S-iridophores as expected (Figure 15I). Our simulations of leo;pfe (Figure 15J’) is fully populated with S-iridophores and randomly distributed melanocytes as expected (Figure 15J). Finally, our in silico representation of leo;shd (Figure 15K’) is fully populated with xanthophores only, as seen in Figure 15K.

The leo;sbr cross mutant phenotype is an emergent property of the model

Previously in Section "Simulation of other mutants" we considered the sbr mutant. The sbr gene encodes Tight Junction Protein 1a (Tjp-1a), which is expressed cell autonomously in dense S-iridophores and causes the shape transition from dense to loose to be delayed (Fadeev et al., 2015). We showed in Figure 7K’–N’ that our model could replicate the sbr pattern development by altering the rate at which S-iridophores attempt to change from loose to dense to be slower than in WT.

Figure 15G—G’ shows a typical example of the leo;sbr phenotype. The adult leo;sbr displays spots which, compared to leo (Figure 15A–A’), are more elongated. Figure 15L displays the flank of an adult leo;sbr. In Figure 15L’, we simulate a mutant that satisfies hypotheses 1–5 as well as our assumptions about sbr to generate an in silico leo;sbr. Consistent with real leo;sbr mutants, our simulation of leo;sbr has melanocyte spots associated with loose S-iridophores that are surrounded by xanthophores and dense S-iridophores. Comparison with our simulation of leo (Figure 15H’) demonstrates that our simulation of leo;sbr also displays more elongated spots than those of our leo simulations.

Robustness of the leo assumptions in generating spots

As a test of robustness, we perform a rigorous robustness analysis by carrying out one hundred repeats of the mutant simulations with perturbed parameter values chosen uniformly at random from the range 0.75–1.25 of their described value as in Section "Simulation of WT pattern". Ten of these randomly sampled repeats are given in Figure Appendix 4—figure 1 for leo. We observe that for all one hundred repeats that small perturbations to the rates still generate consistent spots, demonstrating the robustness of the model.

These results of this section demonstrate the remarkable ability of our model to generate the leo single and double mutant phenotypes under a set of specific proposed changes to the model rules; these proposals can be used to guide experimental exploration of the effects of the leo gene.

Discussion

As a result primarily of beautiful experimental work, zebrafish pigment pattern formation has become the best characterised pigment patterning mechanism (Frohnhöfer et al., 2013; Kelsh et al., 2009; Singh and Nüsslein-Volhard, 2015; Watanabe and Kondo, 2015; Patterson and Parichy, 2019). Zebrafish pigment pattern formation requires three pigment producing cell types: melanocytes, xanthophores and S-iridophores to generate the WT stripes (Frohnhöfer et al., 2013). In an attempt to decipher the mechanisms underlying pattern formation, most previous mathematical models have largely focussed on xanthophores and melanocytes and neglected S-iridophores (Nakamasu et al., 2009; Bullara and De Decker, 2015; Volkening and Sandstede, 2015; Painter et al., 2015; Bloomfield et al., 2011). However, recent studies have shown that S-iridophores are a crucial component, with S-iridophore transitions between dense and loose driving pattern formation. First, the early provision of dense S-iridophores through the horizontal myoseptum orients stripes (Frohnhöfer et al., 2013; Kelsh et al., 2009; Singh and Nüsslein-Volhard, 2015; Watanabe and Kondo, 2015; Patterson and Parichy, 2019). Moreover, S-iridophores sequentially pre-determine stripes and interstripes via respective shape changes in these regions (Frohnhöfer et al., 2013; Fadeev et al., 2015; Krauss et al., 2014).

Here, we have taken a bottom up approach to modelling zebrafish pattern formation with the aim of testing whether the experimentally defined set of biological rules for zebrafish pigment pattern formation might be sufficient to explain both the WT and the diversity of mutant pigment patterns. We used an individual based modelling approach incorporating all five cell-types deemed important for pattern formation in zebrafish. We formalised all respective cell-cell interactions mathematically, with interaction strengths, parametrised, where possible, by the biological literature (see Supplementary file 5). For the less well-studied S-iridophore transitions, we analysed key mutant phenotypes to infer biologically realistic rules for these interactions, aiming to generate assumptions that were the simplest for pattern formation changes seen, but no simpler. We proved our models ability to simulate the distinctive pattern features during developmental stages PB through J+ of each of WT and six mutant patterns that had been used to determine the biological rules. We showed that in each case, our model simulations matched qualitatively the pattern development in real fish at the various developmental stages considered. This is consistent with the proposal that our modelling assumptions were sufficient for pattern development in these cases. As a more rigorous test of the model, we then investigated its ability to successfully simulate the distinctive patterns of three further mutants with defective S-iridophore properties, including two mutants that had not been modelled mathematically before. We showed that in each case, our model also correctly replicated patterns that were qualitatively similar to the corresponding mutant fish at various developmental milestones.

We assessed multiple quantitative features of our simulations against measured data from published studies, focusing on spatial distributions of cell-types, stripe width and melanocyte numbers. We found that in each case our simulations were highly reproducible, and quantitatively matched the biological observations. We conclude that our mathematical modelling approach, built upon the biological literature, provides substantial validation of the sufficiency of that set of biological rules in explaining pattern formation in zebrafish development for WT and many other mutant fish. Furthermore our modelling provides support for the plausibility of the deduced rules for S-iridophore packing transitions during pattern development.

Finally, we demonstrated the capability of our model to give valuable insight into the patterning mechanism and to make testable predictions about the biology.

This paper represents the first demonstration, to the best of our knowledge, of a model being used explicitly to test the impact of the sbr mutation. The sbr gene encodes Tjp1a, a key tight junction protein, and is expressed at much higher levels in S-iridophores in a dense configuration than those of the loose form (Fadeev et al., 2015). Furthermore, double mutant and chimaeric studies show that sbr acts cell-autonomously within the S-iridophores to control adult pigment pattern formation (Fadeev et al., 2015). These authors also show that in sbr mutants the transition from dense to loose S-iridophores is delayed, suggesting that this transition somehow depends upon Tjp1a in dense S-iridophores (Fadeev et al., 2015). Here, we test the patterning impact of this interpretation, by incorporating delayed S-iridophore state transition into our model, and show that this does indeed result in pattern changes consistent with the sbr phenotype. This provides theoretical support for Fadeev and colleagues deductions and deepens the interest in understanding the mechanistic basis for this role for Tjp1a.

Our modelling results demonstrate the applicability of complex models to test hypotheses that are difficult to test experimentally. Previously, it has been hypothesised that S-iridophores contribute to pattern formation by orienting the stripe. In Section "An in silico investigation into important mechanisms for controlling pattern formation", we demonstrate that this is true. We show that simply reorientating the interstripe to a vertical position is not enough to produce vertical bars as our simulations exhibit a labyrinthine pattern instead of vertical bars. Careful observations of our simulations indicate that in addition to the initial condition, growth is important for determining the final pattern. We show that moving the initial position of the interstripe away from the centre of the flank, where growth is centred, also disrupts the patterning. We are also able to show that by enlarging the initial domain so that it is the same width but taller in height we show that we can generate more stripes and interstripes than that are usually observed at stage J+. Finally, by initialising a domain that is already populated with cells in a stripe position that we can replicate a different stripe formation mechanism seen in fish Pomocanthus. Whilst in zebrafish stripe formation is sequential, starting from an initial interstripe and developing bidirectionally from the middle, Pomocanthus develops stripes in-between other stripes. We show that if a striped pattern is fully formed when the fish is still growing (and growth is centred) that this forces new stripes to occur between the old ones, instead of at the periphery.

Our model is the first, to the best of our knowledge, to suggest that S-iridophores play a more important role in melanocyte differentiation than xanthophores. We found that by implementing a stronger signaling capacity of S-iridophores than xanthophores for long range melanocyte differentiation then we obtained a better qualitative and quantitative match for shd and pfe mutant patterns. In particular, a stronger signalling success rate for melanocyte differentiation from S-iridophores in the long range appeared to align subsequent stripes in pfe and WT, resulting in consistency of pseudo-stripe and stripe width respectively. The comparatively reduced signalling success rate for melanocyte differentiation from xanthophores in the long range reduced the number of pseudo-stripes in shd from four, to two at stage J+ which is more consistent with real data. Real shd fish typically exhibit fewer stripes (approx two at stage J+) than the four of WT (Frohnhöfer et al., 2013Figure 4A–B). We also found that this factor was important to produce melanocyte numbers that quantitatively matched data by Frohnhöfer et al when comparing between mutant and WT (Figure 10A) as well as a better stripe width match (Figure 10B). Thus, our study further reinforces findings that S-iridophores play an important role in determining stripe and interstripe width (without S-iridophores, X0 interstripe width in shd is increased) and not just the widely reported role of stripe alignment.

We have further built upon previous mathematical modelling work, by using our model to make predictions about the functions of the leo gene. The gene leo encodes Connexin 39.4, which is required in melanocytes and xanthophores but not S-iridophores (Irion et al., 2014). Connexins play an important role in cell–cell signalling and communication so it has been postulated that the leo gene is involved in the signalling between melanocytes and xanthophores as well as signalling cues to S-iridophores regarding shape-transitions (Irion et al., 2014). Previous to our investigation, there had been one study of the individual cell-cell interactions in leo. This study by demonstrated that unlike WT melanocytes, leo melanocytes are not repelled by xanthophores in the short range (Yamanaka and Kondo, 2014). Using our model, we first demonstrated that this cell-cell interaction alone was not enough to reproduce the leo mutant pattern. Then, in a systematic approach we deduced four hypotheses about cell–cell interactions that might be affected by leo, which, upon implementing in our model, successfully replicated the patterns observed in leo, leo;pfe, leo;nac, leo;shd and leo;sbr. This work provides testable hypotheses about the effect of leo which can now guide future experimental work.

In contrast to most previous mathematical studies of pattern formation, the rules we propose for zebrafish pigment patterns are complex and extensive. For example, successful S-iridophore shape transitions in our model require information from both melanocytes and xanthophores. Many other studies have condensed zebrafish pattern formation to a few simple rules that can often be described by a series of partial differential equations (Kondo, 2017; Painter et al., 2015; Nakamasu et al., 2009), in particular Turing reaction-diffusion-type models posit that combinations of short and long range dynamics between melanocytes and xanthophores generate stripe patterns. Indeed, a lot of the excitement around such models is the ease with which small parameter value changes can sometimes result in diverse patterns, many readily recognisable from nature (Watanabe and Kondo, 2015; Metz et al., 2011; Maini et al., 2012). A major difference between our model and Turing reaction-diffusion models is that small parameter changes in our model do not typically generate qualitatively different patterns, whereas Turing reaction-diffusion models can show substantial pattern changes in response to small alterations to parameters near to bifurcation points (Watanabe and Kondo, 2015; Maini et al., 2012). We suggest that the added complexity of the real system has evolved to make the patterning process robust, with partially redundant mechanisms insulating against the impact of stochastic variation during pattern formation.

Our modelling approach is analogous to that adopted in another recent study of zebrafish pattern formation, one which independently attempts to understand S-iridophore contributions to the process (Volkening and Sandstede, 2018). In a similar fashion, Volkening et al generated an off-lattice model incorporating S-iridophores for which the rules were based upon the experimental literature. Upon implementing these rules they attempted to test the sufficiency of the known biology in describing the patterning process. Importantly, Volkening et al’s modelling approach differs from that adopted here in several ways. First, we use an on-lattice model, whilst Volkening et al use an off-lattice model. Using an on-lattice model allowed us to incorporate volume exclusion by other cells directly. Secondly, we used a continuous-time model, whereas Volkening et al update their model at simulated 24 hr intervals, with all rules implemented simultaneously. By using a continuous-time method, we are able to capture the stochasticity involved in rates of reactions and the ordering of events over time. Finally, we incorporate a hypothesis that S-iridophores contribute more strongly to promoting melanocyte differentiation than xanthophore differentiation, leading us to a better qualitative approximation of shd mutants and a better quantitative estimate of melanocyte numbers in shd, pfe and rse than shown before.

Importantly, the model of Volkening et al. also proves highly capable in generating simulations that accurately mimic WT and various mutant pigment patterns. This observation further strengthens support for the validity of the proposed S-iridophore rules we each postulate. More broadly, we consider that our independent mathematical approaches are mutually reinforcing in reaching the conclusion that the deduced biological rules may be largely sufficient to explain pigment pattern formation in zebrafish.

Further testing of our model should focus on investigating later stages of pigment pattern development and maintenance. To date, our focus has been on the crucial dynamic development between PB and J+ stages when the pattern is evolving; we have not considered pattern maintenance between J+ and adulthood. Work by Frohnhöfer et al., 2013 has demonstrated that in pfe mutants, which up to stage J+ have similar numbers of melanocytes to WTs (approximately 90% of WT numbers), this drops to approximately 50% of WT by adulthood. Whilst our model correctly predicts melanocyte numbers in pfe compared to WT prior to J+, preliminary simulations up to adulthood suggest this sharp decrease in melanocyte numbers shown in real fish are currently not predicted by our model. This clearly indicates that new biological mechanisms concerning maintenance, likely involving the late differentiating L-iridophores, need to be analysed and incorporated into a model that extends to these later developmental stages.

From an analytical perspective, an advantage of our on-lattice model is its amenability to the derivation of a continuum model, although we note that continuum approximations to off-lattice individual-based models can also be derived. Our model therefore opens up the opportunity for future exploratory work using a continuum model for mutants pfe and nac in order to explore whether pattern formation in these cases individually can be described as Turing patterns and to determine parameter ranges for successful pattern formation.

It will also be interesting to investigate the role of growth in pattern formation and maintenance. We observed in our simulations a lower average stripe straightness (0.92) than to that of real WT fish (0.98). We further observed in our simulations that increased height-to-length ratio is correlated with stripes becoming more tortuous (Figure 14C). Stripes in real fish seem, qualitatively, to not show this effect, and we suspect that our model can be further refined here. A search for mechanisms that might increase stripe straightness will be valuable here, and we note that exploration of our model should allow candidate mechanisms to be identified in silico. Further investigation of this feature in our model once extended to adult pattern maintenance may also be important. For example, casual observation suggests that fish that show a particular height to length bias do not show notably tortuous stripes. Therefore, future work will be to understand what preserves the pattern in these cases. We predict that this may be related to the position of growth. Thus, future work will be to fully evaluate the effects of growth on the final pattern formation.

Another significant aspect which deserves attention concerns dorso-ventral pattern differences. In fish this is characterised by having more pigmented cells at the dorsal region than the ventral region. This is certainly true in adult nac (Figure 4C), which has disproportionately more pigmented xanthophores in the dorsal than the ventral region). We have recently noted the subtle impact of dorso-ventral countershading on the WT zebrafish pigment pattern, including in the stripes themselves, and have identified Agouti as a key regulator of this process (Cal et al., 2019). Furthermore, there are clear dorso-ventral asymmetries in some of the adult mutant patterns: e.g, nac mutants exhibit a strong X0 interstripe and a weak ventral interstripe X1D, and completely lack dorsal interstripes. Our model will allow us to explore possible drivers of this asymmetry, which we hypothesise will include Agouti signalling and also differential domain growth.

In conclusion, our on-lattice model, implementing the current biological understanding of adult zebrafish pigment pattern formation, strongly supports the validity of these experimental interpretations, motivating the detailed investigation of their molecular bases. Our model also highlights areas where knowledge is currently incomplete and, importantly, has allowed in silico investigations to identify plausible mechanisms that require experimental testing.

Appendix 1

Implementing the model

This text provides a full description of how a simulation of a WT fish is implemented. All notation used in appendix 1 is summarised in Supplementary files 1 and 2.

Model overview

We model five different cell types; yellow xanthophores (X), unpigmented xanthoblasts (Xb), black melanocytes (M), silver dense S-iridophores (Id) and blue loose (Il) S-iridophores. We account for the three separate hypodermis layers upon which the cells lie as three separate lattice domains, a xanthophore layer represented by matrix 𝑿, a melanocyte layer represented by matrix 𝑴, and an S-iridophore layer represented by a matrix 𝑰. We denote the length and height of domain D{𝑿,𝑴,𝑰} at time t in terms of the numbers of lattice sites as ΠDH(t), ΠDL(t) respectively. At any time, t, a lattice site at row i and column j in 𝑴 denoted 𝑴(i,j,t) can either be occupied by M, i.e. 𝑴(i,j,t)=M or be unoccupied, 𝑴(i,j,t)=0. Similarly 𝑰(i,j,t) can either be occupied by Il, i.e. 𝑰(i,j,t)=Il, be occupied by Id, i.e. 𝑰(i,j,t)=Id, or be unoccupied, 𝑰(i,j,t)=0. Finally, 𝑿(i,j,t) can either be occupied by X, i.e. 𝑿(i,j,t)=X, Xb, i.e. 𝑿(i,j,t)=Xb or be unoccupied, 𝑿(i,j,t)=0. This captures volume exclusion rules on each lattice layer which require that at any given time, any lattice site on the domain can be occupied by at most one cell. If at any time during the simulation an action is chosen which would break this rule, this action is aborted. To account for the different packing densities of the cell types, lattice sites on 𝑿 and 𝑰 are of size 0.02 mm × 0.02 mm, whilst lattice sites in 𝑴 are of size 0.04 mm × 0.04 mm (Parichy et al., 2000a). We denote these as Δ𝑿=Δ𝑰=0.02 mm and Δ𝑴=0.04 mm. Hence the simulated height and length at any given time t is given as ΩH(t)=ΠDL(t)ΔD, ΩL(t)=ΠDH(t)ΔD. An illustration of the three layers and their corresponding cell types are given in Figure 2A in the main text. The lattices correspond to matrices;

𝑿=[00Xb00X00Xb0X0000X0000X0000X0X00Xb00000X0X0XX0X0XXX00Xb0Xb0000X000X0000X00000]
𝑰=[000000Il000000000Il000IlIl00IdIdIdIdIdIdIdIdIdIdIdIdIdIdIdIdIdIdIdIdIdIdIdId00Il0IlIlIl00000000000000000]
𝑴=[M0MM00000M000M00M0]

where Π𝑿H=Π𝑰H=6, Π𝑿L=Π𝑰L=12, Π𝑴H=3, Π𝑴L=6.

Relating size to stage

Zebrafish development is described in the literature with respect to the standard length (SL) of the fish (Parichy et al., 2009). We associate our domain at time t to different zebrafish developmental stages by calculating a simulated SL (ΩSL) from the size of the domain at time t and relating this to the development stage using Supplementary file 3. Since growth is linear with respect to time (in real life and modelled as so) and we initialise the domain 5.7 mm shorter than the real width (SL), the corresponding simulated SL for any simulation at time t can be calculated by:

(2) ΩSL(t)=ΔD×ΠDL(t)+5.7mm,

where D{𝑿,𝑴,𝑰} and ΠDH(t) is the number of sites in the y direction of domain type D at time t. The simulated height ΩH(t)=ΔD×ΠDH(t) directly corresponds to the height at the anterior margin of the anal fin (abbreviated as HAA) of the fish at all times.

Initial conditions (WT only)

The simulation is initialised to represent the zebrafish half way through stage PB at approximately 25dpf. At this stage the fish is approximately 1 mm in height (i.e., HAA is 1 mm) and 7.7 mm in length (i.e. SL is 7.7 mm). We model the full height and a subsection of the full length for convenience. We set ΩH(0)=1mm and ΩL(0)=2mm. Therefore, Π𝑴L(0)=25, Π𝑴L(0)=50, ΠXH(0)=ΠIH(0)=50 and Π𝑿L(0)=Π𝑰L(0)=100.

The initial occupancies for a WT simulation comprise three rows of Id along the middle of the fish, that is, sites on rows 24 to 26 on 𝑰 are fully occupied with Id cells at time t=0. This represents the S-iridophores which differentiate along the horizontal myoseptum between 21 and 25 dpf. We also place NMT(0)=50M uniformly at random on 𝑴 so that the initial occupancy density of melanocytes is

(3) NMT(0)Π𝑴L(0)×Π𝑴H(0)=0.04,

where NCT(t) for any cell type C{M,X,Xb,Id,Il} denotes the total number of cells of that type on the domain at time t. These scattered melanocytes represent the initial larval pattern at low density dispersed across the domain at the beginning of pattern metamorphosis. Similarly, we occupy the domain so that the initial density of xanthoblasts is 0.4. This density is chosen from evidence that larval xanthophores, which de-differentiate around 5 dpf, proliferate and cover much of the entire domain at the start of metamorphosis (Mahalwar et al., 2014). We do this by placing NXb(0)=2000Xb uniformly at random on 𝑿 so

(4) NXbT(0)Π𝑿L(0)×Π𝑿H(0)=0.4.

This is to represent the de-differentiated larval xanthophores which initially appear around 5 dpf, subsequently lose pigment (becoming xanthoblasts) and proliferate between 5 and 25 dpf. For all other cell types i.e. C{X,Il}, NCT(0) = 0.

Model iteration

The model is updated in continuous time according to the Gillespie algorithm (Gillespie, 1977) (see Figure 2B in the main text for an illustration). In our simulation we allow two different event types: continuous and fixed. A continuous event is an event that can happen at any point throughout the simulated time. For example; melanocyte birth or death. A fixed event is an event that occurs once during the simulation, and happens upon meeting predetermined conditions. In WT fish, we define only one fixed event. We assume there is appearance of metamorphic xanthophores X along the horizontal myoseptum at stage SP. This is justified by observations of shd in which delayed appearance of metamorphic X in interstripe X0 were noted. We found this delayed appearance was important for generating pseudo-stripe patterns in rse and shd which have reduced or entirely absent S-iridophores. We model this by occupying all the sites, regardless of prior occupancy, along the middle three rows of 𝑿 with xanthophores (X) at stage SP that is sites 𝑿(i,j,tSP)=X for all iΠXH(tSP)2-1:ΠXH(tSP)2+1 where tSP denotes the time that the given simulation reaches the start of stage SP.

There are 15 continuous event types. These events are summarised in Supplementary file 4. They comprise cell birth, death, movement and shape transitions as well as growth of the domain.

An overview of how the model is updated is given in Figure 2B in the main text and can be described as follows. At any given time t, the model is first assessed for meeting the criteria of a fixed event. If the model meets the criteria, the fixed event occurs, is subsequently marked as complete and the simulation continues. If no fixed time event is to be implemented then one of the fifteen possible continuous time events is attempted. First an exponentially distributed waiting time

(5) τ=1α0log(1u0),

is generated until the next continuous ‘event’ occurs where α0=i=115αi(t) is the sum of all of the event propensities given in Supplementary file 4 and u0 is a uniformly distributed random number in (0,1) (i.e. u0U(0,1)). Next, we generate u1U(0,1) to determine which event occurs at time t+τ. The event i that satisfies

(6) j=1i-1αj(t)α0u1<j=1iαj(t)α0

is chosen to occur and the domain is updated accordingly (provided conditions required for that event to occur are met). Time is also updated. t=t+τ. We continue this process iteratively, checking for fixed events, then subsequently generating a time for the next continuous event to occur. The process repeats until we reach the end of pattern metamorphosis, marked by length condition ΩSL=13.5mm. The algorithm is stochastic in the sense that, within any given simulation, there is variance with regard to the exact rate and order of event occurrence, just as in real fish.

Modeling cell–cell interactions

Continuous events 1–15 given in Supplementary file 4 are mediated by different cell-cell interactions. Cells interact on the fish skin at both short (neighbouring cells) and long (half a stripe width) range, possibly regulated by direct contact, dendrites, or longer extensions (filopodia or airinemes). In our model, uniform disks, with radii on the order of the distance 0.04 mm account for short-range interactions, and an annulus with an outer radius of approximately half an adult stripe width 0.24mm represent long-range dynamics.

We denote SD(r,k) and LD(r,k) where D{𝑿,𝑰,𝑴} is the set of site positions (i,j) such that D(i,j) is in the short (0.04 mm) or long (0.24 mm) distance respectively from a focal site at position D(r,k). Note that SD(r,k) and LD(r,k) are both different for different domain types due to the different lattice site sizes. SD(r,k) and LD(r,k) are visualised for D=𝑴,𝑿 in Figure 3A–H in the main text. Formulae for these sets are given in Supplementary file 2.

To illustrate this see Figure 3A–H in the main text which compares the number of cells in the short (Figure 3A–D) and long (Figure 3E–H) range distance from a central site, D(r,k), on different domain types. In this figure, M are represented as black circles. X are represented as yellow circles. Figure 2A is a visualisation of sites (marked in red) in S𝑴(r,k), where 𝑴(r,k) is the central site marked in grey. S𝑴M(r,k,t) are the sites marked in red in which a melanocyte resides. The number of melanocytes in the short range distance from 𝑴(r,k) at time t is given by NM,MS(r,k,t)=S𝑴M(r,k,t)=2. Figure 2B is a visualisation of the sites considered when calculating NM,XS (formula given in Supplementary file 2). In this example, NM,XS(r,k,t)=9. To compare the number of M and X in the short range distance of 𝑴(r,k), we consider wNM,MS(r,k,t)=4×2=8 which corresponds to the weighted value of melanocytes in the short range distance with NM,XS(r,k,t)=9 the number of xanthophores in the short range. Figure 3C is a visualisation of sites (marked in red) in S𝑿(r,k), where 𝑿(r,k) is the central site marked in grey. S𝑿X(r,k,t) are the sites marked in red within which a xanthophore resides. The number of xanthophores in the short range distance from 𝑿(r,k) at time t is given by NX,XS(r,k,t)=S𝑿X(r,k,t)=7. Figure 2D is a visualisation of the sites considered when calculating NX,MS(r,k,t) (formula given in Supplementary file 2). In this example, NX,MS(r,k,t)=14. To compare the number of M and X in the short range distance of 𝑿(r,k), at time t, we would compare NX,MS(r,k,t)=14 with NX,XS(r,k,t)=7. Figure 3E is a visualisation of sites (marked in red) in L𝑴(r,k), where 𝑴(r,k) is the central site marked in grey. L𝑴M(r,k,t) are the sites marked in red in which a melanocyte resides. The number of melanocytes in the long range distance from 𝑴(r,k) is given by NM,ML(r,k,t)=L𝑴M(r,k,t)=20. Figure 2F is a visualisation of the sites considered when calculating NM,XL(r,k,t) (formula given in Supplementary file 2). In this example, NM,XL(r,k,t)=16. To compare the number of M and X in the long range distance of 𝑴(r,k), we consider wNM,ML(r,k,t)=4×20=80 which corresponds to the weighted value of melanocytes in the long range distance with NM,XS(r,k,t)=16, the number of xanthophores in the long range. Figure 2G is a visualisation of sites (marked in red) in L𝑿(r,k), where 𝑿(r,k) is the central site marked in grey. L𝑿X(r,k) are the sites marked in red in which a xanthophore resides. The number of xanthophores in the long range distance from 𝑿(r,k) is given by NX,XL=L𝑿X(r,k,t)=2. Figure 3H is a visualisation of the sites considered when calculating NX,ML(r,k,t) (formula given in Supplementary file 2). In this example NX,ML(r,k,t)=20. To compare the number of M and X in the long range distance of 𝑿(r,k), we consider NX,ML(r,k,t)=20 with NX,XL(r,k,t)=2.

We denote SDC(r,k,t)SD(r,k) where C lies in layer D to be the sites in SD(r,k) occupied by cell type C{M,X,Id,Il,Xb} at time t. We denote NCi,CjS(r,k,t), (NCi,CjL(r,k,t)) as the number of cells of type Cj in the short (long) range distance 0.04 mm, (0.24 mm) from cell type Ci at Di(r,k). Hence the number of cells of type Ci in the short distance from another cell of the same type at D(r,k) at time t is given by NCi,CiS(r,k,t)=|SDCi(r,k,t)|. The formulae for NCi,CjS(r,k,t), NCi,CjL(r,k,t) where Ci does not equal Cj are more complicated and are given in Equations 7, 8 (and Supplementary file 2 for reference).

(7) NC1,C2S(r,k,t)={SD2C2(r,k,t) if D1=D2orD1 and D2{𝑿,𝑰}.(i,j)SD1(r,k)𝟙𝑴(i2,j2)=M if D1𝑴,D2=𝑴.(i,j)S𝑴(r,k)𝟙D2(2i,2j)=C2+𝟙D2(2i-1,2j)=C2+𝟙D2(2i,2j-1)=C2+𝟙D2(2i-1,2j-1)=C2 if D1=𝑴,D2𝑴.
(8) NC1,C2L(r,k,t)={LD2C2(r,k,t) if D1=D2orD1 and D2{𝑿,𝑰}.(i,j)LD1(r,k)𝟙𝑴(i2,j2)=M if D1𝑴,D2=𝑴.(i,j)L𝑴(r,k)𝟙D2(2i,2j)=C2+𝟙D2(2i-1,2j)=C2+𝟙D2(2i,2j-1)=C2+𝟙D2(2i-1,2j-1)=C2 if D1=𝑴,D2𝑴.

Simply, where domain types D do not have the same lattice site size ΔD, the focal site coordinates (i,j) of the local neighbourhood undergo a transformation for the sites on a domain with a different size. For example, each site 𝑿(i,j), corresponds to a quarter of site 𝑴(i2,j2), this explains case two in Equation 7. Similarly each site 𝑴(i,j) corresponds to four sites on 𝑿, specifically 𝑿(2i,2j), 𝑿(2i,2j-1), 𝑿(2i-1,2j) and 𝑿(2i-1,2j-1). This explains case three in Equations 7 and 8. Note that since M is four times larger than all other cells in our simulation, we provide a weighting system when comparing M with C\M in some cases.

Boundary conditions

Boundary conditions are periodic along the horizontal boundaries and reflecting across the vertical boundaries. We implement periodic boundary conditions along the horizontal axis based on the assumption that the rate at which cells leave along this axis is approximately equal to the rate at which cells enter the domain at the opposite side.

Continuous time events

In this text, we will describe how each of the fifteen events are simulated upon being selected to occur by the Gillespie algorithm. An overview of all continuous time events and their corresponding rates is given in Supplementary file 4.

Movement (continuous time events 1–5)

We implement cell movement so that cells are biased towards cell types they are attracted to and away from cell types they are repelled by. The direction of the cell’s movement is determined using an on-lattice attraction-repulsion mechanism based on a model described by Dini et al., 2018 and is detailed as follows. If a cell is chosen to move, it is able to move in one of eight different possible orientations denoted by O𝑶={No,So,Ea,We,NE,SE,SW,NW}. (See Figure 3(D) in the main text for an illustration of these directions).

The directional neighbours of each cell is given in Figure 2E,F. Figure 2E demonstrates the directional neighbourhoods of a focal cell C{X,Xb,Id,Il} located in position D(i,j) (marked in grey) where D{𝑿,𝑰}. For the cell moving one space on D, the cell considers the occupancy of sites in SkΔXD(i,j) marked in red for D=𝑴 (top) and D{𝑿,𝑰} (bottom). (Q–X) Figure 2F demonstrates the directional neighbourhoods of a focal cell M located in position 𝑴(i,j) (marked in grey). For the cell M moving one space on 𝑴, the cell considers the occupancy of sites in SkΔMD(i,j) marked in red for D{𝑿,𝑰} (top) and D=𝑴 (bottom).

We define the probability that a cell attempts to move in a given direction with orientation O as PO (where O𝑶PO=1). We calculate PO, using bias vector,

(9) v¯=𝑨a¯+𝑹r¯,

where 𝑨 (𝑹) is a matrix whose entries are the number of neighbouring cells within different segments of the attraction (repulsion) range and a¯ (r¯) is a weight of attraction (repulsion) vector.

Matrices 𝑨 and 𝑹 are defined as follows.

(10) 𝑨=[EaC,C1(r,k,t)EaC,C2(r,k,t)EaC,C3(r,k,t)WeC,C1(r,k,t)WeC,C2(r,k,t)WC,C3(r,k,t)NoC,C1(r,k,t)NoC,C2(r,k,t)NoC,C3(r,k,t)SoC,C1(r,k,t)SoC,C2(r,k,t)SoC,C3(r,k,t)NWC,C1(r,k,t)NWC,C2(r,k,t)NWC,C3(r,k,t)NEC,C1(r,k,t)NEC,C2(r,k,t)NEC,C3(r,k,t)SWC,C1(r,k,t)SWC,C2(r,k,t)NEC,C3(r,k,t)SEC,C1(r,k,t)SEC,C2(r,k,t)SEC,C3(r,k,t)],
(11) 𝑹=[WeC,C1(r,k,t)WeC,C2(r,k,t)WeC,C3(r,k,t)EaC,C1(r,k,t)EaC,C2(r,k,t)EaC,C3(r,k,t)SoC,C1(r,k,t)SoC,C2(r,k,t)SoC,C3(r,k,t)NoC,C1(r,k,t)NoC,C2(r,k,t)NoC,C3(r,k,t)NEC,C1(r,k,t)NEC,C2(r,k,t)NEC,C3(r,k,t)NWC,C1(r,k,t)NWC,C2(r,k,t)NWC,C3(r,k,t)SEC,C1(r,k,t)SEC,C2(r,k,t)SEC,C3(r,k,t)SWC,C1(r,k,t)SWC,C2(r,k,t)NEC,C3(r,k,t)]

where for example, OC,C1(r,k,t) denotes the number of cells of type C1{X,M,Id,Il,Xb} within range 0.04 mm, that are orientation O from the focal cell C located in position D(r,k) at time t. This number is calculated differently depending on whether the focal site lies on lattice domain D{𝑴,𝑿,𝑰}. To compute OC,C1(r,k,t) we calculate the following. For any O we define subsets SDO(r,k,t) of SD(r,k,t) (where cell type C lies on domain type D) for each direction O;

(12) SDNo(r,k,t)={(i,j)SD(r,k,t)|j>k},
(13) SDSo(r,k,t)={(i,j)SD(r,k,t)|j<k},
(14) SDE(r,k,t)={(i,j)SD(r,k,t)|(|ri|<ΠDL(t)|ri| and i>r)or (|ri|>ΠDL(t)|ri| and i<r)},
(15) SDW(r,k,t)={(i,j)SD(r,k,t)|(|ri|<ΠDL(t)|ri| and i<r)or (|ri|>ΠDL(t)|ri| and i>r)}.

These are the sets of sites in the short range distance, north, south, east and west of focal site D(r,k). Examples for D=𝑴 and 𝑿 are sites marked in red on the left of Figure 2E–F respectively. By extension we define

(16) SDNE(r,k,t)={(i,j)SD(r,k,t)|jk and (|ri|<ΠDL(t)|ri| and ir)
(17) or (|r-i|>ΠDL(t)-|r-i| and ir)},
(18) SDNW(r,k,t)={(i,j)SD(r,k,t)|jk and (|ri|<ΠDL(t)|ri| and ir)
(19) or (|r-i|>ΠDL(t)-|r-i| and ir)},
(20) SDSE(r,k,t)={(i,j)SD(r,k,t)|jk and (|ri|<ΠDL(t)|ri| and ir)
(21) or (|r-i|>ΠDL(t)-|r-i| and ir)},
(22) SDSW(r,k,t)={(i,j)SD(r,k,t)|jk and (|ri|<ΠDL(t)|ri| and ir)
(23) or (|r-i|>ΠDL(t)-|r-i| and ir)}.

These are the sets of sites in the short range distance, north-east, north-west, south-east and south-west of focal site D(r,k). We denote SDC,O(r,k,t)SDO(r,k,t) to be all the sites in short-range distance in orientation O{No,So,Ea,We,NE,SE,SW,NW} occupied by cell type C. Hence we define

(24) SDC,O(r,k,t)={(i,j)SDO(r,k,t)|D(i,j)=C}.

Therefore to compute the number of cells in each of the eight directions, we compute where C lies on D and C1 lies on D1 by;

(25) OC,C1(r,k,t)={|SD1C,O(r,k,t)| where D=D1orD,D1{X,I}.SD1O(r,k)𝟙M(r2,k2)=M where D=M,D1M.SMO(r,k)𝟙D(2r,2k)=C+𝟙D(2r1,2k)=C+𝟙D(2r,2k1)=C+𝟙D(2r1,2k1)=C where DM,D1=M.

Next we define the weight of the attraction-repulsion vector by:

(26) a¯=[aCC1aCC2aCC3],r¯=[ρCC1ρCC2ρCC3]

where aCC1,aCC2,aCC3 (ρCC1,ρCC2,ρCC3) are the weights of attraction (repulsion) of cells of type C to cells of the type C1, C2 and C3, respectively defined in Supplementary file 5. Vectors a¯ and b¯ are the vectors used to calculate v¯ in Equation 9. To account for the difference in distance between adjacent neighbours and diagonal neighbours we normalise accordingly so that movement diagonally is proportionately less frequent than moving to adjacent squares.

(27) P^W=2v13q=14vq,
(28) P^E=2v23q=14vq,
(29) P^So=2v33q=14vq,
(30) P^No=2v43q=14vq,
(31) P^NE=v53q=58vq,
(32) P^NW=v63q=58vq,
(33) P^SE=v73q=58vq,
(34) P^SW=v83q=58vq.

Finally we normalise PO s.t.;

(35) PO=P^OO𝑶P^O.

To summarise, we simulate movement as follows. Given the event is chosen such that a cell of type C is allocated to move. We choose a random cell of type C. We call its position D(r,k). Next values PO are calculated using Equation 27–33 and Equation 35. A random number uU(0,1) is generated. The cell attempts movement to D(r,k-1) (West) if u[0,PW), D(r,k+1) (East) if u[PW,PW+PE)etc. up to [O𝑶PO-PSW,O𝑶PO]. The cell movement is successful and the cell moves from site D(r,k) to the nearest neighbouring site in direction O, provided this site is empty. Otherwise the movement is aborted. Movement bias is specified by weights of attraction and repulsion (a¯ and r¯). Relevant parameters for aCiCj and ρCiCj for different Ci and Cj are given in Supplementary file 5 . If no such interaction bias is specified then there is no known attraction or repulsion dynamics between those cell types in the literature and so we assume that movement is not biased by this cell type. Note that for cell Xb there are no known short range interactions between Xb and other cell types so we model Xb movement as random, i.e, PO=18, for O𝑶 at all times.

S-iridophore, xanthophore and xanthoblast proliferation (continuous time events 6–7)

Given the event is chosen such that a cell of type C{X,Xb,Id,Il} is determined to proliferate we choose a random cell of type C whose position is given by D(r,k), D{𝑿,𝑰}. Next a random number uU(0,1) is generated. This number is used to determine a neighbouring site into which a mother cell can place a daughter cell. This site is D(r,k+1) if u[0,1/4), D(r+1,k) if u[1/4,1/2), D(r-1,k) if u[1/2,3/4) or D(r,k-1) if u[3/4,1].

For cell types C{Id,Il,Xb}, a proliferation event is successful if the site chosen is empty, and a new cell of the same type is placed into the chosen site, otherwise the event is aborted. In the case C=X, a proliferation event is successful if the site chosen for the daughter cell is unoccupied and

(36) NX,XS(r,k,t)+NX,IdS(r,k,t)>NX,IlS(r,k,t)+NX,MS(r,k,t).

Equation 36 is based on the following assumptions: (1) Dense S-iridophores promote xanthoblasts differentiation into xanthophores in the short range (Patterson and Parichy, 2013). Dense S-iridophores express xanthogenic Colony Stimulating Factor-1 (Csf1) (Patterson and Parichy, 2013) which is essential for xanthophore differentiation, proliferation and survival, allowing unpigmented xanthoblasts near to the dense S-iridophores to mature into xanthophores (Frohnhöfer et al., 2013; Walderich et al., 2016); (2) melanocytes inhibit xanthophore specification in the short range (Nakamasu et al., 2009).

Melanocyte differentiation (continuous time event 8)

If a melanocyte differentiation event is specified then a site 𝑴(r,k) is chosen at random. A differentiation event is successful and a new M appears in this site if the site 𝑴(r,k) is empty and the following is true:

(37) NM,XL(r,k,t)+NM,IdL(r,k,t)>wαNM,ML(r,k,t)+β
(38) and NM,XS(r,k,t)γwNM,MS(r,k,t),
(39) and NM,IdL(r,k,t)κ,

where α=2.5, β=γ=κ=3,w=4. Equation 37 is based on the findings that dense S-iridophores and xanthophores promote the differentiation of melanocytes in the long range. This conclusion is drawn from observations in ablation experiments (Nakamasu et al., 2009) and in the pfe mutant, which retains a high number of M(Frohnhöfer et al., 2013). Equation 38 is based on observations that melanocytes and xanthophores compete in the short range (Nakamasu et al., 2009). Equation 39 is based on observations that in WT fish, melanocyte centers rarely overlap with dense S-iridophores; however, melanocytes frequently settle adjacent to dense S-iridophores, suggesting short range inhibition (Patterson and Parichy, 2013).

We assume there is some melanocyte differentiation into empty space that is independent of cues from other cells. This is from observations of double mutant fish nac;pfe that do not produce S-iridophores or xanthophores but phenotypically display uniformly distributed melanocytes at adulthood (Frohnhöfer et al., 2013). Therefore, alternatively we also allow successful M differentiation if we generate a random number uU(0,1) and we find that u<0.01 in combination with the condition

(40) NM,XS(r,k,t)+NM,IdS(r,k,t)+NM,MS(r,k,t)=0.

Within each attempt, criterion (Equation 37–39) is tried first. If this is not successful, criteria (Equation 40) is tested instead.

Xanthoblast differentiation (continuous time event 9)

If a xanthoblast differentiation event is chosen then an Xb is chosen at random from 𝑿. Suppose the chosen Xb lies in site 𝑿(r,k). The differentiation event is successful and a X replaces the Xb in this site if the site 𝑴(r2,k2) is not occupied by a cell M (as melanocytes and xanthophores are known to compete in the short range [Nakamasu et al., 2009]) and if the following is true.

(41) NX,XS(r,k,t)+NX,IdS(r,k,t)>NX,MS(r,k,t)+NX,IlS(r,k,t).

Equation 41 is based on the following assumptions; (1) Dense S-iridophores promote xanthoblasts differentiation into xanthophores in the short range (Patterson and Parichy, 2013). Dense S-iridophores express Csf1, allowing unpigmented xanthoblasts near to the dense S-iridophores to mature into xanthophores (Frohnhöfer et al., 2013; Walderich et al., 2016). (2) melanocytes inhibit xanthophore specification in the short range (Nakamasu et al., 2009).

Alternatively the differentiation event is successful if:

(42) NX,XS(r,k,t)+NX,MS(r,k,t)=0,

holds and a randomly generated number uU(0,1) is such that u<0.01, or Equation(42) holds and either

(43) t>31dpfandNMT(t)+NXT(t)+NIdT(t)ΠXL(t)×ΠXH(t)<0.2,
(44) t>51dpfandNMT(t)+NXT(t)+NIdT(t)ΠXL(t)×ΠXH(t)<0.4,

is true i.e, the total cell density with all cells combined is less than either 0.2 or 0.4 for different time milestones. We enforce a lower probability (u < 0.01) for the alternative event (Equation 42) since this differentiation event is not influenced by cues from other cells and thus is less likely. Equation 43 and Equation 44 are based on assumptions that if the domain is empty after some key time points, this promotes delayed differentiation of xanthoblasts as is observed in nac (Frohnhöfer et al., 2013). Since melanocytes compete in the short range with xanthophores (Nakamasu et al., 2009) we only allow this to occur when constraints described by Equation(42) are simultaneously held, for consistency.

Iridophore transitions (continuous time event 10–11)

If an Id to Il shape transition is chosen, then an Id is chosen uniformly at random from domain 𝑰 and evaluated for meeting transition criteria. The transition is successful and the chosen Id in position 𝑰(r,k) is replaced with Il in this site if either both

(45) NId,XL(r,k,t)>Lx and NId,XS(r,k,t)<Sx,

or alternatively

(46) NM,MS(r,k,t)Sm,

where Lx,Sx,Sm=12, 9, one respectively. Alternatively if an Il to Id shape transition is chosen, then an Il is chosen at random from domain 𝑰 and evaluated for meeting transition criteria. This transition is successful and the chosen Il in position 𝑰(r,k) is replaced with Id in this site if either

(47) NIl,XL(r,k,t)<Lx2 and NIl,MS(r,k,t)<Sm2,

or alternatively

(48) NIl,XS(r,k,t)>Sx2 and NIl,MS(r,k,t)<Sm2,

where Lx2=16, Sx2=4 and Sm2=1. These conditions are based on observations of the induction set (Frohnhöfer et al., 2013) (for more details see Section "Modelling assumptions" in the main text). The parameters Lx, Sx, Sm, Lx2, Sm2, Sx2 were chosen to give straight stripes with few breaks at stage J+ in WT simulations. However, with some small variations to these parameters,the patterns generated are qualitatively similar.

Melanocyte death (continuous time event 12)

If a melanocyte death event is chosen, then a cell of type M is chosen at random from 𝑴. Suppose the chosen M lies in site 𝑴(r,k). The death event is successful and the melanocyte is removed from site 𝑴(r,k) i.e. 𝐌(r,k) is set to 0, if

(49) wNM,MS(r,k,t)<NM,XS(r,k,t).

Equation 49 is based on findings that xanthophores and melanocytes compete in the short range (Nakamasu et al., 2009). Alternatively, the melanocyte death is successful if a randomly generated number uU(0,1), is such that u<0.001 and both of the following two cases hold:

(50) σNM,XL(r,k,t)  NM,ML(r,k,t),
(51) NM,IlS(r,k,t) < ω,

where σ=ω=3. Equation 50 is based on findings that M appear to inhibit the survival of M in the long range and that X promote the survival of M in the long range (Takahashi and Kondo, 2008). In Equation 51, Il are a proxy for L-iridophores, which promote the survival of M in the short range (Frohnhöfer et al., 2013). We found this equation was important for maintaining melanocytes in mutant pfe in our simulations. Moreover, the two equations Equation 50 and Equation 51 combined were important for maintaining melanocytes in the interstripes in pfe (as seen in real fish) in our simulations. We found it was important that this alternative event had a low probability of success (0.001) but not zero. In pfe there are no xanthophores so Equation 50 enforces melanocyte death where the number of melanocytes in the long range was greater than zero. Therefore, when there was too high a probability of success, in pfe simulations the stripe, interstripe structure that is usually maintained between WT and pfe was broken into randomly dispersed spots of melanocytes distance 0.24 mm apart. On the other hand, when it was too low, we did not observe melanocytes in the interstripe in our pfe simulations, unlike in real pfe.

Xanthoblast pulling event (continuous time event 13)

Suppose a ‘xanthoblast pull melanocyte’ event is chosen at time t, then we choose an Xb uniformly at random from all Xb occupying 𝑿 that meets the following criterion: suppose the chosen Xb resides in site 𝑿(r,k), then the corresponding site on 𝑴, 𝑴(r2,k2) must be empty. We simulate this by randomly choosing Xb on 𝑿, and checking whether 𝑴(r2,k2) is empty. If so, we continue to the next stage. Otherwise we repeat through all Xb on 𝑿 until we either find a suitable Xb. If no Xb satisfy this criteria at time t then we abort the cell-cell pulling event. Given we find a suitable Xb, the chosen Xb will attempt to attach and pull a melanocyte within range (airinemes extend to a length of up to 5–6 cell diameters away or 0.1–0.12 mm [Eom et al., 2015]) to position 𝑴(r2,k2). We simulate this as follows. First, we generate a random position (i,j) uniformly at random from the 28 possible lattice positions euclidean distance 0.1 mm from the site X(r,k) given by P. These sites are shown in Figure 1A–B and their coordinates are given by;

(52)P={(r+5,k),(r5,k),(r,k+5),(r,k5)(53)(r+1,k+5),(r1,k+5),(r+1,k5),(r1,k5),(54)(r+2,k+5),(r2,k+5),(r+2,k5),(r2,k5),(55)(r3,k4),(r+3,k+4),(r3,k+4),(r+3,k4),(56)(r4,k3),(r+4,k+3),(r4,k+3),(r+4,k3),(57)(r+5,k+1),(r5,k+1),(r+5,k1),(r5,k1),(58)(r+5,k+2),(r5,k+2),(r+5,k2),(r5,k2)}.

Next we translate position (i,j) on 𝑿 to its corresponding position on 𝑴. If 𝑴(i2,j2,t)=M, then the melanocyte is pulled from its site 𝑴(i2,j2) into 𝑴(r2,k2)i.e.𝑴(i2,j2,t)=0 and 𝑴(r2,k2,t)=M. The action is deemed complete. Otherwise, if 𝑴(iM,jM,t)=0, i.e, the site is empty, then the process repeats for the same Xb. From here, P becomes P\(i,j) and a new (i,j) pair are generated uniformly at random from P until either an M is found within the range specified or until P is empty (i.e. after 28 tries). At this point, we will have determined that there were no melanocytes close enough to the xanthoblast for the xanthoblast to successfully pull.

Growth (continuous time events 14–15)

Given a growth event is chosen in the horizontal (vertical) direction, from each column (row), a site location is chosen uniformly at random at which insertion of new empty site above or below (randomly chosen) will occur. See Appendix 1—figure 1C-D for an example of a growth event in the vertical direction.

Appendix 1—figure 1
Implementation of xanthoblast-melanocyte communication and domain growth.

(AB) Range for xanthoblast cytonemes in our model. Figure:range_pull The range of xanthoblast cytonemes on 𝑿. The centre lattice site (marked in purple) represents site 𝑿(r,k). The lattice sites labelled in red are those sites 𝑿(i,j) that lie in the range of the cytonemes of an Xb occupying site 𝑿(r,k) (the sites given by set P). Figure:range_pull_mel The range of xanthoblast cytonemes on 𝑴. The sites indicated with in orange or yellow are those sites that can be reached by the range of Xb in 𝑿(i,j). The colour of any given site 𝑴(i,j) reflects the number of corresponding sites on 𝑿(2i,2j), 𝑿(2i-1,2j-1), 𝑿(2i-1,2j), 𝑿(2i,2j-1) that are within range of the cytonemes (marked in red in (A)) and thus the probability of the site being chosen in relation to other coloured sites. Specifically, sites labelled in orange are twice as likely to be chosen than those sites marked in yellow. (CD) An example growth event in the vertical direction. Squares represent sites in a given domain. Yellow circles represent cells occupying these sites. (C) First a site is chosen at random from each column (marked in red). (D) Next a new site is inserted either above or below the first site with probability 12. Each of the chosen sites are marked in dark grey.

To account for the different lattice sizes, each time a growth event in the horizontal (vertical) direction occurs, Π𝑴H(t) (Π𝑴L(t)) increases by one and Π𝑿H(t), Π𝑰H(t) (ΠXL(t), ΠIL(t)) increases by two. This means that when a growth event is chosen, the rules described above occur once in the 𝑴 layer, and twice consecutively in 𝑿 and 𝑰. This is to ensure that the domain size remains consistent between the three layers. That is,

(59) ΠXL(t)×ΔX=ΠIL(t)×ΔI=Π𝑴L(t)×ΔM,
(60) Π𝑿H(t)×ΔX=Π𝑰H(t)×ΔI=Π𝑴H(t)×ΔM,

for all time t.

Appendix 2

Mutant implementation

In this text, we describe how the WT fish simulation is altered in our model to replicate changes known to be present in a variety of mutant fish.

Simulating the shd mutant

The gene shady (shd) encodes zebrafish leukocyte tyrosine kinase (Ltk) which plays a role in S-iridophore specification (Lopes et al., 2008). As a result, strong shd mutants lack S-iridophores. To simulate this defect, we remove all S-iridophores from the initial domain, i.e, we set NIlT(0)=NIdT(0)=0. Since new S-iridophores are only generated by the proliferation of existing S-iridophores, this means that NIlT(t)=NIdT(t)=0 for all time t. (Note: As a result the propensities of events 4, 5, 7, 10 and 11, that is, the movement of Il, movement of Id, proliferation Il and Id, transition of Id to Il and transitions of Il to Id are always zero during this simulation.)

Simulating the pfe mutant

Gene pfeffer (pfe) (alternatively known as salz (sal)) encodes for colony stimulating factor one receptor (csf1ra) which plays a role in xanthophore specification and migration. In strong alleles, adult fish exhibit no detectable xanthophores in the body of embryos. To simulate this defect we remove all xanthophores and xanthoblasts from the initial domain, i.e. we set NXb(0)=NX(0)=0. We also remove the fixed event at stage SP where newly differentiated X cells appear along the horizontal myoseptum. Since new xanthophores are only generated by the proliferation of existing xanthophores, and the differentiation of xanthoblasts this means that NXbT(t)=NXT(t)=0 for all time t. (Note: As a result the propensities of events 2, 3, 6, 9, 15, that is, the movement of Xb, movement of X, proliferation X and Xb, differentiation of Xb to X and Xb pulling of M are always zero during this simulation.)

Simulating the nac mutant

The gene nacre (nac) encodes transcription factor Mitfa (Lister et al., 1999). nac mutants lack melanocytes throughout embryonic and larval development (Lister et al., 1999). To simulate this defect, we remove all melanocytes from the initial conditions, i.e, we set NM(0). We also set the propensity of event 8, the differentiation of new M, to be 0 for all time t. Since this is the only way new M can be produced, NM(t)=0 for all time t. (Note: As a result the propensities of events 1 and 14, that is, the movement of M, and death of M will also always be zero during this simulation.)

Simulating the double mutants: nac;pfe, shd;pfe, shd;nac

To simulate the double mutants we combine the effects of two of the three aforementioned mutation types. For example, to simulate nac;pfe we combine the effects of nac and pfe listed above.

Simulating the rse mutant

Rose (rse), encodes the Endothelin receptor B1a (Krauss et al., 2014) and has been shown to acts cell-autonomously in S-iridophores; homozygous mutants result in a reduction of S-iridophores to approximately 20% of that seen in WT (observed in stage PB and adult fish [Frohnhöfer et al., 2013]). To simulate rse we change the initial conditions so that NId(0) is one fifth of the usual number. In our WT simulations, NId(0) is Π𝑰L(0)×3=300 as Id occupy the three central rows. In rse we set NId(0)=60, and we place these Id uniformly at random within the space of the central three rows, since some S-iridophores still appear along the horizontal myoseptum in rse. We do this by generating a random number r that is uniformly distributed between 24 and 26 (the centre of the horizontal axis is at 25), and a random number k1 that is uniformly distributed between 0 and Π𝑰L(0)=100. If 𝑰(r1,k1) is empty we place a Id cell in this site, otherwise we generate a new position until there are 60 cells of type Id on 𝑰. Furthermore since new S-iridophores are only produced by proliferation of pre-existing S-iridophores, we also reduce the rate of proliferation of Id and Il to one fifth of the usual number. Therefore, we set the propensity of event 7, α7(t)=0.24NIdT(t)+NIlT(t)60×24 for all time t.

Simulating the sbr mutant

Adult sbr mutants exhibit delayed S-iridophore shape transition changes of dense to loose caused by truncations in Tight Junction Protein 1a (Fadeev et al., 2015). To simulate sbr, we reduce the rate at which Id attempts transition to Il to a fortieth of its usual value. Hence the propensity of event 10 becomes α10(t)=NId(t)60×24×40. This reduction acts as a proxy for the delay between receiving a signal and changing shape to loose form.

Simulating the cho mutant

Mutant larvae with mutation cho lack the horizontal myoseptum (Svetic et al., 2007). As a result, dense S-iridophores cannot travel through their usual pathway to generate the initial strip of dense S-iridophores at stage PB seen in WT. Instead loose S-iridophores appear later at stage PR, uniformly across the domain. To simulate the effects of cho, we remove the initial three rows of Id so NId(0)=0, as S-iridophores cannot appear along the horizontal myoseptum. Furthermore, we remove the fixed event of metamorphic xanthophores appearing along the horizontal myoseptum at stage SP. To simulate the appearance of loose S-iridophores at stage PR, we provide a new fixed event when the simulation reaches stage PR. At this point, we place 300 (the same number that usually initially occupies the domain in WT) loose S-iridophores in sites uniformly at random across 𝑰 at stage PR. We do this by generating a random number r between 0 and Π𝑰L(tPR) and a random number k between 0 and Π𝑰H(tPR), where tPR is the time when the simulation first enters stage PR. If 𝑰(r,k) is empty we place a Il cell in this site, otherwise we generate a new r, k and continue until there are 200 cells of type Il on 𝑰.

Simulating the seurat mutant

Homozygous seurat mutants develop fewer adult melanocytes, thus forming irregular spots rather than stripes. This phenotype arises from lesions in the gene encoding Immunoglobulin superfamily member 11 (Igsf11) (Eom et al., 2012) which encodes a cell surface receptor containing two immunoglobulin-like domains which is expressed autonomously by the melanocyte lineage. Igsfl1 promotes the migration and survival of these cells during adult stripe development as well as mediating adhesive interactions in vitro.

To model seurat we reduced the rate at which melanocytes could differentiate to a twentieth of the usual rate. Hence the propensity of event eight becomes α8(t)=0.05×NId(t)2×60×24. This was to reflect the inhibition of migration of melanoblasts across the domain and increased the rate of attempted melanocyte death to one hundred times per day (usually once per day). Hence the propensity of event 12 becomes α12(t)=×NM(t)100×60×24. No other interactions were altered.

Simulating the leo mutant

The gene leo encodes Connexin39.4 (Cx39.4). The leo mutant displays a spotted pattern across the flank of the fish. In Section "An in silico investigation into the function of the leo gene" of the main text we describe how we derive the following hypotheses about the impacts of a mutation in the leo gene;

  • Hypothesis 1: Melanocytes are not repelled by xanthophores

  • Hypothesis 2: Xanthophores do not promote the survival of melanocytes in the long range.

  • Hypothesis 3: Xanthophores do not promote the change of S-iridophores from dense to loose in the long range.

  • Hypothesis 4: Melanocytes lose death signals from local dense S-iridophores and as a result, can differentiate in dense S-iridophore zones.

  • Hypothesis 5: Melanocytes lose directed signalling from S-iridophores and hence in the absense of xanthophores differentiate randomly.

To model hypothesis 1, we change the parameter rmx, the parameter governing repulsion of melanocytes from xanthophores to 0. To incorporate hypothesis two we remove the criteria for successful melanocyte death given in Equation 50 in Appendix 1. To model hypothesis 3, we reduce the rate of successful signalling of xanthophores to iridophores to change to loose form. To do this, if a S-iridophore transition is chosen to occur then a number distributed uniformly at random is generated. If this number is less than 0.5, then normal transition signalling occurs (as described in Appendix 1). If the number is greater than 0.5 then the xanthophores send a signal for S-iridophores to transition to dense, that is, the number NI,XS is changed to five and the number NI,XL is changed to 2. To model hypotheses 4 and 5, we change the melanocyte differentiation success as follows. A melanocyte successfully differentiates into a position on the lattice, if there are no xanthophores on the domain, or if there are xanthophores on the domain and there are three times as many xanthophores in the long range than the number of melanocytes in the long-range.

Appendix 3

Quantitative measures

In this text, we describe in more detail how we take quantitative measurements of our simulations.

Tortuosity

To determine the tortuosity of the X0 stripe, for a specific time point in our simulations the following steps are taken.

  • An occupancy matrix 𝑿^ of zeros and ones is generated such that a matrix entry 𝑿^(i,j)=1 if 𝑿(i,j)=X and 0 otherwise. Sites (i,j) such that 𝑿^(i,j)=1 are yellow and are white otherwise.

  • Representative matrix 𝑿^ is ‘cleaned’ using matlab functions bwmorph(𝑿^,‘clean’) and bwmorph(𝑿^,‘bridge’) consecutively. ‘Clean’ removes any anomalous xanthophores that are not connected (adjacent to) other xanthophores. ‘Bridge’ adds extra xanthophores where there are holes in the population pattern.

  • An algorithm is applied to 𝑿^ to create the outline of the X0 stripe. The algorithm to create the top line (lt) is given below in Algorithm 1. A similar algorithm is used to generate the bottom line (lb).

  • A line L that represents the middle of the stripe is given by lm(i)=lt(i)+lb(i)2 for i=1,Π𝑿L.

  • We smooth L for analysis by applying matlab function smooth. ‘Smooth’ smooths the data in the column vector y using a moving average filter.

  • Finally, we calculate the tortuosity of the line by computing the total length of L divided by the distance (algorithm given below in Algorithm 2).

Algorithm 1. Algorithm to generate a representative line for the top of stripe X0 (lt)

j=ΠXH2
  for i = 1 : ΠXL
    Strike = 0
      if X^(i,j)0 or (𝑿^(i,j)=0 and 𝑿^(i,j+1)0) then
       while strike<2 and j<Π𝑿H do
       j=j+1
        if 𝑿^(i,j)=0 then
       strike=strike+1
        else
       strike=0
       end if
       end while
      j=j-2
     else
    while𝑿^(i,j)=0 and j>Π𝑿H do
   j=j-1
  end while
  end if 
lt(i)=j
 end

Initialise X0 interstripe search in the center.

Initialise the number of consecutive empty sites in a column to be zero.

If we are near to other X, keep moving up to determine the upper bound of the interstripe.

Increment the number of consecutive empty sites in a column by one.



Remove the two consecutively empty sites from the total count.


If we have overshot the interstripe keep moving downwards until we reach the top of the interstripe.
Algorithm 2. Algorithm to calculate the tortuosity (tort) of line L generated from a simulated X0
total = 0
for i = 2 : ΠXL
    total=total+(lm(i)lm(i1))2+1
end
tort=totalΠXL



Compute the length of Lusing the euclidean distance between consecutive points.

To measure the tortuosity of the stripes of real fish, a photograph was taken of the fish at stage J+. Next, matlab function ‘getpts’ was used to mark the outline of the X0 stripe by a set of points [x¯,y¯] where vector x¯ is length k. The tortuosity of this line was then computed using Algorithm 3.

Algorithm 3. Algorithm to calculate the tortuosity (tort) of line L generated from a real fish
total = 0
for i = 2 : k
    total=total+(x(i)x(i1))2+(y(i)y(i1))2
end
tort=total(x(1)x(k))2+(y(i)y(k))2



Compute the length of Lusing the euclidean distance between consecutive points.

X0 interstripe width

To determine the width of the X0 interstripe, used in Section "Quantitative analysis of simulations" of the main text, first, occupancy matrices of 𝑿^ and 𝑰^ of zeros and ones is generated such that a matrix entry 𝑿^(i,j)=1 if 𝑿(i,j)=X and 0 otherwise. Similarly 𝑰^(i,j)=1 if 𝑰(i,j)=Id and 0 otherwise. We then generate lt and lb using either X^ or I^ as described in Section "Necessity of S-iridophore assumptions". From these values, we generate IS (X0 interstripe width) by computing;

(61) IS=i=1Π𝑿Llt(i)-lb(i)Π𝑿L.

Note that the X0 interstripe width is computed using whichever is more appropriate of X on 𝑿 and Id on 𝑰 given the mutation. For example, pfe does not have xanthophores, so we would use the distribution of Id on 𝑰 to determine the width of X0 in this case.

Adapting the pair correlation function (PCF)

PCFs characterise spatial patterns by calculating a numerical value for the deviation from the situation in which the same number of agents are distributed uniformly at random. In this paper, we use the square uniform PCF (Gavagnin et al., 2018) developed for on-lattice systems of agents where distance is measured using the uniform norm. This PCF was originally developed for determining pair correlation between single agents types on a lattice, however, in a technique similar to Dini et al., 2018 we adapt it here so that it can also be used for identifying correlation between two different types of agents (cell types). First we define the PCF. For each distance m the PCF at distance m is given by:

(62) PCF(m)=cC1,C2(m)𝔼[c^C1C2(m)]

where cC1,C2(m) is defined as the number of cells of type C1 we would expect to find at distance m from cells of type C2 using the uniform metric under zero flux boundary conditions. 𝔼[c^C1,C2(m)] if the cells were positioned uniformly at random on the domain. This can be calculated by counting the number of agents of this distance manually from the lattice. To compare cell type M with cells that lie on domains other than 𝑴 we must transform 𝑴 into a matrix 𝑴^ of size Π𝑿H×Π𝑿L where 𝑴^(r,k)=M if 𝑴(r2,k2)=M and 0 otherwise. Hence the set of site positions distance m from each other on any domain D{𝑿,𝑴^,𝑰} under zero flux boundary conditions can be given by

(63) Sm={(r,k),(i,j)(Π𝑿L,Π𝑿H),(Π𝑿L,Π𝑿H)|max{|r-i|,|k-j|}=m}.

Therefore

(64) cC1,C2(m)={(i,j),(r,k)Sm𝟙D1(r,k)=D2(i,j)=C, where C1=C2=C,(i,j),(r,k)Sm𝟙D1(r,k)=C1 and D2(i,j)=C2+𝟙D2(r,k)=C2 and D1(i,j)=C1 otherwise,

where C1 lies on D1, C2 lies on D2 and D1,D2{𝑿,𝑴^,𝑰}. Note that |Sm| is the number of site pairs distance m from one another. |Sm| is computed in Gavagnin et al., 2018 and is given by

(65) |Sm|=4mΠ𝑿HΠ𝑿L-3(Π𝑿H+Π𝑿L)m2+2m3.

To determine 𝔼[c^C1,C2(m)] there are three cases.

Case 1: C1=C2=C{M,X,Id,Il,Xb} lies on domain D{M,X,I} (with volume exclusion on D). For example spatial correlation of M with respect to other M. This is the case discussed in Gavagnin et al., 2018. In this case

(66) E[c^C,C(m)]=P(Two agents of type C are chosen from domain D)|Sm|.
(67) P(Two cells of type C are chosen from domain D)=NCT(NCT1)(ΠXHΠXL)(ΠXHΠXL1).

This is because there are NCT ways to choose a cell of type C and NCT-1 ways to choose a cell of type C given one has been chosen already. There are Π𝑿HΠ𝑿L possible positions for the first cell and then Π𝑿HΠ𝑿L-1 possible positions for the second cell due to volume exclusion on domain D.

Case 2: C1C2, where C1,C2{X,Id,Il,Xb} both lie on domain D{X,I} (with volume exclusion on D). For example spatial correlation of X with Xb. In this case

(68) E[c^dC1,C2(m)]=P( One cell of type C1 and one cell of type C2 are chosen from domain D)2|Sm|,

Notice, |Sm| is multiplied by two here because for each pair (i,j),(r,k)Sm there are two positions C1 and C2 can take and be distance m away from each other. Specifically these are the cases D1(i,j)=C1, D2(r,k)=C2 and D2(i,j)=C2, D1(r,k)=C1. We can compute

(69) P( One cell of type C1 and one cell of type C2are chosen on D)=NC1TNC2T(ΠXHΠXL)(ΠXHΠXL1).

This is because there are NC1T ways to choose a cell of type C1 and NC2T ways to choose a cell of type C2. There are Π𝑿HΠ𝑿L possible positions for the first cell and then Π𝑿HΠ𝑿L-1 possible positions for the second cell due to volume exclusion on domain D.

Case 3: C1, C2 where C1,C2{M,X,Id,Il,Xb} lie on domains D1, D2{M^,X,I} where D1D2, for example the spatial correlation of X and Id. In this case

(70) E[c^C1,C2(m)]=P( One cell of type C1 is chosen on D1 and one cell of type C2 is chosen on D2)2|Sm|,

where

(71) P( One cell of type C1 is chosen on D1 and one cell of type C2 is chosen on D2)=NC1TNC2T(ΠXHΠXL)2.

This is because there are NC1T ways to choose a cell of type C1 and NC2T ways to choose a cell of type C2. There are Π𝑿HΠ𝑿L possible positions for the first cell of type C1 and Π𝑿HΠ𝑿L possible positions for the second cell of type C2 since they lie on different domains (no volume exclusion). A summary of all these cases can be given as follows;

  • C1=C2=C lies on domain D (with volume exclusion on D). For example, spatial correlation of M with respect to other M.

(72) PCF(m)=cC,C(m)(Π𝑿HΠ𝑿L)(Π𝑿HΠ𝑿L-1)(4mΠ𝑿HΠ𝑿L-3(Π𝑿H+Π𝑿L)m2+2m3)(NCT(NCT-1))
  • C1C2, C1,C2 lie on domain D (with volume exclusion on D). For example spatial correlation of X with respect to Xb.

(73) PCF(m)=cC1,C2(m)(Π𝑿HΠ𝑿L)(Π𝑿HΠ𝑿L-1)2(4mΠ𝑿HΠ𝑿L-3(Π𝑿H+Π𝑿L)m2+2m3)(NC1T(NC2T))
  • C1C2, C1,C2 lie on domains D1, D2 where D1D2 for example spatial correlation of X with respect to M.

(74) PCF(m)=cC1,C2(m)(Π𝑿HΠ𝑿L)22(4mΠ𝑿HΠ𝑿L-3(Π𝑿H+Π𝑿L)m2+2m3)(NC1T(NC2T))

Appendix 4

Testing robustness

Due to the abundance of parameters and cell–cell interactions necessary to capture what is known biologically about zebrafish pigment pattern formation, it is not feasible to perform an exhaustive parameter sweep to demonstrate the robustness of the model. Instead, as a test of robustness, we perform a rigorous robustness analysis by carrying out one hundred repeats of all WT and mutant simulations with perturbed parameter values chosen uniformly at random from the range 0.75–1.25 of their described value. The value of each parameter is sampled uniformly from this region, independentally for each parameter and each repeat. Ten of these randomly sampled repeats for shd, nac, pfe and leo are given in Appendix 4—figure 1. We observe, for all one hundred repeats, that small perturbations to the rates still generate consistent patterning, demonstrating the robustness of the model.

Appendix 4—figure 1
Example shd, nac, pfe and leo simulations at stage J+ when the parameters governing the rate of proliferation, movement, differentiation and death are perturbed.

Each rectangle is an example simulation at stage J+ for which each rate parameter is perturbed to 1+x times its normal value. The value x is chosen uniformly at random from the interval [-0.25,0.25].

Appendix 5

Predicting pattern formation for an seurat/sbrcross

In this section, we give an example of how the model can be manipulated to investigate aspects of pattern development and to predict the outcome of mutant pigment patterns.

To the best of our knowledge, the adult pattern of crossed seurat and sbr mutants have not have been published previously. By changing the parameters so that we match both sbr and seurat we can predict pattern formation for a double mutant seurat/sbr. As previously stated (in more detail) in Appendix 2, the mutation in seurat affects melanocyte differentiation (Eom et al., 2012) and the mutation in sbr affects the dense-to-loose S-iridophore transition (Fadeev et al., 2015). To simulate the cross, we simply incorporate the effects of both mutations simultaneously. The results at stage J+ is given in Appendix 5—figure 1. Our model predicts that when both of these mutations occur, dense S-iridophores and associated xanthophores would cover most of the flank of the fish. A few melanocytes associated with loose S-iridophores survive at the very dorsal and ventral regions of the fish.

Appendix 5—figure 1
Prediction for the pattern of a seurat/sbr cross at stage J+.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
    P. zebrafish pigmentation mutations and the processes of neural crest development
    1. RN Kelsh
    2. M Brand
    3. YJ Jiang
    4. CP Heisenberg
    5. S Lin
    6. P Haffter
    7. J Odenthal
    8. MC Mullins
    9. FJ van Eeden
    10. M Furutani-Seiki
    11. M Granato
    12. M Hammerschmidt
    13. DA Kane
    14. RM Warga
    15. D Beuchle
    16. L Vogelsang
    17. C Nüsslein-Volhard
    (1996)
    Development 123:369–389.
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
    Nacre encodes a zebrafish Microphthalmia-Related protein that regulates Neural-Crest-Derived pigment cell fate
    1. JA Lister
    2. CP Robertson
    3. T Lepage
    4. SL Johnson
    5. DW Raible
    (1999)
    Development 126:3757–3767.
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
    An orthologue of the kit-related gene fms is required for development of neural crest-derived xanthophores and a subpopulation of adult melanocytes in the zebrafish, Danio rerio
    1. DM Parichy
    2. DG Ransom
    3. B Paw
    4. LI Zon
    5. SL Johnson
    (2000b)
    Development 127:3031–3044.
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63

Decision letter

  1. Sandeep Krishna
    Reviewing Editor; National Centre for Biological Sciences‐Tata Institute of Fundamental Research, India
  2. Aleksandra M Walczak
    Senior Editor; École Normale Supérieure, France
  3. Sreelaja Nair
    Reviewer; Tata Institute of Fundamental Research, India
  4. Raj K Ladher
    Reviewer; NCBS, India

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This paper constructs and analyses a detailed computational model of zebrafish pigment stripe formation. It highlights the role of iridophores in organizing the melanophores and xanthophores, a facet that has largely been unaccounted for in computational modeling of zebrafish stripe patterns. Of particular interest is the analysis of the schachbrett mutant, which points to the role of tight junctions in the transition of S-iridophore packing.

Decision letter after peer review:

Thank you for submitting your article "A quantitative model for zebrafish pattern formation" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Aleksandra Walczak as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Sreelaja Nair (Reviewer #1); Raj K Ladher (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

The manuscript describes a computational model that attempts to capture the complexities of xanthophore and iridophore interactions to explain zebrafish pigment stripe pattern formation. While the reviewers appreciate that the model is sufficient to reproduce known stripe patterns, it has not been demonstrated that the interactions included are necessary, and all reviewers were of the opinion that the model has not been used to make sufficiently interesting predictions. In order for a revision to be publishable, it must distinguish the work sufficiently from the work of Volkening and Sandstede by including several more predictions (see points 1,2,3 of reviewer 1, and point 3 of reviewer 3). In addition, as pointed out by reviewer 2, it is essential that the existing and added predictions be shown to be robust to changing parameter values within biologically reasonable ranges – tuning parameters one by one, keeping all others fixed, only between low and high values, is not a sufficient exploration of the parameter space to prove such robustness. A more comprehensive exploration of parameter space is also essential to satisfy reviewer 2 that the interactions included in the model are indeed necessary as claimed in the manuscript, and that despite the very large number of parameters the behaviour of the model is constrained enough to make robust predictions. The manuscript should also tone down claims of the novelty of the methods and the superiority of this approach compared to reaction-diffusion approaches – in particular, as reviewer 2 points out, stochasticity is easily included in reaction-diffusion approaches so that cannot count as a reason to favour a lattice-based model.

Reviewer #1:

In the manuscript entitled "A quantitative modeling approach to zebrafish pigment pattern formation", Owen et al. describe in silico modeling strategies to explain zebrafish pigment stripe formation. The authors use on-lattice modeling, in which five pigment cell types that are necessary to form the adult zebrafish pigment stripe pattern are treated as fixed to a mutually exclusive space within the domain in which the pattern will emerge. Onto this basic framework, the authors model known behaviors of the pigment cell types and intercellular interactions to simulate emergence of stripe pattern. Their modeling approach incorporates the role of iridophores in organizing the melanophores and xanthophores, a facet that has largely been unaccounted for in computational modeling of zebrafish pigment stripe pattern formation. The authors validate their model by recapitulating the wild type stripe pattern and test the model by simulating stripe pattern formation in the absence of one or more of the cell types/interactions. These additional simulations very nicely recreate stripe patterns observed in known zebrafish mutants. Overall the study is well carried out and the authors comprehensively validate their modeling strategy, making a strong case for on-lattice agent based modeling rather than a Reaction Diffusion modeling for pigment stripe pattern in zebrafish.

It is not surprising that the model recapitulates known pigment patterns, since the model was built based on biological evidences. On the basis of this the authors state that the current experimental evidences together with their model’s assumptions of S-iridophore behavior is sufficient to explain pigment pattern formation. For a set of biological interactions to be termed as sufficient for formation of a pattern, the model could perhaps be challenged by stating predictions that are currently not easy to test experimentally.

1) For example: As the authors point out, zebrafish do not incorporate additional pigment stripes as they grow. The biological basis for this is not known. In their wild type simulation starting at stage PB, how would a pattern evolve if initially the melanophores were in 3 or 5 stripes instead of the normal 4 stripes?

2) At the start of the simulations at stage PB, the dense iridophores appear along the horizontal myoseptum. Is this spatial location of iridophores relevant to formation of the pattern? Would the pattern that emerges be different if the original metamorphic pattern of dense iridophores at the horizontal myoseptum was displaced dorsally or ventrally?

3) The authors very nicely factor in the layered arrangement of the different cell types along the z-axis, to mirror the biological scenario. It is experimentally tough to switch the order of the layers, but should be possible to do in the simulations? This would be interesting because developmentally the three pigment cell types take distinct dorsal or ventral migratory routes along the horizontal myoseptum to reach their eventual z-layer. The route of migration has been thought to be critical in determining the fate of the pigment cell type (or of the stem cells that resides in the adult) and in establishing the stripe pattern. Could this assumption be partially tested to determine the relevance of the spatial order of the layers in eventual stripe formation?

Overall, I enjoyed the manuscript and the insights mathematical modeling can provide to understand complex phenomena such as pattern formation. My comments stem from limitations that mutants etc pose for understanding what information is sufficient to generate a pattern. The model is based on a finite set of proven interactions, several untested permutations and combinations of these and additional novel interactions are always possible. A major advantage of mathematical modeling is in the predictive zone, some of which experimentalists may find hard/impossible to venture into. The current study falls short of taking that leap, which would have been an interesting and informative exercise.

Reviewer #2:

This manuscript reports the study of pigment patterns in zebrafish embryos claiming a novel bottom-up stochastic model. My understanding is that this is an extensively employed standard approach using the chemical master equation for interacting (and diffusing) species to simulate emergent stochastic patterns. The authors incorporate many detailed and complex interactions between five different cell types and claim that all the interactions are essential to explain the observed patterns. My feeling is that such a study with a large number of control parameters is overkill and does not enhance our understanding of stripe patterns. As such, I do not recommend publishing this study in eLife.

1) In this manuscript, the authors present a quantitative model for simulating the stripe patterns seen in wt and mutant embryos of zebrafish. They develop a stochastic description based on the purported interactions between xanthophores, xanthoblasts, melanocytes and (two kinds of) iridophores while also incorporating motility that is biased by short- and long-ranged interactions between the cell types. These are augmented by cell division, death, differentiation and tissue growth. The resulting detailed numerical simulations are shown to qualitatively reproduce pigment patterns in wt and mutant conditions. Further, the authors quantitatively compare cell density, straightness of stripes and interstripe widths with experimental results. Several other results, such as the pair-correlation-function of cell densities, provide further justifications to the model developed.

2) The main point of this study seems to be the development of a bottom-up stochastic model incorporating as many complex and varied interactions between the cells to simulate pigment patterns. The authors contrast their approach with reaction-diffusion modelling approaches (such as Turing patterns) by saying that such approaches lack the stochasticity observed inherently in pigment patterns.

I would disagree with this viewpoint. The point of reaction-diffusion approaches is minimalistic coarse-grained descriptions of the complex effective interactions between cells (or any other constituents of the patterns) that can capture the essence of the emergent patterns. It is easy to incorporate stochasticity even in approaches that use partial-differential equations.

A complementary approach is to simulate the underlying master equation that governs the chemical interactions, supplemented by a compartment-based approach to diffusion. There are several such studies done in the past using the standard Gillespie algorithm for simulating the chemical master equation. Such studies incorporate the essential chemistry (e.g., activator-inhibitor dynamics) and coupled with diffusion can lead to emergent stripe patterns with inherent stochasticity built-in.

3) The authors, instead, develop a detailed lattice-based model with sufficiently complex interactions between 5 species of cells along with short-ranged local interactions and long-ranged non-local interactions. These interactions bias the hopping rates of the cells along the lattice. The authors speculate on the possible physical origins of these interactions. However, it is not clear if there aren't other possible interactions that could also lead to the same emergent dynamics. And hence the questions arises as to what are the essential features required at this level of modeling and what details are incidental. It should be noted that even though the authors claim that their approach is bottom-up, they are not modelling interactions between molecules. Rather, the interactions incorporated in their modeling are also effective cell-level interactions. In essence, this study has considered a master-equation like approach to simulate zebrafish stripe patterns.

4) The authors claim that all the interactions included in their model are essential for stripe patterns. This is demonstrated by turning off each interaction in-turn and showing that some feature of the pattern is lost. This is an 'all or none' switch. For the set of other parameters chosen, this might indeed lead to the conclusion that all interactions are absolutely needed. One could vary the interaction strengths in a continuous manner and then it is not clear if all the said interactions are absolutely needed for all parameter values.

5) Looking at the appendix and the supplementary material, I feel that the model is sufficiently complicated, and has so many turning parameters, that any range of behaviors is possible. It is not unreasonable to see that the large-scale emergent dynamics of such a complex model is essentially that of Turing systems which could have been simulated with much less complexity.

For the above reasons, I cannot recommend publishing in eLife.

Reviewer #3:

Owen et al. describe the formulation of a model to understand the formation of stripes in zebrafish. They revisit the Turing model that previous studies have described and refine the interaction parameters that were simplified in those models. The model presented recapitulates many of the patterns found in WT and in mutants, and models the formation of these stripes.

The study does capture the complexities of xanthophore and iridophore interactions well. Furthermore, the ability to place weights on the strength of interactions also gives the model flexibility.

The paper is well-written. However, I would like to make some suggestions, listed below.

1) The assumptions and the way the model works could be made more explicit. For this I would suggest that the authors consider incorporating Figure 1, 2 and 3 from Appendix 1 into the main manuscript.

2) I would suggest that the authors emphasise the differences of their model from the Volkening paper published in Aug 2018. I do take issue with the authors characterisation of this paper as very recent, and would suggest that they make greater reference to it in the Introduction.

3) One very interesting piece of data from this study, and the one that does differentiates this study from that of Volkening and Sandstede, is the analysis of the schachbrett mutant. This mutant points to the role of tight junctions in the transition of S-iridophore packing. It would be worth extending this analysis to the compunctiond sbr/leo and sbr/luc mutants described in Fadeev et al. Additionally, there is a suggestion of an interaction between sbr and seurat – I wonder if the model were able to predict what a double sbr/seurat would look like?

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for re-submitting your article "A quantitative model for zebrafish pattern formation" for consideration by eLife. Your article has been re-reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Aleksandra Walczak as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Sreelaja Nair (Reviewer #1); Raj K Ladher (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

1) In the revised manuscript, you write:

"The most commonly used mathematical paradigm for stripe formation takes the form of a Turing reaction diffusion model. In these representations, melanocytes and xanthophores diffuse and interact via a few long and short range 'reactions'. This class of model typically rely on a small number of parameters which, upon being altered, can generate a diverse range of patterns. Simplified models such as these have the benefit that they are often analytically tractable, allowing a deep understanding of the model. However, their main limitation is that, due to the simplicity of the approach, there is often no consistent way to link the parameters with measurable data, making it difficult to relate the model results back to the biology."

This gives the impression that (i) Turing pattern models are always simplified, and analytically tractable, and (ii) Turing models face a difficulty in relating their parameters to measurable data. Both are untenable statements and should be removed – there are several examples of both highly nonlinear Turing-type or reaction-diffusion models for pattern formation that are neither simple nor analytically tractable, and in several cases now people have shown how to relate their parameters to measurable data.

2) The authors write:

"Turing reaction-diffusion type models posit that combinations of short and long range dynamics between melanocytes and xanthophores generate stripe patterns. Indeed, much of the excitement around such models is the ease with which small parameter value changes result in diverse patterns, many readily recognisable from nature. […] A major difference between our model and Turing reaction diffusion models is that small parameter changes in our model do not typically generate qualitatively different patterns, whereas Turing reaction diffusion models can show substantial pattern changes in response to small changes (Budi, Patterson and Parichy, 2011; Yamaguchi, Yoshimoto and Kondo, 2007)."

This gives the impression that large changes in behaviour due to small changes in parameters is a generic feature of Turing models, as opposed to the model in this manuscript. This is not true, and the statements implying this should be removed – it is quite possible for Turing and reaction-diffusion models to produce small changes upon small changes in parameters, and equally it is possible for stochastic lattice models to exhibit large changes in behaviour upon small changes in parameters – it depends on the interactions and nonlinearities included in the model.

3) Finally:

"From an analytical perspective, a significant advantage of our on-lattice model in contrast to off-lattice models is their amenability to the derivation of a continuum model. Our model therefore opens up the opportunity for future exploratory work using a continuum model for mutants pfe and nac in order to explore whether pattern formation in these cases individually can be described as Turing patterns and to determine parameter ranges for successful pattern formation."

Here you argue, and we agree, that the model in the manuscript is in fact simply a discretized version of a reaction-diffusion model (albeit a quite complex one with many dynamical variables and many interactions). Thus, there is no great difference in approach between the discretized stochastic model you analyze and reaction-diffusion models, so please remove all statements that imply a large difference in your approach vs. reaction-diffusion or Turing type models. Further, a continuum limit can easily be constructed for off-lattice models as well, so this is not an advantage of your approach over that of Volkening et al. Please remove this statement, and present the opportunity for building continuum models as one that applies to both your model as well Volkening et al.'s model.

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

Author response

The manuscript describes a computational model that attempts to capture the complexities of xanthophore and iridophore interactions to explain zebrafish pigment stripe pattern formation. While the reviewers appreciate that the model is sufficient to reproduce known stripe patterns, it has not been demonstrated that the interactions included are necessary, and all reviewers were of the opinion that the model has not been used to make sufficiently interesting predictions. In order for a revision to be publishable, it must distinguish the work sufficiently from the work of Volkening and Sandstede by including several more predictions (see points 1,2,3 of reviewer 1, and point 3 of reviewer 3).

We have extended our work substantially to include a number of new predictions as suggested by the reviewers, including:

1) Predicting the potential effects of changing the initial conditions for stripe formation

a) by changing the initial iridophore horizontal striping to be vertical (see subsection “Initial S-iridophore interstripe orientation alone does not determine the orientation of stripes and interstripes and Figure 13A);

b) by displacing the initial horizontal interstripe vertically, so that it appears lower down the initial domain (see subsection “The position of the initial S-iridophore interstripe is important for successful pattern formation” and Figure 13B);

c) by changing the initial domain size while maintaining the usual growth rates (see subsection “Initial domain size contributes to the number of stripes and interstripes” and Figure 13C);

d) by populating the initial domain with pre-formed stripes and observing stripe insertion (see subsection “Initial domain size contributes to the number of stripes and interstripes” and Figure 13D).

2) Formulating and then testing hypotheses about the functions of the leo gene by replicating the range of mutant and cross mutant patterns (see subsection “An in silico investigation into the function of the leo gene” and Figure 14)

3) Predicting the patterning of seurat mutants and double mutants based on the known biology. These extensions are included in the subsection “Model predictions” of the main text along with our original predictions and (for seurat mutants) in the Discussion and Appendix 5.

In addition, as pointed out by reviewer 2, it is essential that the existing and added predictions be shown to be robust to changing parameter values within biologically reasonable ranges – tuning parameters one by one, keeping all others fixed, only between low and high values, is not a sufficient exploration of the parameter space to prove such robustness.

We performed a rigorous robustness analysis by carrying out one hundred repeats for simulations of each of the WT and four of the main mutant types – nac, pfe, shd and leo – with perturbed parameter values chosen uniformly at random from the range between 0.75-1.25 of their described value. The process is described in detail in subsection “Necessity of S-iridophore assumptions”. Twenty (randomly selected) of the repeats for WT are displayed in the new Figure 9 of the manuscript and ten (randomly selected) of the repeats for each of nac, pfe, shd and leo are given in Appendix 4—figure 1. For all one hundred repeats these significant perturbations to the rates still generate consistent patterning, in both WT and all mutant cases, demonstrating the robustness of the model.

A more comprehensive exploration of parameter space is also essential to satisfy reviewer 2 that the interactions included in the model are indeed necessary as claimed in the manuscript, and that despite the very large number of parameters the behaviour of the model is constrained enough to make robust predictions.

Most of the features we have included in the model have been demonstrated through biological experimentation and therefore are incorporated for completeness. A small number of the interactions included in the model are predictions focussed on the role of iridophores. We tested the necessity of these predictions by removing the corresponding model feature and demonstrating that the reduced model fails to replicate either the WT and/or mutant pattern predictions (Figure 11 of the revised manuscript). We further tested our predictions by applying our model to other mutants such as choker and sbr in order to demonstrate that these mutants are also successfully replicated by the model without the need for further assumptions (Figure 7O-Q’ of the revised manuscript). We have also carried out a more comprehensive test of the biologically suggested features demonstrating that the removal of each model feature independently leads to incorrect pattern formation in the model (Figure 12). This acts as a clear demonstration of the necessity of all the interactions included in the model.

The manuscript should also tone down claims of the novelty of the methods and the superiority of this approach compared to reaction-diffusion approaches – in particular, as reviewer 2 points out, stochasticity is easily included in reaction-diffusion approaches so that cannot count as a reason to favour a lattice-based model.

We have removed the suggestion that Turing models cannot incorporate stochasticity and toned-down claims of the superiority of the approaches in comparison to reaction-diffusion methods in the Introduction. Instead we have focussed on the important differences between the Turing reaction-diffusion modelling approach and our more biologically interpretable modelling approach.

Reviewer #1:

[…] It is not surprising that the model recapitulates known pigment patterns, since the model was built based on biological evidences. On the basis of this the authors state that the current experimental evidences together with their model`s assumptions of S-iridophore behavior is sufficient to explain pigment pattern formation. For a set of biological interactions to be termed as sufficient for formation of a pattern, the model could perhaps be challenged by stating predictions that are currently not easy to test experimentally.

1) For example: As the authors point out, zebrafish do not incorporate additional pigment stripes as they grow. The biological basis for this is not known. In their wild type simulation starting at stage PB, how would a pattern evolve if initially the melanophores were in 3 or 5 stripes instead of the normal 4 stripes?

It appears that the reviewer may have slightly misinterpreted our initial conditions. Initially, within our simulations, the melanophores are not in a striped arrangement, instead they appear dynamically across the domain, in a randomly dispersed pattern (Figure 4A’), as reported in real zebrafish (see Figure 4A). Whilst in real WT zebrafish shown in Figure 4A there are some melanocytes in a stripe-orientation that can be observed under the skin, these are much deeper below the skin layer and do not contribute to pattern development. We commented on this in the caption of Figure 4: “White arrows indicate the embryonic pattern of melanocytes in four stripes that are deeper than the skin level and are consequently not included in the model.”

Furthermore, we note that zebrafish do incorporate additional pigment stripes as they grow. However, it is important to note that these new stripes are added at the leading edges of a spreading pattern rather than between pre-formed stripes as has been documented, for instance, in the angelfish genus Pomocanthus.

Nevertheless, in order to investigate the potential for stripe insertion in zebrafish, we simulated a domain which is initially striped (two stripes and one interstripe) at stage PB. The resulting patterns demonstrate, as the reviewer suggested, that in this case new interstripes appear in-between already developed stripes, consistent with the Pomocanthus study we now reference. We have illustrated our results in Figure 13D and provided an associated commentary in the subsection “Stripe insertion can occur on an initially striped domain”.

2) At the start of the simulations at stage PB, the dense iridophores appear along the horizontal myoseptum. Is this spatial location of iridophores relevant to formation of the pattern? Would the pattern that emerges be different if the original metamorphic pattern of dense iridophores at the horizontal myoseptum was displaced dorsally or ventrally?

We have addressed the reviewer’s interesting suggestion by exploring what happens if, hypothetically, the initial iridophore stripe was in a different place. Given how we envisage the role of that initial iridophore stripe in the organisation of the melanophores, we predicted that, if we were to move the position of the initial interstripe, but maintain its horizontal orientation, that the stripes would still form sequentially around this initial interstripe, but in a pattern displaced along the dorso-ventral axis. We tested this hypothesis in the model and confirmed that this is exactly what happens. We have included simulations of pattern development from just such an initial condition in Figure 13B, discussing the similarities and differences to biologically plausible wild-type patterning in subsection “The position of the initial S-iridophore interstripe is important for successful pattern formation”. The stripes are more disrupted than with the wild-type initial condition. This effect is likely due to the altered impact of domain growth on the pattern: In the model, growth is centred at the middle of the domain and so when the initial stripe is not similarly centred, growth disrupts pattern formation.

Another interesting test case inspired by the reviewer’s comments, concerns the implication for pattern formation of altering the orientation of the initial iridophore stripe, so that it is oriented vertically (rather than horizontally as in wild-type fish). Based on the assumed role of the iridophore stripe in orienting the stripe pattern, we predicted that this change would result in vertical barring. In the simulation in Figure 13A and commentary in subsection “Initial S-iridophore interstripe orientation alone does not determine the orientation of stripes and interstripes”, it can be seen that, whilst initially the vertical iridophores drive a vertically oriented banding pattern, this gradually reorganises, so that at the end of the simulation, the pattern is a meshwork of vertical and horizontal striping. This result clearly supports the idea that the initial iridophore pattern is important in establishing stripe/interstripe orientation, but importantly reveals the previously unsuspected importance of domain growth in orienting striping in zebrafish.

3) The authors very nicely factor in the layered arrangement of the different cell types along the z-axis, to mirror the biological scenario. It is experimentally tough to switch the order of the layers, but should be possible to do in the simulations? This would be interesting because developmentally the three pigment cell types take distinct dorsal or ventral migratory routes along the horizontal myoseptum to reach their eventual z-layer. The route of migration has been thought to be critical in determining the fate of the pigment cell type (or of the stem cells that resides in the adult) and in establishing the stripe pattern. Could this assumption be partially tested to determine the relevance of the spatial order of the layers in eventual stripe formation?

The reviewer raises an interesting question: What determines the fate of pigment cells generated from the adult stem cells? To our knowledge there has been little study of the way in which the layered arrangements are built up over juvenile stages. For this reason, our model does not capture the migration pathway to the epidermis, but only the patterning process once the cells have arrived. Cells are modelled as interacting without explicit reference to layering. Thus, the order of the layers does not influence the output of our model. The distinct layers are incorporated in order to aptly represent the role of excluded volume between cells of the same type. That is, melanophores should occupy space available to other melanophores, but not to xanthophores for example, since they occupy a separate z-layer. We have added the sentence ‘We note that in our simulations, the ordering of the layers does not play any role in determining pattern formation.’ to the manuscript in the caption of Figure 2 in order to make this clear.

Reviewer #2:

This manuscript reports the study of pigment patterns in zebrafish embryos claiming a novel bottom-up stochastic model. My understanding is that this is an extensively employed standard approach using the chemical master equation for interacting (and diffusing) species to simulate emergent stochastic patterns. The authors incorporate many detailed and complex interactions between five different cell types and claim that all the interactions are essential to explain the observed patterns. My feeling is that such a study with a large number of control parameters is overkill and does not enhance our understanding of stripe patterns. As such, I do not recommend publishing this study in eLife.

We feel that the reviewer has perhaps misunderstood the motivation for our paper. As they have stated there are many demonstrations in the literature, including those we cite, which show that diverse patterns resembling those in wild-type and mutant zebrafish can be reproduced by a reaction-diffusion master equation (RDME) approach. Our aim, rather than asking if reaction-diffusion mechanisms could generate a pattern, was to ask whether the cell-cell interactions that have already been characterised biologically in the literature are necessary and sufficient to generate the observed biological patterning. Moreover, we investigate whether the resulting model can account for the diverse array of patterns seen in a variety of mutants, if the known biology corresponding to these mutants is incorporated. Though we agree that the number of parameters in this model is large, we argue that this reflects the demonstrated complexity of the underlying biology. Moreover, the iridophore interactions that are incorporated are based on hypotheses derived from biological studies; they have not been chosen for modelling convenience or to replicate particular patterns. Our model then provides a tool to determine whether these hypotheses are consistent with wild-type and mutant pattern formation.

1) In this manuscript, the authors present a quantitative model for simulating the stripe patterns seen in wt and mutant embryos of zebrafish. They develop a stochastic description based on the purported interactions between xanthophores, xanthoblasts, melanocytes and (two kinds of) iridophores while also incorporating motility that is biased by short- and long-ranged interactions between the cell types. These are augmented by cell division, death, differentiation and tissue growth. The resulting detailed numerical simulations are shown to qualitatively reproduce pigment patterns in wt and mutant conditions. Further, the authors quantitatively compare cell density, straightness of stripes and interstripe widths with experimental results. Several other results, such as the pair-correlation-function of cell densities, provide further justifications to the model developed.

2) The main point of this study seems to be the development of a bottom-up stochastic model incorporating as many complex and varied interactions between the cells to simulate pigment patterns. The authors contrast their approach with reaction-diffusion modelling approaches (such as Turing patterns) by saying that such approaches lack the stochasticity observed inherently in pigment patterns.

I would disagree with this viewpoint. The point of reaction-diffusion approaches is minimalistic coarse-grained descriptions of the complex effective interactions between cells (or any other constituents of the patterns) that can capture the essence of the emergent patterns. It is easy to incorporate stochasticity even in approaches that use partial-differential equations.

Whilst we acknowledge that Turing patterns are an excellent tool for inspiring scientists to investigate the mechanisms behind pattern formation in a variety of biological systems, we believe they are not the most appropriate tools for understanding the detailed mechanisms underlying zebrafish pigmentation patterning at the microscale.

A lack of stochasticity in classical partial differential equation models of Turing reaction-diffusion mechanisms is just one reason why we view such models as less appropriate for this inherently stochastic pattern formation process but is of relatively minor concern in comparison to the other drawbacks we perceive these models have. That being said, we acknowledge that reaction-diffusion models can incorporate stochasticity and we have now removed this comment from the manuscript.

More troubling, in our point of view, is the lack of cellular level resolution that such continuum models provide. In an age in which computational power is not an issue and cell-level biological data is available, individual-based models are an appropriate tool for investigating biological hypotheses made at a cellular level.

Even more concerning to us is the lack of detailed insight provided by coarse-grained Turing models. It is precisely the coarse-grained ‘effective’ descriptions of the complex interactions that the reviewer describes, which rob Turing models of zebrafish pigment pattern formation of their predictive power. It is often impossible, in such models, to link parameters specifically to the microscopic individual-level data, especially because many different types of reaction kinetics are capable of generating the same course-grained patterns. Consequently, Turing models are often unable to provide mechanistic insight into the underlying biological processes, something which our detailed, individual-based model is manifestly capable of doing.

Perhaps the starkest illustration of this limitation is that Turing models of zebrafish pigmentation patterns are capable of producing stripes with just two cell species – melanophores and xanthophores – whereas the real biology requires iridophores as well. Early two-species reaction-diffusion models of zebrafish pigmentation pattern formation gave no hint of the importance of the missing cell type or the necessity of the corresponding interactions, whereas detailed individual-based models such as ours, incorporating all the known biology, surely would have. That a model can produce a pattern which looks qualitatively similar to a biological pattern is not sufficient to suggest it can provide biological insight.

In summary, whilst we do not dispute that Turing reaction-diffusion models are useful for suggesting potential pattern formation mechanisms in a variety of biological contexts, and indeed were what first drew one of us, and many others, to become interested in understanding patterning processes many years ago, we do not believe them the most appropriate tool to test the complex individual-cell-based hypotheses that we investigate in this paper.

To address these comments, we have refined our comments about Turing mechanisms to the following:

“They [Turing mechanisms] typically rely on a small number of parameters which, upon being altered, can generate a diverse range of patterns. Simplified models such as these have the benefit that they are often analytically tractable, allowing a deep understanding of the pattern formation mechanism. However, their main limitation is that, due to the simplicity of the approach, there is often no consistent way to link the parameters with measurable data, making it difficult to relate the model results back to the biology.”

A complimentary approach is to simulate the underlying master equation that governs the chemical interactions, supplemented by a compartment-based approach to diffusion. There are several such studies done in the past using the standard Gillespie algorithm for simulating the chemical master equation. Such studies incorporate the essential chemistry (e.g., activator-inhibitor dynamics) and coupled with diffusion can lead to emergent stripe patterns with inherent stochasticity built-in.

Compartment-based reaction-diffusion simulations (often known as a reaction-diffusion master equation (RDME) formalism) are, as the reviewer points out, an alternative way to simulate reaction-diffusion processes which naturally incorporate stochasticity. Indeed, these models are capable of giving stripe formation consistent with the continuum models from which they are derived, but not, crucially, consistent with the biology underlying zebrafish pigmentation pattern formation for the following reasons:

1) In order to produce striped patterns, these RDME models must have multiple cells per compartment (as opposed to the single cell per compartment suggested by the volume-excluding property exhibited by the real cells). Such models typically also rely on reaction kinetics which bear little or no correspondence to the underlying biological interactions between different types of cells.

2) Once antiphase patterns of peaks and troughs have been produced by such models, an appropriately tuned threshold must be chosen in order to decide which regions constitute “stripe” and which “interstripe”. Consequently, the width of the stripes that result from such models is not an inherent property of the model, but an arbitrary choice on behalf of the modeller. The arbitrary threshold beyond which one region is considered stripe and another interstripe also belies the fact that these models result in non-trivial densities of xanthophores in the stripes and non-trivial densities of melanophores in the interstripes. This fact often goes unmentioned in such models but represents a significant departure from the real biology of the system.

Our approach shares some similarity to the RDME approach but incorporates more biological realism. Indeed, we use the Gillespie algorithm to simulate the stochastic events in our individual-based model, but rather than incorporating “the essential chemistry (e.g., activator-inhibitor dynamics)” which corresponds to a phenomenological Turing reaction-diffusion model that the reviewer suggests (NB ‘essential chemistry’ here is a reference to a required feature for the model-type, rather than an ‘essential known feature of the biological mechanism’), our detailed individual-based model begins with the biology. Consequently, we do not incur the limitations associated with Turing models that we have highlighted above.

3) The authors, instead, develop a detailed lattice-based model with sufficiently complex interactions between 5 species of cells along with short-ranged local interactions and long-ranged non-local interactions. These interactions bias the hopping rates of the cells along the lattice. The authors speculate on the possible physical origins of these interactions. However, it is not clear if there aren't other possible interactions that could also lead to the same emergent dynamics. And hence the questions arises as to what are the essential features required at this level of modeling and what details are incidental. It should be noted that even though the authors claim that their approach is bottom-up, they are not modelling interactions between molecules. Rather, the interactions incorporated in their modeling are also effective cell-level interactions. In essence, this study has considered a master-equation like approach to simulate zebrafish stripe patterns.

A bottom-up modelling approach does not require that modelling begins at the level of molecules (or, to go into even more detail, atoms), rather it should be at the lowest level at which pertinent details are available. Given the current state of biological knowledge, the level of cell-cell interactions is the lowest level at which we have sufficient information to build a model capable of recapitulating zebrafish pigmentation patterns. As a result, our model is capable of investigating biological hypotheses about interactions at this level, which higher-level, top-down models, focussed solely on pattern replication irrespective of the details of the underlying biology, would not be. Our model incorporates what is known of the biology and consequently provides a tool to investigate hypotheses about that which is unknown.

For the reasons we highlighted in the response to the reviewer’s previous point, our study is distinct from an RDME approach and avoids the biological unrealism associated with it. We consider fine-grained, biologically informed cell-level interactions as opposed to the coarse-grained interactions employed by the RDME approach.

Whilst it is possible that other cell-cell interaction dynamics might lead to similar emergent dynamics, the cell-cell interactions we have chosen are either biologically evidenced or (in the case of missing iridophore dynamics) previously hypothesised. The origins of the known cell-cell interactions are given in text [1-4]. Including alternative interactions with no basis in biological observation seems illogical.

4) The authors claim that all the interactions included in their model are essential for stripe patterns. This is demonstrated by turning off each interaction in-turn and showing that some feature of the pattern is lost. This is an 'all or none' switch. For the set of other parameters chosen, this might indeed lead to the conclusion that all interactions are absolutely needed. One could vary the interaction strengths in a continuous manner and then it is not clear if all the said interactions are absolutely needed for all parameter values.

An exhaustive continuous parameter sweep over all model parameters is not feasible. Such a parameter sweep would also be inappropriate since the majority of parameters are determined biologically. However, investigating the broad range of patterns that the model can produce is certainly a question of biological interest which we intend to follow up in a subsequent project.

5) Looking at the appendix and the supplementary material, I feel that the model is sufficiently complicated, and has so many turning parameters, that any range of behaviors is possible. It is not unreasonable to see that the large-scale emergent dynamics of such a complex model is essentially that of Turing systems which could have been simulated with much less complexity.

As we have previously stressed, the majority of the model’s parameters are biologically determined leaving the range of behaviours relatively restricted, although still capable of replicating the relevant biological patterns.

Without doubt, the basic striped patterns of wild-type zebrafish could be (and indeed has been) simulated with (simplified two-variable) Turing models. However, the point of our study is not simply to replicate patterns, but to provide mechanistic insight into the biology by which they are generated – something that Turing models are poorly equipped to do.

In order to demonstrate the robustness of our model to changes in the parameters we carried out a rigorous robustness analysis in which we varied the parameter values independently of each other. For 100 repeat simulations each parameter value was assigned at random from the range between 0.75-1.25 of their described value, effectively sampling the high-dimensional parameter space randomly. The process is described in detail in subsection “Necessity of S-iridophore assumptions” and twenty (randomly selected from the total of 100) repeats are displayed in Figure 9. It is clear to see, by looking at the figure, that the wild-type striped pattern is not significantly perturbed by any of the random parameters choice, demonstrating the robustness of the model to parameter changes. Note that we then repeated this robustness analysis for each of four of the key mutants, with similar confirmation of the robustness of the model’s outputs to reasonable parameter value changes. Twenty (randomly selected) of the repeats for WT are displayed in the new Figure 9 of the manuscript and ten (randomly selected) of the repeats for each of nac, pfe, shd and leo are given in Appendix 4—figure 1.

Reviewer #3:

[…] The paper is well-written. However, I would like to make some suggestions, listed below.

1) The assumptions and the way the model works could be made more explicit. For this I would suggest that the authors consider incorporating Figure 1, 2 and 3 from Appendix 1 into the main manuscript.

As suggested, we have added Appendix 1—figures 1, 2 and 3 to the main manuscript, reconstituting them as Figure 2 and Figure 3 to ease understanding. We thank the reviewer for this suggestion which has, we think, strengthened the manuscript by making it easier to understand.

2) I would suggest that the authors emphasise the differences of their model from the Volkening paper published in Aug 2018. I do take issue with the authors characterisation of this paper as very recent, and would suggest that the make greater reference to it in the Introduction.

We have changed ‘very recent’ to simply ‘recent’ where this appears. In the conclusion we had provided a paragraph addressing the similarities and differences between our work and Volkening’s, which we retain. To address the reviewers comment we have also added a sentence about the work of Volkening in the Introduction:

“These findings have paved the way for more detailed modelling, such as that of Volkening et al., 2018, who demonstrated (using an off-lattice individual-based model) the need for understanding S-iridophore behaviour when representing all three cell-types.”

3) One very interesting piece of data from this study, and the one that does differentiates this study from that of Volkening and Sandstede, is the analysis of the schachbrett mutant. This mutant points to the role of tight junctions in the transition of S-iridophore packing. It would be worth extending this analysis to the compunctiond sbr/leo and sbr/luc mutants described in Fadeev et al.

Upon the reviewer’s suggestion we incorporated the known features of leo, leo;sbr, leo;pfe, leo;nac and leo;shd mutants into our model. We found that, upon the incorporation of suitable assumptions we could, in fact, generate all the known patterns of these single and double mutants. We have stated the assumptions that we required in order to generate this range of distinctive patterns as a series of testable hypotheses about the function of the leo gene and have demonstrated the patterns that result when different combinations of these hypotheses are not satisfied in the model. We have included this now in the manuscript (see subsection “The position of the initial S-iridophore interstripe is important for successful pattern formation” and Figure 15). Additionally, we have tested the robustness of this pattern formation in Appendix 4—figure 1.

Additionally, there is a suggestion of an interaction between sbr and seurat – I wonder if the model were able to predict what a double sbr/seurat would look like?

The mutation in seurat affects melanocyte differentiation [3] and the mutation in sbr effects iridophore transition from dense to loose [5]. By changing the parameters to match this known biology we can effectively replicate the patterns associated with both sbr and seurat independently. Consequently, we can model and hence predict pattern formation for a double seurat/sbr mutant by changing the parameters so that the effects of both mutations occur simultaneously. We have added an example of a typical simulation result for a seurat/sbr at stage J+ in Appendix 5—figure 1. Our model predicts, when both of these mutations occur, that by stage J+, dense iridophores and xanthophores would cover most of the flank of the fish. A few melanophores associated with loose iridophores survive at the very dorsal and ventral regions of the fish. We have included a discussion of this double mutant along with the figure in in Appendix 5.

References:

1] Patterson, L. B., and Parichy, D. M. (2013). Interactions with Iridophores and the Tissue Environment Required for Patterning Melanophores and Xanthophores during Zebrafish Adult Pigment Stripe Formation. PLoS Genet., 9(5), 1-14.

2] Takahashi, G., and Kondo, S. (2008). Melanophores in the stripes of adult zebrafish do not have the nature to gather, but disperse when they have the space to move. Pigment Cell Melanoma Res., 21(6), 677–686.

3] Eom, D. S., Inoue, S., Patterson, L. B., Gordon, T. N., Slingwine, R., Kondo, S., Watanabe, M., and Parichy, D. M. (2012). Melanophore Migration and Survival during Zebrafish Adult Pigment Stripe Development Require the Immunoglobulin Superfamily Adhesion Molecule Igsf11. PLoS Genetics, 8(8), 1-16.

4] Singh, A. P., Schach, U., and Nüsslein-Volhard, C. (2014). Proliferation, dispersal and patterned aggregation of iridophores in the skin prefigure striped colouration of zebrafish. Nat. Cell Biol., 16(6), 607–614.

5] Fadeev, A., Krauss, J., Frohnhöfer, H. G., Irion, U., and Nüsslein-Volhard, C. (2015). Tightjunction protein 1a regulates pigment cell organisation during zebrafish colour patterning. eLife, 2015(4), 1–25.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

1) In the revised manuscript, you write:

"The most commonly used mathematical paradigm for stripe formation takes the form of a Turing reaction diffusion model. In these representations, melanocytes and xanthophores diffuse and interact via a few long and short range 'reactions'. This class of model typically rely on a small number of parameters which, upon being altered, can generate a diverse range of patterns. Simplified models such as these have the benefit that they are often analytically tractable, allowing a deep understanding of the model. However, their main limitation is that, due to the simplicity of the approach, there is often no consistent way to link the parameters with measurable data, making it difficult to relate the model results back to the biology."

This gives the impression that (i) Turing pattern models are always simplified, and analytically tractable, and (ii) Turing models face a difficulty in relating their parameters to measurable data. Both are untenable statements and should be removed – there are several examples of both highly nonlinear Turing-type or reaction-diffusion models for pattern formation that are neither simple nor analytically tractable, and in several cases now people have shown how to relate their parameters to measurable data.

We appreciate the reviewer’s considered points and consequently have rectified the potential ambiguity of our statements to reflect the reviewer’s concerns. We have changed the text as follows:

“The most commonly used mathematical paradigm for stripe formation takes the form of a Turing reaction diffusion model. […] However, a potential limitation is that parameters do not always have a clear biological interpretation which, can sometimes make it difficult to link parameters to measurable data.”

2) The authors write:

"Turing reaction-diffusion type models posit that combinations of short and long range dynamics between melanocytes and xanthophores generate stripe patterns. Indeed, much of the excitement around such models is the ease with which small parameter value changes result in diverse patterns, many readily recognisable from nature. […] A major difference between our model and Turing reaction diffusion models is that small parameter changes in our model do not typically generate qualitatively different patterns, whereas Turing reaction diffusion models can show substantial pattern changes in response to small changes (Budi, Patterson and Parichy, 2011; Yamaguchi, Yoshimoto and Kondo, 2007)."

This gives the impression that large changes in behaviour due to small changes in parameters is a generic feature of Turing models, as opposed to the model in this manuscript. This is not true, and the statements implying this should be removed – it is quite possible for Turing and reaction-diffusion models to produce small changes upon small changes in parameters, and equally it is possible for stochastic lattice models to exhibit large changes in behaviour upon small changes in parameters – it depends on the interactions and nonlinearities included in the model.

Whilst in some cases small parameter value changes specifically near bifurcation points can result in diverse patterns for Turing patterns (see Lee et al., 2012; Metz et al., 2011; Watanabe and Kondo, 2015)), we acknowledge that it is not a generic feature of all Turing models. We appreciate that these sentences gave a misleading impression and for this reason we have added ‘can sometimes’ to the first paragraph as follows:

"Turing reaction-diffusion-type models posit that combinations of short and long range dynamics between melanocytes and xanthophores generate stripe patterns.” Indeed, a lot of the excitement around such models is the ease with which small parameter value changes can sometimes result in diverse patterns, many readily recognisable from nature.

Furthermore we have modified the following sentence to be more specific:

“A major difference between our model and Turing reaction-diffusion models is that small parameter changes in our model do not typically generate qualitatively different patterns, whereas Turing reaction-diffusion models can show substantial pattern changes in response to small alterations to parameters near to bifurcation points (Budi, Patterson and Parichy, 2011; Yamaguchi, Yoshimoto and Kondo, 2007)."

3) Finally:

"From an analytical perspective, a significant advantage of our on-lattice model in contrast to off-lattice models is their amenability to the derivation of a continuum model. Our model therefore opens up the opportunity for future exploratory work using a continuum model for mutants pfe and nac in order to explore whether pattern formation in these cases individually can be described as Turing patterns and to

determine parameter ranges for successful pattern formation."

Here you argue, and we agree, that the model in the manuscript is in fact simply a discretized version of a reaction-diffusion model (albeit a quite complex one with many dynamical variables and many interactions). Thus, there is no great difference in approach between the discretized stochastic model you analyze and reaction-diffusion models, so please remove all statements that imply a large difference in your approach vs. reaction-diffusion or Turing type models.

We agree with the reviewer that we can likely derive a continuum reaction-diffusion model from our on-lattice individual model, but we note that this will not necessarily be a Turing reaction-diffusion model. It is for this reason that we express interest in deriving the continuum model, i.e. to assess this. Furthermore, being able to derive a continuum model from an individual-based model does not mean that the individual-based model is simply a discretisation of the resulting continuum model. Undoubtedly, in complex individual-based models which employ non-linear interactions, approximations including (but not limited to) moment closure assumptions will have to be made in order to reach the continuum limit, meaning that the continuum model will exhibit different behaviour to the individual-based model in certain parameter regimes. The advantage of deriving a continuum reaction-diffusion model from an individual-based model, rather than the other way around, is that we know exactly what assumptions have been made in passing from the individual-based model to the continuum model.

Further, a continuum limit can easily be constructed for off-lattice models as well, so this is not an advantage of your approach over that of Volkening et al. Please remove this statement, and present the opportunity for building continuum models as one that applies to both your model as well Volkening et al.'s model.

We acknowledge and address this concern as follows:

"From an analytical perspective, an advantage of our on-lattice model is its amenability to the derivation of a continuum model, although we note that continuum approximations to off-lattice individual-based models can also be derived. Our model therefore opens up the opportunity for future exploratory work using a continuum model for mutants pfe and nac in order to explore whether pattern formation in these cases individually can be described as Turing patterns and to determine parameter ranges for successful pattern formation."

That it is also possible to derive continuum limits from off-lattice models does not diminish the utility of the approach for our on-lattice models. That said, we have moved this statement, so that it fits within our summary of future work instead of our comparison with Volkening et al. as we believe its previous positioning gave a misleading impression of the purpose of the statement. Having checked through the manuscript carefully, we do not see anywhere else where we claim that our model has this advantage over the Volkening model.

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

Article and author information

Author details

  1. Jennifer P Owen

    Department of Biology and Biochemistry and Department of Mathematical Sciences, University of Bath, Claverton Down, Bath, United Kingdom
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology
    For correspondence
    jpo22@bath.ac.uk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8440-6822
  2. Robert N Kelsh

    Department of Biology and Biochemistry and Department of Mathematical Sciences, University of Bath, Claverton Down, Bath, United Kingdom
    Contribution
    Resources, Supervision, Investigation, Methodology
    Contributed equally with
    Christian A Yates
    Competing interests
    No competing interests declared
  3. Christian A Yates

    Department of Biology and Biochemistry and Department of Mathematical Sciences, University of Bath, Claverton Down, Bath, United Kingdom
    Contribution
    Supervision, Investigation, Methodology, Project administration
    Contributed equally with
    Robert N Kelsh
    For correspondence
    c.yates@bath.ac.uk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0461-7297

Funding

Biotechnology and Biological Sciences Research Council (SWBio DTP)

  • Jennifer P Owen

Biotechnology and Biological Sciences Research Council (BB/ L00769X/1 (RNK))

  • Robert N Kelsh

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

Senior Editor

  1. Aleksandra M Walczak, École Normale Supérieure, France

Reviewing Editor

  1. Sandeep Krishna, National Centre for Biological Sciences‐Tata Institute of Fundamental Research, India

Reviewers

  1. Sreelaja Nair, Tata Institute of Fundamental Research, India
  2. Raj K Ladher, NCBS, India

Publication history

  1. Received: October 23, 2019
  2. Accepted: June 21, 2020
  3. Version of Record published: July 27, 2020 (version 1)

Copyright

© 2020, Owen 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

  • 576
    Page views
  • 78
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Computational and Systems Biology
    2. Human Biology and Medicine
    Michael R Garvin et al.
    Short Report Updated

    Neither the disease mechanism nor treatments for COVID-19 are currently known. Here, we present a novel molecular mechanism for COVID-19 that provides therapeutic intervention points that can be addressed with existing FDA-approved pharmaceuticals. The entry point for the virus is ACE2, which is a component of the counteracting hypotensive axis of RAS. Bradykinin is a potent part of the vasopressor system that induces hypotension and vasodilation and is degraded by ACE and enhanced by the angiotensin1-9 produced by ACE2. Here, we perform a new analysis on gene expression data from cells in bronchoalveolar lavage fluid (BALF) from COVID-19 patients that were used to sequence the virus. Comparison with BALF from controls identifies a critical imbalance in RAS represented by decreased expression of ACE in combination with increases in ACE2, renin, angiotensin, key RAS receptors, kinogen and many kallikrein enzymes that activate it, and both bradykinin receptors. This very atypical pattern of the RAS is predicted to elevate bradykinin levels in multiple tissues and systems that will likely cause increases in vascular dilation, vascular permeability and hypotension. These bradykinin-driven outcomes explain many of the symptoms being observed in COVID-19.

    1. Computational and Systems Biology
    2. Neuroscience
    Omer Mazar, Yossi Yovel
    Research Article Updated

    For decades, researchers have speculated how echolocating bats deal with masking by conspecific calls when flying in aggregations. To date, only a few attempts have been made to mathematically quantify the probability of jamming, or its effects. We developed a comprehensive sensorimotor predator-prey simulation, modeling numerous bats foraging in proximity. We used this model to examine the effectiveness of a spectral Jamming Avoidance Response (JAR) as a solution for the masking problem. We found that foraging performance deteriorates when bats forage near conspecifics, however, applying a JAR does not improve insect sensing or capture. Because bats constantly adjust their echolocation to the performed task (even when flying alone), further shifting the signals' frequencies does not mitigate jamming. Our simulations explain how bats can hunt successfully in a group despite competition and despite potential masking. This research demonstrates the advantages of a modeling approach when examining a complex biological system.