Modeling Stripe Formation on Growing Zebrafish Tailfins

Abstract

As zebrafish develop, black and gold stripes form across their skin due to the interactions of brightly colored pigment cells. These characteristic patterns emerge on the growing fish body, as well as on the anal and caudal fins. While wild-type stripes form parallel to a horizontal marker on the body, patterns on the tailfin gradually extend distally outward. Interestingly, several mutations lead to altered body patterns without affecting fin stripes. Through an exploratory modeling approach, our goal is to help better understand these differences between body and fin patterns. By adapting a prior agent-based model of cell interactions on the fish body, we present an in silico study of stripe development on tailfins. Our main result is a demonstration that two cell types can produce stripes on the caudal fin. We highlight several ways that bone rays, growth, and the body–fin interface may be involved in patterning, and we raise questions for future work related to pattern robustness.

Appendices

A Appendix: Parameters

With the exception of the two switch parameters [namely $\zeta$ in Eq. (12) and d in Eqs. (13)–(14)] that we use to implement Mechanism II and IV, we summarize all of the parameters involved in our rules for cell migration, birth, and death in Tables 2, 3, and 4, respectively. We give the simulation-specific values of the switch parameters by figure in “Appendix B.5.” We set $d = 150\,\upmu \hbox {m}$ when we include Mechanism II, and we set $\zeta = 1000$ cells when we test Mechanism IV.

Table 2 Summary of parameters for cell migration in Eqs. (9)–(12) and (18)

Table 3 Summary of our parameters for cell birth [see Eqs. (13)–(14)] and model length scales

Table 4 Summary of our parameters for cell death in Eqs. (15)–(16)

B Appendix: Simulation Conditions

We used MATLAB 9.3, The MathWorks, Inc., Natick, MA, USA, to simulate our model of cell interactions on growing fins. Our code is available from the corresponding author on request. We now describe our model implementation (“Appendix B.1”), our boundary conditions (“Appendix B.2”), our initial conditions (“Appendix B.3”), and our methods for selecting cell-birth locations in tailfin domains and implementing Mechanism V (“Appendix B.4”). In “Appendix B.5,” we summarize the parameters and simulation conditions associated with our simulations in Figs. 5, 6, and 7.

1.1 B.1 Model Implementation

Our simulation begins with an initial condition (see “Appendix B.3”) at $t=18$ dpf. We update cell positions from day t to day $t+1$ in seven steps:

• 1. Set the cycle counter c for migration and birth to $c=1$.

• 2. Increment $N^\text {M}_\text {diff}(t)$ and $N^\text {X}_\text {diff}(t)$ by 20 locations each, so that $N^\text {M}_\text {diff}(t) = n^\text {M}_\text {diff} + 20(t-18)$ locations and $N^\text {X}_\text {diff}(t) = n^\text {X}_\text {diff} + 20(t-18)$ locations.

• 3. Perform one step $\varDelta t_\text {mig,birth}$ of migration [e.g., solve Eqs. (9) and (11) using the forward Euler scheme with time step $\varDelta t_\text {mig,birth}$]. All of our cells migrate simultaneously. We also specify repulsive forces from the discretized fin boundaries at each step $\varDelta t_\text {mig,birth}$ of migration (see “Appendix B.2”).

• 4. Select $N^\text {M}_\text {diff}(t+1)$ and $N^\text {X}_\text {diff}(t+1)$ potential locations for M and X birth on the fin domain (outlined by our boundary curve at time $t+1$), respectively. Evaluate these locations (simultaneously) for cell birth based on Eqs. (13)–(14) and random birth if included (e.g., if $p_\text {M} >0$ and $p_\text {X} >0$). Add the newly born cells to the domain at time $t+c\varDelta t_\text {mig,birth}$.

• 5. If $c \varDelta t_\text {mig,birth} = 1$ day, one day of migration and birth has been completed: go to Step 6. Otherwise, increment the cycle counter c for migration and birth ($c= c+1$) and return to Step 3 with the cell positions at time $t+c\varDelta t_\text {mig,birth}$.

• 6. Evaluate all of the cells for possible death (simultaneously) through Eqs. (15)–(16). Remove any cells that have died from the domain at day $t+1$. We note that the time step for our cell death rules is always $\varDelta t_\text {death} =1$ day.

• 7. If uniform domain growth is included, scale the cell positions at time $t+1$ using the fin domains and bone rays at times $t+1$ and $t+2$ days as we describe in Sect. 3.1. If distal epithelial growth is included, do not scale the cell positions. The result of this process is the updated cell positions at time $t+1$ days.

1.2 B.2 Boundary Conditions

