Synchronization, Oscillator Death, and Frequency Modulation in a Class of Biologically Inspired Coupled Oscillators

The general purpose of this paper is to build up on our understanding of the basic mathematical principles that underlie the emergence of synchronous biological rhythms, in particular, the circadian clock. To do so, we study the role that the coupling strength, coupling type, and noise play in the synchronization of a system of coupled, non-linear oscillators. First, we study a deterministic model based on Van der Pol coupled oscillators, modeling a population of diffusively coupled cells, to find regions in the parameter space for which synchronous oscillations emerge and to provide conditions under which diffusive coupling kills the synchronous oscillation. Second, we study how noise and coupling interact and lead to synchronous oscillations in linearly coupled oscillators, modeling the interaction between various pacemaker populations, each having an endogenous circadian clock. To do so, we use the Fokker-Planck equation associated to the system. We show how coupling can tune the frequency of the emergent synchronous oscillation, which provides a general mechanism to make fast (ultradian) pacemakers slow (circadian) and synchronous via coupling. The basic mechanisms behind the generation of oscillations and the emergence of synchrony that we describe here can be used to guide further studies of coupled oscillations in biophysical non-linear models.


INTRODUCTION
The study of circadian rhythms has been a subject of great interest for a long time. The majority of the first studies were mainly based on observations in plants [1][2][3][4]. The study of circadian rhythms from a mathematical perspective reached a milestone with the work of Kalmus and Wigglesworth, a biologist and a mathematician, respectively, who associated of circadian rhythms to the existence of a limit cycle, using a hydraulic system as analogy. Kalmus and Wigglesworth presented their work entitled "Shock excited system as models for biological rhythms" along with several mathematical models of circadian rhythms at a Symposium on Biological Rhythms carried out at Cold Spring Harbor in the United States of America in 1960 [5][6][7]. Lots of other important works on circadian rhythms were presented in this symposium, but the work of Kalmus and Wigglesworth was key in establishing a better mathematical formalism for the study of circadian rhythms. Although many researchers followed the theoretical path proposed by Kalmus and Wigglesworth, the mathematical study of circadian rhythms was finally established by Arthur Winfree (biologist and mathematician), who introduced topology for the description of several aspects of circadian rhythms. An excellent summary of many of the earlier works can be found in Arthur Winfree's master book entitled "The Geometry of Biological Time" [8]. The number of studies about biological rhythms at large has increased greatly in the last two decades, in part due to new technological advances. Particularly related to circadian rhythms, there is now evidence that there is rhythmic patterns of activity at the molecular [9,10], cellular [11,12], tissue [13,14], and systems levels [15][16][17], and that circadian regulation is involved in jointly regulating activity in all those different levels of biological organization [18][19][20], and also, taking into account interactions with the environment [21] and perturbations induced by behavior [22,23].
Mathematical modeling and experimental characterizations of different properties of circadian rhythms have been combined to produce explanations and hypotheses about the rhythmicity in biological phenomena [24][25][26][27][28]. Of particular interest, the ontogeny of circadian rhythm in the crayfish has been studied by Lara-Aparicio et al. [29] combining theoretical and experimental perspectives, building phenomenological mathematical models that capture a series of experimental results involving the synchronization of electro-retinogram activity in crayfish [30][31][32]. These models couple two van der Pol oscillators [33] represented by state vectors, and include parameters representing the frequency of the oscillators, the radii of the limit cycles, and the first coordinate of the center of the limit cycle. The system is setup such that one oscillator is driving the behavior of the other oscillator, but not vice-versa. One of the main findings with this model is that the driving oscillator induces an Andronov-Hopf bifurcation in the driven oscillator and regulates its frequency.
The model by Aparicio et al. simulates, explains, and has suggested new biological experiments, it is simple enough to allow analytical approaches, and it has provided useful insights about questions referring to the ontogeny of the circadian rhythm in crayfish from the childhood to adult stages. For instance, a hypothesis about the existence of a hormone, which was experimentally detected, was generated from the model. The model also allowed Lara-Aparicio et al. to study synchronization of circadian rhythms with external signals like day and night light cycle [34]. By studying basic principles underlying the generation of oscillations in coupled non-linear systems, these researchers were able to conjecture that circadian rhythms can result from coupling systems of cells, each one oscillating with an ultradian oscillation [35][36][37]. Synchronization among cells emerges naturally as a motivating theme that has been studied through systems of non-linear Equations [38] representing n oscillators with the classical van der Pol non-linear damping for the terms responsible for the oscillatory dynamics.
In the present paper, we extend the work in Lara-Aparicio et al. [29] and Barriga-Montoya et al. [38] by analyzing two qualitative mathematical models, each capturing a different level of organization in the ontogeny of circadian rhythms. Inspired by gap junction coupling between neurons, or similarly, by chemical coupling in self oscillating networks, we study the bifurcation structure in a deterministic model assuming that the coupling between the oscillators is diffusive. The resulting dynamics resemble neuronal activity at the cellular level. Then, using graph theoretical methods and center manifold theory, we show that synchronous oscillations appear via a Hopf bifurcation in a population of pacemaker oscillators. In this case, the bifurcation parameter is thought of as a representative of the developmental stage of the neurons. We further explore the phenomenon of oscillator death: although the single neurons are endogenously oscillating for sufficiently large values of the bifurcation parameter, the population oscillation dies for sufficiently large coupling, which suggests that the weak coupling hypothesis must be satisfied for robust synchronous oscillations to occur. In the case of all-to-all coupling, we provide necessary conditions for oscillator death to occur and leave the derivation of sufficient conditions to a future report.
Frequency modulation is also an important phenomenon that can be studied with these models. In consideration of the results from the first analysis, we shift our focus to the frequency modulations that emerge as a result of the interconnection of various circadian pacemaker populations. In doing so, we estimate the synchronization frequency as a function of the intrinsic frequencies of the oscillators, their coupling strength, and the topology of the network. To do so, we construct a second model that can be thought of as a stochastic, larger-scale extension of the first model we present. In this case, each population is assumed to be an endogenous oscillator and the coupling is assumed to be linear. Using linear stochastic analysis and under the assumption that the population oscillations are synchronized, we derive a scalar Fokker-Plank equation. The model captures an important feature of circadian rhythm ontogeny: the emergence of low frequency (circadian) oscillations from coupled high frequency (ultradian) oscillators [35][36][37][38]. Future work will aim at deriving conditions on the intrinsic frequencies, the coupling strength, and the network topology, that ensure synchronization of the population oscillations.
It is worth noticing that, although the motivation for the present paper is centered around circadian rhythms, the model captures more general phenomena. Our findings include that in diffusively coupled cells resting node dynamics imply global asymptotic stability, oscillating node dynamics imply global synchronization for small coupling, and multistability between oscillator death and global synchronization for large coupling. In stochastic, linearly coupled populations, we describe the dynamical mechanisms through which coupling modulates the frequency of the synchronous oscillation. To the author's knowledge, both phenomena are new from a non-linear collective phenomena perspective. Among other reasons, it is surprising that passive coupling like diffusive coupling could kill an oscillation and create multistability. Similarly, we are not aware of any work providing mechanistic explanations on how coupling can tune a global oscillation frequency.
The paper is organized as follows. In section 2, we present and analyze the model of diffusively coupled oscillators. Two theorems are proved about global asymptotic stability for resting cells and global synchronization for oscillating cells and weak coupling. We then derive sufficient conditions for diffusive coupling-induced oscillator death. In section 3 we present and analyze the model of linearly coupled oscillators. In particular, we derive an explicit formula for the emergent synchronization frequency as a function of the coupling topology and oscillator natural frequencies. Finally we discuss the presented results in section 4.

