A quantitative modelling approach to zebrafish pigment pattern formation
Abstract
Pattern formation is a key aspect of development. Adult zebrafish exhibit a striking striped pattern generated through the selforganisation of three different chromatophores. Numerous investigations have revealed a multitude of individual cellcell interactions important for this selforganisation, but it has remained unclear whether these known biological rules were sufficient to explain pattern formation. To test this, we present an individualbased mathematical model incorporating all the important celltypes 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 ruleset is sufficient to recapitulate wildtype 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 selforganisation 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 selflabelling nature of pigment cells.
The horizontal blue and gold stripes of zebrafish are now one of the beststudied examples of pigment pattern formation, especially at the level of underlying cellular mechanisms (Singh and NüssleinVolhard, 2015; Watanabe and Kondo, 2015; Patterson and Parichy, 2019). Zebrafish are amenable to observational studies, since all development takes place outside the mother and the skin is transparent. This, combined with the availability of multiple key mutants (affecting, for example, celltype differentiation and patterning), and the development of innovative in vivo cell ablation and in vitro cell culture techniques, have provided a unique opportunity to investigate the cellular and molecular basis for pigment pattern formation experimentally (Eom et al., 2015; Budi et al., 2011; Ceinos et al., 2015; Yamanaka and Kondo, 2014; Eom and Parichy, 2017; Hamada et al., 2014; Fadeev et al., 2015; Inoue et al., 2014; Irion et al., 2014; Watanabe et al., 2006; Frohnhöfer et al., 2013; Parichy et al., 2009; Svetic et al., 2007; Mellgren and Johnson, 2006; Parichy et al., 2000b; Iwashita et al., 2006; Hirata et al., 2005; Kelsh et al., 1996; Lister et al., 1999; Parichy et al., 2000a; Maderspacher and NüssleinVolhard, 2003; Walderich et al., 2016; Patterson et al., 2014; McMenamin et al., 2014; Patterson and Parichy, 2013; Krauss et al., 2014; Parichy and Turner, 2003; Eom et al., 2012; Mahalwar et al., 2016; Mahalwar et al., 2014; Asai et al., 1999; Takahashi and Kondo, 2008).
The cellular composition of the stripes and how these become assembled has been 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 pigmentproducing 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 Liridophores, and type Siridophores. Only type Siridophores play a role in stripe formation (type Liridophores appear later and are likely involved in pattern maintenance [Hirata et al., 2003]), and they appear in two different forms. In the light interstripes, Siridophores appear in a ‘dense’ arrangement (dense Siridophores), forming a continuous sheet, whilst in the dark stripes the cells are in a ‘loose’ arrangement and appear more widely spaced (loose Siridophores) (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 Liridophores, melanocytes, loose Siridophores and xanthophores. Similarly, adult light interstripes are made up of layers of dense Siridophores and xanthophores (Figure 1B; Hirata et al., 2003). The final striped pattern is generated by the selforganisation of xanthophores, melanocytes, loose and dense Siridophores 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 postfertilisation (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; Siridophores 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 Siridophores 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 dedifferentiation 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 Siridophores alongside the horizontal myoseptum. In response to the appearance of these Siridophores, 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, Siridophores proliferate rapidly and spread bidirectionally; at the edges of the interstripes they switch to a more scattered (less tightlypacked) form as they continue to spread dorsally and ventrally. Spreading loose Siridophores transition back into dense Siridophores at the locations of the future interstripes X1V and X1D (Figure 1C, stages PR to SP) (Mahalwar et al., 2014). Once Siridophores aggregate at the next interstripe, the process starts again, that is, xanthophores differentiate in response to the dense Siridophores and melanocytes accumulate either side of the new interstripe generating the subsequent stripe. This process of Siridophore aggregation predetermining future interstripe locations and subsequent delamination in future stripe regions repeats until Siridophores 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), longdistance communication, through stabilisation of elongated cellular projections (Eom and Parichy, 2017; Eom et al., 2015), as well as the shape transitions of Siridophores (Fadeev et al., 2015). During this period, there is also simultaneous twodimensional domain growth (Parichy et al., 2009). The pattern is formed by cell–cell interactions of all three pigment producing cell types: melanocytes, xanthophores and Siridophores. 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 Siridophores. The most commonly used mathematical paradigm for stripe formation takes the form of a Turing reactiondiffusion model. In these representations, melanocytes and xanthophores diffuse and interact via a few long and shortrange ‘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 Siridophores (Watanabe and Kondo, 2015; Kondo, 2017; Painter et al., 2015; Bloomfield et al., 2011; Binder and Simpson, 2013; Volkening and Sandstede, 2015; Kondo, 2017; Nakamasu et al., 2009; Moreira and Deutsch, 2005; Bullara and De Decker, 2015; Yamaguchi et al., 2007; Asai et al., 1999). They suggest that the role for iridophores is restricted to simply orienting stripes (Volkening and Sandstede, 2015; Nakamasu et al., 2009; Binder and Simpson, 2013). New biological observations demonstrate that Siridophores play a fundamental role in body stripe formation (Singh and NüssleinVolhard, 2015; Frohnhöfer et al., 2013; Patterson and Parichy, 2013). In particular, it has been shown that without Siridophores, spots of melanocyte aggregates form instead of stripes, which is contrary to what these Turing reactiondiffusion models predict. These findings have paved the way for more detailed modelling, such as that of Volkening and Sandstede, 2018, who demonstrated (using an offlattice individualbased model) the need for understanding Siridophore behaviour when representing all three celltypes. For these reasons, we consider an inclusive modelling approach, incorporating the crucial celltype Siridophores and the full range of interactions depicted above.
Here, in a bottomup 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 agentbased model incorporating all three pigment celltypes and their documented cellular interactions. We use observations of a set of three mutants that each lack an individual celltype, plus the three double mutant combinations lacking pairs of celltypes, to deduce the key rules likely underpinning Siridophore 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 celltypes, and show that these too are successfully matched in silico by our model.
Our work demonstrates that current biological understanding, alongside simple assumptions about Siridophore 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 Siridophores are crucial for stripe formation, suggests a minimum set of their rules, and reveals unexpected subtleties to the phenotypic impact of the wellstudied leo mutant.
Materials and methods
Modelling overview
Request a detailed protocolWe 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 Siridophores 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 Liridophores, 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 celltype 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 redifferentiate 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 Siridophore layer (Hirata et al., 2003; Hirata et al., 2005). For consistency, we model each of the three layers as independent, twodimensional 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 Siridophore domains respectively (Figure 2A). To account for the different packing densities of the cell types, lattice sites within the xanthophore and Siridophore 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 onethird 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 Siridophores along the horizontal myoseptum, scattered melanocytes and dedifferentiated xanthophores (xanthoblasts) scattered across the domain. We model this by populating the central three rows of the Siridophore layer with dense Siridophores, 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 Siridophores 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 selforganisation 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 shortrange 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 longrange 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 longrange 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 Siridophores, loose Siridophores, xanthophores and xanthoblasts, as well as some predictions about Siridophore behaviour which have not been experimentally investigated in the literature. Apart from those involving Siridophores, all the interactions and wherever possible their quantitative properties (strength, frequency etc) come directly from the literature, and are summarised in Figure 4G, 114, 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 Siridophores have not been wellcharacterised experimentally, so we have developed our own predictions based on the literature describing Siridophore behaviour during pattern metamorphosis. It has been shown using clonal cell analysis that during pattern metamorphosis Siridophores 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 Siridophores show a directional bias towards xanthophores in the short range. We propose that dense Siridophores 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 Siridophores show a directional bias away from other loose Siridophores in the short range. We propose that loose Siridophores are repelled by other loose Siridophores as this would facilitate the prompt spreading of loose Siridophores across the stripe regions.
Interestingly, as the Siridophores spread they switch between a loose and dense form, predetermining the positioning of stripes and interstripes consecutively. In the loose form Siridophores are spread and stellate in appearance. In contrast, in the dense form, Siridophores are compact. The transition between the two types appears to be interchangeable. When dense Siridophores initially spread beyond the boundary of the first X0 interstripe, they can later change to loose type (Fadeev et al., 2015). Similarly when loose Siridophores reach an interstripe region, they can aggregate and form dense Siridophores. 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 Siridophores from dense to loose state, suggesting that Tjp1a contributes to regulation of the molecular switch that regulates Siridophore 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 redifferentiating 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 Siridophore 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 celltype. We use these to define the contribution and nature of Siridophore interactions in pattern formation, by considering the outcomes in fish lacking each of the three cell types. Specifically we consider:
Mutants lacking Siridophores: The gene shady (shd) encodes zebrafish leukocyte tyrosine kinase (Ltk) which plays a role in Siridophore specification (Lopes et al., 2008). As a result, strong shd mutants lack all Siridophore types. The resultant adult pattern consists of a widened X0 region of xanthophores, which are flanked dorsally and ventrally by melanocytes organised as spotlike 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 Siridophores and xanthophores with irregular borders, accompanied by spots of dense Siridophores and xanthophores ventrally. The rest of the flank is filled with loose form Siridophores (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 Siridophores. In the regions lacking melanocyte aggregation (the ‘saltandpepper’ region), Siridophores 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 prepattern, but instead is generated by cellcell interaction.
Upon evaluating these mutants, we make the following deductions about Siridophore shape transitions during pattern formation:
Siridophores are initially dense and cannot change shape autonomously. This is based on observations of mutants nac;pfe which only contain Siridophores and in which the adult phenotype consists of dense Siridophores in a coherent sheet across the domain (Frohnhöfer et al., 2013). In contrast, pfe and nac both exhibit loose and dense Siridophores (Frohnhöfer et al., 2013), suggesting that both melanocytes and xanthophores are capable of facilitating Siridophore 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 Siridophores from loose to dense. We propose these interactions for the following reasons. Firstly, in pfe and WT, dense Siridophores are associated with lack of melanocytes, for example within the interstripes, whilst loose Siridophores are associated with melanocytes, for example in the stripe region. Since we predict that melanocytes are required for dense Siridophores transition the simplest assumption is that melanocytes promotes dense to loose transitions in the short range. Since loose Siridophores can reaggregate 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 Siridophores are associated with xanthophores, whilst loose Siridophores are associated with a lack of xanthophores (Frohnhöfer et al., 2013). Since Siridophores initially appear in dense form and become loose for example in nac, when there are xanthophores in a low density local to Siridophores 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, Siridophores 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, 1518. 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 Siridophores 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 Siridophore assumptions").
Comparing simulated fish with real data
Request a detailed protocolIn 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 postfertilisation (dpf) and simulated HAA versus SL, respectively, compared with the averaged data (Parichy et al., 2009). These figures demonstrate that whilst growth rates are variable within simulations (as seen in real fish), the mean of our simulated rates approximately matches that in real fish.
Results
Modelling simulations
Having deduced this minimal ruleset 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 Siridophores 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 Siridophores. 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 ruleset.
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 Siridophores from dense to loose form at the edges of interstripe X0. At stage SP, we observe the spreading of loose Siridophores 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 cellcell 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 Siridophores from the initial conditions and switch off Siridophore 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 pseudostripes as the spots of melanophore aggregates that appear in a stripelike orientation reminiscent of that in WT fish. We describe the pseudostripes in the order they appear as in WT fish. For example, we define pseudo1D and pseudo1V to be the first pseudostripes.
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 sametype 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 pseudostripes of melanocyte spots surrounded by xanthophores. We found that in 100 simulations, 100% of shd stage J+ mutants observed two pseudostripes (1D and 1V) just as in Figure 5H. See Video 2 for a movie of the simulation.
Moreover, pseudostripes varied in how stripelike 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 ‘pseudostripe’ 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 pseudostripe 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 Siridophores 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 Siridophores in the initial X0 interstripe. At stage SP (Figure 5K,K’) we observe the switch of dense Siridophores to loose form and the subsequent spreading of loose Siridophores. 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.
The pfeffer mutant
At stage PB (Figure 5M,M’), we populate the domain with a central stripe of dense Siridophores 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 Siridophores) into spots, surrounded by a sea of dense Siridophores peppered with black melanocytes. In one hundred simulations, the median number of pseudostripes 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 pseudostripes 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 Siridophores 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 Siridophore interactions, our secondary set comprises mutants with Siridophorerelated phenotypes: rose (rse) homozygotes, which show a reduction of Siridophore numbers; schachbrett (sbr) homozygotes, which show a delay in Siridophore 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 Siridophores (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 cellautonomously in Siridophores; homozygous mutants result in a reduction of Siridophores 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 Siridophores 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 Siridophores spread less far from the horizontal myoseptum. At stage SP, (Figure 7J,J’) the stripe boundaries at X0 are poorly defined, and dense Siridophores are still largely associated with the X0 region, with only a few loose Siridophores 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.
The sbr mutant
The sbr gene encodes Tight Junction Protein 1a (Tjp1a), which is expressed cell autonomously in dense Siridophores (but not in loose Siridophores) and truncated in sbr mutants; in adult sbr mutants Siridophore 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 Siridophores transitioning to loose Siridophores at the edge of the interstripe, in sbr, dense Siridophores remain dense and spread over melanocytes dorsally and ventrally bidirectionally . At later stages, some dense Siridophores do switch to loose Siridophores. See Video 6 for a movie of the simulation.
We interpret the sbr mutation as causing a delay in signaling driving the transition of Siridophores from dense to loose Siridophore. 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 Siridophore 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 Siridophores (white cells) and not just with loose Siridophores (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 pseudostripe pattern.
The cho mutant
Homozygous cho mutant larvae lack the horizontal myoseptum (Svetic et al., 2007). As a result, dense Siridophores are prevented from traveling via the horizontal myoseptum to generate the initial stripe of dense Siridophores seen in WT at stage PB. Instead loose Siridophores 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 Siridophores at the PB stage (Figure 7O,O’) and instead place 200 loose Siridophores at random across the Siridophore 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 immunoglobulinlike 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 Siridophores and xanthophores with dense Siridophores, 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 Siridophore 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 Siridophores along the horizontal myoseptum using a Siridophorespecific 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 Siridophores 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 Siridophore 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.
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 hemisegment 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\le 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 onlattice domains and is called the Square Uniform PCF (Gavagnin et al., 2018) adapted for multiple celltypes (see Appendix 3 and Dini et al., 2018). We describe the PCF as homotypic when we are measuring the spatial correlation of one celltype and heterotypic when we consider the spatial correlation between two different celltypes. 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 Siridophores 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 Siridophores occupy the interstripe regions in WT and not the stripe regions at J+, dense Siridophores are spatially correlated at short distances, indicated by a positive value of the PCF at short distances. Conversely, they are anticorrelated 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 Siridophores. 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 Siridophores for WT and cho in Figure 10C. For both cho and WT simulations, the averaged homotypic PCF for dense Siridophores 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 pseudostripes 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 Siridophores for WT and pfe in Figure 10E. For both the WT and pfe simulations, the averaged heterotypic PCF of melanocytes and dense Siridophores 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 Siridophore assumptions
For the less wellstudied Siridophore 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 densetoloose, loosetodense Siridophore transitions first outlined in Figure 4G, 1518 for stripe formation by showing representative images (the model was run 50 times each) of simulations at stage J+ for fish lacking one or more Siridophore transition mechanisms display patterns which diverge from those seen in real fish.
First, we analyse simulations lacking one of the transition types loosetodense or densetoloose (Figure 11 B1–C3). Without a loosetodense transition (Figure 11 B1–B3), note how in all cases (WT, pfe, nac) only one pseudointerstripe is preserved: the initial X0 interstripe. The X0 interstripe is surrounded by a sea of loose Siridophores 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 loosetodense interactions are important for generating subsequent interstripes. Interestingly, without loosetodense 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 densetoloose transition (Figure 11 C1–C3), the Siridophores form a dense sheet over the entire domain with xanthophores and melanophores scattered across the domain. This demonstrates the necessity for Siridophores to be able to transition between dense and loose form.
Next, we consider removing each of the criteria required for an Siridophore 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 longrange 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, Siridophores are consistently promoted to dense type, with a nonzero 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, shortrange xanthophore inhibition is removed from criteria 16. As this is exclusively a xanthophoreiridophore interaction this only effects WT and nac. Without the promotion of xanthophores in the short range, Siridophores change from dense to loose in the interstripes, making the interstripes appear faded. In Figure 11 G1–G3, shortrange melanophore inhibition is removed from criteria 18. Therefore Siridophore 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 Siridophores 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 indepth study, we have incorporated all of the interactions that we have identified from the literature. Consequently, there may be some inbuilt 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 chaserun 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 Siridophores and the number of xanthophores currently on the domain. This is because we assume that dense Siridophores and xanthophores positively influence the rate of melanocyte birth in the long range. In Figure 12E, we change the differentiation rate so it is only dependent on the number of xanthophores currently on the domain and not the number of Siridophores, effectively reducing the rate of differentiation of melanocytes in wildtype fish. The resultant pattern is still striped; however, it is less organised. In Figure 12F–G, we remove the mechanism allowing longrange 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 Siridophores and xanthophores in the longrange and/or loose Siridophores 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 Siridophores (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 Siridophores, Figure 12N). In Figure 12L,M, we change melanocyte and xanthophore movement so that they are no longer influenced by each other (no chaserun 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, wildtype 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 finegrained, celllevel model is the ability to perform in silico experiments that can be directly related to reallife 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 Siridophore interstripe orientation alone does not determine the orientation of stripes and interstripes
Previously it has been hypothesised that the horizontal orientation of the initial Siridophore 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 Siridophore 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 Siridophores (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 (dorsoventrally) down the centre of the domain of width three instead of as a band of width three, dense Siridophores 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 Siridophore 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 Siridophore interstripe is important for successful pattern formation
Another way to better understand the role of the initial Siridophore interstripe is to alter its position in the dorsoventral axis so that it appears more ventrally on the initial domain. We might naively predict, given that it has been hypothesised that dense Siridophores 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 Siridophores 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 adultwidth 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 preestablished 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.
Siridophores 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 Siridophores (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 Siridophore 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 heighttolength 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üssleinVolhard, 2003; Irion et al., 2014). As a result, leo mutants display a leopardlike 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 (Irion et al., 2014). For example, the flank of double mutant leo;nac is covered by xanthophores and dense Siridophores (Figure 15B’’). This is in contrast to nac which also contains large patches of loose Siridophores (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 Siridophores (Watanabe et al., 2006; Maderspacher and NüssleinVolhard, 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 Siridophore 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 invitro 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 Siridophores) 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 melanocytexanthophore interactions. From these observations, we conclude that hypothesis one alone is insufficient for leo pattern formation.
To deduce the other cellcell 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 cellcell interaction changes that could generate the double mutant patterns.
The leo;nac mutant displays an absence of loose Siridophores compared to nac, suggesting that signalling of xanthophores which promote Siridophore 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 longrange 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 Siridophores 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 allencompassing, as so far, none of our hypotheses effect melanocyteiridophore 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 cellcell 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 dorsolaterally, 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 Siridophores and xanthophores. Also, we can replicate leo;nac and leo;shd.
Finally, we attempt to elucidate the melanocyteiridophore cellcell 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 Siridophores, suggesting that leo affects directed differentiation of melanocytes by dense Siridophores. This leads us to two extra hypotheses for the effects of the leo gene:
Hypothesis 4: Melanocytes lose death signals from local dense Siridophores and, as a result, can differentiate in dense Siridophore zones.
Hypothesis 5: Melanocytes lose directed signalling from Siridophores 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 preexisting 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 Siridophores in a sea of dense Siridophores 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 Siridophores as expected (Figure 15I). Our simulations of leo;pfe (Figure 15J’) is fully populated with Siridophores 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 (Tjp1a), which is expressed cell autonomously in dense Siridophores 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 Siridophores 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 Siridophores that are surrounded by xanthophores and dense Siridophores. Comparison with our simulation of leo (Figure 15H’) demonstrates that our simulation of leo;sbr also displays more elongated spots than those of our leo simulations.
Robustness of the leo assumptions in generating spots
As a test of robustness, we perform a rigorous robustness analysis by carrying out one hundred repeats of the mutant simulations with perturbed parameter values chosen uniformly at random from the range 0.75–1.25 of their described value as in Section "Simulation of WT pattern". Ten of these randomly sampled repeats are given in Figure Appendix 4—figure 1 for leo. We observe that for all one hundred repeats that small perturbations to the rates still generate consistent spots, demonstrating the robustness of the model.
These results of this section demonstrate the remarkable ability of our model to generate the leo single and double mutant phenotypes under a set of specific proposed changes to the model rules; these proposals can be used to guide experimental exploration of the effects of the leo gene.
Discussion
As a result primarily of beautiful experimental work, zebrafish pigment pattern formation has become the best characterised pigment patterning mechanism (Frohnhöfer et al., 2013; Kelsh et al., 2009; Singh and NüssleinVolhard, 2015; Watanabe and Kondo, 2015; Patterson and Parichy, 2019). Zebrafish pigment pattern formation requires three pigment producing cell types: melanocytes, xanthophores and Siridophores to generate the WT stripes (Frohnhöfer et al., 2013). In an attempt to decipher the mechanisms underlying pattern formation, most previous mathematical models have largely focussed on xanthophores and melanocytes and neglected Siridophores (Nakamasu et al., 2009; Bullara and De Decker, 2015; Volkening and Sandstede, 2015; Painter et al., 2015; Bloomfield et al., 2011). However, recent studies have shown that Siridophores are a crucial component, with Siridophore transitions between dense and loose driving pattern formation. First, the early provision of dense Siridophores through the horizontal myoseptum orients stripes (Frohnhöfer et al., 2013; Kelsh et al., 2009; Singh and NüssleinVolhard, 2015; Watanabe and Kondo, 2015; Patterson and Parichy, 2019). Moreover, Siridophores sequentially predetermine stripes and interstripes via respective shape changes in these regions (Frohnhöfer et al., 2013; Fadeev et al., 2015; Krauss et al., 2014).
Here, we have taken a bottom up approach to modelling zebrafish pattern formation with the aim of testing whether the experimentally defined set of biological rules for zebrafish pigment pattern formation might be sufficient to explain both the WT and the diversity of mutant pigment patterns. We used an individual based modelling approach incorporating all five celltypes deemed important for pattern formation in zebrafish. We formalised all respective cellcell interactions mathematically, with interaction strengths, parametrised, where possible, by the biological literature (see Supplementary file 5). For the less wellstudied Siridophore 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 Siridophore 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 celltypes, 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 Siridophore 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 Siridophores in a dense configuration than those of the loose form (Fadeev et al., 2015). Furthermore, double mutant and chimaeric studies show that sbr acts cellautonomously within the Siridophores to control adult pigment pattern formation (Fadeev et al., 2015). These authors also show that in sbr mutants the transition from dense to loose Siridophores is delayed, suggesting that this transition somehow depends upon Tjp1a in dense Siridophores (Fadeev et al., 2015). Here, we test the patterning impact of this interpretation, by incorporating delayed Siridophore 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 Siridophores 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 inbetween 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 Siridophores play a more important role in melanocyte differentiation than xanthophores. We found that by implementing a stronger signaling capacity of Siridophores 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 Siridophores in the long range appeared to align subsequent stripes in pfe and WT, resulting in consistency of pseudostripe and stripe width respectively. The comparatively reduced signalling success rate for melanocyte differentiation from xanthophores in the long range reduced the number of pseudostripes 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 Siridophores play an important role in determining stripe and interstripe width (without Siridophores, 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 Siridophores (Irion et al., 2014). Connexins play an important role in cell–cell signalling and communication so it has been postulated that the leo gene is involved in the signalling between melanocytes and xanthophores as well as signalling cues to Siridophores regarding shapetransitions (Irion et al., 2014). Previous to our investigation, there had been one study of the individual cellcell 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 cellcell 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 Siridophore 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 reactiondiffusiontype models posit that combinations of short and long range dynamics between melanocytes and xanthophores generate stripe patterns. Indeed, a lot of the excitement around such models is the ease with which small parameter value changes can sometimes result in diverse patterns, many readily recognisable from nature (Watanabe and Kondo, 2015; Metz et al., 2011; Maini et al., 2012). A major difference between our model and Turing reactiondiffusion models is that small parameter changes in our model do not typically generate qualitatively different patterns, whereas Turing reactiondiffusion models can show substantial pattern changes in response to small alterations to parameters near to bifurcation points (Watanabe and Kondo, 2015; Maini et al., 2012). We suggest that the added complexity of the real system has evolved to make the patterning process robust, with partially redundant mechanisms insulating against the impact of stochastic variation during pattern formation.
Our modelling approach is analogous to that adopted in another recent study of zebrafish pattern formation, one which independently attempts to understand Siridophore contributions to the process (Volkening and Sandstede, 2018). In a similar fashion, Volkening et al generated an offlattice model incorporating Siridophores 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 onlattice model, whilst Volkening et al use an offlattice model. Using an onlattice model allowed us to incorporate volume exclusion by other cells directly. Secondly, we used a continuoustime 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 Siridophores 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 Siridophore 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 Liridophores, need to be analysed and incorporated into a model that extends to these later developmental stages.
From an analytical perspective, an advantage of our onlattice model is its amenability to the derivation of a continuum model, although we note that continuum approximations to offlattice individualbased 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 heighttolength 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 dorsoventral 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 dorsoventral 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 onlattice model, implementing the current biological understanding of adult zebrafish pigment pattern formation, strongly supports the validity of these experimental interpretations, motivating the detailed investigation of their molecular bases. Our model also highlights areas where knowledge is currently incomplete and, importantly, has allowed in silico investigations to identify plausible mechanisms that require experimental testing.
Appendix 1
Implementing the model
This text provides a full description of how a simulation of a WT fish is implemented. All notation used in appendix 1 is summarised in Supplementary files 1 and 2.
Model overview
We model five different cell types; yellow xanthophores ($X$), unpigmented xanthoblasts (${X}^{b}$), black melanocytes ($M$), silver dense Siridophores (${I}^{d}$) and blue loose (${I}^{l}$) Siridophores. 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{\bm{X}}$, a melanocyte layer represented by matrix $\mathit{\bm{M}}$, and an Siridophore layer represented by a matrix $\mathit{\bm{I}}$. We denote the length and height of domain $D\in \{\mathit{\bm{X}},\mathit{\bm{M}},\mathit{\bm{I}}\}$ 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{\bm{M}}$ denoted $\mathit{\bm{M}}(i,j,t)$ can either be occupied by $M$, i.e. $\mathit{\bm{M}}(i,j,t)=M$ or be unoccupied, $\mathit{\bm{M}}(i,j,t)=0$. Similarly $\mathit{\bm{I}}(i,j,t)$ can either be occupied by ${I}^{l}$, i.e. $\mathit{\bm{I}}(i,j,t)={I}^{l}$, be occupied by ${I}^{d}$, i.e. $\mathit{\bm{I}}(i,j,t)={I}^{d}$, or be unoccupied, $\mathit{\bm{I}}(i,j,t)=0$. Finally, $\mathit{\bm{X}}(i,j,t)$ can either be occupied by $X$, i.e. $\mathit{\bm{X}}(i,j,t)=X$, ${X}^{b}$, i.e. $\mathit{\bm{X}}(i,j,t)={X}^{b}$ or be unoccupied, $\mathit{\bm{X}}(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{\bm{X}}$ and $\mathit{\bm{I}}$ are of size 0.02 mm × 0.02 mm, whilst lattice sites in $\mathit{\bm{M}}$ are of size 0.04 mm × 0.04 mm (Parichy et al., 2000a). We denote these as ${\mathrm{\Delta}}_{\mathit{\bm{X}}}={\mathrm{\Delta}}_{\mathit{\bm{I}}}=0.02$ mm and ${\mathrm{\Delta}}_{\mathit{\bm{M}}}=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{\bm{X}}}^{H}={\mathrm{\Pi}}_{\mathit{\bm{I}}}^{H}=6$, ${\mathrm{\Pi}}_{\mathit{\bm{X}}}^{L}={\mathrm{\Pi}}_{\mathit{\bm{I}}}^{L}=12$, ${\mathrm{\Pi}}_{\mathit{\bm{M}}}^{H}=3$, ${\mathrm{\Pi}}_{\mathit{\bm{M}}}^{L}=6$.
Relating size to stage
Zebrafish development is described in the literature with respect to the standard length (SL) of the fish (Parichy et al., 2009). We associate our domain at time $t$ to different zebrafish developmental stages by calculating a simulated SL (${\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{\bm{X}},\mathit{\bm{M}},\mathit{\bm{I}}\}$ 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.
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 ${\mathrm{\Omega}}_{H}(0)=1mm$ and ${\mathrm{\Omega}}_{L}(0)=2mm$. Therefore, ${\mathrm{\Pi}}_{\mathit{\bm{M}}}^{L}(0)$=25, ${\mathrm{\Pi}}_{\mathit{\bm{M}}}^{L}(0)$=50, ${\mathrm{\Pi}}_{X}^{H}(0)={\mathrm{\Pi}}_{I}^{H}(0)=50$ and ${\mathrm{\Pi}}_{\mathit{\bm{X}}}^{L}(0)={\mathrm{\Pi}}_{\mathit{\bm{I}}}^{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{\bm{I}}$ are fully occupied with ${I}^{d}$ cells at time $t=0$. This represents the Siridophores which differentiate along the horizontal myoseptum between 21 and 25 dpf. We also place ${N}_{M}^{T}(0)=50M$ uniformly at random on $\mathit{\bm{M}}$ 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 dedifferentiate 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{X}^{b}$ uniformly at random on $\mathit{\bm{X}}$ so
This is to represent the dedifferentiated 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.
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 pseudostripe patterns in rse and shd which have reduced or entirely absent Siridophores. We model this by occupying all the sites, regardless of prior occupancy, along the middle three rows of $\mathit{\bm{X}}$ with xanthophores ($X$) at stage SP that is sites $\mathit{\bm{X}}(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 ${\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.
Modeling cell–cell interactions
Continuous events 1–15 given in Supplementary file 4 are mediated by different cellcell 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 shortrange interactions, and an annulus with an outer radius of approximately half an adult stripe width 0.24mm represent longrange dynamics.
We denote ${S}_{D}(r,k)$ and ${L}_{D}(r,k)$ where $D\in \{\mathit{\bm{X}},\mathit{\bm{I}},\mathit{\bm{M}}\}$ 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{\bm{M}},\mathit{\bm{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–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{\bm{M}}}(r,k)$, where $\mathit{\bm{M}}(r,k)$ is the central site marked in grey. ${S}_{\mathit{\bm{M}}}^{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{\bm{M}}(r,k)$ at time $t$ is given by ${N}_{M,M}^{S}(r,k,t)={S}_{\mathit{\bm{M}}}^{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{\bm{M}}(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{\bm{X}}}(r,k)$, where $\mathit{\bm{X}}(r,k)$ is the central site marked in grey. ${S}_{\mathit{\bm{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 $\mathit{\bm{X}}(r,k)$ at time $t$ is given by ${N}_{X,X}^{S}(r,k,t)={S}_{\mathit{\bm{X}}}^{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{\bm{X}}(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{\bm{M}}}(r,k)$, where $\mathit{\bm{M}}(r,k)$ is the central site marked in grey. ${L}_{\mathit{\bm{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 $\mathit{\bm{M}}(r,k)$ is given by ${N}_{M,M}^{L}(r,k,t)={L}_{\mathit{\bm{M}}}^{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{\bm{M}}(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{\bm{X}}}(r,k)$, where $\mathit{\bm{X}}(r,k)$ is the central site marked in grey. ${L}_{\mathit{\bm{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 $\mathit{\bm{X}}(r,k)$ is given by ${N}_{X,X}^{L}={L}_{\mathit{\bm{X}}}^{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{\bm{X}}(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{\bm{X}}(i,j)$, corresponds to a quarter of site $\mathit{\bm{M}}(\lceil \frac{i}{2}\rceil ,\lceil \frac{j}{2}\rceil )$, this explains case two in Equation 7. Similarly each site $\mathit{\bm{M}}(i,j)$ corresponds to four sites on $\mathit{\bm{X}}$, specifically $\mathit{\bm{X}}(2i,2j)$, $\mathit{\bm{X}}(2i,2j1)$, $\mathit{\bm{X}}(2i1,2j)$ and $\mathit{\bm{X}}(2i1,2j1)$. 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
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 onlattice attractionrepulsion 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\in \mathit{\bm{O}}=\{No,So,Ea,We,NE,SE,SW,NW\}$. (See Figure 3(D) in the main text for an illustration of these directions).
The directional neighbours of each cell is given in Figure 2E,F. Figure 2E demonstrates the directional neighbourhoods of a focal cell $C\in \{X,{X}^{b},{I}^{d},{I}^{l}\}$ located in position $D(i,j)$ (marked in grey) where $D\in \{\mathit{\bm{X}},\mathit{\bm{I}}\}$. For the cell moving one space on $D$, the cell considers the occupancy of sites in $S_{k}{}_{{\mathrm{\Delta}}_{X}}{}^{D}(i,j)$ marked in red for $D=\mathit{\bm{M}}$ (top) and $D\in \{\mathit{\bm{X}},\mathit{\bm{I}}\}$ (bottom). (Q–X) Figure 2F demonstrates the directional neighbourhoods of a focal cell $M$ located in position $\mathit{\bm{M}}(i,j)$ (marked in grey). For the cell $M$ moving one space on $\mathit{\bm{M}}$, the cell considers the occupancy of sites in $S_{k}{}_{{\mathrm{\Delta}}_{M}}{}^{D}(i,j)$ marked in red for $D\in \{\mathit{\bm{X}},\mathit{\bm{I}}\}$ (top) and $D=\mathit{\bm{M}}$ (bottom).
We define the probability that a cell attempts to move in a given direction with orientation $O$ as ${P}_{O}$ (where $\sum _{O\in \mathit{\bm{O}}}{P}_{O}=1$). We calculate ${P}_{O}$, using bias vector,
where $\mathit{\bm{A}}$ ($\mathit{\bm{R}}$) is a matrix whose entries are the number of neighbouring cells within different segments of the attraction (repulsion) range and $\underset{\xaf}{a}$ ($\underset{\xaf}{r}$) is a weight of attraction (repulsion) vector.
Matrices $\mathit{\bm{A}}$ and $\mathit{\bm{R}}$ are defined as follows.
where for example, ${O}_{C,{C}_{1}}(r,k,t)$ denotes the number of cells of type ${C}_{1}\in \{X,M,{I}^{d},{I}^{l},{X}^{b}\}$ within range 0.04 mm, that are orientation $O$ from the focal cell $C$ located in position $D(r,k)$ at time $t$. This number is calculated differently depending on whether the focal site lies on lattice domain $D\in \{\mathit{\bm{M}},\mathit{\bm{X}},\mathit{\bm{I}}\}$. To compute ${O}_{C,{C}_{1}}(r,k,t)$ we calculate the following. For any $O$ we define subsets ${S}_{D}^{O}(r,k,t)$ of ${S}_{D}(r,k,t)$ (where cell type $C$ lies on domain type $D$) for each direction $O$;
These are the sets of sites in the short range distance, north, south, east and west of focal site $D(r,k)$. Examples for $D=\mathit{\bm{M}}$ and $\mathit{\bm{X}}$ are sites marked in red on the left of Figure 2E–F respectively. By extension we define
These are the sets of sites in the short range distance, northeast, northwest, southeast and southwest of focal site $D(r,k)$. We denote ${S}_{D}^{C,O}(r,k,t)\subseteq {S}_{D}^{O}(r,k,t)$ to be all the sites in shortrange distance in orientation $O\in \{No,So,Ea,We,NE,SE,SW,NW\}$ occupied by cell type $C$. Hence we define
Therefore to compute the number of cells in each of the eight directions, we compute where $C$ lies on $D$ and ${C}_{1}$ lies on ${D}_{1}$ by;
Next we define the weight of the attractionrepulsion vector by:
where ${a}_{C{C}_{1}},{a}_{C{C}_{2}},{a}_{C{C}_{3}}$ (${\rho}_{C{C}_{1}},{\rho}_{C{C}_{2}},{\rho}_{C{C}_{3}}$) are the weights of attraction (repulsion) of cells of type $C$ to cells of the type ${C}_{1}$, ${C}_{2}$ and ${C}_{3}$, respectively defined in Supplementary file 5. Vectors $\underset{\xaf}{a}$ and $\underset{\xaf}{b}$ are the vectors used to calculate $\underset{\xaf}{v}$ in Equation 9. To account for the difference in distance between adjacent neighbours and diagonal neighbours we normalise accordingly so that movement diagonally is proportionately less frequent than moving to adjacent squares.
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)$. Next values ${P}_{O}$ are calculated using Equation 27–33 and Equation 35. A random number $u\sim U(0,1)$ is generated. The cell attempts movement to $D(r,k1)$ (West) if $u\in [0,{P}_{W})$, $D(r,k+1)$ (East) if $u\in [{P}_{W},{P}_{W}+{P}_{E})$etc. up to $[\sum _{O\in \mathit{\bm{O}}}{P}_{O}{P}_{SW},\sum _{O\in \mathit{\bm{O}}}{P}_{O}]$. The cell movement is successful and the cell moves from site $D(r,k)$ to the nearest neighbouring site in direction $O$, provided this site is empty. Otherwise the movement is aborted. Movement bias is specified by weights of attraction and repulsion ($\underset{\xaf}{a}$ and $\underset{\xaf}{r}$). Relevant parameters for ${a}_{{C}_{i}{C}_{j}}$ and ${\rho}_{{C}_{i}{C}_{j}}$ for different ${C}_{i}$ and ${C}_{j}$ are given in Supplementary file 5 . If no such interaction bias is specified then there is no known attraction or repulsion dynamics between those cell types in the literature and so we assume that movement is not biased by this cell type. Note that for cell ${X}^{b}$ there are no known short range interactions between ${X}^{b}$ and other cell types so we model ${X}^{b}$ movement as random, i.e, ${P}_{O}=\frac{1}{8}$, for $O\in \mathit{\bm{O}}$ at all times.
Siridophore, xanthophore and xanthoblast proliferation (continuous time events 6–7)
Given the event is chosen such that a cell of type $C\in \{X,{X}^{b},{I}^{d},{I}^{l}\}$ is determined to proliferate we choose a random cell of type $C$ whose position is given by $D(r,k)$, $D\in \{\mathit{\bm{X}},\mathit{\bm{I}}\}$. Next a random number $u\sim 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\in [0,1/4)$, $D(r+1,k)$ if $u\in [1/4,1/2)$, $D(r1,k)$ if $u\in [1/2,3/4)$ or $D(r,k1)$ if $u\in [3/4,1]$.
For cell types $C\in \{{I}^{d},{I}^{l},{X}^{b}\}$, 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
Equation 36 is based on the following assumptions: (1) Dense Siridophores promote xanthoblasts differentiation into xanthophores in the short range (Patterson and Parichy, 2013). Dense Siridophores express xanthogenic Colony Stimulating Factor1 (Csf1) (Patterson and Parichy, 2013) which is essential for xanthophore differentiation, proliferation and survival, allowing unpigmented xanthoblasts near to the dense Siridophores 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 $\mathit{\bm{M}}(r,k)$ is chosen at random. A differentiation event is successful and a new $M$ appears in this site if the site $\mathit{\bm{M}}(r,k)$ is empty and the following is true:
where $\alpha =2.5$, $\beta =\gamma =\kappa =3,w=4$. Equation 37 is based on the findings that dense Siridophores 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 Siridophores; however, melanocytes frequently settle adjacent to dense Siridophores, 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 Siridophores 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\in U(0,1)$ and we find that $u<0.01$ in combination with the condition
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 $\mathit{\bm{X}}$. Suppose the chosen ${X}^{b}$ lies in site $\mathit{\bm{X}}(r,k)$. The differentiation event is successful and a $X$ replaces the ${X}^{b}$ in this site if the site $\mathit{\bm{M}}(\lceil \frac{r}{2}\rceil ,\lceil \frac{k}{2}\rceil )$ 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.
Equation 41 is based on the following assumptions; (1) Dense Siridophores promote xanthoblasts differentiation into xanthophores in the short range (Patterson and Parichy, 2013). Dense Siridophores express Csf1, allowing unpigmented xanthoblasts near to the dense Siridophores to mature into xanthophores (Frohnhöfer et al., 2013; Walderich et al., 2016). (2) melanocytes inhibit xanthophore specification in the short range (Nakamasu et al., 2009).
Alternatively the differentiation event is successful if:
holds and a randomly generated number $u\in U(0,1)$ is such that $u<0.01$, or Equation(42) holds and either
is true i.e, the total cell density with all cells combined is less than either 0.2 or 0.4 for different time milestones. We enforce a lower probability ($u$ < 0.01) for the alternative event (Equation 42) since this differentiation event is not influenced by cues from other cells and thus is less likely. Equation 43 and Equation 44 are based on assumptions that if the domain is empty after some key time points, this promotes delayed differentiation of xanthoblasts as is observed in nac (Frohnhöfer et al., 2013). Since melanocytes compete in the short range with xanthophores (Nakamasu et al., 2009) we only allow this to occur when constraints described by Equation(42) are simultaneously held, for consistency.
Iridophore transitions (continuous time event 10–11)
If an ${I}^{d}$ to ${I}^{l}$ shape transition is chosen, then an ${I}^{d}$ is chosen uniformly at random from domain $\mathit{\bm{I}}$ and evaluated for meeting transition criteria. The transition is successful and the chosen ${I}^{d}$ in position $\mathit{\bm{I}}(r,k)$ is replaced with ${I}^{l}$ in this site if either both
or alternatively
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 $\mathit{\bm{I}}$ and evaluated for meeting transition criteria. This transition is successful and the chosen ${I}^{l}$ in position $\mathit{\bm{I}}(r,k)$ is replaced with ${I}^{d}$ in this site if either
or alternatively
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.
Melanocyte death (continuous time event 12)
If a melanocyte death event is chosen, then a cell of type $M$ is chosen at random from $\mathit{\bm{M}}$. Suppose the chosen $M$ lies in site $\mathit{\bm{M}}(r,k)$. The death event is successful and the melanocyte is removed from site $\mathit{\bm{M}}(r,k)$ i.e. $\mathbf{\mathbf{M}}\mathit{}\mathrm{(}r\mathrm{,}k\mathrm{)}$ is set to 0, if
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\in U(0,1)$, is such that $u<0.001$ and both of the following two cases hold:
where $\sigma =\omega =3$. Equation 50 is based on findings that $M$ appear to inhibit the survival of $M$ in the long range and that $X$ promote the survival of $M$ in the long range (Takahashi and Kondo, 2008). In Equation 51, ${I}^{l}$ are a proxy for Liridophores, 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 $\mathit{\bm{X}}$ that meets the following criterion: suppose the chosen ${X}^{b}$ resides in site $\mathit{\bm{X}}(r,k)$, then the corresponding site on $\mathit{\bm{M}}$, $\mathit{\bm{M}}(\lceil \frac{r}{2}\rceil ,\lceil \frac{k}{2}\rceil )$ must be empty. We simulate this by randomly choosing ${X}^{b}$ on $\mathit{\bm{X}}$, and checking whether $\mathit{\bm{M}}(\lceil \frac{r}{2}\rceil ,\lceil \frac{k}{2}\rceil )$ is empty. If so, we continue to the next stage. Otherwise we repeat through all ${X}^{b}$ on $\mathit{\bm{X}}$ until we either find a suitable ${X}^{b}$. If no ${X}^{b}$ satisfy this criteria at time $t$ then we abort the cellcell 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 $\mathit{\bm{M}}(\lceil \frac{r}{2}\rceil ,\lceil \frac{k}{2}\rceil )$. 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;
Next we translate position $(i,j)$ on $\mathit{\bm{X}}$ to its corresponding position on $\mathit{\bm{M}}$. If $\mathit{\bm{M}}(\lceil \frac{i}{2}\rceil ,\lceil \frac{j}{2}\rceil ,t)=M$, then the melanocyte is pulled from its site $\mathit{\bm{M}}(\lceil \frac{i}{2}\rceil ,\lceil \frac{j}{2}\rceil )$ into $\mathit{\bm{M}}(\lceil \frac{r}{2}\rceil ,\lceil \frac{k}{2}\rceil )$i.e.$\mathit{\bm{M}}(\lceil \frac{i}{2}\rceil ,\lceil \frac{j}{2}\rceil ,t)=0$ and $\mathit{\bm{M}}(\lceil \frac{r}{2}\rceil ,\lceil \frac{k}{2}\rceil ,t)=M$. The action is deemed complete. Otherwise, if $\mathit{\bm{M}}({i}_{M},{j}_{M},t)=0$, i.e, the site is empty, then the process repeats for the same ${X}^{b}$. From here, $P$ becomes $P\backslash (i,j)$ and a new $(i,j)$ pair are generated uniformly at random from $P$ until either an $M$ is found within the range specified or until $P$ is empty (i.e. after 28 tries). At this point, we will have determined that there were no melanocytes close enough to the xanthoblast for the xanthoblast to successfully pull.
Growth (continuous time events 14–15)
Given a growth event is chosen in the horizontal (vertical) direction, from each column (row), a site location is chosen uniformly at random at which insertion of new empty site above or below (randomly chosen) will occur. See Appendix 1—figure 1CD for an example of a growth event in the vertical direction.
To account for the different lattice sizes, each time a growth event in the horizontal (vertical) direction occurs, ${\mathrm{\Pi}}_{\mathit{\bm{M}}}^{H}(t)$ (${\mathrm{\Pi}}_{\mathit{\bm{M}}}^{L}(t)$) increases by one and ${\mathrm{\Pi}}_{\mathit{\bm{X}}}^{H}(t)$, ${\mathrm{\Pi}}_{\mathit{\bm{I}}}^{H}(t)$ (${\mathrm{\Pi}}_{X}^{L}(t)$, ${\mathrm{\Pi}}_{I}^{L}(t)$) increases by two. This means that when a growth event is chosen, the rules described above occur once in the $\mathit{\bm{M}}$ layer, and twice consecutively in $\mathit{\bm{X}}$ and $\mathit{\bm{I}}$. This is to ensure that the domain size remains consistent between the three layers. That is,
for all time $t$.
Appendix 2
Mutant implementation
In this text, we describe how the WT fish simulation is altered in our model to replicate changes known to be present in a variety of mutant fish.
Simulating the shd mutant
The gene shady (shd) encodes zebrafish leukocyte tyrosine kinase (Ltk) which plays a role in Siridophore specification (Lopes et al., 2008). As a result, strong shd mutants lack Siridophores. To simulate this defect, we remove all Siridophores from the initial domain, i.e, we set ${N}_{{I}^{l}}^{T}(0)={N}_{{I}^{d}}^{T}(0)=0$. Since new Siridophores are only generated by the proliferation of existing Siridophores, this means that ${N}_{{I}^{l}}^{T}(t)={N}_{{I}^{d}}^{T}(t)=0$ for all time $t$. (Note: As a result the propensities of events 4, 5, 7, 10 and 11, that is, the movement of ${I}^{l}$, movement of ${I}^{d}$, proliferation ${I}^{l}$ and ${I}^{d}$, transition of ${I}^{d}$ to ${I}^{l}$ and transitions of ${I}^{l}$ to ${I}^{d}$ are always zero during this simulation.)
Simulating the pfe mutant
Gene pfeffer (pfe) (alternatively known as salz (sal)) encodes for colony stimulating factor one receptor (csf1ra) which plays a role in xanthophore specification and migration. In strong alleles, adult fish exhibit no detectable xanthophores in the body of embryos. To simulate this defect we remove all xanthophores and xanthoblasts from the initial domain, i.e. we set ${N}_{{X}^{b}}(0)={N}_{X}(0)=0$. We also remove the fixed event at stage SP where newly differentiated $X$ cells appear along the horizontal myoseptum. Since new xanthophores are only generated by the proliferation of existing xanthophores, and the differentiation of xanthoblasts this means that ${N}_{{X}^{b}}^{T}(t)={N}_{X}^{T}(t)=0$ for all time $t$. (Note: As a result the propensities of events 2, 3, 6, 9, 15, that is, the movement of ${X}^{b}$, movement of $X$, proliferation $X$ and ${X}^{b}$, differentiation of ${X}^{b}$ to $X$ and ${X}^{b}$ pulling of $M$ are always zero during this simulation.)
Simulating the nac mutant
The gene nacre (nac) encodes transcription factor Mitfa (Lister et al., 1999). nac mutants lack melanocytes throughout embryonic and larval development (Lister et al., 1999). To simulate this defect, we remove all melanocytes from the initial conditions, i.e, we set ${N}_{M}(0)$. We also set the propensity of event 8, the differentiation of new $M$, to be 0 for all time $t$. Since this is the only way new $M$ can be produced, ${N}_{M}(t)=0$ for all time $t$. (Note: As a result the propensities of events 1 and 14, that is, the movement of $M$, and death of $M$ will also always be zero during this simulation.)
Simulating the double mutants: nac;pfe, shd;pfe, shd;nac
To simulate the double mutants we combine the effects of two of the three aforementioned mutation types. For example, to simulate nac;pfe we combine the effects of nac and pfe listed above.
Simulating the rse mutant
Rose (rse), encodes the Endothelin receptor B1a (Krauss et al., 2014) and has been shown to acts cellautonomously in Siridophores; homozygous mutants result in a reduction of Siridophores to approximately 20% of that seen in WT (observed in stage PB and adult fish [Frohnhöfer et al., 2013]). To simulate rse we change the initial conditions so that ${N}_{{I}^{d}}(0)$ is one fifth of the usual number. In our WT simulations, ${N}_{{I}^{d}}(0)$ is ${\mathrm{\Pi}}_{\mathit{\bm{I}}}^{L}(0)\times 3=300$ as ${I}^{d}$ occupy the 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 Siridophores still appear along the horizontal myoseptum in rse. We do this by generating a random number $r\in \mathbb{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 ${\mathrm{\Pi}}_{\mathit{\bm{I}}}^{L}(0)=100$. If $\mathit{\bm{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 $\mathit{\bm{I}}$. Furthermore since new Siridophores are only produced by proliferation of preexisting Siridophores, 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, ${\alpha}_{7}(t)=0.24\frac{{N}_{{I}^{d}}^{T}(t)+{N}_{{I}^{l}}^{T}(t)}{60\times 24}$ for all time $t$.
Simulating the sbr mutant
Adult sbr mutants exhibit delayed Siridophore 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 ${\alpha}_{10}(t)=\frac{{N}_{{I}^{d}}(t)}{60\times 24\times 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 Siridophores cannot travel through their usual pathway to generate the initial strip of dense Siridophores at stage PB seen in WT. Instead loose Siridophores 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 Siridophores 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 Siridophores 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 Siridophores in sites uniformly at random across $\mathit{\bm{I}}$ at stage PR. We do this by generating a random number $r$ between 0 and ${\mathrm{\Pi}}_{\mathit{\bm{I}}}^{L}({t}_{PR})$ and a random number $k$ between 0 and ${\mathrm{\Pi}}_{\mathit{\bm{I}}}^{H}({t}_{PR})$, where ${t}_{PR}$ is the time when the simulation first enters stage PR. If $\mathit{\bm{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 $\mathit{\bm{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 immunoglobulinlike 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 ${\alpha}_{8}(t)=0.05\times \frac{{N}_{{I}^{d}}(t)}{2\times 60\times 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 ${\alpha}_{12}(t)=\times \frac{{N}_{M}(t)}{100\times 60\times 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 Siridophores from dense to loose in the long range.
Hypothesis 4: Melanocytes lose death signals from local dense Siridophores and as a result, can differentiate in dense Siridophore zones.
Hypothesis 5: Melanocytes lose directed signalling from Siridophores 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 Siridophore 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 Siridophores 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 longrange.
Appendix 3
Quantitative measures
In this text, we describe in more detail how we take quantitative measurements of our simulations.
Tortuosity
To determine the tortuosity of the X0 stripe, for a specific time point in our simulations the following steps are taken.
An occupancy matrix $\widehat{\mathit{\bm{X}}}$ of zeros and ones is generated such that a matrix entry $\widehat{\mathit{\bm{X}}}(i,j)=1$ if $\mathit{\bm{X}}(i,j)=X$ and 0 otherwise. Sites $(i,j)$ such that $\widehat{\mathit{\bm{X}}}(i,j)=1$ are yellow and are white otherwise.
Representative matrix $\widehat{\mathit{\bm{X}}}$ is ‘cleaned’ using matlab functions bwmorph($\widehat{\mathit{\bm{X}}}$,‘clean’) and bwmorph($\widehat{\mathit{\bm{X}}}$,‘bridge’) consecutively. ‘Clean’ removes any anomalous xanthophores that are not connected (adjacent to) other xanthophores. ‘Bridge’ adds extra xanthophores where there are holes in the population pattern.
An algorithm is applied to $\widehat{\mathit{\bm{X}}}$ to create the outline of the X0 stripe. The algorithm to create the top line (${l}_{t}$) is given below in Algorithm 1. A similar algorithm is used to generate the bottom line (${l}_{b}$).
A line $L$ that represents the middle of the stripe is given by ${l}_{m}(i)=\frac{{l}_{t}(i)+{l}_{b}(i)}{2}$ for $i=1,\mathrm{\dots}{\mathrm{\Pi}}_{\mathit{\bm{X}}}^{L}$.
We smooth $L$ for analysis by applying matlab function smooth. ‘Smooth’ smooths the data in the column vector $y$ using a moving average filter.
Finally, we calculate the tortuosity of the line by computing the total length of $L$ divided by the distance (algorithm given below in Algorithm 2).
Algorithm 1. Algorithm to generate a representative line for the top of stripe X0 (${l}_{t}$)  

$j=\lceil \frac{{\mathrm{\Pi}}_{\mathit{X}}^{H}}{2}\rceil$ $\mathbf{f}\mathbf{o}\mathbf{r}\text{}i\text{}=\text{}1\text{}:\text{}{\mathrm{\Pi}}_{\mathit{X}}^{L}$ Strike = 0 if $\hat{\mathit{X}}(i,j)\ne 0$ or ($\widehat{\mathit{\bm{X}}}(i,j)=0$ and $\widehat{\mathit{\bm{X}}}(i,j+1)\ne 0$) then while $strike<2$ and $j<{\mathrm{\Pi}}_{\mathit{\bm{X}}}^{H}$ do $j=j+1$ if $\widehat{\mathit{\bm{X}}}(i,j)=0$ then $strike=strike+1\u25b7$ else $strike=0$ end if end while $j=j2\u25b7$ else while$\widehat{\mathit{\bm{X}}}(i,j)=0$ and $j>{\mathrm{\Pi}}_{\mathit{\bm{X}}}^{H}$ do $j=j1$ end while end if ${l}_{t}(i)=j$ end  ▷Initialise X0 interstripe search in the center. ▷Initialise the number of consecutive empty sites in a column to be zero. ▷If we are near to other X, keep moving up to determine the upper bound of the interstripe. ▷Increment the number of consecutive empty sites in a column by one. Remove the two consecutively empty sites from the total count. ▷If we have overshot the interstripe keep moving downwards until we reach the top of the interstripe. 
Algorithm 2. Algorithm to calculate the tortuosity ($tort$) of line $L$ generated from a simulated X0  

total = 0 $\mathbf{f}\mathbf{o}\mathbf{r}\text{}i\text{}=\text{}2\text{}:\text{}{\mathrm{\Pi}}_{\mathit{X}}^{L}$ $total=total+\sqrt{({l}_{m}(i){l}_{m}(i1){)}^{2}+1}$ end $tort=\frac{total}{{\mathrm{\Pi}}_{\mathit{X}}^{L}}$  ▷Compute the length of $L$using the euclidean distance between consecutive points. 
To measure the tortuosity of the stripes of real fish, a photograph was taken of the fish at stage J+. Next, matlab function ‘getpts’ was used to mark the outline of the X0 stripe by a set of points $[\underset{\xaf}{x},\underset{\xaf}{y}]$ where vector $\underset{\xaf}{x}$ is length $k$. The tortuosity of this line was then computed using Algorithm 3.
Algorithm 3. Algorithm to calculate the tortuosity ($tort$) of line $L$ generated from a real fish  

total = 0 $\mathbf{f}\mathbf{o}\mathbf{r}\text{}i\text{}=\text{}2\text{}:\text{}k$ $total=total+\sqrt{(x(i)x(i1){)}^{2}+(y(i)y(i1){)}^{2}}$ end $tort=\frac{total}{\sqrt{(x(1)x(k){)}^{2}+(y(i)y(k){)}^{2}}}$  ▷Compute the length of $L$using the euclidean distance between consecutive points. 
X0 interstripe width
To determine the width of the X0 interstripe, used in Section "Quantitative analysis of simulations" of the main text, first, occupancy matrices of $\widehat{\mathit{\bm{X}}}$ and $\widehat{\mathit{\bm{I}}}$ of zeros and ones is generated such that a matrix entry $\widehat{\mathit{\bm{X}}}(i,j)=1$ if $\mathit{\bm{X}}(i,j)=X$ and 0 otherwise. Similarly $\widehat{\mathit{\bm{I}}}(i,j)=1$ if $\mathit{\bm{I}}(i,j)={I}^{d}$ and 0 otherwise. We then generate ${l}^{t}$ and ${l}^{b}$ using either $\widehat{X}$ or $\widehat{I}$ as described in Section "Necessity of Siridophore assumptions". From these values, we generate $IS$ (X0 interstripe width) by computing;
Note that the X0 interstripe width is computed using whichever is more appropriate of $X$ on $\mathit{\bm{X}}$ and ${I}^{d}$ on $\mathit{\bm{I}}$ given the mutation. For example, pfe does not have xanthophores, so we would use the distribution of ${I}^{d}$ on $\mathit{\bm{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 onlattice 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}^{{C}_{1},{C}_{2}}(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. $\mathbb{E}[{\widehat{c}}^{{C}_{1},{C}_{2}}(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 $\mathit{\bm{M}}$ we must transform $\mathit{\bm{M}}$ into a matrix $\widehat{\mathit{\bm{M}}}$ of size ${\mathrm{\Pi}}_{\mathit{\bm{X}}}^{H}\times {\mathrm{\Pi}}_{\mathit{\bm{X}}}^{L}$ where $\widehat{\mathit{\bm{M}}}(r,k)=M$ if $\mathit{\bm{M}}(\lceil \frac{r}{2}\rceil ,\lceil \frac{k}{2}\rceil )=M$ and 0 otherwise. Hence the set of site positions distance $m$ from each other on any domain $D\in \{\mathit{\bm{X}},\widehat{\mathit{\bm{M}}},\mathit{\bm{I}}\}$ under zero flux boundary conditions can be given by
Therefore
where ${C}_{1}$ lies on ${D}_{1}$, ${C}_{2}$ lies on ${D}_{2}$ and ${D}_{1},{D}_{2}\in \{\mathit{\bm{X}},\widehat{\mathit{\bm{M}}},\mathit{\bm{I}}\}$. Note that ${S}_{m}$ is the number of site pairs distance $m$ from one another. ${S}_{m}$ is computed in Gavagnin et al., 2018 and is given by
To determine $\mathbb{E}[{\widehat{c}}^{{C}_{1},{C}_{2}}(m)]$ there are three cases.
Case 1: ${C}_{\mathrm{1}}\mathrm{=}{C}_{\mathrm{2}}\mathrm{=}C\mathrm{\in}\mathrm{\{}M\mathrm{,}X\mathrm{,}{I}^{d}\mathrm{,}{I}^{l}\mathrm{,}{X}^{b}\mathrm{\}}$ lies on domain $D\mathrm{\in}\mathrm{\{}M\mathrm{,}X\mathrm{,}I\mathrm{\}}$ (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
This is because there are ${N}_{C}^{T}$ ways to choose a cell of type $C$ and ${N}_{C}^{T}1$ ways to choose a cell of type $C$ given one has been chosen already. There are ${\mathrm{\Pi}}_{\mathit{\bm{X}}}^{H}{\mathrm{\Pi}}_{\mathit{\bm{X}}}^{L}$ possible positions for the first cell and then ${\mathrm{\Pi}}_{\mathit{\bm{X}}}^{H}{\mathrm{\Pi}}_{\mathit{\bm{X}}}^{L}1$ possible positions for the second cell due to volume exclusion on domain $D$.
Case 2: ${C}_{\mathrm{1}}\mathrm{\ne}{C}_{\mathrm{2}}$, where ${C}_{\mathrm{1}}\mathrm{,}{C}_{\mathrm{2}}\mathrm{\in}\mathrm{\{}X\mathrm{,}{I}^{d}\mathrm{,}{I}^{l}\mathrm{,}{X}^{b}\mathrm{\}}$ both lie on domain $D\mathrm{\in}\mathrm{\{}X\mathrm{,}I\mathrm{\}}$ (with volume exclusion on $D$). For example spatial correlation of $X$ with ${X}^{b}$. In this case
Notice, ${S}_{m}$ is multiplied by two here because for each pair $(i,j),(r,k)\in {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
This is because there are ${N}_{{C}_{1}}^{T}$ ways to choose a cell of type ${C}_{1}$ and ${N}_{{C}_{2}}^{T}$ ways to choose a cell of type ${C}_{2}$. There are ${\mathrm{\Pi}}_{\mathit{\bm{X}}}^{H}{\mathrm{\Pi}}_{\mathit{\bm{X}}}^{L}$ possible positions for the first cell and then ${\mathrm{\Pi}}_{\mathit{\bm{X}}}^{H}{\mathrm{\Pi}}_{\mathit{\bm{X}}}^{L}1$ possible positions for the second cell due to volume exclusion on domain $D$.
Case 3: ${C}_{\mathrm{1}}$, ${C}_{\mathrm{2}}$ where ${C}_{\mathrm{1}}\mathrm{,}{C}_{\mathrm{2}}\mathrm{\in}\mathrm{\{}M\mathrm{,}X\mathrm{,}{I}^{d}\mathrm{,}{I}^{l}\mathrm{,}{X}^{b}\mathrm{\}}$ lie on domains ${D}_{\mathrm{1}}$, ${D}_{\mathrm{2}}\mathrm{\in}\mathrm{\{}\widehat{M}\mathrm{,}X\mathrm{,}I\mathrm{\}}$ where ${D}_{\mathrm{1}}\mathrm{\ne}{D}_{\mathrm{2}}$, for example the spatial correlation of $X$ and ${I}^{d}$. In this case
where
This is because there are ${N}_{{C}_{1}}^{T}$ ways to choose a cell of type ${C}_{1}$ and ${N}_{{C}_{2}}^{T}$ ways to choose a cell of type ${C}_{2}$. There are ${\mathrm{\Pi}}_{\mathit{\bm{X}}}^{H}{\mathrm{\Pi}}_{\mathit{\bm{X}}}^{L}$ possible positions for the first cell of type ${C}_{1}$ and ${\mathrm{\Pi}}_{\mathit{\bm{X}}}^{H}{\mathrm{\Pi}}_{\mathit{\bm{X}}}^{L}$ 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$.
${C}_{1}\ne {C}_{2}$, ${C}_{1},{C}_{2}$ lie on domain $D$ (with volume exclusion on $D$). For example spatial correlation of $X$ with respect to ${X}^{b}$.
${C}_{1}\ne {C}_{2}$, ${C}_{1},{C}_{2}$ lie on domains ${D}_{1}$, ${D}_{2}$ where ${D}_{1}\ne {D}_{2}$ for example spatial correlation of $X$ with respect to $M$.
Appendix 4
Testing robustness
Due to the abundance of parameters and cell–cell interactions necessary to capture what is known biologically about zebrafish pigment pattern formation, it is not feasible to perform an exhaustive parameter sweep to demonstrate the robustness of the model. Instead, as a test of robustness, we perform a rigorous robustness analysis by carrying out one hundred repeats of all WT and mutant simulations with perturbed parameter values chosen uniformly at random from the range 0.75–1.25 of their described value. The value of each parameter is sampled uniformly from this region, independentally for each parameter and each repeat. Ten of these randomly sampled repeats for shd, nac, pfe and leo are given in Appendix 4—figure 1. We observe, for all one hundred repeats, that small perturbations to the rates still generate consistent patterning, demonstrating the robustness of the model.
Appendix 5
Predicting pattern formation for an seurat/sbrcross
In this section, we give an example of how the model can be manipulated to investigate aspects of pattern development and to predict the outcome of mutant pigment patterns.
To the best of our knowledge, the adult pattern of crossed seurat and sbr mutants have not have been published previously. By changing the parameters so that we match both sbr and seurat we can predict pattern formation for a double mutant seurat/sbr. As previously stated (in more detail) in Appendix 2, the mutation in seurat affects melanocyte differentiation (Eom et al., 2012) and the mutation in sbr affects the densetoloose Siridophore transition (Fadeev et al., 2015). To simulate the cross, we simply incorporate the effects of both mutations simultaneously. The results at stage J+ is given in Appendix 5—figure 1. Our model predicts that when both of these mutations occur, dense Siridophores and associated xanthophores would cover most of the flank of the fish. A few melanocytes associated with loose Siridophores survive at the very dorsal and ventral regions of the fish.
References

1
Zebrafish leopard gene as a component of the putative reactiondiffusion systemMechanisms of Development 89:87–92.https://doi.org/10.1016/S09254773(99)002117
 2

3
How does cellular contact affect differentiation mediated pattern formation?Bulletin of Mathematical Biology 73:1529–1558.https://doi.org/10.1007/s1153801095784
 4
 5
 6

7
Pigment patterns in adult fish result from superimposition of two largely independent pigmentation mechanismsPigment Cell & Melanoma Research 28:196–209.https://doi.org/10.1111/pcmr.12335
 8
 9
 10
 11
 12
 13
 14
 15
 16
 17

18
Exact stochastic simulation of coupled chemical reactionsThe Journal of Physical Chemistry 81:2340–2361.https://doi.org/10.1021/j100540a008
 19

20
Pigment cell organization in the hypodermis of zebrafishDevelopmental Dynamics 227:497–503.https://doi.org/10.1002/dvdy.10334
 21

22
Tetraspanin 3c requirement for pigment cell interactions and boundary formation in zebrafish adult pigment stripesPigment Cell & Melanoma Research 27:190–200.https://doi.org/10.1111/pcmr.12192
 23
 24

25
P. zebrafish pigmentation mutations and the processes of neural crest developmentDevelopment 123:369–389.

26
Stripes and bellyspots  a review of pigment cell morphogenesis in vertebratesSeminars in Cell & Developmental Biology 20:90–104.https://doi.org/10.1016/j.semcdb.2008.10.001

27
Zebrafish adult pigment stem cells are multipotent and form pigment cells by a progressive fate restriction process: clonal analysis identifies shared origin of all pigment cell typesBioEssays : News and Reviews in Molecular, Cellular and Developmental Biology 39:201600234.https://doi.org/10.1002/bies.201600234

28
An updated kernelbased turing model for studying the mechanisms of biological pattern formationJournal of Theoretical Biology 414:120–127.https://doi.org/10.1016/j.jtbi.2016.11.003
 29

30
Black, yellow, or silver: which one leads skin pattern formation?Pigment Cell & Melanoma Research 28:2–4.https://doi.org/10.1111/pcmr.12328
 31

32
Nacre encodes a zebrafish MicrophthalmiaRelated protein that regulates NeuralCrestDerived pigment cell fateDevelopment 126:3757–3767.
 33
 34
 35
 36
 37
 38

39
Pyewacket, a new zebrafish fin pigment pattern mutantPigment Cell Research 19:232–238.https://doi.org/10.1111/j.16000749.2006.00311.x

40
Turing patterns: how the fish got its spotsPigment Cell & Melanoma Research 24:12–14.https://doi.org/10.1111/j.1755148X.2010.00814.x
 41
 42

43
A nonlocal model for contact attraction and repulsion in heterogeneous cell populationsBulletin of Mathematical Biology 77:1132–1165.https://doi.org/10.1007/s115380150080x
 44

45
An orthologue of the kitrelated gene fms is required for development of neural crestderived xanthophores and a subpopulation of adult melanocytes in the zebrafish, Danio rerioDevelopment 127:3031–3044.
 46

47
Normal table of postembryonic zebrafish development: staging by externally visible anatomy of the living fishDevelopmental Dynamics 238:2975–3015.https://doi.org/10.1002/dvdy.22113

48
Zebrafish Puma mutant decouples pigment pattern and somatic metamorphosisDevelopmental Biology 256:242–257.https://doi.org/10.1016/S00121606(03)000150
 49
 50
 51
 52
 53
 54
 55

56
Making digit patterns in the vertebrate limbNature Reviews Molecular Cell Biology 7:45–53.https://doi.org/10.1038/nrm1830

57
Modelling stripe formation in zebrafish: an agentbased approachJournal of the Royal Society Interface 12:20150812.https://doi.org/10.1098/rsif.2015.0812
 58
 59
 60
 61
 62
 63
Decision letter

Sandeep KrishnaReviewing Editor; National Centre for Biological Sciences‐Tata Institute of Fundamental Research, India

Aleksandra M WalczakSenior Editor; École Normale Supérieure, France

Sreelaja NairReviewer; Tata Institute of Fundamental Research, India

Raj K LadherReviewer; NCBS, India
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
This paper constructs and analyses a detailed computational model of zebrafish pigment stripe formation. It highlights the role of iridophores in organizing the melanophores and xanthophores, a facet that has largely been unaccounted for in computational modeling of zebrafish stripe patterns. Of particular interest is the analysis of the schachbrett mutant, which points to the role of tight junctions in the transition of Siridophore packing.
Decision letter after peer review:
Thank you for submitting your article "A quantitative model for zebrafish pattern formation" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Aleksandra Walczak as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Sreelaja Nair (Reviewer #1); Raj K Ladher (Reviewer #3).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
The manuscript describes a computational model that attempts to capture the complexities of xanthophore and iridophore interactions to explain zebrafish pigment stripe pattern formation. While the reviewers appreciate that the model is sufficient to reproduce known stripe patterns, it has not been demonstrated that the interactions included are necessary, and all reviewers were of the opinion that the model has not been used to make sufficiently interesting predictions. In order for a revision to be publishable, it must distinguish the work sufficiently from the work of Volkening and Sandstede by including several more predictions (see points 1,2,3 of reviewer 1, and point 3 of reviewer 3). In addition, as pointed out by reviewer 2, it is essential that the existing and added predictions be shown to be robust to changing parameter values within biologically reasonable ranges – tuning parameters one by one, keeping all others fixed, only between low and high values, is not a sufficient exploration of the parameter space to prove such robustness. A more comprehensive exploration of parameter space is also essential to satisfy reviewer 2 that the interactions included in the model are indeed necessary as claimed in the manuscript, and that despite the very large number of parameters the behaviour of the model is constrained enough to make robust predictions. The manuscript should also tone down claims of the novelty of the methods and the superiority of this approach compared to reactiondiffusion approaches – in particular, as reviewer 2 points out, stochasticity is easily included in reactiondiffusion approaches so that cannot count as a reason to favour a latticebased model.
Reviewer #1:
In the manuscript entitled "A quantitative modeling approach to zebrafish pigment pattern formation", Owen et al. describe in silico modeling strategies to explain zebrafish pigment stripe formation. The authors use onlattice modeling, in which five pigment cell types that are necessary to form the adult zebrafish pigment stripe pattern are treated as fixed to a mutually exclusive space within the domain in which the pattern will emerge. Onto this basic framework, the authors model known behaviors of the pigment cell types and intercellular interactions to simulate emergence of stripe pattern. Their modeling approach incorporates the role of iridophores in organizing the melanophores and xanthophores, a facet that has largely been unaccounted for in computational modeling of zebrafish pigment stripe pattern formation. The authors validate their model by recapitulating the wild type stripe pattern and test the model by simulating stripe pattern formation in the absence of one or more of the cell types/interactions. These additional simulations very nicely recreate stripe patterns observed in known zebrafish mutants. Overall the study is well carried out and the authors comprehensively validate their modeling strategy, making a strong case for onlattice agent based modeling rather than a Reaction Diffusion modeling for pigment stripe pattern in zebrafish.
It is not surprising that the model recapitulates known pigment patterns, since the model was built based on biological evidences. On the basis of this the authors state that the current experimental evidences together with their model’s assumptions of Siridophore behavior is sufficient to explain pigment pattern formation. For a set of biological interactions to be termed as sufficient for formation of a pattern, the model could perhaps be challenged by stating predictions that are currently not easy to test experimentally.
1) For example: As the authors point out, zebrafish do not incorporate additional pigment stripes as they grow. The biological basis for this is not known. In their wild type simulation starting at stage PB, how would a pattern evolve if initially the melanophores were in 3 or 5 stripes instead of the normal 4 stripes?
2) At the start of the simulations at stage PB, the dense iridophores appear along the horizontal myoseptum. Is this spatial location of iridophores relevant to formation of the pattern? Would the pattern that emerges be different if the original metamorphic pattern of dense iridophores at the horizontal myoseptum was displaced dorsally or ventrally?
3) The authors very nicely factor in the layered arrangement of the different cell types along the zaxis, to mirror the biological scenario. It is experimentally tough to switch the order of the layers, but should be possible to do in the simulations? This would be interesting because developmentally the three pigment cell types take distinct dorsal or ventral migratory routes along the horizontal myoseptum to reach their eventual zlayer. The route of migration has been thought to be critical in determining the fate of the pigment cell type (or of the stem cells that resides in the adult) and in establishing the stripe pattern. Could this assumption be partially tested to determine the relevance of the spatial order of the layers in eventual stripe formation?
Overall, I enjoyed the manuscript and the insights mathematical modeling can provide to understand complex phenomena such as pattern formation. My comments stem from limitations that mutants etc pose for understanding what information is sufficient to generate a pattern. The model is based on a finite set of proven interactions, several untested permutations and combinations of these and additional novel interactions are always possible. A major advantage of mathematical modeling is in the predictive zone, some of which experimentalists may find hard/impossible to venture into. The current study falls short of taking that leap, which would have been an interesting and informative exercise.
Reviewer #2:
This manuscript reports the study of pigment patterns in zebrafish embryos claiming a novel bottomup stochastic model. My understanding is that this is an extensively employed standard approach using the chemical master equation for interacting (and diffusing) species to simulate emergent stochastic patterns. The authors incorporate many detailed and complex interactions between five different cell types and claim that all the interactions are essential to explain the observed patterns. My feeling is that such a study with a large number of control parameters is overkill and does not enhance our understanding of stripe patterns. As such, I do not recommend publishing this study in eLife.
1) In this manuscript, the authors present a quantitative model for simulating the stripe patterns seen in wt and mutant embryos of zebrafish. They develop a stochastic description based on the purported interactions between xanthophores, xanthoblasts, melanocytes and (two kinds of) iridophores while also incorporating motility that is biased by short and longranged interactions between the cell types. These are augmented by cell division, death, differentiation and tissue growth. The resulting detailed numerical simulations are shown to qualitatively reproduce pigment patterns in wt and mutant conditions. Further, the authors quantitatively compare cell density, straightness of stripes and interstripe widths with experimental results. Several other results, such as the paircorrelationfunction of cell densities, provide further justifications to the model developed.
2) The main point of this study seems to be the development of a bottomup stochastic model incorporating as many complex and varied interactions between the cells to simulate pigment patterns. The authors contrast their approach with reactiondiffusion modelling approaches (such as Turing patterns) by saying that such approaches lack the stochasticity observed inherently in pigment patterns.
I would disagree with this viewpoint. The point of reactiondiffusion approaches is minimalistic coarsegrained descriptions of the complex effective interactions between cells (or any other constituents of the patterns) that can capture the essence of the emergent patterns. It is easy to incorporate stochasticity even in approaches that use partialdifferential equations.
A complementary approach is to simulate the underlying master equation that governs the chemical interactions, supplemented by a compartmentbased approach to diffusion. There are several such studies done in the past using the standard Gillespie algorithm for simulating the chemical master equation. Such studies incorporate the essential chemistry (e.g., activatorinhibitor dynamics) and coupled with diffusion can lead to emergent stripe patterns with inherent stochasticity builtin.
3) The authors, instead, develop a detailed latticebased model with sufficiently complex interactions between 5 species of cells along with shortranged local interactions and longranged nonlocal interactions. These interactions bias the hopping rates of the cells along the lattice. The authors speculate on the possible physical origins of these interactions. However, it is not clear if there aren't other possible interactions that could also lead to the same emergent dynamics. And hence the questions arises as to what are the essential features required at this level of modeling and what details are incidental. It should be noted that even though the authors claim that their approach is bottomup, they are not modelling interactions between molecules. Rather, the interactions incorporated in their modeling are also effective celllevel interactions. In essence, this study has considered a masterequation like approach to simulate zebrafish stripe patterns.
4) The authors claim that all the interactions included in their model are essential for stripe patterns. This is demonstrated by turning off each interaction inturn and showing that some feature of the pattern is lost. This is an 'all or none' switch. For the set of other parameters chosen, this might indeed lead to the conclusion that all interactions are absolutely needed. One could vary the interaction strengths in a continuous manner and then it is not clear if all the said interactions are absolutely needed for all parameter values.
5) Looking at the appendix and the supplementary material, I feel that the model is sufficiently complicated, and has so many turning parameters, that any range of behaviors is possible. It is not unreasonable to see that the largescale emergent dynamics of such a complex model is essentially that of Turing systems which could have been simulated with much less complexity.
For the above reasons, I cannot recommend publishing in eLife.
Reviewer #3:
Owen et al. describe the formulation of a model to understand the formation of stripes in zebrafish. They revisit the Turing model that previous studies have described and refine the interaction parameters that were simplified in those models. The model presented recapitulates many of the patterns found in WT and in mutants, and models the formation of these stripes.
The study does capture the complexities of xanthophore and iridophore interactions well. Furthermore, the ability to place weights on the strength of interactions also gives the model flexibility.
The paper is wellwritten. However, I would like to make some suggestions, listed below.
1) The assumptions and the way the model works could be made more explicit. For this I would suggest that the authors consider incorporating Figure 1, 2 and 3 from Appendix 1 into the main manuscript.
2) I would suggest that the authors emphasise the differences of their model from the Volkening paper published in Aug 2018. I do take issue with the authors characterisation of this paper as very recent, and would suggest that they make greater reference to it in the Introduction.
3) One very interesting piece of data from this study, and the one that does differentiates this study from that of Volkening and Sandstede, is the analysis of the schachbrett mutant. This mutant points to the role of tight junctions in the transition of Siridophore packing. It would be worth extending this analysis to the compunctiond sbr/leo and sbr/luc mutants described in Fadeev et al. Additionally, there is a suggestion of an interaction between sbr and seurat – I wonder if the model were able to predict what a double sbr/seurat would look like?
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Thank you for resubmitting your article "A quantitative model for zebrafish pattern formation" for consideration by eLife. Your article has been rereviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Aleksandra Walczak as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Sreelaja Nair (Reviewer #1); Raj K Ladher (Reviewer #3).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
We would like to draw your attention to changes in our revision policy that we have made in response to COVID19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.
The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:
1) In the revised manuscript, you write:
"The most commonly used mathematical paradigm for stripe formation takes the form of a Turing reaction diffusion model. In these representations, melanocytes and xanthophores diffuse and interact via a few long and short range 'reactions'. This class of model typically rely on a small number of parameters which, upon being altered, can generate a diverse range of patterns. Simplified models such as these have the benefit that they are often analytically tractable, allowing a deep understanding of the model. However, their main limitation is that, due to the simplicity of the approach, there is often no consistent way to link the parameters with measurable data, making it difficult to relate the model results back to the biology."
This gives the impression that (i) Turing pattern models are always simplified, and analytically tractable, and (ii) Turing models face a difficulty in relating their parameters to measurable data. Both are untenable statements and should be removed – there are several examples of both highly nonlinear Turingtype or reactiondiffusion models for pattern formation that are neither simple nor analytically tractable, and in several cases now people have shown how to relate their parameters to measurable data.
2) The authors write:
"Turing reactiondiffusion type models posit that combinations of short and long range dynamics between melanocytes and xanthophores generate stripe patterns. Indeed, much of the excitement around such models is the ease with which small parameter value changes result in diverse patterns, many readily recognisable from nature. […] A major difference between our model and Turing reaction diffusion models is that small parameter changes in our model do not typically generate qualitatively different patterns, whereas Turing reaction diffusion models can show substantial pattern changes in response to small changes (Budi, Patterson and Parichy, 2011; Yamaguchi, Yoshimoto and Kondo, 2007)."
This gives the impression that large changes in behaviour due to small changes in parameters is a generic feature of Turing models, as opposed to the model in this manuscript. This is not true, and the statements implying this should be removed – it is quite possible for Turing and reactiondiffusion models to produce small changes upon small changes in parameters, and equally it is possible for stochastic lattice models to exhibit large changes in behaviour upon small changes in parameters – it depends on the interactions and nonlinearities included in the model.
3) Finally:
"From an analytical perspective, a significant advantage of our onlattice model in contrast to offlattice models is their amenability to the derivation of a continuum model. Our model therefore opens up the opportunity for future exploratory work using a continuum model for mutants pfe and nac in order to explore whether pattern formation in these cases individually can be described as Turing patterns and to determine parameter ranges for successful pattern formation."
Here you argue, and we agree, that the model in the manuscript is in fact simply a discretized version of a reactiondiffusion model (albeit a quite complex one with many dynamical variables and many interactions). Thus, there is no great difference in approach between the discretized stochastic model you analyze and reactiondiffusion models, so please remove all statements that imply a large difference in your approach vs. reactiondiffusion or Turing type models. Further, a continuum limit can easily be constructed for offlattice models as well, so this is not an advantage of your approach over that of Volkening et al. Please remove this statement, and present the opportunity for building continuum models as one that applies to both your model as well Volkening et al.'s model.
https://doi.org/10.7554/eLife.52998.sa1Author response
The manuscript describes a computational model that attempts to capture the complexities of xanthophore and iridophore interactions to explain zebrafish pigment stripe pattern formation. While the reviewers appreciate that the model is sufficient to reproduce known stripe patterns, it has not been demonstrated that the interactions included are necessary, and all reviewers were of the opinion that the model has not been used to make sufficiently interesting predictions. In order for a revision to be publishable, it must distinguish the work sufficiently from the work of Volkening and Sandstede by including several more predictions (see points 1,2,3 of reviewer 1, and point 3 of reviewer 3).
We have extended our work substantially to include a number of new predictions as suggested by the reviewers, including:
1) Predicting the potential effects of changing the initial conditions for stripe formation
a) by changing the initial iridophore horizontal striping to be vertical (see subsection “Initial Siridophore interstripe orientation alone does not determine the orientation of stripes and interstripes and Figure 13A);
b) by displacing the initial horizontal interstripe vertically, so that it appears lower down the initial domain (see subsection “The position of the initial Siridophore interstripe is important for successful pattern formation” and Figure 13B);
c) by changing the initial domain size while maintaining the usual growth rates (see subsection “Initial domain size contributes to the number of stripes and interstripes” and Figure 13C);
d) by populating the initial domain with preformed stripes and observing stripe insertion (see subsection “Initial domain size contributes to the number of stripes and interstripes” and Figure 13D).
2) Formulating and then testing hypotheses about the functions of the leo gene by replicating the range of mutant and cross mutant patterns (see subsection “An in silico investigation into the function of the leo gene” and Figure 14)
3) Predicting the patterning of seurat mutants and double mutants based on the known biology. These extensions are included in the subsection “Model predictions” of the main text along with our original predictions and (for seurat mutants) in the Discussion and Appendix 5.
In addition, as pointed out by reviewer 2, it is essential that the existing and added predictions be shown to be robust to changing parameter values within biologically reasonable ranges – tuning parameters one by one, keeping all others fixed, only between low and high values, is not a sufficient exploration of the parameter space to prove such robustness.
We performed a rigorous robustness analysis by carrying out one hundred repeats for simulations of each of the WT and four of the main mutant types – nac, pfe, shd and leo – with perturbed parameter values chosen uniformly at random from the range between 0.751.25 of their described value. The process is described in detail in subsection “Necessity of Siridophore assumptions”. Twenty (randomly selected) of the repeats for WT are displayed in the new Figure 9 of the manuscript and ten (randomly selected) of the repeats for each of nac, pfe, shd and leo are given in Appendix 4—figure 1. For all one hundred repeats these significant perturbations to the rates still generate consistent patterning, in both WT and all mutant cases, demonstrating the robustness of the model.
A more comprehensive exploration of parameter space is also essential to satisfy reviewer 2 that the interactions included in the model are indeed necessary as claimed in the manuscript, and that despite the very large number of parameters the behaviour of the model is constrained enough to make robust predictions.
Most of the features we have included in the model have been demonstrated through biological experimentation and therefore are incorporated for completeness. A small number of the interactions included in the model are predictions focussed on the role of iridophores. We tested the necessity of these predictions by removing the corresponding model feature and demonstrating that the reduced model fails to replicate either the WT and/or mutant pattern predictions (Figure 11 of the revised manuscript). We further tested our predictions by applying our model to other mutants such as choker and sbr in order to demonstrate that these mutants are also successfully replicated by the model without the need for further assumptions (Figure 7OQ’ of the revised manuscript). We have also carried out a more comprehensive test of the biologically suggested features demonstrating that the removal of each model feature independently leads to incorrect pattern formation in the model (Figure 12). This acts as a clear demonstration of the necessity of all the interactions included in the model.
The manuscript should also tone down claims of the novelty of the methods and the superiority of this approach compared to reactiondiffusion approaches – in particular, as reviewer 2 points out, stochasticity is easily included in reactiondiffusion approaches so that cannot count as a reason to favour a latticebased model.
We have removed the suggestion that Turing models cannot incorporate stochasticity and toneddown claims of the superiority of the approaches in comparison to reactiondiffusion methods in the Introduction. Instead we have focussed on the important differences between the Turing reactiondiffusion modelling approach and our more biologically interpretable modelling approach.
Reviewer #1:
[…] It is not surprising that the model recapitulates known pigment patterns, since the model was built based on biological evidences. On the basis of this the authors state that the current experimental evidences together with their model`s assumptions of Siridophore behavior is sufficient to explain pigment pattern formation. For a set of biological interactions to be termed as sufficient for formation of a pattern, the model could perhaps be challenged by stating predictions that are currently not easy to test experimentally.
1) For example: As the authors point out, zebrafish do not incorporate additional pigment stripes as they grow. The biological basis for this is not known. In their wild type simulation starting at stage PB, how would a pattern evolve if initially the melanophores were in 3 or 5 stripes instead of the normal 4 stripes?
It appears that the reviewer may have slightly misinterpreted our initial conditions. Initially, within our simulations, the melanophores are not in a striped arrangement, instead they appear dynamically across the domain, in a randomly dispersed pattern (Figure 4A’), as reported in real zebrafish (see Figure 4A). Whilst in real WT zebrafish shown in Figure 4A there are some melanocytes in a stripeorientation that can be observed under the skin, these are much deeper below the skin layer and do not contribute to pattern development. We commented on this in the caption of Figure 4: “White arrows indicate the embryonic pattern of melanocytes in four stripes that are deeper than the skin level and are consequently not included in the model.”
Furthermore, we note that zebrafish do incorporate additional pigment stripes as they grow. However, it is important to note that these new stripes are added at the leading edges of a spreading pattern rather than between preformed stripes as has been documented, for instance, in the angelfish genus Pomocanthus.
Nevertheless, in order to investigate the potential for stripe insertion in zebrafish, we simulated a domain which is initially striped (two stripes and one interstripe) at stage PB. The resulting patterns demonstrate, as the reviewer suggested, that in this case new interstripes appear inbetween already developed stripes, consistent with the Pomocanthus study we now reference. We have illustrated our results in Figure 13D and provided an associated commentary in the subsection “Stripe insertion can occur on an initially striped domain”.
2) At the start of the simulations at stage PB, the dense iridophores appear along the horizontal myoseptum. Is this spatial location of iridophores relevant to formation of the pattern? Would the pattern that emerges be different if the original metamorphic pattern of dense iridophores at the horizontal myoseptum was displaced dorsally or ventrally?
We have addressed the reviewer’s interesting suggestion by exploring what happens if, hypothetically, the initial iridophore stripe was in a different place. Given how we envisage the role of that initial iridophore stripe in the organisation of the melanophores, we predicted that, if we were to move the position of the initial interstripe, but maintain its horizontal orientation, that the stripes would still form sequentially around this initial interstripe, but in a pattern displaced along the dorsoventral axis. We tested this hypothesis in the model and confirmed that this is exactly what happens. We have included simulations of pattern development from just such an initial condition in Figure 13B, discussing the similarities and differences to biologically plausible wildtype patterning in subsection “The position of the initial Siridophore interstripe is important for successful pattern formation”. The stripes are more disrupted than with the wildtype initial condition. This effect is likely due to the altered impact of domain growth on the pattern: In the model, growth is centred at the middle of the domain and so when the initial stripe is not similarly centred, growth disrupts pattern formation.
Another interesting test case inspired by the reviewer’s comments, concerns the implication for pattern formation of altering the orientation of the initial iridophore stripe, so that it is oriented vertically (rather than horizontally as in wildtype fish). Based on the assumed role of the iridophore stripe in orienting the stripe pattern, we predicted that this change would result in vertical barring. In the simulation in Figure 13A and commentary in subsection “Initial Siridophore interstripe orientation alone does not determine the orientation of stripes and interstripes”, it can be seen that, whilst initially the vertical iridophores drive a vertically oriented banding pattern, this gradually reorganises, so that at the end of the simulation, the pattern is a meshwork of vertical and horizontal striping. This result clearly supports the idea that the initial iridophore pattern is important in establishing stripe/interstripe orientation, but importantly reveals the previously unsuspected importance of domain growth in orienting striping in zebrafish.
3) The authors very nicely factor in the layered arrangement of the different cell types along the zaxis, to mirror the biological scenario. It is experimentally tough to switch the order of the layers, but should be possible to do in the simulations? This would be interesting because developmentally the three pigment cell types take distinct dorsal or ventral migratory routes along the horizontal myoseptum to reach their eventual zlayer. The route of migration has been thought to be critical in determining the fate of the pigment cell type (or of the stem cells that resides in the adult) and in establishing the stripe pattern. Could this assumption be partially tested to determine the relevance of the spatial order of the layers in eventual stripe formation?
The reviewer raises an interesting question: What determines the fate of pigment cells generated from the adult stem cells? To our knowledge there has been little study of the way in which the layered arrangements are built up over juvenile stages. For this reason, our model does not capture the migration pathway to the epidermis, but only the patterning process once the cells have arrived. Cells are modelled as interacting without explicit reference to layering. Thus, the order of the layers does not influence the output of our model. The distinct layers are incorporated in order to aptly represent the role of excluded volume between cells of the same type. That is, melanophores should occupy space available to other melanophores, but not to xanthophores for example, since they occupy a separate zlayer. We have added the sentence ‘We note that in our simulations, the ordering of the layers does not play any role in determining pattern formation.’ to the manuscript in the caption of Figure 2 in order to make this clear.
Reviewer #2:
This manuscript reports the study of pigment patterns in zebrafish embryos claiming a novel bottomup stochastic model. My understanding is that this is an extensively employed standard approach using the chemical master equation for interacting (and diffusing) species to simulate emergent stochastic patterns. The authors incorporate many detailed and complex interactions between five different cell types and claim that all the interactions are essential to explain the observed patterns. My feeling is that such a study with a large number of control parameters is overkill and does not enhance our understanding of stripe patterns. As such, I do not recommend publishing this study in eLife.
We feel that the reviewer has perhaps misunderstood the motivation for our paper. As they have stated there are many demonstrations in the literature, including those we cite, which show that diverse patterns resembling those in wildtype and mutant zebrafish can be reproduced by a reactiondiffusion master equation (RDME) approach. Our aim, rather than asking if reactiondiffusion mechanisms could generate a pattern, was to ask whether the cellcell interactions that have already been characterised biologically in the literature are necessary and sufficient to generate the observed biological patterning. Moreover, we investigate whether the resulting model can account for the diverse array of patterns seen in a variety of mutants, if the known biology corresponding to these mutants is incorporated. Though we agree that the number of parameters in this model is large, we argue that this reflects the demonstrated complexity of the underlying biology. Moreover, the iridophore interactions that are incorporated are based on hypotheses derived from biological studies; they have not been chosen for modelling convenience or to replicate particular patterns. Our model then provides a tool to determine whether these hypotheses are consistent with wildtype and mutant pattern formation.
1) In this manuscript, the authors present a quantitative model for simulating the stripe patterns seen in wt and mutant embryos of zebrafish. They develop a stochastic description based on the purported interactions between xanthophores, xanthoblasts, melanocytes and (two kinds of) iridophores while also incorporating motility that is biased by short and longranged interactions between the cell types. These are augmented by cell division, death, differentiation and tissue growth. The resulting detailed numerical simulations are shown to qualitatively reproduce pigment patterns in wt and mutant conditions. Further, the authors quantitatively compare cell density, straightness of stripes and interstripe widths with experimental results. Several other results, such as the paircorrelationfunction of cell densities, provide further justifications to the model developed.
2) The main point of this study seems to be the development of a bottomup stochastic model incorporating as many complex and varied interactions between the cells to simulate pigment patterns. The authors contrast their approach with reactiondiffusion modelling approaches (such as Turing patterns) by saying that such approaches lack the stochasticity observed inherently in pigment patterns.
I would disagree with this viewpoint. The point of reactiondiffusion approaches is minimalistic coarsegrained descriptions of the complex effective interactions between cells (or any other constituents of the patterns) that can capture the essence of the emergent patterns. It is easy to incorporate stochasticity even in approaches that use partialdifferential equations.
Whilst we acknowledge that Turing patterns are an excellent tool for inspiring scientists to investigate the mechanisms behind pattern formation in a variety of biological systems, we believe they are not the most appropriate tools for understanding the detailed mechanisms underlying zebrafish pigmentation patterning at the microscale.
A lack of stochasticity in classical partial differential equation models of Turing reactiondiffusion mechanisms is just one reason why we view such models as less appropriate for this inherently stochastic pattern formation process but is of relatively minor concern in comparison to the other drawbacks we perceive these models have. That being said, we acknowledge that reactiondiffusion models can incorporate stochasticity and we have now removed this comment from the manuscript.
More troubling, in our point of view, is the lack of cellular level resolution that such continuum models provide. In an age in which computational power is not an issue and celllevel biological data is available, individualbased models are an appropriate tool for investigating biological hypotheses made at a cellular level.
Even more concerning to us is the lack of detailed insight provided by coarsegrained Turing models. It is precisely the coarsegrained ‘effective’ descriptions of the complex interactions that the reviewer describes, which rob Turing models of zebrafish pigment pattern formation of their predictive power. It is often impossible, in such models, to link parameters specifically to the microscopic individuallevel data, especially because many different types of reaction kinetics are capable of generating the same coursegrained patterns. Consequently, Turing models are often unable to provide mechanistic insight into the underlying biological processes, something which our detailed, individualbased model is manifestly capable of doing.
Perhaps the starkest illustration of this limitation is that Turing models of zebrafish pigmentation patterns are capable of producing stripes with just two cell species – melanophores and xanthophores – whereas the real biology requires iridophores as well. Early twospecies reactiondiffusion models of zebrafish pigmentation pattern formation gave no hint of the importance of the missing cell type or the necessity of the corresponding interactions, whereas detailed individualbased models such as ours, incorporating all the known biology, surely would have. That a model can produce a pattern which looks qualitatively similar to a biological pattern is not sufficient to suggest it can provide biological insight.
In summary, whilst we do not dispute that Turing reactiondiffusion models are useful for suggesting potential pattern formation mechanisms in a variety of biological contexts, and indeed were what first drew one of us, and many others, to become interested in understanding patterning processes many years ago, we do not believe them the most appropriate tool to test the complex individualcellbased hypotheses that we investigate in this paper.
To address these comments, we have refined our comments about Turing mechanisms to the following:
“They [Turing mechanisms] typically rely on a small number of parameters which, upon being altered, can generate a diverse range of patterns. Simplified models such as these have the benefit that they are often analytically tractable, allowing a deep understanding of the pattern formation mechanism. However, their main limitation is that, due to the simplicity of the approach, there is often no consistent way to link the parameters with measurable data, making it difficult to relate the model results back to the biology.”
A complimentary approach is to simulate the underlying master equation that governs the chemical interactions, supplemented by a compartmentbased approach to diffusion. There are several such studies done in the past using the standard Gillespie algorithm for simulating the chemical master equation. Such studies incorporate the essential chemistry (e.g., activatorinhibitor dynamics) and coupled with diffusion can lead to emergent stripe patterns with inherent stochasticity builtin.
Compartmentbased reactiondiffusion simulations (often known as a reactiondiffusion master equation (RDME) formalism) are, as the reviewer points out, an alternative way to simulate reactiondiffusion processes which naturally incorporate stochasticity. Indeed, these models are capable of giving stripe formation consistent with the continuum models from which they are derived, but not, crucially, consistent with the biology underlying zebrafish pigmentation pattern formation for the following reasons:
1) In order to produce striped patterns, these RDME models must have multiple cells per compartment (as opposed to the single cell per compartment suggested by the volumeexcluding property exhibited by the real cells). Such models typically also rely on reaction kinetics which bear little or no correspondence to the underlying biological interactions between different types of cells.
2) Once antiphase patterns of peaks and troughs have been produced by such models, an appropriately tuned threshold must be chosen in order to decide which regions constitute “stripe” and which “interstripe”. Consequently, the width of the stripes that result from such models is not an inherent property of the model, but an arbitrary choice on behalf of the modeller. The arbitrary threshold beyond which one region is considered stripe and another interstripe also belies the fact that these models result in nontrivial densities of xanthophores in the stripes and nontrivial densities of melanophores in the interstripes. This fact often goes unmentioned in such models but represents a significant departure from the real biology of the system.
Our approach shares some similarity to the RDME approach but incorporates more biological realism. Indeed, we use the Gillespie algorithm to simulate the stochastic events in our individualbased model, but rather than incorporating “the essential chemistry (e.g., activatorinhibitor dynamics)” which corresponds to a phenomenological Turing reactiondiffusion model that the reviewer suggests (NB ‘essential chemistry’ here is a reference to a required feature for the modeltype, rather than an ‘essential known feature of the biological mechanism’), our detailed individualbased model begins with the biology. Consequently, we do not incur the limitations associated with Turing models that we have highlighted above.
3) The authors, instead, develop a detailed latticebased model with sufficiently complex interactions between 5 species of cells along with shortranged local interactions and longranged nonlocal interactions. These interactions bias the hopping rates of the cells along the lattice. The authors speculate on the possible physical origins of these interactions. However, it is not clear if there aren't other possible interactions that could also lead to the same emergent dynamics. And hence the questions arises as to what are the essential features required at this level of modeling and what details are incidental. It should be noted that even though the authors claim that their approach is bottomup, they are not modelling interactions between molecules. Rather, the interactions incorporated in their modeling are also effective celllevel interactions. In essence, this study has considered a masterequation like approach to simulate zebrafish stripe patterns.
A bottomup modelling approach does not require that modelling begins at the level of molecules (or, to go into even more detail, atoms), rather it should be at the lowest level at which pertinent details are available. Given the current state of biological knowledge, the level of cellcell interactions is the lowest level at which we have sufficient information to build a model capable of recapitulating zebrafish pigmentation patterns. As a result, our model is capable of investigating biological hypotheses about interactions at this level, which higherlevel, topdown models, focussed solely on pattern replication irrespective of the details of the underlying biology, would not be. Our model incorporates what is known of the biology and consequently provides a tool to investigate hypotheses about that which is unknown.
For the reasons we highlighted in the response to the reviewer’s previous point, our study is distinct from an RDME approach and avoids the biological unrealism associated with it. We consider finegrained, biologically informed celllevel interactions as opposed to the coarsegrained interactions employed by the RDME approach.
Whilst it is possible that other cellcell interaction dynamics might lead to similar emergent dynamics, the cellcell interactions we have chosen are either biologically evidenced or (in the case of missing iridophore dynamics) previously hypothesised. The origins of the known cellcell interactions are given in text [14]. Including alternative interactions with no basis in biological observation seems illogical.
4) The authors claim that all the interactions included in their model are essential for stripe patterns. This is demonstrated by turning off each interaction inturn and showing that some feature of the pattern is lost. This is an 'all or none' switch. For the set of other parameters chosen, this might indeed lead to the conclusion that all interactions are absolutely needed. One could vary the interaction strengths in a continuous manner and then it is not clear if all the said interactions are absolutely needed for all parameter values.
An exhaustive continuous parameter sweep over all model parameters is not feasible. Such a parameter sweep would also be inappropriate since the majority of parameters are determined biologically. However, investigating the broad range of patterns that the model can produce is certainly a question of biological interest which we intend to follow up in a subsequent project.
5) Looking at the appendix and the supplementary material, I feel that the model is sufficiently complicated, and has so many turning parameters, that any range of behaviors is possible. It is not unreasonable to see that the largescale emergent dynamics of such a complex model is essentially that of Turing systems which could have been simulated with much less complexity.
As we have previously stressed, the majority of the model’s parameters are biologically determined leaving the range of behaviours relatively restricted, although still capable of replicating the relevant biological patterns.
Without doubt, the basic striped patterns of wildtype zebrafish could be (and indeed has been) simulated with (simplified twovariable) Turing models. However, the point of our study is not simply to replicate patterns, but to provide mechanistic insight into the biology by which they are generated – something that Turing models are poorly equipped to do.
In order to demonstrate the robustness of our model to changes in the parameters we carried out a rigorous robustness analysis in which we varied the parameter values independently of each other. For 100 repeat simulations each parameter value was assigned at random from the range between 0.751.25 of their described value, effectively sampling the highdimensional parameter space randomly. The process is described in detail in subsection “Necessity of Siridophore assumptions” and twenty (randomly selected from the total of 100) repeats are displayed in Figure 9. It is clear to see, by looking at the figure, that the wildtype striped pattern is not significantly perturbed by any of the random parameters choice, demonstrating the robustness of the model to parameter changes. Note that we then repeated this robustness analysis for each of four of the key mutants, with similar confirmation of the robustness of the model’s outputs to reasonable parameter value changes. Twenty (randomly selected) of the repeats for WT are displayed in the new Figure 9 of the manuscript and ten (randomly selected) of the repeats for each of nac, pfe, shd and leo are given in Appendix 4—figure 1.
Reviewer #3:
[…] The paper is wellwritten. However, I would like to make some suggestions, listed below.
1) The assumptions and the way the model works could be made more explicit. For this I would suggest that the authors consider incorporating Figure 1, 2 and 3 from Appendix 1 into the main manuscript.
As suggested, we have added Appendix 1—figures 1, 2 and 3 to the main manuscript, reconstituting them as Figure 2 and Figure 3 to ease understanding. We thank the reviewer for this suggestion which has, we think, strengthened the manuscript by making it easier to understand.
2) I would suggest that the authors emphasise the differences of their model from the Volkening paper published in Aug 2018. I do take issue with the authors characterisation of this paper as very recent, and would suggest that the make greater reference to it in the Introduction.
We have changed ‘very recent’ to simply ‘recent’ where this appears. In the conclusion we had provided a paragraph addressing the similarities and differences between our work and Volkening’s, which we retain. To address the reviewers comment we have also added a sentence about the work of Volkening in the Introduction:
“These findings have paved the way for more detailed modelling, such as that of Volkening et al., 2018, who demonstrated (using an offlattice individualbased model) the need for understanding Siridophore behaviour when representing all three celltypes.”
3) One very interesting piece of data from this study, and the one that does differentiates this study from that of Volkening and Sandstede, is the analysis of the schachbrett mutant. This mutant points to the role of tight junctions in the transition of Siridophore packing. It would be worth extending this analysis to the compunctiond sbr/leo and sbr/luc mutants described in Fadeev et al.
Upon the reviewer’s suggestion we incorporated the known features of leo, leo;sbr, leo;pfe, leo;nac and leo;shd mutants into our model. We found that, upon the incorporation of suitable assumptions we could, in fact, generate all the known patterns of these single and double mutants. We have stated the assumptions that we required in order to generate this range of distinctive patterns as a series of testable hypotheses about the function of the leo gene and have demonstrated the patterns that result when different combinations of these hypotheses are not satisfied in the model. We have included this now in the manuscript (see subsection “The position of the initial Siridophore interstripe is important for successful pattern formation” and Figure 15). Additionally, we have tested the robustness of this pattern formation in Appendix 4—figure 1.
Additionally, there is a suggestion of an interaction between sbr and seurat – I wonder if the model were able to predict what a double sbr/seurat would look like?
The mutation in seurat affects melanocyte differentiation [3] and the mutation in sbr effects iridophore transition from dense to loose [5]. By changing the parameters to match this known biology we can effectively replicate the patterns associated with both sbr and seurat independently. Consequently, we can model and hence predict pattern formation for a double seurat/sbr mutant by changing the parameters so that the effects of both mutations occur simultaneously. We have added an example of a typical simulation result for a seurat/sbr at stage J+ in Appendix 5—figure 1. Our model predicts, when both of these mutations occur, that by stage J+, dense iridophores and xanthophores would cover most of the flank of the fish. A few melanophores associated with loose iridophores survive at the very dorsal and ventral regions of the fish. We have included a discussion of this double mutant along with the figure in in Appendix 5.
References:
1] Patterson, L. B., and Parichy, D. M. (2013). Interactions with Iridophores and the Tissue Environment Required for Patterning Melanophores and Xanthophores during Zebrafish Adult Pigment Stripe Formation. PLoS Genet., 9(5), 114.
2] Takahashi, G., and Kondo, S. (2008). Melanophores in the stripes of adult zebrafish do not have the nature to gather, but disperse when they have the space to move. Pigment Cell Melanoma Res., 21(6), 677–686.
3] Eom, D. S., Inoue, S., Patterson, L. B., Gordon, T. N., Slingwine, R., Kondo, S., Watanabe, M., and Parichy, D. M. (2012). Melanophore Migration and Survival during Zebrafish Adult Pigment Stripe Development Require the Immunoglobulin Superfamily Adhesion Molecule Igsf11. PLoS Genetics, 8(8), 116.
4] Singh, A. P., Schach, U., and NüssleinVolhard, C. (2014). Proliferation, dispersal and patterned aggregation of iridophores in the skin prefigure striped colouration of zebrafish. Nat. Cell Biol., 16(6), 607–614.
5] Fadeev, A., Krauss, J., Frohnhöfer, H. G., Irion, U., and NüssleinVolhard, C. (2015). Tightjunction protein 1a regulates pigment cell organisation during zebrafish colour patterning. eLife, 2015(4), 1–25.
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:
1) In the revised manuscript, you write:
"The most commonly used mathematical paradigm for stripe formation takes the form of a Turing reaction diffusion model. In these representations, melanocytes and xanthophores diffuse and interact via a few long and short range 'reactions'. This class of model typically rely on a small number of parameters which, upon being altered, can generate a diverse range of patterns. Simplified models such as these have the benefit that they are often analytically tractable, allowing a deep understanding of the model. However, their main limitation is that, due to the simplicity of the approach, there is often no consistent way to link the parameters with measurable data, making it difficult to relate the model results back to the biology."
This gives the impression that (i) Turing pattern models are always simplified, and analytically tractable, and (ii) Turing models face a difficulty in relating their parameters to measurable data. Both are untenable statements and should be removed – there are several examples of both highly nonlinear Turingtype or reactiondiffusion models for pattern formation that are neither simple nor analytically tractable, and in several cases now people have shown how to relate their parameters to measurable data.
We appreciate the reviewer’s considered points and consequently have rectified the potential ambiguity of our statements to reflect the reviewer’s concerns. We have changed the text as follows:
“The most commonly used mathematical paradigm for stripe formation takes the form of a Turing reaction diffusion model. […] However, a potential limitation is that parameters do not always have a clear biological interpretation which, can sometimes make it difficult to link parameters to measurable data.”
2) The authors write:
"Turing reactiondiffusion type models posit that combinations of short and long range dynamics between melanocytes and xanthophores generate stripe patterns. Indeed, much of the excitement around such models is the ease with which small parameter value changes result in diverse patterns, many readily recognisable from nature. […] A major difference between our model and Turing reaction diffusion models is that small parameter changes in our model do not typically generate qualitatively different patterns, whereas Turing reaction diffusion models can show substantial pattern changes in response to small changes (Budi, Patterson and Parichy, 2011; Yamaguchi, Yoshimoto and Kondo, 2007)."
This gives the impression that large changes in behaviour due to small changes in parameters is a generic feature of Turing models, as opposed to the model in this manuscript. This is not true, and the statements implying this should be removed – it is quite possible for Turing and reactiondiffusion models to produce small changes upon small changes in parameters, and equally it is possible for stochastic lattice models to exhibit large changes in behaviour upon small changes in parameters – it depends on the interactions and nonlinearities included in the model.
Whilst in some cases small parameter value changes specifically near bifurcation points can result in diverse patterns for Turing patterns (see Lee et al., 2012; Metz et al., 2011; Watanabe and Kondo, 2015)), we acknowledge that it is not a generic feature of all Turing models. We appreciate that these sentences gave a misleading impression and for this reason we have added ‘can sometimes’ to the first paragraph as follows:
"Turing reactiondiffusiontype models posit that combinations of short and long range dynamics between melanocytes and xanthophores generate stripe patterns.” Indeed, a lot of the excitement around such models is the ease with which small parameter value changes can sometimes result in diverse patterns, many readily recognisable from nature.
Furthermore we have modified the following sentence to be more specific:
“A major difference between our model and Turing reactiondiffusion models is that small parameter changes in our model do not typically generate qualitatively different patterns, whereas Turing reactiondiffusion models can show substantial pattern changes in response to small alterations to parameters near to bifurcation points (Budi, Patterson and Parichy, 2011; Yamaguchi, Yoshimoto and Kondo, 2007)."
3) Finally:
"From an analytical perspective, a significant advantage of our onlattice model in contrast to offlattice models is their amenability to the derivation of a continuum model. Our model therefore opens up the opportunity for future exploratory work using a continuum model for mutants pfe and nac in order to explore whether pattern formation in these cases individually can be described as Turing patterns and to
determine parameter ranges for successful pattern formation."
Here you argue, and we agree, that the model in the manuscript is in fact simply a discretized version of a reactiondiffusion model (albeit a quite complex one with many dynamical variables and many interactions). Thus, there is no great difference in approach between the discretized stochastic model you analyze and reactiondiffusion models, so please remove all statements that imply a large difference in your approach vs. reactiondiffusion or Turing type models.
We agree with the reviewer that we can likely derive a continuum reactiondiffusion model from our onlattice individual model, but we note that this will not necessarily be a Turing reactiondiffusion model. It is for this reason that we express interest in deriving the continuum model, i.e. to assess this. Furthermore, being able to derive a continuum model from an individualbased model does not mean that the individualbased model is simply a discretisation of the resulting continuum model. Undoubtedly, in complex individualbased models which employ nonlinear interactions, approximations including (but not limited to) moment closure assumptions will have to be made in order to reach the continuum limit, meaning that the continuum model will exhibit different behaviour to the individualbased model in certain parameter regimes. The advantage of deriving a continuum reactiondiffusion model from an individualbased model, rather than the other way around, is that we know exactly what assumptions have been made in passing from the individualbased model to the continuum model.
Further, a continuum limit can easily be constructed for offlattice models as well, so this is not an advantage of your approach over that of Volkening et al. Please remove this statement, and present the opportunity for building continuum models as one that applies to both your model as well Volkening et al.'s model.
We acknowledge and address this concern as follows:
"From an analytical perspective, an advantage of our onlattice model is its amenability to the derivation of a continuum model, although we note that continuum approximations to offlattice individualbased models can also be derived. Our model therefore opens up the opportunity for future exploratory work using a continuum model for mutants pfe and nac in order to explore whether pattern formation in these cases individually can be described as Turing patterns and to determine parameter ranges for successful pattern formation."
That it is also possible to derive continuum limits from offlattice models does not diminish the utility of the approach for our onlattice models. That said, we have moved this statement, so that it fits within our summary of future work instead of our comparison with Volkening et al. as we believe its previous positioning gave a misleading impression of the purpose of the statement. Having checked through the manuscript carefully, we do not see anywhere else where we claim that our model has this advantage over the Volkening model.
https://doi.org/10.7554/eLife.52998.sa2Article and author information
Author details
Funding
Biotechnology and Biological Sciences Research Council (SWBio DTP)
 Jennifer P Owen
Biotechnology and Biological Sciences Research Council (BB/ L00769X/1 (RNK))
 Robert N Kelsh
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Senior Editor
 Aleksandra M Walczak, École Normale Supérieure, France
Reviewing Editor
 Sandeep Krishna, National Centre for Biological Sciences‐Tata Institute of Fundamental Research, India
Reviewers
 Sreelaja Nair, Tata Institute of Fundamental Research, India
 Raj K Ladher, NCBS, India
Publication history
 Received: October 23, 2019
 Accepted: June 21, 2020
 Version of Record published: July 27, 2020 (version 1)
Copyright
© 2020, Owen et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 819
 Page views

 107
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Download citations (links to download the citations from this article in formats compatible with various reference manager tools)
Open citations (links to open the citations from this article in various online reference manager services)
Further reading

 Computational and Systems Biology
To characterize how proteinprotein interaction (PPI) networks change, we quantified the relative PPI abundance of 1.6 million protein pairs in the yeast Saccharomyces cerevisiae across nine growth conditions, with replication, for a total of 44 million measurements. Our multicondition screen identified 13,764 pairwise PPIs, a threefold increase over PPIs identified in one condition. A few ‘immutable’ PPIs are present across all conditions, while most ‘mutable’ PPIs are rarely observed. Immutable PPIs aggregate into highly connected ‘core’ network modules, with most network remodeling occurring within a loosely connected ‘accessory’ module. Mutable PPIs are less likely to coexpress, colocalize, and be explained by simple mass action kinetics, and more likely to contain proteins with intrinsically disordered regions, implying that environmentdependent association and binding is critical to cellular adaptation. Our results show that protein interactomes are larger than previously thought and contain highly dynamic regions that reorganize to drive or respond to cellular changes.

 Computational and Systems Biology
 Structural Biology and Molecular Biophysics
The way in which multidomain proteins fold has been a puzzling question for decades. Until now, the mechanisms and functions of domain interactions involved in multidomain protein folding have been obscure. Here, we develop structurebased models to investigate the folding and DNAbinding processes of the multidomain Yfamily DNA polymerase IV (DPO4). We uncover shifts in folding mechanism among ordered domainwise folding, backtracking folding, and cooperative folding, modulated by interdomain interactions. These lead to "Ushaped' folding kinetics. We characterize the effects of interdomain flexibility on the promotion of DPO4DNA (un)binding, which probably contributes to the ability of DPO4 to bypass DNA lesions, a known biological role of Yfamily polymerases. We suggest that the native topology of DPO4 leads to a tradeoff between fast, stable folding and tight functional DNA binding. Our approach provides an effective way to quantitatively correlate the roles of protein interactions in conformational dynamics at the multidomain level.