# Mathematical modeling of the feather follicle morphogenetic wave in birds

Maxime Estavoyer, Thomas Lepoutre, Marie Manceau

► To cite this version:

Maxime Estavoyer, Thomas Lepoutre, Marie Manceau. Mathematical modeling of the feather follicle morphogenetic wave in birds. Journal of Theoretical Biology, 2025, 612, 10.1016/j.jtbi.2025.112200. hal-04835554v3

HAL Id: hal-04835554

https://inria.hal.science/hal-04835554v3

Submitted on 27 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Highlights

Mathematical modeling of the feather follicle morphogenetic wave
in birds

Maxime Estavoyer, Thomas Lepoutre, Marie Manceau

* Mathematical modeling explores the propagation of feather follicle formation in birds.

* Model analysis predicts the speed of follicle emergence in a competent tissue.

* The predicted speed of competence acquisition depends only on cellular parameters.

* Follicle emergence is constrained by competence acquisition.

* Theoretical and numerical results estimate the overall speed of the wave.
# Mathematical modeling of the feather follicle morphogenetic wave in birds

Maxime Estavoyerª, Thomas Lepoutreª, Marie Manceauᵇ

## Abstract

During the development of the avian skin, feather follicles are produced in a medio-lateral morphogenetic wave that results in their spatial arrangement in typical patterns. This wave involves the timely acquisition of pattern-forming competence followed by a row-by-row production of feather follicles. While several mathematical models combining self-organizing systems accurately reproduced dynamics of feather follicle pattern formation, the events that control timely parameters of wave propagation remain poorly understood. Here, we built on previous modeling work to theoretically calculate the speed at which tissue competence progresses. Using a weakly non-linear analysis, we calculated the speed at which follicles emerge once competence is attained. We produced numerical simulations of our model to predict the respective influences of competence acquisition and follicle emergence on each other and on wave propagation. Our results show that the theoretical speed of follicle emergence is limited by competence acquisition, but that, in turn, competence acquisition is not constrained by follicle emergence. This modeling work provides an approximation of the timely parameters of the morphogenetic wave, and sheds light on the interplay between competence and patterning events in the developing skin.

## Keywords:

Morphogenetic wave, Feather follicle, Tissue competence, Mathematical modeling, Wave speed.

<footnote>ªInria, CNRS, Ecole Centrale de Lyon, INSA Lyon, Université Claude Bernard Lyon 1, Université Jean Monnet, ICJ UMR5208, 69603 Villeurbanne, France,</footnote>

<footnote>ᵇCentre for Interdisciplinary Research in Biology, Collège de France, Université PSL, CNRS, INSERM, France,</footnote>
# 1. Introduction

Many animals display striking geometrical patterns that spatially extend across organs or the entire body. The formation of these patterns is often characterized by the timely emergence of developing organs and appendages in so-called morphogenetic waves, such as observed in the sequential production of shark skin denticles (Cooper et al., 2018), zebra fish and corn snake scales (Evanitsky and Di Talia, 2023; Tzika et al., 2023), mouse and bat teeth (Sadier et al., 2019, 2023), or feathers (Ho et al., 2019; Bailleul et al., 2019). Morphogenetic waves result from the propagation, across the un-patterned tissue, of molecules (i.e., signaling waves) or mechanical cues (i.e., mechanical waves) that provide positional information in a progressive manner. Molecular and/or mechanical cues confer competence to the tissue: above a certain threshold value (e.g., molecular concentration, cellular density), groups of cells become competent to differentiate and shape into patterned elements. Visible patterning ensues, as structures emerge within the competent region in transient or final geometries. By timely pacing competence and patterning events, morphogenetic waves ensure that development occurs in a reproducible manner, thereby contributing to proper organ function and thus, individual fitness. The fundamental importance of morphogenetic waves has led to decades of theoretical and experimental studies of the mechanisms that trigger these waves, the nature of the molecular and mechanical signals that underlie their propagation, and how these signals control competence acquisition and patterning (Baker et al., 2009; Di Talia and Vergassola, 2022). However, the mechanisms that influence the speed of wave propagation, thereby ensuring proper positioning of patterned structures, remain unclear.

To address this question, we took advantage of the sequential formation of feather follicles, which are cutaneous organs in which bird's feathers are implanted. In the dorsum of many avian groups, including finches and poultry birds, feather follicles are spatially distributed in a regular hexagonal pattern where each feather follicle is surrounded by six neighbors (Sengel, 1976; Bailleul et al., 2019). This pattern forms in a morphogenetic wave characterized by two major steps. First, a medial-to-lateral signaling wave of Eda protein produces successive competent regions in the initially naive skin. Each competent region is a thin longitudinal segment characterized by overall thickening of the skin, higher cell density, and an increase in cell shape anisotropy along the antero-posterior axis (Bailleul et al., 2019; Ho
et al., 2019; Curantz et al., 2022). Longitudinal competent segments are also marked by the expression of β-catenin and Dermo-1 transcripts (Ho et al., 2019; Bailleul et al., 2019; Hornik et al., 2005; Dhouailly, 2009). The first competent segments are produced in the skin directly above the spinal cord (e.g., in the chicken embryo), or just lateral to it (e.g., in the quail and finch embryo); following segments are added laterally to the first-formed as the Eda wave progresses. Second, precursors of feather follicles, which are dermal cellular aggregates, individualize within competent segments, in a process involving a crosstalk between the underlying dermis and the overly-ing epidermis (Shyer et al., 2017; Ho et al., 2019; Sengel, 1976; Jiang et al., 2004). This crosstalk relies on self-organizing dynamics of reaction-diffusion operated by a Bmp2/Bmp7 activator-inhibitor couple and chemotaxis medi-ated by FGF signaling, and on mechanical cues, including intrinsic contrac-tile properties of dermal cells and FGF20-mediated epidermis compression (Palmquist et al., 2022; Michon et al., 2008; Painter et al., 2018; Ho et al., 2019; Jiang et al., 1999; Oster et al., 1983; Glover et al., 2017). The in-dividualization of feather follicles in each competent segment results in the sequential production of longitudinal follicle rows. This process repeats later-ally until the morphogenetic wave reaches the ventral boundary of the dorsal tract (See Figure 1A). Together, this work evidenced that feather follicle patterning occurs through a medial-to-lateral morphogenetic wave in which a progressive acquisition of competent segments (i.e., competence wave) is followed by a row-by-row sequence of follicle individualization within each segment (i.e., follicle emergence wave).

Previous theoretical work explored directional and sequential aspects of the feather follicle morphogenetic wave. In the model described in (Bailleul et al., 2019), itself modified from a prior model in (Michon et al., 2008), three equations are combined. Specifically, two reaction-diffusion equations (2)-(3) model the interaction of an activator and its inhibitor that diffuse at long range, which has been widely used to produce periodic patterns visibly resembling those observed in the coat of many animals (Milinkovitch et al., 2023). The third equation involves chemotaxis, which models cell migration towards a chemo-attractant (Hillen and Painter, 2009) and includes a logistic proliferation term (1), which models timely dynamics of cell proliferation in a developing tissue (Fisher, 1937; Kolmogorov, 1937):

$$ \partial_t n = D_n \Delta n - \nabla \cdot (\kappa n \nabla u) + \alpha_n n(1-n), \quad (1) $$
$$ \partial_t u = D_u \Delta u + \frac{\alpha_u n(1 + \omega u^2)}{(\beta_u^2 + u^2)(1 + v)} - \delta_u u, \quad (2) $$

$$ \partial_t v = D_v \Delta v + \alpha_v n u^2 - \delta_v v. \quad (3) $$

The species u and v correspond to the activator and the inhibitor in the reaction-diffusion equations; they diffuse with specific diffusion coefficients $D_u$ and $D_v$. The term n corresponds to the density of cells in the developing skin. The parameter $\kappa$ quantifies the intensity of chemotaxis (represented by the nonlinear term $\kappa \nabla \cdot (n\nabla u)$ in equation (1)); it corresponds to a biased diffusion where cells are attracted towards increasing concentrations of u. The cell proliferation rate is $\alpha_n$.

This model was studied mathematically: it admits, for certain values of parameters, a unique positive equilibrium point denoted as $(1, u_s, v_s)$. This homogeneous state can be Turing unstable, leading to the emergence of dotted arrays closely resembling feather follicle patterns. This result was consistent with a large body of work in which similar reaction-diffusion-chemotaxis models (e.g., (Bailleul et al., 2019; Diez et al., 2024; Tseng et al., 2024; Michon et al., 2008; Painter et al., 2018; Barnafi et al., 2022)) produced dotted motifs in silico. In addition, through changes in initial geometrical conditions, the model intrinsically reproduced species-specific patterning waves (in contrast with other studies in which the wave resulted for an external input (Painter et al., 2018; Ho et al., 2019)) as well as their final patterns. Biological predictions from this modeling work were that (i) the morphogenetic wave is initiated when the profile of initial conditions in the first segment is sharp enough, and (ii) propagates when cell density, which increases through proliferation, reaches a given threshold (Bailleul et al., 2019). Consistent with the latter, drug-induced inhibition of cell proliferation in cultured skin explants modulated the duration of wave propagation (Bailleul et al., 2019). Here, we took advantage of the intrinsic ability of the model to produce a wave to mathematically study its timely characteristics. To do so, we simplified the model described above, and conducted both a linear stability analysis and a weakly non-linear analysis to calculate the speed of follicle emergence within a competent segment. We then numerically studied the competence wave, and showed that its speed depends only on cellular parameters, and is not constrained by the speed of the follicle emergence wave. Conversely, the speed of follicle emergence locks onto the speed of competence acquisition. Finally, we propose an approximation of the speed of the morphogenetic wave in an initially naive tissue. Our results shed light on how competent
segment formation and follicle pattern emergence shape temporal dynamics of the morphogenetic wave.

## 2. Methods

### 2.1. A simplified mathematical model to compute the morphogenetic wave

To study the speed of the competence and follicle emergence waves, we modified the model described above, which was chosen for its ability to intrinsically produce sequential dot formation, but in which the reaction terms of equations (2) and (3) did not have clear biological interpretation. To simplify the reaction terms of equations (2) and (3) and avoid biologically unrealistic interactions, we used a truncation, which is commonly employed in biological pattern modeling (e.g., modeling of lizard or zebra fish coloration (Manukyan et al., 2017; Jahanbakhsh and Milinkovitch, 2022; Nakamasu et al., 2009) or mouse molar, mammalian oral palate ruggae, and shark skin denticle formation (Sadier et al., 2019; Cooper et al., 2018; Economou et al., 2020)). We chose a truncation similar to that introduced in (Nakamasu et al., 2009), as follows:

$$ egin{cases} \partial_t n = D_n \Delta n - \nabla \cdot (\kappa n \nabla u) + \alpha_n n (1-n), \\ \partial_t u = D_u \Delta u + n \chi_M (\gamma_u u - \gamma_v v + c_u) - \delta_u u, \\ \partial_t v = D_v \Delta v + n \chi_M (\eta_u u + c_v) - \delta_v v, \end{cases} $$

with $\chi_M(x) = x \mathbf{1}_{0 \le x \le M}(x) + M \mathbf{1}_{x>M}(x)$, where M is an arbitrary positive constant. With a suitable choice of parameters, this truncated model has the same homogeneous equilibrium and Turing instability as the original system (i.e., the Jacobian matrix around the steady state is the same as for (1)-(3)). More precisely, we selected the constants $\gamma_u$, $\gamma_v$, and $\eta_u$ such that the truncated model has the same unique positive equilibrium point $(1, u_s, v_s)$ and the same Jacobian matrix as the original model. These constants $\gamma_u$, $\gamma_v$ and $\eta_u$ are thus defined as follows:

$$ \gamma_u = \frac{2\alpha_u u_s (\omega\beta_u - 1)}{(\beta_u + u_s^2)^2 (1 + v_s)}, \quad \gamma_v = \frac{\alpha_u (1 + \omega u_s^2)}{(\beta_u + u_s^2) (1 + v_s)^2}, \quad \eta_u = 2\alpha_v u_s. \quad (4) $$

The constants $c_u$ and $c_v$ are defined by $c_u = -\gamma_u u_s + \gamma_v v_s + \delta_u u_s$ and $c_v = -\eta_u u_s + \delta_v v_s$.

Our choice of cut-off function preserves the positivity of $u$ and $v$. If the constant $M$ is assumed to be sufficiently large compared to the amplitudes of
the final stationary pattern, the cutoff function has no impact on the solution of the system. Since our theoretical analysis is based on small amplitude perturbations, we can simplify calculations to the following model without cutoff function:

