## Abstract

We study the formation of auxin peaks in a generic class of concentration-based auxin transport models, posed on static plant tissues. Using standard asymptotic analysis, we prove that, on bounded domains, auxin peaks are not formed via a Turing instability in the active transport parameter, but via simple corrections to the homogeneous steady state. When the active transport is small, the geometry of the tissue encodes the peaks’ amplitude and location: peaks arise where cells have fewer neighbours, that is, at the boundary of the domain. We test our theory and perform numerical bifurcation analysis on two models that are known to generate auxin patterns for biologically plausible parameter values. In the same parameter regimes, we find that realistic tissues are capable of generating a multitude of stationary patterns, with a variable number of auxin peaks, that can be selected by different initial conditions or by quasi-static changes in the active transport parameter. The competition between active transport and production rate determines whether peaks remain localized or cover the entire domain. In particular, changes in the auxin production that are fast with respect to the cellular life cycle affect the auxin peak distribution, switching from localized spots to fully patterned states. We relate the occurrence of localized patterns to a *snaking* bifurcation structure, which is known to arise in a wide variety of nonlinear media, but has not yet been reported in plant models.

## 1. Introduction

The hormone auxin plays a crucial role in plant development [1–4], yet the mechanisms through which it accumulates in certain cells and interacts with cell growth mechanisms remain largely unclear. The patterns formed during the growth of a plant are controlled by the local auxin concentration. For example, it is known that the distribution of auxin maxima in the shoot apex gives rise to the formation of primordia [5–10]. Similarly, the distribution of auxin in the root tip coordinates cell division and cell expansion [11,12]. In models of root hair initiation, intracellular levels and gradients of auxin concentration influence the localization of G proteins, which in turn promote hair formation [13]. In addition, it is known that the distribution of auxin in the leaf primordia mediates vascular patterning [1]. In recent years, many aspects of the molecular basis of these mechanisms have been unravelled, and mathematical models of auxin transport have been proposed to explain growth and development [14–17].

Computer simulations are often used to compare the model output with observed data such as auxin distribution, venation patterns, growth or development. At cellular level, carriers such as PIN-FORMED (PIN) proteins, which are localized in the cell membrane, determine the rate and direction of auxin transport. The coordinated activity of many cells can create peaks of auxin that drive differentiation and growth. Various models that implement and refine these ideas have been proposed [5–7,11,18–20]. Such models differ primarily in the specifics of the transport and the coupling to the cell growth and division, but a common feature is that they generate spatially extended patterns of auxin concentration, which have also been observed experimentally.

Existing transport models can be classified into two main categories, *flux-based* and *concentration-based*, depending on how auxin influences the localization of transport mediators (PINs) to form patterns. In flux-based models, first proposed in reference [21], the polarization depends on the net auxin flux between neighbouring cells: the higher the net flux towards the neighbours, the more PIN will accumulate at the membrane, and changes in the PIN distribution determine changes in auxin fluxes. By contrast, in concentration-based models, it is assumed that the PIN accumulation on the membrane is caused by differences in auxin concentration between neighbouring cells. This type of model was introduced in references [5,7]. For other reviews on flux- and concentration-based models, we refer the reader to references [22–25].

In the models mentioned above, patterns are found by direct numerical simulation, upon choosing control parameters within a plausible biological range. However, there is uncertainty regarding many of the parameter values which are often approximated [26,27], adopted from different systems [28] or estimated with large error margins [29]. Furthermore, the effect of systematic parameter variations on the generated patterns and how this relates to the behaviour of the biological system is unclear: understanding the formation of auxin peaks from a dynamical system standpoint is still an open problem; therefore, a mathematical exploration of the parameter variations may generate new, experimentally testable hypotheses, thereby gaining insights into pattern formation mechanisms [1,9,30].

In spite of the uncertainty regarding experimental parameters, it is believed that active transport is a key player in auxin patterning [27]. Transport models posses an inherent timescale separation: the growth hormone dynamics involve short timescales (of the order of seconds) [31], whereas changes in cellular shapes and proliferation of new cells occur on much slower timescales (hours or days) [32]. In order to determine the distribution of auxin in the plant, it is then possible to concentrate on the fast timescale of the hormone transport, assume a static cell structure and study the plant tissue as a dynamical system, subject to variations in the active transport parameter.

In this paper, we perform such an exploration on concentration-based auxin models, which are studied using standard bifurcation analysis techniques [33]. In particular, we find steady states of the system and explore their dependence upon control parameters, investigating how patterns lose or gain stability in response to parameter changes. The aim is to predict qualitatively the distribution patterns that can occur for a certain parameter range and to understand transitions between different types of patterns.

At present, only a few papers regard transport models as dynamical systems [34–36]: among them, reference [36] stands out for being a *systematic* analysis of flux- and concentration-based models, whose auxin patterns are studied by considering local interactions between cells. In this paper, we take the analysis one step further by finding steady states simultaneously in the whole tissue and by studying the important effects of its boundedness on the auxin patterns.

The main result of our analysis is that, in a large class of concentration-based models posed on finite tissues, peaks do not arise from an instability of the homogeneous steady state, as it was previously reported for unbounded tissues [5,7,37,38]: on regular bounded domains, the geometry of the tissue drives the formation of small auxin peaks, which nucleate without instabilities near the boundary.

Further, we investigate the effects of changes in the active transport coefficient, in the diffusivity coefficient and in the auxin production coefficient for two specific examples: the concentration-based model proposed by Smith *et al.* [5] and a more recent modification studied by Chitwood *et al.* [39]. In these systems, the localized peaks, determined analytically for small values of the active transport coefficient, also persist for moderate and large values of this control parameter. We found that, owing to their boundedness, realistic tissues can select from a multitude of patterned configurations, characterized by a variable number of localized peaks and organized in a characteristic snaking bifurcation diagram.

Snaking bifurcation diagrams are commonly found for localized states arising in (systems of) nonlinear partial differential equations posed in one [40–45], two [46,47] and three [48,49] spatial dimensions, as well as in discrete [50] and non-local systems [51–54]. However, this mechanism has never been reported for auxin models: solutions with localized peaks undergo a series of saddle-node bifurcations, giving rise to a hierarchy of steady states with an increasing number of bumps. A direct consequence of this mechanism is that the resulting patterns are robust against changes in the transport parameter and other control parameters as found, for instance, by Sahlin *et al.* [37]. We argue that this mechanism could be a robust feature in several other types of concentration-based auxin models.

