Spot dynamics in a reaction-diffusion model of plant root hair   initiation

D. Avitabile; V. F. Bre\~na-Medina; M. J. Ward

arXiv:1703.02608·nlin.PS·March 20, 2018·SIAM J. Appl. Math.

Spot dynamics in a reaction-diffusion model of plant root hair initiation

D. Avitabile, V. F. Bre\~na-Medina, M. J. Ward

PDF

TL;DR

This paper models how a plant hormone gradient influences the formation and alignment of root hair initiation sites in Arabidopsis thaliana using a reaction-diffusion system with nonlinear kinetics.

Contribution

It extends a reaction-diffusion model to include spatial heterogeneity, revealing how auxin gradients guide localized ROP activation patterns in plant cells.

Findings

01

Auxin gradient causes spatial alignment of active ROP regions.

02

Localized patterns of ROPs are influenced by nonlinear reaction kinetics.

03

Numerical simulations confirm pattern formation and alignment.

Abstract

We study pattern formation in a 2-D reaction-diffusion (RD) sub-cellular model characterizing the effect of a spatial gradient of a plant hormone distribution on a family of G-proteins associated with root-hair (RH) initiation in the plant cell Arabidopsis thaliana. The activation of these G-proteins, known as the Rho of Plants (ROPs), by the plant hormone auxin, is known to promote certain protuberances on root hair cells, which are crucial for both anchorage and the uptake of nutrients from the soil. Our mathematical model for the activation of ROPs by the auxin gradient is an extension of the model of Payne and Grierson [PLoS ONE, 12(4), (2009)], and consists of a two-component Schnakenberg-type RD system with spatially heterogeneous coefficients on a 2-D domain. The nonlinear kinetics in this RD system model the nonlinear interactions between the active and inactive forms of ROPs.…

Figures36

Click any figure to enlarge with its caption.

Tables1

Table 1. Table 1: Two parameter sets in the original and dimensionless re-scaled variables. The fundamental units of length and time are μ m 𝜇 m \mu\textrm{m} and sec, and concentration rates are measured by an arbitrary datum (con) per time unit; k 20 subscript 𝑘 20 k_{20} is measured by con 2 / s superscript con 2 s \textrm{con}^{2}/\textrm{s} , and diffusion coefficients units are μ m / s 2 𝜇 m superscript s 2 \mu\textrm{m}/\textrm{s}^{2} .

	Parameter Set A		Parameter Set B
Original	Re-scaled	Original	Re-scaled
$D_{1} = 0.075$	$ε^{2} = 1.02 \times 10^{- 4}$	$D_{1} = 0.1$	$ε^{2} = 3.6 \times 10^{- 4}$
$D_{2} = 20$	$D = 0.51$	$D_{2} = 10$	$D = 0.4$
$k_{1} = 0.008$	$τ = 18.75$	$k_{1} = 0.01$	$τ = 11$
$b = 0.008$	$β = 6.25$	$b = 0.01$	$β = 1$
$c = 0.1$		$c = 0.1$
$r = 0.05$		$r = 0.01$
$k_{20} \in [10^{- 3}, 2.934]$	$γ \in [150, 0.051]$	$k_{20} \in [0.01, 1.0]$	$γ \in [11, 0.11]$
$L_{x} = 70$		$L_{x} = 50$
$L_{y} = 29.848$	$d_{y} = 0.4211$	$L_{y} = 20$	$d_{y} = 0.40$

Equations113

(i) α (X) = exp (- \frac{ν X}{L _{x}}), (ii) α (X) = exp (- \frac{ν X}{L _{x}}) sin (\frac{π Y}{L _{y}}), ν = 1.5 .

(i) α (X) = exp (- \frac{ν X}{L _{x}}), (ii) α (X) = exp (- \frac{ν X}{L _{x}}) sin (\frac{π Y}{L _{y}}), ν = 1.5 .

U_{t}

U_{t}

τ V_{t}

ε^{2} \equiv \frac{D _{1}}{L _{x}^{2} ( c + r )}, D \equiv \frac{D _{2}}{L _{x}^{2} k _{1}}, τ \equiv \frac{c + r}{k _{1}}, β \equiv \frac{r}{k _{1}},

ε^{2} \equiv \frac{D _{1}}{L _{x}^{2} ( c + r )}, D \equiv \frac{D _{2}}{L _{x}^{2} k _{1}}, τ \equiv \frac{c + r}{k _{1}}, β \equiv \frac{r}{k _{1}},

γ \equiv \frac{( c + r ) k _{1}^{2}}{k _{20} b ^{2}} .

\int_{Ω} U_{0} (x) d x = \frac{∣Ω∣}{β γ},

\int_{Ω} U_{0} (x) d x = \frac{∣Ω∣}{β γ},

U = ε^{- 2} u, V = ε^{2} v .

U = ε^{- 2} u, V = ε^{2} v .

u_{t}

u_{t}

τ ε^{2} v_{t}

u = u_{0 j} (ρ) + ε u_{1 j} + \dots, v = v_{0 j} (ρ) + ε v_{1 j} + \dots, ρ = ∣ ξ ∣ \equiv ε^{- 1} ∣ x - x_{j} ∣ .

u = u_{0 j} (ρ) + ε u_{1 j} + \dots, v = v_{0 j} (ρ) + ε v_{1 j} + \dots, ρ = ∣ ξ ∣ \equiv ε^{- 1} ∣ x - x_{j} ∣ .

Δ_{ρ} u_{0 j} + α (x_{j}) u_{0 j}^{2} v_{0 j} - u_{0 j} = 0, D_{0} Δ_{ρ} v_{0 j} - τ γ (α (x_{j}) u_{0 j}^{2} v_{0 j} - u_{0 j}) - β γ u_{0 j} = 0,

Δ_{ρ} u_{0 j} + α (x_{j}) u_{0 j}^{2} v_{0 j} - u_{0 j} = 0, D_{0} Δ_{ρ} v_{0 j} - τ γ (α (x_{j}) u_{0 j}^{2} v_{0 j} - u_{0 j}) - β γ u_{0 j} = 0,

α (x) = α (x_{j}) + ε \nabla α (x_{j}) \cdot ξ + \dots, where ξ \equiv ε^{- 1} (x - x_{j}) .

α (x) = α (x_{j}) + ε \nabla α (x_{j}) \cdot ξ + \dots, where ξ \equiv ε^{- 1} (x - x_{j}) .

u_{0 j}^{'} (0) = v_{0 j}^{'} (0) = 0; u_{0 j} \to 0, v_{0 j} \sim S_{j} lo g ρ + χ (S_{j}) + o (1), as ρ \to \infty,

u_{0 j}^{'} (0) = v_{0 j}^{'} (0) = 0; u_{0 j} \to 0, v_{0 j} \sim S_{j} lo g ρ + χ (S_{j}) + o (1), as ρ \to \infty,

ρ \to \infty lim (v_{0 j} - S_{j} lo g ρ) = χ (S_{j}) .

ρ \to \infty lim (v_{0 j} - S_{j} lo g ρ) = χ (S_{j}) .

S_{j} = \frac{β γ}{D _{0}} \int_{0}^{\infty} [\frac{τ}{β} (α (x_{j}) u_{0 j}^{2} v_{0 j} - u_{0 j}) + u_{0 j}] ρ d ρ .

S_{j} = \frac{β γ}{D _{0}} \int_{0}^{\infty} [\frac{τ}{β} (α (x_{j}) u_{0 j}^{2} v_{0 j} - u_{0 j}) + u_{0 j}] ρ d ρ .

u_{0 j} \equiv \frac{D _{0}}{β γ α ( x _{j} )} \lower 0.4 pt \vrule height=13.85553pt,depth=-11.08447pt u_{c}, v_{0 j} \equiv \frac{β γ}{D _{0} α ( x _{j} )} \lower 0.4 pt \vrule height=14.10666pt,depth=-11.28539pt v_{c},

u_{0 j} \equiv \frac{D _{0}}{β γ α ( x _{j} )} \lower 0.4 pt \vrule height=13.85553pt,depth=-11.08447pt u_{c}, v_{0 j} \equiv \frac{β γ}{D _{0} α ( x _{j} )} \lower 0.4 pt \vrule height=14.10666pt,depth=-11.28539pt v_{c},

Δ_{ρ} u_{c} + u_{c}^{2} v_{c} - u_{c} = 0, Δ_{ρ} v_{c} - \frac{τ}{β} (u_{c}^{2} v_{c} - u_{c}) - u_{c} = 0,

Δ_{ρ} u_{c} + u_{c}^{2} v_{c} - u_{c} = 0, Δ_{ρ} v_{c} - \frac{τ}{β} (u_{c}^{2} v_{c} - u_{c}) - u_{c} = 0,

u_{c}^{'} (0) = v_{c}^{'} (0) = 0, u_{c} \to 0, v_{c} \sim S_{c j} lo g ρ + χ_{c} (S_{c j}) + o (1), as ρ \to \infty .

S_{j} \equiv \frac{β γ}{D _{0} α ( x _{j} )} \lower 0.4 pt \vrule height=14.10666pt,depth=-11.28539pt S_{c j}, S_{c j} = \int_{0}^{\infty} [\frac{τ}{β} (u_{c}^{2} v_{c} - u_{c}) + u_{c}] ρ d ρ,

S_{j} \equiv \frac{β γ}{D _{0} α ( x _{j} )} \lower 0.4 pt \vrule height=14.10666pt,depth=-11.28539pt S_{c j}, S_{c j} = \int_{0}^{\infty} [\frac{τ}{β} (u_{c}^{2} v_{c} - u_{c}) + u_{c}] ρ d ρ,

χ (S_{j}) \equiv \frac{β γ}{D _{0} α ( x _{j} )} \lower 0.4 pt \vrule height=14.10666pt,depth=-11.28539pt χ_{c}, where χ_{c} = χ_{c} (S_{c j}) .

- ε^{2} v_{0} + 1 - ε^{- 2} [τ γ (α (x) u_{0}^{2} v_{0} - u_{0}) + β γ u_{0}] = 1 + O (ε^{2}) .

- ε^{2} v_{0} + 1 - ε^{- 2} [τ γ (α (x) u_{0}^{2} v_{0} - u_{0}) + β γ u_{0}] = 1 + O (ε^{2}) .

ε^{- 2} R_{j} (ε^{- 1} (x - x_{j})) \to 2 π (\int_{0}^{\infty} R_{j j} (ρ) ρ d ρ) \updelta (x - x_{j}) = 2 π D_{0} S_{j} \updelta (x - x_{j}),

