Collective Oscillations in Coupled-Cell Systems

We investigate oscillations in coupled systems. The methodology is based on the Hopf bifurcation theorem and a condition extended from the Routh–Hurwitz criterion. Such a condition leads to locating the bifurcation values of the parameters. With such an approach, we analyze a single-cell system modeling the minimal genetic negative feedback loop and the coupled-cell system composed by these single-cell systems. We study the oscillatory properties for these systems and compare these properties between the model with Hill-type repression and the one with protein-sequestration-based repression. As the parameters move from the Hopf bifurcation value for single cells to the one for coupled cells, we compute the eigenvalues of the linearized systems to obtain the magnitude of the collective frequency when the periodic solution of the coupled-cell system is generated. Extending from this information on the parameter values, we further compute and compare the collective frequency for the coupled-cell system and the average frequency of the decoupled individual cells. To compare these scenarios with other biological oscillators, we perform parallel analysis and computations on a segmentation clock model.


Introduction
Biological rhythms play important roles in nature. They include a wide variety of cyclical behaviors which are crucial in organisms, from development to homeostasis. The corresponding periodicity ranges from microseconds for cellular oscillations to seasons or years for recurring movements in ecological systems. Scientists have been continuously making concerted efforts to understand the mechanisms of synchrony, robustness, and sustainability of these oscillations, via mathematical and computational modeling, and experiments (Fall et al. 2002;Forger and Peskin 2003;Gonze 2011;Winfree 1980).
In mammals, a biological clock located in the suprachiasmatic nucleus (SCN) drives remarkably precise circadian rhythm. This master circadian clock is composed of multiple oscillator neurons that are coordinated through molecular regulation (Antle and Silver 2005;Bell-Pedersen et al. 2015). While individual cells oscillate with periods ranging from 20 to 28 hours, the collective rhythm in circadian clock is synchronized through intercellular coupling mediated by neurotransmitters. Experimental evidence shows that intercellular signal vasoactive intestinal polypeptide (VIP) expressed by some of the SCN neurons is such a major neurotransmitter (Aton et al. 2005). A special feature of the master circadian clock within the SCN is that the collective frequency of the synchronized coupled cells is close to the mean of the intrinsic frequencies of individual cells (Aton et al. 2005;Herzog et al. 2004;Honma et al. 2012;Liu et al. 1997;Taylor 2014).
The above-mentioned intricate biological mechanism involves some interesting features and evokes two issues in modeling: When a population of oscillatory cells is coupled, can a collective oscillation be produced? And will a collective oscillation be generated only under sufficiently large coupling strength? Herein, we consider that each individual cell oscillates at its own frequency when uncoupled. We call it a collective oscillation for these cells under coupling if a common frequency of oscillation is attained for all cells. Such common frequency was termed compromise frequency in a pair of phase equations in Strogatz (1994). Therein, such a collective oscillation is generated if the coupling strength is larger than half of the difference between the two individual frequencies. Herein, for generality, we call such common frequency of oscillation the collective frequency in the coupled-cell systems. The second issue is about the closeness between the collective frequency and the mean of the intrinsic frequencies of individual cells.
Denote by ω i the individual frequency of the ith oscillator among the N oscillators, when isolated. Let us call it the average frequency property if the collective frequency of the coupled cells is equal to or close to the average frequency ω Ave , where We note that average frequency is disparate from average period T Ave , and they were sometimes mixed up in the literature. The average period T Ave is given by where T i := 2π/ω i are the individual periods. For the average frequency ω Ave in (1), the corresponding period is which is not equal to T Ave in (2) in general. Certainly, if all ω i are close to each other, say ω i ≈ ω for all i = 1, . . . , N , then the average frequency and the average period are about the same thing, as A mathematical model which well depicts the average frequency property is the coupled phase equations: where θ i is the oscillatory phase of cell i, ω i is the intrinsic frequency of the ith cell, K ji ∈ R is the connection weight, f is an interaction function, and N is the number of cells. System (3) is known as the Kuramoto model (Kuramoto 1984) and has been studied extensively (Chiba 2015;Ha et al. 2016). As pointed out in Liu et al. (1997), any phase-locked solution of (3) oscillates at the mean frequency N i=1 ω i /N , provided that f is an odd function and [K ji ] is a symmetric matrix. This can be seen by adding up all components in (3). However, odd function and symmetric matrix are quite special among all possible interaction functions and connection matrices.
Phase equations are considered from focusing on the collective behavior in terms of the phase of oscillation in a collection of clock cells, instead of concentrating on the internal machinery of cells. However, it is significantly interesting to understand how oscillations are generated in cells. It has been identified that the intracellular transcriptional/translational negative feedback loops between activators and repressors are the key oscillatory mechanism in mammals and other organisms. Therefore, it is appealing to see whether the kinetic models based on such negative feedback loops accommodate the average frequency property.
Investigations on the average frequency property in some kinetic models have been reported in Gonze et al. (2005), Kim (2016), Kim et al. (2014). Recently, in a minimal genetic system constituted solely by a negative feedback loop, two types of gene regulation were studied and compared (Kim 2016;Kim et al. 2014). Therein, a single cell is modeled by where M, P, R are interpreted as the concentrations of a clock gene mRNA, the corresponding protein, and a transcriptional inhibitor, respectively, and the negative feedback loop is realized by the repression function f . The repression considered therein is either of Hill-function type or based on protein sequestration. System (4) is known as the Goodwin model, when the nonlinearity f is a Hill function (Goodwin 1965;Griffith 1968), i.e., f = f 1 , where k H > 0 is the half-saturation constant, and positive integer n is the Hill coefficient.
A graph of f 1 is depicted in Fig. 1a. The Goodwin model has been a prototypical system for accounting the core molecular mechanism associated with generation of self-sustained oscillations. It has been adopted to study circadian clocks (François et al. 2012;Ruoff et al. 1999Ruoff et al. , 2001. Hill functions were largely employed to describe cooperative binding of repressors to the gene promotor in transcription (Keller 1995) or repression based on multisite phosphorylation mechanism (Gonze and Abou-Jaoudé 2013). Hill-function-type repression has been widely adopted in various models for biological rhythms (Invernizzi and Treu 1991; Kim and Forger 2012;Kim 2016;Iwasa 2002, 2005).
Recently, another mechanism of transcriptional repression based on protein sequestration has been proposed and adopted into the negative feedback modeling. In such reaction, repressor protein R binds free activator A into inactive complex. The fraction of activators that are not sequestered is expressed by where k d > 0 is the dissociation constant between the activator and repressor (Buchler and Louis 2008;Buchler and Cross 2009;Kim 2016;Kim et al. 2014). A graph of f 2 is depicted in Fig. 1b. As transcription is in a ratio to this fraction, f 2 (R) was taken as f (R) in system (4) to depict the negative feedback loop. It was reported in Kim and Forger (2012) that tight binding between activators and repressors and balanced stoichiometry are the key for sustained rhythms in a detailed mathematical model on the mammalian circadian clock, and such property is also reflected in the simplified model, system (4) with f = f 2 . Goodwin's model (system (4) with f = f 1 ) was proposed in Goodwin (1965). It was proved in Griffith (1968) that the equilibrium is stable for n ≤ 8, by the Routh-Hurwitz criterion. When n > 8, it was shown that some parameter values can be found at which the equilibrium is unstable. In , for a system similar to system (4), with f = f 1 , it was shown that the equilibrium is stable for n ≤ 8, by the Routh-Hurwitz criterion, and for n > 8, the equilibrium can be unstable for some specially chosen parameter values. In Fall et al. (2002), under the assumption of equal degradation rates, it was justified that the equilibrium is unstable if n ≥ 8, if the degradation rates are small. We note that the steady state depends on the parameters, and so at the bifurcation value, the existence of equilibrium needs to be assured to confirm the occurrence of Hopf bifurcation. This was not addressed in those works. In Woller et al. (2014), under the assumption of equal degradation rates, the steady state at which the linearized system has a pair of purely imaginary eigenvalues was found, where the condition n > 8 is explicitly revealed, yet, the crossing condition was not mentioned. In fact, n > 8 is both a sufficient and necessary condition for the simple Hopf bifurcation. Our formulation considers general repression functions which accommodate both f 1 and f 2 .
Via numerical computations on the coupled-cell system and analyzing the phase response curve, it was asserted in Kim (2016), Kim et al. (2014) that the average frequency property holds at reasonable parameter values for the system with repression based on protein sequestration, whereas the collective frequency is far from the mean if modeled with Hill-type regulation. The mathematical models adopted therein are coupled-cell system comprising subsystems of the form (4) with f = f 1 , f 2 , respectively. It was also reported in Kim (2016), Kim and Forger (2012) that the properties of such coupled-cell system with protein-sequestration-based repression are consistent with data from Drosophila and mammals, whereas the properties of such coupledcell system with Hill-type repression match well with the experimental data from Neurospora. It is therefore interesting to perform a further mathematical study to see and compare the oscillatory properties between the coupled-cell system modeled with f = f 1 and the one modeled with f = f 2 .
Connection between the kinetic models and the phase models has been built in the so-called weakly coupled oscillators theory, cf. Kopell 1984, 1991;Schwemmer and Lewis 2012). However, such connection requires certain assumptions and a transparent correspondence between the dynamics of the kinetic models and its phase equation counterpart remains to be further explored. While weak coupling is a prerequisite in the theory, sufficiently large coupling strength is needed for the existence of phase-locked solution in some phase models. Thus, how weak is weak so that these connection and correspondence are valid appears to be a complicated issue. Recently, it has been reported that phase-amplitude reduction and higher-order reduction, instead of merely phase reduction and linear approximation, are necessary to capture the dynamics in some oscillatory systems, as reviewed in Ermentrout et al. (2019), Kuramoto and Nakao (2019).
The goal of our study is to explore the oscillatory properties of coupled-cell systems modeling biological rhythms. We shall consider the coupled-cell systems which are composed by the single-cell subsystem (4) with f = f 1 , f 2 , respectively. On the one hand, system (4) models the minimal genetic negative feedback loop in single cells, and it is essential to investigate oscillations in the coupled-cell system which comprises subsystems (4). On the other hand, the above-mentioned two issues about the properties of coupled oscillators and circadian clocks deserve a close analysis. The other important concern is about the comparison between the dynamics in the model using Hill-function repression and the one in the model based on protein sequestration repression. Although the mathematics in the kinetic models is in general rather involved, as commented in Baker and Schnell (2009), we take on the challenge to develop efficacious mathematical approaches to analyze the coupled-cell systems to see the dynamical details. In addition, the comparison between indication from the kinetic models and the one from the phase models can be made only after the kinetic models are sufficiently understood.
The first task is to confirm the existence of periodic solutions in the single-cell systems and the coupled-cell systems. Our methodology is based on the Hopf bifurcation theory (Hassard et al. 1981) and a sufficient-and-necessary condition for the simple Hopf bifurcation derived in Liu (1994), which was an extension of the Routh-Hurwitz criterion. As the parameter values vary, we seek for the situation that the determinant D m−1 of the second-to-last Hurwitz matrix H m−1 decreases from positive to zero, for a polynomial of degree m. A pair of complex-conjugate eigenvalues then cross the imaginary axis of the complex plane. This tracing allows to locate the Hopf bifurcation (HB) values.
The criterion in Liu (1994) for detecting Hopf bifurcation has been applied in several works on mathematical biology. For instance, for models on immune responses to persistent viruses, some ODE systems up to five dimension were considered in Liu (1997). In Domijan and Kirkilionis (2009), bistability and oscillations in chemical reaction networks were reported and the criterion was applied to a four-dimensional system of ODEs. In a study on extracellular signal-regulated kinase (Obatake et al. 2019), a polynomial of degree seven was tested to meet the criterion. In those applications, a common predicament is that the linearized system depends on the values of steady states and parameters. Therefore, the basic task is to seek for suitable parameter values, and hence, the steady-state value so that the associated characteristic polynomial meets the criterion. To this end, Newton polytope and symbolic computation using MAPLE programming were employed for such verification in Obatake et al. (2019) and Liu (1994), respectively. On the other hand, numerical bifurcation software, such as AUTO, is also able to locate the HB values and draw the HB curves in the parameter space.
Even with the approach by numerical computation, the condition D m−1 = 0 is also helpful for detecting and confirming the HB values. In this work, we will see from this condition that the Hopf bifurcation does not occur at those parameter values and steady states where f (x) is too small (x is one of the component of the steady state). In addition, the HB value for the single-cell system is larger than the HB value for the coupled-cell system, in the identical-cell case. More information about the Hopf bifurcation and the frequency of bifurcating periodic solution can also be observed from this condition.
While most of the literature about the application of Hopf bifurcation, including the above-mentioned ones, is only concerned with existence of periodic solutions, the Hopf bifurcation theory actually implicates more. It provides a foundation for our concern about the properties for the frequencies of oscillations. At the HB value, the magnitude of the purely imaginary eigenvalue gives the frequency of the emergent periodic solution. As this purely imaginary eigenvalue is a root of the characteristic polynomial, for the identical-cell case, we can link the frequency of oscillation to the parameters and steady states, at the HB values. It is interesting to observe that for the single-cell systems, the repression function f does not play a role in the size of frequency at the HB value. In contrast, γ := − f (x) is a factor to the frequency at the HB value in the coupled-cell systems, and thus, the modeling with f = f 1 or f = f 2 is distinguishable in this respect.
The paper is organized as follows. In Sect. 2, we introduce the Hopf bifurcation theory and the degenerate Routh-Hurwitz criterion to study the existence and stability of periodic solutions. We apply these theories to obtain periodic solutions in the singlecell system (4) with f = f 1 and f 2 , respectively, in Sect. 3. In Sect. 4, we discuss the collective frequency of periodic solutions for the coupled-cell system and compare it with the average of individual frequencies. In addition, we make a comparison between the system with Hill-function repression and the one with protein-sequestration-based repression. To compare with other biological oscillators, in Sect. 5, we conduct similar analysis on a segmentation clock model and make a comparison between these different types of biological oscillators. Some justifications and computations are arranged in Appendices and Supplementary Materials.

