# Synchrony and Oscillatory Dynamics for a 2-D PDE-ODE Model of Diffusion-
Sensing with Small Signaling Compartments

Sarafa A. Iyaniwura and Michael J. Ward

###### Abstract

We analyze a class of cell-bulk coupled PDE-ODE models, motivated by quorum
and diffusion sensing phenomena in microbial systems, that characterize
communication between localized spatially segregated dynamically active
signaling compartments or “cells” that have a permeable boundary. In this
model, the cells are disks of a common radius $\varepsilon\ll 1$ and they are
spatially coupled through a passive extracellular bulk diffusion field with
diffusivity $D$ in a bounded 2-D domain. Each cell secretes a signaling
chemical into the bulk region at a constant rate and receives a feedback of
the bulk chemical from the entire collection of cells. This global feedback,
which activates signaling pathways within the cells, modifies the
intracellular dynamics according to the external environment. The cell
secretion and global feedback are regulated by permeability parameters across
the cell membrane. For arbitrary reaction-kinetics within each cell, the
method of matched asymptotic expansions is used in the limit $\varepsilon\ll
1$ of small cell radius to construct steady-state solutions of the PDE-ODE
model, and to derive a globally coupled nonlinear matrix eigenvalue problem
(GCEP) that characterizes the linear stability properties of the steady-
states. The analysis and computation of the nullspace of the GCEP as
parameters are varied is central to the linear stability analysis. In the
limit of large bulk diffusivity $D={D_{0}/\nu}\gg 1$, where
$\nu\equiv{-1/\log\varepsilon}$, an asymptotic analysis of the PDE-ODE model
leads to a limiting ODE system for the spatial average of the concentration in
the bulk region that is coupled to the intracellular dynamics within the
cells. Results from the linear stability theory and ODE dynamics are
illustrated for Sel’kov reaction-kinetics, where the kinetic parameters are
chosen so that each cell is quiescent when uncoupled from the bulk medium. For
various specific spatial configurations of cells, the linear stability theory
is used to construct phase diagrams in parameter space characterizing where a
switch-like emergence of intracellular oscillations can occur through a Hopf
bifurcation. The effect of the membrane permeability parameters, the reaction-
kinetic parameters, the bulk diffusivity, and the spatial configuration of
cells on both the emergence and synchronization of the oscillatory
intracellular dynamics, as mediated by the bulk diffusion field, is analyzed
in detail. The linear stability theory is validated from full numerical
simulations of the PDE-ODE system, and from the reduced ODE model when $D$ is
large.

Key Words: Green’s function, bulk diffusion, globally coupled eigenvalue
problem (GCEP), Hopf bifurcation, cell-bulk coupling, synchronous
oscillations, diffusion-sensing.

##  1 Introduction

Cell-to-cell communication is an important aspect of microbial systems that is
often achieved through the diffusion of an extracellular signaling molecule,
referred to as an autoinducer, in their common environment (cf. [13], [45],
[9], [55]). This form of bulk-mediated communication involves the secretion
and feedback of signaling molecules from and into the cells, respectively,
which enables each cell to adjust its intracellular dynamics based on the
signals it receives from the autoinducer field that depends on the entire
collection of cells. Specific autoinducers responsible for such an
intercellular communication have been identified in many biological systems
including, cyclic adenosine monophosphate (cAMP) that triggers intracellular
oscillations for a collection of social amoebae Dictyostelium discoideum and
guides them to aggregation in low nutrient environments (cf. [22], [37],
[16]), acetaldehyde (Ace) that leads to glycolytic oscillations in a colony of
yeast cells (cf. [8], [9], [10]), and acylated homoserine lactones (AHLs) that
induces bioluminescence for certain species of squid due to colonies of the
marine bacterium Vibrio fischeri located in the squid’s light organ (cf.
[45]).

In this context of intercellular bulk-mediated communication, quorum-sensing
(QS) refers to the onset of collective dynamics in the cells that occurs when
the cell density increases past a threshold. There are two main categories of
QS systems, for which mathematical models have been developed. The first main
type, which includes yeast cells and social amoeba, involves a switch-like
onset of synchronized oscillatory intracellular dynamics as the cell density
increases (cf. [8], [9], [10], [22], [33], [27]). On the macroscale, triggered
synchronous temporal oscillations also occur in physicochemical systems
involving small catalyst-loaded particles immersed in a BZ chemical mixture
(cf. [47], [46], [49], [48]). The second main type of QS system, for which the
marine bacterium Vibrio fischeri and the human pathogen Pseudomonas aeruginosa
are prototypical examples, is where an increase in the cell density leads to a
sudden transition between bistable steady-states (cf. [12], [52], [32]).

QS systems that involve a switch-like onset of synchronous intracellular
oscillations have most typically been studied for well-stirred systems, where
the bulk diffusivity is taken to be infinite. This well-mixed limit leads to a
globally coupled ODE system, where the global coupling arises from the
spatially homogeneous bulk diffusion field (cf. [22], [10], [46], [49], [30],
[39], [26]). In this more classical ODE setting, well-established mathematical
tools such as ODE bifurcation theory, phase reduction methods and the Kuramoto
order parameter (cf. [41], [36]) can be used to analyze the onset of QS
behavior and predict the degree of synchronization of intracellular
oscillations as the cell density increases.

However, when the bulk diffusivity is finite, diffusion-sensing behavior
associated with the spatial configuration of cells, the spatial gradients of
the bulk signal, reflecting versus absorbing boundary conditions, and the mass
transport properties of the bulk medium, have all been shown experimentally to
play an important role for some QS systems (cf. [15], [11], [50], [24] and the
references therein). In contrast to the study of QS behavior through an ODE
theoretical framework, there have been relatively few theoretical and modeling
studies of how spatial diffusion triggers the onset of collective
intracellular oscillations in spatially segregated, but localized, dynamically
active reaction sites, and these studies have typically been considered in 1-D
spatial contexts (cf. [17], [7], [19], [21], [18], [31], [38], [57], [40].)

The goal of this paper is to study the emergence and synchronization of
intracellular oscillations for the coupled cell-bulk PDE-ODE model of [20] in
a 2-D domain. The formulation of this model was inspired by the 3-D PDE-ODE
cell-based model of [34] with a single intracellular species (see also [35]
and [51]). In contrast to the analysis in [20] that was restricted either to
the well-mixed limit or to the simple case of a ring pattern of identical
cells, our study will focus on analyzing diffusion-sensing behavior, resulting
from a passive bulk diffusion field with finite diffusivity, for various
spatial configurations of a collection of heterogeneous cells.

Figure 1: A schematic diagram showing dynamically active signaling
compartments (in cyan) in an arbitrary bounded 2D domain. The green and red
dot represent the signaling chemicals in the cells, where only the red is
secreted into the extracellular bulk region. On the right: A zoomed-in
illustration of the intracellular concentration of chemicals within each
signaling compartment, the secretion of signaling molecules into the bulk
region, and the feedback of chemical into the cells.

The coupled PDE-ODE model of [20] is formulated as follows: Within $\Omega$ we
assume that there are $m$ dynamically active circular signalling compartments
or “cells” of a common radius $R_{0}$, denoted by $\Omega_{j}$ and centered at
$\boldsymbol{X}_{j}\in\Omega$ for $j=1,\dots,m$. In the bulk, or
extracellular, region $\Omega\setminus\cup_{j=1}^{m}\,\Omega_{j}$, the
concentration $\mathcal{U}(\boldsymbol{X},T)$ of the autoinducer or bulk
signal, which is confined within $\partial\Omega$, is assumed to satisfy the
passive diffusion equation

$\displaystyle\mathcal{U}_{T}=$ $\displaystyle
D_{B}\,\Delta\,\mathcal{U}-k_{B}\,\mathcal{U}\,,\quad
T>0\,,\quad\boldsymbol{X}\in\Omega\setminus\cup_{j=1}^{m}\,\Omega_{j}\,;$
(1.0.1a)
$\displaystyle\,\partial_{n_{\boldsymbol{X}}}\,\mathcal{U}=\,0,\quad\boldsymbol{X}\in\partial\Omega\,;$
$\displaystyle\qquad
D_{B}\,\partial_{n_{\boldsymbol{X}}}\,\mathcal{U}=\beta_{1j}\,\mathcal{U}-\beta_{2j}\,\mu_{j}^{1}\,,\quad\boldsymbol{X}\in\partial\Omega_{j}\,,\qquad
j=1,\ldots,m\,.$ (1.0.1b) Here $D_{B}>0$ and $k_{B}>0$ are the dimensional
bulk diffusivity and rate of degradation of the bulk signal, respectively. In
the Robin boundary condition (1.0.1b) on the cell membrane, $\beta_{1j}>0$ and
$\beta_{2j}>0$ are dimensional parameters modeling the influx and efflux of
chemical into and out of the $j^{\text{th}}$ cell, while
$\partial_{n_{\boldsymbol{X}}}$ denotes the outer normal derivative on the
cell boundary that points into the bulk region. Within each cell we assume
that there are $n$ interacting species represented by the vector
$\boldsymbol{\mu}_{j}=(\mu_{j}^{1},\ldots,\mu_{j}^{n})^{T}$. Assuming that the
cells are sufficiently small so that there are no spatial chemical gradients
within them, the intracellular reaction-kinetics $\boldsymbol{F}_{j}$ for the
$j^{\text{th}}$ cell is coupled to the bulk medium via the integration of the
flux across the cell membrane. In this way, the intracellular dynamics within
the $j^{\text{th}}$ cell is coupled to the bulk signal (1.0.1) by
$\begin{split}\frac{d\boldsymbol{\mu}_{j}}{dT}&=k_{R}\,\mu_{c}\,\boldsymbol{F}_{j}\left(\boldsymbol{\mu}_{j}/\mu_{c}\right)+\boldsymbol{e}_{1}\int_{\partial\Omega_{j}}\,\left(\beta_{1j}\,\mathcal{U}-\beta_{2j}\,\mu_{j}^{1}\right)\,\text{d}S_{\boldsymbol{X}}\,,\qquad
j=1,\ldots,m\,.\end{split}$ (1.0.1c)

Here $\boldsymbol{e}_{1}=(1,0,\ldots,0)^{T}$, $k_{R}>0$ is the dimensional
reaction rate for the intracellular kinetics, and $\mu_{c}>0$ is a typical
value for $\boldsymbol{\mu}_{j}$. In this model, one signaling chemical,
labeled by $\mu_{j}^{1}$ can permeate the cell membrane with an efflux
parameter $\beta_{2j}$ and, by diffusion through the bulk medium, can
communicate with spatially distant cells. The influx permeability parameter
$\beta_{1j}$ controls the global feedback into the $j^{\mbox{th}}$ cell from
the bulk diffusion field, which is determined by the entire collection of
cells. In Fig. 1 we schematically illustrate the cell-bulk coupling for the
case of $n=2$ intracellular species.

We assume that the radius $R_{0}$ of the signaling compartments is small
relative to the domain length-scale $L$, and so we introduce a small parameter
$\varepsilon\equiv R_{0}/L\ll 1$. Then, by non-dimensionalizing the coupled
PDE-ODE model (1.0.1) using the approach of [20], we obtain that the
dimensionless concentration of chemical $U(\boldsymbol{x},t)$ in the bulk
region satisfies

$\displaystyle\tau\frac{\partial U}{\partial t}=$ $\displaystyle\,D\,\Delta
U-\,U\,,\quad
t>0\,,\quad\boldsymbol{x}\in\Omega\setminus\cup_{j=1}^{m}\,\Omega_{\varepsilon_{j}}\,;$
(1.0.2a)
$\displaystyle\partial_{n}\,U=\,0\,,\quad\boldsymbol{x}\in\partial\Omega\,;$
$\displaystyle\qquad\varepsilon
D\,\partial_{n}U=d_{1j}\,U-d_{2j}\,u_{j}^{1}\,,\quad\boldsymbol{x}\in\partial\Omega_{\varepsilon_{j}}\,,\qquad
j=1,\ldots,m\,,$ (1.0.2b) which is coupled to the dimensionless dynamics
within the $j^{\text{th}}$ cell by
$\displaystyle\frac{\text{d}\boldsymbol{u}_{j}}{\text{d}t}$
$\displaystyle=\,\boldsymbol{F}_{j}\left(\boldsymbol{u}_{j}\right)+\frac{\boldsymbol{e}_{1}}{\varepsilon\tau}\int_{\partial\Omega_{\varepsilon_{j}}}\,(d_{1j}\,U-d_{2j}\,u_{j}^{1})\,\,\text{d}s\,,\qquad
j=1,\ldots,m\,,$ (1.0.2c)

where $\boldsymbol{u}_{j}=(u_{j}^{1},\ldots,u_{j}^{n})^{T}$ is the
dimensionless vector representing the $n$ chemical species in the
$j^{\text{th}}$ cell, labeled by
$\Omega_{\varepsilon_{j}}\equiv\\{{\boldsymbol{x}\,|\,\,|\boldsymbol{x}-\boldsymbol{x}_{j}|\leq\varepsilon\\}}$.
We assume that the centers of the cells are well-separated in the sense that
$\mbox{dist}(\boldsymbol{x}_{j},\boldsymbol{x}_{k})={\mathcal{O}}(1)$ for
$j\neq k$ and
$\mbox{dist}(\boldsymbol{x}_{j},\partial\Omega)={\mathcal{O}}(1)$ as
$\varepsilon\to 0$. In this dimensionless formulation (1.0.2), the key
dimensionless parameters are

$D\equiv\frac{D_{B}}{k_{B}L^{2}}\,,\qquad
d_{1j}\equiv\varepsilon\frac{\beta_{1j}}{k_{B}L}={\mathcal{O}}(1)\,,\qquad
d_{2j}\equiv\varepsilon\frac{\beta_{2j}L}{k_{B}}={\mathcal{O}}(1)\,,\qquad\tau\equiv\frac{k_{R}}{k_{B}}\,.$
(1.0.3)

We refer to $D$ and $\tau$ as the effective bulk diffusivity and reaction-time
parameter, respectively. In (1.0.3), the permeability parameters $\beta_{1j}$
and $\beta_{2j}$ are chosen as $\mathcal{O}(\varepsilon^{-1})$ in order to
ensure that there is an ${\mathcal{O}}(1)$ transport across the membrane of
the small cells. The parameter $\tau$ measures the relative rate of the
intracellular dynamics to the time-scale for degradation of the bulk chemical.
When the intracellular reactions proceed slowly, $\tau$ is small, and little
communication between the cells occurs. When the effective bulk diffusivity
$D$ is large, the cells are readily able to communicate through the bulk
medium, and in the well-mixed limit $D\to\infty$ the bulk signal becomes
spatially homogeneous. Alternatively, for smaller values of $D$, only those
cells that are in close spatial proximity should be able to communicate
through the bulk diffusion field.

For arbitrary intracellular reaction kinetics and for an arbitrary spatial
arrangement of cells, in §2 we use strong localized perturbation theory (cf.
[53], [54]) in the limit $\varepsilon\to 0$ to construct steady-state
solutions of (1.0.2) and to derive the linear stability problem for these
steady-states. Unstable eigenvalues of the linearization of a steady-state are
shown to correspond to roots $\lambda$ in $\mbox{Re}(\lambda)>0$ for which a
certain globally coupled nonlinear matrix eigenvalue problem (GCEP)
${\mathcal{M}}(\lambda)\boldsymbol{c}=0$ (see (2.2.13)) has a nontrivial
solution. The $m\times m$ matrix ${\mathcal{M}}$, which couples all the cells
through an eigenvalue-dependent Green’s matrix, depends on the dimensionless
parameters in (1.0.3). The components of the corresponding normalized
eigenvector $\boldsymbol{c}$ determines the relative magnitude of the spatial
gradient of the bulk signal near the cell membranes and it determines the
relative phases and amplitudes of small-scale oscillations within the cells at
the onset of a Hopf bifurcation when $\lambda=i\lambda_{I}$, with
$\lambda_{I}>0$, is a root of $\det{\mathcal{M}}(\lambda)=0$ (see (2.2.17)).

In §3, strong localized perturbation theory is used to reduce the PDE-ODE
model (1.0.2) to an ODE differential algebraic system with global coupling for
the distinguished limit $D={D_{0}/\nu}\gg 1$ of large bulk diffusivity, where
$\nu\equiv{-1/\log\varepsilon}$. In contrast to the simpler ODE system derived
in [20] for the well-mixed limit $D\to\infty$ (i.e.
$D\gg{\mathcal{O}}(\nu^{-1})$), the new derivation in §3, as summarized in
Proposition 2, leads to an ODE system that depends on the scaled bulk
diffusivity parameter $D_{0}$ and it depends weakly on the specific spatial
configuration of cells within the domain.

Our asymptotic theory is applied for the case of Sel’kov reaction kinetics,
which has been used as a conceptual model for glycolysis oscillations in yeast
cells [44]. With Sel’kov kinetics, (1.0.2) has a unique steady-state when
$\varepsilon\ll 1$. As indicated from the experimental and modeling studies of
collective behavior in yeast cells (cf. [9], [8], [10]), individual yeast
cells are typically non-oscillatory when isolated, but readily become
synchronized in a population of such cells. In qualitative agreement with this
observation, the Sel’kov parameters are chosen to be close to the threshold
for the onset of limit-cycle oscillations for an isolated cell. As a result,
the switch-like emergence of intracellular oscillations, resulting from a Hopf
bifurcation and illustrated through various phase diagrams, is inherently due
to the cell-cell interaction, as mediated by the bulk diffusion field. In our
study we will analyze how the onset of intracellular oscillations and
synchronization depends on the membrane permeability parameters, the reaction-
time parameter, the bulk diffusivity, a Sel’kov kinetic parameter, and the
spatial pattern of cells. Diffusion-sensing behavior, whereby cells initiate
oscillations as a result of certain spatial effects such as cell-clustering or
the buildup of large spatial gradients of the autoinducer field near the cell
boundary, are illustrated.

In §4 we illustrate our theory with Sel’kov reaction kinetics for a ring and
center-cell pattern of cells in the unit disk (see Fig. 3), where the ring
cells are taken to have common parameters but where the parameters for the
center cell can be different. For the unit disk, the Green’s matrices needed
in the steady-state and linear stability theory are available analytically.
For a ring and center-cell pattern, the GCEP matrix ${\mathcal{M}}(\lambda)$
has a cyclic sub-block and so in §4.1 we can analytically identify particular
spatial modes for which $\mbox{det}{\mathcal{M}}(\lambda)=0$. By using
arclength continuation in $D$, Hopf bifurcation boundaries in the $(D,\tau)$
parameter space for which $\mbox{det}{\mathcal{M}}(i\lambda_{I})=0$ can be
computed for each of these modes. In open regions of the $(D,\tau)$ parameter
plane, we show how to use a winding number criterion numerically on the roots
of $\mbox{det}{\mathcal{M}}(\lambda)=0$ in $\mbox{Re}(\lambda)>0$, so as to
compute the number of unstable eigenvalues of the linearization of the unique
steady-state. In §4.2 these phase diagrams are shown for the case of two ring
cells. The triggering effect on the emergence of intracellular oscillations
due to a center cell with either different permeability parameters or a
different Sel’kov kinetic parameter is studied in §4.3 and §4.4, respectively.
Diffusion-sensing behavior as a result of changes in the ring radius are
studied in §4.5. The linear stability theory, which predicts the onset of
intracellular oscillations together with the amplitude and phase differences
near the Hopf bifurcation point, is validated through large-scale simulations
of the PDE-ODE model (1.0.2) using FlexPDE [14]. Moreover, we show that in
certain cases the new ODE system in Proposition 2, derived for
$D={D_{0}/\nu}\gg 1$, can still provide a decent agreement with the full
numerical results even when $D={\mathcal{O}}(1)$.

In §5 we study how the onset of intracellular oscillations depends on the
specific spatial arrangement of ten cells in the unit disk that differ only in
their influx permeabilites $d_{1j}$, for $j=1,\ldots,m$. The cell
configurations considered include two clusters of cells (see Fig. 17(b)), two
rings of cells with two isolated cells (see Fig. 23(b)), and arbitrarily
placed cells (see Fig. 17(d)). For computational simplicity, we primarily
focus on the regime $D={D_{0}/\nu}\gg 1$, where matrix perturbation theory can
be used to asymptotically calculate the spectrum of the GCEP matrix
${\mathcal{M}}(\lambda)$, as summarized in Proposition 4. Since this analysis
shows that only one matrix eigenvalue of ${\mathcal{M}}(\lambda)$ can cross
through zero, in §5.1 a numerical root-finding is readily implemented on this
eigenvalue to determine Hopf bifurcation boundaries in the $(D_{0},\tau)$
plane where intracellular oscillations originate for the various spatial
configurations of cells and permeability parameter sets $d_{1j}$ for
$j=1,\ldots,10$. Our linear stability results, validated through FlexPDE
simulations of (1.0.2) and the ODE dynamics of Proposition 2, shows that when
$D={D_{0}/\nu}$ the intracellular dynamics depend sensitively on the influx
permeability set, but only weakly on the cell locations. In §5.1.1 we
implement the linear stability theory based on the GCEP matrix for
$D={\mathcal{O}}(1)$ to predict that small-scale intracellular oscillations
can be highly heterogeneous in terms of amplitude and phase when there are
isolated cells. The linear stability theory is confirmed from full PDE
simulations. Finally, in §6 we briefly summarize our study and discuss some
biological modeling problems that are well-aligned with the cell-based PDE-ODE
framework of (1.0.1).

