Diffusiophoresis-Enhanced Turing Patterns

In his seminal paper, Alan Turing proposed that reaction-diffusion instabilities could give rise to pattern formation in biological systems, a phenomenon now known as Turing patterns [1]. Subsequent research has provided experimental evidence for the existence of Turing patterns in a variety of biological contexts, including zebrafish embryogenesis [2], hair follicle spacing [3], emulsion-based chemical cells [4], patterns found on zebrafish through cell-cell interactions [5], and the development of fingers from early limb-buds [6]. Theoretical advances have led to the development of modifications to Turing’s original model [7], including the widely-used Gierer-Meinhardt [8], Brusselator [9, 10] and cell-cell interaction [5] models. These models and their derivatives have been used to reproduce a wide range of patterns found in nature [11], from the wave-like patterns on marine angelfish [12], to the stripe-like patterns on zebrafish [5, 13], and even patterns on sea shells [14]. However, these models typically rely on diffusive transport mechanisms to generate concentration gradients and have thus far have neglected the role of convective transport. We propose that, in the context of Turing patterns, cells convect and steepen using diffusiophoresis; this mechanism promotes the color sharpening that is observed in nature.

Previous attempts to replicate the sharp gradients observed in natural Turing patterns have prescribed ad-hoc mechanisms including switching systems [7] and reaction pathway coupling [15]. However, a complete physical justification for such processes remains elusive. A more natural mechanism to obtain sharp gradients in concentration comes from operating in the regime of high Péclet numbers, i.e., conditions when the convective transport dominates the diffusive transport. This regime might appear to be difficult to obtain in the absence of fluid flow. However, such a regime is commonly observed in microfluidic experiments with colloidal particles via the process of diffusiophoresis, i.e., transport of colloidal particles in response to solute concentration gradients [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]. In fact, the key feature observed in diffusiophoretic systems is the banding of colloidal particles that creates region of sharp colloidal concentration gradients [16, 29, 19, 23, 24, 28, 27]. The diffusiophoretic transport of colloids, which has been seen to occur with a Péclet number $\textrm{Pe}=O(10)-O(10^{3})$ [16, 23, 28, 27] (the ratio of particle convective flux to diffusive flux), creates these banded structures that can be further enhanced in the presence of acid-base reactions [18]. We note that chemotaxis can be qualitatively similar to diffusiophoresis and has been proposed in the past [30] as an explanation for subtle features arising from time-dependent processes in the development of angelfish. Such an effect has, however, been neglected in more recent state-of-the-art explanations of chromatophore pattern formation [31, 5]. As we will argue, chromatophores are likely diffusiophoretic. Thus, we can recover key features by incorporating colloid physics into models of biological pattern formation models.

Over the past decade, several studies have concluded that diffusiophoresis is a key process that was previously overlooked. For instance, Florea et al. [32] elucidated, both theoretically and experimentally, that diffusiophoresis is responsible for creating an exclusion zone near charged surfaces and in fact stated that diffusiophoresis is “likely to play an important, yet unexplored role” in biological processes. Similarly, Shin et al. [33] underscored the role of diffusiophoresis in removing contaminants from porous, fibrous materials in the presence of surfactant concentration gradients. Therefore, we surmise that diffusiophoresis of chromatophores could play an important role in the process of biological pattern formation. To the best of our knowledge, this effect has not been investigated.

Some of the most readily observable Turing patterns are animal skin patterns. Our central argument is that during biological pattern formation, chromatophores, which are specialized pigment cells known to control the coloration pattern on fish [34], respond diffusiophoretically to physiological reactions. There are different types of chromatophores, including melanophores, xanthophores, erythrophores, iridophores and leucophores [35], which are organized in multiple layers and can interact with each other to create complex, even regenerative or adaptive patterns [36, 37]. We postulate that chromatophores are diffusiophoretic because of three primary reasons. First, the size of a typical chromatophore can range from 1-30 microns, which is typically the range where colloidal particles are diffusiophoretic [38, 39]. Second, recent experimental results suggest that chromatophores are charged and interact with the surrounding medium [40]. Finally, the biological reactions in the surrounding medium can create concentration gradients of small molecules that will interact with the chromatophore; multi-component physiological reactions have been known to control chromatophore aggregation and dispersion [41]. We also emphasize that diffusiophoresis of chromatophore-coated particles has been experimentally reported in literature [40]. In a biological setting, we expect that particle motion is made complex by effects such as cell-cell adhesion and resistance of the porous membrane. Regardless, such effects do not impact the essential physics leading to steepening.

