Modeling the Emergence of Circadian Rhythms in a Clock Neuron Network

Circadian rhythms in pacemaker cells persist for weeks in constant darkness, while in other types of cells the molecular oscillations that underlie circadian rhythms damp rapidly under the same conditions. Although much progress has been made in understanding the biochemical and cellular basis of circadian rhythms, the mechanisms leading to damped or self-sustained oscillations remain largely unknown. There exist many mathematical models that reproduce the circadian rhythms in the case of a single cell of the Drosophila fly. However, not much is known about the mechanisms leading to coherent circadian oscillation in clock neuron networks. In this work we have implemented a model for a network of interacting clock neurons to describe the emergence (or damping) of circadian rhythms in Drosophila fly, in the absence of zeitgebers. Our model consists of an array of pacemakers that interact through the modulation of some parameters by a network feedback. The individual pacemakers are described by a well-known biochemical model for circadian oscillation, to which we have added degradation of PER protein by light and multiplicative noise. The network feedback is the PER protein level averaged over the whole network. In particular, we have investigated the effect of modulation of the parameters associated with (i) the control of net entrance of PER into the nucleus and (ii) the non-photic degradation of PER. Our results indicate that the modulation of PER entrance into the nucleus allows the synchronization of clock neurons, leading to coherent circadian oscillations under constant dark condition. On the other hand, the modulation of non-photic degradation cannot reset the phases of individual clocks subjected to intrinsic biochemical noise.


