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 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 , whereWe note that average frequency is disparate from average period , and they were sometimes mixed up in the literature. The average period is given bywhere are the individual periods. For the average frequency in (1), the corresponding period iswhich is not equal to in (2) in general. Certainly, if all are close to each other, say for all , then the average frequency and the average period are about the same thing, asA mathematical model which well depicts the average frequency property is the coupled phase equations:where is the oscillatory phase of cell i, is the intrinsic frequency of the ith cell, 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 , provided that f is an odd function and 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 bywhere 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., , where is the half-saturation constant, and positive integer n is the Hill coefficient. A graph of 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. 1999, 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; Kurosawa et al. 2002; Kurosawa and Iwasa 2002, 2005).

Fig. 1

The transcription rate function a , the Hill-function-type repression with , , and b , the protein-sequestration-based repression with ,

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 bywhere 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 is depicted in Fig. 1b. As transcription is in a ratio to this fraction, 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 .

Goodwin’s model (system (4) with ) was proposed in Goodwin (1965). It was proved in Griffith (1968) that the equilibrium is stable for , by the Routh–Hurwitz criterion. When , it was shown that some parameter values can be found at which the equilibrium is unstable. In Kurosawa et al. Kurosawa et al. 2002, for a system similar to system (4), with , it was shown that the equilibrium is stable for , by the Routh–Hurwitz criterion, and for , 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 , 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 is explicitly revealed, yet, the crossing condition was not mentioned. In fact, is both a sufficient and necessary condition for the simple Hopf bifurcation. Our formulation considers general repression functions which accommodate both and .

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 , 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 coupled-cell 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 and the one modeled with .

