#  Localized Signaling Compartments in 2-D Coupled by a Bulk Diffusion Field:
Quorum Sensing and Synchronous Oscillations in the Well-Mixed Limit

Sarafa A. Iyaniwura  Michael J. Ward

###### Abstract

We analyze oscillatory instabilities for a coupled PDE-ODE system modeling the
communication between localized spatially segregated dynamically active
signaling compartments that are coupled through a passive extracellular bulk
diffusion field in a bounded 2-D domain. Each signaling compartment is assumed
to secrete a chemical into the extracellular medium (bulk region) and it can
also sense the concentration of this chemical in the region around its
boundary. This feedback from the bulk region, resulting from the entire
collection of cells, in turn modifies the intracellular dynamics within each
cell. In the limit where the signaling compartments are circular disks with a
small common radius $\varepsilon\ll 1$ and where the bulk diffusivity is
asymptotically large, a matched asymptotic analysis is used to reduce the
dimensionless PDE-ODE system into a nonlinear ODE system with global coupling.
For Sel’kov reaction kinetics, this ODE system for the intracellular dynamics
and the spatial average of the bulk diffusion field is then used to
investigate oscillatory instabilities in the dynamics of the cells that are
triggered due to the global coupling. In particular, numerical bifurcation
software on the ODEs is used to study the overall effect of coupling defective
cells (cells that behave differently from the remaining cells) to a group of
identical cells. Moreover, when the number of cells is large, the Kuramoto
order parameter is computed to predict the degree of phase synchronization of
the intracellular dynamics. Quorum sensing behavior, characterized by the
onset of collective behavior in the intracellular dynamics as the number of
cells increases above a threshold, is also studied. Our analysis shows that
the cell population density plays a dual role of triggering and then quenching
synchronous oscillations in the intracellular dynamics.

Key words: Hopf bifurcation, synchronous oscillations, quorum sensing, bulk
diffusion, defector cells, Kuramoto order parameter, global coupling.

##  1 Introduction

Cell-cell communication is an important component of microbiological systems,
which involves the sending and receiving of information from one cell to
another. In a biological system consisting of unicellular organisms, cells
coordinate and communicate with one another to accomplish tasks that cannot be
achieved by a single cell. Those cells that are not in close proximity can
only communicate through the extracellular space between them by both
secreting and then sensing signaling chemicals in the extracellular medium,
called autoinducers, and adjusting their intracellular activities accordingly.
This bulk mediated form of cell-to-cell communication has been observed in
many specific biological systems including, a collection of social amoebae
Dictyostelium discoideum where the secretion of cyclic adenosine monophosphate
(cAMP) by the cells triggers an oscillatory response and helps guide the cells
to aggregation (cf. [10], [21]), colonies of yeast cells in which the exchange
of acetaldehyde (Aca) molecules leads to glycolytic oscillations (cf. [4]),
and the emergence of bioluminescence for the marine bacterium Vibrio fischeri
at large cell densities (cf. [17]). In this context, the term quorum sensing
refers to a process of cell-to-cell communication, mediated by a bulk
diffusion field, that triggers some collective behavior of a group of cells as
the cell density increases.

Quorum sensing phenomenon where the collective dynamics involves a switch-like
emergence of synchronous oscillations as the cell density increases past a
threshold has been observed and studied in several specific biological systems
(cf. [4], [10], [12], [18], [14]) as well as in physicochemical systems
involving small catalyst-loaded particles immersed in a chemical mixture (cf.
[29], [28], [31], [30]). An overview of some universal features for quorum
sensing systems in biology are discussed in [12]. In other cell signalling
problems, spatio-temporal oscillations associated with chemically active, but
spatially localized, sites that are mediated by a bulk diffusion field also
arise in models of certain gene regulatory networks [2].

Figure 1: A schematic diagram for (1.0.1) showing the dynamically active
signaling compartments (in cyan) and the extracellular bulk region (unshaded)
in a bounded 2-D domain. Each circular cell contains two signaling chemicals
(red and green), but only the red chemical is secreted into the extracellular
bulk region. A zoomed-in illustration of the secretion and feedback of
chemicals into the cells is shown on the right.

Our goal herein is to use a theoretical framework based on a coupled PDE-ODE
system to study the onset of synchronous intracellular dynamics in a
collection of small signalling compartments or cells that are coupled through
a 2-D passive bulk diffusion field. We begin by outlining the formulation of
this coupled cell-bulk PDE-ODE model of [9], which was inspired by the models
introduced in [19] (see also [20] and [32]). Let $\Omega\subset\mathbb{R}^{2}$
be a bounded domain that contains $m$ disk-shaped dynamically active
signalling compartments or “cells” of a common radius $R$, denoted by
$\Omega_{j}$ and centered at $\boldsymbol{X}_{j}\in\Omega$ for $j=1,\dots,m$.
We assume that the concentration $\mathcal{U}(\boldsymbol{X},T)$ of the
signaling molecule or autoinducer in the bulk or extracellular region
$\Omega\setminus\cup_{j=1}^{m}\,\Omega_{j}$, which is confined within
$\partial\Omega$, satisfies

$\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_{1,j}\,\mathcal{U}-\beta_{2,j}\,\mu_{j}^{1}\,,\qquad\boldsymbol{X}\in\partial\Omega_{j}\,,\quad
j=1,2,\ldots,m\,.$ (1.0.1b) Here $D_{B}>0$ and $k_{B}>0$ are the bulk
diffusivity and the rate of degradation of the signaling chemical in the bulk
region, respectively. In the Robin boundary condition (1.0.1b) on each cell
boundary, $\beta_{1,j}>0$ and $\beta_{2,j}>0$ are dimensional constants
modeling the permeability of the boundary for the $j^{\text{th}}$ cell, while
$\partial_{n_{\boldsymbol{X}}}$ is the outward normal derivative pointing into
the bulk region. Inside each dynamically active cell, which is assumed to be
well-mixed, we suppose that there are $n$ interacting species, represented by
the vector-field
$\boldsymbol{\mu}_{j}\equiv(\mu_{j}^{1},\ldots,\mu_{j}^{n})^{T}$, that if
isolated from the bulk would interact according to the local reaction kinetics
$\boldsymbol{F}_{j}\left(\boldsymbol{\mu}_{j}/\mu_{c}\right)$. However, owing
to the transport across each cell boundary, the bulk diffusion field of
(1.0.1) is coupled to the local kinetics within the $j^{\text{th}}$ cell 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_{1,j}\,\mathcal{U}-\beta_{2,j}\,\mu_{j}^{1}\right)\,\text{d}S_{\boldsymbol{X}}\,,\qquad
j=1,\ldots,m\,.\end{split}$ (1.0.1c)

Here $\boldsymbol{e}_{1}\equiv(1,0,\ldots,0)^{T}$, $k_{R}>0$ is the typical
reaction rate for the dimensionless local kinetics $\boldsymbol{F}_{j}$, and
$\mu_{c}>0$ is a typical scalar measure of $\boldsymbol{\mu}_{j}$. In this
model, each cell secretes only one signaling chemical, labeled by
$\mu_{j}^{1}$, across its boundary into the bulk medium, and this secretion is
regulated by the efflux permeability parameter $\beta_{2,j}$. However, each
cell receives feedback across its boundary, as controlled by the influx
permeability parameter $\beta_{1,j}$, from the entire collection of cells as
mediated by the bulk concentration field $\mathcal{U}(\boldsymbol{X},T)$ for
$\boldsymbol{X}\in\partial\Omega_{j}$. Figure 1 shows a schematic diagram to
illustrate the cell-bulk coupling in (1.0.1) for the case where there are
$n=2$ intracellular species.

The coupled PDE-ODE model (1.0.1), non-dimensionalized in §2.1 and Appendix A
and written in dimensionless form below in (2.1.1), is analyzed in the limit
of a small common cell radius when the bulk diffusivity is sufficiently large
so that the bulk medium can be considered as “well-mixed”. In this limit,
where the leading-order bulk concentration field is spatially homogeneous, a
matched asymptotic expansion analysis is used in §2.2 to reduce the
dimensionless coupled PDE-ODE system (2.1.1) into an $nm+1$ component
nonlinear ODE system, with global coupling, which is given below in (2.2.18).
A similar reduction, but only for the case of identical cells, was given in §5
of [9]. We remark that our asymptotic reduction, starting from a governing
PDE-ODE coupled system, provides a systematic method for deriving limiting
globally coupled ODE systems characterizing cell-cell communication through a
well-mixed bulk medium. This is in contrast to the more phenemologically-based
ODE models used in previous studies (cf. [4], [22], [25], [18], [7]). Our
limiting ODE system is related to those used in [16] and [15] to model
communication between nonlinear oscillators coupled indirectly through an
external medium.

For the specific case of Sel’kov reaction kinetics [26], the ODE system
(2.2.18) for the intracellular dynamics and the spatial average of the bulk
diffusion field is then used as a conceptual model to investigate the
emergence of oscillatory instabilities in the intracellular dynamics that are
triggered due to the global coupling. In the absence of coupling to the bulk
medium, the Sel’kov kinetic parameters are chosen so that each cell is a
conditional oscillator, in that its dynamics has only a stable steady-state
when uncoupled from the bulk. As a result, in our study the emergence of
intracellular oscillations is inherently due to the cell-cell interaction, as
mediated by the global coupling through the well-mixed bulk region. However,
in qualitative agreement with experimental findings in [4] and [3] involving
yeast cells that are typically non-oscillatory when isolated, the Sel’kov
parameters are chosen relatively close to the threshold for the onset of
limit-cycle oscillations for an isolated cell.

With Sel’kov reaction kinetics, in §3 we compute global branches of steady-
states and periodic solutions for the ODEs (2.2.18) using the numerical
bifurcation software XPPAUT [6] for a small collection of cells. Our goal is
to show how the permeability parameters, the Sel’kov kinetic parameters, and
the number of cells influence the emergence of intracellular oscillations from
a quiescent linearly stable steady-state. In our numerical bifurcation study
in §3 we consider two main scenarios: one where all the cells are identical
with common parameters, and the other a single defective cell, which has
either different permeabilities or kinetic parameters, is coupled to a group
of identical cells. Our numerical results, presented in the form of one- and
two-parameter global bifurcation diagrams, show that for a collection of
identical cells synchronous intracellular oscillations can occur in various
parameter regimes as a result of the global coupling, even when their
individual dynamics are linearly stable. Moreover, we show that small changes
in the permeabilities or the reaction kinetics from a single defective cell
can either extinguish or trigger intracellular oscillations of the entire
group of cells. However, most typically, the range of these parameters that
yields stable periodic solutions are smaller for the cases involving a
defective cell, and this range decreases as the number of identical cells
increases. Our theoretical finding that a single defective cell can trigger
oscillations in the entire group is qualitatively similar to the observations
in [10], based on live-cell imaging of the social bacteria Dictyostelium
discoideum, that stochastic pulses originating from a discrete location play a
critical role in the onset of collective oscillations. Our study of the effect
of a defective cell on the entire group of cells is also related to game-
theoretic concepts, which have recently been applied to quorum sensing
behavior in microbial communities [1].

In §4 we recast the ODE system (2.2.18) in terms of a cell density parameter,
defined by $\rho\equiv{m/|\Omega|}$, in order to study the quorum sensing
behavior associated with a large collection of either identical cells or a
mixture of identical and defective cells. We show that the quorum sensing
transition is characterized by the switch-like emergence of synchronous
intracellular oscillations as $\rho$ crosses through a critical threshold.
This is followed by a switch-like quenching of these oscillations as the cell
density parameter crosses through a second, larger, threshold value. When the
cells are identical, we show that these two threshold values are Hopf
bifurcation points for a cubic polynomial. For a mixture of identical and
defective cells, where a Sel’kov kinetic parameter is randomly chosen from
some range, the quorum sensing and quenching thresholds are identified from a
numerical computation of the Kuramoto order parameter, defined in (4.2.1) (see
[15], [16]). The interval of the cell density parameter where the
synchronization of intracellular dynamics occurs was found, as expected, to
decrease as the number of defective cells increases.

In §5 we summarize some of our main results, and we discuss possible
extensions of this study for investigating some specific biologically-based
cell-cell communication problems that are mediated by a bulk diffusion field.
Related problems that fit within the theoretical framework offered by the
coupled PDE-ODE system (1.0.1) are also discussed.

##  2 The coupled PDE-ODE model: Non-dimensionalization and the well-mixed
limit

In this section we non-dimensionalize the coupled PDE-ODE model (1.0.1) and
study the dimensionless model in the large bulk diffusivity limit in a bounded
2-D domain. In this limit, the bulk region becomes well-mixed and the bulk
concentration field becomes spatially homogeneous. In this well-mixed limit, a
matched asymptotic analysis is used to reduce the dimensionless coupled PDE-
ODE model into a system of ODEs with global coupling. This limiting system is
used to investigate oscillatory instabilities in the intracellular dynamics
and quorum sensing behavior.

###  2.1 Non-dimensionalization