##  2 Asymptotic analysis of the dimensionless coupled model

In this section, strong localized perturbation theory in the limit
$\varepsilon\to 0$ is applied to the dimensionless PDE-ODE model (1.0.2) for
the regime $D=\mathcal{O}(1)$. This theory is used to asymptotically
approximate the steady-state solution of the coupled model, and also to
formulate a globally coupled eigenvalue problem (GCEP) for studying the linear
stability properties of the derived steady-state solution.

###  2.1 Asymptotic construction of the steady-state solution

We construct the steady-state solution for (1.0.2) using matched asymptotics.
In the $j^{\text{th}}$ inner region, defined within an
$\mathcal{O}(\varepsilon)$ neighborhood of the boundary of the $j^{\text{th}}$
cell, we introduce the local variables
$\boldsymbol{y}=\varepsilon^{-1}(\boldsymbol{x}-\boldsymbol{x}_{j})$ and
$U_{j}(\boldsymbol{x})=U_{j}(\varepsilon\boldsymbol{y}+\boldsymbol{x_{j}})$,
where $\rho\equiv|\boldsymbol{y}|$. Upon writing (1.0.2) in terms of the inner
variables, for $\varepsilon\to 0$ the steady-state problem in the
$j^{\text{th}}$ inner region is $\Delta U_{j}=0$ for $\rho\geq 1$, subject to
$D\,\partial_{\rho}U_{j}=d_{1j}U_{j}-d_{2j}u_{j}^{1}$ on $\rho=1$. The
radially symmetric solution to this problem is

$\begin{split}U_{j}(\rho)&=A_{j}\log\rho+\frac{1}{d_{1j}}\left(D\,A_{j}+d_{2j}u_{j}^{1}\right)\,,\qquad
j=1,\ldots,m\,,\end{split}$ (2.1.1)

where $A_{j}$ for $j=1,\ldots,m$ are constants to be determined. Upon
substituting (2.1.1) into the steady-state problem of (1.0.2c), we obtain that
the steady-state intracellular dynamics $\boldsymbol{u}_{j}$ of the
$j^{\text{th}}$ cell satisfies

$\begin{split}\boldsymbol{F}_{j}\left(\boldsymbol{u}_{j}\right)+\frac{2\pi
D}{\tau}\,A_{j}\boldsymbol{e}_{1}\,=\boldsymbol{0}\,,\qquad
j=1,\ldots,m\,.\end{split}$ (2.1.2)

This determines $\boldsymbol{u}_{j}$ in terms of the unknown constant $A_{j}$.
To proceed, we must derive another algebraic system for the constants $A_{j}$
for $j=1,\ldots,m$, which is then coupled to (2.1.2).

By matching the far-field behaviour of the inner solution (2.1.1) to an outer
steady-state solution for (1.0.2a), we obtain that the outer solution must
satisfy a specific singularity behaviour as
$\boldsymbol{x}\to\boldsymbol{x}_{j}$. In this way, the steady-state outer
approximation for the bulk solution satisfies

$\begin{split}\Delta
U-\varphi^{2}\,U&=0\,,\quad\boldsymbol{x}\in\Omega\setminus\\{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{m}\\}\,;\qquad\partial_{n}U=0\,,\quad\boldsymbol{x}\in\partial\Omega\,;\\\
U\sim
A_{j}\log|\boldsymbol{x}-\boldsymbol{x}_{j}|&+\frac{A_{j}}{\nu}+\frac{1}{d_{1j}}(DA_{j}+d_{2j}u_{j}^{1})\,,\quad\text{as}\quad\boldsymbol{x}\to\boldsymbol{x}_{j}\,,\qquad
j=1,\ldots,m\,,\end{split}$ (2.1.3)

where $\nu\equiv{-1/\log\varepsilon}$, with $\varepsilon\ll 1$, and
$\varphi\equiv\sqrt{1/D}$. The pre-specification of the regular part of each
singularity structure in (2.1.3) yields a constraint. Overall these
constraints provide an algebraic system for $A_{j}$ for $j=1,\ldots,m$.

To determine this algebraic system, we next represent the solution to (2.1.3)
as

$\begin{split}U=-2\pi\sum_{i=1}^{m}A_{i}G(\boldsymbol{x};\boldsymbol{x}_{i})\,,\end{split}$
(2.1.4)

where the reduced-wave Green’s function $G(\boldsymbol{x};\boldsymbol{x}_{j})$
satisfies

$\begin{split}\Delta
G\,-\,&\varphi^{2}G=-\delta(\boldsymbol{x}-\boldsymbol{x}_{j})\,,\quad\boldsymbol{x}\in\Omega\,;\qquad\partial_{n}G=0\,,\quad\boldsymbol{x}\in\partial\Omega\,;\\\
G(\boldsymbol{x};\boldsymbol{x}_{j})&\sim-\frac{1}{2\pi}\log|\boldsymbol{x}-\boldsymbol{x}_{j}|+R(\boldsymbol{x}_{j})+o(1)\,,\quad\text{as}\quad\boldsymbol{x}\to\boldsymbol{x}_{j}\,.\end{split}$
(2.1.5)

Here $R_{j}\equiv R(\boldsymbol{x}_{j})$ is the regular part of
$G(\boldsymbol{x};\boldsymbol{x}_{j})$ at $\boldsymbol{x}=\boldsymbol{x}_{j}$.
By expanding (2.1.4) as $\boldsymbol{x}\to\boldsymbol{x}_{j}$, we simply
require that the non-singular terms of the resulting expression agree with
that specified in (2.1.3) for each $j=1,\ldots,m$. This leads to an algebraic
system for the vector
$\boldsymbol{\mathcal{A}}\equiv(A_{1},\ldots,A_{m})^{T}$, which is given in
matrix form as

$\begin{split}\Big{(}I+2\pi\nu\mathcal{G}+\nu
DP_{1}\Big{)}\boldsymbol{\mathcal{A}}=-\nu\,P_{2}\,\boldsymbol{u}^{1}\,,\end{split}$
(2.1.6)

where $\boldsymbol{u}^{1}\equiv(u_{1}^{1},u_{2}^{1},\ldots,u_{m}^{1})^{T}$
denotes the vector of chemicals that is secreted into the bulk region by the
cells. In (2.1.6), $\mathcal{G}$ is the symmetric reduced-wave Green’s
interaction matrix, while $P_{1}$ and $P_{2}$ are $m\times m$ diagonal
matrices, defined by

$\begin{split}\mathcal{G}\equiv\begin{pmatrix}R_{1}&G_{12}&\dots&G_{1m}\\\
G_{21}&R_{2}&\dots&G_{2m}\\\ \vdots&\vdots&\ddots&\vdots\\\
G_{m1}&G_{m2}&\dots&R_{m}\end{pmatrix}\,,\quad
P_{1}\equiv\text{diag}\Big{(}\frac{1}{d_{11}},\ldots,\frac{1}{d_{1m}}\Big{)}\,,\quad
P_{2}\equiv\text{diag}\Big{(}\frac{d_{21}}{d_{11}},\ldots,\frac{d_{2m}}{d_{1m}}\Big{)}\,.\end{split}$
(2.1.7)

Here $G_{ji}=G_{ij}\equiv G(\boldsymbol{x}_{j};\boldsymbol{x}_{i})$ for $i\neq
j$, and $R_{j}\equiv R(\boldsymbol{x}_{j})$ for $j=1,\ldots,m$, are obtained
from the solution to (2.1.5).

Overall, the asymptotic steady-state solution is determined in terms of the
solution $\boldsymbol{\mathcal{A}}$ and $\boldsymbol{u}_{j}$ for
$j=1,\ldots,m$, to the $n\times m$ dimensional nonlinear algebraic system
(NAS) given by (2.1.2) and (2.1.6). This system applies to arbitrary local
reaction kinetics $\boldsymbol{F}_{j}$ and permeabilities parameters
$d_{1j}>0$ and $d_{2j}>0$ for $j=1,\ldots,m$. When the kinetics and
permeability parameters are identical for all the cells, the NAS reduces to
the system given in equations $(2.4)$ and $(2.9)$ of [20].

Depending on the specific reaction kinetics assumed, the solution structure to
(2.1.2) and (2.1.6) as parameters are varied can involve solution
multiplicity, saddle-node points, and other bifurcations. However, to
illustrate our asymptotic theory we will focus on the two-component Sel’kov
reaction kinetics for which (2.1.2) and (2.1.6) has a unique solution.

###  2.2 Linear stability analysis

In the previous subsection, we characterized steady-state solutions of the
coupled model (1.0.2) using strong perturbation theory. Suppose that the NAS
(2.1.2) and (2.1.6) has a solution for a given set of parameters. This then
yields an approximation to the steady-state solution solution
$U_{e}(\boldsymbol{x})$ and $\boldsymbol{u}_{ej}$, for $j=1,\ldots,m$, to
(1.0.2). To determine the linear stability of this steady-state we begin by
introducing the perturbation

$\begin{split}U(\boldsymbol{x},t)=U_{e}(\boldsymbol{x})+e^{\lambda
t}\xi(\boldsymbol{x})\qquad\text{and}\qquad\boldsymbol{u}_{j}(t)&=\boldsymbol{u}_{ej}+e^{\lambda
t}\boldsymbol{\phi}_{j}\,,\quad j=1,\ldots,m\,.\end{split}$ (2.2.1)

where $\lambda$ is the eigenvalue of the linearization, and
$\xi(\boldsymbol{x})$ and
$\boldsymbol{\phi}_{j}\equiv(\phi_{j}^{1},\ldots,\phi_{j}^{n})^{T}$ are the
corresponding eigenfunctions in the bulk region and in the $j^{\text{th}}$
cell, respectively. Upon substituting this perturbation into the PDE-ODE model
(1.0.2), in the bulk region we obtain the linearized problem

$\displaystyle\tau\lambda\xi=$
$\displaystyle\,D\,\Delta\xi-\,\xi\,,\quad\boldsymbol{x}\in\Omega\setminus\cup_{j=1}^{m}\,\Omega_{\varepsilon_{j}}\,;$
(2.2.2a)
$\displaystyle\partial_{n}\,\xi=\,0,\quad\boldsymbol{x}\in\partial\Omega\,;$
$\displaystyle\qquad\varepsilon D\,\partial_{nj}\xi=d_{1j}\,\xi-
d_{2j}\,\phi_{j}^{1}\,,\quad\boldsymbol{x}\in\partial\Omega_{\varepsilon_{j}}\,,\qquad
j=1,\ldots,m\,,$ (2.2.2b) which is coupled to the linearized intracellular
dynamics of the $j^{\text{th}}$ cell given in terms of
$\boldsymbol{\phi}_{j}\equiv(\phi_{j}^{1},\ldots,\phi_{j}^{n})^{T}$, by
$\displaystyle\lambda\boldsymbol{\phi}_{j}$
$\displaystyle=J_{j}\,\boldsymbol{\phi}_{j}+\frac{\boldsymbol{e}_{1}}{\varepsilon\tau}\int_{\partial\Omega_{\varepsilon_{j}}}(d_{1j}\,\xi-
d_{2j}\,\phi_{j}^{1})\,\text{d}s\,,\qquad j=1,\ldots,m\,.$ (2.2.2c)

Here $J_{j}\equiv J(\boldsymbol{u}_{ej})$ is the Jacobian matrix of the local
kinetics $\boldsymbol{F}_{j}$ evaluated at the steady-state
$\boldsymbol{u}_{ej}$.

Next, we use strong localized perturbation theory to analyze the eigenvalue
problem (2.2.2) in the limit $\varepsilon\to 0$. This analysis leads to a
limiting globally coupled eigenvalue problem (GCEP) for $\lambda$ in the form
of a nonlinear matrix eigenvalue problem. This GCEP will be used to
investigate instabilities of the steady-state solution for the PDE-ODE system
(1.0.2).

To derive this GCEP, we first construct an inner region in an
$\mathcal{O}(\varepsilon)$ neighborhood of the $j^{\text{th}}$ cell by
introducing the local variables
$\boldsymbol{y}=\varepsilon^{-1}(\boldsymbol{x}-\boldsymbol{x}_{j})$ and
$\xi_{j}(\boldsymbol{x})\equiv\xi_{j}(\boldsymbol{x}_{j}+\varepsilon\boldsymbol{y})$
with $\rho=|\boldsymbol{y}|$. From (2.2.2), we obtain for $\varepsilon\to 0$
that

$\displaystyle\Delta\,\xi_{j}$ $\displaystyle=0,\quad 1<\rho<\infty;\qquad
D\,\partial_{\rho}\,\xi_{j}=d_{1j}\xi_{j}-d_{2j}\phi_{j}^{1},\quad\text{on}\quad\rho=1,$
(2.2.3)

in the $j^{\text{th}}$ inner region. The radially symmetric solution to
(2.2.3) is

$\begin{split}\xi_{j}&=c_{j}\log\rho+\frac{1}{d_{1j}}\left(D\,c_{j}+d_{2j}\phi_{j}^{1}\right)\,,\qquad
j=1,\ldots,m\,,\end{split}$ (2.2.4)

where $\rho=|\boldsymbol{y}|$ and $c_{j}$ for $j=1,\ldots,m$ are constants to
be determined. Upon substituting (2.2.4) into the linearized intracellular
dynamics of the $j^{\text{th}}$ cell (2.2.2c), we obtain a linear relation
between $\boldsymbol{\phi}_{j}$ and $c_{j}$ given by

$\begin{split}(J_{j}-\lambda I)\boldsymbol{\phi}_{j}=-\frac{2\pi
D}{\tau}\,c_{j}\boldsymbol{e}_{1}\,.\qquad j=1,\ldots,m\,.\end{split}$ (2.2.5)

Next, by analyzing the outer solution in the bulk region, we will derive
another linear system, which will be coupled to (2.2.5). These two systems
will provide the GCEP that is needed to study the linear stability of the
steady-state solution.

To determine this additional linear system, we first match the far-field
behaviour of the inner solution (2.2.4) to the outer solution in order to
obtain the singularity behaviour of the outer solution as
$\boldsymbol{x}\to\boldsymbol{x}_{j}$. This yields in the bulk region that

$\begin{split}\Delta\xi-\varphi_{\lambda}^{2}\,\xi&=0\,,\quad\boldsymbol{x}\in\Omega\setminus\\{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{m}\\}\,;\qquad\partial_{n}\xi=0\,,\quad\boldsymbol{x}\in\partial\Omega\,;\\\
\xi\sim
c_{j}\log|\boldsymbol{x}-\boldsymbol{x}_{j}|&+\frac{c_{j}}{\nu}+\frac{1}{d_{1j}}(Dc_{j}+d_{2j}\phi_{j}^{1})\,,\quad\text{as}\quad\boldsymbol{x}\to\boldsymbol{x}_{j}\,,\qquad
j=1,\ldots,m\,,\end{split}$ (2.2.6)

where $\nu\equiv{-1/\log\varepsilon}$ and
$\varphi_{\lambda}=\sqrt{(1+\tau\lambda)/D}$. The solution to (2.2.6) is
represented as

$\begin{split}\xi(\boldsymbol{x})=-2\pi\sum_{i=1}^{m}c_{i}\,G_{\lambda}(\boldsymbol{x};\boldsymbol{x}_{i})\,,\end{split}$
(2.2.7)

where the eigenvalue-dependent reduced-wave Green’s function
$G_{\lambda}(\boldsymbol{x};\boldsymbol{x}_{j})$ satisfies

$\displaystyle\Delta G_{\lambda}-$
$\displaystyle\varphi_{\lambda}^{2}G_{\lambda}=-\delta(\boldsymbol{x}-\boldsymbol{x}_{j})\,,\quad\boldsymbol{x}\in\Omega\,;\qquad\partial_{n}G_{\lambda}=0\,,\quad\boldsymbol{x}\in\partial\Omega\,;$
(2.2.8a) $\displaystyle G_{\lambda}(\boldsymbol{x};\boldsymbol{x}_{j})$
$\displaystyle\sim-\frac{1}{2\pi}\log|\boldsymbol{x}-\boldsymbol{x}_{j}|+R_{\lambda}(\boldsymbol{x}_{j})+o(1)\,,\quad\text{as}\quad\boldsymbol{x}\to\boldsymbol{x}_{j}\,.$
(2.2.8b)

Here $R_{\lambda j}\equiv R_{\lambda}(\boldsymbol{x}_{j})$ is the regular part
of $G_{\lambda}(\boldsymbol{x};\boldsymbol{x}_{j})$ at
$\boldsymbol{x}=\boldsymbol{x}_{j}$. In (2.2.7) we have specified the
principal branch of $\varphi_{\lambda}$ to ensure that $G_{\lambda}$ is
analytic in $\mbox{Re}(\lambda)>0$ and that $G_{\lambda}$ decays far away from
the cells.

By expanding (2.2.7) as $\boldsymbol{x}\to\boldsymbol{x}_{j}$, we equate the
non-singular terms of the resulting expression with those specified in (2.2.6)
for each $j=1,\ldots,m$. This yields a linear system for
$\boldsymbol{c}\equiv(c_{1},\ldots,c_{m})^{T}$ given in matrix form by

$\begin{split}\big{(}I+2\pi\nu\mathcal{G}_{\lambda}+\nu\,DP_{1}\big{)}\,\boldsymbol{c}=-\nu\,P_{2}\,\boldsymbol{\phi}^{1}\,,\end{split}$
(2.2.9)

where $\boldsymbol{\phi}^{1}\equiv(\phi_{1}^{1},\ldots,\phi_{m}^{1})^{T}$ and
$\mathcal{G}_{\lambda}$ is the eigenvalue-dependent reduced-wave Green’s
matrix whose entries are defined by

$(\mathcal{G}_{\lambda})_{ij}=(\mathcal{G}_{\lambda})_{ji}\equiv
G_{\lambda}(\boldsymbol{x}_{j},\boldsymbol{x}_{i})\quad\text{for}\quad i\neq
j\,,\quad\text{and}\quad\ (\mathcal{G}_{\lambda})_{jj}=R_{\lambda j}\equiv
R_{\lambda}(\boldsymbol{x}_{j})\,.$ (2.2.10)

In (2.2.9), the diagonal matrices $P_{1}$ and $P_{2}$ are defined in (2.1.7).

Next, we will combine the algebraic system (2.2.9) with (2.2.5) in order to
derive the GCEP for $\lambda$ and $\boldsymbol{c}=(c_{1},\ldots,c_{m})^{T}$.
We first use (2.2.5) to determine $\phi_{j}^{1}$ in terms of the constant
$c_{j}$ as $\phi_{j}^{1}=2\pi D\tau^{-1}\boldsymbol{e}_{1}^{T}(\lambda
I-J_{j})^{-1}\boldsymbol{e}_{1}\,c_{j}$ for $j=1,\ldots,m$, provided that
$\lambda$ is not an eigenvalue of $J_{j}$ for any $j=1,\ldots,m$. In matrix
form this yields

$\begin{split}\boldsymbol{\phi}^{1}=\frac{2\pi
D}{\tau}\mathcal{K}\boldsymbol{c}\,,\end{split}$ (2.2.11)

where $\mathcal{K}\equiv\mathcal{K}(\lambda)$ is an $m\times m$ diagonal
matrix
$\mathcal{K}\equiv\mbox{diag}\left(\mathit{K}_{1},\ldots,\mathit{K}_{m}\right)$,
whose entries are given by