Connection between the kinetic models and the phase models has been built in the so-called weakly coupled oscillators theory, cf. (Ermentrout and 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 , 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 of the second-to-last Hurwitz matrix 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 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 is too small ( 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, is a factor to the frequency at the HB value in the coupled-cell systems, and thus, the modeling with or 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 single-cell system (4) with and , 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 ODEswhere , , and are parameters. Let be an equilibrium of system (7) with a fixed . Denote the Jacobian matrix associated with the linearization of system (7) at by . Its characteristic polynomial isThe Hurwitz matrices associated with polynomial are defined asThe Routh–Hurwitz criterion indicates that all roots of (8) have negative real parts if and only if , , cf. (Gantmacher 1959; Hurwitz 1895; Kemperman 1982). Note that , and so the Routh–Hurwitz criterion can also be expressed by , , with .

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.

Lemma 1

has a pair of purely imaginary roots, and all other roots have negative real parts if and only if

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 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 , the determinant of Hurwitz matrices 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 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 equilibrium as a function of , i.e., . The Jacobian matrix associated with the linearization of (7) at then depends on , i.e., and is abbreviated as . The following Hopf bifurcation theorem was stated in Liu (1994), where instead of computing the eigenvalues of , one only needs to compute the determinants of the Hurwitz matrices.

Theorem 1

Consider system (7) whose equilibrium is a function of for near certain , when the other parameters are fixed. Assume that at , the following conditions holdThen, the system undergoes a Hopf bifurcation at , and a small-amplitude periodic solution surrounding 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 , the conditions of the theorem can be clearly expressed byTo 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 , and , and hence , , , and are introduced. Let be the branch of eigenvalue crossing the imaginary axis at , and , with . The following Hopf bifurcation theorem can be found in Hassard et al. (1981).

Theorem 2

Under the conditions of Theorem 1, the periods of the bifurcating periodic solutions of system (7) arewhereFurthermore,

the Hopf bifurcation is supercritical (resp., subcritical) and the periodic solution exists for (resp., ) with near , provided (resp., );

the periodic solution is stable (resp., unstable), provided (resp., );

the period T increases (resp., decreases) as increases, provided (resp., ) and , and decreases (resp., increases) as decreases, provided (resp., ) and .

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 , the Hill-type repression, and as PS model for , 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.

By assuming that the degradation rates are equaland scaling the variables and timeEquation (4) can be transformed into a non-dimensional system:where are derivatives with respect to , andHerein, we retain the same notations , since they have the same forms as (5) and (6). System of equations in form (15) has been studied in Kim (2016), Kim et al. (2014) through numerical simulations.

To analyze oscillatory properties for these models via Hopf bifurcation, herein, we consider another change of variables. Under the same assumption (14), we setThen, system (4) can be transformed intowhere is to serve as the bifurcation parameter.

We call system (16) with the single-cell HT model. In this system, the ratio of rates , the Hill coefficient n, and dissociation constant between the repressor and gene promoter determine the dynamics. We call system (16) with the single-cell PS model, where the rate ratio , activator concentration A, and dissociation constant between the activator and repressor, determine the dynamics. On the one hand, system (16) with 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 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 that is a positive equilibrium of system (16) if and only if , andSuch uniquely exists for or :

Proposition 1

There exists a unique positive equilibrium for system (16) with , for any , , and for system (16) with , for any A, , , where is the unique positive solution to , and , respectively. Moreover, for any fixed , , or fixed A, , is an increasing function of .

Proof

It is obvious that is an equilibrium for system (16) if and only if and satisfies . For , this reads asi.e., satisfies , where . For any , , there exists exactly one positive solution to this equation, due to that q is strictly increasing on , as , and . We also see that is a function of , as q is strictly increasing.

For , we see that for all , , , and for all . Thus, there exists exactly one positive solution to . As is strictly increasing and has an inverse, is thus an increasing function of .

Denote to indicate the dependence on . We note that from (17), can be expressed as a function of , i.e., . This is the inverse expression of . We shall analyze the periodic solutions bifurcating from . The Jacobian matrix associated with the linearization of system (16) at iswhere . The characteristic polynomial isOne 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 ifNotice that depends on which is a function of . Hence, whether the Hopf bifurcation takes place is determined by the existence of solution toIf such exists, then the Hopf bifurcation occurs at equilibrium and . Let us elaborate to confirm such existence for and , respectively. Since will play an important role in the analysis, we denote and .

For ,At , solving (21) with , we obtainThat is, for given , we set as (23). Then, with and computed from (22), equality (18) holds, and hence is an equilibrium of system (16) with .

For ,At , we solve (21) with , and obtainwhere the square root is nonnegative providedHerein, we considerand chooseso that is positive and the corresponding is not large.

Denote . For crossing condition (11), we havewhere if , and if . Note that is always positive for any smooth and strictly decreasing function f, including and . We confirm provided , and for all , under condition (25), at , in Appendix A, where is defined. That is, is both a sufficient and a necessary condition for the simple Hopf bifurcation of system (16) with .

According to Theorem 1, system (16) undergoes a Hopf bifurcation at and , and a small-amplitude periodic solution near emerges as or and is close to . Notice that at ,a constant matrix. The eigenvalues of this matrix are . Therefore, the bifurcating periodic solution has frequency near , for close to , for either or , by Theorem 2. We emphasize that for any repression function f in system (16), when the simple Hopf bifurcation takes place, the bifurcating periodic solution always has frequency about . Let us summarize:

Theorem 3

Assume that if , and (25) holds and if . System (16) undergoes a Hopf bifurcation at , where , and a small-amplitude periodic solution near emerges as or and is close to , with frequency about .

Remark 1

If we multiply each component in system (16) by a factor , or change the time from t to , then the eigenvalues of the linearized system at become . The system still undergoes a Hopf bifurcation at and , and the frequency of the bifurcating periodic solution is near . This property will be used in Sect. 4.2 for a cell-to-cell system which has two different individual frequencies.

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 at the HB value is a combination of the degradation rates (no longer ), but still does not depend on , and hence does not differentiate the form of repression function f. In addition, for , n will be required to be much larger than 8, if the difference among 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).

Example 3.1

Consider system (16) with and parameter values and . The graph of is depicted in Fig. 1a. According to Theorem 3, as , the Hopf bifurcation occurs at , where , computed from (23). Next, we compute to find from (22). Accordingly, a small-amplitude periodic solution emerges as the value of increases through , with frequency about . Furthermore, we compute to find that . The numerics are shown in Supplementary Material I. According to Theorem 2, system (16) with undergoes a supercritical Hopf bifurcation at and . 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 , plotted in Fig. 2. It can be seen that as increases from 0, the frequency decreases from about (i.e., the period is increasing), which is consistent with the assertion of Theorem 2.

Fig. 2

a Frequencies and b amplitudes of oscillations corresponding to various values of , where , , for the single-cell HT model in Example 3.1

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.

Example 3.2

Consider system (16) with , with parameter values and . The graph of is depicted in Fig. 1b. As condition (25) in Theorem 3 is met, the Hopf bifurcation occurs at , with computed from (26). We compute to find from (24), and thus . Therefore, a small-amplitude periodic solution emerges, as the value of passes through , with frequency about . Furthermore, we compute to find that . The numerics are revealed in Supplementary Material I. According to Theorem 2, system (16) with undergoes a supercritical Hopf bifurcation at and . The bifurcating periodic solution is stable and the period increases as increases. The frequencies and amplitudes corresponding to various values of are plotted in Fig. 3. When (), 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.

Fig. 3

a Frequencies and b amplitudes of oscillations corresponding to various values of , where , , for the single-cell PS model in Example 3.2

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 or , system (16) has a unique equilibrium . For each set of parameters with for the case , or satisfying (25) for the case , the Hopf bifurcation occurs at . The Jacobian matrix associated with the linearization of system (16) at with is a constant matrix. Restated, it is always this matrix whenever the Hopf bifurcation takes place in system (16), with or , or other smooth decreasing functions. Subsequently, the frequency of the bifurcating periodic solution near the bifurcation point is about . 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. 2 and 3.

The parameters in Hill-function repression are different from the ones in protein-sequestration-based repression. So how would one choose parameter values in each of these two models to make comparison? In Examples 3.1, 3.2, we have chosen the parameter values , in the single-cell HT model and , in the single-cell PS model. At , the oscillatory wave form in system (16) with is pretty close to the one with , as demonstrated in Fig. 4. This is a consequence of Theorem 3.

Fig. 4

Components , , and of the solutions of system (16) with in Example 3.1 (red), evolved from (0.1, 0.1, 0.1), and in Examples 3.2 (blue), evolved from (0.295, 0.295, 0.295), with . The periods are 3.6367989020 and 3.6367987711 for the HT model and the PS model, respectively (Color figure online)

Remark 2

Recall that system (15) is identical to system (16) with . In Kim et al. (2014), system (15) was considered with , in the HT model and , 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 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 .

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 , 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 (An et al. 2011). 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 is adopted.

When , 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 () and the coupled-cell PS model (). 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 () 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 bywhere parameter represents the coupling strength and is the timescale of intercellular coupling. The discussion herein is also valid for system (29) which comprises N identical cells, as remarked below. Denote , , and . The following proposition about the equilibrium is straightforward.

Proposition 2

For each set of parameter values, system (30) has a unique equilibrium which is homogeneous, i.e., , . For ,and is the unique positive solution to ; for ,and is the unique positive solution to . Moreover, for any fixed , or fixed , and , is a function of .

We show that system (30) has only homogeneous equilibrium 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 . The Jacobian matrix associated with the linearization of system (30) at iswhere . The characteristic polynomial can be factored aswhereThis factorization can be made since system (30) consists of two identical subsystems and the synchronous set is invariant, see Remark 5. Note that , , and . Considering (12) for and , respectively, we see that , whereThat is, if for some , then , and thus the simple Hopf bifurcation does not occur from . Accordingly, has a pair of purely imaginary roots and six roots with negative real parts if and only if , by Lemma 1. From this equality, we findwhere the square root is positive if and only ifHerein, we only consider the case , and hence only the “+” one in (32). Then, providedNote that impliesThus, (34) implies (33), provided . The case can also be discussed.

For crossing condition in (13), we considerat . When , (35) becomeswhere , which is 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 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 .

Theorem 4

Consider and assume that the equilibrium 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 emerges as or and is close to , with frequency aboutwhere if , and if .

Notably, it can be shown that if and only if . Certainly, the collective period of the bifurcating periodic solution is near . We will provide further observation from the assertion of this theorem below.

Implementation and Computation First, we fix the parameter values of n and in the HT model, and A and in the PS model. For fixed values of c and s, we substitute into (32), and express in terms of . Then we substitute into to solve for , and then obtain the equilibrium, as indicated in Proposition 2. Next, we substitute 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 of 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 and still call a HB value. These HB values form the Hopf bifurcation curve in the -plane.

According to Theorem 4, a small-amplitude periodic solution emerges, at or and close to , with frequency about . We are interested in seeing how the collective frequency 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 . Recall that the Hopf bifurcation occurs at in the single-cell systems.

Variations of eigenvalues At ), there are two pairs of purely imaginary eigenvalues for the linearized system at . For the identical-cell case, system (30), these two pairs coincide. It is interesting to see how the two complex-conjugate branches emanating from these two pairs move in complex plane , and which of them reaches the imaginary axis again at . As the magnitudes of the purely imaginary eigenvalues correspond to the frequencies of the individual cells at and coupled-cell at , we can observe the transition of frequency from the variation of eigenvalues from to . 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 with respect to the eigenvalues at . This reflects a transition from single-cell oscillation to coupled-cell oscillation. Moreover, we can also see how the local dynamics around is changing from the variations of eigenvalues. Certainly, there are infinitely many paths from to in -plane. For simplicity and illustration, we take the line segment from to 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 -plane by fixing , and allowing to increase from below to above . We call such segment HB course. We can compute the eigenvalues of to see if there is a stability switch for equilibrium 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 in Example 3.1, and A and in Example 3.2. Recall that the frequency for the isolated () individual cell is approximately , when is close to . In the following example, we fix . 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 , , and . With such parameter values, the single-cell system undergoes a Hopf bifurcation at , as shown in Example 3.1. Through computation, the HB curve is drawn in Fig. 6. For illustration, we take on the curve. At this HB value, we compute to find , and thus confirm that condition (34) in Theorem 4 holds. Note that the line segment lies above the HB curve, as shown in Fig. 6. We demonstrate the synchronous periodic solution at parameter values near in Fig. 5.