Hopf Bifurcation and Degenerate Routh-Hurwitz Criterion
In this section, we review the Hopf bifurcation theory and an extension of the Routh-Hurwitz criterion, which are to be applied to investigate the stable periodic solutions and their periods for the single-cell systems and the coupled-cell systems. Hopf bifurcation theory not only confirms the existence of periodic solution, but also indicates that at the Hopf bifurcation value, the purely imaginary eigenvalue of the linearized system provides the frequency of the bifurcating periodic solution. The Routh-Hurwitz criterion characterizes a polynomial whose roots all have negative real parts. We shall apply a degeneracy of this criterion, established in Liu (1994), that provides a condition for the Hopf bifurcation.
Let us consider an autonomous system of ODEṡ where x = (x 1 , . . . , x m ) ∈ R m , F = (F 1 , . . . , F m ), and p ∈ R k are parameters. Letx be an equilibrium of system (7) with a fixed p. Denote the Jacobian matrix associated with the linearization of system (7) atx by J (x) : The Hurwitz matrices associated with polynomial Δ(λ) are defined as The Routh-Hurwitz criterion indicates that all roots of (8) have negative real parts if and only if D i := det(H i ) > 0, i = 1, 2, . . . , m, cf. (Gantmacher 1959;Hurwitz 1895;Kemperman 1982). Note that D m = b m · D m−1 , and so the Routh-Hurwitz criterion can also be expressed by D i > 0, i = 1, 2, . . . , m − 1, with b m > 0. The following lemma completely characterizes a polynomial which has a pair of purely imaginary roots and negative real parts for all the remaining roots. It can be regarded as a degeneracy of the Routh-Hurwitz criterion. We note that the Routh-Hurwitz criterion has also been formulated as positive determinants of a sequence of matrices different from (9), see (Uspensky 1948). With D i defined as the determinants of such sequence of matrices, this lemma was proved in Liu (1994). By similar arguments, it can be justified that Lemma 1 holds true if we adopt D i , the determinant of Hurwitz matrices H i defined in (9).
To apply the Hopf bifurcation theorem, one first needs to find the parameter values at which a complex-conjugate pair of eigenvalues of J (x) cross the imaginary axis of the complex plane transversally and the real parts of the other eigenvalues remain negative. In a system modeling a biological process, there are usually a number of parameters. We choose one of the parameters as the bifurcation parameter, denoted by α in the following discussions, and fix the other parameters at suitable values. By doing so, we regard the equilibriumx as a function of α, i.e.,x =x(α). The Jacobian matrix associated with the linearization of (7) atx(α) then depends on α, i.e., J (x(α); α) and is abbreviated as J (x; α). The following Hopf bifurcation theorem was stated in Liu (1994), where instead of computing the eigenvalues of J (x; α), one only needs to compute the determinants of the Hurwitz matrices. Theorem 1 Consider system (7) whose equilibriumx =x(α) is a function of α for α near certain α * , when the other parameters are fixed. Assume that at α = α * , the following conditions hold Then, the system undergoes a Hopf bifurcation at x =x(α * ), and a small-amplitude periodic solution surroundingx emerges as α < α * or α > α * and α is close to α * .
The scenario in Theorem 1 was called simple Hopf bifurcation in Liu (1994), as there are more complicated Hopf bifurcations, see (Golubitsky and Schaeffer 1985). The value α * in Theorem 1 is called Hopf bifurcation (HB) value. The advantage of Theorem 1 is that its conditions are transparent and computable, and the existence of periodic solutions is guaranteed. When m = 3, 4, the conditions of the theorem can be clearly expressed by To apply the Hopf bifurcation theory to the ODE systems, we look for values of α * and values of the other parameters, which satisfy the conditions of Theorem 1. This process relies on some mathematical formulations and numerical computations.
Hopf bifurcation theory actually provides more information, including the stability and period of the bifurcating periodic solution. The period (frequency) is especially the focus of the present investigation. Stability of the bifurcating periodic solutions is determined by some higher-order terms of the system. One first transforms the system into normal form and then applies the center manifold theorem to obtain these terms. While this formulation is standard, its computation is cumbersome for large m. We summarize this process in Supplementary Material I, where the terms g 20 , g 11 , g 02 , and g 21 , and hence C 1 (α * ), p 2 , ζ 2 , and T 2 are introduced. Let λ(α) be the branch of eigenvalue crossing the imaginary axis at α = α * , and λ(α * ) = iω * , with ω * > 0. The following Hopf bifurcation theorem can be found in Hassard et al. (1981).
We note that numerical bifurcation software such as AUTO can also detect the Hopf bifurcation and draw the HB curve which comprises the HB values in the parameter space. However, the condition in Theorem 1 allows us to analyze how Hopf bifurcation occurs and locate and confirm the HB values via numerical computations. In addition, from Theorem 2, the frequency of the bifurcating periodic solution is given by the magnitude of the purely imaginary eigenvalue which is the root of the characteristic polynomial. One can thus link the frequency to the parameters and the HB values.

