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

Figure 2

Download asset Open asset

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. 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; $\tau$ is generated until the next possible ‘event’ occurs (based on the rates of all of the possible events). Next, a random number $u_1\in 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+\tau$. 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.

Figure 3

Download asset Open asset

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.

Figure 4

Download asset Open asset

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

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

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

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

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

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

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

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

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

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

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

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

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.

Figure 5

Download asset Open asset

Video 1

Download asset

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.

Figure 6

Download asset Open asset

The shady mutant

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

Video 2

Download asset

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

The nacre mutant

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

Video 3

Download asset

The pfeffer mutant

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

Video 4

Download asset

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

Figure 7

Download asset Open asset

Simulation of other mutants

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

The rse mutant

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

Video 5

Download asset

The sbr mutant

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

Video 6

Download asset

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

The cho mutant

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

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

Video 7

Download asset

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.

Figure 8

Download asset Open asset

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.

Figure 9

Download asset Open asset

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

Melanophore density across WT and mutant fish

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

Figure 10

Download asset Open asset

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.

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.

Figure 13

Download asset Open asset

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

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

Initial domain size contributes to the number of stripes and interstripes

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

Stripe insertion can occur on an initially striped domain

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

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

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

Figure 14

Download asset Open asset

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

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

An in silico investigation into the function of the leo gene

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

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

Figure 15

Download asset Open asset

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

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

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

• Hypothesis 1:Melanocytes are not repelled by xanthophores.

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

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

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

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

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

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

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

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

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

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

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

Necessity of the pre-existing hypothesis about leo

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

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

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

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

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

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

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

Robustness of the leo assumptions in generating spots

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

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

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

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

where ${\mathrm{\Pi }}_{\mathit{𝑿}}^H={\mathrm{\Pi }}_{\mathit{𝑰}}^H=6$, ${\mathrm{\Pi }}_{\mathit{𝑿}}^L={\mathrm{\Pi }}_{\mathit{𝑰}}^L=12$, ${\mathrm{\Pi }}_{\mathit{𝑴}}^H=3$, ${\mathrm{\Pi }}_{\mathit{𝑴}}^L=6$.

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

