A quantitative modelling approach to zebrafish pigment pattern formation

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 cellular composition of the stripes and how these become assembled has been welldescribed. 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., 2003; Figure 1A).
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 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). and xanthophores. Similarly, adult light interstripes are made up of layers of dense S-iridophores and xanthophores ( Figure 1B; Hirata et al., 2003). The final striped pattern is generated by the selforganisation 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 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.

Modelling overview
We build our model with direct reference to the known biology. We model five cell types as individual agents: melanocytes (M), xanthophores (X), xanthoblasts (X b ), the unpigmented precursor cell to xanthophores) and S-iridophores in either dense or loose form (I d , I l , 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).
Agents representing X and X b , M, I d and I l 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. The domain is made up of three layers. Layer X which contains yellow X (yellow circles) and unpigmented X b (clear circles with black outline). Layer I which contains silvery I d (white circles) and blue I l (blue circles). Layer M 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 a 1 ðtÞ; :::a 15 ðtÞ are calculated. (5) Random numbers u 0 ; u 1 2 Uð0; 1Þ are generated. (6) Numbers u 0 ; u 1 are used to determine the time t 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 W SL ¼ 13:5mm. If so, the algorithm finishes. Otherwise the algorithm continues, returning to stage 2. (9) The simulation completes.
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; t is generated until the next possible 'event' occurs (based on the rates of all of the possible events). Next, a random number u 1 2 Uð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 þ 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.

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.
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 (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). 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 dedifferentiating 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 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
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 3; Parichy 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.

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, I l , M) and light interstripes (I d , 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.

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.

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 stripelike 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.
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 Video 1. Simulated development of WT. https://elifesciences.org/articles/52998#video1 Video 2. Simulated development of shd. https://elifesciences.org/articles/52998#video2 Video 3. Simulated development of nac. https://elifesciences.org/articles/52998#video3 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.

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. 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').

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  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 . 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. 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.

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.
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.

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.
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.

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. 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.

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.
The value of SS(x) lies between 0 and 1, since C L. 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 Figure 11. 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. at stage J+ for fish lacking one or more S-iridophore transition mechanisms display patterns which diverge from those seen in real fish.
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-todense interactions are important for generating subsequent interstripes. Interestingly, without looseto-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  Figure 12 continued on next page 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).
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.

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. Figure 13. 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.
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).

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) 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.
Pattern formation is also altered in the double mutants leo;shd, leo;nac and leo;pfe when compared with shd, nac and pfe . 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 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 .
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 . We will refer to this as hypothesis one for the effects of the mutant leo.
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.
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., 2013; Figure 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 . 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 . 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 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 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 continuoustime 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 dorsoventral 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. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. Additional files

Supplementary files
. Supplementary file 1. Summary of all notation used in Appendix 1.
. Supplementary file 2. Summary of all notation used in Appendix 1 regarding short and long range interactions. . Supplementary file 6. Parameters implemented in the model.
. Transparent reporting form