ε^{- 2} R_{j} (ε^{- 1} (x - x_{j})) \to 2 π (\int_{0}^{\infty} R_{j j} (ρ) ρ d ρ) \updelta (x - x_{j}) = 2 π D_{0} S_{j} \updelta (x - x_{j}),

Δ v_{0} + \frac{1}{D _{0}} = 2 π i = 1 \sum N S_{i} \updelta (x - x_{i}), x \in Ω; \partial_{n} v_{0} = 0, x \in \partial Ω .

Δ v_{0} + \frac{1}{D _{0}} = 2 π i = 1 \sum N S_{i} \updelta (x - x_{i}), x \in Ω; \partial_{n} v_{0} = 0, x \in \partial Ω .

v_{0} \sim S_{j} lo g ∣ x - x_{j} ∣ + \frac{S _{j}}{\upsigma} + χ (S_{j}), for j = 1, \dots, N, where \upsigma \equiv - \frac{1}{l o g ε},

Δ G = \frac{1}{d _{y}} - \updelta (x - y), \int_{Ω} G (x; y) d y = 0, x \in Ω, \partial_{n} G = 0, x \in \partial Ω .

Δ G = \frac{1}{d _{y}} - \updelta (x - y), \int_{Ω} G (x; y) d y = 0, x \in Ω, \partial_{n} G = 0, x \in \partial Ω .

G (x; y) \sim - \frac{1}{2 π} lo g ∣ x - y ∣ + R (y; y) + \nabla_{x} R (y; y) \cdot (x - y) + O (∣ x - y ∣^{2}) .

G (x; y) \sim - \frac{1}{2 π} lo g ∣ x - y ∣ + R (y; y) + \nabla_{x} R (y; y) \cdot (x - y) + O (∣ x - y ∣^{2}) .

v_{0} = - 2 π i = 1 \sum N S_{i} G (x; x_{i}) + \overset{v}{ˉ}_{0}, where i = 1 \sum N S_{i} = \frac{d _{y}}{2 π D _{0}} .

v_{0} = - 2 π i = 1 \sum N S_{i} G (x; x_{i}) + \overset{v}{ˉ}_{0}, where i = 1 \sum N S_{i} = \frac{d _{y}}{2 π D _{0}} .

\frac{S _{j}}{\upsigma} + χ (S_{j}) = - 2 π S_{j} R_{j} - 2 π i = 1 \sum N S_{i} G_{j i} + \overset{v}{ˉ}_{0}, for j = 1, \dots N .

\frac{S _{j}}{\upsigma} + χ (S_{j}) = - 2 π S_{j} R_{j} - 2 π i = 1 \sum N S_{i} G_{j i} + \overset{v}{ˉ}_{0}, for j = 1, \dots N .

(G)_{ii} = R_{i}, (G)_{ij} = G_{ij}, \mbox w h er e G_{ij} = G_{j i} .

(G)_{ii} = R_{i}, (G)_{ij} = G_{ij}, \mbox w h er e G_{ij} = G_{j i} .

S + 2 π \upsigma GS = \upsigma (\overset{v}{ˉ}_{0} e - χ), e^{T} S = \frac{d _{y}}{2 π D _{0}} .

S + 2 π \upsigma GS = \upsigma (\overset{v}{ˉ}_{0} e - χ), e^{T} S = \frac{d _{y}}{2 π D _{0}} .

\overset{v}{ˉ}_{0} = \frac{1}{N \upsigma} (\frac{d _{y}}{2 π D _{0}} + 2 π \upsigma e^{T} GS + \upsigma e^{T} χ) .

\overset{v}{ˉ}_{0} = \frac{1}{N \upsigma} (\frac{d _{y}}{2 π D _{0}} + 2 π \upsigma e^{T} GS + \upsigma e^{T} χ) .

S + \upsigma (I - E) (2 π GS + χ) = \frac{d _{y}}{2 π N D _{0}} e .

S + \upsigma (I - E) (2 π GS + χ) = \frac{d _{y}}{2 π N D _{0}} e .

\upalpha \equiv diag (\frac{1}{α ( x _{1} ) \lower 0.4 pt \vrule height=7.5pt,depth=-6.00003pt}, \dots, \frac{1}{α ( x _{N} ) \lower 0.4 pt \vrule height=7.5pt,depth=-6.00003pt}), S_{c} \equiv S_{c 1} ⋮ S_{c N}, χ_{c} \equiv χ_{c} (S_{c 1}) ⋮ χ_{c} (S_{c N}) .

\upalpha \equiv diag (\frac{1}{α ( x _{1} ) \lower 0.4 pt \vrule height=7.5pt,depth=-6.00003pt}, \dots, \frac{1}{α ( x _{N} ) \lower 0.4 pt \vrule height=7.5pt,depth=-6.00003pt}), S_{c} \equiv S_{c 1} ⋮ S_{c N}, χ_{c} \equiv χ_{c} (S_{c 1}) ⋮ χ_{c} (S_{c N}) .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

\headers

Spot dynamics in a reaction–diffusion modelD. Avitabile, V. F. Breña–Medina, M. J. Ward

Spot dynamics in a reaction-diffusion model of plant root hair

initiation

Daniele Avitabile Centre for Mathematical Medicine and Biology, School of Mathematical Sciences, University of Nottingham, University Park, Nottingham, NG7 2RD, , ORCID ID http://orcid.org/0000-0003-3714-7973.

[email protected]

Victor F. Breña–Medina Centro de Ciencias Matemáticas, UNAM, Morelia, Michoacán, México. Present address: Departamento Académico de Matemáticas, ITAM, Ciudad de México, México, . [email protected]

Michael J. Ward Department of Mathematics, UBC, Vancouver, Canada, .

[email protected]

Abstract

We study pattern formation in a 2-D reaction-diffusion (RD) sub-cellular model characterizing the effect of a spatial gradient of a plant hormone distribution on a family of G-proteins associated with root-hair (RH) initiation in the plant cell Arabidopsis thaliana. The activation of these G-proteins, known as the Rho of Plants (ROPs), by the plant hormone auxin, is known to promote certain protuberances on root hair cells, which are crucial for both anchorage and the uptake of nutrients from the soil. Our mathematical model for the activation of ROPs by the auxin gradient is an extension of the model of Payne and Grierson [PLoS ONE, 12(4), (2009)], and consists of a two-component Schnakenberg-type RD system with spatially heterogeneous coefficients on a 2-D domain. The nonlinear kinetics in this RD system model the nonlinear interactions between the active and inactive forms of ROPs. By using a singular perturbation analysis to study 2-D localized spatial patterns of active ROPs, it is shown that the spatial variations in the nonlinear reaction kinetics, due to the auxin gradient, lead to a slow spatial alignment of the localized regions of active ROPs along the longitudinal midline of the plant cell. Numerical bifurcation analysis, together with time-dependent numerical simulations of the RD system are used to illustrate both 2-D localized patterns in the model, and the spatial alignment of localized structures.

1 Introduction

We examine the effect of a spatially-dependent plant hormone distribution on a family of proteins associated with root hair (RH) initiation in a specific plant cell. This process is modeled by a generalized Schnakenberg reaction-diffusion (RD) system on a 2-D domain with both source and loss terms, and with a spatial gradient modeling the spatially inhomogeneous distribution of the plant hormone auxin. This system is an extension of a model proposed by Payne and Grierson in [33], and analyzed in a 1-D context in the companion articles [4, 6]. The new goal of this paper, in comparison with [4, 6], is to analyze 2-D localized spot patterns in the RD system Eq. 1.2, and how these 2-D patterns are influenced by the spatially inhomogeneous auxin distribution.

We now give a brief description of the biology underlying the RD model. In this model, an on-and-off switching process of a small G-protein subfamily, called the Rho of Plants (ROPs), is assumed to occur in a RH cell of the plant Arabidopsis thaliana. ROPs are known to be involved in RH cell morphogenesis at several distinct stages (see [13, 23] for details). Such a biochemical process is believed to be catalyzed by a plant hormone called auxin (cf. [33]). Typically, auxin-transport models are formulated to study polarization events between cells (cf. [14]). However, little is known about specific details of auxin flow within a cell. In [4, 6, 33] a simple transport process is assumed to govern the auxin flux through a RH cell, which postulates that auxin diffuses much faster than ROPs in the cell, owing partially to the in- and out-pump mechanism that injects auxin into RH cells from both RH and non-RH cells (cf. [16, 18]). Recently, in [27], a model based on the one proposed in [33] has been used to describe patch location of ROP GTPases activation along a 2-D root epidermal cell plasma membrane. This model is formulated as a two-stage process, one for ROP dynamics and another for auxin dynamics. The latter process is assumed to be described by a constant auxin production at the source together with passive diffusion and a constant auxin bulk degradation away from the source. In [27] the auxin gradient is included in the ROP finite-element numerical simulation after a steady-state is attained. Since we are primarily interested in the biochemical interactions that promote RHs and key ingredients that geometrical features have on the RHs initiation dynamics, rather than providing a detailed model of auxin transport in the cell, we shall hypothesize specific, but qualitatively reasonable in the sense described below, time-independent spatially varying expressions for the auxin distribution in a 2-D azimuthal projection of a 3-D RH cell wall. In this way, we capture crucial ingredients of the spatially-dependent auxin distribution in the cell. Qualitatively, since a RH cell is an epidermal cell, influx and efflux biochemical gates are distributed along the cell membrane, and the auxin flux is known to be primarily directed from the tip of the cell towards the elongation zone (cf. [16, 21]), leading to a longitudinally spatially decreasing auxin distribution. Our specific form for the auxin distribution will incorporate such a longitudinal spatial dependence. However, as an extension to the auxin model developed in [4, 6], here we also can take into account non-RH lateral contributions of auxin flux in RH cells by considering an auxin distribution with a transverse spatial dependence. This assumption arises as non-RH cells are believed to promote auxin transversal flux into RH cells. However, although auxin is believed to follow pathways either through the cytoplasm or through cell walls of adjacent cells, we are primarily interested in the former pathway that is directed within the cell, and which induces the biochemical switching process (cf. [31]). In this new 2-D analysis of the RD model, we will allow for an auxin distribution that has an arbitrary spatial dependence. However, in illustrating our 2-D theory, and motivated by the biological framework discussed above, we will focus on two specific forms for the auxin distribution, as shown in Fig. 1: a 1-D form that has the spatially decreasing longitudinal dependence of [4, 6], and a modification of this 1-D form that allows for a transverse dependence. We will show that these specific auxin gradients induce an alignment of localized regions of elevated ROP concentration, referred to here as spots. The RD system under consideration exhibits spot-pinning phenomena, in the sense that the position of spots in the stationary patterns is determined by the spatial gradient of auxin. As we shall see below, spots slowly drift until they become aligned with the longitudinal axis of the cell. Their steady-state locations along this cell-midline are at locations determined by the auxin gradient, which effectively “pins” the spots. A similar pinning phenomenon is found for spatially localized structures in other systems: an example is the vortex-pinning phenomenon in the Ginzburg–Landau model of superconductivity in an heterogeneous medium, see [1].

