A design principle of spindle oscillations in mammalian sleep.

Neural oscillations are mainly regulated by molecular mechanisms and network connectivity of neurons. Large-scale simulations of neuronal networks have driven the population-level understanding of neural oscillations. However, cell-intrinsic mechanisms, especially a design principle, of neural oscillations remain largely elusive. Herein, we developed a minimal, Hodgkin-Huxley-type model of groups of neurons to investigate molecular mechanisms underlying spindle oscillation, which is synchronized oscillatory activity predominantly observed during mammalian sleep. We discovered that slowly inactivating potassium channels played an essential role in characterizing the firing pattern. The detailed analysis of the minimal model revealed that leak sodium and potassium channels, which controlled passive properties of the fast variable (i.e., membrane potential), competitively regulated the base value and time constant of the slow variable (i.e., cytosolic calcium concentration). Consequently, we propose a theoretical design principle of spindle oscillations that may explain intracellular mechanisms behind the flexible control over oscillation density and calcium setpoint.

Introduction

Neural oscillations are generated by periodic electrophysiological activity in groups of neurons and are involved in cognitive functions such as memory, attention, and consciousness. Generally, the mechanisms behind these oscillations include (1) network connectivity along with synaptic dynamics (Bartos et al., 2007; Brunel and Wang 2003), and (2) electrical properties of neuronal membranes mainly determined by a combination of ion channels (Hutcheon and Yarom 2000; Llinas, 1988). The involvement of neuronal circuits in generating oscillations has been investigated via classical lesioning experiments or optogenetic manipulations, which have provided profound insights into many aspects of oscillations in the brain (Buzsaki, 2015; Buzsaki and Wang, 2012; Fernandez and Lüthi, 2020; Neske 2015). On the other hand, molecular mechanisms underlying neural oscillations and their functional engagement remain largely elusive. Understanding the molecular configuration of neural oscillations and the properties derived from each component could also be an important step toward controlling various properties of neural oscillations (e.g., amplitude, frequency, and duration) flexibly and independently and investigating their involvement in cognitive functions. Nevertheless, cell-intrinsic molecular mechanisms, especially a design principle, underlying neural oscillations have drawn less attention.

Some previous studies have attempted to determine molecular mechanisms underlying neural oscillations using conductance-based neuron models connected with each other (Bazhenov et al., 2002; Compte et al., 2003; Destexhe et al., 1994, 1999), but the models are heavily dependent on the provided parameter values and thus lacked generalizability. On the other hand, recent studies developing simple computational models of groups of neurons suggest that such models could be a promising way to determine basic molecular configurations of neural oscillations (Alonso and Marder 2019; Rasmussen et al., 2017; Tatsuki et al., 2016; Yoshida et al., 2018). Specifically, a previous study developed an averaged-neuron (AN) model, which is a simplified neuron-network model based on the mean-field approximation of a population of neurons, to investigate the molecular mechanism of slow waves (Tatsuki et al., 2016). Slow waves, or 0.5–4 Hz deflections in electroencephalography (EEG) or local field potential (LFP), are among the most characteristic oscillations during non-rapid eye movement (NREM) sleep in mammals. The neuronal activity underlying slow waves is represented as a slow-wave sleep firing pattern (SWS firing pattern) comprising depolarized “up” states with sustained neuronal firings and hyperpolarized silent “down” states. The AN model contained a constellation of ion channels that were commonly expressed in the cortex and transmitted its output to equivalent neurons (i.e., its output returned to itself as an input). The model recapitulated the SWS firing pattern, identified important channels or ionotropic receptors for the firing pattern, and consequently proposed a new mechanism for the sleep/wake cycle regulation (Shi and Ueda 2018; Tatsuki et al., 2016). Furthermore, in the following study, the authors developed a simplified AN (SAN) model, which suggests a basic configuration of the SWS firing pattern. This simplified model, combined with detailed mathematical analysis enabled by the simplification, suggested that leak potassium channels could play a role in sustaining the SWS firing pattern (Yoshida et al., 2018). Therefore, these studies on the SWS firing pattern demonstrate that a simple but comprehensive neuronal model can reveal molecular mechanisms and design principles underlying neural oscillations.

Sleep spindles, or 10–15 Hz bursts of activity in EEG or LFP, are also typical oscillatory events during NREM sleep in mammals. During sleep spindle generation, the membrane potential of thalamic neurons exhibits waxing and waning and intermittent firings without up or down states (Destexhe et al., 1993; Leresche et al., 1991). This firing pattern is represented as a sleep spindle firing pattern (SS firing pattern) in this study. In vitro electrophysiological studies have demonstrated that the intrinsic properties of neurons in the thalamic reticular nucleus (TRN), as well as the reciprocal connections between thalamocortical relay cells and TRN neurons, play a crucial role in producing sleep spindles (Bal et al., 1995; Krosigk et al., 1993; Steriade et al., 1993). However, recent in vivo electrophysiological experiments or fMRI studies have suggested the existence of non-thalamic origins of sleep spindles outside the thalamus (Bandarabadi et al., 2020; Halassa et al., 2014; Schabus et al., 2007; Xu et al., 2021). Moreover, SS-like firing patterns have also been reported in neural oscillations other than sleep spindles (Amir et al., 2002; Del Negro et al., 1998). These studies underscore the importance of investigating molecular mechanisms, in addition to specific neuron-network connectivity, underlying the SS firing pattern to understand spindle oscillations. Furthermore, sleep spindles have been implicated to be involved in memory consolidation and other cognitive functions (Antony et al., 2019; Eschenko et al., 2006; Fernandez and Lüthi, 2020; Gais et al., 2002). It has been reported that mean spindle density (count per unit time) during NREM sleep is flexibly controlled depending on contexts and positively correlated with cognitive functions (Chatburn et al., 2013; Gais et al., 2002; Reynolds et al., 2018). However, it remains elusive how spindle density is flexibly controlled.

