Pinning of Diffusional Patterns by Non-Uniform Curvature

TL;DR

This paper investigates how non-uniform curvature on surfaces influences the pinning and selection of diffusion-driven patterns, using conformal mapping and perturbation theory, with implications for biological and physical systems.

Contribution

It introduces a novel theoretical framework combining conformal mapping and perturbation theory to analyze pattern pinning by curvature inhomogeneities.

Findings

Curvature inhomogeneities can pin diffusion patterns.

The theory draws an analogy to quantum mechanics in a geometry-dependent potential.

Numerical confirmation supports the theoretical predictions.

Abstract

Diffusion-driven patterns appear on curved surfaces in many settings, initiated by unstable modes of an underlying Laplacian operator. On a flat surface or perfect sphere, the patterns are degenerate, reflecting translational/rotational symmetry. Deformations, e.g. by a bulge or indentation, break symmetry and can pin a pattern. We adapt methods of conformal mapping and perturbation theory to examine how curvature inhomogeneities select and pin patterns, and confirm the results numerically. The theory provides an analogy to quantum mechanics in a geometry-dependent potential and yields intuitive implications for cell membranes, tissues, thin films, and noise-induced quasipatterns.

Figures

Equations

$$\text{differential}{s}^2=(1+R_z^2)\text{differential}{z}^2+R^2\text{differential}{\phi }^2,$$

$$\Delta \varphi$$

$$\text{differential}{s}^2={\Omega }^2(v)\left(\text{differential}{v}^2+R_0^2\text{differential}{\phi }^2\right),$$

$$-{\Delta }^0{\varphi }_k+2\lambda ϵh(v){\varphi }_k=-\lambda {\varphi }_k.$$

$${\varphi }_{sk}^{(0)}$$

$${\lambda }_{sk}^{(0)}$$

$${\lambda }_{0k}^{\pm }$$

$$\displaystyle\partial_{t}\left(\begin{array}[]{c}\Psi_{1}({\bf x},t)\\ \Psi_{2}({\bf x},t)\end{array}\right)=\overbrace{\left(\begin{array}[]{c}R_{1}\left(\Psi_{1},\Psi_{2}\right)\\ R_{2}\left(\Psi_{1},\Psi_{2}\right)\end{array}\right)}^{\text{Reactions}}+\overbrace{\left(\begin{array}[]{cc}\nu_{1}&0\\ 0&\nu_{2}\end{array}\right)\Delta^{G}\left(\begin{array}[]{c}\Psi_{1}\\ \Psi_{2}\end{array}\right)}^{\text{Diffusion}}\,.$$

$$\displaystyle\left(\begin{array}[]{c}\Psi_{1}({\bf x},t)\\ \Psi_{2}({\bf x},t)\end{array}\right)=\left(\begin{array}[]{c}\Psi_{1}\\ \Psi_{2}\end{array}\right)^{*}\!+\displaystyle\sum_{k}\left(\begin{array}[]{c}u_{1k}\\ u_{2k}\end{array}\right)e^{t\omega(k)}\phi_{k}\left({\bf x}\right)\text{ ,}$$

$$\displaystyle\left[\left(\begin{array}[]{cc}R_{1,1}&R_{1,2}\\ R_{2,1}&R_{2,2}\end{array}\right)^{*}\!\!+\left(\begin{array}[]{cc}\nu_{1}\lambda_{k}&0\\ 0&\nu_{2}\lambda_{k}\end{array}\right)\right]\left(\begin{array}[]{c}u_{1k}\\ u_{2k}\end{array}\right)=\omega(k)\left(\begin{array}[]{c}u_{1k}\\ u_{2k}\end{array}\right)\text{.}$$

$${\dot{\Psi }}_1$$

$${\dot{\Psi }}_2$$

$$\left[-{\Delta }^0-2ϵ{\lambda }_{sk}^{(0)}\cos (\bar{p}v)\right]{\varphi }_k=k^2{\varphi }_k.$$

$$\text{differential}{s}^2=\left(R^2+R^{\prime 2}\right)\text{differential}{\theta }^2+R^2{\sin }^2\theta \text{differential}{\phi }^2,$$