To help keep cells in the fin domain, we include wall-like boundary conditions. For each simulated day t, we discretize the associated fin boundary curve at time t into 500 points: $\{\mathbf{F}_i(t)\}_{i=1,\ldots ,500}$. We then specify repulsive forces from these boundary points to our cell agents at each time step of migration $\varDelta t_\text {mig,birth}$:

$$\begin{aligned} \text {boundary force on } i\text {th cell at position } \mathbf{C} _i = -\sum _{j=1}^{500} \triangledown Q^\text {bnd}(\mathbf{F} _j - \mathbf{C} _i) \end{aligned}$$

(17)

where $\mathbf{C} \in \{\mathbf{M}, \mathbf{X}\}$ and

$$\begin{aligned} Q^\text {bnd}(\mathbf{d} ) = R^\text {bnd} e^{-|\mathbf{d} |/r_\text {bnd}}. \end{aligned}$$

(18)

These rules are an approximation of Neumann boundary conditions, and we note that a small number of cells escape from our domain in some simulations. Cells are more likely to escape with increasing time, suggesting that it may be useful for future work to increase the number of boundary agents as the fin grows. When cells escape, we remove them from our final simulated images in post-processing.

1.3 B.3 Initial Condition

The initial condition for our simulations is motivated by images in Parichy et al. (2009). First, we specify a single horizontal strip of M cells (separated $30\,\upmu \hbox {m}$ apart) at the center of our fin domain (e.g., with y-coordinate 0). Second, for X cells, we consider a random distribution of cells that are concentrated more highly toward the proximal edge of the fin. We choose the y-coordinates of these positions by selecting 500 points uniformly at random between the maximum and minimum y-coordinates for the discretized boundary curve that represents our initial domain at 18 dpf (see “Appendix B.2”). We choose the x-coordinates for these points by taking the absolute value of 500 points sampled from a normal distribution with mean 0 and standard deviation $\sigma = 0.25 x_\text {max}(t)$, where $x_\text {max}(t)$ is the maximum x-coordinate for our discretized boundary curve at 18 dpf.

As the penultimate step in setting our initial condition, we remove any M or X cells that fall within $25\,\upmu \hbox {m}$ of our discretized boundary curve in Eq. (18). Finally, if more than 300 of our selected X locations fall in the domain, we use only the first 300 such locations in our initial condition.

For the special case of Fig. 5c, after specifying our initial condition as above, we scale the cell positions to account for one day of uniform epithelial growth (for details on how we implement domain growth, see Sect. 3.1). We use these scaled cell positions as our initial condition for the simulations in Fig. 5c.

1.4 B.4 Selecting Cell Birth Locations

We consider two methods for selecting $N^\text {M}_\text {diff}(t)$ possible locations for M birth:

• Control case (similar to body-model birth) We choose $2\times N^\text {M}_\text {diff}(t)$ locations uniformly at random in a rectangular region surrounding the fin, and we evaluate the first $N^\text {M}_\text {diff}(t)$ of these positions that are inside the fin domain for possible birth simultaneously.

• Mechanism V Under Mechanism V, we first choose $N^\text {M}_\text {diff}(t)$ points uniformly at random from our discretized bone rays, namely $\{\mathbf{B}_i^j\}$ in Eq. (1). For each such point $\mathbf{B} _k$, we then choose a corresponding location to evaluate for M birth by selecting a point uniformly at random in a ball of radius $r \sim \mathcal {N}(0, 2)\,\upmu \hbox {m}$ around $\mathbf{B} _k$.

We always select potential X birth locations in the same way as in the M control case.

Lastly, prior to applying our rules for cell birth to the locations that we randomly selected as we outlined above, we require that these positions are strictly greater than $25\,\upmu \hbox {m}$ away from our discretized boundary curve [see Eq. (17) in “Appendix B.2”]. If a randomly selected location is within $25\,\upmu \hbox {m}$ of a discretized boundary point in $\{\mathbf{F}_i(t)\}_{i=1,\ldots ,500}$, we do not allow cell birth to occur at that location.

1.5 B.5 Instructions for Reproducing Our Figures

We summarize the parameters for our simulated patterns in Figs. 5, 6, and 7 below:

• Figure 5a, Mechanism I (distal epithelial growth)

  • Cell migration: We use the body-model values in Table 2.

  • Cell birth: We use the body-model values in Table 3 with $n^\text {M}_\text {diff}= n^\text {X}_\text {diff} = 600$ locations.

  • Cell death: We use the body-model values in Table 4.

  • Skin growth: We do not scale cell positions with domain growth.

  • Switch parameters: In Eq. (12), we use $\zeta = -1$ cells. In Eqs. (13)–(14), we use $d = -1\,\upmu \hbox {m}$.