The minimum mathematical model that qualitatively reproduces the features of biological patterns consists of two reactive-diffusive species that interact with each other through a nonlinear reaction mechanism [1, 12, 13, 8, 9, 10]. The most natural interpretation in the minimum mathematical model is to assume that these two species are chromatophores. We compare the natural patterns observed in Ornate Boxfish (Aracana Ornata) with the Brusselator model (see Methods for details) assuming a blue and a yellow chromatophore to be our reactive species. As observed in figure 1, the Brusselator model is able to reproduce the patterns of hexagon and stripes seen in the Boxfish, though it is unable to capture the sharpness of color gradients observed in the fish. Besides the size of the fish, in the natural fish patterns, there are two distinct length scales (see Methods). For instance, for the hexagon pattern, there is a clear length scale separation between the edge-thickness, $\ell_{\textrm{t}}$, and the length of the hexagon pattern, $\ell_{\textrm{p}}$. Similarly, for the stripe pattern, there is a distinction between the the thickness of the stripes $\ell_{\textrm{t}}$ and the separation between the stripes $\ell_{\textrm{p}}$. The reaction-diffusion models have an inherent limitation that $\ell_{\textrm{t}}=\ell_{\textrm{p}}$, leading to a diffuse pattern formation. Accordingly, they are unable to recover double spot patterns such as those seen on the Jewel Moray Eel (Muraena lentiginosa) unless multiple nonlinear reaction mechanisms are invoked, i.e., two or more Brusselators are coupled [15], which still yield a diffuse pattern unlike the ones observed in nature. Clearly, we must look for mechanisms beyond the nonlinear reaction-diffusion interactions.

While experiments [31, 5] have soundly implicated chromatophores in pattern-formation in zebrafish skin, recent evidence [5] suggests that a third species, which is more diffusive and is thus a small molecule, is required for a more complete picture of the process. This invites the possibility of robust control of separate length scales, where $\ell_{\textrm{p}}$ is associated with the molecular substance(s) and $\ell_{\textrm{t}}$ is associated with the chromatophores. Furthermore, prior studies on two-component models utilized species diffusion coefficients $O(10^{-10})-O(10^{-9})$ m${}^{2}/$s to obtain patterns that resemble the ones observed in nature [12]. These diffusivity values resemble small molecules instead of microparticles, which have diffusivities $O(10^{-13})$ m${}^{2}/$s. Accordingly, we assert that there are two broad categories of species present in the system: molecular-scale solute species (i.e. morphogens or long-range mediators) and microscopic-scale colloidal species (i.e. chromatophores). This is supported by prior experiments [41].

The transport equations for two-component chromatophore models are primarily dependent on the Damköhler number of the colloids $\textrm{Da}_{N}$, which is the ratio of the diffusion time scale to the reaction time scale, where an increase in $\textrm{Da}_{N}$ corresponds to a decrease in $\ell_{\text{p}}$. However, in our proposed model that builds on the recent experimental findings [5], the reactions are driven by solutes and the chromatophores respond to the concentration gradients of the solutes diffusiophoretically; see Methods for derivation. In this proposed setup, the two dimensionless groups that control the behavior of the patterns are the Damköhler number of the solute, $\textrm{Da}_{C}$, and the Péclet number of the colloids, Pe. Typically, $\textrm{Pe}=O(10)-O(10^{3})$ [16, 24, 26]. Physically, we propose that transport of the solute molecules sets only the pattern type and size $\ell_{\textrm{p}}$, whereas the diffusiophoretic transport of colloids dictates only the thickness $\ell_{\textrm{t}}$.

Our proposed model is able to recover to an considerable degree the natural patterns of hexagons and stripes observed in the Ornate Boxfish as well as the double-spot pattern in the Jewel Moray Eel using the Brusselator model with diffusiophoresis; see figure 1 and Supplementary Videos 1 and 2. While the reaction-diffusion model without diffusiophoresis produces qualitative features of the chromatophore patterns, by including diffusiophoresis we capture the steepening effect. This feature is universal and not exclusive to the Brusselator model. In figure 2 and Supplementary Videos 3 and 4, we perform the same comparison between reaction-diffusion models with and without diffusiophoresis for both the Gierer-Meinhardt model [8] and cell-cell interaction model [5] (see Methods: Numerical simulations for details). In these models, just as with the Brusselator in figure 1, the sharpening effect of diffusiophoresis is apparent both in the finer length scales and the greater magnitude of depletion away from the hotspots (i.e., the color of the two phases is more distinct with the inclusion of diffusiophoresis). Because of its relative abundance and depth of analytical theory in the literature, in this paper we center our analysis on the Brusselator model, but the fundamental physics at play is agnostic to the model.