The paper is structured as follows: in §2, we summarize our working hypotheses and describe the main results of the paper; in §3.1, we present our calculations for a simple one-dimensional tissue, giving a primer on bifurcation analysis for auxin models; in §3.2, we generalize our results to two-dimensional tissues, which are further discussed in §4; finally, we provide a more formal presentation of our general asymptotic results in §5.

## 2. Mathematical formulation and summary of the main results

We begin by giving a generic definition of concentration-based models, and a summary of the main results of the paper. In concentration-based models, cells are identified with an index *i∈*{1, *…* , *n*} and each cell is associated with a set of neighbours , containing neighbours, and a vector of *m* time-dependent state variables *y** _{i}*(

*t*). For instance,

*y**may contain the auxin concentration (*

_{i}*m*= 1) or both auxin and PIN-FORMED1 concentrations (

*m*= 2). Generically, the rate of change of the concentrations is expressed as a balance between production and decay within the cell, diffusion towards neighbouring cells and active transport, hence, we study generic concentration-based auxin models of the form 2.1where

**and**

*π***are the production and decay functions, respectively,**

*δ***is a diffusion matrix,**

*D**T*is the active transport parameter and

*ν**are the active transport functions. In this paper, we concentrate on the fast timescale of hormone transport and hence consider plant organs as static cell structures, so the number of cells*

_{ij}*n*remains constant in time. We make two key assumptions:

### Hypothesis 2.1 (Regular domains)

Cells are arranged in a regular domain, that is, they have all the same shape and size, and they tessellate the tissue. We do not make any assumption on the dimensionality of the domain.

### Hypothesis 2.2 (Active transport functions)

The active transport functions can be expressed as
2.2where the functions *ψ*_{l}, *φ*_{l} depend on the model choices.

Hypothesis 2.2 is a factorization of the active transport functions that is met by several models in the literature [5,7,29,39,55]: active transport between cell *i* and *j* is influenced not only by the respective concentrations *y** _{i}* and

*y**, but also by concentrations in the neighbours of cell*

_{j}*i*. In §2.1, we introduce examples of specific models satisfying hypothesis 2.2 and in §5, we derive explicit expressions for the corresponding functions

*ψ*,

_{l}*φ*.

_{l}It is established in the literature that concentration-based models are capable of reproducing auxin patterns that are found in SAM experiments [5–7,37,39,56], for biologically realistic values of the transport parameters. In this paper, we address the following questions: what type of patterns are generated by the class of models described above? Do they all predict the occurrence of auxin peaks? What are the instabilities that lead to the formation of auxin peaks? Are auxin patterns robust to changes in the control parameters and initial conditions?

These questions have been partially addressed in previous papers [5–7,35–37,39], where analytical results have been obtained only for particular models and only for regular domains without boundaries, where all cells have the same number of neighbours. In such domains, a homogeneous steady state *y** _{i}* =

***, satisfying the balancing condition**

*y***(**

*π*

*y**) =*

_{i}**(**

*δ*

*y**), is known to exist for all values of*

_{i}*T*[5–7,35,37–39]. In domains without boundaries, the formation of peaks has been explained in terms of a Turing bifurcation in the active transport parameter (figure 1): peaks are formed all at once as

*T*is varied, and they derive from an instability of the homogeneous steady state. From a dynamical system viewpoint, however, we expect that boundary conditions and finite domain sizes affect the formation and selection of auxin patterns. The main result of our investigation is that the mechanism for the formation of peaks is radically different in tissues of finite size, and in particular:

### Result 2.3 (Homogeneous steady states)

*In finite domains, concentration-based models* (*2.1*) *satisfying hypotheses 2.1–2.2 support a homogeneous equilibrium y_{i}* =

*y*** for T*=

*0 but, generically, this homogeneous state is not present when T*≠

*0. This result is a direct consequence of the geometry of the domain: an inspection of equation*(

*2.2*)

*shows that the active transport terms*

*ν*_{ij}contain a nested sum over the neighbours of the cell i and, in a finite domain, the number of neighbours varies from cell to cell, namely cells at the boundary of the regular tissue have fewer neighbours than cells in the interior; consequently,*ν*_{ij}−

*ν**=*

_{ji}is generally different from**0**at the boundary, even when**y**_{i}

*y***. The conclusion is that, on finite domains, peaks cannot form with a Turing bifurcation, because the homogeneous steady state exists for T*=

*0, but not for T*≠

*0 (however small). For regular domains where the number of neighbours is the same for all cells, the Turing mechanism is still possible. In addition, a Turing bifurcation is also possible in any domain, provided that T*=

*0 and diffusion is used as a bifurcation parameter. In this regime, however, the active transport is absent, and the tissue is a standard medium with reaction–diffusion mechanisms, which is not a biologically valid hypothesis for auxin models*.

### Result 2.4 (Origin of auxin peaks)

*To understand how peaks are formed from the homogeneous state, we study the case of small active transport coefficients. For* *and in the absence of passive transport,* *D* = *0**, the models above generate steady states with small peaks. Such states take the form*
2.3*where the coefficients ξ_{i} depend on the number of neighbours at distance 2 from the i-th cell^{[1]}, namely*

*Equation*(

*2.3*)

*predicts that peaks are formed as small perturbations of the homogeneous steady state*(

*figure 1*).

*The amplitude of small peaks is proportional to T and*=

*ξ*_{i}. Importantly,*ξ*_{i}*0 in the interior of regular domains, therefore, peaks localize at the boundary and without bifurcations, as opposed to the Turing scenario where they form everywhere owing to an instability of*

**y***: the mechanism on finite domains is purely geometric, as the location of the peaks is determined by the factors*ξ*_{i}.### Result 2.5 (Effect of passive transport)

*For small active transport and at the presence of passive transport, D* ≠

**0**,

*we still obtain solutions with localized peaks, similar to the above-discussed case. The location of the peaks depends again on ξ*

_{1},*…*,

*ξ*.

_{n}### Remark 2.6 (Applicability of analytical results)

