Symmetry-breaking bifurcations for compartmental reaction kinetics coupled by two bulk diffusing species with comparable diffusivities in 2-D

For a 2-D coupled PDE-ODE bulk-cell model, we investigate symmetry-breaking bifurcations that can emerge when two bulk diffusing species are coupled to two-component nonlinear intracellular reactions that are restricted to occur only within a disjoint collection of small circular compartments, or “cells,” of a common small radius that are confined in a bounded 2-D domain. Outside of the union of these cells, the two bulk species with comparable diffusivities and bulk degradation rates diffuse and globally couple the spatially segregated intracellular reactions through Robin boundary conditions across the cell boundaries, which depend on certain membrane reaction rates. In the singular limit of a small common cell radius, we construct steady-state solutions for the bulk-cell model and formulate a nonlinear matrix eigenvalue problem that determines the linear stability properties of the steady-states. For a certain spatial arrangement of cells for which the steady-state and linear stability analysis become highly tractable, we construct a symmetric steady-state solution where the steady-states of the intracellular species are the same for each cell. As regulated by the ratio of the membrane reaction rates on the cell boundaries, we show for various specific prototypical intracellular reactions, and for a specific two-cell arrangement, that our 2-D coupled PDE-ODE model admits symmetry-breaking bifurcations from this symmetric steady-state, leading to linearly stable asymmetric patterns, even when the bulk diffusing species have comparable or possibly equal diffusivities. Overall, our analysis shows that symmetry-breaking bifurcations can occur without the large diffusivity ratio requirement for the bulk diffusing species as is well-known from a Turing stability analysis applied to a spatially uniform steady-state for typical two-component activator-inhibitor systems. Instead, for our theoretical compartmental-reaction diffusion bulk-cell model, our analysis shows that the emergence of stable asymmetric steady-states can be controlled by the ratio of the membrane reaction rates for the two species. Bifurcation theoretic results for symmetric and asymmetric steady-state patterns obtained from our asymptotic theory are confirmed with full numerical PDE simulations.

1. Introduction

A central issue in many chemical and biological systems that involve the coupling of diffusive processes and nonlinear reactions is to determine conditions for which spatio-temporal patterns can form from either a patternless or a pre-patterned state. In a pioneering theoretical study, Turing [1] established that diffusing morphogens with different diffusivities can destablilize a spatially uniform and stable steady-state of the nonlinear reaction kinetics. As applied to two-component activator-inhibitor reaction-diffusion (RD) systems, this Turing stability analysis shows that a sufficiently large diffusivity ratio is typically needed to obtain spatial pattern formation from the destabilization of a spatially uniform state, unless the nonlinear reaction kinetics are finely tuned (cf. Pearson and Horsthemke [2], Baker et al. [3], and Diambra et al. [4]). For certain chemical systems, this large diffusivity ratio requirement needed for pattern formation may be feasible to achieve in situations where one of the chemical species can bind to a substrate, which has the consequence of reducing the effective diffusivity of this species (cf. Lengyel and Epstein [5] and Dulos et al. [6]). However, in many cellular processes related to developmental biology and morphogenesis, the theoretical large diffusivity ratio threshold needed for freely diffusing morphogens to create symmetry-breaking patterns is often unrealistic as different small molecules typically have very comparable diffusivities (cf. Müller et al. [7] and Rauch and Millonas [8]). In Müller et al. [7], various modifications of the simple “freely diffusing” morphogen paradigm such as facilitated diffusion, transient binding, immobilization and transcytosis, among others, have been postulated to play a central role in specific applications of diffusive transport at the cellular level. Qualitatively, the postulated overall effect of these mechanisms is to modify an effective diffusivity ratio of the morphogens, which can, therefore, lead to the emergence of spatial patterns and symmetry-breaking behavior in cellular processes related to developmental biology and early morphogenesis (cf. Sozen et al. [9]).

As a result, one key long-standing theoretical question in RD theory is how to modify the two-component RD paradigm so as to robustly generate stable spatial patterns from a spatially homogeneous state when the time scales for diffusion of the interacting species are comparable. By including an additional non-diffusible component, which roughly models either membrane-bound proteins or an immobile chemically active substrate, it has been shown (cf. Pearson [10], Klika et al. [11], and Korvasová et al. [12]) that this “2+1” extension of the two-component RD framework can yield stable spatial patterns even when the two diffusible species have a common diffusivity. In another direction, which is based on graph-theoretic properties associated with nonlinear reactions between multiple species that are either immobile or freely diffusing, it has been shown that with certain activating and inhibiting feedback relations in the chemical kinetics, spatial patterns can form without the large diffusivity ratio requirement (cf. Marcon et al. [13], Diego et al. [14], and Landge et al. [15]). More recently, the authors in Haas and Goldstein [16] have revealed that in random, multi-component, RD systems the required diffusivity threshold for pattern formation typically decreases as the number of interacting and diffusing species increases.

From a theoretical viewpoint, in specific applications where a large diffusivity ratio is a realistic assumption, it has been shown both analytically and from numerical simulations (cf. Vanag and Epstein [17], Ward [18], Halatek et al. [19], and Halatek and Frey [20]) that two-component RD systems admit a wide range of spatially localized patterns and instabilities that occur in the “far-from-equilibrium” regime, far from where a Turing linear stability analysis will provide any insight into pattern-forming properties.

The goal of this paper is to formulate and quantitatively analyze a new theoretical model in a 2-D setting that robustly leads to pattern formation even when the two diffusing species have a comparable or equal diffusivity. More specifically, we analyze symmetry-breaking pattern formation for a 2-D PDE-ODE bulk-cell RD model in which spatially segregated localized reaction compartments, referred to as “cells,” are coupled to a two-component linear bulk diffusion field with constant bulk degradation rates. In the cells, which are assumed to have a common radius that is small compared to the domain length-scale and the inter-cell distances, two-component intracellular activator-inhibitor reaction kinetics are specified. The intracelluar species undergo an exchange with the two bulk species across the cell boundaries, as mediated by membrane reaction rates in a Robin boundary condition that is specified on each cell boundary. The two extracellular diffusing bulk species, with comparable diffusivities and degradation rates, provide the mechanism that couples the nonlinear intracellular reactions that occur in the union of the spatially segregated cells. We refer to this modeling framework as a compartmental-reaction diffusion system.

The numerical implementation of our theoretical analysis for this model for various specific intracellular reaction kinetics reveals that it is the ratio of the reaction rate of the inhibitor component to that of the activator component on the compartment boundaries that plays a central role in the initiation of symmetry-breaking bifurcations of a symmetric steady-state. The magnitude of this ratio ultimately controls whether linearly stable asymmetric steady-states for the bulk-cell model can occur even when the bulk diffusivities are comparable or equal. The bifurcation threshold condition for this key membrane reaction rate ratio parameter is distinct from the usual large diffusivity ratio threshold that is required for pattern formation from a spatially uniform state for typical two-component activator-inhibitor RD systems (cf. Maini et al. [21] and Krause et al. [22]). We emphasize that our linear stability analysis predicting symmetry-breaking bifurcations for the bulk-cell model, as regulated by the membrane reaction rate ratio, is significantly more challenging than performing a simple Turing stability analysis [1] since it is based on the linearization of the bulk-cell model around a spatially non-uniform symmetric steady-state. In our previous 1-D study [23], where nonlinear reactions were restricted either to domain boundaries or at lattice site on a 1-D periodic chain, it has been shown for some specific nonlinear kinetics that symmetry-breaking bifurcations can occur from a symmetric steady-state when the ratio of membrane reaction rates exceeds a threshold.

We remark that our 2-D study, and related 1-D analysis in Pelz and Ward [23], is largely inspired by the agent-based numerical computations in Rauch and Millonas [8] where it was shown that nonlinear kinetic reactions restricted to lattice sites on a 2-D lattice can generate stable Turing-type spatial patterns when coupled through a spatially discretized two-component bulk diffusion field in which the two diffusible species have a comparable diffusivity.

In a broader context, the study of novel pattern-forming properties associated with compartmentalized reactions interacting through a passive bulk diffusion field originates from the 1-D analysis in Gomez-Marin et al. [24] for the FitzHugh-Nagumo model and the bulk-membrane analysis of Levine and Rappel [25] in disk-shaped domains. In a 1-D context, and with one bulk diffusing species, this compartmental-reaction diffusion system modeling paradigm has been shown to lead to triggered oscillatory instabilities for various reaction kinetics involving conditional oscillators (cf. Gou et al. [26] Gou and Ward [27], and Gou et al. [28]). Amplitude equations characterizing the local branching behavior for these triggered oscillations have been derived in Paquin-Lefebvre et al. [29] using a weakly nonlinear analysis. Applications of this framework have been used to model intracellular polarization and oscillations in fission yeast (cf. Xu and Bressloff [30] and Xu and Jilkine [31]). In a 2-D domain, similar bulk-cell models, but with only one diffusing bulk species, have been formulated and used to model quorum-sensing behavior (cf. Gou and Ward [32], Iyaniwura and Ward [33], Ridgway et al. [34], and Gomez et al. [35]). With regards to bulk-membrane RD models in a multi-spatial dimensional context, where nonlinear kinetics are restricted to the membrane, the associated pattern-forming properties have been studied both theoretically (cf. Rätz [36], Elliott et al. [37], Madzvamuse et al. [38], Madzvamuse and Chung [39], and Paquin-Lefebvre et al. [40]), and for some specific biological applications (cf. Cusseddu et al. [41], Rätz and Röger [42, 43], Stolerman et al. [44], and Paquin-Lefebvre et al. [45]).

The outline of this paper is as follows. In Section 2 we formulate our bulk-cell model and use a singular perturbation approach in the limit of a small common cell radius to derive a nonlinear algebraic system characterizing all steady-state solutions of the model. In Section 3 we show that the discrete eigenvalues of the linearization of the bulk-cell model around a steady-state solution are determined by a root-finding condition on a nonlinear matrix eigenvalue problem. For a certain type of spatial configuration of the cells, the bulk-cell model is shown to admit a symmetric steady-state solution in which the steady-states of the intracellular reactions are identical. The possibility of symmetry-breaking bifurcations along this symmetric steady-state solution branch, leading to the existence of linearly stable asymmetric patterns, are analyzed by applying solution path continuation software to our bifurcation-theoretic analytical results. For a certain two-cell configuration in the unit disk, and for either Rauch and Millonas [8], Gierer and Meinhardt [46], or FitzHugh and Nagumo [24] intracellular reactions, we show in Section 4 that it is the magnitude of the ratio of the reaction rates for the two bulk species on the cell membranes that controls whether linearly stable asymmetric patterns can bifurcate from the symmetric steady-state. Our theoretical predictions of symmetry-breaking behavior, leading to stable asymmetric steady-states even when the two bulk species have comparable or equal diffusivities, are confirmed from full PDE numerical simulations. For a closely-spaced arrangement of cells as is typical in biological tissues, and where our asymptotic theory no longer applies, the PDE numerical simulations shown in Section 4.4 illustrate that symmetry-breaking bifurcations can still be controlled by the reaction rate ratio on the cell boundaries. In particular, our numerical results suggest that such bifurcations occur with a smaller membrane reaction-rate ratio than for the situation where the cells are more spatially segregated. In Section 5 we discuss our theoretical results in a wider context, and suggest a few open directions.

2. Compartmental-reaction diffusion system in 2-D

2.1. Model formulation

We consider a bounded 2-D domain with length scale L, denoted by Ω^L ⊂ ℝ², that contains m disconnected circular compartments ${\Omega }_j^L$, for j ∈ {1, …, m}, referred to as “cells”. We will assume that these cells have a common radius that is small in comparison with the length scale L of the domain. The bulk or extracellular medium is the region ${\Omega }^L\backslash {\bigcup }_{j=1}^m{\Omega }_j^L$.

In the bulk we assume that there are two extracellularly diffusing and degrading chemical species with concentrations U and V. These messenger molecules are synthesized on the “cell” membranes through the interaction with two corresponding intracellular species M_j and H_j. With the molecule counts 𝔘, 𝔙, 𝔐_j and $ℌ$_j corresponding to respectively U, V, M_j and H_j, the chemical equations are

$\begin{array}{ll}𝔘\underset{[{\beta }_U]}{\overset{[{\beta }_U]}{\rightleftharpoons }}{𝔐}_j,𝔙\underset{[{\beta }_V]}{\overset{[{\beta }_V]}{\rightleftharpoons }}{ℌ}_j.& (2.1)\end{array}$

Here we made the assumption that the exponential forward reaction rates equal the backward reaction rates and that all compartments are identical in that they have common membrane reaction rates. The intra-compartmental species, in turn, are produced by certain reaction kinetics, denoted by f(M, H) and g(M, H), that are assumed to be identical in each compartment.

More precisely, in dimensional variables, our bulk-cell coupled model is