Here, to elucidate molecular mechanisms of spindle oscillation and derive its functional implications, we investigated the SS firing pattern using the AN model. We demonstrated that the AN model could sufficiently recapitulate the SS firing pattern. To extract a basic configuration of the SS firing pattern, we subsequently developed a minimal model for the SS firing pattern based on the AN model. Based on the configuration of the minimal model and analysis of current, we discovered that slowly inactivating potassium channels played an important role in characterizing the SS firing pattern. Furthermore, dynamical system analysis revealed that, generally, the subcritical Andronov-Hopf bifurcation and fold limit cycle bifurcation underlay the SS firing pattern, which explained the waxing and waning subthreshold oscillations of spindle oscillations. Based on this mathematical structure of the SS firing pattern, we discovered that conductance of leak sodium and potassium channels, which controlled passive properties of the fast variable (i.e., membrane potential), competitively controlled the base value and time constant of the slow variable (i.e., cytosolic calcium concentration). Overall, we proposed a theoretical design principle of spindle oscillation by leveraging a minimal model of groups of neurons, and discovered that the balance of background inward and outward currents could flexibly control spindle density along with intracellular calcium concentration.

Results

The averaged-neuron model recapitulates the SS firing pattern

To elucidate the basic intracellular mechanism of spindle oscillation, we constructed an AN model consisting of 13 components (nine types of ion channels, three types of ionotropic receptors, and one type of exchangers/pumps) with 10 variables (V, , , , , , , , , and [Ca2+]) described by ordinary differential equations (Figure 1A) (Tatsuki et al., 2016). We tested whether the model recapitulated the SS firing pattern, which is putatively the intrinsic firing pattern in groups of neurons during sleep spindle generation. The AN model recapitulated the SS firing pattern with the features corresponding to the previous studies such as spikes with large afterhyperpolarizations (AHPs) and silent phases without up or down states (Figure 1B). Close inspection revealed that there were the waxing and waning oscillations before and after the bursting phases, respectively, and the rhythmic activity of the membrane potential (i.e., alternating bursting and silent phases) was accompanied by [Ca2+] oscillation (Figure 1C). The slower oscillation of [Ca2+] was mediated by the influx of Ca2+ through voltage-gated Ca2+ channels and NMDA receptors and by exiting through Ca2+ pumps in this model (see STAR Methods). This means that the SS firing pattern consists of a fast spiking subsystem and a slower subsystem of [Ca2+]. Therefore, the AN model reproduced the waveform of the SS firing pattern, and its oscillation density corresponded to the frequency of the intracellular calcium concentration. However, we need to investigate the general properties of the SS firing pattern in the AN model, which are not specific to a certain parameter set to obtain valid insights for spindle oscillation.

Figure 1

The configuration of the averaged-neuron (AN) model and the reduced-AN (RAN) model for the sleep spindle (SS) firing pattern

(A) Schematic illustration of the AN model. It is composed of extrinsic components in dendrites and intrinsic components in soma.

(B) Membrane potential in the slow-wave sleep (SWS) (top) and SS (bottom) firing patterns. While the intermittent firing pattern is shared between the two firing patterns, the active up and silent down states are not observed in the SS firing pattern.

(C) Detailed membrane potential (top) and [Ca2+] (bottom) dynamics in the SS firing pattern. The waxing and waning subthreshold oscillations and afterhyperpolarizations (AHPs) are unique to the SS firing pattern.

(D) SD of the standardized V (left) and the average of the standardized [Ca2+] (right) of the 1,112 parameter sets collected by random parameter search for the AN model.

(E) Distinct distributions of several parameters between the SS and the SWS parameter sets.

(F) Factor loadings of each component for PC 1 in the AN model.

(G) Distributions of the parameters whose factor loadings are highlighted in (F). The components colored yellow or magenta are implied to be characterizing the SS firing pattern.

(H) Proportions of parameter sets that bifurcate to other firing patterns after KO of each component. The components colored magenta are implied to be indispensable for generating the SS firing pattern.

(I) Table shows channels, receptors, and exchangers/pumps included in each model.

(J) Hit rates for the SS firing pattern in each model. Details are listed in Table S6.

(K) Schematic illustration of the RAN model

To this end, we collected parameter sets that recapitulated the firing pattern from more than parameter sets which were generated randomly from exponential distributions based on automatic and visual inspections based on the waveforms (see STAR Methods). The success rate of the parameter search was 2.0 . This success rate suggests that the SS firing pattern is sufficiently stable and valid in the AN model, comparing the success rate for the SWS firing pattern (3.9 ). Moreover, all the collected SS firing patterns were homogeneous in terms of the dynamics of the membrane potential and [Ca2+] oscillation (Figure 1D). Hence, we demonstrated that the AN model validly recapitulated the SS firing pattern and collected over a thousand SS parameter sets for further analysis.

In terms of intermittent firings accompanied by [Ca2+] oscillation, the SS and SWS firing patterns resemble each other. Thus, to characterize the SS firing pattern, we compared it with the SWS firing pattern. We first compared the distributions of each parameter in the collected parameter sets between the two firing patterns (Figures 1E and S1A). In particular, the conductance of leak channels (), voltage-gated Na+ channels (), Hodgkin-Huxley-type K+ channels (), slowly inactivating channels (), persistent channels (), and the time constant of pumps () displayed distributions distinct from those of the SWS firing pattern (Figure 1E). This result indicates that these components can play a role in characterizing SS firing patterns. To examine the contribution of each component to the differentiation between the SS and SWS firing patterns, we conducted a primary component analysis (PCA) and factor analysis. The parameter sets for the SS firing pattern and those for the SWS firing pattern were clustered separately along the primary component (PC) 1 axis (Figure 1F). The factor loadings of PC 1 revealed that , , , , (yellow bars), and especially (magenta bar) contributed to PC 1 (Figure 1G). Moreover, , , and of the SS firing pattern tended to be larger than those of the SWS firing pattern, while Na+ and tended to be smaller (Figure 1G). Therefore, these components, especially , make the SS parameter sets distinct from those of the SWS. We subsequently compared the oscillation parameters (average length of bursting/silent phases, average number of spikes per bursting phase, and average oscillation density) between the SS and SWS firing patterns. The SS firing pattern was characterized by fewer spikes per burst, longer bursting phases, shorter silent phases, and larger oscillation density than the SWS firing pattern (Figure S1B). These results imply that the difference between the two firing patterns is not just a difference in the fast subsystem but also comes from the fundamental mechanisms underlying them. Based on these analyses of the AN model, we discovered that the AN model validly recapitulated the SS firing pattern which was distinct from the SWS firing pattern, but it is still difficult to describe and interpret the detailed molecular mechanism because the model is biologically and mathematically too complex to conduct more elaborate analyses.

