article: The suprachiasmatic nucleus (SCN) of the hypothalamus is a multicellular system that drives daily rhythms in mammalian behavior and physiology. Although the gene regulatory network that produces daily oscillations within individual neurons is well characterized, less is known about the electrophysiology of the SCN cells and how firing rate correlates with circadian gene expression. We developed a firing rate code model to incorporate known electrophysiological properties of SCN pacemaker cells, including circadian dependent changes in membrane voltage and ion conductances. Calcium dynamics were included in the model as the putative link between electrical firing and gene expression. Individual ion currents exhibited oscillatory patterns matching experimental data both in current levels and phase relationships. VIP and GABA neurotransmitters, which encode synaptic signals across the SCN, were found to play critical roles in daily oscillations of membrane excitability and gene expression. Blocking various mechanisms of intracellular calcium accumulation by simulated pharmacological agents (nimodipine, IP3- and ryanodine-blockers) reproduced experimentally observed trends in firing rate dynamics and core-clock gene transcription. The intracellular calcium concentration was shown to regulate diverse circadian processes such as firing frequency, gene expression and system periodicity. The model predicted a direct relationship between firing frequency and gene expression amplitudes, demonstrated the importance of intracellular pathways for single cell behavior and provided a novel multiscale framework which captured characteristics of the SCN at both the electrophysiological and gene regulatory levels.
  In mammals many physiological and behavioral responses are subject to internal time-keeping mechanisms or biological clocks. Daily rhythms are generated by an internal, self-sustained oscillator located in the suprachiasmatic nucleus (SCN) of the hypothalamus. The SCN produces autonomous 24h cycles in gene expression and firing frequency from the synchronization of multiple individual oscillatory signals across the network [1]. The single cell gene regulatory mechanism involves a number of interlocking positive and negative feedback loops in which the Period (Per) gene occupies the central position [2]. The circadian modulation of neural firing affects a number of electrophysiological properties of the cell membrane which also fluctuate over the course of the day [3]. In vitro studies of SCN slices and cultures have demonstrated diurnal modulation of neural firing [4], resting potential [5] and membrane resistance [6], as well as daily oscillations in a number of ionic currents that include the fast delayed rectifier potassium [7], L-type calcium [6] and the large-conductance Ca2+-activated potassium [8], [9] channels. Although individual SCN neurons contain molecular feedback loops that drive such rhythms, membrane excitability and synaptic transmission also play significant roles in generating daily oscillations. Experimental studies in Drosophila have demonstrated the dependence of core clock oscillations on electrical activity, as electrical silencing resulted in abolishment of circadian oscillations of the free-running molecular clock [10]. In mammalian organisms a direct association between membrane excitability and core-clock rhythms has also been reported in multiple studies, providing evidence for a positive correlation between Per gene transcription and neural spike frequency output [11]–[13]. For example, activation of GABAA receptors via muscimol enhanced inhibitory postsynaptic currents (IPSCs) leads to lower firing rates [14], [15] and suppression of Per1 mRNA [12]. Another example involves mice deficient in vasoactive intenstinal peptide (VIP) receptors known to display lower amplitude oscillations of both core clock genes [16] and neural firing [17]. The mechanisms by which the single cells produce synchronized rhythms in neural firing, gene expression and neuropeptide secretion are postulated to involve intracellular second messengers [18]. A candidate second messenger that regulates diverse cellular processes is intracellular calcium. Cytosolic calcium is known to oscillate over the course of the day preceding rhythms in multiple-unit-activity (MUA) recordings by a mean phase of ∼4. 5 hr [4]. Variations in intracellular calcium concentrations have been demonstrated to induce Per1 gene expression by activating the Ca2+/calmodulin dependent kinase, which in turn phosphorylates the cAMP-response-element binding (CREB) protein [19]. Reduced Ca2+concentrations have been shown to abolish daily Per1 mRNA oscillations in SCN slices [20]. Cytosolic Ca2+rhythms also affect neural firing frequency, as dampening of Ca2+ oscillations via blockade of calcium release from ryanodine-sensitive pools results in decreased firing activity [4], [6]. To our knowledge, detailed cell models with molecular descriptions of gene expression and neural firing coupled by intracellular signaling pathways are not currently available for any circadian system. In a recent study (Sim and Forger 2007), a Hodgkin-Huxley type model of SCN neurons was developed and shown to reproduce a significant amount of experimentally observed electrophysiological behavior on a millisecond timescale. While this study facilitated the formulation of our electrophysiology model by providing guidelines for the mathematical representation of a number of relevant ionic currents, our modeling study was distinct due to its focus on the circadian timescale. In addition to incorporating single-cell electrophysiological properties, our model has accounted for circadian rhythmicity by coupling electrophysiology to daily oscillations in core-clock gene expression and calcium dynamics. The objective of the present study was to model couplings between the circadian gene-regulatory pathway, cellular electrophysiology and cytosolic calcium dynamics to evaluate the role of extracellular synaptic stimuli on firing rate behavior over a circadian timescale. The role of distinct intracellular pathways as well as the directionality of information flow along the network nodes was evaluated by analyzing single cell model behavior following the introduction of various external stimuli. Calcium dynamics, adapted from a published model [21], included the contributions of IP3- and ryanodine stores as well as the flux of Ca2+ in and out of the cell membrane. Our model has demonstrated the dependence of membrane excitability on synaptic input conveyed by VIP and GABA, and predicted reduced neural firing and Per mRNA oscillation amplitudes as well as shorter circadian periods with reduced cytosolic Ca2+ concentrations. These results suggest a dual effect of signaling pathways instigated by VIP and calcium that potentially operate as coupling agents between the gene regulatory network and the electrophysiology of SCN neurons.
  In this work we developed a firing rate code model to capture the circadian fluctuations of relevant ion channels as well as 24 hour trends in core-clock gene expression. Individual currents (IK, INa, ICa, IKCa, Iex and Iinhib) were assumed to interact with the circadian gene regulatory network via signaling pathways that included VIP and Ca2+ contributions (Fig. 1). Our model produced daily oscillations in constant darkness with 23. 6 hour periodicity in the ionic and synaptic currents (Figs. 2A–2F), intracellular calcium concentration (Fig. 2G), Per mRNA expression (Fig. 2H) and neural firing rate (Fig. 2I). These rhythmic profiles constitute the nominal output of our model and will be referred to as the “control. ” To calculate the phase relationships of individual rhythmic signals we regarded the cytosolic calcium peak at CT 1. 5 [4] as our reference point and computed the phase differences between the Ca2+ peak and the peaks of the various rhythmic profiles. This calibration method was used to produce a well tuned model, where circadian components exhibited rhythmic behavior generally matching experimental data both in their relative phase relationships and oscillation peaks as summarized in Table 1. The circadian trend of the sodium current remains unknown [22], but our model predicted circadian oscillations of INa with a peak during the subjective night at CT 14. 3 (Fig. 2B). The effects on calcium dynamics on circadian behavior was first studied by clamping the membrane voltage at a hyperpolarized level while maintaining constant calcium concentrations. As observed experimentally [17], the model neuron produced arrhythmic behavior (results not shown). Next we blocked various mechanisms of intracellular calcium accumulation to investigate their effects on firing rate and rhythmic behavior. In an effort to reproduce experimentally observed trends we specifically eliminated IP3-stores, ryanodine stores and L-type Ca2+ currents. Ikeda et. al (2003) [4] showed that IP3-stores did not contribute to calcium oscillations and action potential dynamics. In our model, IP3-blockade was simulated by zeroing the rate of calcium release from InsP3-sensitive stores (v1) and observed to reproduce this data (results not shown), demonstrating circadian dynamics of Per gene expression were not affected by cytosolic calcium stores. Ryanodine stores, responsible for Ca2+ release into the cytosol, are known to play an important role in the 24h oscillations of the intracellular calcium concentration and action potential frequency [4]. Ryanodine blockade was implemented by zeroing the term that accounts for calcium release from the stores into the cytosol (v3). Elimination of this term resulted in a decrease in the peak of the neural firing rate by 13% compared to the control during the subjective day (Figs. 3A, 3B), consistent with Ikeda et. al (2003) [4] who reported a 22±8% reduction. Ryanodine blockade had a minimal effect on the firing frequency trough during the subjective night in agreement with experimental findings [4] (Figs. 3A, 3B). Decreases of 11% and 45% in the intracellular calcium trough (subjective night) and peak (subjective day), respectively, were observed compared to the control (Figs. 3B, 3C) consistent with Ikeda et. al (2003) [4]. The reduction in intracellular calcium levels had an effect on ICa, which was predicted to decrease by 20% during the subjective day. Furthermore, our simulations produced a 17% decrease in the Per mRNA amplitude (Fig. 3B, 3D). This trend is consistent with Lundkvist et. al (2005) [20] who demonstrated abolishment of Per gene expression rhythms with decreasing intracellular calcium concentrations. L-type Ca2+ currents experimentally blocked via nimodipine application have been shown to affect firing rates, but not intracellular calcium levels, within a single cell [4], [22]. The effects of nimodipine were implemented by setting the L-type Ca2+ conductance (gCa) to zero. Experimental studies involving nimodipine application were carried out over a maximum period of 5 hours [4], whereas our simulations involved constitutive nimodipine application over the course of the day. Under these conditions the model predicted a slight period decrease of the core-oscillator period to 22. 2 h. Simulations of nimodipine application produced a 45% decrease in the firing frequency peak compared to the control during the subjective day, while minimal effects were observed on the minimum firing rate during the subjective night (Figs. 3A, 3E) consistent with experimental studies [4], [22], [23]. Decrease in the intracellular calcium concentration peak by 10% was also observed (Figs. 3C, 3E) in agreement with experiments reporting a 4±6% reduction [4]. The Ca2+-activated potassium current (IKCa) decreased by 23% and 38% during the subjective day and night, respectively (Fig. 3E), in agreement with the literature where 30–50% reductions have been reported [22]. The model produced a reduction of the Per mRNA amplitude by 30% (Figs. 3D, 3E), similar to that reported by Lundkvist et. al (2005) [20]. We simulated autocrine response of the GABA neurotransmitter to investigate the role of inhibitory postsynaptic currents (IPSCs) on single cell neural firing and circadian behavior. Our initial objective was to simulate the effects of incrementally decreasing GABA concentrations by imposing step reductions in the mean cytosolic GABA level (Eq. 20). The effect of GABA reduction in our system during the subjective day and night was evaluated by computing the percent changes in firing frequency peak and trough versus the control. IPSC levels exhibited dose dependent reductions (Fig. 4A) leading to increased membrane excitability and therefore increased neural firing rate (Fig. 4B, 4E) as GABA concentrations decreased, in agreement with the literature [15]. Consistent with experimental data by Gribkoff et. al (2003) [14] complete blockade of GABA was seen to increase the firing rate peak and trough by 11% and 57%, respectively (Fig. 4E). Additional simulations involved incremental increase of the GABA concentration binding to the cell surface, implemented by including an additive term in the mathematical expression of GABA (Eq. 20). Our model produced a gradual increase in IPSC levels accompanied by a dose-dependent decrease in neural firing as increasing GABA concentrations were applied (Figs. 4C, 4D), consistent with experimental findings [14]. GABA application was shown to have a minimal effect on the neural firing rate during the subjective day (less than 2%) while significantly decreasing the firing frequency by a maximum of 28% during the subjective night (Fig 4E). This dependence of cellular response on the circadian time of GABA administration has also been demonstrated in experimental studies [14]. Our simulations produced reduced firing frequencies during the subjective night comparable to experimental data. The smaller responses produced during the subjective day, however, did not agree with available data [14]. We simulated autocrine response of the VIP neurotransmitter to investigate the effects of the VIP signaling pathway on single cell behavior. Initially, we zeroed the VIP concentration responsible for the GABA oscillations (Eq. 20) and CREB activation (Eq. 30). Complete VIP blockade reduced firing rate amplitudes by 66% (Figs. 5A, 5C) in agreement with Brown et. al (2007) who reported a 52% average reduction [17]. Our simulations showed an acute decrease in Per mRNA levels that oscillated with much lower amplitudes compared to the control (73% decrease; Fig. 5B, 5C), as observed experimentally by Maywood et. al (2006) [16]. Because VIP blockade affects GABA release a 57% reduction in IPSC amplitude was also observed (Fig. 5C), consistent with data from Itri et. al (2004) [24]. VIP elimination resulted in a period decrease to 22. 2 h, consistent with Brown et. al (2007) [17] who measured a 22. 9±1. 9 h period across VIP−/− SCN populations. Additional simulations involved increase of the VIP concentration binding to the cell surface to investigate the effects of constitutive VIP application throughout the circadian cycle. Simulation of VIP application was implemented by adding 1nM to the VIP concentration responsible for GABA release (Eq. 20) and CREB activation (Eq. 30), leading to saturation of the VPAC2 receptors. To our knowledge, experimental data on constant VIP administration throughout a 24 hour period are not currently available. Our simulations showed a 25% increase in the circadian amplitude of the firing frequency (Fig 5A, 5D). Because Per gene expression is the final target of the VIP signaling cascade [25], the imposed VPAC2 receptor saturation resulted in increased Per mRNA amplitudes of approximately 54% (Figs. 5B, 5D). IPSC amplitudes were observed to increase by 110% (Figs. 5D), comparable to experimental studies [26]. The simultaneous increase in neural spiking and IPSCs can be attributed to the dominant effects of Per gene dynamics within the network. Constitutive VIP application decreased the predicted period to 22. 5h. These model predictions can be tested experimentally by applying constant VIP concentration throughout a 24 hour period while conducting bioluminescence recordings to measure Per gene activation and utilizing multielectrode arrays on highly dispersed SCN cultures to measure firing activity of single cells. We varied the intracellular calcium concentration to investigate its effects on the rhythmic output of the circadian clock. Changes in effective Ca2+ levels were achieved by scaling the output of Eq. 12, responsible for the circadian evolution of intracellular calcium, by multiplying with a scaling factor ranging from 0. 5 to 1. 5. Hence mean levels of calcium were varied by ±50% of their nominal value. Because our model was constructed under the assumption that cytosolic calcium instigates a signaling cascade with Per gene transcription as the final product [19] incrementally increasing Ca2+ concentrations had a positive effect on Per mRNA amplitudes (Figs. 6A, 6C). A similar trend was observed for neural firing, as increasing intracellular calcium increased firing frequency amplitudes (Figs. 6B, 6C). The calcium concentration also affected the periodicity of the model system. Increased Ca2+ levels produced longer periods of the core oscillator, reaching a maximum of 25. 6h for a 50% Ca2+ increase (Fig. 6B). Our model predictions can be compared with experimental data on SCN explants, which show abolished mean Per mRNA rhythms when averaged over the entire population as a function of increasing concentrations of Ca2+ buffer [20]. Our simulations suggest that the observed elimination of collective Per gene expression rhythm across the population can be attributed to reductions in the amplitude of individual Per mRNA signals accompanied by a period decrease on the single cell level. This hypothesis requires further experimental studies for validation. Incrementally increasing current levels were applied on the cell membrane of our neuron model to investigate the effects of extracellular electrical stimuli on single cell behavior. Current application was implemented by including an additive term in the mathematical expression of I* (Eq. 3). As expected, a positive, linear correlation of firing rate with inward current levels (I) was observed (Fig 7A). Because VIP is released as a function of firing rate (Eq. 29), VIP concentrations were also seen to increase with electrical stimuli (Fig. 7A). The direct relationship between mean firing rate and mean VIP concentration over the course of a circadian cycle is shown in Fig 7B. Increasing electrical stimuli did not significantly contribute to mean intracellular calcium levels (results not shown). Neural firing was predicted to correlate with core-clock gene transcriptional activity as demonstrated by Quintero et. al (2003) [11]. Mean Per mRNA levels were predicted to increase as a function of the mean neuronal spiking frequency and ultimately obtained their maximum after a threshold in firing rate had been reached (Fig 7B). Our model demonstrated a positive relationship between electrical firing and core-clock gene activity consistent with the literature [11], [12]. Simulations of increasing electrical stimuli suggest correlations between Per mRNA levels and the firing frequency potentially mediated via a VIP-instigated signaling pathway. Our model postulates an underlying intracellular network where VIP, released due to elevated neuronal spiking, binds on the cell surface initiating a signaling mechanism that leads to Period gene activation.
  We developed a multiscale mathematical model to investigate the association between the circadian gene regulatory pathway, electrophysiology and cytosolic calcium and to evaluate SCN single cell behavior over a circadian time scale. Initially, we investigated the effects of calcium dynamics on cell behavior by simulating the blockade of L-type Ca2+ channels, IP3- and ryanodine stores. The model successfully replicated a number of experimental observations. The effect of synaptic input on membrane excitability was explored by assuming autocrine response of the VIP and GABA neurotransmitters. Simulations of variable GABA concentrations and VIP blockade were conducted to compare the model output with experimental data and test model validity. Increasing the GABA concentration was shown to have an inhibitory effect on neural firing, matching experimental data. Several experimental studies have reported excitatory responses in a subset of SCN neurons depending on the circadian timing of GABA administration [27], [28]. These effects have not been included in the present model but will be considered in our future work. Simulations of VIP blockade produced decreased amplitudes in firing frequency, Per mRNA and IPSCs compared to the control, all in agreement with the literature. We further tested the effects of constitutive VIP application, for which experimental data are not currently available. Our model predicted increased amplitudes in Per mRNA, firing rate and IPSCs compared to the control as VPAC2 receptors became saturated. Our model postulates that calcium plays an important role in the coordination of neural firing and core clock gene expression. To test this hypothesis, we gradually altered the intracellular Ca2+ concentration and determined the effect on single-cell rhythmic behavior. Amplitude reductions in Per mRNA and neural firing accompanied by a period decrease were observed as the cytosolic Ca2+ concentration was gradually reduced. Experimental data showing the role of cytosolic Ca2+ on circadian behavior are currently available only for SCN explants [20]. These data showed elimination of the mean Per mRNA rhythm, averaged over the entire population, as cytosolic Ca2+ was gradually buffered. Thus, our simulations suggest that reduction in the amplitude of individual Per mRNA signals accompanied by a period decrease and deregulation of the phase relationships between the various circadian components may be manifested experimentally as the elimination of the collective Per gene expression rhythm across the population. This hypothesis can be tested experimentally by adding a buffer (e. g. the intracellular chelator BAPTA-AM) to reduce cytosolic Ca2+ while conducting bioluminescence recordings to measure Per gene activity of a single cell and utilizing multielectrode arrays on highly dispersed SCN cultures to measure neural firing activity. Our model demonstrated a positive correlation of core-clock gene expression with the neural spike frequency consistent with the literature [11], [12]. Incrementally increasing electrical stimuli applied on the cell membrane of our neuron model affected the mean levels of the firing rate, Per gene expression and VIP release. Furthermore positive relationships between VIP release and core-clock gene transcription with firing rate were observed. We hypothesized an underlying intracellular network where VIP, released due to elevated neuronal spiking, binds on the cell surface initiating a signaling cascade that leads to Per gene activation. Our single cell model was shown to replicate a number of experimentally observed trends, as well as to provide predictions concerning intracellular couplings. One of the key features of the model is the ability to test effects of blockers, neurotransmitters or extracellular stimuli for prolonged periods of time, which can be challenging in an experimental setup. Long-term application of ryanodine blockers or intracellular calcium buffers (BAPTA-AM), for example, is known to disrupt intracellular Ca2+ from physiological levels [4], [20]. A number of vital cell processes depend on calcium activity, rendering the interpretation of such experiments in terms of a single effect on the circadian network difficult. The model may therefore provide predictions that assist in the development of carefully designed experiments to test these hypotheses.
  The core oscillator utilized in our model originates from a previous study [29] and consists of 16 ordinary differential equations in time that describe intertwined negative and positive regulatory transcriptional loops. Transcription of the Per and Cry genes is activated by a heterodimer formed from the CLOCK and BMAL1 proteins. This activation is rhythmically suppressed and reestablished by a complex of the PER and CRY proteins, which blocks the activity of CLOCK/BMAL1 dimer and negatively autoregulates transcription of the Per and Cry genes. Our model did not include the loop involving Rev-Erbα since it was not required for sustained circadian oscillations. Nominal parameters values were mostly obtained from the original reference, with the exception of the parameter k1 that determines the transport rate of the PER/CRY complex from the cytosol to the nucleus, KAP involved in Per activation due to elevated BMAL1 concentrations and vsp0 that denotes the basal value of Per transcription rate. These values were modified as part of the model tuning process (see Table 2). In this work we utilized a modified integrate-and-fire model [30] that takes into account the contributions of the relevant ion channels and includes the effects of extracellular synaptic stimuli. Following the methodology of Liu and Wang (2001) [31] we constructed a SCN membrane model that included sodium (INa), potassium (IK), calcium (ICa) and calcium-activated potassium (IKCa) currents, as proposed by Brown et al (2007) [3]. The contributions of inhibitory and excitatory input signals, known to influence membrane excitability, were also incorporated in our model. Individual ionic and synaptic currents (Ir) were modeled as: (1) where V represents the membrane voltage, gr is the conductance and Er is the reversal potential of current r. Our single SCN model neuron was described by a modified integrate-and-fire model [31]–[34]: (2) where Cm denotes the membrane capacitance and gex, gGABA, gNa, gCa, gΚ, gΚCa, and gL are the conductances for the excitatory, inhibitory, sodium, calcium, potassium, calcium-activated potassium and leakage currents, respectively, and Eex, EGABA, ENa, ECa, EΚ, and EL are the corresponding reversal potentials. Activation and inactivation variables associated with the relevant conductances were not included in this model (Eq. 2) as they evolve on the millisecond timescale rather than the circadian timescale considered in the study. By defining: (3) (4) (5) Eq. (2) can be reformulated to yield [30]: (6) The integrate-and-fire model does not produce the form and shape of an action potential. Rather neural spikes can be characterized by a “firing time” tf, defined by the criterion: (7) where θ is the firing threshold. Immediately after tf the potential is reset to a new value Vreset<θ. To calculate the trajectory of the membrane potential after the occurrence of a spike at time tf, Eq. 6 was integrated with the initial condition V (t (1) ) = Vreset. Because time variations in I*, R*and τm were much faster than the circadian timescale they were considered constant for each simulated 10 minute time step. Therefore integration of Eq. 6 yields: (8) The membrane potential described by Eq. 8 approaches the asymptotic value V (∞) = R*I* as t→ ∞. Therefore since R*I*>θ the membrane potential reaches the threshold θ at time t (2): (9) The time interval t (2) −t (1) constitutes the firing period T′. Thus the firing rate (fr) can be calculated as the inverse of T′ [30], yielding the firing rate code model used in this study: (10) The firing rate (eq. 10) fluctuates due to circadian variations in conductances and reversal potentials of the various currents, which are included within the terms I*, R* and τm (eqs. 3–5). Membrane properties such as the resting potential and resistance display sustained circadian rhythms [5], [6]. Our model included oscillations of the membrane resting potential (Vrest) by utilizing a modified version of the Goldman-Hodgkin-Katz equation derived by Piek (1975) [46] that takes into account both monovalent ions and divalent ions, such as Ca2+: (23) (24) (25) (26) where R denotes the gas constant, F is the Faraday constant, PCa, PK, PNa, and PCl are the membrane permeabilities of Ca2+, K+, Na+ and Cl−, respectively, Kin and Nain represent the K+ and Na+ concentrations within the cytosol, whereas Kex, Caex and Naex are the K+, Ca2+ and Na+ concentrations in the extracellular space. Values for PCa, PNa, PCl, Kex, Clex and Naex were chosen to match experimental measurements from the literature [39], [47], whereas Kin and Nain were computed by inversion of the Nernst equation. PK values were modeled to vary over the course of the day in agreement with Kuhlman et. al. [5], who demonstrated circadian rhythmicity in K+ currents underlying the membrane potential oscillations: (27) where vPK is the maximum value, KPK is the saturation constant and npk is the cooperativity coefficient of the PK oscillations. Eq. 27 was not intended to imply a mechanistic relationship between PK and BC. We modeled the membrane potential to oscillate over the course of the day and to peak during the subjective day. The membrane resistance (R) oscillated in a constant phase relationship with the resting potential and peaked during the subjective day as shown by experimental studies [5], [6]: (28) where VR represents the maximum value and KR the saturation constant of the membrane resistance oscillations. Core clock gene transcription displays self-sustained circadian rhythms that are likely modulated via VIP/VPAC2 activation [25] and fluctuations in intracellular calcium dynamics [19] (Fig. 1). We utilized a revised signaling transduction mechanism from our previous study [48] to capture the effects of these two components on gene regulation. In this study, VIP oscillations were assumed to depend on the neural spike frequency as well as the rate governing the depletion of the neurotransmitter from the synaptic cleft: (29) where vVIP is the maximum rate, KVIP the saturation constant and nVIP the cooperative coefficient of VIP release, while kdVIP denotes the rate constant and ndVIP the cooperativity coefficient of VIP depletion. The VIP concentration binding on the membrane surface of our model neuron was assumed to be equal to the VIP concentration released by the neuron (Eq. 20), consistent with an autocrine response. The intracellular calcium concentration also displayed rhythmic variations over the course of the day (Eqs. 12–13). The transduction mechanism involving the VPAC2 receptor and Ca2+ likely includes the activation of protein kinases, which in turn phosphorylate CREB, leading to core clock gene activation. Protein kinase activity was modeled as: (30) where vk is the rate of kinase activity, VMK and Vβ denote the maximum rates of Ca2+- and VIP-induced protein kinase activation, respectively, and KMK and Kβ represent the saturation constants of Ca2+- and VIP-induced protein kinase activation, respectively. Circadian fluctuations in the phosphorylated CREB fraction, as well as the dynamics of the Per gene activation, are described in detail in our original study [48]. Parameter values altered from the original reference are displayed in Table 2. The complete single cell model was formulated within MATLAB (The MathWorks, Natick, MA) and consisted of twenty ordinary differential equations (ODEs). Sixteen ODEs described the circadian evolution of the gene transcriptional loop (for details refer to [29]), two ODEs described intracellular calcium rhythms (Eqs. 12–13 modified from [21]), and the remaining two ODEs described the VIP concentration (Eq. 29) and the phosphorylated CREB concentration (for details refer to [48]). The model was integrated numerically using the differential-algebraic equation solver ode23 with a 10 minute time step to ensure accurate solutions with reasonable computational cost. Nominal parameter values utilized in our model are listed in Table 2, with parameters directly obtained from the literature accompanied by the corresponding reference. The tuning of the remaining parameters is discussed in detail below. The initial conditions utilized in simulations of the 20 ordinary differential equations characterizing our model system are listed in Table 3. These values were chosen to produce individual rhythmic profiles that oscillated within a reasonable range, consistent with experimental data.
