The circadian rhythm induced by the heterogeneous network structure of the suprachiasmatic nucleus

The circadian rhythms of physiological and behavioral activity are universal in living things, reflecting the period of the earth's axial rotation. In mammals, the circadian rhythms are regulated by a main clock in the suprachiasmatic nucleus (SCN) which is situated above the optic chiasm.^1–4 The SCN is composed of about twenty thousand of nonidentical neuronal oscillators expressing distinct intrinsic periods centered around 24 h.^5,6 These neurons are synchronized to constitute an SCN network with a coherent output.^3,4,7 The period of the SCN network output is dependent on environmental clues. Under constant darkness, the period (free running period) is around 24 h;^8,9 and under an external light-dark cycle, the period of SCN output is precisely equal to the period of the external light-dark cycle.¹⁰ 

The SCN network is heterogeneous.¹¹ It can be divided into two functional distinct subgroups, including a ventrolateral part (VL) and a dorsomedial part (DM).^12–14 The VL contains about 25% SCN neurons which receive photic input from the retina, and the DM comprises the remaining 75% SCN neurons. Both subgroups take part in the regulation of overt circadian rhythms. The periods vary in different regions of the SCN, with the DM running faster than the VL.¹⁵ The SCN neurons are coupled through peptide neurotransmitters which differ between the VL and the DM. In the VL, the neurons express vasoactive intestinal polypeptide (VIP)¹⁶ and arginine vasopressin (AVP) in the DM.¹⁷ The VL and the DM are coupled through Gamma aminobutyric acid (GABA) which is present throughout the SCN.¹⁸ The synapses from the VL to the DM are evidently denser than those from the DM to the VL.^3,19 Accordingly, the coupling is asymmetrical between the VL and the DM, i.e., the VL dominates the DM. Therefore, the minority of the SCN neurons in the VL are coupled with a large number of SCN neurons. Thus far, the detailed network structure of the SCN has not been discovered.

Much theoretical work has been motivated by a desire to understand how the heterogeneity in the SCN affects the circadian rhythm. The heterogeneity in the coupling strength or the neuronal periods not only leads to reduction of the free running period^20–23 but also affects entrainment ability of the SCN to an external cycle.^24–27 In addition, several papers investigated the influence of the SCN network structure on the circadian rhythm.^28,29 Noticeable difference was not found in the synchronization degree of neuronal oscillators or the SCN amplitude between an all-to-all network structure and an small world network structure with only 10% of the total possible links.²⁸ The amplitude of the circadian rhythm in a scale free network is larger than in a completely random network or a regular network.²⁹ 

The Goodwin model has been widely used to describe the SCN network, in which one neuronal oscillator is driven by a transcription-translation negative feedback loop.^{20,21,24–27,30–33} Most previous studies assumed that the SCN is an all-to-all network for simplicity. In the present study, we performed research on the effect of the heterogeneity in the network structure on the circadian rhythm of the SCN based on the Goodwin model. In addition to the all-to-all network, we examined the effects of a Newman-Watts (NW) small world network,³⁴ an Erdös-Rényi (ER) random network,³⁵ and a Barabási-Albert (BA) scale free network,³⁶ which range from low level to high level of heterogeneity in the network structure. We found that the heterogeneous network structure induces the circadian rhythm of the SCN in the case of weak cellular coupling where the rhythm is absent in the all-to-all network structure.

A network is composed of two components, nodes and links. The node here is the neuronal oscillator which is driven by a negative feedback loop consisting of a gene mRNA, its protein, and an inhibitor in the Goodwin model.^{20,21,24–27,30–33} The link is defined as the coupling through transmitters between the neuronal oscillators. The topology of the SCN network is described by a matrix L: if there is an undirected link between the nodes i and j, $L(i,j)=L(j,i)=1$⁠; if there is no link between the nodes i and j, $L(i,j)=L(j,i)=0$⁠. We here let the Goodwin model composed of N neuronal oscillators under constant darkness be read as