GLOBAL SYNCHRONIZATION AND OSCILLATOR DEATH IN DIFFUSIVELY COUPLED OSCILLATORS
We regard a network of N coupled oscillators as a directed graph G with N vertices, with a network topology codified by a matrix , where a ij ≥ 0 represent connection weights. If oscillator i receives signals from oscillator j, then the graphical representation of G has an arrow from j to i, and a ij > 0. If a ij > 0, the signal received by oscillator i from oscillator j is a ij (x j − x i ). Assume that the dynamics for each oscillator satisfies the following coupled oscillator dynamicṡ where µ is the global coupling strength, which uniformly scales the coupling weights a ij .
In the absence of coupling the oscillator dynamics reduces to the modified van der Pol equation These equations define a simple dynamical system that can transition between global asymptotic stability and almost global convergence to a hyperbolic limit cycle through variations of the control parameter λ ∈ R, providing a simple model for various biological systems that exhibit the same transition, in particular neurons and molecular oscillators. For λ < 0 the non- is always positive, which leads to damped oscillations. For λ > 0 the dissipation coefficient becomes negative close to the origin, which leads to sustained oscillations through a Hopf bifurcation. A generic trajectory belonging to the family of periodic orbits born at the Hopf bifurcation has the form √ λ(cos(ω i t + θ 0 ), sin(ω i t + θ 0 )), where the θ 0 is the initial phase.
In the presence of coupling, Equations (1) represent a generic network of diffusively coupled non-linear oscillators. As mentioned earlier, the diffusive form of the coupling can be thought of as gap junction coupling in a neuronal population, or diffusive coupling between non-linear molecular oscillators. In both interpretations, the graph topology is necessarily undirected, that is, a ij = a ji . However, the mathematical results presented in this section hold under more general assumptions and we present them in the general form.