Construction of the reduced AN model specific to the SS firing pattern

The complexity of the AN model still makes it difficult to identify the important components for generating sleep spindles and to perform detailed mathematical analyses. Therefore, a minimum model, particularly one comprising fewer components with fewer than three variables, is needed. Indeed, a previous study constructed SAN model that contained six components and three variables (V, , and [Ca2+]), and revealed the important role of leak K+ channels in the SWS firing pattern through bifurcation analysis and dynamical system analysis (Yoshida et al., 2018). Therefore, we simplified the AN model for the SS firing pattern by eliminating redundant components. To screen components that could not be replaced, we first knocked out each component in each parameter set. We subsequently calculated the percentage of parameter sets in which the SS firing pattern turned into any other firing pattern after knocking out one of its components. This revealed that knocking out leak channels, slowly inactivating K+ channels, calcium-dependent K+ channels, persistent Na+ channels, and Ca2+ pumps (magenta bars) abolished the SS firing pattern in almost all of the parameter sets (Figure 1H), indicating that these components are irreplaceable in the SS firing pattern. Consequently, we constructed a model with only these five components (Model 1) and conducted a random parameter search (Figure 1I). We generated more than parameter sets, none of which recapitulated the SS firing pattern (Figure 1J). It was assumed that Model 1 did not contain any channels or receptors that mediate Ca2+ influx, failing to generate the intermittency of the SS firing pattern without slow [Ca2+] rhythm. Therefore, we added either NMDA receptors () or voltage-gated Ca2+ channels () to Model 1 and named them Model 2 and Model 3, respectively. Voltage-gated Ca2+ channels mediated Ca2+ influx in the models, and knocking out these components abolished the SS firing pattern in more than half of the parameter sets in the screening (Figures 1H and 1I). The success rates of the random parameter search for Models 2 and 3 were 6.4 and 5.0 , respectively (Figure 1J). Although Model 2 was slightly more efficient in producing the SS firing pattern, it contained five variables (V, , , , and [Ca2+]), whereas Model 3 contained three variables (V, , and [Ca2+]) (Figure 1I). Therefore, we selected Model 3 as a simplified AN model for the SS firing pattern and named it the reduced averaged-neuron (RAN) model (Figure 1K). The RAN model recapitulated the SS firing pattern as in the AN model (Figures S1C-S1E). Notably, the composition of the ion channels in the RAN model was similar to that in the single thalamocortical neuron model which recapitulated spindle oscillations (Destexhe et al., 1993), and the subtypes ion channels that were reported to be important for generating sleep spindles were included in the model (Astori et al., 2011; Espinosa et al., 2008; Kim and McCormick 1998; Wimmer et al., 2012). This means that the simple molecular configuration of the RAN model covers previously reported models for spindle oscillation. Therefore, we constructed the RAN model as a simplified model for the SS firing pattern (Figure 1K).

Molecular mechanism of rhythmic firing without up or down states based on the RAN model

To investigate the SS firing pattern using the RAN model, we first focused on the molecular configuration of the model. The RAN model consisted of Ca2+-related channels ( and ), leak channels (), Na+ channels (), and K+ channels (). Notably, the only difference between the RAN and SAN (i.e., the simplified AN model for the SWS firing pattern) models was voltage-dependent, delayed rectifier K+ channels: slowly inactivating K+ channels in the RAN model and Hodgkin-Huxley-type K+ channels in the SAN model (Figures 1I and S2A). Therefore, it was speculated that both models could recapitulate both firing patterns, but they could not. This suggests that the different kinetics of the delayed rectifier K+ channels (i.e., activation at lower voltages and/or slower time constants) was critical for distinguishing the SS firing pattern from the SWS firing pattern (Figure S2A). Subsequently, to investigate the contribution of each conductance for producing the SS firing pattern, we compared the parameter distributions of the SS and SWS parameter sets. and of the SS firing pattern in the RAN model were generally larger than those of the SWS firing pattern in the SAN model (Figure S2B). These distributions are consistent with those in the AN model (Figure 1E). Therefore, the kinetics and conductance of slowly inactivating K+ channels and conductances of leak channels and persistent Na+ channels are key for characterizing the SS firing pattern.

To understand how the unique kinetics or conductance contributed to the SS firing pattern, we examined inward and outward currents through each channel. Because leak current is mainly mediated by inward Na+ currents and outward K+ currents through different types of channels, we divided leak channels into Na+ leak channels and K+ leak channels (see STAR Methods). We chose a representative parameter set based on PCA and kernel density estimation (KDE) and calculated the proportion of each channel's inward or outward currents over time. The current analysis revealed that (1) spikes and AHPs without sustaining depolarization were generated by fast, transient and subsequent delayed-rectifying , respectively, and that (2) outward currents during silent phases were dominated by and which resulted in silent phases without profound hyperpolarization (Figure 2A). These observations were supported by other parameter sets (Figures 2B and 2C) and unique to the SS firing pattern (Figures S2C and S2D). Therefore, the characteristic channel kinetics and conductances contributed to the intermittent firings without up or down states, which is one of the major features of the SS firing pattern. Moreover, the only calcium-dependent current, , mediated the transitions between bursting and silent phases in synchronization with [Ca2+] oscillation, despite its comparably small occupancy in the outward current (Figures 2A–2C). Consequently, these results partly explain the waveform of the SS firing pattern composed of the current through each channel. However, it remains unclear how small changes in drive the dynamic transitions between phases, and how subthreshold waxing and waning oscillations appear.

Figure 2