Fig. 6

For Example 4.1.1(i): a HB curve comprising HB values plotted by green curve on -plane; is a HB value for the single-cell system, is a HB value for the coupled-cell system. b Line segment in red line, the HB course through in blue line, from solid square to hollow square (Color figure online)

Fig. 5

Synchronous periodic solution in Example 4.1.1(i)

Variation of eigenvalues in along the segment is plotted by the curve in Fig. 7a, c. There are two branches of complex-conjugate eigenvalues and , and four negative real eigenvalues. We denote by and the eigenvalues and at point P. When , and are purely imaginary, with . As the parameters c and vary along , the -branch moves to the right complex plane and makes a turn downward to reach the imaginary axis at parameter value at . Notice that lies below , i.e., the collective frequency at is smaller than the individual frequency at . At the meantime, -branch moves to the left complex plane. To summarize, at , there are a pair of purely imaginary eigenvalues and six eigenvalues with negative real parts, including . They are , , , , , and .

Fig. 7

Variations of two complex eigenvalues as moves along segment (red) and the HB course (blue) in Fig. 6: a -branch along reaches the imaginary axis at , where . b -branch along the course, from solid square to hollow square, crosses the imaginary axis at parameter value . c -branches fall on the left complex plane along both and the course. A supercritical Hopf bifurcation occurs at (Color figure online)