In Fig. 2, a sketch of an idealized 3-D RH cell is shown. In this figure, heavily dashed gray lines represent the RH cell membrane. The auxin flux is represented by black arrows, which schematically illustrates a longitudinal decreasing auxin pathway through adjacent cells. As ROP activation is assumed to occur near the cell wall and transversal curvature of the domain is not included in our model, we consider a projection onto a 2-D rectangular domain of this cell, as is also shown in Fig. 2.

In [6] we considered the 1-D version of the RD model given in Eq. 1.2, where the RH cell is assumed to be slim and flat, and where the auxin distribution depends only on the longitudinal spatial direction. There it was shown that 1-D localized steady-state patterns, representing active-ROP proteins, slowly drift towards the interior of the domain and asymmetrically get pinned to specific spatial locations by the auxin gradient. This pinning effect, induced by the spatial distribution of auxin, has been experimentally observed in, for instance, [13, 15]. Biologically, although multiple swelling formation in transgenic RH cells may occur, not all of these will undergo the transition to tip growth (cf. [22]). In [6] a linear stability analysis showed that multiple-active-ROP 1-D spikes may be linearly stable on $\mathcal{O}(1)$ time-scales to 1-D perturbations, with the stability properties depending on an inversely proportional relationship between length and auxin catalytic strength. This dynamic phenomenon is a consequence of an auxin catalytic process sustaining cell wall softening as the RH cell grows. For further discussion of this correlation between growth and biochemical catalysis see [6].

In [4] the linear stability of 1-D stripe patterns of elevated ROP concentration to 2-D spatially transverse perturbations was analyzed. There it was shown that single interior stripes and multi-stripe patterns are unstable under small transverse perturbations. The loss of stability of the stripe pattern was shown to lead to more intricate patterns consisting of spots, or a boundary stripe together with multiple spots (see Figure 9 of [4]). Moreover, it was shown that a single boundary stripe can also lose stability and lead to spot formation in the small diffusivity ratio limit.

Our main new goal in this paper is to analyze these 2-D spot patterns, which results from the breakup of a stripe, and characterize their slow evolution towards a true steady-state of the RD system. In particular, we will examine how specific forms of the auxin gradient lead to a steady-state spot pattern where the spots are aligned with the direction of the auxin gradient.

In the RD model below, $U(\mathbf{X},T)$ and $V(\mathbf{X},T)$ denote active- and inactive-ROP densities, respectively, at time $T>0$ and point $\mathbf{X}=(X,Y)\in\tilde{\Omega}\subset\mathbb{R}^{2}$ , where $\tilde{\Omega}$ is the rectangular domain $\tilde{\Omega}=[0,L_{x}]\times[0,L_{y}]$ . The on-and-off switching interaction is assumed to take place on the cell membrane, which is coupled through a standard Fickian-type diffusive process for both densities (cf. [6, 33]). Although RH cells are flanked by non-RH cells (cf. [17]), there is in general no exchange of ROPs between them. Therefore, no-flux boundary conditions are assumed on $\partial\tilde{\Omega}$ . On the other hand, we will consider two specific forms for the spatially-dependent dimensionless auxin distribution $\alpha$ : (i) a steady monotone decreasing longitudinal gradient, which is constant in the transverse direction, and (ii) a steady monotone decreasing longitudinal gradient, which decreases on either side of the midline $Y={L_{y}/2}$ ; see the left- and right-hand panels in Fig. 1, respectively. These two forms are modelled by

[TABLE]

As formulated in [4] and [6], the basic dimensionless RD model is

[TABLE]

with homogeneous Neumann boundary conditions, and where $\Omega=[0,1]\times[0,d_{y}]$ . Here $\Delta$ is the 2-D Laplacian, and the domain aspect ratio $d_{y}\equiv{L_{y}/L_{x}}$ is assumed to satisfy $0<d_{y}<1$ , which represents a cell as is shown in Fig. 2. The other dimensionless parameters in Eq. 1.2 are defined in terms of the original parameters, as in [4, 6], by

[TABLE]

In the original dimensional model of [6], $D_{1}\ll D_{2}$ are the diffusivities for $U$ and $V$ , $b$ is the rate of production of inactive ROP (source parameter), $c$ is the rate constant for deactivation, $r$ is the rate constant describing active ROPs being used up in cell softening and subsequent hair formation (loss parameter), and the activation step is assumed to be proportional to $k_{1}V+k_{20}\alpha(\mathbf{x})U^{2}V$ . The activation and overall auxin level within the cell, initiating an autocatalytic reaction, is represented by $k_{1}>0$ and $k_{20}>0$ respectively. The parameter $k_{20}$ , arising in the bifurcation parameter $\gamma$ of Eq. 1.3b, will play a key role for pattern formation (see [6, 33] for more details on the model formulation).

We will examine the role that geometrical features, such as the 2-D domain and the auxin-gradient, play on the dynamics of localized regions, referred to as spots or patches, where the active-ROP has an elevated concentration. This analysis is an extension to the analysis performed in [4], where it was shown that a 1-D localized stripe pattern of active-ROP in a 2-D domain will, generically, exhibit breakup into spatially localized spots. To analyze the subsequent dynamics and stability of these localized spots in the presence of the auxin gradient, we will extend the hybrid asymptotic-numerical methodology developed in [25, 35] for prototypical RD systems with spatially homogeneous coefficients. This analysis will lead to a novel finite-dimensional dynamical system characterizing slow spot evolution. Although in our numerical simulations we will focus on the two specific forms for the auxin gradient in Eq. 1.1, our analysis applies to an arbitrary gradient.

The outline of the paper is as follows. In Section 2 we introduce the basic scaling for Eq. 1.2, and we use singular perturbation methods to construct an $N$ -spot quasi steady-state solution. In Section 3 we derive a differential algebraic system (DAE) of ODEs coupled to a nonlinear algebraic system, which describes the slow dynamics of the centers of a collection of localized spots. In particular, we explore the role that the auxin gradient (i) in Eq. 1.1 has on the ultimate spatial alignment of the localized spots that result from the breakup of an interior stripe. For the auxin gradient (i) of Eq. 1.1, in Section 4 we study the onset of spot self-replication instabilities and other ${\mathcal{O}}(1)$ time-scale instabilities of quasi steady-state solutions. In Section 5 we examine localized spot patterns for the more biologically realistic model auxin (ii) of Eq. 1.1. To illustrate the theory, throughout the paper we perform various numerical simulations and numerical bifurcation analyses using the parameter sets in Table 1, which are qualitatively close to the biological parameters (see [6, 33]). Finally, a brief discussion is given in Section 6.

2 Asymptotic Regime for Spot Formation

We shall assume a shallow, oblong cell, which will be modelled as a 2-D flat rectangular domain, as shown in Fig. 2. In this domain, the biochemical interactions leading to an RH initiation process are assumed to be governed by the RD model Eq. 1.2.

For a RD system, a spatial pattern in 2-D consisting of spots is understood as a collection of structures where, at least, one of its solution components is highly localized. These structures typically evolve in such a manner that the spatial locations of the spots change slowly in time. However, depending on the parameter regime, these localized patterns can also exhibit fast ${\mathcal{O}}(1)$ time-scale instabilities, leading either to spot creation or destruction. For prototypical RD systems such as the Gierer–Meinhardt, Gray–Scott, Schnakenberg, and Brusselator models, the slow dynamics of localized solutions and the possibility of self-replication and competition instabilities, leading either to spot creation or destruction, respectively, have been studied using a hybrid asymptotic-numerical approach in [26, 35]. In addition, in the large inhibitor diffusion limit, a leading-order-in- ${-1/\log\varepsilon}$ analysis, shows that the linear stability problem for localized spot patterns characterizing competition instabilities can be reduced to the study of classes of nonlocal eigenvalue problems (NLEPs) (cf. [37, 38]). In a 1-D context, rigorous results for the spectral properties of “far-from-equilibrium” periodic 1-D patterns have been obtained recently in [11].

Our previous studies in [6] of 1-D spike patterns and in [4] of 1-D stripe patterns have extended these previous NLEP analyses of prototypical RD systems to the more complicated, and biologically realistic, plant root-hair system Eq. 1.2. In the analysis below we extend the 2-D hybrid asymptotic-numerical methodology developed in [26, 35] for prototypical RD systems to Eq. 1.2. The primary new feature in Eq. 1.2, not considered in these previous works, is to analyze the effect of a spatial gradient in the right-hand sides of Eq. 1.2, which represents variations in the auxin distribution. This spatial gradient is central to the spatial alignment of the localized structures.

To initiate our hybrid approach for Eq. 1.2, we first need to rescale the RD model Eq. 1.2 in $\Omega$ in order to determine the parameter regime where localized spots exist. To do so, we integrate the steady-state version of the RD model Eq. 1.2 over the domain to show, as in Proposition 2.2 of [6] for the 1-D problem, that for any steady-state solution $U_{0}(\mathbf{x})$ of Eq. 1.2 we must have

[TABLE]

where $|\Omega|=d_{y}$ is the area of $\Omega$ . Since the parameters $d_{y}$ , $\beta$ and $\gamma$ are $\mathcal{O}(1)$ , the constraint Eq. 2.4 implies that the average of $U$ across the whole domain is also ${\mathcal{O}}(1)$ . In other words, if we seek localized solutions such that $U\to 0$ away from localized ${\mathcal{O}}(\varepsilon)$ regions near a collection of spots, then from Eq. 2.4 we must have that $U=\mathcal{O}\left(\varepsilon^{-2}\right)$ near each spot.

2.1 A Multiple-Spot Quasi Steady-State Pattern

In order to derive a finite-dimensional dynamical system for the slow dynamics of active-ROP spots, we must first construct a quasi steady-state pattern describing ROP aggregation in $N$ small distinct spatial regions.