$\begin{array}{l}\text{bulk}\\ \{\begin{array}{ll}{\partial }_{T}U=D_U{\Delta }_{X}U-{\kappa }_{U}U,& \text{X}\in {\Omega }^L\backslash {\displaystyle {\bigcup }_{j=1}^m}{\Omega }_j^L,\\ {\partial }_{T}V=D_V{\Delta }_{X}V-{\kappa }_{V}V,& \text{X}\in {\Omega }^L\backslash {\displaystyle {\bigcup }_{j=1}^m}{\Omega }_j^L,\\ {\partial }_{{\tilde{n}}_X}U={\partial }_{{\tilde{n}}_X}V=0,& \text{X}\in \partial {\Omega }^L,\text{(Neumann condition)}\end{array}& (\mathrm{2.2a})\end{array}$

$\begin{array}{l}\text{reaction fluxes}\\ \{\begin{array}{ll}{D}_U{\partial }_{n_{j,X}}U={\beta }_{U,1}U-{\beta }_{U,2}{M}_j,& \text{X}\in \partial {\Omega }_j^L,\text{(Robin condition)}\\ D_V{\partial }_{n_{j,X}}V={\beta }_{V,1}V-{\beta }_{V,2}{H}_j,& \text{X}\in \partial {\Omega }_j^L,\end{array}& (\mathrm{2.2b})\end{array}$

$\begin{array}{l}\text{compartments}\\ \{\begin{array}{l}\frac{d}{dT}{M}_j={\kappa }_R{\mu }_{c}f\left(\frac{1}{{\mu }_c}{M}_j,\frac{1}{{\mu }_c}{H}_j\right)+{\int }_{\partial {\Omega }_j^L}({\beta }_{U,1}U-{\beta }_{U,2}{M}_j)d{S}_X,\\ \frac{d}{dT}{H}_j={\kappa }_R{\mu }_{c}g\left(\frac{1}{{\mu }_c}{M}_j,\frac{1}{{\mu }_c}{H}_j\right)+{\int }_{\partial {\Omega }_j^L}({\beta }_{V,1}V-{\beta }_{V,2}{H}_j)d{S}_X\\ \text{(reaction kinetics)},\end{array}& (\mathrm{2.2c})\end{array}$

with j ∈ {1, …, m} and where n_j,X is the outward unit normal vector to ${\Omega }_j^L$ while $\tilde{n}$_X is the outward unit normal vector to Ω^L. The diffusivities (diffusion coefficients) for U and V are D_U and D_V, and U and V are degrading in the bulk with exponential rate constants κ_U and κ_V, respectively. The exponential reaction rates on the compartment boundaries are β_U and β_V with corresponding rates β_U,1 and β_V,1 per area times length and β_U,2 and β_V,2 per length and time units, and μ_c is a normalizing constant for the intracellular species. Lastly, κ_R is a dimensional reaction rate for the intracellular reactions.

In Appendix A we non-dimensionalize (Equation 2.2) to obtain the dimensionless PDE-ODE model

$\begin{array}{l}\text{bulk}\\ \{\begin{array}{ll}{\partial }_{t}u=D_u\Delta u-{\sigma }_{u}u,& \text{x}\in \Omega \backslash {\displaystyle {\bigcup }_{j=1}^m}{\Omega }_j,\\ {\partial }_{t}v=D_v\Delta v-{\sigma }_{v}v,& \text{x}\in \Omega \backslash {\displaystyle {\bigcup }_{j=1}^m}{\Omega }_j,\\ {\partial }_{\tilde{n}}u={\partial }_{\tilde{n}}v=0,& \text{x}\in \partial \Omega ,\end{array}& (\mathrm{2.3a})\end{array}$

$\begin{array}{l}\text{reaction fluxes}\\ \{\begin{array}{ll}\epsilon D_u{\partial }_{n_j}u=d_1^{u}u-d_2^u{\mu }_j,& \text{x}\in \partial {\Omega }_j,\\ \epsilon D_v{\partial }_{n_j}v=d_1^{v}v-d_2^v{\eta }_j,& \text{x}\in \partial {\Omega }_j,\end{array}& (\mathrm{2.3b})\end{array}$

$\begin{array}{l}\text{compartments}\\ \{\begin{array}{l}\frac{d{\mu }_j}{dt}=f({\mu }_j,{\eta }_j)+\frac{1}{\epsilon }{\int }_{\partial {\Omega }_j}(d_1^{u}u-d_2^u{\mu }_j)dS,\\ \frac{d{\eta }_j}{dt}=g({\mu }_j,{\eta }_j)+\frac{1}{\epsilon }{\int }_{\partial {\Omega }_j}(d_1^{v}v-d_2^v{\eta }_j)dS,\end{array}& (\mathrm{2.3c})\end{array}$

for j ∈ {1, …, m}. Here n and $\tilde{n}$ are the outward unit normal vectors to Ω_j and Ω, respectively, and we have dropped the label “x" for Δ and the outward unit normal vectors. In Equation (2.3), the compartments are disks of a common radius ε ≪ 1 centered at x_j ∈ Ω, i.e. Ω_j ≡ {x||x − x_j| ≤ ε}. We will refer to $d_1^u$, $d_1^v$, $d_2^u$, and $d_2^v$ as dimensionless membrane reaction rates. An illustration of the bulk-cell model is shown in Figure 1.

FIGURE 1

Figure 1. A 2-D bounded domain with four diffusion-coupled circular cells of a common radius. In the j^th cell, activator-inhibitor reaction kinetics occur for the activator μ_j and inhibitor η_j. Across the cell membrane, there is an exchange between the intracellular and bulk species. In the bulk, u and v diffuse and undergo degradation. The cells are assumed to be small but are drawn larger here for illustration only.

We will use strong localized perturbation theory [18] to construct the steady-state solutions of Equation (2.3) and to analyze their linear stability properties in the asymptotic limit ε ≪ 1 and under the assumption that m circular cells are well-separated in the sense that the cell centers satisfy $|{\text{x}}_i-{\text{x}}_j|={O}(1)$, for i, j ∈ {1, …, m} and i ≠ j.

2.2. Asymptotic construction of the steady-states

Our main goal is to construct a symmetric steady-state solution for Equation (2.3) in which the concentration of each species is the same inside and in the local vicinity of each compartment. We will show below that even when the bulk diffusing species have comparable diffusivities this symmetric steady-state is unstable to symmetry-breaking perturbations that occur beyond a pitchfork bifurcation point associated with the membrane reaction rate ratio $\rho \equiv d_1^v/d_1^u=d_2^v/d_2^u$. This leads to the existence of linearly stable asymmetric steady-state solutions to Equation (2.3).

In the absence of diffusion, the ODE system for the intra-compartmental species is decoupled from the bulk medium and reduces to

$\begin{array}{ll}\stackrel{.}{\mu }(t)=f(\mu ,\eta ),\stackrel{.}{\eta }(t)=g(\mu ,\eta ).& (2.4)\end{array}$

Let (μ_e, η_e) be an equilibrium point for Equation (2.4) and label F(μ, η) ≡ (f(μ, η), g(μ, η)). For a specific parameter set, the linear stability property of the equilibrium state is characterized by whether the eigenvalues λ of the Jacobian matrix DF(μ_e, η_e) have positive (unstable, exponentially growing perturbations) or negative (stable, exponentially decaying perturbations) real parts Re(λ). However, when there is bulk diffusion and the compartments are coupled through the bulk, the steady-state solution in the compartments depends on the bulk diffusivities, the membrane reaction rates, and the spatial configuration of the cells.

We now use the method of matched asymptotic expansions to construct steady-state solutions for Equation (2.3). In the j^th inner region, defined within an ${O}(\epsilon )$ neighborhood of the boundary of the j^th cell, we introduce the local variables ${\text{y}}_j={\epsilon }^{-1}(\text{x}-{\text{x}}_j)$, u_j(x) = u(εy_j + x_j), and v_j(x) = v(εy_j + x_j), where p_j ≡ |y_j|. Upon writing the steady-state of Equation (2.3a) in terms of the inner variables, for ε → 0 the steady-state problem in the j^th inner region is Δu_j = 0 and Δv_j = 0, for p_j ≥ 1, subject to $D_u{\partial }_{p_j}{u}_j=d_1^u{u}_j-d_2^u{\mu }_j$ and $D_v{\partial }_{p_j}{v}_j=d_1^v{v}_j-d_2^v{\eta }_j$ on p_j = 1. The radially symmetric solutions to these problems are

$\begin{array}{l}{u}_j(p_j)=A_j^{u}log{p}_j+\frac{1}{d_1^u}\left(D_u{A}_j^u+d_2^u{\mu }_j\right),\\ v_j(p_j)=A_j^{v}log{p}_j+\frac{1}{d_1^v}\left(D_v{A}_j^v+d_2^v{\eta }_j\right),& (2.5)\end{array}$

for j ∈ {1, …, m}, where $A_j^u$ and $A_j^v$ for j = 1, …, m are constants to be determined. Upon substituting (Equation 2.5) into the steady-state problem of Equation (2.3c), we obtain for the j^th cell that

$\begin{array}{l}f({\mu }_j,{\eta }_j)+2\pi D_u{A}_j^u=0,g({\mu }_j,{\eta }_j)+2\pi D_v{A}_j^v=0,\\ j\in \left\{1,\dots ,m\right\}.& (2.6)\end{array}$

Next, we must determine $A_j^u$ and $A_j^v$ by matching the far-field behavior of the inner solutions (Equation 2.5) to the outer solutions defined in the bulk region.

In the limit ε → 0, in the bulk region the compartments formally shrink to points and from the far-field behavior of Equation (2.5), when written in outer variables, we obtain that the steady-state bulk species U satisfies

$\begin{array}{ll}\Delta U-{\omega }_u^{2}U=0,\text{x}\in \Omega \backslash \left\{{\text{x}}_1,\dots ,{\text{x}}_m\right\};{\partial }_{n}U=0,\text{x}\in \partial \Omega ;& (\mathrm{2.7a})\end{array}$

$\begin{array}{l}U~A_j^{u}log|\text{x}-{\text{x}}_j|+\frac{A_j^u}{\nu }+\frac{1}{d_1^u}(D_u{A}_j^u+d_2^u{\mu }_j),\\ \text{as\hspace{1em}}\text{x}\to {\text{x}}_j,j\in \left\{1,\dots ,m\right\},& (\mathrm{2.7b})\end{array}$

where ν ≡ −1/log ε ≪ 1 and ${\omega }_u\equiv \sqrt{{\sigma }_u/D_u}$. Similarly, with ${\omega }_v\equiv \sqrt{{\sigma }_v/D_v}$, for the bulk species V we have that

$\begin{array}{l}\Delta V-{\omega }_v^{2}V=0,\text{x}\in \Omega \backslash \left\{{\text{x}}_1,\dots ,{\text{x}}_m\right\}\\ {\partial }_{n}V=0,\text{x}\in \partial \Omega ;& (\mathrm{2.8a})\end{array}$

$\begin{array}{l}V~A_j^{v}log|\text{x}-{\text{x}}_j|+\frac{A_j^v}{\nu }+\frac{1}{d_1^v}(D_v{A}_j^v+d_2^v{\eta }_j)\\ \text{as\hspace{1em}}\text{x}\to {\text{x}}_j,j\in \left\{1,\dots ,m\right\}.& (\mathrm{2.8b})\end{array}$

To represent solutions to (2.7) and (2.8), we introduce the reduced-wave Green's function G_ω(x, x_j) that satisfies

$\begin{array}{ll}\Delta G_{\omega }-{\omega }^2{G}_{\omega }=-\delta (\text{x}-{\text{x}}_j),\text{x}\in \Omega ;{\partial }_n{G}_{\omega }=0,\text{x}\in \partial \Omega ;& (\mathrm{2.9a})\end{array}$

$\begin{array}{ll}{G}_{\omega }~-\frac{1}{2\pi }log|\text{x}-{\text{x}}_j|+R_{\omega }({\text{x}}_j)+o(1),\text{as\hspace{1em}}\text{x}\to {\text{x}}_j.& (\mathrm{2.9b})\end{array}$

Here R_ω(x_j) is the regular, or non-singular, part of the singularity at x = x_j. The solutions to (2.7) and (2.8) are represented as

$\begin{array}{ll}U(\text{x})=-2\pi {\displaystyle \sum _{i=1}^m}{A}_i^u{G}_{{\omega }_u}(\text{x};{\text{x}}_i),V(\text{x})=-2\pi {\displaystyle \sum _{i=1}^m}{A}_i^v{G}_{{\omega }_v}(\text{x};{\text{x}}_i).& (2.10)\end{array}$

The pre-specification of the regular part of each singularity condition in Equations (2.7), (2.8) yields a constraint. These constraints provide algebraic systems for $A_j^u$ and $A_j^v$ for j ∈ {1, …, m}. By expanding (Equation 2.10) as x → x_j, we enforce that that the non-singular terms in the resulting expression agree with the conditions that are required in Equations (2.7b), (2.8b) for each j ∈ {1, …, m}. This leads to linear algebraic systems for ${{A}}^u\equiv {(A_1^u,\dots ,A_m^u)}^T$ and ${{A}}^v\equiv {(A_1^v,\dots ,A_m^v)}^T$, given in matrix form by

$\begin{array}{l}\left(\left(1+\frac{\nu D_u}{d_1^u}\right)I+2\pi \nu {{G}}_{{\omega }_u}\right){{A}}^u=-\frac{\nu d_2^u}{d_1^u}\mu ,\\ \left(\left(1+\frac{\nu D_v}{d_1^v}\right)I+2\pi \nu {{G}}_{{\omega }_v}\right){{A}}^v=-\frac{\nu d_2^v}{d_1^v}\eta ,& (2.11)\end{array}$

where $\mu \equiv {({\mu }_1,\dots ,{\mu }_m)}^T$ and $\eta \equiv {({\eta }_1,\dots ,{\eta }_m)}^T$. In Equation (2.11), ${{G}}_{\omega }$ with either ω = ω_u or ω = ω_v is the symmetric reduced-wave Greens' interaction matrix defined by

$\begin{array}{ll}{{G}}_{\omega }\equiv \left(\begin{array}{cccc}{R}_{\omega 1}& G_{\omega 12}& \dots & G_{\omega 1m}\\ G_{\omega 21}& R_{\omega 2}& \dots & G_{\omega 2m}\\ ⋮& ⋮& \ddots & ⋮\\ G_{\omega m1}& G_{\omega m2}& \dots & R_{\omega m}\end{array}\right).& (2.12)\end{array}$

Here G_ωji = G_ωij ≡ G_ω(x_j; x_i) for i ≠ j, and R_ωj ≡ R_ω(x_j) for j ∈ {1, …, m}, are obtained from the solution to Equation (2.9).

To determine a nonlinear algebraic system that characterizes our steady-state solution, we solve (Equation 2.11) for ${{A}}^v$ and ${{A}}^u$, and substitute the resulting expressions into Equation (2.6). In this way, we obtain a 2m dimensional nonlinear algebraic system for μ_j and η_j, for j = 1, …, m, given by

$\begin{array}{l}f({\mu }_j,{\eta }_j)-{\text{e}}_j^T{\Theta }_u\mu =0,g({\mu }_j,{\eta }_j)-{\text{e}}_j^T{\Theta }_v\eta =0,\\ \text{ for\hspace{1em}}j\in \left\{1,\dots ,m\right\},& (\mathrm{2.13a})\end{array}$

Where ${\text{e}}_j\equiv {(0,\dots ,0,1,0,\dots ,0)}^T$ is the unit vector in the j^th direction. In Equation (2.13a), Θ_u and Θ_v are defined by

$\begin{array}{l}{\Theta }_u\equiv 2\pi \nu D_u\frac{d_2^u}{d_1^u}{\left[\left(1+\frac{\nu D_u}{d_1^u}\right)I+2\pi \nu {{G}}_{{\omega }_u}\right]}^{-1},\\ {\Theta }_v\equiv 2\pi \nu D_v\frac{d_2^v}{d_1^v}{\left[\left(1+\frac{\nu D_v}{d_1^v}\right)I+2\pi \nu {{G}}_{{\omega }_v}\right]}^{-1}.& (\mathrm{2.13b})\end{array}$

We can simplify our steady-state analysis for the special case where g(μ, η) is linear and inhibiting in η, with the form

$\begin{array}{ll}g(\mu ,\eta )=g_1(\mu )-g_2\eta ,& (2.14)\end{array}$

where g₂ ≥ 0 is a constant. This specific form applies to Rauch and Millonas [8], FitzHugh and Nagumo [24], and Gierer and Meinhardt [46], reaction kinetics, and is relevant for the illustrations of the theory given in Section 4. In this case, we obtain from the second equation in Equation (2.13a) that

$\begin{array}{ll}\eta ={\left[g_{2}I+{\Theta }_v\right]}^{-1}{\text{g}}_1\text{where\hspace{1em}}{\text{g}}_1\equiv {(g_1({\mu }_1),\dots ,g_1({\mu }_m))}^T.& (2.15)\end{array}$

Then, from the first equation in Equation (2.13a) we obtain an m dimensional nonlinear algebraic system for $\mu ={({\mu }_1,\dots ,{\mu }_m)}^T$ given by

$\begin{array}{ll}f\left({\mu }_j,{\text{e}}_j^T{(g_{2}I+{\Theta }_v)}^{-1}{\text{g}}_1\right)-{\text{e}}_j^T{\Theta }_u\mu =0,j\in \left\{1,\dots ,m\right\}.& (2.16)\end{array}$

Next, we define a symmetric cell arrangement for which the steady-state analysis can be further simplified.

Definition 2.1. A symmetric cell arrangement is defined by the condition that the symmetric Green's matrix ${{G}}_{\omega }$ satisfies the following two properties:

• Property 1: e ≡ (1, …, 1)^T is an eigenvector of ${{G}}_{\omega }$ for all ω > 0:

• Property 2: The eigenspace of ${{G}}_{\omega }$ orthogonal to e is independent of ω.

These two properties certainly hold when ${{G}}_{\omega }$ is a circulant matrix. In particular, ${{G}}_{\omega }$ is a circulant matrix when m small cells are equidistantly spaced on a ring that is concentric within a circular domain Ω. Such an arrangement of cells is called a ring pattern.

For a symmetric cell arrangement, ${{G}}_{{\omega }_u}$ and ${{G}}_{{\omega }_v}$ have a common eigenspace, and so we can seek a symmetric solution to Equation (2.13) of the form

$\begin{array}{ll}\mu ={\mu }_c\text{e},\eta ={\eta }_c\text{e},{{A}}^u=A_c^u\text{e},{{A}}^v=A_c^v\text{e},& (2.17)\end{array}$

where the scalars μ_c, η_c, $A_c^u$, and $A_c^v$ are to be found. Upon substituting (Equation 2.17) into Equation (2.13), we obtain that μ_c and η_c satisfy the nonlinear algebraic system

$\begin{array}{ll}f({\mu }_c,{\eta }_c)-{\alpha }_u{\mu }_c=0,g({\mu }_c,{\eta }_c)-{\alpha }_v{\eta }_c=0,& (2.18)\end{array}$

where α_u and α_v, denoting the eigenvalues of Θ_u and Θ_v for the eigenvector e, respectively, are defined by

$\begin{array}{ll}{\alpha }_u\equiv \frac{2\pi \nu D_u{d}_2^u/d_1^u}{1+\nu D_u/d_1^u+2\pi \nu {\kappa }_u},{\alpha }_v\equiv \frac{2\pi \nu D_v{d}_2^v/d_1^v}{1+\nu D_v/d_1^v+2\pi \nu {\kappa }_v}.& (\mathrm{2.19a})\end{array}$

Here κ_u and κ_v are the eigenvalues of the Green's matrices for the eigenvector e, given by

$\begin{array}{ll}{{G}}_{{\omega }_u}\text{e}={\kappa }_u\text{e},{{G}}_{{\omega }_v}\text{e}={\kappa }_v\text{e}.& (\mathrm{2.19b})\end{array}$

Moreover, if g(μ, η) has the specific form in Equation (2.14), we obtain from Equation (2.16) that for a symmetric pattern of cells, there is a symmetric steady-state solution whenever there is a root μ_c to the scalar nonlinear algebraic equation

$\begin{array}{ll}f\left({\mu }_c,\frac{g_1({\mu }_c)}{g_2+{\alpha }_v}\right)-{\alpha }_u{\mu }_c=0.& (2.20)\end{array}$

In summary, for a symmetric pattern of cells, the asymptotic construction of a symmetric steady-state solution for Equation (2.3) is reduced to the much simpler problem of determining a solution to the two-dimensional nonlinear algebraic problem (Equation 2.18) for general reaction kinetics, or to Equation (2.20) when g has the specific form in Equation (2.14). In these algebraic problems, the eigenvalues κ_u and κ_v, as needed in Equation (2.19a), are the constant row sums of the Green's matrices for the two bulk species. The bulk diffusivities, the membrane reaction rates, and the spatial configuration of the cells all influence α_u and α_v.

2.3. Symmetry-breaking bifurcations

To detect any symmetry-breaking pitchfork bifurcation points along the symmetric steady-state solution branch we can perform a linear stability analysis of Equation (2.3) around the steady-state solution and seek λ = 0 eigenvalue crossings. An equivalent, but simpler, approach to detect zero-eigenvalue crossings for the linearized problem is to determine bifurcation points associated with the linearization of the nonlinear algebraic system (Equation 2.13) around a symmetric steady-state.

To do so, we introduce the perturbations

$\begin{array}{l}\mu ={\mu }_c\text{e}+\tilde{\mu },\eta ={\eta }_c\text{e}+\tilde{\eta },\\ {{A}}^u=A_c^u\text{e}+{\tilde{{A}}}^u,{{A}}^v=A_c^v\text{e}+{\tilde{{A}}}^v,& (2.21)\end{array}$

into Equation (2.13) and linearize the resulting system. In this way, we obtain that a symmetry-breaking bifurcation occurs if and only if there is a non-trivial solution to the 2m × 2m homogeneous linear system

$\begin{array}{ll}\left(\begin{array}{cc}{f}_{\mu }^{c}I-{\Theta }_u& f_{\eta }^{c}I\\ g_{\mu }^{c}I& g_{\eta }^{c}I-{\Theta }_v\end{array}\right)\left(\begin{array}{c}\tilde{\mu }\\ \tilde{\eta }\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),& (2.22)\end{array}$

at some point along the symmetric solution branch given by Equation (2.18). In Equation (2.22) we have labeled $f_{\mu }^c$ by $f_{\mu }^c\equiv {\partial }_{\mu }f(\mu ,\eta )$ when evaluated at μ = μ_c and η = η_c, while I is the m × m identity matrix. For the special case where g has the specific form in Equation (2.14), we can solve (Equation 2.22) for $\tilde{\eta }$ and reduce (Equation 2.22) to the m-dimensional homogeneous linear system

$\begin{array}{ll}\left(f_{\mu }^{c}I+f_{\eta }^c{g}_1^{\prime }({\mu }_c){(g_{2}I+{\Theta }_v)}^{-1}-{\Theta }_u\right)\tilde{\mu }=0.& (2.23)\end{array}$

Next, by Property 2 for a symmetric cell arrangement, it follows that ${{G}}_{{\omega }_u}$ and ${{G}}_{{\omega }_v}$ have a common orthogonal eigenspace ${{Q}}^{\perp }\equiv \mathrm{\text{span}}\left\{{\text{q}}_2,\dots ,{\text{q}}_m\right\}$, where ${\text{q}}_j^T\text{e}=0$ for j ∈ {2, …, m} and ${\text{q}}_i^T{\text{q}}_j=0$ for i ≠ j. The eigenvalues of ${{G}}_{{\omega }_u}$ and ${{G}}_{{\omega }_v}$ in this common eigenspace are labeled by

$\begin{array}{ll}{{G}}_{{\omega }_u}{\text{q}}_j={\kappa }_{u,j}^{\perp }{\text{q}}_j,{{G}}_{{\omega }_v}{\text{q}}_j={\kappa }_{v,j}^{\perp }{\text{q}}_j,j\in \left\{2,\dots ,m\right\},& (2.24)\end{array}$

so that

$\begin{array}{ll}{\Theta }_u{\text{q}}_j={\alpha }_{u,j}^{\perp }{\text{q}}_j,{\Theta }_v{\text{q}}_j={\alpha }_{v,j}^{\perp }{\text{q}}_j,j\in \left\{2,\dots ,m\right\},& (2.25)\end{array}$

with

$\begin{array}{ll}{\alpha }_{u,j}^{\perp }\equiv \frac{2\pi \nu D_u{d}_2^u/d_1^u}{1+\nu D_u/d_1^u+2\pi \nu {\kappa }_{u,j}^{\perp }},{\alpha }_{v,j}^{\perp }\equiv \frac{2\pi \nu D_v{d}_2^v/d_1^v}{1+\nu D_v/d_1^v+2\pi \nu {\kappa }_{v,j}^{\perp }}.& (2.26)\end{array}$

By setting $\tilde{\mu }={\tilde{\mu }}_c{\text{q}}_j$ and $\tilde{\eta }={\tilde{\eta }}_c{\text{q}}_j$ in Equation (2.22), we conclude that a symmetry-breaking bifurcation occurs for the j^th mode with j ∈ {2, …, m} whenever

$\begin{array}{ll}\left(\begin{array}{cc}{f}_{\mu }^c-{\alpha }_{u,j}^{\perp }& f_{\eta }^c\\ g_{\mu }^c& g_{\eta }^c-{\alpha }_{v,j}^{\perp }\end{array}\right)\left(\begin{array}{c}{\tilde{\mu }}_c\\ {\tilde{\eta }}_c\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),& (2.27)\end{array}$

has a nontrivial solution. This is equivalent to the condition that

$\begin{array}{ll}\left(f_{\mu }^c-{\alpha }_{u,j}^{\perp }\right)\left(g_{\eta }^c-{\alpha }_{v,j}^{\perp }\right)-f_{\eta }^c{g}_{\mu }^c=0,j\in \left\{2,\dots ,m\right\},& (2.28)\end{array}$

is satisfied at some point along the symmetric solution branch defined by the solution to (2.18).

Finally, for the special case where g has the specific form in Equation (2.14), we obtain that there is a symmetry-breaking bifurcation for the j^th mode, with j ∈ {2, …, m}, when there is a root to the scalar problem

$\begin{array}{ll}{f}_{\mu }^c+\frac{f_{\eta }^c{g}_1^{\prime }({\mu }_c)}{g_2+{\alpha }_{v,j}^{\perp }}-{\alpha }_{u,j}^{\perp }=0,& (2.29)\end{array}$

whenever η_c = g₁(μ_c)/(g₂ + α_v) where μ_c satisfies (Equation 2.20).

In the examples shown in Section 4 we will use $\rho \equiv d_1^v/d_1^u=d_2^v/d_2^u$ as the bifurcation parameter to detect whether symmetry-breaking bifurcations can occur along the symmetric solution branch.

2.4. A symmetric cell arrangement with two cells

Consider a symmetric cell arrangement with two cells, i.e. m = 2, for the special case where g has the form in Equation (2.14). Then, to determine all steady-state solutions we need only solve the nonlinear algebraic system (Equation 2.16) for μ₁ and μ₂. The symmetric steady-state solution, for which μ_c ≡ μ₁ = μ₂, is obtained by solving the scalar problem (Equation 2.20). To detect whether symmetry-breaking bifurcations can occur, we note that ${\text{q}}_2={(1,-1)}^T$ spans the common eigenspace of ${{G}}_{{\omega }_u}$ and ${{G}}_{{\omega }_v}$ orthogonal to e, and that κ_u = R_ωu1 − G_ωu12 and κ_v = R_ωv1 − G_ωv12 are the associated eigenvalues for q₂. This yields that the root-finding condition (Equation 2.29) becomes

$\begin{array}{ll}{f}_{\mu }^c+\frac{f_{\eta }^c{g}_1^{\prime }({\mu }_c)}{g_2+{\alpha }_{v,2}^{\perp }}-{\alpha }_{u,2}^{\perp }=0,& (2.30)\end{array}$

where in terms of the entries of the Green's matrices we have

$\begin{array}{l}{\alpha }_{u,2}^{\perp }\equiv \frac{2\pi \nu D_u{d}_2^u/d_1^u}{1+\nu D_u/d_1^u+2\pi \nu \left[R_{{\omega }_{u}1}-G_{{\omega }_{u}12}\right]},\\ {\alpha }_{v,2}^{\perp }\equiv \frac{2\pi \nu D_v{d}_2^v/d_1^v}{1+\nu D_v/d_1^v+2\pi \nu \left[R_{{\omega }_{v}1}-G_{{\omega }_{v}12}\right]}.& (2.31)\end{array}$

To detect any pitchfork bifurcation points on the symmetric steady-state branch parameterized by $\rho =d_1^v/d_1^u=d_2^v/d_2^u$ we numerically solve (Equation 2.20) together with Equation (2.30). In Section 4 we illustrate this approach for certain reaction kinetics when Ω is the unit disk. The advantage of considering a disk-shaped confining domain is that the reduced-wave Green's function is known analytically by using separation of variables (see Appendix B). We remark that it would also be readily feasible to illustrate our asymptotic theory for a rectangular-shaped confining domain, since the reduced-wave Green's function is also available analytically for such a domain.

3. The linear stability analysis

In this section, we formulate the linear stability problem for the steady-state solutions constructed in Section 2.2. We denote the bulk steady-state solutions of Section 2.2 by u_e(x) and v_e(x), and the steady-state vector of intracellular steady-states by ${\mu }_e={({\mu }_{e1},\dots ,{\mu }_{em})}^T$ and ${\eta }_e={({\eta }_{e1},\dots ,{\eta }_{em})}^T$.

To formulate the linear stability problem, we first introduce the perturbations

$\begin{array}{l}u(t,\text{x})=u_e(\text{x})+e^{\lambda t}\varphi (\text{x}),v(t,\text{x})=v_e(\text{x})+e^{\lambda t}\psi (\text{x}),\\ {\mu }_j(t)={\mu }_{ej}+e^{\lambda t}{\xi }_j,{\eta }_j(t)={\eta }_{ej}+e^{\lambda t}{\zeta }_j,\text{for\hspace{1em}}j\in \left\{1,\dots ,m\right\},\end{array}$

into Equation (2.3) and linearize the resulting system. This yields the eigenvalue problem

$\begin{array}{lll}\text{bulk}& \{\begin{array}{ll}\Delta \varphi -{\Omega }_u^2\varphi =0,& \text{x}\in \Omega \backslash {\displaystyle {\bigcup }_{j=1}^m}{\Omega }_j,\\ \Delta \psi -{\Omega }_v^2\psi =0,& \text{x}\in \Omega \backslash {\displaystyle {\bigcup }_{j=1}^m}{\Omega }_j,\\ {\partial }_{\tilde{n}}\varphi ={\partial }_{\tilde{n}}\psi =0,& \text{x}\in \partial \Omega ,\end{array}& (\mathrm{3.1a})\end{array}$

$\begin{array}{lll}\text{reaction fluxes}& \{\begin{array}{cc}\epsilon D_u{\partial }_{n_j}\varphi =d_1^u\varphi -d_2^u{\xi }_j,& \text{x}\in \partial {\Omega }_j,\\ \epsilon D_v{\partial }_{n_j}\psi =d_1^v\psi -d_2^v{\zeta }_j,& \text{x}\in \partial {\Omega }_j,\end{array}& (\mathrm{3.1b})\end{array}$

$\begin{array}{lll}\text{compartments}& \{\begin{array}{l}(\lambda I-J_j)\left(\begin{array}{c}{\xi }_j\\ {\zeta }_j\end{array}\right)={\epsilon }^{-1}\left(\begin{array}{c}{\int }_{\partial {\Omega }_j}(d_1^u\varphi -d_2^u{\xi }_j)dS\\ {\int }_{\partial {\Omega }_j}(d_1^v\psi -d_2^v{\zeta }_j)dS\end{array}\right),\\ j\in \left\{1,\dots ,m\right\}.\end{array}& (\mathrm{3.1c})\end{array}$

Here the Jacobian matrix J_j of the intracellular kinetics, as well as Ω_u and Ω_v are defined by

$\begin{array}{l}{J}_j\equiv {\left(\begin{array}{cc}{\partial }_{\mu }f(\mu ,\eta )& {\partial }_{\eta }f(\mu ,\eta )\\ {\partial }_{\mu }g(\mu ,\eta )& {\partial }_{\eta }g(\mu ,\eta )\end{array}\right)|}_{\mu ={\mu }_{ej},\eta ={\eta }_{ej}},\\ {\Omega }_u\equiv \sqrt{\frac{\lambda +{\sigma }_u}{D_u}},{\Omega }_v\equiv \sqrt{\frac{\lambda +{\sigma }_v}{D_v}}.& (3.2)\end{array}$

We now use strong localized perturbation theory [18] to analyze (Equation 3.1) in the limit ε → 0. In this way we will derive a nonlinear matrix eigenvalue problem, referred to as the globally coupled eigenvalue problem (GCEP), for the discrete eigenvalues λ of the linearization. This GCEP will be used to investigate various instabilities of the steady-state solutions constructed in Section 2.2.

In the ${O}(\epsilon )$ inner region near the j^th cell we introduce the local variables ${\text{y}}_j={\epsilon }^{-1}(\text{x}-{\text{x}}_j)$, ϕ_j(x) ≡ ϕ(x_j + εy_j) and ψ_j(x) ≡ ψ(x_j + εy_j), with p_j = |y_j|. Upon writing (Equation 3.1a) in terms of the inner variables, for ε → 0 we obtain in the j^th inner region that Δϕ_j = 0 and Δψ_j = 0, for p_j ≥ 1, subject to $D_u{\partial }_{p_j}{\varphi }_j=d_1^u{\varphi }_j-d_2^u{\xi }_j$ and $D_v{\partial }_{p_j}{\psi }_j=d_1^v{\psi }_j-d_2^v{\zeta }_j$ on p_j = 1. The radially symmetric solutions to these problems are

$\begin{array}{l}{\varphi }_j(p_j)=c_j^{u}log{p}_j+\frac{1}{d_1^u}\left(D_u{c}_j^u+d_2^u{\xi }_j\right),\\ {\psi }_j(p_j)=c_j^{v}log{p}_j+\frac{1}{d_1^v}\left(D_v{c}_j^v+d_2^v{\zeta }_j\right),& (3.3)\end{array}$

for j ∈ {1, …, m}, where $c_j^u$ and $c_j^v$ for j ∈ {1, …, m} are constants to be determined. Upon substituting (Equation 3.3) into Equation (3.1c) we obtain, in terms of the Jacobian J_j of Equation (3.2), that

$\begin{array}{ll}(\lambda I-J_j)\left(\begin{array}{c}{\xi }_j\\ {\zeta }_j\end{array}\right)=\left(\begin{array}{c}2\pi D_u{c}_j^u\\ 2\pi D_v{c}_j^v\end{array}\right),\text{for\hspace{1em}}j\in \left\{1,\dots ,m\right\}.& (3.4)\end{array}$

To determine $c_j^u$ and $c_j^v$ we must match the far-field behavior of the inner solutions (Equation 3.3) to the outer solutions defined in the bulk region. Similar to the analysis of the steady-state solution, we obtain that

$\begin{array}{l}\Delta \varphi -{\Omega }_u^2\varphi =0,\text{x}\in \Omega \backslash \left\{{\text{x}}_1,\dots ,{\text{x}}_m\right\};{\partial }_n\varphi =0,\text{x}\in \partial \Omega ;\\ U~c_j^{u}log|\text{x}-{\text{x}}_j|+\frac{c_j^u}{\nu }+\frac{1}{d_1^u}(D_u{c}_j^u+d_2^u{\xi }_j),& (\mathrm{3.5a})\end{array}$

$\begin{array}{ll}\text{as\hspace{1em}}\text{x}\to {\text{x}}_j,j\in \left\{1,\dots ,m\right\},& (\mathrm{3.5b})\end{array}$

where ν ≡ −1/log ε ≪ 1. Similarly, for the perturbation of the other bulk species we obtain

$\begin{array}{l}\Delta \psi -{\Omega }_v^2\psi =0,\text{x}\in \Omega \backslash \left\{{\text{x}}_1,\dots ,{\text{x}}_m\right\};{\partial }_n\psi =0,\text{x}\in \partial \Omega ;\\ \psi ~c_j^{v}log|\text{x}-{\text{x}}_j|+\frac{c_j^v}{\nu }+\frac{1}{d_1^v}(D_v{c}_j^v+d_2^v{\zeta }_j),& (\mathrm{3.6a})\end{array}$

$\begin{array}{ll}\text{as\hspace{1em}}\text{x}\to {\text{x}}_j,j\in \left\{1,\dots ,m\right\}.& (\mathrm{3.6b})\end{array}$

The solutions to Equations (3.5), (3.6) are represented as

$\begin{array}{ll}\varphi (x)=-2\pi {\displaystyle \sum _{i=1}^m}{c}_i^u{G}_{u,\lambda }(\text{x};{\text{x}}_i),\psi (x)=-2\pi {\displaystyle \sum _{i=1}^m}{c}_i^v{G}_{v,\lambda }(\text{x};{\text{x}}_i),& (3.7)\end{array}$

where, to simplify the notation and emphasize the dependence on the eigenvalue parameter λ, we have defined

$\begin{array}{ll}{G}_{u,\lambda }(\text{x};{\text{x}}_j)\equiv G_{{\Omega }_u}(\text{x};{\text{x}}_j),G_{v,\lambda }(\text{x};{\text{x}}_j)\equiv G_{{\Omega }_v}(\text{x};{\text{x}}_j),& (3.8)\end{array}$

where G_ω(x; x_j) is defined by the solution to Equation (2.9). Upon letting x → x_j in Equation (3.7) and ensuring that the singularity conditions in Equations (3.5b), (3.6b) are satisfied, we obtain a linear algebraic system for the vectors ${\text{c}}^u\equiv {(c_1^u,\dots ,c_m^u)}^T$ and ${\text{c}}^v\equiv {(c_1^v,\dots ,c_m^v)}^T$, given in matrix form by

$\begin{array}{l}\left(\left(1+\frac{\nu D_u}{d_1^u}\right)I+2\pi \nu {{G}}_{u,\lambda }\right){\text{c}}^u=-\frac{\nu d_2^u}{d_1^u}\xi ,\\ \left(\left(1+\frac{\nu D_v}{d_1^v}\right)I+2\pi \nu {{G}}_{v,\lambda }\right){\text{c}}^v=-\frac{\nu d_2^v}{d_1^v}\zeta ,& (3.9)\end{array}$

where $\xi \equiv {({\xi }_1,\dots ,{\xi }_m)}^T$ and $\zeta \equiv {({\zeta }_1,\dots ,{\zeta }_m)}^T$. In Equation (3), ${{G}}_{u,\lambda }$ and ${{G}}_{v,\lambda }$ denote the reduced-wave Green's matrix given in Equation (2.12) with either ω = Ω_u or ω = Ω_v, respectively. Here Ω_u and Ω_v are defined in terms of λ by Equation (3.2).

Assuming that λ is not an eigenvalue of J_j for any j ∈ {1, …, m}, we obtain upon inverting (Equation 3.4) and writing the system in matrix form that

$\begin{array}{ll}\xi =2\pi D_u{\text{K}}_{11}{\text{c}}^u+2\pi D_v{\text{K}}_{12}{\text{c}}^v,\zeta =2\pi D_u{\text{K}}_{21}{\text{c}}^u+2\pi D_v{\text{K}}_{22}{\text{c}}^v,& (3.10)\end{array}$

where $\xi ={({\xi }_1,\dots ,{\xi }_m)}^T$ and $\zeta ={({\zeta }_1,\dots ,{\zeta }_m)}^T$. Here K₁₁, K₁₂, K₂₁, and K₂₂ are the diagonal matrices defined by

$\begin{array}{l}{\text{K}}_{11}\equiv \text{diag}(K_{11j}),{\text{K}}_{12}\equiv \text{diag}(K_{12j}),\\ {\text{K}}_{21}\equiv \text{diag}(K_{21j}),{\text{K}}_{22}\equiv \text{diag}(K_{22j}),& (\mathrm{3.11a})\end{array}$

with diagonal entries given by

$\begin{array}{l}{K}_{11j}\equiv {\text{e}}_1^T{(\lambda I-J_j)}^{-1}{\text{e}}_1,K_{12j}\equiv {\text{e}}_1^T{(\lambda I-J_j)}^{-1}{\text{e}}_2,\\ K_{21j}\equiv {\text{e}}_2^T{(\lambda I-J_j)}^{-1}{\text{e}}_1,K_{22j}\equiv {\text{e}}_2^T{(\lambda I-J_j)}^{-1}{\text{e}}_2,& (\mathrm{3.11b})\end{array}$

where ${\text{e}}_1={(1,0)}^T$ and ${\text{e}}_2={(0,1)}^T$.

Then, upon substituting (Equation 3) into Equation (3.10), we obtain the 2m × 2m homogeneous algebraic system, which we write in block matrix form as

$\begin{array}{ll}{M}\left(\lambda \right)\left(\begin{array}{c}{\text{c}}^u\\ {\text{c}}^v\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),\text{where\hspace{1em}}{M}(\lambda )\equiv \left(\begin{array}{cc}{{M}}_u(\lambda )& {{H}}_u(\lambda )\\ {{M}}_v(\lambda )& {{H}}_v(\lambda )\end{array}\right),& (\mathrm{3.12a})\end{array}$

with

$\begin{array}{l}{{M}}_u(\lambda )\equiv \left(1+\frac{\nu D_u}{d_1^u}\right)I+2\pi \nu D_u\frac{d_2^u}{d_1^u}{\text{K}}_{11}+2\pi \nu {{G}}_{u,\lambda },\\ {{H}}_u(\lambda )\equiv 2\pi \nu D_v\frac{d_2^u}{d_1^u}{\text{K}}_{12},& (\mathrm{3.12b})\end{array}$

$\begin{array}{l}{{H}}_v(\lambda )\equiv 2\pi \nu D_u\frac{d_2^v}{d_1^v}{\text{K}}_{21},\\ {{M}}_v(\lambda )\equiv \left(1+\frac{\nu D_v}{d_1^v}\right)I+2\pi \nu D_v\frac{d_2^v}{d_1^v}{\text{K}}_{22}+2\pi \nu {{G}}_{v,\lambda }.& (\mathrm{3.12c})\end{array}$

The nonlinear matrix eigenvalue problem (Equation 3.12) is referred to as the globally coupled eigenvalue problem (GCEP). The GCEP has a nontrivial solution (c^u, c^v)T ≠ (0, 0)^T, if and only λ satisfies $det{M}(\lambda )=0$. The set $\Lambda ({M})$, defined by

$\begin{array}{ll}\Lambda ({M})\equiv \left\{\lambda |det{M}(\lambda )=0\right\},& (3.13)\end{array}$

is the union of all such roots. Any element $\lambda \in \Lambda ({M})$ satisfying Re(λ) > 0 provides an approximation, valid as ε → 0, for an unstable discrete eigenvalue of the linearized problem (Equation 3.1). However, if for all $\lambda \in \Lambda ({M})$ we have Re(λ) < 0, then the steady-state solution is linearly stable.

When there is a large number of cells m, the determination of the discrete eigenvalues comprising $\Lambda ({M})$ is in general a very challenging numerical problem. A survey of nonlinear matrix eigenvalue problems and available solution strategies that apply to only certain classes of matrices is given in Güttel and Tisseur [47] and Betcke et al. [48]. Specific applications of nonlinear matrix problems in simpler contexts where ${M}(\lambda )$ is either a polynomial or a rational function of λ are discussed in Betcke et al. [49]. Since for our problem, ${M}(\lambda )$ is not symmetric and has a complicated dependence on the eigenvalue parameter through the Green's matrices and from the diagonal K matrices of Equation (3.11), these previously developed numerical strategies are not applicable for computing the set $\Lambda ({M})$ in Equation (3.13) for a steady-state solution with an arbitrary collection of cells.

For a symmetric steady-state solution corresponding to a symmetric cell arrangement, we now verify that the condition $det({M}(0))=0$ in Equation (3.12) is equivalent to the zero-eigenvalue crossing condition derived in Equation (2.22), which was based on a linearization of the nonlinear algebraic system around the symmetric steady-state. When λ = 0, we have

$\begin{array}{l}2\pi D_u{\text{c}}^u=-{\Theta }_u\xi ,2\pi D_v{\text{c}}^v=-{\Theta }_v\zeta ,\end{array}$

where Θ_u and Θ_v are defined in Equation (2.13b). Since K₁₁, K₁₂, K₂₂, and K₂₁ are all multiples of the identity for a symmetric steady-state, we obtain from Equation (3.12) together with Equation (3) that when λ = 0 we have in block matrix form

$\begin{array}{ll}\left(\begin{array}{c}\xi \\ \zeta \end{array}\right)+\left(\begin{array}{cc}{K}_{11}I& K_{12}I\\ K_{12}I& K_{22}I\end{array}\right)\left(\begin{array}{cc}{\Theta }_u& 0\\ 0& {\Theta }_v\end{array}\right)\left(\begin{array}{c}\xi \\ \zeta \end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).& (3.14)\end{array}$

Here we have re-defined the scalars K₁₁, K₁₂, K₂₁, and K₂₂ by K₁₁ = K₁₁I, K₁₂ = K₁₂I, K₂₁ = K₂₁I, and K₂₂ = K₂₂I. When λ = 0, we calculate that

$\begin{array}{l}{K}=\left(\begin{array}{cc}{K}_{11}I& K_{12}I\\ K_{21}I& K_{22}I\end{array}\right)=-J_c^{-1},\text{where\hspace{1em}}{J}_c\equiv \left(\begin{array}{cc}{f}_{\mu }^{c}I& f_{\eta }^{c}I\\ g_{\mu }^{c}I& g_{\eta }^{c}I\end{array}\right).\end{array}$

Finally, upon multiplying (Equation 3.14) by J_c, and using $J_c{K}=-I$, we readily obtain that

$\begin{array}{ll}\left(\begin{array}{cc}{f}_{\mu }^{c}I-{\Theta }_u& f_{\eta }^{c}I\\ g_{\mu }^{c}I& g_{\eta }^{c}I-{\Theta }_v\end{array}\right)\left(\begin{array}{c}\xi \\ \zeta \end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right),& (3.15)\end{array}$

which is precisely the same as in Equation (2.22).

3.1. Re-formulation of the linear stability problem

A simpler formulation of the linear stability problem that applies to both symmetric and asymmetric steady-state solutions can be done when g has the specific form in Equation (2.14). In this situation, we can write (Equation 3.4) in the form

$\begin{array}{l}{\zeta }_j=\frac{g_1^{\prime }({\mu }_j)}{\lambda +g_2}{\xi }_j+\frac{2\pi D_v}{\lambda +g_2}{c}_j^v,\left(\lambda -f_{\mu }({\mu }_j,{\eta }_j)\right){\xi }_j-f_{\eta }({\mu }_j,{\eta }_j){\zeta }_j=2\pi D_u{c}_j^u.\end{array}$

Then, upon relating $c_j^v$ and $c_j^u$ to ζ_j and ξ_j by using Equation (3), we obtain in matrix form that

$\begin{array}{ll}\zeta =\frac{1}{\lambda +g_2}{\Lambda }_2\xi -\frac{1}{\lambda +g_2}{\Theta }_{v,\lambda }\zeta ,{\Lambda }_3\xi -{\Lambda }_4\zeta =-{\Theta }_{u,\lambda }\xi ,& (3.16)\end{array}$

Where Θ_u,λ,Θ_v,λ, and the diagonal matrices Λ₂, Λ₃, and Λ₄ are defined by

$\begin{array}{l}{\Theta }_{u,\lambda }\equiv 2\pi \nu D_u\frac{d_2^u}{d_1^u}{\left[\left(1+\frac{\nu D_u}{d_1^u}\right)I+2\pi \nu {{G}}_{u,\lambda }\right]}^{-1},\\ {\Theta }_{v,\lambda }\equiv 2\pi \nu D_v\frac{d_2^v}{d_1^v}{\left[\left(1+\frac{\nu D_v}{d_1^v}\right)I+2\pi \nu {{G}}_{v,\lambda }\right]}^{-1},& (\mathrm{3.17a})\end{array}$

$\begin{array}{l}{\Lambda }_2\equiv \text{diag}(g_1^{\prime }({\mu }_j)),{\Lambda }_3\equiv \text{diag}(\lambda -f_{\mu }({\mu }_j,{\eta }_j)),\\ {\Lambda }_4\equiv \text{diag}(f_{\eta }({\mu }_j,{\eta }_j)).& (\mathrm{3.17b})\end{array}$

Upon eliminating ζ in Equation (3.16), we obtain the following nonlinear eigenvalue problem for the case where g has the form in Equation (2.14):

$\begin{array}{ll}{N}(\lambda )\xi =0,\text{where }{N}(\lambda )\equiv {\Lambda }_3-{\Lambda }_4{\left[(\lambda +g_2)I+{\Theta }_{v,\lambda }\right]}^{-1}{\Lambda }_2+{\Theta }_{u,\lambda }.& (3.18)\end{array}$

Observe that setting $det({N}(\lambda ))=0$ involves root-finding on the determinant of a matrix of size m × m rather than that for the larger 2m × 2m matrix, as needed for Equation (3.13).

The characterization (Equation 3.18) is particularly useful for determining the linear stability properties of a symmetric steady-state for a symmetric cell arrangement when g has the form in (2.14). For a symmetric steady-state with λ = 0, we obtain from Equation (3.17) that ${\Lambda }_2=g_1^{\prime }({\mu }_c)I$, ${\Lambda }_3=-f_{\mu }^{c}I$, and ${\Lambda }_4=f_{\eta }^{c}I$. From Equation (3.18), this yields that

$\begin{array}{l}{N}(0)=-f_{\mu }^{c}I-f_{\eta }^c{g}_1^{\prime }({\mu }_c){\left[g_{2}I+{\Theta }_{v,0}\right]}^{-1}+{\Theta }_{u,0}.\end{array}$

Since Θ_u,λ = Θ_u and Θ_v,λ = Θ_v when λ = 0, where Θ_u and Θ_v were defined in Equation (2.13b), we obtain that the condition $det({N}(0))=0$ is equivalent to the formulation (Equation 2.23) derived in Section 2.3, which was based on linearizing the nonlinear algebraic system around the symmetric steady-state solution.

For a symmetric steady-state solution of a symmetric cell arrangement, one key advantage of the re-formulation (Equation 3.18) is that the discrete eigenvalues of the linearization (Equation 3.1) can be determined by finding the union of the roots of m scalar problems. This is done by using $det({N}(\lambda ))=\prod _{j=1}^m{\sigma }_j(\lambda )$, where σ_j(λ) for j = {1, …, m} are the eigenvalues of ${N}(\lambda )$. More specifically, since ${{G}}_{u,\lambda }$ and ${{G}}_{v,\lambda }$ have the common eigenspace

$\begin{array}{l}{{G}}_{u,\lambda }\text{e}={\kappa }_{u,\lambda }\text{e},{{G}}_{v,\lambda }\text{e}={\kappa }_{v,\lambda }\text{e};{{G}}_{u,\lambda }{\text{q}}_j={\kappa }_{u,\lambda j}^{\perp }{\text{q}}_j,\\ {{G}}_{v,\lambda }{\text{q}}_j={\kappa }_{v,\lambda j}^{\perp }{\text{q}}_j,j\in \left\{2,\dots ,m\right\},& (\mathrm{3.19a})\end{array}$

the eigenvalue σ₁(λ) corresponding to e and the eigenvalues σ_j(λ) corresponding to q_j, for j ∈ {2, …, m} are readily calculated. A simple calculation yields that

$\begin{array}{ll}{\sigma }_1(\lambda )=\lambda -f_{\mu }^c-\frac{f_{\eta }^c{g}_1^{\prime }({\mu }_c)}{\lambda +g_2+{\alpha }_{v,\lambda }}+{\alpha }_{u,\lambda },& (\mathrm{3.20a})\end{array}$

$\begin{array}{l}{\sigma }_j(\lambda )=\lambda -f_{\mu }^c-\frac{f_{\eta }^c{g}_1^{\prime }({\mu }_c)}{\lambda +g_2+{\alpha }_{v,\lambda j}^{\perp }}+{\alpha }_{u,\lambda j}^{\perp },\\ j\in \left\{2,\dots ,m\right\},& (\mathrm{3.20b})\end{array}$

where we have defined

$\begin{array}{ll}{\alpha }_{u,\lambda }\equiv \frac{2\pi \nu D_u{d}_2^u/d_1^u}{1+\nu D_u/d_1^u+2\pi \nu {\kappa }_{u,\lambda }},{\alpha }_{v,\lambda }\equiv \frac{2\pi \nu D_v{d}_2^v/d_1^v}{1+\nu D_v/d_1^v+2\pi \nu {\kappa }_{v,\lambda }},& (\mathrm{3.20c})\end{array}$

$\begin{array}{l}{\alpha }_{u,\lambda j}^{\perp }\equiv \frac{2\pi \nu D_u{d}_2^u/d_1^u}{1+\nu D_u/d_1^u+2\pi \nu {\kappa }_{u,\lambda j}^{\perp }},{\alpha }_{v,\lambda j}^{\perp }\equiv \frac{2\pi \nu D_v{d}_2^v/d_1^v}{1+\nu D_v/d_1^v+2\pi \nu {\kappa }_{v,\lambda j}^{\perp }},\\ j\in \left\{2,\dots ,m\right\}.& (\mathrm{3.20d})\end{array}$

With this re-formulation, for a symmetric steady-state of a symmetric cell arrangement, and with g of the form in Equation (2.14), the set of discrete eigenvalues of the linearization, $\Lambda ({M})$, in Equation (3.13) can be written conveniently as

$\begin{array}{ll}\Lambda ({M})\equiv \left\{\lambda |{\sigma }_1(\lambda )=0\right\}\cup {\displaystyle \bigcup _{j=2}^m}\left\{\lambda |{\sigma }_j(\lambda )=0\right\}.& (3.21)\end{array}$

In summary, to determine the linear stability properties of this symmetric steady-state solution we need only solve m scalar root-finding problems and determine whether there are any roots in Re(λ) > 0. This is considerably more tractable numerically than performing a root-finding based on the determinant of the GCEP in Equation (3.18).

4. Illustrations of the theory: A ring pattern of cells

In this section we illustrate the steady-state and linear stability theory developed in Sections 2.2, 3 for a ring pattern of cells inside unit disk Ω, for which the Green's function is given analytically in Appendix B. We will show that symmetry-breaking bifurcations can occur for the Rauch and Millonas [8], Gomez-Marin et al. [24], and Gierer and Meinhardt [46] reaction kinetics. The theoretical prediction of stable asymmetric patterns will be confirmed through full time-dependent numerical simulations of Equation (2.3) computed using FlexPDE [50].

For a ring pattern of m cells in the unit disk with ring radius r, with 0 < r < 1, the cell centers are located at

$\begin{array}{ll}{\text{x}}_k=r\left(cos\left(\frac{2\pi (k-1)}{m}\right),sin\left(\frac{2\pi (k-1)}{m}\right)\right),k\in \left\{1,\dots ,m\right\}.& (4.1)\end{array}$

For a ring pattern of cells, all Green matrices are symmetric and circulant and have the common eigenspace

$\begin{array}{l}{\text{v}}_k={(1,Z_k,Z_k^2,...,Z_k^{m-1})}^T\text{with\hspace{1em}}{Z}_k\equiv \text{exp}\left(\frac{2\pi i(k-1)}{m}\right)\\ \text{and\hspace{1em}}k\in \left\{1,...,m\right\},\end{array}$

which are a basis of ℂ^m. In Equation (B.2) of Appendix B we summarize how to obtain the matrix spectrum of a symmetric and circulant matrix that has a real-valued basis for the eigenspace.

In our illustrations of the theory below, we will assume for simplicity that the membrane reaction rates satisfy

$\begin{array}{ll}{d}_1^u=d_2^u\equiv d_u{d}_1^v=d_2^v\equiv d_v,\text{with\hspace{1em}}\rho \equiv \frac{d_v}{d_u},& (4.2)\end{array}$

and that g(μ, η) has the specific form in Equation (2.14). We will focus on a two-cell ring pattern in the unit disk, as shown schematically in Figure 2, for three specific reaction kinetics.

FIGURE 2

Figure 2. A schematic plot of a ring pattern in the unit disk with two cells. The bifurcation parameter for symmetry-breaking is ρ, while the diffusivities satisfy D_u = D_v.

To numerically implement our asymptotic theory, the steady-state solution branches are computed from Equation (2.16) with m = 2 using the parameter continuation software MatCont [51], while the symmetric solution branch is obtained from Equation (2.20). Symmetry-breaking bifurcation points in ρ along the symmetric branch are identified by numerically solving (Equation 2.30) together with Equation (2.20). Finally, to determine the linear stability properties of the symmetric branch we need only determine if there exists a λ with Re(λ) > 0 in the set $\Lambda ({M})$ given in Equation (3.21). For m = 2, this is done by calculating all the roots of σ₁(λ) = 0 and σ₂(λ) = 0 by using Equation (3.20) and the explicit expressions for the eigenvalues of the Green's matrices as can be obtained from Appendix B.

4.1. Gierer-Meinhardt reaction kinetics

We consider a prototypical Gierer-Meinhardt model (GM), where the nonlinear reaction kinetics are confined within the compartments. The original GM model, introduced in Gierer and Meinhardt [46] and Gierer [52] to model pattern formation in biological morphogenesis, has the form

$\begin{array}{l}{\partial }_{t}u=D_u\Delta u-{\sigma }_{u}u+{\varrho }_0(x)+c_u{\varrho }_u(x)\frac{u^2}{v},\\ {\partial }_{t}v=D_v\Delta v-{\sigma }_{v}v+c_v{\varrho }_v(x)u^2.\end{array}$

This model disregards that in biological tissues morphogen-producing reactions mostly occur intracellularly and on the membranes of cells. For simplicity, we illustrate our compartmental-reaction diffusion theory for the specific case where ϱ₀ ≡ 0, c_uϱ_u(x) ≡ 1, and c_vϱ_v ≡ 1. When there is no bulk diffusion, the compartments are uncoupled from the bulk and we impose the intracellular reaction kinetics

$\begin{array}{ll}\stackrel{.}{\mu }(t)=f(\mu ,\eta )\equiv \frac{{\mu }^2}{\eta },\stackrel{.}{\eta }(t)=g(\mu ,\eta )\equiv {\mu }^2.& (4.3)\end{array}$

The uncoupled equilibrium for Equation (4.3) given by μ_e = 0, and where η_e is an arbitrary constant, is non-hyberbolic in all directions.

To apply the bulk-cell steady-state analysis of Section 2 for a two-cell ring pattern, we first identify that g(μ, η) = g₁(μ) − g₂η, where $g_1={\mu }^2$ and g₂ = 0. For m = 2, we conclude from Equation (2.16) that all steady-states of the bulk-cell system are associated with the nonlinear algebraic problem

$\begin{array}{ll}\begin{array}{c}f({\mu }_{e1},{\text{e}}_1^T{\Theta }_v^{-1}{({({\mu }_{e1})}^2,{({\mu }_{e2})}^2)}^T)-{\text{e}}_1^T{\Theta }_u{({\mu }_{e1},{\mu }_{e2})}^T=0\\ f({\mu }_{e2},{\text{e}}_2^T{\Theta }_v^{-1}{({({\mu }_{e1})}^2,{({\mu }_{e2})}^2)}^T)-{\text{e}}_2^T{\Theta }_u{({\mu }_{e1},{\mu }_{e2})}^T=0.\end{array}& (4.4)\end{array}$

The symmetric equilibrium (μ_e, η_e), which satisfies (Equation 2.20), is readily calculated as

$\begin{array}{ll}{\mu }_e=\frac{{\alpha }_v}{{\alpha }_u},{\eta }_e=\frac{{\alpha }_v}{{\alpha }_u^2},& (4.5)\end{array}$

where α_u and α_v are defined in Equation (2.19). By combining (Equation 4.5) with Equation (2.30), we conclude that a symmetry-breaking bifurcation from the symmetric steady-state occurs whenever the condition

$\begin{array}{ll}\frac{{\alpha }_v}{{\alpha }_{v,2}^{\perp }}+\frac{{\alpha }_{u,2}^{\perp }}{2{\alpha }_u}-1=0,& (4.6)\end{array}$

is satisfied at some point along the symmetric solution branch. Here ${\alpha }_{u,2}^{\perp }$ and ${\alpha }_{v,2}^{\perp }$ were defined in (2.31).

In Figure 3, Left we plot the bifurcation diagram of solutions to Equation (4.4) for a parameter set where D_v = D_u and with the other parameter values as in the figure caption. We observe that a supercritical symmetry-breaking pitchfork bifurcation from the symmetric branch occurs at the critical value ρ = ρ_p ≈ 9.79168. In Figure 4 we show full PDE results for Equation (2.3) computed with FlexPDE [50] for values of ρ on either side of this theoretically predicted bifurcation value. In the left panels of this figure, we observe that when ρ = 5 < ρ_p, an initial perturbation of the symmetric steady-state converges back to the symmetric steady-state as time increases. In contrast, when ρ = 15 > ρ_p, we observe from the Figure 4, Right that, for an initial condition near the symmetric steady-state, the time-dependent PDE solution converges as time increases to the asymmetric steady-state predicted in Figure 3, Left.

FIGURE 3

Figure 3. Left: 3-D Bifurcation diagram, computed from Equation (4.4), showing symmetric and asymmetric steady-states of a two-cell ring pattern with ring radius r = 0.5 and GM kinetics (Equation 4.3). Asymmetric steady-states emerge at the supercritical pitchfork bifurcation point ρ = ρ_p ≈ 9.79168 along the symmetric branch. Right: The pitchfork bifurcation value of ρ increases rapidly as the ring radius r, and consequently the distance between the cells, increases. The dots are the values computed from Equation (4.6), while the curve is the interpolation by the plotting function in Bezanson et al. [53]. Parameters: D_u = D_v = 5, σ_u = σ_v = 0.6, d_u = 0.09, and ε = 0.03.

FIGURE 4

Figure 4. Full numerical PDE simulation results of Equation (2.3) with FlexPDE [50] for GM kinetics Equation (4.3). The bottom panels show the concentration of u. Left: convergence to the symmetric branch for ρ = 5 before the supercritical pitchfork point ρ_p = 9.79168, for an initial condition close to the symmetric branch. Right: convergence to the asymmetric branch selected by the eigenperturbation direction ${\text{q}}_2={(1,-1)}^T$ for ρ = 15 and starting near the symmetric branch. Parameters: D_u = D_v = 5, σ_u = σ_v = 0.6, d_u = 0.09, ε = 0.03, and r = 0.5.

In the Figure 3, Right we show that the pitchfork bifurcation value for the emergence of asymmetric steady-states increases substantially when the ring radius for the two-cell pattern increases. As a result, we conclude that for cells that are farther apart, a larger value of ρ is needed to create a stable asymmetric pattern.

We now show that by varying the membrane reaction rate d_u, which necessarily varies the membrane reaction rate to the v-species according to d_v = ρd_u, the steady-state solution branches with GM kinetics (Equation 4.3) computed from Equation (4.4) can exhibit a hysteresis structure for low d_u. The numerical results of Figure 5 show such a hysteretic bifurcation structure between the asymmetric and symmetric solution branches for two values of d_u. We observe that as d_u decreases the extent of the hysteresis increases. The range where hysteresis occurs is given by the separation ρ_p − ρ_s between the pitchfork point ρ_p and the secondary fold bifurcation point ρ_s along each asymmetric branch. Numerical results for this range for a parameter set where hysteresis occurs when d_u < 0.09 is given in Table 1.

FIGURE 5

Figure 5. 3-D Bifurcation diagram, computed from Equation (4.4), for symmetric and asymmetric steady-states of a two-cell ring pattern with ring radius r = 0.5 and two different values of d_u with GM kinetics (4.3). Left: d_u = 0.08. Right: d_u = 0.05. For these values of d_u, the steady-states exhibit hysteresis, i.e., a subcritical pitchfork bifurcation occurs from the symmetric equilibrium branch, with the emerging unstable asymmetric equilibrium branches regaining stability at a secondary fold point. Observe that the extent of the hysteresis increases when d_u decreases. Parameters: D_u = D_v = 5, σ_u = σ_v = 0.6, ε = 0.03, and r = 0.5.

TABLE 1

Table 1. Numerical values (rounded to 5th decimal place) of the subcritical (or supercritical) pitchfork bifurcation point ρ_p, the fold bifurcation point ρ_s, and the associated values for the symmetric μ_e and one of the asymmetric (μ_e1, μ_e2) solution branches.

For the parameter set with d_u = 0.05, which corresponds to the bifurcation diagram shown in the Figure 5, Right, the full time-dependent computations of Equation (2.3) with FlexPDE [50], as shown in Figure 6, illustrate that for an initial condition near the symmetric steady-state branch, and with ρ either satisfying ρ < ρ_s or ρ_s < ρ < ρ_p, the time-dependent solution converges to the stable symmetric steady-state solution. However, as shown in Figure 7, Left, for an initial condition near the asymmetric branch when ρ is in the hysteresis region ρ_s < ρ < ρ_p, the time-dependent solution converges to the asymmetric steady-state. Moreover, if ρ > ρ_p, the Figure 7, Right shows that for an initial condition near the unstable symmetric steady-state the time-dependent solution converges to the asymmetric steady-state solution.

FIGURE 6

Figure 6. Full numerical simulation results of Equation (2.3) with FlexPDE [50] for GM kinetics (Equation 4.3). Left: for an initial condition near the symmetric branch we observe convergence to the symmetric branch when ρ = 4, which is before the hysteresis region bounded by the fold point ρ_s ≈ 6.27945 and the subcritical pitchfork point ρ_p ≈ 7.70971. Right: convergence to the symmetric branch for ρ = 7.2, which lies on the range ρ_s < ρ < ρ_p, when starting near the symmetric branch. The bottom panels show the concentration of u. Parameters: D_u = D_v = 5, σ_u = σ_v = 0.6, d_u = 0.05, ε = 0.03, r = 0.5.

FIGURE 7

Figure 7. Full numerical simulation results of Equation (2.3) with FlexPDE [50] for GM kinetics (Equation 4.3). Left: convergence to the asymmetric branch for an initial condition near this branch when ρ = 7.2 lies in the hysteresis region between the fold point ρ_s ≈ 6.27945 and the subcritical pitchfork point ρ_p ≈ 7.70971. Right: convergence to an asymmetric steady-state as selected by a small initial perturbation of the symmetric solution in the direction q₂ = (1, −1) when ρ = 15 > ρ_p. The bottom panels show the concentration of u. Parameters: D_u = D_v = 5, σ_u = σ_v = 0.6, d_u = 0.05, ε = 0.03, and r = 0.5.

To determine the linear stability properties of the symmetric steady-state solution branch as ρ is varied in Figures 3, 5 we must determine the eigenvalues λ in the set (3.21). This is done by numerically computing the largest roots to σ₁(λ) = 0 and to σ₂(λ) = 0, as defined in Equation (3.20). In Figure 8 we plot these roots vs. ρ for two values of d_u. From this figure, we observe that in-phase perturbations of the symmetric steady-state solution branch, as determined by the roots of σ₁(λ) = 0, are always linearly stable. In contrast, anti-phase perturbations of the symmetric steady-state, as characterized by the roots of σ₂(λ) = 0, are linearly stable only for ρ < ρ_p, where ρ_p is the symmetry-breaking threshold. For ρ > ρ_p, the symmetric steady-state solution branch is unstable to an anti-phase eigenperturbation ${\text{q}}_2={(1,-1)}^T$.

FIGURE 8

Figure 8. Plots of the numerically computed largest roots of σ₁(λ) = 0 and σ₂(λ) = 0 vs. ρ, as defined in Equation (3.20), that determine the linear stability properties to either in-phase e = (1, 1)^T or anti-phase ${\text{q}}_2={(1,-1)}^T$ eigenperturbations of the symmetric steady-state solution, respectively. Left: for d_u = 0.09 we have ρ_p ≈ 9.79168. Right: for d_u = 0.05 we have ρ_p ≈ 7.70971. Observe that in-phase eigenperturbations are always linearly stable, whereas anti-phase eigenperturbations are linearly stable only on the range ρ < ρ_p before the pitchfork point ρ_p. Parameters: D_u = D_v = 5, σ_u = σ_v = 0.6, ε = 0.03, r = 0.5.

Next, to study the linear stability properties of the asymmetric steady-state solution branches we must determine whether the root-finding condition $det({N}(\lambda ))=0$ in Equation (3.18) yields an eigenvalue with Reλ > 0. The numerical results shown in Figure 9 for d_u = 0.05 (corresponding to in Figure 5, Right), establishes that the asymmetric branch on the subcritical range ρ_s < ρ < ρ_p, which emanates from the symmetric steady-state branch, is unstable. However, as observed from Figure 9, the upper portion of the asymmetric branch for ρ > ρ_s is linearly stable.

FIGURE 9

Figure 9. Zero-crossings of $det({N}(\lambda ))=0$, as defined in Equation (3.18), determine the linear stability properties of an asymmetric steady-state solution with two cells. On the range ρ_p < ρ < ρ_s, before the secondary fold point along the asymmetric branch, we observe that λ > 0. This implies that the subcritical portion of the asymmetric steady-state branch between the pitchfork point and the fold point is unstable. Further along past the fold point the asymmetric branch regains stability.
Parameters: D_u = D_v = 5, σ_u = σ_v = 0.6, d_u = 0.05, ε = 0.03, r = 0.5.

We now explore how the pitchfork bifurcation point depends on decreasing values of the diffusion coefficient ratio D_v/D_u when d_u = 0.09. When this ratio is unity, there was no hysteresis between the symmetric and asymmetric steady-state solution branches (see Table 1). We remark that a similar numerical experiment was performed in Section 2.3 of Pelz and Ward [23] for a 1-D compartmental-reaction diffusion model with GM kinetics with dynamically reactive boundaries. In our 2-D setting, we observe from the numerical results in Table 2 that a symmetry-breaking bifurcation can occur on the range D_v/D_u < 1, but only up until some minimum diffusion ratio threshold at which the pitchfork bifurcation point given by the root of Equation (4.6) no longer exists. In addition, we observe that reducing the diffusion ratio threshold D_v/D_u below unity for fixed d_u = 0.09 does not introduce new hysteresis behavior, and the symmetry-breaking bifurcation remains supercritical.

TABLE 2

Table 2. Decreasing the ratio D_v/D_u does not trigger hysteresis when d_u = 0.09, but rather there is a minimum threshold of the diffusivity ratio where the symmetry-breaking pitchfork bifurcation point exists.

Next, we set d_u = 0.08 where hysteresis occurs when D_v/D_u = 1, and we vary this diffusion ratio to determine whether hysteresis can be eliminated. Our numerical results in Figure 10 and Table 3 indicates that varying D_v/D_u does not eliminate the hysteresis between the symmetric and asymmetric steady-state branches. However, the extent of the hysteresis decreases as the ratio D_v/D_u increases.

FIGURE 10

Figure 10. Effect of the diffusivity ratio D_v/D_u of the two bulk species on the extent of the hysteresis when d_u = 0.08, as measured by the distance between the fold bifurcation points and the subcritical pitchfork bifurcation point (Left) and by the distance of the two asymmetric equilibria μ_e1 and μ_e2 from each other (Right). The diffusivity D_u = 5 is fixed and the remaining parameters are as in Table 2. The dots are the numerically computed values using MatCont [51] that are interpolated by the plotting function in Julia [53].

TABLE 3

Table 3. Increasing the diffusivity ratio D_v/D_u when d_u = 0.08 does not eliminate hysteresis, as the symmetry-breaking bifurcation point remains subcritical.

To obtain some analytical insight into the disappearance of the pitchfork point as shown in Tables 2, 3 when the diffusivity ratio D_v/D_u decreases below a threshold, in Figure 11 we plot the function $F_{\alpha }(\rho )\equiv {\alpha }_v/{\alpha }_{v,2}^{\perp }+{\alpha }_{u,2}^{\perp }/(2{\alpha }_u)-1$ vs. ρ, representing the left-hand side of the pitchfork bifurcation condition (Equation 4.6), for several values of D_v/D_u, and for either d_u = 0.08 (Figure 11, Left) or d_u = 0.09 (Figure 11, Right). From Equation (4.6) a root of F_α(ρ) = 0 corresponds to a symmetry-breaking bifurcation point along the symmetric solution branch. From Figure 11 we observe that the asymptote of F_α(ρ) as ρ → ∞ is positive when D_v/D_u is below a threshold, which eliminates the possibility of a pitchfork bifurcation point.

FIGURE 11

Figure 11. Effect of the diffusivity ratio D_v/D_u on the existence of the pitchfork point when d_u = 0.08 (Left) and d_u = 0.09 (Right). The numerical results show that the asymptote of $F_{\alpha }(\rho )\equiv {\alpha }_v/{\alpha }_{v,2}^{\perp }+{\alpha }_{u,2}^{\perp }/(2{\alpha }_u)-1$ is positive for smaller values of D_v/D_u, as suggested by Tables 2, 3. Therefore, when D_v/D_u falls below a threshold, the pitchfork bifurcation condition F_α(ρ) = 0, which is equivalent to Equation (4.6), no longer holds for any ρ > 0.

4.2. Rauch-Millonas reaction kinetics

Next, we consider the activator-inhibitor system proposed in Rauch and Millonas [8] to universally model two-species signal transduction reaction kinetics between cells. This Rauch-Millonas (RM) intracellular kinetics of Rauch and Millonas [8] is given by

$\begin{array}{ll}\begin{array}{ccc}\stackrel{.}{\mu }=c_u-q_u\mu +\frac{a_1^u\mu }{b_1^u+\mu }-\frac{a_2^u\mu \eta }{b_2^u+\mu }& \equiv & f(\mu ,\eta )\\ \stackrel{.}{\eta }=c_v+w_v\mu -q_v\eta \hfill & \equiv & g(\mu ,\eta ).\end{array}& (4.7)\end{array}$

Since g has the form in Equation (2.14), we identify that g₁(μ) = c_v + w_vμ and g₂ = q_v. We will choose a parameter set for which the reaction kinetics when uncoupled from the bulk has a unique linearly stable steady-state.

From Equation (2.16), all steady-states of the bulk-cell model for a two-cell ring pattern are associated with the nonlinear algebraic system

$\begin{array}{l}f({\mu }_{e1},{\text{e}}_1^T{(q_{v}I+{\Theta }_v)}^{-1}{(c_v+w_v{\mu }_{e1},c_v+w_v{\mu }_{e2})}^T)\\ -{\text{e}}_1^T{\Theta }_u{({\mu }_{e1},{\mu }_{e2})}^T=0\\ f({\mu }_{e2},{\text{e}}_2^T{(q_{v}I+{\Theta }_v)}^{-1}{(c_v+w_v{\mu }_{e1},c_v+w_v{\mu }_{e2})}^T)\\ -{\text{e}}_2^T{\Theta }_u{({\mu }_{e1},{\mu }_{e2})}^T=0.& (4.8)\end{array}$

By using Equation (2.20), the symmetric steady-state solution branch is obtained from the solution μ_e to

$\begin{array}{ll}{c}_u-q_u{\mu }_e+\frac{a_1^u{\mu }_e}{b_1^u+{\mu }_e}-\frac{a_2^u{\mu }_e}{b_2^u+{\mu }_e}\frac{(c_v+w_v{\mu }_e)}{q_v+{\alpha }_v}-{\alpha }_u{\mu }_e=0,& (4.9)\end{array}$

where α_u and α_v are given in Equation (2.19). Symmetry-breaking bifurcation points are obtained by solving the zero-eigenvalue crossing condition (Equation 2.30) together with Equation (4.9). By solving for w_v = w_v(μ_e) in Equation (2.30), we calculate

$\begin{array}{l}{w}_v({\mu }_e)\\ =\frac{-q_u+\frac{a_1^u}{b_1^u+{\mu }_e}-\frac{a_1^u{\mu }_e}{{(b_1^u+{\mu }_e)}^2}-\frac{a_2^u}{b_2^u+{\mu }_e}\frac{c_v}{q_v+{\alpha }_v}+\frac{a_2^u{\mu }_e}{{(b_2^u+{\mu }_e)}^2}\frac{c_v}{q_v+{\alpha }_v}-{\alpha }_{u,2}^{\perp }}{\frac{a_2^u}{b_2^u+{\mu }_e}\frac{{\mu }_e}{q_v+{\alpha }_v}-\frac{a_2^u{\mu }_e}{{(b_2^u+{\mu }_e)}^2}\frac{{\mu }_e}{q_v+{\alpha }_v}+\frac{a_2^u{\mu }_e}{b_2^u+{\mu }_e}\frac{1}{q_v+{\alpha }_{v,2}^{\perp }}},& (4.10)\end{array}$

where ${\alpha }_{u,2}^{\perp }$ and ${\alpha }_{v,2}^{\perp }$ are defined in Equation (2.31). By using Equation (4.10) to eliminate ω_v in Equation (4.9) we obtain a nonlinear algebraic equation that determines any symmetry-breaking bifurcation value for μ_e along the symmetric solution branch. The corresponding bifurcation value for w_v is obtained from Equation (4.10).

For the parameter set given in the figure caption, we show in Figure 12, Left that, for the fixed value ρ = 15, there is a degenerate w_v-pitchfork bubble, which is characterized by the emergence of asymmetric steady-state solutions at two values of w_v. From the Figure 12, Right we observe that in terms of ρ, and at a fixed w_v, the symmetry-breaking bifurcation is supercritical in ρ. For the parameter set in Figure 12, Right, we observe from Figure 13 that the eigenvalue λ determined by the root-finding condition σ₂(λ) = 0, with σ₂ given in Equation (3.20), crosses through zero at the ρ-pitchfork bifurcation point along the symmetric steady-state branch. As a result, when $w_v=w_v^{P,2}\approx 7.08723$, the symmetric steady-state solution is linearly stable for ρ < ρ_p = 15, and is unstable on the range ρ > ρ_p = 15 to eigenperturbations in the direction of ${\text{q}}_2={(1,-1)}^T$.

FIGURE 12

Figure 12. 3-D Bifurcation diagram, computed from Equation (4.8) using MatCont [51], for symmetric and asymmetric steady-states of a two-cell ring pattern with ring radius r = 0.5 and with RM kinetics Equation (4.7). Left: 3-D Plot of (μ_e1, μ_e2) vs. the kinetic parameter w_v in Equation (4.7) at a fixed ρ = d_v/d_u = 15, showing that asymmetric steady-states occur inside the degenerate pitchfork bubble delimited by $w_v^{P,1}\approx 6.88285$ and $w_v^{P,2}\approx 7.08723$. Note that the bubble lobes are stretching into decreasing w_v and that there exists hysteresis at $w_v^{P,2}$. Right: In terms of ρ, a supercritical pitchfork bifurcation from the symmetric branch occurs when $w_v=w_v^{P,2}\approx 7.08723$. The asymmetric branches are linearly stable past this bifurcation threshold in ρ. Parameters: $D_u=D_v=1,$ ${\sigma }_u={\sigma }_v=0.1,$ $\epsilon =0.03,$ $c_u=c_v=1,$ $q_u=1/100,$ $q_v=7,$ $a_1^u=600,$ $a_2^u=6,$ $b_1^u=100,$ $b_2^u=1/10$, and d_u = 0.14.

FIGURE 13

Figure 13. For the two-cell system with RM kinetics (Equation 4.7) and parameters as in the caption of Figure 12 with $w_v=w_v^{P,2}$, we plot the eigenvalue λ, satisfying σ₂(λ) = 0 in Equation (3.20), vs. ρ that determines the linear stability of the symmetric steady-state solution branch to eigenperturbations of the form ${\text{q}}_2={(1,-1)}^T$. We observe that the symmetric steady-state branch is unstable only for ρ > ρ_p = 15. There is no root to σ₁(λ) = 0 on this range.

4.3. FitzHugh-Nagumo reaction kinetics

Finally, we consider a ring pattern for the bulk-cell system (Equation 2.3) with two cells and with FitzHugh-Nagumo (FN) intracellular reaction kinetics [24]. The uncoupled intracellular kinetics are

$\begin{array}{ll}\stackrel{.}{\mu }(t)=\mu -q{(\mu -2)}^3+4-\eta \equiv f(\mu ,\eta ),\stackrel{.}{\eta }(t)=\delta z\mu -\delta \eta \equiv g(\mu ,\eta ),& (4.11)\end{array}$

with q > 0, δ > 0 and z > 0. Since g has the form in Equation (2.14), we identify that g₁(μ) = δz and g₂ = δ.

We will choose a parameter set for which there is a unique linearly stable steady-state of the intra-compartmental dynamics (Equation 4.11). From Equation (2.16), all steady-states of the bulk-cell model for a two-cell ring pattern are obtained from the nonlinear algebraic problem

$\begin{array}{ll}\begin{array}{c}f({\mu }_{e1},\delta z{\text{e}}_1^T{(\delta I+{\Theta }_v)}^{-1}{({\mu }_{e1},{\mu }_{e2})}^T)-{\text{e}}_1^T{\Theta }_u{({\mu }_{e1},{\mu }_{e2})}^T=0\\ f({\mu }_{e2},\delta z{\text{e}}_2^T{(\delta I+{\Theta }_v)}^{-1}{({\mu }_{e1},{\mu }_{e2})}^T)-{\text{e}}_2^T{\Theta }_u{({\mu }_{e1},{\mu }_{e2})}^T=0.\end{array}& (4.12)\end{array}$

The symmetric steady-state solution branch, as characterized by Equation (2.20), is obtained from the root μ_e of the cubic equation

$\begin{array}{ll}{\mu }_e-q{({\mu }_e-2)}^3+4-\frac{\delta z{\mu }_e}{\delta +{\alpha }_v}-{\alpha }_u{\mu }_e=0,& (4.13)\end{array}$

where α_u and α_v are given in Equation (2.19). The symmetry-breaking bifurcation condition (Equation 2.30) along the symmetric steady-state solution branch yields that

$\begin{array}{l}1-3q{({\mu }_e-2)}^2-\frac{\delta z}{\delta +{\alpha }_{v,2}^{\perp }}-{\alpha }_{u,2}^{\perp }=0\\ \iff z({\mu }_e)=\frac{\delta +{\alpha }_{v,2}^{\perp }}{\delta }\left(1-3q{({\mu }_e-2)}^2-{\alpha }_{u,2}^{\perp }\right),\end{array}$

where ${\alpha }_{u,2}^{\perp }$ and ${\alpha }_{v,2}^{\perp }$ are defined in Equation (2.31). We substitute z(μ_e) into Equation (4.13), and solve the resulting equation numerically for μ_e. For ρ = 150, and with the parameters as in the caption of Figure 13, we obtain that there are two supercritical pitchfork bifurcation points z_P,1 and z_P,2 on the symmetric steady-state branch. The linearly stable asymmetric steady-state branches that exist on the range z_P,1 < z < z_P,2 between the two pitchfork points are shown in the Figure 14, Left. When z = z_P,2, we observe from the bifurcation diagram in the Figure 14, Right, together with the eigenvalue computations in Figure 15, that the symmetry-breaking bifurcation is supercritical in terms of ρ.

FIGURE 14

Figure 14. 3-D Bifurcation diagram, computed from Equation (4.12) using MatCont [51], for symmetric and asymmetric steady-states of a two-cell ring pattern with ring radius r = 0.5 and with FN kinetics (Equation 4.11). Left: 3-D Plot of (μ_e1, μ_e2) showing that asymmetric steady-states occur inside the supercritical pitchfork bubble delimited by z_P,1 ≈ 36.75458 and z_P,2 ≈ 41.26889 when ρ = d_v/d_u = 150. Right: Supercritical pitchfork bifurcation from the symmetric branch occurs at ρ_p = 150 when z = z_P,2. Linearly stable asymmetric branches exist past this threshold in ρ. Parameters: D_u = 1, D_v = 4, σ_u = σ_v = 1, ε = 0.03, r = 0.5, q = 1, δ = 0.1, and d_u = 0.04.

FIGURE 15

Figure 15. For the two-cell system with FN kinetics (Equation 4.11) and parameters as in the caption of Figure 14 with z = z_P,2, we plot the numerically computed eigenvalue λ, satisfying σ₂(λ) = 0 in Equation (3.20), vs. ρ. Since λ > 0 only on the range ρ > ρ_p = 150, we conclude that the symmetric steady-state solution is linearly stable to anti-phase eigenperturbations of the form ${\text{q}}_2={(1,-1)}^T$ only when ρ < ρ_p = 150. Moreover, since the root to σ₁(λ) = 0 satisfies λ < 0, we conclude that the symmetric steady-state branch is always linearly stable to in-phase eigenperturbations of the symmetric steady-state.

Next, we illustrate that the bulk-cell model with FN kinetics can also exhibit oscillatory instabilities for in-phase perturbations of the symmetric steady-state. In Figure 16 we plot the bifurcation diagram of μ_e1 vs. z for the same parameter set as in the caption of Figure 14 except that the bulk degradation rates have been decreased slightly to σ_u = σ_v = 0.9. We observe that there are now two Hopf bifurcation values z_H,1 and z_H,2 of z along the symmetric steady-state branch for the in-phase mode that lie within the interval delimited by the two pitchfork bifurcation points. In the Figure 16, Right, we plot the real and imaginary parts of the complex-valued root of σ₁(λ) = 0, as computed from Equation (3.20), which shows that Re(λ) > 0 and Im(λ) ≠ 0 when z_H,1 < z < z_H,2. This leads to the possibility of a synchronous oscillatory instability. As a result, on the range z_H,1 < z < z_H,2, the symmetric steady-state solution branch is unstable to both anti-phase and in-phase perturbations. However, as seen from the Figure 16, Right, where we also plot the growth rate λ for the anti-phase mode as obtained by setting σ₂(λ) = 0 in Equation (3.20), the anti-phase instability has a larger growth rate than the in-phase instability.

FIGURE 16

Figure 16. Left: Bifurcation diagram of μ_e1 vs. z, computed from Equation (4.12) using MatCont [51], for symmetric and asymmetric steady-states of a two-cell ring pattern with the same parameters as in the Figure 14, Left except that now the degradation rates are decreased slightly to σ_u = σ_v = 0.9. For this parameter set, there are Hopf bifurcation points associated with in-phase perturbations of the symmetric steady-state that emerge at z = z_H,1 ≈ 34.65328 and z = z_H,2 ≈ 38.02834 between the two symmetry-breaking pitchfork bifurcation points located at z = z_P,1 ≈ 33.41022 and z = z_P,2 ≈ 38.80742. Right: the root λ of σ₂(λ) = 0 vs. z (pink curve), as computed from Equation (3.20), shows that the symmetric steady-state solution branch is unstable to anti-phase perturbations on the range z_P,1 < z < z_P,2. The plotted real and imaginary parts of the complex-valued root λ_h ≡ λ_r + iλ_i to σ₁(λ_h) = 0, from Equation (3.20), shows that Re(λ_h) > 0 on the range z_H,1 < z < z_H,2. On this range of z, a synchronous oscillatory instability of the symmetric steady-state solution can also occur, but it has a smaller growth rate than that for the anti-phase mode.

4.4. Numerical experiments with closely-spaced cells: GM kinetics

We now briefly explore, from full PDE simulations of Equation (2.3), symmetry-breaking behavior leading to stable asymmetric patterns that can occur for closely spaced cells when the ratio ρ = d_v/d_u is increased. For realistic modeling of pattern formation properties of biological tissues one needs to consider the situation where cells are closely spaced in the sense that the cell radii are either comparable to the distance between the cells, or that there are only narrow gaps between cells. Although the asymptotic theory of Sections 2, 3 is no longer valid for such closely spaced cell arrangements, the FlexPDE simulations of Equation (2.3) shown below reveal a similar qualitative solution behavior as we have analyzed for spatially segregated cells. More specifically, although we no longer have an analytical theory to predict a bifurcation diagram of all steady-state solutions, our full PDE numerical results suggest that stable symmetric steady-states occur only when ρ is below some threshold. When ρ exceeds some symmetry-breaking threshold, stable asymmetric steady-states will be the preferred state. Our numerical results suggest that the critical threshold of ρ that is needed to establish this symmetry-breaking behavior for closely spaced cells is smaller than that needed for spatially segregated cells, if in fact such a threshold exists.

To illustrate this, in Figure 17 we take two closely spaced cells centered near the origin that have a minimum separation of 0.002. The degradation rates, cell radius, and the value of d_u used for Figure 17 are the same as in Table 3, where bifurcation values were given for the two-cell arrangement at different ratios of D_v/D_u with D_u = 5 and for a ring radius r = 0.5. In the cell arrangement in Figure 17, the only difference is that the two cells are now much more closely spaced than in Table 3 and we fix D_v/D_u = 0.3 and D_u = 5. For these parameter values, we observe from Table 3 that no symmetry-breaking bifurcations are possible for this diffusivity ratio when the ring radius is r = 0.5. However, as suggested from the results shown in Figure 17, when the cells are closely spaced there is a symmetry-breaking pitchfork bifurcation point that occurs on the range 3 < ρ < 8. We remark that, rather surprisingly, if we use the symmetry-breaking bifurcation condition (Equation 4.6) from the asymptotic theory for this case of two-closely spaced cells it predicts that ρ_p ≈ 6.53, which lies within the range 3 < ρ < 8. However, we emphasize that the asymptotic theory is not valid for closely-spaced cells.

FIGURE 17

Figure 17. Full numerical simulation results of Equation (2.3) with FlexPDE [50] for GM kinetics (Equation 4.3) with two closely spaced cells centered on a ring of radius r = 0.031 and with minimum cell separation of 0.002. The other parameters are the same as in Table 3. The only difference here is that the cells are now much more closely spaced. The bottom two panels show the concentration of u. Left: convergence to a stable symmetric steady-state solution when ρ = 3. Right: convergence to a stable asymmetric steady-state soluton for ρ = 8 when starting with a symmetric initial condition. Parameters: D_u = 5, D_v = 1.5, σ_u = σ_v = 0.6, d_u = 0.08, and ε = 0.03.

Finally, in Figure 18 we show that the stable asymmetric steady-state patterns can also occur for three closely-spaced cells when ρ exceeds some threshold.

FIGURE 18

Figure 18. Full numerical simulation results of Equation (2.3) with FlexPDE [50] for GM kinetics (Equation 4.3) with three closely spaced cells with cell centers located on the vertices of an equilateral triangle centered at the origin. The ring radius is $r=0.2\sqrt{3}/6$. Left: convergence to a symmetric steady-state when ρ = 5. Right: convergence to an asymmetric steady-state when ρ = 10. Parameters: D_u = D_v = 5, σ_u = σ_v = 0.6, d_u = 0.05, ε = 0.099. The cells have a minimum separation of 0.002. The bottom two panels show the concentration of v. By plotting v rather than u, the bottom right panel clearly shows the asymmetry in the three cells.

5. Discussion

We have analyzed symmetry-breaking behavior associated with the PDE-ODE bulk-cell model (Equation 2.3) where identical two-component intra-compartmental reactions occur only within a disjoint collection of small circular compartments, or “cells,” of a common radius within a bounded 2-D domain. In the bulk, or extra-cellular, medium two bulk species with comparable diffusivities and bulk degradation rates diffuse and globally couple the spatially segregated intracellular reactions. The bulk species are coupled to the intracellular species through an exchange across the compartment boundaries, as modeled by a Robin boundary condition that depends on certain membrane reaction rates. In the limit of a small cell radius, we have used a singular perturbation methodology to derive a nonlinear algebraic system (Equation 2.13) characterizing all the steady-states for the bulk-cell model (Equation 2.3). Moreover, the linear stability properties of the steady-state solutions of the bulk-cell model were shown to be determined by the nonlinear matrix eigenvalue problem (Equation 3.12) of size 2m × 2m, where m is the number of compartments. A root-finding condition on the determinant of this matrix yields the discrete eigenvalues of the linearization (Equation 3.1) around an arbitrary steady-state solution, as defined by the set (Equation 3.13).

We have shown that the steady-state and linear stability theory simplifies considerably for a symmetric cell arrangement, as characterized by Definition (Equation 2.1), and when one of the intracellular species has a linear dependence of the form Equation (2.14). In this more restricted scenario, we have shown that a symmetric steady-state solution, in which the steady-states of the intracellular species are the same for each cell, will exist if the scalar nonlinear algebraic (Equation 2.20) has a solution. We emphasize that since our bulk-cell model does not admit spatially homogeneous steady-state solutions that can be analyzed by a simple Turing-type linear stability approach [1], this symmetric steady-state solution of the bulk-cell model (2.3) represents the base state in our construction. Instabilities and bifurcations associated with this base state are challenging to analyze owing to the fact that the base state is not spatially uniform. Asymmetric steady-state solutions, as determined from Equation (2.16), were shown to bifurcate from the symmetric steady-state solution branch whenever the algebraic criterion (Equation 2.29) is satisfied at some point on the symmetric branch. For a symmetric cell arrangement, the linear stability properties of the symmetric and asymmetric steady-state solution branches are characterized by Equation (3.21) and the roots of the nonlinear matrix eigenvalue problem (Equation 3.18), respectively.

We have implemented our steady-state and linear stability theory for a specific symmetric cell arrangement in which two cells are equally spaced on a ring concentric within a unit disk, and where we have specified either Gierer-Meinhardt, Rauch-Millonas, or FitzHugh-Nagumo intracellular reactions, which all have the simplified form in Equation (2.14). By using parameter continuation numerical software [51] to implement the asymptotic theory, we have shown that the symmetric steady-state solution branch can undergo symmetry-breaking pitchfork bifurcations, leading to linearly stable asymmetric patterns, even when the two bulk diffusing species have identical diffusivities and degradation rates. Overall, we have shown that it is the magnitude of the ratio of the reaction rates for the two bulk species to the cell membranes that determines whether stable asymmetric patterns can occur. This membrane reaction rate ratio threshold condition for the emergence of symmetry-breaking bifurcations is in marked contrast to the well-known large diffusivity ratio threshold condition for pattern formation from a spatially uniform state that is typically derived by a Turing stability analysis for two-component activator-inhibitor RD systems. For FitzHugh-Nagumo and Rauch-Millonas kinetics we have also shown that stable asymmetric patterns can also emerge from a symmetric steady-state pattern at a fixed, but large, membrane reaction rate ratio when a control parameter in the intracellular kinetics is varied. Our theoretical predictions of symmetry-breaking behavior leading to linearly stable asymmetric patterns for a symmetric two-cell arrangement were confirmed through full time-dependent PDE computations of Equation (2.3).

We now briefly relate our theoretical results to some qualitative behavior that has been suggested in chemical and biological applications. Firstly, compartmental-reaction diffusion models of the form Equation (2.3) could potentially be useful for theoretically modeling the collective behavior that occurs for a microemulsion consisting of Belousov-Zhabotinsky (BZ) chemical reactants that are confined within small aqueous droplets that is dispersed in oil [54] (see also Epstein and Xu [55] and Budroni et al. [56]). In this experimental set-up, polar BZ reactants and a catalyst are confined within small immobile droplets, while two non-polar intermediate species generated during the reaction can be transported across the droplet boundaries. These intermediate species diffuse across the domain, with comparable diffusivities, and provide the mechanism for inter-drop coupling [54]. The recent experimental study of Budroni et al. [56] has suggested that it is the relative magnitude of the membrane reaction rates of these intermediates on the droplet boundaries that plays a key role for determining pattern-forming properties for BZ microemulsions. Secondly, with regards to the transport of biological morphogens, it has been suggested in Müller et al. [7] that a differential reaction rate ratio on the cell boundaries for two morphogen species with comparable diffusivities can yield the large effective diffusivity ratio that is needed for pattern formation and symmetry-breaking in tissues. This membrane attachment mechanism, which reduces the effective diffusivity of one of the morphogens and is referred to in Müller et al. [7] as a binding-mediated hindrance diffusion process, may be relevant in many biological applications. Moreover, detailed intracellular mechanisms in biological cells, such as signaling pathways and gene expression rate constants, may also play a pivotal role in large-scale pattern-forming properties of biological tissues [7]. By way of qualitiative comparison, our theoretical analysis of the 2-D bulk-cell model (Equation 2.3) for a very simple 2-cell pattern has revealed that a large membrane reaction rate ratio, together possibly with changes in a parameter in the intracellular kinetics, can trigger the emergence of stable asymmetric steady-state patterns that bifurcate from a symmetry steady-state. However, owing to the complexity of the analysis needed for Equation (2.3), where certain Green's matrices were found to be central to the analysis, it appears rather intractable analytically to isolate via a simple scaling analysis an effective diffusivity for the bulk species that incorporates the membrane reaction rates.

Although we have only applied our theoretical framework to a simple two-cell arrangement, it is rather straightforward to numerically implement the steady-state and linear stability theory for a symmetric cell arrangement with a much larger number of cells. For this scenario, the symmetric steady-state solutions are again determined by the scalar nonlinear algebraic (Equation 2.20) and the linear stability properties of this steady-state are readily studied by computing the union of all the roots of the scalar problems σ_j(λ) = 0, for j ∈ {1, …, m}, in Equation (3.20) that comprise the set (Equation 3.21) that approximates the discrete eigenvalues of the linearization (Equation 3.1) of Equation (2.3) around the steady-state. However, for an arbitrary spatial arrangement of a large number of cells, one key impediment for implementing the linear stability theory for steady-state solutions is with regards to numerically computing the eigenvalues λ from a root-finding condition on the determinant of the full 2m × 2m GCEP matrix ${M}(\lambda )$ in Equation (3.12). This matrix is non-Hermitian, is not sparse, and has an intricate dependence on λ through the Green's matrices. In contrast to the availability of efficient numerical solution strategies for nonlinear matrix eigenvalue problems with special structure, as was discussed in Güttel and Tisseur [47], Betcke et al. [48], and Betcke et al. [49], it appears to a significant open challenge to develop efficient numerical methods to determine all such eigenvalues λ for which ${M}(\lambda )$ is a singular matrix when m is large. Recall that if there are any such eigenvalues in Re(λ) > 0, the steady-state for Equation (2.3) is unstable.

A few other open problems related to our analysis are as follows. Firstly, it would be interesting to analyze symmetry-breaking bifurcation for Equation (2.3) on ℝ² where identical cells of small radii are centered at the lattice points of an arbitrary Bravais lattice. In this periodic setting, it should be possible to analyze symmetry-breaking bifurcations of a periodic steady-state solution by using Floquet-Bloch theory, combined with the explicit analytical formulae for the reduced-wave Bloch Green's function as derived in Iyaniwura et al. [57]. Secondly, it would be interesting to develop an extension of our asymptotic approach to treat closely-spaced cell configurations that are more relevant to modeling pattern-forming properties in biological tissues. Our numerical results shown in Section 4.4 have suggested that only a smaller membrane reaction rate ratio is needed to initiate symmetry-breaking behavior for closely-spaced cells than for arrangements with more spatially segregated cells. To theoretically analyze pattern-forming properties of the bulk-cell model with closely-spaced cells, an extension of the approach developed in Iyaniwura and Ward [58] to analyze the mean first passage time for a cluster of small traps may be fruitful. Thirdly, it would be interesting to formulate and analyze a related bulk-cell model where the chemical reactions occur on the boundaries of a collection of small compartments, rather than in the interior of the compartments. In this scenario, chemical species produced on the membrane can then detach and diffuse in the bulk medium. Such an extension is relevant for analyzing collective behavior that occurs for dynamically reactive solid pellets that are chemically coated and are coupled through a bulk diffusion field (cf. Taylor et al. [59], Taylor et al. [60], Tinsley et al. [61], and Tinsley et al. [62]). Finally, it would be worthwhile to extend our 2-D analysis to a 3-D setting. For a 3-D bounded domain that contains a collection of small spherical compartments, the analysis would be rather different than in 2-D since the free space Green's function has a rapid decay at infinity instead of a logarithmic growth. This would suggest that the inter-cell coupling effect would be, in general, much weaker in 3-D than in 2-D.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

MP and MW analysis, conception, and writing of this article was done 50–50. All authors contributed to the article and approved the submitted version.

Funding

MW was supported by NSERC grant 81541.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2023.1110497/full#supplementary-material

References

1. Turing AM. The chemical basis of morphogenesis. Phil Trans Roy Soc B. (1952) 237:37–72. doi: 10.1098/rstb.1952.0012

CrossRef Full Text

2. Pearson J, Horsthemke W. Turing instabilities with nearly equal diffusivities. J Chem Phys. (1989) 90:1588. doi: 10.1063/1.456051

CrossRef Full Text | Google Scholar

3. Baker RE, Gaffney EA, Maini PK. Partial differential equations for self-organization in cellular and developmental biology. Nonlinearity. (2008) 21:R251. doi: 10.1088/0951-7715/21/11/R05

CrossRef Full Text | Google Scholar

4. Diambra L, Senthivel VR, Menendez DB, Isalan M. Cooperativity to increase Turing pattern space for synthetic biology. AVS Synthetic Biol. (2015) 4:177–86. doi: 10.1021/sb500233u

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Lengyel I, Epstein IR. Modeling of Turing structures in the chlorite-iodide-malonic acid-starch reaction system. Science. (1991) 251:650. doi: 10.1126/science.251.4994.650

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Dulos E, Boissonade J, Perraud JJ, Rudovics B, De Kepper P. Chemical morphogenesis: Turing patterns in an experimental chemical system. Acta Biotheor. (1996) 44:249. doi: 10.1007/BF00046531

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Müller P, Rogers KW, You SR, Brand M, Schier AF. Morphogen transport. Development. (2013) 140:1621–39. doi: 10.1242/dev.083519

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Rauch EM, Millonas MM. The role of trans-membrane signal transduction in Turing-type cellular pattern formation. J Theor Biol. (2004) 226:401–7. doi: 10.1016/j.jtbi.2003.09.018

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Sozen B, Cornwall-Scoones J, Zernicka-Goetz M. The dynamics of morphogenesis in stem cell-based embryology: novel insights for symmetry breaking. Development. (2021) 474:82–90. doi: 10.1016/j.ydbio.2020.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Pearson J. Pattern formation in a (2+1)-species activator-inhibitor immobilizer system. Physica A. (1992) 188:178–89. doi: 10.1016/0378-4371(92)90264-Q

CrossRef Full Text | Google Scholar

11. Klika V, Baker RE, Headon D, Gaffney EA. The influence of receptor-mediated interactions on reaction-diffusion mechanisms of cellular self-organization. Bull Math Bio. (2012) 74:935–57. doi: 10.1007/s11538-011-9699-4

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Korvasová K, A GE, K MP, A FM, Klika V. Investigating the Turing conditions for diffusion-driven instability in the presence of a binding immobile substrate. J Theor Biol. (2015) 367:286–95. doi: 10.1016/j.jtbi.2014.11.024

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Marcon L, Diego X, Sharpe J, Müller P. High throughput mathematical analysis identifies Turing networks for patterning with equal diffusing signals. eLife. (2016) 5:e14022. doi: 10.7554/eLife.14022

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Diego X, Marcon L, Müller P, Sharpe J. Key features of Turing systems are determined purely by network topology. Phys Rev X. (2018) 8:021071. doi: 10.1103/PhysRevX.8.021071

CrossRef Full Text | Google Scholar

15. Landge A, Jordan BM, Diego X, Müller P. Pattern formation mechanisms of self-organizing reaction-diffusion systems. Dev Biol. (2020) 460:2–11. doi: 10.1016/j.ydbio.2019.10.031

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Haas P, Goldstein R. Turing's diffusive threshold in random reaction-diffusion systems. Phys Rev Lett. (2021) 126:238101. doi: 10.1103/PhysRevLett.126.238101

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Vanag VK, Epstein IR. Localized patterns in reaction-diffusion systems. Chaos. (2007) 17:037110. doi: 10.1063/1.2752494

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Ward MJ. Spots, traps, and patches: asymptotic analysis of localized solutions to some linear and nonlinear diffusive processes. Nonlinearity. (2018) 31:R189. doi: 10.1088/1361-6544/aabe4b

CrossRef Full Text | Google Scholar

19. Halatek J, Brauns F, Frey E. Self-organization principles of intracellular pattern formation. Phil Trans R Soc B. (2018) 373:20170107. doi: 10.1098/rstb.2017.0107

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Halatek J, Frey E. Rethinking pattern formation in reaction-diffusion systems. Nat Phys. (2018) 14:507. doi: 10.1038/s41567-017-0040-5

CrossRef Full Text | Google Scholar

21. Maini P, Painter K, Chau HNP. Spatial pattern formation in chemical and biological systems. J Chem Soc Faraday Trans. (1997) 93:3601–3610. doi: 10.1039/a702602a

CrossRef Full Text | Google Scholar

22. Krause A, Gafney EA, Maini P, Klika V. Modern perspectives on near-equilibrium analysis of Turing systems. Phil Trans R Soc A. (2021) 379:20200268. doi: 10.1098/rsta.2020.0268

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Pelz M, Ward MJ. The emergence of spatial patterns for compartmental reaction kinetics coupled by two bulk diffusing species with comparable diffusivities. Phil Trans Roy Soc A. (2022).

Google Scholar

24. Gomez-Marin A, Garcia-Ojalvo J, Sancho JM. Self-sustained spatiotemporal oscillations induced by membrane-bulk coupling. Phys Rev Lett. (2007) 98:168303. doi: 10.1103/PhysRevLett.98.168303

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Levine H, Rappel WJ. Membrane-bound Turing patterns. Phys Rev E. (2005) 72:061912. doi: 10.1103/PhysRevE.72.061912

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Gou J, Li YX, Nagata W, Ward MJ. Synchronized oscillatory dynamics for a 1-D model of membrane kinetics coupled by linear bulk diffusion. SIAM J Appl Dyn Sys. (2015) 14:2096–137. doi: 10.1137/15M1039122

CrossRef Full Text | Google Scholar

27. Gou J, Ward MJ. Oscillatory dynamics for a coupled membrane-bulk diffusion model with Fitzhugh-Nagumo kinetics. SIAM J Appl Math. (2016) 76:776–804. doi: 10.1137/15M1028297

CrossRef Full Text | Google Scholar

28. Gou J, Chiang WY, Lai PY, Ward MJ, Li YX. A theory of synchrony by coupling through a diffusive chemical signal. Physica D. (2017) 339:1–17. doi: 10.1016/j.physd.2016.08.004

CrossRef Full Text | Google Scholar

29. Paquin-Lefebvre F, Nagata W, Ward MJ. Weakly nonlinear theory for oscillatory dynamics in a one-dimensional PDE-ODE model of membrane dynamics coupled by a bulk diffusion field. SIAM J Appl Math. (2020) 80:1520–45. doi: 10.1137/19M1304908

CrossRef Full Text | Google Scholar

30. Xu B, Bressloff P. A PDE-DDE model for cell polarization in fission yeast. SIAM J Appl Math. (2016) 76:1844–70. doi: 10.1137/16M1065458

CrossRef Full Text | Google Scholar

31. Xu B, Jilkine A. Modeling the dynamics of Cdc42 oscillation in fission yeast. Biophysical J. (2018) 114:711–22. doi: 10.1016/j.bpj.2017.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Gou J, Ward MJ. An asymptotic analysis of a 2-D model of dynamically active compartments coupled by bulk diffusion. J Nonlin Sci. (2016) 26:979–1029. doi: 10.1007/s00332-016-9296-7

CrossRef Full Text | Google Scholar

33. Iyaniwura S, Ward MJ. Synchrony and oscillatory dynamics for a 2-D PDE-ODE model of diffusion-mediated communication between small signalling compartments. SIAM J Appl Dyn Sys. (2021) 20:438–99. doi: 10.1137/20M1353666

CrossRef Full Text | Google Scholar

34. Ridgway W, Ward MJ, Wetton BT. Quorum-sensing induced transitions between bistable steady-states for a cell-bulk ODE-PDE model with Lux intracellular kinetics. J Math Bio. (2021) 84:1–2. doi: 10.1007/s00285-021-01705-z

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Gomez D, Iyaniwura S, Paquin-Lefebvre F, Ward MJ. Pattern forming systems coupling linear bulk diffusion to dynamically active membranes or cells. Phil Trans Roy Soc A. (2021) 379:20200276. doi: 10.1098/rsta.2020.0276

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Rätz A. Turing-type instabilities in bulk-surface reaction-diffusion systems. J Comp Appl Math. (2015) 289:142–52. doi: 10.1016/j.cam.2015.02.050

CrossRef Full Text | Google Scholar

37. Elliott CM, Ranner T, Venkataraman C. Coupled bulk-surface free boundary problems arising from a mathematical model of receptor-ligand dynamics. SIAM J Math Anal. (2017) 49:360–397. doi: 10.1137/15M1050811

CrossRef Full Text | Google Scholar

38. Madzvamuse A, Chung AHW, Venkataraman C. Stability analysis and simulations of coupled bulk-surface reaction-diffusion systems. Proc R Soc A. (2015) 471:20140546. doi: 10.1098/rspa.2014.0546

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Madzvamuse A, Chung AHW. The bulk-surface finite element method for reaction-diffusion systems on stationary volumes. Finite Elem Anal Design. (2016) 108:9–21. doi: 10.1016/j.finel.2015.09.002

CrossRef Full Text | Google Scholar

40. Paquin-Lefebvre F, Nagata W, Ward MJ. Pattern formation and oscillatory dynamics in a two-dimensional coupled bulk-surface reaction-diffusion system. SIAM J Appl Dyn Syst. (2019) 18:1334–90. doi: 10.1137/18M1213737

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Cusseddu D, Edelstein-Keshet L, Mackenzie JA, Portet S, Madzvamuse A. A coupled bulk-surface model for cell polarisation. J Theor Biol. (2019) 481:119–35. doi: 10.1016/j.jtbi.2018.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Rätz A, Röger M. Turing instabilities in a mathematical model for signaling networks. J Math Biol. (2012) 65:1215–44. doi: 10.1007/s00285-011-0495-4

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Rätz A, Röger M. Symmetry breaking in a bulk-surface reaction-diffusion model for signalling networks. Nonlinearity. (2014) 27:1805. doi: 10.1088/0951-7715/27/8/1805

CrossRef Full Text | Google Scholar

44. Stolerman LM, Getz M, Llewellyn Smith S, Holst M, Rangamani P. Stability analysis of a bulk-surface reaction model for membrane protein clustering. Bull Math Bio. (2020) 82:2. doi: 10.1007/s11538-020-00703-4

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Paquin-Lefebvre F, Xu B, DiPietro KL, Lindsay AE, Jilkine A. Pattern formation in a coupled membrane-bulk reaction-diffusion model for intracellular polarization and oscillations. J Theor Biol. (2020) 497:110242. doi: 10.1016/j.jtbi.2020.110242

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Gierer A, Meinhardt H. A theory of biological pattern formation. Kybernetik. (1972) 12:30–39. doi: 10.1007/BF00289234

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Güttel S, Tisseur F. The nonlinear eigenvalue problem. Acta Numerica. (2017) 26:1–94. doi: 10.1017/S0962492917000034

CrossRef Full Text | Google Scholar

48. Betcke T, Highan NG, Mehrmann V, Schröder C, F T. NLEVP: a collection of nonlinear eigenvalue problems. ACM Trans Math Software. (2013) 39:7.1–7.28. doi: 10.1145/2427023.2427024

CrossRef Full Text | Google Scholar

49. Betcke T, Highan NG, Mehrmann V, Porzio GMN, Schröder C, F T. NLEVP: a collection of nonlinear eigenvalue problems. In: Users' Guide. MIMS EPring 2011117. Manchester: The University of Manchester, UK (2019).

Google Scholar

50. FlexPDE P. (2015). Solutions Inc. Available online at: http://www.pdesolutions.com.

Google Scholar

51. Dhooge A, Govaerts W, Kuznetsov YA. MatCont: a MATLAB package for numerical bifurcation analysis of ODEs. ACM Trans Math Software. (2003) 29:141–64. doi: 10.1145/779359.779362

CrossRef Full Text | Google Scholar

52. Gierer A. Generation of biological patterns and form: some physical, mathematical, and logical aspects. Progr Biophys Mol Biol. (1981) 37:1–47. doi: 10.1016/0079-6107(82)90019-0

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Bezanson J, Karpinski S, Shah VB, Edelman A. Julia: a fast dynamic language for technical computing. arXiv preprint arXiv:12095145. (2012). doi: 10.48550/arXiv.1209.5145

CrossRef Full Text | Google Scholar

54. Tompkins N, Li N, Girabwe C, Fraden S. Testing Turing's theory of morphogenesis in chemical cells. Pro Natl Acad Sci USA. (2014) 111:4397–402. doi: 10.1073/pnas.1322005111

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Epstein I, Xu B. Reaction-diffusion processes at the nano- and microscales. Nat Technol. (2016) 11:312–9. doi: 10.1038/nnano.2016.41

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Budroni MA, Torbensen K, Ristori S, A AH, Rossi F. Membrane structure drives synchronization patterns in arrays of diffusively coupled self-oscillating droplets. J Phys Chem Lett. (2020) 11:2014–20. doi: 10.1021/acs.jpclett.0c00072

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Iyaniwura S, Gou J, Ward MJ. Synchronous oscillations for a coupled cell-bulk PDE-ODE model with localized cells on ℝ². J Eng Math. (2021) 127:24. doi: 10.1007/s10665-021-10113-7

CrossRef Full Text | Google Scholar

58. Iyaniwura S, Ward MJ. Asymptotic analysis for the mean first passage time in finite or spatially periodic 2-D domains with a cluster of small traps. ANZIAM J. (2021) 63:1–22. doi: 10.21914/anziamj.v63.15976

CrossRef Full Text | Google Scholar

59. Taylor AF, Tinsley M, Showalter K. Insights into collective cell behavior from populations of coupled chemical oscillators. Phys Chemistry Chem Phys. (2015) 17:20047–55. doi: 10.1039/C5CP01964H

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Taylor AF, Tinsley M, Wang F, Huang Z, Showalter K. Dynamical quorum sensing and synchronization in large populations of chemical oscillators. Science. (2009) 323:614–017. doi: 10.1126/science.1166253

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Tinsley MR, Taylor AF, Huang Z, Wang F, Showalter K. Dynamical quorum sensing and synchronization in collections of excitable and oscillatory catalytic particles. Physica D. (2010) 239:785–90. doi: 10.1016/j.physd.2009.08.001

CrossRef Full Text | Google Scholar

62. Tinsley MR, Taylor AF, Huang Z, Showalter K. Emergence of collective behavior in groups of excitable catalyst-loaded particles: spatiotemporal dynamical quorum sensing. Phys Rev Lett. (2009) 102:158301. doi: 10.1103/PhysRevLett.102.158301

PubMed Abstract | CrossRef Full Text | Google Scholar