$$\text{differential}{s}^2={\Omega }^2\left(\text{differential}{\Theta }^2+{\sin }^2\Theta \text{differential}{\phi }^2\right).$$

$$\left[-{\Delta }^0+k^2\left(R_0^2-R^2\frac{{\sin }^2\theta }{{\sin }^2\Theta }\right)\right]\Phi =k^2{R}_0^2\Phi ,$$

$$\left[-{\Delta }^0-2ϵk^2{R}_0^{2}h\right]\Phi =k^2\Phi .$$

$$⟨\text{absolutevalue}{\mathbf{\Psi }(x)}^2⟩\propto \sum _k\frac{{\cos }^2(\bar{k}x)}{\det \mathbf{M}({\lambda }_{0k}^{+})}+\frac{{\sin }^2(\bar{k}x)}{\text{det}\mathbf{M}({\lambda }_{0k}^{-})}+\text{order}{ϵ}^2.$$

Full text

In 1952, Turing coined the term “morphogen” in a seminal paper that showed how combining diffusion with generalized reactions can create spatial and temporal patterns even though separately each leads to uniform, static concentrations Turing (1952). Such reaction-diffusion (RD) patterns are but one example of diffusion-driven instabilities that have since been studied on scales ranging from microns in cells Hecht et al. (2010) and active fluids Bois et al. (2011); Giomi and Marchetti (2012); Pearce et al. (2018); Dürre et al. (2018); Marchetti et al. (2013); Kumar et al. (2014), millimeters in hydrodynamics Almarcha et al. (2010), centimeters in zoology Murray (1981), to meters in ecology Butler and Goldenfeld (2011). Diffusion-driven patterns are known to determine morphology Kondo and Miura (2010) in model organisms like zebrafish Nakamasu et al. (2009); Yamanaka and Kondo (2014) and complex organs like the eye Sasai (2013). Recent theoretical progress in patterning Biancalani et al. (2017); Corson et al. (2017) encourages further study.

Substrate curvature plays a role in pattern formation in many systems, including cell membranes Isaac et al. (2013) and thin films Turner et al. (2010). The importance of surface curvature on collective behavior has recently been explored in liquid crystals Fialho et al. (2017), flocking Shankar et al. (2017), and wave propagation Miller et al. (2018). Closer to our work, the geometric dependence of pattern formation has recently been studied in various models of protein Thalmeier et al. (2016); Vandin et al. (2016) and molecular bonding Yu et al. (2017).

Recent studies of Turing patterns have explored the effects of curvature on highly symmetric shapes such as spheres, cylinders, toroids Nampoothiri (2016); Sánchez-Garduño et al. (2018), and ellipsoids Nampoothiri and Medhi (2017) where the Laplacian is known in closed form. Inhomogeneities in curvature, such as protrusions or cavities, reduce such symmetries and can pin or modify the patterns. To understand how nonuniform curvature can entrain and modify patterns, we study perturbations to the Laplacian, and its eigenmodes. To our knowledge, the intimate link between pinning of patterns and the spectrum of the Laplacian has not been pursued.

We follow a two prong strategy. First, we identify the onset of instabilities, by linearizing evolution equations expressed in terms of the appropriate ‘modes’, e.g., Fourier, cylindrical, or spherical harmonics. Modes with the largest positive real part grow fastest and are harbingers of the final patterns molded by nonlinearities. The modes are eigenfunctions of the diffusion (Laplacian) operator on the relevant manifold. Symmetries of the manifold, reflected in degeneracies of the eigenfunctions, must be broken in the final patterns. Previous work on the Laplacian on Riemannian manifolds focused on its determinant Kac (1966); Schoen and Yau (1994); Rosenberg (1997); Braun et al. (2008) and short-time behavior appropriate to field theory David (1987); Faraudo (2002); Castro-Villarreal (2010). We focus instead on how non-uniform curvature breaks degeneracies, pinning eigenfunctions to inhomogeneities. To do this, we utilize conformal mappings and perturbation theory. There is no guarantee, however, that patterns resulting from non-linear evolution are similarly entrained, so the second step in our study explores patterns with simulations. We implemented the Thomas-Murray RD equations Murray (1981) on COMSOL Multiphysics® COMSOL . We conclude with suggestions for experiments.

