A model of fast synchronized firing patterns in rat brain cells controlling daily rhythms grown on electrode arrays

In this work, we use (i) molecular processes originally modeled by Bhalla [18] describing the VIP-activated cAMP/PKA signal transduction [19], (ii) experimental data demonstrating functioning of cyclic nucleotide-gated (CNG) channels in SCN neurons [20], [21], and (iii) experimental data that characterizes VPAC₂ receptor desensitization and internalization [22]–[25]. Signaling cascades described in the model are shown schematically in Fig. 1. A detailed description of the model can be found in Appendix (Text S1). Equations were integrated using MATLAB ode15s solver and analyzed using MatCont software [26].

We summarized the key data collected in experiments into the following model assumptions:

To study the effect of FOFR on synchronization of circadian rhythms in the SCN neuronal network we incorporated a model of the mammalian circadian clock developed by Leloup and Goldbeter [34] into our model of SCN neuron developed to reproduce ∼30-min oscillations of firing rate. The interaction of VPAC₂ receptors with circadian oscillations was described by introducing an additional term into the Eq.1 of Leloup and Goldbeter's model [34] (see also Eq.1 in [19]) which described dependence of Per gene expression on the CREB activity (see Text S1, Eq. 36, for details):

In this work we considered three models of circadian regulation of SCN neurons firing (Eqs. 49–52 in): Text S1

The degree of the circadian activity synchronization was estimated in two different ways: first, as a standard deviation of either phases of circadian oscillations of firing rate or Per gene product concentration in a population of neurons and second, by introducing the synchronization index (SI) [35]:

The value of SI ranged from 0 (no synchronization) to 1. For each cell, the phaseswere calculated as positions of the center of mass of the cell activity during the successive 24-h periods.

VPAC₂ is a Gs protein-coupled receptor activated by the endogenous peptide VIP. As VPAC₂ activates the adenylate cyclase (AC) - cAMP signaling pathways [36], simulation of VIP application should produce an increase of cAMP concentration in the cytosol.

For low VIP concentrations (<∼1 nM), the time course of simulated cAMP accumulation represented an exponential rise to some steady-state level, while for higher concentrations (>∼1 nM) this accumulation demonstrated a transient peak within 5–10 min after the onset of VIP application followed by a slow decrease to a steady-state level. In both cases, the subsequent washout of VIP returned cAMP concentration to a basal level (Fig. 2A, upper plot). Simultaneously, in both cases simulated VIP treatment resulted in slowly developing internalization of VPAC₂ receptors reaching a plateau (Fig. 2A, lower plot). Washout of VIP resulted in recovery of the VPAC₂ receptors concentration in the plasma membrane to a basal level. A half-time of VPAC₂ receptors internalization at 1000 nM of VIP was ∼4.6 min, and the number of receptors in the plasma membrane reached 85% of its initial value within 1 h after VIP washout. The maximal level of internalization VIP was 76%, and the recovery half-time was ∼20 min in this case. Obviously, a transient peak with a subsequent plateau in cAMP concentration was due to internalization of VPAC₂ receptors, forming in this manner a negative feedback, which stabilized intracellular cAMP concentration. It is noteworthy that the results of these simulations were in agreement with experimental data (Fig. 10 in [25]).

We have also modeled the dependence of AC activity, estimated as a concentration of AC bound with G_s subunit, on the VIP concentration before and after VPAC₂ receptors desensitization and internalization. To do this we have simulated 30 min application of 1 µM of VIP. Then, after simulated 5-min washout, the AC activity in a response to the application of different testing VIP concentrations was calculated and plotted versus these testing VIP concentrations. AC activity was measured in two minutes after VIP application and was expressed as a percent of the highest possible concentration of G_αs-AC complexes. Figure 2B demonstrates that the prolonged activation of VPAC₂ receptors resulted in ∼4-fold decrease in the maximal value of AC activity without a significant shift of AC half-maximal activity. The observed dose-dependence curves (Fig. 2B) were also similar to the experimental data (Fig. 13 in Murthy et al. [25]) indicating the values of experimentally measurable parameters produced by our model are close to experimental data.