summary: Circadian rhythms are ∼24 hour cycles in biochemical, physiological and behavioral processes observed in a diverse range of organisms including Cyanobacteria, Neurospora, Drosophila, mice and humans. In mammals, the dominant circadian pacemaker that drives daily rhythms is located in the suprachiasmatic nucleus (SCN) of the hypothalamus. The SCN is composed of a highly connected network of ∼20,000 neurons. Within each individual SCN neuron core clock genes and proteins interact through intertwined regulatory loops to generate circadian oscillations on the molecular level. These neurons express daily rhythmicity in their firing frequency and other electrophysiological properties. The mechanisms by which the core clock produces synchronized rhythms in neural firing and gene expression are postulated to involve intracellular calcium, a second messenger that regulates many cellular processes. The interaction between the various clock components however remains unknown. In this paper, we present a single cell model that incorporates the circadian gene regulatory pathway, cellular electrophysiological properties, and cytosolic calcium dynamics. Our results suggest a possible system architecture that accounts for the robustness of the circadian clock at the single cell level. Our simulations predict a dual role for intracellular pathways instigated by intracellular calcium and VIP: maintaining the periodicity and amplitude of the core clock genes as well as the firing frequency oscillations.
keywords: computational biology/systems biology
  computational biology/computational neuroscience
title: A Multiscale Model to Investigate Circadian Rhythmicity of Pacemaker Neurons in the Suprachiasmatic Nucleus