We assume that the signaling compartments are circular with equal radius
$R_{0}$, which is small relative to the length-scale $L$ of the domain. As
such, we introduce a small scaling parameter $\varepsilon\equiv R_{0}/L\ll 1$.
With this assumption, in Appendix A the coupled PDE-ODE model (1.0.1) is non-
dimensionalized such that the dimensionless spatio-temporal bulk concentration
field $U(\boldsymbol{x},t)$ 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}}\,;$
(2.1.1a)
$\displaystyle\partial_{n}\,U=\,0,\quad\boldsymbol{x}\in\partial\Omega\,,$
$\displaystyle\qquad\varepsilon
D\,\partial_{n_{j}}U=d_{1,j}\,U-d_{2,j}\,u_{j}^{1}\,,\quad\boldsymbol{x}\in\partial\Omega_{\varepsilon_{j}}\,,\quad
j=1,\ldots,m\,,$ (2.1.1b) where $\Omega_{\varepsilon_{j}}$ is a circular
region of radius $\varepsilon\ll 1$ centered at $\boldsymbol{x}_{j}\in\Omega$,
representing the $j^{\text{th}}$ cell, $D$ is the effective diffusion
coefficient in the bulk region, $\partial_{n}$ is the normal derivative
pointing outward on the boundary of the domain $\Omega$, and
$\partial_{n_{j}}$ is the normal derivative on the boundary of the
$j^{\text{th}}$ cell pointing out of the cell into the bulk medium. It is
assumed that the cells are well-separated from each other in the sense that
dist$(\boldsymbol{x}_{i},\boldsymbol{x}_{j})=\mathcal{O}(1)$ for $i\neq j$,
and that the center of each cell is at an $\mathcal{O}(1)$ distance from the
boundary of the domain $\Omega$, i.e.
dist$(\boldsymbol{x}_{j},\partial\Omega)=\mathcal{O}(1)$ for $j=1,\dots,m$ as
$\varepsilon\to 0$. The dimensionless bulk concentration $U(\boldsymbol{x},t)$
is coupled to the local kinetics within the $j^{\text{th}}$ cell by
$\begin{split}\frac{d\boldsymbol{u}_{j}}{dt}&=\,\boldsymbol{F}_{j}\left(\boldsymbol{u}_{j}\right)+\frac{\boldsymbol{e}_{1}}{\varepsilon\tau}\int_{\partial\Omega_{\varepsilon_{j}}}\,(d_{1,j}\,U-d_{2,j}\,u_{j}^{1})\,\,\text{d}s\,,\quad
j=1,\ldots,m\,,\end{split}$ (2.1.1c)

where $\boldsymbol{u}_{j}\equiv(u_{j}^{1},\dots,u_{j}^{n})^{T}$ is a
dimensionless vector representing the $n$ chemical species in the
$j^{\text{th}}$ cell, $\boldsymbol{e}_{1}\equiv(1,0,\dots,0)^{T}$, while
$d_{1,j}>0$ and $d_{2,j}>0$ are dimensionless permeability parameters
modelling the influx and efflux on the $j^{\text{th}}$ cell boundary,
respectively. The intracellular, or local, kinetics within the $j^{\text{th}}$
cell is specified by the vector-field
$\boldsymbol{F}_{j}\left(\boldsymbol{u}_{j}\right)$. In (2.1.1), the
dimensionless parameters are defined in terms of the parameters of (1.0.1) by

$\begin{split}\tau\equiv\frac{k_{R}}{k_{B}}\,,\qquad
D\equiv\frac{D_{B}}{k_{B}L^{2}}\,,\qquad
d_{1,j}\equiv\frac{\varepsilon\beta_{1,j}}{k_{B}L}\,,\qquad
d_{2,j}\equiv\frac{\varepsilon\beta_{2,j}L}{k_{B}}\,.\end{split}$ (2.1.2)

With this non-dimensionalization, $\tau$ is large if the intracellular
reactions occur quickly with respect to the time it takes the autoinducer
signal to decay in the bulk medium. The dimensionless permeabilities $d_{1,j}$
and $d_{2,j}$ are taken to be $\mathcal{O}(1)$, which implies that
$\beta_{1,j}$ and $\beta_{2,j}$ are $\mathcal{O}(\varepsilon^{-1})$ since the
length-scale $L$ of the domain and the rate $k_{B}$ of chemical decay in the
bulk region are $\mathcal{O}(1)$. In our analysis below, this scaling of the
permeabilities is necessary in order to have significant secretion and
feedback of chemicals into each signaling compartment in the limit
$\varepsilon\to 0$.

###  2.2 Asymptotic analysis in the well-mixed limit

In this subsection, we study the dimensionless coupled PDE-ODE model (2.1.1)
in the well-mixed limit $D\gg\mathcal{O}(\nu^{-1})\gg\mathcal{O}(1)$, where
$\nu\equiv-1/\log\varepsilon$. In this limit, the method of matched asymptotic
expansions is used to reduce the coupled model (2.1.1) into a system of ODEs
with global coupling.

For $D\gg\mathcal{O}(\nu^{-1})$, we expand $U(\boldsymbol{x},t)$ in the
(outer) bulk region at ${\mathcal{O}}(1)$ distances from the cells as

$\begin{split}U(\boldsymbol{x},t)=U_{0}(\boldsymbol{x},t)+\frac{1}{D}\,U_{1}(\boldsymbol{x},t)+\cdots\end{split}$
(2.2.1)

Upon substituting (2.2.1) into (2.1.1), and collecting terms in powers of $D$,
we obtain the leading-order problem

$\begin{split}\Delta
U_{0}=0\,,\quad&\boldsymbol{x}\in\Omega\setminus\\{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{m}\\}\,;\qquad\partial_{n}U_{0}=0,\quad\boldsymbol{x}\in\partial\Omega\,,\end{split}$
(2.2.2)

which has the solution $U_{0}\equiv U_{0}(t)$. The next order problem for
$U_{1}$ in the outer bulk region is

$\begin{split}\Delta U_{1}=U_{0}+\tau
U_{0t}\,,\quad&\boldsymbol{x}\in\Omega\setminus\\{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{m}\\}\,;\qquad\partial_{n}U_{1}=0\,,\quad\boldsymbol{x}\in\partial\Omega\,,\end{split}$
(2.2.3)

which must be supplemented by singularity conditions for $U_{1}$ as
$\boldsymbol{x}\to\boldsymbol{x}_{j}$, for each $j=1,\ldots,m$.

To determine these conditions, we consider the inner region defined within an
$\mathcal{O}(\varepsilon)$ neighborhood of the $j^{\text{th}}$ cell. In this
region, we introduce the inner variables
$\boldsymbol{y}=\varepsilon^{-1}(\boldsymbol{x}-\boldsymbol{x}_{j})$ and
$U(\boldsymbol{x},t)=\mathcal{U}(\boldsymbol{x}_{j}+\varepsilon\boldsymbol{y},t)$,
with $\rho=|\boldsymbol{y}|$. Upon writing (2.1.1a) and (2.1.1b) in terms of
the inner variables, we obtain

$\displaystyle\tau\,\mathcal{U}_{t}$
$\displaystyle=\frac{D}{\varepsilon^{2}}\,\Delta_{\boldsymbol{y}}\,\mathcal{U}-\,\mathcal{U}\,,\quad
1<\rho<\infty\,;\qquad
D\,\partial_{\rho}\,\mathcal{U}=d_{1,j}\,\mathcal{U}-d_{2,j}\,u_{j}^{1}\,,\,\,\,\text{on}\,\,\,\rho=1\,,\qquad
j=1,\ldots,m\,.$ (2.2.4)

Since $D\gg\mathcal{O}(1)$, we expand the inner solution $\mathcal{U}(\rho,t)$
in powers of $D$ as

$\begin{split}\mathcal{U}(\rho,t)=\mathcal{U}_{0}(\rho,t)+\frac{1}{D}\,\mathcal{U}_{1}(\rho,t)+\cdots\,.\end{split}$
(2.2.5)

Upon substituting (2.2.5) into (2.2.4), and collecting terms in powers of $D$,
we derive the leading-order inner problem

$\begin{split}\Delta_{\rho}\,\mathcal{U}_{0}\equiv\frac{1}{\rho}\frac{\partial}{\partial\rho}\left(\rho\frac{\partial{\mathcal{U}}_{0}}{\partial\rho}\right)=0\,,&\quad
1<\rho<\infty\,;\qquad\partial_{\rho}\,\mathcal{U}_{0}=0\,,\quad\text{on}\quad\rho=1\,,\end{split}$
(2.2.6)

near the $j^{\text{th}}$ cell. The solution to this problem which matches the
outer solution as $\rho\to\infty$ is simply $\mathcal{U}_{0}\equiv U_{0}(t)$.

We proceed to consider the next order problem in the inner region given by

$\begin{split}\Delta_{\rho}\,\mathcal{U}_{1}&=0\,,\quad 1<\rho<\infty\,;\qquad
D\,\partial_{\rho}\,\mathcal{U}_{1}=\,d_{1,j}\,\mathcal{U}_{0}-d_{2,j}\,u_{j}^{1}\,,\quad\text{on}\quad\rho=1\,,\qquad
j=1,\ldots,m\,.\end{split}$ (2.2.7)

Allowing logarithmic growth as $\rho\to\infty$, the radially symmetric
solution to this problem is

$\begin{split}\mathcal{U}_{1,j}=(d_{1,j}\,\mathcal{U}_{0}-d_{2,j}\,u_{j}^{1})\log\rho+c_{j}\,,\qquad
j=1,\ldots,m\,,\end{split}$ (2.2.8)

where $c_{j}$ for $j=1,\ldots,m$ are unknown constants to be determined. We
write this solution in terms of the outer variables and substitute the
resulting expression into the inner expansion (2.2.5). For $\rho\to\infty$,
this yields that

$\begin{split}\mathcal{U}\sim
U_{0}+\frac{1}{D}\left[(d_{1,j}\,U_{0}-d_{2,j}u_{j}^{1})\,\log|\boldsymbol{x}-\boldsymbol{x}_{j}|+\frac{1}{\nu}(d_{1,j}\,U_{0}-d_{2,j}u_{j}^{1})\,+c_{j}\right]+\cdots,\quad
j=1,\ldots,m\,.\end{split}$ (2.2.9)

Since $D\gg\mathcal{O}(\nu^{-1})$, where $\nu={-1/\log\varepsilon}$, we have
$1/{D\nu}\ll 1$, which implies that $U_{0}$ is the largest term in the far-
field behavior (2.2.9) of the inner solution. This estimate establishes the
range of $D$ for which our analysis is valid. Upon matching the far-field
behaviour of the inner solution given in (2.2.9) to the outer solution, we
obtain the singularity behaviour of the outer solution $U_{1}$ as
$\boldsymbol{x}\to\boldsymbol{x}_{j}$. In this way, the complete outer problem
for $U_{1}$ is given by

$\displaystyle\Delta U_{1}$ $\displaystyle=U_{0}+\tau
U_{0t}\,,\quad\boldsymbol{x}\in\Omega\setminus\\{\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{m}\\}\,;\qquad\partial_{n}U_{1}=0\,,\quad\boldsymbol{x}\in\partial\Omega\,;$
(2.2.10a) $\displaystyle\vspace{0.5cm}U_{1}$
$\displaystyle\sim(d_{1,j}\,U_{0}-d_{2,j}\,u_{j}^{1})\,\log|\boldsymbol{x}-\boldsymbol{x}_{j}|\,,\quad\text{as}\quad\boldsymbol{x}\to\boldsymbol{x}_{j}\,,\quad
j=1,\ldots,m\,.$ (2.2.10b)

By representing the singularity behaviour for $U_{1}$ given in (2.2.10b) in
terms of a delta function for the PDE (2.2.10a), the problem for $U_{1}$ is
equivalently written as

$\begin{split}\Delta U_{1}&=(\tau
U_{0t}+U_{0})+2\pi\sum_{j=1}^{m}(d_{1,j}\,U_{0}-d_{2,j}\,u_{j}^{1})\,\delta(\boldsymbol{x}-\boldsymbol{x}_{j})\,,\quad\boldsymbol{x}\in\Omega\,;\qquad\partial_{n}U_{1}=0\,,\quad\boldsymbol{x}\in\partial\Omega\,.\end{split}$
(2.2.11)

Upon integrating (2.2.11) over the domain $\Omega$ and then applying the
divergence theorem, we obtain

$\begin{split}U_{0}^{\prime}=-\frac{1}{\tau}U_{0}-\frac{2\pi}{\tau|\Omega|}\sum_{j=1}^{m}(d_{1,j}\,U_{0}-d_{2,j}u_{j}^{1})\,,\end{split}$
(2.2.12)

where $|\Omega|$ is the area of $\Omega$. The chemical concentration in the
bulk region is described at leading-order by this ODE, where the first term on
the right hand side of the equation represents the decay of chemical in the
bulk region, and the summation averages the net flux of chemical through the
boundaries of all the signaling compartments.

To represent the solution to (2.2.11) we introduce the unique Neumann Green’s
function $G(\boldsymbol{x};\boldsymbol{x}_{j})$, which satisfies

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

where $R_{j}\equiv R(\boldsymbol{x}_{j})$ is its regular part at
$\boldsymbol{x}=\boldsymbol{x}_{j}$ for $j=1,\ldots,m$. Without loss of
generality, we impose $\int_{\Omega}U_{1}\,\text{d}\boldsymbol{x}=0$, so that
the spatial average of $U$ in the bulk region to order ${\mathcal{O}}(\nu)$ is
simply $U_{0}$. Then, the solution to (2.2.11) is

$U_{1}=-2\pi\sum_{i=1}^{m}\,(d_{1,i}\,U_{0}-d_{2,i}u_{i}^{1})\,G(\boldsymbol{x};\boldsymbol{x}_{i})\,.$
(2.2.14)

Upon substituting (2.2.14) into the outer expansion (2.2.1), we obtain the
two-term asymptotic expansion

$\begin{split}U\sim
U_{0}-\frac{2\pi}{D}\sum_{i=1}^{m}\,(d_{1,i}\,U_{0}-d_{2,i}\,u_{i}^{1})\,G(\boldsymbol{x};\boldsymbol{x}_{i})+\cdots\,,\end{split}$
(2.2.15)

for the outer solution in the bulk region. Then, we expand (2.2.14) as
$\boldsymbol{x}\to\boldsymbol{x}_{j}$ and use the singularity behaviour of the
Neumann Green’s function given in (2.2.13b) to determine the constants $c_{j}$
for $j=1,\ldots,m$ in (2.2.8). Upon matching the local behavior of (2.2.15) to
(2.2.9), we obtain that