Global Synchronization
We start by recalling some basic graph-theoretical definitions and facts. The graph G is said to be strongly connected if for each pair of nodes in G, there exists a directed path between them. G is balanced if N j=1 a ij = N j=1 a ji for all i. The Laplacian matrix L = L ij : i, j = 1, ..., n for the graph G is such that L ij = −a ij if i = j, and L ii = N j=1 a ij . Note that the vector of ones is always a right null eigenvector of L and zero is always an eigenvalue of L (L1 N = 0). It can be shown that a graph is strongly connected if and only if zero is a simple eigenvalue of the Laplacian matrix [39]. Obviously, symmetric graphs (i.e., satisfying a ij = a ji ) are balanced, but the converse is not true. Consider, for instance, a directed ring.
The global behavior of the system (1) for λ < 0 and µ ≥ 0, for a network with connectivity represented by a generic balanced graph is characterized by the next theorem (Figure 1). Theorem 2.1. Assume that the graph G is balanced and that ω i = ω j = ω for all i, j = 1, . . . , N. If λ < 0, then the origin is globally asymptotically stable and locally exponentially stable for all µ ≥ 0.
Proof. We consider the quadratic Lyapunov function The derivative of V along the trajectories of the system (1) giveṡ which after substitution of the derivativesẋ i andẏ i , can be bounded from above as follows To show that the last term is negative definite, we use the balanced interconnection hypothesis. Let a ij x i .
As a consequence, Then the global part of statement follows by LaSalle invariance principle [40]. The local part follows by observing that for λ < 0 the linearization of model (1) at the origin is non-singular and therefore asymptotic stability of the origin implies that all eigenvalues have negative real part. Remark. Because exponentially stability implies robustness to small perturbations, Theorem 2.1 remains true for small heterogeneity in the natural frequencies.
Note that, for λ < 0, the system in Equation (1) exhibits exponentially damped oscillations toward the origin (Figure 1, top panels).
Next we show that, at λ = 0 and identical natural frequencies model (1) undergoes a supercritical Hopf bifurcation inside the consensus space provided the graph is strongly connected. The linearization of the system (1) is given by where I N is the N-dimensional identity matrix and L is the network Laplacian defined in section 2.1. Let 1 N be the Ndimensional vector of ones. Given a (complex) vector ν = (v, w) in the consensus space C, that is, v = a1 N and w = b1 N for some a, b ∈ C, the eigenvalue problem for the Jacobian matrix Equation (4), restricted to the consensus space, takes the form where ξ is a (complex) eigenvalue. In the third equality we used the fact that 1 N is a right null eigenvector of L. The last equality shows that ν ∈ C is an eigenvector of J with eigenvalue ξ if and only if its components a and b satisfy We can now easily solve Equation (5) to obtain the eigenvalues/eigenvectors pairs For λ = 0 we got two purely imaginary eigenvalues which correspond to a supercritical Hopf bifurcation of model (1) inside the consensus space, as summarized in the following theorem and illustrated in (Figure 1, bottom panels).