Results 2.3–2.5 are valid for generic models of the form (2.1), provided they satisfy hypotheses 2.1–2.2. In particular, these results are valid for regular cellular arrays in any spatial dimension. The main implication of this result is that a wide class of concentration-based models are able to generate spontaneously auxin peaks in various geometries, irrespective of the model specifics. A similar derivation can be done for irregular domains, albeit the coefficients *ξ*_{i} depend in this case on the local cellular volumes as well as the number of neighbours.

The above-described analytical theory, which is presented in more detail in §5, is valid only for small values of the active transport coefficient and does not explain the formation of peaks in the interior of the domain [6,37,57]. While it is difficult to make general analytical predictions for larger values of *T*, it is possible to explore the solution landscape of specific models via numerical methods. We investigated regular domains in two models by Smith *et al.* [5] and Chitwood *et al.* [39] (henceforth called the Smith model and the Chitwood model, respectively) which satisfy hypotheses 2.1–2.2 and are described in detail in §2.1. Our computations confirm the theoretical results 2.3–2.5 and provide numerical evidence for the following conclusions:

### Result 2.7 (Formation of stable auxin spots in the interior)

*As active transport increases, the two models by Smith and Chitwood predict the formation of peaks in the interior. Peaks are formed progressively, from the boundary towards the interior via saddle-node bifurcations, with a characteristic snaking bifurcation diagram. From a biological perspective, this means that the tissue can form peaks that are robust with respect to changes in the control parameters. In addition, the tissue is capable of selecting from a variety of auxin patterns, with a variable number of spots, depending on the initial conditions and on the value of the auxin transport parameter T* (*figure 1*).

### Result 2.8 (Robustness of the snaking mechanism)

*The scenario above is robust to perturbations in secondary parameters, that is, snaking is found for a range of D values, even when the parameters in the production function π are changed. This suggests that solutions with localized peaks at the boundary should be observable in experiments for which the two models are applicable, if the active transport is inhibited with respect to passive transport.*

Furthermore, in numerical calculations, we are able to violate hypotheses 2.1 and 2.2 and see how this affects our results. An important conclusion is the following:

### Result 2.9 (Irregular domains)

*When the Smith and the Chitwood models are posed on irregular domains, the bifurcation structure presented above persists. Auxin peaks are formed progressively via saddle-node bifurcations, albeit they may, in principle, nucleate in the interior of the domain before the boundary is filled. Furthermore, this scenario is robust to changes in secondary parameters, so result 2.8 is still valid on irregular domains for both models*.

### 2.1. Two models of auxin transport

We now present two models that will be analysed in detail using the asymptotic theory of §5 and numerical simulations. As a first example, we consider the Smith model [5], which features two state variables per cell, namely the indole-3-acetic acid (IAA) concentration, *a _{i}*(

*t*), and the PIN-FORMED1 (PIN1) amount,

*p*(

_{i}*t*). The model features IAA production, decay, active and passive transport terms, whereas for PIN1, only production and decay are included. This results in the following set of coupled nonlinear ODEs 2.4 2.5for

*i*= 1,

*…*,

*n*. In this model,

*D*is a diffusion coefficient,

*V*is the cellular volume,

_{i}*l*=

_{ij}*S*/(

_{ij}*W*+

_{i}*W*) is the ratio between the contact area

_{j}*S*of the adjacent cells

_{ij}*i*and

*j*, and the sum of the corresponding cellular wall thicknesses

*W*and

_{i}*W*. In addition,

_{j}*T*is the active transport coefficient and

*P*is the number of PIN1 proteins on the cellular membrane of cell

_{ij}*i*facing cell

*j*, 2.6The Smith model posed on a regular domain satisfies hypotheses 2.1 and 2.2 in §2, and the reader can find explicit expressions for the functions

*φ*,

_{l}*ψ*in §5.1.2. More details on the model and simulations of realistic phyllotactic patterns can be found in [5].

_{l}The second concentration-based transport model that will be studied below is the more recent Chitwood model [39]. This modification of the Smith model is able to produce stable spiral phyllotactic patterns once cell division is included. The system also features two variables per cell, the IAA concentration and the PIN1 amount, and it is given by the following set of coupled nonlinear ODEs
2.7and
2.8for *i* = 1, *…* , *n*, where *P _{ij}* are given by (2.6), and the only new parameter,

*c*

_{2}, controls the exponential transport. The Chitwood model posed on a regular domain also satisfies hypotheses 2.1 and 2.2, as shown in §5.1.2.

In the reminder of the paper, we shall fix most parameters in the Smith and the Chitwood models, and we will examine variations in the active transport parameter *T*, auxin diffusion coefficient *D* and auxin production coefficient *ρ*_{IAA}. In table 1, we report a brief description of parameters for both models, together with characteristic values and units, which are taken from [5,39].

### Remark 2.10 (Comparison with experimental parameters)

The model equations (2.4)–(2.8) show that the transport parameters *T* and *D* are scaled by the cellular volumes *V*_{i}: comparisons with experimental parameters and other computer simulations available in the literature should be based on the ratios *T/V*_{i} and *D/V*_{i}. Through these ratios, we implicitly specify *T/D*, so as to account for competition between active and passive transports.

### Remark 2.11 (Tissue types)

We model three plant tissues: identical cubic cells arranged on a line of finite length (here and henceforth, one-dimensional regular), identical hexagonal prismic cells tessellating a finite square (two-dimensional regular) and irregular prismic cells tessellating an almost-circular domain (two-dimensional irregular, taken from Merks *et al.* [58]). We stress that the nomenclature one-dimensional and two-dimensional refers to the domain, not the cells, which are assumed to have consistently assigned volumes and contact areas. We note that cellular volumes and contact areas may change between different domains (see also remark (2.10)).

## 3. Results

### 3.1. A primer on the formation of auxin peaks—a one-dimensional regular tissue

Here, we present in detail results 2.3–2.8 for the Smith model posed on a one-dimensional regular tissue, which is represented in figure 2. The tissue consists of a file of *n* identical cubic cells with volume *V* = 1 µm^{3} and *l _{ij}* = 1 µm. We assume that there exist cells to the left of cell 1 and that the net proximal flux is zero; hence, we prescribe Neumann boundary conditions at

*i*= 1. We impose that there are no cells to the right of cell

*n*, so we use a free boundary condition

^{2}at

*i*=

*n*by setting . In this geometry, the tissue has a physical boundary only at