$\begin{split}\mathit{K}_{j}=\boldsymbol{e}_{1}^{T}(\lambda
I-J_{j})^{-1}\boldsymbol{e}_{1}=\frac{\boldsymbol{e}_{1}^{T}N_{j}\boldsymbol{e}_{1}}{\det(\lambda
I-J_{j})}=\frac{(N_{j})_{11}}{\det(\lambda I-J_{j})}\,.\end{split}$ (2.2.12a)
Here $N_{j}$ is the $n\times n$ matrix of cofactors of $(\lambda I-J_{j})$,
while $(N_{j})_{11}$ is its entry in the first row and first column given by
$(N_{j})_{11}\equiv(N_{j}(\lambda))_{11}=\det\begin{pmatrix}\lambda-\left.\frac{\partial
F_{j}^{2}}{\partial
u_{2}}\right|_{\boldsymbol{u}=\boldsymbol{u}_{e,j}}&\dots&\left.-\frac{\partial
F_{j}^{2}}{\partial u_{n}}\right|_{\boldsymbol{u}=\boldsymbol{u}_{e,j}}\\\
\vdots&\ddots&\vdots\\\ \left.-\frac{\partial F_{j}^{n}}{\partial
u_{2}}\right|_{\boldsymbol{u}=\boldsymbol{u}_{e,j}}&\dots&\lambda-\left.\frac{\partial
F_{j}^{n}}{\partial
u_{n}}\right|_{\boldsymbol{u}=\boldsymbol{u}_{e,j}}\end{pmatrix}\,.$ (2.2.12b)

The functions $F_{j}^{2},\ldots,F_{j}^{n}$ are the components of the local
reaction kinetics $\boldsymbol{F_{j}}\equiv(F_{j}^{1},\ldots,F_{j}^{n})^{T}$
of the $j^{\text{th}}$ cell.

Upon substituting (2.2.11) into (2.2.9), we obtain the $m$-dimensional
homogeneous algebraic system

$\begin{split}\mathcal{M}(\lambda)\boldsymbol{c}&=\boldsymbol{0}\,,\qquad\mbox{where}\qquad\mathcal{M}(\lambda)\equiv
I+2\pi\nu\mathcal{G}_{\lambda}+\nu\,D\,P_{1}+\frac{2\pi\nu
D}{\tau}P_{2}\mathcal{K}\,,\end{split}$ (2.2.13a) where
$\nu\equiv{-1/\log\varepsilon}$. The homogeneous system (2.2.13a), where
${\mathcal{M}}$ is a symmetric but non-Hermitian matrix, is referred to as the
globally coupled eigenvalue problem (GCEP). The GCEP is a nonlinear matrix
eigenvalue problem for $\lambda$, and it has a nontrivial solution
$\boldsymbol{c}\neq\boldsymbol{0}$ if and only $\lambda$ satisfies
$\det\mathcal{M}(\lambda)=0$. We label the set $\Lambda({\mathcal{M}})$ as the
union of all such roots, i.e.
$\begin{split}\Lambda({\mathcal{M}})\equiv\\{{\lambda\,|\,\det\mathcal{M}(\lambda)=0\\}}\,.\end{split}$
(2.2.13b)

The parameters in the GCEP (2.2.13a) are the bulk diffusivity $D$, the
reaction-timescale $\tau$ and the permeabilities $d_{1j}$ and $d_{2j}$, for
$j=1,\ldots,m$, which are encoded in the matrices $P_{1}$ and $P_{2}$ given in
(2.1.7). Moreover, in (2.2.13a), the effect of the spatial configuration of
the centers $\\{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{m}\\}\in\Omega$ of
the cells and the domain shape $\Omega$ arises from both the eigenvalue-
dependent reduced wave Green’s interaction matrix $\mathcal{G}_{\lambda}$ and
the steady-state solution, which determines ${\mathcal{K}}$ in (2.2.12).

A recent survey of nonlinear matrix eigenvalue problems and available solution
strategies for certain classes of matrices is given in [23] and [4]. A range
of applications of such problems, but in simpler contexts where
${\mathcal{M}}(\lambda)$ is either a polynomial or rational function of
$\lambda$, are discussed in [3].

Any element $\lambda\in\Lambda({\mathcal{M}})$ satisfying
$\mbox{Re}(\lambda)>0$ provides an approximation, valid as $\varepsilon\to 0$,
for an unstable discrete eigenvalue of the linearized system (2.2.2). This
leads to the following criterion regarding the linear stability of the steady-
state solution.

######  Proposition 1.

For $\varepsilon\to 0$, a steady-state solution to (1.0.2) is linearly stable
when there are no roots to $\mbox{det}{\mathcal{M}}(\lambda)=0$ in
$\mbox{Re}(\lambda)>0$, i.e. for all $\lambda\in\Lambda({\mathcal{M}})$ we
have $\mbox{Re}(\lambda)<0$. Moreover, if ${\mathcal{A}}_{e}$ and
$\boldsymbol{u}_{ej}$ for $j=1,\ldots,m$ is a non-degenerate solution to the
NAS (2.1.2) and (2.1.6), for which $J_{j}$ is non-singular, then $\lambda=0$
is not a root of $\mbox{det}{\mathcal{M}}(\lambda)=0$.

###### Proof.

For $\varepsilon\to 0$, any discrete eigenvalue $\lambda$ of the linearization
(2.2.1) corresponds to a non-trivial solution to (2.2.2). Since, for
$\varepsilon\to 0$, these discrete eigenvalues comprise the set
$\Lambda({\mathcal{M}})$ in (2.2.13b), the steady-state is linearly stable if
all $\lambda\in\Lambda({\mathcal{M}})$ satisfy $\mbox{Re}(\lambda)<0$. Next,
suppose that ${\mathcal{A}}_{e}$ and $\boldsymbol{u}_{ej}$ for $j=1,\ldots,m$
is a non-degenerate solution to the NAS (2.1.2) and (2.1.6). Introducing the
perturbation ${\mathcal{A}}={\mathcal{A}}_{e}+\boldsymbol{\psi}$ and
$\boldsymbol{u}_{j}=\boldsymbol{u}_{ej}+\boldsymbol{v}_{j}$, where
$|\boldsymbol{v}_{j}|\ll 1$ and $|\boldsymbol{\psi}|\ll 1$, we linearize
(2.1.2) and (2.1.6) to obtain

$\big{(}I+2\pi\nu\mathcal{G}+\nu\,DP_{1}\big{)}\,\boldsymbol{\psi}=-\nu\,P_{2}\,\boldsymbol{v}^{1}\,,\qquad
J_{j}\boldsymbol{v}_{j}=-\frac{2\pi
D}{\tau}\,\psi_{j}\boldsymbol{e}_{1}\,.\qquad j=1,\ldots,m\,,$ (2.2.14)

where $\boldsymbol{v}^{1}\equiv(v_{1}^{1},\ldots,v_{m}^{1})^{T}$,
$\boldsymbol{\psi}=(\psi_{1},\ldots,\psi_{m})^{T}$ and
$\boldsymbol{v}_{j}=(v_{j}^{1},\ldots,v_{j}^{n})^{T}$. Assuming that $J_{j}$
is invertible, (2.2.14) yields that

$\mathcal{J}\boldsymbol{\psi}=\boldsymbol{0}\,,\quad\mbox{where}\quad\mathcal{J}\equiv
I+2\pi\nu\mathcal{G}+\nu\,D\,P_{1}+\frac{2\pi\nu
D}{\tau}P_{2}\mathcal{K}_{0}\quad\mbox{and}\quad{\mathcal{K}}_{0}\equiv-\mbox{diag}\left(\boldsymbol{e}_{1}^{T}J_{1}^{-1}\boldsymbol{e}_{1},\ldots,\boldsymbol{e}_{1}^{T}J_{m}^{-1}\boldsymbol{e}_{1}\right)\,.$
(2.2.15)

We conclude that $\det{\mathcal{J}}\neq 0$ owing to the fact that
${\mathcal{A}}_{e}$ and $\boldsymbol{u}_{ej}$ for $j=1,\ldots,m$ is assumed to
be a non-degenerate solution of the NAS (2.1.2) and (2.1.6). Finally, since it
is readily verified that $\mathcal{J}={\mathcal{M}}(0)$, where
${\mathcal{M}}(\lambda)$ is the GCEP matrix in (2.2.13a), we conclude that
$\lambda=0$ is not a root of (2.2.13b). ∎

Proposition 1 implies that branches of non-degenerate solutions to the NAS
(2.1.2) and (2.1.6) obtained as a parameter is varied cannot lose stability
through a zero-eigenvalue crossing of the GCEP (2.2.13a). This simple
observation is the motivation for analzying whether stability can be lost
through Hopf bifurcations associated with the linearization.

In terms of the eigenvectors $\boldsymbol{c}$ and eigenvalues
$\lambda\in\Lambda({\mathcal{M}})$ of the GCEP (2.2.13), we obtain from
(2.2.1), (2.2.7) and (2.2.5) that the linearization around the bulk and
intracellular steady-state solutions is, for $\varepsilon\to 0$, given by the
superposition

$U(\boldsymbol{x},t)\sim
U_{e}-\sum_{\lambda\in\Lambda({\mathcal{M}})}e^{\lambda
t}\left(\sum_{j=1}^{m}c_{j}G_{\lambda}(\boldsymbol{x};\boldsymbol{x}_{j})\right)\,;\qquad\boldsymbol{u}_{j}(t)\sim\boldsymbol{u}_{ej}+\frac{2\pi
D}{\tau}\sum_{\lambda\in\Lambda({\mathcal{M}})}e^{\lambda t}c_{j}\left(\lambda
I-J_{j}\right)^{-1}\boldsymbol{e}_{1}\,,\quad j=1,\ldots,m\,,$ (2.2.16)

where each $\boldsymbol{c}=(c_{1},\ldots,c_{m})^{T}$ depends on the particular
eigenvalue $\lambda$. To relate the diffusive flux into the $j$-th cell to the
components of $\boldsymbol{c}$ we use the the inner solutions (2.1.1) and
(2.2.4). To determine the effect on the intracellular component that can be
transported across the membrane we calculate
$\boldsymbol{e}_{1}^{T}\boldsymbol{u}_{j}\equiv u_{j}^{1}$ in (2.2.16). This
yields that

$D\partial_{\rho}U|_{\rho=1}\sim
D\left(A_{j}+\sum_{\lambda\in\Lambda({\mathcal{M}})}c_{j}e^{\lambda
t}\right)\,,\qquad u_{j}^{1}\sim u_{ej}^{1}+\frac{2\pi
D}{\tau}\sum_{\lambda\in\Lambda({\mathcal{M}})}\left({\mathcal{K}}\boldsymbol{c}\right)_{j}e^{\lambda
t}\,,\qquad j=1,\ldots,m\,,$ (2.2.17)

where
${\mathcal{K}}\equiv\mbox{diag}\left(\mathit{K}_{1},\ldots,\mathit{K}_{m}\right)$,
with $\mathit{K}_{j}$ defined in (2.2.12). As evident from (2.2.17), and
discussed for various examples in §4.2 and §5.1, the modulus
$|({\mathcal{K}}\boldsymbol{c})_{j}|$ and argument
$\mbox{arg}({\mathcal{K}}\boldsymbol{c})_{j}$ of the components of a complex-
valued ${\mathcal{K}}\boldsymbol{c}$, resulting from pure imaginary
eigenvalues with $\mbox{Re}(\lambda)=0$ and $\mbox{Im}(\lambda)\neq 0$,
determines the relative amplitudes and phase differences of small scale
intracellular oscillations near a Hopf bifurcation of the steady-state
solution.

Our linear stability theory for steady-state solutions of PDE-ODE model
(1.0.2) can be applied to any configuration of cells in an arbitrary 2-D
bounded domain and for arbitrary local reaction kinetics. However, in our
illustrations of the theory below in section 4 and section 5, we will consider
the two-component Sel’kov model, which is used in simple models of glycolysis
(cf. [37], [44]). From (1.0.2c), for an isolated cell with no influx from the
bulk, the intracellular dynamics within the $j^{\text{th}}$ cell that accounts
for efflux from the cell boundary is given by
${\text{d}\boldsymbol{u}_{j}/\text{d}t}=\boldsymbol{F}_{j}\left(\boldsymbol{u}_{j}\right)-{2\pi
d_{2j}u_{j}^{1}\boldsymbol{e}_{1}/\tau}$, where
$\boldsymbol{u}_{j}=(u_{j}^{1},u_{j}^{2})^{T}$ and the Sel’kov kinetics
$\boldsymbol{F}_{j}(v,w)=(F_{j1}(v,w),F_{j2}(v,w))^{T}$ are defined by

$\begin{split}F_{j1}(v,w)=\alpha_{j}w+wv^{2}-v\,,\qquad
F_{j2}(v,w)=\zeta_{j}\left[\mu_{j}-\left(\alpha_{j}w+wv^{2}\right)\right]\,.\end{split}$
(2.2.18)

The steady-state for this isolated cell is given by
$u_{j}^{1}={\mu_{j}/\chi_{j}}$ and
$u_{j}^{2}={\mu_{j}/\left(\alpha_{j}+(u_{j}^{1})^{2}\right)}$, where
$\chi_{j}\equiv 1+{2\pi d_{2j}/\tau}$. The determinant and trace of the
Jacobian $J_{ej}$ for this isolated cell with boundary efflux is

$\mbox{det}(J_{ej})=\zeta_{j}\chi_{j}\left(\alpha_{j}+\frac{\mu_{j}^{2}}{\chi_{j}^{2}}\right)\,,\qquad\mbox{tr}(J_{ej})=\frac{2\mu_{j}^{2}}{\chi_{j}}\left(\alpha_{j}+\frac{\mu_{j}^{2}}{\chi_{j}^{2}}\right)^{-1}-\chi_{j}-\zeta_{j}\left(\alpha_{j}+\frac{\mu_{j}^{2}}{\chi_{j}^{2}}\right)\,.$
(2.2.19)

Since $\mbox{det}(J_{ej})>0$, the steady-state for this isolated cell is
linearly stable only if $\mbox{tr}(J_{ej})<0$. We will choose Sel’kov kinetic
parameters $\alpha_{j}$, $\mu_{j}$ and $\zeta_{j}$ so that an isolated cell
with zero boundary efflux (i.e. $d_{2j}=0$) is linearly stable, but with
parameters rather close to the stability threshold. A set of such parameters
is shown in Fig. 2(a) where we verify that $\mbox{tr}(J_{ej})<0$ on
$0.7<\alpha_{j}<1.0$ when $\mu_{j}=2$ and $\zeta_{j}=0.15$. The Hopf
bifurcation (HB) boundary in the $\alpha_{j}$ versus $\mu_{j}$ parameter
plane, as obtained by setting $\text{tr}(J_{ej})=0$ is

$\alpha_{j}=-\frac{\mu_{j}^{2}}{\chi_{j}^{2}}+\frac{1}{2\zeta_{j}}\left[-\chi_{j}+\sqrt{\chi_{j}^{2}+\frac{8\zeta_{j}\mu_{j}^{2}}{\chi_{j}}}\right]\,,\qquad\mbox{where}\quad\chi_{j}\equiv
1+\frac{2\pi d_{2j}}{\tau}\,.$ (2.2.20)

For an isolated cell, a simple application of the Poincare-Bendixson theorem
shows that when the steady-state is unstable the cell will have limit cycle
oscillations. When there is no boundary efflux, i.e. $d_{2j}=0$, this
parameter range of periodic solutions is given by the green-shaded region in
Fig. 2(b). However, no time-periodic oscillations with Sel’kov kinetics are
possible when the steady-state is linearly stable. In Fig. 2(c) we show how
the HB boundary for an isolated cell depends on the boundary efflux parameter
$d_{2j}$. As expected, for the fixed value $\mu_{j}=2$, the interval in
$\alpha_{j}$ where oscillations are possible is decreased when there is an
efflux out of the cell boundary. The shifting of the HB boundaries to the
right in Fig. 2(c) indicates that a greater rate $\mu_{j}$ of production of
$u_{j}^{2}$ is needed to ensure oscillations when there is a boundary efflux.
The interval in $\mu_{j}$ where oscillations are possible, at least for some
range of $\alpha_{j}>0$, is $0<\mu_{j}<{\chi_{j}^{3/2}/\sqrt{\zeta_{j}}}$.

For the baseline parameter set $\mu_{j}=2$ and $\zeta_{j}=0.15$, in section 4
and section 5 we will show that the inter-cell coupling via the bulk diffusion
field can be sufficient to trigger an oscillatory instability in the cells
through a Hopf bifurcation.

(a) $\mbox{tr}(J_{ej})$ versus $\alpha_{j}$ for isolated cell

(b) Instability region: isolated cell

(c) HB boundaries: isolated cell with efflux

Figure 2: Left panel: $\mbox{tr}(J_{ej})$, from (2.2.19) versus $\alpha_{j}$,
for the steady-state of the Sel’kov kinetics (2.2.18) for an isolated cell,
with $\mu_{j}=2$ and $\zeta_{j}=0.15$. This steady-state is linearly stable
but the parameters are close to the stability threshold. Middle panel: Green-
shaded region of instability where $\mbox{tr}(J_{ej})>0$ in the $\alpha_{j}$
versus $\mu_{j}$ plane for the steady-state of an isolated cell with no
boundary efflux and $\zeta_{j}=0.15$. Within this region, a time-periodic
solution (limit cycle) occurs for an isolated cell. The HB boundary is given
by (2.2.20) with $d_{2j}=0$. In the unshaded region the steady-state is
linearly stable. Right panel: HB boundaries for an isolated cell (see
(2.2.20)) with $\zeta_{j}=0.15$, $\tau=0.5$ and four boundary efflux
parameters. A larger production rate $\mu_{j}$ is needed to support
oscillations.

For Sel’kov local reaction kinetics (2.2.18) it is readily shown that, with an
arbitrary arrangement of cells, there is a unique solution to the NAS (2.1.2)
and (2.1.6) given by

$u_{ej}^{1}=\mu_{j}+\frac{2\pi D}{\tau}A_{j}\,,\qquad
u_{ej}^{2}=\frac{\mu_{j}}{\alpha_{j}+\left(u_{ej}^{1}\right)^{2}}\,,\qquad
j=1,\ldots,m\,,$ (2.2.21a) where ${\mathcal{A}}=(A_{1},\ldots,A_{m})^{T}$
satisfies the linear algebraic system
$\left(I+2\pi\nu\mathcal{G}+\nu\,D\,P_{1}+\frac{2\pi\nu
D}{\tau}P_{2}\right){\mathcal{A}}=-\nu
P_{2}\,\boldsymbol{\mu}\,,\qquad\mbox{with}\qquad\boldsymbol{\mu}\equiv(\mu_{1},\ldots,\mu_{m})^{T}\,.$
(2.2.21b)

For $\nu\ll 1$ sufficiently small, the matrix in (2.2.21b) is invertible,
yielding a unique solution for ${\mathcal{A}}$. For Sel’kov kinetics we
conclude that steady-state solutions of the PDE-ODE model (1.0.2) are always
non-degenerate as $\varepsilon\to 0$. As such, since by Proposition 1
stability cannot be lost via a zero-eigenvalue crossing, in section 4 and
section 5 we will focus on analyzing instabilities of the steady-state arising
from Hopf bifurcations.

##  3 A limiting ODE system with global coupling: $D=\mathcal{O}(\nu^{-1})$

In this section, we use a singular perturbation approach to reduce the
dimensionless coupled PDE-ODE model (1.0.2) to an ODE system that is valid for
the limiting regime $D=\mathcal{O}(\nu^{-1})\gg\mathcal{O}(1)$, where
$\nu\equiv{-1/\log\varepsilon}$ and $\varepsilon\ll 1$. This ODE system
depends weakly on the spatial configuration of the cells and on the scaled
diffusivity $D_{0}={\mathcal{O}}(1)$, defined by $D={D_{0}/\nu}$.

Consider a collection of $m$ small cells centered at the points
$\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{m}$ in a 2-D bounded domain
$\Omega$. Define $\Omega_{p}\equiv\cup_{j=1}^{m}\Omega_{\varepsilon_{j}}$ as
the region formed by the union of the cells, and the average bulk
concentration $\overline{U}=\overline{U}(t)$ by

$\overline{U}=\frac{1}{|\Omega\setminus\Omega_{p}|}\int_{\Omega\setminus\Omega_{p}}U\,\text{d}\boldsymbol{x}\,.$
(3.0.1)

Our goal is to derive an ODE for $\overline{U}\equiv\overline{U}(t;\nu)$,
accurate to ${\mathcal{O}}(\nu)$, which is coupled to the intracellular
dynamics of the cells, as given in (1.0.2b). Upon multiplying (1.0.2a) by
$1/|\Omega\setminus\Omega_{p}|$ and using the divergence theorem, we obtain

$\tau\overline{U}_{t}+\overline{U}=\frac{2\pi}{|\Omega\setminus\Omega_{p}|}\sum_{j=1}^{m}\left(d_{2j}u_{j}^{1}-\frac{d_{1j}}{2\pi\varepsilon}\int_{\partial\Omega_{\varepsilon_{j}}}U\,\text{d}s\right)\,,$
(3.0.2)