Variation of eigenvalues in 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 close to , by Theorem 4. Our numerical simulation shows that for slightly above , the collective frequency is for .

(ii) Consider coupled-cell PS model (30) with parameter values , , and . With such parameter values, the single-cell system undergoes a Hopf bifurcation at , as shown in Example 3.2. We compute and draw the HB curve in Fig. 8. For illustration, we take in the curve. We compute to find , and then confirm the condition (34) in Theorem 4. Note that this line segment lies below the HB curve, as shown in Fig. 8. Variations of eigenvalues along are depicted by the curves in Fig. 9a, c. When , and are purely imaginary, with . As the parameters c and vary along , -branch moves to the left complex plane, and makes a turn downward to reach the imaginary axis at parameter value . As lies below , we see that the collective frequency at is smaller than the individual frequency at . Meanwhile, -branch moves to the left complex plane. At , the eigenvalues are computed as , , , , , and .

Fig. 8

For Example 4.1.1(ii): a HB curve comprising HB values plotted by green curve on -plane; is a HB value for the single-cell system, is a HB value for the coupled-cell system. b Line segment in red line; the HB course through in blue line, from solid square to hollow square (Color figure online)

Fig. 9

Variations of two complex eigenvalues as moves along segment (red) and the course (blue) in Fig. 8: a -branch along reaches the imaginary axis at parameter value . b -branch along the course, from solid square to hollow square, crosses the imaginary axis at . c -branches fall on the left complex plane along both and the course. A supercritical Hopf bifurcation occurs at (Color figure online)

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 close to . Our numerical simulation reveals that for slightly above , the collective frequency is at .

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 , resemble the ones for in Figs. 7a, c and 9a, c. The variation of eigenvalues along each associated HB course is also similar to the one for 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, -branch moves to the right half plane of and then makes a turn downward to reach the imaginary axis at , whereas -branch remains on the left complex plane. In the PS model, both -branch and -branch move into the left complex plane. Then, -branch makes a turn downward and reaches the imaginary axis at , and -branch continues to stay in the left complex plane. The basic reason for the different movement of -branch is that line segment 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 between these two models. For the HT model, the dimension of the unstable manifold of increases from zero to two (the dimension of the stable manifold from four to six) when moves from along , before reaching . On the other hand, for the PS model, the dimension of the stable manifold of increases from four to eight when moves from , before reaching . 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 is smaller than the individual frequencies at in both models. Indeed, the purely imaginary eigenvalues lie below the purely imaginary eigenvalues in both models.