Small changes in mediates the transition between phases in the SS firing pattern

(A) Membrane potential proportions of each channel's inward and outward current of the representative parameter set in the RAN model over time.

(B) Proportion of current through leak K+ channel () (left), slowly inactivating K+ channel () (middle), and calcium-dependent K+ channel () (right) of 1,166 parameter sets obtained by the random parameter searches in the RAN model (top) and their averages with standard deviations (bottom).

(C) Distributions of the proportion of outward current in bursting (left) and silent (right) phases of the SS and SWS firing patterns in the RAN and SAN models, respectively.

(D) A trajectory over time (red line) and nullclines of V and of the representative parameter set of the SS firing pattern in the phase space. The coiled part of the trajectory represents bursting phases, while the linear part represents silent phases.

(E) Phase planes and stream plot of the representative parameter set on the V- plane when [Ca2+] is low (67 μM) and high (90 μM). Nullclines determine the streams in the phase plane and therefore the trajectory over time. Each concentration corresponds to the transition from a silent to a bursting phase and a bursting phase to a silent phase, respectively.

(F) Fixed points (that is, intersections of V and nullclines) on the V-[Ca2+] plane and the trajectory of the SS firing pattern (left), and a bifurcation diagram representing the linear stability of the fixed points with the arrows indicating the movement of the trajectory (right). The transition from a silent to a bursting phase occurs via the subcritical Andronov-Hopf bifurcation. On the other hand, the transition from a bursting phase to a silent phase occurs theoretically via the fold (saddle-node) limit cycle bifurcation

Mathematical mechanism of subthreshold waxing and waning oscillation based on the RAN model

To determine the hidden properties of the SS firing pattern, we conducted a dynamical system analysis. We first drew the trajectory of the SS firing pattern over time in the phase space in which each point represented a state of the system at a certain moment. The trajectory for the representative parameter set revealed that the coiled and straight parts corresponded to the bursting and silent phases, respectively (Figure 2D). To uncover how this system evolved over time along the trajectory, we plotted nullclines of V and in the phase space and conducted a linear stability analysis for the set of fixed points (i.e., intersections of the nullclines) (Figure 2D). At lower [Ca2+] ([Ca2+] 67 μM), the trajectory diverged from the unstable equilibrium and converged to the stable limit cycle, whereas at higher [Ca2+] ([Ca2+] 90 μM), the trajectory converged to the stable equilibrium point (Figure 2E). To illustrate how the stability of the fixed points controls the entire system, we plotted the equilibrium points with their stability on the V-[Ca2+] plane (see STAR Methods). This bifurcation diagram revealed that the transition from the silent phase to the subsequent bursting phase was driven by the transition from the stable focus to the unstable focus, meaning the subcritical Andronov-Hopf bifurcation. On the other hand, the transition from a bursting phase to the subsequent silent phase was driven by the transition from the stable limit cycle to the stable focus, implicating the fold (saddle-node) limit cycle bifurcation (Figure 2F). This type of oscillation is mathematically called “sub-Hopf/fold cycle” oscillation or “elliptic” oscillation (Izhikevich 2000; Rinzel 1987). On the other hand, when we conducted the same analysis on the SWS firing pattern, the transition from a silent phase to the following bursting phase and that from a bursting phase to the following silent phase corresponded to the fold bifurcation and the saddle homoclinic orbit bifurcation respectively, and this type of oscillation is called “fold/homoclinic” oscillation or “square-wave” oscillation (Figures S2E–S2G) (Izhikevich 2000; Rinzel 1987). In the sub-Hopf/fold cycle oscillation, the transitions between the phases occur gradually via delayed stability loss in the subcritical Andronov-Hopf bifurcation and slow convergence to the stable focus in the fold limit cycle bifurcation (Figure 2F, right). This graduality in the phase transitions resulted in subthreshold waxing and waning oscillations (Figure 2F, left). Moreover, the trajectory during the bursting phases needed to produce large AHPs for each spike to trace the limit cycle around the unstable focus (Figure 2F). Consequently, firing in the SS firing pattern depended on and , which generated bursting phases without up states. Furthermore, small changes in [Ca2+] and can sufficiently drive the phase transitions because these bifurcations are caused by tiny shifts of the V nullcline (Figures 2D–2F). Therefore, the mathematical mechanism behind the SS firing pattern explains the subthreshold waxing and waning of spindle oscillations and dynamical shifts between phases with little change in in the SS firing pattern.

Bidirectional control of oscillation density and calcium setpoint by the balance between background inward and outward currents based on the RAN model

Although it has been reported that oscillation density of the sub-Hopf/fold cycle oscillation was flexibly controlled in previous studies, which used highly simplified computational models (Burke et al., 2012; Ying and Qin-Sheng 2011), it has not been investigated whether Hodgkin-Huxley-based models could reproduce the control of the oscillatory rhythm by a specific parameter, and moreover, which molecular components were responsible. To test this, we examined the oscillation density of the SS firing pattern (i.e., the density of spindle oscillations) in response to changes in each conductance or time constant in the RAN model. It was demonstrated that in most cases and were negatively correlated with oscillation density (denoted as “decreasing”), while and were positively correlated (denoted as “increasing”) (Figures 3A–3C, S3C, and S3D). On the other hand, there were some parameter sets in which oscillation density was controlled in the opposite direction in response to the parameter shifts (Figures 3B, S3A, S3D, and S3E). These results confirmed that oscillation density could be controlled by the channel conductance, which is not directly involved in the generation slower subsystem of the SS firing pattern (i.e., [Ca2+] oscillation). Furthermore, the mean intracellular calcium concentration (i.e., setpoint) was controlled along with oscillation density: the larger and (or the smaller and ) were, the lower [Ca2+] was, in the major parameter sets (Figures 3A and S3C). In the minor parameter sets, on the other hand, the opposite phenomena were observed (Figures S3A and S3E). Notably, these minor subgroups were related to higher [Ca2+] compared to the major subgroups (Figures S3B and S3F). Therefore, we discovered that the conductance of certain ion channels, particularly those involved in determining the base value of the fast variable V, controlled the frequency and setpoint of the calcium oscillation described by the slow variable [Ca2+] in two ways: the major one in which oscillation density and calcium setpoint were positively correlated and the minor one in which they were negatively correlated. To investigate the controlling mechanism in more detail, we examined oscillation density in response to changes in two out of four parameters (, , , and ). The representative parameter sets were chosen from the common part of the parameter sets of the two channels for the major and minor subgroups (Figure S3G). This analysis revealed that the conductance of channels carrying current in the same direction (e.g., and ) cooperatively controlled oscillation density and the calcium setpoint, while the conductance of channels carrying current in the opposite directions (e.g., and ) competitively controlled them (Figures 3D and S3H). These results suggest that the control over oscillation density and the calcium setpoint is not mediated by specific properties of each channel, but rather by the balance between background inward and outward currents. Therefore, the passive properties of the fast variable were responsible for the flexible control of the frequency and setpoint of the slow variable in the SS firing pattern.