*i*=

*n*, as illustrated in figure 2, so each cell has two neighbours, except cell

*n*, which has only one neighbour. This geometry (but with other boundary conditions) has been studied in other auxin-patterning simulations: it was used in reference [35] as an approximation of a root tissue and in reference [59] to model part of a leaf, between the midvein and the margin. In this paper, we use this geometry only as a primer to illustrate how the asymptotic and numerical calculations are used to predict the formation of auxin spots in concentration-based models. While the focus is on the Smith model, several of the results presented in this section are valid in more general models and spatial configurations (we refer the reader to remark 2.6 and the whole §2 for further comments on the generality of these results).

#### 3.1.1. From homogeneous to patterned solutions

The Smith model posed on a regular domain satisfies hypotheses 2.1 and 2.2 (as shown in §5.1.2), so we can apply the asymptotic theory in §5 (in particular, lemma 5.3). On an unbounded array (or on a bounded array with periodic boundary conditions), every cell has the same number of neighbours; therefore, the model admits the following homogeneous steady state
3.1for *i* = 1, *…* *n*. However, the one-dimensional regular domain considered here is finite, so the homogeneous solution exists only for *T* = 0 μm^{3} h^{−1}, as for positive *T*, the sums in the transport term in (2.4) do not vanish in general (result 2.3). This is in accordance with experimental data presented in reference [5] (figure 2), where auxin peaks were suppressed by inhibiting auxin transport on a tissue of finite size. For *D* = 0 µm^{2} h^{−1} and , lemma 5.3 shows the existence of patterns in the form of small deviations from the homogeneous steady state (see §5.1.2 for a full derivation).

Equations (3.2) predict that, for infinitesimal values of the transport coefficient *T*, the perturbed solution coincides with the homogeneous solution, except for a small peak at cell *n −* 1 and a small dip at cell *n* (see also bold curve in figure 2). In other words, peaks are present where the number of neighbours differs from the number of neighbours in the unbounded domain, that is, where the sum in the active transport in equation (2.4) is non-zero (result 2.4). The above-mentioned analysis applies also in the limit of slow dynamics of the PIN1 proteins: upon assuming *p* constant and homogeneous in the tissue, we find that *a _{i}* is as in equation (3.2), with

*p** replaced by

*p*(see §5.1.1 for details).

We have thus established that, in regular domains, a small auxin transport coefficient *T* elicits low auxin peaks. This correlation was previously reported in numerical experiments on various models [7,37], and we now provide a mathematical explanation of this phenomenon.

Asymptotic calculations can also be carried out in the presence of diffusion (*D* ≠ 0), leading to a linear system for the perturbations. In §5.2, we present a derivation for generic models in generic tissues, which is then specialized for the Smith model as an example. Quantitative results of this calculation are shown in figure 2, where we plot approximate steady states for the Smith model towards the boundary *i* = *n* for *T* = 3 × 10^{−}^{5} μm^{3} h^{−1} and various values of the diffusion coefficient. We note that, in the regime of small active transport and comparatively much bigger diffusion coefficient, a peak is still present at the boundary. Inspecting the solid line (*D* = 0 μm^{2} h^{−1}) and the dashed lines (*D* = 0.06 μm^{2} h^{−1} and *D* = 0.18 μm^{2} h^{−1}), we see that the peaks decrease in amplitude and are more spread out, as expected (result 2.5). As a concluding remark, we point out that these results are not influenced by the no-flux boundary conditions specified at *i* = 1: the only physical boundary is at *i* = *n*, where the number of neighbours differs from the interior.

#### 3.1.2. Forming high peaks in the interior via snaking

We now turn to the more interesting question of how the tissue develops high auxin peaks (result 2.6) which are observed in experiments. Once again, we illustrate our findings in the one-dimensional regular case and generalize in the following sections.

In realistic simulations, the transport coefficient *T* is not necessarily small [5,7,27]; therefore, it is interesting to explore the solution landscape when *T* is increased at the presence of diffusion. This is done using numerical bifurcation analysis, that is, equilibria of systems (2.1) are followed in parameter space using the Newton–Raphson method and pseudo-arclength continuation [33]. Linear stability is then inferred computing the spectrum of the Jacobian at the steady state.

In figure 3, we show a branch of solutions of the Smith model for the one-dimensional regular domain obtained with the parameter set in table 1. We start the computation from the homogeneous solution at *T* = 0 and follow the pattern for increasing values of *T*. As *T* changes, we plot the 2-norm of the auxin vector, **|| a||**

_{2}, which is a measure of the spatial extent of the solution (the lower

**||**

*a*||_{2}, the more localized the pattern) and denote stable (unstable) branches with solid (dashed) or thick (thin) lines.

The low peak found close to the boundary persists for increasing values of *T* and grows steadily until we meet a first turning point (TP_{1}). Before exploring the diagram in full, we compare our numerical findings with the analytical predictions of the asymptotic theory, valid for small *T*. The analytic asymptotic profile (3.2) gives a relative error ||** a** − (

***+**

*a**T*)||

*α*_{2}/||

**||**

*a*_{2}less than 0.4% for

*T*≤ 0.2 μm

^{3}h

^{−1}, after which higher-order terms become predominant. This is shown in figure 4 where we compare a branch of approximate solutions (magenta) to a branch of solutions to the full nonlinear problem (blue) for

*D*= 1 μm

^{2}h

^{−1}and small values of

*T*.

As we ascend the bifurcation diagram in figure 3, new peaks are formed on the left side of the existing peaks, that is, towards the interior of the domain, until the whole domain is filled with peaks. We note that peaks are regularly spaced, as it was also observed in references [5,7,60].

This bifurcation diagram resembles the one found for reaction–diffusion PDEs posed on the real line [42–44] except that here peaks are formed at the boundary rather than at the core of the domain. When peaks fill the entire domain, the branch enters an unstable irregular regime without snaking (not shown). Branches of solutions with peaks covering the entire domain are also present (not shown) and are partially discussed in §3.1.4.

### Remark 3.1 (Biological interpretation of snaking)