Remark 3

The intercellular coupling is relatively fast, and we took in the above examples. We observe from system (30) that and may be regarded as being in quasi-steady state, i.e., , , and thus the first equation is approximatelyRecalling the Hopf bifurcation in the single-cell system (16), we see that the occurrence of the Hopf bifurcation in the coupled system (30) is approximately along the curve where in Example 4.1.1(i), and in Example 4.1.1(ii). This explains why the HB curves in Figs. 6 and 8 are almost straight 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

It is seen from (32) that the HB value for the coupled-cell system is smaller than the HB value for the single-cell system (when ). 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.

We observe from (32) that as . That is, the Hopf bifurcation does not occur at the flat parts of the repression functions and (small values of and ), see Fig. 10. On the other hand, the Hopf bifurcation occurs when , the magnitude of , 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 and , and thus, condition (34) is automatically met.

Fig. 10

The graphs of and , for the parameter values in Example 4.1.1, with , , , . The green diamond represents , the magenta diamond represents when in the HT model, and the magenta circle represents when in the PS model (Color figure online)

From Theorem  4, we see that frequency decreases as increases. That is, in contrast to the single-cell case, the frequency of the periodic solution near the HB value distinguishes between and , and hence the repression function from . From the graphs of and in Fig. 10, we see that while has larger maximum value, 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 on which is determined by the parameter values. For the HT model, we may compare the point where attains its maximum value with the component of the equilibrium at , that is, Their distance can be measured from When the coupling strength c is small, in the coupled system is close to in (36), and the quantity in (37) estimates how far the point that takes on is away from , where the maximum value of is attained. Note that for , lies on the right hand side of , whereas it is reverse for . In addition, we see that the gap between and can be enlarged, by choosing larger and n. The corresponding frequencies will then further distinguish the repression function from . The analysis for the graphs of and is arranged in Appendix C. From the expression for frequency , 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 of , , and the corresponding frequencies , along the HB curves for the cases with , and .

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 and the average frequency 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 for the single-cell 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 (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 (resp., ), the stable periodic solutions for the single-cell 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 , we compute the average frequency from the frequencies of individual cells in Examples 3.1, 3.2. For each of and its corresponding , with increasing values of , by solving the ODE system numerically, we observe the periodic solutions of coupled-cell 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 at near , 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 individual frequencies (when ), and are close to the average frequency when c is small.

Fig. 11

Collective frequencies of the periodic solutions for the coupled-cell systems and average frequencies for the periodic solutions of the single-cell systems, corresponding to increasing value of , at , in a the HT model, b the PS model, with the values of the other parameters given in Examples 4.1.1 and 4.1.2, respectively

Remark 5

In system (29), the synchronous set is invariant. Indeed, , , will constitute a solution of (29) if (M(t), P(t), R(t), V(t)) is a solution ofAccordingly, we can study the dynamics of system (29) restricted to the synchronous set . 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 at is exactly . Recall that the bifurcation analysis in this section unfolds from the roots of . Thus, Theorem 4 and the analysis herein remains valid for the restricted system. And the bifurcating periodic solution for the restricted system extends, with its copies , 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 , 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 , system (39) reduces to two decoupled subsystemsandComponents can be obtained by integrating the fourth and eighth components of system (39) once are solved. According to the discussions in Sect. 3, both of systems (40) and (41) undergo a Hopf bifurcation at , with respective frequencyThe positive equilibrium of system (39) satisfies , , and , , andThe existence of such heterogenous equilibrium is shown in Appendix E. The Jacobian matrix associated with the linearization of system (39) at iswhere , and . We denote the characteristic polynomial bywhere the coefficients depend on , and , which in turn depend on and . 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 . We can also observe that as , for almost all parameter values. That is, the Hopf bifurcation does not occur at the flat parts of and . This is due to that always appears together with or in .

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 and and which satisfy 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 .

Example 4.2

(i) Consider the coupled-cell HT model (39) with parameter values , , , and . The HB curve is plotted in Fig. 12. For illustration, we take one of the HB values . We demonstrate the periodic solution at parameter values near in Fig. 13, which appears to be phase-locked (asynchronous). Variations of the eigenvalues along the line segment are shown by the curves in Fig. 14a, c. There are two branches of complex-conjugate eigenvalues and , and four negative real eigenvalues. We denote by and the eigenvalues of and at point P. When , and are purely imaginary, with , . As c and vary along , the -branch moves to the left complex plane, and makes a turn downward to reach the imaginary axis at . Observe that lies below , i.e., the magnitude of is smaller than that of . On the other hand, along , -branch moves to and stays in the left complex plane. That is, at , there are a pair of purely imaginary eigenvalues and six eigenvalues with negative real parts, including , namely , , , , , .

Fig. 12

For Example 4.2.1(i): The HB curve plotted in green curve on -plane; is a HB value for the single-cell system. is a HB value for the coupled-cell system; line segment in red line, the HB course in blue line (Color figure online)

Fig. 13

Phase-locked periodic solution, with phase difference approximately 0.310757 (the difference of t between the peaks of and ), in Example 4.2(i)

Fig. 14

Variations of two complex eigenvalues as moves along segment (red) and the HB course (blue) in Fig. 12: a Along , -branch reaches the imaginary axis again at . b Along the course, from solid square to hollow square, -branch crosses the imaginary axis at . (c) -branches fall on the left complex plane along both of and the course. A supercritical Hopf bifurcation occurs when . (Color figure online)

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 becomes unstable and a stable periodic solution emerges, with frequency close to .

(ii) Consider the coupled-cell PS model (39) with parameter values , , , and . We compute and draw the HB curve in Fig. 15. We take the HB value . Variations of eigenvalues along the segment are shown by the curves in Fig. 16a, c. The scenario is similar to the one in Fig. 14a, c.

Fig. 15

For Example 4.2(ii): The HB curve plotted by green curve on -plane; is a HB value for the single-cell system, is a HB value for the coupled-cell system; in red line, the HB course in blue line (Color figure online)

Fig. 16

Variations of two complex eigenvalues as moves along segment (red) and the HB course (blue) in Fig. 15: a -branch along reaches the imaginary axis again at . b -branch along the course, from solid square to hollow square, crosses the imaginary axis at . c The branches of fall on the left complex plane along both of and the course. A supercritical Hopf bifurcation occurs when . (Color figure online)

Variations of eigenvalues along the HB courses are shown in Fig. 16b, c. As increases from below to above , the stable equilibrium becomes unstable and a stable periodic solution emerges, with frequency about .

The scenarios at the other Hopf bifurcation values in Fig. 12 (resp. Fig. 15) are qualitatively similar to the ones in Fig. 14 (resp. Fig. 16). For each of in the HT model and 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 , with , the frequencies of the first cell and second cell, respectively, when isolated. The leftmost points of the plots indicate the values for , and for each , which are different.

Fig. 18

Average frequency with respect to at , and collective frequency with respect to at fixed c in a the HT model with , and b the PS model with , where and are the frequencies for the two single cells corresponding to the value of

Remark 6

We note that the HB curve is not connected to point in Figs. 12 and 15, since it is not known if the Hopf bifurcation can occur at arbitrarily small coupling strength c.

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 and . Thus, for an larger than for the HT model, it is even 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.

For coupled system (39) consisting of two nonidentical cells with individual frequencies and when isolated, the average frequency is . To indicate dependence on , we denote , and , where is the bifurcation value for the single cell. Note that . In Examples 4.2(i) (HT model) and 4.2(ii) (PS model), we have seen that both lie below , with , and , . That is, when the collective periodic solution is formed at a HB value , the collective frequency drops below the larger individual frequency at (decoupled) and near . However, on the one hand, the average frequency at and near is about . On the other hand, the comparison of the collective frequency (at ) and the average frequency (at ) should be made at the same parameter values, including . Notably, our notation means , .

Our strategy of comparison starts with comparing , , with , , and , where is a HB value. That is, we compare the eigenvalues at with the ones at . For example, for the case , if lies below on the imaginary axis of , i.e., , and if decreases as increases (with fixed ), then we are sure that at slightly larger than (with fixed ), should the periodic solution exist up to such . An example for such a scenario is the HB value in Example 4.2(ii) (not shown in Fig. 15); the movement of is shown in Fig. 17. If lies above , i.e., , we can still compute numerically for increasing values of , starting from , and compare it with as passes over .

Fig. 17

Variations of complex eigenvalue a along . b along the associated HB course. c along (dotted curve) and the course (solid curve) in Fig. 15, where is at ; ,

With parameter values in Example 4.2 and three different coupling strengths c, it is indicated in Fig. 18 that at larger coupling strength , the collective frequency is always smaller than the average frequency , 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 is close to , 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.

Example 4.3.1

For parameters of the individual cells, we take , as in Examples 3.1, 4.1.1(i), 4.2(i) (HT model), and , as in Examples 3.2, 4.1.1(ii), 4.2(ii) (PS model). In addition, we take the same time scaling factor and as in the above examples. Recall that Hopf bifurcation occurs at in Example 3.1 (single-cell HT model), and at in Example 3.2 (single-cell PS model). We fix so that the single-cell HT model and the single-cell PS model have analogous oscillatory wave forms, as in Fig. 4.

Via solving the ODE system (39) numerically, we find the phase-locked periodic solutions for c starting from in the coupled-cell HT model and from 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 , the collective frequency is less than and close to the average frequency in the coupled-cell PS model, whereas the difference between and 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 .

Fig. 19

Relation between the collective frequency and the average of individual frequencies for a range of coupling strength c, when a , the average frequency is 1.77086089 for , and 1.77086095 for , b , the average frequency is 1.774609880 for , and 1.774625695 for

In Example 4.3.1, may be too far from the values of in Examples 4.2.1, 4.2.2. Let us consider the value of which is closer to .

Example 4.3.2

Consider the same parameter values in Example 4.3.1, except that here we choose in both models (larger than and all in Example 4.2) and vary the coupling strength c. We compute numerically and observe the phase-locked periodic solutions starting from in both models. The results are shown in Fig. 19b. At , 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 . At this value of c, for the coupled-cell 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 .

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 protein-sequestration-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 , and two identical cells if . Synchronous oscillations and traveling waves for the N-cell system, expanded from system (44) with , were studied in Chen et al. (2018), Liao and Shih (2012), Liao et al. (2012), Uriu et al. (2009, 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 , 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 , and in delay models in Chen et al. (2018), Liao and Shih (2012). It was shown that the collective frequency decreases as the coupling strength increases. Herein, we investigate the periodic solution of system (44) which comprises two nonidentical subsystems, i.e., when .

The existence of positive equilibrium 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 can be computed, as shown in Supplementary Materials II.

The characteristic polynomial iswhere the coefficients can be computed. To apply the Hopf bifurcation theorem to system (44), we look for the value of and the values for the other parameters, which satisfy the conditions of Theorem 1 with . We follow the computation steps in Sect. 4.2.

We are interested in seeing how the collective frequency of oscillation in system (44) varies with the coupling strength and parameter . We adopt the parameter values in Example S for single cell (in Supplementary Material II), where denotes the HB value, and the frequency for the isolated () individual cell is approximately in Theorem 5 (in Supplementary Material II). The HB values are denoted by .

Example 5.1

Consider coupled system (44) with parameter values in Example S and . When , the system is decoupled, and each of the two subsystems has a periodic solution with respective frequency about and . The HB curve consisting of the HB values is drawn in Fig. 20. For illustration, we take the line segment , with , , where and . The variation of eigenvalues along is depicted in Fig. 21a, c. There are two branches of complex-conjugate eigenvalues and , and four negative real eigenvalues. We denote by and the eigenvalues of and at point P. When , and are purely imaginary, with . As the parameters and vary along , -branch moves to the left complex plane, and makes a turn upward to reach the imaginary axis at . That is, lies above . On the other hand, along the route, -branch moves to and stays in the left complex plane. To summarize, at , there is a pair of purely imaginary eigenvalues and the remaining eigenvalues have negative real parts, including . These eigenvalues are , , , , , and .

Fig. 20

For Example 5.1: The HB curve plotted by green curve on -plane; is a HB value for the single-cell system, is a HB value for the coupled-cell system, the segment in red line, the HB course in blue line (Color figure online)

Fig. 21

Variations of two complex eigenvalues as moves along segment (red), and the HB course (blue) in Fig. 20: a -branch reaches the imaginary axis at ; -branch along the course, from solid square to hollow square, crosses the imaginary axis at parameter value . b -branches enter into and stay in the left complex plane along and the course. A supercritical Hopf bifurcation occurs at . (Color figure online)

Variation of eigenvalues along the HB course is depicted in Fig. 21a, b. We observe from numerical computations that periodic solutions emerge at for each , and close to , with frequency about . The scenarios at other Hopf bifurcation points in Fig. 20 are similar to the ones for in Fig. 21.

For , respectively, we increase further over , while holding fixed at each , and perform numerical computation on system (44) to observe the periodic solutions and see how their frequencies vary with . The result is shown in Fig. 22. The leftmost points of the plots represent the collective frequencies which are approximately , corresponding to slightly larger than , respectively, and the average frequency about corresponding to slightly larger than . It can be seen that, in each case, both the collective frequency and the average frequency decrease as increases.

Fig. 22

Average frequency with respect to for , at , and collective frequency with respect to for , at , respectively

Next, let us see how coupling strength 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 , and coupled-cell system (44) with undergoes a Hopf bifurcation at . We take so that the stable periodic solutions exist in both coupled-cell system (44) and single-cell systems (54), (55). Next, we increase the coupling strength from and observe the relation between the collective frequency and the average frequency. The relations are depicted in Fig. 23. When , the collective frequency , i.e., the collective period , 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.

Fig. 23

Collective frequency versus the average of individual frequencies, for increasing values of , with

We remark that for some other parameter values, the periodic solutions may not persist for extending well over in system (44). From the theory, the Hopf bifurcation occurs, and hence, the stable periodic solutions are confirmed to exist, at slightly over in coupled-cell system (44), and slightly over in single-cell systems (54) and (55). In our computations, the values of are all smaller than those of . If the periodic solution of system (44) does not exist for at least larger than , 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 for the single-cell 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 on the line segment from to , where is a HB value of the coupled-cell system at coupling strength . 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 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 phase-locked 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 in AUTO computation). Our numerical computation shows that a supercritical Hopf bifurcation still occurs for , but the bifurcating periodic solution has very small amplitude. The difference of individual frequencies for these two decoupled cells is about 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.