$\begin{split}c_{j}=-\left(\frac{1}{\nu}+2\pi\,R_{j}\right)(d_{1,j}\,U_{0}-d_{2,j}u_{j}^{1})\,-2\pi\sum_{i\neq
j}^{m}\,(d_{1,i}\,U_{0}-d_{2,i}u_{i}^{1})\,G(\boldsymbol{x}_{j},\boldsymbol{x}_{i})\,,\qquad
j=1,\ldots,m\,.\end{split}$ (2.2.16)

In this way, (2.2.9) and (2.2.16) provide a two-term asymptotic expansion of
the inner solution $\mathcal{U}(\rho,t)$, which is valid within an
$\mathcal{O}(\varepsilon)$ neigbourhood of the $j^{\text{th}}$ cell. This
inner solution is given by

$\begin{split}\mathcal{U}=U_{0}+\frac{1}{D}\left[(d_{1,j}\,U_{0}-d_{2,j}u_{j}^{1})\,\log\rho+\frac{1}{\nu}(d_{1,j}\,U_{0}-d_{2}u_{j}^{1})\,+c_{j}\right]+\cdots\,,\qquad
j=1,\ldots,m\,.\end{split}$ (2.2.17)

Lastly, we substitute the leading-order inner solution $U_{0}(t)$ into
(2.1.1c) and evaluate the integral over the boundary of the $j^{\text{th}}$
cell. This provides the global coupling of the intracellular dynamics to the
bulk solution.

In summary, in the well-mixed limit $D\gg{\mathcal{O}}(\nu^{-1})$, where
$\nu\equiv{-1/\log\varepsilon}$ and $\varepsilon\ll 1$, the coupled PDE-ODE
model (2.1.1) reduces to the following $nm+1$ dimensional nonlinear ODE system
where the intracellular dynamics for $\boldsymbol{u}_{j}$ are globally coupled
through the spatially uniform bulk diffusion field $U_{0}(t)$:

$\frac{dU_{0}}{dt}=-\frac{U_{0}}{\tau}-\frac{2\pi}{\tau|\Omega|}\sum_{j=1}^{m}(d_{1,j}\,U_{0}-d_{2,j}u_{j}^{1})\,;\qquad\frac{d\boldsymbol{u}_{j}}{dt}=\,\boldsymbol{F}_{j}\left(\boldsymbol{u}_{j}\right)+\frac{2\pi\boldsymbol{e}_{1}}{\tau}\,(d_{1,j}\,U_{0}-d_{2,j}\,u_{j}^{1})\,,\quad
j=1,\ldots,m\,.$ (2.2.18)

Here $\boldsymbol{e}_{1}=(1,0,...,0)^{T}$, $m$ is the total number of
signaling compartments/cells, $|\Omega|$ is the area of $\Omega$,
$\boldsymbol{u}_{j}=(u_{j}^{1},\dots,u_{j}^{n})^{T}$ is an $n$-dimensional
vector whose entries are the $n$ chemical species in the $j^{\text{th}}$ cell,
and $\boldsymbol{F}_{j}(\boldsymbol{u}_{j})$ is the vector-field describing
the local reaction kinetics within the $j^{\text{th}}$ cell.

##  3 Global bifurcation diagrams

In this section, we compute global bifurcation diagrams for (2.2.18) using the
numerical bifurcation software XPPAUT [6] for the case where the intracellular
species interact according to Sel’kov reaction kinetics. Sel’kov kinetics
provides a two-component system of ODEs, which was originally used to
qualitatively model oscillations in the glycolytic pathway (cf. [26]). For
convenience of notation, we write the intracellular variable as
$\mathbf{u}=(v,w)^{T}$, and the local Sel’kov reaction kinetics
$\boldsymbol{F}(\mathbf{u})$ as $\boldsymbol{F}(v,w)=(F(v,w),G(v,w))^{T}$,
where

$\begin{split}F(v,w)\equiv\alpha w+wv^{2}-v\,,\qquad
G(v,w)\equiv\epsilon_{0}(\mu-(\alpha w+wv^{2}))\,.\end{split}$ (3.0.1)

Here $\alpha>0$, $\mu>0$ and $\epsilon_{0}>0$ are the Sel’kov reaction
parameters. The steady-state for Sel’kov dynamics ${dv/dt}=F(v,w)$ and
${dw/dt}=G(v,w)$ is $v_{e}=\mu$ and $w_{e}={\mu/(\alpha+\mu^{2})}$, while the
trace and determinant of the Jacobian matrix $J_{e}$ for this steady-state is
readily calculated as

$\det(J_{e})=\epsilon_{0}(\alpha+\mu^{2})\,,\qquad\text{tr}(J_{e})=2w_{e}v_{e}-1-\det(J_{e})=\frac{2\mu^{2}}{\alpha+\mu^{2}}-1-\epsilon_{0}(\alpha+\mu^{2})\,.$
(3.0.2)

Since $\det(J_{e})>0$, we conclude for the Sel’kov reaction kinetics that a
single isolated cell, decoupled from the bulk medium, has a steady-state that
is linearly stable if and only if $\text{tr}(J_{e})<0$. The Hopf bifurcation
boundary in the $\alpha$ versus $\mu$ parameter plane, as obtained by setting
$\text{tr}(J_{e})=0$ is given by

$\alpha=-\mu^{2}-\frac{1}{2\epsilon_{0}}+\frac{\sqrt{1+8\epsilon_{0}\mu^{2}}}{2\epsilon_{0}}\,.$
(3.0.3)

By constructing a trapping region and appealing to the Poincare-Bendixson
theorem, it is well-known that when the steady-state is unstable there is a
time-periodic solution, characterized by a limit cycle, for the isolated cell.
In contrast, when the steady-state is linearly stable there are no time-
periodic oscillations with Sel’kov kinetics.

Figure 2: The region of instability (green-shaded) in the $\alpha$ versus
$\mu$ parameter plane for the steady-state of a single isolated cell,
decoupled from the bulk region, that undergoes Sel’kov kinetics (3.0.1) with
$\epsilon_{0}=0.15$. Within this region, where $\text{tr}(J_{e})>0$ in
(3.0.2), there is a time-periodic solution (limit cycle) for an isolated cell.
On the boundary of this region, given by (3.0.3), is where Hopf bifurcations
occur. In the unshaded region the steady-state is linearly stable.

In our study of the ODE system (2.2.18), the Sel’kov kinetic parameters are
chosen so that the steady-state of the intracellular dynamics within each cell
is linearly stable, and hence no intracellular oscillations occur, when there
is no coupling to the bulk diffusion field. When $\epsilon_{0}=0.15$, this
corresponds to choosing a pair $(\mu,\alpha)$ not within the green-shaded
region in Figure 2. Our goal is to investigate intracellular oscillations that
are triggered as a result of the global coupling through the bulk medium. Our
results below are presented in the form of one- and two-parameter global
bifurcation diagrams computed from (2.2.18) using XPPAUT [6]. When the cells
are all identical, the bifurcation parameters are common to all of them.
However, for the cases involving a defective cell, the bifurcation parameter
is that of the defective cell only, which is identified as the cell with index
$j=1$.

####  3.0.1 Instability triggered by the rate of influx and efflux of the bulk
chemical

We begin by studying the effect of the permeability parameters $d_{1}$ and
$d_{2}$ on the intracellular dynamics. Recall that $d_{1,j}$ is the rate of
chemical feedback (influx) from the bulk into the $j^{\text{th}}$ cell, while
$d_{2,j}$ is the rate of secretion of chemical (efflux) by the cell into the
bulk region. Our goal below is to investigate synchronous oscillations in the
cells that are triggered by these parameters. One particular focus is to study
how the intracellular dynamics of a collection of identical cells is altered
due to changes in either the efflux or influx rate of one specific cell. The
numerical bifurcation diagrams shown below exhibit how the permeabilities of
such a single “defective” cell can significantly alter the parameter region
where collective intracellular oscillations of the entire group of cells can
occur.

In Figures 3(a) and 3(b) global branches of steady-state and periodic
solutions for (2.2.18) are shown as the permeability parameter $d_{1}$ for
$m=3$ and $m=8$ identical cells, respectively, is varied. In contrast, in
Figures 3(c) and 3(d) we plot similar bifurcation diagrams but for the case
where a single defective cell, having a modified permeability parameter
$d_{1,1}$, is coupled to the remaining identical cells of Figures 3(a) and
3(b). In both sets of figures the parameters used are $\alpha=0.9$,
$\epsilon_{0}=0.15$, $\mu=2$, $\tau=0.5$ and $d_{2}=0.2$. In Figures 3(c) and
3(d), $d_{1}=0.8$ is used for the identical cells. We observe from these
figures that the range of $d_{1}$ for which periodic solutions occur for
(2.2.18) is bounded for this parameter set. Since $d_{1}$ is the rate of
chemical feedback into the cells, a small value of $d_{1}$ implies little
feedback of the bulk chemical into the cells and, hence, a weak coupling
between the cells. As a result, since each cell has a linearly stable steady-
state and no time-periodic oscillations when uncoupled from the bulk, one
would expect no synchronous oscillations when $d_{1}$ is too small. This
intuition is confirmed in Figures 3(a)–3(b). Alternatively, when $d_{1}$ is
large relative to $d_{2}$ (rate of secretion of chemical into the
extracellular medium), the secreted chemical is almost immediately fed back
into the cells, which leads to a lower chemical concentration in the bulk
region. As a result, the coupling between the cells is again weak, and one
would expect that the steady-state is linearly stable and that no
intracellular oscillations occur. This intuition is again confirmed in Figures
3(a)–3(b). Later in our numerical study (see Figure 5 below), we use two-
parameter bifurcation diagrams to show that some specific ratios of $d_{1}$
and $d_{2}$ are required in order to trigger oscillatory instabilities of the
steady-state. For $m=3$ identical cells (see Figure 3(a)), the range of
$d_{1}$ that yields synchronous oscillations is $0.7613\leq d_{1}\leq 2.241$,
while for $m=8$ (see Figure 3(b)) that range is $0.2517\leq d_{1}\leq 0.8403$.
This indicates that the range of $d_{1}$ for which stable periodic solutions
occur shrinks as the number of cells increases, and as a result it becomes
increasingly difficult to synchronize the dynamics of the cells using the
permeability parameter $d_{1}$ when $d_{2}$ is fixed. Figure 3(c) shows the
result for two identical cells coupled to a defective cell, with intracellular
oscillations occurring for $0.6277\leq d_{1,1}\leq 1.134$. In Figure 3(d),
where seven identical cells are coupled to a defective cell, intracellular
oscillations occur only when $0.639\leq d_{1,1}\leq 0.9248$. Similar to the
case of identical cells, the range of $d_{1,1}$ for which linearly stable
periodic solutions exist is larger when $m=3$ as compared to when $m=8$. Most
importantly, however, Figures 3(c) and 3(d) clearly show how rather small
changes in the influx rate of a single defective cell, relative to the common
value of the group, can act as a switch that quenches, or extinguishes,
intracellular oscillations for the entire collection of cells.

(a) Three identical cells.

(b) Eight identical cells

(c) Two identical cells and a defective cell.

(d) Seven identical cells and a defective cell

Figure 3: Global bifurcation diagrams for $u_{1e}$ versus $d_{1}$, as computed
from the ODEs (2.2.18), showing steady-states and global branches of periodic
solutions. The thin red lines and bold black lines correspond to linearly
stable and unstable steady-state solution branches, respectively. The green
loops indicate linearly stable branches of periodic solutions. The parameters
are $\tau=0.5$, $\alpha=0.9$, $\epsilon_{0}=0.15$, and $\mu=2$, with
$d_{2,j}\equiv d_{2}=0.2$ for $j=1,\ldots,m$ and $d_{1,j}=0.8$ for
$j=2,\dots,m$. (a) and (b) are for $m=3$ and $m=8$ identical cells, with the
range of $d_{1}$ yielding synchronous oscillations is computed as $0.7613\leq
d_{1}\leq 2.241$ and $0.2517\leq d_{1}\leq 0.8403$, respectively. (c) is for
$m=3$ cells, where two of cells are identical with the remaining one
defective. The Hopf bifurcation points are $d_{1,1}=0.6277$ and
$d_{1,1}=1.134$. (d) is for $m=8$ cells (seven identical and a defective),
with Hopf bifurcations occurring at $d_{1,1}=0.639$ and $d_{1,1}=0.9248$. In
(c) and (d), the bifurcation parameter is that of the defective cell only.

The global bifurcation diagrams where the permeability parameter $d_{2}$ (the
efflux rate of secretion of chemicals into the bulk region by the cells) is
the bifurcation parameter are shown in Figure 4. The other parameter values
are the same as those in Figure 3. Similar to the results for $d_{1}$ in
Figure 3, the range of $d_{2}$ for which stable periodic solutions occur is
bounded above and below for all the scenarios considered. Linearly stable
steady-state solutions exist for small values of $d_{2}$ because the cells do
not secrete enough chemical required to trigger oscillations. On the other
hand, when the rate of secretion of chemicals is large relative to the rate of
feedback of chemicals into the cells, there is insufficient chemical
concentration within the cells to trigger oscillations. For $m=3$ and $m=8$
identical cells, we find that linearly stable synchronous periodic solutions
occur on the ranges $0.07624\leq d_{2}\leq 0.2283$ and $0.1908\leq d_{2}\leq
0.5229$, respectively. However, unlike the results in Figures 3(a) and 3(b),
the range of $d_{2}$ that predicts synchronous oscillations becomes larger as
the number of identical cells increases. For the case involving a defective
cell, when $m=3$ (two identical cells and one defective cell) the Hopf
bifurcation points are $d_{2,1}=0.1471$ and $d_{2,1}=0.2402$, while for $m=8$
(seven identical and one defective cell) the Hopf bifurcation points are
$d_{2,1}=0.1743$ and $d_{2,1}=0.2407$. Since the range of $d_{2,1}$ that
predicts stable periodic solutions is smaller when $m=8$ relative to when
$m=3$, this shows that when there is a defector cell it is more difficult to
achieve synchronous oscillations as the overall population of cells increases.
Comparing the identical cells cases (Figures 4(a) and 4(b)) with the cases
involving a defective cell (Figures 4(c) and 4(d)), we observe that the range
of $d_{2}$ that predicts linearly stable periodic solutions is larger, as
expected, when the cells are identical. This indicates that cell-cell
communication is more difficult to achieve when a single cell is secreting
chemical at a different rate as compared to the other cells. More
specifically, we observe that intracellular oscillations of the group can be
extinguished by a defective cell that reduces its efflux rate. From comparing
the results in Figures 3 and 4, we notice that the range of $d_{1}$ that
synchronizes the intracellular dynamics of the cells is larger than that of
$d_{2}$ for all the cases considered. This shows, for this parameter set, that
it is easier to trigger communication among the cells using the rate of
chemical feedback or influx from the bulk into the cells, rather than the rate
of secretion of chemicals into the bulk.