We first consider a cylinder with axially symmetric deformations described in cylindrical polar coordinates as $\rho=R(z)$. The surface line element is

[TABLE]

where $R_{z}$ denotes the derivative with respect to $z$. The Laplacian of a scalar $\phi$ is

[TABLE]

where $g_{ab}$ is the metric of the underlying geometry; $g$ is the determinant of the metric. Conformal mapping simplifies analysis through mapping to a flat geometry. We introduce a conformal axial coordinate $v$, such that the line element acquires the conformally flat form

[TABLE]

where $R_{0}$ is the asymptotic radius and $\Omega$ is the conformal factor. In the conformal coordinates, $\sqrt{g}=\Omega^{2}R_{0}$, and thus the Laplacian on the deformed geometry takes the simple form $\Delta^{G}=\Omega^{-2}\Delta^{0}$, where $\Delta^{0}$ is the Laplacian in the conformally flat coordinates. Since the behavior of $\Delta^{0}$ is well understood, and $\Omega$ is determined by the equality of Eqs. 1 and 3, this conformal mapping provides a tractable method of understanding $\Delta^{G}$.

Solutions for $v$ and $\Omega$ for arbitrary $R(z)$ are in the supplement. To develop perturbation theory, we set $R=R_{0}(1+\epsilon h(z))$, in which case $v\approx z$ and $\Omega\approx 1+2\epsilon h(v)$ to lowest order, such that the eigenfunctions $\phi_{k}$ of $\Delta^{G}$, with eigenvalues $\lambda$, satisfy

[TABLE]

In analogy with quantum mechanics, one can interpret the deformation as giving rise to a potential in the conformal coordinates, whose magnitude is dependent on the eigenvalue $\lambda$. (This differs from da Costa’s geometric potential, which comes from confining a particle to a surface da Costa (1981).) For physical phenomena described by the Laplacian, this mapping can be interpreted as the replacement of the diffusion operator $-D\Delta$ on a deformed geometry, with a spatially-dependent diffusion coefficient $\tilde{D}(v)\equiv D/\Omega^{2}$ in the conformally-related homogeneous geometry. This provides an intuitive picture of how diffusion is modified by curvature.

The undistorted cylinder of length $L$ with periodic boundary conditions has eigenfunctions and eigenvalues

[TABLE]

with integers $s$ and $k$. Consider an axially-symmetric bump on the cylinder, shaped like a Gaussian of standard deviation $\sigma$ and height $v_{0}$. We apply Rayleigh-Schrödinger perturbation theory to Eq. 4 to calculate eigenvalue corrections to first order in $\epsilon\equiv v_{0}/R$, which simplify for the case $s=0$ to

[TABLE]

The positive/negative sign is set by the mode: positive for cosine (symmetric) modes, and negative for sine (antisymmetric) modes. The sign of $\epsilon$ depends on the orientation of the ridge-like deformation: positive for a bulge and negative for a constriction. Thus, to first order, a ridge breaks the degeneracy in the eigenvalues of a Laplacian on the cylinder, causing the eigenvalues of sine (cosine) modes to become more negative in the case of an outward (inward) ridge. This correction is largest for $2\pi/\bar{k}=4\pi\sigma$, when the wavelength is approximately twice the width of the ridge. This is a general trait of deformations, explored later for a rippled cylinder. The ridge can be positioned anywhere along the cylinder, and the modes will shift with the ridge.

Modifications to the Laplacian spectrum affect any physical system involving diffusion. For example, consider a two-component RD system:

[TABLE]

After linearization, small deviations around a uniform stable fixed point of the reactions, $(\Psi_{1},\Psi_{2})^{*}$, evolve as

[TABLE]

with $\omega(k)$ satisfying eigenvalue equation, $\mathbf{M}(\lambda)=\mathbf{R}+\lambda\boldsymbol{\nu}$,