where four variables, i.e., the gene mRNA X_i, the protein Y_i, the inhibitor Z_i, and the transmitter V_i, form the neuronal oscillator (node) −i. The transmitter V_i is produced by the gene mRNA X_i, whose mean value from neighbors is represented by F_i. The neuronal oscillators are coupled through the mean value F_i, with the coupling strength (absorbing ability) g. p_i is the degree of node i, which is equal to the number of node i's neighbors.²⁹ Let the average node degree of the network is $m≡1N∑j=1Npi$⁠. If the SCN is an all-to-all network, i.e., $∀1≤i,j≤N, L(i,j)=1$⁠, the current model is reduced to the original Goodwin model.²⁰ With a note, the neuronal oscillator i is linked to itself here. When N = 1, the transmitter V_i is produced and absorbed by the neuron i itself, thus the neuron oscillates with its intrinsic period. In order to model intrinsic periods of the neurons varying from 22 h to 28 h, the factors η_i are introduced,^20,31 i.e., the neuron in the DM with larger η_i shows smaller intrinsic period. The factors η_i satisfy a normal distribution with the mean unity and the standard deviation δ. As reported in Ref. 22, the standard deviation δ plays a role in the SCN network amplitude. When the coupling strength g ∈ [0.77, 0.80], the maximal amplitude is achieved with δ ≥ 0.1. Whereas g < 0.77, the presence of zero amplitude is independent of the value of δ. In the present study, δ is set as 0.1.

In the present study, four types of networks are considered ranging from low level to high level of heterogeneity in the network structure, including the all-to-all network (homogeneous network), the NW network, the ER network, and the BA network (most heterogeneous network). In the BA network, the distribution of node degree satisfies a power law, i.e., few nodes (hubs) are linked with a large number of nodes and most nodes possess small node degrees.³⁶ In the ER network, each pair of nodes is linked with certain possibility q_e.³⁵ In the NW network, the link between two nodes is from two aspects.³⁴ First if the physical distance between two nodes is smaller than a predefined distance, there is a link. Second, each pair of nodes is linked with certain possibility q_n as that in the ER network. In the NW network, the numbers of links from the former and the latter aspects are selected as the same for the present simulations. We compare the effects of network structure among the ER, NW, and BA networks, when the average node degree m is chosen equally for these three networks. The average node degree m is equal to q_eN and 2q_nN for the ER network and the NW network, respectively. For the BA network, the construction of network starts with a very small number of nodes. During the growth of the network, in each time step, a new incoming node is added with M links. These links are connected to the existing nodes with possibility determined by their node degree. For the BA network, the average node degree m is equal to 2M.

The values of other parameters are set as in Ref. 31: α₁ = 6.8355 nM/h, k₁ = 2.7266 nM, n = 5.6645, α₂ = 8.4297 nM/h, k₂ = 0.2910 nM, k₃ = 0.1177/h, α₄ = 1.0841 nM/h, k₄ = 8.1343 nM, k₅ = 0.3352/h, α₆ = 4.6645 nM/h, k₆ = 9.9849/h, k₇ = 0.2282/h, α₈ = 3.5216 nM/h, k₈ = 7.4519 nM, α_c = 6.7924 nM/h, and k_c = 4.8283 nM.

V_i and $Π=1N∑j=1NVj$ are taken into account as markers for the evolution of individual neuronal oscillators and the evolution of the SCN network, respectively. The strength of circadian rhythms is determined by the SCN amplitude and the synchronization degree. When the synchronization degree R is low, the evolution of Π versus time is arrhythmic. The SCN amplitude ρ is defined as

The synchronization degree R between the neuronal oscillators is measured over time as in Refs. 20 and 31 

where $〈…〉$ represents average over time. R = 0 means fully out of synchronization and 1 means perfect synchronization. We used the fourth-order Runge-Kutta method for numerical simulation with time increments of 0.01 h. In order to avoid the influence of transients, the initial 5 000 000 time steps (50 000 h) were neglected. The number of neurons N was set as 1000 for the numerical simulations. The initial conditions for each variable were chosen randomly from a uniform distribution in the range (0–1) in the Goodwin model.

The illustrative examples for the effects of heterogeneity of network topology are presented in Fig. 1. As reported in Refs. 21, 22, and 31, when the coupling strength g = 0.75, the individual neurons go to the steady state, and the amplitude of individual neurons is zero in the all-to-all network (a). In the BA network, the circadian rhythm is induced and the amplitude of the neuronal oscillators is larger than zero, and the phase of the neuronal oscillators becomes slightly different (b). The zero amplitude of individual neurons leads to zero amplitude of the SCN (ρ = 0) in the all-to-all network, while in the BA network the amplitude of the SCN is remarkable (ρ = 0.22) (c).