(a) Three identical cells

(b) Eight identical cells

(c) Two identical cells and a defector cell.

(d) Seven identical cells and a defector cell.

Figure 4: Global bifurcation diagrams for $u_{1e}$ versus $d_{2}$, as computed
from the ODEs (2.2.18), showing steady-states and global branches of periodic
solutions. The labeling of the curves is the same as in Figure 3. The
parameters are $\tau=0.5$, $\alpha=0.9$, $\epsilon_{0}=0.15$, and $\mu=2$,
with $d_{1,j}\equiv d_{1}=0.8$ for $j=1,\ldots,m$ and $d_{2,j}=0.2$ for
$j=2,\dots,m$. (a) and (b) are for $m=3$ and $m=8$ identical cells, and the
range of $d_{2}$ yielding synchronous oscillations is computed as $0.07624\leq
d_{2}\leq 0.2283$ and $0.1908\leq d_{2}\leq 0.5229$, respectively. (c) is for
$m=3$ cells, where two of cells are identical, and the remaining one is
defective. The Hopf bifurcation points are $d_{2,1}=0.1471$ and
$d_{2,1}=0.2402$. (d) is for $m=8$ cells (seven identical and one defective),
with Hopf bifurcations occurring at $d_{2,1}=0.1743$ and $d_{2,1}=0.2407$. In
(c) and (d), the bifurcation parameter is that of the defective cell only.

(a) Three identical cells

(b) Eight identical cells

(c) Two identical cells and a defective cell

(d) Seven identical cells and a defective cell

Figure 5: Two-parameter Hopf bifurcation diagrams in the plane of
permeabilities, as computed from the ODEs (2.2.18), showing regions of stable
periodic solutions (green) and regions of linearly stable steady-state
solutions (unshaded). The Sel’kov parameters are $\alpha=0.9$, $\mu=2$,
$\epsilon_{0}=0.15$, and we fix $\tau=0.5$. The results in (a) and (b) are for
$m=3$ and $m=8$ identical cells, respectively, having common permeability
parameters $d_{1}$ and $d_{2}$. (c) shows the result for $m=3$ cells (two
identical and a defective cell), while the result in (d) is for $m=8$ cells
(seven identical and a defective cell). For (c) and (d), the permeability
parameters $d_{1,1}$ and $d_{2,1}$ are for the defective cell only, while they
are fixed as $d_{1,j}=0.8$ and $d_{2,j}=0.2$ for $j=2,\dots,m$, for the
remaining identical cells.

Next, we present two-parameter global bifurcation diagrams in terms of the
permeability parameters $d_{1}$ and $d_{2}$. Figures 5(a) and 5(b) are for
$m=3$ and $m=8$ identical cells, respectively, while Figures 5(c) and 5(d) are
for the cases involving a defective cell. The remaining parameters are the
same as those used in Figures 3 and 4. In Figures 5(a)–5(d), stable periodic
solutions are predicted in the green wedge-shaped regions, while linearly
stable steady-state solutions occur in the unshaded regions. From Figure 5 the
region of intracellular oscillations are bounded by two straight black lines,
consisting of Hopf bifurcation values for the steady-state solution that
correspond to the one-parameter bifurcation diagrams in Figures 3 and 4. In
other words, the one-parameter bifurcation diagrams in Figures 3 and 4
represent to taking vertical and horizontal slices through the parameter plane
shown in Figures 5, respectively. We observe that as the number of identical
cells increases there is a larger region in the permeability parameter space
that lead to linearly stable synchronous intracellular oscillations (see
Figure 5(b)). Moreover, the permeability parameter space that yields
oscillations is relatively larger when the cells are all identical (Figure
5(a) and 5(b)), as compared to when a single defective cell has permeability
parameters that are different from those of the group (Figure 5(c) and 5(d)).
However, this is not the case for the scenario involving a defective cell. For
this case, as the number of remaining identical cells increases, it becomes
increasingly difficult to trigger intracellular oscillations of the group
through changes in the permeability parameters of a single defective cell (see
Figures 5(c) and 5(d)).

In Figure 14 of Appendix B we show a gallery of time-dependent solutions to
the ODE system (2.2.18) for permeability parameter values sampled from the
bifurcation diagrams in Figures 3–5.

####  3.0.2 Instability triggered by local reaction kinetics of the cells

In the previous subsection, we studied how intracellular oscillations are
triggered by modifying the rate of influx ($d_{1}$) and efflux ($d_{2}$) of
signaling chemicals across the cell boundaries. In this subsection, we use the
ODE system (2.2.18) to study intracellular oscillations that are triggered by
changes in the parameters of the local reaction kinetics. More specifically,
in (3.0.1), the Sel’kov parameters $\mu>0$ and $\alpha>0$ are used as
bifurcation parameters, while $\epsilon_{0}=0.15$ is fixed. In the results
below, the permeabilities were fixed at $d_{1,j}=0.8$ and $d_{2,j}=0.2$ for
$j=1,\dots,m$.

(a) Three identical cells

(b) Eight identical cells

(c) Two identical cells and a defective cell

(d) Seven identical cells and a defective cell

Figure 6: Global bifurcation diagrams for $u_{1e}$ versus $\alpha$, as
computed from the ODEs (2.2.18), showing steady-states and global branches of
periodic solutions. The thin red lines and bold black lines correspond to
linearly stable and unstable steady-state solution branches, respectively,
while the green curves are linearly stable branches of periodic solutions. The
parameters are $\tau=0.5$, $\epsilon_{0}=0.15$, $\mu=2$, with $d_{1,j}=0.8$
and $d_{2,j}=0.2$ for $j=1,\ldots,m$. (a) and (b) are for $m=3$ and $m=8$
identical cells, respectively. When $m=3$, the bifurcation point is
$\alpha_{c}=0.9182$, while $\alpha_{c}=0.9042$ for $m=8$. (c) is for two
identical cells and a defective cell, with a Hopf bifurcation occurring at
$\alpha_{c}=0.9562$, and (d) is for $m=8$, where seven of the cells are
identical and the remaining one is defective. The Hopf bifurcation occurs at
$\alpha_{c}=0.9342$. The bifurcation parameter in (c) and (d) is that of the
defective cell only, and $\alpha=0.9$ is fixed for the identical cells.

Figure 6 shows the stability of the steady-state solutions of the ODEs
(2.2.18) as $\alpha>0$ is varied. We observe from this figure that there is
only one Hopf bifurcation point, $\alpha_{c}$ for $\alpha>0$. For any value of
$\alpha\leq\alpha_{c}$, the ODE system has linearly stable periodic solutions
with increasing amplitude as $\alpha$ decreases. Linearly stable steady-state
solutions exist for $\alpha>\alpha_{c}$. Figures 6(a) and 6(b) are for $m=3$
and $m=8$ identical cells, respectively. The amplitudes of the predicted
oscillations for $m=3$ and $m=8$ are roughly similar, and the Hopf bifurcation
threshold for $m=3$ at $\alpha_{c}=0.9182$ is a little larger than for $m=8$
where $\alpha_{c}=0.9042$. Observe from Figure 2 that for a single isolated
cell, uncoupled from the bulk, a Hopf bifurcation occurs when $\alpha\approx
0.694$ as computed from (3.0.3). For the defective cell case (Figures 6(c) and
6(d)), when $m=3$ (two identical cells with $\alpha=0.9$ and a defective
cell), the Hopf bifurcation occurs at $\alpha_{c}=0.9562$, while
$\alpha_{c}=0.9342$ when $m=8$ cells (seven identical cells with $\alpha=0.9$
and a defective cell). As expected, similar to the case of identical cells,
the Hopf bifurcation value of $\alpha$ for the defective cell is smaller than,
and therefore closer to the stability boundary shown in Figure 2 (see equation
(3.0.3)), as the number of cells increases. Most importantly, since for
$\alpha=0.9$ there would be intracellular oscillations in the identical cells
when the defective cell is absent, we observe from Figures 6(c) and 6(d) that
when there is a defective cell then an increase in its kinetic parameter
$\alpha$ above $\alpha_{c}$ will extinguish the intracellular oscillations of
the entire group of cells. In contrast, from Figures 7(a) and Figures 7(b) we
observe that a defective cell can trigger collective oscillations in a group
of identical cells when the kinetic parameter of the defective cell decreases
below a threshold. In Figures 7(a) and Figures 7(b) we fix $\alpha=0.93$ so
that that no oscillations are possible for this group of identical cells.
However, by adding a defective cell, we observe that intracellular
oscillations for the entire group are triggered when the kinetic parameter for
the defector satisfies $\alpha<\alpha_{c}\approx 0.8953$ for $m=3$ and
$\alpha<\alpha_{c}\approx 0.7404$ for $m=8$. Since $\alpha_{c}>0.694$, the
defective cell is still a conditional oscillator, and it would have a linearly
stable steady-state when uncoupled from the bulk.

(a) Two identical cells and a defective cell

(b) Seven identical cells and a defective cell

Figure 7: Global bifurcation diagram for $u_{1e}$ versus the Sel’kov parameter
$\alpha$ for the defective cell. Same caption and parameter values as in
Figure 6 except that now $\alpha=0.93$ for the identical cells. When
$\alpha=0.93$, there are no intracellular oscillations for $m=3$ or $m=8$ (see
Figures 6(a) and 6(b)). As $\alpha$ decreases below the thresholds
$\alpha_{c}\approx 0.8953$ for $m=3$ and $\alpha_{c}\approx 0.7404$ for $m=8$,
the defector cell triggers intracellular oscillations for the entire group of
cells. On the range $0.694<\alpha<\alpha_{c}$, the defector cell still has a
linearly stable steady-state if uncoupled from the bulk.

The global bifurcation diagrams with respect to the Sel’kov parameter $\mu>0$
are shown in Figure 8. For an isolated cell, uncoupled from the bulk, a
horizontal slice at a sufficiently small value of $\alpha$ through the
instability region in Figure 2 shows that there are two Hopf bifurcation
values of $\mu$ for an isolated cell. The results in Figure 8 show how this
instability region of an isolated cell is modified due to coupling with the
bulk and in the presence of a defective cell. Similar to the results shown
above when $\alpha$ was the bifurcation parameter, the range of $\mu$ that
yields synchronous oscillations shrinks as the number of identical cells
increases (see Figures 8(a) and 8(b)). For $m=3$ identical cells, Hopf
bifurcations occur at $\mu_{1}=1.946$ and $\mu_{2}=2.987$, and when $m=8$ they
are at $\mu_{1}=1.845$ and $\mu_{2}=2.05$. For the case of two identical cells
(with $\mu=2$) and a defective cell, the Hopf bifurcation points are
$\mu_{1}=1.842$ and $\mu_{2}=3.312$ (Figure 8(c)). Figure 8(d) shows the
result for $m=8$ (seven identical cells with $\mu=2$ and a defective cell),
which has three Hopf bifurcation points. This result shows that linearly
stable periodic solutions exist for $0<\mu<0.3164$ and for $1.396<\mu<2.26$.
By comparing the results in Figures 8(c) and 8(d) we notice that the range of
$\mu$ that predicts linearly stable periodic solutions is smaller when $m=8$
than when $m=3$, even though synchronous oscillations are predicted in two
different intervals when $m=8$. This agrees with our earlier observation that
an increase in the population of identical cells makes synchronization more
difficult when they are coupled to a defective cell.

(a) Three identical cells

(b) Eight identical cells

(c) Two identical cells and a defective cell

(d) Seven identical cells and a defective cell

Figure 8: Global bifurcation diagrams for $u_{1e}$ versus $\mu$, computed from
the ODEs (2.2.18), showing steady-states and global branches of periodic
solutions. The labeling of the curves is the same as in Figure 6. The
parameters are $\tau=0.5$, $\epsilon_{0}=0.15$, $\alpha=0.9$, with
$d_{1,j}=0.8$ and $d_{2,j}=0.2$ for $j=1,\ldots,m$. (a) and (b) are for $m=3$
and $m=8$ identical cells respectively. The Hopf bifurcation points are
$\mu_{1}=1.946$ and $\mu_{2}=2.987$ for $m=3$, and $\mu_{1}=1.845$ and
$\mu_{2}=2.05$ for $m=8$. (c) is for two identical cells (with $\mu=2$) and a
defective cell, with Hopf bifurcation points at $\mu_{1}=1.842$ and
$\mu_{2}=3.312$. (d) is for $m=8$ cells, where seven of the cells are
identical (with $\mu=2$) and the remaining one is defective. There are three
Hopf bifurcation points for this case: $\mu_{1}=0.3164,\mu_{2}=1.396$ and
$\mu_{3}=2.26$. The bifurcation parameter in (c) and (d) is that of the
defective cell only.

(a) Three identical cells.

(b) Eight identical cells.

(c) Two identical cells and a defective cell.