Data availability
All mathematical modelling assumptions/ methods have been provided in the supplementary material. Code relating to this paper have been made available on github, at https://github.com/ JenniferOwen/Zebrafish-stripe-model (copy archived at https://github.com/elifesciences-publications/Zebrafish-stripe-model).
shorter than the real width (SL), the corresponding simulated SL for any simulation at time t can be calculated by: where D 2 fX; M; Ig and P H D ðtÞ is the number of sites in the y direction of domain type D at time t. The simulated height W H ðtÞ ¼ D D Â P H D ð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 W H ð0Þ ¼ 1mm and W L ð0Þ ¼ 2mm. Therefore, P L M ð0Þ=25, P L M ð0Þ=50, P H X ð0Þ ¼ P H I ð0Þ ¼ 50 and P L X ð0Þ ¼ P L I ð0Þ ¼ 100. The initial occupancies for a WT simulation comprise three rows of I d along the middle of the fish, that is, sites on rows 24 to 26 on I are fully occupied with I d cells at time t ¼ 0. This represents the S-iridophores which differentiate along the horizontal myoseptum between 21 and 25 dpf. We also place N T M ð0Þ ¼ 50 uniformly at random on M so that the initial occupancy density of melanocytes is where N T C ðtÞ for any cell type C 2 fM; X; X b ; I d ; I l g 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 N X b ð0Þ ¼ 2000 uniformly at random on X so 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 2 fX; I l g, N T C ð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 X with xanthophores (X) at stage SP that is sites Xði; j; t SP Þ ¼ X for all i 2 d P H X ðtSPÞ 2 e À 1 : d 2 e þ 1 where t SP 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 is generated until the next continuous 'event' occurs where a 0 ¼ P 15 i¼1 a i ðtÞ is the sum of all of the event propensities given in Supplementary file 4 and u 0 is a uniformly distributed random number in ð0; 1Þ (i.e. u 0~U ð0; 1Þ). Next, we generate u 1~U ð0; 1Þ to determine which event occurs at time t þ t . The event i that satisfies 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 þ 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 W 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 S D ðr; kÞ and L D ðr; kÞ where D 2 fX; I; Mg 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 S D ðr; kÞ and L D ðr; kÞ are both different for different domain types due to the different lattice site sizes. S D ðr; kÞ and L D ðr; kÞ are visualised for D ¼ M; X in Figure 3A-H in the main text. Formulae for these sets are given in Supplementary file 2.
To illustrate this see Figure 3A-  Figure 2B is a visualisation of the sites considered when calculating N S M;X (formula given in Supplementary file 2). In this example, N S M;X ðr; k; tÞ ¼ 9. To compare the number of M and X in the short range distance of Mðr; kÞ, we consider wN S M;M ðr; k; tÞ ¼ 4 Â 2 ¼ 8 which corresponds to the weighted value of melanocytes in the short range distance with N S M;X ðr; k; tÞ ¼ 9 the number of xanthophores in the short range. Figure 3C is a visualisation of sites (marked in red) in S X ðr; kÞ, where Xðr; kÞ is the central site marked in grey. S X 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 Xðr; kÞ at time t is given by N S X;X ðr; k; tÞ ¼ S X X ðr; k; tÞ ¼ 7. Figure 2D is a visualisation of the sites considered when calculating N S X;M ðr; k; tÞ (formula given in Supplementary file 2). In this example, N S X;M ðr; k; tÞ ¼ 14. To compare the number of M and X in the short range distance of Xðr; kÞ, at time t, we would compare N S X;M ðr; k; tÞ=14 with N S X;X ðr; k; tÞ ¼ 7. Figure 3E is a visualisation of sites (marked in red) in L M ðr; kÞ, where Mðr; kÞ is the central site marked in grey. L M 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 Mðr; kÞ is given by N L M;M ðr; k; tÞ ¼ L M M ðr; k; tÞ ¼ 20. Figure 2F is a visualisation of the sites considered when calculating N L M;X ðr; k; tÞ (formula given in Supplementary file 2). In this example, N L M;X ðr; k; tÞ ¼ 16. To compare the number of M and X in the long range distance of Mðr; kÞ, we consider wN L M;M ðr; k; tÞ ¼ 4 Â 20 ¼ 80 which corresponds to the weighted value of melanocytes in the long range distance with N S M;X ðr; k; tÞ ¼ 16, the number of xanthophores in the long range. Figure 2G is a visualisation of sites (marked in red) in L X ðr; kÞ, where Xðr; kÞ is the central site marked in grey. L X X ðr; kÞ are the sites marked in red in which a xanthophore resides. The number of xanthophores in the long range distance from Xðr; kÞ is given by N L X;X ¼ L X X ðr; k; tÞ ¼ 2. Figure 3H is a visualisation of the sites considered when calculating N L X;M ðr; k; tÞ (formula given in Supplementary file 2). In this example N L X;M ðr; k; tÞ ¼ 20. To compare the number of M and X in the long range distance of Xðr; kÞ, we consider N L X;M ðr; k; tÞ ¼ 20 with N L X;X ðr; k; tÞ ¼ 2. We denote S C D ðr; k; tÞ & S D ðr; kÞ where C lies in layer D to be the sites in S D ðr; kÞ occupied by cell type C 2 fM; X; I d ; I l ; X b g at time t. We denote N S Ci;Cj ðr; k; tÞ, (N L Ci;Cj ðr; k; tÞ) as the number of cells of type C j in the short (long) range distance 0.04 mm, (0.24 mm) from cell type C i at D i ðr; kÞ. Hence the number of cells of type C i in the short distance from another cell of the same type at Dðr; kÞ at time t is given by N S Ci;Ci ðr; k; tÞ ¼ jS Ci D ðr; k; tÞj. The formulae for N S Ci;Cj ðr; k; tÞ, N L Ci;Cj ðr; k; tÞ where C i does not equal C j are more complicated and are given in N L C1;C2 ðr; k; tÞ ¼ L C2 D2 ðr; k; tÞif D 1 ¼ D 2 or D 1 andD 2 2 fX; Ig: P ð Þ¼C2 þ ::: Simply, where domain types D do not have the same lattice site size D 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 Xði; jÞ, corresponds to a quarter of site Mðd i 2 e; d j 2 eÞ, this explains case two in Equation 7. Similarly each site Mði; jÞ corresponds to four sites on X, specifically Xð2i; 2jÞ, Xð2i; 2j À 1Þ, Xð2i À 1; 2jÞ and Xð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 CM 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 2 O ¼ fNo; So; Ea; We; NE; SE; SW; NWg. (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 where A (R) 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 A and R are defined as follows.
Finally we normalise P O s.t.; 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Þ. S-iridophore, xanthophore and xanthoblast proliferation (continuous time events 6-7) Given the event is chosen such that a cell of type C 2 fX; X b ; I d ; I l g is determined to proliferate we choose a random cell of type C whose position is given by Dðr; kÞ, D 2 fX; Ig. Next a random number u~Uð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 2 ½0; 1=4Þ, Dðr þ 1; kÞ if u 2 ½1=4; 1=2Þ, Dðr À 1; kÞ if u 2 ½1=2; 3=4Þ or Dðr; k À 1Þ if u 2 ½3=4; 1.
For cell types C 2 fI d ; I l ; X b g, 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 N S X;X ðr; k; tÞ þ N S X;I d ðr; k; tÞ>N S X;I l ðr; k; tÞ þ N S X;M ð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 Mðr; kÞ is chosen at random. A differentiation event is successful and a new M appears in this site if the site Mðr; kÞ is empty and the following is true: where a ¼ 2: 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 u 2 Uð0; 1Þ and we find that u<0:01 in combination with the condition N S M;X ðr; k; tÞ þ N S M;I d ðr; k; tÞ þ N S M;M ð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 X b is chosen at random from X. Suppose the chosen X b lies in site Xðr; kÞ. The differentiation event is successful and a X replaces the X b in this site if the site Mðd r 2 e; d k 2 eÞ 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. N S X;X ðr; k; tÞ þ N S X;I d ðr; k; tÞ>N S X;M ðr; k; tÞ þ N S X;I l ð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).
Alternatively the differentiation event is successful if: N S X;X ðr; k; tÞ þ N S X;M ðr; k; tÞ ¼ 0; holds and a randomly generated number u 2 Uð0; 1Þ is such that u<0:01, or Equation (42) where L x ; S x ; S m =12, 9, one respectively. Alternatively if an I l to I d shape transition is chosen, then an I l is chosen at random from domain I and evaluated for meeting transition criteria. This transition is successful and the chosen I l in position Iðr; kÞ is replaced with I d in this site if either N L I l ;X ðr; k; tÞ<L x2 and N S I l ;M ðr; k; tÞ<S m2 ; or alternatively N S I l ;X ðr; k; tÞ>S x2 and N S I l ;M ðr; k; tÞ<S m2 ; where L x2 ¼ 16, S x2 ¼ 4 and S m2 ¼ 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 L x , S x , S m , L x2 , S m2 , S x2 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.
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 u 2 Uð0; 1Þ, is such that u<0:001 and both of the following two cases hold: where 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, I l 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 X b uniformly at random from all X b occupying X that meets the following criterion: suppose the chosen X b resides in site Xðr; kÞ, then the corresponding site on M, Mðd r 2 e; d k 2 eÞ must be empty. We simulate this by randomly choosing X b on X, and checking whether Mðd r 2 e; d k 2 eÞ is empty. If so, we continue to the next stage. Otherwise we repeat through all X b on X until we either find a suitable X b . If no X b satisfy this criteria at time t then we abort the cell-cell pulling event. Given we find a suitable X b , the chosen X b 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 Mðd r 2 e; d k 2 eÞ. 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; three central rows. In rse we set N I d ð0Þ=60, and we place these I d 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 2 Z that is uniformly distributed between 24 and 26 (the centre of the horizontal axis is at 25), and a random number k 1 that is uniformly distributed between 0 and P L I ð0Þ ¼ 100. If Iðr 1 ; k 1 Þ is empty we place a I d cell in this site, otherwise we generate a new position until there are 60 cells of type I d on I. Furthermore since new S-iridophores are only produced by proliferation of preexisting S-iridophores, we also reduce the rate of proliferation of I d and I l to one fifth of the usual number. Therefore, we set the propensity of event 7, a 7 ðtÞ ¼ 0: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 I d attempts transition to I l to a fortieth of its usual value. Hence the propensity of event 10 becomes a 10 ð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 I d so N I d ð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 I at stage PR. We do this by generating a random number r between 0 and P L I ðt PR Þ and a random number k between 0 and P H I ðt PR Þ, where t PR is the time when the simulation first enters stage PR. If Iðr; kÞ is empty we place a I l cell in this site, otherwise we generate a new r, k and continue until there are 200 cells of type I l on I.

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 a 8 ðtÞ ¼ 0:05 Â N I d ð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 a 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 r mx , 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 N I ; X S is changed to five and the number N I ; X L 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.
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.

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 ofX andÎ of zeros and ones is generated such that a matrix entryXði; jÞ ¼ 1 if Xði; jÞ ¼ X and 0 otherwise. SimilarlyÎði; jÞ ¼ 1 if Iði; jÞ ¼ I d and 0 otherwise. We then generate l t and l b using eitherX orÎ as described in Section "Necessity of S-iridophore assumptions". From these values, we generate IS (X0 interstripe width) by computing; IS ¼ X P L X i¼1 l t ðiÞ À l b ðiÞ P L X : Note that the X0 interstripe width is computed using whichever is more appropriate of X on X and I d on I given the mutation. For example, pfe does not have xanthophores, so we would use the distribution of I d on I 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: where c C1;C2 ðmÞ is defined as the number of cells of type C 1 we would expect to find at distance m from cells of type C 2 using the uniform metric under zero flux boundary conditions. E½ĉ 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 M we must transform M into a matrixM of size P H X Â P L X whereMðr; kÞ ¼ M if Mðd r 2 e; d k 2 eÞ ¼ M and 0 otherwise. Hence the set of site positions distance m from each other on any domain D 2 fX;M; Ig under zero flux boundary conditions can be given by S m ¼ fðr; kÞ; ði; jÞ 2 ðP L X ; P H X Þ; ðP L X ; P H X Þ j maxfjr À ij; jk À jjg ¼ mg: Therefore c C1;C2 ðmÞ ¼ P ði;jÞ;ðr;kÞ2Sm 1 D1ðr;kÞ¼D2ði;jÞ¼C ; whereC 1 ¼ C 2 ¼ C; P ði;jÞ;ðr;kÞ2Sm 1 D1ðr;kÞ¼C1andD2ði;jÞ¼C2 þ 1 D2ðr;kÞ¼C2andD1ði;jÞ¼C1 otherwise; ( where C 1 lies on D 1 , C 2 lies on D 2 and D 1 ; D 2 2 fX;M; Ig. Note that jS m j is the number of site pairs distance m from one another. jS m j is computed in Gavagnin et al., 2018 and is given by To determine E½ĉ C1;C2 ðmÞ there are three cases. Case 1: C 1 ¼ C 2 ¼ C 2 fM; X; I d ; I l ; X b g lies on domain D 2 fM; X; Ig (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 E½ĉ C;C ðmÞ ¼ PðTwo agents of type C are chosen from domain DÞjS m j: PðTwo cells of type C are chosen from domain DÞ ¼ N T C ðN T C À 1Þ ðP H X P L X ÞðP H X P L X À 1Þ : This is because there are N T C ways to choose a cell of type C and N T C À 1 ways to choose a cell of type C given one has been chosen already. There are P H X P L X possible positions for the first cell and then P H X P L X À 1 possible positions for the second cell due to volume exclusion on domain D.
Case 2: C 1 6 ¼ C 2 , where C 1 ; C 2 2 fX; I d ; I l ; X b g both lie on domain D 2 fX; Ig (with volume exclusion on D). For example spatial correlation of X with X b . In this case E½ĉ C1;C2 d ðmÞ ¼ Pð One cell of type C 1 and one cell of type C 2 are chosen from domain DÞ2jS m j; (68) Notice, jS m j is multiplied by two here because for each pair ði; jÞ; ðr; kÞ 2 S m there are two positions C 1 and C 2 can take and be distance m away from each other. Specifically these are the cases D 1 ði; jÞ ¼ C 1 , D 2 ðr; kÞ ¼ C 2 and D 2 ði; jÞ ¼ C 2 , D 1 ðr; kÞ ¼ C 1 . We can compute Pð One cell of type C 1 and one cell of type C 2 are chosen on DÞ ¼ N T C1 N T

C2
ðP H X P L X ÞðP H X P L X À 1Þ : (69) This is because there are N T C1 ways to choose a cell of type C 1 and N T C2 ways to choose a cell of type C 2 . There are P H X P L X possible positions for the first cell and then P H X P L X À 1 possible positions for the second cell due to volume exclusion on domain D.
Case 3: C 1 , C 2 where C 1 ; C 2 2 fM; X; I d ; I l ; X b g lie on domains D 1 , D 2 2 fM; X; Ig where D 1 6 ¼ D 2 , for example the spatial correlation of X and I d . In this case E½ĉ C1;C2 ðmÞ ¼ Pð One cell of type C 1 is chosen on D 1 and one cell of type C 2 is chosen on D 2 Þ2jS m j; where Pð One cell of type C 1 is chosen on D 1 and one cell of type C 2 is chosen on D 2 Þ ¼ N T C1 N T

C2
ðP H X P L X Þ 2 : (71) This is because there are N T C1 ways to choose a cell of type C 1 and N T C2 ways to choose a cell of type C 2 . There are P H X P L X possible positions for the first cell of type C 1 and P H X P L X possible positions for the second cell of type C 2 since they lie on different domains (no volume exclusion). A summary of all these cases can be given as follows; . C 1 ¼ C 2 ¼ C lies on domain D (with volume exclusion on D). For example, spatial correlation of M with respect to other M.