Introduction
Most living organisms present rhythmic phenomena whose periods range from few milliseconds to years. Circadian oscillation is an important example of this kind of phenomenon. Several recent observations have suggested that the circadian rhythm at the molecular level results from a gene regulatory network [1]. In Drosophila melanogaster the expression of the circadian clock genes occurs within approximately 150 clock neurons. The molecular mechanisms of these fundamental oscillators consist mainly of two interlocked transcriptional feedback loops involving per, tim, clk, vri and pdp1 genes [2,3]. In one loop, the CLK protein activates per expression, and the PER-TIM complex then inhibits the activity of clk. Furthermore, clk is regulated negatively by the VRI protein and positively by the protein PDP1. In the second loop, vri and pdp1 are directly activated by the CLK protein. After its synthesis, the PER protein is phosphorylated at several residues. This leads to a time delay between the rise of mRNAs and that of the PER acting as transcriptional repressor for the clk gene. Thus alternating protein production, gene repression, and protein degradation may lead to self-sustained oscillations. The circuit is further complicated by the influence of the kinase doubletime on degradation and transport of the PER protein [4]. Moreover, the degradation rate of the TIM protein is indirectly controlled by light, which enables entrainment with the environment [5]. Over the past years, a number of deterministic and stochastic models for the circadian clock have been proposed [6][7][8][9][10][11]. They differ largely in the detail of the specific oscillator and, consequently, in their complexity.
Recent observations have revealed a more complex organization in the fly brain. They indicate that the neurons (*150) are organized into functional groups. Each group contributes to the control of the rhythmic behavior in an orchestrated manner. In Drosophila, the clock neurons have been divided into two major groups: the lateral neurons (LN) and the dorsal neurons (DN). The lateral neurons are subdivided into three subgroups: the large and small ventrolateral neurons (LN V s), and the dorsolateral neurons (LN D s). There is ample evidence indicating that the small LN V s are especially important as circadian pacemakers for locomotor rhythms in constant darkness (DD) [12,13]. It is known that the molecular oscillations that underlie circadian rhythms in the small LN V s, persist for weeks in constant darkness [14,15], while in other Drosophila tissues these oscillations damp rapidly in the absence of light [16,17]. This is an indication that the coherent circadian output results from synchronized electrical activity of the LN group [18,19] and suggests that the neuronal interactions of the various groups of clock neurons in the fly brain play a critical role in the circadian behavioral rhythms.
Today there is evidence suggesting that the interneuron synchronization in Drosophila is achieved through the pigment-dispersing factor (PDF) neuropeptide, which is specifically expressed in the ventral lateral neuron group (LN V ) [18,[20][21][22]. Pdf 01 mutant flies gradually lose behavioral rhythms in constant darkness, even though the molecular oscillations persist within individual LN V neurons [21,22]. This loss of oscillatory behavior in pdf 01 mutant flies has been interpreted as a loss of neuron synchrony due to the noisy nature of the molecular clock machinery. Therefore, in constant darkness, PDF is required to maintain the circadian rhythm of a group of neurons, though it is not required to maintain circadian rhythm in a single neuron [21,22]. Thus, in the presence of zeitgebers the phases of individual clocks are reset, but in their absence, the molecular fluctuations disperse the phase of individual clocks relative to each other, leading to the oscillation damping. Despite this progress towards understanding the genesis of circadian rhythms, the precise action of the synchronizing agent over the clock machinery is essentially unknown at the present time.
The mathematical models mentioned above were developed to explain cell-autonomous oscillations and could not explain some synchronization aspects of the cellular clocks in constant darkness. Models that take into account interactions between the clock cells themselves have been proposed [23][24][25][26][27][28]. More recently, detailed models for a clock neuron network in mammals have also been proposed [29,30]. These multiscale models provide details of the suprachiasmatic nucleus organization at both the gene regulatory and electrophysiological levels.
In this paper, we present a model of clock network to assess the putative mechanisms through which the clock neurons in the fly brain are coordinated to produce a circadian coherent output even under constant darkness condition. The model consists of an array of connected circadian pacemakers. The individual pacemakers are described by a well-known biochemical model for circadian oscillation [6], to which we have added degradation of PER protein by light and multiplicative noise. Each clock neuron has some of its parameters modulated by the PER protein level of all the other clock neurons. This corresponds to a fully connected network, without self-interactions. In particular we have investigated the effect of modulation of the parameters associated with (i) the control of net entrance of PER into the nucleus; and (ii) the non-photic degradation of PER.

The Model
To model the synchronization of the oscillations in the small LN V s group, we assume that, at the level of individual clock neurons, the circadian clock is represented by the core negative feedback loop established by PER. For the individual clock neurons we have adapted a model originally proposed by Goldbeter [6] that explicitly includes (i) transcription: the gene is transcribed into mRNA (M p ) in the absence of phosphorylated PER in the nucleus (P N ), assuming, that the repression is cooperative; (ii) translation: a portion of this mRNA is degraded, and another portion is translated into PER protein (P 0 ) in the cytoplasm; (iii) phosphorylation: PER protein is phosphorylated in a reversible way twice (from P 0 to P 1 and from P 1 to P 2 ); (iv) degradation: the fully phosphorylated PER (P 2 ) is degraded by the default molecular machinery following a Michaelis-Menten rate expression; (v) transport: the entrance of PER into the nucleus is assumed to be a reversible first-order process. These processes correspond to the model introduced by Goldbeter [6] for autonomous oscillations. In order to simulate the light effect we have added a PER degradation process due to light exposure: the fully phosphorylated PER (P 2 ) is linearly degraded by light exposure at the rate k light . Figure 1 displays the model schematic diagram.
The stochastic biochemical processes underlying the molecular machinery are subject to noise or fluctuations [7,31]. For this reason, all the above mentioned processes are affected by molecular noise, this noise term being assumed to be multiplicative. The temporal evolution of the above-mentioned chemical species is then governed by the following set of differential equations: Figure 1. Sketch of the model (adapted from [7]). A portion of per mRNA (M P ) is translated to PER protein in the cytosol (P 0 ), where it undergoes two phosphorylations. Part of the fully phosphorylated PER protein (P 2 ) enters into the nucleus (dashed box) and the remains are degraded in the cytosol either induced by light or basally. Per activity is inhibited by P 2 forming a negative feedback loop. The Input Network box represents hypothetical regulation of the P 2 entrance and/or degradation. doi:10.1371/journal.pone.0033912.g001 Table 1. Parameter values used in the simulations.
This cellular model incorporates enough details about the PER loop, thus allowing us to test some of the hypotheses regarding the emergence of circadian rhythms in a network of noisy neurons. It should be remarked that the model described by Eqs. (1) is reduced to the original Goldbeter model [6] when g~0 and k light~0 . The additional term representing the degradation of PER by light is null during the night. In our simulations this is accomplished by setting the parameter k light equal to zero. This additional degradation term allows the synchronization of the clocks by an external clue (the light). However, under DD conditions, when there is no light-induced degradation, each clock neuron oscillates almost independently of each other. The molecular fluctuations mentioned above could disperse the rhythmic phase of the clock neurons relative to each other, leading to the loss of synchrony between the oscillator units [7]. In this way the overall oscillatory output of the network is damped. Neuronal interactions may also be required to maintain proper molecular oscillations in constant darkness. Our model assumes that the clock neurons form a network, and that the interneuron communication between the clock neurons is through the current value of cytoplasmic P 2 concentration. The input for each clock neuron is the averaged P 2 concentration over all other clock neurons with equal weight. So the input signal over each clock neuron is S N~N where N is the number of clock neurons in the network, which was set to be N~50 in all cases, unless otherwise specified. This corresponds to a fully interconnected network, without self-interactions. Unlike to the models [23,25] that use a heterogeneous population of clock neurons, the network of our model is formed by identical cells. In this case, the dispersion in the oscillations' period and phase arises from the molecular fluctuations rather than from a heterogeneous clock population. As the precise mechanism for maintaining circadian regime under free-running condition remains essentially unknown, we will test the following alternative hypotheses: (i) The input signal S N modulates the net rate of the fully phosphorylated PER entrance into the nucleus. This is implemented by increasing and/or decreasing the parameters k 1 and k 2 : The input signal S N modulates the degradation rate of the fully phosphorylated PER. This is implemented by increasing or decreasing the degradation constant v d ?v d +aS N .

Numerical simulations
In our simulations we have considered a network of N~50 clock neurons, where each clock neuron obeys Eqs. (1). They are connected through the network input that modulates the parameters k 1 , k 2 and v d , as explained above.
The numerical integration was performed by using the Runge-Kutta method, with the integration time step set to 0.1 min. To obtain fluctuations similar to those exhibited by a stochastic version [7] of the model we have used g~2:0. Other parameter values used in all simulations are given in Table 1. The initial condition of the i-th clock neuron corresponds to the chemical state of a referential clock neuron at the time t~t i , where t i is a Gaussian distributed random number with mean 0 h and standard deviation of 1 h. Thus the initial condition of the network corresponds to an ensemble of pacemakers whose phases are Gaussian distributed around 0 h. Using this initial condition, the network is subjected to 9 days of LD (12 h of light and 12 h of darkness), followed by 7 days under DD condition (24 h of darkness) or free-running condition. In our simulation the absence of light (darkness) is accomplished by setting k light~0 .
For each set of parameter values we determine the amplitude, the period and the phase of the P N level. This is accomplished by a nonlinear fitting of a cosine function of the form where b is the baseline, A the amplitude, T is the period, and t s the phase shift. Thus, the phase shift corresponds to the time interval from day onset to the P N peak. The fitting is done for several fixed values of T ranging from 16 h to 32 h, with a resolution of 15 min. For each clock neuron, we select the period T corresponding to the best fitting.
It should be remarked that there are two procedures for obtaining the average period and the average phase of the network: (i) averaging the period and phase obtained for each individual clock neuron over the clock neuron population; (ii) first averaging the P N level over the clock neuron population, and then fitting the cosine function for the average level SP N T. These procedures for obtaining the average phase and period lead to the same results when the clock neurons are well synchronized. However, under DD condition, there is a loss of synchronization that causes a decrease in the amplitude, which that does not allow the precise determination of the average value of the phase and period using the procedure (ii) above. For this reason we have adopted procedure (i) to determine the period and phase of the clock network.
The overall degree of synchronization over a specified time interval is analyzed by computing the parameter R This parameter R has been computed using the concentration of PER in the nucleus, i.e., x~P N . Furthermore, the parameter R, the period T, and phase shift t s were obtained by averaging over the last 5 days, independently of the case (LD condition or DD condition). Figure 2 displays the temporal behavior of P N of the noninteracting neuron population for 6 days under LD condition followed by 6 days under DD condition (v d~3 :0 and k light~0 :65, other parameter values being given in Table 1). We can see that in the presence of the zeitgeber the clock neurons are synchronized, but in free running (DD condition) there is a loss of synchronization. Nevertheless, the molecular oscillations of each clock neuron are little affected by the absence of the zeitgeber.  high values of k light and low values of v d where the parameter R reaches 0.75. In the bottom-left panel (LD condition) we can see that v d must be smaller than 2.85 to obtain circadian oscillation; for higher values of v d the period is shorter than 24 hs. At the topright panel one can observe any degree of synchronized behavior, and that oscillations in the range 2:7vv d v3:0 have periods longer than 24 hs. Now we will test the hypotheses regarding the effects of interneuron interaction on synchronization:

Results
(i) We assume that the input signal S N increases the rate of entrance of P 2 into the nucleus, this being accomplished by replacing k 1 by k 1 za 1 S N in Eqs. (1). Figure 4 Figure 4 shows a small region around v d~2 :91 and k light~0 :65 that satisfies these requirements when k 1 is increased by the interneuron communication. Figure 5A   (indicated by the white/black bars) followed by seven days under DD condition. The resulting oscillations preserve the 24 h period under both LD and DD conditions, and there is no phase shift between the LD and DD conditions. A decrease in the oscillation amplitude is observed under the DD condition due to a lower degree of interneuron synchronization.
In addition, we have also determined the period and phase of the oscillations of the individual clock neurons by fitting the time course of P N to the function (2). Figure 5B shows the histogram of phase shift (B1 for LD and B2 for DD), and Fig. 5C shows the period histogram of the cell population (C1 for LD condition, and C2 for DD condition). One can observe that our model suggests that the dispersion both in the oscillation period and in the phase shift contributes to a decrease in the amplitude of oscillation of the average P N . Figure 6 displays the temporal behavior of P N of 10 interacting neurons (from a population of 50 clock neurons) in 6 days during LD condition followed by 6 days during DD condition for v d~2 :91 and k light~0 :65. We can see that in free-running (DD) condition the synchronization of the clock neurons is similar to that obtained in the presence of the zeitgebers.
We have also done numerical simulations for two other values of a 1 : 0:4 and 1:2. For a 1~0 :4, we found that under DD condition there exists circadian oscillation for 2:79vv d v2:82, although the degree of synchronization is poor, Rv0:2 (see Fig. 7). For a 1~1 :2 (see Fig. 8) we found that, for DD condition, there exists circadian oscillation for 0:5vk light v0:85 and v d^2 :8 or v d^3 :0, the degree of synchronization being 0.2 and 0.4, respectively. However, in these parameter regions the oscillation has a period T smaller than 24 h under the LD condition.
We have also analyzed the effect of increasing k 1 and decreasing k 2 at the same time. Numerical simulations using a 1~0 :8 and a 2~{ 0:4 show that, for both LD and DD conditions, the degree of synchronization R is high in the k light interval ½0:5,0:65 in the whole range of the parameter v d used in the simulations (top panels of Fig. 9). Circadian oscillations are obtained for both LD and DD conditions in the region around the segment that links the points (0.4,2.95) and (0.65,2.80) in the (k light ,v d )-space (bottom panels of Fig. 9). The top panels of Fig. 9 suggest that this type of modulation produces better synchronization, and that the parameter regions where the system presents circadian oscillations are larger than in the case a 1~0 :8 displayed in Fig. 4 for which k 2 is not modulated by the input network feedback S N . Figure 10A depicts the time course (double plot actogram) of P N averaged over the clock neuron population when v d~2 :95 and k light~0 :45. We can observe that the interneuron communication, mediated by the simultaneous modulation of k 1 and k 2 , produces oscillations that preserves the period and amplitude, but introduces a phase shift of 6 h in advance. Fig. 10B shows the histogram of the phase shift, and Fig. 10C the periods of the cells. The period and phase of each clock neuron were obtained by fitting the time course of P N to the function (2).  Table 1). The 7 days under LD condition are indicated by the white-black bars. B1, B2: Histograms of phase shift. C1, C2: period obtained by fitting Eq. (2) to the time course of P N concentration at each clock neuron. B1, C1 correspond to LD condition, and B2, C2 to DD condition. doi:10.1371/journal.pone.0033912.g010 (ii) Now we assume that the input network feedback S N increases the non-photic degradation rate v d , this being accomplished by replacing v d by v d z0:8S N in Eqs. (1). In this case the positive modulation of the v d parameter does not lead to synchronization in the whole region of parameters v d and k light used in the simulation (Fig. 11). Furthermore, the oscillation period T is substantially larger than 24 h. However, a decrease in the non-photic degradation v d can lead to synchronization. Figure 12 displays the results obtained when v d is decreased according to v d {0:3S N , showing that in free running there is synchronization and circadian oscillation. However, for the set of parameter values studied here the associated oscillations have periods longer than 25 h.

Discussion
During the past years many important studies have elucidated the mechanisms underlying circadian oscillations at the molecular level. These mechanisms are common to several types of cells. However, the persistence of the oscillation at tissue level is tissue specific, even though circadian rhythms at the single cell level persist for weeks in constant darkness. It is known that in the ventral lateral neuron group the neuropeptide PDF is required to maintain the circadian rhythm of this group under DD condition. However, it is not known how this neuropeptide acts to promote the synchronization of the clock neurons. In order to uncover the putative mechanisms (or discard wrong hypotheses) that allow us to explain the observed behavior in Drosophila, we have developed a mathematical model for a network of interacting clock neurons, where some parameters of each clock neuron are modulated by the PER protein level averaged over the whole network. Unlike other clock network models [25], which use a heterogenous population of neurons, the network in our model is formed by identical cells. The dispersion in the period and phase of the oscillations arises from the molecular fluctuations rather than from clocks with different parameters.
In contrast with [29,30], our simple model has not considered a neurotransmitter hypothesis. The literature is scarce regarding the mechanistic details connecting the clock neuron state to the release of neurotransmitters. So our model directly uses the doublephosphorylated PER as synchronizer agent. This assumption would not adequately represent the real world. Nevertheless, it admits a more realistic interpretation: a particular state of the clock (the high concentration of PER in our simple model) promotes PDF releases.
With our model we have studied the degree of synchronization of oscillators, and the period and phase of the oscillation, at the network level, in the space of parameters for some alternative hypotheses. In particular, we have investigated the effect on the synchronization when (i) a 1 is varied so as to modulate the Figure 11. Positive modulation of non-photic degradation does not lead to synchronization. R (top) and period (bottom), in the (k light ,v d )-space, obtained by averaging over the cell population in the last 5 days, when the non-photic degradation v d is replaced by v d z0:8S N . Panels on the left correspond to LD condition, and panels on the right correspond to DD condition (k light~0 ). doi:10.1371/journal.pone.0033912.g011 parameter k 1 that controls the net entrance of PER into the nucleus; and (ii) a is varied so as to modulate the parameter v d that controls the non-photic degradation of PER.
Our results indicate that for a 1~0 :8, the modulation of PER entrance into the nucleus allows the synchronization of clock neurons leading to coherent circadian oscillations under constant darkness condition. These results were obtained for a network with N~50, a number sensibly greater than the actual number of lateral neurons in the fly brain. However, similar circadian oscillations have been also observed for networks with a different number of clock neurons (N~10, 25, 80, 150 and 200, data not shown). This lack of dependence on N is likely due to the fact that S N is normalized by N. We have also tested other values of the parameter a 1 : for lower modulation (a 1~0 :4) there is a poor synchronization of the clock neurons in free running, and for higher modulation (a 1~1 :2) the synchronization increases but the oscillations have period shorter than 24 h. In contrast, the modulation of the non-photic degradation cannot reset the phases of individual clocks subjected to intrinsic biochemical noise. The present network model is able to display circadian oscillation even in the absence of external zeitgebers when the input feedback signal S N (defined as the average of P 2 over the entire network) increases the P 2 entrance rate into the nucleus. There exists a small region in the parameter space that supports circadian oscillation both under the LD and DD conditions simultaneously, the phase shift between these conditions being small. We have also observed that there exist several ways to reach a high degree of neuronal synchronization, but at the expense of noncircadian oscillations.
Our results indicate that mechanisms based on a positive feedback acting over the rate of entrance of the phosphorylated PER into the nucleus could be essential for maintaining the circadian oscillation under free-running condition. This fact suggests a putative way of action for the neuropeptide PDF, which could be acting as an agent that promotes the entrance of the phosphorylated PER into the nucleus.