(d) Seven identical cells and a defective cell.

Figure 9: Two-parameter bifurcation diagrams in the $\alpha$ versus $\mu$
parameter plane, computed from the ODEs (2.2.18), showing regions of linearly
stable periodic solutions (shaded in green) and linearly stable steady-state
solutions (unshaded). The black dashed curve in each figure is the Hopf
bifurcation boundary shown in Figure 2 for a single isolated cell, as given by
(3.0.3). The parameters are $\epsilon_{0}=0.15$, $\tau=0.15$ and with
$d_{1,j}=0.8$ and $d_{2,j}=0.2$ for $j=1,\dots,m$. (a) and (b) are for $m=3$
and $m=8$ identical cells, respectively, while (c) and (d) are for $m=3$ (two
identical cells and a defective cell) and $m=8$ (seven identical cells and a
defective cell), respectively. The bifurcation parameter in (c) and (d) is for
the defective cell only, with $\mu=2$ and $\alpha=0.9$ fixed for the identical
cells.

In Figure 9, we present two-parameter bifurcation diagrams for the Sel’kov
parameters $\alpha$ and $\mu$. In this figure, the green-shaded region is
where linearly stable periodic solutions occur, while the unshaded region
predicts linearly stable steady-state solutions. The results in Figures 6 and
8 correspond to vertical and horizontal slices through the two-parameter plane
in Figure 9. The black dashed curve in Figure 9(a-d) is the Hopf bifurcation
boundary of a single isolated cell, uncoupled from the bulk, as given in
(3.0.3). In each subfigure, for any pair $(\mu,\alpha)$ in a green-shaded
region that lies above this dashed curve, we conclude that the cell-cell
coupling from the bulk medium is the trigger for intracellular oscillations
that would otherwise not occur when the cells are decoupled from the bulk. In
contrast, and most notably for three cells, we observe from Figures 8(a) and
8(c) that there are narrow parameter regimes where the black dashed curve is
above the green region. These thin regimes correspond to where an oscillatory
instability of an isolated cell is quenched by cell-cell coupling through the
bulk medium.

When the cells are all identical, we observe by comparing Figure 8(a) and
Figure 8(b) that the region where synchronous oscillations occur for $m=3$
cells is slightly larger than that for $m=8$ cells. This feature is in
contrast to our observation from Figures 5(a) and 5(b) when the permeabilities
were used as bifurcation parameters, where the parameter region of
intracellular oscillations get larger as the number of cells increases. For
the defective cell case, the existence of an additional regime of periodic
solutions, as observed in Figure 8(d), is apparent in Figure 9(d). This new
regime predicts linearly stable periodic solutions for large values of
$\alpha$, provided that $\mu>0$ is small enough, which is not apparent in the
one-parameter bifurcation diagram for $\alpha$ presented in Figure 6(d). We
notice from the result in Figure 9(b) that synchronous oscillations are not
predicted as $\mu$ approaches zero when the cells are all identical. This
indicates that the stable periodic solutions observed in Figure 9(d) as
$\mu\to 0$ are triggered by the existence of a defective cell.

In Figure 15 of Appendix B we show a gallery of time-dependent solutions to
the ODE system (2.2.18) for different Sel’kov kinetic parameters $\alpha$ and
$\mu$ sampled from the bifurcation diagrams in Figures 6, 8 and 9.

####  3.0.3 Instabilities due to increasing the number of cells

In this subsection, we determine how parameter regimes for linearly stable
periodic solutions change as the number of cells increases. Our findings are
presented as two-parameter Hopf bifurcation diagrams, where the vertical axis
is one of our main bifurcation parameters and the horizontal axis is the
number of cells, the latter of which is treated for convenience as a
continuous rather than a discrete variable. In Figure 10 and Figure 11,
linearly stable periodic solutions are predicted in the green-shaded regions,
while linearly stable solutions are predicted in the unshaded regions. The
parameter values are in the caption of these figures. Since synchronous
oscillations occur as the cell population increases, these results can also be
used to study quorum sensing behavior for a fixed value of each parameter.
Starting with the case of identical cells, Figure 10(a) shows the result for
the permeability parameter $d_{1}$, which agrees with the result in Figure
$18$ of [9]. It predicts that linearly stable periodic solutions can be
triggered in the system for $m\leq 35$. For each population size $m$, and for
fixed $d_{2}$, there is a range of $d_{1}$ for which synchronous oscillations
exist, and this range widens as the number of cells decreases. This trend is
qualitatively different than that shown in Figure 10(b) with regards to
$d_{2}$ at a fixed value of $d_{1}$. For this permeability parameter, linearly
stable periodic solutions exist for all values of $m$ in some finite interval
of $d_{2}$, and this range widens as the number of cells increases. Since
$d_{2}$ is the efflux rate out of a cell, the observation that the Hopf
bifurcation threshold for $d_{2}$ increases as $m$ increases, indicates that
more signaling chemical is needed to synchronize the intraceulluar dynamics of
a larger number of cells. Figure 10(c) shows the corresponding result for the
Sel’kov parameter $\alpha>0$. The black dashed horizontal line in this figure
at $\alpha\approx 0.694$ is the Hopf bifurcation threshold of a single
isolated cell when uncoupled from the bulk. We observe from this figure that
linearly stable periodic solutions can occur for any number of cells, even
when each cell has a linearly stable steady-state when uncoupled from the
bulk. However, as seen in Figure 10(d), this is not the case for $\mu>0$,
where linearly stable synchronous solutions are predicted only for $m\leq 8$
on some range of $\mu$ that decreases as the number of cells increases. No
oscillation is predicted for any value of $\mu$ when $m\geq 9$.

(a) $d_{1}$ versus $m$.

(b) $d_{2}$ versus $m$.

(c) $\alpha$ versus $m$.

(d) $\mu$ versus $m$.

Figure 10: Two-parameter bifurcation diagrams computed from the ODEs (2.2.18)
showing regions of instability with respect to certain parameters and the
number of identical cells. Linearly stable periodic solutions are predicted in
the green-shaded regions, while the unshaded regions correspond to linearly
stable steady-state solutions. The parameters are $d_{1,j}=0.8$ and
$d_{2,j}=0.2$ for $j=1,\dots,m$, with $\mu=2$, $\alpha=0.9$,
$\epsilon_{0}=0.15$ and $\tau=0.5$ (except when used as a bifurcation
parameter). (a) and (b) are for the permeability parameters $d_{1}$ and
$d_{2}$, respectively, (c) shows the result for $\alpha$. The black dashed
horizontal line at $\alpha\approx 0.694$ is the Hopf bifurcation threshold of
a single isolated cell when uncoupled from the bulk. (d) is for $\mu$.

Similar results, but predicting qualitatively different outcomes, are
presented in Figure 11 for the case of a single defective cell coupled to a
group of identical cells. In these results, the bifurcation parameters
$d_{1}$, $d_{2}$, $\alpha$ and $\mu$ are for the defective cell only, and $m$
is the population of identical cells. From Figure 11(a), for $2\leq m\leq 7$,
there is a finite range of permeabilities $d_{1,1}$ for which intracellular
oscillations can occur. On this range of $m$, variations in $d_{1,1}$ can
either trigger or extinguish intracellular oscillations. No oscillations can
occur for $m\geq 8$. From Figure 11(b) similar conclusions hold for the second
permeability parameter $d_{2,1}$. We observe from Figures 11(a) and 11(b) that
variations in the cell population $m$ also play a dual role of triggering and
quenching intracellular oscillations.

(a) $d_{1,1}$ versus $m$.

(b) $d_{1,2}$ versus $m$.

(c) $\alpha$ versus $m$.

(d) $\mu$ versus $m$.

Figure 11: Two-parameter bifurcation diagrams computed from the ODEs (2.2.18)
showing regions of instability with respect to certain parameters and the
number of cells. Linearly stable periodic solutions are predicted in the
green-shaded regions while the unshaded regions correspond to linearly stable
steady-state solutions. The parameters are $d_{1,j}=0.8$ and $d_{2,j}=0.2$ for
$j=1,\dots,m$, with $\mu=2$, $\alpha=0.9$ $,\epsilon_{0}=0.15$ and $\tau=0.5$
(except when used as the bifurcation parameter). The cells are not all
identical, and the bifurcation parameters $d_{1,1}$, $d_{2,1}$, $\alpha$ and
$\mu$ are for the defective cell only. $m$ is the population of identical
cells. (a) and (b) are for the permeability parameters $d_{1,1}$ and
$d_{2,1}$, respectively, (c) is for $\alpha$. The black dashed horizontal line
at $\alpha\approx 0.694$ is the Hopf bifurcation threshold of a single
isolated cell when uncoupled from the bulk. (d) is for $\mu$.

In Figure 11(c) we show that when the number of identical cells with a common
Sel’kov parameter $\alpha=0.9$ satisfies $m\leq 18$, intracellular
oscillations for the group of cells can still occur when a defective cell is
included, provided that the Sel’kov parameter for the defective cell is below
some threshold. Without the defective cell, if there are more than seven
identical cells with a common Sel’kov parameter $\alpha=0.9$ there are no
oscillations, as $m\geq 7$ and $\alpha=0.9$ is not within the green-shaded
region of Figure 10(c). However, on the range $10\leq m\leq 18$, we observe
from Figure 11(c) that to maintain intracellular oscillations the defective
cell is no longer a conditional oscillator, as its required value of $\alpha$
is below the horizontal black dashed line in Figure 11(c). With regards to
$\mu$, Figure 11(d) shows the qualitatively different feature that there are
two regimes of $\mu$ for which intracellular oscillations occur for different
ranges of $m$. As $m$ increases, there is initially only one regime in $\mu$
for which periodic solutions exists. However, as the cell population
increases, a new regime emerges when $m=5$. These two regimes co-exist until
$m=8$, after which the first regime vanishes, and the second one continues
until $m=9$. There are no intracellular oscillations in the system for $m>9$.

##  4 Quorum sensing and phase synchronization: The Kuramoto order parameter

In this section, we use the ODE system (2.2.18) to study quorum sensing (QS)
behavior, characterized by the emergence of collective cell dynamics that
arises from a quiescent state as the cell population density exceeds a
threshold. Since the ratio of the area occupied by $m$ circular cells of a
common radius $\varepsilon$ to the overall domain area is
${m\pi\varepsilon^{2}/|\Omega|}$, our measure of the cell density is taken as
$\rho\equiv{m/|\Omega|}$. With our Sel’kov kinetics, we will define collective
dynamics as the emergence of synchronous oscillations in the intracellular
dynamics as $\rho$ exceeds a critical value.

In terms of the parameter $\rho\equiv m/|\Omega|$, the ODE system (2.2.18)
becomes

$\begin{split}\frac{dU_{0}}{dt}&=-\frac{1}{\tau}U_{0}-\frac{\rho}{m}\sum_{j=1}^{m}(k_{1,j}\,U_{0}-k_{2,j}u_{j}^{1})\,,\quad\frac{d\boldsymbol{u}_{j}}{dt}=\,\boldsymbol{F}_{j}\left(\boldsymbol{u}_{j}\right)+\boldsymbol{e}_{1}\,(k_{1,j}\,U-k_{2,j}\,u_{j}^{1})\,,\qquad
j=1,\ldots,m\,,\end{split}$ (4.0.1)

where we have defined $k_{1,j}\equiv(2\pi d_{1,j})/\tau$ and
$k_{2,j}\equiv(2\pi d_{2,j})/\tau$. For the case of identical cells, we will
perform a Hopf bifurcation analysis on this system in order to identify the
range of instability of the steady-state as $\rho$ is varied. For the case of
non-identical cells, in §4.2 the Kuramoto order parameter (cf. [23], [16],
[15]) will be used to measure the degree of phase synchrony of the
intracellular oscillations.

###  4.1 Hopf bifurcation analysis for identical cells

For the case where the cells are all identical we label $k_{1}\equiv k_{1,j}$
and $k_{2}\equiv k_{2,j}$, and introduce the common local variables
$\boldsymbol{u}=\boldsymbol{u}_{j}$ and
$\boldsymbol{F}(\boldsymbol{u})\equiv\boldsymbol{F}_{j}(\boldsymbol{u}_{j})$
for $j=1,\dots,m$. For this identical cell case, (4.0.1) reduces to

$\begin{split}\frac{dU_{0}}{dt}&=-\frac{1}{\tau}U_{0}-\rho\,(k_{1}\,U_{0}-k_{2}u^{1})\,,\qquad\frac{d\boldsymbol{u}}{dt}=\,\boldsymbol{F}\left(\boldsymbol{u}\right)+\boldsymbol{e}_{1}\,(k_{1}\,U-k_{2}\,u^{1})\,,\end{split}$
(4.1.1)

where $U_{0}\equiv U_{0}(t)$ is the leading-order concentration of the
signaling chemical in the bulk region and
$\boldsymbol{u}=(u^{1},\ldots,u^{n})^{T}$ represents the $n$ interacting
chemical species in the cells. Upon substituting the Sel’kov reaction kinetics
(3.0.1) into (4.1.1), for which $n=2$, we readily calculate that the steady-
state solution $(U_{0e},u_{e}^{1},u_{e}^{2})$ is

$\begin{split}U_{0e}=\frac{\tau\mu\rho k_{2}}{(1+k_{2}+\tau\rho
k_{1})}\,,\quad u_{e}^{1}=\frac{\mu(1+\tau\rho k_{1})}{(1+k_{2}+\tau\rho
k_{1})}\,,\quad
u_{e}^{2}=\frac{\mu}{\left(\alpha+(u_{e}^{1})^{2}\right)}\,,\quad\text{where}\quad\boldsymbol{u}_{e}=(u_{e}^{1},u_{e}^{2})^{T}\,.\end{split}$
(4.1.2)