To do so, we seek a quasi steady-state solution where $U$ is localized within an $\mathcal{O}\left(\varepsilon\right)$ region near each spot location $\mathbf{x}_{j}$ for $j=1,\ldots,N$ . From the conservation law Eq. 2.4 we get that $U=\mathcal{O}\left(\varepsilon^{-2}\right)$ and, as a consequence, $V=\mathcal{O}(\varepsilon^{2})$ near each spot. In this way, if we replace $U=\varepsilon^{-2}U_{j}$ near each spot, and define $\boldsymbol{\xi}=\varepsilon^{-1}\left(\mathbf{x}-\mathbf{x}_{j}\right)$ , we obtain from Eq. 2.4 that $\sum_{j=1}^{N}\int_{\mathbb{R}^{2}}U_{j}\;d\boldsymbol{\xi}\sim{d_{y}/(\beta\gamma)}$ . This scaling law motivates the introduction of the new variables $u$ and $v$ defined by

[TABLE]

In terms of these new variables, Eq. 1.2 with $\partial_{n}u=\partial_{n}v=0$ on $\mathbf{x}\in\partial\Omega$ becomes

[TABLE]

where $D=\varepsilon^{-2}D_{0}$ comes from balancing terms in Eq. 2.6b.

In the inner region near the $j$ -th spot, the leading-order terms in the inner expansion are locally radially symmetric. We expand this inner solution as

[TABLE]

Substituting Eq. 2.7 into the steady-state problem for Eq. 2.6, we obtain the following leading-order radially symmetric problem on $0<\rho<\infty$ :

[TABLE]

where $\Delta_{\rho}\equiv\partial_{\rho\rho}+\rho^{-1}\partial_{\rho}$ is the Laplacian operator in polar co-ordinates. As a remark, in our derivation of the spot dynamics in Section 3 we will need to retain the next term in the Taylor series of the auxin distribution, representing the auxin gradient in the regions where active-ROP is localized. This is given by

[TABLE]

These higher-order terms are key to deriving a dynamical system for spot dynamics.

We consider Eq. 2.8 together with the boundary conditions

[TABLE]

where $S_{j}$ is called the source parameter. The logarithmic far-field behavior for $v_{0j}$ arises since, owing to the fact that there is no $-v_{0j}$ term in the second equation of (2.8), we cannot assume that $v_{0j}$ is bounded as $\rho\to\infty$ . The correction term $\chi_{j}=\chi(S_{j})$ in the far-field behavior is determined from

[TABLE]

By integrating the equation for $v_{0j}$ in Eq. 2.8, and using the limiting behavior $v_{0j}\sim S_{j}\log\rho$ as $\rho\to\infty$ , we obtain the identity

[TABLE]

In Eq. 2.8, we introduce the change of variables

[TABLE]

to obtain what we refer to as the shape canonical core problem (SCCP), which consists of

[TABLE]

Upon substituting (2.13) into the identity (2.12), we obtain

[TABLE]

To numerically compute solutions to Eq. 2.14 we use the MATLAB code BVP4C and solve the resulting BVP on the large, but finite, interval $[0,\rho_{0}]$ , for a range of values of the parameter $S_{cj}$ . The constant $\chi_{c}\left(S_{cj}\right)$ is identified from $v_{c}\left(\rho_{0}\right)-S_{c}\log\rho_{0}=\chi_{c}$ , where $v_{c}$ is the numerically computed solution, and so it depends only on $S_{cj}$ and the specified ratio ${\tau/\beta}$ in Eq. 2.14. We solve the boundary-value problem on $0<\rho_{0}<12$ with a tolerance of $10^{-4}$ . Increasing $\rho_{0}$ did not change the computational results shown in Fig. 3. The results of our computation are shown in Fig. 3 where we plot $\chi_{c}\left(S_{cj}\right)$ versus $S_{cj}$ . In Fig. 3 we also plot $u_{c}$ and $v_{c}$ at a few values of $S_{cj}$ . We observe that $u_{c}(\rho)$ has a volcano-shaped profile when $S_{cj}=5$ , which suggests, by means of the identity Eq. 2.12, that this parameter value should be near where spot self-replication will occur (see Section 4.1).

We remark that the SCCP Eq. 2.14 depends only on the ratio ${\tau/\beta}$ , which characterizes the deactivation and removal rates of active-ROP (see Eq. 1.3a). This implies that the SCCP Eq. 2.14 only describes the aggregation process. On the other hand, notice that Eq. 2.13 and Eq. 2.15 reveal the role that the auxin plays on the activation process. Since $\alpha(\mathbf{x})$ decreases longitudinally for any of the two forms in Eq. 1.1, the scaling Eq. 2.13 is such that the further from the left boundary active-ROP localizes, the larger will be the solution amplitude. This feature will be confirmed by numerical simulations in Section 2.2. In addition, the scaling of the source-parameter in Eq. 2.15 indicates that the auxin controls the inactive-ROP distribution along the cell. That is, while $S_{cj}$ characterizes how $v$ interacts with $u$ in the region where a spot occurs, the parameter $S_{j}$ will be determined by the inhomogeneous distribution of source points, which are also governed by the auxin gradient. In other words, the gradient controls and catalyzes the switching process, which leads to localized elevations of concentration of active-ROP.

Next, we examine the outer region away from the spots centered at $\mathbf{x}_{1},\ldots,\mathbf{x}_{N}$ , which are locations where $u$ is localized. To leading-order, the steady-state problem in Eq. 2.6a yields that $u_{0}\sim{\varepsilon^{4}v_{0}/(\beta\gamma)}$ as obtained from balancing terms. Therefore, from this order in the reaction terms of Eq. 2.6b, we obtain in the outer region that

[TABLE]

Next, we calculate the effect of the localized spots on the outer equation for $v$ . Upon assuming $u\sim u_{j}$ and $v\sim v_{j}$ near the $j$ -th spot, we approximate the singular terms in the sense of distributions as Dirac masses. Indeed, labelling the reactive terms by $\mathcal{R}_{j}\equiv\tau\gamma\left(\alpha\left(\mathbf{x}\right)u_{j}^{2}v_{j}-u_{j}\right)+\beta\gamma u_{j}$ , we use (2.12) and the Mean Value Theorem to calculate, in the sense of distributions, that

[TABLE]

as $\varepsilon\to 0$ , where $\mathcal{R}_{jj}(\rho)\equiv\tau\gamma\left(\alpha\left(\mathbf{x}_{j}\right)u_{0j}^{2}v_{0j}-u_{0j}\right)+\beta\gamma u_{0j}$ .

Combining Eq. 2.16, Eq. 2.17 and Eq. 2.6b we obtain that the leading-order outer solution $v_{0}$ , in the region away from ${\mathcal{O}}(\varepsilon)$ neighborhoods of the spots, satisfies

[TABLE]

We now proceed to solve Eq. 2.18 for $v_{0}$ and derive a nonlinear algebraic system for the source parameters $S_{c1},\ldots,S_{cN}$ . To do so, we introduce the unique Neumann Green’s function $G(\mathbf{x};\mathbf{y})$ satisfying (cf. [26])

[TABLE]

This Green’s function satisfies the reciprocity relation $G\left(\mathbf{x};\mathbf{y}\right)=G\left(\mathbf{y};\mathbf{x}\right)$ , and can be globally decomposed as $G\left(\mathbf{x};\mathbf{y}\right)=-1/(2\pi)\log\left|\mathbf{x}-\mathbf{y}\right|+R\left(\mathbf{x};\mathbf{y}\right)$ , where $R\left(\mathbf{x};\mathbf{y}\right)$ is the smooth regular part of $G\left(\mathbf{x};\mathbf{y}\right)$ . For $\mathbf{x}\to\mathbf{y}$ , this Green’s function has the local behavior

[TABLE]

An explicit infinite series solution for $G$ and $R$ in the rectangle $\Omega$ is available (cf. [26]).

In terms of this Green’s function, the solution to Eq. 2.18a is

[TABLE]

Here $\bar{v}_{0}$ is a constant to be found. The condition on $S_{i}$ in Eq. 2.19 arises from the Divergence Theorem applied to Eq. 2.18. We then expand $v_{0}$ as $\mathbf{x}\to\mathbf{x}_{j}$ and enforce the required singularity behavior Eq. 2.18b. This yields that

[TABLE]