Single-Cell HT Model and PS Model
In this section, we consider system (4), abbreviated as HT model for f = f 1 , the Hill-type repression, and as PS model for f = f 2 , the protein-sequestration-based repression. We shall apply Theorem 1 in Sect. 2 to locate the parameter values at which the Hopf bifurcation takes place in system (4). Near these bifurcation values, Theorem 2 provides the approximate periods for the bifurcating stable periodic solutions.
To analyze oscillatory properties for these models via Hopf bifurcation, herein, we consider another change of variables. Under the same assumption (14), we set Then, system (4) can be transformed into ⎧ ⎪ ⎨ ⎪ ⎩ẋ where α := α 1 α 2 α 3 /β 3 is to serve as the bifurcation parameter.
We call system (16) with f = f 1 the single-cell HT model. In this system, the ratio of rates α, the Hill coefficient n, and dissociation constant k H between the repressor and gene promoter determine the dynamics. We call system (16) with f = f 2 the singlecell PS model, where the rate ratio α, activator concentration A, and dissociation constant k d between the activator and repressor, determine the dynamics. On the one hand, system (16) with α = 1 reduces to system (15). On the other hand, dynamical properties of system (16) certainly correspond to the kinetics in system (4). System in the form (16) with f = f 1 has also been studied in Woller et al. (2014). We now take α as the bifurcation parameter and apply Theorems 1 and 2 to investigate the periodic solutions of system (16).
Note thatx = (x 1 ,x 2 ,x 3 ) is a positive equilibrium of system (16) if and only if x 1 =x 2 =x 3 =x > 0, and Suchx uniquely exists for f = f 1 or f = f 2 : Proposition 1 There exists a unique positive equilibriumx = (x,x,x) for system (16) with f = f 1 , for any n ≥ 1, k H , α > 0, and for system (16) with f = f 2 , for any A, k d , α > 0, wherex is the unique positive solution to α f 1 (x) =x, and α f 2 (x) =x, respectively. Moreover, for any fixed n ≥ 1, k H > 0, or fixed A, k d > 0,x is an increasing function of α.
Proof It is obvious that (x 1 ,x 2 ,x 3 ) is an equilibrium for system (16) if and only if For f = f 1 , this reads as i.e.,x satisfies q(ξ ) = αk n H , where q(ξ ) := ξ n+1 + k n H ξ . For any n ≥ 1, k H , α > 0, there exists exactly one positive solution to this equation, due to that q is strictly increasing on [0, ∞), q(ξ ) → ∞ as ξ → ∞, and q(0) = 0 < αk n H . We also see that x is a function of α, as q is strictly increasing. For As ξ/ f 2 (ξ ) is strictly increasing and has an inverse,x is thus an increasing function of α.
Denotex =x(α) to indicate the dependence on α. We note that from (17), α can be expressed as a function ofx, i.e., α =x/ f (x). This is the inverse expression of x(α). We shall analyze the periodic solutions bifurcating fromx. The Jacobian matrix associated with the linearization of system (16) atx is One can factorize this cubic polynomial and find all its roots. Certainly, we can also apply Theorem 1. From (10), we see that Δ(λ) has a pair of purely imaginary roots and a negative root if and only if Notice that γ depends onx which is a function of α. Hence, whether the Hopf bifurcation takes place is determined by the existence of solutionx to If suchx exists, then the Hopf bifurcation occurs at equilibrium (x,x,x) and α = α 0 = −8/ f (x). Let us elaborate to confirm such existence for f = f 1 and f = f 2 , respectively. Since γ will play an important role in the analysis, we denote γ 1 := − f 1 (x) and γ 2 := − f 2 (x).
Remark 1 (i) If we multiply each component in system (16) by a factor σ , or change the time from t to t/σ , then the eigenvalues of the linearized system atx become λ 1,2 = ±i √ 3σ, λ 3 = −3σ . The system still undergoes a Hopf bifurcation at α = α 0 and x =x, and the frequency of the bifurcating periodic solution is near σ √ 3. This property will be used in Sect. 4.2 for a cell-to-cell system which has two different individual frequencies.
(ii) For the general situation of degradation rates, i.e., when assumption (14) does not hold, assertions similar to the ones in Theorem 3 still hold. In particular, the frequency ω 0 at the HB value is a combination of the degradation rates (no longer √ 3), but still does not depend on γ , and hence does not differentiate the form of repression function f . In addition, for f = f 1 , n will be required to be much larger than 8, if the difference among β 1 , β 2 , β 3 is large. This was mentioned in (Fall et al. 2002) and is shown precisely in Shih and Yang (2021).
The following examples illustrate Theorems 2 and 3. The parameter values considered here are mostly taken from Kim et al. (2014). (16) with f = f 1 and parameter values n = 11 and k H = 0.136. The graph of f 1 is depicted in Fig. 1a. According to Theorem 3, as n > 8, the Hopf bifurcation occurs atx = (x,x,x), wherex ≈ 0.148684, computed from (23). Next, we compute to find γ 1 ≈ 14.674227 from (22). Accordingly, a smallamplitude periodic solution emerges as the value of α increases through α 0 = 8/γ 1 ≈ 0.545174, with frequency about √ 3 ≈ 1.732051. Furthermore, we compute to find that p 2 > 0, ζ 2 < 0, T 2 > 0. The numerics are shown in Supplementary Material I. According to Theorem 2, system (16) with f = f 1 undergoes a supercritical Hopf bifurcation at α = α 0 andx = (x,x,x). The bifurcating periodic solution is stable and the period increases as α increases. We numerically solve system (16) and compute the frequencies and amplitudes of the periodic solutions corresponding to various values of μ := α − α 0 , plotted in Fig. 2. It can be seen that as μ increases from 0, the frequency decreases from about √ 3 (i.e., the period is increasing), which is consistent with the assertion of Theorem 2.

Example 3.1 Consider system
Next, we adopt the parameter values in Kim et al. (2014) and illustrate that the numerically computed oscillations therein appear to be a continuation of the periodic solutions generated by the Hopf bifurcation. (16) with f = f 2 , with parameter values A = 0.0659 and k d = 0.00001. The graph of f 2 is depicted in Fig. 1b. As condition (25) in Theorem 3 is met, the Hopf bifurcation occurs atx = (x,x,x), withx ≈ 0.058730 computed from (26). We compute to find γ 1 ≈ 14.986634 from (24), and thus α 0 = 8/γ 1 ≈ 0.533809. Therefore, a small-amplitude periodic solution emerges, as the value of α passes through α 0 , with frequency about √ 3. Furthermore, we compute to find that p 2 > 0, ζ 2 < 0, T 2 > 0. The numerics are revealed in Supplementary Material I. According to Theorem 2, system (16) with f = f 2 undergoes a supercritical Hopf bifurcation at α = α 0 andx = (x,x,x). The bifurcating periodic solution is stable and the period increases as α increases. The frequencies and amplitudes corresponding to various values of μ := α − α 0 are plotted in Fig. 3. When α = 1 (μ = 0.466191), system (16) becomes (15), and the parameter values herein are exactly the ones adopted in Kim et al. (2014). The system with such parameter values generates an oscillation with frequency about 1.664346.