From (4.1.2) we observe that both the steady-state bulk concentration $U_{0e}$
and $u_{e}^{1}$ are increasing functions of the cell density parameter $\rho$.
Moreover, as expected, $U_{0e}$ is an increasing function of the cell efflux
rate $k_{2}$, but decreases as the influx rate $k_{1}$ into the cells
increases. As the influx rate $k_{1}$ into the cells from the bulk medium
increases, the intracellular level $u_{e}^{1}$ also increases. In the large
cell density limit we obtain from (4.1.2) that

$u_{e}^{1}\to\mu\,,\qquad u_{e}^{2}\to\frac{\mu}{\alpha+\mu^{2}}\,,\qquad
U_{0e}\to\frac{\mu k_{2}}{k_{1}}\,,\qquad\mbox{as}\quad\rho\to\infty\,.$
(4.1.3)

These limiting values for the intracellular species are the same as those for
a single cell that is uncoupled from the bulk medium, as summarized in §3.

To determine the linear stability of this steady-state, we introduce the
perturbation

$\begin{split}U&=U_{0e}+e^{\lambda
t}\eta\,,\qquad\boldsymbol{u}=\boldsymbol{u}_{e}+e^{\lambda
t}\boldsymbol{\phi}\,,\end{split}$ (4.1.4)

where $|\boldsymbol{\phi}|\ll 1$ and $\eta\ll 1$. Upon substituting (4.1.4)
into (4.1.1), we obtain the linearized system

$\begin{split}\lambda\eta&=-\frac{\eta}{\tau}-\rho(k_{1}\eta-
k_{2}\phi_{1})\,,\qquad\lambda\boldsymbol{\phi}=J_{e}\boldsymbol{\phi}+\boldsymbol{e}_{1}(k_{1}\eta-
k_{2}\phi_{1})\,,\end{split}$ (4.1.5)

where $J_{e}$ is the Jacobian matrix of the local kinetics
$\boldsymbol{F}(\boldsymbol{u})=(F(u^{1},u^{2}),G(u^{1},u^{2}))^{T}$ evaluated
at the steady-state solution $\boldsymbol{u}_{e}=(u_{e}^{1},u_{e}^{2})^{T}$,
and $\boldsymbol{\phi}=(\phi_{1},\phi_{2})^{T}$. For convenience of notation,
we let $(v,w)\equiv(u_{e}^{1},u_{e}^{2})$, and write the linearized system
(4.1.5) in matrix form as

$\begin{split}\mathcal{H}(\lambda)\boldsymbol{\Psi}=\boldsymbol{0}\,,\end{split}$
(4.1.6)

where $\boldsymbol{\Psi}=(\eta,\phi_{1},\phi_{2})^{T}$ and
$\mathcal{H}(\lambda)$ is the $3\times 3$ matrix defined by

$\begin{split}\mathcal{H}(\lambda)=\begin{pmatrix}-\left(\frac{1}{\tau}+\rho
k_{1}\right)-\lambda&\rho k_{2}&0\\\ \\\
k_{1}&(F_{v}^{e}-k_{2}-\lambda)&F_{w}^{e}\\\ \\\
0&G_{v}^{e}&(G_{w}^{e}-\lambda)\end{pmatrix}\,.\end{split}$ (4.1.7)

Here $F_{v}^{e},F_{w}^{e},G_{v}^{e}$, and $G_{w}^{e}$ denote the partial
derivatives of $F$ or $G$ with respect to $v$ or $w$ evaluated at the steady-
state $\boldsymbol{u}_{e}$. A simple calculation shows that there is a
nontrivial solution to (4.1.6) if only if $\lambda$ is a root of the cubic

$\begin{split}\lambda^{3}+p_{1}\lambda^{2}+p_{2}\lambda+p_{3}=0\,,\end{split}$
(4.1.8a) whose coefficients are given by
$\begin{split}p_{1}\equiv\left(k_{2}+\frac{1}{\tau}+\rho
k_{1}\right)&-\text{tr}(J_{e})\,,\qquad
p_{2}\equiv(1+k_{2})\det(J_{e})-\left(\frac{1}{\tau}+\rho
k_{1}\right)\text{tr}(J_{e})+\frac{k_{2}}{\tau}\,,\\\
p_{3}&\equiv\left(\frac{1}{\tau}+\rho
k_{1}+\frac{k_{2}}{\tau}\right)\det(J_{e})\,.\end{split}$ (4.1.8b)

Here $\det(J_{e})=F_{v}^{e}G^{e}_{w}-F_{w}^{e}G^{e}_{v}>0$ and
$\text{tr}(J_{e})=F_{v}^{e}+G_{w}^{e}$ are the determinant and trace of
$J_{e}$ for Sel’kov reaction kinetics, for which we readily calculate that
$\det(J_{e})=-G_{w}^{e}=\epsilon_{0}(\alpha+v_{e}^{2})$. By the Routh-Hurwitz
criterion for cubic functions, the three roots of the polynomial (4.1.8a)
satisfy $\mbox{Re}(\lambda)<0$ if and only if the following conditions hold

$\begin{split}p_{1}>0\,,\qquad p_{3}>0\,,\qquad\text{and}\quad
p_{1}p_{2}>p_{3}\,.\end{split}$ (4.1.9)

From (4.1.8b) we observe that $p_{3}>0$ always holds. Since we are interested
in finding Hopf bifurcation points, we need to consider a special cubic
polynomial whose roots satisfy $\lambda_{1}=a<0$, $\lambda_{2,3}=\pm i\omega$,
for which

$\begin{split}(\lambda-a)(\lambda-i\omega)(\lambda+i\omega)=\lambda^{3}-a\lambda^{2}+\omega^{2}\lambda-a\omega^{2}=0\,.\end{split}$
(4.1.10)

Upon comparing the polynomials (4.1.8a) and (4.1.10), we conclude that a Hopf
bifurcation threshold must satisfy

$\begin{split}p_{1}>0\,,\quad p_{3}>0\,,\quad\text{and}\quad
p_{1}p_{2}=p_{3}\,.\end{split}$ (4.1.11)

For the parameter values in the caption of Figure 12, we use (4.1.11) to
numerically compute two Hopf bifurcation points at $\rho_{1}=0.8010$ and
$\rho_{2}=2.6750$. These results agree with those computed in Figure 12 using
XPPAUT [6].

(a) One-parameter Hopf Bifurcation for $u_{1e}$ vs $\rho$

(b) Hopf bifurcation diagrams for $m$ vs $|\Omega|$

Figure 12: (a) Global bifurcation diagrams for $u_{1e}$ versus the cell
population density parameter $\rho\equiv{m/|\Omega|}$, as computed from the
ODEs (4.1.1), showing steady-states and global branches of periodic solutions
for the case of identical cells. The red and black lines correspond to
linearly stable and unstable steady-state solution branches, respectively. The
green loop indicates a linearly stable branch of periodic solutions. The Hopf
bifurcation points are at $\rho_{1}=0.8010$ and $\rho_{2}=2.6750$. (b) A two-
parameter bifurcation diagram for the number of cells $m$ and the domain area
$|\Omega|$. The green wedge-shaped region predicts stable periodic solutions
while the unshaded regions correspond to linearly stable steady-states. In (a)
and (b) the Sel’kov parameters are $\mu=2,\alpha=0.9$ and $\epsilon_{0}=0.15$,
with $\tau=0.5,k_{1}=10.0531$ and $k_{2}=2.5133$, which corresponds to the
permeabilities $d_{1}=0.8$ and $d_{2}=0.2$.

Figure 12(a) shows a one-parameter Hopf bifurcation diagram for the ODE system
(4.1.1) with respect to the cell population density $\rho$ parameter. Observe
from this figure that the range of $\rho$ that predicts stable periodic
solutions is bounded. This indicates that the cell population density plays a
dual role of both triggering and then quenching intracellular oscillations,
which agrees with our earlier observation in Figure 11. Our results in Figure
12(a) show that there are no synchronous oscillations for $\rho<0.8010$ and
for $\rho>2.6750$, and that linearly stable periodic solutions exist only on
the range $0.8010\leq\rho\leq 2.6750$. As $\rho$ increases above
$\rho=0.8010$, synchronous oscillations are triggered through quorum sensing
behavior. These two Hopf bifurcation values for $\rho$ generate the green
wedge-shaped region shown in Figure 12(b) in the $m$ versus $|\Omega|$
parameter plane where synchronous intracellular oscillations occur. For any
vertical slice through the bifurcation diagram in Figure 12(b), corresponding
to a domain of fixed area $|\Omega|$, there is a quorum sensing threshold for
the number of cells required to trigger synchronous oscillations. Moreover,
there is an additional threshold on the number of cells for which the
oscillations are quenched. This clearly shows that the cells adjust their
intracellular dynamics in accordance with the number of cells in the
population. When the domain area is small, fewer cells are required to trigger
oscillations, while as the domain area increases the number of cells required
for quorum sensing also increases. On the other hand, the quenching of
oscillations predicted by Figure 12(b) results from having too many cells in a
given domain area. As the domain area increases, the number of cells required
to quench oscillations also increases, as expected intuitively. The
qualitative mechanism underlying the quenching of oscillations is that when
$\rho\gg 1$ the collection of cells behaves like a single “effective” cell,
which has no intracellular oscillations and a unique linearly stable steady-
state as result of the Sel’kov parameters chosen. Mathematically, as
$\rho\to\infty$, we have $v_{e}\to\mu$ and $w_{e}\sim{\mu/(\alpha+\mu^{2})}$
from (4.1.3), and so from (3.0.2) and the fact that the Sel’kov kinetic
parameters were chosen so that a single isolated cell has a linearly steady-
state, we must have $\text{tr}(J_{e})<0$ for $\rho\gg 1$. Therefore, from
(4.1.8b) we obtain that $p_{1}>0$ for $\rho$ sufficiently large. Moreover,
since from (4.1.8b) we have

$p_{1}p_{2}\sim\rho^{2}k_{1}^{2}|\text{tr}(J_{e})|={\mathcal{O}}(\rho^{2})\,,\qquad
p_{3}\sim\rho
k_{1}\left(\det{J_{e}}\right)={\mathcal{O}}(\rho)\,,\qquad\mbox{as}\quad\rho\to\infty\,,$
(4.1.12)

we conclude that $p_{1}p_{2}>p_{3}$ for $\rho$ sufficiently large. Therefore,
when $\rho$ is sufficiently large the Routh-Hurwitz conditions in (4.1.9) are
satisfied. As a result, there must be a critical value $\rho_{m}$ of the cell
density parameter $\rho$, with $\rho_{m}$ sufficiently large, for which the
steady-state of the ODE system (4.1.1) is linearly stable for all
$\rho>\rho_{m}$.

We emphasize that the two bifurcation diagrams in Figure 12 are related. A
vertical slice of the one-parameter Hopf bifurcation diagram in Figure 12(a)
at a population density of $\rho=\rho_{s}$ corresponds to a straight line with
slope $\rho_{s}$ (since $\rho=m/|\Omega|$) passing through the origin in the
two-parameter bifurcation diagram in Figure 12(b). The Hopf bifurcation points
$\rho_{1}=0.8010$ and $\rho_{2}=2.6750$ in Figure 12(a) are the slopes of the
Hopf bifurcation boundaries in Figure 12(b). Any straight line passing through
the origin and contained between these boundaries corresponds to a vertical
slice of the one-parameter bifurcation diagram in Figure 12(a) at a $\rho$
value between $\rho_{1}$ and $\rho_{2}$.

###  4.2 Phase synchronization and the order parameter

Next, we study quorum sensing using the ODE system (4.0.1) together with the
Kuramoto order parameter. The Kuramoto order parameter was originally used to
measure the degree of phase synchrony of a system of coupled oscillators [13].
We will use the version of the order parameter presented in [15] and [16],
given by

$\displaystyle
R=\left\langle\left|N^{-1}\sum_{j=1}^{N}\,\,\exp[i\theta_{j}(t)]-\left\langle
N^{-1}\sum_{j=1}^{N}\,\,\exp[i\theta_{j}(t)]\right\rangle\right|\right\rangle\,,$
(4.2.1)

which was used in [15] and [16] to measure the degree of phase synchrony of
indirectly coupled chaotic oscillators. Here $N$ is the number of oscillators,
$\theta_{j}(t)$ is the instantaneous phase of the $j^{\text{th}}$ oscillator,
and $\langle\dots\rangle$ represents average over time. The value of $R$
indicates the level of phase synchronization of the oscillators, and it varies
between 0 and 1, inclusive. When $R=1$, the oscillators are perfectly in
phase, while the value $R=0$ indicates that they are perfectly out of phase.
Two oscillators are said to be in phase if their frequency and phase angles
are the same at each cycle as they oscillate. Our goal is to use the system of
ODEs (4.0.1) together with the order parameter (4.2.1) to measure the degree
of phase synchronization of the intracellular dynamics within the signaling
compartments as the population density $\rho$ is varied. In our numerical
examples below, the cell population will be taken to be a mixture of identical
and defective cells, where the heterogeneous cells have one different Sel’kov
parameter.