$$ \left\{ \begin{aligned} \partial_t n &= D_n \Delta n - \nabla \cdot (\kappa n \nabla u) + \alpha_n n (1-n), \\ \partial_t u &= D_u \Delta u + n(\gamma_u u - \gamma_v v + c_u) - \delta_u u, \\ \partial_t v &= D_v \Delta v + n(\eta_u u + c_v) - \delta_v v. \end{aligned} \right. \qquad (5) $$

This simpler model, (5), provides the opportunity to obtain theoretical results. Because the morphogenetic wave is unidirectional (i.e., from the medial to the lateral part of the developing embryo), we studied model (5) in a one-dimensional space, thereby representing the dorso-medial axis of the spinal tract. Moreover, because the morphogenetic wave progresses bilaterally, we considered a space between the medial-most and lateral-most limits of the dorsal tract, namely a domain $\Omega = [0, \ell]$, where $\ell$ is a positive constant. Two mathematical steps in numerical simulations of the model represent the two main biological steps of the morphogenetic wave (Figure 1B). Specifically, for the competence wave, the increase in cell density is represented by a traveling wave connecting the point $(0, 0, 0)$, which represents the naive, un-patterned tissue, to the unique positive homogeneous equilibrium point $(1, u_s, v_s)$, which represents the last-formed competent segment. Setting initial cell density as $n = 0$ does not exactly reflect the biology of the naïve tissue, which contains cells, but this mathematical approximation simplifies the model for further analyses and does not drastically modify the result of simulations (see below). The follicle emergence wave is represented by dynamic oscillations in cell density that are caused by the Turing-instability of the point $(1, u_s, v_s)$. In this part of the curve, each peak represents a feather follicle (i.e., a locally dense cellular aggregate separated by regions of lower cell density).

## 2.2. Mathematical analyses of the model

Because the competence wave is intrinsic to our model, mathematical analyses allow exploring the events that shape the follicle emergence wave when competence is attained. We performed these analyses using a three-pronged methodology.
<figure><img src="image_2.png" /><figcaption>Figure 1: **The morphogenetic wave involves competence acquisition followed by follicle emergence.** A. Dorsal views of Japanese quail embryos and corresponding schematics of regions indicated with black squares marked for β-catenin transcripts from stage HH29 (left panel), when the first competent segments become visible (lined in green on the left half of the embryo) to HH30 (middle panel), when follicles individualize within the first segment to form a longitudinal row (in yellow), and HH32, when two follicle rows have been produced (in orange and yellow). The red arrow marks the direction of the morphogenetic wave; the blue dotted line marks the dorsal midline. B. Transverse sections and corresponding schematics of regions indicated with black squares in A, and marked with DAPI to visualize cell nuclei, show that skin histology changes as the wave progresses. The skin thickens in the competent area (in green) compared to the naive skin (in white). Follicles emerge sequentially as cells aggregate in the skin (in orange and yellow). C. Numerical simulation of the model (5) in which cell density *n* is represented along the medio-lateral axis (x) illustrate that during wave propagation (i.e., from left to right; red arrow), cell density, initially low in the naive tissue (right part of the curve, in white), gradually increases to a given threshold, which marks the acquisition of competence (visible here as a transition zone, in green). This is followed by the emergence of follicles, represented by regularly distributed peaks of high cell density (left part of the curve, in yellow).</figcaption></figure>
## Linear analysis

First, we performed a linear analysis to determine under which parameter conditions the positive equilibrium point is unstable and thus, produces patterns. For the equilibrium (0, 0, 0), which is always unstable, the linear analysis has been performed previously (Bailleul et al., 2019). We completed this work by studying the positive equilibrium point $(1, u_s, v_s)$ and the occurrence of a Turing bifurcation (Appendix A). The equilibrium point $(1, u_s, v_s)$ is stable without chemotaxis and diffusion. We determined the critical value of the chemotaxis term $\kappa_c$ at which the equilibrium point $(1, u_s, v_s)$ becomes unstable for all $\kappa > \kappa_c$. While chemotaxis is not required for pattern formation, it greatly facilitates the emergence of patterns: an increase in the chemotaxis value leads to a significant increase in the zone of parameter values where the positive equilibrium point $(1, u_s, v_s)$ exhibits a Turing-type instability (Fraga Delfino Kunz et al., 2024). In addition, increased chemotaxis causes a decrease in the regularity of stationary pattern, such as modifications in the number of follicles, fusions between adjacent follicles, or the emergence of ectopic follicles (Fraga Delfino Kunz et al., 2024). Thus, to obtain regular stationary patterns that mimic the observed arrangement of follicles in birds that display a morphogenetic wave, we opted for chemotaxis coefficients close to the critical parameter $\kappa \gtrsim \kappa_c$.

## Weakly non-linear analysis

Second, to estimate values of cell density that theoretically represent a formed feather follicle (i.e., the amplitude of peaks) and the speed of the follicle emergence wave, we included non-linear terms in the analysis. Solutions of a linear analysis indeed grow exponentially according to linear theory and do not have clear biological significance. Weakly non linear analyses allow deriving explicit quantitative information linking system parameters to pattern amplitude, thereby providing a powerful comparison between theoretical system expectations and biological observations (Cross and Hohenberg, 1993; Wollkind et al., 1994; Giunta et al., 2021).

We studied the system near the Turing bifurcation threshold. At this key bifurcation value, the pattern evolves on a slow time scale such that, using the multiple scales method, one can derive an evolutionary equation for pattern amplitude. In the case of one-dimensional space, a real Ginzburg-Landau type amplitude equation is obtained that theoretically admits traveling wave type solutions for which critical wave speed and wave steepness can be theoretically calculated. Similar analysis has been conducted for var-
ious models related to the system (5), such as non-linear reaction-diffusion systems (Bozzini et al., 2015; Klika et al., 2024), with cross-diffusion terms (Gambino et al., 2012; Tulumello et al., 2014; Gambino et al., 2013b; Giunta et al., 2021; Han and Dai, 2017), and for volume-filling chemotaxis model (Ma et al., 2016; Han et al., 2016, 2017).

The weakly non-linear analysis is detailed in Appendix B; we provide here a summary. We first translated the uniform steady state to the origin

$$ (\tilde{n}, \tilde{u}, \tilde{v}) = (n - 1, u - u_s, v - v_s). $$

To simplify reading, we omitted the tildes and referred to this vector as

$$ \mathbf{w} := (\tilde{n}, \tilde{u}, \tilde{v})^T = (n, u, v)^T. $$

The system (5) can thus be rewritten as the combination of a linear term and of a non-linear term

$$ \partial_t \mathbf{w} = \mathcal{L}\mathbf{w} + \mathcal{B}(\mathbf{w}, \mathbf{w}), \qquad (6) $$

where the linear part $\mathcal{L}$ is defined by

$$ \begin{gathered} \mathcal{L} = J + D^{\kappa} \partial_{xx}, \quad \text{with} \quad J := \begin{pmatrix} -\alpha_n & 0 & 0 \\ \delta_u u_s & \gamma_u - \delta_u & -\gamma_v \\ \delta_v v_s & \eta_u & -\delta_v \end{pmatrix}, \\ D^{\kappa} := \begin{pmatrix} D_n & -\kappa & 0 \\ 0 & D_u & 0 \\ 0 & 0 & D_v \end{pmatrix}, \end{gathered} $$

and the non-linear part by

$$ \mathcal{B}(\mathbf{w}, \mathbf{w}') := \frac{1}{2} \begin{pmatrix} -\kappa \partial_x (n' \partial_x u + n \partial_x u') - 2\alpha_n n' n \\ \gamma_u (n'u + nu') - \gamma_v (n'v + nv') \\ \eta_u (n'u + un') \end{pmatrix}. $$

To describe the spatial dynamics of the pattern, we introduced a slow time scale as follows

$$ t = \mathcal{T}(T_1, T_2, T_3, \dots), \quad T_i = \epsilon^i t, \quad i = 1, 2, \dots, $$

and a space scale by defining the characteristic length scale of the spatial modulation ($X = \epsilon x$), where $\epsilon \ll 1$ represents the distance from the Turing
bifurcation point as detailed mathematically in the Appendix B. Considering these changes in scale, we obtained

$$ \left\{ \begin{aligned} \partial_x &\rightarrow \partial_x + \epsilon \partial_X, \\ \partial_{xx} &\rightarrow \partial_{xx} + 2\epsilon \partial_{xX} + \epsilon^2 \partial_{XX}, \\ \partial_t &\rightarrow \epsilon \partial_{T_1} + \epsilon^2 \partial_{T_2} + \epsilon^3 \partial_{T_3} + \dots \end{aligned} \right. \qquad (7) $$

We expanded $\kappa$ and $\mathbf{w}$ in the small parameter as

$$ \left\{ \begin{aligned} \kappa &= \kappa_c + \epsilon\kappa_1 + \epsilon^2\kappa_2 + \epsilon^3\kappa_3 + \dots, \\ \mathbf{w} &= \epsilon\mathbf{w}_1 + \epsilon^2\mathbf{w}_2 + \epsilon^3\mathbf{w}_3 + \dots. \end{aligned} \right. \qquad (8) $$

where $\mathbf{w}_i = (n_i, u_i, v_i)^T$. Using the relations (7) and (8), we expanded the linear and non-linear parts in terms of powers of $\epsilon$, which yields for fixed $f$

$$ \mathcal{L}(f) = \mathcal{L}^0(f) + \epsilon\mathcal{L}^1(f) + \epsilon^2\mathcal{L}^2(f) + \epsilon^3\mathcal{L}^3(f) + \dots, \qquad (9) $$

and,

$$ \mathcal{B}(f,g) = \mathcal{B}^0(f,g) + \epsilon\mathcal{B}^1(f,g) + \epsilon^2\mathcal{B}^2(f,g) + \dots. \qquad (10) $$

where $(\mathcal{L}^i)_i$ and $(\mathcal{B}^i)_i$ are operators specified in Appendix B. We expanded each term of equation (6) to identify the terms of each power of $\epsilon$, which results in

$$ \begin{aligned} \partial_t \mathbf{w} &= 0, &+\epsilon.0, &+\epsilon^2 \partial_{T_1} \mathbf{w}_1, &+\epsilon^3 (\partial_{T_1} \mathbf{w}_2 + \partial_{T_2} \mathbf{w}_1) + O(\epsilon^3), \\ \mathcal{L}(\mathbf{w}) &= 0, &+\epsilon \mathcal{L}^0(\mathbf{w}_1), &+\epsilon^2 (\mathcal{L}^0(\mathbf{w}_2) + \mathcal{L}^1(\mathbf{w}_1)), &+\epsilon^3 (\mathcal{L}^0(\mathbf{w}_3) + \mathcal{L}^1(\mathbf{w}_2) + \mathcal{L}^2(\mathbf{w}_1)) + O(\epsilon^3), \\ \mathcal{B}(\mathbf{w}, \mathbf{w}) &= 0, &+\epsilon.0, &+\epsilon^2 \mathcal{B}^0(\mathbf{w}_1, \mathbf{w}_1), &+\epsilon^3 (2\mathcal{B}^0(\mathbf{w}_1, \mathbf{w}_2) + \mathcal{B}^1(\mathbf{w}_1, \mathbf{w}_1)) + O(\epsilon^3). \end{aligned} $$

To ensure a matching up to order 3, it was thus necessary to solve for $\mathbf{w}_{1,2,3}$ the system

$$ 0 = \mathcal{L}^0(\mathbf{w}_1) + \qquad 0, \qquad (11) $$

$$ \partial_{T_1} \mathbf{w}_1 = \mathcal{L}^0(\mathbf{w}_2) + \mathcal{L}^1(\mathbf{w}_1) + \mathcal{B}^0(\mathbf{w}_1, \mathbf{w}_1), \qquad (12) $$

$$ (\partial_{T_1} \mathbf{w}_2 + \partial_{T_2} \mathbf{w}_1) = \mathcal{L}^0(\mathbf{w}_3) + \mathcal{L}^1(\mathbf{w}_2) + \mathcal{L}^2(\mathbf{w}_1) + 2\mathcal{B}^0(\mathbf{w}_1, \mathbf{w}_2) + \mathcal{B}^1(\mathbf{w}_1, \mathbf{w}_1). \qquad (13) $$
At $O(\epsilon)$, we recovered the linear problem $\mathcal{L}^0(\mathbf{w}_1) = 0$ whose solution satisfying the Neumann boundary conditions is

$$ \mathbf{w}_1 = A(T_1, T_2, \dots; X)\rho \cos(k_c x), \quad (14) $$

where $k_c$ corresponds to the first and unique unstable mode and $A(T_1, T_2, \dots; X)$ is the spatial amplitude of follicles. The vector $\rho := (1, \rho_u, \rho_v)$ is derived in Appendix B and given by the relation (B.1).

At $O(\epsilon^2)$, we combined formula (14) and the Fredholm alternative to determine the following equalities

$$ T_1 = 0, \quad \kappa_1 = 0, \quad (15) $$

as well as the definition of the vector $\mathbf{w}_2$, defined as follows

$$ \mathbf{w}_2 = A^2 \mathbf{w}_{20} + A^2 \mathbf{w}_{22} \cos(2k_c x) + \frac{\partial A}{\partial X} \mathbf{w}_{21} \sin(k_c x), \quad (16) $$

where $\mathbf{w}_{2i}$, $i = 0, 1, 2$ are the solutions of the system (B.8).

Using the resolvability condition and equality (16) at $O(\epsilon^3)$, we found the following real Ginzburg-Landau equation for amplitude $A$

$$ \partial_{T_2} A = \nu \partial_{XX} A + \sigma A - LA^3. \quad (17) $$

where the parameters $\nu$, $\sigma$ and $L$ are given by the equations (B.9), (B.10) and (B.11), respectively. Since the growth rate coefficient $\sigma$ is always positive, the dynamics of the real Ginzburg-Landau equation (17) could be divided into two cases depending on the sign of the constant $L$, namely the supercritical case, when $L$ is positive, and the subcritical case, when $L$ is negative.

When $L < 0$, the real Ginzburg-Landau equation (17) was unable to capture the amplitude of cell density, and it was necessary to increase the weakly non-linear expansion to order five. However, in this case, the expected amplitude $\epsilon A$ is of order $O(1)$ (Giunta et al., 2021). Thus, the approximation proposed by the weakly non-linear analysis did not generally match the asymptotic solutions. This lack of precision has been observed in many studies (e.g., (Gambino et al., 2012; Han et al., 2017; Giunta et al., 2021)). In addition, the theoretical existence of a positive traveling wave for the complex equation obtained at order 5 upon adding the slow modulation in space has not been established (Han et al., 2016). We therefore chose to focus on the case where $L$ is positive.
The weakly non-linear analysis allowed drawing predictions on the stationary pattern and the speed of the follicle emergence wave. Specifically, the asymptotic expression for the spatial pattern of (5) when $L > 0$ is given by

$$ \begin{pmatrix} n \\ u \\ v \end{pmatrix} = \begin{pmatrix} 1 \\ u_s \\ v_s \end{pmatrix} + \epsilon A(\epsilon x, \epsilon^2 t) \cos(k_c x) \begin{pmatrix} 1 \\ \rho_u \\ \rho_v \end{pmatrix} + O(\epsilon^2), \quad (18) $$

where the amplitude $A$ is the solution of the Eq. (17). The Eq. (17) admits two non-negative points equilibrium, namely the unstable point $A = 0$ and the stable positive point $A_\infty := \sqrt{\sigma/L}$, which is globally asymptotically stable, leading to $\lim_{T_2 \to +\infty} A(X, T_2) = A_\infty$ for all $X$. Thus, for the final stationary pattern, the weakly non-linear analysis at order 3 resulted in

$$ \begin{pmatrix} n(x) \\ u(x) \\ v(x) \end{pmatrix} = \begin{pmatrix} 1 \\ u_s \\ v_s \end{pmatrix} + \epsilon A_\infty \rho \cos(k_c x) + \epsilon^2 A_\infty^2 (\mathbf{w}_{20} + \mathbf{w}_{22} \cos(2k_c x)) + O(\epsilon^3). \quad (19) $$

For positive $L$, equation (17) is of monostable KPP type and its associated propagation speed (in rescaled variables) is given by $c_{\text{WNL}}^* := 2\sqrt{\nu\sigma}$ (Kolmogorov, 1937; Aronson and Weinberger, 1975). The solution propagates as a traveling wave front, with a steepness denoted $\mu$, given by $\mu = \sqrt{\sigma/\nu}$. To determine the speed of the follicle emergence wave in a competent area at the correct scale (recall in (17) $X = \epsilon x, T_2 = \epsilon^2 t$), we rescaled to $c_{\text{WNL}}^*$ by $\epsilon$, resulting in

$$ c_{\text{comp}\to\text{FE}} = \epsilon c_{\text{WNL}}^* = 2\epsilon \sqrt{\nu\sigma}. $$

### 3. Results

#### 3.1. Mathematical analyses of the model predict the speed of the follicle emergence wave in a competent region

*Results of the weakly non-linear analysis:* The weakly non-linear analysis provided predictions for three major characteristics of the morphogenetic wave near the Turing bifurcation (i.e. for $\epsilon \ll 1$) (Figure 2).

The amplitude of cell density at each peak $A_{\text{CD}}$ mathematically corresponds to the distance between maximal cell density in the patterned region and cell density in the un-patterened region of the curve. Biologically, this value represents the density of cellular aggregates within emerging follicles,
in the two-dimensional plane of the early developing skin (feather follicles later develop as 3D outgrowths, but this occurs after the propagation of the patterning wave).

$$A_{CD} := \epsilon \sqrt{\sigma/L}. \qquad (20)$$

The speed of the follicle emergence wave in a competent region, which we define as $c_{comp→FE}$, is mathematically given by the propagation speed of the pattern envelope (by analogy, in the next part of the article, the speed of the competence wave is defined as $c_{naive→comp}$). Biologically, $c_{comp→FE}$ represents how fast follicles emerge in a competent skin region.

$$c_{comp→FE} := 2\epsilon\sqrt{\nu\sigma}. \qquad (21)$$

The steepness of the wave front $\mu$ (in red) mathematically corresponds to the exponential decay rate of the wavefront tail. Biologically, it represents the sharpness of the transition between a patterned and an un-patterned, competent region (e.g., the steeper the curve, the less follicles will emerge at the same time along the one-directional patterning space). Indirectly, the sharpness of the wavefront thus indicates to what extent follicles emerge in a sequential fashion along the dorso-ventral axis.

$$\mu = \sqrt{\sigma/\nu}. \qquad (22)$$

The parameters $\sigma$, $\nu$, $L$ are given, respectively, by equations (B.9), (B.10), and (B.11). The parameter $\epsilon$ represents the relative distance to the Turing bifurcation, mathematically defined as $\epsilon := \sqrt{(\kappa - \kappa_c)/\kappa_c}$, where $\kappa_c$ corresponds to the critical value at which the Turing instability occurs. Since $\epsilon$ is influenced by chemotaxis, an increase in chemotactic sensitivity enhances cell aggregation in a given follicle, leading to higher values of cell density amplitude $A_{CD}$ and an increase in the speed of the follicle emergence wave in a competent region $c_{comp→FE}$.

*Numerical validation:* Third, we validated the weakly non-linear analysis using numerical simulations of the model (equation (5)). To test predictions on the amplitude of cell density that characterizes feather follicles, we compared the approximations of the weakly non-linear analysis with numerical simulations of model (5) solutions at stationary pattern state. We found that when $\epsilon$ is small enough, the approximation (19) predicted by the weakly non-linear analysis is accurate (Figure D.9A). Similarly, for $\epsilon = 0.3$, the analysis remained close to the numerical solution (Figure D.9B). We computed
<figure><img src="image_3.png" /><figcaption>Figure 2: **Results of the weakly non-linear analysis allow calculating characteristics of the follicle emergence wave.** Results of the weakly non-linear analysis allow predicting the amplitude of cell density at each peak $A_{CD}$ (in black), the speed of the follicle emergence wave $c_{comp→FE}$ (in blue) and the steepness of the wave front $\mu$ (in red).</figcaption></figure>

the bifurcation diagram of model (5) using the Matlab package `pde2path` (Uecker et al., 2012; Uecker, 2021). The code was adapted from the tutorial (Uecker, 2019). We found that the bifurcation diagram validates both numerical simulations (see Figure D.9A-B) and the precision of the weakly non-linear analysis (Figure D.9C).

To test predictions on the speed of the follicle emergence wave in a competent area, we compared the result of the weakly non-linear analysis with numerical solution of the system (5) over time. We chose an initial amplitude of cell density given by $A_0(X) = A_\infty \cdot 1_{|X|\le K}$, with $K$ being a positive constant, because with these initial values, the speed of the morphogenetic wave corresponds to the follicle emergence speed $c_{comp→FE}$ (Kolmogorov, 1937). For a small $\epsilon$ ($\epsilon = 0.05$), we found that the approximation obtained by the weakly non-linear analysis closely matches numerical results (Figure D.10). For $\epsilon = 0.2$, as expected, the weakly non-linear analysis lost precision. However, the approximation of the speed of the follicle emergence wave remained precise. In particular, we noticed an overshoot phenomenon: when the inhomogeneity appeared, it tended to increase quickly up to a given maximum, then decrease toward a stationary amplitude close to the amplitude predicted by the weakly non-linear analysis, $A_{CD}$. In other words, far enough from the Turing bifurcation, the amplitude of the wavefront became non-monotonic (Figure D.11). Together, these simulations show that close to the Turing bifurcation, approximations of the weakly non-linear analysis accurately reproduce solutions of the model (5) for the amplitude of follicle cell density
and the speed of the follicle emergence wave in a competent region.

## 3.2. Chemical and cellular parameters influence the follicle emergence wave

To predict which biological cues modulate timely dynamics of the morphogenetic wave, we tested the impact of various parameters of the model (5) on the parameters $\sigma$, $L$ and $\nu$, which are defined explicitly in our system (17) and influence the stationary amplitude of follicle cell density $A_{CD}$, the propagation speed of the follicle emergence wave $c_{comp \to FE}$ and the steepness of the wave front $\mu$. For each parameter variation, we recalculated $\kappa_c$.

**Amplitude of follicle cell density ($A_{CD}$):** We found that when the Turing bifurcation is supercritical, the final amplitude of follicle cell density is given by $A_{CD} = \epsilon\sqrt{\sigma/L}$, with $\sigma$ and $L$ defined respectively by (B.10) and (B.11). We simulated the impact of each parameter of the system (5) on $A_{CD}$ for a fixed $\epsilon$ (Figure D.12) and found that an increase in the parameters $D_u$, $\alpha_n$, $\delta_u$, $\delta_v$, and $\gamma_v$ leads to an increase in the amplitude of feather follicle cell density $A_{CD}$. Conversely, an increase in the parameter $\gamma_u$ reduced this amplitude. The parameters $D_v$ and $\eta_u$ correlated positively or negatively depending on their values. The diffusion coefficient of dermal cells $D_n$ had a negligible effect on cell density amplitude $A_{CD}$.

**Propagation speed of the follicle emergence wave in a competent area ($c_{comp \to FE}$):** We found that increasing the diffusion coefficient $D_n$ leads to faster propagation speed, with an initial exponential increase that then tends to zero. Similarly, increasing $D_u$, $\delta_u$, $\gamma_v$, and $\eta_u$ increased propagation speed while increasing $D_v$, $\delta_v$ and $\gamma_u$ decreased it. Thus, in the theoretical context of our model, both molecular and cellular parameters modulate the speed of the follicle emergence wave $c_{comp \to FE}$. Interestingly, however, the impact of the proliferation rate $\alpha_n$ on the speed of the follicle emergence wave $c_{comp \to FE}$ was minor. Experimentally reducing cell proliferation using drug treatments on cultured skin explant lengthens the duration of the morphogenetic wave (Bailleul et al., 2019; Curantz and Manceau, 2021). Our simulation results suggest that this experimental observation is not due to an effect of cell proliferation on the speed of the follicle emergence wave, but rather to an effect on competence acquisition.

**Steepness of the follicle emergence wave in a competent area ($\mu$):** We found that the parameters $D_v$, $\delta_u$, $\delta_v$, $\gamma_u$, $\gamma_v$ and $\eta_u$ had a similar influence than on the speed of the wave $c_{comp \to FE}$. However, varying the parameter $D_n$ had the opposite effect: when it increased, the wave front became steeper.
The steepness of the wave front $\mu$ was proportional to the proliferation pa- rameter $\alpha_n$, and dynamic with respect to the parameter $D_u$, with an initial increase followed by a decrease (Figure D.14). These results suggest that high values of cell proliferation are necessary to produce a sharp transition between patterned and competent regions, and thus, strict sequentiality in the emergence of follicles along the dorso-ventral axis.

### 3.3. The speed of the competence wave depends only on cellular parameters

We thus far studied our system when the tissue is already competent, that is, near the unstable equilibrium point $(1, u_s, v_s)$ and with initial con- ditions of high cell density. To predict the behavior of the morphogenetic wave in a naive tissue, we studied properties of our model with initial condi- tions of cell density lower than in the competent region. To do so, we used initial data with finite support (i.e., cell density equals zero outside of the competence interval). The transition between the naive region and the com- petent region was defined as half the value of cell density in the competent region (i.e., above this value, the tissue is considered competent). We thus mathematically defined the naive region by $\{x: n(t, x) < 1/2, x \in \Omega\}$. Nu- merically, our system (5) admits traveling wave solutions connecting the null point $(0, 0, 0)$ to the equilibrium point $(1, u_s, v_s)$. Intuitively, in conditions of very low cell density, molecular factors are not produced and/or maintained ($u$ and $v$ do not vary) and therefore, chemotaxis has no effect (i.e., $\kappa = 0$). The competence wave may thus be reduced to the Fisher-KPP equation:

$$ \partial_t n = D_n \Delta n + \alpha_n n(1 - n). \qquad (23) $$

For this equation, traveling wave type solutions have the following speed $c_{\text{F-KPP}}^* := 2\sqrt{D_n\alpha_n}$. For initial data with finite support, the selected speed of the wave of the associated Fisher-KPP equation (23) corresponds to the (minimal) Fisher-KPP speed (Kolmogorov, 1937). We thus produced nu- merical simulations of the morphogenetic wave using initial data with finite support and compared them with results of the Fisher-KPP equation.

We first studied the competence wave in isolation, independently of the formation of follicles. To do so, we produced numerical simulations around the equilibrium point $(1, u_s, v_s)$ when it is stable (i.e., when feather follicles can not emerge from Turing instability). We found that numerical simula- tions of our system produce a competence wave that has a speed and front very similar to those obtained with simulations of the Fisher-KPP equation
(23) (Figure 3). This validates our mathematical approximations of the state of a naive tissue and suggests that alone, the speed of the competence wave $c_{naive→comp}$ depends only on cell density:

$$ \text{Speed of the competence wave: } c_{naive→comp} := 2\sqrt{\alpha_n D_n}. \qquad (24) $$

### 3.4. The speed of follicle emergence is limited by the speed of the competence wave

We next tested the influence of follicle emergence on competence acquisition. To do so, We produced numerical simulations of our model in conditions of initial data with finite support in which we varied the chemotaxis parameter. As seen in the weakly non-linear analysis, this modifies the amplitude of follicle cell density $A_{CD}$, consistent with the observation that *in vivo*, increased chemotaxis causes the formation of denser follicles (Ho et al., 2019). We found that under these conditions, the competence wave is not affected by changes in the chemotaxis parameter and follows the behavior of the Fisher-KPP solution (Figure 4). Thus, follicle emergence does not impact the speed of the competence wave, suggesting that in living tissues, the molecular and cellular events that trigger and shape follicle formation are independent from, or are active later than, those responsible for competence acquisition.

### 3.5. The emergence of follicles in the wake of the competence wave is influenced by competence acquisition

In poultry birds and finches, follicles individualize in a longitudinal row concomitantly to competence acquisition in the neighboring lateral segment (i.e., follicle emergence and competence acquisition are spatially coupled). We thus tested numerically the production of follicles in the wake of competence. In addition, *in vivo*, each competent segment is a thin longitudinal band. We thus numerically assessed the spatial extent of the competence wave. To do so, we simulated differences in the respective speeds of the competence and follicle emergence waves.

In the first scenario, the speed of the follicle emergence wave in an initially competent tissue $c_{comp→FE}$ is greater than that of the competence wave $c_{naive→comp}$. We found that with such conditions, model solutions are more complex: around the unstable equilibrium point $(1, u_s, v_s)$, heterogeneous patterns form in the wake of the competence wave, displaying differences
<figure><img src="image_4.png" /><figcaption>Figure 3: In isolation, the *competence* wave progresses at a speed similar to solutions of the Fisher-KPP equation. Snapshots of a numerical simulation of our model (5) (black line) at different times show that when the equilibrium point (1, *u*<sub>s</sub>, *v*<sub>s</sub>) is stable (i.e., in the absence of follicle emergence), the *competence* wave (in green) progresses at a speed and with a front similar to those produced by numerical simulations of the Fisher-KPP equation (23) (red dashed line). Parameter values: *D*<sub>n</sub> = 7e-5, *D*<sub>u</sub> = 6e-3, *D*<sub>v</sub> = 0.13, *γ*<sub>u</sub> = 14.5389, *γ*<sub>v</sub> = 0.099932, *η*<sub>u</sub> = 3931, *c*<sub>u</sub> = 2.9672, *c*<sub>v</sub> = -1104, *α*<sub>n</sub> = 0.001 and *κ* = 4.84e-5 < *κ*<sub>c</sub>.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Case 2 – Intermediate chemotaxis</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Case 3 – High chemotaxis</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: The speed of the competence wave is unaffected by follicle formation.</figcaption></figure>

Three independent numerical simulations of model (5) (black lines) show that varying the chemotaxis parameter $κ$ such that either no follicles emerge (Case 1; $κ = 2.1e - 5 < κ_c$), low-amplitude follicles emerge (Case 2; $κ = 2.17e - 5 \gtrsim κ_c$) or high-amplitude follicles emerge (Case 3; $κ = 2.42e - 5 > κ_c$) does not modify the competence wave, which follows the behavior of the Fisher-KPP solution (red dashed lines). Parameter values: $D_n = 3e-5$, $D_u = 6e-3$, $D_v = 0.13$, $γ_u = 14.5389$, $γ_v = 0.099932$, $η_u = 3931$, $c_u = 2.9672$, $c_v = -1104$, $α_n = 0.001$.
in the amplitude of follicle cell density between the naive region, the com-
petent region (in which follicles may start to form), and the patterned re-
gion (in which follicles are formed and well-established). Specifically, when
$c_{\text{comp}\to\text{FE}} > c_{\text{naive}\to\text{comp}}$, the spread of peaks locked with the competence
wave, which was itself very close to the speed of the solution of the associ-
ated Fisher-KPP equation (23) (Figure 5). This phenomenon occurred for
various sets of parameter values. The size of the competent region remained
small and constant throughout the simulation (Figure D.15). Thus, when
$c_{\text{comp}\to\text{FE}} > c_{\text{naive}\to\text{comp}}$, the overall speed of the morphogenetic wave, which
we can define as $c_{\text{morph}}$, corresponds to the speed of the competence wave:
$c_{\text{morph}} = c_{\text{naive}\to\text{comp}} < c_{\text{comp}\to\text{FE}}$, and the competent region remains limited
in size. From a biological perspective, these results are expected, as follicles
can not form without competence. They also predict that the speed of com-
petence acquisition is not influenced by follicle formation, consistent with
conclusions above (see Figure 4).

In a second scenario, the speed of the competence wave $c_{\text{naive}\to\text{comp}}$ is greater than the speed of the follicle emergence wave $c_{\text{comp}\to\text{FE}}$. In this case, we found that the overall speed of the morphogenetic wave $c_{\text{morph}}$ was similar to predictions from the weakly non-linear analysis, and that much like in tests above, the competence wave connecting $(0,0,0)$ to $(1, u_s, v_s)$ was very close to the Fisher-KPP wave. In addition, the weakly non-linear analysis described precisely the emergence of follicles in the wake of the competence wave. Because the speed of the competence wave was higher than follicle emergence, the competent region grew continuously, corresponding mathematically to the formation of a “terrace”. The latter had no impact on the speed of follicle emergence (Figure 6). Together, these results predict that increasing the speed of competence acquisition does not modify the speed of follicle emergence but increases the extent of the competent region. *In vivo*, large competent regions dynamic in size are never observed: competent segments are thin, and dorso-ventral sequentiality is strict. This suggests that the second scenario described here may not accurately reflect biological reality.

3.6. *Comparing theoretical calculations and numerical results allows formu-
lating a formula for the speed of the morphogenetic wave*

To compute the overall speed of the morphogenetic wave, we plotted the
theoretical speeds of the competence and follicle emergence waves, respec-
tively predicted by the Fisher-KPP and weakly non-linear analysis, relative to
<figure><img src="image_8.png" /><figcaption>Figure 5: **When it is faster, the follicle emergence wave locks with the competence wave.** Snapshots of a numerical simulation of the system (5) (black line) at three different time points (t=7500, 15000 and 30000) when the point (1, $u_s$, $v_s$) is Turing unstable and when $c_{naive→comp} < c_{comp→FE}$ show that the speed of the follicle emergence wave (in the yellow region) corresponds to the speed of the competence wave connecting $n=0$ to $n=1$ (from the white to the dorsal limit of the green region). The competence wave is similar to solutions of the Fisher-KPP equation (23) (red dashed line). Parameter values: $D_n = 7e-5$, $D_u = 6e-3$, $D_v = 0.13$, $\gamma_u = 14.5389$, $\gamma_v = 0.099932$, $\eta_u = 3931$, $c_u = 2.9672$, $c_v = -1104$, $\alpha_n = 0.001$ and $\kappa = 5.5e-5 > \kappa_c$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6: **When it is faster, the competence wave creates larger follicle-free competent regions.** Snapshots of a numerical simulation of the system (5) (black line) at three different time points ($t = 24000, 42500$ and $82500$) when the point $(1, u_s, v_s)$ is Turing unstable and when $c_{naive\tocomp} \ge c_{comp\toFE}$ show that the speed of the follicle emergence wave (in the yellow region) is slower than in the reverse situation, and fits the theoretical envelope calculated using the equation (18) given by the weakly non-linear analysis (blue dotted line). The competence wave, similar to solutions of the Fisher-KPP equation (23) (red dashed line), extends over a larger and expanding competent region (in green). Parameter values: $D_n = 1.5e-5$, $D_u = 6e-3$, $D_v = 0.13$, $\gamma_u = 14.5389$, $\gamma_v = 0.099932$, $\eta_u = 3931$, $c_u = 2.9672$, $c_v = -1104$, $\alpha_n = 0.001$ and $\kappa = 4.9e-5 > \kappa_c$.</figcaption></figure>
the proliferation of skin cells $\alpha_n$, which influences both, and compared those with results obtained in numerical simulations. When proliferation is low, and competence acquisition is therefore slower than follicle emergence, the morphogenetic wave progresses at the speed of the competence wave (follicle emergence locks with the progression of competence): $c_{morph} = c_{naive \to comp}$ (Figure 7, left part). As proliferation increases, the speed of both competence and follicle emergence waves increases. However, the speed of competence acquisition grows significantly faster with proliferation than the speed of follicle emergence. When proliferation becomes sufficiently large, competence acquisition becomes faster than follicle emergence, and in this case, the speed of the morphogenetic wave is constrained by the speed of follicle emergence, which is the same as the speed predicted by the weakly non-linear analysis for follicle emergence in an already competent region: $c_{morph} \approx c_{comp \to FE}$ (Figure 7, right part). These results hold for different values of the parameters $\alpha_n$ and $D_n$ (Figure D.16).

In both cases, we obtain an accurate approximation of the speed of the morphogenetic wave, given by:

$$ \text{Speed of the morphogenetic wave: } c_{morph} \approx \text{minimum}(c_{naive \to comp}, c_{comp \to FE}). \qquad (25) $$

## 4. Discussion

In this study, we produce a simpler model modified from (Bailleul et al., 2019) that provides the opportunity to perform mathematical validation of theoretical results and propose an approximation of the speed of the follicle morphogenetic wave under two conditions, namely when competence acquisition is faster than follicle emergence, and vice versa. We found that in conditions of low cell proliferation, when follicles emerge faster than the acquisition of competence, they catch up and synchronize with the progression of competence, without impacting it. In that case, the speed of the morphogenetic wave aligns with that of competence acquisition, which in vivo, would depend only on cell density, and is mathematically given by $c_{comp} = 2\sqrt{D_n\alpha_n}$, with $D_n$ being the diffusion coefficient and $\alpha_n$ the dermal cell proliferation rate. Conversely, when follicle emergence is slower than competence, it aligns with the speed predicted by the weakly non-linear analysis (i.e., the speed of follicle emergence in an already competent region). In this scenario, competence acquisition does not directly impact follicle emergence, and the overall
<figure><img src="image_10.png" /><figcaption>Figure 7: Morphogenetic wave speed transitions from competence acquisition speed to follicle emergence speed depending on proliferation rate. Plotting the theoretical speeds of the competence wave (i.e., corresponding to the Fisher KPP equation, in red) and the follicle emergence wave in a competent region (corresponding to the results of the weakly non-linear analysis, in blue) relative to cell proliferation rate allows comparing theoretical results with the behavior of the morphogenetic wave observed in numerical observation (in black). When $α_n$ is low, the speed of the morphogenetic wave corresponds to the speed of the competence wave. When $α_n$ reaches a threshold value, the speed of the morphogenetic wave corresponds to the speed of the follicle emergence wave in a competent region.</figcaption></figure>
speed of the morphogenetic wave approximates with the speed of follicle
emergence. *In vivo*, the speed of the morphogenetic wave would thus likely
depend on a combination of increased cell density and other factors, such as
follicle-forming factors. Together, this work provides groundwork for further
mathematical and experimental analyses.

Mathematically, our results demonstrate that when competence acquisi-
tion is faster than follicle emergence, studying the speed of morphogenetic
wave is akin to studying the speed of follicle emergence within a competent
zone (since the speed of the competence wave does not affect the speed of
follicle emergence). Here, we chose to do the latter using a weakly non-linear
analysis, which allows obtaining a precise approximation of the envelope,
whose shape represents the timely evolution of pattern emergence. Specifi-
cally, in our analysis, it is necessary for the parameter *L* to be positive (i.e.,
supercritical case) and to be sufficiently close to the critical bifurcation. The
envelope of equation (17) admits traveling wave solutions that have a prop-
agation speed $c_{\text{comp}\to\text{FE}} = 2\epsilon\sqrt{\sigma\nu}$. From this analysis, we derived the speed
of the morphogenetic wave as given by $c_{\text{morph}} \approx \min(2\sqrt{\alpha_n D_n}, 2\epsilon\sqrt{\sigma\nu})$.

It is important to note several limitations to that prediction. First, while
our study focuses primarily on the asymptotic speed of the competence wave,
such a regime may not be fully achieved in small domains (i.e., over very
short biological distances). Second, though we obtained a good approxima-
tion based on previous results Bailleul et al. (2019), it is challenging to link
the predicted speed of the patterning wave with real biological values of wave
speed due to (i) the variability in wave duration observed across individuals
and avian species Bailleul et al. (2019), and (ii) the difficulty of linking exper-
imental data with continuous PDE-based models. This prediction however
provides a framework to test how the speed of the patterning wave is modi-
fied upon experimental perturbation. Third, to allow for analytical progress,
our model simplifies several biological processes: the underlying molecular
process is reduced to minimal interactions, the cell movement term assumes
a linear chemotactic response, and cell proliferation follows logistic growth.
The latter has been experimentally verified Bailleul et al. (2019), but not the
first two. Fourth, we approximate cell density to zero ($n = 0$) in the naive
tissue, which does not reflect biological reality.

Future work based on our methodology may involve extending analyses
to a two dimensional setting. In this condition, new phenomena appear,
such as rhombi, hexagons or mixed-mode patterns, which can also be pre-
dicted by a weakly non-linear analysis (Bozzini et al., 2015; Gambino et al.,
2013a; Cruywagen et al., 1997; Gambino et al., 2022; Ma et al., 2019). It would also be interesting to integrate other mathematical methods. For example, an alternative methodology to a weakly non-linear analysis would be a well-described marginal stability criterion (Dee and Langer, 1983; Van Saarloos, 2003; Holzer and Scheel, 2014; Ben-Jacob et al., 1985; Van Saarloos and Hohenberg, 1992; Myerscough and Murray, 1992). For a system of two reaction-diffusion equations, the application of this method, close to a supercritical bifurcation, predicts the same propagation speed as a weakly non-linear analysis (Klika et al., 2024). However, the weakly non-linear analysis provides insights into the amplitude of the pattern and the transient profile of the wave that are not captured by a marginal stability method. It is also possible to consider modulated front solutions, in which a front connects an unstable homogeneous stationary state to a spatially periodic pattern. This has been initially studied for the Swift-Hohenberg equation with cubic nonlinearity (Collet and Eckmann, 1986, 1990; Eckmann and Wayne, 1991), before having been extended to other models such as the non-local Fisher-KPP equation (Faye and Holzer, 2015). Another useful tool may be the theory of a “triggered/quenched front”, which combines a pattern-forming front with a wave of external competence (Goh and Scheel, 2014, 2016; Avery et al., 2019; Goh and Scheel, 2018, 2023). Altogether, this mathematical work provides a simple and computable framework to study periodic pattern formation in biological systems other than the feather follicle pattern, as it can be easily adapted to other periodic patterns observed in numerous species and at different stages of embryonic development (e.g., (Sudderick and Glover, 2024)).

Biologically, our work provides a number of testable predictions. First, it suggests that the progression of the competence wave is independent from the speed of follicle emergence. This suggests that, provided biological information and material necessary for follicle formation are available in the tissue, the speed of the competence wave largely sets the overall speed of the morphogenetic wave. In other words, the mechanisms driving the progression of the morphogenetic wave in a naive tissue are independent from the molecular and/or mechanical cues that drive follicle emergence within a confined patterning space. This hypothesis is intuitively compatible with a molecular basis of competence acquisition, consistent with the observation that the signaling molecule Eda creates a traveling wave concomitant with competence acquisition during feather follicle formation (but also in other systems such as rodent teeth). In any case, experimentally testing this hy-
pothesis would imply perturbing competence acquisition without affecting follicle emergence, which may be possible using inhibitors or enhancers of Eda activity as in (Ho et al., 2019). While differences in the extent of the patterned region where observed upon treatment, the effect of perturbations on the respective sizes of the competent and patterned regions was not assessed. It would be interesting to perform similar experiments and study the effect of Eda signaling perturbation on cellular characteristics that mark the un-patterened, competent region, and that are easily measurable with simple immuno-histological stains, such as an increase in tissue thickness, cell density, and cell elongation (Curantz and Manceau, 2021). This result also indicates that the speed of follicle emergence in the wake of the competence wave (and thus the emergence of a visible pattern) is constrained by other or additional factors than those providing competence. Cellular changes such as increased cell density and cell elongation (Ho et al., 2019; Bailleul et al., 2019) may indeed not be sufficient for a follicle to emerge, consistent with experimental observation showing that follicle formation necessitates local interactions between skin cells and morphogens of the FGF and BMP families. Testing this hypothesis may be achieved by locally perturbing cell density and/or cell shape anisotropy in a competent region as in (Curantz et al., 2022) and assessing the effect of perturbations on the speed of follicle emergence.

Second, results indicate that cell proliferation is one of the major elements shaping the speed of the morphogenetic wave. Indeed, both when competence is acquired faster than follicle emergence, and *vice versa*, results of numerical simulations align with previous observations that experimentally inhibiting cell proliferation increases the duration of pattern emergence. Much like above, repeating these experiments and assessing the effect of changes in cell proliferation on the respective extents of the competent and follicle-forming regions would allow discriminating whether *in vivo*, the naive skin tissue is in a natural state in which competence acquisition is faster or slower than follicle emergence (i.e., as in the first case, one would observe a larger competent region in front of forming follicles, while in the second case, follicle formation would "lock" with competence progression).

## Acknowledgements

We thank Coline Hermine, Richard Bailleul, and Jonathan Touboul for helpful discussions. MM is a member of PSL Q-Life. This work was funded
by the ANR Grant "PLUME", under grant agreement ANR-21-CE13-0040.

## Appendix A. Linear analysis

We investigated the possibility of pattern emergence for the system (5) with a linear stability analysis. We found that the equilibrium point $(1, u_s, v_s)$ may be stable for the kinetic part of the system (5) but is unstable for the full reaction-diffusion system. The linearized system in the neighborhood of $(1, u_s, v_s)$ is:

$$ w_t = Jw + D^\kappa \nabla^2 w, \quad \text{with } w = \begin{pmatrix} n-1 \\ u-u_s \\ v-v_s \end{pmatrix}, \qquad (A.1) $$

where $J$ and $D^\kappa$ are given by

$$ J := \begin{pmatrix} -\alpha_n & 0 & 0 \\ \delta_u u_s & \gamma_u - \delta_u & -\gamma_v \\ \delta_v v_s & \eta_u & -\delta_v \end{pmatrix}, \quad D^\kappa := \begin{pmatrix} D_n & -\kappa & 0 \\ 0 & D_u & 0 \\ 0 & 0 & D_v \end{pmatrix}. $$

The solution of (A.1) with the form $w = \bar{w}e^{ikx+\lambda t}$ led to the following dispersion relation, which characterizes the relation between the eigenvalue $\lambda$ and the wavenumber $k$

$$ 0 = \lambda^3 + a_2(k^2)\lambda^2 + a_1(k^2)\lambda + a_0(k^2), \qquad (A.2) $$

where

$$ a_2(k^2) = \delta_u - \gamma_u + k^2 D_u + \delta_v + k^2 D_v + \alpha_n + k^2 D_n, $$

$$ a_1(k^2) = (\alpha_n + k^2 D_n) [\delta_u - \gamma_u + k^2 D_u + \delta_v + k^2 D_v] \\ + [\gamma_v \eta_u + (\gamma_u - \delta_u - k^2 D_u)(-\delta_v - k^2 D_v)] - k^2 \kappa \delta_u u_s, $$

and $a_0(k^2) = -\det(A(k^2))$, where $A(k^2) := J - k^2 D$. The Routh-Hurwitz criterion provides three conditions for ensuring that the equilibrium is stable (i.e. eigenvalues have strictly negative real parts), which occurs if and only if the following conditions hold

$$ a_0(k^2) > 0, a_2(k^2) > 0 \quad \text{and} \quad a_2(k^2)a_1(k^2) - a_0(k^2) > 0. \qquad (A.3) $$
Turing instability emerged when the real parts of two eigenvalues re-
mained negative and that of the remaining one passed through zero. With-
out loss of generality, we assumed that $\lambda_1, \lambda_2$ and $\lambda_3$ represent the roots of
the characteristic equation and accordingly obtained

$$
\lambda_1 + \lambda_2 + \lambda_3 = -a_2(k^2), \quad \lambda_1\lambda_2 + \lambda_2\lambda_3 + \lambda_3\lambda_1 = a_1(k^2), \quad \lambda_1\lambda_2\lambda_3 = -a_0(k^2),
$$

$$
\text{and, } -(\lambda_1 + \lambda_2)(\lambda_2 + \lambda_3)(\lambda_3 + \lambda_1) = a_2(k)a_1(k^2) - a_0(k^2).
$$

We assumed there is a Turing instability and $k_c$ denotes the critical wave number $k = k_c$, where, without loss of generality

$$
\lambda_1|_{k^2=k_c^2} = 0, \quad \mathrm{Re} (\lambda_2)|_{k^2=k_c^2} < 0, \quad \mathrm{and} \quad \mathrm{Re} (\lambda_3)|_{k^2=k_c^2} < 0.
$$

Thus, $a_0(k^2) = 0$ at the critical wave number $k = k_c$.

In addition, conditions for Turing instability led to $a_2(k_c^2) > 0$, $a_1(k_c^2) = \lambda_2\lambda_3 > 0$ and $a_2(k_c^2)a_1(k_c^2) - a_0(k_c^2) = a_1(k_c^2)a_2(k_c^2) > 0$. Thus, the equilibrium point $(1, u_s, v_s)$ remained stable when $\det(A(k^2)) < 0$ held for all $k \ge 0$ and became Turing unstable when $\det(A(k^2)) > 0$ held for at least one $k > 0$.

The determinant of the matrix $A(k^2)$ had the following form

$$
\det(A(k^2)) = k^6 A_6 + k^4 A_4(\kappa) + k^2 A_2(\kappa) + A_0,
$$

with

$$
A_6 = -D_n D_u D_v < 0,
$$

$$
A_4(\kappa) = -\alpha_n D_u D_v - D_n D_u \delta_v + D_n D_v (\gamma_u - \delta_u) + \kappa \delta_u u_s D_v,
$$

$$
A_2(\kappa) = \alpha_n (D_v(\gamma_u - \delta_u) - D_u\delta_v) - D_n(\gamma_v\eta_u - \delta_v(\gamma_u - \delta_u)) - \kappa [\gamma_v\delta_v v_s - \delta_u u_s \delta_v],
$$

$$
A_0 = -\alpha_n [\gamma_v \eta_u - \delta_v (\gamma_u - \delta_u)].
$$

The condition for marginal stability at some $k = k_c$ was

$$
\max [\det (A(k_c^2))] = 0.
$$

The minimum of $a_0$ was attained when

$$
k_{\max}^2(\kappa) = \frac{A_4(\kappa) + \sqrt{A_4(\kappa)^2 + 3A_2(\kappa)|A_6|}}{3|A_6|},
$$
<figure><img src="image_11.png" /><figcaption>Figure A.8: Left panel: Plot of $a_0(k)$. Right panel: Growth rate of the kth mode, with a band of growing modes.</figcaption></figure>

which required $A_4(\kappa)^2+3A_2(\kappa)|A_0| > 0$ and $A_4(\kappa)+\sqrt{A_4(\kappa)^2+3A_2(\kappa)|A_0|} > 0$. The graph of $a_0(k)$ depended on $\kappa$, which plays the role of the bifurcation parameter (Figure A.8). Bifurcation could occur at the critical value $\kappa_c$ defined by the root of the function

$$ \kappa \mapsto k_{\max}^6(\kappa)A_6 + k_{\max}^4(\kappa)A_4(\kappa) + k_{\max}^2(\kappa)A_2(\kappa) + A_0, $$

with $k_c$ is given by $k_c^2 := k_{\max}^2(\kappa_c)$. This function does not necessarily admits a unique root, but when it is the case, the root is approximated numerically. When $\kappa < \kappa_c$, the equilibrium point $(1, u_s, v_s)$ was stable for all $k$ (Figure A.8), while for $\kappa > \kappa_c$, the system admitted a pattern-forming stationary instability. The unstable wave numbers stayed comprised between the roots of $a_0(k)$, denoted by $k_1$ and $k_2$ (Figure A.8).

## Appendix B. Weakly non-linear analysis

We explicitly provided the definitions of the operators ($\mathcal{L}^i$) and ($\mathcal{B}^i$) for $i = 0, 1, 2$ introduced in the equations (9) and (10),

$$ \mathcal{L}^0(\mathbf{w}) := (J + D^{\kappa_c} \partial_{xx})\mathbf{w}, \quad \mathcal{L}^1(\mathbf{w}) := 2D^{\kappa_c} \partial_{xX}\mathbf{w} - M_{\kappa_1} \partial_{xx}\mathbf{w}, $$
$$ \mathcal{L}^2(\mathbf{w}) := D^{\kappa_c} \partial_{XX} \mathbf{w} - M_{\kappa_2} \partial_{xx} \mathbf{w} - 2M_{\kappa_1} \partial_{X} \mathbf{w}, \quad \text{with } M_{\kappa_i} := \begin{pmatrix} 0 & \kappa_i & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, $$

$$ \mathcal{B}^0(\mathbf{w}, \mathbf{w}') := \frac{1}{2} \begin{pmatrix} -2\alpha_n n'n \\ \gamma_u(n'u + nu') - \gamma_v(n'v + nv') \\ \eta_u(n'u + un') \end{pmatrix} - \frac{\kappa_c}{2} \begin{pmatrix} \partial_x(n'\partial_x u + n\partial_x u') \\ 0 \\ 0 \end{pmatrix}, $$

and

$$ \mathcal{B}^i(\mathbf{w}, \mathbf{w}') := -\frac{\kappa_i}{2} \begin{pmatrix} \partial_x(n'\partial_x u + n\partial_x u') \\ 0 \\ 0 \end{pmatrix}, \quad \text{for } i = 1, 2. $$

The vector $\rho$ introduced in the equality (14) was defined up to a constant and could be normalized as

$$ \rho = \begin{pmatrix} 1 \\ M_y \\ M_z \end{pmatrix} \in \ker(K - k_c^2 D^{\kappa_c}), \qquad (B.1) $$

with

$$ M_y = \frac{\alpha_n + k_c^2 D_n}{k_c^2 \kappa_c}, \quad M_z = \frac{\delta_u u_s + M_y[\gamma_u - \delta_u - k_c^2 D_u]}{\gamma_v}. $$

Equations (12), and (13) could be rewritten as

$$ O(\epsilon^2): \qquad \mathcal{L}^0(\mathbf{w}_2) = \mathbf{F}, \qquad (B.2) $$

$$ O(\epsilon^3): \qquad \mathcal{L}^0(\mathbf{w}_3) = \mathbf{G}, \qquad (B.3) $$

with

$$ \mathbf{F} := \partial_{T_1} \mathbf{w}_1 - \mathcal{L}^1(\mathbf{w}_1) - \mathcal{B}^0(\mathbf{w}_1, \mathbf{w}_1), \qquad (B.4) $$

and

$$ \mathbf{G} := \partial_{T_2} \mathbf{w}_1 + \partial_{T_1} \mathbf{w}_2 - \mathcal{L}^1(\mathbf{w}_2) - \mathcal{L}^2(\mathbf{w}_1) - \mathcal{B}^1(\mathbf{w}_1, \mathbf{w}_1) - 2\mathcal{B}^0(\mathbf{w}_1, \mathbf{w}_2). \quad (B.5) $$

Substituting (14) into (B.4) led to
$$ 
\begin{aligned}
\mathbf{F} = \cos(k_c x) \rho \left[ \frac{\partial A}{\partial T_1} - k_c^2 A M_{\kappa_1} \right] + 2k_c D^{\kappa_c} \rho \frac{\partial A}{\partial X} \sin(k_c x) \\
- \frac{A^2 \cos(2k_c x)}{2} \begin{pmatrix} 2\kappa_c k_c^2 M_y - \alpha_n \\ \gamma_u M_y - \gamma_v M_z \\ \eta_u M_y \end{pmatrix} - \frac{A^2}{2} \begin{pmatrix} -\alpha_n \\ \gamma_u M_y - \gamma_v M_z \\ \eta_u M_y \end{pmatrix}.
\end{aligned}
$$

By the Fredholm alternative, Equation (B.2) admits a solution, if and only if,

$$ \langle \mathbf{F}, \mathbf{W}^* \rangle = 0, \qquad (B.6) $$

where $\mathbf{W}^* = \psi \cos(k_c x) \in \ker(L^*)$ with $L^*$ the adjoint operator of $\mathcal{L}^0$ and $\langle \cdot, \cdot \rangle$ the scalar product in $L^2(0, 2\pi/k_c)$. Thanks to the Neumann boundary condition, we could construct for (5) a solution in the whole domain using reflection and periodic extension. The vector $\psi$ was defined up to a constant,

$$ 
\begin{aligned}
\psi &:= \begin{pmatrix} 1 \\ \text{with } M_y^* \\ M_z^* \end{pmatrix}, \quad M_y^* = -M_z^* \frac{\delta_v + k_c^2 D_v}{\gamma_v}, \\
M_z^* &= \frac{\gamma_v k_c^2 \kappa_c}{(\delta_v + k_c^2 D_v)(\gamma_u - \delta_u - k_c^2 D_u) - \eta_u \gamma_v}.
\end{aligned}
$$

The condition (B.6) gave

$$ \partial_{T_1} A = \gamma A, \quad \gamma = \frac{k_c^2 \kappa_1 M_y}{1 + M_y M_y^* + M_z M_z^*}. \qquad (B.7) $$

Because the above equation does correctly predict the evolution of amplitude, we chose $T_1 = 0$ and $\kappa_1 = 0$ such that the solvability condition (B.6) holds. Thus, the solution of (B.2) with boundary condition is

$$ \mathbf{w}_2 = A^2 \mathbf{w}_{20} + A^2 \mathbf{w}_{22} \cos(2k_c x) + \frac{\partial A}{\partial X} \mathbf{w}_{21} \sin(k_c x), $$
where vectors $\mathbf{w}_{2i}$, $i = 0, 1, 2$ are the solutions of the following linear systems

$$ egin{cases} J\mathbf{w}_{20} = \frac{1}{2} \begin{pmatrix} -\alpha_n \\ \gamma_u M_y - \gamma_v M_z \\ \eta_u M_y \end{pmatrix}, \\ (J - k_c^2 D^{\kappa_c})\mathbf{w}_{21} = 2k_c D^{\chi_c} \rho, \\ (J - 4k_c^2 D^{\kappa_c})\mathbf{w}_{22} = -\frac{1}{2} \begin{pmatrix} 2\kappa_c k_c^2 M_y - \alpha_n \\ \gamma_u M_y - \gamma_v M_z \\ \eta_u M_y \end{pmatrix}. \end{cases} \quad (B.8) $$

To understand the evolution of the feather follicle amplitude $A$, we ought to consider the third order coefficient equation (B.3). The solvability condition $\langle \mathbf{G}, \mathbf{W}^* \rangle = 0$ led to

$$ \partial_{T_2} A = \nu \partial_{XX} A + \sigma A - LA^3, $$

where parameters $\nu$, $\sigma$ and $L$ were given by

$$ \nu := \frac{\langle 2k_c D^{\kappa_c} \mathbf{w}_{21} + D^{\kappa_c} \rho, \psi \rangle}{\langle \rho, \psi \rangle}, \quad (B.9) $$

$$ \sigma := \frac{\langle \mathbf{G}_{11}, \psi \rangle}{\langle \rho, \psi \rangle}, \quad (B.10) $$

$$ L := \frac{\langle \mathbf{G}_{13}, \psi \rangle}{\langle \rho, \psi \rangle}, \quad (B.11) $$

with

$$ \mathbf{G}_{11} = -k_c^2 M_{\kappa_2} \rho, \text{ and } \mathbf{G}_{13} := -\kappa_c k_c^2 \begin{pmatrix} M_y W_{20}^{(1)} + W_{22}^{(2)} - \frac{M_y}{2} W_{22}^{(1)} \\ 0 \\ 0 \end{pmatrix} - \mathbf{V}_{20} - \mathbf{V}_{22}, $$

where $W_{2j}^{(i)}$ corresponded to the $i$-th component of $\mathbf{w}_{2j}$. The two vectors $\mathbf{V}_{20}$ and $\mathbf{V}_{22}$ were defined by

$$ \mathbf{V}_{20} := \begin{pmatrix} -2\alpha_n W_{20}^{(1)} \\ \gamma_u [W_{20}^{(2)} + M_y W_{20}^{(1)}] - \gamma_v [W_{20}^{(1)} M_z + W_{20}^{(3)}] \\ \eta_u [W_{20}^{(2)} + M_y W_{20}^{(1)}] \end{pmatrix}, $$
$$ \mathbf{V}_{22} := \frac{1}{2} \left( \begin{array}{c} -2\alpha_n W_{22}^{(1)} \\ \gamma_u \left[ W_{22}^{(2)} + M_y W_{22}^{(1)} \right] - \gamma_v \left[ W_{22}^{(1)} M_z + W_{22}^{(3)} \right] \\ \eta_u \left[ W_{22}^{(2)} + M_y W_{22}^{(1)} \right] \end{array} \right). $$

## Appendix C. Numerical simulations

*Numerical scheme:* Numerical simulations were performed using a finite difference scheme. A semi-implicit time discretization was employed, with the diffusion term treated implicitly, while the reaction and chemotaxis terms were discretized explicitly. The corresponding codes are available at the following address: link.

*Parameter values:* The parameter values used in this article are those introduced in Bailleul et al. (2019). In that study, the parameters were adjusted to reproduce the spatiotemporal patterning of follicle formation in embryos. These values are reported in Table C.1. For the simplified model (5), we selected the constants $\gamma_u$, $\gamma_v$, $\eta_u$, $c_u$, and $c_v$ such that the truncated model has the same unique positive equilibrium point $(1, u_s, v_s)$ and the same Jacobian matrix as the original model (1)-(3) (see equation (4)). These values are listed in Table C.2.

<table><caption>Table C.1: Parameter values of the model (1)-(3), introduced in Bailleul et al. (2019), used in our numerical simulations, and their biological interpretation.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Biological interpretation</th></tr></thead><tbody><tr><td>$D_n$</td><td>Variable</td><td>Cell diffusion rate</td></tr><tr><td>$D_u$</td><td>$6 \times 10^{-3}$</td><td>Activator's diffusion rate</td></tr><tr><td>$D_v$</td><td>$0.13$</td><td>Inhibitor diffusion rate</td></tr><tr><td>$\kappa$</td><td>Variable</td><td>Chemotaxis sensitivity</td></tr><tr><td>$\delta_u$</td><td>$10$</td><td>Activator degradation rate</td></tr><tr><td>$\delta_v$</td><td>$20$</td><td>Inhibitor degradation rate</td></tr><tr><td>$\alpha_u$</td><td>$100$</td><td>Activator production rate by cells</td></tr><tr><td>$\omega$</td><td>$20$</td><td>Activator autocatalysis sensitivity</td></tr><tr><td>$\beta_u$</td><td>$2$</td><td>Activator saturation threshold</td></tr><tr><td>$\alpha_v$</td><td>$3500$</td><td>Inhibitor production rate by cells</td></tr><tr><td>$\alpha_n$</td><td>$1 \times 10^{-3}$</td><td>Proliferation rate</td></tr></tbody></table>

<table><caption>Table C.2: Parameter values of the simplified model (5).</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>$D_n$</td><td>Variable</td></tr><tr><td>$D_u$</td><td>$6 \times 10^{-3}$</td></tr><tr><td>$D_v$</td><td>$0.13$</td></tr><tr><td>$\kappa$</td><td>Variable</td></tr><tr><td>$\delta_u$</td><td>$10$</td></tr><tr><td>$\delta_v$</td><td>$20$</td></tr><tr><td>$\alpha_n$</td><td>$1 \times 10^{-3}$</td></tr><tr><td>$\gamma_u$</td><td>$14.54$</td></tr><tr><td>$\gamma_v$</td><td>$0.1$</td></tr><tr><td>$\eta_u$</td><td>$3931$</td></tr><tr><td>$c_u$</td><td>$2.97$</td></tr><tr><td>$c_v$</td><td>$-1104$</td></tr></tbody></table>

*Continuation and bifurcation analysis:* In Figure D.9, we also use the Matlab package `pde2path` to compute the stationary amplitude of the patterns (Uecker et al., 2012; Uecker, 2021). The corresponding codes are available at the following address: link.
Appendix D. Supplementary Figures

<figure><img src="image_12.png" /><figcaption>Figure D.9: Numerical simulations of the model at stationary state validate the weakly non-linear analysis. A-B. Comparison between the approximation of the weakly non-linear analysis (19); dashed red line) and the numerical solution of (5) (black line) for $\epsilon = 0.1$ (A) or $\epsilon = 0.3$ (B). C. Comparison between the bifurcation diagram obtained with pde2path (black line) and the prediction obtained using the weakly non-linear analysis (19); red line). The equilibrium point becomes unstable following a super-critical Turing bifurcation. Results of A ($\epsilon = 0.1$) and B ($\epsilon = 0.3$) are respectively shown on the diagram by a circle and a square. Parameter values: $D_n = 0.1$, $D_u = 6e-3$, $D_v = 0.13$, $\gamma_u = 15.539$, $\gamma_v = 0.01$, $\eta_u = 3931$, $c_u = 2.9672$, $c_v = -1104$ and $\alpha_n = 0.5$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure D.10: **Numerical simulations in the transient regime, close to the Turbing bifurcation ($\epsilon = 0.05$), validate the weakly non-linear analysis.** Comparison between the envelope approximation predicted by the weakly non-linear analysis ((17); blue dashed line) and the numerical solution of (5) (black solid line) at different times. Parameter values are given by: $D_n = 7e-5$, $D_u = 6e-3$, $D_v = 0.13$, $\gamma_u = 14.5389$, $\gamma_v = 0.099932$, $\eta_u = 3931$, $c_u = 2.9672$, $c_v = -1104$, $\alpha_n = 0.001$ and $\epsilon = 0.05$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure D.11: Numerical simulations in the transient regime, further from the Turing bifurcation ($\epsilon = 0.2$), confirm the accuracy of the weakly non-linear analysis on key aspects. Comparison between the envelope approximation from the weakly non-linear analysis ((17); blue dashed line) and the numerical solution of (5) (black solid line) at different times. The weakly non-linear approximation accurately predicts the propagation speed of the envelope and the amplitude of the patterns in its wake. However, at early stages of pattern formation, newly formed follicle structures appear denser in the numerical solution than in the prediction before gradually relaxing to the expected values. Parameter values are the same as in the previous figure D.10 except for $\epsilon = 0.2$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure D.12: **Final predicted amplitude of the patterns, $A_{CD}$, depends on the parameter values.** Each panel shows the effect of varying a single parameter while keeping $\epsilon$ fixed. We observe that $D_u$, $\alpha_n$, $\delta_u$, $\delta_v$, and $\gamma_v$ increase $A_{CD}$, whereas $\gamma_u$ decreases it. The influence of $D_v$ and $\eta_u$ depends on their values, while $D_n$ has a negligible effect. All other parameters are fixed at their values from Fig. D.9.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure D.13: **Predicted speed of the follicle emergence wave in a competent region, $c_{\text{comp→FE}}$, depends on the parameter values.** Each panel shows the effect of varying a single parameter while keeping $\epsilon$ fixed. We observed that increasing $D_n$, $D_u$, $\delta_u$, $\gamma_v$, and $\eta_u$ increases $c_{\text{comp→FE}}$, while increasing $D_v$, $\delta_v$, and $\gamma_u$ decreases it. In contrast, the proliferation rate $\alpha_n$ has only a minor impact on $c_{\text{comp→FE}}$. The parameters that are not varied are fixed at their values from Figure D.9.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure D.14: The predicted steepness of the follicle emergence wave in a competent region $\mu$ depends on parameter values. Each panel shows the effect of varying one parameter while keeping $\epsilon$ fixed. We observed that $D_v$, $\delta_u$, $\delta_v$, $\gamma_u$, $\gamma_v$, and $\eta_u$ influence $\mu$ similarly to their effect on wave speed $c_{comp\to FE}$. However, increasing $D_n$ has the opposite effect, leading to a steeper wave front. The steepness $\mu$ is proportional to the proliferation parameter $\alpha_n$ and varies non-monotonically with $D_u$, showing an initial increase followed by a decrease. The parameters that are not varied are fixed at their values from Figure D.9.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure D.15: The predicted morphogenetic wave speed closely matches the numerical solution over time in both scenarios: $c_{\text{comp}\to\text{FE}} > c_{\text{naive}\to\text{comp}}$ (Top) and $c_{\text{comp}\to\text{FE}} < c_{\text{naive}\to\text{comp}}$ (Bottom). Kymograph of $n$ from the numerical solution of system (5). The dashed red and blue lines correspond to the theoretical speeds $c_{\text{naive}\to\text{comp}}$ and $c_{\text{comp}\to\text{FE}}$, respectively. In the top panel, the speed of follicle emergence locks onto the competence wave, leading to near-instantaneous emergence of patterns upon competence acquisition. Conversely, in the bottom panel, the speed of follicle emergence predicted by the weakly non-linear analysis is slower than the speed of the competence wave. As a result, a growing region of competent but follicle-free tissue emerges over time (between the red and blue lines). Parameter values: $D_n = 1.5e-5$, $D_u = 6e-3$, $D_v = 0.13$, $\gamma_u = 14.5389$, $\gamma_v = 0.099932$, $\eta_u = 3931$, $c_u = 2.9672$, $c_v = -1104$, $\alpha_n = 0.001$, and $\kappa = 1.18e-5$ (Top), $\kappa = 1.14e-5$ (Bottom).</figcaption></figure>
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Figure D.16: **Numerical simulations validate the formula (25) for the speed of the morphogenetic wave, for different values of $D_n$ and $\alpha_n$.** Lines represent the theoretical predictions for the speed of the competence wave $c_{naive\tocomp}$ (in red), the speed of follicle emergence in a competent region $c_{comp\toFE}$ (in blue), and the speed of the morphogenetic wave $c_{morph}$ (in black). Black crosses represent the speed of follicle emergence estimated from numerical simulations of system (5). Colored dots show numerically computed speeds from traveling wave solutions of the Fisher-KPP equation (23) (in red) and the envelope equation (17) (in blue), using the same simulation parameters (i.e., time and space discretization) as black crosses. The close fit between black crosses and the explicit prediction (25) confirms the validity of the formula. Other parameter values are the same as in Figure D.9.</figcaption></figure>

## References

* Aronson, D.G., Weinberger, H.F., 1975. Nonlinear diffusion in population genetics, combustion, and nerve pulse propagation, in: Goldstein, J.A. (Ed.), Partial Differential Equations and Related Topics. Springer Berlin Heidelberg
berg, Berlin, Heidelberg. volume 446, pp. 5–49. doi:10.1007/BFb0070595. series Title: Lecture Notes in Mathematics.

Avery, M., Goh, R., Goodloe, O., Milewski, A., Scheel, A., 2019. Growing Stripes, with and without Wrinkles. SIAM Journal on Applied Dynamical Systems 18, 1078–1117. doi:10.1137/18M1221989.

Bailleul, R., Curantz, C., Desmarquet-Trin Dinh, C., Hidalgo, M., Touboul, J., Manceau, M., 2019. Symmetry breaking in the embryonic skin triggers directional and sequential plumage patterning. PLOS Biology 17, e3000448. doi:10.1371/journal.pbio.3000448.

Baker, R.E., Schnell, S., Maini, P.K., 2009. Waves and patterning in developmental biology: vertebrate segmentation and feather bud formation as case studies. The International Journal of Developmental Biology 53, 783–794. doi:10.1387/ijdb.072493rb.

Barnafi, N.A., De Oliveira Vilaca, L.M., Milinkovitch, M.C., Ruiz-Baier, R., 2022. Coupling Chemotaxis and Growth Poromechanics for the Modelling of Feather Primordia Patterning. Mathematics 10, 4096. doi:10.3390/math10214096.

Ben-Jacob, E., Brand, H., Dee, G., Kramer, L., Langer, J., 1985. Pattern propagation in nonlinear dissipative systems. Physica D: Nonlinear Phenomena 14, 348–364. doi:10.1016/0167-2789(85)90094-6.

Bozzini, B., Gambino, G., Lacitignola, D., Lupo, S., Sammartino, M., Sgura, I., 2015. Weakly nonlinear analysis of Turing patterns in a morphochemical model for metal growth. Computers & Mathematics with Applications 70, 1948–1969. doi:10.1016/j.camwa.2015.08.019.

Collet, P., Eckmann, J.P., 1986. The existence of dendritic fronts. Communications in Mathematical Physics 107, 39–92. doi:10.1007/BF01206953.

Collet, P., Eckmann, J.P., 1990. Instabilities and Fronts in Extended Systems: Princeton University Press. doi:10.1515/9781400861026.

Cooper, R.L., Thiery, A.P., Fletcher, A.G., Delbarre, D.J., Rasch, L.J., Fraser, G.J., 2018. An ancient Turing-like patterning mechanism regulates skin denticle development in sharks. Science Advances 4, eaau5484. doi:10.1126/sciadv.aau5484.
Cross, M.C., Hohenberg, P.C., 1993. Pattern formation outside of equilibrium. Reviews of Modern Physics 65, 851-1112. doi:10.1103/RevModPhys.65.851.

Cruywagen, G.C., Murray, J.D., Maini, P.K., 1997. Biological pattern formation on two-dimensional spatial domains: a nonlinear bifurcation analysis. SIAM Journal on Applied Mathematics 57, 1485-1509. Publisher: SIAM.

Curantz, C., Bailleul, R., Castro-Scherianz, M., Hidalgo, M., Durande, M., Graner, F., Manceau, M., 2022. Cell shape anisotropy contributes to self-organized feather pattern fidelity in birds. PLOS Biology 20, e3001807. doi:10.1371/journal.pbio.3001807.

Curantz, C., Manceau, M., 2021. Trends and variation in vertebrate patterns as outcomes of self-organization. Current Opinion in Genetics & Development 69, 147-153. doi:10.1016/j.gde.2021.05.001.

Dee, G., Langer, J.S., 1983. Propagating Pattern Selection. Physical Review Letters 50, 383-386. doi:10.1103/PhysRevLett.50.383.

Dhouailly, D., 2009. A new scenario for the evolutionary origin of hair, feather, and avian scales. Journal of Anatomy 214, 587-606. doi:10.1111/j.1469-7580.2008.01041.x.

Di Talia, S., Vergassola, M., 2022. Waves in Embryonic Development. Annual Review of Biophysics 51, 327-353. doi:10.1146/annurev-biophys-111521-102500.

Diez, A., Krause, A.L., Maini, P.K., Gaffney, E.A., Seirin-Lee, S., 2024. Turing Pattern Formation in Reaction-Cross-Diffusion Systems with a Bi-layer Geometry. Bulletin of Mathematical Biology 86, 13. doi:10.1007/s11538-023-01237-1.

Eckmann, J.P., Wayne, C.E., 1991. Propagating fronts and the center manifold theorem. Communications in Mathematical Physics 136, 285-307. doi:10.1007/BF02100026.

Economou, A.D., Monk, N.A.M., Green, J.B.A., 2020. Perturbation analysis of a multi-morphogen turing reaction-diffusion stripe patterning system reveals key regulatory interactions. Development, dev.190553doi:10.1242/dev.190553.
Evanitsky, M.N., Di Talia, S., 2023. An active traveling wave of Eda/NF-kappaB signaling controls the timing and hexagonal pattern of skin appendages in zebrafish. Development 150, dev201866. doi:10.1242/dev.201866.

Faye, G., Holzer, M., 2015. Modulated traveling fronts for a nonlocal Fisher-KPP equation: A dynamical systems approach. Journal of Differential Equations 258, 2257-2289. doi:10.1016/j.jde.2014.12.006.

Fisher, R.A., 1937. The wave of advance of advantageous genes. Annals of Eugenics 7, 355-369. doi:10.1111/j.1469-1809.1937.tb02153.x.

Fraga Delfino Kunz, C., Gerisch, A., Glover, J., Headon, D., Painter, K.J., Matthäus, F., 2024. Novel Aspects in Pattern Formation Arise from Coupling Turing Reaction-Diffusion and Chemotaxis. Bulletin of Mathematical Biology 86, 4. doi:10.1007/s11538-023-01225-5.

Gambino, G., Giunta, V., Lombardo, M.C., Rubino, G., 2022. Cross-diffusion effects on stationary pattern formation in the FitzHugh-Nagumo model. Discrete and Continuous Dynamical Systems - B 27, 7783. doi:10.3934/dcdsb.2022063.

Gambino, G., Lombardo, M., Sammartino, M., 2012. Turing instability and traveling fronts for a nonlinear reaction-diffusion system with cross-diffusion. Mathematics and Computers in Simulation 82, 1112-1132. doi:10.1016/j.matcom.2011.11.004.

Gambino, G., Lombardo, M., Sammartino, M., 2013a. Pattern formation driven by cross-diffusion in a 2D domain. Nonlinear Analysis: Real World Applications 14, 1755-1779. doi:10.1016/j.nonrwa.2012.11.009.

Gambino, G., Lombardo, M.C., Sammartino, M., Sciacca, V., 2013b. Turbing pattern formation in the Brusselator system with nonlinear diffusion. Physical Review E 88, 042925. doi:10.1103/PhysRevE.88.042925.

Giunta, V., Lombardo, M.C., Sammartino, M., 2021. Pattern Formation and Transition to Chaos in a Chemotaxis Model of Acute Inflammation. SIAM Journal on Applied Dynamical Systems 20, 1844-1881. doi:10.1137/20M1358104.
Glover, J.D., Wells, K.L., Matthäus, F., Painter, K.J., Ho, W., Riddell, J., Johansson, J.A., Ford, M.J., Jahoda, C.A.B., Klika, V., Mort, R.L., Headon, D.J., 2017. Hierarchical patterning modes orchestrate hair follicle morphogenesis. PLOS Biology 15, e2002117. doi:10.1371/journal.pbio.2002117.

Goh, R., Scheel, A., 2014. Triggered fronts in the complex Ginzburg Landau equation. Journal of Nonlinear Science 24, 117–144. Publisher: Springer.

Goh, R., Scheel, A., 2016. Pattern formation in the wake of triggered pushed fronts. Nonlinearity 29, 2196. Publisher: IOP Publishing.

Goh, R., Scheel, A., 2018. Pattern-forming fronts in a Swift–Hohenberg equation with directional quenching — parallel and oblique stripes. Journal of the London Mathematical Society 98, 104–128. doi:10.1112/jlms.12122.

Goh, R., Scheel, A., 2023. Growing patterns. Nonlinearity 36, R1–R51. doi:10.1088/1361-6544/acf265.

Han, R., Dai, B., 2017. Cross-diffusion induced Turing instability and amplitude equation for a toxic-phytoplankton-zooplankton model with non-monotonic functional response. International Journal of Bifurcation and Chaos 27, 1750088. Publisher: World Scientific.

Han, Y., Li, Z., Tao, J., Ma, M., 2017. Pattern formation for a volume-filling chemotaxis model with logistic growth. Journal of Mathematical Analysis and Applications 448, 885–907. Publisher: Elsevier.

Han, Y., Li, Z., Zhang, S., Ma, M., 2016. Wavefront invasion for a volume-filling chemotaxis model with logistic growth. Computers & Mathematics with Applications 71, 471–478. Publisher: Elsevier.

Hillen, T., Painter, K.J., 2009. A user's guide to PDE models for chemotaxis. Journal of Mathematical Biology 58, 183–217. doi:10.1007/s00285-008-0201-3.

Ho, W.K.W., Freem, L., Zhao, D., Painter, K.J., Woolley, T.E., Gaffney, E.A., McGrew, M.J., Tzika, A., Milinkovitch, M.C., Schneider, P., Drusko, A., Matthäus, F., Glover, J.D., Wells, K.L., Johansson, J.A., Davey, M.G.,
Sang, H.M., Clinton, M., Headon, D.J., 2019. Feather arrays are patterned by interacting signalling and cell density waves. PLOS Biology 17, e3000132. doi:10.1371/journal.pbio.3000132.

Holzer, M., Scheel, A., 2014. Criteria for Pointwise Growth and Their Role in Invasion Processes. Journal of Nonlinear Science 24, 661-709. doi:10.1007/s00332-014-9202-0.

Hornik, C., Krishan, K., Yusuf, F., Scaal, M., Brand-Saberi, B., 2005. cDermo-1 misexpression induces dense dermis, feathers, and scales. Developmental Biology 277, 42-50. doi:10.1016/j.ydbio.2004.08.050.

Jahanbakhsh, E., Milinkovitch, M.C., 2022. Modeling convergent scale-by-scale skin color patterning in multiple species of lizards. Current Biology 32, 5069-5082.e13. doi:10.1016/j.cub.2022.10.044.

Jiang, T.X., Jung, H.S., Widelitz, R.B., Chuong, C.M., 1999. Self-organization of periodic patterns by dissociated feather mesenchymal cells and the regulation of size, number and spacing of primordia. Development 126, 4997-5009. doi:10.1242/dev.126.22.4997.

Jiang, T.X., Widelitz, R.B., Shen, W.M., Will, P., Wu, D.Y., Lin, C.M., Jung, H.S., Chuong, C.M., 2004. Integument pattern formation involves genetic and epigenetic controls: feather arrays simulated by digital hormone models. The International Journal of Developmental Biology 48, 117-135. doi:10.1387/ijdb.15272377.

Klika, V., Gaffney, E.A., Maini, P.K., 2024. On the speed of propagation in Turing patterns for reaction-diffusion systems. ArXiv:2403.09247 [nlin].

Kolmogorov, A.N., 1937. A study of the equation of diffusion with increase in the quantity of matter, and its application to a biological problem. Moscow University Bulletin of Mathematics 1, 1-25.

Ma, M., Gao, M., Carretero-González, R., 2019. Pattern formation for a two-dimensional reaction-diffusion model with chemotaxis. Journal of Mathematical Analysis and Applications 475, 1883-1909. doi:10.1016/j.jmaa.2019.03.060.
Ma, M., Gao, M., Tong, C., Han, Y., 2016. Chemotaxis-driven pattern formation for a reaction-diffusion-chemotaxis model with volume-filling effect. *Computers & Mathematics with Applications* 72, 1320-1340. doi:10.1016/j.camwa.2016.06.039.

Manukyan, L., Montandon, S.A., Fofonjka, A., Smirnov, S., Milinkovitch, M.C., 2017. A living mesoscopic cellular automaton made of skin scales. *Nature* 544, 173-179. doi:10.1038/nature22031.

Michon, F., Forest, L., Collomb, E., Demongeot, J., Dhouailly, D., 2008. BMP2 and BMP7 play antagonistic roles in feather induction. *Development* 135, 2797-2805. doi:10.1242/dev.018341.

Milinkovitch, M.C., Jahanbakhsh, E., Zakany, S., 2023. The Unreasonable Effectiveness of Reaction Diffusion in Vertebrate Skin Color Patterning. *Annual Review of Cell and Developmental Biology* 39, 145-174. doi:10.1146/annurev-cellbio-120319-024414.

Myerscough, M.R., Murray, J.D., 1992. Analysis of propagating pattern in a chemotaxis system. *Bulletin of Mathematical Biology* 54, 77-94. doi:10.1007/BF02458621.

Nakamasu, A., Takahashi, G., Kanbe, A., Kondo, S., 2009. Interactions between zebrafish pigment cells responsible for the generation of Turing patterns. *Proceedings of the National Academy of Sciences* 106, 8429-8434. doi:10.1073/pnas.0808622106.

Oster, G.F., Murray, J.D., Harris, A.K., 1983. Mechanical aspects of mesenchymal morphogenesis. *Development* 78, 83-125. doi:10.1242/dev.78.1.83.

Painter, K.J., Ho, W., Headon, D.J., 2018. A chemotaxis model of feather primordia pattern formation during avian development. *Journal of Theoretical Biology* 437, 225-238. doi:10.1016/j.jtbi.2017.10.026.

Palmquist, K.H., Tiemann, S.F., Ezzeddine, F.L., Yang, S., Pfeifer, C.R., Erzberger, A., Rodrigues, A.R., Shyer, A.E., 2022. Reciprocal cell-ECM dynamics generate supracellular fluidity underlying spontaneous follicle patterning. *Cell* 185, 1960-1973.e11. doi:10.1016/j.cell.2022.04.023.
Sadier, A., Anthwal, N., Krause, A.L., Dessalles, R., Lake, M., Bentolila, L.A., Haase, R., Nieves, N.A., Santana, S.E., Sears, K.E., 2023. Bat teeth illuminate the diversification of mammalian tooth classes. Nature Communications 14, 4687. doi:10.1038/s41467-023-40158-4.

Sadier, A., Twarogowska, M., Steklikova, K., Hayden, L., Lambert, A., Schneider, P., Laudet, V., Hovorakova, M., Calvez, V., Pantalacci, S., 2019. Modeling Edar expression reveals the hidden dynamics of tooth signaling center patterning. PLOS Biology 17, e3000064. doi:10.1371/journal.pbio.3000064.

Sengel, P., 1976. Morphogenesis of skin. Number 3 in Developmental and cell biology series, Cambridge Univ. Press, Cambridge.

Shyer, A.E., Rodrigues, A.R., Schroeder, G.G., Kassianidou, E., Kumar, S., Harland, R.M., 2017. Emergent cellular self-organization and mechanosensation initiate follicle pattern in the avian skin. Science 357, 811-815. doi:10.1126/science.aai7868.

Sudderick, Z.R., Glover, J.D., 2024. Periodic pattern formation during embryonic development. Biochemical Society Transactions 52, 75-88. doi:10.1042/BST20230197.

Tseng, C.C., Woolley, T.E., Jiang, T.X., Wu, P., Maini, P.K., Widelitz, R.B., Chuong, C.M., 2024. Gap junctions in Turing-type periodic feather pattern formation. PLOS Biology 22, e3002636. doi:10.1371/journal.pbio.3002636.

Tulumello, E., Lombardo, M.C., Sammartino, M., 2014. Cross-Diffusion Driven Instability in a Predator-Prey System with Cross-Diffusion. Acta Applicandae Mathematicae 132, 621-633. doi:10.1007/s10440-014-9935-7.

Tzika, A.C., Ullate-Agote, A., Zakany, S., Kummrow, M., Milinkovitch, M.C., 2023. Somitic positional information guides self-organized patterning of snake scales. Science Advances 9, eadf8834. doi:10.1126/sciadv.adf8834.

Uecker, H., 2019. Pattern formation with pde2path – a tutorial. doi:10.48550/ARXIV.1908.05211. version Number: 3.
Uecker, H., 2021. Numerical continuation and bifurcation in nonlinear PDEs. Number 174 in Other Titles, Society for Industrial and Applied Mathematics, Philadelphia.

Uecker, H., Wetzel, D., Rademacher, J.D.M., 2012. pde2path - A Matlab package for continuation and bifurcation in 2D elliptic systems. ArXiv:1208.3112 [math].

Van Saarloos, W., 2003. Front propagation into unstable states. Physics Reports 386, 29-222. doi:10.1016/j.physrep.2003.08.001.

Van Saarloos, W., Hohenberg, P., 1992. Fronts, pulses, sources and sinks in generalized complex Ginzburg-Landau equations. Physica D: Nonlinear Phenomena 56, 303-367. doi:10.1016/0167-2789(92)90175-M.

Wollkind, D.J., Manoranjan, V.S., Zhang, L., 1994. Weakly Nonlinear Stability Analyses of Prototype Reaction-Diffusion Model Equations. SIAM Review 36, 176-214. doi:10.1137/1036052.