Example 3.2 Consider system
A comparison and summary Let us summarize the above discussions, and make a comparison between the single-cell HT model and PS model. For either f = f 1 or f = f 2 , system (16) has a unique equilibriumx = (x,x,x). For each set of parameters with n > 8 for the case f = f 1 , or satisfying (25) for the case f = f 2 , the Hopf bifurcation occurs at (α 0 ,x). The Jacobian matrix associated with the linearization of system (16) atx with α = α 0 is a constant matrix. Restated, it is always this matrix whenever the Hopf bifurcation takes place in system (16), with f = f 1 or f 2 , or other smooth decreasing functions. Subsequently, the frequency of the bifurcating periodic solution near the bifurcation point is about √ 3. As α increases, the frequency decreases and the amplitude increases in both models, but the slopes are slightly different between these two models. The frequency drops faster in the single-cell PS model, whereas the amplitude increases faster in the single-cell HT model, see Figs (16) with f = f 1 is pretty close to the one with f = f 2 , as demonstrated in Fig. 4. This is a consequence of Theorem 3.
Remark 2 Recall that system (15) is identical to system (16) with α = 1. In Kim et al. (2014), system (15) was considered with n = 11, k H = 0.04 in the HT model and A = 0.0659, k d = 0.00001 in the PS model. While the amplitudes in the two models with these parameter values are close to each other, our computation shows that the period for the single-cell HT model is 3.9302912, whereas the period for the single-cell PS model is 3.7751695. Therefore, we take k H = 0.136 in the HT model in Example 3.1, and we can tune α so that the oscillatory wave forms in these two models are similar at α = 0.60139.

Coupled-Cell HT Model and PS Model
The synchronous or coherent rhythmicity for a collection of clock neurons is mediated by the neurotransmitters. It has been identified that vasoactive intestinal polypeptide (VIP) is the key synchronizing agent in several experiments (Aton et al. 2005;To et al. 2007). The collective period of oscillation is generated through intercellular coupling which is determined by the concentrations of neurotransmitters in extracellular medium. Modeling with such intercellular signaling via VIP, the following coupled-cell system has been considered in Kim et al. (2014) to investigate the average frequency property: where i = 1, 2, . . . , N , and parameters c and s are the coupling strength and the timescale of intercellular coupling, respectively. In this model, each cell releases VIP into the extracellular space at a rate proportional to the activity of the promoter f (R). The fourth component of the system indicates the release of V (VIP) by each cell into the extracellular space at a rate proportional to f (R) which describes the activity of the promoter. Note that the experimental data indicate that the intercellular coupling strength occurs much faster, compared to the intracellular feedback loop . This model is also based on the fact that the release of VIP is fast with respect to the 24 hours timescale, and thus the mean field expression ( N j=1 V j )/N is adopted. When c = 0, the cells are decoupled, and we consider the parameter values at which each cell oscillates at the same frequency, as these N subsystems are identical. As c increases from 0, the cells become coupled. However, whether a collective periodic solution of the coupled system is then generated requires a justification.
In this section, we consider the coupled-cell HT model ( f = f 1 ) and the coupledcell PS model ( f = f 2 ). We apply Theorem 1 to confirm the existence of periodic solution for these coupled-cell systems. We shall investigate how coupling strength affects the collective oscillations, and whether and how average frequency property holds in these two coupled systems. Referring to the frequency of periodic solution indicated in Theorem 2, we trace the eigenvalues of the linearized systems from the HB value for the single-cell subsystems (c = 0) to the HB value for the coupled-cell system, and see which eigenvalue branch reaches the imaginary axis of the complex plane. Also, we compare the variations of collective frequency in these two models as the coupling strength c increases. In Sect. 4.1, we discuss the coupled system comprising two identical cells. This analysis is also valid for the coupled system comprising N identical cells, i.e., system (29). The coupled system comprising two nonidentical cells will be addressed in Sect. 4.2.
We note that synchronization in a more complicated model with Michaelian kinetics for the degradations and coupling through the mean field, has been reported in Gonze et al. (2005).

Identical Cells
In this section, we consider coupled system (29) which consists of two identical cells, expressed by where parameter c > 0 represents the coupling strength and s > 0 is the timescale of intercellular coupling. The discussion herein is also valid for system (29) which comprises N identical cells, as remarked below. Denote x = (x 1 , x 2 , x 3 , x 4 ), y = (y 1 , y 2 , y 3 , y 4 ), and X = (x, y). The following proposition about the equilibrium is straightforward.
Proposition 2 For each set of parameter values, system (30) has a unique equilibrium X which is homogeneous, i.e., We show that system (30) has only homogeneous equilibrium X = (x,x) in Appendix B(I). Now we take α as the bifurcation parameter, while holding all other parameters fixed at suitable values, and apply Theorem 1 to confirm the existence of periodic solution bifurcating from X. The Jacobian matrix associated with the linearization of system (30) at X is can be factored as where This factorization can be made since system (30) consists of two identical subsystems and the synchronous set That is, if D − imaginary roots and six roots with negative real parts if and only if D + 3 (α * ) = 0, by Lemma 1. From this equality, we find where the square root is positive if and only if Herein, we only consider the case s > 1, and hence only the "+" one in (32). Then, Note that s > 1 implies Thus, (34) implies (33), provided s > 1. The case 0 < s ≤ 1 can also be discussed. For crossing condition in (13), we consider , which is (s + 3)(s 2 + 3) multiplied to the term for determining the crossing condition in the decoupled single-cell case (27), discussed in Sect. 3, see Appendix A. Hence, by continuity, (35) is nonzero for small coupling strength c, under the condition in Theorem 3. We confirm the existence of homogeneous equilibrium X = (x,x) at the bifurcation value α = α * for system (30) with sufficiently small coupling strength c in Appendix B(II). By substituting α * in (32) into the fourth-degree polynomial Δ + (λ), we can actually find its purely imaginary roots ±iω * c .
Theorem 4 Consider s > 1 and assume that the equilibrium X = (x,x) exists at α = α * defined in (32), and that (34) and the crossing condition hold. System (30) undergoes a Hopf bifurcation at α * and a small-amplitude periodic solution near X emerges as α < α * or α > α * and α is close to α * , with frequency ω c about Notably, it can be shown that α * > 0 if and only if ω * c > 0. Certainly, the collective period of the bifurcating periodic solution is near T * c = 2π/ω * c . We will provide further observation from the assertion of this theorem below.
Implementation and Computation First, we fix the parameter values of n and k H in the HT model, and A and k d in the PS model. For fixed values of c and s, we (32), and express α * in terms ofx. Then we substitute α = α * into (α + c) f (x) =x to solve forx, and then obtain the equilibrium, as indicated in Proposition 2. Next, we substitutex into γ , and then compute α * in (32) to confirm that it is positive or check condition (34). Note that for the PS model, there are two values ofx and we choose the smaller one to have smaller value of α * . By examining the crossing condition, we confirm that the assertion of Theorem 4 holds. We denote such c by c * and still call (c * , α * ) a HB value. These HB values (c * , α * ) form the Hopf bifurcation curve in the (c, α)-plane.
According to Theorem 4, a small-amplitude periodic solution emerges, at α > α * or α < α * and α close to α * , with frequency ω c about ω * c . We are interested in seeing how the collective frequency ω c of oscillation varies with the coupling strength c and parameter α in system (30). Let us pick one of the HB values and denote it by P 1 = P 1 (c * , α * ). Recall that the Hopf bifurcation occurs at α = α 0 in the single-cell systems.
Variations of eigenvalues At (c, α) = P 0 (0, α 0 ), there are two pairs of purely imaginary eigenvalues for the linearized system at X(α 0 ). For the identical-cell case, system (30), these two pairs coincide. It is interesting to see how the two complexconjugate branches emanating from these two pairs move in complex plane C, and which of them reaches the imaginary axis again at P 1 (c * , α * ). As the magnitudes of the purely imaginary eigenvalues correspond to the frequencies of the individual cells at P 0 and coupled-cell at P 1 , we can observe the transition of frequency from the variation of eigenvalues from P 0 to P 1 . More precisely, it is interesting to see, as c increases from 0 so that the coupling becomes effective and a collective periodic solution is formed, the relative positions of the eigenvalues at P 1 with respect to the eigenvalues at P 0 . This reflects a transition from single-cell oscillation to coupled-cell oscillation. Moreover, we can also see how the local dynamics around X is changing from the variations of eigenvalues. Certainly, there are infinitely many paths from P 0 to P 1 in (c, α)-plane. For simplicity and illustration, we take the line segment from P 0 to P 1 to trace the variations of eigenvalues.
To confirm that the Hopf bifurcation is supercritical so that the emergent periodic solution is stable and see whether it occurs for α > α * or α < α * , we take another line segment in (c, α)-plane by fixing c = c * , and allowing α to increase from below α * to above α * . We call such segment HB course. We can compute the eigenvalues of J (X; α) to see if there is a stability switch for equilibrium X when α crosses α * , along such HB courses.
Let us illustrate the implementation of these ideas by the following examples. We take the same parameter values of n and k H in Example 3.1, and A and k d in Example 3.2. Recall that the frequency for the isolated (c = 0) individual cell is approximately ω 0 = √ 3, when α is close to α 0 . In the following example, we fix s = 20. These numerical results are carried out by MATLAB programming and are consistent with the computations by software AUTO.
Example 4.1.1 (i) Consider coupled-cell HT model (30) with parameter values n = 11, k H = 0.136, and s = 20. With such parameter values, the single-cell system undergoes a Hopf bifurcation at α 0 ≈ 0.54174, as shown in Example 3.1. Through computation, the HB curve is drawn in Fig. 6. For illustration, we take P 1 (c * = 0.05, α * = 0.4742) on the curve. At this HB value, we compute to find γ 1 ≈ 15.05911, and thus confirm that condition (34) in Theorem 4 holds. Note that the line segment P 0 P 1 lies above the HB curve, as shown in Fig. 6. We demonstrate the synchronous periodic solution at parameter values near P 1 in Fig. 5.
Variation of eigenvalues in C along the HB course in Fig. 6 is plotted by the curve in Fig. 7b, c. We thus confirm that a periodic solution emerges as α passes slightly over α * , with frequency ω c close to ω * c ≈ 1.721229, by Theorem 4. Our numerical simulation shows that for α slightly above α * , the collective frequency is ω c ≈ 1.71999 for (c, α) = (0.05, 0.475).
Variations of eigenvalues along the HB course are shown by the curves in Fig. 9b, c. According to Theorem 4, a periodic solution emerges for α > α * and close to α * , with frequency ω c close to ω * c ≈ 1.721278. Our numerical simulation reveals that for α slightly above α * = 0.4766, the collective frequency is ω c ≈ 1.72115 at (c, α) = (0.05, 0.498).
We carry out similar computation for the system at the other HB values in Figs. 6 and 8. The scenarios, the variations of eigenvalues along line segments from P 0 , resemble the ones for P 0 P 1 in Figs. 7a, c and 9a, c. The variation of eigenvalues along each associated HB course is also similar to the one for P 1 in Figs. 7b, c and 9b, c.
Let us make a comparison about variations of eigenvalues along the line segments between the coupled-cell HT model and PS model. In the HT model, λ 1 -branch moves to the right half plane of C and then makes a turn downward to reach the imaginary axis at P 1 , whereas λ 2 -branch remains on the left complex plane. In the PS model, both λ 1 -branch and λ 2 -branch move into the left complex plane. Then, λ 1 -branch makes a turn downward and reaches the imaginary axis at P 1 , and λ 2 -branch continues to stay in the left complex plane. The basic reason for the different movement of λ 1 -branch is that line segment P 0 P 1 is above the HB curve for the HT model, whereas it lies below the HB curve in the PS model, as indicated in Figs. 6 and 8, respectively. Accordingly, there is a difference on the transition of the local dynamics around equilibrium X between these two models. For the HT model, the dimension of the unstable manifold of X increases from zero to two (the dimension of the stable manifold from four to six) when (c, α) moves from P 0 along P 0 P 1 , before reaching P 1 . On the other hand, for the PS model, the dimension of the stable manifold of X increases from four to eight when (c, α) moves from P 0 , before reaching P 1 . Other than that distinction, both models undergo a supercritical Hopf bifurcation at α = α * . As the stable periodic solutions of the coupled-cell systems emerge, the collective frequency at each (c, α) = (c * , α * ) is smaller than the individual frequencies at (c, α) = (0, α 0 ) in both models. Indeed, the purely imaginary eigenvalues λ 1 (P 1 ) lie below the purely imaginary eigenvalues λ 1 (P 0 ) in both models.