To measure the degree of synchronization of the cells, we first solved the
coupled system of ODEs (4.0.1) numerically using ODE45 in MATLAB [27]. The
solution over the transient period was discarded, after which we extracted the
time evolution of the same chemical species in all the cells. Only one of the
intracellular chemical species was used for this analysis because both of them
started oscillating and synchronizing at the same time. A single mode Fourier
series expansion was fitted to each of the selected solutions, and their
instantaneous phase $\theta(t)$ was computed from the coefficients of the
series. At each time point, the phase of each oscillator was mapped to the
complex number $e^{i\theta}$ (where $i=\sqrt{-1}$) on the unit circle, and an
average of these points was computed using
$z=N^{-1}\sum_{j=1}^{N}e^{i\theta_{j}}$. For better intuition of the meaning
of this average, take each point on the unit circle to be a unit vector. Then
$z$ is a vector whose direction is the average direction of all the vectors on
the unit circle. The length of this average vector is bounded between 0 and 1,
which indicates the degree of phase synchrony of the oscillators at that time
point. For instance, if two oscillators are perfectly synchronized at a
specific time point, their phases at that point will be the same and they will
both be mapped to the same point on the unit circle. The average of these two
points/vectors gives another unit vector in the same direction. The fact that
their average is a unit vector verifies that the two oscillators are in
perfect synchrony at this point. On the other hand, if the two oscillators are
completely out of phase, their phases will be mapped to opposite ends of the
unit circle, and their average will be the zero vector, so that $|z|=0$. The
more synchronized oscillators are at a time point, the closer the
corresponding mapping of their phases are on the unit circle, and as a result
the closer the value of $|z|$ is to 1. Thus, $|z|$ indicates the level of
synchrony of the oscillators at each time point. Upon computing the
instantaneous averages at each time point, we computed their average over a
specified time interval (after the system has reached a quasi steady-state).
Lastly, the modulus of the difference between the instantaneous averages $z$
and the time-average $\langle z\rangle$ was computed for each time point, and
then the corresponding result was averaged over time to calculate the order
parameter $R$, as given in (4.2.1) (cf. [16], [23]). In our computations
below, we set $R=0$ when the cells are in a quiescent state and whenr
oscillations have an amplitude less than $1\times 10^{-4}$.

(a) $1000$ identical cells

(b) $600$ identical and $400$ defective cells

(c) $100$ identical and $900$ defective cells

(d) $1000$ identical cells

(e) $600$ identical and $400$ defective cells

(f) $100$ identical and $900$ defective cells

Figure 13: The degree of phase synchrony and amplitudes of the oscillations
predicted in $1000$ coupled cellular compartments measured with the order
parameter $R$ of (4.2.1), as computed from the ODE system (4.0.1). (a), (b),
and (c) show the computed values of $R$ with respect to $\rho$, for $1000$
identical cells, $600$ identical and $400$ defective cells, and $100$
identical and $900$ defective cells, respectively. (d), (e), and (f) show the
corresponding average amplitudes of the oscillations in the cells (black
curve) and in the bulk region (red curve). The Sel’kov parameter $\alpha=0.9$
is used for the identical cells, while for the defective cells it is chosen
uniformly from the range $0.9220\leq\alpha\leq 0.9520$. The other parameters
are $\mu=2,\epsilon_{0}=0.15,\tau=0.5,k_{1}=10.0531$ and $k_{2}=2.5133$.

In terms of the cell population density $\rho$, in Figure 13 we plot the
numerically computed order parameter $R$ and the corresponding average
amplitudes of predicted oscillations, as computed from the ODEs (4.0.1), for
$1000$ cells. Since the cell population is fixed at $m=1000$, increasing the
population density $\rho$ corresponds to decreasing the area $|\Omega|$ of the
2-D bounded domain $\Omega$. The results in Figure 13 are for the Sel’kov
kinetic parameters $\mu=2$, $\epsilon_{0}=0.15$, the reaction-time parameter
$\tau=0.5$, and the permeabilities $k_{1}=10.0531$ and $k_{2}=2.5133$. The
remaining Sel’kov parameter value $\alpha$ was fixed at $\alpha=0.9$ for the
identical cells, and was chosen uniformly from the interval
$0.9220\leq\alpha\leq 0.9520$ for the defective cells. As seen from Figure 2,
each cell has Sel’kov kinetic parameters that are chosen so that intracellular
steady-state is linearly stable when there is no coupling to the bulk medium.
From Figure 13 we observe that there are unique critical values of $\rho$ at
which the order parameter $R$ switches from $0$ to $1$ (perfectly out of phase
to perfect synchrony) and from $1$ to $0$ (perfect synchrony to perfectly out
of phase). These transitions occur in a switch-like manner, in the sense that
when intracellular oscillations emerge they all synchronize simultaneously.
When $R$ switches from $0$ to $1$, the synchronous intracellular oscillations
that are triggered can be interpreted as a quorum sensing transition in the
sense that it occurs for $1000$ cells only when the domain area $|\Omega|$ is
small enough. Similarly, the cells desynchronize simultaneously as $\rho$
increases past a further threshold, leading to a linearly stable steady state.
For identical cells, this linearly stable steady-state when $\rho$ exceeds a
threshold was predicted theoretically in §4.1. Having two critical values of
$\rho$, labeled by $\rho_{1}$ and $\rho_{2}$, at which the dynamics of the
cells change state agrees with our earlier observation in Figure 12 that cell
population density plays a dual role of initiating and then quenching
intracelluar oscillations. When the cells are all identical, the critical
population densities are $\rho_{1}=0.8$ and $\rho_{2}=2.77$, when 600 of the
cells are identical (remaining 400 are defective) we have $\rho_{1}=0.9$ and
$\rho_{2}=2.34$ , and when only $100$ of the cells are identical (remaining
900 are defective) we have $\rho_{1}=1.14$ and $\rho_{2}=1.76$. For any $\rho$
in $\rho_{1}\leq\rho<\rho_{2}$, we have $R=1$ and perfect synchrony in the
intracellular dynamics. On the other hand, when $\rho<\rho_{1}$, we have $R=0$
so that the cells are either perfectly out of phase or there are no
oscillations. This occurs when there are too few cells relative to the area of
the domain, as there is insufficient secreted chemical from the cells to
trigger synchronous oscillations. Similarly, $R=0$ when $\rho>\rho_{2}$, since
in this case the domain area is sufficiently small so much that the cells
behave like a single effective cell, which has a linearly stable steady state.

In Figure 13 we observe that an increase in the relative proportion of
defective cells in the population significantly alters the critical values
$\rho_{1}$ and $\rho_{2}$ for which synchronous oscillations are triggered and
quenched, respectively. As the number of defective cells increases, the range
of $\rho$ that predicts perfect synchrony shrinks (see Figures 13(a), 13(b),
and 13(c)), and so does the amplitudes of the predicted oscillations (see
Figures 13(d), 13(e), and 13(f)). This shows that synchronous intracellular
dynamics, resulting from effective communication between the cells, becomes
more difficult to sustain as the relative number of defective cells increases.

##  5 Discussion

We have studied the emergence of oscillatory intracellular dynamics for a
collection of $m$ small dynamically active disk-shaped signaling compartments,
referred to as “cells”, that are coupled together through a scalar passive
bulk diffusion field in a bounded 2-D domain. In our theoretical model,
originating from [19] (see also [20] and [32]) and later extended and
formalized in [9], each cell secretes across its boundary only one signaling
chemical into the common bulk medium, while it receives feedback across its
boundary from the entire collection of cells. This boundary efflux and influx
of signaling material is regulated by permeability parameters. The resulting
PDE-ODE model (1.0.1), given in dimensionless form in (2.1.1), was studied in
the limit of a small common cell radius $\varepsilon\ll 1$ and in the limit
where the bulk diffusion coefficient $D$ is asymptotically large, on the range
$D\gg\mathcal{O}(\nu^{-1})\gg\mathcal{O}(1)$ where
$\nu\equiv{-1/\log\varepsilon}\ll 1$. In this limit, where the bulk region
becomes well-mixed and the bulk concentration becomes spatially homogeneous, a
matched asymptotic analysis was used to reduce the dimensionless PDE-ODE model
into the system of ODEs (2.2.18) with global coupling. A similar reduction,
but for identical cells, was done in [9].

For the specific case of Sel’kov reaction kinetics, we have used this ODE
system to provide a detailed study of the emergence of intracellular
oscillations in a collection of cells that is due to the cell-cell coupling
through the bulk medium, and which otherwise would not occur for isolated
cells. As such we have assumed that the Sel’kov kinetic parameters are chosen
so that each cell is a conditional oscillator, and so its dynamics has only a
stable steady-state when uncoupled from the bulk. For identical cells, the
emergence of intracellular oscillations is found to arise via a Hopf
bifurcation of the globally coupled system as parameters are varied.

One key focus has been to provide a detailed bifurcation study of the role of
cell membrane permeabilites and the Sel’kov kinetic parameters on the
emergence of intracellular oscillations. Our numerical bifurcation study using
XPPAUT [6] has considered two main scenarios: one where all the cells are
identical, and the other where there is a defector cell that is coupled to a
group identical cells. A cell is considered defective if its Sel’kov kinetic
parameters or permeabilities are different from those of the group of
identical cells. Our results were presented in the form of one- and two-
parameter global bifurcation diagrams of steady-states and periodic solution
branches of the ODEs (2.2.18). Some highlights of the effect of a defector
cell on a small group of identical cells include the following: Small changes
in the influx rate of a single defective cell, relative to the common influx
rate of the group, can act as a switch that extinguishes intracellular
oscillations for the entire collection of cells (see Figures 3(c) and 3(d)).
Intracellular oscillations of an entire group of cells can be extinguished by
a single defective cell that reduces its efflux rate (see Figures 4(c) and
4(d)). Intracellular oscillations of an entire group of identical cells can be
extinguished (see Figure 6) or triggered (see Figure 7) through increases or
decreases, respectively, in the Sel’kov kinetic parameter $\alpha$ of one
defective cell. A more detailed discussion of the effect of a defector cell on
a small group of identical cells is given in §3.

Our second key focus was to use the ODE system (2.2.18) to study quorum
sensing behavior associated with a large collection of either identical cells
or a mixture of identical and defective cells. In terms of a cell density
parameter, defined by $\rho\equiv{m/|\Omega|}$, we showed that quorum sensing
is characterized by the switch-like emergence of synchronous intracellular
oscillations as $\rho$ crosses through a critical threshold $\rho_{1}$ (see
Figures 12 and 13). Moreover, we showed that further increases in the cell
density parameter eventually leads to a switch-like quenching of the
oscillations when a second threshold density $\rho_{2}$ is reached. For the
case of identical cells, we showed that these two thresholds $\rho_{1}$ and
$\rho_{2}$ are Hopf bifurcation points associated with a cubic polynomial (see
(4.1.8)), while the existence of the second threshold follows from (4.1.12).
The qualitative explanation underlying the quenching threshold is that for
large $\rho$ the collection of cells behaves like a single “effective” cell,
which has no intracellular oscillations since it is a conditional Sel’kov
oscillator. For two different mixtures of identical and defective cells, for
which the Sel’kov kinetic parameter $\alpha$ was uniformly distributed in some
interval, the Kuramoto order parameter in (4.2.1) was computed numerically to
establish intervals in the cell density parameter $\rho$ where synchronous
intracellular oscillations occur.

There are several possible extensions of our modeling framework and analysis
that should be considered. One main direction would be to analyze (2.1.1) for
the finite bulk diffusion regime $D={\mathcal{O}}(1)$ for an arbitrary
collection of non-identical small cells in a bounded 2-D domain. In
particular, it would be interesting to investigate the effect of various
spatial factors, such as diffusion sensing and the role of cell-clustering
(cf. [8]), on the emergence of intracellular oscillations for a collection of
cells whose reaction kinetics are conditional Sel’kov oscillators. In this
context, the spatial diffusion gradient and the overall spatial configuration
of cells should determine which cells can effectively communicate and
synchronize their dynamics. In particular, is a chimera-type state involving
distinct intracellular oscillations for different cell-clusters within the
entire group possible in this parameter regime?

Our PDE-ODE system (2.1.1) with Sel’kov reaction kinetics provides a simple
conceptual model that predicts the switch-like emergence of intracellular
oscillations with increasing cell density, known as quorum sensing, and this
feature is qualitatively similar to that seen in many biological and chemical
settings (cf. [4], [10], [12], [18], [14], [29], [28], [31], [30]). As a
second main direction, it would be worthwhile to investigate triggered
intracellular oscillations and quorum sensing behavior when the Sel’kov
kinetics in our theoretical PDE-ODE system (2.1.1) are replaced with specific
biologically-motivated reaction kinetics, such as the detailed model of [34]
and [11] for the glycolytic pathway in yeast cells and the model of [18] for
bacterial communication of Escherichia coli cells. In other settings, such as
in the modeling of bioluminescence emitted by the marine bacterium Vibrio
fischeri, the quorum sensing behavior of microbial cells due to an autoinducer
bulk field does not involve the triggering of intracellular oscillations, but
instead is believed to involve the sudden transition to a new stable steady-
state (the bioluminescent state) as the cell population increases past a
threshold (cf. [17]). It would be interesting to incorporate the intracellular
signaling pathways of [17] into our PDE-ODE system (2.1.1) to show that quorum
sensing in this context results from a saddle-node structure of equilibria and
the sudden transition between stable states as the number of cells increases.
Other modeling approaches to analyze quorum sensing behavior associated with
sudden transitions between stable steady-states as the cell density increases
past a threshold are given in [33] and [5].

Finally, other cell-cell communication problems mediated by a bulk diffusion
field involve two-bulk diffusing species that mediate interactions between
small signalling compartments. In this direction, it would be interesting to
reframe the class of agent-based models considered in [24] for Turing pattern
formation into an extended version of our PDE-ODE system (2.1.1) where there
are now two-bulk diffusing species, consisting of an activator and an
inhibitor field. For this scenario, it would be interesting to investigate, by
using the permeabilities as bifurcation parameters, whether a Turing
bifurcation can occur even when the activator and inhibitor bulk diffusion
fields have very similar diffusion coefficients, as is typical in real-world
biological systems.

## Acknowledgements

Michael Ward gratefully acknowledges the financial support from the NSERC
Discovery grant program.

##  Appendix A  Non-dimensionalization of the coupled PDE-ODE model

In this appendix, the dimensional coupled PDE-ODE model (1.0.1) is non-
dimensionalized. With $[z]$ denoting the dimension of $z$, the dimensional
units of the variables in this model are

