For the linear model, d u was increased from 0 to 0.2 within limits of the reaction term along the x -axis (Figures 4A, D–H). It can change the distance between the equilibrium point of u and the upper limitation [23] without shifting the maximum of the limit, and the parameters p v w and c w v were decreased linearly from 1 to 0 (Figures 4B,D–H) for investigation of the effect of gates on each cell. For short-range effects, i u v and i v u were decreased linearly from 1 to 0.6 in Figure 4C. Accordingly, the arbitrary parameter set generating stripes (wild type) was placed in the right top of the phase plane (Figures 4D–I). Partial differential equations (PDEs) were calculated with 20,000 and 40,000 iterations with d t = 0.1 and d t = 0.05 in fields sized x l = 56.25 , y l = 225 with d x = 0.75 , and x l =   y l = 200 with d s = 0.5 for Figures 4A–C and Figures 4D–H, respectively. These conditions satisfy each CLF condition. In Figure 4I, to investigate the simultaneous gap-junction effects with hemi-channels s d ( x , y ) = 1 − 0.002 max { x , y } , u d ( x ) = 0.0004 x , and v d ( y ) = 1 − 0.002 y were utilized in the field sized x = y = 300 ( d x = 0.75 ) and d t = 0.1 for 20,000 iterations.

For the nonlinear improved model with nonlinear terms, the parameters d u , i u v , and i v u were decreased linearly from 1 to 0.6. c wv and p vw were decreased linearly from 0.5 to 0.1 in Figures 5A–C and decreased by s d ( x , y ) = 1 − 0.002 max { x , y } simultaneously in Figure 5D. Accordingly, the arbitrary parameter set generating stripes (wild type) was placed in the right top of the phase plane. PDEs were calculated in fields sized x l = 10 , y l = 40 in Figures 5A–C and x l = y l = 50 with d s = 0.25 in Figure 5D. Then, after 500,000 iterations calculated with d t = 0.01 , we obtained the result.

Calculations were performed in the language Full BASIC ver. 8.1 with no-diffusion boundary conditions with difference calculus; then, results are shown as density plots of u . Parameters utilized in this study are summarized in Table 1.

Color pattern complexity and overall tone were quantified from binarized images using ImageJ as described in [24]. Briefly, the pattern simplicity score (PSS) is defined as the area weighted mean isoperimetric quotient of the contours extracted from each image. The overall color tone (OCT) of a pattern is defined and calculated as the ratio of white pixels in the binarized image. Analyzed images were prepared by the quaternary connection of a numerical result ( 100 ×   100 individual fields with periodic boundary conditions) of u in each parameter.

Previous laser ablation experiments reveal that the presence of mutual interactions between two types of pigment cells are necessary to generate Turing patterns (Figure 1E). Briefly, the density of melanophore existing and newly generated in a stripe decreases when the xanthophores in adjacent interstripes are ablated. On the other hand, that of xanthophore was not drastically changed. Then, two types of pigment cell inhibit each other at a one-cell distance even though the inhibition from xanthophore is inapparent until melanophores in adjacent stripes are eliminated. A mathematical model is derived from these apparent interactions in [12] although the details of the relationship between the experimental results and the model are not described. This model is based on the following set of RD equations: { ∂ u ∂ t = D u ∇ 2 u + f ( v , w ) − d u u ∂ v ∂ t = D v ∇ 2 v + g ( u , w ) − d v v ∂ w ∂ t = D w ∇ 2 w + h ( u , v ) − d w w . (1)

Here the alternative distribution of two types of pigment cells (Figure 1C) is expressed by two factors ( U , V ) of the three components. Then, u and v are each volume (viability). The numerical simulation of this model results in a Turing pattern in which u and v are distributed with antiphase, and a concentration of third factor W ( w ) presents peaks synchronized with u (Figure 2). It should be noted that cell divisions of differentiated melanophores contribute only minimally to the pigment-pattern formation in fish (Figure 1F). Therefore, the number of melanophores is changed by 1) the supply of new cells from randomly scattered precursor cells, 2) the death of existing pigment cells, or 3) the migration from a position close to the skin surface [25]. In the case of melanophores, it is known that cell movements and cell deaths are complementary to each other [26]. Even though they are inhibited, xanthophores are found in the stripe region, where they exist with a pale color [27–29]. Xanthophores do not move actively in vivo as may be the case for iridophores. As detailed in Figure 1F, motilities of the cells are approximated by small diffusion coefficients ( D ) . The rapidly diffusing factor W is assumed to have a large diffusion coefficient 1 in the outer region of the cells based on the results of electro-physiological experiments [30–32].