[TABLE]

While the diffusion and stable reaction matrices separately posses negative eigenvalues, Turing showed that their sum can have positive eigenvalues ($\omega_{+}(k)>0$) signaling finite-wavelength instabilities Turing (1952). The possibly degenerate modes, $k^{*}$, with largest eigenvalue evolve to the final pattern. On a cylinder with a ridge, the degeneracy between sine and cosine modes is broken: an outward (inward) ridge leads to sine (cosine) growing faster.

This linear analysis indicates only the onset of instabilities. To understand the patterns formed after nonlinearities stabilize the dynamics, we conducted finite element simulations of the Thomas-Murray model Murray (1981):

[TABLE]

We used COMSOL Multiphysics®COMSOL , which approximates the Laplacian by finite differences on a mesh and computes the fully non-linear reaction terms. (Supplement includes .mph file.) It should be noted that, although periodic boundary conditions are enforced in the numerical calculations, the computational methods using by COMSOL have a tendency to pin the resulting Turing patterns in a specific configuartion, even on a flat cylinder. Despite this, as shown in Fig. 1, density maxima of Turing patterns on a deformed cylinder are entrained by a ridge. An inward (outward) ridge selects the sine (cosine) mode.

This conformal mapping approach enables studies of other geometries, including bumps on spheres (below), drums (supplement), and rippled cylinders. An axially symmetric rippled cylinder, $h(z)=\cos\left[(2\pi p/L)z\right]\equiv\cos(\bar{p}z)$, gives a Schrödinger-like equation at $\order{\epsilon}$ (cf. Eq. 4):

[TABLE]

This Schrödinger equation describes a particle moving in a weak periodic potential, whose properties are well understood in the context of solid-state physics Kittel (2004). At leading order, this perturbation gives rise to a broken degeneracy (band-gap) at $k=p/2$ with magnitude $4\epsilon\lambda_{sk}^{(0)}$. Our analysis predicts that when an RD system governs surface concentrations, the effective diffusion rate will increase in troughs and slow down on ridges. Hence, diffusion is enhanced (diminished) where Gaussian curvature is negative (positive). This agrees with the short-time analysis of diffusion on Riemannian manifolds, where the leading order correction to diffusion is proportional to Gaussian curvature David (1987); Faraudo (2002); Castro-Villarreal (2010). For Turing patterns, steady-state regions of high concentration switch sharply from ridges to troughs as the most unstable wavelength is dialed past twice the ripple wavelength, see Fig. 2. Note that Eq. LABEL:periodic predicts a broken degeneracy only when the most unstable wavelength is commensurate with the ripple wavelength. However, in numerical simulations (presumably due to higher-order effects), we observe pinning for a range of wavelengths close to commensurability, although the patterns become unpinned for sufficiently incommensurate wavelengths.

An axially-symmetric distorted sphere, described by $R(\theta)$, has line element

[TABLE]

where prime denotes derivatives. Mapping to conformally flat coordinates:

[TABLE]

The Laplacian eigenvalue equation in conformal coordinates becomes (see supplement)

[TABLE]

where $\Delta^{0}$ is the Laplacian on a round sphere. Setting $R=R_{0}\left(1+\epsilon h(\theta)\right)$, we expand in powers of $\epsilon$. To $\order{\epsilon}$, $\theta=\Theta$, and Eq. (32) reduces to

[TABLE]

Using Equation (33), we study how deformations modify the Laplacian eigenfunctions $Y_{\ell}^{m}(\theta,\phi)$. The perturbation depends only on $\theta$, breaking rotational symmetry by fixing polar orientation of eigenmodes while preserving azimuthal symmetry. This can potentially entrain Turing patterns, although there are competing influences from nonlinear effects and incommensurability of length scales as seen in Fig. 3. In fact, in numerical simulations, we often found that the initial patterning predicted by linear stability analysis would stabilize to more complex patterns. In contrast, on a cylinder, initial patterning tends to persist indefinitely. We attribute this to differing degrees of degeneracy. The spherical eigenmodes, $Y_{\ell}^{m}$, have $2\ell+1$ degeneracy. Thus, near the most unstable eigenmode, there are $2\ell$ additional unstable modes that can contribute to non-linear pattern formation. In contrast, a cylinder presents only a twofold degeneracy.

