A simple model for circadian timing by mammals

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).


Introduction
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 gener-ated.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 A simple model for circadian timing by mammals www.bjournal.com.br 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 mm 3 .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.

Model adopted
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: (Equation 1) Figure 1. Figure 1. Figure 1. Figure 1. Figure 1.Schematic representation of the neural circuit adopted to simulate the circadian timing system.E1 = sensorial element; E2 = suprachiasmatic nucleus (SCN); E3 = suprasegmental afferences; E4 = segmental afferences.RHT = retinohypothalamic tract.where i 0 is the initial luminous intensity on the retina, ω is the width at half-height of the Gauss curve, t c is the instant of maximal luminosity and t is a given instant within the 24h 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: (Equation 2) where f 0 represents the firing frequency in the dark phase and k f 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 r ij 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: (Equation 4) 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 γ value: γ = s Θ / Θ med (s Θ = variance of phases, Θ -= 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 poten-     A simple model for circadian timing by mammals www.bjournal.com.brtials (APs) whose excitatory and inhibitory postsynaptic effects are given by: (Equation 5) (Equation 6) where E j and I j are the characteristic amplitudes of the excitatory and inhibitory postsynaptic potentials, respectively, w SEj and w SIj are the synaptic weight characteristic of the jth excitatory and inhibitory synapses, respectively, n is the total amount of synapses, and t i , 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 P R is the resting potential, P o is the synaptic reversal potential, τ P is the time constant of the period after potential hyperpolarization, T R is the threshold potential at rest, T o is the excitation threshold after AP, τ H is the decay time constant for the relative refractory period, and t p 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][19][20][21][22].Some of these values were obtained from experimental studies.

Results and Discussion
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 longrange coupling, forming a spherical lattice.The kvalues (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 kvalues 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 highintensity (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.Figure 3. Figure 3. Figure 3. Figure 3. Figure 3. Input frequency of the oscillator system in E2 resulting from the summation of simultaneous inhibitory and excitatory postsynaptic effects on E2 along 20 h.      4. Mean output frequency of the circadian timing system (CTS) (or output of the oscillator system) along 20 h.The frequency ranged from 0 Hz to a maximal value of 3.5 Hz within the time interval from 6:00 to 12:00 h (midday).After 12:00 h, the output frequency of the CTS fell until it reached 0 Hz at 18:00 h.
Figure 2.Figure 2.Figure 2.Figure 2.Figure 2. Variation in the output frequency of E1 with time along 96 h.
Figure 2.Figure 2.Figure 2.Figure 2.Figure 2. Variation in the output frequency of E1 with time along 96 h.
Figure 2.Figure 2.Figure 2.Figure 2.Figure 2. Variation in the output frequency of E1 with time along 96 h.

Figure 4 .
Figure 4.Figure 4.Figure 4.Figure 4.Figure4.Mean output frequency of the circadian timing system (CTS) (or output of the oscillator system) along 20 h.The frequency ranged from 0 Hz to a maximal value of 3.5 Hz within the time interval from 6:00 to 12:00 h (midday).After 12:00 h, the output frequency of the CTS fell until it reached 0 Hz at 18:00 h.
Figure 4.Figure 4.Figure 4.Figure 4.Figure4.Mean output frequency of the circadian timing system (CTS) (or output of the oscillator system) along 20 h.The frequency ranged from 0 Hz to a maximal value of 3.5 Hz within the time interval from 6:00 to 12:00 h (midday).After 12:00 h, the output frequency of the CTS fell until it reached 0 Hz at 18:00 h.
Figure 4.Figure 4.Figure 4.Figure 4.Figure4.Mean output frequency of the circadian timing system (CTS) (or output of the oscillator system) along 20 h.The frequency ranged from 0 Hz to a maximal value of 3.5 Hz within the time interval from 6:00 to 12:00 h (midday).After 12:00 h, the output frequency of the CTS fell until it reached 0 Hz at 18:00 h.
Figure 4.Figure 4.Figure 4.Figure 4.Figure4.Mean output frequency of the circadian timing system (CTS) (or output of the oscillator system) along 20 h.The frequency ranged from 0 Hz to a maximal value of 3.5 Hz within the time interval from 6:00 to 12:00 h (midday).After 12:00 h, the output frequency of the CTS fell until it reached 0 Hz at 18:00 h.

Figure
Figure 4.Figure 4.Figure 4.Figure 4.Figure4.Mean output frequency of the circadian timing system (CTS) (or output of the oscillator system) along 20 h.The frequency ranged from 0 Hz to a maximal value of 3.5 Hz within the time interval from 6:00 to 12:00 h (midday).After 12:00 h, the output frequency of the CTS fell until it reached 0 Hz at 18:00 h.