Simple Reaction-Diffusion Modeling Predicts Inconspicuous Neighborhood-Dependent Color Subclustering of Lizard Scales

Turing proposed a reaction-diffusion (RD) process as the chemical basis of morphogenesis. Despite the elegance of this model, its relevance for the precise description of morphogenesis in real organisms is largely disputed. Here, we show that a simple RD system, predicting the cellular-automaton-like patterning of ocellated lizards’ skin into green and black labyrinthine chains of scales, additionally predicts unsuspected subtle color subclustering that correlates with the colors of the scales’ neighbors. Hyperspectral imaging indicates that color subclustering is present in real lizards, confirming the numerical model nontrivial prediction. In addition, extensive histological analyses show that melanophores’ spatial distribution correlates with scale neighborhood, confirming that color subclustering is associated to the underlying microscopic system of chromatophore interactions. We then show that the observed subclustering is efficiently captured by RD models, irrespective of their form, discretization, and spatial dimensionality. We also show that sets of values can be identified in the 12-dimensional RD parameter space to yield the correct direction of correlation (i.e., observed in real lizards) between green-scale blackness and their neighborhood configuration, hence instructing the mathematical model. More generally, our results show that subtle mesoscopic properties of biological dynamical systems, as well as some of the underlying microscopic features, are quantitatively captured by simple RD models without integrating the unmanageable profusion of variables at lower scales.

Popular Summary

Turing posited that morphological patterns in living organisms—a zebra’s stripes, a human’s digits—could emerge from the interplay between chemical reactions and particle diffusion. This reaction-diffusion mathematical model is still viewed today with suspicion by most developmental biologists. They regard it as an artificial mathematical toy because it does not incorporate the many cellular and molecular variables that underlie chemical-based patterning. Here, we show that reaction-diffusion modeling is highly effective in quantitatively predicting the skin color patterning process.

When superposed on the geometry of reptilian scaled skin, a simple system of reaction-diffusion equations not only describes the cellular automaton—like scale-by-scale skin color patterning observed in some lizards, but also predicts unsuspected and subtle color subclustering that correlates with the colors of the scales’ neighbors. Hyperspectral imaging and histological analyses indicate that this color subclustering, observed in numerical simulations, is present in real lizards and correlates with the behavior of colored cells in the skin. We also show that the quantitative correlation observed in real lizards can be leveraged to produce an improved mathematical model.

These results are important in advancing the field of biological pattern formation by raising the bar on how to move beyond simple qualitative comparisons between theory and experiments, and by showing that reaction-diffusion models can predict unobserved features of these systems. More generally, our study shows that biology, despite its “messy” nature, with its unmanageable profusion of cellular and molecular variables, can be efficiently and quantitatively investigated mathematically, including with simple phenomenological models.

Figure 1

Discrete reaction-diffusion (dRD) simulations indicate neighborhood-dependent scale color subclustering. (a) Signs of interactions between chromatophores [17]; short-range compounds $U$ and $V$ represent melanophores (M) and xanthophores (X) in the dRD model; $W$ is the long-range factor. (b),(c) Enlarged detail of a lattice of 3850 hexagons with (b) juvenilelike initial condition and (c) steady-state pattern [$t>1800$; levels of green indicate unnormalized (left) and normalized (right) values of $V$]. (d) Heat maps: evolution of densities $U$ (top), $V$ (middle), and $W$ (bottom) for all hexagons. Subclustering of $(U,V,W)$ is apparent at late times in both the green and black clusters. Dynamical “lineages” of time evolution are visible at earlier times [insets (i) and (ii)]. Time trajectory of each hexagon is dynamically colored according to its nearest-neighbor configuration (bottom panel, $V_{\max }=10$). (e) Clustering and subclustering of $(U,V,W)$ of all hexagons at steady state ($t>1800$); values subcluster according to the number of black ($n_B$) or green ($n_G$) nearest neighbors. Coloring of data points as in (d). Top: green cluster (gc, $U<2.5$, $V>4.5$) and black cluster (bc, $U\ge 2.5$, $V\le 4.5$); the 3D graph is provided as Supplemental Material, File S1 (34). Middle: enlargement of the black cluster in $(U,V)$ projection; inset illustrates subsubclustering (for $n_B=2$). Bottom: enlargement of the green cluster in $(U,V)$ projection. White crosses indicate analytical approximation of steady-state subcluster locations (see Appendix pp1-s2).