Mechanisms producing periodical cAMP oscillations have been studied and modeled in many cellular systems [17], [37]–[39]. In the presented model, oscillations of external VIP concentration as well as oscillations of other parameters controlling its release, e.g. electrical firing rate, are expected to be observed due to the existence of the positive autocrine effect of VIP. Indeed, our simulations demonstrated that the model exhibited oscillations of these and other parameters with a period that fall within the range of experimentally observed values for the period of FOFR in cultured SCN neurons (20–60 min, [15]).

An example of such oscillations for internal cAMP and external VIP concentrations as well as concentration of membrane VPAC₂ receptors and firing rate (model parameters are described in Table 1 in Text S1) is shown in Fig. 3A. A proposed mechanism of autocrine control of observed 30-min oscillations of firing rate suggests that binding of external VIP to VPAC₂ receptor activates G_s protein (Fig. 1F). The activated α-subunit of G_s protein dissociates from its respective βγ-subunits and activates the production of cAMP by adenylyl cyclase (AC) (Fig. 1B). cAMP activates cation CNG channels, which depolarize the SCN neurons (Fig. 1A). Depolarization of model neuron evokes action-potential (AP) firing that, in turn, induces VIP secretion (Fig. 1A). This sequence provided a positive feedback loop for the mechanism of 30-min oscillations. Simultaneously, cascades of events interrupting the positive feedback loop were incorporated into our model. First, four cAMP molecules sequentially bind to each protein kinase A (PKA) receptor subunits and release two activated catalytic subunits (Fig. 1D). Then, the PKA activates cAMP phosphodiesterase (PDE), which transforms cAMP to AMP (Fig. 1C). Second, the same PKA evokes desensitization and internalization of VPAC₂ receptors via phosphorylation of G protein-coupled receptor kinase (GRK) (Fig. 1F). Recovery from desensitization and internalization, which are the slowest processes in the model, determines the 30-min period of firing rate oscillations.

Analysis of the oscillatory mechanism in the system of 47 equations is complicated, but it can be simplified taking into account that the recovery from desensitization and internalization of VPAC₂ receptors are the slowest processes in the system. In the limit case of very slow rate of receptors internalization, evolution of the system on the VIP - VPAC₂ plane takes place along an orbit lying close to the VIP nullcline (i.e. a line of equilibrium VIP concentration for a fixed number of active VPAC₂ receptors). This nullcline has a characteristic N-shape (a red line in Fig. 3B), which reflects the fact that if the number of VPAC₂ receptors in the system is fixed then the system has a region of bistability. This bistability arises as a property of VIP-cAMP positive feedback loop: at low VIP and cAMP levels this feedback does not work - both variables evolved rather independently and their concentrations remained low; at higher levels strong positive interaction raised concentrations of both molecules closely to some high stationary levels restricted by the saturation rate of cAMP production. These interactions are illustrated in Fig. 3C, where nullclines of the system with a fixed number of VPAC₂ receptors are shown in the cAMP – VIP plane. These nullclines have two or three points of intersection depending on the number of VPAC₂ receptor molecules in the plasma membrane. In the particular case shown in Fig. 3C there are 3 points of intersection and thus, the system has three fixed points two of which are stable (shown by filled circles). In the ideal case of slow receptor internalization, relaxation oscillations would arise in the system when the VPAC₂ nullcline (a green line in Fig. 3B) intersects the VIP nullcline (a red line in Fig 3B) in its descending phase (i.e., along the region shown in yellow in Fig. 3B). In this case, the system evolves repeatedly along the two positive slope regions of the VIP nullcline with rapid transitions between them (blue line in Fig. 3B). For experimentally observed rates of receptor desensitization and internalization, the set of intersection points between the nullclines that yield stable oscillations was estimated by continuous shifting of the VPAC₂ nullcline (green line in Fig. 3D) to the right and monitoring the oscillation amplitude. Shifting of the VPAC₂ nullcline was implemented by changing initial concentration of GRK in the model. As it can be seen from Fig. 3D, the obtained set of intersection points (yellow region in Fig. 3D) mostly coincided with the descending phase of the VIP nullcline. The stable limit cycle also arised in this case and it is shown in Fig. 3D (blue line).