Figure 3

Balance of background inward and outward currents controls oscillation density and calcium setpoint of the SS firing pattern

(A) Relationships between channel conductance (leak K+ channel (left) and leak Na+ channel (right)) and oscillation density of the SS firing pattern (top). The corresponding representative membrane potential and [Ca2+] at each are shown in the bottom panels. Parameter sets are the representatives chosen from major subgroups of parameter sets shown in (B).

(B) The “decreasing” and “increasing” parameter sets are determined by changes in oscillation density when each parameter is modified 2.5. Parameter sets that do not change their oscillation density more than 2.5 are not shown in this figure.

(C) Proportion of the decreasing and the increasing parameter sets among 1,166 parameter sets obtained by the random parameter search.

(D) Changes in oscillation frequencies when two sets of channel conductance are modulate at the same time. and modulate oscillation density competitively (left), while and , and and modulated cooperatively (middle and right, respectively).

(E and F) Relative changes in current through each channel in each phase in response to the shifts in (E) and (F) values in each representative parameter set of the major subgroup. Changes in are larger than another current in both cases.

(G and H) Trajectories and bifurcation diagrams in response to the parameter shifts of representative parameter sets in major subgroups of (G) and (H), corresponding to (A). In major subgroups, when outward current increases (G, left) or inward current decreases (H, right), the bifurcation structures change, followed by decreases in calcium setpoint and oscillation density

To investigate the underlying mechanism behind these controls, we examined the changes in the current proportions in response to changes in the background inward/outward current balance. It was shown that the proportion of was significantly controlled by the changes in background currents in a manner that compensates them (Figures 3E, 3F, S4A, S4B, S4E, S4F, S4I, and S4J). This implies that the changes in calcium setpoint are forced to compensate for the changes in background currents by controlling . These results support the uniqueness of the flexible control of oscillation density and calcium setpoint in the SS firing pattern because of its relatively small dependence on in the outward current (Figure 2A). We subsequently examined the changes in the phase space and bifurcation diagram in response to changes in the background inward/outward current balance. In the major subgroups, the decrease in outward background current (or the increase in background inward current) resulted in a higher calcium setpoint, faster delayed stability loss, and shorter interval between the two bifurcations, whereas the increase in outward background current (or the decrease in background inward current) resulted in the opposite change (Figures 3G, 3H, S4G, and S4H). On the other hand, in the minor subgroups, the increased calcium setpoint resulted in faster delayed stability loss and a much longer interval between the two bifurcations, meaning a longer duration of each phase in total (Figures S4C, S4D, S4K, and S4L). These results indicate that the changes in background currents dynamically control the mathematical system of the SS firing pattern, resulting in a flexible control of oscillation density and calcium setpoint. Consequently, in the SS firing pattern, both the small dependence on and flexible mathematical structure dominated by the stability of a single fixed point enabled the dynamical control in response to the changes in background currents.

Bidirectional control of oscillation density and the calcium setpoint by balance between background inward and outward currents based on the AN model

To assess whether the molecular configuration and mathematical structure of the SS firing pattern in the RAN model met the necessary conditions in the AN model, we analyzed the current during the SS firing pattern and the responses against the changes in the conductance in the AN model. The analysis of current for the representative parameter set and the whole parameter sets revealed that (1) and contributed to the spikes predominantly during bursting phases; (2) outward current during silent phases was mainly dominated by and , although also participated a little; (3) changes in [Ca2+] had little effect on the composition of current during each phase because of the small occupancy of (Figures 4A, 4B, and S5A–S5C). These results indicate that the intracellular mechanism required for producing intermittent firings without up or down states was preserved in the AN model. Regarding the control of oscillation density and calcium setpoint of the SS firing pattern by channel conductance, , , , and regulated oscillation density and the calcium setpoint in the same way for both major and minor subgroups as in the RAN model (Figures 4C–4E and S6A–S6F), although the minor subgroups for and were relatively larger in the AN model than in the RAN model (Figures 4E and S6D). Furthermore, channels through which ions flowed in the same direction controlled oscillation density cooperatively, while channels through which ions flowed in the opposite direction controlled competitively (Figures 4F, S6G, and S6H). Therefore, these results indicate that the characteristics of the SS firing pattern derived from its mathematical structure in the RAN model are generally preserved. Consequently, even under perturbation from various sources of current, the molecular and mathematical configurations and properties of the SS firing pattern, which were found in the simplified model, were well preserved in the AN model.

Figure 4

The balance between inward and outward currents through leak channels controls oscillation density of the SS firing pattern based on the AN model

(A) Proportion of (left), (middle), and (right) of 1,112 parameter sets obtained by the random parameter searches in the AN model (top) and their averages with standard deviations (bottom).

(B) Distribution of proportions of outward current in bursting (top) and silent (bottom) phases of the SS and SWS firing pattern in the AN model.

(C) Relations between channel conductance (leak K+ channel (left) and leak Na+ channel (right)) and oscillation density of the SS firing pattern based on the AN model (top). The corresponding representative membrane potential and [Ca2+] at each are shown in the bottom panels. Parameter sets are the representatives chosen from subgroups of parameter sets shown in (D).