Figure 2

Geometry perturbations increase dispersion of subclusters. Steady-state subclustering of the black and green clusters in four different geometries. (a) Weakly perturbed hexagonal lattice, (b) mildly perturbed hexagonal lattice, (c) strongly perturbed hexagonal lattice, and (d) Voronoi diagram of skin-scale positions on an actual ocellated lizard. From left to right, we show (i) a detail of the lattice with steady-state colors, (ii) the distribution of edge lengths $L$ (inset: distribution of the number of edges $Nn$ of the polygons), (iii) the black cluster in $(U,V)$ projection, (iv) the green cluster in $(U,V)$ projection, and (v) for the $n_G=2$ green subcluster: the correlation between the steady-state value $V_k$ and the total length ${\mathcal{L}}_k$ of edges in Eq. (1) that each polygon shares with neighbors of the same color. This correlation between $(U,V,W)$ and ${\mathcal{L}}_k$ explains overlapping between subclusters in strongly irregular geometries (c),(d). Note that the 3D distribution graphs in $(U,V,W)$ space of all scales are provided as Supplemental Material, Files S2–S5 [34] for (a)–(d).

Figure 3

Neighborhood-dependent scale color subclustering in real lizards. (a) Dorsal skin of an ocellated lizard (individual TL66) with a sample of 224 skin scales marked in red. Scales “a” and “b” are examples of neighborhood configurations $n_G=4$ (“a”) and $n_G=1$ (“b”). (b) Top: spectral modes associated to the first (solid line) and second (dashed line) principal components (PC) accounting for 90.6% of the dataset total variance. Bottom: deviation from the mean ($r_{\lambda ,i}-{\overline{r}}_{\lambda }$) of relative reflectance spectra (red solid line) and its projection on the first two PCs (blue dashed line) for scales “a” and “b”. (c) Top: coordinates of first two PCs for relative reflectance data, colored according to nearest-neighbor configurations $n_G=1$, 2, 3, 4; scales “a” and “b” are indicated. Bottom: Bayesian credibility ellipses for clusters’ averages of the same data ($p=0.9$ for inner ellipses, $p=0.99$ for outer ellipses). The first PC strongly correlates with neighborhood configuration. (d) Left: box plot and analysis of variance ($F$-test) for first PC data; clusters’ averages are significantly different. Right: $p$ values for Tukey-Kramer honestly significant difference (Tukey HSD) tests of pairwise comparisons between clusters’ averages. (e) Top: coordinates of first two PCs [dimensional reduction of $(U,V,W)$ space] for simulated data [dRD on real lizard Voronoi diagram of Fig. 2]; data points with different neighborhoods are sampled with the same frequencies as in real data (c). Bottom: Bayesian credibility ellipses represented as in (c). The origin of the inverse correlation with nearest-neighbor configurations of simulated (e) vs observed (c) data is analyzed in Sec. 4.

Figure 4

Histological analysis of chromatophore spatial distribution in lizard skin scales. (a) Example of histological section (HE staining) through scales “a” and “b”. (b) Segmentation (black, melanophores; blue, iridophores; red, epidermis) and alignment of sections for reconstruction of skin-scale chromatophore distribution. Left: a detail of raw image of section and its color segmentation. Top right: alignment of segmented sections for 3D reconstruction of the skin scale; a large proportion of each section in front of no. 44 is made transparent. The region of interest (“inside ellipse”) used for extraction of chromatophore distribution is indicated in light red. Bottom right: top view of the same alignment of sections, highlighting the construction of the light red ellipse used for data extraction. Stars indicate missing sections. (c) Spatial distribution of chromatophores (distance from epidermis of each labeled voxel) extracted from light red region. The 10th, 50th, and 90th percentiles are indicated.

Figure 5