$\displaystyle[\mathcal{U}]$
$\displaystyle=\frac{\text{moles}}{(\text{length})^{2}}\,,\qquad[\mu_{c}]=\text{moles}\,,\qquad[D_{B}]=\frac{(\text{length})^{2}}{\text{time}}\,,\qquad[k_{R}]=[k_{B}]=\frac{1}{\text{time}}\,,$
(A.0.1a) $\displaystyle\quad[\boldsymbol{\mu}_{j}]$
$\displaystyle=\text{moles}\,,\qquad[\beta_{1,j}]=\frac{\text{length}}{\text{time}}\,,\qquad[\beta_{2,j}]=\frac{1}{\text{length}\times\text{time}}\quad\text{for}\quad
j=1,\ldots,m\,.$ (A.0.1b)

We assume that the cells are circular with a common radius $R_{0}$, which is
small relative to the length-scale of the domain $L$, and so we define a small
parameter $\varepsilon={R_{0}/L}\ll 1$. We then introduce dimensionless
variables defined by

$\displaystyle
U=\frac{L^{2}}{\mu_{c}}\,\mathcal{U}\,,\qquad\boldsymbol{x}=\frac{\boldsymbol{X}}{L}\,,\qquad
t=k_{R}\,T\,,\qquad\boldsymbol{u}_{j}=\frac{\boldsymbol{\mu}_{j}}{\mu_{c}}\,.$
(A.0.2)

Observe that the dimensionless time-scale that is chosen is the time-scale for
the reaction kinetics of the cells. By introducing (A.0.2) into (1.0.1), we
obtain that the dimensionless bulk concentration satisfies

$\displaystyle\frac{k_{R}}{k_{B}}\,U_{t}=D\Delta U-\,U\,,\qquad$
$\displaystyle\boldsymbol{x}\in\Omega\setminus\cup_{j=1}^{m}\Omega_{\varepsilon_{j}}\,;$
(A.0.3a)
$\displaystyle\varepsilon\,D\,\partial_{n}U=\frac{\varepsilon\,\beta_{1,j}}{L\,k_{B}}\,\,U-\frac{\varepsilon\,\beta_{2,j}L}{k_{B}}\,\,u_{j}^{1}$
$\displaystyle\qquad\text{on}\quad\partial\Omega_{\varepsilon_{j}}\,;\qquad\partial_{n}U=0\quad\text{on}\quad\partial\Omega\,,$
(A.0.3b) where $D\equiv D_{B}/(L^{2}\,k_{B})$ is the dimensionless bulk
diffusivity, $\Omega_{\varepsilon_{j}}$ for $j=1,\dots,m$ are disks of a
common radius $\varepsilon$, and $\Delta$ is the Laplace operator with respect
to $\boldsymbol{x}$. This PDE for the dimensionless bulk concentration is
coupled to the intracellular dynamics of the $j^{\text{th}}$ cell using the
dimensionless integro-differential equation given by
$\displaystyle\frac{\text{d}\boldsymbol{u}_{j}}{\text{d}t}=\boldsymbol{F}_{j}(\boldsymbol{u}_{j})+\frac{\boldsymbol{e}_{1}}{\varepsilon}\,\frac{k_{B}}{k_{R}}\,\int_{\partial\Omega_{\varepsilon_{j}}}\left(\frac{\varepsilon\,\beta_{1,j}}{L\,k_{B}}\,U-\frac{\varepsilon\,\beta_{2,j}L}{k_{B}}\,u_{j}^{1}\right)\,\,\text{d}S_{x}\,,$
(A.0.3c)

where we used $\text{d}S_{\boldsymbol{X}}=L\,\text{d}S_{\boldsymbol{x}}$. In
terms of the dimensionless variables

$\displaystyle\tau\equiv\frac{k_{R}}{k_{B}}\,,\qquad
D\equiv\frac{D_{B}}{L^{2}\,k_{B}}\,,\qquad
d_{1,j}\equiv\frac{\varepsilon\,\beta_{1,j}}{L\,k_{B}}\,,\quad\text{and}\quad
d_{2,j}\equiv\frac{\varepsilon\,\beta_{2,j}L}{k_{B}}\,,$ (A.0.4)

the dimensionless form of the coupled PDE-ODE model (1.0.1) is given by

$\displaystyle\tau\,U_{t}=D\Delta U-\,U\,,\qquad$
$\displaystyle\boldsymbol{x}\in\Omega\setminus\cup_{j=1}^{m}\Omega_{\varepsilon_{j}}\,;$
(A.0.5a)
$\displaystyle\varepsilon\,D\,\partial_{n}U=d_{1,j}\,\,U-d_{2,j}\,\,u_{j}^{1}$
$\displaystyle\qquad\text{on}\quad\partial\Omega_{\varepsilon_{j}}\,;\qquad\partial_{n}U=0\quad\text{on}\quad\partial\Omega\,,$
(A.0.5b) in the bulk region. This PDE is coupled to the intracellular dynamics
of the $j^{\text{th}}$ cell through
$\displaystyle\frac{\text{d}\boldsymbol{u}_{j}}{\text{d}t}=\boldsymbol{F}_{j}(\boldsymbol{u}_{j})+\frac{\boldsymbol{e}_{1}}{\varepsilon\,\tau}\,\,\int_{\partial\Omega_{\varepsilon_{j}}}\left(d_{1,j}\,U-d_{2,j}\,u_{j}^{1}\right)\,\,\text{d}S_{x}\,.$
(A.0.5c)

##  Appendix B  Gallery of Numerical Simulations

Here, we present some numerical simulations of the reduced system of ODEs
(2.2.18).

(a) 2 identical,1 defective

(b) 2 identical,1 defective

(c) 7 identical,1 defective

(d) 7 identical,1 defective

(e) 2 identical,1 defective

(f) 2 identical,1 defective

(g) 7 identical,1 defective

(h) 7 identical,1 defective

(i) 2 identical,1 defective

(j) 2 identical,1 defective

(k) 7 identical,1 defective

(l) 7 identical,1 defective

Figure 14: Numerical simulations of the system of ODEs (2.2.18) for different
values of the permeability parameters $d_{1}$ and $d_{2}$ sampled from the
bifurcation diagrams in Figure 3 (a-d), Figure 4 (e-f), and Figure 5 (i-l).

(a) 2 identical,1 defective

(b) 2 identical,1 defective

(c) 7 identical,1 defective

(d) 7 identical,1 defective

(e) 2 identical,1 defective

(f) 2 identical,1 defective

(g) 7 identical,1 defective

(h) 7 identical,1 defective

(i) 2 identical,1 defective

(j) 2 identical,1 defective

(k) 7 identical,1 defective

(l) 7 identical,1 defective

Figure 15: Numerical simulations of the system of ODEs (2.2.18) for different
values of the Sel’kov parameters $\alpha$ and $\mu$ sampled from the
bifurcation diagrams in Figure 6 (a-d), Figure 8 (e-f), and Figure 9 (i-l).

## References

  * [1] K. L. Asfahl and M. Schuster.  Social interactions in bacterial cell-cell signaling.  FEMS Microbiology Reviews, 41(1):92–107, 2017.
  * [2] M. Chaplain, M. Ptashnyk, and M. Sturrock.  Hopf bifurcation in a gene regulatory network model: Molecular movement causes oscillations.  Math. Mod. Meth. Appl. Sci., 25(6):1179–1215, 2015.
  * [3] S. Danø, M. F. Madsen, and P. G. Sorensøn.  Quantitative characterization of cell synchronization in yeast.  Proceedings of the National Academy of Sciences, 104(31):12732–12736, 2007.
  * [4] S. De Monte, F. d’Ovidio, S. Danø, and P. G. Sørensen.  Dynamical quorum sensing: Population density encoded in cellular dynamics.  Proceedings of the National Academy of Sciences, 104(47):18377–18381, 2007.
  * [5] J. D. Dockery and J. P. Keener.  A mathematical model for quorum sensing in pseudomonas aeruginosa.  Bull Math Biol., 63(1):95–116, 2001.
  * [6] B. Ermentrout and Analyzing Simulating.  Animating dynamical systems: A guide to xppaut for researchers and students (software, environments, tools).  Society for Industrial and Applied Mathematics, 2002.
  * [7] K. Fujimoto and Sawai S.  A design principle of group-level decision making in cell populations.  PLoS Computational Biology, 9(6):e1003110, 2013.
  * [8] M. Gao, H. Zheng, Y. Ren, R. Lou, F. Wu, X. Yu, W. Liu, and X. Ma.  A crucial role for spatial distribution in bacterial quorum sensing.  Scientific reports, 6:34695, 2016.
  * [9] J. Gou and M. J. Ward.  An asymptotic analysis of a 2-D model of dynamically active compartments coupled by bulk diffusion.  Journal of Nonlinear Science, 26(4):979–1029, 2016.
  * [10] T. Gregor, K. Fujimoto, N. Masaki, and S. Sawai.  The onset of collective behavior in social amoebae.  Science, 328(5981):1021–1025, 2010.
  * [11] M. A. Henson, Müller D., and M. Reuss.  Cell population modelling of yeast glycolytic oscillations.  Biochem J., 368:433–446, 2002.
  * [12] K. Kamino, K. Fujimoto, and Sawai S.  Collective oscillations in developing cells: Insights from simple systems.  Develop. Growth Differ., 53:503–517, 2011.
  * [13] Y. Kuramoto.  Self-entrainment of a population of coupled non-linear oscillators.  In International symposium on mathematical problems in theoretical physics, pages 420–422. Springer, 1975.
  * [14] E. J. Leaman, B. Q. Geuther, and B. Behkam.  Quantitative investigation of the role of intra-/intercellular dynamics in bacterial quorum sensing.  ACS synthetic biology, 7(4):1030–1042, 2018.
  * [15] B. W. Li, X. Z. Cao, and C. Fu.  Quorum sensing in populations of spatially extended chaotic oscillators coupled indirectly via a heterogeneous environment.  Journal of Nonlinear Science, 27(6):1667–1686, 2017.
  * [16] B. W. Li, C. Fu, H. Zhang, and X. Wang.  Synchronization and quorum sensing in an ensemble of indirectly coupled chaotic oscillators.  Physical Review E, 86(4):046207, 2012.
  * [17] P. Melke, P. Sahlin, A. Levchenko, and H. Jonsson.  A cell-based model for quorum sensing in heterogeneous bacterial colonies.  PLoS Computational Biology, 6(6):e1000819, 2010.
  * [18] P. Mina, M. di Benardo, N. J. Savery, and K. Tsaneva-Atanasova.  Modeling emergence of oscillations in communicating bacteria: a structured approach from one to many cells.  J. Royal Society Interface, 10:20120612, 2012.
  * [19] J. Müller, C. Kuttler, B. A. Hense, M. Rothballer, and A. Hartmann.  Cell–cell communication by quorum sensing and dimension-reduction.  Journal of Mathematical Biology, 53(4):672–702, 2006.
  * [20] J. Müller and H. Uecker.  Approximating the dynamics of communicating cells in a diffusive medium by ODEs — homogenization with localization.  Journal of Mathematical Biology, 67(5):1023–1065, 2013.
  * [21] V. Nanjundiah.  Cyclic AMP oscillations in dictyostelium discoideum: models and observations.  Biophysical Chemistry, 72(1-2):1–8, 1998.
  * [22] J. Noorbakhsh, D. Schwab, A. Sgro, T. Gregor, and P. Mehta.  Modeling oscillations and spiral waves in Dictyostelium populations.  Phys. Rev. E., 91:062711, 2015.
  * [23] Nykamp D. Q.  The idea of synchrony of phase oscillators.  From Math Insight: http://mathinsight.org/synchrony_phase_oscillator_idea.
  * [24] E. M. Rauch and M. Millonas.  The role of trans-membrane signal transduction in turing-type cellular pattern formation.  J. Theor. Biol., 226:401–407, 2004.
  * [25] D. J. Schwab, A. Baetica, and P. Mehta.  Dynamical quorum-sensing in oscillators coupled through an external medium.  Physica D, 241(21):1782–1788, 2012.
  * [26] E. E. Sel’Kov.  Self-oscillations in glycolysis 1. A simple kinetic model.  European Journal of Biochemistry, 4(1):79–86, 1968.
  * [27] L. F. Shampine and M. W. Reichelt.  The Matlab ODE suite.  SIAM journal on scientific computing, 18(1):1–22, 1997.
  * [28] A. F. Taylor, M. Tinsley, and K. Showalter.  Insights into collective cell behavior from populations of coupled chemical oscillators.  Phys. Chemistry Chem Phys., 17(31):20047–20055, 2015.
  * [29] A. F. Taylor, M. Tinsley, F. Wang, Z. Huang, and K. Showalter.  Dynamical quorum sensing and synchronization in large populations of chemical oscillators.  Science, 323(5914):614–6017, 2009.
  * [30] M. R. Tinsley, A. F. Taylor, Z. Huang, and K. Showalter.  Emergence of collective behavior in groups of excitable catalyst-loaded particles: Spatiotemporal dynamical quorum sensing.  Phys. Rev. Lett., 102:158301, 2009.
  * [31] M. R. Tinsley, A. F. Taylor, Z. Huang, F. Wang, and K. Showalter.  Dynamical quorum sensing and synchronization in collections of excitable and oscillatory catalytic particles.  Physica D, 239(11):785–790, 2010.
  * [32] H. Uecke, J. Müller, and B. A. Hense.  Individual-based model for quorum sensing with background flow.  Bulletin of mathematical biology, 76(7):1727–1746, 2014.
  * [33] J. P. Ward, J. R. King, A. J. Koerber, P. Williams, J. M. Croft, and R. E. Sockett.  Mathematical modelling of quorum sensing in bacteria.  Mathematical Medicine and Biology, 18(3):263–292, 2001.
  * [34] J. Wolf and R. Heinrich.  Effect of cellular interaction on glycolytic oscillations in yeast: a theoretical investigation.  Biochem J., 345:321–334, 2000.