To assess the robustness of oscillations with respect to the variation of model parameters, we performed two series of computational experiments. Firstly, we systematically varied initial concentrations of all molecules in order to define the value range in which oscillations of cAMP concentration with amplitude higher than 10% of average cAMP level would be observed for longer than 15 hours of model evolution. Secondly, we defined a range of all kinetic rates and parameters added to Hao's et al model [19] in which oscillations with the above mentioned characteristics occured. The varied parameters were: the rates for activation of GRK by PKA_active (k₆₁), dephosphorylation of GRK (k₆₂), phosphorylation of VPAC₂ by GRK (k₆₃), dephosphorylation of VPAC2 (k₆₄), phosphorylated VPAC₂ receptor internalization (k₄₂), insertion of new VPAC2 receptors (k₇₀). The other varied parameters were: maximal conductivity of CNG channels (g_{CNG, 0}), steady state VIP concentration (VIP₀), threshold for AP generation (I_threshold) and the ratio of the rates of binding with Gs of phosphorylated to dephosphorylated VPAC₂ receptors (DE).

Results of these computational experiments are presented in Fig. 4A. Only for 10 out of 47 molecules described in the model variations of their initial concentrations in the range (where is the default initial value) resulted in disappearance of oscillations. These were AC, Receptor-G-protein complex (R.GDP.G_asbg), GDP bound Gs trimeric complex (GDP.G_asbg), cAMP phosphodiesterase (cAMP-PDE), VPAC₂ Receptors (R), PKA heterotetramer with two molecules of regulatory subunits and two molecules of catalytic subunits (R2C2), PKA heterotetramer reforming intermediate (R₂C), Cytoplasmic PKA inhibitor (PKA-inhibitor_C) and ATP. Parameter ranges (blue bars; Fig. 4A) in which oscillations were preserved during variation of parameters are shown in a logarithmic scale for the 10 most critical initial concentrations described above and for all added parameters. The oscillations were most sensitive to the following parameters (listed here in order of their influence on the oscillations): concentrations of PDE, ATP and G_βγ and conductance of CNG channels. For these parameters, the ratio of upper to lower boundaries of their values at which oscillations were still observed was 2.4, 2.6, 3.0 and 3.15, respectively. At the same time, this ratio was larger than 19 for other parameters. Thus, FOFR could be observed in our model for a wide range of model parameters, indicating model robustness.

To study in detail the parameter dependence of system dynamics, we have calculated positions of fixed points and stable limit cycles in the system while changing 6 most critical model parameters. Figure 4 shows dependencies of equilibrium positions (a red line) on concentrations of the key molecular players – VIP₀ (Fig. 4B), CNG-channels (Fig. 4C), AC (Fig. 4D), VPAC₂ (Fig. 4F), and GRK (Fig. 4G), as well as characteristic time of recovery from internalization for VPAC₂ receptors, τ_D (Fig. 4E; τ_D was set to an inverse rate of recovery from internalization, k₇₀). Fig. 4 also shows amplitudes of stable limit cycle (green lines). Dependencies of the oscillation periods on each model parameter are shown in the insets. All the bifurcation diagrams have a characteristic structure: when VIP-induced depolarization exceeds a firing threshold, oscillations appear abruptly having low frequency and high relative amplitude. Subsequently, they disappear in a fashion observed in the canard explosion scenario [40]: a fast transition from relaxation oscillations to relatively high frequency low amplitude oscillations. It can be concluded from these bifurcation diagrams that the oscillations with properties that fell within the experimentally observed range could be observed for a wide range of the key model parameters. The parameters differ in their potency to influence the period of oscillations. As it can be predicted, the period of oscillations is determined by the time course of receptors internalization (Fig. 4E), but it is also strongly correlated with the number of VPAC₂ receptors (Fig. 4F) and the concentration of GRK (Fig. 4G).

An important property of FOFR observed in dispersed and cultured SCN neurons, which we expected to have a physiological significance, was synchronization of electrical activity between cells located in distant parts of the 1-mm recording area of a multielectrode array dish [15]. In many cases, such synchronization of FOFR was observed without obvious correlation of AP firing in studied neurons. This fact allowed us to exclude a principal role of excitatory synaptic connections in this phenomenon and to suggest that synchronization is provided by diffusible VIP in the intercellular milieu [14], [41]. In order to investigate whether a heterogeneous population of model neurons could oscillate in synchrony, we developed a model of a neuronal network with substantial random scattering of initial parameters between individual neurons. Initially, a short period of neurons behavior was simulated supposing that there is no VIP exchange between them. In this case each cell sensed only autocrine VIP. After this short period VIP exchange between the neurons within the network was allowed, i.e. VIP concentration sensed by all the cells in the network was set to the same value so that the cells could oscillate in synchrony.