Remark 3
The intercellular coupling is relatively fast, and we took s = 20 in the above examples. We observe from system (30) that x 4 and y 4 may be regarded as being in quasi-steady state, i.e., x 4 (t) ≈ f (x 3 (t)), y 4 (t) ≈ f (y 3 (t)), and thus the first equation is approximatelyẋ Recalling the Hopf bifurcation in the single-cell system (16), we see that the occurrence of the Hopf bifurcation in the coupled system (30)   P 0 is a HB value for the single-cell system, P 1 is a HB value for the coupled-cell system. b Line segment P 0 P 1 in red line; the HB course through P 1 in blue line, from solid square to hollow square (Color figure online) Fig. 9 Variations of two complex eigenvalues as (c, α) moves along segment P 0 P 1 (red) and the course (blue) in Fig. 8: a λ 1 -branch along P 0 P 1 reaches the imaginary axis at parameter value P 1 (0.05, 0.4766). b λ 1 -branch along the course, from solid square to hollow square, crosses the imaginary axis at P 1 . c λ 2branches fall on the left complex plane along both P 0 P 1 and the course. A supercritical Hopf bifurcation occurs at α * ≈ 0.4766 (Color figure online) lines and that the frequency decreases with increasing coupling strength is linked to the frequency decrease with increasing α in the single-cell system.
Remark 4 (i) It is seen from (32) that the HB value α * for the coupled-cell system is smaller than the HB value α 0 for the single-cell system (when c = 0). Based on the Hopf bifurcation theorem, this order provides a potential parameter regime where oscillations exist in both single-cell and coupled-cell systems, as well as a potential regime where oscillations exist in the coupled-cell system, but not in the single-cell system.
(ii) We observe from (32) that α * → ∞ as γ → 0. That is, the Hopf bifurcation does not occur at the flat parts of the repression functions f 1 and f 2 (small values of γ 1 and γ 2 ), see Fig. 10. On the other hand, the Hopf bifurcation occurs when γ , the magnitude of f (x), is lower than the bound in (34). This condition was actually derived to guarantee that α * is positive. This bound is large when c is small. In our Example 4.1.1, the bound is above the maximal values of γ 1 and γ 2 , and thus, condition (34) is automatically met. (iii) From Theorem 4, we see that frequency ω * c decreases as γ increases. That is, in contrast to the single-cell case, the frequency of the periodic solution near the HB value distinguishes between γ = γ 1 and γ = γ 2 , and hence the repression function f 1 from f 2 . From the graphs of γ 1 and γ 2 in Fig. 10, we see that while γ 1 has larger maximum value, γ 2 has wider range of large values. Figure 10 is depicted with the data in Example 4.1.1. It is indirect to track the value of γ , as it depends onx which is determined by the parameter values. For the HT model, we may compare the point x M where γ 1 attains its maximum value with the componentx of the equilibrium at c = 0, that is, Their distance can be measured from When the coupling strength c is small,x in the coupled system is close tox in (36), and the quantity in (37) estimates how far the pointx that γ 1 takes on is away from x M , where the maximum value of γ 1 is attained. Note that for n < 17,x lies on the right hand side of x M , whereas it is reverse for n > 17. In addition, we see that the gap between x M andx can be enlarged, by choosing larger k H and n. The corresponding frequencies ω * c will then further distinguish the repression function f 1 from f 2 . The analysis for the graphs of γ 1 and γ 2 is arranged in Appendix C. From the expression for frequency ω * c , we also see that large coupling strength c enhances the effect from γ , whereas large intercellular coupling time scale s suppresses the factor from γ . In Appendix D, we provide more numerics for the values ofx, γ 1 , γ 2 , and the corresponding frequencies ω * c , along the HB curves for the cases with s = 10, and n = 18.
What are discussed in Remark 4 are about the scenario at the HB points and HB values. To extend the understanding of the frequency, we perform further numerical simulations and follow the continuation from the HB values. We are interested in comparing the collective frequency ω c and the average frequency ω Ave of the individual cells. For system (30) comprising identical cells, the average frequency coincides with the individual frequency. To this end, we need to choose parameter values at which the periodic solutions exist for both of the single-cell and coupled-cell systems. In addition, to make comparison, we consider the same values of α in both single-cell and coupled-cell systems. As mentioned in Remark 4(i), the HB value α 0 for the singlecell system is larger than the HB value α * for the coupled-cell system. This can also be seen in Figs. 6 and 8. For α slightly larger than α 0 (resp., α * ), the stable periodic solutions emerge in the single-cell systems (resp., coupled-cell systems), so that we can compute the average frequency (resp., collective frequency). For those values of α which deviate farther from α 0 (resp., α * ), the stable periodic solutions for the singlecell system (resp., coupled-cell systems) may persist, following the continuation of bifurcating periodic solutions from the Hopf bifurcation.
The following example illustrates a comparison between the average frequencies and the collective frequencies.
Example 4.1.2 For increasing values of α > α 0 , we compute the average frequency from the frequencies of individual cells in Examples 3.1, 3.2. For each of c * = 0.01, 0.05, 0.1 and its corresponding α * , with increasing values of α > α * , by solving the ODE system numerically, we observe the periodic solutions of coupledcell system (30) at these parameter values, and identify the frequency for each of these periodic solutions. The computed collective frequencies of these periodic solutions are plotted in Fig. 11. It can be seen that, in both HT and PS models, the frequency decreases as α increases, and it drops faster in the PS model. Notably, the leftmost points of the plots in Fig. 11 correspond to the average frequency about √ 3 at α near α 0 , and the collective frequencies at α near each α * , respectively. Moreover, for each fixed α, the collective frequency decreases as the coupling strength increases. These computations all show that the collective frequency is smaller than the average of