Theorem 2.2.
For almost all balanced, strongly connected interconnection topologies the following holds. For all µ > 0, the system (1) undergoes a supercritical Hopf bifurcation at λ = 0 with center manifold given by the consensus space C. Moreover, the family of periodic solutions born at the Hopf bifurcation are exponentially asymptotically stable and correspond to synchronous oscillations of the oscillator network.
Proof. By Theorem 2.1 the origin is locally exponentially stable for λ < 0. We further observe that, if the interconnection topology is strongly connected, then zero is a simple eigenvalue of L and therefore no either eigenvalue of J satisfies the same eigenvalue problem defined by Equation (5). It then follows by the center manifold theorem [41] and Equation (6), that the system (1) possesses a two-dimensional center manifold W c that is tangent to the consensus space C, for λ = 0. Moreover, this center manifold is exponentially attractive. By the Hopf bifurcation theorem [42], Equation (6) also implies that the system Equation (1) undergoes a supercritical Hopf bifurcation inside W c when λ crosses zero from negative to positive. By direct substitution inside the model equations, we see that along a generic member of the family of periodic orbits born at the Hopf bifurcation, oscillators are synchronously oscillating with each oscillator orbit given by (2). Remark 2.3. Because Hopf bifurcation is codimension-zero (in the sense of [43]), it is persistent under small perturbations, which ensures that Theorem 2.2 remains true for small heterogeneity in the natural frequencies.

Oscillator Death and Multi-Stability for Stronger Coupling
We now explore the phenomenon of "oscillator death, " induced by strong coupling in model (1). We restrict our attention to the all-to-all coupling case, i.e., a ij = 1 for all i = j. For λ > 0 and µ sufficiently small the synchronous oscillations born at Hopf bifurcation (Theorem 2.2) attract all trajectories. However, the system is multistable, as can be noted from the fact that increasing µ leads to the appearance of a family of steady states that attract some of the trajectories, but the synchronous periodic orbits remain locally exponentially stable (Figure 2). Indeed, depending on the initial conditions, only some trajectories converge to the synchronous oscillations. In the following we will provide geometric insights, without formal proof, about the mechanisms underlying oscillator death and multi-stability in model (1). We start by observing that the oscillator death state is characterized by the presence of two dead oscillator clusters.
Inside each cluster, oscillators converge to the same steady state. To analyze the appearance of oscillator death steady-states, we can simplify the model by assuming that (x i , y i ) = (x 1 , y 1 ) for all i = 1, . . . , N 1 and (x i , y i ) = (x 2 , y 2 ) for all i = 1, . . . , N 2 , where N 1 , N 2 < N, N 1 + N 2 = N, are the cluster sizes. The pairs (x 1 , y 1 ) and (x 2 , y 2 ) define the cluster states.
The cluster state dynamics can be easily derived and reaḋ Each cluster state dynamics has the forṁ where N j is the other cluster size and x j the other cluster state. A sufficient condition for the appearance of multiple steady states is that there must exist values of x j for which the model (8) has multiple steady-states. This condition can easily be verified by analyzing the dependence of the intersection of the nullclines of model (8) as a function of x j and the parameters µ and λ (Figure 3). If the coupling strength is too small the origin is the only steady state (Figure 3). This steady state is unstable and all trajectories are attracted toward the synchronous periodic orbit. However, new steady states appear for larger values of µ. The critical value of µ for which the new steady-states appear can be found by computing the slope of the nullclines at the origin. The slope of the x-nullcline is evidently µN j . The slope of the xnullcline can be computed by implicit differentiation and is given by ω 2 λ . Multiple steady-states appear if µ > ω 2 N j λ (Figure 2).