A typical example of model experiment that was conducted with a network of interacting SCN neurons is shown in Fig. 5. The model included 50 cells that oscillated independently until VIP exchange was enabled. Each initial concentration in the model was randomly selected from an interval , and each rate constant was also randomly selected from an interval , where and are default values for initial concentrations and kinetic rates for each molecule and chemical transition in the model (see Text S1). Oscillations were initially observed in 70% of the cells (35 of 50). A period of oscillations varied from 14 min to 36 min (mean±SD = 23.0±5.3 min) and their amplitudes varied from 0.1 Hz to 9 Hz (4.8±2.3 Hz). Examples of firing rate changes for silent and spontaneously active cells are shown on the left from a vertical red line in Fig. 5A, and the respective distribution of oscillation phases is shown in Fig. 5B. After 20 hours of the independent evolution of each cell (D = 0), VIP exchange was enabled (D = 0.5 after a vertical red line in Fig. 5A). As a result, oscillations of different neurons in the network became highly synchronized and a distribution of their amplitudes gradually became narrower (Fig. 5A, to the right of the vertical red line, and C). The oscillation period varied in different cells in a very narrow range from 16.7 min to 16.9 min (16.77±0.06 min) and the range of the respective oscillation amplitudes was from 1.2 Hz to 3.1 Hz (1.75±0.45 Hz). Distribution of oscillation phases was markedly narrowed as a result of neuronal network synchronization (standard deviation 24 sec).

As was shown in [9], [42], [43], less than 20% of SCN cells exhibit circadian firing rate oscillations in low-density cultures, while almost all cells oscillate when cultured at high density. Therefore, we further addressed a similar question of whether the 20% of cells that exhibited FOFR evoked by autocrine VIP activation without VIP exchange could evoke persistent synchronous oscillations in other cells when intercellular exchange of VIP was enabled. In this set of experiments a network of 20 cells was modeled, with 4 cells having default parameters (see Text S1) and other 16 cells having CNG conductance twice lower than default, while all other parameters were set to their default values. Initially there was no intercellular VIP exchange, and only these 4 cells generated oscillations of firing rate (Fig. 6A, D = 0). When VIP exchange was enabled (D = 0.5), all 20 cells started to oscillate synchronously (Fig. 6B) with a mean period of 28.19±0.01 min, amplitude of 4.20±0.03 Hz, phase SD of 47 sec. Thus, our modeling indicates that VIP exchange in the SCN network having 20% of intrinsically oscillating neurons results in the synchronous oscillations of entire neuronal network. This result was robust with respect to both CNG-channel conductance and moderate random perturbations of parameters. Figures 6C, D demonstrate that 4 cells having default parameters were capable to induce synchronous oscillations in the other 16 cells having 20% of default CNG conductivity while all other parameters of the model were selected for each of these 16 cells randomly and uniformly from an interval where is a default parameter value. It is known [44] that the amplitude of network oscillations inversely depends on the strength of coupling between oscillators (so called oscillation death phenomenon). Correspondingly, in our model increase of coupling strength resulted in expected rise of synchrony while the amplitude of FOFR decreased. For weakly coupled oscillators (D = 0.3, Fig. 6C, E) the phase SD was 49 sec, amplitude mean±SD – 2.5±1.5 Hz, period of oscillations – 26.6 min. For strongly coupled oscillators (D = 6, Fig. 6D, F), the degree of synchrony was characterized by phase SD of 18 sec, amplitude mean±SD of 1.0±0.1 Hz, period of oscillations – 27.77±0.03 min.

In accordance with our conception, existence of synchronous FOFR in a cultured network of SCN neurons suggests that their firing activity can be promptly and strongly coupled to the activation of VPAC₂ receptors. Such coupling may be important for the synchronization of the circadian activity of SCN network, since the phase of circadian electrical activity for each SCN neuron could be determined in this case by the VIP level in extracellular milieu rather than exclusively from the gene expression pathway [14], [35]. In order to study the importance of this coupling in a computational experiment, we added the model of the mammalian circadian clock developed by Leloup and Goldbeter [34] to our model of the SCN neuronal network with FOFR. The interaction of VPAC₂ receptors with a circadian genes network was modeled as described in Materials and Methods (Fig. 1). We considered three models of circadian regulation of neuronal firing (see ‘Modeling the circadian synchronization in the SCN network with FOFR’ in Materials and Methods).