To further emphasize the distinction between the proposed and the reaction-diffusion mechanisms, in figure S1, we compare the hexagon pattern for different values of $\textrm{Da}_{N}$ in the reaction-diffusion model and different values of Pe in the proposed model. As is clear from the results, increasing $\textrm{Da}_{N}$ simply decreases $\ell_{\textrm{p}}$ without changing the relative thickness of concentration gradients, i.e., $\frac{\ell_{\textrm{t}}}{\ell_{\textrm{p}}}$ is constant. In contrast, an increase in Pe does not modify $\ell_{\textrm{p}}$ but reduces $\ell_{\textrm{t}}$, and better resembles the patterns observed in nature. In our simulations, we found strongest qualitative agreement with natural patterns when setting Pe in the range $10-50$. This is consistent with the effective Pe number observed in recent theoretical and experimental findings [28, 27]. We acknowledge that we neglect inter-particle interactions and invoke the assumption of dilute suspension, which in some limits may cause our simulations to under-predict values of $\ell_{\textrm{t}}$.

Next, we focus on quantifying the pattern formation of chromatophores via diffusiophoresis. For simplicity, we focus on the Brusselator model for solute reactions (figure S2 left panel), though our analysis is readily extended to other reaction-diffusion frameworks [13]. In the Brusselator model [9, 10], two abundant reactants with concentrations $A$ and $B$ produce solutes with concentrations $C_{1}$ and $C_{2}$, where concentrations are scaled by $c^{*}$. An initially homogeneous state of $C_{1}$ and $C_{2}$ can spontaneously come into formation when perturbations of permitted wavelengths grow into patterns such as hexagons or stripes; the resulting gradients $\nabla C_{1}$ and $\nabla C_{2}$ generate the diffusiophoretic velocity $V_{\text{DP}_{j}}=M_{j1}\nabla C_{1}+M_{j2}\nabla C_{2}$ of the colloid (i.e. chromatophore) for diffusiophoretic mobility $M_{i}$ corresponding to $C_{i}$, where $M_{ji}$ is the dimensionless diffusiophoretic mobility that induces the velocity in the $j^{\textrm{th}}$ colloid due to the $i^{\textrm{th}}$ solute. The results in figure S2 highlight that two different diffusiophoretic colloids with different mobilities, denoted with white and red dots, mirror the pattern of the solute to form hexagons and stripes, but produce sharper and more well-defined features than the solute alone. Only concentration hotspots of $C_{1}$ are shown (in green) since the concentration distribution of $C_{2}$ has the same spatial structure as $C_{1}$. Here, we further reduce the parameter space by focusing on a system with one colloid only and therefore we drop the $j^{\textrm{th}}$ subscript in $M_{ji}$ and refer to it as $M_{i}$.

By employing the framework of amplitude equations, we analytically predict solute and relative colloid concentrations for different patterns produced with the Brusselator model, including an analytical expression for the steepening ratio $\frac{\lambda_{N}}{\lambda_{C}}$; see Methods for details. Theoretical investigations of the Brusselator model [42] have revealed a complex parameter space where different Turing patterns form. The perturbations are sinusoidal in nature and have amplitudes that obey complex partial differential equations that encode the stability and form of these structures. It is worth noting that these equations are not deterministic, and merely describe the shape and stability of all possible structures [11] such that in a simulation, any stable pattern can be forced, and in experiments, we expect to the fastest-growing wavelength to dominate.

One of the dominant patterns in the Brusselator model is hexagons. Simulating stable hexagon patterns, in figure 3(a) we plot the dependence of solute (left) and colloidal (right) concentrations with distance away from the center of the hexagon (averaged over all hexagons in the simulation box). The other dominant pattern in the Brusselator model is stripes, which we investigate in figure 3(b). The solute (left) and colloidal (right) concentrations here are plotted along the periodic direction.

There are two length scales in each of the patterns described in figure 3. In hexagons, it is the size of the hexagons and the edge-thickness of the hexagons. For stripes, it is the spacing between the stripes and the thickness of the stripes. The size of the hexagons and the spacing between the stripes is set by $\lambda_{C}$, i.e., the wavelength of the solute patterns, and the thickness of hexagons and stripes is set by $\lambda_{N}$. Crucially, the ratio of $\frac{\lambda_{N}}{\lambda_{C}}$ is set by the Péclet number of the colloid for both the patterns. To explore this effect in more detail, we mathematically derive (see Methods for derivation) the Péclet number as

$$\text{Pe}=\bigg{\rvert}\frac{\alpha\left(M_{1}-M_{2}\frac{\eta(1+A\eta)}{A}\right)}{\mathcal{D}_{N}}\bigg{\rvert},$$

(1)

where $\alpha$ is the amplitude of the perturbations to $C_{1}$ and $\eta=\mathcal{D}_{C_{2}}^{-1/2}$. We note that equation (1) is valid for both hexagon and stripe patterns. In the Brusselator model, the non-dimensional solute concentration local to hexagon or stripe structures as a function of distance $\mathscr{R}$ from the local maximum is given by