As evident from Eq. 33, eigenfunction modifications are sensitive to the sign of $h$. An inward (outward) bump causes an increase (decrease) in diffusion. Non-linear effects in the Thomas-Murray model are known to stabilize spotted patterns instead of the single most unstable spherical harmonic, see Ref. Barrio et al. (1999). However, during the initial development of Turing patterns, the most unstable $Y_{\ell}^{m}$ is visible, and the mode selected by a deformed sphere differs from that of an uniform sphere, see Fig. 3.

Including extrinsic noise in RD systems leads to temporally fluctuating quasi-patterns for a broader range of parameters than those required for a Turing instability  Lugo and McKane (2008); Butler and Goldenfeld (2011); Schumacher et al. (2013); Biancalani et al. (2017). We have investigated how curvature influences such noise-induced transient patterns. We add uncorrelated white noise $\eta(\mathbf{x},t)$ of zero mean and variance $\mathcal{D}$ to Eq. 13. In the linearly stable regime, this noise leads to eigenfunction fluctuations with a power spectrum $P(\lambda)$ proportional to $\mathcal{D}\left[\text{det}\mathbf{M}(\lambda)\right]^{-1}$. On a uniform surface, time averaged fluctuations are translationally invariant, $\langle\absolutevalue{\mathbf{\Psi}(x)}^{2}\rangle=\mathbf{\Psi}_{0}^{2}$. Deformations break this symmetry, leading to average fluctuations on a cylinder of the form

[TABLE]

On spheres, the spherical harmonics replace the Fourier modes. Figure 4 shows numerical verification that such deformations fix the phase of fluctuations and create non-uniform time-averaged intensities. Such behavior is expected to hold generically for quasi-patterns induced by intrinsic noise Lugo and McKane (2008); Butler and Goldenfeld (2011); Schumacher et al. (2013); Biancalani et al. (2017), where the power spectrum describing fluctuations is qualitatively similar.

If these effects can be realized in manufactured or experimental systems they could enable systematic manipulation of patterns. For example, collagen vitrigel (CV) for corneal endothelial regenerative treatment Yoshida et al. (2014) could be molded in a hemispherical or eggcrate geometry and used to support zebrafish chromatophores, which were recently cultured for in vitro studies Yamanaka and Kondo (2014) to explore RD models Nakamasu et al. (2009). Mice hair follicle patterning is also driven by an RD system involving WNT growth factor Sick et al. (2006), and might be cultured on molded CV using established in vitro techniques Philpott et al. (1990). Zebrafish stripes and hair follicles have spacing on the order of hundred microns. Three orders of magnitude smaller, cytoskeletal suspensions exhibit patterns Bois et al. (2011); Giomi and Marchetti (2012); Pearce et al. (2018) that may realize our results. For example, a pinched sphere might be constructed in the recently studied encapsulations of actomyosin in giant unilamellar vesicles Dürre et al. (2018). A recent review of such active fluids suggests that coupling to RD systems looks like a fruitful direction for research Marchetti et al. (2013).

Even in a system as simple as a cylinder with a ridge, geometry can dramatically affect pattern formation. Our approach to the Laplacian on curved surfaces opens a new route to analytical understanding of patterns in real systems, taking advantage of intuition and tools from quantum mechanics. Future directions include pattern formation in the presence of advection Kumar et al. (2014), and on time-varying shapes. Our approach can also be applied to scenarios involving interactions between the surface and boundary physics. In particular, it would be interesting to examine the effect of nonuniform surface curvature when additional reaction-diffusion processes take place at the boundary Halatek and Frey (2012). Our mathematical methods apply to any process described by the Laplacian, and may also find application in soap films Osserman (2002), or Marangoni flows from surface tension gradients Pereira et al. (2007).

JRF acknowledges support from the Hertz Foundation. JG acknowledges support from CONACyT grant 180901. MK acknowledges support from DMR-1708280.