Circadian rhythms in mammals are controlled by the neurons located in the suprachiasmatic nucleus of the hypothalamus. In physiological conditions, the system of neurons is very efficiently entrained by the 24 h light–dark cycle. Most of the studies carried out so far emphasize the crucial role of the periodicity imposed by the light–dark cycle in neuronal synchronization. Nevertheless, heterogeneity as a natural and permanent ingredient of these cellular interactions seemingly plays a major role in these biochemical processes. In this paper, we use a model that considers the neurons of the suprachiasmatic nucleus as chemically coupled modified Goodwin oscillators, and introduce non-negligible heterogeneity in the periods of all neurons in the form of quenched noise. The system response to the light–dark cycle periodicity is studied as a function of the interneuronal coupling strength, external forcing amplitude and neuronal heterogeneity. Our results indicate that the right amount of heterogeneity helps the extended system to respond globally in a more coherent way to the external forcing. Our proposed mechanism for neuronal synchronization under external periodic forcing is based on heterogeneity-induced oscillator death, damped oscillators being more entrainable by the external forcing than the self-oscillating neurons with different periods.
Circadian rhythms are light–dark-dependent cycles of roughly 24 h present in the biochemical and physiological processes of many living entities . In mammals, the main mediator between the light–dark periodicity and the biological rhythms is formed by two interconnected suprachiasmatic nuclei (SCN), located in the hypothalamus. These nuclei form the so-called ‘circadian pacemaker’ and contain about 10 000 neurons each [1,2].
The main property of the SCN is that their activity displays self-sustained oscillations in synchrony with the external forcing imposed by the light–dark cycle. The exact mechanism leading to this behaviour has been the subject of intense research. It has been shown that, when taken individually, neurons produce oscillations with a constant period ranging from 20 to 28 h [3,4]. The oscillatory behaviour originates in a regulatory circuit with a negative feedback loop. The relevant question is how this individual oscillatory behaviour translates into common, global, oscillations of the SCN activity synchronized with the external light stimulus.
It has been shown that the origin of the oscillatory activity of the circadian pacemaker at the global level resides on the interaction between the SCN neurons. Coupling between cells in the SCN is achieved partly by neurotransmitters [3,5] and it is by means of those neurotransmitters that external forcing by light influences the neuronal synchronization. For example, the vasoactive intestinal polypeptide (VIP) has been shown to be necessary in mediating both the periodicity and the internal synchrony of mammalian clock neurons [6–8]. Therefore, a model of coupled and forced neurons appears quite naturally as responsible for the circadian rhythms. Along these lines, an interesting mechanism has been put forward recently by Gonze et al.  and by Bernard et al. . They proposed that synchronization to the external forcing is facilitated by the fact that interneuronal coupling transforms SCN into damped oscillators that can then be easily entrained.
In this paper, we show that the presence of some level of heterogeneity or dispersion in the intrinsic periods of the oscillators [11,12] can improve the response of the coupled neuronal system to the external light–dark forcing. The proposed mechanism for the improvement of the neuronal synchronization under external periodic forcing bears some similarities to the one proposed in Gonze et al.  and Bernard et al.  in the sense that the oscillators are brought to a regime of oscillator death [13,14], but in our case, this regime is induced by the presence of heterogeneity. Once this regime has been reached, the damped oscillators are more entrainable by the external forcing than the self-oscillating neurons with different periods, or the synchronized oscillatory state that appears in the strong coupling regime but with a period larger than the individual neuronal periods.
To be more specific, we will assume that the periods of the individual neurons are random variables drawn from a normal distribution. We will then analyse the global response of the system to the light–dark cycle periodicity as a function of the interneuronal coupling strength, external forcing amplitude and neuronal heterogeneity. We show that the presence of the right amount of dispersion in the periods of the neurons can indeed enhance the synchronization to the external forcing.
Period dispersion arises as a consequence of the cellular heterogeneity at the biochemical level, which is an experimentally well-observed fact [3,15]. It can act in either physiological or pathological conditions. An example of the latter is the diversification of antigenic baggage present in tumour cells that makes them more difficult to be recognized and captured by the defence mechanisms and therefore more prone to migrate and develop metastasis . Our results show that some level of disorder can be of help when synchronizing neuronal activity to the external forcing. Although counterintuitive, it has been unambiguously shown that the addition of various forms of disorder can improve the order in the output of a large variety of nonlinear systems. For example, the mechanism of stochastic resonance [17,18] shows that the response of a bistable system to a weak signal can be optimally amplified by the presence of an intermediate level of dynamical noise. Stochastic resonance is not a rare phenomenon; it has been repeatedly shown to be relevant in physical and biological systems described by nonlinear dynamical equations [17,18]. In large systems with many coupled elements, noise is responsible for a large variety of ordering effects, such as pattern formation, phase transitions, phase separation, spatio-temporal stochastic resonance, noise-sustained structures, doubly stochastic resonance, among many others . All these examples have in common that some sort of order at the macroscopic level appears only in the presence of the right amount of noise or disorder at the microscopic level. Furthermore, it has been proved that noise may play a constructive role in nonlinear systems, by enhancing coherent (periodic) behaviour near bifurcations and phase transitions [20,21]. In this paper, we introduce non-negligible random heterogeneity into the periods of all neurons, the so-called quenched noise. Numerical simulations suggest (data not shown) that the results are valid as well when the quenched noise is introduced into the model parameters. A different approach is the consideration of intracellular stochastic variability owing to low molecule numbers  or both variability and heterogeneity.
Close to our work is the study by Ueda et al. , where the effect of fluctuations in neuron parameter values is assessed and it is shown that the coupled system is relatively robust to noise. Previous theoretical studies have addressed the effect of noise on genetic oscillators [24–26], and some have proposed an ordering influence of noise on circadian clocks at the single cell level in cases where neither light intensity nor coupling strength by themselves can synchronize the system. Collective phenomena induced by heterogeneity in autonomous, non-forced systems have also been discussed in the literature. For example, de Vries & Sherman  and Cartwright  have shown that collective bursting or firing can appear in excitable systems and a general theory of the role of heterogeneity in those systems has been developed by Tessone et al. . In this paper, we refer to the collective response in systems of nonlinear oscillators subjected to the action of an external forcing representing the dark–light cycle.
The paper is organized as follows. In §2, we will describe in detail the model of circadian oscillators and the methods we use. It is a coupled extension of the original Goodwin oscillator  as developed by Gonze et al. . In §3, we analyse the system response to the periodic external forcing, as a function of the external forcing amplitude, coupling strength and neuronal diversity or heterogeneity. By simulating numerically the governing differential equations, we identify the range of these parameters for which the extended system oscillates in synchrony and entrained to the external light period. Section 3.1 describes the mechanism through which the neuronal heterogeneity favours the synchronization with the external forcing and analyses the combined influence of the coupling strength, neuronal heterogeneity and light amplitude on the stability of the linearized system of coupled oscillators. We show that a mean variable in this model exhibits a transition from a rhythmic to an arrhythmic dynamics (the so-called oscillator death [13,14]). Concluding remarks are found in §4.
2. Model and methods
2.1. The circadian pacemaker
As stated in the introduction, our aim is to consider the role that the heterogeneity in the population of neurons plays in the global response of the SCN to an external oscillating stimulus. To this end, we consider an ensemble of coupled neurons subject to a periodic forcing. Each of the neurons, when uncoupled from the others and from the external stimulus, acts as an oscillator with an intrinsic period. Heterogeneity is considered insofar as the individual periods are not identical, but show some degree of dispersion around a mean value. For each one of the neurons in the SCN we use a four-variable model proposed by Gonze et al. , which is based originally on the Goodwin oscillator , to describe circadian oscillations in single cells. The variables of the model are as follows. The clock gene mRNA (X) produces a clock protein (Y), which activates a transcriptional inhibitor (Z) and this in turn inhibits the transcription of the clock gene, closing a negative feedback loop. The mRNA X also excites the production of neurotransmitter V, which in the coupled system will be then responsible for an additional positive feedback loop. In order to overcome the high Hill coefficients required for self-oscillations, Gonze et al. replaced the linear degradation by nonlinear Michaelis–Menten terms. This leads to the system of equations 2.1 2.2 2.3 2.4 which, depending on parameters, might produce oscillations in a stable limit cycle. Using the values ν1 = 0.7 nM h−1, ν2 = ν4 = ν6 = 0.35 nM h−1, ν8 = 1 nM h−1, K1 = K2 = K4 = K6 = K8 = 1 nM, k3 = k5 = 0.7 h−1, k7 = 0.35 h−1, the period of the limit cycle oscillations is T = 23.5 h.
For the complete model, we take N neuronal oscillators, each one of them described by four variables (Xi, Yi, Zi, Vi), i = 1, … , N, satisfying the above evolution equations. Heterogeneity in the intrinsic periods is introduced by multiplying the left-hand side of each one of the equations (2.1)–(2.4) by a scale factor τi. Hence, the intrinsic period Ti of the isolated neuron i is τiT. The numbers τi are independently taken from a normal random distribution of mean 1 and standard deviation σ. Since the periods must be positive, in the numerical simulations, we have explicitly checked that, for the values of σ considered later, τi never takes a negative value, which would be unacceptable. The standard deviation σ will be taken as a measure of the diversity. A value of σ = 0.1 for example corresponds to a standard deviation of 10 per cent in the individual periods of the uncoupled neurons, close to the observed variation of periods between 20 and 28 h.
Two additional factors influence the dynamics of single cell oscillations: forcing by light and intercellular coupling. Both are assumed to act independently from the negative feedback loop and are added as independent terms in the transcription rate of X . Light is incorporated through a periodic time-dependent function L(t), which can be realized in various forms. In the majority of the presented results, and unless stated otherwise, we have used a sinusoidal signal, L(t) = (L0/2)(1 + sin ωt). In some cases, for comparison and to simulate different day lengths, we have used a square wave In both ways, the signal oscillates between the values L(t) = 0 and L(t) = L0 with a period 2π/ω = 24 h.
Coupling between the neurons is assumed to depend on the concentration F of the synchronizing factor (the neurotransmitter) in the extracellular medium, which builds up by contributions from all neurons. Under a fast transmission hypothesis, the extracellular concentration is assumed to equilibrate to the average, mean-field, cellular neurotransmitter concentration, F = (1/N)∑i=1N Vi. The resulting model is 2.5 2.6 2.7 2.8 2.9 with νc = 0.4 nM h−1, Kc = 1 nM.
There is experimental evidence supporting the assumption of a chemical (rather than electrical) mechanism of intercell communication among SCN neurons as a synchronization factor and, in fact, mechanisms other than neurotransmitters or electrical coupling for the SCN communication have been suggested (e.g. by Pol & Dudek ). Furthermore, more realistic modelling that takes into account all variables known to participate of the negative feedback loop has been introduced. These models may include up to 10 variables and corresponding equations for each single cell .
It seems, however, that in order to get an understanding of the SCN dynamics, a sufficient tool is the four-variable model described above. In fact, the synchronization of damped oscillators is independent from the particular intracellular model used and as discussed by Bernard et al. , this system, the model developed by Leloup & Goldbeter , and other simple negative feedback oscillators have similar synchronization properties. In this paper, we have decided to use the simpler four-variable model although most of our results are also valid in the more complex 10-variable model.
A model close to equations (2.5)–(2.9) has been used by Ullner et al. , where the authors investigate how the interplay between fluctuations of constant light and intercellular coupling affects the dynamics of the collective rhythm in a large ensemble of non-identical, globally coupled oscillators. In their case, however, an inverse dependence of the cell–cell coupling strength on the light intensity was implemented, in such a way that the larger the light intensity the weaker the coupling.
2.2. Measures of synchrony and entrainment
Owing to the effect both of coupling and of forcing, the neurons might synchronize their oscillations. There are several possible measures of how good this synchronization is. In this paper, the interneuronal synchronization will be quantified by the parameter of synchrony ρ, defined as 2.10 where 〈⋯〉 denotes a time average in the long-time asymptotic state. The parameter ρ varies between a value close to 0 (no synchronization) and 1 (perfect synchronization, with all neurons in phase, Vi(t) = Vj(t), ∀i,j). It is important to note that even if the neurons synchronize perfectly their oscillations, the period of those oscillations does not necessarily coincide with the mean period T of the individual oscillators or with the period 2π/ω of the external forcing. In fact, in the unforced (no light) case, the period of the common oscillations (for the set of parameters given before and a dispersion of σ = 0.05 and coupling K = 0.5) is approximately equal to 26.5 h whereas the period of the forcing is 2π/ω = 24 h and the mean period of the individual uncoupled oscillators is T = 23.5 h .
Besides the previous measure of synchronization among the oscillators, we are also concerned about the quality of the global response of the neuronal ensemble to the external forcing L(t). A suitable measure of this response can be defined using the average gene concentration, 2.11 and computing the so-called spectral amplification factor R , 2.12 R is nothing but the normalized amplitude of the Fourier component at the forcing frequency ω of the time series X(t). We will show that, under some circumstances, the response R will increase with the intrinsic diversity σ and that the period of the oscillations at the global level coincides with that of the external forcing, these being the main results of this paper.
The synchronization properties of the set of circadian oscillators are influenced by the amplitude of the external forcing L0, the coupling strength K and the diversity in the individual periods σ. The role of the first two has been studied in Bernard et al. , Gonze et al.  and Becker-Weimann et al. . In this section, we focus on the heterogeneity of neuronal periods and analyse the combined influence of L0, K and σ on the different parameters quantifying interneuronal synchronization and response to the forcing.
Figure 1 shows colour plots of the parameter of synchrony ρ as a function of the diversity σ and the light intensity L0, for different values of the coupling strength K. High values of the light intensity L0 favour interneuronal synchrony. Also in agreement with its intuitive disordering role, high neuronal diversity leads to a low synchrony parameter ρ in several parts of the diagrams. However, there is a region of values of L0∈ [0,Lmax] for which there is a non-monotonic dependence of the synchrony order parameter with respect to the diversity. This can be seen more clearly in figures 2a and 3a where we plot ρ as a function of diversity σ for fixed values of K = 0.6 and L0 = 0.005. ρ first decreases by increasing σ within the interval 0 ≤ σ≤ 0.05, but then it develops to a maximum. The range of values of L0 for which this non-monotonic behaviour is observed depends on the coupling constant K: the larger is K, the larger is Lmax.
As stated before, the fact that neurons synchronize among themselves does not mean that they synchronize to the forcing by light. To study this point, we have computed the individual periods Ti, i = 1, … , N, of the oscillators in the ensemble. In those cases in which the concentrations do not oscillate with exact periodicity, we still define the period as the average time between maxima of the dynamical variables. In figure 4, we plot the mean value as a function of σ and L0 for different values of K. As the dispersion in Ti is small, it turns out that is close to the period of the average variable X(t).
Although, by construction, individual neurons have periods that fluctuate around T = 23.5 h, it turns out that the period of the resulting synchronized oscillations that occur in the unforced but coupled (L0 = 0, K > 0) case increases with increasing coupling K. For example, h for K = 0.6, mostly independent of the value of σ. As the forcing sets in, at low values of the coupling strength, the mean period is now h for all values of L0 and σ. As the coupling between neurons increases, larger values of L0 and/or σ are needed in order for the mean period to coincide with that of the external forcing. An important feature that emerges from these plots is that for low light intensity it is possible to achieve a mean period of 24 h by increasing the neuronal diversity. For example, in the areas at the left of the different panels of figure 4, or figures 2b and 3b corresponding to the case K = 0.6, while identical neurons have periods of ≈30 h, increasing σ induces an adjustment of the period to 24 h. The transition between h and h is rather sharp, specially for large K. This is a clear manifestation that diversity indeed is able to improve the response to the external forcing. The same conclusion about the constructive role of diversity can be reached by looking at the measure of response R (figures 5, 2c and 3c). These figures show that there is a region in parameter space in which the system response to the periodic light forcing displays a maximum value as a function of diversity σ. This indicates that it is possible to improve neuronal synchronization to the daily-varying light input by taking σ close to an optimal value. Too small or too large a diversity will not yield an optimal response. This is a clear manifestation that diversity indeed is able to improve the response to the external forcing.
A complementary perspective on this constructive role of diversity is attained looking at spectral amplification factor, R, from equation (2.12). This is a normalized measure of the amplitude of the oscillation of the neuronal system at the frequency of the daily forcing. Figures 5 and 2c show that there is a region in parameter space in which the system response to the periodic light forcing displays a maximum value as a function of diversity σ. In fact this maximum is very large when compared with the R value at zero diversity, so that one can say that one of the most noticeable effects of a non-vanishing neuronal diversity is to give the system the capacity to respond efficiently to the 24 h forcing in situations of small or no response at this frequency in the absence of diversity (the non-diverse neuronal ensemble could be oscillating at a different frequency, as revealed by high values of ρ). In summary, it is possible to largely improve neuronal synchronization to the daily-varying light input by taking σ close to an optimal value. Too small or too large a diversity will not yield an optimal response at this frequency, although the response is generally larger than for zero diversity.
An external signal of square wave form and with different day lengths leads to similar results. As can be seen in figure 3, the response R to the external signal passes through a maximum at an intermediate value of diversity. The mean period and the synchrony parameter behave as in the case with a pure sinusoidal as the driving force. Furthermore, the qualitative result is independent of the chosen day length.
3.1. Diversity-induced oscillator death
Why does an increase in the diversity of the oscillators lead to an improved response to the external forcing? We argue that the main effect of the increase of the diversity is to take the oscillators into a regime of oscillator death [13,14] in which they can be easily entrained by the varying part of the forcing. To understand this mechanism, we first split the forcing into a constant (the mean) and a time varying part: L(t) = L0/2 + (L0/2)sin(ωt). Taking only the constant part, L(t) = L0/2, figure 6a–c shows that the oscillators go from self-sustained oscillations to oscillator death, i.e. the amplitude of the self-sustained oscillations decreases, as σ increases. Once oscillators are damped, they would respond quasi-linearly to periodic forcing, at least if this forcing is not too large, and linear oscillators always become synchronized to the external forcing, independently of their internal frequency. This is consistent with what is seen in figure 6d–f, where the neurons in the case of low heterogeneity oscillate synchronously with each other, but their common period is larger than the one of the light forcing. Only when diversity brings the neurons to oscillator death can all of them be entrained to the period of the forcing signal. The mechanism is related to the one discussed by Gonze et al.  and Bernard et al. , but here we stress that neuron heterogeneity, as opposed to internal neuron parameters and couplings, is enough to damp the collective neuron oscillations and bring the system to a non-oscillating state where it can be more easily entrained. It is interesting to note that it has been shown experimentally for fruitflies that only a subset of the pacemaker neurons sustain cyclic gene expression after changing the laboratory light conditions to constant darkness, whereas the oscillations of the other pacemaker neurons are damped out . Although this does not reveal the mechanism by which the oscillations die out it suggests that some of the circadian oscillators do indeed work in the damped regime, at least in Drosophila.
An alternative way of checking this mechanism based on diversity-induced oscillator death is by analysing the stability of the steady state of the system of equations (2.5)–(2.9) when considering a constant forcing L(t) = L0/2. The numerical calculation of the fixed point of the dynamics is greatly simplified by the fact that the concentrations of the biochemical variables are the same for each one of the N neurons irrespectively of their specific value of τi. The system (2.5)–(2.9) is linearized around this steady state and the eigenvalues of the stability matrix computed for several realizations of diversity parameters τi. In each case, the positive or negative character of the real part of the eigenvalue with the largest real part indicates the instability or stability, respectively, of the fixed point solution. In figure 7, we plot the mean of that maximum real part of the eigenvalues averaged over various realizations of the time scales τi, for N = 200 coupled neurons, as a function of L0 and σ, and different values of the coupling K (see figure 2d). In every diagram, we can see that low diversity or low forcing yield an unstable steady state (yellow region). This is where self-sustained oscillations are observed. A thick black line in the contour plots indicates a zero real part. The relevance of this line separating positive from negative maximum average eigenvalues is more apparent when we note that it also delimits regions of interest in figures 1, 4 and 5.
In summary, increasing the diversity or the (constant) forcing term decreases (on average) the maximum eigenvalue of the coupled system and thus a Hopf bifurcation can be crossed backwards, so that self-oscillations disappear. When applying the periodic external forcing on the system formed by self-sustained neurons, coherence with the external frequency is difficult to achieve because there is the competing effect of mutual neuron synchronization to a different frequency. However, when the periodic external forcing is applied on the system of damped neurons, they all synchronize to the external forcing, and thus with each other since this is the only dynamical regime available to forced damped oscillators (if forcing is not too strong to excite further resonances). Increased coupling strength increases the range of unentrained self-oscillations.
Oscillator death by diversity is not particular to this system. Mirollo & Strogatz  analyse a large system of limit-cycle oscillators with mean field coupling and randomly distributed frequencies. They proved that when the coupling is sufficiently strong and the distribution of frequencies has a sufficiently large variance, the system undergoes ‘amplitude death’. In their approach, the oscillators pull each other off their limit cycles, which are translated into a stable equilibrium point for the coupled system. Thus, this mechanism suggests that the quenched noise we introduced in the system ‘pushes apart’ the limit cycles of the different neurons, so that their competition enlarges the range of parameters where fixed point behaviour is stable.
A qualitative argument explaining the diversity-induced oscillator death in our system of coupled neurons goes as follows. We know from Gonze et al.  that a single oscillator can switch from a limit cycle to a stable steady state by adding a constant mean field (the term containing F in equation (2.5) but with time-independent F) of sufficient strength to equation (2.1). A constant light forcing term has the same effect (see the zero coupling case in figure 7). Furthermore, we have observed that the amplitude of the oscillations decreases with rising diversity (compare figure 6), but the mean does not change. In a system with low diversity, we have large oscillations of F around that mean value. If this value, taken as a constant, determines a stable steady state, then we argue that the large oscillations lead the system into unstable regions, whereas by increasing σ the amplitude is decreased and the concentrations do not leave the neighbourhood of the stable fixed point, thus finding themselves damped all the time. This is a possible mechanism for the diversity-induced oscillator death phenomenon.
4. Concluding remarks
In this work, we have analysed the role of diversity in favouring the entrainment of a system of coupled circadian oscillators. We introduce non-negligible heterogeneity in the periods of all neurons in the form of quenched noise. This is achieved by rescaling the individual neuronal periods by a scaling factor drawn from a normal distribution. The system response to the light–dark cycle periodicity is studied as a function of the interneuronal coupling strength, external forcing amplitude and neuronal heterogeneity.
Most of the cases of order induced by heterogeneity or noise carried out so far [17,18,21,29,33,36,37] emphasize the fact that the diversity directly improves oscillator synchronization. In our case, the mechanism is rather different. Diversity does not improve system synchronization directly. This is achieved indirectly, leading first to a diversity-induced stabilization of the fixed points of the neurons forming the system. Once steady concentrations are asymptotically stable, it is much better entrainable by the external forcing, so that the damped neurons adapt easily to the external forcing (and then, in addition, they appear as synchronized between them). The synchronization arises, therefore, not as a result of a direct diversity-induced neuronal synchronization but indirectly, as a result of the diversity-induced oscillator death. Our results indicate, therefore, that the right amount of heterogeneity helps the extended system to respond globally in a more coherent way to the external forcing. In addition to the robustness of the results against the use of non-sinusoidal forcing we have checked that resonances in the responses to the external forcing and matching of the circadian period to the external forcing appear in more complex models, such as the 10-variable model of Bernard et al.  with diversity in the time scales τi, or the four-variable model of Gonze et al.  with heterogeneity in all the reaction rate parameters νi. We expect that a similar behaviour will be found in models of non-mammalian clocks like those of Drosophila , Arabidopsis , Neurospora  or cyanobacteria .
Of course, it is an open question whether the observed diversity in the periods of the neurons of the SCN has been tuned by evolution in order to display a maximum response to the 24 h dark–light natural cycle. A detailed experimental check of our predictions would require one to be able to vary the amount of diversity. In this sense, we suggest the possibility of using chimeric organisms  to introduce diversity in a controlled way.
We acknowledge financial support from the EU NoE BioSim, contract LSHB-CT-2004-005137, and MEC (Spain) and FEDER (EU) through project FIS2007-60327. N.K. is supported by a grant from the Govern Balear.
One contribution of 16 to a Theme Issue ‘Advancing systems medicine and therapeutics through biosimulation’.
- Received June 16, 2010.
- Accepted September 20, 2010.
- This Journal is © 2010 The Royal Society