The diagram in figure 3 makes a plausible biological prediction for the formation of high peaks. A tissue composed of a string of cells, in the presence of passive diffusion, selects auxin patterns depending on the value of the active transport. Our analysis of the Smith model predicts that there exist two main regimes: for small *T*, there is a single low auxin peak at the boundary. As *T* becomes larger, we enter a regime where the tissue can select from a large variety of auxin patterns. If, for instance *T* ≈ 2.2 µm^{3} h^{−1}, the tissue is able to select patterns 1, 2 or 3, which have a variable number of peaks. Therefore, the pattern selected in experiments depends highly upon the initial conditions of the system, similar to what was reported by Jönsson & Krupinski [56]. As we shall see, a slanted version of the snaking bifurcation diagram is also present in two-dimensional domains (both regular and irregular and for both the Smith and the Chitwood models).

#### 3.1.3. Instabilities on the snaking branch

To understand whether a solution pattern is stable to small perturbations, we study in detail the eigenvalues of the Jacobian matrix (result 2.7). The Jacobian matrix for the spatially extended system (2.4) and (2.5) is sparse with a characteristic block structure determined by transport and diffusion terms (we refer the reader to reference [59] for a detailed description) and, for relatively small systems such as this one, eigenvalues are computed with dense linear algebra routines.

In this example, the solution with one small peak at the boundary becomes unstable at a Hopf bifurcation (HP_{1}) at *T* ≈ 2.1 µm^{3} h^{−1}, closely followed by other oscillatory instabilities and a saddle-node bifurcation (TP_{1}) at *T* ≈ 2.2 µm^{3} h^{−1}, after which the solution remains unstable. On the snaking branch, we find that saddle-node and Hopf bifurcations alternate regularly, as documented in figure 3: saddle-node bifurcations align at *T* ≈ 1.9 μm^{3} h^{−1} and *T* ≈ 2.5 µm^{3} h^{−1}, whereas Hopf bifurcations depart from each other as patterns become less localized. In this parameter setting, stable portions of the branch are delimited by two Hopf bifurcations, which, to the best of our knowledge, has not been reported before for snaking in reaction–diffusion systems. It should be noted, however, that Burke & Dawes [61] found Hopf bifurcations at the bottom of the snaking branch for an extended Swift–Hohenberg equation, which may lead to a bifurcation structure similar to the one in figure 3 if secondary parameters are varied.

In figure 5, we show spectra of solutions at selected points on the branch. Overall, these spectra resemble those found in discretized advection–diffusion PDEs, with largely negative real eigenvalues associated with diffusion terms of the governing equations. In this context, however, increasing the number of cells does not alter the cell spacing; hence, the spectrum does not grow in the negative real direction for larger system sizes.

We monitored spectra of localized solutions as the snaking branch was ascended (figure 5): immediately after the Hopf bifurcation HP_{1}, multiple eigenvalues cross the imaginary axis, therefore, several oscillatory instabilities exist between HP_{1} and TP_{1} (HP_{2} and TP_{2}, etc.). In the bottom panel of figure 5, we show that the Hopf eigenfunction at HP_{1} has a maximum near the boundary and the one for HP_{2} is also spatially localized. We expect that branches of time-periodic (possibly spatially localized) solutions emerge from the Hopf bifurcations. We have not observed stable small-amplitude oscillations in direct numerical simulations, but we report the existence of stable periodic states in which a temporal oscillation of the peak at *i* = *n* initiates a wave of auxin moving towards the boundary at *i* = 1, with long oscillation periods.

In figure 6, we show such a periodic solution obtained via time simulation in the neighbourhood of HP_{1} (which is also visible in figure 4). We set and use as the initial condition a steady state (with one peak at the boundary), obtained for . In the resulting periodic state, auxin peaks are dynamically formed from the tip towards the interior of the domain: we point out that the period of oscillations (about 377 h) is much greater than the period associated with the unstable Hopf eigenvalues. In addition, on such long timescales, it is reasonable to assume that new cells are formed, so the geometry of the problem should also change.

It was recently shown by Farcot and Yuan that, in one-dimensional flux-based models with no-flux boundary conditions, active transport is sufficient to elicit auxin oscillations [35]. In the concentration-based model considered here, oscillatory states in regular one-dimensional arrays are also found in a regime where active transport dominates over diffusion.

#### 3.1.4. Changes in the auxin production parameter

We conclude this primer on the one-dimensional tissue by investigating the robustness of the snaking scenario described above (result 2.8). In reference [59], it was shown that the auxin production parameter *ρ*_{IAA} has a significant influence on the solution profiles; therefore, it is interesting to study how changes in this parameter affect the bifurcation structure. We repeated the numerical continuation of the Smith model for 20 values of *ρ*_{IAA} in the interval [0.3 µM h^{−1}, 1.5 µM h^{−1}]. For low values of *ρ*_{IAA}, both the oscillatory instability HP_{1} and the saddle-node TP_{1} move to the right and give rise to snaking bifurcation diagrams with increasingly wider stable segments (figure 7). As a consequence, in biological experiments where the auxin production was kept to a low value, the tissue would support multiple auxin patterns for a wider range of active transport coefficients. In the limit decay dominates over production, hence, large peaks cannot be sustained, and indeed, we find that the solution with a single small peak at the tip persists for very large values of *T*. This is in line with previous papers by Sahlin *et al.* [37] and De Reuille *et al.* [6] where it was postulated that auxin patterning demands a minimal level of auxin production within the tissue.

On the other hand, increasing *ρ*_{IAA} causes the snaking diagram to shrink and then disappear for *ρ*_{IAA} ≥ 1.2 µM h^{−1}. In figure 7, we show a fully stable branch for *ρ*_{IAA} = 1.5 µM h^{−1}. On this branch, peaks develop at once from the small-amplitude solution, without turning points. We mention, however, that for *ρ*_{IAA} between 1.2 and 1.3 µM h^{−1}, Hopf bifurcations are found along the non-snaking branch (not shown), similar to what is found for the infinite domain [59].

As snaking branches distort, several types of secondary instabilities and collisions with neighbouring branches occur. In particular, we report codimension-2 Bogdanov–Takens bifurcations originating from the collision between TP_{2} and HP_{2}, (TP_{4} and HP_{4}, TP_{6} and HP_{6}, etc.) when *ρ*_{IAA} is varied. The existence of these codimension-2 bifurcations could also be envisaged from the spectra in figure 5. These instabilities, as the ones reported in the previous section, indicate that the tissue is capable of sustaining oscillations and dynamical auxin patterning, as well as steady states with multiple peaks. Dynamic states with spatio-temporal coherence (such as the one reported in figure 6) are interesting from a biological standpoint [20], as they occur for the biologically plausible parameter values reported in table 1. However, we could not find them in two-dimensional domains with realistic parameter values; hence, we do not discuss them further in this paper.

