iSciencell

OPEN ACCESS

ArticleA design principle of spindle oscillations inmammalian sleep
Minimal model of spindle oscillations
Na+, K+

Leak channel

gL+

Persistent Na channel

gNaP

Ca2+

Voltage-gated Ca2+ channel

gCa

K+

Slowly inactivating K+ channel gKS

K+

Ca2+-dependent K+ channel

gKCa

Ca2+

Ca2+-pump/exchanger

TCa

Na

+

Tetsuya Yamada,Shoi Shi, Hiroki R.
Uedauedah-tky@umin.ac.jp

Highlights
A minimal, HodgkinHuxley-type model ofspindle oscillations isdeveloped

Mathematical analysis of spindle oscillations
Bifurcation diagram

Phase space

The property of delayed
rectifier K+ channels
characterizes spindle
oscillations

V nullcline
mKS nullcline

Limit cycle (Vmin)mKS
V

Limit cycle (Vmax)[Ca2+]

Stable focus
Unstable focus
Unstable node
SaddleBalance of leak currents controls spindle oscillations
Outward current

The combination ofbifurcations specifies
spindle oscillations/ Inward current

Spindle oscillations arecontrolled by the balanceof inward and outward
currents

IKCa, Calcium concentration

ExtracellularNa+

Sensitivity of KCa channel
against Ca2+ increase
Oscillation density

Intracellular

K+

Ca2+

Leak channel
Ca

2+

pump

Ca2+-dependent K+ channel
Voltage-dependent Ca2+ channel

Yamada et al., iScience 25,103873
March 18, 2022 ª 2022 The
Authors.
https://doi.org/10.1016/
j.isci.2022.103873

iSciencellOPEN ACCESS

ArticleA design principle of spindle
oscillations in mammalian sleep
Tetsuya Yamada,1,2,4 Shoi Shi,2,3,4 and Hiroki R. Ueda2,3,5,*
SUMMARY

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, cellintrinsic mechanisms, especially a design principle, of neural oscillations remain
largely elusive. Herein, we developed a minimal, Hodgkin-Huxley-type modelof 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 intracellularmechanisms 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 ingenerating 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 Lü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.

1Faculty of Medicine, The

University of Tokyo, 113-0033
Tokyo, Japan
2Department of Systems

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 computationalmodels 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

Pharmacology, Graduate
School of Medicine, The
University of Tokyo, 113-0033
Tokyo, Japan
3Laboratory for Synthetic

Biology, RIKEN Center for
Biosystems Dynamics
Research, 565-0871 Osaka,Japan
4These authors contributed

equally
5Lead contact*Correspondence:uedah-tky@umin.ac.jp
https://doi.org/10.1016/j.isci.
2022.103873

iScience 25, 103873, March 18, 2022 ª 2022 The Authors.
This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1ll

OPEN ACCESS

