Synchrony and Oscillatory Dynamics for a 2-D PDE-ODE Model of Diffusion-Sensing with Small Signaling Compartments
SSynchrony and Oscillatory Dynamics for a 2-D PDE-ODE Model ofDiffusion-Sensing with Small Signaling Compartments
Sarafa A. Iyaniwura and Michael J. WardJuly 20, 2020
Abstract
We analyze a class of cell-bulk coupled PDE-ODE models, motivated by quorum and diffusion sensing phenomena inmicrobial systems, that characterize communication between localized spatially segregated dynamically active signalingcompartments or “cells” that have a permeable boundary. In this model, the cells are disks of a common radius ε (cid:28) andthey are spatially coupled through a passive extracellular bulk diffusion field with diffusivity D in a bounded 2-D domain.Each cell secretes a signaling chemical into the bulk region at a constant rate and receives a feedback of the bulk chemicalfrom the entire collection of cells. This global feedback, which activates signaling pathways within the cells, modifiesthe intracellular dynamics according to the external environment. The cell secretion and global feedback are regulated bypermeability parameters across the cell membrane. For arbitrary reaction-kinetics within each cell, the method of matchedasymptotic expansions is used in the limit ε (cid:28) of small cell radius to construct steady-state solutions of the PDE-ODEmodel, and to derive a globally coupled nonlinear matrix eigenvalue problem (GCEP) that characterizes the linear stabilityproperties of the steady-states. The analysis and computation of the nullspace of the GCEP as parameters are variedis central to the linear stability analysis. In the limit of large bulk diffusivity D = D /ν (cid:29) , where ν ≡ − / log ε , anasymptotic analysis of the PDE-ODE model leads to a limiting ODE system for the spatial average of the concentrationin the bulk region that is coupled to the intracellular dynamics within the cells. Results from the linear stability theoryand ODE dynamics are illustrated for Sel’kov reaction-kinetics, where the kinetic parameters are chosen so that each cellis quiescent when uncoupled from the bulk medium. For various specific spatial configurations of cells, the linear stabilitytheory is used to construct phase diagrams in parameter space characterizing where a switch-like emergence of intracellularoscillations can occur through a Hopf bifurcation. The effect of the membrane permeability parameters, the reaction-kineticparameters, the bulk diffusivity, and the spatial configuration of cells on both the emergence and synchronization of theoscillatory intracellular dynamics, as mediated by the bulk diffusion field, is analyzed in detail. The linear stability theoryis validated from full numerical simulations of the PDE-ODE system, and from the reduced ODE model when D is large. Key Words:
Green’s function, bulk diffusion, globally coupled eigenvalue problem (GCEP), Hopf bifurcation, cell-bulkcoupling, synchronous oscillations, diffusion-sensing.
Cell-to-cell communication is an important aspect of microbial systems that is often achieved through the diffusion of anextracellular signaling molecule, referred to as an autoinducer, in their common environment (cf. [13], [45], [9], [55]).This form of bulk-mediated communication involves the secretion and feedback of signaling molecules from and into thecells, respectively, which enables each cell to adjust its intracellular dynamics based on the signals it receives from theautoinducer field that depends on the entire collection of cells. Specific autoinducers responsible for such an intercellularcommunication have been identified in many biological systems including, cyclic adenosine monophosphate (cAMP) thattriggers intracellular oscillations for a collection of social amoebae
Dictyostelium discoideum and guides them to aggregationin low nutrient environments (cf. [22], [37], [16]), acetaldehyde (Ace) that leads to glycolytic oscillations in a colony of yeastcells (cf. [8], [9], [10]), and acylated homoserine lactones (AHLs) that induces bioluminescence for certain species of squiddue to colonies of the marine bacterium
Vibrio fischeri located in the squid’s light organ (cf. [45]).In this context of intercellular bulk-mediated communication, quorum-sensing (QS) refers to the onset of collective dynamicsin the cells that occurs when the cell density increases past a threshold. There are two main categories of QS systems, forwhich mathematical models have been developed. The first main type, which includes yeast cells and social amoeba, involvesa switch-like onset of synchronized oscillatory intracellular dynamics as the cell density increases (cf. [8], [9], [10], [22], [33],[27]). On the macroscale, triggered synchronous temporal oscillations also occur in physicochemical systems involving smallcatalyst-loaded particles immersed in a BZ chemical mixture (cf. [47], [46], [49], [48]). The second main type of QS system,for which the marine bacterium
Vibrio fischeri and the human pathogen
Pseudomonas aeruginosa are prototypical examples,is where an increase in the cell density leads to a sudden transition between bistable steady-states (cf. [12], [52], [32]).QS systems that involve a switch-like onset of synchronous intracellular oscillations have most typically been studied forwell-stirred systems, where the bulk diffusivity is taken to be infinite. This well-mixed limit leads to a globally coupled ODE1 a r X i v : . [ n li n . PS ] J u l ystem, where the global coupling arises from the spatially homogeneous bulk diffusion field (cf. [22], [10], [46], [49], [30],[39], [26]). In this more classical ODE setting, well-established mathematical tools such as ODE bifurcation theory, phasereduction methods and the Kuramoto order parameter (cf. [41], [36]) can be used to analyze the onset of QS behavior andpredict the degree of synchronization of intracellular oscillations as the cell density increases.However, when the bulk diffusivity is finite, diffusion-sensing behavior associated with the spatial configuration of cells,the spatial gradients of the bulk signal, reflecting versus absorbing boundary conditions, and the mass transport propertiesof the bulk medium, have all been shown experimentally to play an important role for some QS systems (cf. [15], [11], [50],[24] and the references therein). In contrast to the study of QS behavior through an ODE theoretical framework, therehave been relatively few theoretical and modeling studies of how spatial diffusion triggers the onset of collective intracellularoscillations in spatially segregated, but localized, dynamically active reaction sites, and these studies have typically beenconsidered in 1-D spatial contexts (cf. [17], [7], [19], [21], [18], [31], [38], [57], [40].)The goal of this paper is to study the emergence and synchronization of intracellular oscillations for the coupled cell-bulkPDE-ODE model of [20] in a 2-D domain. The formulation of this model was inspired by the 3-D PDE-ODE cell-basedmodel of [34] with a single intracellular species (see also [35] and [51]). In contrast to the analysis in [20] that was restrictedeither to the well-mixed limit or to the simple case of a ring pattern of identical cells, our study will focus on analyzingdiffusion-sensing behavior, resulting from a passive bulk diffusion field with finite diffusivity, for various spatial configurationsof a collection of heterogeneous cells. Y Figure 1: A schematic diagram showing dynamically active signaling compartments (in cyan) in an arbitrary bounded 2D domain. Thegreen and red dot represent the signaling chemicals in the cells, where only the red is secreted into the extracellular bulk region. Onthe right: A zoomed-in illustration of the intracellular concentration of chemicals within each signaling compartment, the secretion ofsignaling molecules into the bulk region, and the feedback of chemical into the cells.
The coupled PDE-ODE model of [20] is formulated as follows: Within Ω we assume that there are m dynamically activecircular signalling compartments or “cells” of a common radius R , denoted by Ω j and centered at XXX j ∈ Ω for j = 1 , . . . , m .In the bulk, or extracellular, region Ω \ ∪ mj =1 Ω j , the concentration U ( XXX, T ) of the autoinducer or bulk signal, which isconfined within ∂ Ω , is assumed to satisfy the passive diffusion equation U T = D B ∆ U − k B U , T > , XXX ∈ Ω \ ∪ mj =1 Ω j ; (1.0.1a) ∂ n XXX U = 0 , XXX ∈ ∂ Ω ; D B ∂ n XXX U = β j U − β j µ j , XXX ∈ ∂ Ω j , j = 1 , . . . , m . (1.0.1b)Here D B > and k B > are the dimensional bulk diffusivity and rate of degradation of the bulk signal, respectively. Inthe Robin boundary condition (1.0.1b) on the cell membrane, β j > and β j > are dimensional parameters modeling theinflux and efflux of chemical into and out of the j th cell, while ∂ n XXX denotes the outer normal derivative on the cell boundarythat points into the bulk region. Within each cell we assume that there are n interacting species represented by the vector µµµ j = ( µ j , . . . , µ nj ) T . Assuming that the cells are sufficiently small so that there are no spatial chemical gradients within them,the intracellular reaction-kinetics FFF j for the j th cell is coupled to the bulk medium via the integration of the flux across thecell membrane. In this way, the intracellular dynamics within the j th cell is coupled to the bulk signal (1.0.1) by dµµµ j dT = k R µ c FFF j ( µµµ j /µ c ) + eee (cid:90) ∂ Ω j (cid:0) β j U − β j µ j (cid:1) d S XXX , j = 1 , . . . , m . (1.0.1c)Here eee = (1 , , . . . , T , k R > is the dimensional reaction rate for the intracellular kinetics, and µ c > is a typical valuefor µµµ j . In this model, one signaling chemical, labeled by µ j can permeate the cell membrane with an efflux parameter β j and, by diffusion through the bulk medium, can communicate with spatially distant cells. The influx permeability parameter β j controls the global feedback into the j th cell from the bulk diffusion field, which is determined by the entire collection ofcells. In Fig. 1 we schematically illustrate the cell-bulk coupling for the case of n = 2 intracellular species.We assume that the radius R of the signaling compartments is small relative to the domain length-scale L , and so weintroduce a small parameter ε ≡ R /L (cid:28) . Then, by non-dimensionalizing the coupled PDE-ODE model (1.0.1) using the2pproach of [20], we obtain that the dimensionless concentration of chemical U ( xxx, t ) in the bulk region satisfies τ ∂U∂t = D ∆ U − U , t > , xxx ∈ Ω \ ∪ mj =1 Ω ε j ; (1.0.2a) ∂ n U = 0 , xxx ∈ ∂ Ω ; εD ∂ n U = d j U − d j u j , xxx ∈ ∂ Ω ε j , j = 1 , . . . , m , (1.0.2b)which is coupled to the dimensionless dynamics within the j th cell byd uuu j d t = FFF j ( uuu j ) + eee ετ (cid:90) ∂ Ω εj ( d j U − d j u j ) d s , j = 1 , . . . , m , (1.0.2c)where uuu j = ( u j , . . . , u nj ) T is the dimensionless vector representing the n chemical species in the j th cell, labeled by Ω ε j ≡{ xxx | | xxx − xxx j | ≤ ε } . We assume that the centers of the cells are well-separated in the sense that dist ( xxx j , xxx k ) = O (1) for j (cid:54) = k and dist ( xxx j , ∂ Ω) = O (1) as ε → . In this dimensionless formulation (1.0.2), the key dimensionless parameters are D ≡ D B k B L , d j ≡ ε β j k B L = O (1) , d j ≡ ε β j Lk B = O (1) , τ ≡ k R k B . (1.0.3)We refer to D and τ as the effective bulk diffusivity and reaction-time parameter, respectively. In (1.0.3), the permeabilityparameters β j and β j are chosen as O ( ε − ) in order to ensure that there is an O (1) transport across the membrane of thesmall cells. The parameter τ measures the relative rate of the intracellular dynamics to the time-scale for degradation ofthe bulk chemical. When the intracellular reactions proceed slowly, τ is small, and little communication between the cellsoccurs. When the effective bulk diffusivity D is large, the cells are readily able to communicate through the bulk medium,and in the well-mixed limit D → ∞ the bulk signal becomes spatially homogeneous. Alternatively, for smaller values of D ,only those cells that are in close spatial proximity should be able to communicate through the bulk diffusion field.For arbitrary intracellular reaction kinetics and for an arbitrary spatial arrangement of cells, in §2 we use strong localizedperturbation theory (cf. [53], [54]) in the limit ε → to construct steady-state solutions of (1.0.2) and to derive the linearstability problem for these steady-states. Unstable eigenvalues of the linearization of a steady-state are shown to correspondto roots λ in Re ( λ ) > for which a certain globally coupled nonlinear matrix eigenvalue problem (GCEP) M ( λ ) ccc = 0 (see(2.2.13)) has a nontrivial solution. The m × m matrix M , which couples all the cells through an eigenvalue-dependent Green’smatrix, depends on the dimensionless parameters in (1.0.3). The components of the corresponding normalized eigenvector ccc determines the relative magnitude of the spatial gradient of the bulk signal near the cell membranes and it determines therelative phases and amplitudes of small-scale oscillations within the cells at the onset of a Hopf bifurcation when λ = iλ I ,with λ I > , is a root of det M ( λ ) = 0 (see (2.2.17)).In §3, strong localized perturbation theory is used to reduce the PDE-ODE model (1.0.2) to an ODE differential algebraicsystem with global coupling for the distinguished limit D = D /ν (cid:29) of large bulk diffusivity, where ν ≡ − / log ε . Incontrast to the simpler ODE system derived in [20] for the well-mixed limit D → ∞ (i.e. D (cid:29) O ( ν − ) ), the new derivationin §3, as summarized in Proposition 2, leads to an ODE system that depends on the scaled bulk diffusivity parameter D and it depends weakly on the specific spatial configuration of cells within the domain.Our asymptotic theory is applied for the case of Sel’kov reaction kinetics, which has been used as a conceptual model forglycolysis oscillations in yeast cells [44]. With Sel’kov kinetics, (1.0.2) has a unique steady-state when ε (cid:28) . As indicatedfrom the experimental and modeling studies of collective behavior in yeast cells (cf. [9], [8], [10]), individual yeast cellsare typically non-oscillatory when isolated, but readily become synchronized in a population of such cells. In qualitativeagreement with this observation, the Sel’kov parameters are chosen to be close to the threshold for the onset of limit-cycleoscillations for an isolated cell. As a result, the switch-like emergence of intracellular oscillations, resulting from a Hopfbifurcation and illustrated through various phase diagrams, is inherently due to the cell-cell interaction, as mediated by thebulk diffusion field. In our study we will analyze how the onset of intracellular oscillations and synchronization depends onthe membrane permeability parameters, the reaction-time parameter, the bulk diffusivity, a Sel’kov kinetic parameter, andthe spatial pattern of cells. Diffusion-sensing behavior, whereby cells initiate oscillations as a result of certain spatial effectssuch as cell-clustering or the buildup of large spatial gradients of the autoinducer field near the cell boundary, are illustrated.In §4 we illustrate our theory with Sel’kov reaction kinetics for a ring and center-cell pattern of cells in the unit disk (seeFig. 3), where the ring cells are taken to have common parameters but where the parameters for the center cell can be different.For the unit disk, the Green’s matrices needed in the steady-state and linear stability theory are available analytically. Fora ring and center-cell pattern, the GCEP matrix M ( λ ) has a cyclic sub-block and so in §4.1 we can analytically identifyparticular spatial modes for which det M ( λ ) = 0 . By using arclength continuation in D , Hopf bifurcation boundaries in the ( D, τ ) parameter space for which det M ( iλ I ) = 0 can be computed for each of these modes. In open regions of the ( D, τ ) parameter plane, we show how to use a winding number criterion numerically on the roots of det M ( λ ) = 0 in Re ( λ ) > ,so as to compute the number of unstable eigenvalues of the linearization of the unique steady-state. In §4.2 these phasediagrams are shown for the case of two ring cells. The triggering effect on the emergence of intracellular oscillations dueto a center cell with either different permeability parameters or a different Sel’kov kinetic parameter is studied in §4.3 and§4.4, respectively. Diffusion-sensing behavior as a result of changes in the ring radius are studied in §4.5. The linear stability3heory, which predicts the onset of intracellular oscillations together with the amplitude and phase differences near the Hopfbifurcation point, is validated through large-scale simulations of the PDE-ODE model (1.0.2) using FlexPDE [14]. Moreover,we show that in certain cases the new ODE system in Proposition 2, derived for D = D /ν (cid:29) , can still provide a decentagreement with the full numerical results even when D = O (1) .In §5 we study how the onset of intracellular oscillations depends on the specific spatial arrangement of ten cells in the unitdisk that differ only in their influx permeabilites d j , for j = 1 , . . . , m . The cell configurations considered include two clustersof cells (see Fig. 17b), two rings of cells with two isolated cells (see Fig. 23b), and arbitrarily placed cells (see Fig. 17d). Forcomputational simplicity, we primarily focus on the regime D = D /ν (cid:29) , where matrix perturbation theory can be used toasymptotically calculate the spectrum of the GCEP matrix M ( λ ) , as summarized in Proposition 4. Since this analysis showsthat only one matrix eigenvalue of M ( λ ) can cross through zero, in §5.1 a numerical root-finding is readily implemented onthis eigenvalue to determine Hopf bifurcation boundaries in the ( D , τ ) plane where intracellular oscillations originate forthe various spatial configurations of cells and permeability parameter sets d j for j = 1 , . . . , . Our linear stability results,validated through FlexPDE simulations of (1.0.2) and the ODE dynamics of Proposition 2, shows that when D = D /ν the intracellular dynamics depend sensitively on the influx permeability set, but only weakly on the cell locations. In §5.1.1we implement the linear stability theory based on the GCEP matrix for D = O (1) to predict that small-scale intracellularoscillations can be highly heterogeneous in terms of amplitude and phase when there are isolated cells. The linear stabilitytheory is confirmed from full PDE simulations. Finally, in §6 we briefly summarize our study and discuss some biologicalmodeling problems that are well-aligned with the cell-based PDE-ODE framework of (1.0.1). In this section, strong localized perturbation theory in the limit ε → is applied to the dimensionless PDE-ODE model(1.0.2) for the regime D = O (1) . This theory is used to asymptotically approximate the steady-state solution of the coupledmodel, and also to formulate a globally coupled eigenvalue problem (GCEP) for studying the linear stability properties ofthe derived steady-state solution. We construct the steady-state solution for (1.0.2) using matched asymptotics. In the j th inner region, defined within an O ( ε ) neighborhood of the boundary of the j th cell, we introduce the local variables yyy = ε − ( xxx − xxx j ) and U j ( xxx ) = U j ( εyyy + x j x j x j ) ,where ρ ≡ | yyy | . Upon writing (1.0.2) in terms of the inner variables, for ε → the steady-state problem in the j th inner regionis ∆ U j = 0 for ρ ≥ , subject to D ∂ ρ U j = d j U j − d j u j on ρ = 1 . The radially symmetric solution to this problem is U j ( ρ ) = A j log ρ + 1 d j (cid:0) D A j + d j u j (cid:1) , j = 1 , . . . , m , (2.1.1)where A j for j = 1 , . . . , m are constants to be determined. Upon substituting (2.1.1) into the steady-state problem of (1.0.2c),we obtain that the steady-state intracellular dynamics uuu j of the j th cell satisfies FFF j ( uuu j ) + 2 πDτ A j eee = 000 , j = 1 , . . . , m . (2.1.2)This determines uuu j in terms of the unknown constant A j . To proceed, we must derive another algebraic system for theconstants A j for j = 1 , . . . , m , which is then coupled to (2.1.2).By matching the far-field behaviour of the inner solution (2.1.1) to an outer steady-state solution for (1.0.2a), we obtain thatthe outer solution must satisfy a specific singularity behaviour as xxx → xxx j . In this way, the steady-state outer approximationfor the bulk solution satisfies ∆ U − ϕ U = 0 , xxx ∈ Ω \ { xxx , . . . , xxx m } ; ∂ n U = 0 , xxx ∈ ∂ Ω ; U ∼ A j log | xxx − xxx j | + A j ν + 1 d j ( DA j + d j u j ) , as xxx → xxx j , j = 1 , . . . , m , (2.1.3)where ν ≡ − / log ε , with ε (cid:28) , and ϕ ≡ (cid:112) /D . The pre-specification of the regular part of each singularity structure in(2.1.3) yields a constraint. Overall these constraints provide an algebraic system for A j for j = 1 , . . . , m .To determine this algebraic system, we next represent the solution to (2.1.3) as U = − π m (cid:88) i =1 A i G ( xxx ; xxx i ) , (2.1.4)4here the reduced-wave Green’s function G ( xxx ; xxx j ) satisfies ∆ G − ϕ G = − δ ( xxx − xxx j ) , xxx ∈ Ω ; ∂ n G = 0 , xxx ∈ ∂ Ω ; G ( xxx ; xxx j ) ∼ − π log | xxx − xxx j | + R ( xxx j ) + o (1) , as xxx → xxx j . (2.1.5)Here R j ≡ R ( xxx j ) is the regular part of G ( xxx ; xxx j ) at xxx = xxx j . By expanding (2.1.4) as xxx → xxx j , we simply require that thenon-singular terms of the resulting expression agree with that specified in (2.1.3) for each j = 1 , . . . , m . This leads to analgebraic system for the vector AAA ≡ ( A , . . . , A m ) T , which is given in matrix form as (cid:16) I + 2 πν G + νDP (cid:17) AAA = − ν P uuu , (2.1.6)where uuu ≡ ( u , u , . . . , u m ) T denotes the vector of chemicals that is secreted into the bulk region by the cells. In (2.1.6), G is the symmetric reduced-wave Green’s interaction matrix, while P and P are m × m diagonal matrices, defined by G ≡ R G . . . G m G R . . . G m ... ... . . . ... G m G m . . . R m , P ≡ diag (cid:16) d , . . . , d m (cid:17) , P ≡ diag (cid:16) d d , . . . , d m d m (cid:17) . (2.1.7)Here G ji = G ij ≡ G ( xxx j ; xxx i ) for i (cid:54) = j , and R j ≡ R ( xxx j ) for j = 1 , . . . , m , are obtained from the solution to (2.1.5).Overall, the asymptotic steady-state solution is determined in terms of the solution AAA and uuu j for j = 1 , . . . , m , to the n × m dimensional nonlinear algebraic system (NAS) given by (2.1.2) and (2.1.6). This system applies to arbitrary localreaction kinetics FFF j and permeabilities parameters d j > and d j > for j = 1 , . . . , m . When the kinetics and permeabilityparameters are identical for all the cells, the NAS reduces to the system given in equations (2 . and (2 . of [20].Depending on the specific reaction kinetics assumed, the solution structure to (2.1.2) and (2.1.6) as parameters are variedcan involve solution multiplicity, saddle-node points, and other bifurcations. However, to illustrate our asymptotic theorywe will focus on the two-component Sel’kov reaction kinetics for which (2.1.2) and (2.1.6) has a unique solution. In the previous subsection, we characterized steady-state solutions of the coupled model (1.0.2) using strong perturbationtheory. Suppose that the NAS (2.1.2) and (2.1.6) has a solution for a given set of parameters. This then yields an approxi-mation to the steady-state solution solution U e ( xxx ) and uuu ej , for j = 1 , . . . , m , to (1.0.2). To determine the linear stability ofthis steady-state we begin by introducing the perturbation U ( xxx, t ) = U e ( xxx ) + e λt ξ ( xxx ) and uuu j ( t ) = uuu ej + e λt φφφ j , j = 1 , . . . , m . (2.2.1)where λ is the eigenvalue of the linearization, and ξ ( xxx ) and φφφ j ≡ ( φ j , . . . , φ nj ) T are the corresponding eigenfunctions in thebulk region and in the j th cell, respectively. Upon substituting this perturbation into the PDE-ODE model (1.0.2), in thebulk region we obtain the linearized problem τ λξ = D ∆ ξ − ξ , xxx ∈ Ω \ ∪ mj =1 Ω ε j ; (2.2.2a) ∂ n ξ = 0 , xxx ∈ ∂ Ω ; εD ∂ nj ξ = d j ξ − d j φ j , xxx ∈ ∂ Ω ε j , j = 1 , . . . , m , (2.2.2b)which is coupled to the linearized intracellular dynamics of the j th cell given in terms of φφφ j ≡ ( φ j , . . . , φ nj ) T , by λφφφ j = J j φφφ j + eee ετ (cid:90) ∂ Ω εj ( d j ξ − d j φ j ) d s , j = 1 , . . . , m . (2.2.2c)Here J j ≡ J ( uuu ej ) is the Jacobian matrix of the local kinetics FFF j evaluated at the steady-state uuu ej .Next, we use strong localized perturbation theory to analyze the eigenvalue problem (2.2.2) in the limit ε → . Thisanalysis leads to a limiting globally coupled eigenvalue problem (GCEP) for λ in the form of a nonlinear matrix eigenvalueproblem. This GCEP will be used to investigate instabilities of the steady-state solution for the PDE-ODE system (1.0.2).To derive this GCEP, we first construct an inner region in an O ( ε ) neighborhood of the j th cell by introducing the localvariables yyy = ε − ( xxx − xxx j ) and ξ j ( xxx ) ≡ ξ j ( xxx j + εyyy ) with ρ = | yyy | . From (2.2.2), we obtain for ε → that ∆ ξ j = 0 , < ρ < ∞ ; D ∂ ρ ξ j = d j ξ j − d j φ j , on ρ = 1 , (2.2.3)5n the j th inner region. The radially symmetric solution to (2.2.3) is ξ j = c j log ρ + 1 d j (cid:0) D c j + d j φ j (cid:1) , j = 1 , . . . , m , (2.2.4)where ρ = | yyy | and c j for j = 1 , . . . , m are constants to be determined. Upon substituting (2.2.4) into the linearizedintracellular dynamics of the j th cell (2.2.2c), we obtain a linear relation between φφφ j and c j given by ( J j − λI ) φφφ j = − πDτ c j eee . j = 1 , . . . , m . (2.2.5)Next, by analyzing the outer solution in the bulk region, we will derive another linear system, which will be coupled to(2.2.5). These two systems will provide the GCEP that is needed to study the linear stability of the steady-state solution.To determine this additional linear system, we first match the far-field behaviour of the inner solution (2.2.4) to the outersolution in order to obtain the singularity behaviour of the outer solution as xxx → xxx j . This yields in the bulk region that ∆ ξ − ϕ λ ξ = 0 , xxx ∈ Ω \ { xxx , . . . , xxx m } ; ∂ n ξ = 0 , xxx ∈ ∂ Ω ; ξ ∼ c j log | xxx − xxx j | + c j ν + 1 d j ( Dc j + d j φ j ) , as xxx → xxx j , j = 1 , . . . , m , (2.2.6)where ν ≡ − / log ε and ϕ λ = (cid:112) (1 + τ λ ) /D . The solution to (2.2.6) is represented as ξ ( xxx ) = − π m (cid:88) i =1 c i G λ ( xxx ; xxx i ) , (2.2.7)where the eigenvalue-dependent reduced-wave Green’s function G λ ( xxx ; xxx j ) satisfies ∆ G λ − ϕ λ G λ = − δ ( xxx − xxx j ) , xxx ∈ Ω ; ∂ n G λ = 0 , xxx ∈ ∂ Ω ; (2.2.8a) G λ ( xxx ; xxx j ) ∼ − π log | xxx − xxx j | + R λ ( xxx j ) + o (1) , as xxx → xxx j . (2.2.8b)Here R λj ≡ R λ ( xxx j ) is the regular part of G λ ( xxx ; xxx j ) at xxx = xxx j . In (2.2.7) we have specified the principal branch of ϕ λ toensure that G λ is analytic in Re ( λ ) > and that G λ decays far away from the cells.By expanding (2.2.7) as xxx → xxx j , we equate the non-singular terms of the resulting expression with those specified in (2.2.6)for each j = 1 , . . . , m . This yields a linear system for ccc ≡ ( c , . . . , c m ) T given in matrix form by (cid:0) I + 2 πν G λ + ν DP (cid:1) ccc = − ν P φφφ , (2.2.9)where φφφ ≡ ( φ , . . . , φ m ) T and G λ is the eigenvalue-dependent reduced-wave Green’s matrix whose entries are defined by ( G λ ) ij = ( G λ ) ji ≡ G λ ( xxx j , xxx i ) for i (cid:54) = j , and ( G λ ) jj = R λj ≡ R λ ( xxx j ) . (2.2.10)In (2.2.9), the diagonal matrices P and P are defined in (2.1.7).Next, we will combine the algebraic system (2.2.9) with (2.2.5) in order to derive the GCEP for λ and ccc = ( c , . . . , c m ) T .We first use (2.2.5) to determine φ j in terms of the constant c j as φ j = 2 πDτ − eee T ( λI − J j ) − eee c j for j = 1 , . . . , m , providedthat λ is not an eigenvalue of J j for any j = 1 , . . . , m . In matrix form this yields φφφ = 2 πDτ K ccc , (2.2.11)where K ≡ K ( λ ) is an m × m diagonal matrix K ≡ diag ( K , . . . , K m ) , whose entries are given by K j = eee T ( λI − J j ) − eee = eee T N j eee det( λI − J j ) = ( N j ) det( λI − J j ) . (2.2.12a)Here N j is the n × n matrix of cofactors of ( λI − J j ) , while ( N j ) is its entry in the first row and first column given by ( N j ) ≡ ( N j ( λ )) = det λ − ∂F j ∂u (cid:12)(cid:12)(cid:12) uuu = uuu e,j . . . − ∂F j ∂u n (cid:12)(cid:12)(cid:12) uuu = uuu e,j ... . . . ... − ∂F nj ∂u (cid:12)(cid:12)(cid:12) uuu = uuu e,j . . . λ − ∂F nj ∂u n (cid:12)(cid:12)(cid:12) uuu = uuu e,j . (2.2.12b)6he functions F j , . . . , F nj are the components of the local reaction kinetics F j F j F j ≡ ( F j , . . . , F nj ) T of the j th cell.Upon substituting (2.2.11) into (2.2.9), we obtain the m -dimensional homogeneous algebraic system M ( λ ) ccc = 000 , where M ( λ ) ≡ I + 2 πν G λ + ν D P + 2 πνDτ P K , (2.2.13a)where ν ≡ − / log ε . The homogeneous system (2.2.13a), where M is a symmetric but non-Hermitian matrix, is referred toas the globally coupled eigenvalue problem (GCEP). The GCEP is a nonlinear matrix eigenvalue problem for λ , and it hasa nontrivial solution ccc (cid:54) = 000 if and only λ satisfies det M ( λ ) = 0 . We label the set Λ( M ) as the union of all such roots, i.e. Λ( M ) ≡ { λ | det M ( λ ) = 0 } . (2.2.13b)The parameters in the GCEP (2.2.13a) are the bulk diffusivity D , the reaction-timescale τ and the permeabilities d j and d j ,for j = 1 , . . . , m , which are encoded in the matrices P and P given in (2.1.7). Moreover, in (2.2.13a), the effect of the spatialconfiguration of the centers { xxx , . . . , xxx m } ∈ Ω of the cells and the domain shape Ω arises from both the eigenvalue-dependentreduced wave Green’s interaction matrix G λ and the steady-state solution, which determines K in (2.2.12).A recent survey of nonlinear matrix eigenvalue problems and available solution strategies for certain classes of matrices isgiven in [23] and [4]. A range of applications of such problems, but in simpler contexts where M ( λ ) is either a polynomialor rational function of λ , are discussed in [3].Any element λ ∈ Λ( M ) satisfying Re ( λ ) > provides an approximation, valid as ε → , for an unstable discrete eigenvalueof the linearized system (2.2.2). This leads to the following criterion regarding the linear stability of the steady-state solution. Proposition 1.
For ε → , a steady-state solution to (1.0.2) is linearly stable when there are no roots to det M ( λ ) = 0 inRe ( λ ) > , i.e. for all λ ∈ Λ( M ) we have Re ( λ ) < . Moreover, if A e and uuu ej for j = 1 , . . . , m is a non-degenerate solutionto the NAS (2.1.2) and (2.1.6), for which J j is non-singular, then λ = 0 is not a root of det M ( λ ) = 0 .Proof. For ε → , any discrete eigenvalue λ of the linearization (2.2.1) corresponds to a non-trivial solution to (2.2.2). Since,for ε → , these discrete eigenvalues comprise the set Λ( M ) in (2.2.13b), the steady-state is linearly stable if all λ ∈ Λ( M ) satisfy Re ( λ ) < . Next, suppose that A e and uuu ej for j = 1 , . . . , m is a non-degenerate solution to the NAS (2.1.2) and(2.1.6). Introducing the perturbation A = A e + ψψψ and uuu j = uuu ej + vvv j , where | vvv j | (cid:28) and | ψψψ | (cid:28) , we linearize (2.1.2) and(2.1.6) to obtain (cid:0) I + 2 πν G + ν DP (cid:1) ψψψ = − ν P vvv , J j vvv j = − πDτ ψ j eee . j = 1 , . . . , m , (2.2.14)where vvv ≡ ( v , . . . , v m ) T , ψψψ = ( ψ , . . . , ψ m ) T and vvv j = ( v j , . . . , v nj ) T . Assuming that J j is invertible, (2.2.14) yields that J ψψψ = 000 , where J ≡ I + 2 πν G + ν D P + 2 πνDτ P K and K ≡ − diag (cid:0) eee T J − eee , . . . , eee T J − m eee (cid:1) . (2.2.15)We conclude that det J (cid:54) = 0 owing to the fact that A e and uuu ej for j = 1 , . . . , m is assumed to be a non-degenerate solutionof the NAS (2.1.2) and (2.1.6). Finally, since it is readily verified that J = M (0) , where M ( λ ) is the GCEP matrix in(2.2.13a), we conclude that λ = 0 is not a root of (2.2.13b).Proposition 1 implies that branches of non-degenerate solutions to the NAS (2.1.2) and (2.1.6) obtained as a parameteris varied cannot lose stability through a zero-eigenvalue crossing of the GCEP (2.2.13a). This simple observation is themotivation for analzying whether stability can be lost through Hopf bifurcations associated with the linearization.In terms of the eigenvectors ccc and eigenvalues λ ∈ Λ( M ) of the GCEP (2.2.13), we obtain from (2.2.1), (2.2.7) and (2.2.5)that the linearization around the bulk and intracellular steady-state solutions is, for ε → , given by the superposition U ( xxx, t ) ∼ U e − (cid:88) λ ∈ Λ( M ) e λt m (cid:88) j =1 c j G λ ( xxx ; xxx j ) ; uuu j ( t ) ∼ uuu ej + 2 πDτ (cid:88) λ ∈ Λ( M ) e λt c j ( λI − J j ) − eee , j = 1 , . . . , m , (2.2.16)where each ccc = ( c , . . . , c m ) T depends on the particular eigenvalue λ . To relate the diffusive flux into the j -th cell to thecomponents of ccc we use the the inner solutions (2.1.1) and (2.2.4). To determine the effect on the intracellular componentthat can be transported across the membrane we calculate eee T uuu j ≡ u j in (2.2.16). This yields that D∂ ρ U | ρ =1 ∼ D A j + (cid:88) λ ∈ Λ( M ) c j e λt , u j ∼ u ej + 2 πDτ (cid:88) λ ∈ Λ( M ) ( K ccc ) j e λt , j = 1 , . . . , m , (2.2.17)where K ≡ diag ( K , . . . , K m ) , with K j defined in (2.2.12). As evident from (2.2.17), and discussed for various examples in§4.2 and §5.1, the modulus | ( K ccc ) j | and argument arg ( K ccc ) j of the components of a complex-valued K ccc , resulting from pure7maginary eigenvalues with Re ( λ ) = 0 and Im ( λ ) (cid:54) = 0 , determines the relative amplitudes and phase differences of small scaleintracellular oscillations near a Hopf bifurcation of the steady-state solution.Our linear stability theory for steady-state solutions of PDE-ODE model (1.0.2) can be applied to any configuration ofcells in an arbitrary 2-D bounded domain and for arbitrary local reaction kinetics. However, in our illustrations of the theorybelow in §4 and §5, we will consider the two-component Sel’kov model, which is used in simple models of glycolysis (cf. [37],[44]). From (1.0.2c), for an isolated cell with no influx from the bulk, the intracellular dynamics within the j th cell thataccounts for efflux from the cell boundary is given by d uuu j / d t = FFF j ( uuu j ) − πd j u j eee /τ , where uuu j = ( u j , u j ) T and the Sel’kovkinetics FFF j ( v, w ) = ( F j ( v, w ) , F j ( v, w )) T are defined by F j ( v, w ) = α j w + wv − v , F j ( v, w ) = ζ j (cid:2) µ j − (cid:0) α j w + wv (cid:1)(cid:3) . (2.2.18)The steady-state for this isolated cell is given by u j = µ j /χ j and u j = µ j / (cid:0) α j + ( u j ) (cid:1) , where χ j ≡ πd j /τ . Thedeterminant and trace of the Jacobian J ej for this isolated cell with boundary efflux isdet ( J ej ) = ζ j χ j (cid:32) α j + µ j χ j (cid:33) , tr ( J ej ) = 2 µ j χ j (cid:32) α j + µ j χ j (cid:33) − − χ j − ζ j (cid:32) α j + µ j χ j (cid:33) . (2.2.19)Since det ( J ej ) > , the steady-state for this isolated cell is linearly stable only if tr ( J ej ) < . We will choose Sel’kov kineticparameters α j , µ j and ζ j so that an isolated cell with zero boundary efflux (i.e. d j = 0 ) is linearly stable, but with parametersrather close to the stability threshold. A set of such parameters is shown in Fig. 2a where we verify that tr ( J ej ) < on . < α j < . when µ j = 2 and ζ j = 0 . . The Hopf bifurcation (HB) boundary in the α j versus µ j parameter plane, asobtained by setting tr ( J ej ) = 0 is α j = − µ j χ j + 12 ζ j − χ j + (cid:115) χ j + 8 ζ j µ j χ j , where χ j ≡ πd j τ . (2.2.20)For an isolated cell, a simple application of the Poincare-Bendixson theorem shows that when the steady-state is unstablethe cell will have limit cycle oscillations. When there is no boundary efflux, i.e. d j = 0 , this parameter range of periodicsolutions is given by the green-shaded region in Fig. 2b. However, no time-periodic oscillations with Sel’kov kinetics arepossible when the steady-state is linearly stable. In Fig. 2c we show how the HB boundary for an isolated cell depends on theboundary efflux parameter d j . As expected, for the fixed value µ j = 2 , the interval in α j where oscillations are possible isdecreased when there is an efflux out of the cell boundary. The shifting of the HB boundaries to the right in Fig. 2c indicatesthat a greater rate µ j of production of u j is needed to ensure oscillations when there is a boundary efflux. The interval in µ j where oscillations are possible, at least for some range of α j > , is < µ j < χ / j / (cid:112) ζ j .For the baseline parameter set µ j = 2 and ζ j = 0 . , in §4 and §5 we will show that the inter-cell coupling via the bulkdiffusion field can be sufficient to trigger an oscillatory instability in the cells through a Hopf bifurcation. (a) tr ( J ej ) versus α j for isolated cell (b) Instability region: isolated cell (c) HB boundaries: isolated cell with effluxFigure 2: Left panel: tr ( J ej ) , from (2.2.19) versus α j , for the steady-state of the Sel’kov kinetics (2.2.18) for an isolated cell, with µ j = 2 and ζ j = 0 . . This steady-state is linearly stable but the parameters are close to the stability threshold. Middle panel: Green-shadedregion of instability where tr ( J ej ) > in the α j versus µ j plane for the steady-state of an isolated cell with no boundary efflux and ζ j = 0 . . Within this region, a time-periodic solution (limit cycle) occurs for an isolated cell. The HB boundary is given by (2.2.20)with d j = 0 . In the unshaded region the steady-state is linearly stable. Right panel: HB boundaries for an isolated cell (see (2.2.20))with ζ j = 0 . , τ = 0 . and four boundary efflux parameters. A larger production rate µ j is needed to support oscillations. For Sel’kov local reaction kinetics (2.2.18) it is readily shown that, with an arbitrary arrangement of cells, there is a uniquesolution to the NAS (2.1.2) and (2.1.6) given by u ej = µ j + 2 πDτ A j , u ej = µ j α j + (cid:0) u ej (cid:1) , j = 1 , . . . , m , (2.2.21a)8here A = ( A , . . . , A m ) T satisfies the linear algebraic system (cid:18) I + 2 πν G + ν D P + 2 πνDτ P (cid:19) A = − νP µµµ , with µµµ ≡ ( µ , . . . , µ m ) T . (2.2.21b)For ν (cid:28) sufficiently small, the matrix in (2.2.21b) is invertible, yielding a unique solution for A . For Sel’kov kinetics weconclude that steady-state solutions of the PDE-ODE model (1.0.2) are always non-degenerate as ε → . As such, since byProposition 1 stability cannot be lost via a zero-eigenvalue crossing, in §4 and §5 we will focus on analyzing instabilities ofthe steady-state arising from Hopf bifurcations. D = O ( ν − ) In this section, we use a singular perturbation approach to reduce the dimensionless coupled PDE-ODE model (1.0.2) to anODE system that is valid for the limiting regime D = O ( ν − ) (cid:29) O (1) , where ν ≡ − / log ε and ε (cid:28) . This ODE systemdepends weakly on the spatial configuration of the cells and on the scaled diffusivity D = O (1) , defined by D = D /ν .Consider a collection of m small cells centered at the points xxx , . . . , xxx m in a 2-D bounded domain Ω . Define Ω p ≡ ∪ mj =1 Ω ε j as the region formed by the union of the cells, and the average bulk concentration U = U ( t ) by U = 1 | Ω \ Ω p | (cid:90) Ω \ Ω p U d xxx . (3.0.1)Our goal is to derive an ODE for U ≡ U ( t ; ν ) , accurate to O ( ν ) , which is coupled to the intracellular dynamics of the cells,as given in (1.0.2b). Upon multiplying (1.0.2a) by / | Ω \ Ω p | and using the divergence theorem, we obtain τ U t + U = 2 π | Ω \ Ω p | m (cid:88) j =1 (cid:32) d j u j − d j πε (cid:90) ∂ Ω εj U d s (cid:33) , (3.0.2)where we used | ∂ Ω ε j | = 2 πε for the perimeter of each cell. Since | Ω p | = mπε , we estimate | Ω \ Ω p | = | Ω | − O ( ε ) , so that | Ω \ Ω p | → | Ω | as ε → in (3.0.2). Next, by evaluating the integral in (1.0.2b), and using | ∂ Ω ε j | = 2 πε , we obtaind uuu j d t = FFF j ( uuu j ) − πeee τ (cid:32) d j u j − d j πε (cid:90) ∂ Ω εj U d s (cid:33) , j = 1 , . . . , m . (3.0.3)In (3.0.2) and (3.0.3), we must estimate the bulk concentration U ≡ U ( xxx, t ) on the j th cell boundary ∂ Ω ε j .For D = O ( ν − ) we introduce the scaling D = D ν , where ν ≡ − ε , ε (cid:28) and D = O (1) . (3.0.4)Upon substituting (3.0.4) into (1.0.2), we obtain in the bulk region that τ U t = D ν ∆ U − U , t > , xxx ∈ Ω \ ∪ mj =1 Ω ε j ; ∂ n U = 0 , xxx ∈ ∂ Ω ε ; ε D ν ∂ n U = d j U − d j u j on ∂ Ω ε j , j = 1 , . . . , m . (3.0.5)In the inner region at an O ( ε ) neighborhood of the j th cell, we introduce the inner variables yyy = ε − ( xxx − xxx j ) and U ( xxx ) = V j ( xxx j + εyyy ; ν ) , with ρ = | yyy | . Writing (3.0.5) in terms of these inner variables, we obtain for ε → that ∆ V j = 0 , ρ ≥ ∂ ρ V j = νD (cid:16) d j V j − d j u j (cid:17) , on ρ = 1 , j = 1 , . . . , m , (3.0.6)which has the radially symmetric solution V j = ν b j log ρ + U j , where b j ≡ D (cid:16) d j U j − d j u j (cid:17) , j = 1 , . . . , m , (3.0.7)where V j | ρ =1 = U j is to be determined. By writing (3.0.7) in the outer xxx variable, and using | yyy | = ε − | xxx − xxx j | , we obtainthe following asymptotic matching condition for the outer solution in the bulk region: U ∼ ν b j log | xxx − xxx j | + (cid:16) d j D + 1 (cid:17) U j − d j D u j , as xxx → xxx j . (3.0.8)9ince V j ( ρ ) = U j on ρ = 1 from (3.0.7), we have (cid:82) ∂ Ω εj U d s = ε (cid:82) π V j | ρ =1 d θ = 2 πε U j . Then, from (3.0.2) and (3.0.3), andrecalling (3.0.7) for b j , we obtain for ε → that τ U t + U = 2 π | Ω | m (cid:88) j =1 (cid:0) d j u j − d j U j (cid:1) = − πD | Ω | m (cid:88) j =1 b j , d uuu j d t − FFF j ( uuu j ) = − πeee τ (cid:0) d j u j − d j U j (cid:1) = 2 πD τ eee b j , j = 1 , . . . , m . (3.0.9)To complete the derivation of the ODE system we must obtain an algebraic system for b j for j = 1 , . . . , m from the analysisof the outer solution. From (3.0.5) and (3.0.8), and relating U j to b j using (3.0.7), the outer problem for U ( xxx, t ) is τ U t = D ν ∆ U − U , t > , xxx ∈ Ω \ ∪ mj =1 Ω ε j ; ∂ n U = 0 , xxx ∈ Ω ; U ∼ νb j log | xxx − xxx j | + b j (cid:18) D d j (cid:19) + d j d j u j , as xxx → xxx j . (3.0.10)We then expand U ( xxx, t ) as U ( xxx, t ) = U + νD U ( xxx, t ) + . . . , where (cid:90) Ω U dxxx = 0 . (3.0.11)The zero average constraint on U ensures that U is the spatial average of U to terms of order O ( ν ) . Upon substituting(3.0.11) into (3.0.10), we obtain in the sense of distributions that U satisfies ∆ U = τ U t + U + 2 πD m (cid:88) i =1 b i δ ( xxx − xxx i ) , xxx ∈ Ω ; ∂ n U = 0 , xxx ∈ ∂ Ω ; U ∼ b j D log | xxx − xxx j | − D ν U + D ν (cid:20) b j (cid:16) D d j (cid:17) + d j d j u j (cid:21) , as xxx → xxx j . (3.0.12)The divergence theorem applied to (3.0.12) yields the ODE given in (3.0.9) while the linear system for b j , for j = 1 , . . . , m ,is obtained from the constraints involved with specifying the form of the regular part of the singularity behavior in (3.0.12).The solution to (3.0.12), with (cid:82) Ω U dxxx = 0 is written as U = − πD m (cid:88) i =1 b i G ( xxx ; xxx i ) , (3.0.13)where G ( xxx ; xxx j ) is the unique Neumann Green’s function satisfying ∆ G = 1 | Ω | − δ ( xxx − xxx j ) , xxx ∈ Ω ; ∂ n G = 0 , xxx ∈ ∂ Ω ; (3.0.14a) G ( xxx ; xxx j ) ∼ − π log | xxx − xxx j | + R j + o (1) , as xxx → xxx j , and (cid:90) Ω G d xxx = 0 . (3.0.14b)Here R j is the regular part of G at xxx = xxx j . By expanding (3.0.13) as xxx → xxx j , we enforce that the nonsingular part of theresulting expression agrees with that in (3.0.12). This yields that b j (cid:18) D d j (cid:19) + 2 πν b j R j + m (cid:88) i (cid:54) = j b i G ji = U − d j d j u j , j = 1 , . . . , m , (3.0.15)where G ji = G ( xxx j ; xxx i ) . This linear system for b , . . . , b m is then coupled to the ODEs given in (3.0.9). Upon writing thisODE system in matrix form we summarize the result as follows: Proposition 2.
Let ε → and assume that D = D /ν (cid:29) where D = O (1) and ν = − / log ε (cid:28) . Then, the PDE-ODEsystem (1.0.2) reduces to the following nm + 1 dimensional ODE system for U ≈ | Ω | − (cid:82) Ω U dxxx and the intracellular species:dd t U = − τ U − πD τ | Ω | eee T bbb ; d uuu j d t = FFF j ( uuu j ) + 2 πD eee τ b j , j = 1 , . . . , m , (3.0.16a) where eee ≡ (1 , . . . , T , eee ≡ (1 , , . . . , T and bbb ≡ ( b , . . . , b m ) T . In (3.0.16a), bbb is the solution to the linear system (cid:0) I + D P + 2 πν G (cid:1) bbb = U eee − P uuu , (3.0.16b)10 here uuu ≡ ( u , . . . , u m ) T and P and P are the diagonal matrices defined in terms of the permeabilities by (2.1.7). In(3.0.16b), G is the Neumann Green’s matrix with matrix entries ( G ) ij = ( G ) ji = G ( xxx i , xxx j ) , i (cid:54) = j and ( G ) jj = R j . (3.0.17) For ν (cid:28) , this ODE system is equivalent up to O ( ν ) terms todd t U = − τ U − πτ | Ω | (cid:2) Ueee T C eee − eee T C P uuu (cid:3) + O ( ν ) , d uuu j d t = FFF j ( uuu j ) + 2 πeee τ (cid:104) U ( C eee ) j − (cid:0) C P uuu (cid:1) j (cid:105) , j = 1 , . . . , m , (3.0.18a) where the matrix C is defined in terms of G and a diagonal matrix P by C ≡ P − πνD PG P , P ≡ diag (cid:18) D d d + D , . . . , D d m d m + D (cid:19) . (3.0.18b) In the well-mixed limit for which D → ∞ , (3.0.16b) yields D b j ∼ U d j − d j u j , so that (3.0.16a) reduces to U t = − τ U − πτ | Ω | m (cid:88) j =1 ( U d j − d j u j ) ; d uuu j d t = FFF j ( uuu j ) + 2 πeee τ ( U d j − d j u j ) , j = 1 , . . . , m . (3.0.19)To derive (3.0.18) from (3.0.16), we approximate the solution bbb to (3.0.16b) up to terms of order O ( ν ) . By inverting thediagonal matrix I + D P , we obtain from (3.0.16b) that bbb = 1 D (cid:18) I + 2 πνD PG (cid:19) − (cid:0) U P eee − P P uuu (cid:1) ∼ D (cid:18) I − πνD PG (cid:19) (cid:0) U P eee − P P uuu (cid:1) = 1 D (cid:0) U C eee − C P uuu (cid:1) + O ( ν ) , (3.0.20)where C and P are given in (3.0.18b). The ODE system (3.0.18) results from substituting (3.0.20) into (3.0.16).The ODE systems (3.0.16), or alternatively (3.0.18), for the regime D = O ( ν − ) are accurate up to and including termsof order O ( ν ) and show how the intracellular species are globally coupled through the spatial average of the bulk field.Since these ODE systems depend on the scaled diffusivity parameter D and include the effect of the spatial configuration xxx , . . . , xxx m of the cells through the Neumann Green’s matrix, these ODE systems can account for both diffusion-sensing andquorum-sensing behavior (see §4 and §5). In contrast, the limiting well-mixed ODE system (3.0.19), originally derived in[20] for the simpler case of identical cells, depends only on the number m of cells. As a result, the well-mixed ODE dynamicsis independent of the diffusivity and the spatial configuration of the cells.In our numerical experiments in §4 and §5 using the ODE system (3.0.16) the domain Ω is the unit disk. For the unitdisk, the Neumann Green’s function G ( xxx ; xxx j ) and its regular part R j , satisfying (3.0.14), are (see equation (4.3) of [28]) G ( xxx ; xxx j ) = − π log | xxx − xxx j | − π log (cid:0) | xxx | | xxx j | + 1 − xxx · xxx j (cid:1) + ( | xxx | + | xxx j | )4 π − π ,R j = − π log (cid:0) − | xxx j | (cid:1) + | xxx j | π − π . (3.0.21)For an arbitrary cell pattern { xxx , . . . , xxx m } , (3.0.21) is used to evaluate the Neumann Green’s matrix G as needed in (3.0.16). With Sel’kov reaction kinetics, we apply the theory developed in §2 to a ring and center cell configuration in the unit disk.This pattern is characterized by m − ≥ equally spaced cells on a concentric ring within the unit disk, and with one atthe center of the disk (see Fig. 3). For this pattern, the GCEP (2.2.13) will be used to obtain tractable nonlinear algebraicequations that can be solved numerically to compute HB boundaries in the τ versus D parameter plane. In addition, awinding number criterion is developed to count the number of unstable eigenvalues in open regions of this parameter plane.Some of our examples will show that rather small changes in either the permeabilities or reaction kinetic parameters of thecenter cell can significantly alter the region in parameter space where oscillations occurs. Consider a ring and center cell pattern of m cells where the cells on the ring of radius r have identical parameters, butwhere the center cell has possibly different permeabilities or Sel’kov kinetic parameters (see Fig. 3). The cell centers are at xxx j = r (cid:18) cos (cid:18) π ( j − m − (cid:19) , sin (cid:18) π ( j − m − (cid:19)(cid:19) , j = 1 , . . . , m − xxx m = , (4.1.1)11 igure 3: Schematic showing a ring and center cell configuration of m = 9 (left) and m = 3 (right) cells in the unit disk. The ringcells are identical and equally spaced on a concentric ring within the disk (in cyan). The center cell (in red), possibly has differentparameters. The red dots represent the signaling molecules secreted in the bulk region by the cells. where < r < . For this pattern, the reduced-wave Green’s matrix G in (2.1.7) can be partitioned as G = g m G m − ... g m g m . . . g m R m ; g m ≡ G ( xxx j , xxx m ) = G ( xxx m , xxx j ) , j = 1 . . . , m − , R m = R ( xxx m ) . (4.1.2)Here G m − is the ( m − × ( m − symmetric matrix block representing the interaction between the cells on the ring. Sincethis block is also cyclic it has the eigenpair G m − eee = ω eee , with eee = (1 , . . . , T ∈ R m − , and ω ≡ R + m − (cid:88) j =2 G ( xxx , xxx j ) . (4.1.3)In (4.1.2) there is a common interaction, represented by g m , between each ring cell and the center cell owing to the rotationalsymmetry and the fact that the ring cells are all equidistant from the center cell.For the identical ring cells, we label their permeabilities as d = d j and d = d j for j = 1 , . . . , m − and their commonSel’kov kinetic parameters as µ = µ j , α = α j and ζ = ζ j for j = 1 , . . . , m − . Since the unique steady-state solution to(2.2.21b) has the form A = ( A c , . . . , A c , A m ) T , we readily find that A c and A m satisfies the × linear system (cid:18) νDd + 2 πνω + 2 πνDτ d d (cid:19) A c + 2 πνg m A m = − νµ d d ; (4.1.4a) (cid:18) νDd m + 2 πνR m + 2 πνDτ d m d m (cid:19) A m + 2 πν g m ( m − A c = − νµ m d m d m , (4.1.4b)where ω is the eigenvalue of G m − in (4.1.3). In terms of A = ( A c , . . . , A c , A m ) T , the steady-state for the intracellularspecies as obtained from (2.2.21a) is u ej = (cid:40) u e ≡ µ + πDτ A c , j = 1 , . . . , m − ,u em ≡ µ m + πDτ A m , j = m , ; u ej = (cid:40) u e ≡ µ α + ( u e ) , j = 1 , . . . , m − ,u em ≡ µ m α m + ( u em ) , j = m . (4.1.5)Next, we determine the GCEP for the ring and center cell pattern using (2.2.13a). For this pattern, the GCEP matrix M ( λ ) in (2.2.13a) is written as M ( λ ) = 2 πν G λ + M , (4.1.6a)where G λ is the eigenvalue-dependent Green’s matrix, as defined in (2.2.10), and where the diagonal M is defined by M = M . . . M M m , where M ≡ νDd + πνDτ d d K c ,M m ≡ νDd m + πνDτ d m d m K m . (4.1.6b)Here K c = K c ( λ ) and K m = K m ( λ ) are the entries of the m × m diagonal matrix K = diag ( K c , . . . , K c , K m ) defined in(2.2.12a). For the Sel’kov kinetics given in (2.2.18), (2.2.12a) yields K c ≡ λ + det ( J ) λ − tr ( J ) λ + det ( J ) , K m ≡ λ + det ( J m ) λ − tr ( J m ) λ + det ( J m ) , (4.1.7a)12here the trace and determinant of the Jacobians of the intracellular dynamics for the identical ring cells and the center cellare given in terms of the steady-state values in (4.1.5) bydet ( J ) = ζ (cid:0) α + ( u e ) (cid:1) , det ( J m ) = ζ m (cid:0) α m + ( u em ) (cid:1) , tr ( J ) = (cid:104) u e ) µ − (cid:0) α + ( u e ) (cid:1) − ζ (cid:0) α + ( u e ) (cid:1) (cid:105) α + ( u e ) , tr ( J m ) = (cid:104) u em ) µ m − (cid:0) α m + ( u em ) (cid:1) − ζ m (cid:0) α m + ( u em ) (cid:1) (cid:105) α m + ( u em ) . (4.1.7b)From Proposition 1, discrete eigenvalues λ associated with the ring and center cell pattern are roots of det M ( λ ) = 0 . Aconvenient way to implement this determinant root-finding problem numerically is to use the special structure of M ( λ ) inorder to determine explicit formulae for its matrix spectrum M ( λ ) ccc j = σ j ccc j , for j = 1 , . . . , m , where σ j = σ j ( λ ) . Then, weneed only numerically solve the scalar root-finding problems σ j ( λ ) = 0 for λ for each j = 1 , . . . , m .To do so, we use the convenient fact that G λ can be partitioned, similar to that in (4.1.2), as G λ = g λm G λ ( m − ... g λm g λm . . . g λm R λm ; g λm ≡ G λ ( xxx j , xxx m ) = G λ ( xxx m , xxx j ) , j = 1 . . . , m − , R λm ≡ R λ ( xxx m ) , (4.1.8)where G λ is the eigenvalue-dependent reduced-wave Green’s function with regular part R λ satisfying (2.2.8). In (4.1.8), the ( m − × ( m − matrix block G λ ( m − , representing cell interactions on the ring, is symmetric and cyclic. As a result, ithas the well-defined eigenspace G λ ( m − vvv j = ω λj vvv j , j = 1 , . . . , m − eee T vvv j = 0 , j = 1 , . . . , m − , vvv m − = eee ≡ (1 , . . . , T ∈ R m − . (4.1.9)By using this special matrix structure, it readily follows that the GCEP matrix M ( λ ) in (4.1.6) admits m − anti-phasemodes , characterized by eigenvectors of the form ccc j = ( vvv j , T ∈ R m where vvv j ∈ R m − are those eigenvectors of G λ ( m − in(4.1.9), which satisfy eee T vvv j = 0 for j = 1 , . . . , m − . With this choice, (4.1.6) becomes πν G λ ( m − vvv j + M vvv j πνg λm eee T vvv j = σ j vvv j , for j = 1 , . . . , m − . (4.1.10)Since eee T vvv j = 0 for j = 1 , . . . , m − , we obtain from (4.1.10), (4.1.9) and (4.1.6b) that m − eigenpairs of M ( λ ) are σ j = 2 πνω λj + M = 2 πνω λj + 1 + νDd + 2 πνDτ d d K c , ccc j = ( vvv j , T , for j = 1 , . . . , m − , (4.1.11)where K c = K c ( λ ) is defined in (4.1.7a). We remark that since K c depends on the steady-state values A c and A m , asobtained from the linear system (4.1.4), this term depends on the permeabilities and local kinetics of the center cell. Discreteeigenvalues λ for the anti-phase modes are union of the zeroes of σ j ( λ ) = 0 for j = 1 , . . . , m − .For the remaining two eigenpairs of the GCEP matrix M ( λ ) the associated eigenvector ccc has the form ccc = ( eee, γ ) T forsome scalar γ to be determined and eee = (1 , . . . , T ∈ R m − . This eigenvector is referred to as the in-phase mode since anyinstability associated with this mode has the same phase for the cells on the ring. With this choice, (4.1.6) reduces to (cid:18) ω λ ( m − + M / (2 πν ) g λm ( m − g λm R λm + M m / (2 πν ) (cid:19) (cid:18) γ (cid:19) = σ πν (cid:18) γ (cid:19) , (4.1.12)where G λ ( m − eee = ω λ ( m − eee from (4.1.9).Upon eliminating σ from the × matrix problem (4.1.12) we obtain that γ ± are the roots of the quadratic equation γ + 12 πνg λm (cid:104) ( M − M m ) + 2 πν ( ω λ ( m − − R λm ) (cid:105) γ − ( m −
1) = 0 , (4.1.13)given by γ ± = − β λ ± (cid:113) β λ + 4( m − , where β λ ≡ πνg λm (cid:104) ( M − M m ) + 2 πν ( ω λ ( m − − R λm ) (cid:105) . (4.1.14)Since γ + γ − = − ( m − > , but with γ ± possibly complex-valued, we confirm that the two possible in-phase modes ccc ± = (1 , . . . , , γ ± ) T are orthogonal. The two eigenvalues σ = σ ± ( λ ) , given by σ ± ≡ πν (cid:0) ω λ ( m − + γ ± g λm (cid:1) + M , can bewritten as σ ± ( λ ) = ( h + h )2 ± (cid:113) ( h − h ) + 16 π ν ( m − g λm , where h ≡ πν ω λ ( m − + M , h ≡ πν R λm + M m . (4.1.15)13ere M = M ( λ ) and M m = M m ( λ ) are given in (4.1.6b) and ω λ ( m − is defined by G λ ( m − eee = ω λ ( m − eee from (4.1.9).Alternatively, rather than solving (4.1.12) for σ ± , and then setting σ ± ( λ ) = 0 by using a root-finder for λ , we can moredirectly conclude that M ( λ ) ccc = for ccc = ( eee, γ ) T if and only if the determinant of the × matrix in (4.1.12) vanishes. Inthis way, a discrete eigenvalue λ of the GCEP (2.2.13) for the in-phase mode, which satisfies det ( M ( λ )) = 0 , is a root of H c ( λ ) = 0 defined by H c ( λ ) = (cid:18) ω λ ( m − + M πν (cid:19) (cid:18) R λm + M m πν (cid:19) − ( m − g λm . (4.1.16)For any root of (4.1.16), the corresponding eigenvector ccc of M ( λ ) is ccc = (1 , . . . , , γ ) T , γ = − g λm (cid:18) ω λ ( m − + M πν (cid:19) . (4.1.17)It is readily verified using (4.1.15) that if λ (cid:63) satisfies either σ + ( λ (cid:63) ) = 0 or σ − ( λ (cid:63) ) = 0 , then we must have H c ( λ (cid:63) ) = 0 .In contrast, if λ (cid:63) satisfies H c ( λ (cid:63) ) = 0 , then we can only conclude that either σ + ( λ (cid:63) ) = 0 or σ − ( λ (cid:63) ) = 0 . Therefore, inimplementing a root-finding strategy based on the single scalar equation (4.1.16) instead of the two scalar equations (4.1.15)care must be taken to identify all possible roots of H c ( λ ) = 0 for the same parameter set.We summarize our result for eigenvalues λ of the GCEP (2.2.13) for a ring and center cell pattern as follows: Proposition 3.
Consider a ring and center hole pattern of m ≥ cells in the unit disk with cell centers at (4.1.1). The set Λ( M ) as obtained from the GCEP (2.2.13), and which approximates as ε → all the discrete eigenvalues of the linearizationof the PDE-ODE system (1.0.2) around the steady-state solution, is Λ( M ) ≡ (cid:110) λ (cid:12)(cid:12) m − (cid:91) j =1 { σ j ( λ ) = 0 } , (cid:91) { σ ± ( λ ) = 0 } (cid:111) . (4.1.18) Here σ j ( λ ) , for j = 1 , . . . , m − , for the anti-phase modes are defined in (4.1.11), while σ ± ( λ ) for the in-phase modes aredefined in (4.1.15). As shown in Remark 1 of Appendix A, due to mode degeneracy of G λ ( m − , there are ( m − / distinctanti-phase modes if m is odd and ( m − / distinct anti-phase modes if m is even. For the unit disk, where an infinite series representation of the solution to (2.2.8) is available, explicit formulae for theeigenvalues ω λj of G λ ( m − , as needed in (4.1.11) and (4.1.15), are given in Appendix A. The result in Remark 1 of AppendixA regarding mode degeneracy results from the fact that G λ ( m − is both symmetric and cyclic. In Appendix A we also showhow to readily calculate the quantities in (4.1.2) and (4.1.3), which are needed in (4.1.4) for determining the steady-state.As indicated by Proposition 1, stability boundaries in the τ versus D parameter space for the steady-state under Sel’kovkinetics are determined by HB boundaries where λ = iλ I ∈ Λ( M ) . To determine Hopf bifurcation boundaries for theanti-phase modes we set Re ( σ j ( iλ I )) = 0 and Im ( σ j ( iλ I )) = 0 for j = 1 , . . . , m − in (4.1.11). Upon separating real andimaginary parts in (4.1.7a) we obtain the following nonlinear algebraic system for each j = 1 , . . . , m − : πν Re ( ω λj ) + (cid:18) ν Dd (cid:19) + 2 πνDτ d d Re ( K c ( iλ I )) = 0 , Im ( ω λj ) + Dτ d d Im ( K c ( iλ I )) = 0 , (4.1.19a)whereIm ( K c ( iλ I )) = λ I ( det ( J ) − λ I ) + λ I det ( J ) tr ( J )( det ( J ) − λ I ) + ( λ I tr ( J )) , Re ( K c ( iλ I )) = det ( J )( det ( J ) − λ I ) − λ I tr ( J )( det ( J ) − λ I ) + ( λ I tr ( J )) . (4.1.19b)Here J is the Jacobian of the Sel’kov kinetics for the ring cells with determinant and trace given in (4.1.7b). Similarly,the Hopf bifurcation boundaries for the in-phase modes are obtained by setting σ ± ( iλ I ) = 0 in (4.1.15), which yields thenonlinear algebraic system Re ( σ ± ( iλ I )) = 0 , Im ( σ ± ( iλ I )) = 0 , (4.1.20)or equivalently H c ( iλ I ) = 0 from (4.1.16). We now apply the theory developed in §4.1 to a population of m = 3 cells, where two of the cells are equally spaced ona concentric ring of radius r within the unit disk, with the remaining one centered at the origin (see Fig. 3). For thisconfiguration, the eigenvalues of the × GCEP matrix M ( λ ) are given in (4.1.11) and (4.1.15) for a single anti-phasemode and the two in-phase modes, respectively. To compute the HB boundaries in the τ versus D parameter plane for thesemodes, we solve (4.1.19) and (4.1.20) numerically by implementing the psuedo-arclength continuation algorithm TEST_CON(cf. [1]) with respect to D , while using Newton’s method to compute τ and λ I at each point on the solution path. Such acontinuation scheme in D is needed owing to the possibility of fold points along the HB boundary.14o determine regions of instability in open sets of the τ versus D parameter plane we use (4.1.18) of Proposition 3, togetherwith the winding number criterion of complex analysis, to identify the number N of eigenvalues λ ∈ Λ( M ) with Re ( λ ) > .To do so, we first define F ( λ ) ≡ det( M ( λ )) , where M ( λ ) is the GCEP matrix in (4.1.6). Provided that there are no zeroesor poles on the imaginary axis, N is the number of zeroes of F ( λ ) = 0 in Re ( λ ) > , which from the argument principle is N = 12 π (cid:2) arg F ( λ ) (cid:3) Γ + P . (4.2.1)Here P is the number of poles of F ( λ ) in Re ( λ ) > , while (cid:2) arg F ( λ ) (cid:3) Γ denotes the change in the argument of F ( λ ) overthe closed, counter-clockwise oriented contour Γ . This contour Γ is the limit as R → ∞ of the union of the imaginaryaxis Γ I = iλ I , for | λ I | ≤ R , and the semi-circle Γ R , defined by | λ | = R with | arg ( λ ) | ≤ π/ . To count the number ofpoles of F ( λ ) in Re ( λ ) > , we must examine the analyticity properties of the GCEP matrix M ( λ ) in (4.1.6). Since theentries of the Green’s matrix G λ are analytic in Re ( λ ) > , any singularity of F ( λ ) must arise from the diagonal matrix K ( λ ) ≡ diag ( K c , K c , K m ) of (4.1.6b), which is given explicitly in (4.1.7) in terms of the Jacobians J and J m of the Sel’kovkinetics for the identical ring cells and the center cell, respectively. Since det ( J ) > , (4.1.7a) yields that K c has a complexconjugate pair of poles in Re ( λ ) > only if tr ( J ) > . Since this term involves two rows of M ( λ ) , P must be incrementedby four whenever tr ( J ) > . Similarly, since det ( J m ) > , P is increased by two when tr ( J m ) > for the center cell.With P determined in this way, we numerically compute the number of unstable eigenvalues N of the linearization of thesteady-state by evaluating (4.2.1) for each point ( D, τ ) in the τ versus D plane. For each such point, we numerically constructthe closed contour Γ for some value R (chosen so that Γ encloses all the poles of F ) in the complex λ -plane. As the closedcurve Γ is traversed in a counter-clockwise direction, the closed image curve F ( λ ) = F R + i F I is evaluated numerically inthe complex F -plane. The winding number, denoting the number of times F encloses/winds around the origin, is computednumerically from the algorithm of [2], and this is used to calculate (cid:2) arg F ( λ ) (cid:3) Γ . If the orientation of F ( λ ) = F R + i F I aroundthe origin is in the counter-clockwise direction, then (cid:2) arg F ( λ ) (cid:3) Γ is positive; otherwise, it is negative. In our computations,we chose R = 1 . (since any pole of F is close to the imaginary axis of the λ -plane), and we discretized the closed contour Γ = Γ I ∪ Γ R into subintervals, with Γ I having 800 subintervals while the semi-circle Γ R had subintervals. In the functionevaluation, the identity F ( λ ) = F ( λ ) was used to halve the computational effort. D (a) HB boundaries: Identical cells D (b) Number of unstable eigenvalues of the GCEPFigure 4: Left panel: HB boundaries in the τ versus D plane for a ring and center hole pattern of m = 3 identical cells with ring radius r = 0 . , parameters as in (4.3.1) and permeabilities d = 0 . and d = 0 . . The dashed and heavy solid curves are for the in-phasemodes computed from (4.1.20) with (+) and ( − ) , respectively. The thin solid curve is for the anti-phase mode computed from (4.1.19).Each mode is linearly unstable within its respective lobe. Linearly stable steady-state solutions exist outside the union of the lobes. FullPDE simulations of (1.0.2) are shown in Figs. 6 and 7 at the red and blue dots, respectively. Right panel: The regions of instabilitycomputed from the winding number (4.2.1). Blue region: in-phase ‘+’ mode is unstable with 2 roots of F ( λ ) = det( M ( λ )) = 0 inRe ( λ ) > . Green region: anti-phase mode is also unstable, yielding 4 roots. Magenta region: all modes are unstable, and there are 6roots. Plot on right is a zoom of the one on the left. The HB boundaries in the left panel are superimposed in this figure. In our results below, except when otherwise stated, the Sel’kov parameters α , µ and ζ , and permeabilities d and d for theidentical cells on the ring, and the common cell radius ε are α = 0 . , µ = 2 , ζ = 0 . , d = 0 . , d = 0 . , ε = 0 . . (4.3.1)From Fig. 2, we conclude that each cell, when isolated, has no intracellular oscillations. The permeabilities d and d forthe center cell will be stated in the figure captions below.Fig. 4a shows the computed HB boundaries in the τ versus D plane for a ring of radius r = 0 . when the cells areall identical. We observe that one of the in-phase lobes is open/unbounded, which predicts the existence of intracellular15scillations even for large D . In Fig. 4b, we show the corresponding regions of instability in the τ versus D parameter plane.In generating this figure, we pixelated the τ versus D plane with the uniform spacing ∆ τ = ∆ D = 0 . , and at each discretepoint ( D, τ ) used our winding number algorithm to count the number of roots N of det M ( λ ) = 0 in Re ( λ ) > . In Fig. 4b,each point in the blue-shaded region has two unstable eigenvalues for the GCEP, and they correspond to the in-phase ‘+’mode. The green-shaded region contains four unstable eigenvalues, two of which are for the anti-phase mode while the othertwo for the in-phase ‘+’ mode. Finally, in the magenta-shaded region there are six unstable eigenvalues of the linearization,with two such eigenvalues associated with each of the three possible modes of instability (in-phase ± and anti-phase). TheHB boundaries in the left panel of Fig. 4 are superimposed on these instability regions.The F ( λ ) = F R + i F I curve in the complex F -plane is shown in Fig. 5 for a specific point in each of the three instabilityregions in Fig. 4b. These plots show how F ( λ ) winds around the origin ( F R , F I ) = (0 , (shown with a green dot) as λ traverses Γ in the counterclockwise direction. For the point ( D, τ ) = (0 . , . in the magenta-shaded region, we observefrom the left panel of Fig. 5 that (cid:2) arg F ( λ ) (cid:3) Γ = 0 . At this point, F ( λ ) has two poles, one of order four and the other of order2, so that P = 6 . As such, (4.2.1) yields that there are 6 roots (counting multiplicity) to F ( λ ) = 0 in Re ( λ ) > . At thepoint ( D, τ ) = (0 . , . in the anti-phase mode instability region (green-shaded region in Fig. 4b), F ( λ ) winds round theorigin twice in the clockwise direction as shown in the middle panel of Fig. 5, which yields (cid:2) arg F ( λ ) (cid:3) Γ = − π . Since P = 6 at this point, (4.2.1) yields that F ( λ ) = 0 has four roots (counting multiplicity) in Re ( λ ) > . In the right panel of Fig. 5,we present a similar result for the point ( D, τ ) = (0 . , . in the blue-shaded region in Fig. 4b. At this point, we calculate (cid:2) arg F ( λ ) (cid:3) Γ = − π and that F has a pole of order four and a pole of order two in Re ( λ ) > . As such, (4.2.1) yields that F ( λ ) = 0 has two roots in Re ( λ ) > . -0.2 0 0.2 0.4 0.6-0.4-0.3-0.2-0.100.10.20.30.4 -0.01 0 0.01-505 10 -3 -1 -0.5 0 0.5 1 1.5-1.5-1-0.500.511.5 -8 -6 -4 -2 010 -3 -0.0100.01 -4 -2 0 2 4-5-2.502.55 Figure 5: The closed curve F ( λ ) = det( M ( λ )) plotted for a specific points in each of the three regions of instability in Fig. 4b. Theclosed contour Γ = Γ I ∪ Γ R in the λ -plane was constructed with R = 1 . , and discretized with subintervals on Γ I and subintervalson Γ R . The green dot represents the origin ( F R , F I ) = (0 , , the black dot indicates the starting point of the curve, and the blue arrowsshow the direction of the curve. The inserts detail the behavior near the origin. Left panel: for ( D, τ ) = (0 . , . in the magenta-shaded region of Fig. 4b, we have (cid:2) arg F ( λ ) (cid:3) Γ = 0 . Middle panel: for ( D, τ ) = (0 . , . in the green-shaded region in Fig. 4b, wehave (cid:2) arg F ( λ ) (cid:3) Γ = − π . Right panel: for ( D, τ ) = (0 . , . in the blue-shaded region of Fig. 4b, we have (cid:2) arg F ( λ ) (cid:3) Γ = − π . The real and imaginary parts of the normalized eigenvector ccc of the GCEP matrix M in (4.1.6) is given in Table 1for selected points on the HB boundaries in Fig. 4a. Recall from (2.2.17) that the magnitude of the components of theeigenvector ccc measure the diffusive flux at the boundary of each cell, while ˜ ccc ≡ K ccc predicts the relative amplitude and phaseshifts of the intracellular oscillations within the cells at the Hopf bifurcation point. For our ring and center-cell pattern K = diag ( K c , K c , K m ) , where K c and K m are given in (4.1.7a).mode ( D, τ ) j (cid:16) Re ( c j ) , Im ( c j ) (cid:17) θ j ( rad ) (cid:16) Re (˜ c j ) , Im (˜ c j ) (cid:17) (0 . ,
0) 0 ( − . , . In-phase ( + ) (1 . , . (0 . ,
0) 0 ( − . , . (dashed curve) 3 (0 . , . . − . , . ( − . , − . .
15 (0 . , − . In-phase ( − ) (0 . , . ( − . , − . .
15 (0 . , − . (heavy solid) 3 (0 . ,
0) 0 ( − . , . (0 . ,
0) 0 ( − . , . Anti-phase (0 . , . ( − . , π (0 . , − . (thin solid) 3 (0 ,
0) 0 (0 , Table 1: Real and imaginary parts of the eigenvector ccc of the GCEP matrix M ( λ ) in (4.1.6), together with ˜ ccc ≡ K ccc , as computed fora few points on the HB boundaries in Fig. 4a for three identical cells. The second to last column shows the phase shifts measured interms of the angle each component of the vector ccc makes with the positive real axis in anticlockwise direction. In Fig. 6 we show full numerical simulations of the coupled PDE-ODE model (1.0.2) obtained using the commercial PDE16oftware package FlexPDE [14] for τ = 0 . and D = 1 , which corresponds to the red dot in the phase diagram of Fig. 4a. Weobserve from the results in this figure that the intracellular dynamics of the cells are synchronized with a very slight phaseshift, which agrees with the prediction by the eigenvector K ccc in the first three rows of Table 1 from the linearized theory.Although the cells have identical parameters, the center and ring cells have slightly different dynamics owing to the fact thatthe full Green’s matrix is not cyclic for a ring and center cell pattern. Our numerical computations of det M ( λ ) = 0 for ( D, τ ) = (1 , . using the GCEP matrix in (4.1.6) yields that Re ( λ ) ≈ . , Im ( λ ) ≈ . , Re ( ccc ) ≈ (0 . , . , . and Im ( ccc ) ≈ (0 , , . . Observe that the eigenvector is rather close to that on the nearby point on the HB boundary, asgiven in the first three rows of Table 1. The prediction from linearized theory is that the period of oscillations is approximately π/ Im ( λ ) ≈ . , which is rather close to the period observed in the full PDE simulations of Fig. 6. A similar full numericalresult is presented in Fig. 7 for τ = 1 . and D = 1 , corresponding to the blue dot in Fig. 4a. At this pair ( D, τ ) , our phasediagram predicts no intracellular oscillations. This is confirmed from the full numerical results shown in Fig. 7.
470 480 490 50011.21.41.61.82
Ring cellsCenter cell
470 480 490 5000.40.50.60.70.8
Ring cellsCenter cell
Figure 6: Full PDE simulations of (1.0.2), computed with FlexPDE [14], for τ = 0 . and D = 1 for three identical cells correspondingto the red dot in Fig. 4a. Left panel: surface plot at time t = 400 . Middle panel: intracellular species u versus t . Right panel:intracellular species u versus t . The blue and red curve is for the cells on the ring and the center cell, respectively. Ring cellsCenter cell
Ring cellsCenter cell
Figure 7: Same caption as in Fig. 6 except that now τ = 1 . and D = 1 , corresponding to the blue dot in Fig. 4a. There are nointracellular oscillations and the steady-state is linearly stable. In the left panel of Fig. 8 we show the computed HB boundaries in the τ versus D plane for the same parameters as inFig. 4a, except that the influx rate into the center cell is reduced to d = 0 . (keeping d = 0 . ). With this lower influxrate, Fig. 8 shows that the in-phase lobes are now closed, so that there are no longer intracellular oscillations when D is large.Qualitatively, when D is large, the bulk chemical diffuses quickly in the entire disk and there is insufficient feedback of itinto the center cell, owing to the smaller value of d , to sustain intracellular oscillations. The winding number algorithm of(4.2.1) can be used, with the same result as in Fig. 4, to determine the number of unstable eigenvalues of the GCEP withinthe lobes (not shown). In the right panel of Fig. 8 we give the real and imaginary parts of the normalized eigenvector ccc ofthe GCEP matrix in (4.1.6) at a few selected points on the HB boundaries shown in the left panel of Fig. 8.In Fig. 9, we show full FlexPDE [14] numerical simulations of (1.0.2) for τ = 0 . and D = 0 . , which corresponds to the reddot in the left panel of Fig. 8. Our numerical computations of det M ( λ ) = 0 for ( D, τ ) = (0 . , . using the GCEP matrix in(4.1.6) yields that Re ( λ ) ≈ . , Im ( λ ) ≈ . , Re ( K ccc ) ≈ ( − . , − . , − . and Im ( K ccc ) ≈ (0 . , . , . .As such, our linearized theory predicts that the center cell will have larger amplitude oscillations near onset than the ringcells, and there will be a ≈ ◦ phase shift between the oscillations. Our FlexPDE numerical results in the middle andright panels of Fig. 9 show that the prediction of the linearized theory does extend to the fully nonlinear regime in that thedefective center cell has larger amplitude oscillations than do the identical ring cells, and there is a slight phase shift in theoscillations. From the surface plot shown in the left panel of Fig. 9, we observe that the coupling between the cells mediatedby the bulk medium is rather weak. Moreover, the rather large concentration of the bulk chemical close to the center cell, withflux measured by the modulus of the third component of Re ( ccc ) ≈ (0 . , . , . and Im ( ccc ) ≈ (0 . , . , . , is dueto its smaller rate d = 0 . of influx into the center cell than for the identical ring cells centered at ( ± . , . Paradoxically,17 .2 0.4 0.6 0.8 1 D mode ( D, τ ) j Re( c j ) Im( c j ) θ j ( rad ) .
345 0 .
171 0 . In-phase ( + ) (0 . , . .
345 0 .
171 0 . (dashed curve) 3 .
839 0 0 − .
446 0 .
256 2 . In-phase ( − ) (0 . , . − .
446 0 .
256 2 . (heavy solid) 3 .
686 0 0 .
707 0 0
Anti-phase (0 . , . − .
707 0 π (thin solid) 3 Figure 8: Left panel: same caption as in Fig. 4a except that now the center cell is a defector with permeabilities d = 0 . and d = 0 . ,corresponding to a reduced influx into the center cell. Full PDE simulations are shown in Figs. 9 and 10 at the red and blue dots,respectively. Right panel: same caption as in Table 1 except that now d = 0 . and d = 0 . . The real and imaginary parts of theeigenvector ccc of the GCEP matrix M ( λ ) in (4.1.6) are computed for a few points on the HB boundaries shown in the left panel.
450 460 470 480 490 50011.522.5
Ring cellsCenter cell
450 460 470 480 490 5000.40.60.81
Ring cellsCenter cell
Figure 9: Full PDE simulations of (1.0.2), computed with FlexPDE [14], for τ = 0 . and D = 0 . , corresponding to the red dot in theleft panel of Fig. 8. The center cell is a defector with permeabilites d = 0 . and d = 0 . . Left panel: surface plot at time t = 400 .Middle panel: intracellular species u versus t . Right panel: intracellular species u versus t .
200 300 400 500 600 7000123
Ring cells
200 300 400 500 600 7000123
Center cell
Figure 10: Same caption as in Fig. 9, except that the FlexPDE [14] simulation of (1.0.2) is done at τ = 0 . and D = 0 . , correspondingto the blue dot in the left panel of Fig. 8. At this point, both in-phase modes and the anti-phase mode are unstable. The beating-behaviorobserved arises from the fact that these three modes all have comparable frequencies. D mode ( D, τ ) j Re( c j ) Im( c j ) θ j ( rad ) .
707 0 0
In-phase ( + ) (0 . , . .
707 0 0 (dashed curve) 3 . − . . − . . . In-phase ( − ) (1 . , . − . . . (heavy solid) 3 .
986 0 0 .
707 0 0
Anti-phase (0 . , . − .
707 0 π (thin solid) 3 Figure 11: Same caption as in Fig. 8 except that now the center cell is a defector with permeabilities d = 0 . and d = 0 . ,corresponding to a reduced influx and a larger efflux out of the center cell. Full PDE simulations are shown in Fig. 12 at the blue dotin the left panel. however, this large buildup of the bulk signal near the center cell counteracts the relatively smaller rate of influx into thecenter cell, and has the effect of triggering a larger amplitude oscillation in the center cell than for the ring cells.In Fig. 10, we show full FlexPDE [14] numerical simulations of (1.0.2) for τ = 0 . and D = 0 . , corresponding to the bluedot in the left panel of Fig. 8. As seen from Fig. 8, this point is located within the region of instability that is common to allthree modes of instability. By solving det M ( λ ) = 0 for ( D, τ ) = (0 . , . numerically, the eigenvalues λ and eigenvectors ccc for the three modes are:in-phase (+): λ ≈ . . i , Re ( c ) ≈ (0 . , . , . , Im ( c ) ≈ (0 , , − . , in-phase (-): λ ≈ . . i , Re ( c ) ≈ (0 . , . , − . , Im ( c ) ≈ (0 , , − . , anti-phase: λ ≈ . . i , Re ( c ) ≈ (0 . , − . , , Im ( c ) = (0 , , . (4.3.2)We observe that this linearized theory predicts three distinct unstable oscillatory modes with roughly similar frequenciesIm ( λ ) and growth rates Re ( λ ) . The beating-type intracellular oscillations observed in Fig. 10 is likely related to the well-known linear phenomenon of superimposing two or more single-mode oscillations with comparable frequencies.
450 460 470 480 490 50011.52
Ring cellsCenter cell
450 460 470 480 490 5000.30.40.50.60.70.8
Ring cellsCentre cell
Figure 12: Full PDE simulations of (1.0.2), computed with FlexPDE [14], for τ = 3 and D = 0 . , corresponding to the blue dot in theleft panel of Fig. 11. The center cell is a defector with permeabilites d = 0 . and d = 0 . . Left panel: surface plot at time t = 400 .Observe the buildup of the bulk chemical near the center cell in comparison to the two ring cells. Middle panel: intracellular species u versus t . Right panel: intracellular species u versus t . As predicted by the second row in the table in the right panel of Fig. 11, weconfirm that the center cell has larger amplitude oscillations than do the ring cells. In the left panel of Fig. 11 we plot the HB boundaries for the case where the center cell has permeabilities d = 0 . and d = 0 . . This corresponds to a larger secretion or efflux rate for the center cell, while the feedback it receives from the ringcells is reduced. Unlike the HB boundaries shown in Figs. 4 and 8 where the instability regions for the modes are nested withineach other, we observe from the insert in the left panel of Fig. 11 that only the in-phase (+) mode and the anti-phase modeoverlap. Moreover, since the other in-phase ( − ) lobe is unbounded in D , intracellular oscillations always occur in some rangeof τ as D increases. In the right panel of Fig. 11, we give the real and imaginary parts of the normalized eigenvector ccc of theGCEP matrix (4.1.6) at a few points on the HB boundaries. From the second row of this table, the linearized theory suggeststhat the amplitude of intracellular oscillations associated with the dominant in-phase ( − ) instability lobe will be much largerin the center cell than in the identical ring cells at the HB point, and that there will be a significant phase shift in theoscillations between the center cell and the ring cells. More precisely, at the point ( D, τ ) = (0 . , interior to the instabilitylobe, we solve det M ( λ ) = 0 numerically to obtain that Re ( λ ) ≈ . , Im ( λ ) ≈ . , Re ( ccc ) ≈ ( − . , − . , . ,19m ( ccc ) ≈ (0 . , . , , Re ( K ccc ) ≈ ( − . , − . , − . , and Im ( K ccc ) ≈ (0 . , . , . . From K ccc and Im ( λ ) thisindicates that the center cell will have much larger oscillations near onset than the ring cells, with a period of oscillations of ≈ and with a phase shift of ≈ ◦ between the ring and center cell oscillations. From the FlexPDE simulations of (1.0.2)when τ = 3 and D = 0 . shown in Fig. 12, corresponding to the blue dot in the left panel of Fig. 11, we observe that thesepredictions of the linearized theory are roughly satisfied. Moreover, since the efflux rate out of the center cell is larger, whilethe influx rate is smaller, the rather small bulk diffusivity D = 0 . should qualitatively lead to a buildup of the bulk chemicalnear the center cell at certain times in the oscillation. This feature is observed in the surface plot in the left panel of Fig. 12,and provides a clear example of the diffusion-sensing behavior, as regulated by the permeability parameters.Next, we examine whether the ODE system (3.0.16) of Proposition 2, derived under the assumption of large bulk diffusivity D = O ( ν − ) (cid:29) , can still be used to reliably approximate the intracellular dynamics observed in the full FlexPDEsimulation results of (1.0.2), performed for O (1) values of D , in Figs. 6, 7, 9 and 12. In Fig. 13 we show that the numericalresults computed from the ODE system (3.0.16) compare surprisingly well with the full PDE simulations with respect to theamplitude and period of intracellular oscillations (first, third and fourth rows of Fig. 13) and the prediction of a linearly stablesteady-state (second row of Fig. 13). In using the ODE system (3.0.16) we calculated D as D = Dν , where ν = − / log ε with ε = 0 . . We remark that the beating-type oscillations observed in Fig. 10 for the very small value D = 0 . arenot captured by the ODE system (3.0.16). Moreover, we emphasize that the simpler ODE system corresponding to thewell-mixed limit D → ∞ as given in (3.0.19), and which was used in [20] and [26] for studying quorum-sensing behavior,does not reliably approximate the intracellular oscillations for the values of bulk diffusivity given in Figs. 6, 7, 9, 9, and 12. Next, we study how the HB boundaries in the ( D, τ ) plane are altered by varying a Sel’kov kinetic parameter of the centercell. In Fig. 14 we plot the HB boundaries for a ring and center cell pattern of m = 3 cells, where the identical cells onthe ring have parameters as in (4.3.1), while the Sel’kov parameter α for the defective center cell is either α = 0 . (redcurves), α = 0 . (blue curves), or α = 0 . (same as in the left panel of Fig. 4). The HB boundaries in the right panelof Fig. 14 show a zoom for smaller values of D than the figure in the left panel. In these figures, the dashed and the heavysolid curves are for the in-phase modes computed from (4.1.20) with (+) and ( − ) , respectively, while the thin solid curve isfor the anti-phase mode computed from (4.1.19). We observe that the HB boundary for the anti-phase mode is independentof the Sel’kov kinetic parameter α for the center cell. This follows from the facts that the steady-state solution to (4.1.4)for the ring cells and its Jacobian J in (4.1.7b) do not depend on α . As a result, the anti-phase HB boundary, computedfrom (4.1.19), is independent of α .From the left panel of Fig. 14 we observe that as α increases the parameter regions in the ( D, τ ) plane where intracellularoscillations occur decreases. In particular, the in-phase (+) instability lobe that was open for α = 0 . becomes closedwhen α = 0 . , thereby precluding the possibility of intracellular oscillations when D is sufficiently large. To interpretthis result, we observe from the middle panel of Fig. 2 that, as α increases, the center cell becomes less activated, withthe parameters drifting further from the HB boundary for the uncoupled cell. As a result, it becomes more difficult totrigger in-phase intracellular oscillations for a group of coupled cells as α increases. Overall, Fig. 14 does show that rathersmall increases or decreases in a parameter value of the nonlinear Sel’kov kinetics for one specific cell can either extinguishor trigger intracellular oscillations for an entire group of cells. The corresponding eigenvectors for selected points on thebifurcation diagrams in Fig. 14 are shown in Table 2 for α = 0 . and α = 0 . . α = 0 . α = 0 . mode j ( D, τ ) ( Re c j , Im c j ) θ j ( rad ) ( D, τ ) ( Re c j , Im c j ) θ j ( rad ) (0 . , − . .
14 (0 . , + ) 2 (0 . , . . , − . .
14 (0 . , . . , (0 . ,
0) 0 (0 . , − . . − (0 . , . .
23 ( − . , . . In-phase ( − ) 2 (0 . , . − (0 . , . .
23 (0 . , . − . , . . (dashed curve) 3 (0 . ,
0) 0 (0 . ,
0) 0 (0 . ,
0) 0 (0 . , (0 . , . − (0 . , π (0 . , . − ( . , π (thin solid) 3 (0 ,
0) 0 (0 ,
0) 0
Table 2: Real and imaginary parts of the eigenvector ccc of the GCEP matrix M ( λ ) in (4.1.6), computed for a few points on the HBboundaries in Fig. 14. The cells on the ring are identical and located at ( ± . , with parameters as given in (4.3.1). The center cellhas the same parameters, with the exception that it has a Sel’kov kinetic parameter α different than those on the ring. Middle threecolumns: α = 0 . . Last three columns: α = 0 . .
70 480 490 50011.21.41.61.82
470 480 490 5000.40.50.60.70.8
450 460 470 480 490 50011.522.5
450 460 470 480 490 5000.40.60.81
450 460 470 480 490 50011.52
450 460 470 480 490 5000.30.40.50.60.70.8
Figure 13: Numerical results for intracellular dynamics versus time computed from the ODE system (3.0.16) of Proposition 2 cor-responding to the PDE simulations shown in Fig. 6 (first row), Fig. 7, (second row), Fig. 9 (third row) and Fig. 12 (fourth row),respectively. In each case, although D is not large, results from the ODE system (3.0.16) are seen to compare surprisingly well with thefull FlexPDE simulations of the PDE-ODE system (1.0.2). First row: ( D, τ ) = (1 . , . (red dot in Fig. 4a). Compare with Fig. 6.Second row: ( D, τ ) = (1 . , . (blue dot in Fig. 4a). Compare with Fig. 7. Third row: ( D, τ ) = (0 . , . (red dot in the left panel ofFig. 8). Compare with Fig. 9. Fourth row: ( D, τ ) = (0 . , (blue dot in the left panel of Fig. 11). Compare with Fig. 12. Althoughthere is a phase shift between the ODE and full PDE results for intracellular oscillations due to different initial conditions used, theODE system (3.0.16) captures well the amplitude and period of intracellular oscillations observed in the full PDE simulations. For a ring and center cell pattern, with cells centered at ( ± r , and (0 , , we study how the HB boundaries in the ( D, τ ) plane depend on the ring radius r . In the first row of Fig. 15 we show the HB boundaries for the in-phase and anti-phase modes when r = 0 . , r = 0 . (same as Fig. 4a), and r = 0 . for the case where the cells are all identical withpermeabilities d = 0 . and d = 0 . . Similar plots for the three values of r are shown in the second row of Fig. 15 for the21
10 20 30 40 D D Figure 14: HB boundaries in the τ versus D plane for a ring and center cell pattern of m = 3 cells, where the Sel’kov kinetic parameter α for the center cell is varied. The ring cells are centered at ( ± , and, with the exception of the Sel’kov kinetic parameter for thecenter cell, all three cells have parameters as in (4.3.1). Left panel: α = 0 . (red curves), α = 0 . (black curves), α = 0 . (bluecurves). Right panel: a zoomed-in version of the left panel for smaller values of D . In both panels, the dashed and heavy solid curvesare for the in-phase modes computed from (4.1.20), while the thin solid curve is for the anti-phase mode computed from (4.1.19). Eachmode is linearly unstable in their respective lobes, while linearly stable steady-state solutions exists outside the union of the lobes. Theanti-phase HB boundaries are independent of α , and so are plotted on top of each other. case where the center cell is now defective with d = 0 . and d = 0 . (Fig. 15e for r = 0 . is the same as in Fig. 8).When the cells are all identical, we observe from the top row in Fig. 15 that the regions of instability for the anti-phasemode and the in-phase ( − ) mode shrink noticeably as r decreases. In contrast, the instability region for the in-phase (+) mode is relatively insensitive to changes in r . To qualitatively explain this observation, we examine the eigenvector of theGCEP matrix given in Table 1 at points on the HB boundaries for the three modes. From this table, for the in-phase ( − ) mode only the two ring cells oscillate in phase while the center cell has larger amplitude oscillations that are roughly ◦ out of phase. For the anti-phase ( − ) mode, the center cell is quiescent while the ring cells oscillate ◦ out of phase. Finally,for the dominant in-phase (+) mode, the three identical cells are synchronized with very similar amplitudes and phases.Intuitively, as the ring radius r decreases the cells become more clustered and, as a result, intracellular dynamics can stillbe synchronized even for smaller values of the bulk diffusivity D . As a result, when r is small, the anti-phase and in-phase ( − ) instability lobes, in which the center and ring cells are not synchronized, should exist only for very small values of D (see Fig. 15a). When D is very small, communication between cells that are close can still be rather weak. Notice that for r = 0 . , where the cells are farther apart, the anti-phase lobe in Fig. 15c exists for larger values of D than when r = 0 . or r = 0 . . For r = 0 . , the ring cells can oscillate, maintaining a quiescent center cell, provided that the bulk chemicalis not washed away, i.e. D is not too large, owing to the fact that each of the two ring cells are now relatively close to theirimages across the domain boundary due to the reflecting boundary condition imposed. However, when r is large, the ringcells are far from each other and so, within the anti-phase instability lobe, their oscillations are ◦ out of phase.The second row of Fig. 15 shows similar results for the case where the center cell is defective, with permeabilities d = 0 . and d = 0 . , corresponding to a reduced rate of influx into the center cell. Similar to the results for identical cellspresented in the first row of Fig. 15, as r decreases the region of instability for the in-phase ( − ) and the anti-phase modesshrink. Moreover, we observe that the maximum extent in D of the in-phase (+) lobes, in which all the cells are essentiallysynchronized in amplitude and phase, decreases as r increases. This is because when the cells are farther apart, a relativelysmaller value of the bulk diffusivity can lead to a washing out of the bulk signal, which thereby weakens the communicationbetween the three cells and precludes synchronization. As the number of cells increases it becomes increasingly more challenging numerically to determine the stability boundariesin parameter space from the root-finding condition det M ( λ ) = 0 given an arbitrary cell configuration xxx j with arbitrarypermeability parameters d j and d j for j = 1 , . . . , m . In this section we provide a simpler approach to determine thestability boundaries from the GCEP (2.2.13) for the large bulk diffusion parameter regime where D = O ( ν − ) . To isolateonly the effect of different cell configurations, such as shown in Fig. 16, as well as the effect of different cell permeabilities,in our analysis below we will assume that the Sel’kov kinetic parameters are all cell-independent, i.e. that α j = α , µ j = µ and ζ j = ζ for j = 1 , . . . , m . Our analysis is easily extended to remove this assumption.We first recall from (2.2.21b) of §2.2 that, for ε → , the steady-state solution with Sel’kov kinetics is determined in termsof the solution AAA = ( A , . . . , A m ) T to the linear system (cid:18) I + 2 πν G + ν D P + 2 πνDτ P (cid:19) A = − µνP eee ; P = diag (cid:16) d , . . . , d m (cid:17) P = diag (cid:16) d d , . . . , d m d m (cid:17) , (5.0.1a)22 D (a) r = 0 . D (b) r = 0 . D (c) r = 0 . D (d) r = 0 . D (e) r = 0 . D (f) r = 0 . Figure 15: HB boundaries in the τ versus D plane for a ring and center hole configuration of m = 3 cells, for three different ringradii r as indicated. The Sel’kov parameters and cell radius are as in (4.3.1). Top row: identical cells with permeabilities d = 0 . and d = 0 . . Bottom row: center cell is defective with d = 0 . and d = 0 . . The dashed and heavy solid curves are the HBboundaries for the in-phase (+) and ( − ) modes, respectively, as computed from (4.1.20). The thin solid curve is the HB boundary forthe anti-phase mode computed from (4.1.19). Each mode is unstable in their respective lobes, while linearly stable steady-state solutionsoccur outside the union of the lobes.(a) Randomly placed cells (b) Two clusters of cellsFigure 16: Schematic showing different cell configurations in the unit disk. The cells are represented by smaller disks (in cyan) and thediffusing bulk species corresponds to the red dots. Left panel: randomly placed cells. Right panel: two spatially clustered groups of cells. where eee = (1 , . . . , T and G is the Green’s interaction matrix defined in (2.1.7). In terms of the solution AAA to (5.0.1a), weobtain from (2.2.21a) that the steady-state for the intracellular species uuu e = ( u ej , u ej ) T in each cell is u ej = µ + 2 πDτ A j , and u ej = µα + ( u ej ) . (5.0.1b)Since this simple steady-state construction is accurate to all orders in ν for any D > , we will not seek an approximationto it valid for the regime D = O ( ν − ) . Therefore, in the GCEP matrix M ( λ ) given in (2.2.13a), the diagonal matrix K ( λ ) ,as defined in (2.2.12a), will be evaluated at the solution to (5.0.1). With Sel’kov kinetics, this yields K = K ( λ ) ≡ diag (cid:16) K , . . . , K m (cid:17) , where K j = λ + det ( J j ) λ − λ tr ( J j ) + det ( J j ) , (5.0.2a)where J j is the Jacobian matrix for the local Sel-kov kinetics at the j th cell, for whichdet ( J j ) = ζ (cid:0) α + ( u ej ) (cid:1) , tr ( J j ) = (cid:104) µu ej − (cid:0) α + ( u ej ) (cid:1) − ζ (cid:0) α + ( u ej ) (cid:1) (cid:105) α + ( u ej ) . (5.0.2b)23ur goal in this subsection of determining a more tractable formula for determining the roots of det M ( λ ) = 0 when m islarge and in the limit D = D /ν , where ν = − / log ε and D = O (1) , is based on approximating the eigenvalue-dependentGreen’s matrix G λ in (2.2.13a), while retaining the full K matrix, as defined in (5.0.2).To do so, we readily calculate from (2.2.8) that when D = D /ν (cid:29) we have (see [20]) G λ = D ν | Ω | (1 + τ λ ) + G + O ( ν ) . (5.0.3)Here G is the Neumann Green’s matrix of (3.0.17) given in terms of the Neumann Green’s function G of (3.0.14), which isknown analytically for the the unit disk in (3.0.21). By using (5.0.3) in (2.2.13a), we find that the GCEP matrix reduces to M ( λ ) = M ∞ ( λ ) + O ( ν ) , where M ∞ ≡ B + γE + 2 πν G . (5.0.4a)In (5.0.4a) the scalar γ , the rank-one matrix E and the diagonal matrix B are defined by γ = γ ( λ ) ≡ πmD (1 + τ λ ) | Ω | , E ≡ meee eee T , B = B ( λ ) ≡ I + D P + 2 πD τ P K ( λ ) . (5.0.4b)Neglecting terms of order O ( ν ) , in the limit D = D /ν (cid:29) we conclude that the GCEP (2.2.13) reduces to finding valuesof λ for which there is a nontrivial solution ccc to M ∞ ( λ ) ccc = 000 , which occurs iff det M ∞ ( λ ) = 0 . (5.0.4c)We remark that in our asymptotic reduction (5.0.4) of the GCEP (2.2.13), the cell configuration enters through the diagonalmatrix K , as defined in (5.0.2), and at order ν from the Neumann Green’s matrix G in (5.0.4a).To analyze (5.0.4) we first observe from (5.0.4b) and (5.0.2a) that B = diag ( b , . . . , b m ) , where b j = ( d j + D ) d j (cid:104) η j τ K j (cid:105) , η j ≡ πd j D d j + D , j = 1 , . . . , m . (5.0.5)The key advantage of our asymptotic reduction is that by using the matrix structure in (5.0.4) we can readily calculatea tractable analytical expression for det M ∞ ( λ ) , which can be used in implementing the root-finding condition in (5.0.4c)numerically. Our first result for (5.0.4) in this direction is as follows: Proposition 4.
Suppose that b j (cid:54) = 0 for j = 1 , . . . , m so that B is invertible. Then,det M ∞ ( λ ) = det ( B + γE + 2 πν G ) = (cid:32) m (cid:89) i =1 b i (cid:33) κ ε (1 + O ( ν )) , (5.0.6) where κ ε ≡ γmeee T vvv + 2 πν (cid:18) vvv T G vvv eee T vvv (cid:19) + O ( ν ) , with vvv T ≡ (1 /b , . . . , /b m ) . (5.0.7) Proof.
Since B is a diagonal matrix and invertible by assumption, we havedet ( B + γE + 2 πν G ) = det ( B ) det (cid:0) I + γ B − E + 2 πν B − G (cid:1) , (5.0.8)where det B = (cid:81) mi =1 b i . The second term in (5.0.8) is calculated by multiplying all the eigenvalues κ ε of (cid:0) I + γ B − E + 2 πν B − G (cid:1) vvv ε = κ ε vvv ε . (5.0.9)To determine O ( ν ) accurate expressions for these eigenvalues we expand an eigenpair κ ε and vvv ε as κ ε = κ + ν ˜ κ + · · · , vvv ε = vvv + ν ˜ vvv + · · · . (5.0.10)Upon substituting (5.0.10) into (5.0.9) we equate powers of ν to obtain the leading-order problem ( I + γ B − E ) vvv = κ vvv , (5.0.11)and the following problem at O ( ν ) : ( I + γ B − E − κ I ) ˜ vvv = ˜ κ vvv − π B − G vvv . (5.0.12)For the leading order problem (5.0.11), we first observe that E qqq j = 000 for j = 2 , . . . , m , where span { qqq , . . . , qqq m } is the m − dimensional subspace orthogonal to eee . We readily conclude that (5.0.11) has the eigenpairs κ j = 1 and vvv j = qqq j , for j = 2 , . . . , m . (5.0.13)24o determine the O ( ν ) correction to these eigenvalues, we substitute (5.0.13) into (5.0.12) to obtain γ B − E ˜ vvv = ˜ κqqq j − π B − G qqq j . (5.0.14)Since B qqq j , with B T = B , is a left null-vector of B − E , the solvability condition for (5.0.14) is qqq Tj B (cid:0) ˜ κqqq j − π B − G qqq j (cid:1) = 0 ,which determines ˜ κ for each eigenpair. In this way, we obtain that m − eigenvalues of (5.0.9) are κ jε = 1 + 2 πν qqq Tj G qqq j qqq Tj B qqq j + O ( ν ) , where qqq Tj eee = 0 , for j = 2 , . . . , m . (5.0.15)To find the remaining eigenpair of the leading order problem we write (5.0.11) as vvv (1 − κ ) = − γm (cid:0) eee T vvv (cid:1) (1 /b , . . . , /b m ) T .This has the (unnormalized) nontrivial solution vvv = (1 /b , . . . , /b m ) T , iff κ = 1 + γmeee T vvv = 1 + γm m (cid:88) i =1 b i . (5.0.16)To determine the O ( ν ) correction to this eigenvalue we observe, for this eigenpair, that eee is a left null-vector of the matrixin (5.0.12) since, by using (5.0.16) for κ , we calculate eee T (cid:0) I + γ B − E − κ I (cid:1) = eee T (cid:32) γm (cid:32) m (cid:88) i =1 b i (cid:33) − κ (cid:33) = 000 . (5.0.17)Therefore, upon left-multiplying (5.0.12) by eee T , the solvability condition for (5.0.12) is eee T (cid:0) ˜ κvvv − π B − G vvv (cid:1) = 0 , whichyields ˜ κ = 2 πvvv T G vvv / (cid:0) eee T vvv (cid:1) where we used eee T B − = vvv T . We conclude that a two-term expansion κ ε = κ + ˜ κν for thisremaining eigenvalue of (5.0.9) is as given in (5.0.7). Finally, by multiplying κ ε with the other eigenvalues given in (5.0.15)we obtain det (cid:0) I + γ B − E + 2 πν B − G (cid:1) = κ ε (1 + O ( ν )) . In view of (5.0.8) this completes the derivation of (5.0.6).The key assumption in Proposition 4 is that b j (cid:54) = 0 for any j = 1 , . . . , m . By using (5.0.5) and (5.0.2a) for b j and K j ,respectively, we conclude that b j = 0 if and only if λ is a root of the quadratic equation Q j ( λ ) = 0 , where Q j ( λ ) ≡ λ − (cid:16) tr ( J j ) − η j τ (cid:17) λ + det ( J j ) (cid:104) η j τ (cid:105) , where η j ≡ πd j D d j + D , (5.0.18)with det ( J j ) and tr ( J j ) as given in (5.0.2b). With this criterion and together with Proposition 4 we readily formulate asimple scalar root-finding problem to identify values of λ for which the reduced GCEP (5.0.4) has a nontrivial solution. Proposition 5.
Suppose that λ = λ (cid:63) is a root of Q s ( λ ) = 0 , where Q s ( λ ) ≡ γmeee T vvv + 2 πν (cid:18) vvv T G vvv eee T vvv (cid:19) with vvv T ≡ (1 /b , . . . , /b m ) . (5.0.19) Here γ = γ ( λ ) and b j = b j ( λ ) for j = 1 , . . . , m are given in (5.0.4b) and (5.0.5), respectively. Suppose that λ (cid:63) satisfies λ (cid:63) / ∈ m (cid:91) j =1 { λ j ± } , where Q j ( λ j ± ) = 0 , (5.0.20) and Q j ( λ ) is the quadratic defined in (5.0.18). Then, det M ∞ ( λ (cid:63) ) = 0 and the corresponding (unnormalized) nontrivialsolution ccc to the reduced GCEP (5.0.4c) is ccc = (cid:18) b ( λ (cid:63) ) , . . . , b m ( λ (cid:63) ) (cid:19) T + O ( ν ) (5.0.21)In this way, we can use Proposition 5 to determine Hopf bifurcation boundaries in the τ versus D parameter plane bysimply letting λ = iλ I , with λ I > , and settingRe [ Q s ( iλ I )] = 0 , and Im [ Q s ( iλ I )] = 0 , (5.0.22)while ensuring that the condition (5.0.20) holds with λ (cid:63) = iλ I . A sufficient condition for (5.0.20) to hold along solutions of(5.0.22) as parameters are varied is that tr ( J j ) (cid:54) = η j /τ for all j = 1 , . . . , m .A simple analytically tractable special case of Proposition 5 is when the permeabilities are all identical, i.e. d j = d c , d j = d c , and when the cell configuration { xxx , . . . , xxx m } is such that eee is an eigenvector of the reduced-wave Green’s matrix25 , and consequently the Neumann Green’s matrix G . Therefore, G eee = βeee for some eigenvalue β . A ring patterns of cellsconcentric within the unit disk with common cell permeabilities and Sel’kov parameters is an example of such a cell pattern.For this case, (5.0.1) admits the solution A = A c eee and so J j = J c for j = 1 , . . . , m . Since b j = b c for j = 1 , . . . , m , we cantake vvv = b − c (1 , . . . , T and readily obtain that the root-finding condition Q s ( λ ) = 0 in (5.0.19) reduces to γb c + 2 πνmb c eee T G eee = 0 , where G eee = βeee , (5.0.23)while from (5.0.18) we have Q j ( λ ) = Q c ( λ ) for j = 1 , . . . , m , where Q c ( λ ) ≡ λ − (cid:16) tr ( J c ) − η c τ (cid:17) λ + det ( J c ) (cid:104) η c τ (cid:105) , where η c ≡ πd c D d c + D , (5.0.24)By using (5.0.5) and (5.0.2a), we obtain after a little algebra that (5.0.23) reduces to λ + det ( J c ) λ − λ tr ( J c ) + det J c = − τ πd c (cid:20) d c D (1 + 2 πνβ ) + 2 πmd c | Ω | (1 + τ λ ) (cid:21) , (5.0.25)which is a cubic equation in λ . For this special cell pattern, we conclude that if λ = λ (cid:63) is a root of (5.0.25) for which Q c ( λ (cid:63) ) (cid:54) = 0 , then det M ∞ ( λ (cid:63) ) = 0 . The corresponding eigenvector of M ∞ ( λ (cid:63) ) is the in-phase mode ccc = eee .Next, we will show how to determine roots of the reduced GCEP (5.0.4c) in the case where B ( λ ) = diag ( b ( λ ) , . . . , b m ( λ )) T in (5.0.4b) is not invertible. We first observe that a nontrivial ccc to the leading order problem in (5.0.4c) (with ν = 0 ), existsif and only if there is a λ = λ (cid:63) at which at least two b j cross through zero simultaneously. After relabelling the indices asnecessary, this occurs without loss of generality when there is a λ (cid:63) and an integer J with ≤ J ≤ m for which b ( λ (cid:63) ) = . . . = b J ( λ (cid:63) ) = 0 , with b j ( λ (cid:63) ) (cid:54) = 0 for j = J + 1 , . . . , m if J < m . (5.0.26)Then, for the leading-order problem in (5.0.4c) at λ = λ (cid:63) the nontrivial solutions ccc are ( B ( λ (cid:63) ) + γ ( λ (cid:63) ) E ) ccc = 000 , = ⇒ ccc ∈ C ⊥ ≡ (cid:26) ccc (cid:12)(cid:12)(cid:12) eee T ccc = 0 , ccc = (cid:18) ccc J (cid:19) , ccc J ∈ R J , ∈ R m − J (cid:111) . (5.0.27)Next, we introduce an orthonormal basis for the subspace C ⊥ of dimension J − and we decompose ccc as ccc = ω vvv + . . . + ω J − vvv J − , where C ⊥ ≡ span { vvv , . . . , vvv J − } , vvv Tj vvv i = δ ij , (5.0.28)where δ ij is the Kronecker symbol, and ω j for j = 1 , . . . , J − are scalar coefficients to be found.To determine the O ( ν ) correction to the leading-order eigenvalue λ (cid:63) of the GCEP and identify the constants ω j , we lookfor a nontrivial solution to M ∞ ( λ ) ccc = 000 , as defined in (5.0.4), of the form λ = λ (cid:63) + O ( ν ) . We substitute the expansion λ = λ (cid:63) + ν ˜ λ + · · · , ccc = ccc + νccc + · · · , (5.0.29)into (5.0.4) and then equate O ( ν ) terms to obtain ( B ( λ (cid:63) ) + γ ( λ (cid:63) ) E ) ccc = − ˜ λγ (cid:48) ( λ (cid:63) ) Eccc − ˜ λ B (cid:48) ( λ (cid:63) ) ccc − π G ccc , (5.0.30)where we observe that Eccc = 000 . The solvability condition for (5.0.30) is that the right-hand side of (5.0.30) is orthogonalto each vvv j for j = 1 , . . . , J − . In this way, we readily obtain that ˜ λ and ωωω ≡ ( ω , . . . , ω J − ) T are eigenpairs of the J − dimensional symmetric generalized matrix eigenvalue problem V T G V ωωω = − ˜ λ π V T B (cid:48) ( λ (cid:63) ) V ωωω , V ≡ (cid:0) vvv , . . . , vvv J − (cid:1) . (5.0.31)Here V is the m × J − matrix whose columns provide an orthonormal basis for C ⊥ . In summary, for any such ˜ λ satisfying(5.0.31) a two-term expansion for a root of the reduced GCEP det M ∞ ( λ ) = 0 is λ = λ (cid:63) + ν ˜ λ where λ (cid:63) satisfies (5.0.26).We now illustrate this theory for the special case where J = 2 . This analysis, given below, will be shown in §5.1.1 tobe relevant for analyzing anti-phase instabilities associated with the cell configuration of Fig. 23b where a pair of isolatedidentical cells is spatially segregated from two symmetric ring clusters. Suppose that cells 1 and 2 have common permeabilities d c ≡ d = d , d c ≡ d = d and that they have the same intracellular steady-states. Then, from (5.0.2) we obtaintr J c ≡ tr J = tr J , det J c ≡ det J = det J , and so (5.0.5) and (5.0.2a) yields that b c ( λ ) ≡ b ( λ ) = b ( λ ) . Then, in (5.0.31) wecan take V = (cid:0) / √ , − / √ (cid:1) T , and readily calculate that ˜ λ = − πR c /b (cid:48) c ( λ (cid:63) ) , where b c ( λ (cid:63) ) = 0 . Here R c ≡ R ( xxx ) = R ( xxx ) is the common value of the regular part of the Neumann Green’s function at the centers xxx and xxx of the two cells. For26 (cid:28) , a simple perturbation argument shows that λ ∼ λ (cid:63) + ν ˜ λ can be identified as the root to b c ( λ ) + 2 πνR c = 0 . Finally,by using (5.0.5) for b c , together with (5.0.2), we conclude for ν (cid:28) that λ must be a root of the quadratic λ − λ (cid:18) tr J c − η c τ f (cid:19) + (cid:18) η c τ f (cid:19) det J c = 0 , where η c ≡ πd c D d c + D , f ≡ πνd c d c + D R c . (5.0.32)Since det J c > , an anti-phase instability occurs for cells 1 and 2, with the other cells remaining quiescent, only whentr ( J c ) − η c τ f > . (5.0.33)Although this criterion gives a region in the ( D , τ ) parameter plane, in §5.1.1 we will only implement it for the cellconfiguration in Fig. 23b at a fixed value of τ in order to determine a threshold in D . m = 10 cells in the unit disk We now apply our simplified theory for the regime D = O ( ν − ) to a population of m = 10 cells in the unit disk. Differentspatial configurations of these cells are considered, and we will focus on three different scenarios: (a) all cells are identical, (b)some groups of cells are identical and, (c) none of the cells are identical. For all the examples considered in this subsection,the cells have a common radius ε = 0 . and common Sel’kov kinetic parameters as given in (4.3.1). The heterogeneity inthe cells is introduced through the cell locations and their permeability parameters d j , j = 1 , . . . , m , which specifies therate of feedback of the bulk signal into the cells. The secretion rate is fixed at d = 0 . for all the cells. For each spatialconfiguration of cells and permeability parameter set d , . . . , d m we will compute the HB bifurcation boundary in the τ versus D plane by solving (5.0.22) numerically using Newton’s method with arclength continuation in D .In Fig. 17 we plot the HB boundaries (left panel) in the τ versus D parameter plane together with the cell pattern(right panel) for a pattern with two clusters of cells (top row) and a pattern with an arbitrary arrangement of cells (bottomrow). The precise locations and influx rate d j for the cells for these two specific configurations are given in Table 4 ofAppendix B. By using a numerical winding number computation we have verified that there are exactly two roots of (5.0.19),corresponding to two unstable eigenvalues of the GCEP matrix, inside the lobes spanned by the HB boundaries. (a) HB boundaries (b) Two spatially segregated clusters of cells (c) HB boundaries (d) Arbitrary (random) arrangement of cellsFigure 17: HB boundaries in the τ versus D plane for m = 10 cells with two groups/clusters of cells (top row) and arbitrarily placedcells (bottom row). The dashed curve corresponds to identical cells where d = 0 . for each cell. The thin solid curve is for when thereare two groups of identical cells with d = 0 . for the first group and d = 0 . for the second group (permeability set II). The heavysolid curve is for non-identical cells with d uniformly selected from the interval . ≤ d ≤ . (permeability set III). The locationsand influx rates d j for the cells are given in Table 4 of Appendix B. The remaining parameters are as given in (4.3.1). By comparing the HB boundaries in Fig. 17a and Fig. 17c we observe from the dashed lines in these figures that when thecells are all identical (with d = 0 . for each cell) the unbounded region of instability in the ( D , τ ) parameter plane is very27imilar for the two cell patterns. Therefore, when D is large and when there are a sufficiently large number of identical cells,this observation suggests that the parameter region where intracellular oscillations occur should depend only weakly on thespatial configuration of cells. In addition, by comparing the heavy solid curves in Fig. 17a and Fig. 17c, we conclude that theHB boundaries for cell patterns where d is uniformly selected from the interval . ≤ d ≤ . (permeability set III) alsodepend only weakly on the spatial arrangement of the cells, but that there are no intracellular oscillations when D is large.The influx rates d j for the cells are given in Table 4 of Appendix B. In this case, the bounded instability region in the ( D , τ ) plane exists only within a rather narrow range of τ , which measures the ratio of the time-scale for the reaction kinetics tothe decay rate of the bulk signal. From the thin solid curve in Fig. 17a, corresponding to where we assign the different influxrates d = 0 . and d = 0 . to all the cells in the two different groups (permeability set II in Table 4 of Appendix B), weobserve that the parameter region where intracellular oscillations occur is rather small in D but has a larger extent in τ .Even when D is small, identical cells within a cluster are in close proximity and so are able to synchronize their activitiesand generate oscillations within their respective group. However, as D increases, the bulk signal diffuses rapidly away fromeach of the two clusters and this signal no longer has a sufficient spatial gradient to coordinate synchronized oscillationsbetween the two spatially segregated groups of cells. For an arbitrary arrangement of cells, we observe from the thin solidcurve in Fig. 17c that when cells are assigned either d = 0 . or d = 0 . in such a way that two neighboring cells are notidentical, intracellular oscillations only occur when D is significantly smaller than for the case when the cells with the sameinflux rates are clustered. Qualitatively, this indicates that identical cells within a group of more closely spaced cells canmore readily synchronize their activity.For the case of two clustered groups of cells (see Fig. 17b), we will validate our linear stability predictions with full numericalresults computed from the coupled PDE-ODE model (1.0.2) using FlexPDE [14] at the indicated points in the HB phasediagram in Fig. 17a and for specific permeability sets. In Table 3 we show the real and imaginary parts of the componentsof the normalized eigenvectors ccc and K ccc , together with the unique complex conjugate pair of unstable eigenvalues of the × GCEP matrix M ( λ ) in (2.2.13a), as computed from the root-finding condition det M ( λ ) = 0 , at the red, blue, andgreen dots shown in the phase diagram in Fig. 17a. In Table 4 of Appendix B we indicate the specific permeability set forthe influx rate used at these three pairs of ( D , τ ) . From (2.2.17), the components of ccc measure the diffusive flux into thecells, while ˜ ccc ≡ K ccc , with K given in (5.0.2), predicts the relative amplitude and phase shifts of the intracellular oscillations.In the top row of Fig. 18 we show FlexPDE simulation results of (1.0.2) corresponding to the red dot at ( D , τ ) = (5 , . in the phase diagram of Fig. 17a for the case where the cells are all identical with d j = 0 . for j = 1 , . . . , . The eigenpair ccc and K ccc of the GCEP matrix is given in the top third of Table 3. As predicted by the eigenvector K ccc of the linearizedtheory, the intracellular oscillations observed in the full simulations are nearly synchronized both in amplitude and phase.In the bottom row of Fig. 18 we show, for this parameter set, that results computed from the ODE system (3.0.16). Weobserve that the amplitude and period of intracellular oscillations predicted by the ODEs (3.0.16) compare very favorablywith corresponding FlexPDE results computed from (1.0.2). To explain this very favorable comparison, we observe from thesurface plot in Fig. 18a that the bulk signal is roughly spatially uniform when D = 5 . Recall that the asymptotic analysisfor the derivation of the ODE system (3.0.16) in §3 relies on a nearly spatially uniform bulk signal.In the top row of Fig. 19 we show full FlexPDE results computed from (1.0.2) at the blue and red stars in Fig. 17acorresponding to ( D , τ ) = (5 , . (top left panel) and ( D , τ ) = (5 , . (top right panel), respectively. In this case, wherethe cells are all identical with d = 0 . , the linear stability analysis predicts that the steady-state is linearly stable and that nosustained intracellular oscillations should occur. This prediction is confirmed from the FlexPDE simulations. In the bottompanels of Fig. 19 we show that the corresponding results predicted by the ODE system (3.0.16) compare very favorably withthe FlexPDE results with regards to the long time limiting behavior. At the red star point in Fig. 17a we calculate fromthe root-finding condition det M ( λ ) = 0 , where M is given in (2.2.13a), that the eigenvalue nearest the origin in Re ( λ ) < is λ ≈ − .
01 + 0 . i . This predicts a monotone, non-oscillatory, decay to the steady-state. This feature is observedin the right panels of Fig. 19. Alternatively, at the blue star point in Fig. 19 we calculate that the nearest eigenvalue inRe ( λ ) < for the GCEP matrix at ( D , τ ) = (5 . , . is λ ≈ − .
030 + 0 . i . This eigenvalue predicts an oscillatory decayto the steady-state, and is confirmed by the results in the left panels of Fig. 19. To further explain this transition betweenmonotone and oscillatory decay to the steady-state, in Fig. 20 we plot the real and imaginary parts of the eigenvalue withthe largest real part of the GCEP matrix along the vertical slice D = 5 for . ≤ τ ≤ . in the phase diagram in Fig. 17a,as obtained by numerically solving det M ( λ ) = 0 . This figure shows that the imaginary part of this eigenvalue becomessmall when τ decreases below the lower HB boundary.In the top and middle rows of Fig. 21 we show FlexPDE simulation results corresponding to the blue dot in Fig. 17a where ( D , τ ) = (0 . , . . For this case, the influx rate for the cells in the right and left clusters were assigned as d = 0 . and d = 0 . , respectively (parameter set II in Table 4 of Appendix B). The corresponding normalized eigenpair K ccc obtainedfrom the GCEP matrix, as given in the middle third of Table 3, predicts that the cells will synchronize their oscillationswithin their respective groups, but that there will be a slight phase difference in the intracellular oscillations between thetwo groups. The results from the full FlexPDE simulations in the middle row of Fig. 21 confirm these predictions from thelinearized theory. However, from comparing the middle and bottom rows of Fig. 21, we observe that results from the ODEsystem (3.0.16) do not compare as favorably with FlexPDE simulations as they do in Fig. 18 when D = 5 . This pooreragreement is likely due to the fact that there is a noticeable spatial gradient in the bulk signal at the lower value D = 0 . ,28 a) Surface plot at t = 400
650 660 670 680 690 7000.20.30.40.50.6 (b) U at xxx = (0 , . (FlexPDE)
650 660 670 680 690 70011.52650 660 670 680 690 7000.40.6 (c) u and u (FlexPDE)
650 660 670 680 690 7000.20.30.40.50.6 (d) ¯ U (ODEs)
650 660 670 680 690 70011.52 (e) u (ODEs)
650 660 670 680 690 7000.40.6 (f) u (ODEs)Figure 18: Top row: FlexPDE numerical results computed from the PDE-ODE model (1.0.2) for m = 10 identical cells, arrangedin two clusters (see Fig. 17b), at the red dot in Fig. 17a where ( D , τ ) = (5 . , . . The cells have identical influx rates d j = 0 . for j = 1 , . . . , m (set I) and almost identical intracellular dynamics. The cell locations are in Table 4 of Appendix B. Lower row:corresponding results for ¯ U , u , and u , as computed from the ODE system (3.0.16). The eigenvector and eigenvalue for the GCEPmatrix for the linearization are in the top third of Table 3. The results from the ODEs compare well with the FlexPDE simulations. (a) U , u , u at ( D , τ ) = (5 . , . (FlexPDE) (b) U , u , u at ( D , τ ) = (5 . , . (FlexPDE) (c) ¯ U , u , u at ( D , τ ) = (5 . , . (ODEs) (d) ¯ U , u , u at ( D , τ ) = (5 . , . (ODEs)Figure 19: Top row: FlexPDE numerical results computed from the PDE-ODE model (1.0.2) for m = 10 identical cells, arrangedin two clusters (see Fig. 17b) at the blue star in Fig. 17a where ( D , τ ) = (5 . , . (left panel) and at the red star in Fig. 17awhere ( D , τ ) = (5 . , . (right panel). The cells have identical influx rates d j = 0 . for j = 1 , . . . , (set I) and almost identicalintracellular dynamics. The cell locations are in Table 4 of Appendix B. Lower row: corresponding results for ¯ U , u , and u , ascomputed from the ODE system (3.0.16). As expected, there are no sustained oscillations and the steady-state is stable. ( D , τ ) ; D = D /ν ( Re λ, Im λ ) Cell j (cid:16) Re ( c j ) , Im ( c j ) (cid:17) θ j ( rad ) (cid:16) Re (˜ c j ) , Im (˜ c j ) (cid:17) (0 . , . . − . , . (0 . , . . − . , . (0 . , . . − . , . (0 . , − . .
284 ( − . , . Identical cells (5 , .
3) (0 . , . (0 . , − . .
28 ( − . , . (Set I) D ≈ . (0 . , . . − . , . Red dot 7 (0 . , − . .
28 ( − . , . (0 . , . . − . , . (0 . , − . .
28 ( − . , . (0 . , − . .
28 ( − . , . (0 . , . . − . , . (0 . , − . .
17 ( − . , . (0 . , . . − . , . (0 . , − . .
10 ( − . , . Two groups (0 . , .
35) (0 . , . (0 . , − . .
06 ( − . , . (Set II) D ≈ . (0 . , − . .
23 ( − . , . Blue dot 7 (0 . , − . .
19 ( − . , . (0 . , − . .
23 ( − . , . (0 . , − . .
16 ( − . , . (0 . , − . .
17 ( − . , . (0 . , . . − . , . (0 . , . .
175 ( − . , − . (0 . , . . − . , − . (0 . , . .
331 ( − . , − . Random (1 . , .
3) (0 . , . (0 . , . .
30 ( − . , − . (Set III) D ≈ . (0 . , . .
307 ( − . , − . Green dot 7 (0 . , . .
630 ( − . , − . (0 . , . . − . , − . (0 . , . .
208 ( − . , − . ( − . , . .
12 ( − . , − . Table 3: Real and imaginary parts of the normalized eigenvector ccc of the GCEP matrix (2.2.13a), together with ˜ ccc ≡ K ccc , computed atthe red, blue, and green dot shown in the phase diagram in Fig. 17a for a two-cluster arrangement of cells shown in Fig. 17b. The celllocations and permeability sets for the influx rate are given in Table 4 of Appendix B. The other parameters are as given in (4.3.1).The third column gives the unique unstable eigenvalue in Re ( λ ) > of the GCEP matrix at these three pairs of ( D , τ ) . Figure 20: Real (blue curve; left y axis) and imaginary (red curve; right y axis) parts of the eigenvalue of the GCEP matrix withthe largest real part versus τ along the vertical slice D = 5 for . ≤ τ ≤ . in the phase diagram in Fig. 17a. The eigenvalue iscomputed using root-finding on det M ( λ ) = 0 , where M ( λ ) is given in (2.2.13a). Observe the two HB values in τ where Re ( λ ) = 0 (dotted black line). When τ = 0 . , we have Im ( λ ) ≈ . . However, as τ decreases below the lower HB point, Im ( λ ) decreases to zero. as observed from the surface plot of Fig. 21a. Moreover, Fig. 21a shows that the bulk signal can concentrate, as expected,in the left cluster owing to the lower rate of chemical influx into the cells in this cluster.Next, we give FlexPDE results computed from (1.0.2) corresponding to the green dot in Fig. 17a where ( D , τ ) = (1 . , . ,corresponding to the two cluster cell pattern of Fig. 17b). For this case, the influx rates are selected randomly from theinterval . ≤ d ≤ . and are given in Table 4 under permeability set III. In the top and middle row of Fig. 22 we showthe FlexPDE results for the bulk solution at the point (0 , . as well as the intracellular dynamics in cells 1, 2, 6 and 7.30 a) Surface plot at t = 400
650 660 670 680 690 7000.20.30.40.50.60.7
Bulk concentration (b) U at xxx = (0 , . (FlexPDE)
650 660 670 680 690 70011.52
Left Cluster
650 660 670 680 690 7000.60.8 (c) u , u Left cluster (FlexPDE)
650 660 670 680 690 70012
Right Cluster
650 660 670 680 690 7000.40.50.6 (d) u , u Right cluster (FlexPDE)
650 660 670 680 690 7000.20.30.40.50.60.7 (e) ¯ U (ODEs)
650 660 670 680 690 70011.52650 660 670 680 690 7000.60.8 (f) u , u Left cluster (ODEs)
650 660 670 680 690 70011.52650 660 670 680 690 7000.40.50.6 (g) u , u Right cluster (ODEs)Figure 21: Top and middle row: FlexPDE numerical results computed from the PDE-ODE model (1.0.2) for m = 10 cells, arrangedin two clusters (see Fig. 17b) at the blue dot in Fig. 17a where ( D , τ ) = (0 . , . . The cells in the left cluster have identical influxrates d = 0 . while the cells in the right cluster have d = 0 . . The cell locations are in Table 4 of Appendix B. Within each clusterthere is very similar intracellular dynamics. Lower row: corresponding results for ¯ U , u , and u , as computed from the ODE system(3.0.16). The eigenvector and eigenvalue for the GCEP matrix for the linearization are in the middle third of Table 3. From Table 4, cells 1 and 2 are from the left cluster, while 6 and 7 belong to the right cluster. The corresponding normalizedeigenpair K ccc of the GCEP matrix, as given in the bottom third of Table 3, predicts that the amplitudes of the intracellulardynamics will differ from cell to cell due to the different cell influx rates, but that there will only be a small phase shift in theintracellular oscillations. Since | ( K ccc ) | < | ( K ccc ) | (see Table 3), the linearized theory predicts larger amplitude oscillations for u in cell 7 than in cell 1, which is expected since cell 7 has a larger influx rate than does cell 1 (see Table 4). This is confirmedin the FlexPDE simulations in Fig. 22. In the bottom row of Fig. 22 we observe that the corresponding results from the ODEsystem (3.0.16) compare moderately well with the FlexPDE results regarding the amplitude and period of oscillations in thebulk medium and within the cells. We remark that if, instead, we used the simpler ODE system (3.0.19), corresponding tothe well-mixed limit D → ∞ , the ODE results would be in very poor agreement with the full PDE simulations. In this subsection, we analyze in detail the role of the spatial configuration of the cells on the triggering of intracellularoscillations. Specifically, we consider the cell configuration shown in Fig. 23b, with cell centers given in Table 5 of Appendix B.This pattern consists of two spatially segregated rings of cells together with two cells that are spatially isolated from therings. For this symmetric pattern we will consider five permeability parameter sets for the cell influx rate, as given in Table 5,which lead to solution behavior that can be interpreted both qualitatively and from our linear stability analysis.In Fig. 23a we plot the HB boundaries in the τ versus D parameter plane for permeability sets II − V of Table 5, asobtained by solving for the roots of (5.0.22) numerically. There is no HB boundary in this parameter plane from (5.0.22) for31 a) surface plot at t = 400
650 660 670 680 690 7000.350.40.450.50.550.6
Bulk concentration (b) U at xxx = (0 . . (FlexPDE)
650 660 670 680 690 70011.5
Cell 1
650 660 670 680 690 7000.80.91 (c) u , u Cell 1 (FlexPDE)
650 660 670 680 690 70011.52
Cell 2
650 660 670 680 690 7000.60.70.8 (d) u , u Cell 2 (FlexPDE)
650 660 670 680 690 7001.52
Cell 6
650 660 670 680 690 7000.50.60.7 (e) u , u Cell 6 (FlexPDE)
650 660 670 680 690 7001.52
Cell 7
650 660 670 680 690 7000.40.50.6 (f) u , u Cell 7 (FlexPDE)
650 660 670 680 690 7000.350.40.450.50.550.6 (g) ¯ U (ODEs)
650 660 670 680 690 70011.52 (h) u , Cells 1,2,6,7 (ODEs)
650 660 670 680 690 7000.40.60.81 (i) u , Cells 1,2,6,7 (ODEs)Figure 22: Top and middle row: FlexPDE numerical results computed from the PDE-ODE model (1.0.2) for m = 10 cells, arrangedin two clusters (see Fig. 17b) at the green dot in Fig. 17a where ( D , τ ) = (1 . , . . The influx rates are selected uniformly from theinterval . ≤ d ≤ . (see Table 4 for set III). The cell locations are in Table 4 of Appendix B, with Cell 1 at (0 . , − . , Cell 2at (0 . , . , Cell 6 at ( − . , . , and Cell 7 at ( − . , . . Each cell has different intracellular dynamics owing to the differentinflux rates. Lower row: corresponding results for ¯ U , u , and u , as computed from the ODE system (3.0.16). The eigenvector andeigenvalue for the GCEP matrix for the linearization are in the lower third of Table 3. permeability set I, where d = 0 . for the ring cells and d = 0 . for the isolated cells. As verified from a numerical windingnumber computation, only within the lobes spanned by the HB boundaries is the steady-state unstable to an oscillatoryinstability. From the dotted and dashed-dotted curves in Fig. 23a corresponding to where all the cells have the commoninflux rates d = 0 . or d = 0 . , respectively, we observe that the instability lobe is unbounded in D . Therefore, when thecells are all identical, intracellular oscillations occur within some finite band of the reaction-time parameter τ for all bulkdiffusivities. For identical cells with d = 0 . , the HB boundaries in Fig. 23a are very similar to that shown in Fig. 17 forthe arbitrary cell pattern and the two-cluster pattern. For identical cells with the larger cell influx rate d = 0 . we observefrom Fig. 23a that the lower threshold in τ where oscillatory instabilities first occur is smaller than when the common influxrate is d = 0 . . To explain this qualitatively, we observe that τ decreases as the bulk decay rate k B increases (see (1.0.3)).With a larger bulk decay rate, the bulk signal becomes less diffuse and has larger spatial gradients, which leads to a buildupof the bulk signal, near the cells. As a result, when there is a larger cell influx rate, the bulk signal can more readily enterthe cell to trigger intracellular oscillations.Since from Fig. 23 the steady-state is always linearly stable for permeability set I when D = D /ν , this motivates studyingthe full GCEP matrix M ( λ ) in (2.2.13) that is valid for D = O (1) . In the left and right panels of Fig. 24 we plot spectralinformation versus D , for fixed τ = 0 . , obtained from the roots of det M ( λ ) = 0 for the two-ring and isolated cell pattern forpermeability sets I (left panel) and II (right panel). In the left y -axes of Fig. 24 we plot the real parts of the three eigenvaluesof the GCEP matrix that have the largest real parts. These three eigenvalue branches are associated with different spatialmodes for the cells as indicated in the figure legends and discussed below. In the right y -axes of Fig. 24 we plot the sum32 .2 0.4 0.6 0.8 100.20.40.60.811.21.41.6 (a) HB boundaries (b) Two rings of cells with two isolated cellsFigure 23: HB boundaries (left panel) for m = 10 cells for a pattern with two rings of cells and two isolated cells (right panel). The celllocations and permeability parameter sets for the influx rates d j are given in Table 5 of Appendix B. Four out of five of the permeabilitysets in Table 5 give HB boundaries in the ( D , τ ) plane. The thin solid curve (set II) has d = 0 . for the ring cells and d = 0 . forthe isolated cells. The heavy solid curve has d = 0 . for the isolated cells, d = 0 . for one group of ring cells, and d = 0 . for thecells on the other ring (set III). The dash-dotted curve has d = 0 . for all cells (set IV), while the dotted curve has d = 0 . for allcells (set V). There is no HB boundary for set I where d = 0 . and d = 0 . for the isolated and ring cells, respectively. The remainingparameters are given in (4.3.1). (a) set I: d = 0 . (rings), d = 0 . (isolated) (b) set II: d = 0 . (isolated), d = 0 . (rings)Figure 24: Spectral information from the eigenvalues λ and normalized eigenvectors ˜ ccc = K ccc , with (cid:80) i | ( K ccc ) i | = 1 , obtained from theGCEP matrix (2.2.13a) and (5.0.2a) for K , as computed using root-finding on det M ( λ ) = 0 . The cell configuration is the two-ring andtwo-isolated cell pattern of Fig. 23b for permeability set I (left panel) and permeability set II (right panel) for the fixed value τ = 0 . .Left y -axis: Real parts (blue curves) of the three eigenvalues of the GCEP matrix with the largest real parts versus D . Right y -axis:The sum | ˜ c | + | ˜ c | (red curves) of the first two components of the normalized eigenvector K ccc , which measures the relative amplitudeof oscillations in the two isolated cells in comparison to the cells on the ring. The linetype (solid, dashed, dot-dashed) of the blue andred curves correspond to the same eigenpair of the GCEP matrix. | ˜ c | + | ˜ c | (red curves) of the first two components of the normalized eigenvector ˜ ccc ≡ K ccc , which measures the relativeamplitude of oscillations in the two isolated cells in comparison to the cells on the ring (see (2.2.17)). The linetype (solid,dashed, dot-dashed) of the blue and red curves in Fig. 24 correspond to the same eigenpair of the GCEP matrix.For permeability set I, where the influx rate d = 0 . on the isolated cells is lower than that on the ring cells ( d = 0 . ), weobserve from the solid and dashed lines in Fig. 24a that there can be two unstable modes when D is small enough. On therange D < . , there is an unstable mode given by the blue dashed curve where anti-phase oscillations occur for the twoisolated cells (cells 1 and 2), which has a much higher amplitude than for the ring cells. On the two rings, cells 4 and 6 aswell as cells 8 and 10 oscillate out of phase. The other ring cells (cells 3, 5, 7, and 9) are essentially quiescent for this mode.Due to the low influx rate into the isolated cells, the isolated cells communicate imperfectly with their images across thedomain boundary when D is small, which leads to the anti-phase instability. Observe that as D increases above D ≈ . this mode becomes stable and the amplitude in the isolated cells decreases (decreasing dashed red curve) since the bulksignal near the isolated cells is washed away. The other possible unstable mode, given by the blue solid curve in Fig. 24a,is one for which the two isolated cells are in-phase. The ring cells, which have smaller amplitude oscillations than do theisolated cells, are all roughly in phase and have only a slight phase difference with the isolated cells. This in-phase modeis unstable only for D < . . We remark that this low stability threshold value of D is consistent with the observationsin Fig. 23a that the steady-state is always linearly stable for D = O ( ν − ) (cid:29) . Finally, the dashed-dotted blue curve in33ig. 24a corresponds to a linearly stable mode where the isolated cells are quiescent, with each ring having cell oscillationsthat are (roughly) synchronized in amplitude and phase. However, for this mode there is a large phase shift between theoscillations in the two rings clusters. With the larger influx rate for the ring cells, there is effective communication betweenthe two spatially segregated rings as D increases, which precludes any anti-phase instability between the two ring clusters. (a) Surface plot at t = 400
650 660 670 680 690 7000.110.120.130.140.15 (b) U at xxx = (0 , (FlexPDE)
650 660 670 680 690 70012
Cells 1 and 2
650 660 670 680 690 7000.40.60.8 (c) u , u isolated cells 1 and 2 (FlexPDE)
650 660 670 680 690 7001.81.851.9
Cells 3 and 4
650 660 670 680 690 7000.4550.460.4650.47 (d) u , u cells 3 and 4 (left ring)
650 660 670 680 690 7001.851.9
Cells 7 and 8
650 660 670 680 690 7000.460.47 (e) u , u cells 7 and 8 (right ring)Figure 25: FlexPDE numerical results computed from the PDE-ODE model (1.0.2) for m = 10 cells, arranged in two rings with twoisolated cells (see Fig. 23b) when D = 0 . and τ = 0 . for the permeability set I and cell locations given in Table 5. Other parametersare given in (4.3.1). The linear stability predictions are given in Fig. 24a.(a) Surface plot at t = 400
650 660 670 680 690 7000.20.220.240.260.280.3 (b) U at xxx = (0 , (FlexPDE)
650 660 670 680 690 70011.52
Cells 1 and 2
650 660 670 680 690 7000.60.70.8 (c) u , u isolated cells 1 and 2 (FlexPDE)
650 660 670 680 690 7001.71.81.9
Cells 3 and 4
650 660 670 680 690 7000.450.5 (d) u , u cells 3 and 4 (left ring)
650 660 670 680 690 7001.71.81.9
Cells 7 and 8
650 660 670 680 690 7000.450.5 (e) u , u cells 7 and 8 (right ring)Figure 26: Same caption as in Fig. 25 except that now D is increased to D = 0 . . Observe that there is large in-phase signalinggradient near the isolated cells in the surface plot (a), and that the isolated cells have now synchronized their dynamics. The oscillationsin the ring cells are more synchronized than when D = 0 . , but still have smaller amplitude for u than do the isolated cells. As a test of the predictions of the linear stability theory, as summarized by Fig. 24a, in Fig. 25 and Fig. 26 we show34lexPDE results for (1.0.2) for permeability set I when D = 0 . and D = 0 . , respectively. For D = 0 . where both thein-phase and anti-phase modes (solid and dashed curves in Fig. 24a) are unstable with comparable growth rates, we observefrom Figs. 25c–25e that, as predicted, the intracellular dynamics in cells 1 and 2 are not in-phase and that the amplitudeof oscillations in the ring cells is much smaller than in the isolated cells. The strong anti-phase bulk signaling gradient atthe isolated cells, as shown in Fig. 25a, is also consistent with the linear stability theory. In contrast, when D = 0 . ,the in-phase mode is the dominant instability as seen from Fig. 24a. For this larger value of D , we observe from Fig. 26aand Fig. 26c that the bulk signaling gradient and the intracellular oscillations are now in-phase at the two isolated cells.Moreover, by comparing Figs. 25 and 26, we observe that the oscillations within the ring cells are more synchronized andhave a larger amplitude, while the isolated cells have a smaller amplitude, when D = 0 . as compared to when D = 0 . .These observations are all consistent with the linear stability predictions shown in Fig. 24a. (a) Surface plot at t = 400
650 660 670 680 690 7000.150.20.250.30.35 (b) U at xxx = (0 , (FlexPDE)
650 660 670 680 690 7001.61.82
Cells 1 and 2
650 660 670 680 690 7000.450.5 (c) u , u isolated cells 1 and 2 (FlexPDE)
650 660 670 680 690 70011.52
Cells 3 and 4
650 660 670 680 690 7000.40.6 (d) u , u cells 3 and 4 (left ring)
650 660 670 680 690 70011.52
Cells 7 and 8
650 660 670 680 690 7000.40.6 (e) u , u cells 7 and 8 (right ring)Figure 27: FlexPDE numerical results computed from the PDE-ODE model (1.0.2) for m = 10 cells, arranged in two rings with twoisolated cells (see Fig. 23b) when D = 0 . and τ = 0 . for the permeability set II and cell locations given in Table 5. Other parametersare given in (4.3.1). The linear stability predictions are given in Fig. 24b. In Fig. 24b we show the corresponding spectral plot, as computed from the GCEP (2.2.13), for permeability set II wherethe isolated cells have a higher influx rate d = 0 . than do the ring cells d = 0 . . In contrast to the case for permeabilityset I, we now observe that in addition to the usual in-phase mode that is unstable for D < . , the only other possibleunstable mode corresponds to an anti-phase instability between the two ring clusters. Since the influx rate into the ring cellshas been decreased, this unstable mode is due to a relatively poorer communication between the two rings when the bulkdiffusivity is small. This anti-phase cluster mode is unstable for D < . . In contrast, the anti-phase mode for the isolatedcells, which was unstable for permeability set I when D is small, is now always linearly stable for set II. Another key featurefrom Fig. 24b is that, for the in-phase mode, the intracellular oscillations in the isolated cells are much smaller when D issmall than that in Fig. 24b. However, we observe from the increasing solid red curve in Fig. 24b that the amplitude of theoscillations in the isolated cells grows as D increases. The interpretation of this observation is that the secretion from thering cells can more easily diffuse across the domain to the isolated cells, where it is readily absorbed, as D is increased. InFig. 27 we show results from FlexPDE simulations of (1.0.2) for permeability set II when D = 0 . . The strong signallinggradient near the ring cells where the influx rate is small, the smaller amplitude oscillations in the isolated cells than in therings, and the slight phase shift between the oscillations in the two ring clusters (cells 3 and 4 versus cells 7 and 8) are allconsistent with the predictions in Fig. 24b of the linear stability theory.Finally, to obtain some analytical insight into the possibility of an unstable mode consisting of anti-phase oscillationsin the isolated cells together with quiescent ring cells, we implement the degenerate perturbation theory on the reducedGCEP matrix (5.0.4) from the D = D /ν regime. By using the simple criterion in (5.0.33), based on the quadratic equation(5.0.32), we calculate for permeability set I that this anti-phase mode for the isolated cells is unstable when τ = 0 . onlywhen D < . . This corresponds to D < . using D = D /ν with ν = − / log ε and ε = 0 . . Although this simplecriterion (5.0.33) predicts the existence of an unstable anti-phase mode for the isolated cells, the threshold value is not soaccurate (see Fig. 24a) owing to the fact that D is not asymptotically large. In contrast, for permeability set II, we calculatefrom the criterion (5.0.33) that the anti-phase mode for the isolated cells is always stable, as is consistent with Fig. 24b.35 Discussion and Outlook
We have analyzed diffusion-sensing behavior for the coupled cell-bulk PDE-ODE system (1.0.2) that was used to modelthe switch-like onset of intracellular oscillations for a collection of heterogeneous cells, as mediated through a passive bulkdiffusion field with finite bulk diffusivity D . For the case of Sel’kov reaction kinetics, we have studied how the onset ofintracellular oscillations depends on the spatial configuration of cells, on the membrane permeability parameters, and onthe triggering effect of a single “defective” cell with a different kinetic parameter. In contrast, in [26] QS behavior resultingfrom increase in cell density were studied using only the ODE system (3.0.19) with global coupling, which pertains to thewell-mixed regime of infinite bulk diffusivity.There are several open theoretical and computational challenges that should be explored for (1.0.2). One key numericalissue concerns developing efficient and well-conditioned numerical techniques to implement the full linear stability theorybased on the root-finding condition det ( M ( λ )) = 0 for the GCEP given in (2.2.13) for a large number, i.e. m ≥ , ofrandomly located cells with arbitrary permeabilities. The solution strategies for such nonlinear matrix eigenvalue problemsare typically restricted to matrices wih special structure, such as Hermitian matrices, matrices with low-rank dependenceon λ , or matrices that are quadratic or rational in the eigenvalue parameter λ (cf. [23] and [4]). In contrast, the GCEPmatrix in (2.2.13) is not Hermitian when λ is complex-valued, and is not of low-rank owing to its dependence on the fulleigenvalue-dependent Green’s matrix G λ .From a mathematical viewpoint, a second open direction would be to develop a weakly nonlinear theory for the Hopfbifurcations of steady-state solutions of (1.0.2). The numerical results shown in §4 and 5 have suggested that the intracellularoscillations are supercritical for the Sel’kov kinetics, and that their relative amplitude and phases within the cells are well-predicted by the eigenvector of the GCEP matrix M ( λ ) in (2.2.13). Also of interest would be to explore whether ODEphase-reduction techniques, such as surveyed in [41] and [36], can be extended to the PDE-ODE setting of (1.0.2).From a modeling perspective, the PDE-ODE system (1.0.1) is well-suited for investigating triggered intracellular oscillationsand collective dynamics associated with specific microbial systems for which detailed models of the signaling pathways areavailable, the autoinducer is known, and where membrane permeabilities can be estimated from biological data. Suchcalibrated biological models are readily incorporated into the theoretical framework (1.0.1). In particular, it would beworthwhile to analyze (1.0.1) with the detailed model of [56] and [25] for the glycolytic pathway in yeast cells and the modelof [33] for bacterial communication of Escherichia coli cells. For non-oscillatory QS systems, it would be interesting to usethe intracellular Lux-signaling pathways of [32] into our PDE-ODE model (1.0.1) to analyze sudden jumps between smalland large amplitude steady-states for bistable QS systems as the cell density increases past a saddle-node bifurcation value.Moreover, (1.0.1) can be used to study the effect of spatial diffusion on the parallel QS signaling pathways associated withthe marine bacterium
V. harveyi , which were recently modeled in [5] through an ODE well-mixed limit. Our analysis in §3has shown that the dimensionless PDE-ODE system (1.0.2) can be reduced to a limiting ODE system only when the bulkdiffusivity satisfies D = O ( − log ε ) , where ε is the dimensionless ratio of the cell radius to the length scale of the domain.The classic ODE system for the well-mixed limit corresponds to the regime D (cid:29) O ( − log ε ) .Finally, it would be worthwhile to extend the 2-D model (1.0.1) to allow for two passive bulk diffusion fields, and toconsider 3-D domains, where the cell-cell interaction is weaker than in 2-D owing to the rapid decay of the 3-D free-spaceGreen’s function. By including two bulk species, and with a re-formulation of the class of lattice-based models introducedin [42] into the framework of (1.0.1), it should be possible to show analytically that changes in the permeability parametersof a collection of cells can induce a Turing or transcritical bifurcation to a patterned steady-state for an activator-inhibitorsystem even when the two bulk diffusing species have very similar diffusivities. In a spatially homogeneous medium withoutcells, Turing bifurcations only occur for activator inhibitor systems when there is a sufficiently large diffusivity ratio, whichis not typical for most biological systems. The PDE-ODE framework of (1.0.1), where the small signaling compartmentsare modeled explicitly, also provides an alternative approach for studying large-scale self-organized structures that havepreviously been modeled through PDE-ODE lattice dynamics in which the dynamically active “cells” are treated as pointsources restricted to lattice sites and where a discrete Laplacian averaging over nearby lattice points replaces the Laplaceoperator (cf. [6], [43], [29]). Such lattice PDE-ODE systems have been used to model traveling waves of yeast activation dueto substrate addition of glucose at a localized source (cf. [43]), and the emergence of spiral waves resulting from the couplingof Fitzhugh-Nagumo or Rossler cell kinetics to discrete lattice-based diffusion (cf. [29], [6]). Acknowledgements
Michael Ward gratefully acknowledges the financial support from the NSERC Discovery grant program. We are grateful toJustin Tzou (Macquarie U.) for his initial help with the FlexPDE simulations.36 ppendices
A The reduced-wave Green’s function for the unit disk
In this appendix we provide explicit expressions for the matrix spectrum of (4.1.9) corresponding to the matrix block of theeigenvalue-dependent reduced-wave Green’s matrix in (4.1.8) representing the interactions of m − ≥ cells on a ring ofradius r , with < r < , concentric within the unit disk. For notational convenience we define N by N = m − . The N cells on the ring of radius r are centered at xxx k = r (cos ψ k , sin ψ k ) T , where ψ k ≡ π ( k − /N for k = 1 , . . . , N . In the unitdisk, the reduced-wave Green’s function G λ ( xxx ; ξξξ ) and its regular part, satisfying (2.2.8), can be calculated using separationof variables as (see equations (6.10) and (6.11) of [20]) G λ ( xxx ; ξξξ ) = 12 π K ( ϕ λ | xxx − ξξξ | ) − π ∞ (cid:88) n =0 β n cos ( n ( ψ − ψ )) K (cid:48) n ( ϕ λ ) I (cid:48) n ( ϕ λ ) I n ( ϕ λ r ) I n ( ϕ λ ρ ) , (A.1a) R λ ( ξξξ ) = 12 π (cid:20) log (cid:16) √ D (cid:17) − γ e −
12 log (1 + τ λ ) (cid:21) − π ∞ (cid:88) n =0 β n K (cid:48) n ( ϕ λ ) I (cid:48) n ( ϕ λ ) [ I n ( ϕ λ ρ )] , (A.1b)where ϕ λ ≡ (cid:112) (1 + τ λ ) /D is the principal branch of ϕ λ , γ e = 0 . is Euler’s constant, and I n ( z ) and K n ( z ) are themodified Bessel functions of the first and second kind of order n , respectively. In (A.1), β ≡ , β n ≡ for n ≥ , while xxx ≡ r (cos ψ, sin ψ ) T , and ξξξ ≡ ρ (cos ψ , sin ψ ) T .For a ring pattern, the N × N symmetric Green’s matrix G λN is also cyclic and can be generated by a cyclic permutationof its first row aaa λ ≡ ( a λ, , . . . , a λ,N ) T , which is defined term-wise by a λ, ≡ R λ ( xxx ) ; a λ,k = G λ ( xxx k ; xxx ) , k = 2 , . . . , N . (A.2)The matrix spectrum G λN vvv j = ω λj vvv j in (4.1.9) with N = m − is calculated as in §6 of [20]. The in-phase eigenpair is ω λN = N (cid:88) k =1 a λ,k , vvv N = (1 , . . . , T , (A.3a)while the other eigenvalues, corresponding to the anti-phase modes for which vvv Tj vvv N = 0 for j = 1 , . . . , N − , are ω λj = N − (cid:88) k =0 cos (cid:18) πjkN (cid:19) a λ,k +1 , j = 1 , . . . , N − . (A.3b)Since ω λj = ω λ,N − j for j = 1 , . . . , (cid:100) N/ (cid:101) − , there are (cid:100) N/ (cid:101) − pairs of degenerate eigenvalues for G λN . Here the ceilingfunction (cid:100) x (cid:101) is defined as the smallest integer not less than x . When N is even, there is an eigenvalue of multiplicity onegiven by ω λ, N = (cid:80) N − k =0 ( − k a λ,k +1 . The other eigenvectors for j = 1 , . . . , (cid:100) N/ − (cid:101) are vvv j = (cid:18) , cos (cid:18) πjN (cid:19) , . . . , cos (cid:18) πj ( N − N (cid:19) (cid:19) T , vvv N − j = (cid:18) , sin (cid:18) πjN (cid:19) , . . . , sin (cid:18) πj ( N − N (cid:19) (cid:19) T . (A.3c)Finally, when N is even, there is an additional eigenvector vvv N = (1 , − , . . . , − T . By using the explicit formulae for G λ andits regular part from (A.1), the eigenvalues ω λj of the Green’s matrix G λN are then easily computed from (A.3). Remark 1.
The symmetric and cyclic matrix G λN has N/ distinct anti-phase modes if N is even and ( N − / distinctanti-phase modes if N is odd. Finally, in order to calculate the quantities in (4.1.2) and (4.1.3), which are needed in (4.1.4) for determining the steady-state solution, we need only set λ = 0 in (A.1) and (A.3a). B Cell locations and permeability parameters
The cell locations and permeability parameter sets for the influx rate d j , for j = 1 , . . . , m , are given in the next two tablesfor the cell patterns studied in §5. References [1] Continuation Test: a MATLAB library which defines test functions for continuation codes. http://people.math.sc.edu/Burkardt/m_src/test_con/test_con.html . Accessed: 2020-02-26.37wo clusters of cells Arbitrary cell locationsCell i x i y i d i (II) d i (III) x i y i d i (II) d i (III) . − . . . . − . . . . . . . − . − . . . . − . . . . − . . . . . . . − . . . . . − . . . . . . . − . . . . − . . . . − . . . . − . − . . . − . − . . . . . . . − . . . . . . . . − . − . . . − . − . . . Table 4: Locations of the cell centers and influx permeability rates for parameter sets II and III corresponding to the two clusterarrangement of cells in Fig. 17b (columns 2–5) and the arbitrary arrangement of cells in Fig. 17d (columns 6–9). The permeabilityparameter set I (not shown) is for identical cells where d j = 0 . for j = 1 , . . . , m . Two rings with two isolated cellsCell i x i y i d i (I) d i (II) d i (III) d i (IV) d i (V) . − . . . . . .
32 0 . . . . . . . − . . . . . . . − . . . . . . . − . . . . . . . − . − . . . . . .
37 0 . . . . . . .
38 0 . . . . . . .
39 0 . . . . . . .
310 0 . − . . . . . . Table 5: Locations of the cell centers and influx permeability rates for parameter sets I–V corresponding to the pattern of two rings ofcells together with two isolated cells, as shown in Fig. 23. [2] D. Alciatore and R. Miranda. A winding number and point-in-polygon algorithm.
Glaxo Virtual Anatomy ProjectResearch Report, Department of Mechanical Engineering, Colorado State University , 1995.[3] T. Betcke, N. G. Highan, V. Mehrmann, G. M. N. Porzio, C. Schröder, and Tisseur F. NLEVP: A collection ofnonlinear eigenvalue problems. users’ guide.
MIMS EPring 2011.117, Manchester Institue for Mathematical Sciences,The University of Manchester, UK , page 10 pages, updated 2019.[4] T. Betcke, N. G. Highan, V. Mehrmann, C. Schröder, and Tisseur F. NLEVP: A collection of nonlinear eigenvalueproblems.
ACM Trans. Math. Software , 39(2):7.1–7.28, 2013.[5] P. Bressloff. Ultrasensitivity and noise amplification in a model of v. harveyi quorum sensing.
Phys. Rev. E , 93:062418,2016.[6] X. Z. Cao, H. Yuan, and B. W. Li. Selection of spatiotemporal patterns in arrays of spatially distributed oscillatorsindirectly coupled via a diffusive environment.
Chaos , 29:043104, 2019.[7] M. A. J. Chaplain, M. Ptashnyk, and M. Sturrock. Hopf bifurcation in a gene regulatory network model: Molecularmovement causes oscillations.
Math. Mod. Meth. Appl. Sci. , 25(6):1179–1215, 2015.[8] S. Danø, M. F. Madsen, and P. G. Sørensen. Quantitative characterization of cell synchronization in yeast.
Proceedingsof the National Academy of Sciences , 104(31):12732–12736, 2007.[9] S. Danø, P. G. Sørensen, and F. Hynne. Sustained oscillations in living cells.
Nature , 402(6759):320–322, 1999.[10] S. De Monte, F. d’Ovidio, S. Danø, and P. G. Sørensen. Dynamical quorum sensing: Population density encoded incellular dynamics.
Proceedings of the National Academy of Sciences , 104(47):18377–18381, 2007.[11] G. E. Dilanji, J. B. Langebrake, P. De Leenheer, and S. J. Hagen. Quorum activation at a distance: Spatiotemporalpatterns of gene regulation from diffusion of an autoinducer signal.
J. Am. Chem. Soc. , 6:34695, 2016.[12] J. D. Dockery and J. P. Keener. A mathematical model for quorum sensing in pseudomonas aeruginosa . Bull MathBiol. , 63(1):95–116, 2001. 3813] G. M. Dunny and B. Leonard. Cell-cell communication in gram-positive bacteria.
Annual review of microbiology ,51(1):527–564, 1997.[14] PDE FlexPDE. Solutions inc. , 2015.[15] M. Gao, H. Zheng, Y. Ren, R. Lou, F. Wu, W. Yu, X. Liu, and X. Ma. A crucial role for spatial distribution in bacterialquorum sensing.
Scientific Reports , 6(34695), 2016.[16] A. Goldbeter.
Biochemical oscillations and cellular rhythms: the molecular bases of periodic and chaotic behaviour .Cambridge university press, 1997.[17] A. Gomez-Marin, J. Garcia-Ojalvo, and J. M. Sancho. Self-sustained spatiotemporal oscillations induced by membrane-bulk coupling.
Phys. Rev. Lett. , 98:168303, Apr 2007.[18] J. Gou, W.-Y. Chiang, P.-Y. Lai, M. J. Ward, and Y.-X. Li. A theory of synchrony by coupling through a diffusivechemical signal.
Physica D , 339:1–17, 2017.[19] J. Gou, Y. X. Li, W. Nagata, and M. J. Ward. Synchronized oscillatory dynamics for a 1-D model of membrane kineticscoupled by linear bulk diffusion.
SIAM J. Appl. Dyn. Sys. , 14(4):2096–2137, 2015.[20] J. Gou and M. J. Ward. An asymptotic analysis of a 2-d model of dynamically active compartments coupled by bulkdiffusion.
Journal of Nonlinear Science , 26(4):979–1029, 2016.[21] J. Gou and M. J. Ward. Oscillatory dynamics for a coupled membrane-bulk diffusion model with Fitzhugh-Nagumokinetics.
SIAM J. Appl. Math. , 76(2):776–804, 2016.[22] T. Gregor, K. Fujimoto, N. Masaki, and S. Sawai. The onset of collective behavior in social amoebae.
Science ,328(5981):1021–1025, 2010.[23] S. Güttel and F. Tisseur. The nonlinear eigenvalue problem.
Acta Numerica , 26(1):1–94, 2017.[24] B. A. Hense, C. Kuttler, J. Müller, M. Rothballer, A. Hartmann, and J. U. Kreft. Does efficiency sensing unify diffusionand quorum sensing?
Nature Reviews. Microbiology , 5:230–239, 2007.[25] M. A. Henson, Müller D., and M. Reuss. Cell population modelling of yeast glycolytic oscillations.
Biochem J. , 368:433–446, 2002.[26] S. Iyaniwura and M. J. Ward. Localized signaling compartments in 2-d coupled by a bulk diffusion field: Quorumsensing and synchronous oscillations in the well-mixed limit. to appear, Europ. J. Appl. Math. , 2020.[27] K. Kamino, K. Fujimoto, and Sawai S. Collective oscillations in developing cells: Insights from simple systems.
Develop.Growth Differ. , 53:503–517, 2011.[28] T. Kolokolnikov, M. S Titcombe, and M. J. Ward. Optimizing the fundamental Neumann eigenvalue for the Laplacianin a domain with small traps.
Europ. J. Appl. Math. , 16(2):161–200, 2005.[29] B. W. Li, X. Z. Cao, and C. Fu. Quorum sensing in populations of spatially extended chaotic oscillators coupledindirectly via a heterogeneous environment.
Journal of NonLinear Science , 27(6):1667–1686, 2017.[30] B. W. Li, C. Fu, H. Zhang, and X. Wang. Synchronization and quorum sensing in an ensemble of indirectly coupledchaotic oscillators.
Physical Review E , 86(4):046207, 2012.[31] C. K. Macnamara and M. A. J. Chaplain. Spatio-temporal models of synthetic genetic oscillations.
Math. Biosc. Eng. ,14:249–262, 2017.[32] P. Melke, P. Sahlin, A. Levchenko, and H. Jonsson. A cell-based model for quorum sensing in heterogeneous bacterialcolonies.
PLoS Computational Biology , 6(6):e1000819, 2010.[33] P. Mina, M. di Benardo, N. J. Savery, and K. Tsaneva-Atanasova. Modeling emergence of oscillations in communicatingbacteria: a structured approach from one to many cells.
J. Royal Society Interface , 10:20120612, 2012.[34] J. Müller, C. Kuttler, B. A. Hense, M. Rothballer, and A. Hartmann. Cell–cell communication by quorum sensing anddimension-reduction.
Journal of mathematical biology , 53(4):672–702, 2006.[35] J. Müller and H. Uecker. Approximating the dynamics of communicating cells in a diffusive medium byodes—homogenization with localization.
Journal of mathematical biology , 67(5):1023–1065, 2013.[36] H. Nakao. Phase reduction approach to synchronisation of nonlinear oscillators.
Contemporary Physics , 57(2):188–214,2015. 3937] V. Nanjundiah. Cyclic AMP oscillations in Dictyostelium Discoideum: Models and observations.
Biophysical chemistry ,72(1-2):1–8, 1998.[38] F. Naqib, T. Quail, L. Musa, H. Vulpe, J. Nadeau, J. Lei, and L. Glass. Tunable oscillations and chaotic dynamics insystems with localized synthesis.
Phys. Rev. E , 85:046210, Apr 2012.[39] J. Noorbakhsh, D. J. Schwab, A. E. Sgro, T. Gregor, and P. Mehta. Modeling oscillations and spiral waves in dictyostelium populations.
Phys Rev. E. , 91(6):062711, 2015.[40] F. Paquin-Lefebvre, W. Nagata, and M. J Ward. Pattern formation and oscillatory dynamics in a two-dimensionalcoupled bulk-surface reaction-diffusion system.
SIAM J. Appl. Math. , 80(3):1520–1545, 2020.[41] B. Pietras and A. Daffertshofer. Network dynamics of coupled oscillators and phase reduction techniques.
PhysicsReports , 819:1–105, 2019.[42] 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.[43] J. Schütze and J. Wolf. Spatio-temporal dynamics of glycolysis in cell layers: A mathematical model.
Biosystems ,99(2):104–108, 2010.[44] E. E. Sel’kov. Self-oscillations in glycolysis 1. a simple kinetic model.
European Journal of Biochemistry , 4(1):79–86,1968.[45] M. E. Taga and B. L. Bassler. Chemical communication among bacteria.
Proceedings of the National Academy ofSciences , 100(suppl 2):14549–14554, 2003.[46] A. F. Taylor, M. Tinsley, and K. Showalter. Insights into collective cell behavior from populations of coupled chemicaloscillators.
Phys. Chemistry Chem Phys. , 17(31):20047–20055, 2015.[47] A. F. Taylor, M. Tinsley, F. Wang, Z. Huang, and K. Showalter. Dynamical quorum sensing and synchronization inlarge populations of chemical oscillators.
Science , 323(5914):614–6017, 2009.[48] M. R. Tinsley, A. F. Taylor, Z. Huang, and K. Showalter. Emergence of collective behavior in groups of excitablecatalyst-loaded particles: Spatiotemporal dynamical quorum sensing.
Phys. Rev. Lett. , 102:158301, 2009.[49] M. R. Tinsley, A. F. Taylor, Z. Huang, F. Wang, and K. Showalter. Dynamical quorum sensing and synchronization incollections of excitable and oscillatory catalytic particles.
Physica D , 239(11):785–790, 2010.[50] A. Trovato, F. Seno, M. Zanardo, S. Alberghini, A. Tondello, and A. Squartini. Quorum vs. diffusion sensing: Aquantitative analysis of the relevance of absorbing or reflecting boundaries.
FEMS Microbiology Letters , 352:198–203,2014.[51] H. Uecker, J. Müller, and B. A. Hense. Individual-based model for quorum sensing with background flow.
Bulletin ofmathematical biology , 76(7):1727–1746, 2014.[52] J. P. Ward, J. R. King, A. J. Koerber, P. Williams, J. M. Croft, and R. E. Sockett. Mathematical modelling of quorumsensing in bacteria.
Mathematical Medicine and Biology , 18(3):263–292, 2001.[53] M. J Ward. Spots, traps, and patches: Asymptotic analysis of localized solutions to some linear and nonlinear diffusivesystems.
Nonlinearity , 31(8):R189, 2018.[54] M. J. Ward, W. D. Henshaw, and J. B. Keller. Summing logarithmic expansions for singularly perturbed eigenvalueproblems.
SIAM J. Appl. Math. , 53(3):799–828, 1993.[55] M. Whiteley, S. P. Diggle, and E. P. Greenberg. Bacterial quorum sensing: the progress and promise of an emergingreserch areas.
Nature , 551:7680, November 15 2017.[56] J. Wolf and R. Heinrich. Effect of cellular interaction on glycolytic oscillations in yeast: a theoretical investigation.
Biochem J. , 345:321–334, 2000.[57] B. Xu and P. Bressloff. A PDE-DDE model for cell polarization in fission yeast.