will constitute a solution of (29) if (M(t), P(t), R(t), V (t)) is a solution of
Accordingly, we can study the dynamics of system (29) restricted to the synchronous set S. That is, we can consider the N -cell analog of system (30) on its synchronous set, which is a four-dimensional system like (38). The characteristic polynomial for the linearization of the restricted system atx is exactly Δ + (λ). Recall that the bifurcation analysis in this section unfolds from the roots of Δ + (λ) = 0. Thus, Theorem 4 and the analysis herein remains valid for the restricted system. And the bifurcating periodic solution (x 1 (t), x 2 (t), x 3 (t), x 4 (t)) for the restricted system extends, with its copies (x 1 (t), x 2 (t), x 3 (t), x 4 (t)) N , to a synchronous periodic solution for the N -cell system.

Two Nonidentical Cells
In this section, we consider the coupled-cell system consisting of two nonidentical cells, by adding a scaling factor σ into the second cell: When c = 0, two cells are decoupled. We consider the parameter values at which each cell oscillates at its own frequency. As c increases from 0, the cells become coupled, yet it is not assured whether a collective periodic solution is then generated. The periodic solution, if exists, is no longer synchronous, and becomes phase-locked. Indeed, as system (39) comprises nonidentical subsystems, the synchronous set is no longer invariant and the solutions are asynchronous in general. Therefore, for different components of solutions to attain (or compromise to reach) the same period, it is quite natural that there exists a phase difference between the corresponding components of a collective periodic solution.
We again choose α as the bifurcation parameter and employ Theorem 1 to analyze the existence of periodic solution for system (39). Note that when c = 0, system (39) reduces to two decoupled subsystems and Components x 4 (t), y 4 (t) can be obtained by integrating the fourth and eighth components of system (39) once x 3 (t), y 3 (t) are solved. According to the discussions in Sect. 3, both of systems (40) and (41) undergo a Hopf bifurcation at α = α 0 , with respective frequency The positive equilibrium X = (x,ȳ) = (x 1 ,x 2 ,x 3 ,x 4 ,ȳ 1 ,ȳ 2 ,ȳ 3 ,ȳ 4 ) of system (39) satisfiesx 1 =x 2 =x 3 , f (x 3 ) =x 4 , andȳ 1 =ȳ 2 =ȳ 3 , f (ȳ 3 ) =ȳ 4 , and The existence of such heterogenous equilibrium is shown in Appendix E. The Jacobian matrix associated with the linearization of system (39) at X is We denote the characteristic polynomial Δ(λ) := det(λI − J (X; α)) by where the coefficients b i depend on s, σ, c, and γ,γ , which in turn depend onx 3 and y 3 . In contrast to the discussion on systems comprising two identical cells, herein, we could not factorize this polynomial Δ(λ). To apply the Hopf bifurcation theory, we look for the values of α * and the values for the other parameters, which satisfy the conditions of Theorem 1 with m = 8. We can also observe that α * → ∞ as γ → 0, for almost all parameter values. That is, the Hopf bifurcation does not occur at the flat parts of f 1 and f 2 . This is due to that α always appears together with γ orγ in J (X; α).
The following example is parallel to the discussions in Example 4.1.1. The computation process is similar, but slightly different, as we do not have an explicit expression of α * like (32) for the identical-cell case, nor a concrete formula like the one in Theorem 4 to resort to. We solve numerically for α * andx 3 andȳ 3 which satisfy D 7 (α * ) = 0 and (42). Then, we confirm that the other conditions of Theorem 1 are met. Again, for parameters of the individual cells, we take the same values as in Examples 3.1 and 3.2. For the scaling factor in the second cell, we take σ = 1.05. (39) with parameter values n = 11, k H = 0.136, s = 20, and σ = 1.05. The HB curve is plotted in Fig. 12. For illustration, we take one of the HB values P 1 (c * = 0.05, α * = 0.5165). We demonstrate the periodic solution at parameter values near P 1 in Fig. 13, which appears to be phase-locked (asynchronous). Variations of the eigenvalues along the line segment  Fig. 14a, c. There are two branches of complexconjugate eigenvalues λ M , λ M and λ m , λ m , and four negative real eigenvalues. We denote by λ M (P) and λ m (P) the eigenvalues of λ M and λ m at point P. When c = 0, λ M (P 0 ) and λ m (P 0 ) are purely imaginary, with λ m (P 0 ) = √ 3 i, λ M (P 0 ) = 1.05 √ 3 i. As c and α vary along P 0 P 1 , the λ M -branch moves to the left complex plane, and makes a turn downward to reach the imaginary axis at P 1 . Observe that λ M (P 1 ) = iω * c lies below λ M (P 0 ), i.e., the magnitude of λ M (P 1 ) is smaller than that of λ M (P 0 ). On the other hand, along P 0 P 1 , λ m -branch moves to and stays in the left complex plane. That is, at P 1 , there are a pair of purely imaginary eigenvalues λ M (P 1 ), λ M (P 1 ) and six eigenvalues with negative real parts, including λ m (P 1 ), namely ±1.801918i, −0.018114 ± 1.706965, −2.974896, −3.141083, −19.997793, −20.

Example 4.2 (i) Consider the coupled-cell HT model
Variations of the eigenvalues along the HB course are shown by the curves in Fig. 14b, c. As α increases from below α * to above α * , the stable equilibrium X becomes unstable and a stable periodic solution emerges, with frequency ω c close to ω * c ≈ 1.801918. (ii) Consider the coupled-cell PS model (39) with parameter values A = 0.0659, k d = 0.00001, s = 20, and σ = 1.05. We compute and draw the HB curve in Fig. 15. We take the HB value P 1 (c * = 0.038, α * = 0.5043). Variations of eigenvalues along the segment P 0 P 1 are shown by the curves in Fig. 16a, c. The scenario is similar to the one in Fig. 14a, c.
Variations of eigenvalues along the HB courses are shown in Fig. 16b, c. As α increases from below α * to above α * , the stable equilibrium X becomes unstable and a stable periodic solution emerges, with frequency ω c about ω * c ≈ 1.806613. The scenarios at the other Hopf bifurcation values (c * , α * ) in Fig. 12 (resp. Fig. 15) are qualitatively similar to the ones in Fig. 14 (resp. Fig. 16). For each of c * = 0.05, 0.075, 0.1 in the HT model and c * = 0.04, 0.05, 0.1 in the PS model, there corresponds an α * . We further trace the frequency of oscillation as α increases over each α * by solving system (39) numerically. The results are plotted in Fig. 18, where the average frequency is ω Ave = (ω 1 + ω 2 )/2, with ω 1 , ω 2 the frequencies of the first cell and second cell, respectively, when isolated. The leftmost points of the plots indicate the values α 0 for c = 0, and α * for each c * , which are different.
Remark 6 (i) We note that the HB curve is not connected to point P 0 in Figs. 12 and 15, since it is not known if the Hopf bifurcation can occur at arbitrarily small coupling strength c. (ii) We observe from Figs. 12 and 15 that the HB values α * for the HT model are larger than the ones for the PS model, at the same coupling strength c. This difference is also shown in Fig. 18, at the leftmost endpoints of these plots at coupling strength c = 0.05 and c = 0.1. Thus, for an α larger than α * for the HT model, it is even Variations of two complex eigenvalues as (c, α) moves along segment P 0 P 1 (red) and the HB course (blue) in Fig. 15: a λ M -branch along P 0 P 1 reaches the imaginary axis again at P 1 (0.038, 0.5043). b λ M -branch along the course, from solid square to hollow square, crosses the imaginary axis at P 1 . c The branches of λ m fall on the left complex plane along both of P 0 P 1 and the course. A supercritical Hopf bifurcation occurs when α * ≈ 0.5043. λ Ave (P 0 ) := [λ M (P 0 ) + λ m (P 0 )]/2 (Color figure online) larger than α * for the PS model. This could be a factor that the frequency for the PS model drops faster than the HT model as α increases over α * .