SYNCHRONIZATION AND FREQUENCY MODULATION IN LINEARLY COUPLED OSCILLATORS
In this section we present an alternative approach to study synchronization under the influence of noise using the Fokker-Planck equation (FPE). The modeling in this section can be thought of in the context of interacting populations of oscillators. Related work in the context of populations of synchronized neurons can be found in the work by Jiao et al. [44]. Introducing noise in the equation is natural from the biological perspective. However, even in the absence of noise, the introduction of random perturbations allows the extraction of information about the deterministic system. This is done by letting the perturbation amplitudes go to zero. To investigate the dependence of the synchronization frequency on the frequencies of the coupled oscillators, and the different coupling parameters, we assume  synchronization, which reduces the FPE equation to an equation in two variables. We divide this section in two parts. First, we consider a simple deterministic system in which the effect of coupling can be understood. In the second part, we randomly perturb a more general version of the previous model to show that the FP equation provides an approximation for the syncrhonization frequency, and obtain some insights on the effect of noise.
Let us then consider a system similar to the one already studiedẍ Notice that x(t) = sin(ωt) is still a solution for the uncoupled system (µ = 0), independently of the value of ν. This system has the advantage of allowing direct calculations around the limit cycle, which can be written explicitly.
If we linearize the equation and take ν << 1, we can neglect the contribution of the disipative term. The resulting linear system isẍ a ij x j , for i = 1, . . . , N. If all the a ij = 1, corresponding to a fully connected network of oscillators, the previous system can be reduced to the equation assuming that a synchronized regime is established. This assumption may not always be biologically realistic, but it allows us to obtain the common synchronization frequency in a simple way. Later on we consider the general case and recover this formula as a particular case. This provides an estimate for the synchronization frequency of Moreover, this reduction suggests that synchronized oscillatory behavior takes place for sufficiently small µ. That is, when Otherwise, exponentially large growth can be expected. Notice that unless the a ij are equal, the previous reasoning is not consistent and no conclusion can be drawn. We claim that introducing random perturbations and using the FP equation allows us to circumvent this difficulty and analyze the general case. This is the content of what follows. First of all, we write the system in the forṁ Notice that in the linearized regime, anologous to the reasoning for small ν in the previous example, we might naturally assume that Perturbing the equation with Brownian noise we have where the W ij are uncorrelated Brownian motions for i, j ∈ {1, ..., N}. The probability density, u(x 1 , ..., x n , y 1 , ..., y n , t) of the system being in the state x 1 , ..., x n , y 1 , ..., y n at time t satisfies the FP Equation where F is the vector field determined by the right hand side of the stochastic system. Looking for stationary solutions, i.e., u t = 0 and explicitly substituting F in terms of x and y, the equation becomes If we use the synchronization condition x 1 = ... = x n and y 1 = ... = y n , we obtain the equation we can write the Equation (11) as Assuming ε is small, it is reasonable to expect that the probability u will concentrate around the characteristic curves of the first order equation that are solutions to the systeṁ x = ny, By taking the scalar product with the vector (x/n, y/b) we obtain the relationẋ x n +ẏ y b = 0, or equivalently, which defines the characteristic curves as ellipses. In turn, interpreting these as curves in the phase portrait, the resulting solutions would correspond to periodic trajectories with frequency which provides an estimate for the synchronization frequency in terms of the original frequencies and the coupling parameters. Importantly, it shows that the synchronization frequency decreases with the coupling strength. In particular, formula Equation (12) can be used to study how coupled ultradian oscillations can give rise to circadian oscillations (Figure 4).

DISCUSSION AND SUMMARY
We have described, through basic geometrical analysis, the relationship between the dissipation coefficient, an intrinsic property of the oscillators we study, and the coupling strength µ in a network of diffusively coupled non-linear oscillators. Our analysis predicts the emergence of sustained oscillations for increasing values of the parameter λ in the system, only for a limited range of coupling strengths. It is reasonable to conjecture from this result, that there is a functional limit in the coupling strength for oscillating tissues in nature above which the tissue oscillations dies. To the best of our knowledge, it is the first time that diffusive coupling has been shown to be able to induce such oscillator death. We have also derived an estimation for the synchronization frequency of a linearly coupled network of non-linear oscillators in terms of the oscillator natural frequencies and the coupling parameters [Equation (12)]. The presented results are indeed local, that is, the synchronous oscillation is only locally asymptotically stable. The oscillators are not synchronized at the beginning of the simulations, but their spread in state space is very small to ensure convergence to the synchronous oscillation. For a larger spread, a more complex behavior is observed.
We believe that these results constitute predictions that, although possibly difficult to test experimentally, would be worth verifying in light of the existing evidence about the joint frequency modulation of activity between different tissues during the day [23,45].
The results we have presented thus far emphasize the importance of simple mathematical models in understanding situations where synchronization of multiple oscillating populations appears. The results presented here may help to shed light on both physiological and pathological phenomena involving synchronization of oscillators in different tissues (Parkinson's disease [46,47], epilepsy [48,49]). The other way around, it is also of potential importance to unravel mechanisms underlying the disappearance of coordinated oscillatory regimes. In a future publication, we plan to formally justify our estimations, and further, integrate the analysis of oscillations in the cellular and network levels of biological organization, to build up on our understanding of coupling oscillators at the tissue level. Two important extensions of the current models that we are studying are, the full characterization in higher codimension of the bifurcation structures of the system (1), and also, replacements of the van der Pol dynamics with biophysical models of excitable cells [50]. This last extension may prove useful to explain possible compensatory mechanisms that take place during the beginning of a pathology [51].