where $D\in \{\mathit{𝑿},\mathit{𝑴},\mathit{𝑰}\}$ and ${\mathrm{\Pi }}_D^H(t)$ is the number of sites in the $y$ direction of domain type $D$ at time $t$. The simulated height ${\mathrm{\Omega }}_H(t)={\mathrm{\Delta }}_D\times {\mathrm{\Pi }}_D^H(t)$ directly corresponds to the height at the anterior margin of the anal fin (abbreviated as HAA) of the fish at all times.

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 ${\displaystyle {\mathrm{\Omega }}_H(0)=1mm}$ and ${\displaystyle {\mathrm{\Omega }}_L(0)=2mm}$. Therefore, ${\mathrm{\Pi }}_{\mathit{𝑴}}^L(0)$=25, ${\mathrm{\Pi }}_{\mathit{𝑴}}^L(0)$=50, ${\mathrm{\Pi }}_X^H(0)={\mathrm{\Pi }}_I^H(0)=50$ and ${\mathrm{\Pi }}_{\mathit{𝑿}}^L(0)={\mathrm{\Pi }}_{\mathit{𝑰}}^L(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 $\mathit{𝑰}$ 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_M^T(0)=50M$ uniformly at random on $\mathit{𝑴}$ so that the initial occupancy density of melanocytes is

where $N_C^T(t)$ for any cell type $C\in \{M,X,X^b,I^d,I^l\}$ 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)=2000X^b$ uniformly at random on $\mathit{𝑿}$ 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\in \{X,I^l\}$, $N_C^T(0)$ = 0.

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 $\mathit{𝑿}$ with xanthophores ($X$) at stage SP that is sites $\mathit{𝑿}(i,j,t_{SP})=X$ for all $i\in \lceil \frac{{\mathrm{\Pi }}_X^H(t_{SP})}{2}\rceil -1:\lceil \frac{{\mathrm{\Pi }}_X^H(t_{SP})}{2}\rceil +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 ${\alpha }_0={\sum }_{i=1}^{15}{\alpha }_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\sim U(0,1)$). Next, we generate $u_1\sim U(0,1)$ to determine which event occurs at time $t+\tau$. 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+\tau$. 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 ${\displaystyle {\mathrm{\Omega }}_{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.

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\in \{\mathit{𝑿},\mathit{𝑰},\mathit{𝑴}\}$ 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=\mathit{𝑴},\mathit{𝑿}$ in Figure 3A–H in the main text. Formulae for these sets are given in Supplementary file 2.

To illustrate this see Figure 3A–H in the main text which compares the number of cells in the short (Figure 3A–D) and long (Figure 3E–H) range distance from a central site, $D(r,k)$, on different domain types. In this figure, $M$ are represented as black circles. $X$ are represented as yellow circles. Figure 2A is a visualisation of sites (marked in red) in $S_{\mathit{𝑴}}(r,k)$, where $\mathit{𝑴}(r,k)$ is the central site marked in grey. $S_{\mathit{𝑴}}^M(r,k,t)$ are the sites marked in red in which a melanocyte resides. The number of melanocytes in the short range distance from $\mathit{𝑴}(r,k)$ at time $t$ is given by $N_{M,M}^S(r,k,t)=S_{\mathit{𝑴}}^M(r,k,t)=2$. Figure 2B is a visualisation of the sites considered when calculating $N_{M,X}^S$ (formula given in Supplementary file 2). In this example, $N_{M,X}^S(r,k,t)=9$. To compare the number of $M$ and $X$ in the short range distance of $\mathit{𝑴}(r,k)$, we consider $w{N}_{M,M}^S(r,k,t)=4\times 2=8$ which corresponds to the weighted value of melanocytes in the short range distance with $N_{M,X}^S(r,k,t)=9$ the number of xanthophores in the short range. Figure 3C is a visualisation of sites (marked in red) in $S_{\mathit{𝑿}}(r,k)$, where $\mathit{𝑿}(r,k)$ is the central site marked in grey. $S_{\mathit{𝑿}}^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 $\mathit{𝑿}(r,k)$ at time $t$ is given by $N_{X,X}^S(r,k,t)=S_{\mathit{𝑿}}^X(r,k,t)=7$. Figure 2D is a visualisation of the sites considered when calculating $N_{X,M}^S(r,k,t)$ (formula given in Supplementary file 2). In this example, $N_{X,M}^S(r,k,t)=14$. To compare the number of $M$ and $X$ in the short range distance of $\mathit{𝑿}(r,k)$, at time $t$, we would compare $N_{X,M}^S(r,k,t)$=14 with $N_{X,X}^S(r,k,t)=7$. Figure 3E is a visualisation of sites (marked in red) in $L_{\mathit{𝑴}}(r,k)$, where $\mathit{𝑴}(r,k)$ is the central site marked in grey. $L_{\mathit{𝑴}}^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 $\mathit{𝑴}(r,k)$ is given by $N_{M,M}^L(r,k,t)=L_{\mathit{𝑴}}^M(r,k,t)=20$. Figure 2F is a visualisation of the sites considered when calculating $N_{M,X}^L(r,k,t)$ (formula given in Supplementary file 2). In this example, $N_{M,X}^L(r,k,t)=16$. To compare the number of $M$ and $X$ in the long range distance of $\mathit{𝑴}(r,k)$, we consider $w{N}_{M,M}^L(r,k,t)=4\times 20=80$ which corresponds to the weighted value of melanocytes in the long range distance with $N_{M,X}^S(r,k,t)=16$, the number of xanthophores in the long range. Figure 2G is a visualisation of sites (marked in red) in $L_{\mathit{𝑿}}(r,k)$, where $\mathit{𝑿}(r,k)$ is the central site marked in grey. $L_{\mathit{𝑿}}^X(r,k)$ are the sites marked in red in which a xanthophore resides. The number of xanthophores in the long range distance from $\mathit{𝑿}(r,k)$ is given by $N_{X,X}^L=L_{\mathit{𝑿}}^X(r,k,t)=2$. Figure 3H is a visualisation of the sites considered when calculating $N_{X,M}^L(r,k,t)$ (formula given in Supplementary file 2). In this example $N_{X,M}^L(r,k,t)=20$. To compare the number of $M$ and $X$ in the long range distance of $\mathit{𝑿}(r,k)$, we consider $N_{X,M}^L(r,k,t)=20$ with $N_{X,X}^L(r,k,t)=2$.

We denote $S_D^C(r,k,t)\subset S_D(r,k)$ where $C$ lies in layer $D$ to be the sites in $S_D(r,k)$ occupied by cell type $C\in \{M,X,I^d,I^l,X^b\}$ at time $t$. We denote $N_{C_i,C_j}^S(r,k,t)$, ($N_{C_i,C_j}^L(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_{C_i,C_i}^S(r,k,t)=|S_D^{C_i}(r,k,t)|$. The formulae for $N_{C_i,C_j}^S(r,k,t)$, $N_{C_i,C_j}^L(r,k,t)$ where $C_i$ does not equal $C_j$ are more complicated and are given in Equations 7, 8 (and Supplementary file 2 for reference).

Simply, where domain types $D$ do not have the same lattice site size ${\mathrm{\Delta }}_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 $\mathit{𝑿}(i,j)$, corresponds to a quarter of site $\mathit{𝑴}(\lceil \frac{i}{2}\rceil ,\lceil \frac{j}{2}\rceil )$, this explains case two in Equation 7. Similarly each site $\mathit{𝑴}(i,j)$ corresponds to four sites on $\mathit{𝑿}$, specifically $\mathit{𝑿}(2i,2j)$, $\mathit{𝑿}(2i,2j-1)$, $\mathit{𝑿}(2i-1,2j)$ and $\mathit{𝑿}(2i-1,2j-1)$. This explains case three in Equations 7 and 8. Note that since $M$ is four times larger than all other cells in our simulation, we provide a weighting system when comparing $M$ with $C\mathrm{\backslash }M$ in some cases.

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.

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.