Collective Frequency Versus Average Frequency
With the discussions in Sects. 4.1 and 4.2, we like to observe further how the collective frequency of coupled-cell systems is compared to the average frequency of individual cells when isolated. For coupled-cell system (30) comprising two identical cells, the average frequency is the individual frequency. We have seen in Example 4.1.2 and Fig. 11 that the collective frequencies are lower than the average frequencies for both HT model and PS model. In addition, the collective frequencies are close to the average frequencies when the coupling strength c is small. This result also holds for coupled-cell system (29) with N identical cells, as mentioned in Remark 5.
With parameter values in Example 4.2 and three different coupling strengths c, it is indicated in Fig. 18 that at larger coupling strength c = 0.1, the collective frequency ω c is always smaller than the average frequency ω Ave , for α in a range where periodic solutions exist for both single-cell and coupled-cell HT model and PS model. In addition, the coupling strength c is smaller in the PS model than in the HT model, at which the collective frequency stays close to the average frequency in a range of α.
In the following examples, we fix α and allow coupling strength c to increase. If it happens that ω c is close to ω Ave , we are interested to see whether it holds only for certain parameter values and the coupling strength. If it holds only for one of the HT model and PS model, we like to see the scenario for the other model. Via solving the ODE system (39) numerically, we find the phase-locked periodic solutions for c starting from c = 0.035 in the coupled-cell HT model and from c = 0.037 in the coupled-cell PS model. Then, we observe the relation between the collective frequencies and the average of the individual frequencies. The results are shown in Fig. 19a. At the coupling strength c ≈ 0.037, the collective frequency ω c is less than and close to the average frequency ω Ave in the coupled-cell PS model, whereas the difference between ω c and ω Ave is substantial in the coupled-cell HT model. If we further increase the coupling strength c in the coupled-cell HT model, the collective frequency tends to the average frequency when c ≈ 0.075.

Example 4.3.2
Consider the same parameter values in Example 4.3.1, except that here we choose α = 0.56 in both models (larger than α 0 and all α * in Example 4.2) and vary the coupling strength c. We compute numerically and observe the phase-locked periodic solutions starting from c = 0.035 in both models. The results are shown in Fig. 19b. At c = 0.035, the collective frequency is greater than the average frequency in both models. For the coupled-cell PS model, the collective frequency matches the average frequency at coupling strength c ≈ 0.04. At this value of c, for the coupledcell HT model, the collective frequency is above the average frequency. If we still increase the coupling strength in the coupled-cell HT model, the collective frequency matches the average frequency when c ≈ 0.0776.
The above examples indicate that starting with similar oscillatory wave forms for individual cells in the HT model and in the PS model, the collective frequency of the coupled-cell system decreases as the coupling strength increases in both models. Moreover, the collective frequency drops faster in the PS model than in the HT model, as the coupling strength increases. In both models, there exist coupling strengths c such that the collective frequency equals the average frequency of individual cells. For the HT model, such strength c is larger. Our result in this regard is basically consistent with the one reported in Kim et al. (2014). It was reported in Kim et al. (2014) that proteinsequestration-based repression is more suitable for modeling circadian rhythms of mammals. One of the reasons is that the average frequency property holds for the PS model at suitable coupling strength, whereas this property holds for the HT model at coupling strength unreasonably large. Herein, we have proposed a methodology to examine closely the parameter values at which the average frequency property holds. Regarding which repression mechanism is more suitable for the modeling, among other concerns such as robustness, pertinence of parameter values in the models is certainly crucial.