In order to build the network, the value of the average node degree is required, yet the average node degree of the SCN network has not been found. For the NW, ER, and BA networks, the average node degree was selected from 3 to 51 in the present work. Figure 2 shows the relationship of the parameters for the collective behaviors to the average node degree m in the four kinds of networks when the coupling strength is g = 0.75. Panel (a) shows that in the all-to-all network, the average node degree m is equal to a constant N. Accordingly, we used a dashed line to present the SCN amplitude ρ. In the NW, ER, and BA networks, the amplitude ρ is larger than 0 for certain average node degree m, i.e., m = 3, m = 3 or 5, and m ≥ 3, and the maximal amplitude is achieved when the average node degree is m_c = 3, 5, and 11, respectively. Interestingly, when the average node degree m is larger than m_c, with the increase of m, the amplitude ρ is decreased in the NW, ER, and BA networks, i.e., more links of the network result in smaller amplitude. In the BA network, the amplitude ρ is increased to 0.22 when m ≥ 5, which is visibly larger than the amplitude ρ in the all-to-all, ER, and NW networks. Panel (b) shows that the synchronization degree R between the neuronal oscillators is also larger in the BA network than that in the other networks when m ≥ 5. With a note, when the amplitude ρ is zero, it is impossible to calculate the synchronization degree R.

Further, we investigated the dependence of the parameters for the collective behaviors on the deviation δ in the four kinds of networks when the coupling strength is g = 0.75 (Figs. 2(c) and 2(d)). From panels (a) and (b), we observed that there is no noticeable difference in the amplitude ρ or the synchronization degree R when the average node degree m ≥ 7 in each kind of network. Accordingly, we selected m = 7 for the numerical simulation. In panel (c), the relationship of the amplitude ρ with the deviation δ is not monotonous, and the peak is around δ = 0.1 in the BA network, around 0.14 in the ER network, and around 0.18 in the NW network. At the peak, the amplitude ρ is 0.23 in the BA network, which is evidently larger than 0.09 in the ER network and 0.06 in the NW network. With the increase of deviation δ, the synchronization degree R is reduced in each kind of network (d).

The effects of the average node degree m is examined in the four kinds of networks when the coupling strength g is larger than 0.75 (Fig. 3). With the increase of coupling strength g, the amplitude ρ is larger than zero in the all-to-all network, and the difference of the amplitude ρ among the networks is reduced ((a), (c), and (e)). For each kind of network, when m ≥ 7, the evolution of the amplitude ρ versus m is stable ((a), (c), and (e)), as well as the evolution of the synchronization degree R is stable ((b), (d), and (f)).

In Figs. 2 and 3, we observed that the rhythm is absent in the all-to-all network but the rhythm is induced in the heterogeneous networks with the coupling g = 0.75. When the coupling g ≥ 0.8, in either the all-to-all network or the heterogeneous networks, the SCN shows robust rhythm. In this subsection, an explanation was given for the appearance of the rhythm with strong coupling (g ≥ 0.8) and the emergence of the rhythm induced by the heterogeneity of network topology when the coupling strength is weak (g = 0.75). The absence of a rhythm in the all-to-all network corresponds to the states of the neurons being stable equilibrium points; and the presence of the rhythm in the other types of networks with a certain average node degree corresponds to the states being limit cycles with the coupling strength g = 0.75. In order to examine the existence of stable equilibrium points or limit cycles, there are three steps, i.e., finding the Jacobian matrix of Eq. (1), searching the equilibrium points of Eq. (1), and examining the stability of the equilibrium points.^22,33,37 For simplicity, let the number of the neurons be N = 3, and let “a,” “b,” and “c” represent these three neurons where “a” is in the VL, and “b” and “c” are in the DM. In order to model the neurons in the DM running faster, the factors η_i are set as η_a = 1 − 2θ and η_b = η_c = 1 + θ. Hence, the heterogeneous degree, $δ=(2θ)2+θ2+θ2N=2θ$⁠.

The Jacobian matrix of Eq. (1) with three coupled neurons is shown in Eq. (4)

where

Next, we sought for the equilibrium points. Letting $Ẋ=0, Ẏ=0, Ż=0$⁠, and $V̇=0$⁠, Eq. (1) with N = 3 can be rewritten as

Two network topologies of the SCN are taken into account as shown in Fig. 4, an all-to-all network, i.e., a homogeneous network (a) and a heterogeneous network (b).