where we used $|\partial\Omega_{\varepsilon_{j}}|=2\pi\varepsilon$ for the
perimeter of each cell. Since $|\Omega_{p}|=m\pi\varepsilon^{2}$, we estimate
$|\Omega\setminus\Omega_{p}|=|\Omega|-\mathcal{O}(\varepsilon^{2})$, so that
$|\Omega\setminus\Omega_{p}|\to|\Omega|$ as $\varepsilon\to 0$ in (3.0.2).
Next, by evaluating the integral in (1.0.2b), and using
$|\partial\Omega_{\varepsilon_{j}}|=2\pi\varepsilon$, we obtain

$\frac{\text{d}\boldsymbol{u}_{j}}{\text{d}t}=\boldsymbol{F}_{j}(\boldsymbol{u}_{j})-\frac{2\pi\boldsymbol{e}_{1}}{\,\tau}\left(d_{2j}u_{j}^{1}-\frac{d_{1j}}{2\pi\varepsilon}\int_{\partial\Omega_{\varepsilon_{j}}}U\,\text{d}s\right)\,,\qquad
j=1,\ldots,m\,.$ (3.0.3)

In (3.0.2) and (3.0.3), we must estimate the bulk concentration $U\equiv
U(\boldsymbol{x},t)$ on the $j^{\text{th}}$ cell boundary
$\partial\Omega_{\varepsilon_{j}}$.

For $D={\mathcal{O}}(\nu^{-1})$ we introduce the scaling

$D=\frac{D_{0}}{\nu},\quad\text{where}\quad\nu\equiv\frac{-1}{\log\varepsilon}\,,\quad\varepsilon\ll
1\quad\text{and}\quad D_{0}=\mathcal{O}(1)\,.$ (3.0.4)

Upon substituting (3.0.4) into (1.0.2), we obtain in the bulk region that

$\begin{split}\tau U_{t}&=\frac{D_{0}}{\nu}\,\Delta U-\,U\,,\qquad
t>0\,,\quad\boldsymbol{x}\in\Omega\setminus\cup_{j=1}^{m}\Omega_{\varepsilon_{j}}\,;\\\
\partial_{n}\,U=\,0,\quad\boldsymbol{x}\in\partial\Omega_{\varepsilon}\,;&\qquad\varepsilon\frac{D_{0}}{\nu}\,\partial_{n}U=d_{1j}\,U-d_{2j}u_{j}^{1}\quad\text{on}\quad\partial\Omega_{\varepsilon_{j}}\,,\qquad
j=1,\ldots,m\,.\end{split}$ (3.0.5)

In the inner region at an $\mathcal{O}(\varepsilon)$ neighborhood of the
$j^{\text{th}}$ cell, we introduce the inner variables
$\boldsymbol{y}=\varepsilon^{-1}(\boldsymbol{x}-\boldsymbol{x}_{j})$ and
$U(\boldsymbol{x})=V_{j}(\boldsymbol{x}_{j}+\varepsilon\boldsymbol{y};\nu)$,
with $\rho=|\boldsymbol{y}|$. Writing (3.0.5) in terms of these inner
variables, we obtain for $\varepsilon\to 0$ that

$\begin{split}\Delta V_{j}=0,\quad\rho\geq
1\,;\qquad\,\partial_{\rho}\,V_{j}=\frac{\nu}{D_{0}}\Big{(}\,d_{1j}\,V_{j}-d_{2j}\,u_{j}^{1}\Big{)}\,,\quad\text{on}\quad\rho=1\,,\qquad
j=1,\ldots,m\,,\end{split}$ (3.0.6)

which has the radially symmetric solution

$\begin{split}V_{j}=\nu\,b_{j}\log\rho+U^{0}_{j}\,,\qquad\text{where}\qquad
b_{j}\equiv\frac{1}{D_{0}}\Big{(}d_{1j}\,U^{0}_{j}-d_{2j}\,u_{j}^{1}\Big{)}\,,\qquad
j=1,\ldots,m\,,\end{split}$ (3.0.7)

where $V_{j}|_{\rho=1}=U_{j}^{0}$ is to be determined. By writing (3.0.7) in
the outer $\boldsymbol{x}$ variable, and using
$|\boldsymbol{y}|=\varepsilon^{-1}|\boldsymbol{x}-\boldsymbol{x}_{j}|$, we
obtain the following asymptotic matching condition for the outer solution in
the bulk region:

$\begin{split}U\sim\nu\,b_{j}\log|\boldsymbol{x}-\boldsymbol{x}_{j}|+\Big{(}\frac{d_{1j}}{D_{0}}+1\Big{)}\,U^{0}_{j}-\frac{d_{2j}}{D_{0}}\,u_{j}^{1}\,,\qquad\text{as}\qquad\boldsymbol{x}\to\boldsymbol{x}_{j}\,.\end{split}$
(3.0.8)

Since $V_{j}(\rho)=U_{j}^{0}$ on $\rho=1$ from (3.0.7), we have
$\int_{\partial\Omega_{\varepsilon
j}}U\,\text{d}s=\varepsilon\int_{0}^{2\pi}V_{j}|_{\rho=1}\,\text{d}\theta=2\pi\varepsilon\,U_{j}^{0}$.
Then, from (3.0.2) and (3.0.3), and recalling (3.0.7) for $b_{j}$, we obtain
for $\varepsilon\to 0$ that

$\begin{split}\tau\overline{U}_{t}+\overline{U}&=\frac{2\pi}{|\Omega|}\sum_{j=1}^{m}\left(d_{2j}u_{j}^{1}-d_{1j}U_{j}^{0}\right)=-\frac{2\pi
D_{0}}{|\Omega|}\sum_{j=1}^{m}b_{j}\,,\\\
\frac{\text{d}\boldsymbol{u}_{j}}{\text{d}t}-\boldsymbol{F}_{j}(\boldsymbol{u}_{j})&=-\frac{2\pi\boldsymbol{e}_{1}}{\,\tau}\left(d_{2j}u_{j}^{1}-d_{1j}U_{j}^{0}\right)=\frac{2\pi
D_{0}}{\tau}\boldsymbol{e}_{1}b_{j}\,,\qquad j=1,\ldots,m\,.\end{split}$
(3.0.9)

To complete the derivation of the ODE system we must obtain an algebraic
system for $b_{j}$ for $j=1,\ldots,m$ from the analysis of the outer solution.
From (3.0.5) and (3.0.8), and relating $U_{j}^{0}$ to $b_{j}$ using (3.0.7),
the outer problem for $U(\boldsymbol{x},t)$ is

$\begin{split}\tau U_{t}&=\frac{D_{0}}{\nu}\,\Delta U-\,U\,,\quad
t>0\,,\quad\boldsymbol{x}\in\Omega\setminus\cup_{j=1}^{m}\Omega_{\varepsilon_{j}}\,;\qquad\partial_{n}\,U=\,0\,,\quad\boldsymbol{x}\in\Omega\,;\\\
\qquad U&\sim\nu
b_{j}\log|\boldsymbol{x}-\boldsymbol{x}_{j}|+b_{j}\left(1+\frac{D_{0}}{d_{1j}}\right)+\frac{d_{2j}}{d_{1j}}u_{j}^{1}\,,\qquad\text{as}\qquad\boldsymbol{x}\to\boldsymbol{x}_{j}\,.\end{split}$
(3.0.10)

We then expand $U(\boldsymbol{x},t)$ as

$U(\boldsymbol{x},t)=\overline{U}+\frac{\nu}{D_{0}}\,U_{1}(\boldsymbol{x},t)+\ldots\,,\qquad\mbox{where}\quad\int_{\Omega}U_{1}\,d\boldsymbol{x}=0\,.$
(3.0.11)

The zero average constraint on $U_{1}$ ensures that $\overline{U}$ is the
spatial average of $U$ to terms of order ${\mathcal{O}}(\nu)$. Upon
substituting (3.0.11) into (3.0.10), we obtain in the sense of distributions
that $U_{1}$ satisfies

$\begin{split}\Delta U_{1}&=\tau\,\overline{U}_{t}+\overline{U}+2\pi
D_{0}\sum_{i=1}^{m}b_{i}\delta(\boldsymbol{x}-\boldsymbol{x}_{i})\,,\quad\boldsymbol{x}\in\Omega\,;\qquad\partial_{n}U_{1}=0\,,\quad\boldsymbol{x}\in\partial\Omega\,;\\\
U_{1}&\sim
b_{j}D_{0}\log|\boldsymbol{x}-\boldsymbol{x}_{j}|-\frac{D_{0}}{\nu}\overline{U}+\frac{D_{0}}{\nu}\left[b_{j}\Big{(}1+\frac{D_{0}}{d_{1j}}\Big{)}+\frac{d_{2j}}{d_{1j}}\,u_{j}^{1}\right]\,,\qquad\text{as}\qquad\boldsymbol{x}\to\boldsymbol{x}_{j}\,.\end{split}$
(3.0.12)

The divergence theorem applied to (3.0.12) yields the ODE given in (3.0.9)
while the linear system for $b_{j}$, for $j=1,\ldots,m$, is obtained from the
constraints involved with specifying the form of the regular part of the
singularity behavior in (3.0.12).

The solution to (3.0.12), with $\int_{\Omega}U_{1}\,d\boldsymbol{x}=0$ is
written as

$U_{1}=-2\pi
D_{0}\sum_{i=1}^{m}b_{i}\,G_{0}(\boldsymbol{x};\boldsymbol{x}_{i})\,,$
(3.0.13)

where $G_{0}(\boldsymbol{x};\boldsymbol{x}_{j})$ is the unique Neumann Green’s
function satisfying

$\displaystyle\Delta G_{0}$
$\displaystyle=\frac{1}{|\Omega|}-\delta(\boldsymbol{x}-\boldsymbol{x}_{j})\,,\quad\boldsymbol{x}\in\Omega\,;\qquad\partial_{n}G_{0}=0\,,\quad\boldsymbol{x}\in\partial\Omega\,;$
(3.0.14a) $\displaystyle G_{0}(\boldsymbol{x};\boldsymbol{x}_{j})$
$\displaystyle\sim-\frac{1}{2\pi}\log|\boldsymbol{x}-\boldsymbol{x}_{j}|+R_{0j}+o(1)\,,\quad\text{as}\quad\boldsymbol{x}\to\boldsymbol{x}_{j}\,,\quad\text{and}\quad\int_{\Omega}\,G_{0}\,\text{d}\boldsymbol{x}=0\,.$
(3.0.14b)

Here $R_{0j}$ is the regular part of $G_{0}$ at
$\boldsymbol{x}=\boldsymbol{x}_{j}$. By expanding (3.0.13) as
$\boldsymbol{x}\to\boldsymbol{x}_{j}$, we enforce that the nonsingular part of
the resulting expression agrees with that in (3.0.12). This yields that

$b_{j}\left(1+\frac{D_{0}}{d_{1j}}\right)+2\pi\nu\left(b_{j}R_{0j}+\sum_{i\neq
j}^{m}b_{i}G_{0ji}\right)=\overline{U}-\frac{d_{2j}}{d_{1j}}u_{j}^{1}\,,\qquad
j=1,\ldots,m\,,$ (3.0.15)

where $G_{0ji}=G_{0}(\boldsymbol{x}_{j};\boldsymbol{x}_{i})$. This linear
system for $b_{1},\ldots,b_{m}$ is then coupled to the ODEs given in (3.0.9).
Upon writing this ODE system in matrix form we summarize the result as
follows:

######  Proposition 2.

Let $\varepsilon\to 0$ and assume that $D={D_{0}/\nu}\gg 1$ where
$D_{0}={\mathcal{O}}(1)$ and $\nu={-1/\log\varepsilon}\ll 1$. Then, the PDE-
ODE system (1.0.2) reduces to the following $nm+1$ dimensional ODE system for
$\overline{U}\approx|\Omega|^{-1}\int_{\Omega}U\,d\boldsymbol{x}$ and the
intracellular species:

$\begin{split}\frac{\text{d}}{\text{d}t}\overline{U}=-\frac{1}{\tau}\overline{U}-\frac{2\pi
D_{0}}{\tau|\Omega|}\boldsymbol{e}^{T}\boldsymbol{b}\,;\qquad\frac{\text{d}\boldsymbol{u}_{j}}{\text{d}t}=\boldsymbol{F}_{j}(\boldsymbol{u}_{j})+\frac{2\pi
D_{0}\boldsymbol{e}_{1}}{\,\tau}b_{j}\,,\qquad j=1,\ldots,m\,,\end{split}$
(3.0.16a) where $\boldsymbol{e}\equiv(1,\ldots,1)^{T}$,
$\boldsymbol{e}_{1}\equiv(1,0,\ldots,0)^{T}$ and
$\boldsymbol{b}\equiv(b_{1},\ldots,b_{m})^{T}$. In (3.0.16a), $\boldsymbol{b}$
is the solution to the linear system
$\big{(}I+D_{0}P_{1}+2\pi\nu\,\mathcal{G}_{0}\big{)}\boldsymbol{b}=\overline{U}\,\boldsymbol{e}-P_{2}\,\boldsymbol{u}^{1}\,,$
(3.0.16b)

where $\boldsymbol{u}^{1}\equiv(u^{1}_{1},\ldots,u^{1}_{m})^{T}$ and $P_{1}$
and $P_{2}$ are the diagonal matrices defined in terms of the permeabilities
by (2.1.7). In (3.0.16b), ${\mathcal{G}}_{0}$ is the Neumann Green’s matrix
with matrix entries

$(\mathcal{G}_{0})_{ij}=(\mathcal{G}_{0})_{ji}=G_{0}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})\,,\quad
i\neq j\qquad\text{and}\quad(G_{0})_{jj}=R_{0j}\,.$ (3.0.17)

For $\nu\ll 1$, this ODE system is equivalent up to ${\mathcal{O}}(\nu)$ terms
to

$\begin{split}\frac{\text{d}}{\text{d}t}\overline{U}&=-\frac{1}{\tau}\overline{U}-\frac{2\pi}{\tau|\Omega|}\left[\overline{U}\boldsymbol{e}^{T}{\mathcal{C}}\boldsymbol{e}-\boldsymbol{e}^{T}{\mathcal{C}}P_{2}\,\boldsymbol{u}^{1}\right]+{\mathcal{O}}(\nu^{2})\,,\\\
\frac{\text{d}\boldsymbol{u}_{j}}{\text{d}t}&=\boldsymbol{F}_{j}(\boldsymbol{u}_{j})+\frac{2\pi\boldsymbol{e}_{1}}{\,\tau}\left[\overline{U}\left({\mathcal{C}}\boldsymbol{e}\right)_{j}-\left({\mathcal{C}}P_{2}\,\boldsymbol{u}^{1}\right)_{j}\right]\,,\qquad
j=1,\ldots,m\,,\end{split}$ (3.0.18a) where the matrix ${\mathcal{C}}$ is
defined in terms of ${\mathcal{G}}_{0}$ and a diagonal matrix ${\mathcal{P}}$
by
${\mathcal{C}}\equiv{\mathcal{P}}-\frac{2\pi\nu}{D_{0}}{\mathcal{P}}{\mathcal{G}}_{0}{\mathcal{P}}\,,\qquad{\mathcal{P}}\equiv\mbox{diag}\left(\frac{D_{0}d_{11}}{d_{11}+D_{0}}\,,\ldots,\frac{D_{0}d_{1m}}{d_{1m}+D_{0}}\right)\,.$
(3.0.18b)

In the well-mixed limit for which $D_{0}\to\infty$, (3.0.16b) yields
$D_{0}b_{j}\sim\overline{U}d_{1j}-d_{2j}u_{j}^{1}$, so that (3.0.16a) reduces
to

$\overline{U}_{t}=-\frac{1}{\tau}\overline{U}-\frac{2\pi}{\tau|\Omega|}\sum_{j=1}^{m}(\overline{U}d_{1j}-d_{2j}u_{j}^{1})\,;\qquad\frac{\text{d}\boldsymbol{u}_{j}}{\text{d}t}=\boldsymbol{F}_{j}(\boldsymbol{u}_{j})+\frac{2\pi\boldsymbol{e}_{1}}{\,\tau}(\overline{U}d_{1j}-d_{2j}u_{j}^{1})\,,\qquad
j=1,\ldots,m\,.$ (3.0.19)

To derive (3.0.18) from (3.0.16), we approximate the solution $\boldsymbol{b}$
to (3.0.16b) up to terms of order ${\mathcal{O}}(\nu)$. By inverting the
diagonal matrix $I+D_{0}P_{1}$, we obtain from (3.0.16b) that

$\boldsymbol{b}=\frac{1}{D_{0}}\left(I+\frac{2\pi\nu}{D_{0}}{\mathcal{P}}{\mathcal{G}}_{0}\right)^{-1}\left(\overline{U}{\mathcal{P}}\boldsymbol{e}-{\mathcal{P}}P_{2}\boldsymbol{u}^{1}\right)\sim\frac{1}{D_{0}}\left(I-\frac{2\pi\nu}{D_{0}}{\mathcal{P}}{\mathcal{G}}_{0}\right)\left(\overline{U}{\mathcal{P}}\boldsymbol{e}-{\mathcal{P}}P_{2}\boldsymbol{u}^{1}\right)=\frac{1}{D_{0}}\left(\overline{U}{\mathcal{C}}\boldsymbol{e}-{\mathcal{C}}P_{2}\boldsymbol{u}^{1}\right)+{\mathcal{O}}(\nu^{2})\,,$
(3.0.20)

where ${\mathcal{C}}$ and ${\mathcal{P}}$ are given in (3.0.18b). The ODE
system (3.0.18) results from substituting (3.0.20) into (3.0.16).

The ODE systems (3.0.16), or alternatively (3.0.18), for the regime
$D={\mathcal{O}}(\nu^{-1})$ are accurate up to and including terms of order
${\mathcal{O}}(\nu)$ and show how the intracellular species are globally
coupled through the spatial average of the bulk field. Since these ODE systems
depend on the scaled diffusivity parameter $D_{0}$ and include the effect of
the spatial configuration $\boldsymbol{x}_{1}\,,\ldots,\boldsymbol{x}_{m}$ of
the cells through the Neumann Green’s matrix, these ODE systems can account
for both diffusion-sensing and quorum-sensing behavior (see section 4 and
section 5). In contrast, the limiting well-mixed ODE system (3.0.19),
originally derived in [20] for the simpler case of identical cells, depends
only on the number $m$ of cells. As a result, the well-mixed ODE dynamics is
independent of the diffusivity and the spatial configuration of the cells.

In our numerical experiments in section 4 and section 5 using the ODE system
(3.0.16) the domain $\Omega$ is the unit disk. For the unit disk, the Neumann
Green’s function $G_{0}(\boldsymbol{x};\boldsymbol{x}_{j})$ and its regular
part $R_{0j}$, satisfying (3.0.14), are (see equation (4.3) of [28])

$\begin{split}G_{0}(\boldsymbol{x};\boldsymbol{x}_{j})&=-\frac{1}{2\pi}\log|\boldsymbol{x}-\boldsymbol{x}_{j}|-\frac{1}{4\pi}\log\left(|\boldsymbol{x}|^{2}|\boldsymbol{x}_{j}|^{2}+1-2\boldsymbol{x}\cdot\boldsymbol{x}_{j}\right)+\frac{(|\boldsymbol{x}|^{2}+|\boldsymbol{x}_{j}|^{2})}{4\pi}-\frac{3}{8\pi},\\\
R_{0j}&=-\frac{1}{2\pi}\log\left(1-|\boldsymbol{x}_{j}|^{2}\right)+\frac{|\boldsymbol{x}_{j}|^{2}}{2\pi}-\frac{3}{8\pi}\,.\end{split}$
(3.0.21)

For an arbitrary cell pattern
$\\{{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{m}\\}}$, (3.0.21) is used to
evaluate the Neumann Green’s matrix ${\mathcal{G}}_{0}$ as needed in (3.0.16).

##  4 A ring and center cell pattern

With Sel’kov reaction kinetics, we apply the theory developed in section 2 to
a ring and center cell configuration in the unit disk. This pattern is
characterized by $m-1\geq 2$ equally spaced cells on a concentric ring within
the unit disk, and with one at the center of the disk (see Fig. 3). For this
pattern, the GCEP (2.2.13) will be used to obtain tractable nonlinear
algebraic equations that can be solved numerically to compute HB boundaries in
the $\tau$ versus $D$ parameter plane. In addition, a winding number criterion
is developed to count the number of unstable eigenvalues in open regions of
this parameter plane. Some of our examples will show that rather small changes
in either the permeabilities or reaction kinetic parameters of the center cell
can significantly alter the region in parameter space where oscillations
occurs.