In the reaction, rather than it should be called an interaction, formulae, the dimension-less parameters are chosen arbitrarily from the sets that bring diffusion-driven instability. They are positive constants as shown in Table 1. Each formula is composed of a set of linear terms with upper and lower limits as utilized in the two-component system mentioned by [1] as follows: { f ( v , w ) = ( 0 ; − i u v v − i u w w + s u < 0 − i u v v − i u w w + s u ; 0 < − i u v v − i u w w + s u < f max f max ; − i u v v − i u w w + s u > f max ) g ( u , w ) = ( 0 ; − i v u u + p v w w + s v < 0 − i v u u + p v w w + s v ; 0 < − i v u u + p v w w + s v < g max g max ; − i v u u + p v w w + s v > g max ) h ( u , v ) = ( 0 ; p w u u − c w v v + s w < 0 p w u u − c w v v + s w ; 0 < p w u u − c w v v + s w < h max ; p w u u − c w v v + s w > h max h max ) . (2)

These equations simply indicate negative or positive interactions among two cells and molecules by the coefficients with different signs. They are derived from the interaction network obtained by the experimental results in [12]. The cells are assumed to be inhibiting mutually ( − i u v , − i v u ) , and then W is assumed to be produced by U ( p w u ) and consumed by V ( − c w v ) . Accordingly, W is assumed to inhibit pigment cells of the producer ( − i u w ) and produce (or preserve) the consumer ( p v w ) . The self-coupling parameters − d u , − d v , − d w corresponded to degradation (or death) coefficients, whereas s ’s represent constants related to the supply (also called “support sustainability”) of each component. The producer U activates V but then inhibits itself at long range via W . By consuming W , V also indirectly inhibits itself but then is activating U by double inhibition at long range. As a result, U and V exhibit no difference in apparent interactions, making it difficult to identify which factor corresponds to which cell type. Besides melanophores, xanthophores also show self-inhibition at long range. In laser ablation experiments, melanophore elimination in adjacent stripes causes pale-colored xanthophores in an interstripe region. Therefore, the pale color reflects xanthophore inhibition even without a change in cell number.

The three-component model in (1) is somewhat complex though it can be roughly regarded as a combination of two-component systems originally suggested by [33] as follows. For ease of mathematical analyses, I use the following two-component system that shares a component with high-diffusivity. { ∂ x ∂ t = f x x + f y y + D x ∇ 2 x ∂ y ∂ t = g x x + g y y + D y ∇ 2 y . (3)

There are two different cases that bring diffusion-driven instabilities, i.e., activator–inhibitor type (Figure 3A) or activator–substrate–depletion type (Figure 3B), characterized by different signs of the parameters [34, 35]. Both conditions are included in the three-component model sharing with the high-diffusible component (Figure 3C). Therefore, each parameter in the model can be regarded as part of a two-component system. The first and second reaction Eq. 2 include mutual inhibitions ( − i u v v , − i v u u ) , each of which corresponds to the respective self-activation ( f x > 0 for x ) though it cannot be realized without each partner. Recent experiments reveal that xanthophores are generated from division of other xanthophores [29]. Therefore, at least one self-reaction parameter; − d for u or v in Eq. 2 cannot be a degradation coefficient; i.e., it might be a nonnegative parameter.

Variation in patterns observed in most zebrafish mutants is explained by changes in the types and sizes of patterns. The former is defined by pattern selection [23, 36] and manifests as a general variation of the Turing pattern from spots to stripes to reverse spots. The latter is dictated by the characteristic wavelength of the pattern [37]. The following sections analytically describe the effect of parameters in the model (1)–(2) on these characteristics.

Next, correspondence between this model and molecular functions is assumed (Figure 1G). In zebrafish, leopard mutants are known to have an aberrant pigment pattern, whereby stripes are changed to spots [6] (Figures 1A,B). The gene responsible for the leopard mutation is a connexin cx41.8, which encodes a four trans-membrane connexin protein. Additionally, mutation of connexin cx39.4 results in wavy stripes. Recent analyses of connexin activity reveal different functions associated with hemi-channels and gap junctions [30, 31]. Hemi-channels are open to the extracellular environment, whereas gap junctions form connections between cells to allow the exchange of small molecules (Figure 1D). Accordingly, connexins may be involved in both long- and short-range interactions. These channels may function as gates for the transport of molecule W across the membrane (Figure 1G). Accordingly, producer U produces W , which then diffuses outside the cell into the extracellular space via some kind of gate. Therefore, gate defects affect survival of the producer by preventing the release of harmful W (i.e., d u of − d u is increased). I consider two other possibilities, i.e., the effect of the gate defect on D w and/or i u w . If D w is changed, W outside of the cells is also e affected. A u -dependent decrease in D w might be more appropriate for the assumption of enclosed W though it gives difficulties in the analytical approach and both finally result in an increase in the death of U . Then, the degree of harmful effect on U by the same concentration of w does not change. Therefore, D w and i u w are not affected. w has peaks with the peaks of producer U, and gates on cells assume passive effects on W movements. Defects of gate on a U cell will enclose W into the cell, then the production of W from U, i.e., p w u , other than d _u will also be affected. If p w u is decreased, v cells will be affected to decrease. However, the supply of w from other than u cells to the system will absorb the negative effects on v, even though the s _w in Table 1 is 0.0. Then the increase in d _u itself has the effect to decrease the concentration of W. Consumer V is assumed to incorporate beneficial W into the cell across the gates and then to consume it, indicating that the gate defects decrease the parameters c v w of − c v w and p w v . Both parameters are assumed to be related to intracellular events; therefore, higher w is needed for the same rate of consumption and V production compared with the case of intact hemi-channels. At the same time, parameters for mutual inhibitions i u v and i v u seem to be decreased by the leopard mutation [38] through gap junctions (combination of gates) at short range. They are simultaneously affected, linking with the hemi-channel defect on the corresponding cell. In Figure 1G, these parameters linked to different gates are indicated by different hatchings.

First, the independent effects of the gate on each cell were analyzed numerically. When an arbitrary stripe is set as a starting point, only gates that open to the outside on U and V cells were removed along the x- and y-axes. Numerical simulation of this model yielded a spot pattern of u in the case in which gates on U cells have defects as expected by the sign of the parameter (Figure 4A). Reversed u spots are yielded in V cell defects (Figure 4B) though the change is not strong because it includes opposite effects on V volume. PSS increases in both cases, and then OCT are decreased and increased by respective defects. That is, the asymmetry of changes in pattern selection can be observed by the removals of gates on respective cell types.

Defects on short-range inhibition do not have drastic effects on pattern selection though the pattern does finally disappear (Figure 4C). On the 2-D plane, the stripe region is recessive together with the defect in short-range effects by gap junctions (Figures 4D–H,J) though the tendency to shift the pattern selection is not changed. In Figure 4I, short-range effects are simultaneously decreased by respective x or y values that link with U or V defects as shown in the right panel in Figure 4J. This also keeps the same tendency to shift the pattern selection with individual cases.

The results mostly consist of the positive effects of connexin additions to WKO that increase the respective pigment cells [8]. However, experimental eliminations of the gate on each wild-type pigment cell lead to an increase in the rate of respective pigment cells; melanophore defects generate net (or rather wavy stripe) patterns of melanophores, and the xanthophore defects result in dot patterns of melanophores (i.e., net patterns of xanthophore). Therefore, the simulation results are partly inconsistent with experimental results in [8] in which the effects of connexins on each or both pigment cells are investigated in detail (Figure 4K).

From the electro-physiological experiments in [30]; each type of connexin can be considered to have different strengths of the (hemi-) channel functions on the two types of cells. Deduced patterns of respective transgenic fish are indicated by small letters on the phase plane in Figure 4I. Even though the differences in the strengths of channels are taken into account, the removal of hemi-channels on wild-type xanthophores tends to increase the proportion of xanthophores; then a faint increase of proportion in melanophore is brought in the case of connexin on melanophores (Figure 4K). If cell U is a melanophore, other than that the gray-framed patterns shown in Figure 4K do not seem to correspond to Figure 4I, all of the experimentally obtained patterns exist in the simulation.

From the analyses of wavelength mentioned above, defects to the gates on cell U (increasing in d u of − d u ) or cell V (decreasing in p v w and c w v ) cause a decrease or an increase in pattern size, respectively. In numerical simulations, the pattern size tends to be decreased and increased by the hemi-channel defects on both U and V cells as expected (Figures 3A,B). The characteristic thinning of u stripes and widening of v interstripes are observed in the simulation of U -cell defects though it is not explained by the analyses. From the characteristic thinning of the u stripe, each cell corresponding to one of two short-range factors may be deduced. However, the thinning is observed on melanophore stripe in the case of a defect in the xanthophore. Inconsistent with the U defect, the melanophore defect tends to result in wide melanophore stripes in the experiment.

To describe the detailed behavior of the components in the system, the model is changed to a model including nonlinear terms as shown in (5). { ∂ u ∂ t = D ∇ 2 u − i u v u v − i u w u w + d u u + s u 1 + u + v ∂ v ∂ t = D ∇ 2 v − i v u v u + p v w v w − d v v + s v 1 + u + v ∂ w ∂ t = ∇ 2 w + p w u u − c w v w v − d w w + s w . (5) The qualitative relationships for pattern formation shown in Figure 1F are not changed from model (1) and (2). That is, mutual inhibitions between U and V ( − i u v u v ,   − i v u v u ) are assumed; then, substance W with high diffusivity is assumed to inhibit the producer U ( − i u w u w ) and to promote consumer V ( p v w v w ) . Considering the interactions between different types of cells and between cells and molecules, the inhibitions were given by multiplication terms (e.g., the inhibition of U by W was described as u w and so on). These multiple terms are based on the description of the second order reaction in the chemical kinetic equation or dimer reactions by [39]; it enables only limited reaction by the contacts between the components. W is needed for maintenance of cell V , so this reaction is also given by a multiplication term of their volume. On the other hand, production of W by U only occurs u -dependently, and degradation (and death) is also w -dependent. Therefore, those terms are linearly related to each component. This model as an example of possible improvements also generated Turing patterns (Figure 5). These improvements can identify the functions on existing cells or newly differentiating cells. Then, new generations of pigment cells occur only with eliminations of existing cells [40]. Therefore, mature cells would inhibit the new generation of cells. The sign of d u seemed preferably to be positive for the starting point of pattern selection corresponding to zebrafish, i.e., the start from stripe. Even though the sign was opposite to the linear model, total u change may become minus with relation to other components (i.e., − i u v v − i u w w + d u < 0 in the deformed reaction terms ( ( − i u v v − i u w w + d u )   u + s u / ( 1 + u + v ) ) . Therefore, the self-productivity can be small enough to agree with the low proliferation rate of melanophores. Similarly, concerns about the self-productivity of xanthophores mentioned above, ∂ v / ∂ t in (Eq. 5) already have self-productivity by the multiplication term v w regardless of the sign of d v (i.e., − i v u u + p v w w − d v > 0 in deformed reaction terms ( − i v u u + p v w w − d v )   v + s v / ( 1 + u + v ) ). The sign of d v can be inverted though the change is not expected to substantially affect pattern characteristics.

Next, numerical calculations of the nonlinear model are performed to consider the condition in which respective gates on U and V cells and both-gate defects are added (Figures 5A–C). Again, numerical results consistent with the linear model can be obtained from an arbitrary parameter set generating a stripe pattern. A biased pattern shift can also be obtained by simulations, partly corresponding to connexin-mutation experiments. When gates on the U cell are deleted, it results in a u dot pattern (Figure 5A). Simultaneous decreases in p v w and c w v by increasing the gate defect on V cells generate a net of u though not so drastic (Figure 5B). A defect in the gap junction by decreasing i v u and i u v has a minimal effect though the stripe region became recessive with a combination of defects on each outer gate (Figures 5C,D). The thinning of the u population is not clear because of the thin stripe at the start.