(D) The “decreasing” and “increasing” parameter sets are determined by changes in oscillation density when each parameter is modified 2.5.

(E) Proportion of the decreasing and the increasing parameter sets among 1,166 parameter sets obtained by the random parameter search.

(F) Changes in oscillation density when two sets of channel conductance are modified at the same time.

(G) Schematic illustration of the mechanism in which the balance between inward and outward currents through leak channels control oscillation density of the SS firing pattern

Discussion

Molecular configuration of spindle oscillations

In this study, we first demonstrated that the AN model validly recapitulated the SS firing pattern. We subsequently developed the RAN model, which is a minimal model for the SS firing pattern based on the AN model, composed of six components and three variables. Notably, the composition of the ion channels in the RAN model was consistent with previous computational and experimental studies (Astori et al., 2011; Destexhe et al., 1993; Espinosa et al., 2008; Kim and McCormick 1998; Wimmer et al., 2012). The RAN model proposed the basic molecular configuration for generating the characteristic waveform of the SS firing pattern (i.e., intermittent firings without up or down states, and subthreshold waxing and waning oscillations). First, silent phases without down states resulted from a lower dependence on compared to and . This composition of outward currents made silent phases stable and less susceptible to changes in [Ca2+]. On the other hand, bursting phases without up states attributed to spikes are predominantly generated by voltage-dependent currents and , resulting in large AHPs. From a mathematical perspective, spikes with large AHPs correspond to the limit cycle around the unstable equilibrium. Second, subthreshold waxing and waning oscillations were caused by delayed stability loss in the subcritical Andronov-Hopf bifurcation and slow convergence to the stable focus in fold-cycle bifurcation. Furthermore, this molecular mechanism revealed that the oscillation density and the setpoint of intracellular calcium concentration (i.e., the slow variable) were controlled bidirectionally by the balance between background inward and outward currents, or the passive properties of the fast variable V. The insights obtained from the RAN model were ascertained in the original model, which contained various types of ion channels and ionotropic receptors. Therefore, we developed a minimal model of groups of neurons recapitulating spindle oscillations and proposed a theoretical design principle of spindle oscillations and its properties to control calcium setpoint and oscillation density.

Based on these findings, we suggest a simple, molecular circuit that controls the oscillation density and calcium setpoint of the SS firing pattern by balancing background currents (Figure 4G). Here, we focus on the major subgroups of the parameter sets, considering their dominance and high calcium concentrations in the minor subgroups. When background outward current increases or inward current decreases, the outward decreases to compensate for the change in background currents by lowering the calcium setpoint. The decrease in the calcium setpoint results in a weakening of the reactivity of calcium-dependent K+ channels against [Ca2+] based on the formulation. Hence, the whole system becomes robust to changes in [Ca2+], lengthening both bursting and silent phases and decreasing the oscillation density. When background outward current decreases or inward current increases, the opposite phenomenon occurs. Consequently, we propose a theoretical design principle of spindle oscillations, based on which we discover a molecular circuit controlling the oscillation density and calcium setpoint of sleep spindles.

Possible role of the control of spindle density and calcium setpoint

The oscillation density of the SS firing pattern in our model corresponds to the density of sleep spindles under the physiological conditions (Destexhe et al., 1993; Leresche et al., 1991). Therefore, the molecular mechanism proposed in this study could underlie the flexible control of spindle density after memory tasks and the adaptable nesting by slow oscillations ( 1 Hz) and delta waves (0.5–4 Hz), both of which are deeply involved in memory consolidation during sleep (Chatburn et al., 2013; Gais et al., 2002; Kim et al., 2019; Latchoumane et al., 2017; Maingret et al., 2016; Reynolds et al., 2018; Silversmith et al., 2020; Sirota et al., 2003). Furthermore, in pathophysiological contexts, patients with schizophrenia or early-stage Alzheimer disease exhibit a lower density of sleep spindles (Manoach et al., 2016), and it has been suggested that sleep spindles and their density can be diagnostic indicators and therapeutic markers and targets. Thus, the proposed molecular mechanism controlling spindle density may have far-reaching biological and therapeutic implications. On the other hand, it was previously reported that optogenetic activation and inhibition of TRN inhibitory neurons induced an increase and decrease in spindle density, respectively (Halassa et al., 2011; Thankachan et al., 2019). In the next step, therefore, it would be important to combine the insights for the network connectivity, particularly the thalamocortical circuit, and the neuron-intrinsic properties provided in this study for further understanding of sleep spindles.

In addition, the model revealed that intracellular calcium concentration changed with the density of spindle oscillations. Previous studies have demonstrated that intracellular calcium plays an essential role in generating sleep spindles and maintaining their periodicity by controlling hyperpolarization-activated cation current (Bal and McCormick 1996; Luthi and McCormick, 1998a, 1988b). In the present study, our models endorsed and expanded these findings; we connected the rhythmogenesis of sleep spindles by calcium-related components and the flexible control of spindle density (Destexhe et al., 1993; Luthi and McCormick, 1998b). Notably, the models demonstrated that density and calcium setpoint could be controlled by the balance between inward and outward currents at the resting membrane potential, which was not specific to the specific current such as . We believe that the proposed mechanism could be helpful in understanding the molecular mechanisms of the relationship between spindle density and cognitive functions such as memory consolidation.

Mathematical structure of spindle oscillations