Segmentation Clock Model
It is interesting to see the scenario of collective frequencies for other biological oscillators and compare it with the ones in Sect. 4. In this section, we perform parallel analysis and computation to a mathematical model on segmentation clock genes in zebrafish (Chen et al. 2018;Herrgen et al. 2012;Uriu et al. 2010). There are delay models and ODE models depicting the somitogenesis of zebrafish. We consider the following ODE system: System (44) depicts the interaction of two nonidentical cells if σ = 1, and two identical cells if σ = 1. Synchronous oscillations and traveling waves for the N -cell system, expanded from system (44) with σ = 1, were studied in Chen et al. (2018), , , Uriu et al. (2009Uriu et al. ( , 2010. All parameters are positive, and their meanings can be found in those papers. In particular, Hill coefficients n and h are positive integers. When the coupling strength ν c = 0, system (44) reduces to two decoupled subsystems, each of four dimension. We arrange the discussion of periodic solutions of single-cell systems in Supplementary Materials II. The Hopf bifurcation theory has been applied to investigate synchronous oscillations in system (44) with σ = 1, and in delay models in Chen et al. (2018), . It was shown that the collective frequency decreases as the coupling strength ν c increases. Herein, we investigate the periodic solution of system (44) which comprises two nonidentical subsystems, i.e., when σ = 1.
The existence of positive equilibrium X = (x,ȳ) = (x 1 ,x 2 ,x 3 ,x 4 ,ȳ 1 ,ȳ 2 ,ȳ 3 ,ȳ 4 ) for system (44) can be argued by the implicit function theorem, and is similar to the one in Appendix E. The Jacobian matrix associated with the linearization of system (44) at X can be computed, as shown in Supplementary Materials II.
The characteristic polynomial Δ(λ) := det(λI − J (X; ν 1 )) is where the coefficients b i can be computed. To apply the Hopf bifurcation theorem to system (44), we look for the value of ν 1 = ν * 1 and the values for the other parameters, which satisfy the conditions of Theorem 1 with m = 8. We follow the computation steps in Sect. 4.2.
We are interested in seeing how the collective frequency ω c of oscillation in system (44) varies with the coupling strength ν c and parameter ν 1 . We adopt the parameter values in Example S for single cell (in Supplementary Material II), where ν 1 denotes the HB value, and the frequency for the isolated (ν c = 0) individual cell is approximately ω 0 in Theorem 5 (in Supplementary Material II). The HB values are denoted by (ν * c , ν * 1 ).
Variation of eigenvalues along the HB course is depicted in Fig. 21a, b. We observe from numerical computations that periodic solutions emerge at ν 1 > ν * 1 for each ν * 1 , and ν 1 close to ν * 1 , with frequency about ω * c ≈ 0.322124. The scenarios at other Hopf bifurcation points (ν * c , ν * 1 ) in Fig. 20 are similar to the ones for P 1 in Fig. 21. For ν c = 0.006, 0.008, 0.01, respectively, we increase ν 1 further over ν * 1 , while holding ν c fixed at each ν * c , and perform numerical computation on system (44) to observe the periodic solutions and see how their frequencies vary with ν 1 . The result is shown in Fig. 22. The leftmost points of the plots represent the collective frequencies which are approximately ω * c , corresponding to ν 1 slightly larger than ν * 1 , respectively, and the average frequency about ω * Ave = (ω * 1 + ω * 2 )/2 ≈ (0.305785 + 1.05 · 0.305785)/2 corresponding to ν 1 slightly larger than ν 1 . It can be seen that, in each case, both the collective frequency and the average frequency decrease as ν 1 increases. Next, let us see how coupling strength ν c affects the collective frequency and compare it with the average of the individual frequencies.
Example 5.2 With parameter values in Example S, single-cell systems (54) and (55) undergo a Hopf bifurcation at ν 1 = ν 1 ≈ 0.074725, and coupled-cell system (44) with ν c = 0.006 undergoes a Hopf bifurcation at ν 1 = ν * 1 ≈ 0.06507. We take ν 1 = 0.076 so that the stable periodic solutions exist in both coupled-cell system (44) and singlecell systems (54), (55). Next, we increase the coupling strength from ν c = 0.006 and observe the relation between the collective frequency and the average frequency.

Fig. 22
Average frequency ω Ave = (ω 1 + ω 2 )/2 with respect to ν 1 for ν 1 > ν 1 , at ν c = 0, and collective frequency ω c with respect to ν 1 for ν 1 > ν * 1 , at ν c = 0.006, 0.008, 0.01, respectively The relations are depicted in Fig. 23. When ν c = 0.1, the collective frequency ω c ≈ 0.209346, i.e., the collective period T c ≈ 30.013387, falls within the range [25,35] minutes of biological interest. We observe that as the coupling strength increases, the collective frequency decreases and deviate farther from the average frequency of individual cells. This can also be seen in Fig. 22. The scenario is different from the one in coupled-cell HT and PS models, discussed in Sect. 4.
We remark that for some other parameter values, the periodic solutions may not persist for ν 1 extending well over ν * 1 in system (44). From the theory, the Hopf bifurcation occurs, and hence, the stable periodic solutions are confirmed to exist, at ν 1 slightly over ν * 1 in coupled-cell system (44), and ν 1 slightly over ν 1 in single-cell systems (54) and (55). In our computations, the values of ν * 1 are all smaller than those of ν 1 . If the periodic solution of system (44) does not exist for ν 1 at least larger than ν 1 , then we will not be able to compare the collective frequency and the average frequency at the same parameter values. Nevertheless, the present approach still leads us to locate the parameter values at which the existence of periodic solutions for the single-cell systems and the ones for the coupled-cell systems is assured.

Discussions and Conclusions
We have investigated the stable periodic solutions and their properties in some single-cell systems and coupled-cell systems. The methodology is based on the Hopf bifurcation theory and numerical simulations on the considered systems at parameter values unfolding from where the bifurcations occur. In addition, vanishing determinant for one of the Hurwitz matrices in the Routh-Hurwitz criterion leads to a condition which determines the parameter values where the Hopf bifurcation takes place. The condition not only detects and confirms the occurrence of Hopf bifurcation, but also allows us to analyze the bifurcation. With this approach, we further explored how the frequency of oscillation in the coupled-cell systems varies with respect to a system parameter and the coupling strength.
We applied this approach to investigate oscillations, represented by stable periodic solutions, and their frequencies, in a system modeling minimal genetic negative feedback loop. We analyzed the oscillatory properties and compared these properties between Hill-type repression and protein-sequestration-based repression. Taking α as the bifurcation parameter, we computed the bifurcation value α 0 for the singlecell systems. Then, we located the parameter values at which the oscillatory wave forms for the two systems, each with one of these two repressions, resemble each other. Then, for the coupled system comprising such single cells, we investigated how (at what parameter values) the collective periodic solutions emerge. We computed the eigenvalues of the linearized system along parameters (c, α) on the line segment from (0, α 0 ) to (c * , α * ), where α * is a HB value of the coupled-cell system at coupling strength c = c * . The purely imaginary eigenvalues provide the magnitudes of the collective frequencies at α = α * . This exhibits the relative magnitude between the average frequency of individual cells at α near α 0 and the collective frequency of coupled-cell system at α near α * . Unfolding from such information, we further computed to observe how the collective frequency of oscillation varies with the coupling strength and α. Moreover, we extended the computation to compare the average frequency with the collective frequency at the same value of α, and discussed the average frequency property. The collective behaviors for the coupled systems with two types of repression were compared. We observed that the collective frequency in the coupled-cell system with protein-sequestration-based repression approaches the average frequency at smaller coupling strength than the one with Hill-type repression.
For the system comprising two identical cells, we were able to express explicitly the HB values in terms of the component of the equilibrium, and the collective frequency at the HB values. This expression reveals the role played by the repression function on the frequency property. The influence from the Hill-type repression or from the protein-sequestration-based repression is indicated in the expression by the slope of the repression function. It can also be seen that the effect in the collective frequency at the HB values from this slope is enhanced by large coupling strength and suppressed by large timescale of the intercellular coupling. For further oscillatory properties at parameter values beyond the HB values, we resorted to numerical computations on the ODE systems. For the system comprising two nonidentical cells, an eight-dimensional ODE system, we still used the criterion in Liu (1994) to detect the Hopf bifurcation, but this process was carried out through numerical computations, as described in Sect. 4.2.
To compare with other biological oscillators, we performed a similar analysis and computation to a segmentation clock model. We observed different variations of eigenvalues along the parameter values connecting the HB value of the single-cell system and the HB value of the coupled-cell system. It appears in all our computations that the collective frequency of such coupled-cell system does not match the average frequency of individual cells.
Another significant issue is that our results enable a comparison between the oscillatory properties indicated in the kinetic models which take into account gene regulations and the ones obtained from the phase equations. As mentioned in Introduction, the average frequency property in the coupled phase Eq. (3) holds under the assumption of odd interaction function and symmetric connection. In our discussions in Sects. 4 and 5, the average frequency property only holds for specific parameter values or does not hold, even for small coupling strengths, see Figs. 18 and 22. Therefore, whether and under what assumption can the coupled phase equations accommodate the oscillatory properties of the coupled-cell systems remains an issue to be further examined.
In a sense of coordinating two individual oscillators with different frequencies, it was conceived that the coupling strength c has to attain a threshold (coupling strength is strong enough) to generate a collective periodic solution and the cells oscillate at a compromise frequency. The notion of compromise frequency was discussed in a coupled phase equation in Strogatz (1994). It was shown therein that the phaselocked solution exists only if the total coupling strength is larger than the difference of the individual frequencies, in a two-component phase equation. On the contrary, in Sect. 4.2, we saw that the HB curves in Figs. 12 and 15 can extend to small values of coupling strength c (as small as 10 −8 in AUTO computation). Our numerical computation shows that a supercritical Hopf bifurcation still occurs for c = 10 −4 , but the bifurcating periodic solution has very small amplitude. The difference of individual frequencies for these two decoupled cells is about 0.05 · √ 3 which is much larger than the coupling strength.
It has been our goal to develop an efficacious mathematical and computational approach to tackle the problems about how the oscillations of the individual cells with different intrinsic frequencies compromise to generate a collective periodic solution through intercellular coupling, and how the collective frequency of oscillation changes with the coupling strength and other parameters. The Hopf bifurcation theory provides a theoretical basis for confirming the existence of stable periodic solutions so that the findings on oscillations, their frequencies, and relevant dynamics are not merely based on numerical simulations. Combining effective numerical computation with the Hopf bifurcation theory, as demonstrated in this paper, one expands the capacity to observe the correspondence between oscillations and parameter values. Some of the examples in this work are based on the HB theory, whereas the other examples are extended from the HB theory via numerical simulations. Although the later are phenomenological findings, the present approach provides a way to explore the nature of single-cell systems and coupled-cell systems and make a close observation and comparison on the kinetics induced by Hill-type repression and the ones by protein-sequestration-based repression. The kinetics in terms of the properties imbedded in these mathematical models on biological processes can then be understood more thoroughly. We hope that the present approach has also contributed toward the selection or tuning of parameter values in the modeling.
C. The graphs of γ 1 and γ 2 : We compute to find γ 1 = γ 1 (x) = k n H nx n−1 (k n H +x n ) 2 . We drop the bar on x and regard γ 1 as a function of x. Then, There is a critical point x M := k H n−1 n+1 1 n , where γ 1 attains a local maximum, by the first derivative test. In addition, γ 1 (0) = 0, and lim x→∞ γ 1 (x) = 0.
We compute to find γ 2 = γ 2 (x) = Dropping the bar on x, then Thus, γ 2 is decreasing for x > 0, concave downward for 0 < x < A − k d and concave upward for x > A − k d . In addition, γ 2 (0) = 1/(A + k d ) and lim x→∞ γ 2 (x) = 0. The graphs of γ 1 and γ 2 for the data in Example 4.1.1 are drawn in Fig. 10. D. Numerics ofx, γ 1 , γ 2 , and ω * c , along the HB curves:   frequency ω * c at the HB values along the HB curve are plotted in Fig. 28. The graphs of γ 1 and γ 2 are drawn in Fig. 29. Note thatx now lies on the left hand side of x M .

Supplementary Materials
I. The terms determining the properties of Hopf bifurcation in Theorem 2: Let α denote one of the parameters in system (7). While holding the other parameters at fixed values, we write the system as dx/dt = f(x; α). After translation to the equilibriumx, we put it into the form:ẋ = J (x; α)x +f(x; α), where x = (x 1 , x 2 , . . . , x m ), J (x; α)x is the linear term, andf = ( f 1 , f 2 , . . . , f m ) is the nonlinear term. At α = α * , J (x; α * ) has a pair of purely imaginary eigenvalues λ(α * ) = ±iω * . Let us make a transformation so that the linear part is in normal form.
These values either coincide with the ones in the example in Uriu et al. (2010) or fall within the designated ranges therein. We look for a value ν 1 which satisfies (57). First, we express ν 1 in (57) in terms of γ 1 which in turn is a function ofx 3 . We substitute ν 1 = ν 1 into system (54) to find the equilibrium. After computation, the equilibrium isx = (x 1 ,x 2 ,x 3 ,x 4 ), withx 1 ≈ 0.026935,x 2 ≈ 0.629151,x 3 ≈ 0.541554,x 4 ≈ 1.207095. With thisx, we then compute to obtain the bifurcation value ν 1 ≈ 0.074726. Therefore, with further computation, we can confirm that a periodic solution emerges at ν 1 > ν 1 and ν 1 close to ν 1 , with frequency approximately ω 0 = 0.305785. We illustrate such result by computing numerically system (54) with ν 1 = 0.076 and the other parameter values as above. We observe the periodic solution with frequency 0.305142, shown in Fig. 30.