Melanophore distribution and comparison with reflectance results. (a) Top: median distance to epidermis and volume fraction of voxels labeled as “melanophore” for each of the 41 skin scales (individual TL66) selected for histology. Bottom: Bayesian credibility ellipses for cluster averages of the same data, grouped according to nearest-neighbor configuration ($p=0.9$ for inner ellipses and dashed ellipses, $p=0.99$ for outer ellipses). (b) Box plot, analysis of variance ($F$-test, cluster averages are significantly different) and $p$ values for Tukey-Kramer pairwise comparison tests between cluster averages for melanophore distance data. See Supplemental Material, Fig. S10 [34], for melanophore volume fraction analysis. (c) PCA of relative reflectance data [as in Fig. 3] for this subset of 41 skin scales and (d) corresponding box plot, analysis of variance ($F$-test) and Tukey-Kramer pairwise comparison for first PC of relative reflectance data. The significance level of clustering is similar to that of melanophore distance in (b) and melanophore volume fraction in Supplemental Material Fig. S10 [34]. (e) Correlation between reflectance data (first PC) and melanophore volume fraction (left) or median distance to epidermis (right). Scales “a” and “b” are indicated in (a), (c), and (e).

Figure 6

Scanning of the dRD system for correct correlation between green-scale blackness and neighborhood configuration. (a) Left: green subcluster correlation direction $(\delta U_{\mathrm{rel}}^{(G)},\delta V_{\mathrm{rel}}^{(G)})$ in $U$ and $V$ directions for sets of parameter values associated to natural nearest-neighbor statistics; right: enlargement close to origin. Red and black data points are the 160 and 1194 sets giving correct and incorrect correlation directions, respectively. Note that 58 sets have undetermined correlation. No set simultaneously gives both $\delta U_{\mathrm{rel}}^{(G)}>0$ and $\delta V_{\mathrm{rel}}^{(G)}<0$. (b) An example of good parameter set ($c_1=-0.22810$, $c_2=-0.02610$, $c_3=0.63861$, $c_4=-0.17365$, $c_5=0.21950$, $c_6=0.83467$, $c_7=0.04327$, $c_8=-0.08629$, $c_9=0.12975$, $c_u=0.02467$, $c_v=0.20014$, $c_w=0.08955$). Left: steady-state pattern. Middle: corresponding nearest-neighbor statistics (solid lines) compared to those observed in an adult lizard (dotted lines) for both green (top) and black (bottom) scales. Right: $U$ and $V$ values of black (top) and green (bottom) subclusters, colored according to the number of isochromatic nearest neighbors. The green cluster (gc) has correct blackness ($U$ direction) correlation. (C) Parameter sets with natural nearest-neighbor statistics, colored according to correct (red) or incorrect (black) blackness correlation direction. Production rates $c_3$ and $c_6$ (left), decay rate $c_v$ and interaction $c_1$ (right) influence blackness correlation. (d) Lower panel: components of the Fisher linear discriminant direction ${\mathbf{v}}_{\mathrm{FLD}}$ best separating the sets with correct and incorrect blackness correlation directions; significant (bootstrapping, 10 000 random resamplings) nonzero coordinates are annotated with stars (*, $p\in [0.05,0.01]$; **, $p\in [0.01,0.001]$; ***, $p<0.001$). $F$-test value for the discrimination of correct vs incorrect correlations in the FLD direction is indicated. Upper panels: projections along ${\mathit{v}}_{\mathrm{FLD}}$ for parameters $c_v$ and $c_5$ that, respectively, most and least influence blackness correlation [correct (incorrect) directions are labeled in red (black)]. (e) ${\mathit{v}}_{\mathrm{FLD}}$ in (d) identifies the interaction (positive and negative arrows), production (baby icon), and decay (crossbones) terms whose absolute value should be increased (green) or decreased (red) to obtain correct correlation direction; i.e., all colored parameters should see their (nonabsolute) value increased, as visible in the bottom panel of (d). The blue square in (a), (c), and (d) indicates the set used in Figs. 1 and 2 giving incorrect blackness correlation of the green subclusters.