Two series of experiments with 10-cell networks were conducted to demonstrate that VPAC₂/CNG-channel coupling in the SCN network could help to achieve rapid and precise synchronization of electrical firing activity even if synchronization of circadian genes expression is slow and imperfect.

In the first series of experiments, all cells had default parameters for both circadian clock and FOFR models. Circadian oscillations in each cell had a uniformly distributed phase shift. Initially, cells oscillated independently (D = 0) up to the 91 hour, when coupling between oscillators was introduced (D = 0.5). The typical results of such experiments are shown in Fig. 7, where each row represents a different model of circadian regulation of firing activity.

The existence of FOFR in cultured SCN neuronal network suggests the existence of a strong and fast coupling between receptors for biologically active substance(s) and action potential generation ([15], Figure S1). This coupling previously was shown in direct electrophysiological experiments with SCN neurons in slice preparations [20]. In these experiments, application of VIP induced depolarization of neuronal membrane within several seconds. VIP is considered as a main paracrine signaling molecule in the SCN and it is known that its interaction with VPAC2 receptors on the neuronal surface induces elevation of cAMP concentration. Thus, the model developed in our work could be considered as a member of the class of biochemical oscillator models based on the cAMP signaling. Other examples of such models include oscillations of cAMP level in Dictyostelium discoideum amoebae [37], GnRH neurons of hypothalamus [17] and pancreatic insulin-releasing β cells [38], [39]. All these models have a similar structure with a positive feedback loop which consists of the G-protein coupled receptors for some paracrine substance, AC and some link between cAMP level and the rate of exocytosis of this substance. Another important component of these models is a negative feedback loop which breaks this positive feedback at high cAMP levels. This loop could be organized in a variety of ways. In GnRH neurons, this loop is organized via switching from Gs to Gi activity at high hormone levels [17], ]; in D. discoideum it is organized through blocking of AC by PKA [37] and desensitization of receptors for cAMP [45], [46] and in pancreatic cells it is organized through several mechanisms, including periodic activation of PDE by Ca2+-calmodulin kinase or PKA and calcium-dependent inhibition of AC [39]. In our model we have hypothesized that a negative feedback in a single SCN neuron is based on the well-known phenomenon of VPAC₂ receptor desensitization and internalization [22]–[25]. We have shown that reversible internalization of receptors with properties that fall within the experimentally observed range is enough for the emergence of oscillations. However, we could not exclude that some other mechanisms not yet well characterized for SCN neurons, e.g., G_i and G_q subunit activation at high VIP levels, could also contribute to feedback generation [47], [48].

There were several reasons to introduce CNG channels in our model of FOFR. First, cAMP-signaling plays an important role in SCN pacemaking [39], [42]. Moreover, CNG-like channels were observed in acutely isolated SCN neurons [15], [49] and their activation by db-cGMP, membrane-permeable derivate of cGMP, resulted in the long-lasting firing [15]. Besides, an application of diltiazem, a blocker of CNG-channels, produced a reversible inhibition of firing in SCN neurons (N.I. Kononenko, unpublished data). Finally, cAMP with corresponding CNG-channels are key players of pulsatile electrical firing in a well-studied mechanism of 30-min oscillations in GnRH neurons [50]. In order to verify incorporation of CNG-channels in our model, one could study effects of their putative blockers on experimentally observed FOFR in MED.

We have also investigated the robustness of proposed mechanism of FOFR. The obtained results have shown that oscillations were highly robust with respect to the variations of all model parameters (Fig. 4). Therefore, the occurrence of FOFR in our model is a direct consequence of the model properties rather than of any specific choice of parameters.

It has been shown earlier that separate neurons in distant parts of the 1-mm MED recording area exhibited synchronized FOFR ([15], Figure S1). Our model successfully mimicked synchronization of 30-min oscillations of firing rate in the SCN network containing 20–50 cells (Fig. 5, 6). Notably, reliable synchronized FOFR in separate cells were observed when model parameters of each cell were randomly chosen in a broad range. Even if most of the cells were initially silent, as it is observed in dispersed and cultured SCN neurons [9], a small group of oscillating model neurons (∼20%) were capable to involve the silent cells into FOFR (Fig. 5, 6).