### 3.2. Two-dimensional domains

We now move to more realistic geometries and study two-dimensional domains with approximately square and circular boundaries, on which we prescribe free boundary conditions. Here, we revisit results 2.3–2.9 in the two-dimensional setting for both the Smith and the Chitwood models, so we refer the reader to the general summary in §2 and the one-dimensional primer in §3.1.

The methods described in §3.1 apply straightforwardly to the two-dimensional case. In the first example, we consider the Smith model on a grid of 50 × 50 hexagonal prismic cells with *l _{ij}* = 1 μm and . Cells have six neighbours in the interior, three neighbours at the left and right edges and four neighbours at the top and bottom edges. In this domain, corners are not all equal (figure 8), and we chose this configuration intentionally, to illustrate the influence of the number of neighbouring cells on the emerging patterns in two-dimensional domains.

Figure 9 shows the values of the geometric pre-factors ** ξ** for two corners of the domain: our asymptotic analysis for

*D*= 0 µm

^{2}h

^{−1}(see the discussion in §3.1.1 and the generic derivation in §5) predicts the formation of peaks at the boundaries, with the highest peak at the top-left and bottom-right corners. Numerical computations for positive

*D*show that these peaks persist and become prominent for increasing

*T*(figure 8, pattern 2). As in the one-dimensional case, patterns are arranged on a snaking bifurcation branch, even though in two-dimensions the snaking is slanted. Peaks arise initially in all four corners, then, new spots are formed along the left and right edges (where cells have fewer neighbours), and then, for sufficiently large values of

*ρ*

_{IAA}, along the top and bottom edges. In contrast to the one-dimensional case, we have not found oscillatory bifurcations in this region of parameter space, so we conclude that stable portions of the branch are now delimited by turning points (figure 8). From a biological perspective, this means that patterns with peaks at the boundary are more likely to be observable in experiments, as they are stable in a wider region of parameter space.

In a second example, we consider the Chitwood model posed on the two-dimensional regular domain, using the parameters of table 1. Remarkably, the resulting bifurcation diagram (not shown) is analogous to the one found in figure 8 for the Smith model. As a further confirmation, we tested the robustness of the snaking mechanism to changes in the auxin production coefficient, as it was done for the one-dimensional case in figure 7. We set *ρ*_{IAA} = 2 µM h^{−1} and show in figure 10 the corresponding bifurcation diagram and solution patterns. While active transport remains responsible for the selection of peaks towards the boundary, the interplay with auxin production allows the formation of a ring of spots at the boundary as opposed to single spots at the corners (see pattern 1 in figure 10). After the first turning point, spots are formed in pairs (pattern 2), until they fill a full row (pattern 3) and the whole domain (not shown). As in §3.1.4, the auxin production coefficient has a large influence on the resulting peaks. For this two-dimensional regular domain, we also scanned several values of the auxin production coefficient and confirmed that comparatively high values of *ρ*_{IAA} induce the formation of peaks all-at-once (similar to top panel of figure 7); hence, a fully patterned tissue is possible without a Turing bifurcation.

In the remaining two examples, 742 irregular prismic cells cover an almost-circular domain in a realistic tissue with free boundary conditions (figure 11, whose geometry has been extracted from reference [58]). Even though the asymptotic analysis is not valid for irregular arrays, we expect results to be qualitatively similar if cellular volumes and contact areas do not vary greatly from cell to cell. In these examples, the number of neighbours varies over the domain; however, the cells at the boundaries have predominantly fewer neighbours and this is where peaks are formed initially. On the *x*-axis of the bifurcation diagram, we now use the scaled parameter , where denotes the average in the tissue (in the cases analysed thus far, , so the numerical values of *T* and coincided). We note that in this domain and As in the two-dimensional regular cases, stable portions of the branch are enclosed between consecutive saddle-node bifurcations, and there are no oscillatory instabilities on the stable branches. Figure 11 shows the results of the Smith model: as usual, peaks are formed initially at the boundary and then fill the interior (see pattern 4), and the slanted snaking ensures the existence of stable solutions with localized peaks in a wide regime of the parameter *T*.

When we pose the Chitwood model to the same irregular domain, the results are strikingly similar, as seen in figure 12: peaks form initially at the top-left quadrant of the circular domain (see patterns 1 and 2 in figures 11 and 12), confirming that the geometry of the tissue drives the location of the spots. An inspection of the fully patterned tissues (pattern 4 in figures 11 and 12) reveals that model parameters and functional forms for the active transport functions have an influence on the size and structure of the peaks. Variations in *ρ*_{IAA} also confirmed the trend seen in figure 7 (not shown).

## 4. Discussion

In this paper, we investigated the origin of auxin peaks in generic concentration-based models and proposed a robust mechanism for their formation over short timescales, using a combination of asymptotic and numerical bifurcation analysis.

The asymptotic calculations, valid for a class of models with identical cells and weak active transport, show that peaks emerge as boundary corrections to the homogeneous steady state: the peak amplitude depends on the local geometry and is higher in regions where cells have fewer neighbours, that is, next to the boundary. Crucially, this is a direct consequence of the mathematical structure of the models considered here (hypotheses 2.1 and 2.2): because the active transport depends on the number of neighbours at distance 2 from the *i*th cell via the geometric coefficients *ξ _{i}*, then deviations from the homogeneous state will always appear at the boundaries, where the number of neighbours is different from the interior of the tissue. This mechanism is different from (and not in contrast with) the Turing bifurcation scenario reported in previous studies on unbounded domains: on finite tissues, peaks do not emerge from instabilities of the flat state, but they simply morph from it for low values of

*T*. The most immediate consequence of our mathematical analysis is that, in concentration-based models, active transport and geometry concur to promote localization of auxin peaks [10]. The results are immediately applicable to models of generic plant tissues, or to other models of morphogen-transporter systems, provided they satisfy hypotheses 2.1 and 2.2.