Figure 3: Schematic showing a ring and center cell configuration of $m=9$
(left) and $m=3$ (right) cells in the unit disk. The ring cells are identical
and equally spaced on a concentric ring within the disk (in cyan). The center
cell (in red), possibly has different parameters. The red dots represent the
signaling molecules secreted in the bulk region by the cells.

###  4.1 Analysis of the GCEP

Consider a ring and center cell pattern of $m$ cells where the cells on the
ring of radius $r_{0}$ have identical parameters, but where the center cell
has possibly different permeabilities or Sel’kov kinetic parameters (see Fig.
3). The cell centers are at

$\boldsymbol{x}_{j}=r_{0}\left(\cos\left(\frac{2\pi(j-1)}{(m-1)}\right),\,\,\sin\left(\frac{2\pi(j-1)}{(m-1)}\right)\right)\,,\quad
j=1\,,\ldots,m-1\,;\qquad\boldsymbol{x}_{m}={\bf 0}\,,$ (4.1.1)

where $0<r_{0}<1$. For this pattern, the reduced-wave Green’s matrix
$\mathcal{G}$ in (2.1.7) can be partitioned as

$\mathcal{G}=\left(\begin{array}[]{ccc:c}&&&g_{m}\\\
&\mathcal{G}_{m-1}&&\vdots\\\ &&&g_{m}\\\
\hdashline[2pt/2pt]g_{m}&\dots&g_{m}&R_{m}\end{array}\right)\,;\qquad
g_{m}\equiv
G(\boldsymbol{x}_{j},\boldsymbol{x}_{m})=G(\boldsymbol{x}_{m},\boldsymbol{x}_{j})\,,\quad
j=1\,\ldots,m-1\,,\quad R_{m}=R(\boldsymbol{x}_{m})\,.$ (4.1.2)

Here ${\mathcal{G}}_{m-1}$ is the $(m-1)\times(m-1)$ symmetric matrix block
representing the interaction between the cells on the ring. Since this block
is also cyclic it has the eigenpair

${\mathcal{G}}_{m-1}\boldsymbol{e}=\omega_{1}\boldsymbol{e}\,,\qquad\mbox{with}\qquad\boldsymbol{e}=(1,\dots,1)^{T}\in{\mathbb{R}}^{m-1}\,,\quad\mbox{and}\quad\omega_{1}\equiv
R_{1}+\sum_{j=2}^{m-1}G(\boldsymbol{x}_{1},\boldsymbol{x}_{j})\,.$ (4.1.3)

In (4.1.2) there is a common interaction, represented by $g_{m}$, between each
ring cell and the center cell owing to the rotational symmetry and the fact
that the ring cells are all equidistant from the center cell.

For the identical ring cells, we label their permeabilities as $d_{1}=d_{1j}$
and $d_{2}=d_{2j}$ for $j=1,\dots,m-1$ and their common Sel’kov kinetic
parameters as $\mu_{1}=\mu_{j}$, $\alpha_{1}=\alpha_{j}$ and
$\zeta_{1}=\zeta_{j}$ for $j=1,\ldots,m-1$. Since the unique steady-state
solution to (2.2.21b) has the form
${\mathcal{A}}=(A_{c},\ldots,A_{c},A_{m})^{T}$, we readily find that $A_{c}$
and $A_{m}$ satisfies the $2\times 2$ linear system

$\displaystyle\left(1+\frac{\nu D}{d_{1}}+2\pi\nu\omega_{1}+\frac{2\pi\nu
D}{\tau}\frac{d_{2}}{d_{1}}\right)$ $\displaystyle A_{c}+2\pi\nu
g_{m}\,A_{m}=-\nu\mu_{1}\frac{d_{2}}{d_{1}}\,;$ (4.1.4a)
$\displaystyle\left(1+\frac{\nu D}{d_{1m}}+2\pi\nu R_{m}+\frac{2\pi\nu
D}{\tau}\frac{d_{2m}}{d_{1m}}\right)$ $\displaystyle
A_{m}+2\pi\nu\,g_{m}\,(m-1)A_{c}=-\nu\mu_{m}\frac{d_{2m}}{d_{1m}}\,,$ (4.1.4b)

where $\omega_{1}$ is the eigenvalue of ${\mathcal{G}}_{m-1}$ in (4.1.3). In
terms of ${\mathcal{A}}=(A_{c},\ldots,A_{c},A_{m})^{T}$, the steady-state for
the intracellular species as obtained from (2.2.21a) is

$\begin{split}u_{ej}^{1}=\begin{cases}u_{e1}^{1}\equiv\mu_{1}+\frac{2\pi
D}{\tau}\,A_{c}\,,&j=1\,,\ldots,m-1\,,\\\ u_{em}^{1}\equiv\mu_{m}+\frac{2\pi
D}{\tau}\,A_{m}\,,&j=m\,,\end{cases}\,;\qquad
u_{ej}^{2}=\begin{cases}u_{e1}^{2}\equiv\frac{\mu_{1}}{\alpha_{1}\,\,+\,\,(u_{e1}^{1})^{2}}\,,&j=1\,,\ldots,m-1\,,\\\
u_{em}^{2}\equiv\frac{\mu_{m}}{\alpha_{m}\,\,+\,\,(u_{em}^{1})^{2}}\,,&j=m\,.\end{cases}\end{split}$
(4.1.5)

Next, we determine the GCEP for the ring and center cell pattern using
(2.2.13a). For this pattern, the GCEP matrix $\mathcal{M}(\lambda)$ in
(2.2.13a) is written as

$\mathcal{M}(\lambda)=2\pi\nu\,\mathcal{G}_{\lambda}+\mathcal{M}_{0}\,,$
(4.1.6a) where $\mathcal{G}_{\lambda}$ is the eigenvalue-dependent Green’s
matrix, as defined in (2.2.10), and where the diagonal $\mathcal{M}_{0}$ is
defined by $\begin{split}\mathcal{M}_{0}=\begin{pmatrix}M_{0}&&&&\\\
&\ddots&&&\\\ &&&M_{0}&\\\
&&&&M_{m}\end{pmatrix}\,,\qquad\mbox{where}\qquad\begin{cases}M_{0}\equiv
1+\frac{\nu D}{d_{1}}+\frac{2\pi\nu
D}{\tau}\frac{d_{2}}{d_{1}}\mathit{K}_{c}\,,\\\ \\\ M_{m}\equiv 1+\frac{\nu
D}{d_{1m}}+\frac{2\pi\nu
D}{\tau}\frac{d_{2m}}{d_{1m}}\mathit{K}_{m}\,.\end{cases}\end{split}$ (4.1.6b)

Here $\mathit{K}_{c}=\mathit{K}_{c}(\lambda)$ and
$\mathit{K}_{m}=\mathit{K}_{m}(\lambda)$ are the entries of the $m\times m$
diagonal matrix
$\mathcal{K}=\text{diag}(\mathit{K}_{c},\dots,\mathit{K}_{c},\mathit{K}_{m})$
defined in (2.2.12a). For the Sel’kov kinetics given in (2.2.18), (2.2.12a)
yields

$\mathit{K}_{c}\equiv\frac{\lambda+\mbox{det}(J_{1})}{\lambda^{2}-\mbox{tr}(J_{1})\lambda+\mbox{det}(J_{1})}\,,\qquad\mathit{K}_{m}\equiv\frac{\lambda+\mbox{det}(J_{m})}{\lambda^{2}-\mbox{tr}(J_{m})\lambda+\mbox{det}(J_{m})}\,,$
(4.1.7a) where the trace and determinant of the Jacobians of the intracellular
dynamics for the identical ring cells and the center cell are given in terms
of the steady-state values in (4.1.5) by
$\begin{split}\mbox{det}(J_{1})=\zeta_{1}\left(\alpha_{1}+(u_{e1}^{1})^{2}\right)\,,\qquad\mbox{det}(J_{m})&=\zeta_{m}\left(\alpha_{m}+(u_{em}^{1})^{2}\right)\,,\\\
\mbox{tr}(J_{1})=\frac{\left[2(u_{e1}^{1})\mu_{1}-\left(\alpha_{1}+(u_{e1}^{1})^{2}\right)-\zeta_{1}\left(\alpha_{1}+(u_{e1}^{1})^{2}\right)^{2}\right]}{\alpha_{1}+(u_{e1}^{1})^{2}}\,,\quad\mbox{tr}(J_{m})&=\frac{\left[2(u_{em}^{1})\mu_{m}-\left(\alpha_{m}+(u_{em}^{1})^{2}\right)-\zeta_{m}\left(\alpha_{m}+(u_{em}^{1})^{2}\right)^{2}\right]}{\alpha_{m}+(u_{em}^{1})^{2}}\,.\end{split}$
(4.1.7b)

From Proposition 1, discrete eigenvalues $\lambda$ associated with the ring
and center cell pattern are roots of $\mbox{det}{\mathcal{M}}(\lambda)=0$. A
convenient way to implement this determinant root-finding problem numerically
is to use the special structure of ${\mathcal{M}}(\lambda)$ in order to
determine explicit formulae for its matrix spectrum
${\mathcal{M}}(\lambda)\boldsymbol{c}_{j}=\sigma_{j}\boldsymbol{c}_{j}$, for
$j=1,\ldots,m$, where $\sigma_{j}=\sigma_{j}(\lambda)$. Then, we need only
numerically solve the scalar root-finding problems $\sigma_{j}(\lambda)=0$ for
$\lambda$ for each $j=1,\ldots,m$.

To do so, we use the convenient fact that ${\mathcal{G}}_{\lambda}$ can be
partitioned, similar to that in (4.1.2), as

$\mathcal{G}_{\lambda}=\left(\begin{array}[]{ccc:c}&&&g_{\lambda m}\\\
&\mathcal{G}_{\lambda(m-1)}&&\vdots\\\ &&&g_{\lambda m}\\\
\hdashline[2pt/2pt]g_{\lambda m}&\dots&g_{\lambda m}&R_{\lambda
m}\end{array}\right)\,;\qquad g_{\lambda m}\equiv
G_{\lambda}(\boldsymbol{x}_{j},\boldsymbol{x}_{m})=G_{\lambda}(\boldsymbol{x}_{m},\boldsymbol{x}_{j})\,,\quad
j=1\,\ldots,m-1\,,\quad R_{\lambda m}\equiv
R_{\lambda}(\boldsymbol{x}_{m})\,,$ (4.1.8)

where $G_{\lambda}$ is the eigenvalue-dependent reduced-wave Green’s function
with regular part $R_{\lambda}$ satisfying (2.2.8). In (4.1.8), the
$(m-1)\times(m-1)$ matrix block $\mathcal{G}_{\lambda(m-1)}$, representing
cell interactions on the ring, is symmetric and cyclic. As a result, it has
the well-defined eigenspace

$\mathcal{G}_{\lambda(m-1)}\,\boldsymbol{v}_{j}=\omega_{\lambda
j}\,\boldsymbol{v}_{j},\quad
j=1,\dots,m-1\,;\qquad\boldsymbol{e}^{T}\boldsymbol{v}_{j}=0\,,\quad
j=1,\ldots,m-2\,,\quad\boldsymbol{v}_{m-1}=\boldsymbol{e}\equiv(1,\ldots,1)^{T}\in{\mathbb{R}}^{m-1}\,.$
(4.1.9)

By using this special matrix structure, it readily follows that the GCEP
matrix ${\mathcal{M}}(\lambda)$ in (4.1.6) admits $m-2$ anti-phase modes,
characterized by eigenvectors of the form
$\boldsymbol{c}_{j}=(\boldsymbol{v}_{j},0)^{T}\in{\mathbb{R}}^{m}$ where
$\boldsymbol{v}_{j}\in{\mathbb{R}}^{m-1}$ are those eigenvectors of
${\mathcal{G}}_{\lambda(m-1)}$ in (4.1.9), which satisfy
$\boldsymbol{e}^{T}\boldsymbol{v}_{j}=0$ for $j=1,\ldots,m-2$. With this
choice, (4.1.6) becomes

$\begin{pmatrix}2\pi\nu\mathcal{G}_{\lambda(m-1)}\boldsymbol{v}_{j}+M_{0}\boldsymbol{v}_{j}\\\
\\\ 2\pi\nu g_{\lambda
m}\,\boldsymbol{e}^{T}\boldsymbol{v}_{j}\end{pmatrix}=\begin{pmatrix}\sigma_{j}\boldsymbol{v}_{j}\\\
\\\ 0\end{pmatrix}\,,\quad\mbox{for}\quad j=1,\ldots,m-2\,.$ (4.1.10)

Since $\boldsymbol{e}^{T}\boldsymbol{v}_{j}=0$ for $j=1,\ldots,m-2$, we obtain
from (4.1.10), (4.1.9) and (4.1.6b) that $m-2$ eigenpairs of
${\mathcal{M}}(\lambda)$ are

$\sigma_{j}=2\pi\nu\omega_{\lambda j}+M_{0}=2\pi\nu\omega_{\lambda
j}+1+\frac{\nu D}{d_{1}}+\frac{2\pi\nu
D}{\tau}\frac{d_{2}}{d_{1}}\mathit{K}_{c}\,,\quad\boldsymbol{c}_{j}=(\boldsymbol{v}_{j},0)^{T}\,,\quad\mbox{for}\quad
j=1,\dots,m-2\,,$ (4.1.11)

where $\mathit{K}_{c}=\mathit{K}_{c}(\lambda)$ is defined in (4.1.7a). We
remark that since $\mathit{K}_{c}$ depends on the steady-state values $A_{c}$
and $A_{m}$, as obtained from the linear system (4.1.4), this term depends on
the permeabilities and local kinetics of the center cell. Discrete eigenvalues
$\lambda$ for the anti-phase modes are union of the zeroes of
$\sigma_{j}(\lambda)=0$ for $j=1,\ldots,m-2$.

For the remaining two eigenpairs of the GCEP matrix ${\mathcal{M}}(\lambda)$
the associated eigenvector $\boldsymbol{c}$ has the form
$\boldsymbol{c}=(\boldsymbol{e},\gamma)^{T}$ for some scalar $\gamma$ to be
determined and $\boldsymbol{e}=(1,\ldots,1)^{T}\in{\mathbb{R}}^{m-1}$. This
eigenvector is referred to as the in-phase mode since any instability
associated with this mode has the same phase for the cells on the ring. With
this choice, (4.1.6) reduces to

$\begin{pmatrix}\omega_{\lambda(m-1)}+{M_{0}/(2\pi\nu)}&g_{\lambda m}\\\
(m-1)g_{\lambda m}&R_{\lambda
m}+{M_{m}/(2\pi\nu)}\end{pmatrix}\begin{pmatrix}1\\\
\gamma\end{pmatrix}=\frac{\sigma}{2\pi\nu}\begin{pmatrix}1\\\
\gamma\end{pmatrix}\,,$ (4.1.12)

where
${\mathcal{G}}_{\lambda(m-1)}\boldsymbol{e}=\omega_{\lambda(m-1)}\boldsymbol{e}$
from (4.1.9).

Upon eliminating $\sigma$ from the $2\times 2$ matrix problem (4.1.12) we
obtain that $\gamma_{\pm}$ are the roots of the quadratic equation

$\gamma^{2}+\frac{1}{2\pi\nu g_{\lambda
m}}\Big{[}(M_{0}-M_{m})+2\pi\nu(\omega_{\lambda(m-1)}-R_{\lambda
m})\Big{]}\gamma-(m-1)=0\,,$ (4.1.13)

given by

$\begin{split}\gamma_{\pm}=-\frac{\beta_{\lambda}}{2}\pm\frac{1}{2}\sqrt{\beta^{2}_{\lambda}+4(m-1)}\,,\qquad\mbox{where}\qquad\beta_{\lambda}\equiv\frac{1}{2\pi\nu
g_{\lambda m}}\Big{[}(M_{0}-M_{m})+2\pi\nu(\omega_{\lambda(m-1)}-R_{\lambda
m})\Big{]}\,.\end{split}$ (4.1.14)

Since $\gamma_{+}\gamma_{-}=-(m-1)>0$, but with $\gamma_{\pm}$ possibly
complex-valued, we confirm that the two possible in-phase modes
$\boldsymbol{c}_{\pm}=(1,\dots,1,\gamma_{\pm})^{T}$ are orthogonal. The two
eigenvalues $\sigma=\sigma_{\pm}(\lambda)$, given by $\sigma_{\pm}\equiv
2\pi\nu\big{(}\omega_{\lambda(m-1)}+\gamma_{\pm}\,g_{\lambda m}\big{)}+M_{0}$,
can be written as

$\sigma_{\pm}(\lambda)=\frac{(h_{1}+h_{2})}{2}\pm\frac{1}{2}\sqrt{\left(h_{1}-h_{2}\right)^{2}+16\pi^{2}\nu^{2}(m-1)g_{\lambda
m}^{2}}\,,\quad\mbox{where}\quad h_{1}\equiv
2\pi\nu\,\omega_{\lambda(m-1)}+M_{0}\,,\quad h_{2}\equiv 2\pi\nu\,R_{\lambda
m}+M_{m}\,.$ (4.1.15)

Here $M_{0}=M_{0}(\lambda)$ and $M_{m}=M_{m}(\lambda)$ are given in (4.1.6b)
and $\omega_{\lambda(m-1)}$ is defined by
$\mathcal{G}_{\lambda(m-1)}\boldsymbol{e}=\omega_{\lambda(m-1)}\boldsymbol{e}$
from (4.1.9).

Alternatively, rather than solving (4.1.12) for $\sigma_{\pm}$, and then
setting $\sigma_{\pm}(\lambda)=0$ by using a root-finder for $\lambda$, we can
more directly conclude that ${\mathcal{M}}(\lambda)\boldsymbol{c}=\bf{0}$ for
$\boldsymbol{c}=(\boldsymbol{e},\gamma)^{T}$ if and only if the determinant of
the $2\times 2$ matrix in (4.1.12) vanishes. In this way, a discrete
eigenvalue $\lambda$ of the GCEP (2.2.13) for the in-phase mode, which
satisfies $\mbox{det}({\mathcal{M}}(\lambda))=0$, is a root of
${\mathcal{H}}_{c}(\lambda)=0$ defined by

${\mathcal{H}}_{c}(\lambda)=\left(\omega_{\lambda(m-1)}+\frac{M_{0}}{2\pi\nu}\right)\left(R_{\lambda
m}+\frac{M_{m}}{2\pi\nu}\right)-(m-1)g_{\lambda m}^{2}\,.$ (4.1.16)

For any root of (4.1.16), the corresponding eigenvector $\boldsymbol{c}$ of
${\mathcal{M}}(\lambda)$ is

$\boldsymbol{c}=(1,\ldots,1,\gamma)^{T}\,,\qquad\gamma=-\frac{1}{g_{\lambda
m}}\left(\omega_{\lambda(m-1)}+\frac{M_{0}}{2\pi\nu}\right)\,.$ (4.1.17)

It is readily verified using (4.1.15) that if $\lambda^{\star}$ satisfies
either $\sigma_{+}(\lambda^{\star})=0$ or $\sigma_{-}(\lambda^{\star})=0$,
then we must have ${\mathcal{H}}_{c}(\lambda^{\star})=0$. In contrast, if
$\lambda^{\star}$ satisfies ${\mathcal{H}}_{c}(\lambda^{\star})=0$, then we
can only conclude that either $\sigma_{+}(\lambda^{\star})=0$ or
$\sigma_{-}(\lambda^{\star})=0$. Therefore, in implementing a root-finding
strategy based on the single scalar equation (4.1.16) instead of the two
scalar equations (4.1.15) care must be taken to identify all possible roots of
${\mathcal{H}}_{c}(\lambda)=0$ for the same parameter set.

We summarize our result for eigenvalues $\lambda$ of the GCEP (2.2.13) for a
ring and center cell pattern as follows:

######  Proposition 3.

Consider a ring and center hole pattern of $m\geq 3$ cells in the unit disk
with cell centers at (4.1.1). The set $\Lambda({\mathcal{M}})$ as obtained
from the GCEP (2.2.13), and which approximates as $\varepsilon\to 0$ all the
discrete eigenvalues of the linearization of the PDE-ODE system (1.0.2) around
the steady-state solution, is

$\Lambda({\mathcal{M}})\equiv\Big{\\{}\lambda\,\,\big{|}\,\,\bigcup\limits_{j=1}^{m-2}\\{{\,\sigma_{j}(\lambda)=0\,\\}}\,,\,\,\bigcup\,\\{{\sigma_{\pm}(\lambda)=0\\}}\Big{\\}}\,.$
(4.1.18)