The fast oscillations of electrical activity similar to modeled in this work have not been observed in the intact SCN yet. The reason for this could be the differences in tonic activity state and connectivity density between SCN cell culture, SCN in vitro and intact SCN, similar to the dependence of the amplitudes and other properties of globally synchronized oscillations of firing activity on the same parameters observed in cultures of neocortical neurons [51], [52]. Another possible explanation is that FOFR only occurs within a small subnetwork of the SCN. Nevertheless, the existence of these oscillations in the SCN neuronal cultures raises the possibility that the membrane potential of a particular SCN neuron during circadian oscillations can be controlled by the level of external neuromodulator(s) through coupling of their respective receptors and specific membrane conductance, rather than by the internal signal from the circadian genetic oscillator of the neuron, as it was assumed in other models [35], [53]–[55]. An intercellular exchange of these neuromediator(s) can efficiently amplify and average influences of the genetic oscillations of a large neuronal populations thus far inducing synchronous circadian oscillations of membrane potential in each neuron. In line with this assumption, recent experiments with SCN neurons cultured in MED have shown that synchronous FOFR arise together with synchronization of circadian peaks of electrical firing ([15], Fig. B, C, D1 and D2 in Figure S1). The synchronization of circadian electrical peaks in separate SCN cells is now considered to be regulated by Per gene expression only [14], but it is plausible to suggest that FOFR contributes to this process. Although this interpretation of experimental results ([15], Fig. D1 in Figure S1) should be done with caution, we attempted to model this putative process and found that FOFR could influence synchronization of circadian firing activity peaks.

We have shown that the Leloup-­Goldbeter model of circadian oscillators [34], with intercellular VIP exchange without FOFR, produces synchronization of circadian firing within 2.5–5 days after enabling VIP intercellular exchange (Fig. 7A1 and Fig. A in Figure S2). We have shown that the Leloup-Goldbeter model of circadian oscillators [8], with intercellular VIP exchange without FOFR, produces synchronization of circadian firing within 2.5–5 days after enabling VIP intercellular exchange (Fig. 7A1 and Fig. A in Figure S2). These results are close to those obtained in To's et al model [14] in which Per gene expression regulates VIP secretion. Incorporation of FOFR in the model including Leloup and Goldbeter circadian oscillator resulted in practically immediate synchronization of circadian peaks of electrical firing (Fig. 7B1, Fig. C in Figure S2). This effect can be easily observed in the model of SCN network with FOFR (Fig. 7 B1,B2, Fig. Fig. A, B in Figure S2) and still exists in the networks without FOFR although only for the narrow range of parameters (e.g., in the network near saddle-node bifurcation, Fig. 7 C1,C2, see ‘Modifications of the default parameter set’ and ‘Sensitivity of the synchronization of the circadian oscillations of firing rate to the model parameters’ in Text S1). Nevertheless, it may be assumed that even a small subnetwork of cells with FOFR can improve synchronization of circadian firing activity in the whole network. This assumption is in line with prediction that a small group of strongly connected circadian oscillators could underlie temporal precision of the electrical activity of the whole SCN network [43], [50]. Alternatively, synchronization can be achieved via slow signaling from the membrane to the genetic network in the neuronal nucleus and back to the membrane [35], [53]–[55]. It is natural to suggest that in this case the precision and rapidness of synchronization highly depends on the heterogeneity of signaling in these feedback loops and in circadian genetic networks in the population of neurons.

Several predictions based upon the FOFR mechanism's contribution to synchronization of circadian firing activity peaks can be experimentally tested: (i) reactivation of the paracrine signaling in the SCN network with partially desynchronized genetic oscillations should immediately improve synchronization; (ii) block of paracrine signaling [29] should reduce synchronization of circadian oscillations of firing activity in separate SCN neurons; (iii) synchronization of membrane potential oscillations should be stronger in comparison with synchronization of genetic networks. Experimental testing of these suggestions and further elaboration of the presented model of synchronization within SCN network are necessary for complete understanding of circadian clock mechanisms.