On-line version ISSN 1414-431X
Braz J Med Biol Res vol.42 no.1 Ribeirão Preto Jan. 2009
A simple model for circadian timing by mammals
1Departamento de Ciências Fisiológicas, Universidade do Estado do Rio de Janeiro, Rio de Janeiro, RJ, Brasil
2Departamento de Física, Universidade Federal Rural do Rio de Janeiro, Rio de Janeiro, RJ, Brasil
Circadian timing is structured in such a way as to receive information from the external and internal environments, and its function is the timing organization of the physiological and behavioral processes in a circadian pattern. In mammals, the circadian timing system consists of a group of structures, which includes the suprachiasmatic nucleus (SCN), the intergeniculate leaflet and the pineal gland. Neuron groups working as a biological pacemaker are found in the SCN, forming a biological master clock. We present here a simple model for the circadian timing system of mammals, which is able to reproduce two fundamental characteristics of biological rhythms: the endogenous generation of pulses and synchronization with the light-dark cycle. In this model, the biological pacemaker of the SCN was modeled as a set of 1000 homogeneously distributed coupled oscillators with long-range coupling forming a spherical lattice. The characteristics of the oscillator set were defined taking into account the Kuramoto's oscillator dynamics, but we used a new method for estimating the equilibrium order parameter. Simultaneous activities of the excitatory and inhibitory synapses on the elements of the circadian timing circuit at each instant were modeled by specific equations for synaptic events. All simulation programs were written in Fortran 77, compiled and run on PC DOS computers. Our model exhibited responses in agreement with physiological patterns. The values of output frequency of the oscillator system (maximal value of 3.9 Hz) were of the order of magnitude of the firing frequencies recorded in suprachiasmatic neurons of rodents in vivo and in vitro (from 1.8 to 5.4 Hz).
Key words: Circadian cycle; Coupling oscillators; Kuramoto's model
Over the last few years, behaviors with periodic oscillations in biological processes have attracted the interest of many areas of science (1) and it has been demonstrated that light has non-visual systemic effects in healthy humans (2). Physico-mathematical modeling and computational simulation have proven to be very useful for the understanding of the factors that influence these oscillations.
Biological rhythms present frequencies ranging from a millisecond to several years and are genetically determined, evolutionally preserved and endogenously generated. Periodic oscillations ranging from 20 to 28 h are called circadian rhythms and are generated by a circadian timing system (CTS), which is organized in such a way as to receive information from the external environment (mainly from the light-dark cycle) and the internal one (associated with the physiological reactions occurring in the organism). The CTS function is driving most daily rhythms in physiology and behavior in a circadian pattern in order to allow the expression of a "maximal adaptive behavior" (3). Body temperature and metabolic rate have strong circadian patterns (4). The CTS controls food processing and energy homeostasis by regulating the expression and/or activity of enzymes involved in cholesterol, amino acid, lipid, glycogen, and glucose metabolism. In addition, many hormones involved in metabolism, such as insulin, glucagon, adiponectin, corticosterone, leptin, and ghrelin, exhibit circadian oscillation (5).
The circadian rhythms possess two fundamental characteristics: endogenous generation and synchronization with the light-dark cycle (6). The endogenous generation is probably caused by the sequenced pattern of the organic reactions at the cellular and molecular level of biological organization, allowing the organic metabolism to occur in phases that permit its manifestation in later stages at the same frequency, thus functioning as a timer mechanism as they occur (7).
In mammals, the CTS is a hierarchical system consisting of a group of structures including the suprachiasmatic nucleus (SCN), the intergeniculate leaflet (IGL), the pineal gland, and other structures (8). The SCN is a paired structure found in the anterior hypothalamus on each side of the third ventriculus, lying on the optic chiasm. In the rat, the SCN has an approximately ovoid shape, with a volume of approximately 0.066 mm3. From the neurochemical point of view, the SCN has various distinct neuronal subpopulations characterized by the presence of different neuroactive substances found in the terminals coming from its afference sources (9).
The SCN receives the afferences from the retinohypothalamic tract, geniculohypothalamic tract and raphe (10). The retinohypothalamic tract is formed by fibers from the ganglion cells of the retina and is necessary and sufficient for the maintenance of synchronization with the light-dark cycle (11). The geniculohypothalamic tract originates in the IGL and plays an important role in the synchronization of the circadian rhythm (12).
Biological time keeping is coordinated by a hypothalamic pacemaker located in the SCN, known as a central biological pacemaker or master clock (4). It is synchronized with solar time by direct retinal afferents. Analytical action spectra in rodents, primates, and humans have identified 446-484 nm (predominantly the blue part of the spectrum) as the most potent wavelength region for neuroendocrine, circadian, and neurobehavioral responses (2). Most suprachiasmatic neurons form a biological clock that drives circadian rhythms in vivo and in vitro (13). Individual SCN neurons generate time because they are coupled with each other and their intrinsic oscillator consists of a series of interlinked autoregulatory transcriptional/post-translational feedback loops. Thus, the clock is an intracellular transcriptional mechanism sharing the same molecular components in SCN neurons (9). Clock activity in the SCN neurons is expressed as activity-dependent electric signals (14).
Important recent findings have suggested clock-like activities in many neural and non-neural tissues. Certain areas of the brain, including other hypothalamic nuclei, the amygdala, olfactory bulb, and lateral habenula, express circadian rhythms in core clock gene expression, hormone output and electrical activity (3), as is also the case for peripheral cells such as the liver, intestine, and retina (5).
We present here a simple model for the CTS using a set of 1000 coupled oscillators homogeneously distributed as a spherical lattice with long-range coupling to represent the biological pacemaker of SCN. The characteristics of this oscillator set were defined based on Kuramoto's oscillator dynamics (15), but we used a new method for estimating the equilibrium order parameter (16). Our model has been able to reproduce the two important properties of the mammalian CTS, i.e., the endogenous generation of pulses and the response to the light-dark cycle.
The neural circuit model adopted to simulate the CTS is presented in Figure 1, and its basic elements as well as the equations used are described later.
The E1 element or sensorial element represents the group of retina photoreceptors. The signals generated in element E1 are sent to element E2 mainly via the retinohypothalamic tract. To determine E1 output, we must consider the temporal variation of luminous intensity on the retina as a Gaussian function given by:
where i0 is the initial luminous intensity on the retina, ω is the width at half-height of the Gauss curve, tc is the instant of maximal luminosity and t is a given instant within the 24-h interval. The maximum function corresponds to the instant of highest light intensity on the retina and highest activity in E1.
Based on experimental rat models, we have considered in our modeling of the mammalian CTS a subjective day beginning at 6:00 am with a luminous incidence of 400 lux (luminous intensity) on the retina. From this instant, the luminous intensity increases (according to Equation 1) until midday, when it reaches a maximal value of 800 lux. From midday, the luminous intensity decreases to 400 lux at 18:00 h. To simulate phototransduction, we have used the following equation:
where f0 represents the firing frequency in the dark phase and kf is the inverse of mean luminous intensity.
Figure 2 shows the variation in output frequency of E1 with time along 96 h. From 6:00 to 12:00 h, the output frequency gradually increases from 4 Hz to a maximal value of 62 Hz occurring at the middle of the light phase.
The E2 element represents a population of excitatory and inhibitory neurons connected with each other in the SCN. In E2, the neurons have similar patterns of synaptic connections and dynamic behavior and represent a neural oscillator circuit. This element has been modeled as a set of 1000 coupled oscillators homogeneously distributed as a 10 x 10 x 10 spherical lattice with long-range coupling. The synchronism of the oscillator system was studied using the following equation:
where Θ is the phase and F(Θi) = f(Θi(t + 1)) represents the coupling function, which is a harmonic oscillator function in our model, k is coupling strength, and rij is distance between oscillators i and j. For α = 0, the coupling is uniform, and α → ∞ corresponds to the nearest neighbor coupling. We have assumed that F(Θi) = A sin(Θ), for Θ(t) = Θ0 + ωt for -π / 2 ≤ ωt ≤ π / 2.
Studies using coupled oscillators require the estimation of a certain normalized parameter, order parameter, to characterize the degree of synchronization in the asymptotic regime. According to Kuramoto (15), the order parameter related to phase can be written as:
where δ is the order parameter and Φj is the phase of oscillator j.
In a previous study, Cruz and Cortez (16) presented a new method for estimating the equilibrium order parameter (δ). This is a method that adopts the statistical analysis of the phases of oscillators (Θ) to study the average behavior of oscillators and cluster formation. In this method, δ is given by δ = exp(-γ), where γ is Pearson's variation coefficient of different Θ values.
Thus, the Cruz-Cortez (16) method uses the evaluation of the mean value and the variance of the Θ-distribution within the oscillator system to estimate the g value: g = sΘ / Θmean (sΘ = variance of phases, Θmean = arithmetic mean value of phases). The variation coefficient is usually given as percentage and becomes null (γ = 0) for equal Θ values, with the following correspondence being normally used: γ< 10, low variation of values; 20<γ ≤ 10, intermediate variation of values; 30 < γ ≤ 20, high variation of values; 30 ≤ γ, very disperse values. Similarly, the following classification can be used for the order parameter: δ = 1 for a coherent system or completely synchronized state, 0.9 < δ < 1 for low variation of Θ-values or an almost synchronized state, and δ ≤ 0.9 for incoherent oscillations or a desynchronized state. To determine the phase mean value, we used 50 phase measures. At t = 0, we assume that oscillations are incoherent and δ = 0.
In Figure 1 one can see the E2 element receiving afferences directly from elements E1, E3 and E4. E2 input has been considered to consist of trains of action potentials (APs) whose excitatory and inhibitory postsynaptic effects are given by:
where Ej and Ij are the characteristic amplitudes of the excitatory and inhibitory postsynaptic potentials, respectively, wSEj and wSIj are the synaptic weight characteristic of the jth excitatory and inhibitory synapses, respectively, n is the total amount of synapses, and ti, the time of arrival of the ith presynaptic AP plus synaptic delay. The total amount of synaptic inputs is represented by p and q, and τE and τI are decay time constants.
Since simultaneous inhibitory and excitatory postsynaptic potential generate by widely distributed synapses can summate, postsynaptic neurons fire when a sufficiently large depolarizing effect is provided and the plasma membrane potential, V(t), becomes larger than the threshold potential, T(t). In our model, these potentials are given by Equation 7:
where PR is the resting potential, Po is the synaptic reversal potential, τP is the time constant of the period after potential hyperpolarization, TR is the threshold potential at rest, To is the excitation threshold after AP, τH is the decay time constant for the relative refractory period, and tp is the absolute refractory period. Here, we adopted the same synaptic conditions as used by Dalcin et al. (17).
The E3 element represents several suprasegmental structures (IGL, thalamic paraventricular nucleus, hypothalamus), and the E4 element represents several segmental structures (mainly raphe nuclei). Both send signals to E2. Because of the lack of available information about segmental and suprasegmental afferences, we decided to adopt random routines to represent E3 and E4 outputs.
The programs used were written in Fortran 77, compiled and run on PC DOS computers. Numerical values of parameters found in the literature were applied to these programs (Table 1; Refs. 18-22). Some of these values were obtained from experimental studies.
Our model for the CTS of mammals incorporates the excitatory and inhibitory properties of the neuronal membrane and takes into account a set of known synaptic characteristics, including certain mechanisms of temporal and spatial bioelectric responses, such as synaptic afterfiring and the adaptive and facilitation processes shown in Equation 7. The circuit inputs model the depolarizing signals generated by activity of the E1 element (Figure 1), and the output frequency of E1 as a function of luminous intensity is given by Equation 1.
The E2 element represents the SCN and has been programmed to process signals coming from E1, E3 (suprasegmental afference) and E4 (segmental afference). E3 and E4 send excitatory and inhibitory random signals, respectively, to E2, with their frequencies being maintained close to 30 Hz. This value is within the physiological scale of central nervous activity (23). The hyperpolarizing effect produced by E3 has been considered important to avoid synaptic fatigue during the light phase of the day. Figure 3 shows typical simultaneous activities of elements E4 and E5 over a period of 20 h.
The biological pacemaker in E2 has been simulated as a set of 1000 coupled oscillators with long-range coupling, forming a spherical lattice. The k-values (Equation 3) have been estimated based on the result of synaptic effects (Equations 5-7) on element E2 from elements E1, E3 and E4. For phases of higher and lower light intensities, the k-values were 0.8 and 0.1, respectively. Figure 3 shows the variation of the input frequency of the oscillator system in E2, illustrating the summation of simultaneous inhibitory and excitatory postsynaptic effects on E2 due to widely distributed synapses.
Figure 4 presents the mean output frequency of CTS (or output of the oscillator system), which ranges from 0 Hz to a maximal value of 3.9 Hz within the time interval from 6:00 to 12:00 h (midday). After 12:00 h, the output frequency of CTS falls until it reaches 0 Hz at 18:00 h. These results agree with experimental data reported in the literature. Suprachiasmatic neurons presenting spontaneous activity have been found by several investigators in vivo and in vitro experiments, and firing frequencies recorded in some of these experiments in rodents have varied from 1.8 to 5.4 Hz (13,18,24). Discharges of about 5.4 Hz and two troughs of about 1.8 and 2.4 Hz were recorded by Zlomanczuk et al. (24) from SCN neurons of hamsters, which exhibited a split in overt behavior. The troughs coincided with the projected time of wheel-running activity. When this activity was suppressed with high-intensity (500 lux) constant light, the SCN firing profile failed to exhibit a daily rhythm, with firing frequencies of about 5.8 Hz throughout the 24-h period.
In conclusion, the values of output frequency of the oscillator system have been of the order of magnitude of the firing frequencies recorded in SCN neurons in vivo and in vitro. The model equations were able to describe the behavior of elements of the adopted circuit. In addition, the simulation environment allows model expansion to include additional neuronal geometry and/or cell populations in the model.
[View larger version of this image (87 K JPG file)]
[View larger version of this image (73 K JPG file)]
[View larger version of this image (80 K JPG file)]
[View larger version of this image (73 K JPG file)]
|Table 1. Values of the biological parameters used.|
[View larger version of this table (129 K JPG file)]
1. Ichikawa T. Mutual coupling among insect neurosecretory cells with an ultradian firing rhythm. Neurosci Lett 2001; 299: 73-76. [ Links ]
2. Hanifin JP, Brainard GC. Photoreception for circadian, neuroendocrine, and neurobehavioral regulation. J Physiol Anthropol 2007; 26: 87-94. [ Links ]
3. Moore RY. The suprachiasmatic nucleus and the circadian timing system. In: Klein DC, Moore RY, Reppert SM (Editors), Suprachiasmatic nucleus: the mind's clock. New York: Oxford UP; 1991. p 13-15. [ Links ]
4. Mortola JP. Hypoxia and circadian patterns. Respir Physiol Neurobiol 2007; 158: 274-279. [ Links ]
5. Froy O. The relationship between nutrition and circadian rhythms in mammals. Front Neuroendocrinol 2007; 28: 61-71. [ Links ]
6. Morin LP. The circadian visual system. Brain Res Brain Res Rev 1994; 19: 102-127. [ Links ]
7. Cipolla-Neto J, Afeche SC, Menna-Barreto L, Marques N, Benedito-Silva AA, Fortunato G, et al. Lack of similarity between the effect of lesions of the suprachiasmatic nucleus and subparaventricular hypothalamic zone on behavioral circadian rhythms. Braz J Med Biol Res 1988; 21: 653-654. [ Links ]
8. Bob P, Fedor-Freybergh P. Melatonin, consciousness, and traumatic stress. J Pineal Res 2008; 44: 341-347. [ Links ]
9. van den Pol AN, Tsujimoto KL. Neurotransmitters of the hypothalamic suprachiasmatic nucleus: immunocytochemical analysis of 25 neuronal antigens. Neuroscience 1985; 15: 1049-1086. [ Links ]
10. Cassone VM, Speh JC, Card JP, Moore RY. Comparative anatomy of the mammalian hypothalamic suprachiasmatic nucleus. J Biol Rhythms 1988; 3: 71-91. [ Links ]
11. Nishino H. Activity of neurons of the hypothalamic suprachiasmatic nuclei and circadian rhythms: the role of the suprachiasmatic nuclei in rhythmic activities of neurons in the lateral hypothalamic area and ventromedian nuclei and the control of the pineal gland. Nippon Yakurigaku Zasshi 1976; 72: 941-954. [ Links ]
12. Wickland C, Turek FW. Lesions of the thalamic intergeniculate leaflet block activity-induced phase shifts in the circadian activity rhythm of the golden hamster. Brain Res 1994; 660: 293-300. [ Links ]
13. Aujard F, Herzog ED, Block GD. Circadian rhythms in firing rate of individual suprachiasmatic nucleus neurons from adult and middle-aged mice. Neuroscience 2001; 106: 255-261. [ Links ]
14. Okamura H. Suprachiasmatic nucleus clock time in the mammalian circadian system. Cold Spring Harb Symp Quant Biol 2007; 72: 551-556. [ Links ]
15. Kuramoto Y. Chemical oscillation, waves and turbulence. Berlin: Springer-Verlag; 1984. [ Links ]
16. Cruz FAO, Cortez CM. Computer simulation of a central pattern generator via Kuramoto model. Physica A 2005; 353: 258-270. [ Links ]
17. Dalcin BL, Cruz FAO, Cortez CM, Passos EP. Computer modeling of a spinal reflex circuit. Braz J Physics 2005; 35: 987-994. [ Links ]
18. Walsh IB, van den Berg RJ, Marani E, Rietveld WJ. Spontaneous and stimulated firing in cultured rat suprachiasmatic neurons. Brain Res 1992; 588: 120-131. [ Links ]
19. Kudela P, Franaszczuk PJ, Bergey GK. A simple computer model of excitable synaptically connected neurons. Biol Cybern 1997; 77: 71-77. [ Links ]
20. Kim YI, Dudek FE. Membrane properties of rat suprachiasmatic nucleus neurons receiving optic nerve input. J Physiol 1993; 464: 229-243. [ Links ]
21. Carp JS. Physiological properties of primate lumbar motoneurons. J Neurophysiol 1992; 68: 1121-1132. [ Links ]
22. Mendell LM, Henneman E. Terminals of single Ia fibers: location, density, and distribution within a pool of 300 homonymous motoneurons. J Neurophysiol 1971; 34: 171-187. [ Links ]
23. Baker SN, Olivier E, Lemon RN. Coherent oscillations in monkey motor cortex and hand muscle EMG show task-dependent modulation. J Physiol 1997; 501 (Part 1): 225-241. [ Links ]
24. Zlomanczuk P, Margraf RR, Lynch GR. In vitro electrical activity in the suprachiasmatic nucleus following splitting and masking of wheel-running behavior. Brain Res 1991; 559: 94-99. [ Links ]
Address for correspondence: C.M. Cortez, Departamento de Ciências Fisiológicas, UERJ, Av. Prof. Manuel de Abreu, 444, 5° andar, 20555-170 Rio de Janeiro, RJ, Brasil. Fax: +55-21 2238-8101. E-mail: firstname.lastname@example.org
Presented at the IV Miguel R. Covian Symposium, Ribeirão Preto, SP, Brazil, May 23-25, 2008. Research supported by FAPERJ. Received June 14, 2008. Accepted January 21, 2009.