constellation of ion channels that were commonly expressed in the cortex and transmitted its output toequivalent 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 patterndemonstrate that a simple but comprehensive neuronal model can reveal molecular mechanisms anddesign 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 thethalamic 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). Thesestudies underscore the importance of investigating molecular mechanisms, in addition to specificneuron-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 Lü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, weinvestigated 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 theconfiguration 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, dynamicalsystem 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 oscillationsof 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 fastvariable (i.e., membrane potential), competitively controlled the base value and time constant of theslow variable (i.e., cytosolic calcium concentration). Overall, we proposed a theoretical design principleof spindle oscillation by leveraging a minimal model of groups of neurons, and discovered that the balanceof 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, hNa , nK , hA , mKS , sAMPA , xNMDA , sNMDA , sGABA , and [Ca2+]) described byordinary 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 downstates (Figure 1B). Close inspection revealed that there were the waxing and waning oscillations before andafter the bursting phases, respectively, and the rhythmic activity of the membrane potential (i.e., alternating

2

iScience 25, 103873, March 18, 2022

iScience
ArticleiScience
Articlell

OPEN ACCESS

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.

iScience 25, 103873, March 18, 2022

3

ll

OPEN ACCESS

iScience
ArticleFigure 1. Continued(C) Detailed membrane potential (top) and [Ca2+] (bottom) dynamics in the SS firing pattern. The waxing and waning subthreshold oscillations andafterhyperpolarizations (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 beindispensable 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

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 patternconsists 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.
To this end, we collected 1; 139 parameter sets that recapitulated the firing pattern from more than50; 000; 000 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 3 103 %. 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 3 103 %).
Moreover, all the collected SS firing patterns were homogeneous in terms of the dynamics of themembrane 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 firstcompared 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 (gL ), voltage-gated Na+ channels(gNa ), Hodgkin-Huxley-type K+ channels (gK ), slowly inactivating K + channels (gKS ), persistent Na + channels(gNaP ), and the time constant of Ca2 + pumps (tCa ) 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 gL , gNa , gK , gNaP , tCa(yellow bars), and especially gKS (magenta bar) contributed to PC 1 (Figure 1G). Moreover, gNaP , tCa , andgKS of the SS firing pattern tended to be larger than those of the SWS firing pattern, while Na+ and gK
tended to be smaller (Figure 1G). Therefore, these components, especially gKS , make the SS parameter
sets distinct from those of the SWS. We subsequently compared the oscillation parameters (average lengthof 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 differencein the fast subsystem but also comes from the fundamental mechanisms underlying them. Based on theseanalyses of the AN model, we discovered that the AN model validly recapitulated the SS firing patternwhich 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 moreelaborate analyses.

4

iScience 25, 103873, March 18, 2022

iSciencell

ArticleOPEN ACCESS

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, nK , and [Ca2+]), and revealed the important role of leak K+ channels in the SWS firing pattern through bifurcation analysis anddynamical 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 firstknocked out each component in each parameter set. We subsequently calculated the percentage ofparameter sets in which the SS firing pattern turned into any other firing pattern after knocking out oneof 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 theSS 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 10; 000; 000
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 intermittencyof the SS firing pattern without slow [Ca2+] rhythm. Therefore, we added either NMDA receptors (gNMDA ) or
voltage-gated Ca2+ channels (gCa ) 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 3 103 % and 5.0 3
103 %, respectively (Figure 1J). Although Model 2 was slightly more efficient in producing the SS firing
pattern, it contained five variables (V, mKS , xNMDA , sNMDA , and [Ca2+]), whereas Model 3 contained three variables (V, mKS , 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 ionchannels 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 beimportant 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 configurationof 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 (gCa and gKCa ), leak channels (gL ), Na+ channels (gNaP ), and K+ channels (gKS ). 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 patternfrom 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. gL and gNaP of the SS firing pattern in the RAN model were generally larger than thoseof 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 basedon 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

iScience 25, 103873, March 18, 2022

5

ll

OPEN ACCESS

iScience
ArticleFigure 2. Small changes in IKCa 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 (IKL ) (left), slowly inactivating K+ channel (IKS ) (middle), and calcium-dependent K+ channel (IKCa ) (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 mKS of the representative parameter set of the SS firing pattern in the phase space. The coiled partof the trajectory represents bursting phases, while the linear part represents silent phases.

6

iScience 25, 103873, March 18, 2022

iSciencell

ArticleOPEN ACCESS

Figure 2. Continued(E) Phase planes and stream plot of the representative parameter set on the V-mKS plane when [Ca2+] is low (67 mM) and high (90 mM). 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 anda bursting phase to a silent phase, respectively.
(F) Fixed points (that is, intersections of V and mKS 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

depolarization were generated by fast, transient INaP and subsequent delayed-rectifying IKS , respectively,
and that (2) outward currents during silent phases were dominated by IKS and IKL which resulted in silent
phases without profound hyperpolarization (Figure 2A). These observations were supported by otherparameter 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, IKCa , 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 througheach channel. However, it remains unclear how small changes in IKCa drive the dynamic transitions between
phases, and how subthreshold waxing and waning oscillations appear.

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 thatthe 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 mKS in thephase 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 mM), the trajectory diverged from the unstable equilibrium
and converged to the stable limit cycle, whereas at higher [Ca2+] ([Ca2+] = 90 mM), the trajectory converged
to the stable equilibrium point (Figure 2E). To illustrate how the stability of the fixed points controls theentire 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 subsequentsilent 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 gNaP and gKS , which generated bursting phases without up states.
Furthermore, small changes in [Ca2+] and IKCa 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 IKCa 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;

iScience 25, 103873, March 18, 2022

7

ll

OPEN ACCESS

iScience
ArticleFigure 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 gKL are shown in the bottom panels. Parameter sets are the representatives chosenfrom 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.

8

iScience 25, 103873, March 18, 2022

iSciencell

ArticleOPEN ACCESS

Figure 3. Continued(D) Changes in oscillation frequencies when two sets of channel conductance are modulate at the same time. gKL and gNaL modulate oscillation density
competitively (left), while gKL and gKS , and gNaL and gNaP 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 gKL (E) and gNaL (F) values in each representative
parameter set of the major subgroup. Changes in IKCa 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 gKL (G) and gNaL
(H), corresponding to (A). In major subgroups, when outward current increases (G, left) or inward current decreases (H, right), the bifurcation structureschange, followed by decreases in calcium setpoint and oscillation density

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 molecularcomponents 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 gKL and gKS were negatively correlated with oscillation density (denoted as ‘‘decreasing’’), while gNaL and gNaP 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 gKL and gKS (or the smaller gNaL and gNaP ) were, thelower [Ca2+] was, in the major parameter sets (Figures 3A and S3C). In the minor parameter sets, on theother 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 thebase value of the fast variable V, controlled the frequency and setpoint of the calcium oscillation describedby the slow variable [Ca2+] in two ways: the major one in which oscillation density and calcium setpoint werepositively 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 offour parameters (gNaL , gNaP , gKL , and gKS ). 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., gKL and gKS )cooperatively controlled oscillation density and the calcium setpoint, while the conductance of channels
carrying current in the opposite directions (e.g., gKS and gNaP ) competitively controlled them (Figures 3Dand 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.
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 thatthe proportion of IKCa was significantly controlled by the changes in background currents in a manner thatcompensates 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 IKCa . 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 IKCa 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, fasterdelayed 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 thechanges 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 IKCa and flexible mathematical structure dominated by the stabilityof a single fixed point enabled the dynamical control in response to the changes in background currents.

iScience 25, 103873, March 18, 2022

9ll

OPEN ACCESS

10

iScience 25, 103873, March 18, 2022

iScience
ArticleiSciencell

ArticleOPEN ACCESS

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 IKL (left), IKS (middle), and IKCa (right) of 1,112 parameter sets obtained by the random parameter searches in the AN model (top) and theiraverages 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 gKL are shown in the bottom panels. Parameter sets are therepresentatives 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

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) INaP and IKS contributed
to the spikes predominantly during bursting phases; (2) outward current during silent phases was mainly
dominated by IKS and IKL , although IAR 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 IKCa (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, gNaL , gKL , gNaP , and gKS 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 gNaP and gKS were relatively
larger in the AN model than in the RAN model (Figures 4E and S6D). Furthermore, channels through whichions flowed in the same direction controlled oscillation density cooperatively, while channels through whichions flowed in the opposite direction controlled competitively (Figures 4F, S6G, and S6H). Therefore, theseresults 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.

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 ANmodel, 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 IKCa compared to
IKL and IKS . 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 INaP and IKS , resulting in large AHPs. From a mathematicalperspective, 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, thismolecular 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

iScience 25, 103873, March 18, 2022

11

ll

OPEN ACCESS

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 IKCa 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 amolecular 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 molecularmechanism 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 biologicaland therapeutic implications. On the other hand, it was previously reported that optogenetic activationand 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 Ih (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 calciumrelated 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 thebalance between inward and outward currents at the resting membrane potential, which was not specificto the specific current such as Ih . 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 toclassify 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

12

iScience 25, 103873, March 18, 2022

iScience
ArticleiSciencell

ArticleOPEN ACCESS

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 uniquecharacteristics of the sub-Hopf/fold cycle oscillation, such as the SS firing pattern, is the flexibility of thesystem. 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 arecontrolled by the position of a set of fixed points (Figures 3G and 3H). Thus, we characterized a theoreticaldesign 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 orproperties 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 afast 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, therelationships 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 ofthe 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 validationfrom experimental perspectives. Given that each type of channel in the AN model has many subtypes andvariants 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 forspindle 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 comparedto other electrophysiological experiments (Dunlop et al., 2008). Furthermore, the SWS firing pattern wasrecently 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 byprevious 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 (Ih ) 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

iScience 25, 103873, March 18, 2022

13

ll

OPEN ACCESS

important for generating and controlling sleep spindles (Fernandez and Lüthi, 2020; Luthi and McCormick,
1998a, 1988b). Although we incorporated limited types of molecular components in our models to focus ona 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 tomore 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 theircomplexity, 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
Detailed methods are provided in the online version of this paper and include the following:d KEY RESOURCES TABLE
d RESOURCE AVAILABILITY

B Lead contact
B Materials availability
B Data and code availabilityd METHOD DETAILSB Original model (AN model)
B Simplified-AN (SAN) and Reduced-AN (RAN) models
B Parameter search
B PCA and factor analysis
B Normalization of the oscillation density
B Knock-out (KO) analysis
B Selection of a representative parameter set
B Analysis of current
B Analysis of the oscillation density in response to parameter shifts
B Phase space
B Linear stability analysis
B Bifurcation diagram
d QUANTIFICATION AND STATISTICAL ANALYSIS

SUPPLEMENTAL INFORMATION
Supplemental information can be found online at https://doi.org/10.1016/j.isci.2022.103873.

ACKNOWLEDGMENTS
The authors would like to thank all the laboratory members at RIKEN Center for Biosystems Dynamics
Research and the University of Tokyo, particularly F. Tatsuki and T. Katsumata for discussions on mathematical modeling. This work was supported by grants from the Japan Agency for Medical Research and Development – Core Research for Evolutional Science and Technology (AMED-CREST) (JP17gm0610006,
H.R.U.), Brain Mapping by Integrated Neurotechnologies for Disease Studies (Brain/MINDS)(JP17DM0207049, H.R.U.), Basic Science and Platform Technology Program for Innovative Biological Medicine (JP17AM0301025, H.R.U.), KAKENHI Grant-in-Aid from JSPS (Scientific Research S 18H05270, H.R.U.),
Young Scientist (19K16487 and 21K15136, S.S.), Transformative Research Areas A (20H05894 and
20H05903, S.S.), Human Frontier Science Program (HFSP) Research Grant Program (HFSP RGP0019/
2018, H.R.U.), Exploratory Research for Advanced Technology (ERATO) grant (JPMJER1904, H.R.U.) fromthe Japan Science and Technology Agency (JST), Moonshot R and D – MILLENNIA Program Grant(JPMJMS2023–25, S.S.) from JST, and an intramural Grant-in-Aid from the RIKEN Center for Biosystems Dynamics Research to H.R.U.

14

iScience 25, 103873, March 18, 2022

iScience
ArticleiSciencell

ArticleOPEN ACCESSAUTHOR CONTRIBUTIONS
T.Y., S.S., and H.R.U. conceptualized the study; T.Y. and S.S. designed the experiments; T.Y. performed thesimulations, analyzed the data, and prepared the figures; T.Y., S.S., and H.R.U. interpreted the results; T.Y.,
S.S., and H.R.U. wrote the manuscript; S.S. and H.R.U. supervised T.Y.

DECLARATION OF INTERESTS
H.R.U. is a member of the iScience’s editorial advisory board.

Received: October 14, 2021
Revised: January 7, 2022
Accepted: February 1, 2022
Published: March 18, 2022
REFERENCES
Alonso, L.M., and Marder, E. (2019). Visualizationof currents in neural models with similar behavior
and different conductance densities. Elife 8,
e42722.

Burke, J., Desroches, M., Barry, A.M., Kaper, T.J.,
and Kramer, M.A. (2012). A showcase of torus
canards in neuronal bursters. J. Math. Neurosci. 2,
1–30.

Amir, R., Michaelis, M., and Devor, M. (2002).
Burst discharge in primary sensory neurons:
triggered by subthreshold oscillations,maintained by depolarizing afterpotentials.
J. Neurosci. 22, 1187–1198.

Buzsaki, G. (2015). Hippocampal sharp waveripple: a cognitive biomarker for episodic
memory and planning. Hippocampus 25, 1073–
1188.

Antony, J.W., Schönauer, M., Staresina, B.P., andCairney, S.A. (2019). Sleep spindles and memory
reprocessing. Trends Neurosci. 42, 1–3.
Bal, T., von Krosigk, M., and McCormick, D.A.
(1995). Role of the ferret perigeniculate nucleus in
the generation of synchronized oscillations
in vitro. J. Physiol. 483, 665–685.
Astori, S., Wimmer, R.D., Prosser, H.M., Corti, C.,
Corsi, M., Liaudet, N., Volterra, A., Franken, P.,
Adelman, J.P., and Luthi, A. (2011). The CaV3.3
calcium channel is the major sleep spindle
pacemaker in thalamus. Proc. Natl. Acad. Sci. U S
A 108, 13823–13828.
Bal, T., and McCormick, D.A. (1996). What stops
synchronized thalamocortical oscillations?
Neuron 17, 297–308.
Bandarabadi, M., Herrera, C.G., Gent, T.C.,
Bassetti, C., Schindler, K., and Adamantidis, A.R.
(2020). A role for spindles in the onset of rapid eye
movement sleep. Nat. Commun. 11, 5247.
Bartos, M., Vida, I., and Jonas, P. (2007). Synaptic
mechanisms of synchronized gamma oscillationsin inhibitory interneuron networks. Nat. Rev.
Neurosci. 8, 45–56.
Bazhenov, M., Timofeev, I., Steriade, M., andSejnowski, T.J. (2002). Model of thalamocorticalslow-wave sleep oscillations and transitions toactivated states. J. Neurosci. 22, 8691–8704.
Borisuk, M.T., and Tyson, J.J. (1998). Bifurcation
analysis of a model of mitotic control in frog eggs.
J. Theor. Biol. 195, 69–85.
Brunel, N., and Wang, X. (2003). What determines
the frequency of fast network oscillations withirregular neural discharges? I. Synaptic dynamics
and excitation-inhibition balance.
J. Neurophysiol. 90, 415–430.

Buzsaki, G., and Wang, X. (2012). Mechanisms of
gamma oscillations. Annu. Rev. Neurosci. 35,
203–225.
Chatburn, A., Coussens, S., Lushington, K.,
Kennedy, D., Baumert, M., and Kohler, M. (2013).
Sleep spindle activity and cognitive performancein healthy children. Sleep 36, 237–243.
Compte, A., Sanchez-Vives, M.V., McCormick,D.A., and Wang, X. (2003). Cellular and network
mechanisms of slow oscillatory activity (< 1 Hz)
and wave propagations in a cortical networkmodel. J. Neurophysiol. 89, 2707–2725.
Del Negro, C.A., Hsiao, C., Chandler, S.H., andGarfinkel, A. (1998). Evidence for a novel burstingmechanism in rodent trigeminal neurons.
Biophys. J. 75, 174–182.
Destexhe, A., Babloyantz, A., and Sejnowski, T.J.
(1993). Ionic mechanisms for intrinsic slow
oscillations in thalamic relay neurons. Biophys. J.
65, 1538–1552.
Destexhe, A., Contreras, D., Sejnowski, T.J., andSteriade, M. (1994). A model of spindle
rhythmicity in the isolated thalamic reticular
nucleus. J. Neurophysiol. 72, 803–818.
Destexhe, A., Contreras, D., and Steriade, M.
(1999). Cortically induced coherence of athalamic- generated oscillation. Neuroscience 92,
427–443.
DiFrancesco, D., and Tortora, P. (1991). Direct
activation of cardiac pacemaker channels byintracellular cyclic AMP. Nature 351, 145–147.

learning or after retrieval in rats. J. Neurosci. 26,
12914–12920.
Espinosa, F., Torres-Vega, M.A., Marks, G.A., andJoho, R.H. (2008). Ablation of Kv3.1 and Kv3.3
potassium channels disrupts thalamocorticaloscillations in vitro and in vivo. J. Neurosci. 28,
5570–5581.
Fernandez, L.M., and Lüthi, A. (2020). Sleep
spindles: mechanisms and functions. Physiol. Rev.
100, 805–868.
Fussmann, G.F., Ellner, S.P., Shertzer, K.W., andHairston, N.G., Jr. (2000). Crossing the Hopf
bifurcation in a live predator-prey system. Science
290, 1358–1360.
Gais, S., Molle, M., Helms, K., and Born, J. (2002).
Learning-dependent increases in sleep spindle
density. J. Neurosci. 22, 6830–6834.
Halassa, M.M., Chen, Z., Wimmer, R.D., Brunetti,
P.M., Zhao, S., Zikopoulos, B., Wang, F., Brown,
E.N., and Wilson, M.A. (2014). State-dependentarchitecture of thalamic reticular subnetworks.
Cell 158, 808–821.
Halassa, M.M., Siegle, J.H., Ritt, J.T., Ting, J.T.,
Feng, G., and Moore, C.I. (2011). Selective opticaldrive of thalamic reticular nucleus generates
thalamic bursts and cortical spindles. Nat.
Neurosci. 14, 1118–1120.
Hutcheon, B., and Yarom, Y. (2000). Resonance,
oscillation and the intrinsic frequency preferencesof neurons. Trends Neurosci. 23, 216–222.
Izhikevich, E.M. (2000). Neural excitability, spiking
and bursting. Int. J. Bifurcat. Chaos 10, 1171–
1266.
Izhikevich, E.M. (2007). Dynamical Systems in
Neuroscience: The Geometry of Excitability andBursting (MIT Press).
Jolley, C.C., Ode, K.L., and Ueda, H.R. (2012). A
design principle for a posttranslational
biochemical oscillator. Cell Rep. 2, 938–950.

Dunlop, J., Bowlby, M., Peri, R., Vasilyev, D., andArias, R. (2008). High-throughput electrophysiology: an emerging paradigm for ion-channel
screening and physiology. Nat. Rev. Drug Discov.
7, 358–368.

Kim, J., Gulati, T., and Ganguly, K. (2019).
Competing roles of slow oscillations and delta
waves in memory consolidation versus forgetting.
Cell 179, 514–526.e13.

Eschenko, O., Molle, M., Born, J., and Sara, S.J.
(2006). Elevated sleep spindle density after

Kim, U., and McCormick, D.A. (1998).
Functional and ionic properties of a slow

iScience 25, 103873, March 18, 2022

15

iSciencell

ArticleOPEN ACCESS

afterhyperpo- larization in ferret perigeniculate
neurons in vitro. J. Neurophysiol. 80, 1222–1235.

reveals epigenetic landscape. Proc. Natl. Acad.
Sci. U S A 111, E5643–E5650.

modulates neural synchrony. J. Neurosci. 40,
4673–4684.

Thankachan, S., Katsuki, F., McKenna, J.T., Yang,
C., Shukla, C., Deisseroth, K., Uygun, D.S.,
Strecker, R.E., Brown, R.E., McNally, J.M., et al.
(2019). Thalamic reticular nucleus parvalbumin
neurons regulate sleep spindles andelectrophysiological aspects of schizophrenia inmice. Sci. Rep. 9, 1–16.

Neske, G.T. (2015). The slow oscillation in corticaland thalamic networks: mechanisms andfunctions. Front. Neural Circuit. 9, 88.

Sirota, A., Csicsvari, J., Buhl, D., and Buzsaki, G.
(2003). Communication between neocortex andhippocampus during sleep in rodents. Proc. Natl.
Acad. Sci. U S A 100, 2065–2069.

von Krosigk, M., Bal, T., and McCormick, D.A.
(1993). Cellular mechanisms of a synchronizedoscillation in the thalamus. Science 261, 361–364.

Rasmussen, R., Jensen, M.H., and Heltberg, M.L.
(2017). Chaotic dynamics mediate brain statetransitions, driven by changes in extracellular ion
concentrations. Cell Syst. 5, 591–603.e4.

Latchoumane, C.V., Ngo, H.-V.V., Born, J., andShin, H. (2017). Thalamic spindles promotememory formation during sleep through triple
phase-locking of cortical, thalamic, and hippocampal rhythms. Neuron 95, 424–435.e6.
Leresche, N., Lightowler, S., Soltesz, I., JassikGerschenfeld, D., and Crunelli, V. (1991). Lowfrequency oscillatory activities intrinsic to rat andcat thalamocortical cells. J. Physiol. 441, 155–174.
Llinas, R.R. (1988). The intrinsic
electrophysiological properties of mammalian
neurons: insights into central nervous system
function. Science 242, 1654–1664.
Luthi, A., and McCormick, D.A. (1998a). H-current:
properties of a neuronal and network pacemaker.
Neuron 21, 9–12.
Luthi, A., and McCormick, D.A. (1998b).
Periodicity of thalamic synchronized oscillations:the role of Ca2+-mediated upregulation of Ih.
Neuron 20, 553–563.
Maingret, N., Girardeau, G., Todorova, R.,
Goutierre, M., and Zugaro, M. (2016).
Hippocampo- cortical coupling mediates
memory consolidation during sleep. Nat.
Neurosci. 19, 959–964.

Prinz, A.A., Abbott, L.F., and Marder, E. (2004).
The dynamic clamp comes of age. Trends
Neurosci. 27, 218–224.

Reynolds, C.M., Short, M.A., and Gradisar, M.
(2018). Sleep spindles and cognitive performanceacross adolescence: a meta-analytic review.
J. Adolescence 66, 55–70.
Rinzel, J. (1987). A formal classification of burstingmechanisms in excitable systems. In Lecture
Notes in Biomathematics (Springer), pp. 267–281.
Robinson, R.B., and Siegelbaum, S.A. (2003).
Hyperpolarization-activated cation currents: frommolecules to physiological function. Annu. Rev.
Physiol. 65, 453–480.
Saberi-Moghadam, S., Simi, A., Setareh, H.,
Mikhail, C., and Tafti, M. (2018). In vitro cortical
network firing is homeostatically regulated: amodel for sleep regulation. Sci. Rep. 8, 1–14.
Schabus, M., Dang-Vu, T.T., Albouy, G., Balteau,E., Boly, M., Carrier, J., Darsaud, A., Deguel- dre,C., Desseilles, M., Gais, S., et al. (2007).
Hemodynamic cerebral correlates of sleep
spindles during human non-rapid eye movement
sleep. Proc. Natl. Acad. Sci. U S A 104, 13164–
13169.
Sharp, A.A., O’Neil, M.B., Abbott, L.F., andMarder, E. (1993). The dynamic clamp: artificial
conductances in biological neurons. Trends
Neurosci. 16, 389–394.

Manoach, D.S., Pan, J.Q., Purcell, S.M., andStickgold, R. (2016). Reduced sleep spindles inschizophrenia: a treatable endophenotype thatlinks risk genes to impaired cognition? Biol.
Psychiat. 80, 599–608.

Shi, S., and Ueda, H.R. (2018). Ca -Dependent
hyperpolarization pathways in sleep homeostasis
and mental disorders. BioEssays 40, 1700105.

Marco, E., Karp, R.L., Guo, G., Robson, P., Hart,
A.H., Trippa, L., and Yuan, G. (2014). Bifurcation
analysis of single-cell gene expression data

Silversmith, D.B., Lemke, S.M., Egert, D., Berke,
J.D., and Ganguly, K. (2020). The degree ofnesting between spindles and slow oscillations

16

iScience 25, 103873, March 18, 2022

2+

Steriade, M., McCormick, D.A., and Sejnowski,
T.J. (1993). Thalamocortical oscillations in thesleeping and aroused brain. Science 262,
679–685.
Sugai, S.S., Ode, K.L., and Ueda, H.R. (2017). A
design principle for an autonomous posttranslational pattern formation. Cell Rep. 19,
863–874.
Tatsuki, F., Sunagawa, G.A., Shi, S., Susaki, E.A.,
Yukinaga, H., Perrin, D., Sumiyama, K., UkaiTadenuma, M., Fujishima, H., Ohno, R., et al.
(2016). Involvement of Ca2+-dependent hyperpolarization in sleep duration in mammals. Neuron
90, 70–85.
Weimann, J.M., and Marder, E. (1994). Switchingneurons are integral members of multiple
oscillatory networks. Curr. Biol. 4, 896–902.
Wimmer, R.D., Astori, S., Bond, C.T., Rovo, Z.,
Chatton, J., Adelman, J.P., Franken, P., and Luthi,
A. (2012). Sustaining sleep spindles throughenhanced SK2-channel activity consolidates
sleep and elevates arousal threshold. J. Neurosci.
32, 13917–13928.
Xu, W., De Carvalho, F., Clarke, A.K., andJackson, A. (2021). Communication from thecerebellum to the neocortex during sleep
spindles. Prog. Neurobiol. 199, 101940.
Yamaguchi, H.Q., Ode, K.L., and Ueda, H.R.
(2021). A design principle for posttranslational
chaotic oscillators. iScience 24, 101946.
Ying, J., and Qin-Sheng, B. (2011). subHopf/foldcycle bursting in the hindmarsh-rose neuronalmodel with periodic stimulation. Chin. Phys. Lett.
28, 090201.
Yoshida, K., Shi, S., Ukai-Tadenuma, M.,
Fujishima, H., Ohno, R., and Ueda, H.R. (2018).
Leak potassium channels regulate sleep duration.
Proc. Natl. Acad. Sci. U S A 115, E9459–E9468.

iSciencell

ArticleOPEN ACCESS

STAR+METHODSKEY RESOURCES TABLEREAGENT or RESOURCESOURCEIDENTIFIER

Python 3.7.3

Python Software Foundation

https://www.python.org/

Scripts for simulations and analyses

This paper

https://github.com/DSPsleeporg/an_spindle

Software and algorithms

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.

Data and code availabilityd No raw data was collected for this study. All results were produced from deterministic simulations andcan be reproduced by using the same parameters.
d All of the code used for simulations and analyses can be found at https://github.com/DSPsleeporg/

an_spindled For any additional questions or information, please contact the lead contact.

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.

CA

  




 


 
dV
=  A IL V + INa V ; hNa + IK V ; Nk + IA V + IKS V ; mKS + ICa V
dt
 

 
 
+ IKCa V ; Ca2 + + INaP V + IAR V  INMDA ðV ; sNMDA ; xNMDA ÞIAMPA ðV ; sAMPA Þ  IGABA ðV ; sGABA Þ
dhNa= 4ðah ðV Þð1  hNa Þ  bh ðV ÞhNa Þ
dtdnK= 4ðan ðV Þð1  nK Þ  bn ðV ÞnK Þdt
dhA hAN ðV Þ  hA
=
dt
t hA
dmKS mKSN ðV Þ  mKS=
dt
t mKS ðV Þ
dsAMPAsNMDA= 3:48f ðV Þ dt
tsNMDA
dsNMDAsNMDA= 0:5xNMDA ð1  sNMDA Þ 
dt
tsNMDA
dxNMDA
xNMDA= 3:48f ðV Þ dt
txNMDA
dsGABAsGABA= f ðV Þ 
dt
tGABA

iScience 25, 103873, March 18, 2022

17

iSciencell

ArticleOPEN ACCESS 2+ 
d Ca2 +Ca=  aCa ðAICa ðV Þ + INMDA ðV ; sNMDA ; xNMDA ÞÞ dttCa
, 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 IX (X: any types of ion channels or ionotropic receptors in themodel) denotes the current through X. Each function is listed below.

IL ðV Þ = gL ðV  VL Þ
INa ðV ; hNa Þ = gNa m3NaN ðV ÞðV  VNa Þ
am ðV ÞmNaN ðV Þ =am ðV Þ + bm ðV Þ
V + 33am ðV Þ =+ 33
1  exp V 10

V + 50bm ðV Þ = 4exp 10V + 50ah ðV Þ = 0:07exp 
101bh ðV Þ =+ 20
1 + exp V 10
IK ðV ; nK Þ = gK n4K ðV  VK Þ
0:01ðV + 34Þan ðV Þ =+ 34
1  exp V 10

V + 44bn ðV Þ = 0:125exp 
25
3IA ðV ; hA Þ = gA mAN ðV ÞhA ðV  VK Þ1mAN ðV Þ =+ 50
1 + exp V 20
hAN ðV Þ =1
1 + exp V +680

IKS ðV ; mKS Þ = gKS mKS ðV  VK Þ
1mKSN ðV Þ =+ 34
1 + exp V 6:5


tmKS ðV Þ =
exp

8



+ 55
+ exp V 30

+ 55
V 30

ICa ðV Þ = gCa m2CaN ðV ÞðV  VCa Þ1mCaN ðV Þ =1 + exp V +920
 



IKCa V ; Ca2 + = gKCa mKCaN Ca2 + ðV  VK Þ1
mKCaN Ca2 + =
!3:5
1+

KD

½Ca2 + 

INaPN ðV Þ = gNaP m3NaPN ðV ÞðV  VNa Þ
1mNaPN ðV Þ =55:7
1 + exp  V +7:7
IAR ðV Þ = gAR hARN ðV ÞðV  VK Þ

18

iScience 25, 103873, March 18, 2022

iSciencell

ArticleOPEN ACCESShARN ðV Þ =f ðV Þ =1
1 + exp V +4751


1 + exp V 20
2

IAMPA ðV ; sAMPA Þ = gAMPA sAMPA ðV  VAMPA Þ
INMDA ðV ; sNMDA ; xNMDA Þ = gNMDA sNMDA ðV  VNMDA Þ
IGABA ðV ; sGABA Þ = gGABA sGABA ðV  VGABA Þ
Please note that intrinsic (non-synaptic) ionic currents (e.g. ICa ) 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 t˛½0; 10. Sampling frequency was set to 1,000 Hz forparameter searches and 10,000 Hz for detailed analyses.
In some analyses, the conductance of leak channels gL was divided into the conductance of leak K + channels (gKL ) and the conductance of leak Na+ channels (gNaL ) as follows:

VL  VNaL
VK  VNaLVL  VK
gNaL = gL 3
VNaL  VK
gKL = gL 3

These definitions satisfy the following conditions:IL ðV Þ = gL ðV  VL Þ = gKL ðV  VK Þ + gNaL ðV  VNaL Þ
gL = gKL + gNaL
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 gNa ;gKS ;gA ;gAR ;gAMPA ;gNMDA , and gGABA 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 gNa ; gK ; gA ; gAR ; gAMPA ; gNMDA , and gGABA 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 (gL , gNa , gK , gA , gKS , gCa , gKCa , gNaP , gAR , gAMPA , gNMDA ,gGABA , and tCa ). All parameter sets were generated so that the logarithm of each parameter has a uniform
distribution. Let Uða; bÞ be a distribution whose probability density function ðf ðxÞÞ is as follows.

f ðxÞ =

8
>
<
>
:

1
ðlogb  logaÞx

ða%x%bÞ

0

ðx<a; x>bÞThe values of gL , gNa , gK , gA , gKS , gCa , gKCa , gNaP , and gAR (mS/cm2) were generated according to Uð0:01;
100Þ, gAMPA , gNMDA , and gGABA (mS) were generated according to Uð0:002; 20Þ, and tCa (msec) was generated according to Uð10; 1000Þ.
For each parameter set, the differential equations were solved by using integrate.odeint method in SciPy
1.1.0. Let V ðtÞ be a membrane potential obtained by solving the differential equation, and vi was defined asV ð5:000 + 0:001iÞ ði = 0; 1; .; 4999Þ. The lists of vi were first analyzed by the discrete Fourier transform after
normalization and detrending. The frequency whose Fourier power was the largest was defined as the peakfrequency. The sets Ml and ml was defined as follows.


Ml = j˛Z 1%j%4998; vj Rvj1 ; vj Rvj + 1

ml = j˛Z 1%j%4998; vj %vj1 ; vj %vj + 1

iScience 25, 103873, March 18, 2022

19

iSciencell

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


P = j ˛ Ml vj >  40
B = j ˛ ml vj <  80
P ðjPjRjBjÞ
B ðjBjRjPjÞ

S=

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

jX
max
8
>
Sj
S
=
>
>
>
>
j
>
<
X
Sj =
sjk
>
>
k
>


>
>
>
1%k% Sj  1: sjk  sjðk1Þ %50
sj1  sðj1ÞjSj1 j >50 ð1%j%jmax  1Þno8
B = l˛Z 1 + p%l%5000  p; sj1  p%l%sjjSj j + p
>
>
< k
!
jmax
1
When Sj R3, Bk , defined as follows:
X
>
>
:
p=
sj  sj1 = jSj jð1%k%kmax Þ1

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
jvi j ði = 0; 1; .; 4999Þ 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 thenumber 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 Uða; bÞ. The normalized gnorm was expressed as follows.

gnorm =

log g  log a
log b  log a

The PCA was conducted for gnorm . Let li be the contribution ratio and hci be the coefficient of component c˛C of PCi (1%i%jCj). The principal component score ti and the factor loading fci were calculated as follows.

ti =

X
c˛C

fci =

hci gðcÞnorm
pﬃﬃﬃ
li hci

Normalization of the oscillation density
The membrane potential, [Ca2+], and IC were normalized based on the rhythms of altering bursting and silent phases. First, differential equations were solved in t˛½0; 10n; n = 1. When kmax <7, n was multiplied bytwo, and the differential equations were solved again. This procedure was recursively repeated until n = 4.
The parameter sets which could not reach kmax R7 were excluded from the following analyses. Let bki be theith element of Bk , and the normalization was executed for each Nk defined as follows.


Nk = j ˛ Z bkjBk j % j % bðk + 1ÞjBk + 1 j

ð1 % k % 6Þ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, thedifferential equations were solved, and subsequently the firing pattern was classified based on the criteria

20

iScience 25, 103873, March 18, 2022

iSciencell

ArticleOPEN ACCESS

above mentioned. The percentage of the parameter set in which the SS firing pattern was not abolishedafter 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 PEX and outward currents PIX was defined as follows:

PEX = ðIX + jIX jÞ =X

ðIc + jIc jÞc˛CX
PIX = ðIX  jIX jÞ = ðIc  jIc jÞc˛C

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

10
bk jB j
6
Xk
1X
1
@
BIX =
IX ð5:000 + 0:001tÞ A
6 k = 1 bkjBk j  bk1 t = b
k1!
sk j S j
6
Xk
1X
1
SIX =
IX ð5:000 + 0:001tÞ
6 k = 1 skjSk j  sk1 t = s
k1

The average proportion of channel X’s outward current PSEX and inward current PSIX during the silent
phases, and the proportion of channel X’s outward current PBEX and inward current PBIX during bursting
phases were defined as follows.

PBEX = ðBIX + jBIX jÞ =XðBIc + jBIc jÞc˛C
XPBIX = ðBIX  jBIX jÞ = ðBIc  jBIc jÞc˛C

PSEX = ðSIX + jSIX jÞ =PSIX = ðSIX  jSIX jÞ =XðSIc + jSIc jÞc˛CXðSIc  jSIc jÞc˛CAnalysis 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 g0 as
follows.

g0 = g 3 ð1 + 0:001f Þðf ˛ Z;  100 % f % 100Þ

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

g01 = g1 3 ð1 + 0:005f1 Þ
g02 = g2 3 ð1 + 0:005f2 Þðf ˛ Z;  20 % f1 % 20Þðf ˛ Z;  20 % f2 % 20Þ

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

fo = 1 =6
1Xbðk + 1Þ1  bk1
6 k =1!

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

iScience 25, 103873, March 18, 2022

21

iSciencell

ArticleOPEN ACCESS 2+ 
d Ca2 +Ca=  aCa AICa ðV Þ 
=0
dt
tCa
The nullclines of V and mKS were described as follows.



 
dV
1
= 
IL ðV Þ + IKS ðV ; mKS Þ + ICa ðV Þ + IKCa V ; Ca2 + fixed + INaP ðV Þ = 0
dt
CdmKS mKSN ðV Þ  mKS==0
dt
t mKS ðV Þ
[Ca2+]fixed values in figures are described on the top of each figure.

Linear stability analysis
The function fs and gs were defined as follows.

fs ðV ; mKS Þ = 



 
1
IL ðV Þ + IKS ðV ; mKS Þ + ICa ðV Þ + IKCa V ; Ca2 + fixed + INaP ðV ÞC
mKSN ðV Þ  mKS
gs ðV ; mKS Þ =t mKS ðV Þ

When [Ca2+] was fixed at [Ca2+]fixed, the fixed point ðV  ; mKS Þ was calculated by solving fs ðV ; mKS Þ = gs ðV ;mKS Þ = 0. These nonlinear equations were solved optimize.fsolve method in SciPy 1.1.0. Subsequently, theJacobian matrix J was described as follows.

2
vfs  V ; mKS
6
6 vV
6
6
J=6
6
6
6

4 vgs  
V ; mKSvVvfs  V ; mKSvmKS
vgs  V ; mKS
vmKS

3
7
7
7
7
7
7
7
7
5

The stability of the fixed point was categorized as follows: (1) saddle (detðJÞ<0); (2) stable focus (detðJÞ> 0,trðJÞ>0, and trðJÞ2  4detðJÞ<0); (3) stable node (detðJÞ>0, trðJÞ>0, and trðJÞ2  4detðJÞ>0); (4)unstable focus (detðJÞ>0, trðJÞ<0, and trðJÞ2  4detðJÞ<0); (5) unstable node (detðJÞ>0, trðJÞ<0, andtrðJÞ2  4detðJÞ>0).

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

8>
fs ðV ; mKS Þ = 0>>><
gs ðV ; mKS Þ = 0
>
>


ba
>
>
: Ca2 + fixed = a +i1000

ði˛Z; 0%i%1000Þ

The stability of the solutions were determined based on the linear stability analysis. The solutions wereplotted 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.

22

iScience 25, 103873, March 18, 2022

iScience, Volume 25

Supplemental information

A design principle of spindle
oscillations in mammalian sleep
Tetsuya Yamada, Shoi Shi, and Hiroki R. Ueda

Supplemental Figures

Figure S1. Basic properties of the SS firing pattern in the AN and RAN model,
related to Figure 1(A) Distributions of the SS and SWS parameters in the AN model which are not shown in Figure
1E.
(B) Properties of intermittent firings of the SS and SWS firing pattern in the AN model.
(C) The typical SS firing pattern in the RAN model. The properties that are observed in the ANmodel are conserved in the RAN model.
(D) Standard deviation of the standardized V (left) and the average of the standardized [Ca2+ ](right) of the 1, 166 parameter sets collected by random parameter search for the RAN model.
(E) Properties of intermittent firings of the SS and SWS firing pattern in the RAN model.

Figure S2. Current analysis and mathematical analysis for the simplified models,related to Figure 2(A) Activation (left) and inactivation (right) voltage-clamp assay for the two types of delayedrectifier, voltage-gated K+ channels in the model.
(B) Distributions of the SS and SWS parameters in the RAN and SAN model respectively.
(C) Proportion of the current through leak Na+ channel (INaL ), voltage-gated Ca2+ channel(ICa ), and persistent Na+ channel (INaP ) of 1,166 parameter sets obtained by the random parameter searches in the RAN model.
(D) Distributions of the proportion of inward current during bursting and silent phases of the SS
and SWS firing pattern in the RAN and SAN model respectively.
(E) A trajectory over time and nullclines of V and nK of the representative parameter set of the
SWS firing pattern in the phase space. The coiled part of the trajectory represents bursting
phases, while the linear part represents silent phases.
(F) Phase planes and stream plot of the representative parameter set when the concentrationof intracellular calcium is low (5.0 µM) and high (6.0 µM). Each concentration corresponds to thetransition from a silent to a bursting phase and a bursting phase to a silent phase respectively.
(G) Fixed points (i.e., intersections of V and nKS nullclines) on the V -[Ca2+ ] plane (left) andtheir linear stability (right). The transition from a silent to a bursting phase occurs via the mergeof the stable focus and the saddle, meaning the fold (saddle-node) bifurcation. On the other
hand, the transition from a bursting phase to a silent phase occurs by the trajectory passing
through the saddle, meaning saddle homoclinic orbit bifurcation.

Figure S3. Balance of background inward/outward current controls oscillation
density and calcium setpoint of the SS firing pattern in two manners in the RAN
model, related to Figure 3(A) Relations between channel conductances (leak K+ channel (left) and leak Na+ channel(right)) and oscillation density of the SS firing pattern in the RAN model. Parameter sets are
the representatives chosen from minor subgroups of parameter sets shown in Figure 3B.
(B) Mean [Ca2+ ] for the decreasing and increasing parameter sets for each subgroup shown in
Figure 3A. The mean [Ca2+ ] in the major subgroup is smaller than that in the minor subgroupin both channels.
(C-F) The same analysis for the gKS and gNaP , corresponding to Figures 3A and 3B; S3A andS3B.
(G) Venn diagrams for the major subgroups of each channel (i.e., decreasing subgroups of theK+ channels and increasing subgroups of the Na+ channels) (left) and the minor subgroups ofeach channel.
(H) Changes in oscillation frequencies when two sets of channel conductances are changed
at the same time. As in Figure 3D, K+ channels and Na+ channels control oscillation density
competitively.

Figure S4. Changes in currents and mathematical structures in response toshifts in conductance values, related to Figure 3(A, B) Relative changes in current through each channel during each phase in response to theshifts in gKL (E) and gNaL values in each representative parameter set of the minor subgroup.
(C, D) Trajectories and bifurcation diagrams in response to the parameter shifts of representative parameter sets in minor subgroups of gKL (G) and gNaL (H), corresponding to (A). In minor
subgroups, when outward current increases (C, left) or inward current decreases (D, right), thebifurcation structures change, followed by increases in [Ca2+ ] setpoints and oscillation density.
When outward current decreases (C, right) or inward current increases (D, left), vice versa. In
such minor subgroups, the distance between the subcritical Andronov-Hopf bifurcation and thefold limit cycle bifurcation significantly increase.
(E-L) The same analysis for the gKS and gNaP , corresponding to Figures 3E-3H; S3A-S3D.

Figure S5. Current analysis for the SS firing pattern in the AN model, related to
Figure 4(A) Membrane potential, proportions of each channel’s inward and outward current, and proportion of inward and outward synaptic current of the representative parameter set in the ANmodel.
(B) Proportion of the current through leak Na+ channel (INaL ), voltage-gated Ca2+ channel(ICa ), and persistent Na+ channel (INaP ) of 1, 112 parameter sets obtained by the random parameter searches in the AN model.
(C) Distributions of the proportion of inward current during bursting and silent phases of the SS
and SWS firing pattern in the AN model.

Figure S6. Balance of background inward/outward current controls oscillation
density and calcium setpoint of the SS firing pattern in two manners in the ANmodel, related to Figure 4(A) Relations between channel conductances (leak K+ channel (left) and leak Na+ channel(right)) and oscillation density of the SS firing pattern in the AN model. Parameter sets are therepresentatives chosen from minor subgroups of parameter sets shown in Figure 4D.
(B) Mean [Ca2+ ] for the decreasing and increasing parameter sets for each subgroup shown in
Figure 4C. The mean [Ca2+ ] in the major subgroup is smaller than that in the minor subgroupin both channels.
(C-F) The same analysis for the gKS and gNaP , corresponding to Figures 4C and 4D; S6A andS6B.
(G) Venn diagrams for the major subgroups of each channel (i.e., decreasing subgroups of theK+ channels and increasing subgroups of the Na+ channels) (left) and the minor subgroups ofeach channel.
(H) Changes in oscillation frequencies when two sets of channel conductances are changed
at the same time. As in Figure 4F, K+ channels and Na+ channels control oscillation density
competitively.

Supplemental Tables

Table S1. Parameter values in the models, related to STAR Methods
C
A
VL
VNaL
VNa
VK
τhA
VCa
KD
VAMPA
VNMDA
VGABA
αCa
τAMPA
τsNMDA
τxNMDA
τGABA
gL
gNaL
gKL
gNa
gK
gA
gKS
gNaP
gAR
gCa
gKCa
gAMPA
gNMDA
gGABA
τCa

Membrane capacitanceArea of neuron
Leak reversal potential
Reversal potential of leak Na+ channels
Na+ reversal potential
K+ reversal potential
Time-constant of hA
Ca2+ reversal potential
Dissociation constant of Ca2+ -dependent K+
channels
Reversal potential of AMPA receptors
Reversal potential of NMDA receptors
Reversal potential of GABAA receptors
Coefficient of Ca2+ influx
Time constnant of sAMPA
Time constant of sNMDATime constant of xNMDATime constant of sGABA
Conductance of leak channels
Conductance of leak K+ channels
Conductance of leak Na+ channels
Conductance of voltage-gated Na+ channels
Conductance of voltage-gated K+ channels
Conductance of fast A-type K+ channels
Conductance of slowly inactivating K+ channels
Conductance of persistent Na+ channels
Conductance of inwardly rectifying K+ channels
Conductance of voltage-gated Ca2+ channels
Conductance of fast Ca2+ -dependent K+ channels
Conductance of AMPA receptors
Conductance of NMDA receptors
Conductance of GABAA receptors
Time constant of Ca2+ efflux

1 µF/cm2
0.02 mm2−60.95 mV
0 mV
55 mV
−100 mV
15 msec
120 mV
30 µM
0 mV
0 mV
−70 mV
0.5 µM/(nA · ms)
2 msec
100 msec
2 msec10 msec
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗
∗

gL [mS/cm2 ]gNa [mS/cm2 ]gK [mS/cm2 ]gA [mS/cm2 ]gKS [mS/cm2 ]gNaP [mS/cm2 ]gAR [mS/cm2 ]gCa [mS/cm2 ]gKCa [mS/cm2 ]gAMPA [µS]gNMDA [µS]gGABA [µS]
τCa [msec]

Figure 1B(SWS)
0.03573
12.2438
2.61868
1.79259
0.0350135
0.0717984
0.0166454
0.0256867
2.34906
0.513425
0.00434132
0.00252916
121.403

Figure 1B(SS)
1.073449
4.934444
0.122135
0.013062
56.409369
12.034643
0.174262
0.192476
0.245811
0.859253
0.048610
0.515353
828.725007

Figure
S1C
8.222564
00097.524247
11.690555
00.766501
0.459510000565.240040
1.40603000019.138263
6.63643801.914677
0.296167000884.719189

Figure 2A

0.076208019.258326
000.697291
00.084111
14.093700
000709.874820

Figure S2E

Table S2. Channel conductance in the models (1), related to STAR Methods

gL [mS/cm2 ]gNa [mS/cm2 ]gK [mS/cm2 ]gA [mS/cm2 ]gKS [mS/cm2 ]gNaP [mS/cm2 ]gAR [mS/cm2 ]gCa [mS/cm2 ]gKCa [mS/cm2 ]gAMPA [µS]gNMDA [µS]gGABA [µS]
τCa [msec]

Figure 3A(left)
0.16702400018.009547
3.601729
00.437785
0.615421000611.720170

Figure 3A(right)1.50793200033.860581
4.197253
00.511429
0.604971000692.292870

Figure S3A(left)
0.27677600038.062753
9.200761
02.525039
0.056078000591.128972

Figure S3A(right)0.61399800069.305695
21.476862
00.964186
0.212878000593.107158

FigureS3C (left)
4.512591
00057.452328
6.952845
00.783477
0.350631000580.502638

FigureS3C (right)
0.16810200024.998703
5.35043300.240259
5.158678000470.158046

Figure S3E(left)
0.86890700076.655055
25.006008
01.178154
0.271071000605.562460

Table S3. Channel conductance in the models (2), related to STAR Methods
Figure S3E(right)1.72966900082.533728
18.148719
03.467343
0.338148000455.608044

gL [mS/cm2 ]gNa [mS/cm2 ]gK [mS/cm2 ]gA [mS/cm2 ]gKS [mS/cm2 ]gNaP [mS/cm2 ]gAR [mS/cm2 ]gCa [mS/cm2 ]gKCa [mS/cm2 ]gAMPA [µS]gNMDA [µS]gGABA [µS]
τCa [msec]

Figure 4C(left)
0.417727
1.583492
0.3011300.151947
11.331956
1.603688
0.016290
0.170668
20.433243
0.098587
2.065959
0.053074
588.964631

Figure 4C(right)1.870606
2.628010
3.342107
0.023800
99.700384
16.441094
1.143548
1.293300
0.397588
0.835785
0.160560
0.517354
243.598190
0.223534
0.012612
1.509700
0.092186
31.319209
7.291272
0.043571
0.550217
0.237943
0.263511
0.032661
0.005638
851.362910

Figure S5A

Figure S5A(left)
4.656359
11.755551
0.255884
0.155798
81.681239
19.976798
0.019561
0.155672
0.976006
0.002020
0.041042
0.167732
334.598243

Figure S6A(right)4.656359
11.755551
0.255884
0.155798
81.681239
19.976798
0.019561
0.155672
0.976006
0.002020
0.041042
0.167732
334.598243

FigureS6C (left)
1.396760
3.312551
1.410942
0.206154
75.115235
11.182473
0.384665
0.056088
16.918662
0.007833
0.134228
0.095727
571.740400

FigureS6C (right)
2.745233
0.881557
0.316212
0.017476
64.192433
7.541997
0.248578
0.421745
1.822533
0.336310
0.016536
0.006981
823.122356

Table S4. Channel conductance in the models (3), related to STAR Methods
Figure S6E(left)
0.080461
0.118620
0.471621
0.033450
8.221155
3.639074
0.023680
0.878336
0.182650
0.004408
0.091796
0.038631
138.631890

Figure S6E(right)1.256775
0.864784
0.268804
4.827961
41.349411
6.677935
0.048519
4.236638
0.493004
0.005055
2.096988
0.021540
36.494085

Table S5. Initial values for the simulations, related to STAR Methods
V
hNa
nK
hA
mKS
sAMPA
xNMDA
sNMDA
sGABA
[Ca2+ ]

−45 mV
0.045
0.54
0.045
0.34
0.01
0.01
0.01
0.01
1.0 µM

Table S6. Hit rate of the SS firing pattern in each model, related to Figure 1
Components in the model
gCa , gNMDA , τCa , gKCa , gL , gNa , gNaP , gAMPA , gK ,gKS , gA , gAR , gGABA (Full model)τCa , gKCa , gL , gNaP , gKS (Model 1)
gNMDA , τCa , gKCa , gL , gNaP , gKS (Model 2)
gCa , τCa , gKCa , gL , gNaP , gKS (RAN model)
gCa , τCa , gKCa , gL , gNaP , gK (SAN model)

All

SS

Hit rate [%]55, 715, 098

1, 1122.00 × 10−3

10, 385, 607
10, 869, 316
23, 146, 942
62, 808, 489

0696
1, 1660

0.00
6.40 × 10−3
5.04 × 10−30.00