Here $G_{ji}\equiv G\left(\mathbf{x}_{j};\mathbf{x}_{i}\right)$ and $R_{j}\equiv R\left(\mathbf{x}_{j};\mathbf{x}_{j}\right)$ . The system Eq. 2.20, together with the constraint in Eq. 2.19, defines a nonlinear algebraic system for the $N+1$ unknowns $S_{1},\ldots,S_{N}$ and $\bar{v}_{0}$ . To write this system in matrix form, we define $\mathbf{e}\equiv\left(1,\ldots,1\right)^{T}$ , $\mathbf{S}\equiv(S_{1},\ldots,S_{N})^{T}$ , $\boldsymbol{\chi}\equiv(\chi(S_{1}),\ldots,\chi(S_{N})^{T}$ , $E\equiv{\mathbf{e}\mathbf{e}^{T}/N}$ , and the symmetric Green’s matrix $\mathbf{G}$ with entries

[TABLE]

Then, Eq. 2.20, together with the constraint in Eq. 2.19, can be written as

[TABLE]

Left multiplying the first equation in Eq. 2.21 by $\mathbf{e}^{T}$ , and using the second equation in Eq. 2.21, we find an expression for $\bar{v}_{0}$

[TABLE]

Substituting $\bar{v}_{0}$ back into the first equation of Eq. 2.21, we obtain a nonlinear system for $S_{1},\ldots,S_{N}$ given by

[TABLE]

In terms of $\mathbf{S}$ , $\bar{v}_{0}$ is given in Eq. 2.22. To relate $S_{j}$ and $\chi_{j}$ with the SCCP Eq. 2.14, we define

[TABLE]

Then, by using Eq. 2.15, we find

[TABLE]

Finally, upon substituting Eq. 2.24 into Eq. 2.23 we obtain our main result regarding quasi steady-state spot patterns.

Proposition 2.1.

Let $0<\varepsilon\ll 1$ and assume $D_{0}={\mathcal{O}}(1)$ in Eq. 2.6 so that $D=\mathcal{O}\left(\varepsilon^{-2}\right)$ in Eq. 1.2. Suppose that the nonlinear algebraic system for the source parameters, given by

[TABLE]

has a solution for the given spatial configuration $\mathbf{x}_{1},\ldots,\mathbf{x}_{N}$ of spots. Then, there is a quasi steady-state solution for Eq. 1.2 with $U=\varepsilon^{-2}u$ and $V=\varepsilon^{2}v$ . In the outer region away from the spots, we have $u={\mathcal{O}}(\varepsilon^{4})$ and $v\sim v_{0}$ , where $v_{0}$ is given asymptotically in terms of $\mathbf{S}=\upomega\boldsymbol{\upalpha}\mathbf{S}_{c}$ by Eq. 2.19. In contrast, in the $j$ -th inner region near $\mathbf{x}_{j}$ , we have

[TABLE]

*where $u_{c}$ and $v_{c}$ is the solution of the SCCP Eq. 2.14 in terms of $S_{cj}$ satisfying Eq. 2.25a. *

2.2 Spots and Baby Droplets

Proposition 2.1 provides an asymptotic characterization of a $N$ -spot quasi steady-state solution. As has been previously studied in the previous paper [4], localized spots can be generated by the breakup of a straight homoclinic stripe. The numerical simulations in Fig. 4, for the 1-D auxin gradient in Eq. 1.1, illustrate how small spots are formed from the breakup of an interior homoclinic stripe and how they evolve slowly in time to a steady-state two-spot pattern aligned with the direction of the auxin gradient. The initial condition for these simulations is obtained by trivially extending a 1-D spike in the transversal direction to obtain a homoclinic stripe, which is a steady-state of Eq. 1.2, and then perturbing it with a transverse perturbation of small amplitude. As shown in Fig. 4–Fig. 4, the amplitude of the left-most spot decreases as the spot drifts towards the left domain boundary at $X=0$ , while the amplitude of the right-most spot grows as it drifts towards the right domain boundary at $X=L_{x}$ . The spot profiles are described by Proposition 2.1 for each fixed re-scaled variables $\mathbf{x}_{j}$ . Our realizations of the asymptotic theory shown below in Fig. 7 confirm the numerical results in Fig. 4 that the spots become aligned with the direction of the auxin gradient.

In order to explore solutions that bifurcate from a multi-spot steady state, we adapted the MATLAB code of [2, 34] and performed a numerical continuation for the RD system Eq. 1.2 in the bifurcation parameter $k_{20}$ , starting from the solution in Fig. 4 ( all other parameters as in Parameter Set A in Table 1). The results are shown in Fig. 5. In the bifurcation diagram depicted in Fig. 5 the linearly stable branch is plotted as a solid curve and unstable ones by light gray-dashed curves. We found fold bifurcations and a pitchfork bifurcation, represented by a filled circle and a filled star, respectively. To the right of the pitchfork bifurcation, the symmetric solution branch destabilises; the emerging branch of asymmetric states is fully unstable, and features a sequence of fold bifurcations. Overall, these results suggest that a snaking bifurcation structure occurs in 2-D domains, similar to what was found in [4] for the corresponding 1-D system. This intricate branch features a set of fold points (not explicitly labelled in Fig. 5), and around each of them a new spot emerges. A pattern consisting of a spot and a peanut structure is found in the solution branch labelled as P. This solution is shown in Fig. 5.

We found that solutions with peanut structures or baby droplets are unstable steady-states in this region of parameter space, and so they do not persist in time; see a sample of these in Fig. 5. In addition, as can be seen in Fig. 5, there are stable multi-spot solutions when $k_{20}$ is small enough. Overall, this suggests that certain $\mathcal{O}(1)$ time-scale instabilities play an important role in transitioning between steady-state patterns (cf. [8, 20]). Further, in Section 4, peanut-splitting instabilities are addressed.

3 The Slow Dynamics of Spots

In Section 2.1 we constructed an $N$ -spot quasi steady-state solution for an arbitrary, but “frozen”, spatial configuration $\mathbf{x}_{1},\ldots,\mathbf{x}_{N}$ of localized patches of active-ROP. In this section, we derive a differential-algebraic system of ODEs (DAE) for the slow time-evolution of $\mathbf{x}_{1},\ldots,\mathbf{x}_{N}$ towards their equilibrium locations in $\Omega$ . This system will consist of a constraint that involves the slow time evolution of the source parameters $S_{c1},\ldots,S_{cN}$ . The derivation of this DAE system is an extension of the 1-D analysis of [6].

To derive the DAE system, we must extend the analysis in Section 2.1 by retaining the ${\mathcal{O}}(\varepsilon)$ terms in Eq. 2.7 and Eq. 2.9. By allowing the spot locations to depend on a distinguished limit for the slow time scale $\eta=\varepsilon^{2}t$ as $\mathbf{x}_{j}=\mathbf{x}_{j}\left(\eta\right)$ , we obtain in the $j$ -th inner region that the corrections terms $u_{1j}$ and $v_{1j}$ to the SCCP Eq. 2.14 satisfy

[TABLE]

Here $\nabla_{\mathbf{x}}R_{j}\equiv\nabla_{\mathbf{x}}R\left(\mathbf{x}_{j};\mathbf{x}_{j}\right)$ and $\nabla_{\mathbf{x}}G_{j}\equiv\nabla_{\mathbf{x}}G\left(\mathbf{x}_{j};\mathbf{x}_{j}\right)$ for $j=1,\ldots,N$ .

To impose a solvability condition on Eq. 3.27, which will lead to ODEs for the spot locations, we first rewrite Eq. 3.27 in terms of the canonical variables Eq. 2.13 of the SCCP Eq. 2.14. Using this transformation, the right-hand side $\mathbf{F}$ of Eq. 3.27g becomes

[TABLE]

In addition, the matrix $\mathbf{M}$ on the left-hand side of Eq. 3.27g becomes

[TABLE]

Together with Eq. 3.28 this suggests that we define new variables $\hat{u}_{1}$ and $\hat{v}_{1}$ by

[TABLE]

In terms of these new variables, and upon substituting Eq. 3.28 and Eq. 3.32 into Eq. 3.27, we obtain that $\hat{\mathbf{Q}}_{1}\equiv(\hat{u}_{1},\hat{v}_{1})^{T}$ in $\boldsymbol{\xi}\in\mathbb{R}^{2}$ satisfies

[TABLE]

where $u_{c}$ and $v_{c}$ satisfy the SCCP Eq. 2.14, which is parametrized by $S_{cj}$ .

In Eq. 3.33k, the source parameters $S_{cj}$ depend indirectly on the auxin distribution through the nonlinear algebraic system Eq. 2.25, as summarized in Proposition 2.1. We further observe from the second term on the right-hand side of Eq. 3.33g that its coefficient $\nabla\alpha_{j}/\alpha_{j}$ is independent of the magnitude of the spatially inhomogeneous distribution of auxin, but instead depends on the direction of the gradient.

3.1 Solvability Condition

We now derive a DAE system for the spot dynamics by applying a solvability condition to Eq. 3.33. We label $\boldsymbol{\xi}=(\xi_{1},\xi_{2})^{T}$ and ${\mathbf{u}}_{c}\equiv(u_{c},v_{c})^{T}$ . Differentiating the core problem Eq. 2.14a with respect to the $\xi_{i}$ , we get

[TABLE]

This demonstrates that the dimension of the nullspace of ${\mathcal{L}}$ in Eq. 3.33, and hence its adjoint ${\mathcal{L}}^{\star}$ , is at least two-dimensional. Numerically, for our two parameter sets in Table 1, we have checked that the nullspace of ${\mathcal{L}}$ is exactly two-dimensional, provided that $S_{j}$ does not coincide with the spot self-replication threshold given in Section 4.1.

There are two independent nontrivial solutions to the homogeneous adjoint problem ${\mathcal{L}}^{\star}\mathbf{\Psi}\equiv\Delta\mathbf{\Psi}+\mathbf{M}_{c}^{T}\mathbf{\Psi}=\mathbf{0}$ given by $\mathbf{\Psi}_{i}\equiv{\mathbf{P}}(\rho){\xi_{i}/\rho}$ for $i=1,2$ , where $\mathbf{P}$ satisfies

[TABLE]

Here the condition at infinity in Eq. 3.34 is used as a normalization condition for $\mathbf{P}$ .

Next, to derive our solvability condition we use Green’s identity over a large disk $\Omega_{\rho_{0}}$ of radius $|\boldsymbol{\xi}|=\rho_{0}\gg 1$ to obtain for $i=1,2$ that

[TABLE]

With $\mathbf{\Psi}_{i}\equiv\mathbf{P}(\rho){y_{i}/\rho}$ , for $i\in\{{1,2\}}$ , we first calculate the left-hand side (LHS) of (3.35) using Eq. 3.33g to obtain with $\mathbf{P}\equiv(P_{1}(\rho),P_{2}(\rho))^{T}$ and $\mathbf{x}_{j}^{\prime}\equiv(x_{j1}^{\prime},x_{j2}^{\prime})^{T}$ that

[TABLE]

where $\mathbf{x}^{\prime}_{j}=d\mathbf{x}_{j}/d\eta$ . In deriving the result above, we used the identity $\int_{\Omega_{\rho_{0}}}\xi_{i}\xi_{k}f(\rho)\,d\boldsymbol{\xi}=\delta_{ik}\pi\int_{0}^{\rho_{0}}\rho^{3}f(\rho)\,d\rho$ , for any radially symmetric function $f(\rho)$ , where $\delta_{ik}$ is the Kronecker delta.

Next, we calculate the right-hand side (RHS) of Eq. 3.35 using the far-field behaviors of $\hat{\mathbf{Q}}_{1}$ and $\mathbf{P}$ as $\rho\to\infty$ . We derive that

[TABLE]

In the last passage, we used $dS=\rho_{0}d\theta$ where $\theta$ is the polar angle. Finally, we equate LHS and RHS for $i=1,2$ , and write the resulting expression in vector form:

[TABLE]

We summarize our main, formally-derived, asymptotic result for slow spot dynamics as follows:

Proposition 3.1.

Under the same assumptions as Proposition 2.1, and assuming that the $N$ -spot quasi steady-state solution is stable on an ${\mathcal{O}}(1)$ time-scale, the slow dynamics on the long time-scale $\eta=\varepsilon^{2}t$ of this quasi steady-state spot pattern consists of the constraints Eq. 2.25 coupled to the following ODEs for $j=1,\ldots,N$ :

[TABLE]

*In Eq. 3.37a, the source parameters $S_{cj}$ satisfy the nonlinear algebraic system Eq. 2.25, which depends on the instantaneous spatial configuration of the $N$ spots. Overall, this leads to a DAE system characterizing slow spot dynamics. *

Proposition 3.1 describes the slow dynamics of a collection of $N$ localized spots under an arbitrary, but smooth, spatially-dependent auxin gradient. It is an extension of the 1-D analysis of spike evolution, considered in [6]. The dynamics in Eq. 3.37a, shows that the spot locations depend on the gradient of the Green’s function, which depends on the domain $\Omega$ , as well as the spatial gradient of the auxin distribution. In particular, the spot dynamics depends only indirectly on the magnitude of the auxin distribution $\alpha(\mathbf{x}_{j})$ through the source parameter $S_{cj}$ . The auxin gradient $\nabla\alpha_{j}$ , however, is essential to determining the true steady-state spatial configuration of spots. In addition, the spatial interaction between the spots arises from the terms in the finite sum of Eq. 3.37a, mediated by the Green’s function. Since the Green’s function and its regular part can be found analytically for a rectangular domain (cf. [26]), we can readily use Eq. 3.37 to numerically track the slow time-evolution of a collection of spots for a specified auxin gradient.

Before illustrating results from the DAE dynamics, we must determine $n_{1}$ and $n_{2}$ as a function of $S_{cj}$ for a prescribed ratio ${\tau/\beta}$ . This ratio is associated with the linear terms in the kinetics of the original dimensional system Eq. 1.2, which are related to the deactivation of ROPs and production of other biochemical complexes which promote cell wall softening (cf. [6, 33]). To determine $n_{1}$ and $n_{2}$ , we first solve the adjoint problem Eq. 3.34 numerically using the MATLAB routine BVP4C. This is done by enforcing the local behavior that $\mathbf{P}={\mathcal{O}}(\rho)$ as $\rho\to 0$ and by imposing the far-field behavior for $\mathbf{P}$ , given in Eq. 3.34, at $\rho=\rho_{0}=12$ . In Fig. 6 we plot $P_{1}$ and $P_{2}$ for three values of $S_{cj}$ , where ${\tau/\beta}=3$ , and Parameter Set A in Table 1 was used. For each of the three values of $S_{cj}$ , we observe that the far-field behavior $P_{2}(\rho)\sim{1/\rho}$ and $P_{1}(\rho)\sim\left(\tau/\beta-1\right)/\rho$ as $\rho\to\infty$ , which is readily derived from Eq. 3.34, is indeed satisfied. Upon performing the required quadratures in Eq. 3.37b, in Fig. 6 we plot $n_{1}$ and $n_{2}$ versus $S_{cj}$ . These numerical results show that $n_{1}>0$ and $n_{2}>0$ for $S_{cj}>0$ , which will ensure existence of stable fixed points of the DAE dynamics. These stable fixed points correspond to realizable steady-state spot configurations for our two specific forms for the auxin gradient.

3.2 Comparison Between Theory and Asymptotics for Slow Spot Dynamics

In this subsection we compare predictions from our asymptotic theory for slow spot dynamics with corresponding full numerical results computed from Eq. 2.6 using a spatial mesh with 500 and 140 gridpoints in the $x$ and $y$ directions, respectively. For the time-stepping a modified Rosenbrock formula of order 2 was used.

The procedure to obtain numerical results from our asymptotic DAE system in Proposition 3.1 is as follows. This DAE system is solved numerically by using Newton’s method to solve the nonlinear algebraic system Eq. 2.25 together with a Runge–Kutta ODE solver to evolve the dynamics in Eq. 3.37a. The solvability integrals $n_{1}(S_{cj})$ and $n_{2}(S_{cj})$ in the dynamics Eq. 3.37a and the function $\boldsymbol{\chi}_{c}$ in the nonlinear algebraic system Eq. 2.25 are pre-computed at $200$ grid points in $S_{cj}$ and a cubic spline interpolation is fitted to the discretely sampled functions to compute them at arbitrary values of $S_{cj}$ . For the rectangular domain, explicit expressions for the Green’s functions $G$ and $R$ , together with their gradients, as required in the DAE system, are calculated from the expressions in §4 of [26].

To compare our results for slow spot dynamics we use Parameter Set A of Table 1 and take $k_{20}=0.3$ . For the auxin gradient we took the monotone 1-D gradient in type (i) in Eq. 1.1. By a small perturbation of the unstable 1-D stripe solution, our full numerical computations of Eq. 2.6 lead to the creation of two localized spots. Numerical values for the centers of the two spots are calculated (see the caption of Fig. 7) and these values are used as the initial conditions for our numerical solution of the DAE system in Proposition 3.1. In Fig. 7 we compare our full numerical results for the $x$ and $y$ coordinates for the spot trajectories, as computed from the RD system Eq. 2.6, and those computed from the corresponding DAE system. The key distinguishing feature in these results is that the spots become aligned to the longitudinal midline of the domain. We observe that the $y$ -components of the spot trajectories are predicted rather accurately over long time intervals by the DAE system. However, although the $x$ -components of the trajectories are initially close, they deviate somewhat as $t$ increases. Since we only have a formal asymptotic theory, with no error bounds, and because it is computationally expensive to perform more refined numerical simulations over very-long time scales from the full RD system, we are not able to precisely identify why the agreement in the $y$ -coordinate is better than for the $x$ -coordinate. However, overall, the asymptotic theory accurately identifies the time-scale over which the two localized spots become aligned with the auxin gradient along the mid-line of the cell.

4 Fast $\mathcal{O}(1)$ Time-Scale Instabilities

We now briefly examine the stability properties of the $N$ -spot quasi-equilibrium solution of Proposition 2.1 to ${\mathcal{O}}(1)$ time-scale instabilities, which are fast relative to the slow dynamics. We first consider spot self-replication instabilities associated with non-radially symmetric perturbations near each spot.

4.1 The Self-Replication Threshold

Since the speed of the slow spot drift is ${\mathcal{O}}(\varepsilon^{2})\ll 1$ , in our stability analysis below we will assume that the spots are “frozen” at some configuration $\mathbf{x}_{1},\ldots,\mathbf{x}_{N}$ . We will consider the possibility of instabilities that are locally non-radially symmetric near each spot. In the inner region near the $j$ -th spot at $\mathbf{x}_{j}$ , where $\boldsymbol{\xi}=\varepsilon^{-1}(\mathbf{x}-\mathbf{x}_{j})$ , we linearize Eq. 2.6 around the leading-order core solution $u_{0j}$ , $v_{0j}$ , satisfying Eq. 2.8, by writing $u=u_{0j}+e^{\lambda t}\phi_{0j}$ and $v=v_{0j}+e^{\lambda t}\psi_{0j}$ . From Eq. 2.6, we obtain to leading order that

[TABLE]

where $\Delta_{\boldsymbol{\xi}}$ is the Laplacian in the local variable $\boldsymbol{\xi}$ . Then, upon relating $u_{0j}$ , $v_{0j}$ to the SCCP by using (2.13), and defining $\psi_{0j}\equiv{\beta\gamma\tilde{\psi}_{0j}/D_{0}}$ , the system above reduces to

[TABLE]

where $u_{c}$ , $v_{c}$ is the solution to the SCCP Eq. 2.14.

We then look for an ${\mathcal{O}}(1)$ time-scale instability associated with the local angular integer mode $m\geq 1$ by introducing the new variables $\Phi_{0}(\rho)$ and $\Psi_{0}(\rho)$ defined by

[TABLE]

and $\boldsymbol{\xi}^{T}=\rho(\cos\theta,\sin\theta)$ . Substituting Eq. 4.39 into Eq. 4.38, we obtain the eigenvalue problem:

[TABLE]

Here we have defined ${\cal L}_{m}\Upsilon\equiv\partial_{\rho\rho}\Upsilon+\rho^{-1}\partial_{\rho}\Upsilon-m^{2}\rho^{-2}\Upsilon$ . We impose the usual regularity condition for $\Phi_{0}$ and $\Psi_{0}$ at $\rho=0$ . The appropriate far-field boundary conditions for Eq. 4.40 is discussed below.

Since the eigenvalue problem Eq. 4.40 is difficult to study analytically, we solve it numerically for various integer values of $m$ . We denote $\lambda_{\textrm{max}}$ to be the eigenvalue of Eq. 4.40 with the largest real part. Since $u_{c}$ and $v_{c}$ depend on $S_{cj}$ from the SCCP Eq. 2.14, we have implicitly that $\lambda_{\textrm{max}}=\lambda_{\textrm{max}}(S_{cj},m)$ . To determine the onset of any instabilities, for each $m$ we compute the smallest threshold value $S_{cj}=\Sigma_{m}$ where $\mbox{Re}(\lambda_{\textrm{max}}(\Sigma_{m},m))=0$ . In our computations, we only consider $m=2,3,4,\ldots$ , since $\lambda_{\textrm{max}}=0$ for any value of $S_{cj}$ for the translational mode $m=1$ . Any such instability for $m=1$ is reflected in instability in the DAE system Eq. 3.37a.

For $m\geq 2$ we impose the far-field behavior that $\Phi_{0}$ decays exponentially as $\rho\to\infty$ while $\Psi_{0}\sim{\mathcal{O}}(\rho^{-m})\to 0$ as $\rho\to\infty$ . With this limiting behavior, Eq. 4.40 is discretized with centered differences on a large but finite domain. We then determine $\lambda_{\textrm{max}}(S_{cj},m)$ by computing the eigenvalues of the discretized eigenvalue problem, in matrix form. For $m\geq 2$ our computations show that $\lambda_{\textrm{max}}(S_{cj},m)$ is real and that $\lambda_{\textrm{max}}(S_{cj},m)>0$ if and only if $S_{cj}>\Sigma_{m}$ . In our computations we took $400$ meshpoints on the interval $0\leq\rho<15$ . For the ratio ${\tau/\beta}=3$ , corresponding to Parameter Set A of Table 1, the results for the threshold values $\Sigma_{m}$ for $m=2,3,4$ given in Fig. 8 are insensitive to increasing either the domain length or the number of grid points. Our main conclusion is that as $S_{cj}$ is increased, the solution profile of the $j$ -th spot first becomes unstable to a non-radially symmetric peanut-splitting mode, corresponding to $m=2$ , as $S_{cj}$ increases above the threshold $\Sigma_{2}\approx 4.16$ when ${\tau/\beta}=3$ . As a remark, if we take ${\tau/\beta}=11$ , corresponding to Parameter Set B of Table 1, we compute instead that $\Sigma_{2}\approx 3.96$ .

To illustrate this peanut-splitting threshold, we consider a pattern with a single localized spot. Using Eq. 2.19 we find $S_{1}={d_{y}/(2\pi D_{0})}$ . Then, Eq. 2.15 yields the following expression for the source parameter for the SCCP Eq. 2.14

[TABLE]

When ${\tau/\beta}=3$ , corresponding to Parameter Set A of Table 1, we predict that the one-spot pattern will first undergo a shape-deformation, due to the $m=2$ peanut-splitting mode, when $S_{c1}$ exceeds the threshold $\Sigma_{2}\approx 4.16$ . Since $\gamma$ is inversely proportional to the parameter $k_{20}$ , as shown in Eq. 1.3b, which determines the auxin level, we predict that as the auxin level is increased above a threshold a peanut structure will develop from a solitary spot. As shown below in Section 4.2 for an auxin distribution of the type (i) in Eq. 1.1, this linear peanut-splitting instability triggers a nonlinear spot replication event.

In addition, in Fig. 8 we plot the source parameters corresponding to the two-spot steady-state aligned along the midline $Y=L_{y}/2$ as the parameter $k_{20}$ is varied. These asymptotic results, which asymptotically characterize states in branch A of Fig. 5, are computed from the steady-state of the DAE system Eq. 3.37. Notice that the spot closest to the left-hand boundary loses stability at $k_{20}^{*}\approx 0.851$ , when $S_{c1}$ meets the peanut-splitting threshold $\Sigma_{2}$ . On the other hand, $S_{c2}$ values corresponding to the spot closest to the right-hand boundary are always below $\Sigma_{2}$ for the parameter range considered. Hence, no shape change for this spot is observed. This result rather accurately predicts the pitchfork bifurcation in Fig. 5 portrayed by a filled star at $k_{20}^{*}\approx 0.902$ .

4.2 Numerical Illustrations of ${\mathcal{O}}(1)$ Time-Scale

Instabilities

In Fig. 5 we showed a steady-state solution consisting of a spot and a peanut structure, which is obtained from numerical continuation in the parameter $k_{20}$ . This solution belongs to the unstable branch labelled by P in Fig. 5, which is confined between two fold bifurcations (not explicitly shown). We take a solution from this branch as an initial condition in Fig. 9 for a time-dependent numerical simulation of the full RD model Eq. 1.2-Eq. 2.6 using a centered finite-difference scheme. Our numerical results in Fig. 9 show that this initial state evolves dynamically into a stable two-spot solution. This is a result of the overlapping of the stable branch A with an unstable one in the bifurcation diagram.

Other initial conditions and parameter ranges also exhibit ${\mathcal{O}}(1)$ time-scale instabilities. To illustrate these instabilities, and see whether a self-replication spot process is triggered by a linear peanut-splitting instability, we perform a direct numerical simulation of the full RD model Eq. 1.2 taking as initial condition the unstable steady-state solution labelled by B in Fig. 5. This steady-state, shown in Fig. 10, consists of two spots, with one having four small droplets associated with it. This particular steady-state is chosen since both competition and self-replication instabilities can be seen in the time evolution of this initial condition in the panels of Fig. 10. Firstly, we observe in Fig. 10 that the small droplets for the left-most spot in Fig. 10 are annihilated in $\mathcal{O}(1)$ times. Then, the right-most spot, which has a relatively low source parameter owing to its distance away from the left-hand boundary (see Eq. 2.15a), is annihilated on an ${\mathcal{O}}(1)$ time-scale triggered by a linear competition instability. The solitary remaining spot at $(X_{1},Y_{1})=(28,14.9240)$ , shown in Fig. 10, now has a source parameter $S_{c1}=9.1403$ , given in Eq. 4.41 of Section 4.1, that exceeds the peanut-splitting threshold $\Sigma_{2}\approx 4.16$ . This linear shape-deforming instability gives rise to the peanut structure shown in Fig. 10, and the direction of spot-splitting is aligned with the direction of the monotone decreasing auxin gradient, which depends only on the longitudinal direction. In Fig. 10 and Fig. 10 we show that this linear peanut-splitting instability has triggered a nonlinear spot-splitting event that creates a second localized spot. Finally, these two spots drift away from each other and their resulting slow dynamics is characterized by the DAE dynamics in Proposition 3.1.

Competition and self-replication instabilities, such as illustrated above, are two types of fast ${\mathcal{O}}(1)$ time-scale instabilities that commonly occur for localized spot patterns in singularly perturbed RD systems. Although it is beyond the scope of this paper to give a detailed analysis of a competition instability for our plant RD model Eq. 1.2, based on analogies with other RD models (cf. [9, 25, 26, 35]), this instability typically occurs when the source parameter of a particular spot is below some threshold or, equivalently, when there are too many spots that can be supported in the domain for the given substrate diffusivity. In essence, a competition instability is an overcrowding instability. Alternatively, as we have unravelled in Section 4.1, the self-replication instability is an undercrowding instability and is triggered when the source parameter of a particular spot exceeds a threshold, or equivalently when there are too few spots in the domain. For the standard Schnakenberg model with spatially homogeneous coefficients, in [26] it was shown through a center-manifold type calculation that the direction of splitting of a spot is always perpendicular to the direction of its motion. However, with an auxin gradient the direction of spot-splitting is no longer always perpendicular to the direction of motion of the spot. In our plant RD model, the auxin gradient not only enhances robustness of solutions, which results in the overlapping of solution branches and is suggestive that a strong slanted homoclinic snaking mechanism can occur (see the next section), but it also controls the location of steady-state spots (see Section 2.2 and Section 3).

5 More Robust 2D Patches and Auxin Transport

We now consider a more biologically plausible model for the auxin transport, which initiates the localization of ROP. This model allows for both a longitudinal and transverse spatial dependence of auxin. The original ROP model, derived in [6, 33] and analyzed in [4, 6] and the sections above, depends crucially on the spatial gradient of auxin. The key assumption we have made above was to assume a decreasing auxin distribution along the $X$ -direction. Indeed, such an auxin gradient controls the $x$ -coordinate spot location in such a way that the larger the overall auxin parameter $k_{20}$ , the more spots are likely to occur. However, as was discussed in Section 4, there are instabilities that occur when an extra spatial dimension is present. In other words, when an RH cell is modelled as a two-dimensional flat and oblong cell, certain pattern formation attributes become relevant that are not present in a 1-D setting. In particular, 1-D stripes generically break up into localized spots, which are then subject to possible secondary instabilities. We now explore the effect that a 2-D spatial distribution of auxin has on such localized ROP spots.

5.1 A 2-D Auxin Gradient

Next, we perform a numerical bifurcation analysis of the RD system Eq. 1.2 for an auxin distribution of the form given in (ii) of Eq. 1.1. Such a distribution represents a decreasing concentration of auxin in the $X$ -direction, as is biologically expected, but with a greater longitudinal concentration of auxin along the midline $Y={L_{y}/2}$ of the flat rectangular cell than at the edges $Y=0,L_{y}$ .

To perform a numerical bifurcation study we discretize Eq. 1.2 using centered finite differences and we adapt the 2-D continuation code written in MATLAB given in [34]. We compute branches of steady-state solutions of this system using pseudo arclength continuation, and the stability of these solutions is computed a posteriori using MATLAB eigenvalue routines. The resulting bifurcation diagram is shown in Fig. 11. We observe that there are similarities between this bifurcation diagram and the one for the 1-D problem shown in Figure 6 of [6] in that only fold-point bifurcations were found for the branches depicted. As similar to bifurcation diagrams associated with homoclinic snaking theory on a finite domain (cf. [3, 5, 7, 14, 30]), we observe a key distinguishing feature: solutions belong to a single branch of steady states, undergoing a sequence of fold bifurcations and, in some cases, a change in stability. In large intervals of $k_{20}$ we observe multistability which, in turn, indicates that hysteretic transitions between solutions with varying number of spots can occur.

To explore fine details of this bifurcation behavior, the bifurcation diagram in Fig. 11 has been split into two parts, with the top and bottom panels for lower and higher values of $k_{20}$ , respectively. From Fig. 11 we observe that a boundary spot emerges for very low values of the overall auxin rate. This corresponds to branch 1 of Fig. 11. As $k_{20}$ is increased, stability is lost at a fold-point bifurcation which gives rise to branch 2 of Fig. 11, which gathers a family of single interior spot steady-states, as shown in Fig. 11. As $k_{20}$ increases further, this interior spot solution persists until stability is lost at a fold-point bifurcation. Branch 3 in Fig. 11 consists of stable steady-states of one interior spot and one boundary spot, as shown in Fig. 11. Furthermore, Fig. 11 shows that steady-states consisting of two interior spots occur at even larger values of $k_{20}$ . This corresponds to branch 4 of Fig. 11. Finally, at much larger values of $k_{20}$ , Fig. 11 shows that an additional spot is formed at the domain boundary. This qualitative behavior associated with increasing $k_{20}$ continues, and leads to a creation-annihilation cascade similar to that observed for 1-D pulses and stripes in [6] and [4], respectively. In other words, overlapping stable branches yield steady-states consisting of interior and/or boundary spots that appear or disappear as $k_{20}$ is either slowly increased or decreased, respectively.

5.2 Instabilities with a 2-D Spatially-Dependent Auxin Gradient

Similar to the 1-D studies in [4, 6], the auxin gradient controls the location and number of 2-D localized regions of active ROP. As the level of auxin increases in the cell, an increasing number of active ROPs are formed, and their spatial locations are controlled by the spatial gradient of the auxin distribution. Moreover, in analogy with the theory of homoclinic snaking, overlapping of stable solution branches occurs, and this leads to a wide range of different observable steady-states in the RD system. These signature features of the bifurcation structure suggest that homoclinic snaking can occur for Eq. 1.2, similar to that observed in [5] for a RD system with spatially homogeneous coefficients. In passing, we note that the snaking observed in our system is slanted. Slanted snaking has been reported in systems with conserved quantities (see for instance [10, 36]). The system under consideration does not have a conserved quantity, and the slanting is caused by the auxin gradient. In other words, the gradient gives rise to a strongly slanted homoclinic snaking behavior.

Transitions from unstable to stable branches are determined through fold bifurcations, and are controlled by $\mathcal{O}(1)$ time-scale instabilities. To illustrate this behavior, we perform a direct numerical simulation of the full RD model Eq. 1.2 taking as initial condition an unstable steady-state solution consisting of two rather closely spaced spots, as shown in Fig. 12. The time-evolution of this initial condition is shown in the different panels of Fig. 12. We first observe from Fig. 12 that these two spots begin to merge at the left domain boundary. As this merging process continues, a new-born stripe emerges near the right-hand boundary, as shown in Fig. 12. This interior stripe is weakest at the top and bottom boundaries, due to the relatively low levels of auxin in these regions. In Fig. 12, the interior stripe is observed to give rise to a transversally aligned peanut structure, which remains centered at the $Y$ -midline for the same reasons as above at the top and bottom boundaries, whilst another peanut structure in the longitudinal direction occurs near the left-hand boundary. The right-most peanut structure slowly collapses to a solitary spot, while the left-most form undergoes a breakup instability, yielding two distinct spots, as shown in Fig. 12. Finally, Fig. 12 shows a steady-state solution with three localized spots that are spatially aligned along the $Y$ -midline.

In the instabilities and dynamics described above, the 2-D auxin gradient apparently play a very important role, and leads to ROP activation in the cell primarily where the auxin distribution is higher, rather than at the transversal boundaries. In particular, although a homoclinic-type transverse interior stripe is formed, it quickly collapses near the top and bottom domain boundaries where the auxin gradient is weakest. This leads to a peanut structure centered near the $Y$ -midline, which then does not undergo a self-replication process but instead leads to the merging, or aggregation, of the peanut structure into a solitary spot. The transverse component of the 2-D auxin distribution is essential to this behavior.

We remark that the asymptotic analysis in Section 2 and Section 3 for the existence and slow dynamics of 2-D quasi steady-state localized spot solutions can also readily be implemented for the 2-D auxin gradient type (ii) of Eq. 1.1. The self-replication threshold in Section 4.1 also applies to a 2-D spatial gradient for auxin. However, for an auxin gradient that depends on both the longitudinal and transverse directions, it is difficult to analytically construct a 1-D quasi steady-state stripe solution such as that done in [6]. As a result, for a 2-D auxin gradient, it is challenging to obtain a dispersion relation predicting the linear stability properties of a homoclinic stripe (see for instance [12, 24]).

6 Discussion

We have used a combination of asymptotic analysis, numerical bifurcation theory, and full time-dependent simulations of the governing RD system to study localized 2-D patterns of active ROP patches, referred to as spots, in a two-component RD model for the initiation of RH in the plant cell Arabidopsis thaliana. This 2-D study complements our previous analyses in [4, 6] of corresponding 1-D localized patterns. In our RD model, where the 2-D domain is a projection of a 3-D cell membrane and cell wall, as shown in Fig. 2, the amount of available plant hormone auxin has been shown to play a key role in both the formation and number of 2-D localized ROP patches observed, while the spatial gradient of the auxin distribution is shown to lead to the alignment of these patches along the transversal midline of the cell. The interaction between localized active ROP patches and our specific 2-D domain geometry is mediated by the inactive ROP concentration, which has a long range diffusion. This long-range spatial interaction can trigger the formation of localized patches of ROP, which are then bound to the cell membrane. In addition, as is illustrated in Fig. 7, our results in Proposition 3.1 suggest that once ROP patches consisting of localized multi-spots are quickly formed, their location is controlled by the auxin gradient and RH cell shape, where the latter seems to play a more relevant part than in a less realistic 1-D setting (cf. [6]). In this way, a sensitive interplay between overall auxin levels and geometrical properties is crucial to regulate separation of active-ROP patches, which promotes the onset of RHs in mutants. This may qualitatively describe findings on the rhd6 mutation of Arabidopsis thaliana reported in [29].

From a mathematical viewpoint, our hybrid asymptotic-numerical analysis of the existence and slow dynamics of localized active ROP patches, in the absence of any ${\mathcal{O}}(1)$ time-scale instabilities, extends the previous analyses of [9, 24, 25, 35] of spot patterns for prototypical RD systems by allowing for a spatially-dependent coefficient in the nonlinear kinetics, representing the spatially heterogeneous auxin distribution. In particular, we have derived an explicit DAE system for the slow dynamics of active ROP patches that includes typical Green’s interaction terms that appear in other models (cf. [9, 24, 25, 35]), but that now includes a new term in the dynamics associated with the spatial gradient of the auxin distribution. This new term contributing to spot dynamics is shown to lead to a spot-pinning effect (cf. [32]) whereby localized active ROP patches become aligned along the transverse midline of the RH cell. We have also determined a specific criterion that predicts when an active ROP patch will become linearly unstable to a shape deformation. This criterion, formulated in terms of the source parameter $S_{j}$ for an individual spot encodes the domain geometry and parameters associated with the inactive ROP density. In addition, the effect of auxin on instability and shape of active-ROP patches is captured by this parameter, and the nonlinear algebraic system for the collection of source parameters encodes the spatial interaction between patches via a Green’s matrix. When $S_{j}$ exceeds a threshold, the predicted linear instability of peanut-splitting type is shown to lead to a nonlinear spot self-replication event.

Although our asymptotic analysis of localized active ROP patches has been developed for an arbitrary spatial distribution of auxin, we have focused our study to two specific, biologically reasonable, forms for the auxin distribution given in Eq. 1.1. The first form assumes that the auxin concentration is monotone decreasing in the longitudinal direction, as considered in [6] and [4] and motivated experimentally, with no spatially transverse component. The second form is again monotone decreasing in the longitudinal direction, but allows for a spatially transverse component. This second form is motivated by the fact that our 2-D model is a projection of a 3-D cell membrane and cell wall. As a result, it is biologically plausible that some auxin can diffuse out of the transverse boundaries of our 2-D projected domain, leading to lower auxin levels near the transverse boundaries than along the midline of the cell.

We performed a numerical bifurcation analysis on the full RD model for both specific forms for auxin to study how the level of auxin determines both the number of ROP patches that will appear and the increasing complexity of the spatial pattern. This study showed similar features as in the 1-D studies in [4] in that there are stable solution branches that overlap (see Fig. 5 and Fig. 11), which allows for hysteretic transitions between various 2-D pattern types through fold-point bifurcations. In addition, owing to a homoclinic snaking type structure, a creation-annihilation cascade yields 2-D spatial patterns where both spots, baby droplets, and peanut structures can occur. Starting from unstable steady-states, our full time-dependent simulations have shown that $\mathcal{O}(1)$ time-scale competition and self-replication instabilities can occur leading to transitions between various spatial patterns. For the auxin distribution model with a transverse spatial dependence, we showed numerically that the active ROP patches are more confined to the cell midline and that there are no active ROP homoclinic stripes.

As an extension to our analysis, it would be worthwhile to consider a more realistic domain geometry rather than the projected 2-D geometry considered herein. In particular, it would be interesting to study the effect of transversal curvature on the dynamics and stability of localized active ROP patches. In addition, in [4, 6] and here, transport dynamics of ROPs have been modeled under the assumption of a standard Brownian diffusion process. Even though the experimental observations of pinning in [13, 15] seem to reproduced, at least qualitatively, by our model, it would be interesting to extend our model to account for, possibly, more realistic diffusive processes. In particular, it would be interesting to consider a hyperbolic-type diffusive process having a finite propagation speed for signals, as governed by the Maxwell–Cattaneo transfer law (cf. [28, 39]), or to allow for anomalous diffusion (cf. [19]).

Acknowledgments

D. A., V. F. B.–M. and M. J. W. were supported by EPSRC grant EP/P510993/1 (United Kingdom), UNAM–PAPIIT grant IA104316–RA104316 (Mexico) and NSERC Discovery Grant 81541 (Canada), respectively.

Bibliography39

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] E. S. A. Aftalion and S. Serfaty , Pinning Phenomena In The Ginzburg-Landau Model Of Superconductivity , J. Math. Pures Appl., 80 (2001), pp. 339–372.
2[2] D. Avitabile , Numerical Computation Of Coherent Structures In Spatially-Extended Systems . Second International Conference on Mathematical Neuroscience, Antibes Juan-les-Pins, 2016, https://icmns 2016.inria.fr/pre-conference-workshop/ .
3[3] D. Avitabile and H. Schmidt , Snakes And Ladders In An Inhomogeneous Neural Field Model , Physica D: Nonlinear Phenomena, 294 (2015), pp. 24–36.
4[4] V. Brena–Medina, D. Avitabile, A. Champneys, and M. Ward , Stripe To Spot Transition In A Plant Root Hair Initiation Model , SIAM J. Appl. Math., 75 (2015), pp. 1090–1119.
5[5] V. Brena–Medina and A. R. Champneys , Subcritical Turing Bifurcation And The Morphogenesis Of Localized Structures , Phys. Rev. E, 90 (2014), p. 032923.
6[6] V. Brena–Medina, A. R. Champneys, C. Grierson, and M. J. Ward , Mathematical Modelling Of Plant Root Hair Initiation; Dynamics Of Localized Patches , SIAM J. Appl. Dyn. Sys., 13 (2014), pp. 210–248.
7[7] J. Burke and E. Knobloch , Localized States In The Generalized Swift-Hohenberg Equation , Phys. Rev. E, 73 (2006).
8[8] W. Chen and M. Ward , The Stability And Dynamics Of Localized Spot Patterns In The Two-Dimensional Gray–Scott Model , SIAM J. Appl. Dyn. Sys., 10 (2011), pp. 582–666.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Spot dynamics in a reaction-diffusion model of plant root hair

Abstract

1 Introduction

2 Asymptotic Regime for Spot Formation

2.1 A Multiple-Spot Quasi Steady-State Pattern

Proposition 2.1**.**

2.2 Spots and Baby Droplets

3 The Slow Dynamics of Spots

3.1 Solvability Condition

Proposition 3.1**.**

3.2 Comparison Between Theory and Asymptotics for Slow Spot Dynamics

4 Fast O(1)\mathcal{O}(1)O(1) Time-Scale Instabilities

4.1 The Self-Replication Threshold

4.2 Numerical Illustrations of O(1){\mathcal{O}}(1)O(1) Time-Scale

5 More Robust 2D Patches and Auxin Transport

5.1 A 2-D Auxin Gradient

5.2 Instabilities with a 2-D Spatially-Dependent Auxin Gradient

6 Discussion

Acknowledgments

Proposition 2.1.

Proposition 3.1.

4 Fast $\mathcal{O}(1)$ Time-Scale Instabilities

4.2 Numerical Illustrations of ${\mathcal{O}}(1)$ Time-Scale