Here $\sigma_{j}(\lambda)$, for $j=1,\ldots,m-2$, for the anti-phase modes are
defined in (4.1.11), while $\sigma_{\pm}(\lambda)$ for the in-phase modes are
defined in (4.1.15). As shown in Remark 1 of Appendix A, due to mode
degeneracy of ${\mathcal{G}}_{\lambda(m-1)}$, there are ${(m-1)/2}$ distinct
anti-phase modes if $m$ is odd and ${(m-2)/2}$ distinct anti-phase modes if
$m$ is even.

For the unit disk, where an infinite series representation of the solution to
(2.2.8) is available, explicit formulae for the eigenvalues $\omega_{\lambda
j}$ of ${\mathcal{G}}_{\lambda(m-1)}$, as needed in (4.1.11) and (4.1.15), are
given in Appendix A. The result in Remark 1 of Appendix A regarding mode
degeneracy results from the fact that ${\mathcal{G}}_{\lambda(m-1)}$ is both
symmetric and cyclic. In Appendix A we also show how to readily calculate the
quantities in (4.1.2) and (4.1.3), which are needed in (4.1.4) for determining
the steady-state.

As indicated by Proposition 1, stability boundaries in the $\tau$ versus $D$
parameter space for the steady-state under Sel’kov kinetics are determined by
HB boundaries where $\lambda=i\lambda_{I}\in\Lambda({\mathcal{M}})$. To
determine Hopf bifurcation boundaries for the anti-phase modes we set
$\mbox{Re}(\sigma_{j}(i\lambda_{I}))=0$ and
$\mbox{Im}(\sigma_{j}(i\lambda_{I}))=0$ for $j=1,\ldots,m-2$ in (4.1.11). Upon
separating real and imaginary parts in (4.1.7a) we obtain the following
nonlinear algebraic system for each $j=1,\ldots,m-2$:

$2\pi\nu\,\mbox{Re}(\omega_{\lambda
j})+\left(1+\nu\frac{D}{d_{1}}\right)+\frac{2\pi\nu
D}{\tau}\frac{d_{2}}{d_{1}}\,\mbox{Re}(\mathit{K}_{c}(i\lambda_{I}))=0\,,\qquad\mbox{Im}(\omega_{\lambda
j})+\frac{D}{\tau}\frac{d_{2}}{d_{1}}\,\mbox{Im}(\mathit{K}_{c}(i\lambda_{I}))=0\,,$
(4.1.19a) where
$\mbox{Im}(\mathit{K}_{c}(i\lambda_{I}))=\frac{\lambda_{I}(\mbox{det}(J_{1})-\lambda_{I}^{2})\,+\lambda_{I}\mbox{det}(J_{1})\mbox{tr}(J_{1})}{(\mbox{det}(J_{1})-\lambda_{I}^{2})^{2}+(\lambda_{I}\,\mbox{tr}(J_{1}))^{2}}\,,\quad\mbox{Re}(\mathit{K}_{c}(i\lambda_{I}))=\frac{\mbox{det}(J_{1})(\mbox{det}(J_{1})-\lambda_{I}^{2})-\lambda_{I}^{2}\,\mbox{tr}(J_{1})}{(\mbox{det}(J_{1})-\lambda_{I}^{2})^{2}+(\lambda_{I}\,\mbox{tr}(J_{1}))^{2}}\,.$
(4.1.19b)

Here $J_{1}$ is the Jacobian of the Sel’kov kinetics for the ring cells with
determinant and trace given in (4.1.7b). Similarly, the Hopf bifurcation
boundaries for the in-phase modes are obtained by setting
$\sigma_{\pm}(i\lambda_{I})=0$ in (4.1.15), which yields the nonlinear
algebraic system

$\mbox{Re}(\sigma_{\pm}(i\lambda_{I}))=0\,,\qquad\mbox{Im}(\sigma_{\pm}(i\lambda_{I}))=0\,,$
(4.1.20)

or equivalently ${\mathcal{H}}_{c}(i\lambda_{I})=0$ from (4.1.16).

###  4.2 Example: Two cells on a ring with a center cell

We now apply the theory developed in section 4.1 to a population of $m=3$
cells, where two of the cells are equally spaced on a concentric ring of
radius $r_{0}$ within the unit disk, with the remaining one centered at the
origin (see Fig. 3). For this configuration, the eigenvalues of the $3\times
3$ GCEP matrix $\mathcal{M}(\lambda)$ are given in (4.1.11) and (4.1.15) for a
single anti-phase mode and the two in-phase modes, respectively. To compute
the HB boundaries in the $\tau$ versus $D$ parameter plane for these modes, we
solve (4.1.19) and (4.1.20) numerically by implementing the psuedo-arclength
continuation algorithm TEST_CON (cf. [1]) with respect to $D$, while using
Newton’s method to compute $\tau$ and $\lambda_{I}$ at each point on the
solution path. Such a continuation scheme in $D$ is needed owing to the
possibility of fold points along the HB boundary.

To determine regions of instability in open sets of the $\tau$ versus $D$
parameter plane we use (4.1.18) of Proposition 3, together with the winding
number criterion of complex analysis, to identify the number ${\mathcal{N}}$
of eigenvalues $\lambda\in\Lambda({\mathcal{M}})$ with $\mbox{Re}(\lambda)>0$.
To do so, we first define
$\mathcal{F}(\lambda)\equiv\det(\mathcal{M}(\lambda))$, where
${\mathcal{M}}(\lambda)$ is the GCEP matrix in (4.1.6). Provided that there
are no zeroes or poles on the imaginary axis, $\mathcal{N}$ is the number of
zeroes of $\mathcal{F}(\lambda)=0$ in $\mbox{Re}(\lambda)>0$, which from the
argument principle is

$\mathcal{N}=\frac{1}{2\pi}\big{[}\text{arg}\,\mathcal{F}(\lambda)\big{]}_{\Gamma}+{\mathcal{P}}\,.$
(4.2.1)

Here ${\mathcal{P}}$ is the number of poles of $\mathcal{F}(\lambda)$ in
$\mbox{Re}(\lambda)>0$, while
$\big{[}\text{arg}\,\mathcal{F}(\lambda)\big{]}_{\Gamma}$ denotes the change
in the argument of $\mathcal{F}(\lambda)$ over the closed, counter-clockwise
oriented contour $\Gamma$. This contour $\Gamma$ is the limit as
$\mathcal{R}\to\infty$ of the union of the imaginary axis
$\Gamma_{I}=i\lambda_{I}$, for $|\lambda_{I}|\leq\mathcal{R}$, and the semi-
circle $\Gamma_{\mathcal{R}}$, defined by $|\lambda|=\mathcal{R}$ with
$|\text{arg}(\lambda)|\leq\pi/2$. To count the number of poles of
$\mathcal{F}(\lambda)$ in $\mbox{Re}(\lambda)>0$, we must examine the
analyticity properties of the GCEP matrix $\mathcal{M}(\lambda)$ in (4.1.6).
Since the entries of the Green’s matrix $\mathcal{G}_{\lambda}$ are analytic
in $\mbox{Re}(\lambda)>0$, any singularity of $\mathcal{F}(\lambda)$ must
arise from the diagonal matrix
${\mathcal{K}}(\lambda)\equiv\mbox{diag}(\mathit{K}_{c},\mathit{K}_{c},\mathit{K}_{m})$
of (4.1.6b), which is given explicitly in (4.1.7) in terms of the Jacobians
$J_{1}$ and $J_{m}$ of the Sel’kov kinetics for the identical ring cells and
the center cell, respectively. Since $\mbox{det}(J_{1})>0$, (4.1.7a) yields
that $\mathit{K}_{c}$ has a complex conjugate pair of poles in
$\mbox{Re}(\lambda)>0$ only if $\mbox{tr}(J_{1})>0$. Since this term involves
two rows of ${\mathcal{M}}(\lambda)$, ${\mathcal{P}}$ must be incremented by
four whenever $\mbox{tr}(J_{1})>0$. Similarly, since $\mbox{det}(J_{m})>0$,
${\mathcal{P}}$ is increased by two when $\mbox{tr}(J_{m})>0$ for the center
cell.

With ${\mathcal{P}}$ determined in this way, we numerically compute the number
of unstable eigenvalues $\mathcal{N}$ of the linearization of the steady-state
by evaluating (4.2.1) for each point $(D,\tau)$ in the $\tau$ versus $D$
plane. For each such point, we numerically construct the closed contour
$\Gamma$ for some value $\mathcal{R}$ (chosen so that $\Gamma$ encloses all
the poles of $\mathcal{F}$) in the complex $\lambda$-plane. As the closed
curve $\Gamma$ is traversed in a counter-clockwise direction, the closed image
curve $\mathcal{F}(\lambda)=\mathcal{F}_{R}+i\mathcal{F}_{I}$ is evaluated
numerically in the complex $\mathcal{F}$-plane. The winding number, denoting
the number of times $\mathcal{F}$ encloses/winds around the origin, is
computed numerically from the algorithm of [2], and this is used to calculate
$\big{[}\text{arg}\,\mathcal{F}(\lambda)\big{]}_{\Gamma}$. If the orientation
of ${\mathcal{F}}(\lambda)={\mathcal{F}}_{R}+i{\mathcal{F}}_{I}$ around the
origin is in the counter-clockwise direction, then
$\big{[}\text{arg}\,\mathcal{F}(\lambda)\big{]}_{\Gamma}$ is positive;
otherwise, it is negative. In our computations, we chose $\mathcal{R}=1.5$
(since any pole of $\mathcal{F}$ is close to the imaginary axis of the
$\lambda$-plane), and we discretized the closed contour
$\Gamma=\Gamma_{I}\cup\Gamma_{\mathcal{R}}$ into subintervals, with
$\Gamma_{I}$ having 800 subintervals while the semi-circle
$\Gamma_{{\mathcal{R}}}$ had $50$ subintervals. In the function evaluation,
the identity
${\mathcal{F}}(\overline{\lambda})=\overline{{\mathcal{F}}(\lambda)}$ was used
to halve the computational effort.

(a) HB boundaries: Identical cells

(b) Number of unstable eigenvalues of the GCEP

Figure 4: Left panel: HB boundaries in the $\tau$ versus $D$ plane for a ring
and center hole pattern of $m=3$ identical cells with ring radius $r_{0}=0.5$,
parameters as in (4.3.1) and permeabilities $d_{1}=0.8$ and $d_{2}=0.2$. The
dashed and heavy solid curves are for the in-phase modes computed from
(4.1.20) with $(+)$ and $(-)$, respectively. The thin solid curve is for the
anti-phase mode computed from (4.1.19). Each mode is linearly unstable within
its respective lobe. Linearly stable steady-state solutions exist outside the
union of the lobes. Full PDE simulations of (1.0.2) are shown in Figs. 6 and 7
at the red and blue dots, respectively. Right panel: The regions of
instability computed from the winding number (4.2.1). Blue region: in-phase
‘+’ mode is unstable with 2 roots of
$\mathcal{F}(\lambda)=\det(\mathcal{M}(\lambda))=0$ in $\mbox{Re}(\lambda)>0$.
Green region: anti-phase mode is also unstable, yielding 4 roots. Magenta
region: all modes are unstable, and there are 6 roots. Plot on right is a zoom
of the one on the left. The HB boundaries in the left panel are superimposed
in this figure.

###  4.3 A defective center cell: different permeabilities

In our results below, except when otherwise stated, the Sel’kov parameters
$\alpha$, $\mu$ and $\zeta$, and permeabilities $d_{1}$ and $d_{2}$ for the
identical cells on the ring, and the common cell radius $\varepsilon$ are

$\alpha=0.9\,,\qquad\mu=2\,,\qquad\zeta=0.15\,,\qquad d_{1}=0.8\,,\qquad
d_{2}=0.2\,,\qquad\varepsilon=0.05\,.$ (4.3.1)

From Fig. 2, we conclude that each cell, when isolated, has no intracellular
oscillations. The permeabilities $d_{13}$ and $d_{23}$ for the center cell
will be stated in the figure captions below.

Fig. 4(a) shows the computed HB boundaries in the $\tau$ versus $D$ plane for
a ring of radius $r_{0}=0.5$ when the cells are all identical. We observe that
one of the in-phase lobes is open/unbounded, which predicts the existence of
intracellular oscillations even for large $D$. In Fig. 4(b), we show the
corresponding regions of instability in the $\tau$ versus $D$ parameter plane.
In generating this figure, we pixelated the $\tau$ versus $D$ plane with the
uniform spacing $\Delta\tau=\Delta D=0.005$, and at each discrete point
$(D,\tau)$ used our winding number algorithm to count the number of roots
$\mathcal{N}$ of $\mbox{det}{\mathcal{M}}(\lambda)=0$ in
$\mbox{Re}(\lambda)>0$. In Fig. 4(b), each point in the blue-shaded region has
two unstable eigenvalues for the GCEP, and they correspond to the in-phase ‘+’
mode. The green-shaded region contains four unstable eigenvalues, two of which
are for the anti-phase mode while the other two for the in-phase ‘+’ mode.
Finally, in the magenta-shaded region there are six unstable eigenvalues of
the linearization, with two such eigenvalues associated with each of the three
possible modes of instability (in-phase $\pm$ and anti-phase). The HB
boundaries in the left panel of Fig. 4 are superimposed on these instability
regions.

The $\mathcal{F}(\lambda)={\mathcal{F}}_{R}+i{\mathcal{F}}_{I}$ curve in the
complex $\mathcal{F}$-plane is shown in Fig. 5 for a specific point in each of
the three instability regions in Fig. 4(b). These plots show how
${\mathcal{F}}(\lambda)$ winds around the origin
$({\mathcal{F}}_{R},{\mathcal{F}}_{I})=(0,0)$ (shown with a green dot) as
$\lambda$ traverses $\Gamma$ in the counterclockwise direction. For the point
$(D,\tau)=(0.05,0.15)$ in the magenta-shaded region, we observe from the left
panel of Fig. 5 that
$\big{[}\text{arg}\,\mathcal{F}(\lambda)\big{]}_{\Gamma}=0$. At this point,
${\mathcal{F}}(\lambda)$ has two poles, one of order four and the other of
order 2, so that ${\mathcal{P}}=6$. As such, (4.2.1) yields that there are 6
roots (counting multiplicity) to $\mathcal{F}(\lambda)=0$ in
$\mbox{Re}(\lambda)>0$. At the point $(D,\tau)=(0.15,0.35)$ in the anti-phase
mode instability region (green-shaded region in Fig. 4(b)),
${\mathcal{F}}(\lambda)$ winds round the origin twice in the clockwise
direction as shown in the middle panel of Fig. 5, which yields
$\big{[}\text{arg}\,\mathcal{F}(\lambda)\big{]}_{\Gamma}=-4\pi$. Since
${\mathcal{P}}=6$ at this point, (4.2.1) yields that $\mathcal{F}(\lambda)=0$
has four roots (counting multiplicity) in $\mbox{Re}(\lambda)>0$. In the right
panel of Fig. 5, we present a similar result for the point
$(D,\tau)=(0.4,0.5)$ in the blue-shaded region in Fig. 4(b). At this point, we
calculate $\big{[}\text{arg}\,\mathcal{F}(\lambda)\big{]}_{\Gamma}=-8\pi$ and
that ${\mathcal{F}}$ has a pole of order four and a pole of order two in
$\mbox{Re}(\lambda)>0$. As such, (4.2.1) yields that $\mathcal{F}(\lambda)=0$
has two roots in $\mbox{Re}(\lambda)>0$.

Figure 5: The closed curve $\mathcal{F(\lambda)}=\det(\mathcal{M}(\lambda))$
plotted for a specific points in each of the three regions of instability in
Fig. 4(b). The closed contour $\Gamma=\Gamma_{I}\cup\Gamma_{{\mathcal{R}}}$ in
the $\lambda$-plane was constructed with $\mathcal{R}=1.5$, and discretized
with $800$ subintervals on $\Gamma_{I}$ and $50$ subintervals on
$\Gamma_{{\mathcal{R}}}$. The green dot represents the origin
$({\mathcal{F}}_{R},{\mathcal{F}}_{I})=(0,0)$, the black dot indicates the
starting point of the curve, and the blue arrows show the direction of the
curve. The inserts detail the behavior near the origin. Left panel: for
$(D,\tau)=(0.05,0.15)$ in the magenta-shaded region of Fig. 4(b), we have
$\big{[}\text{arg}\,\mathcal{F}(\lambda)\big{]}_{\Gamma}=0$. Middle panel: for
$(D,\tau)=(0.15,0.35)$ in the green-shaded region in Fig. 4(b), we have
$\big{[}\text{arg}\,\mathcal{F}(\lambda)\big{]}_{\Gamma}=-4\pi$. Right panel:
for $(D,\tau)=(0.4,0.5)$ in the blue-shaded region of Fig. 4(b), we have
$\big{[}\text{arg}\,\mathcal{F}(\lambda)\big{]}_{\Gamma}=-8\pi$.

The real and imaginary parts of the normalized eigenvector $\boldsymbol{c}$ of
the GCEP matrix ${\mathcal{M}}$ in (4.1.6) is given in Table 1 for selected
points on the HB boundaries in Fig. 4(a). Recall from (2.2.17) that the
magnitude of the components of the eigenvector $\boldsymbol{c}$ measure the
diffusive flux at the boundary of each cell, while
$\tilde{\boldsymbol{c}}\equiv{\mathcal{K}}\boldsymbol{c}$ predicts the
relative amplitude and phase shifts of the intracellular oscillations within
the cells at the Hopf bifurcation point. For our ring and center-cell pattern
${\mathcal{K}}=\mbox{diag}\left(\mathit{K}_{c},\mathit{K}_{c},\mathit{K}_{m}\right)$,
where $\mathit{K}_{c}$ and $\mathit{K}_{m}$ are given in (4.1.7a).

mode | $(D,\tau)$ | j | $\Big{(}\mbox{Re}(c_{j}),\mbox{Im}(c_{j})\Big{)}$ | $\theta_{j}\,(\text{rad})$ | $\Big{(}\mbox{Re}(\tilde{c}_{j}),\mbox{Im}(\tilde{c}_{j})\Big{)}$
---|---|---|---|---|---
|  | 1 | $(0.579,0)$ | $0$ | $(-0.569,0.067)$
In-phase ($+$) | $(1.021,0.262)$ | 2 | $(0.579,0)$ | $0$ | $(-0.569,0.067)$
(dashed curve) |  | 3 | $(0.575,0.0144)$ | $0.0251$ | $(-0.582,0.0625)$
|  | 1 | $(-0.412,-0.004)$ | $3.15$ | $(0.395,-0.00429)$
In-phase ($-$) | $(0.0857,0.199)$ | 2 | $(-0.412,-0.004)$ | $3.15$ | $(0.395,-0.00429)$
(heavy solid) |  | 3 | $(0.813,0)$ | $0$ | $(-0.829,0.0239)$
|  | 1 | $(0.707,0)$ | $0$ | $(-0.706,0.042)$
Anti-phase | $(0.211,0.365)$ | 2 | $(-0.707,0)$ | $\pi$ | $(0.706,-0.042)$
(thin solid) |  | 3 | $(0,0)$ | $0$ | $(0,0)$
Table 1: Real and imaginary parts of the eigenvector $\boldsymbol{c}$ of the
GCEP matrix ${\mathcal{M}}(\lambda)$ in (4.1.6), together with
$\tilde{\boldsymbol{c}}\equiv{\mathcal{K}}\boldsymbol{c}$, as computed for a
few points on the HB boundaries in Fig. 4(a) for three identical cells. The
second to last column shows the phase shifts measured in terms of the angle
each component of the vector $\boldsymbol{c}$ makes with the positive real
axis in anticlockwise direction.

In Fig. 6 we show full numerical simulations of the coupled PDE-ODE model
(1.0.2) obtained using the commercial PDE software package FlexPDE [14] for
$\tau=0.55$ and $D=1$, which corresponds to the red dot in the phase diagram
of Fig. 4(a). We observe from the results in this figure that the
intracellular dynamics of the cells are synchronized with a very slight phase
shift, which agrees with the prediction by the eigenvector
${\mathcal{K}}\boldsymbol{c}$ in the first three rows of Table 1 from the
linearized theory. Although the cells have identical parameters, the center
and ring cells have slightly different dynamics owing to the fact that the
full Green’s matrix is not cyclic for a ring and center cell pattern. Our
numerical computations of $\det{\mathcal{M}}(\lambda)=0$ for
$(D,\tau)=(1,0.55)$ using the GCEP matrix in (4.1.6) yields that
$\mbox{Re}(\lambda)\approx 0.0143$, $\mbox{Im}(\lambda)\approx 0.762$,
$\mbox{Re}(\boldsymbol{c})\approx(0.588,0.588,0.556)$ and
$\mbox{Im}(\boldsymbol{c})\approx(0,0,0.0178)$. Observe that the eigenvector
is rather close to that on the nearby point on the HB boundary, as given in
the first three rows of Table 1. The prediction from linearized theory is that
the period of oscillations is approximately ${2\pi/\mbox{Im}(\lambda)}\approx
8.25$, which is rather close to the period observed in the full PDE
simulations of Fig. 6. A similar full numerical result is presented in Fig. 7
for $\tau=1.1$ and $D=1$, corresponding to the blue dot in Fig. 4(a). At this
pair $(D,\tau)$, our phase diagram predicts no intracellular oscillations.
This is confirmed from the full numerical results shown in Fig. 7.