In irregular domains, a similar asymptotic analysis can be carried out, but the peak selection mechanism in this case also depends on the cellular volumes and contact areas, so we cannot exclude *a priori* that peaks will form in the interior as well as on the boundary. A statistical characterization of the peaks location in relation to the variance of the cellular array is possible and should be considered in future studies. The two models by Smith [5] and Chitwood [39] fit in the above-discussed framework and, for these systems, we have provided numerical evidence that peaks persist for moderate and large values of the active transport rate *T*: large-amplitude peaks are arranged on a snaking branch, which becomes slanted in two-dimensional tissues. A major implication of the numerical findings of figures 11 and 12 is that localized auxin peaks are observable: a slow experimental sweep in the active transport from low to high values should reveal an initial localization of the peaks, which then progressively de-localize and fill the domain. Conversely, if the active transport is kept constant, then the tissue should be able to select from more than one pattern, depending on the initial condition.

We speculate that multistability could have a functional role: we found that the tissue can switch between localized to non-localized states by either changing transport, or auxin production and this could be at the origin of the multiple number of primordia arising at the SAM.

Because snaking is now recognized as the footprint of localization in a wide variety of nonlinear media, we expect that bifurcation diagrams similar to the ones shown in this paper for the Smith and the Chitwood models might also arise generically. It is known that localized states occur also at the presence of external forcings [62]; hence, an interesting avenue of research would be to test the robustness of the snaking mechanism when the system is subject to a time-dependent stimulus. We expect to find a similar mechanism (with or without external stimuli) in other transport systems and also in transport models applied to other plant tissues, and particularly root tissues, for which coexisting multiple steady states have been reported but not yet analysed [63]. One class of models that are of interest and can be analysed via numerical bifurcation analysis are flux-based auxin models, which have not been studied in this paper. While numerical continuation is readily applicable to such models, they are likely to require a separate analytic treatment, as some of them do not possess the factorization presented in hypothesis 2.2.

Importantly, we find that the bifurcation scenario is influenced by the auxin production rate, because the selectable configurations depend sensitively on the balance between auxin production and active transport. The results in figure 7 (also confirmed by the two-dimensional irregular calculations of figures 11 and 12) support the conclusion that if the auxin production rate was decreased quasi-statically, either actively or passively, the organism would be able to switch from fully patterned states to configurations with few peaks at the boundary. We find that our theoretical predictions are supported by the experimental results in reference [60] (figure 3*a*,*c*): depending on the amount of auxin in the system, the resulting auxin pattern varies from one localized peak (giving an isolated primordium) to solutions with peaks at the boundary (ring-shaped primordium).