$$C_{i}(\mathscr{R})=\sigma_{i}\alpha\cos\left(2\pi\frac{\mathscr{R}}{\lambda_{C}}\right)+\text{const},$$

(2)

where we note that the one-dimensional form of equation (2) is exact in the direction of maximal gradient for stripes and approximate for radially-averaged hexagons. From the gradient of equation (2), we calculate the diffusiophoretic velocity

$$V_{\text{DP}}=-\frac{2\pi}{\lambda_{C}}\alpha\left(M_{1}-\frac{\eta(1+A\eta)}{A}M_{2}\right)\sin\left(2\pi\frac{\mathscr{R}}{\lambda_{C}}\right)$$

(3)

and the form of the steady-state colloid concentration

$$N(\mathscr{R})\propto\exp\left(\text{Pe}\cos\left(2\pi\frac{\mathscr{R}}{\lambda_{C}}\right)\right).$$

(4)

We define the colloid length-scale $\lambda_{N}$ as the exponential decay distance, giving the analytical expression

$$\cos\left(2\pi\frac{\lambda_{N}}{\lambda_{C}}\right)=1-\frac{1}{\text{Pe}};$$

(5)

see Methods: Analytical model for derivations. In figure 3, we show that our analytical model is in good quantitative agreement with simulation data extracted from a large number of both hexagons and stripes for solute and colloidal concentrations. We note that the analytical model has no discontinuities, and the colloid hotspots are predicted on a continuous and periodic curve.

It is clear from equation (5) that $\frac{\lambda_{N}}{\lambda_{C}}$ can vary significantly through changes in Pe, revealing a control mechanism for the sharpness of Turing patterns involving colloidal particles. We further elucidate this control mechanism in figure 4 by comparing numerical predictions of $\frac{\lambda_{N}}{\lambda_{C}}$ vs. Pe to the analytical theory for the Brusselator model (equation 5). We obtain an approximately linear relationship between the two variables, consistent with our theoretical predictions, and show a collapse for two orders of magnitude of Pe. In the limit of low Pe, the steepening effect is diminished and the colloid length scale approaches the solute length scale. If Pe is exactly zero, however, there would be no mechanism for colloidal formation and thus $\frac{\lambda_{N}}{\lambda_{C}}$ would be undefined.

In addition to the Brusselator model, we include in figure 4 simulation data points for the Gierer-Meinhardt and the cell-cell interaction models. In the absence of an analytical prediction for the form of the solute field perturbations in these models, we estimate $\textrm{Pe}=\frac{M}{\mathcal{D}_{N}}$, noting that we observe perturbation amplitudes of $\mathcal{O}(1)$ in our simulations. A more accurate representation of Péclet number would be the general form of equation (1), $\textrm{Pe}=\sum_{i}\alpha_{i}\frac{M_{i}}{\mathcal{D}_{N}}$, where the sum is across each solute species and the coefficients $\alpha_{i}$ are dependent upon the specific solute reaction-diffusion model. However, the analytical theory necessary to obtain such coefficients does not exist, to the best of our knowledge, for the Gierer-Meinhardt [8] and cell-cell [5] interaction models. Nonetheless, the data points from our simulations of these models are well-represented by a best-fit line, suggesting the collapse to a master curve of $\frac{\lambda_{N}}{\lambda_{C}}$ vs. Pe for models other than the Brusselator.

The findings reported in this article open up numerous possibilities for future research and applications. For instance, further exploration of biological systems could reveal more instances where diffusiophoresis contributes to pattern formation such as embryo morphogenesis, which could subsequently lead to a better understanding of how these patterns emerge, evolve, and adapt in different species. In recent years, an embryonic framework has gained attention for understanding cancer ecosystems as chaotic expressions of the complex interactions between cellular and chemical components [43], and Turing patterns are theorized to play an important role in tumorigenesis [44]. Unveiling diffusiophoresis as a mechanism in cellular organization processes can be crucial to future cancer research. Additionally, a more comprehensive understanding of the interaction between different species of chromatophores and solutes might be achieved by exploring various reaction rates, diffusivities, and mobilities. Going forward, research might consider complex colloidal interactions [45], dense packing effects [46], anisotropy [47], particle self-propulsion [48], and hydrodynamic interactions, which could reveal important control mechanisms of pattern formation in various biological systems. Experimental studies can observe the dynamics of the formation of Turing patterns with colloidal particles in controlled conditions. Beyond biophysics, understanding the role of diffusiophoresis in pattern formation has potential applications in engineering and materials science. By harnessing these mechanisms, researchers could develop novel methods for fabricating materials with precise micro-scale patterns. These materials could find applications in areas such as photonics, lab-on-a-chip devices, and biotechnology. In conclusion, our findings on the role of diffusiophoresis in pattern formation provide a foundation for future research and have the potential to impact a wide range of fields.