Figure 6: Full PDE simulations of (1.0.2), computed with FlexPDE [14], for
$\tau=0.55$ and $D=1$ for three identical cells corresponding to the red dot
in Fig. 4(a). Left panel: surface plot at time $t=400$. Middle panel:
intracellular species $u_{1}$ versus $t$. Right panel: intracellular species
$u_{2}$ versus $t$. The blue and red curve is for the cells on the ring and
the center cell, respectively.

Figure 7: Same caption as in Fig. 6 except that now $\tau=1.1$ and $D=1$,
corresponding to the blue dot in Fig. 4(a). There are no intracellular
oscillations and the steady-state is linearly stable.

mode | $(D,\tau)$ | j | Re(${c}_{j}$) | Im(${c}_{j}$) | $\theta_{j}\,(\text{rad})$
---|---|---|---|---|---
|  | 1 | $0.345$ | $0.171$ | $0.460$
In-phase ($+$) | $(0.820,0.584)$ | 2 | $0.345$ | $0.171$ | $0.460$
(dashed curve) |  | 3 | $0.839$ | $0$ | $0$
|  | 1 | $-0.446$ | $0.256$ | $2.62$
In-phase ($-$) | $(0.0775,0.216)$ | 2 | $-0.446$ | $0.256$ | $2.62$
(heavy solid) |  | 3 | $0.686$ | $0$ | $0$
|  | 1 | $0.707$ | $0$ | $0$
Anti-phase | $(0.129,0.325)$ | 2 | $-0.707$ | $0$ | $\pi$
(thin solid) |  | 3 | $0$ | $0$ | $0$

Figure 8: Left panel: same caption as in Fig. 4(a) except that now the center
cell is a defector with permeabilities $d_{13}=0.4$ and $d_{23}=0.2$,
corresponding to a reduced influx into the center cell. Full PDE simulations
are shown in Figs. 9 and 10 at the red and blue dots, respectively. Right
panel: same caption as in Table 1 except that now $d_{13}=0.4$ and
$d_{23}=0.2$. The real and imaginary parts of the eigenvector $\boldsymbol{c}$
of the GCEP matrix ${\mathcal{M}}(\lambda)$ in (4.1.6) are computed for a few
points on the HB boundaries shown in the left panel.

Figure 9: Full PDE simulations of (1.0.2), computed with FlexPDE [14], for
$\tau=0.6$ and $D=0.4$, corresponding to the red dot in the left panel of Fig.
8. The center cell is a defector with permeabilites $d_{13}=0.4$ and
$d_{23}=0.2$. Left panel: surface plot at time $t=400$. Middle panel:
intracellular species $u_{1}$ versus $t$. Right panel: intracellular species
$u_{2}$ versus $t$.

Figure 10: Same caption as in Fig. 9, except that the FlexPDE [14] simulation
of (1.0.2) is done at $\tau=0.2$ and $D=0.05$, corresponding to the blue dot
in the left panel of Fig. 8. At this point, both in-phase modes and the anti-
phase mode are unstable. The beating-behavior observed arises from the fact
that these three modes all have comparable frequencies.

In the left panel of Fig. 8 we show the computed HB boundaries in the $\tau$
versus $D$ plane for the same parameters as in Fig. 4(a), except that the
influx rate into the center cell is reduced to $d_{13}=0.4$ (keeping
$d_{23}=0.2$). With this lower influx rate, Fig. 8 shows that the in-phase
lobes are now closed, so that there are no longer intracellular oscillations
when $D$ is large. Qualitatively, when $D$ is large, the bulk chemical
diffuses quickly in the entire disk and there is insufficient feedback of it
into the center cell, owing to the smaller value of $d_{13}$, to sustain
intracellular oscillations. The winding number algorithm of (4.2.1) can be
used, with the same result as in Fig. 4, to determine the number of unstable
eigenvalues of the GCEP within the lobes (not shown). In the right panel of
Fig. 8 we give the real and imaginary parts of the normalized eigenvector
$\boldsymbol{c}$ of the GCEP matrix in (4.1.6) at a few selected points on the
HB boundaries shown in the left panel of Fig. 8.

In Fig. 9, we show full FlexPDE [14] numerical simulations of (1.0.2) for
$\tau=0.6$ and $D=0.4$, which corresponds to the red dot in the left panel of
Fig. 8. Our numerical computations of $\det{\mathcal{M}}(\lambda)=0$ for
$(D,\tau)=(0.4,0.6)$ using the GCEP matrix in (4.1.6) yields that
$\mbox{Re}(\lambda)\approx 0.0115$, $\mbox{Im}(\lambda)\approx 0.778$,
$\mbox{Re}({\mathcal{K}}\boldsymbol{c})\approx(-0.551,-0.551,-0.624)$ and
$\mbox{Im}({\mathcal{K}}\boldsymbol{c})\approx(0.0197,0.0197,0.0569)$. As
such, our linearized theory predicts that the center cell will have larger
amplitude oscillations near onset than the ring cells, and there will be a
$\approx 22^{\circ}$ phase shift between the oscillations. Our FlexPDE
numerical results in the middle and right panels of Fig. 9 show that the
prediction of the linearized theory does extend to the fully nonlinear regime
in that the defective center cell has larger amplitude oscillations than do
the identical ring cells, and there is a slight phase shift in the
oscillations. From the surface plot shown in the left panel of Fig. 9, we
observe that the coupling between the cells mediated by the bulk medium is
rather weak. Moreover, the rather large concentration of the bulk chemical
close to the center cell, with flux measured by the modulus of the third
component of $\mbox{Re}(\boldsymbol{c})\approx(0.267,0.267,0.914)$ and
$\mbox{Im}(\boldsymbol{c})\approx(0.107,0.107,0.0)$, is due to its smaller
rate $d_{13}=0.4$ of influx into the center cell than for the identical ring
cells centered at $(\pm 0.5,0)$. Paradoxically, however, this large buildup of
the bulk signal near the center cell counteracts the relatively smaller rate
of influx into the center cell, and has the effect of triggering a larger
amplitude oscillation in the center cell than for the ring cells.

mode | $(D,\tau)$ | j | Re(${c}_{j}$) | Im(${c}_{j}$) | $\theta_{j}\,(\text{rad})$
---|---|---|---|---|---
|  | 1 | $0.707$ | $0$ | $0$
In-phase ($+$) | $(0.0931,0.165)$ | 2 | $0.707$ | $0$ | $0$
(dashed curve) |  | 3 | $0.0245$ | $-0.0314$ | $5.375$
|  | 1 | $-0.0638$ | $0.0985$ | $2.146$
In-phase ($-$) | $(1.241,3.577)$ | 2 | $-0.0638$ | $0.0985$ | $2.146$
(heavy solid) |  | 3 | $0.986$ | $0$ | $0$
|  | 1 | $0.707$ | $0$ | $0$
Anti-phase | $(0.0836,0.159)$ | 2 | $-0.707$ | $0$ | $\pi$
(thin solid) |  | 3 | $0$ | $0$ | $0$

Figure 11: Same caption as in Fig. 8 except that now the center cell is a
defector with permeabilities $d_{13}=0.2$ and $d_{23}=0.4$, corresponding to a
reduced influx and a larger efflux out of the center cell. Full PDE
simulations are shown in Fig. 12 at the blue dot in the left panel.

In Fig. 10, we show full FlexPDE [14] numerical simulations of (1.0.2) for
$\tau=0.2$ and $D=0.05$, corresponding to the blue dot in the left panel of
Fig. 8. As seen from Fig. 8, this point is located within the region of
instability that is common to all three modes of instability. By solving
$\det{\mathcal{M}}(\lambda)=0$ for $(D,\tau)=(0.05,0.2)$ numerically, the
eigenvalues $\lambda$ and eigenvectors $\boldsymbol{c}$ for the three modes
are:

$\begin{split}\mbox{in-phase (+):}\quad\lambda&\approx
0.0350+0.779\,i\,,\quad\mbox{Re}(c)\approx(0.152,0.152,0.829)\,,\quad\mbox{Im}(c)\approx(0,0,-0.517)\,,\\\
\mbox{in-phase (-):}\quad\lambda&\approx
0.0144+0.813\,i\,,\quad\mbox{Re}(c)\approx(0.624,0.624,-0.434)\,,\quad\mbox{Im}(c)\approx(0,0,-0.180)\,,\\\
\mbox{anti-phase:}\quad\lambda&\approx
0.0189+0.812\,i\,,\quad\mbox{Re}(c)\approx(0.707,-0.707,0)\,,\quad\mbox{Im}(c)=(0,0,0)\,.\\\
\end{split}$ (4.3.2)

We observe that this linearized theory predicts three distinct unstable
oscillatory modes with roughly similar frequencies $\mbox{Im}(\lambda)$ and
growth rates $\mbox{Re}(\lambda)$. The beating-type intracellular oscillations
observed in Fig. 10 is likely related to the well-known linear phenomenon of
superimposing two or more single-mode oscillations with comparable
frequencies.

Figure 12: Full PDE simulations of (1.0.2), computed with FlexPDE [14], for
$\tau=3$ and $D=0.5$, corresponding to the blue dot in the left panel of Fig.
11. The center cell is a defector with permeabilites $d_{13}=0.2$ and
$d_{23}=0.4$. Left panel: surface plot at time $t=400$. Observe the buildup of
the bulk chemical near the center cell in comparison to the two ring cells.
Middle panel: intracellular species $u_{1}$ versus $t$. Right panel:
intracellular species $u_{2}$ versus $t$. As predicted by the second row in
the table in the right panel of Fig. 11, we confirm that the center cell has
larger amplitude oscillations than do the ring cells.

In the left panel of Fig. 11 we plot the HB boundaries for the case where the
center cell has permeabilities $d_{13}=0.2$ and $d_{23}=0.4$. This corresponds
to a larger secretion or efflux rate for the center cell, while the feedback
it receives from the ring cells is reduced. Unlike the HB boundaries shown in
Figs. 4 and 8 where the instability regions for the modes are nested within
each other, we observe from the insert in the left panel of Fig. 11 that only
the in-phase $(+)$ mode and the anti-phase mode overlap. Moreover, since the
other in-phase $(-)$ lobe is unbounded in $D$, intracellular oscillations
always occur in some range of $\tau$ as $D$ increases. In the right panel of
Fig. 11, we give the real and imaginary parts of the normalized eigenvector
$\boldsymbol{c}$ of the GCEP matrix (4.1.6) at a few points on the HB
boundaries. From the second row of this table, the linearized theory suggests
that the amplitude of intracellular oscillations associated with the dominant
in-phase $(-)$ instability lobe will be much larger in the center cell than in
the identical ring cells at the HB point, and that there will be a significant
phase shift in the oscillations between the center cell and the ring cells.
More precisely, at the point $(D,\tau)=(0.5,3)$ interior to the instability
lobe, we solve $\mbox{det}{\mathcal{M}}(\lambda)=0$ numerically to obtain that
$\mbox{Re}(\lambda)\approx 0.00665$, $\mbox{Im}(\lambda)\approx 0.781$,
$\mbox{Re}(\boldsymbol{c})\approx(-0.0534,-0.0534,0.9908)$,
$\mbox{Im}(\boldsymbol{c})\approx(0.0792,0.0792,0)$,
$\mbox{Re}({\mathcal{K}}\boldsymbol{c})\approx(-0.177,-0.177,-0.926)$, and
$\mbox{Im}({\mathcal{K}}\boldsymbol{c})\approx(0.187,0.187,0.106)$. From
${\mathcal{K}}\boldsymbol{c}$ and $\mbox{Im}(\lambda)$ this indicates that the
center cell will have much larger oscillations near onset than the ring cells,
with a period of oscillations of $\approx 8$ and with a phase shift of
$\approx 40^{\circ}$ between the ring and center cell oscillations. From the
FlexPDE simulations of (1.0.2) when $\tau=3$ and $D=0.5$ shown in Fig. 12,
corresponding to the blue dot in the left panel of Fig. 11, we observe that
these predictions of the linearized theory are roughly satisfied. Moreover,
since the efflux rate out of the center cell is larger, while the influx rate
is smaller, the rather small bulk diffusivity $D=0.5$ should qualitatively
lead to a buildup of the bulk chemical near the center cell at certain times
in the oscillation. This feature is observed in the surface plot in the left
panel of Fig. 12, and provides a clear example of the diffusion-sensing
behavior, as regulated by the permeability parameters.

Figure 13: Numerical results for intracellular dynamics versus time computed
from the ODE system (3.0.16) of Proposition 2 corresponding to the PDE
simulations shown in Fig. 6 (first row), Fig. 7, (second row), Fig. 9 (third
row) and Fig. 12 (fourth row), respectively. In each case, although $D$ is not
large, results from the ODE system (3.0.16) are seen to compare surprisingly
well with the full FlexPDE simulations of the PDE-ODE system (1.0.2). First
row: $(D,\tau)=(1.0,0.55)$ (red dot in Fig. 4(a)). Compare with Fig. 6. Second
row: $(D,\tau)=(1.0,1.1)$ (blue dot in Fig. 4(a)). Compare with Fig. 7. Third
row: $(D,\tau)=(0.4,0.6)$ (red dot in the left panel of Fig. 8). Compare with
Fig. 9. Fourth row: $(D,\tau)=(0.5,3)$ (blue dot in the left panel of Fig.
11). Compare with Fig. 12. Although there is a phase shift between the ODE and
full PDE results for intracellular oscillations due to different initial
conditions used, the ODE system (3.0.16) captures well the amplitude and
period of intracellular oscillations observed in the full PDE simulations.

Next, we examine whether the ODE system (3.0.16) of Proposition 2, derived
under the assumption of large bulk diffusivity $D={\mathcal{O}}(\nu^{-1})\gg
1$, can still be used to reliably approximate the intracellular dynamics
observed in the full FlexPDE simulation results of (1.0.2), performed for
${\mathcal{O}}(1)$ values of $D$, in Figs. 6, 7, 9 and 12. In Fig. 13 we show
that the numerical results computed from the ODE system (3.0.16) compare
surprisingly well with the full PDE simulations with respect to the amplitude
and period of intracellular oscillations (first, third and fourth rows of Fig.
13) and the prediction of a linearly stable steady-state (second row of Fig.
13). In using the ODE system (3.0.16) we calculated $D_{0}$ as $D_{0}=D\nu$,
where $\nu={-1/\log\varepsilon}$ with $\varepsilon=0.05$. We remark that the
beating-type oscillations observed in Fig. 10 for the very small value
$D=0.05$ are not captured by the ODE system (3.0.16). Moreover, we emphasize
that the simpler ODE system corresponding to the well-mixed limit $D\to\infty$
as given in (3.0.19), and which was used in [20] and [26] for studying quorum-
sensing behavior, does not reliably approximate the intracellular oscillations
for the values of bulk diffusivity given in Figs. 6, 7, 9, 9, and 12.

###  4.4 A defective center cell: different Sel’kov kinetics

Next, we study how the HB boundaries in the $(D,\tau)$ plane are altered by
varying a Sel’kov kinetic parameter of the center cell. In Fig. 14 we plot the
HB boundaries for a ring and center cell pattern of $m=3$ cells, where the
identical cells on the ring have parameters as in (4.3.1), while the Sel’kov
parameter $\alpha_{3}$ for the defective center cell is either
$\alpha_{3}=0.86$ (red curves), $\alpha_{3}=0.96$ (blue curves), or
$\alpha_{3}=0.9$ (same as in the left panel of Fig. 4). The HB boundaries in
the right panel of Fig. 14 show a zoom for smaller values of $D$ than the
figure in the left panel. In these figures, the dashed and the heavy solid
curves are for the in-phase modes computed from (4.1.20) with $(+)$ and $(-)$,
respectively, while the thin solid curve is for the anti-phase mode computed
from (4.1.19). We observe that the HB boundary for the anti-phase mode is
independent of the Sel’kov kinetic parameter $\alpha_{3}$ for the center cell.
This follows from the facts that the steady-state solution to (4.1.4) for the
ring cells and its Jacobian $J_{1}$ in (4.1.7b) do not depend on $\alpha_{3}$.
As a result, the anti-phase HB boundary, computed from (4.1.19), is
independent of $\alpha_{3}$.

Figure 14: HB boundaries in the $\tau$ versus $D$ plane for a ring and center
cell pattern of $m=3$ cells, where the Sel’kov kinetic parameter $\alpha_{3}$
for the center cell is varied. The ring cells are centered at $(\pm 5,0)$ and,
with the exception of the Sel’kov kinetic parameter for the center cell, all
three cells have parameters as in (4.3.1). Left panel: $\alpha_{3}=0.86$ (red
curves), $\alpha_{3}=0.9$ (black curves), $\alpha_{3}=0.96$ (blue curves).
Right panel: a zoomed-in version of the left panel for smaller values of $D$.
In both panels, the dashed and heavy solid curves are for the in-phase modes
computed from (4.1.20), while the thin solid curve is for the anti-phase mode
computed from (4.1.19). Each mode is linearly unstable in their respective
lobes, while linearly stable steady-state solutions exists outside the union
of the lobes. The anti-phase HB boundaries are independent of $\alpha_{3}$,
and so are plotted on top of each other.

From the left panel of Fig. 14 we observe that as $\alpha_{3}$ increases the
parameter regions in the $(D,\tau)$ plane where intracellular oscillations
occur decreases. In particular, the in-phase $(+)$ instability lobe that was
open for $\alpha_{3}=0.86$ becomes closed when $\alpha_{3}=0.96$, thereby
precluding the possibility of intracellular oscillations when $D$ is
sufficiently large. To interpret this result, we observe from the middle panel
of Fig. 2 that, as $\alpha_{3}$ increases, the center cell becomes less
activated, with the parameters drifting further from the HB boundary for the
uncoupled cell. As a result, it becomes more difficult to trigger in-phase
intracellular oscillations for a group of coupled cells as $\alpha_{3}$
increases. Overall, Fig. 14 does show that rather small increases or decreases
in a parameter value of the nonlinear Sel’kov kinetics for one specific cell
can either extinguish or trigger intracellular oscillations for an entire
group of cells. The corresponding eigenvectors for selected points on the
bifurcation diagrams in Fig. 14 are shown in Table 2 for $\alpha_{3}=0.86$ and
$\alpha_{3}=0.96$.

|  | $\qquad\alpha_{3}=0.86$ | $\qquad\alpha_{3}=0.96$
---|---|---|---
mode | j | $(D,\tau)$ | $(\mbox{Re}c_{j},\mbox{Im}c_{j})$ | $\theta_{j}\,(\text{rad})$ | $(D,\tau)$ | $(\mbox{Re}c_{j},\mbox{Im}c_{j})$ | $\theta_{j}\,(\text{rad})$
| 1 |  | $(0.539,-0.080)$ | $6.14$ |  | $(0.639,0)$ | 0
In-phase ($+$) | 2 | $(0.233,0.744)$ | $(0.539,-0.080)$ | $6.14$ | $(0.475,0.672)$ | $(0.639,0)$ | 0
(heavy solid) | 3 |  | $(0.638,0)$ | $0$ |  | $(0.427,-0.0118)$ | $6.256$
| 1 |  | $-(0.461,0.0415)$ | $3.23$ |  | $(-0.254,0.0428)$ | $2.974$
In-phase ($-$) | 2 | $(0.108,0.284)$ | $-(0.461,0.0415)$ | $3.23$ | $(0.0492,0.168)$ | $(-0.254,0.0428)$ | $2.974$
(dashed curve) | 3 |  | $(0.756,0)$ | $0$ |  | $(0.932,0)$ | $0$
| 1 |  | $(0.707,0)$ | $0$ |  | $(0.707,0)$ | 0
Anti-phase | 2 | $(0.211,0.365)$ | $-(0.707,0)$ | $\pi$ | $(0.211,0.365)$ | $-(.707,0)$ | $\pi$
(thin solid) | 3 |  | $(0,0)$ | $0$ |  | $(0,0)$ | $0$
Table 2: Real and imaginary parts of the eigenvector $\boldsymbol{c}$ of the
GCEP matrix ${\mathcal{M}}(\lambda)$ in (4.1.6), computed for a few points on
the HB boundaries in Fig. 14. The cells on the ring are identical and located
at $(\pm 0.5,0)$ with parameters as given in (4.3.1). The center cell has the
same parameters, with the exception that it has a Sel’kov kinetic parameter
$\alpha_{3}$ different than those on the ring. Middle three columns:
$\alpha_{3}=0.86$. Last three columns: $\alpha_{3}=0.96$.