• Figure 5b, Mechanisms I and III (distal epithelial growth with alignment cues from the body)

  • Cell migration: We use the body-model values in Table 2 with one exception: $\varDelta t_\text {mig,birth} = 0.25$ days (more frequent cell birth and migration than we use in Fig. 5a).

  • Cell birth: We use the body-model values in Table 3 with three exceptions: $\varDelta t_\text {mig,birth} = 0.25$ days, $p_\text {M} = 0$, and $p_\text {X} = 0$ (no random birth). We use $n^\text {M}_\text {diff}= n^\text {X}_\text {diff} = 600$ locations.

  • Cell death: We use the body-model values in Table 4.

  • Skin growth: We do not scale cell positions with domain growth.

  • Switch parameters: In Eq. (12), we use $\zeta = -1$ cells. In Eqs. (13)–(14), we use $d = 150\,\upmu \hbox {m}$ (special body–fin interface dynamics; see Fig. 4g).

• Figure 5c, Mechanisms II and III (uniform epithelial growth with alignment cues from the body)

  • Cell migration: We use the body-model values in Table 2.

  • Cell birth: We use the body-model values in Table 3 with $n^\text {M}_\text {diff}= 1200$ and $n^\text {X}_\text {diff} =600$ locations.

  • Cell death: We use the body-model values in Table 4.

  • Skin growth: We scale cell positions along the bone rays with fin growth as we describe in Sect. 3.1.

  • Switch parameters: In Eq. (12), we use $\zeta = -1$ cells. In Eqs. (13)–(14), we use $d = 150\,\upmu \hbox {m}$.

• Figure 6, Mechanisms I, III, and IV (M migration along bone rays with alignment cues from the body and distal epithelial growth)

  • Cell migration: We use the fin-specific values in Table 2.

  • Cell birth: We use the fin-specific values in Table 3.

  • Cell death: We use the fin-specific values in Table 4.

  • Skin growth: We do not scale cell positions with domain growth.

  • Switch parameters: In Eq. (12), we use $\zeta =1000$ cells (so that M movement is always projected along the bones). In Eqs. (13)–(14), we use $d = 150\,\upmu \hbox {m}$.

  • Additional note in Fig. 6b: When we calculate the nearest-neighbor distances between cells, we only consider M–M and X–X distances that are less than $200\,\upmu \hbox {m}$ (this ensures that cells that have escaped our fin domain or appear at low density do not affect our measurements of cell–cell distances in developing stripes). When we calculate M–X distances at stripe–interstripe boundaries, we only consider measurements that are less than $110\,\upmu \hbox {m}$. We made this choice so that our M–X distances measure stripe–interstripe separation (in comparison, Takahashi and Kondo (2008) showed that M and X cells are roughly $82\,\upmu \hbox {m}$ apart at stripe–interstripe boundaries).

• Figure 7, Mechanisms I, III, and V (M birth in association with bone rays, with alignment cues from the body and distal epithelial growth)

  • Cell migration: We use the fin-specific values in Table 2.

  • Cell birth: We use the fin-specific values in Table 3. Additionally, we select locations for cell birth using the coordinates of the discretized bone rays in Eq. (1) (see “Appendix B.4” for details).

  • Cell death: We use the fin-specific values in Table 4.

  • Skin growth: We do not scale cell positions with domain growth.

  • Switch parameters: In Eq. (12), we use $\zeta = -1$ cells. In Eqs. (13)–(14), we use $d = 150\,\upmu \hbox {m}$.

  • Additional note in Fig. 7b: To calculate the distances that cells move per day, we consider the differences in their locations between consecutive days. In particular, the distance the ith M cell moves in one day is $||\mathbf{M} _i(t) - \mathbf{M} _i (t + \varDelta t_\text {mig,birth})|| + ||\mathbf{M} _i(t+\varDelta t_\text {mig,birth}) - \mathbf{M} _i (t + 2\varDelta t_\text {mig,birth})|| +||\mathbf{M} _i(t +2\varDelta t_\text {mig,birth}) - \mathbf{M} _i (t + 1)||$, since $\varDelta t_\text {mig,birth} = 1/3$ days in this simulation. If a new cell is born at position $\mathbf{M} _j$ at, for example, time $t+\varDelta t_\text {mig,birth}$, then we define the distance that cell agent moved between day t and day $t+1$ as just $||\mathbf{M} _j(t+\varDelta t_\text {mig,birth}) - \mathbf{M} _j (t + 2\varDelta t_\text {mig,birth})|| +||\mathbf{M} _j(t +2\varDelta t_\text {mig,birth}) - \mathbf{M} _j (t + 1)||$.