A bifurcation in the dynamical system occurs when a continuous change in a parameter value causes a qualitative or topological change in the system, as has been implied in many biological systems (e.g., cell cycles, developmental processes, and predator-prey systems) (Borisuk and Tyson 1998; Fussmann et al., 2000; Marco et al., 2014). These previously proposed design principles of biological systems have helped to understand the complex behavior of the system (Jolley et al., 2012; Sugai et al., 2017; Yamaguchi et al., 2021). Intermittent firings alternating between near-steady state and trains of rapid spikes are also described by the bifurcation theory. The type of bifurcations that drive the transition between phases has been used to classify the firing patterns into different categories (Izhikevich 2000, 2007; Rinzel 1987), but their connections to molecular properties have been largely elusive. In this study, using the simplified model, we revealed that the spindle oscillation was classified into the sub-Hopf/fold cycle oscillation or elliptic bursting, whereas the slow oscillation was classified as the fold/homoclinic oscillation or square-wave bursting. The sub-Hopf/fold cycle oscillation, along with the fold/homoclinic oscillation, is one of the mathematical structures for generating intermittent firings (Rinzel 1987). They have also been characterized in several electrophysiological experiments (Amir et al., 2002; Del Negro et al., 1998). Notably, one of the unique characteristics of the sub-Hopf/fold cycle oscillation, such as the SS firing pattern, is the flexibility of the system. Both the subcritical Hopf bifurcation and the fold limit cycle bifurcation are controlled by the stability of the fixed point, while both the fold bifurcation and the saddle homoclinic orbit bifurcation are controlled by the position of a set of fixed points (Figures 3G and 3H). Thus, we characterized a theoretical design principle of spindle oscillations, which explained their unique features by the underlying mathematical mechanism. The proposed design principle of sleep spindles is not dependent on specific circuits or properties of neurons; thus, it could be positioned as a simplified and general view of the potential mechanisms of sleep spindles. Furthermore, considering its validity and reasonability as a mathematical structure, the basic molecular configuration of the SS firing pattern and its flexibility characterized in this study might underlie neural oscillation other than sleep spindles, especially in regions where the flexible changes in calcium setpoint or oscillation density are required. Indeed, the stomatogastric ganglion of the crab Cancer borealis dynamically changes its oscillation density depending on resting potential, implementing a fast pyloric rhythm and a slower gastric mill rhythm (Weimann and Marder 1994). Therefore, the proposed design principle of spindle oscillations explained the generation of a variety of rhythmic firings and calcium dynamics from a mathematical perspective, which could be applied to other oscillatory systems.

The AN model as a tool for investigating the basic configuration of collective neuronal activity

Here, we leveraged the AN model and RAN model to determine a design principle for generating sleep spindles. This result supports the idea that simple computational models of groups of neurons could be a promising approach for extracting the basic molecular configuration of neural oscillations. However, the relationships between oscillations with different rhythms, such as slow and spindle oscillations, remain to be investigated in our models. In the present study, we demonstrated that there were high similarities between the minimum models for these two oscillations and the only difference was in the property of the delayed rectifier potassium channel. We consider that it would be important as a next step to investigate the interplay and similarity of these two oscillations using these simplified models, which could lead to a better understanding of the molecular mechanisms of the oscillations controlling cognitive processes (Kim et al., 2019; Latchoumane et al., 2017; Maingret et al., 2016). Furthermore, as in previous studies (Rasmussen et al., 2017; Tatsuki et al., 2016; Yoshida et al., 2018), the proposed mechanism behind sleep spindles should be validated in subsequent studies. In particular, we need further validation from experimental perspectives. Given that each type of channel in the AN model has many subtypes and variants that can affect their conductance, we need a comprehensive genetic screening system for collective firing patterns occurring in groups of neurons. However, such a comprehensive genetic screening for spindle phenotypes in animals is not practical because it is costly. One promising approach is the combination of primary neuronal cultures and a microelectrode array (MEA) system. The key features of the MEA system with neuronal cultures are that they can record and stimulate neurons at multiple sites simultaneously, reduce the effects of specific neuronal connectivity, and achieve higher throughput compared to other electrophysiological experiments (Dunlop et al., 2008). Furthermore, the SWS firing pattern was recently recorded in that system (Saberi-Moghadam et al., 2018). Therefore, it can be a suitable platform for comprehensive genetic screening to validate the hypothesis obtained from the AN model. For finer control over the conductance of a specific channel than genetic manipulations, a dynamic clamp would be suitable. In the dynamic clamp, we give a neuron real-time feedback and thus can control the conductance of the specific channels as we like (Prinz et al., 2004; Sharp et al., 1993). Therefore, by controlling the conductance of leak channels, we would be able to test the mechanism behind spindle oscillations suggested in this study. Taken together, simple models of groups of neurons such as the AN models, combined with validation by these experiments, could be a systematic approach for elucidating collective cellular activities.

Limitations of the study

In this study, we developed a minimal model of spindle oscillations by simplifying the AN model which contained a variety of ion channels, receptors, and exchangers/pumps expressed in many types of neurons. However, there are some components that are not included in our models but have been suggested by previous studies to be important for generating sleep spindles. Hyperpolarization-activated cyclic nucleotide-gated (HCN) channels are one of such components. The current through HCN channels () is known as pacemaker current which plays a role in controlling cardiac and neuronal rhythmicity (DiFrancesco and Tortora 1991; Robinson and Siegelbaum 2003). As discussed above, they have also been implicated to be important for generating and controlling sleep spindles (Fernandez and Lüthi, 2020; Luthi and McCormick, 1998a, 1988b). Although we incorporated limited types of molecular components in our models to focus on a design principle of spindle oscillations, the involvement of other components, such as HCN channels, should be considered carefully.

In addition, although our Hodgkin-Huxley-type models made it possible to connect molecular properties and electrical activity of neurons, the complexity of the models should be taken into account. Compared to more computationally efficient models proposed in previous studies (Izhikevich 2000, 2007; Rinzel 1987), there are various solutions to generate and control spindle oscillations in our models due to their complexity, for example, as demonstrated in Figures S3G and S6G. These different types of solutions sometimes make it difficult to interpret the results and could possibly lead to inconsistent conclusions. We performed simulations and analysis for thousands of parameter sets to ensure the generalizability of our model, but it would also be important to develop one of the simplest models and compare it with our models.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by Prof. Hiroki R. Ueda with e-mail: uedah-tky@umin.ac.jp.

Materials availability

This study did not generate new materials.

Method details

Original model (AN model)

To analyze the spindle firing pattern, a computational model was constructed based on the previous study (Tatsuki et al., 2016). The differential equations for this model are described below., where C is the membrane capacitance, A is the area of a single neuron, V is the membrane potential, [Ca2+] is the intracellular calcium concentration, and (X: any types of ion channels or ionotropic receptors in the model) denotes the current through X. Each function is listed below.