In addition, because the parameters of the concentration-based models considered here are scaled by cellular volumes, we expect tissues with different cell sizes to behave similarly: in tissues with larger *V*, the snaking limits are expected to occur for larger values of *T*, so as to keep the ratio *T*/*V* constant (and a similar reasoning is valid for the passive transport parameter *D*).

The above-reported conclusions are naturally limited to experiments that are well approximated by concentration-based models, and for the plausible biological parameter values selected in the original papers by Smith *et al.* [5] and Chitwood *et al.* [39]. We note that an experimental validation of the predictions presented here requires the ability to detect changes in the auxin distribution during development. An experimental technique that could help testing the predictions of these two models is the one recently proposed by Brunoud *et al.* [31], which allows visualization of auxin with high spatio-temporal resolution. We note that it would be possible to also apply numerical bifurcation analysis to a modified model that accounted for markers’ dynamics and auxin–sensor interactions.

A further desirable property of the experimental set-up would be the ability to stimulate auxin peaks locally (thereby changing initial conditions) and test whether the tissue settles to a new equilibrium. In view of the large uncertainty regarding the parameter values of the models, we expect our predictions to agree only qualitatively with experimental results. Finally, by comparing the solution landscape found with the bifurcation analysis for different models and the experimental results, we can determine the biological plausibility of each model and exclude models.

## 5. Material and methods

Here, we present an analytical framework to construct steady-state solutions featuring localized auxin peaks in generic concentration-based models.

### 5.1. Asymptotic derivation of peak solutions

We begin by giving a generic definition of concentration-based models in the absence of diffusion: as mentioned above, several examples from literature can be cast in this form.

### Definition 5.1 (Concentration-based model without diffusion)

A concentration-based model without diffusion is a set of *m**n* ODEs of the form
5.1for *i* = 1, …, *n*, where , are the production and decay functions, respectively, , is the (non-negative) active transport parameter, {1, … , n} are vertices of a static undirected graph G, is the set of neighbours of cell i, containing elements and are the active transport functions. We will assume ** π**,

**and**

*δ*

*ν*_{ij}to be smooth vector fields depending on a set of control parameters , but we omit this dependence for simplicity and write, for instance,

**(**

*π*

*y*_{i}) instead of

**(**

*π*

*y*_{i};

**).**

*p*### Remark 5.2.

Specific examples of concentration-based models in this form can be found in §5.1.1 and §5.1.2.

We now prove the following result.

### Lemma 5.3.

*Let us consider the concentration-based model* (*5.1*) *and let us suppose that there exist vector-valued functions* *and* *such that*
5.2*for all* *i*,*j* = 1, …, *n*, *where* *and* *denote the standard Hadamard product and division between vectors. Further, let* *be such that* ** π**(

***) =**

*y***(**

*δ****),**

*y***(**

*ψ****,**

*y****) ≠**

*y***0**

*and*

*(*

*φ****) ≠ 0**

*y**, then*

(

*1*)*If T*=*0 or all cells have the same number of neighbours,**, then the homogeneous solution**is a steady state for the concentration-based model*.(

*2*)*If**and cells have different numbers of neighbours and the Jacobian matrix*′(*π***y***)−′(*δ***y***) is non-singular, then a inhomogeneous steady state (to leading order) is given by

*where the coefficients ξ_{i} depend on the local properties of the cellular array, namely*

### Proof.

If *T* = 0, the statement is clearly true, so henceforth, we assume *T* ≠ 0. Because ** π**(

***)**

*y**−*

**(**

*δ****) =**

*y***0**, the right-hand side of (5.1) vanishes for all

*i*if . On the other hand, if not all cells have the same number of neighbours and

*T*is small, then we may seek for asymptotic steady states in the form for

*i*= 1,

*…*,

*n*and . A Taylor expansion of the right-hand side around gives, to leading order, 5.3

In order to find an expression for *η** _{i}*, we evaluate the sums in (5.3):
5.4and combining (5.3) with (5.4), we obtain the assert. ▪

### Remark 5.4 (Small-amplitude peak solutions)

In finite regular arrays, cells in the interior all have the same number of neighbours, so we can use these properties to give formal definitions of interior and boundary sets
*In passing, we note that* *contains in general more cells than the physical boundary. Lemma 5.3 shows that, to leading order, steady states for small T deviate from the homogeneous solutions only in* *, where auxin peaks and dips are proportional to T and ξ_{i}, namely*
5.5

### Remark 5.5

(**Irregular domains**)**.** Lemma 2.1 cannot be applied in general if the domain is irregular, that is, if cells have different volumes and contact lengths: if, say, the active transport function **v**_{ij} depends explicitly on the cellular volume V_{i} and the V_{i} are not all equal, then it is not possible to express **v**_{ij} as in (5.2). However, the theory can be extended to the case of irregular domains. We do not report this derivation here, but note that it features, as expected, cellular volumes and contact areas.

#### 5.1.1. One-dimensional domain and one component per cell

As an example, we consider the Smith model [5] with constant fixed PIN1 amount, posed on a one-dimensional array of identical cells with volume *V*, Neumann boundary conditions at *i* = 1 and free boundary conditions at *i* = *n*. The Neumann boundary conditions are obtained by considering ghost cells [59], therefore, boundary and interior sets are given by and , respectively. Furthermore, we denote by *p* the fixed PIN concentration and apply lemma 5.3 with *m* = 1, *y _{i}* =

*a*and where all parameters are assumed to be strictly positive. By balancing production and decay terms, we find the positive homogeneous state

_{i}In the absence of active transport, *a** is a stable steady state of the model, because *π*′(*a**)*−δ*′(*a**) < 0. For , ** a*** is not a steady state, because is non-empty and

*ξ*

_{n−}_{1}= −1/2,

*ξ*= 1/2, hence we obtain, to leading order

_{n}

#### 5.1.2. One-dimensional domain and two components per cell

The Smith model [5] features two ODEs per cell. If we pose this model to a one-dimensional array of identical cells, then the graph *G* associated with the nodes is the same as in our previous example, hence, , and **ξ** are unchanged. We can now apply lemma 5.3 with *m* = 2, *y** _{i}* = (

*a*,

_{i}*p*)

_{i}^{T}and Balancing production and decay terms, we find a homogeneous strictly positive steady state for

*T*= 0 which is stable, because

Because the parameters are assumed to be positive with the exception of *ρ*_{PIN0} which is non-negative (see also table 1), we do not have a zero eigenvalue.

The inverse of ** π**′(

***)**

*y**−*

**′(**

*δ****) can be computed explicitly, and for we obtain to leading order**

*y*

Another example of a concentration-based transport model that features two ODEs per cell is the Chitwood model [39]. The model is given by

For the same domain as above and with *m* = 2 and *y** _{i}* = (

*a*,

_{i}*p*)

_{i}^{T}, we can apply lemma 5.3.

The homogeneous steady state for *T* = 0 is the same as for the Smith model and for we obtain to leading order

#### 5.1.3. Two-dimensional domain of identical hexagonal cells

Lemma 5.3 also applies when the Smith or the Chitwood models are posed on a two-dimensional array of identical hexagonal cells. In this case, the computations for ** y*** are identical to the previous example, and the asymptotic derivation is also straightforward. Instead of writing the full expressions for , and

**, we refer the reader to figure 9, where the values of**

*ξ***are shown for two corners of the domain. As claimed in §3.2, the highest peaks occur at the top-left and right-bottom corners.**

*ξ*### 5.2. Models with diffusion

We can extend the definition of concentration-based models to the case where diffusion is present.

### Definition 5.6 (Concentration-based model with diffusion)

A concentration-based model with diffusion is a set of *mn* ODEs of the form
5.6for *i* = 1, …, *n*, where is a diagonal diffusion matrix and all other quantities are as in definition 5.1.

### Remark 5.7.

Reasoning as in lemma 5.3, we obtain where *η*_{i} satisfy
is block-diagonal with blocks ** π**′(

***) −**

*y***′(**

*δ****), is the Laplacian matrix associated with the graph G with Neumann boundary conditions and denotes the Kronecker product between matrices. The operator is negative semi-definite and it has a zero eigenvalue, corresponding to a constant eigenvector. However, summing this matrix to**

*y***(**

*J***) =**

*y***′(**

*π****)−**

*y***′(**

*δ****), makes the resulting linear operator non-singular.**

*y*In the presence of diffusion, we cannot directly apply the formula (5.5), even for regular cellular arrays: owing to diffusion, cells in will also deviate from the homogeneous state; hence, peaks and dips are not necessarily formed within , but may occur in interior cells that are close to the boundary. First-order corrections for these cases can be computed analytically using Chebyshev polynomials [64] or numerically using linear algebra routines. Even though we report below an example of this calculation, we point out that in practice this is not necessary, because the numerical bifurcation software gives access to the full nonlinear solution and to its linear stability.

#### 5.2.1. One-dimensional domain with diffusion and two components per cell

We return to the Smith model with *m* = 2, posed on a row of identical cells, and we now add diffusion only in the auxin component.

The expressions for *a** and *p** are unchanged from §5.1.2, as is the first-order approximation of *p _{i}* (because there is no diffusion for

*p*). Expressions for the first-order approximations in

*a*are more involved: proceeding as explained above for generic models with diffusion, we obtain, to leading order

_{i}Solving this linear equation above led to the approximate solution profiles in figure 2 and the red solution branch in figure 4. The same derivation and figures can be obtained for the Chitwood model (not shown).

## Funding statement

D.D. acknowledges financial support from the Department of Mathematics and Computer Science of the University of Antwerp. D.A. acknowledges the University of Nottingham Research Development Fund, supported by the Engineering and Physical Sciences Research Council (EPSRC).

## Acknowledgements

We acknowledge fruitful discussions with Gerrit T.S. Beemster, Jan Broeckhove, Dirk De Vos, Abdiravuf Dzhurakhalov, Etienne Farcot and Przemyslaw Klosiewicz.

## Footnotes

- Received December 28, 2014.
- Accepted March 19, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.