Abstract
In many micro- and macro-scale systems, collective dynamics occurs from the coupling of small spatially segregated, but dynamically active, units through a bulk diffusion field. This bulk diffusion field, which is both produced and sensed by the active units, can trigger and then synchronize oscillatory dynamics associated with the individual units. In this context, we analyze diffusion-induced synchrony for a class of cell-bulk ODE–PDE system in \({\mathbb {R}}^2\) that has two spatially segregated dynamically active circular cells of small radius. By using strong localized perturbation theory in the limit of small cell radius, we calculate the steady-state solution and formulate the linear stability problem. For Sel’kov intracellular reaction kinetics, we analyze how the effect of bulk diffusion can trigger, via a Hopf bifurcation, either in-phase or anti-phase intracellular oscillations that would otherwise not occur for cells that are uncoupled from the bulk medium. Phase diagrams in parameter space where these oscillations occur are presented, and the theoretical results from the linear stability theory are validated from full numerical simulations of the ODE–PDE system. In addition, the two-cell case is extended to study the onset of synchronous oscillatory instabilities associated with an infinite hexagonal arrangement of small identical cells in \({\mathbb {R}}^2\) with Sel’kov intracellular kinetics. This analysis for the hexagonal cell pattern relies on determining a new, computationally efficient, explicit formula for the regular part of a certain periodic reduced-wave Green’s function.
Similar content being viewed by others
References
Stankovski T, Pereira T, McClintock PVE, Stefanovska A (2017) Coupling functions: universal insights into dynamical interaction mechanisms. Rev Mod Phys 89:045001
Müller J, Kuttler C, Hense BA, Rothballer M, Hartmann A (2006) Cell–cell communication by quorum sensing and dimension-reduction. J Math Biol 53(4):672–702
Gou J, Ward MJ (2016) Oscillatory dynamics for a coupled membrane-bulk diffusion model with Fitzhugh–Nagumo kinetics. SIAM J Appl Math 76(2):776–804
Iyaniwura S, Ward MJ (2020) Synchrony and oscillatory dynamics for a 2-D PDE-ODE model of diffusion-sensing with small signaling compartments. submitted. SIAM J Appl Dynam Syst 4
Danø S, Madsen MF, Sørensen PG (2007) Quantitative characterization of cell synchronization in yeast. Proc Natl Acad Sci 104(31):12732–12736
Danø S, Sørensen PG, Hynne F (1999) Sustained oscillations in living cells. Nature 402(6759):320–322
De Monte S, d’Ovidio F, Danø S, Sørensen PG (2007) Dynamical quorum sensing: population density encoded in cellular dynamics. Proc Natl Acad Sci 104(47):18377–18381
Gao M, Zheng H, Ren Y, Lou R, Wu F, Yu W, Liu X, Ma X (2016) A crucial role for spatial distribution in bacterial quorum sensing. Sci. Rep. 6: 34695
Gregor T, Fujimoto K, Masaki N, Sawai S (2010) The onset of collective behavior in social amoebae. Science 328(5981):1021–1025
Marenda M, Zanardo M, Trovato A, Seno F, Squartini A (2016) Modeling quorum sensing trade-offs between bacterial cell density and system extension from open boundaries. Sci Rep 6: 39142
Trovato A, Seno F, Zanardo M, Alberghini S, Tondello A, Squartini A (2014) Quorum vs. diffusion sensing: A quantitative analysis of the relevance of absorbing or reflecting boundaries. FEMS Microbiol Lett 352(2):198–203
Taylor AF, Tinsley M, Showalter K (2015) Insights into collective cell behavior from populations of coupled chemical oscillators. Phys Chem Chem Phys 17(31):20047–20055
Tinsley MR, Taylor AF, Huang Z, Showalter K (2009) Emergence of collective behavior in groups of excitable catalyst-loaded particles: spatiotemporal dynamical quorum sensing. Phys Rev Lett 102:158301
Tinsley MR, Taylor AF, Huang Z, Wang F, Showalter K (2010) Dynamical quorum sensing and synchronization in collections of excitable and oscillatory catalytic particles. Physica D 239(11):785–790
Yamamoto SY, Surko CM, Maple MB (1995) Spatial coupling in heterogeneous catalysis. J Chem Phys 103(18):8209
FlexPDE. PDE Solutions inc, 2015
Linton CM (2010) Lattice sums for the Helmholtz equation. SIAM Rev 52(4):630–674
Moroz A (2006) Quasi-periodic Green’s functions of the Helmholtz and Laplace equations. J Phys A: Math Gen 39(36):11247–11282
Beylkin G, Kurcz C, Monzón L (2008) Fast algorithms for Helmholtz Green’s functions. Proc R Soc A 464:3301–3326
Iron D, Rumsey J, Ward MJ (2015) A hybrid asymptotic-numerical method for stability thresholds of periodic patterns of localized spots for reaction-diffusion systems in \({\mathbb{R}}^2\). Eur J Appl Math 26(3):325–353
Gou J, Ward MJ (2016) An asymptotic analysis of a 2-D model of dynamically active compartments coupled by bulk diffusion. J Nonlinear Sci 26(4):979–1029
Ward MJ (2018) Spots, traps, and patches: asymptotic analysis of localized solutions to some linear and nonlinear diffusive systems. Nonlinearity 31(8):R189
Ward MJ, Henshaw WD, Keller JB (1993) Summing logarithmic expansions for singularly perturbed eigenvalue problems. SIAM J Appl Math 53(3):799–828
Abramotitz M, Stegun I (1965) Handbook of mathematical functions, 9th edn. Dover Publications, New York
Betcke T, Highan NG, Mehrmann V, Schröder C, NLEVP Tisseur F (2013) A collection of nonlinear eigenvalue problems. ACM Trans Math Softw 39(2):7.1–7.28
Güttel S, Tisseur F (2017) The nonlinear eigenvalue problem. Acta Numer 26(1):1–94
Betcke T, Highan NG, Mehrmann V, Porzio GMN, Schröder C, NLEVP Tisseur F (2019) A collection of nonlinear eigenvalue problems. users’ guide. MIMS EPring 2011.117, Manchester Institue for Mathematical Sciences, The University of Manchester, UK, p 10, updated
Sel’kov EE (1968) Self-oscillations in glycolysis 1. A simple kinetic model. Eur J Biochem 4(1):79–86
Continuation Test: a MATLAB library which defines test functions for continuation codes. Accessed: 2020-02-26
Alciatore D, Miranda R (1995) A winding number and point-in-polygon algorithm. Glaxo Virtual Anatomy Project Research Report, Department of Mechanical Engineering, Colorado State University
Iron D, Rumsey J, Ward MJ, Wei JC (2014) Logarithmic expansions and the stability of periodic patterns of localized spots for reaction-diffusion systems in \({\mathbb{R}}^2\). J Nonlinear Sci 24(5):564–627
Kuchment P (1993) Floquet theory for partial differential equations. Birkhauser, Basel
Trefethon N (2018) Series solution of Laplace problems. ANZIAM J 60:1–26
Uecker H, Müller J, Hense BA (2014) Individual-based model for quorum sensing with background flow. Bull Math Biol 76:1727–1746
Piessens R (2018) The Hankel transform. In Poularikas AD (ed.), Transforms and applications handbook chapter 9. 3rd ed. CRC Press, Boca Raton
Chen X, Oshita Y (2007) An application of the modular function in nonlocal variational problems. Arch Rat Mech Anal 186(1):109–132
Acknowledgements
Michael Ward gratefully acknowledges the financial support from the NSERC Discovery Grant Program.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendices
The quasi-periodic reduced-wave Green’s function
We develop a rapidly converging expansion for the Bloch Green’s function \(G_\mathrm{{b}}(\mathbf{x})\) and its regular part \(R_\mathrm{{b}}({\pmb k};z)\) as defined by (4.6). To do so, we extend the analysis in §5.1 of [20] for the case \(z=1\), as was motivated by the methodology in [19], to allow for complex-valued z with \(\text{ Re }(z)>0\). Then, by setting \(z=1+\tau \lambda \) we obtain \(R_\mathrm{{b}}({\pmb k};1+\tau \lambda )\), as needed in (4.5). Moreover, we identify that \(R_p\equiv R_\mathrm{{b}}(\pmb {0};1)\) in (4.3) and (4.4).
We begin by representing the solution to (4.6) as the sum of free-space Green’s functions
which ensures that the quasi-periodicity condition in (4.6) is satisfied. Then, to analyze (A.1) we use the Poisson summation formula which converts a sum over the hexagonal lattice \(\Lambda \) to a sum over the reciprocal lattice \(\Lambda ^{\star }\) of (4.2). In the notation of [19], we have (see Proposition 2.1 of [19])
Here \(|\Omega |\) is the area of the primitive cell of the lattice, while \({\hat{f}}\) is the Fourier transform of f, defined in \({\mathbb {R}}^2\) by
Upon applying (A.2) to (A.1) we obtain that the sum over the reciprocal lattice consists of free-space Green’s functions in the Fourier domain. By taking the Fourier transform of the PDE \(\Delta G_{\mathrm{free}} - z D^{-1} G_{\mathrm{free}}=-\delta (\mathbf{x})\) for the free-space Green’s function, we obtain that \({\hat{G}}_{\mathrm{free}} (\pmb p) = {\hat{G}}_{\mathrm{free}} (|\pmb p|)\), where
Then, from (A.2) and (A.1), and by using \(|\Omega |=1\), we obtain that
In order to obtain a rapidly converging infinite series representation for \(G_\mathrm{{b}}(\mathbf{x})\), we introduce the decomposition
where the function \(\alpha (2\pi \pmb d-\pmb k,\eta )\) is defined by
Here \(\eta >0\) is a real-valued cut-off parameter, which is specified below. Since \(\text{ Re }(z)>0\), we readily calculate that
With this choice for \(\alpha \), the sum over \(\pmb d\in \Lambda ^{\star }\) in (A.5) of the first set of terms in (A.6), as labeled by
converges absolutely when \(\text{ Re }(z)>0\).
Next, we calculate the lattice sum in (A.5) over the second set of terms in (A.6) by using the inverse transform (A.3) after first writing \((1-\alpha )\,{\hat{G}}_{\mathrm{free}}\) as an integral. To do so, we define \(\rho \equiv |2\pi \pmb d-\pmb k|\), so that from (A.7) and (A.4)
In recognizing the middle term in (A.9) as the integral \({\hat{\chi }}(\rho )\) we used the easily verified definite integral
Next, we take the inverse Fourier transform of (A.9). To do so, we use two key facts. Firstly, the inverse Fourier transform of a radially symmetric function \({\hat{f}}(\rho )\) is the inverse Hankel transform of order zero (cf. [35]), so that \(f(r)=(2\pi )^{-1}\int _0^\infty {\hat{f}}(\rho )\,J_0(\rho r)\,\rho \, \mathrm{{d}}\rho \). Secondly, from [35], we recall the well-known inverse Hankel transform
In this way, we calculate using the definition of \({\hat{\chi }}(\rho )\) in (A.9) that
In the notation of [19], we then define \(F_{\mathrm{sing}}(|\mathbf{x}|)\) as
Therefore, by using the Poisson summation formula (A.2) to calculate lattice sum in (A.5) over the second set of terms in (A.6), we obtain
In summary, the Bloch Green’s function for the reduced-wave operator in the spatial domain, satisfying (4.6), is \(G_\mathrm{{b}}(\mathbf{x})=G_{\mathrm{fourier}}(\mathbf{x}) + G_{\mathrm{spatial}}(\mathbf{x})\), representing the sum of (A.8) and (A.12). In this way, we have
where \(F_{\mathrm{sing}}(|\mathbf{x}|)\) is defined in (A.11). In terms of the Ewald cut-off parameter \(\eta \), we observe from (A.8) that \(G_{\mathrm{fourier}}(\mathbf{x}) \rightarrow 0\) as \(\eta \rightarrow 0\), while from (A.12) and (A.11) we get that \(G_{\mathrm{spatial}}(\mathbf{x}) \rightarrow 0\) as \(\eta \rightarrow \infty \).
The next step in the analysis is to isolate the logarithmic singularity in (A.13) as \(\mathbf{x}\rightarrow \mathbf{0}\), so as to identify the regular part \(R_\mathrm{{b}}({\pmb k};z)\) defined by
From (A.8) we observe that \(G_{\mathrm{fourier}}(\mathbf{0})\) is finite for all \(|2\pi \pmb d -\pmb k|\) and \(0<\eta <\infty \) when \(\text{ Re }(z)>0\). Likewise, for the last set of terms in (A.13), we can take the limit \(\mathbf{x}\rightarrow \mathbf{0}\) to conclude that \(\sum _{\genfrac{}{}{0.0pt}{}{\pmb l\in \Lambda }{\pmb l\ne \mathbf{0}}}\mathrm{{e}}^{\mathrm{{i}}\pmb k\cdot \pmb l}\, F_{\mathrm{sing}}(|\pmb l|)\) is also finite. As a result, the logarithmic singularity in \(G_\mathrm{{b}}\) as \(\mathbf{x}\rightarrow 0\) must arise from the middle term, \(F_{\mathrm{sing}}(r)\), in (A.13), where \(r\equiv |\mathbf{x}|\).
To analyze \(F_{\mathrm{sing}}(r)\) when \(\text{ Re }(z)>0\), we write \(z=|z|\mathrm{{e}}^{\mathrm{{i}}\theta }\), where \(\theta \equiv \text{ Arg }z\) satisfies \(|\theta |<{\pi /2}\). In the integrand of (A.11) we introduce the new variable t by
so that (A.11) becomes
Here \(\theta \equiv \text{ Arg }\,z\) and we have specified the principal branch for \(\sqrt{z}\). We then split the integration range in (A.16) as
To calculate the first term in (A.17), we use a contour integration over the box-shaped contour \(-b\le \text{ Re }(t)\le b\) and \(0\le \text{ Im }(t)<{\pi /2}\) in the complex t-plane. Since there are no residues within the contour, and we have exponential decay on the sides of the box as \(b\rightarrow \infty \), we can replace \(t-{\mathrm{{i}}\theta /2}\) by t in (A.17) and then use symmetry to obtain
where the last equality in (A.18) follows from the substitution \(\xi =\cosh (t)\). The last integral in (A.18) is identified by using the well-known integral representation \( K_{0}(\omega ) = \int _{1}^{\infty } {\mathrm{{e}}^{-\omega \xi }/\sqrt{\xi ^2-1}} \, \mathrm{{d}}\xi \) for the modified Bessel function of the second kind of order zero, which is valid for \(|\text{ Arg }\,\omega |<{\pi /2}\). In this way, we obtain from (A.18) and (A.17) that
Next, we provide a more tractable representation for J(r; z). We introduce a new variable \(\xi \) defined by
so that \(\xi =1\) when \(t=\beta \) and \(\xi \rightarrow +\infty \) as \(t\rightarrow -\infty \). In terms of \(\xi \), the exponential in the integral J is
By using (A.21) and \(\mathrm{{d}}t=-{\mathrm{{d}}\xi /\xi }\) within the integral J in (A.19), we obtain that (A.19) transforms to
Finally, we use (A.22) to extract the local behavior of \(F_{\mathrm{sing}}(r)\) as \(r\rightarrow 0\). As \(r\rightarrow 0\), we calculate for J(r; z) that
where \(E_1(w)\equiv \int _{1}^{\infty } \xi ^{-1} \mathrm{{e}}^{-w\xi }\, \mathrm{{d}}\xi \), for \(\text{ Re }(w)>0\), is the well-known exponential integral. In addition, in (A.22) we use \(K_0(w) = -\log {w}+\log {2}-\gamma _e+o(1)\) as \(|w|\rightarrow 0\), where \(\gamma _e\) is Euler’s constant. In this way, (A.22) and (A.23) yield
Finally, by substituting (A.13) in (A.14), and using (A.24), we obtain the result for \(R_\mathrm{{b}}({\pmb k};z)\) given in (4.7).
Formulation on the Wigner–Seitz cell
In this appendix, we provide a key result for \(R_\mathrm{{b}}({\pmb k};z)\). However, before doing so, we first must obtain a more refined description of the fundamental Wigner–Seitz (FWS) cell, as was discussed in §2.2 of [31]. For a general lattice, there are eight nearest neighbor lattice points to \(\mathbf{x}=0\) given by the set
For each (vector) point \({\pmb P}_i\in P\), for \(i=1,\ldots ,8\), the Bragg line \(L_i\) is defined as the line that crosses the point \({{\pmb P}_i/2}\) orthogonally to \({\pmb P}_i\). The unit outer normal to \(L_i\) is labeled by \({\pmb \eta }_i\equiv {{\pmb P}_i/|{\pmb P}_i|}\). The convex hull generated by these Bragg lines is the FWS cell \(\Omega \). Specifically, for the hexagonal lattice (4.1) its boundary \(\partial \Omega \) is the union of exactly six Bragg lines, and the centers of the Bragg lines generating \(\partial \Omega \) are re-indexed as \({{\pmb P}_i/2}\) for \(i=1,\ldots ,6\). The boundary \(\partial \Omega \) of \(\Omega \) is the union of the re-indexed Bragg lines \(L_i\) for \(i=1,\ldots ,6\), and is parameterized segment-wise as
where \(2t_i\) is the length of \(L_i\). Here \({\pmb \eta }_{i}^{\perp }\) is the direction perpendicular to \({\pmb P}_i\), and is, therefore, tangent to \(L_i\). From this construction, Bragg lines on \(\partial \Omega \) must come in pairs. In particular, if \({\pmb P}\) is a neighbor of \(0\) and the Bragg line crossing \({{\pmb P}/2}\) lies on \(\partial \Omega \), it follows by symmetry that the Bragg line crossing \({-{\pmb P}/2}\) must also lie on \(\partial \Omega \).
Next, we reformulate the PDE (4.6) for the reduced-wave Bloch Green’s function \(G_\mathrm{{b}}\) on \({\mathbb {R}}^2\) to an equivalent PDE on the FWS \(\Omega \). This is done by imposing a boundary operator \(\mathcal{{P}}_k\) on \(\partial \Omega \) that incorporates the quasi-periodic condition in (4.6). This equivalent PDE is
where \({{\pmb k}/(2\pi )}\in \Omega _\mathrm{{b}}\). In (B.3), the boundary operator is defined by
Here \(L_{i}\) and \(L_{-i}\) denote two parallel Bragg lines on opposite sides of \(\partial \Omega \) for \(i=1,\ldots ,3\), while \(\mathbf{x}_{i1}\in L_i\) and \(\mathbf{x}_{i2}\in L_{-i}\) are any two opposing points on these Bragg lines. In this notation, periodic boundary conditions refer to \(G_\mathrm{{b}}\in \mathcal{{P}}_{0}\).
With this reformulation of (4.6) to the PDE (B.3) on the FWS cell \(\Omega \), with unit area \(|\Omega |=1\), we now derive a result is needed in Sect. 4.2. Firstly, we obtain for the periodic problem, with \(G_\mathrm{{b}}\in \mathcal{{P}}_{0}\) on \(\partial \Omega \), that \(G_\mathrm{{b}}={{\mathcal {O}}}(D)\) when \(D\gg 1\). By expanding in powers of D, we readily derive for \({\pmb k}={\pmb 0}\) that
where \(G_{p0}(\mathbf{x})\) is the periodic source-neutral Green’s function satisfying
and normalized by \(\int _{\Omega } G_{p0} \, \mathrm{{d}}\mathbf{x}= 0\). From Theorem 1 of [36], \(R_{p0}\approx -0.21027\) is given explicitly by
By taking the regular part of (B.5) we obtain the two-term result for \(R_\mathrm{{b}}({\pmb 0};1+\tau \lambda )\) in (4.12), which is valid for \(D\gg 1\).
Rights and permissions
About this article
Cite this article
Iyaniwura, S.A., Gou, J. & Ward, M.J. Synchronous oscillations for a coupled cell-bulk ODE–PDE model with localized cells on \({\mathbb {R}}^2\). J Eng Math 127, 18 (2021). https://doi.org/10.1007/s10665-021-10113-7
Received:
Accepted:
Published:
DOI: https://doi.org/10.1007/s10665-021-10113-7