Please note that intrinsic (non-synaptic) ionic currents (e.g. ) should be multiplied by 10 to adjust its unit to nanoampere (nA) if the numerical values listed in Table S1 are directly used in the numerical simulation. The constant values are listed in Tables S1–S4 and the initial values of these differential equations are listed in Table S5. All differential equations were solved in . Sampling frequency was set to 1,000 Hz for parameter searches and 10,000 Hz for detailed analyses.

In some analyses, the conductance of leak channels was divided into the conductance of leak K + channels () and the conductance of leak Na+ channels () as follows:

These definitions satisfy the following conditions:

Simplified-AN (SAN) and Reduced-AN (RAN) models

The SAN model, a minimal model for the SWS firing pattern, is described on the condition that , and are set to zero in the original model according to the previous study (Yoshida et al., 2018). On the other hand, the RAN model, a minimal model for the SS firing pattern, is described on the condition that , and are set to zero in the original model.

Parameter search

The method of the parameter searches was based on the previous studies (Tatsuki et al., 2016; Yoshida et al., 2018). A parameter set includes 13 values (, , , , , , , , , , , , and ). All parameter sets were generated so that the logarithm of each parameter has a uniform distribution. Let be a distribution whose probability density function is as follows.

The values of , , , , , , , , and (mS/cm2) were generated according to , , , and (μS) were generated according to , and (msec) was generated according to .

For each parameter set, the differential equations were solved by using integrate.odeint method in SciPy 1.1.0. Let be a membrane potential obtained by solving the differential equation, and was defined as . The lists of were first analyzed by the discrete Fourier transform after normalization and detrending. The frequency whose Fourier power was the largest was defined as the peak frequency. The sets and was defined as follows.

A set of time points when spikes occurred, S, was defined as follows.

Empirically, P was referred to for the SWS firing pattern while B for the SS firing pattern. The time points of spike events were divided into subgroups as follows.

When , , defined as follows: was denoted as a bursting phase, while the other was denoted as silent phases. Based on these spike metrices, the firing pattern was classified as follows. First, the parameter sets that at least one of the value was bigger than 200 mV were excluded. Then the solutions were classified into four categories (Resting, SWS/SS, Awake, and SWS/SS with few spikes) based on the peak frequency and the number of spikes. The category SWS/SS was further divided into SWS and SS depending on the minimum membrane potentials during bursting and silent phases. If the minimum membrane potential during bursting phases were smaller than that during silent phases, the firing pattern was classified as SS, and if not, the firing pattern was classified as SWS. All the SWS and SS firing patterns were confirmed manually. The results of the parameter searches are listed in Table S6.

PCA and factor analysis

To conduct dimension reduction, each parameter value g was converted into logarithmic scale and normalized based on . The normalized was expressed as follows.

The PCA was conducted for . Let be the contribution ratio and be the coefficient of component of PCi (). The principal component score and the factor loading were calculated as follows.

Normalization of the oscillation density

The membrane potential, [Ca2+], and were normalized based on the rhythms of altering bursting and silent phases. First, differential equations were solved in . When , n was multiplied by two, and the differential equations were solved again. This procedure was recursively repeated until . The parameter sets which could not reach were excluded from the following analyses. Let be the ith element of , and the normalization was executed for each defined as follows.

Knock-out (KO) analysis

In the KO analysis for channels or ionotropic receptors, the conductance of the target channel was divided by 1,000. For Ca2+ pumps, the time constant was multiplied by 1,000. After KO of the target component, the differential equations were solved, and subsequently the firing pattern was classified based on the criteria above mentioned. The percentage of the parameter set in which the SS firing pattern was not abolished after KO was calculated for each channels, receptors, and exchangers/pumps.

Selection of a representative parameter set

PCA was conducted to the parameter sets as mentioned above. Subsequently, the probability density function was calculated for PC1 and PC2 based on KDE. The bandwidth of KDE was estimated based on Scott's method. The representative parameter set was finally chosen manually from the three parameter sets whose probability density was the highest.

Analysis of current

The proportion of currents were calculated as in the previous study. The proportion of channel X's inward current and outward currents was defined as follows:

Let be the sets of time points of silent phases (i.e., ), and be the ith element of the . The average current during the bursting phases and during the silent phases were calculated as follows.

The average proportion of channel X's outward current and inward current during the silent phases, and the proportion of channel X's outward current and inward current during bursting phases were defined as follows.

Analysis of the oscillation density in response to parameter shifts

When shifting a single parameter value for each parameter set, the parameter value g was shifted to as follows.

When shifting two parameter values at the same time, the parameter values and were shifted to and as follows.

After modifying the parameter value(s), the differential equations were solved. Subsequently, the density of the SS firing pattern was defined as follows.

Phase space

[Ca2+] was fixed based on the following equation.

The nullclines of V and were described as follows.

[Ca2+]fixed values in figures are described on the top of each figure.

Linear stability analysis

The function and were defined as follows.

When [Ca2+] was fixed at [Ca2+]fixed, the fixed point was calculated by solving . These nonlinear equations were solved optimize.fsolve method in SciPy 1.1.0. Subsequently, the Jacobian matrix J was described as follows.

The stability of the fixed point was categorized as follows: (1) saddle (); (2) stable focus (, , and ); (3) stable node (, , and ); (4) unstable focus (, , and ); (5) unstable node (, , and ).

Bifurcation diagram

To draw a bifurcation diagram, [Ca2+] [] was determined arbitrarily based on the solution of differential equations. Subsequently, the equations expressed as follows were solved.

The stability of the solutions were determined based on the linear stability analysis. The solutions were plotted on V-[Ca2+] plane and colored based on their stability.

Quantification and statistical analysis

Figures represent averaged or representative results of multiple independent simulations. The figure legends provide details concerning the simulations.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Software and algorithms
Python 3.7.3	Python Software Foundation	https://www.python.org/
Scripts for simulations and analyses	This paper	https://github.com/DSPsleeporg/an_spindle