###  4.5 Diffusion-sensing: effect of cell locations

For a ring and center cell pattern, with cells centered at $(\pm r_{0},0)$ and
$(0,0)$, we study how the HB boundaries in the $(D,\tau)$ plane depend on the
ring radius $r_{0}$. In the first row of Fig. 15 we show the HB boundaries for
the in-phase and anti-phase modes when $r_{0}=0.25$, $r_{0}=0.5$ (same as Fig.
4(a)), and $r_{0}=0.75$ for the case where the cells are all identical with
permeabilities $d_{1}=0.8$ and $d_{2}=0.4$. Similar plots for the three values
of $r_{0}$ are shown in the second row of Fig. 15 for the case where the
center cell is now defective with $d_{13}=0.4$ and $d_{23}=0.2$ (Fig. 15(e)
for $r_{0}=0.5$ is the same as in Fig. 8).

When the cells are all identical, we observe from the top row in Fig. 15 that
the regions of instability for the anti-phase mode and the in-phase $(-)$ mode
shrink noticeably as $r_{0}$ decreases. In contrast, the instability region
for the in-phase $(+)$ mode is relatively insensitive to changes in $r_{0}$.
To qualitatively explain this observation, we examine the eigenvector of the
GCEP matrix given in Table 1 at points on the HB boundaries for the three
modes. From this table, for the in-phase $(-)$ mode only the two ring cells
oscillate in phase while the center cell has larger amplitude oscillations
that are roughly $180^{\circ}$ out of phase. For the anti-phase $(-)$ mode,
the center cell is quiescent while the ring cells oscillate $180^{\circ}$ out
of phase. Finally, for the dominant in-phase $(+)$ mode, the three identical
cells are synchronized with very similar amplitudes and phases. Intuitively,
as the ring radius $r_{0}$ decreases the cells become more clustered and, as a
result, intracellular dynamics can still be synchronized even for smaller
values of the bulk diffusivity $D$. As a result, when $r_{0}$ is small, the
anti-phase and in-phase $(-)$ instability lobes, in which the center and ring
cells are not synchronized, should exist only for very small values of $D$
(see Fig. 15(a)). When $D$ is very small, communication between cells that are
close can still be rather weak. Notice that for $r_{0}=0.75$, where the cells
are farther apart, the anti-phase lobe in Fig. 15(c) exists for larger values
of $D$ than when $r_{0}=0.25$ or $r_{0}=0.5$. For $r_{0}=0.75$, the ring cells
can oscillate, maintaining a quiescent center cell, provided that the bulk
chemical is not washed away, i.e. $D$ is not too large, owing to the fact that
each of the two ring cells are now relatively close to their images across the
domain boundary due to the reflecting boundary condition imposed. However,
when $r_{0}$ is large, the ring cells are far from each other and so, within
the anti-phase instability lobe, their oscillations are $180^{\circ}$ out of
phase.

(a) $r_{0}=0.25$

(b) $r_{0}=0.5$

(c) $r_{0}=0.75$

(d) $r_{0}=0.25$

(e) $r_{0}=0.5$

(f) $r_{0}=0.75$

Figure 15: HB boundaries in the $\tau$ versus $D$ plane for a ring and center
hole configuration of $m=3$ cells, for three different ring radii $r_{0}$ as
indicated. The Sel’kov parameters and cell radius are as in (4.3.1). Top row:
identical cells with permeabilities $d_{1}=0.8$ and $d_{2}=0.2$. Bottom row:
center cell is defective with $d_{13}=0.4$ and $d_{23}=0.2$. The dashed and
heavy solid curves are the HB boundaries for the in-phase $(+)$ and $(-)$
modes, respectively, as computed from (4.1.20). The thin solid curve is the HB
boundary for the anti-phase mode computed from (4.1.19). Each mode is unstable
in their respective lobes, while linearly stable steady-state solutions occur
outside the union of the lobes.

The second row of Fig. 15 shows similar results for the case where the center
cell is defective, with permeabilities $d_{13}=0.4$ and $d_{23}=0.2$,
corresponding to a reduced rate of influx into the center cell. Similar to the
results for identical cells presented in the first row of Fig. 15, as $r_{0}$
decreases the region of instability for the in-phase $(-)$ and the anti-phase
modes shrink. Moreover, we observe that the maximum extent in $D$ of the in-
phase $(+)$ lobes, in which all the cells are essentially synchronized in
amplitude and phase, decreases as $r_{0}$ increases. This is because when the
cells are farther apart, a relatively smaller value of the bulk diffusivity
can lead to a washing out of the bulk signal, which thereby weakens the
communication between the three cells and precludes synchronization.

##  5 Asymptotic analysis for a large population of cells

As the number of cells increases it becomes increasingly more challenging
numerically to determine the stability boundaries in parameter space from the
root-finding condition $\mbox{det}{\mathcal{M}}(\lambda)=0$ given an arbitrary
cell configuration $\boldsymbol{x}_{j}$ with arbitrary permeability parameters
$d_{1j}$ and $d_{2j}$ for $j=1,\ldots,m$. In this section we provide a simpler
approach to determine the stability boundaries from the GCEP (2.2.13) for the
large bulk diffusion parameter regime where $D={\mathcal{O}}(\nu^{-1})$. To
isolate only the effect of different cell configurations, such as shown in
Fig. 16, as well as the effect of different cell permeabilities, in our
analysis below we will assume that the Sel’kov kinetic parameters are all
cell-independent, i.e. that $\alpha_{j}=\alpha$, $\mu_{j}=\mu$ and
$\zeta_{j}=\zeta$ for $j=1,\ldots,m$. Our analysis is easily extended to
remove this assumption.

(a) Randomly placed cells

(b) Two clusters of cells

Figure 16: Schematic showing different cell configurations in the unit disk.
The cells are represented by smaller disks (in cyan) and the diffusing bulk
species corresponds to the red dots. Left panel: randomly placed cells. Right
panel: two spatially clustered groups of cells.

We first recall from (2.2.21b) of section 2.2 that, for $\varepsilon\to 0$,
the steady-state solution with Sel’kov kinetics is determined in terms of the
solution $\boldsymbol{\mathcal{A}}=(A_{1},\dots,A_{m})^{T}$ to the linear
system

$\left(I+2\pi\nu\mathcal{G}+\nu\,D\,P_{1}+\frac{2\pi\nu
D}{\tau}P_{2}\right){\mathcal{A}}=-\mu\nu P_{2}\,\boldsymbol{e}\,;\qquad
P_{1}=\mbox{diag}\Big{(}\frac{1}{d_{11}},\ldots,\frac{1}{d_{1m}}\Big{)}\quad
P_{2}=\mbox{diag}\Big{(}\frac{d_{21}}{d_{11}},\ldots,\frac{d_{2m}}{d_{1m}}\Big{)}\,,$
(5.0.1a) where $\boldsymbol{e}=(1,\dots,1)^{T}$ and $\mathcal{G}$ is the
Green’s interaction matrix defined in (2.1.7). In terms of the solution
$\boldsymbol{\mathcal{A}}$ to (5.0.1a), we obtain from (2.2.21a) that the
steady-state for the intracellular species
$\boldsymbol{u}_{e}=(u^{1}_{ej},u^{2}_{ej})^{T}$ in each cell is
$u^{1}_{ej}=\mu+\frac{2\pi D}{\tau}A_{j}\,,\qquad\mbox{and}\qquad
u^{2}_{ej}=\frac{\mu}{\alpha+(u^{1}_{ej})^{2}}\,.$ (5.0.1b)

Since this simple steady-state construction is accurate to all orders in $\nu$
for any $D>0$, we will not seek an approximation to it valid for the regime
$D={\mathcal{O}}(\nu^{-1})$. Therefore, in the GCEP matrix
${\mathcal{M}}(\lambda)$ given in (2.2.13a), the diagonal matrix
${\mathcal{K}}(\lambda)$, as defined in (2.2.12a), will be evaluated at the
solution to (5.0.1). With Sel’kov kinetics, this yields

$\mathcal{K}=\mathcal{K}(\lambda)\equiv\mbox{diag}\Big{(}\mathit{K}_{1},\dots,\mathit{K}_{m}\Big{)}\,,\qquad\mbox{where}\qquad\mathit{K}_{j}=\frac{\lambda+\mbox{det}(J_{j})}{\lambda^{2}-\lambda\,\mbox{tr}(J_{j})+\mbox{det}(J_{j})}\,,$
(5.0.2a) where $J_{j}$ is the Jacobian matrix for the local Sel-kov kinetics
at the $j^{\text{th}}$ cell, for which
$\mbox{det}(J_{j})=\zeta\left(\alpha+(u_{ej}^{1})^{2}\right)\,,\qquad\mbox{tr}(J_{j})=\frac{\left[2\mu
u_{ej}^{1}-\left(\alpha+(u_{ej}^{1})^{2}\right)-\zeta\left(\alpha+(u_{ej}^{1})^{2}\right)^{2}\right]}{\alpha+(u_{ej}^{1})^{2}}\,.$
(5.0.2b)

Our goal in this subsection of determining a more tractable formula for
determining the roots of $\det\mathcal{M}(\lambda)=0$ when $m$ is large and in
the limit $D={D_{0}/\nu}$, where $\nu={-1/\log\varepsilon}$ and
$D_{0}={\mathcal{O}}(1)$, is based on approximating the eigenvalue-dependent
Green’s matrix ${\mathcal{G}}_{\lambda}$ in (2.2.13a), while retaining the
full $\mathcal{K}$ matrix, as defined in (5.0.2).

To do so, we readily calculate from (2.2.8) that when $D={D_{0}/\nu}\gg 1$ we
have (see [20])

$\mathcal{G}_{\lambda}=\frac{D_{0}}{\nu|\Omega|(1+\tau\lambda)}+\mathcal{G}_{0}+\mathcal{O}(\nu)\,.$
(5.0.3)

Here $\mathcal{G}_{0}$ is the Neumann Green’s matrix of (3.0.17) given in
terms of the Neumann Green’s function $G_{0}$ of (3.0.14), which is known
analytically for the the unit disk in (3.0.21). By using (5.0.3) in (2.2.13a),
we find that the GCEP matrix reduces to

$\mathcal{M}(\lambda)={\mathcal{M}}_{\infty}(\lambda)+\mathcal{O}(\nu^{2})\,,\qquad\mbox{where}\qquad{\mathcal{M}}_{\infty}\equiv{\mathcal{B}}+\gamma
E+2\pi\nu{\mathcal{G}}_{0}\,.$ (5.0.4a) In (5.0.4a) the scalar $\gamma$, the
rank-one matrix $E$ and the diagonal matrix ${\mathcal{B}}$ are defined by
$\gamma=\gamma(\lambda)\equiv\frac{2\pi
mD_{0}}{(1+\tau\lambda)|\Omega|}\,,\qquad
E\equiv\frac{1}{m}\boldsymbol{e}\,\boldsymbol{e}^{T}\,,\qquad{\mathcal{B}}={\mathcal{B}}(\lambda)\equiv
I+D_{0}P_{1}+\frac{2\pi D_{0}}{\tau}P_{2}\,\mathcal{K}(\lambda)\,.$ (5.0.4b)
Neglecting terms of order ${\mathcal{O}}(\nu^{2})$, in the limit
$D={D_{0}/\nu}\gg 1$ we conclude that the GCEP (2.2.13) reduces to finding
values of $\lambda$ for which there is a nontrivial solution $\boldsymbol{c}$
to
$\mathcal{M}_{\infty}(\lambda)\boldsymbol{c}=\boldsymbol{0}\,,\qquad\mbox{which
occurs iff}\qquad\det\mathcal{M}_{\infty}(\lambda)=0\,.$ (5.0.4c)

We remark that in our asymptotic reduction (5.0.4) of the GCEP (2.2.13), the
cell configuration enters through the diagonal matrix ${\mathcal{K}}$, as
defined in (5.0.2), and at order $\nu$ from the Neumann Green’s matrix
${\mathcal{G}}_{0}$ in (5.0.4a).

To analyze (5.0.4) we first observe from (5.0.4b) and (5.0.2a) that

${\mathcal{B}}=\mbox{diag}(b_{1},\ldots,b_{m})\,,\qquad\mbox{where}\qquad
b_{j}=\frac{(d_{1j}+D_{0})}{d_{1j}}\left[1+\frac{\eta_{j}}{\tau}\mathit{K}_{j}\right]\,,\qquad\eta_{j}\equiv\frac{2\pi
d_{2j}D_{0}}{d_{1j}+D_{0}}\,,\quad j=1,\ldots,m\,.$ (5.0.5)

The key advantage of our asymptotic reduction is that by using the matrix
structure in (5.0.4) we can readily calculate a tractable analytical
expression for $\mbox{det}{{\mathcal{M}}_{\infty}(\lambda)}$, which can be
used in implementing the root-finding condition in (5.0.4c) numerically. Our
first result for (5.0.4) in this direction is as follows:

######  Proposition 4.

Suppose that $b_{j}\neq 0$ for $j=1\,,\ldots,m$ so that ${\mathcal{B}}$ is
invertible. Then,

$\mbox{det}{\mathcal{M}}_{\infty}(\lambda)=\mbox{det}\left({\mathcal{B}}+\gamma
E+2\pi\nu{\mathcal{G}}_{0}\right)=\left(\prod_{i=1}^{m}b_{i}\right)\kappa_{1\varepsilon}\left(1+{\mathcal{O}}(\nu)\right)\,,$
(5.0.6)

where

$\kappa_{1\varepsilon}\equiv
1+\frac{\gamma}{m}\boldsymbol{e}^{T}\boldsymbol{v}_{1}+2\pi\nu\,\left(\frac{\boldsymbol{v}^{T}_{1}\,\mathcal{G}_{0}\,\boldsymbol{v}_{1}}{\boldsymbol{e}^{T}\,\boldsymbol{v}_{1}}\right)+\mathcal{O}(\nu^{2})\,,\qquad\mbox{with}\qquad\boldsymbol{v}_{1}^{T}\equiv\left({1/b_{1}},\ldots,{1/b_{m}}\right)\,.$
(5.0.7)

###### Proof.

Since ${\mathcal{B}}$ is a diagonal matrix and invertible by assumption, we
have

$\mbox{det}\left({\mathcal{B}}+\gamma
E+2\pi\nu{\mathcal{G}}_{0}\right)=\mbox{det}({\mathcal{B}})\det\left(I+\gamma{\mathcal{B}}^{-1}E+2\pi\nu{\mathcal{B}}^{-1}{\mathcal{G}}_{0}\right)\,,$
(5.0.8)

where $\mbox{det}{\mathcal{B}}=\prod_{i=1}^{m}b_{i}$. The second term in
(5.0.8) is calculated by multiplying all the eigenvalues
$\kappa_{\varepsilon}$ of

$\left(I+\gamma{\mathcal{B}}^{-1}E+2\pi\nu{\mathcal{B}}^{-1}{\mathcal{G}}_{0}\right)\,\boldsymbol{v}_{\varepsilon}=\kappa_{\varepsilon}\boldsymbol{v}_{\varepsilon}\,.$
(5.0.9)

To determine ${\mathcal{O}}(\nu)$ accurate expressions for these eigenvalues
we expand an eigenpair $\kappa_{\varepsilon}$ and
$\boldsymbol{v}_{\varepsilon}$ as

$\kappa_{\varepsilon}=\kappa_{1}+\nu\,\tilde{\kappa}+\cdots\,,\qquad\boldsymbol{v}_{\varepsilon}=\boldsymbol{v}_{1}+\nu\,\tilde{\boldsymbol{v}}+\cdots\,.$
(5.0.10)

Upon substituting (5.0.10) into (5.0.9) we equate powers of $\nu$ to obtain
the leading-order problem

$(I+\gamma{\mathcal{B}}^{-1}E)\,\boldsymbol{v}_{1}=\kappa_{1}\,\boldsymbol{v}_{1}\,,$
(5.0.11)

and the following problem at $\mathcal{O}(\nu)$:

$(I+\gamma{\mathcal{B}}^{-1}E-\kappa_{1}I)\,\tilde{\boldsymbol{v}}=\tilde{\kappa}\,\boldsymbol{v}_{1}-2\pi{\mathcal{B}}^{-1}{\mathcal{G}}_{0}\,\boldsymbol{v}_{1}\,.$
(5.0.12)

For the leading order problem (5.0.11), we first observe that
$E\,\boldsymbol{q}_{j}=\boldsymbol{0}$ for $j=2,\dots,m$, where
$\mbox{span}\\{{\boldsymbol{q}_{2},\ldots,\boldsymbol{q}_{m}\\}}$ is the $m-1$
dimensional subspace orthogonal to $\boldsymbol{e}$. We readily conclude that
(5.0.11) has the eigenpairs

$\kappa_{1j}=1\quad\mbox{and}\quad\boldsymbol{v}_{1j}=\boldsymbol{q}_{j}\,,\qquad\mbox{for}\quad
j=2,\ldots,m\,.$ (5.0.13)

To determine the ${\mathcal{O}}(\nu)$ correction to these eigenvalues, we
substitute (5.0.13) into (5.0.12) to obtain

$\gamma{\mathcal{B}}^{-1}E\,\tilde{\boldsymbol{v}}=\tilde{\kappa}\boldsymbol{q}_{j}-2\pi{\mathcal{B}}^{-1}{\mathcal{G}}_{0}\,\boldsymbol{q}_{j}.$
(5.0.14)

Since ${\mathcal{B}}\,\boldsymbol{q}_{j}$, with
${\mathcal{B}}^{T}={\mathcal{B}}$, is a left null-vector of
${\mathcal{B}}^{-1}E$, the solvability condition for (5.0.14) is
$\boldsymbol{q}_{j}^{T}{\mathcal{B}}\left(\tilde{\kappa}\boldsymbol{q}_{j}-2\pi{\mathcal{B}}^{-1}{\mathcal{G}}_{0}\,\boldsymbol{q}_{j}\right)=0$,
which determines $\tilde{\kappa}$ for each eigenpair. In this way, we obtain
that $m-1$ eigenvalues of (5.0.9) are

$\kappa_{j\varepsilon}=1+2\pi\nu\,\frac{\boldsymbol{q}_{j}^{T}\,\mathcal{G}_{0}\,\boldsymbol{q}_{j}}{\boldsymbol{q}_{j}^{T}\,{\mathcal{B}}\,\boldsymbol{q}_{j}}+\mathcal{O}(\nu^{2})\,,\qquad\mbox{where}\quad\boldsymbol{q}_{j}^{T}\boldsymbol{e}=0\,,\quad\mbox{for}\quad
j=2,\dots,m\,.$ (5.0.15)

To find the remaining eigenpair of the leading order problem we write (5.0.11)
as
$\boldsymbol{v}_{1}(1-\kappa_{1})=-\frac{\gamma}{m}\left(\boldsymbol{e}^{T}\boldsymbol{v}_{1}\right)\left({1/b_{1}},\ldots,{1/b_{m}}\right)^{T}$.
This has the (unnormalized) nontrivial solution

$\boldsymbol{v}_{1}=\left({1/b_{1}},\ldots,{1/b_{m}}\right)^{T}\,,\qquad\mbox{iff}\qquad\kappa_{1}=1+\frac{\gamma}{m}\boldsymbol{e}^{T}\boldsymbol{v}_{1}=1+\frac{\gamma}{m}\sum_{i=1}^{m}\frac{1}{b_{i}}\,.$
(5.0.16)

To determine the $\mathcal{O}(\nu)$ correction to this eigenvalue we observe,
for this eigenpair, that $\boldsymbol{e}$ is a left null-vector of the matrix
in (5.0.12) since, by using (5.0.16) for $\kappa_{1}$, we calculate

$\boldsymbol{e}^{T}\left(I+\gamma{\mathcal{B}}^{-1}E-\kappa_{1}I\right)=\boldsymbol{e}^{T}\left(1+\frac{\gamma}{m}\left(\sum_{i=1}^{m}\frac{1}{b_{i}}\right)-\kappa_{1}\right)=\boldsymbol{0}\,.$
(5.0.17)

Therefore, upon left-multiplying (5.0.12) by $\boldsymbol{e}^{T}$, the
solvability condition for (5.0.12) is
$\boldsymbol{e}^{T}\left(\tilde{\kappa}\boldsymbol{v}_{1}-2\pi{\mathcal{B}}^{-1}{\mathcal{G}}_{0}\boldsymbol{v}_{1}\right)=0$,
which yields