After Eq. (5) has been established and the topologies have been defined, the equilibrium points versus the coupling strength g are found by the method of numerical simulation. The relationships of the equilibrium points represented by the marker V with the coupling strength g are shown in Fig. 5. Because neurons b and c are identical for both network structures, we used neuron b to represent neuron c. With the increase of g, the markers V_a and V_b go up in both network structures. There is no evident difference in the neuron a or b between these two network structures.

Next to the equilibrium points having been numerically found, the stability of these equilibrium points is examined. If the equilibrium point is stable, the zero amplitude is to appear; while the unstability of the equilibrium point is accompanied by a limit cycle (oscillation). The Hopf bifurcation locates at the transition from the stable equilibrium point to the unstable equilibrium point. The (un)stability is determined by the maximal value λ_max of the eigenvalues' real parts of the Jacobian matrix of Eq. (1). If λ_max is smaller than 0, the equilibrium point is stable, which corresponds to zero amplitude of the neurons; if λ_max is no smaller than 0, the equilibrium point is unstable, which corresponds to the limit cycle. The relationship of λ_max with g is exhibited in Fig. 5(b). The dashed line λ_max = 0 separates the investigated plane into two regions, i.e., the region of “limit cycle” and the region of “stable equilibrium point.” For the all-to-all network, λ_max is below λ_max = 0 when g < 0.8 and λ_max is above λ_max = 0 when g ≥ 0.8, thus the Hopf bifurcation happens around g = 0.8. Accordingly, the limit cycle (the rhythm) appears when g ≥ 0.8; while for the heterogeneous network, λ_max < 0 when g < 0.75 and λ_max > 0 when g ≥ 0.75. Consequently, the Hopf bifurcation happens around g = 0.75, and the limit cycle (the rhythm) appears when g ≥ 0.75. Therefore, we have proven that the heterogeneous network topology induces the circadian rhythm (limit cycle) when the coupling strength is g = 0.75 and the heterogeneous degree is δ = 0.1.

In the present study, we investigated the role of the heterogeneity of network structure in the circadian rhythm of the SCN reflected by the SCN amplitude and the synchronization degree between the neuronal oscillators. With the increase of the heterogeneity, the all-to-all network (homogeneous network), the NW small world network, the ER random network, and the BA scale free network (most heterogeneous network) are considered. Thus far, the detailed topology of the SCN network has not been revealed. Nevertheless, the topology is suggested to be heterogeneous, since the VL comprising 25% neurons dominates the DM containing the rest 75% neurons.^12,13 Physically, the synapses from the VL to the DM are dense but from the DM to the VL are sparse.^3,19 For reference, in the other regions of brain, the complex topology characterized by small-world properties has been found in the architecture of synaptic connections.^38,39

Most theoretical studies assumed that the SCN is an all-to-all network, and only several works considered the heterogeneous network structure of the SCN. In Ref. 28, the synchronization degree or the SCN amplitude do not differ evidently between the all-to-all network structure and the small world network structure. In Ref. 29, the SCN amplitude was found to be larger in the BA network than that in the regular network and the random network, yet the all-to-all network is not taken into account. In this work, similar results were observed with large coupling strength as in Ref. 29. When the coupling strength is weak, we observed that the circadian rhythm is absent in the all-to-all network where each pair of nodes (neuronal oscillators) is linked. Surprisingly, we found that the robust circadian rhythm is induced by a heterogeneous network structure, such as the BA network which occupies much less links than the all-to-all network. Additionally, the maximal amplitude of the circadian rhythm was achieved with a certain average node degree in each of ER, NW, or BA network, i.e., increasing the number of the total links of each network would lead to the reduction of the SCN amplitude.

Our results provide some insights for the experimental findings. First, in aging,⁴⁰ the circadian rhythm is weakened due to the reduced coupling. The weak circadian rhythm may be strengthened by increasing the heterogeneity of the network structure in aging. The SCN network structure is suggested to be plastic, which is assumed to be altered by the external photopic period.³⁰ We imply that applying proper light signal may improve the heterogeneity of the network structure in order to strengthen the circadian rhythm. Second, our finding provides an alternative explanation for the heterogeneous network structure of the SCN that is beneficial for enhancing or inducing circadian rhythm.

This work was supported by the National Science Foundation of China under Grant Nos. 11505114 and 10975099, the Program for Professor of Special Appointment (Orientational Scholar) at Shanghai Institutions of Higher Learning under Grant Nos. QD2015016 and D-USST02, and the Shanghai project for construction of discipline peaks.