Leak potassium channels regulate sleep duration
Kensuke Yoshidaa,1, Shoi Shia,b,c,1, Maki Ukai-Tadenumab,c, Hiroshi Fujishimab, Rei-ichiro Ohnoa, and Hiroki R. Uedaa,b,c,2
aDepartment of Systems Pharmacology, Graduate School of Medicine, The University of Tokyo, 113-0033 Tokyo, Japan; bLaboratory for Synthetic Biology,
RIKEN Center for Biosystems Dynamics Research, 565-0871 Osaka, Japan; and cInternational Research Center for Neurointelligence, The University of Tokyo
Institutes for Advanced Study, The University of Tokyo, 113-0033 Tokyo, Japan

Edited by Joseph S. Takahashi, Howard Hughes Medical Institute and University of Texas Southwestern Medical Center, Dallas, TX, and approved August 21,
2018 (received for review April 15, 2018)

|

|

|

sleep computational model slow-wave sleep firing pattern
leak potassium channel Ca2+-dependent hyperpolarization pathway|

S

leep is a physiological phenomenon widely conserved among
species (e.g., nematodes, jellyfish, flies, and mammals). The
sleep/wake cycle occurs with a time scale of hours in primatesand one of minutes in rodents. Although the daily sleep duration
varies among species, it is conserved within a species (1), suggesting that a genetic mechanism is involved in its regulation.
Studies on the cellular and molecular nature of sleep/wake cycle
regulation have revealed some sleep/wake-promoting circuits
and secreted molecules (2, 3). Indeed, the synaptic projection
from sleep/wake-promoting nuclei in the brainstem to the cerebral cortex has been found to regulate the sleep/wake cycle (3).
However, it is still elusive how the daily sleep duration and thesleep/wake cycle are genetically regulated.
Whereas the regulation of sleep duration has relatively slow
dynamics (i.e., from minutes to hours), the EEG during eachstate of the sleep/wake cycle exhibits faster dynamics (i.e., frommilliseconds to seconds). For example, nonrapid eye movement(NREM) sleep is accompanied by a high-amplitude, lowfrequency (0.5–4 Hz) pattern of EEG, whereas other states(wakefulness or REM sleep) are associated with low-amplitude
high frequencies, which are respectively generated by the synchronous or irregular firing of cortical neurons (4–7). The
membrane potential of cortical neurons exhibits a depolarized
bursting pattern (awake firing pattern) during wakefulness or
REM sleep, whereas it alternates between depolarized “up”
states and hyperpolarized “down” states during NREM sleep.
This firing pattern during NREM sleep is represented as theslow-wave sleep (SWS) firing pattern in this study. Several
computational studies have attempted to reveal the mechanism
underlying these firing patterns (8–12). These studies led to theidea that the transition from the up state to the down state occurs in cortical neurons as a result of their regulation by neuralnetworks (e.g., the thalamocortical network). However, although
the role of synaptic connections in the sleep/wake cycle has been
revealed, the intracellular mechanisms of sleep/wake cycle regulation are still elusive.

www.pnas.org/cgi/doi/10.1073/pnas.1806486115

A recent study proposed a new mechanism for the sleep/wakecycle regulation by focusing on intracellular mechanisms of the
SWS firing pattern (13). The study succeeded in computationally
recapitulating the SWS firing pattern in an averaged-neuron(AN) model, whereby a cortical neuron transmitted its output
to equivalent neurons (i.e., its output returns to itself as input),suggesting that an intrinsic mechanism, the regulation of a Ca2+dependent hyperpolarization pathway, is involved in generating
the SWS firing pattern (13). Notably, the components of the Ca2+dependent hyperpolarization pathway, which are predicted to beimportant for generating the SWS firing pattern, are also involved
in regulating sleep duration (13), suggesting that the regulation of
the SWS firing pattern (fast dynamics) and sleep duration (slowdynamics) could share a common mechanism (13, 14). In otherwords, the changes in cortical firing patterns might not just accompany the sleep/wake cycle but rather cause it. Therefore, determining the mechanism by which the SWS firing pattern isgenerated could reveal clues about the regulation mechanism of
the sleep/wake cycle.
Although the AN model succeeded in revealing the role of theCa2+-dependent hyperpolarization pathway in generating the
SWS firing pattern, its detailed mathematical structure was still
unknown. Because the major difference between the SWS andawake firing patterns is the presence or absence of the downstate, it is important to elucidate how the transition from the up
to the down state occurs. This question can be addressed by mathematical analyses, for example, to analyze the currents responsible
SignificanceTo address molecular mechanisms regulating sleep duration, a
simple computational model of a cortical neuron [simplified
averaged neuron (SAN) model], which recapitulates the electrophysiological characteristics of slow-wave sleep (SWS) andwakefulness, is developed in this study. Comprehensive bifurcation and detailed mathematical analyses predicted thatleak K+ channels play a role in generating the cortical electrophysiological characteristics of SWS, leading to a hypothesisthat leak K+ channels play a role in the regulation of sleep
duration. We comprehensively generated and analyzed 14
knockout mice of the leak K+ channel family, which demonstrated that impairment of the leak K+ channel (Kcnk9) decreases sleep duration. The results confirm the validity of the
SAN model and suggest a molecular mechanism regulating
sleep duration.
Author contributions: K.Y., S.S., and H.R.U. designed research; K.Y., S.S., M.U.-T., H.F.,
R.-i.O., and H.R.U. performed research; K.Y., S.S., and H.R.U. contributed new reagents/
analytic tools; K.Y. and S.S. analyzed data; and K.Y., S.S., and H.R.U. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
Published under the PNAS license.
1

K.Y. and S.S. contributed equally to this work.

2

To whom correspondence should be addressed. Email: uedah-tky@umin.ac.jp.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.
1073/pnas.1806486115/-/DCSupplemental.
Published online September 17, 2018.

PNAS | vol. 115 | no. 40 | E9459–E9468

NEUROSCIENCE

A primary goal of sleep research is to understand the molecularbasis of sleep. Although some sleep/wake-promoting circuits andsecreted substances have been identified, the detailed molecularmechanisms underlying the regulation of sleep duration havebeen elusive. Here, to address these mechanisms, we developed asimple computational model of a cortical neuron with five channels
and a pump, which recapitulates the cortical electrophysiologicalcharacteristics of slow-wave sleep (SWS) and wakefulness. Comprehensive bifurcation and detailed mathematical analyses predictedthat leak K+ channels play a role in generating the electrophysiological characteristics of SWS, leading to a hypothesis that leak K+
channels play a role in the regulation of sleep duration. To test
this hypothesis experimentally, we comprehensively generated
and analyzed 14 KO mice, and found that impairment of the leakK+ channel (Kcnk9) decreased sleep duration. Based on theseresults, we hypothesize that leak K+ channels regulate sleep duration in mammals.

for the transition from the up to the down state, but the ANmodel, which contains 13 components, is too complicated to
conduct such analyses. Studies of other rhythmic bursting patterns
that are also characterized by up and down states (e.g., the preBötzinger complex or hippocampal pyramidal neurons) have indicated that approximately four to nine channels are sufficient togenerate them (15–17), raising the possibility that some of the
13 components in the AN model are redundant for generating the
SWS firing pattern. Indeed, another study indicated that NMDA
receptors and voltage-gated Ca2+ channels have redundant roles
and that GABA receptors are minimally involved in generating
the SWS firing pattern (13). These reports suggest that the ANmodel might be simplified to conduct mathematical analyses.
Therefore, to reveal the mechanism for generating the down state
and the SWS firing pattern, here we developed a simplified modeland used it to perform detailed mathematical analyses.
Results
Construction of the Simplified AN Model, Which Recapitulates the
SWS Firing Pattern. To analyze the SWS firing pattern occurringin the cortical neurons during SWS (6), we simplified the ANmodel (13) by eliminating its components: channels, receptors,
and pump (Fig. 1A). We first sought to identify which components were redundant for generating the SWS firing pattern byconducting random parameter searches and calculating the hit
rates (i.e., the percentage of parameter sets that showed the SWS
firing pattern among all randomly generated parameter sets)
under KO of the components (i.e., setting the conductance tozero). We compared the hit rate of each model, and when the hit
rate of a model is more than 100 times lower than the other, we
assumed that the model is highly unlikely and inappropriate. If
eliminating components had only a minor effect on the hit rate,
those components were considered redundant for the generationof the SWS firing pattern and were eliminated to simplify the model.
Because the voltage-gated Ca2+ channels and the NMDA receptors play redundant roles, the elimination of the NMDA receptors was reported to only slightly affect the hit rate (13).
Therefore, the NMDA receptors were eliminated first frommodel 0 to generate model 1. We then classified the channels inmodel 1 into five categories: Ca2+-related channels (gCa, gKCa),
leak channels (gL), Na+ channels (gNa, gNaP, gAMPA), K+ channels(gK, gKS, gA, gAR), and Cl− channels (gGABA). Because both of theCa2+-related channels (gCa, gKCa) are reported to be essential for
generating the SWS firing pattern under NMDA receptor KO(13), we conducted the random parameter searches and compared the hit rates under the KO of each category except for the
Ca2+-related channels. The parameter searches under KO ofleak channels exhibited a much lower hit rate (∼1.9 × 10−5%)
than model 1 (∼4.9 × 10−3%; Fig. 1B), indicating that the leakchannels are considered essential for generating the SWS firing
pattern and should therefore be preserved in the model. The hit
rate under KO of all three Na+ channels was zero (Fig. 1C),
indicating that at least one of the Na+ channels is necessary for
generating the SWS firing pattern. Among the three Na+ channels, the parameter searches under KO of voltage-gated Na+channels exhibited a higher hit rate (∼1.2 × 10−3%) than thoseunder KO of persistent Na+ channels (∼4.6 × 10−4%) or of
AMPA receptors (∼4.5 × 10−4%; Fig. 1C), indicating that thevoltage-gated Na+ channels are less important for generating the
SWS firing pattern than the other two. Because the persistentNa+ channel is an intrinsic channel that was found to be important for generating the rhythmic bursting in other studies (15,16, 18, 19), the persistent Na+ channels were chosen to be preserved in the model. In addition, the parameter searches under
KO of all four K+ channels (gK, gKS, gA, gAR) or the GABA receptors exhibited hit rates (∼2.3 × 10−3% or ∼6.7 × 10−3%, respectively) in the same range as model 1 (Fig. 1 D and E),
indicating that the K+ channels and the GABA receptors play
E9460 | www.pnas.org/cgi/doi/10.1073/pnas.1806486115

limited roles in generating the SWS firing pattern. Based onthese results, the voltage-gated Na+ channels, AMPA receptors,
four K+ channels, and GABA receptors were considered candidates for elimination from model 1.
However, the hit rate under KO of the voltage-gated Na+channels, AMPA receptors, four K+ channels (gK, gKS, gA, gAR),
and GABA receptors (model 2) was zero (Fig. 1F), indicatingthat model 2 lacks some components needed to generate the
SWS firing pattern. To find the sufficient model, we generated
seven models by adding each eliminated component back tomodel 2 and conducted the random parameter searches. The
models consisting of model 2 plus gNa, model 2 plus gAMPA,
model 2 plus gK, model 2 plus gA, and model 2 plus gGABA
exhibited nearly normal hit rates (∼1.6 × 10−3%, ∼3.5 × 10−4%,∼1.2 × 10−3%, ∼3.9 × 10−4%, and ∼2.6 × 10−3%, respectively;
Fig. 1G); hence, these five models were further considered as
candidates for the simplified model.
As the simplified model should exhibit the important characteristics of model 1, in which down-regulation of the Ca2+dependent hyperpolarization pathway causes the transition fromthe SWS to awake firing pattern (SI Appendix, Fig. S1A), we
conducted a bifurcation analysis (i.e., a gradual change in parameter values) in the five models. In the bifurcation analysis of
1,008 parameter sets in the model 2 plus gK model, the Ca2+dependent hyperpolarization pathway was preserved, whereby
down-regulating the conductance of voltage-gated Ca2+ channels
or Ca2+-dependent K+ channels or reducing the time constant of
Ca2+ efflux caused the transition from the SWS to awake firing
pattern in 188 (18.7%), 816 (81.0%), and 815 (80.9%) parameter
sets, respectively (Fig. 1H). However, this transition derived fromthe Ca2+-dependent hyperpolarization pathway was not observedin the other four models (Fig. 1 I and J and SI Appendix, Fig. S1
B and C). These results indicated that the model 2 plus gK model
preserved more of the characteristics of model 1 than the otherfour models. Therefore, we chose the model 2 plus gK model as
the simplified AN (SAN) model (Fig. 1K and SI Appendix, SI
Materials and Methods, SAN Model), which contains six components and three variables in differential equations, as opposed to
the 13 components and 10 variables in differential equations in
the AN model. The SAN model exhibited a similar membranepotential and intracellular Ca2+ concentration as model 1 (Fig. 1
L and M and SI Appendix, Fig. S1 D–H), suggesting that the SAN
model preserved the characteristics of model 1.
Leak K+ Channels Play a Role in Generating the SWS Firing Pattern
Based on the SAN Model. The histogram of the 1,008 parameter

sets with the SWS firing patterns in the SAN model showed bimodality of the gK and gCa values, implying that these parameter
sets could be divided into two clusters (SI Appendix, Fig. S2A). In
fact, the logarithmic plot of gK against gCa exhibited two clusters(Fig. 2A): cluster 1 with 480 parameter sets and cluster 2 with
528 parameter sets. These clusters were also obtained by principal component analysis (SI Appendix, Fig. S2B), supporting
their validity. We then determined the representative parameter
set in each cluster by performing the density estimation (SIAppendix, Table S2). The representative parameter set exhibited
the SWS firing pattern in each cluster (Fig. 2B). We analyzed
cluster 1 and then confirmed that cluster 1 and cluster 2 represent the same results. As the major difference between the SWS
and awake firing patterns is whether the down state exists, we
focused on the currents of channels during the down state of the
SWS firing pattern. The analysis of current with the representative parameter sets revealed that the leak K+ channels andCa2+-dependent K+ channels were two major sources of outward
current in the down state (Fig. 2C and SI Appendix, Fig. S2C).
These results suggested that the leak K+ channels, in addition to
the Ca2+-dependent K+ channels, which were previously studied(13), are important for generating the down state in the SWS
Yoshida et al.

A
Voltage-gated Ca channel
NMDA receptor
2+
Ca -pump/exchanger
2++
Ca -dependent K channel

ΤCa
gKCa
gL
gNa
gNaP
gAMPA

+

Na, K Leak channel++
Na
Voltage-gated Na channel++
Persistent Na channel
Na+
Na
AMPA receptor++
K
Voltage-gated K channel++
K
Slowly inactivating K channel++
K
Fast A-type K channel++
K
Inwardly rectifying K channel
Cl
GABAA receptor

Hit rate (%)

0.004
0.003

0.003
0.002
0.001

2

1

el
od

M

el
od
M
40
20

0

Model 2M

SAN Model
Voltage-gated Ca channel

2+

Ca -pump/exchanger
2++
Ca -dependent K channel

2+

Ca
+
K

+

+

Na, K Leak channel+
Persistent Na+ channel
Na++
K
Voltage-gated K channel

-2

g Ca
Τ Ca
g KCa
gL
g NaP

Parameter sets

2+

2+

Ca

gK

L

Membrane Potential
n = 1,008
-1
0
1

2

-2

Intracellular Calciumn = 1,008
-1
0
1

2

Parameter sets

K

MembranePotential(mV)

SAN Model
-20

10

8

2+

Standardized Ca

Intracellular
Calcium
(µM)

0

13

Standardized
MembranePotential

3

-80

0

-3
0

Time (sec)

10

0

3

0

-3
6
0
Bursting cycles

6Fig. 1. Construction of the SAN model by simplifying the AN model. (A) Schematic illustration of the AN model (Left). Table shows the channels included ineach model (Right). (B–E) Hit rates of model 1 or model 1 under the KO of leak channels (B), Na+ channels (C), K+ channels (D), or GABA receptors (E). (F and
G) Hit rates of model 1 or model 2 (F) or model 1 or model 2 plus each ion channel condition (G). (H–J) Summary of the bifurcation analysis in each condition:model 2 plus gK (H), model 2 plus gNa (I), or model 2 plus gGABA model (J). Conductance and time constants related to the Ca2+-dependent hyperpolarization
pathway are in red. (K) Schematic illustration of the SAN model. (L) Typical patterns of the membrane potential and intracellular Ca2+ concentration in the
SAN model. (M) The standardized membrane potential and intercellular Ca2+ concentration of the 1,008 parameter sets obtained by parameter searches in
the SAN model. The line and shadow (Lower) indicate the mean and SD, respectively.

Yoshida et al.

PNAS | vol. 115 | no. 40 | E9461

NEUROSCIENCE

0 0 0

G
AB
A

0

g

a

0 0 0

60

g
L
g
N
g aP

20

80

a

40

100

KC

60

Model 2
+gGABA

ΤC
ga

Bifurcation from SWS to Awake (%)

80

KC

-

J

100

0

0

0.000
+

Model 2
+gNa

ga
L
g
N
a
gP

0

g

G
AB
A

AR

0.004

+g

A

+g

KS

+g

+g

K

PA

+g

+g

AM

a

20
0

I

a

0

1
N

el

+g

od
M

0.005

0.000
gGABA

+
+
+
-

+
+
+

C

0

+
+
+

g

0.001

40

K

0.002

60

ga
L
g
N
a
gP

0.003

80

a

0.004

100

KC

Bifurcation from SWS to Awake (%)

0.005

+
+
+

Bifurcation from SWS to Awake (%)

H

+
+
+
+
Model 2+ gK
(SAN Model)

ΤC
ga

0.000
gK
gKS
gA
gAR

+
+
-

ΤC
ga

+
+

C

+
+

-

G

Hit rate (%)

F 0.005

0.001
0

0.000
g Na
+
g NaP +g AMPA +

0.000

+
+
+
+
+
+
-

0.002

0.001

-

+

0.010
0.005

1.

0.000
gL

+
+
+
+
+
-

a

0.002

9×

0.001

0.003

SAN Model

0.006

N

0.002

0.004

Model 2

E 0.007
Hit rate (%)

0.003

ΤCa
gKCa
gL
gNa
gNaP
gAMPA
gK
gKS
gA
gAR
gGABA

Model 1(NMDAR KO)
+
+
+
+
+
+
+
+
+
+
+
+

0.015
Hit rate (%)

Hit rate (%)

0.004

gCa
gNMDA

D

C 0.005

1 0 -5Hit rate (%)

B 0.005

gK
gKS
gA
gAR
gGABA

Model 0(Original)
+
+
+
+
+
+
+
+
+
+
+
+
+

C

+

gCa
gNMDA

2+

2+

Ca
2+
Ca
2+
Ca
+
K

A

B

C

D

E

F

H

G

I

Fig. 2. Leak K+ channels play a role in generating the SWS firing pattern based on the SAN model. (A) Scatter plot of the gK and gCa values for all
1,008 parameter sets obtained by the parameter searches. All parameter sets were divided into two clusters, cluster 1 (blue) and cluster 2 (orange), according
to the gK and gCa values. (B) Membrane potential and intracellular Ca2+ concentration of the representative parameter set in each cluster over time. (C)Membrane potential, proportion of each channel’s inward current, and proportion of each channel’s outward current of the representative parameter set ineach cluster over time. (D) Bifurcation analysis of the representative parameter set in each cluster. The definitions and examples of SWS, awake, resting, andslow-wave with few spikes are described in SI Appendix, Fig. S2E and Table S5, respectively. (E) Two-dimensional bifurcation analysis of the representative
parameter set in cluster 1. (F) Three- (Upper) and two-dimensional (Lower) plots of the trajectory of the representative parameter set in cluster 1. (G) Phase
planes of the representative parameter set in cluster 1 in cases in which the Ca2+ concentration is fixed at 7 μM, 8.5 μM, 10.6 μM, or 12.5 μM. The red, green,
and blue lines indicate the trajectory of the differential equation, V nullcline, and nK nullcline, respectively. (H) Membrane potential and phase plane of theparameter set whose gKCa was lower (Left) or whose gKCa and gKL were lower (Right) than the representative parameter set in cluster 1. (I) Summary of thebifurcation analysis in cluster 1. Conductance of the leak K+ channels is in red.

E9462 | www.pnas.org/cgi/doi/10.1073/pnas.1806486115

Yoshida et al.

Leak K+ Channels Play a Role in Generating the Down State of the
SWS Firing Pattern Based on Detailed Mathematical Analyses. To

reveal the role of the leak K+ channels in generating the SWS
firing pattern, we focused on their role in the transition from theup to the down state during the SWS firing pattern by analyzingcluster 1 and then confirmed that cluster 1 and cluster 2 represent the same results. First, we plotted the trajectory of the SWS
firing pattern in the phase space. Whereas the AN model has
10 variables, the SAN model has only three variables (V, themembrane potential; nK, the dimensionless quantity associatedwith the activation of the voltage-gated K+ channels; and [Ca2+],
the intracellular calcium concentration) in the differential
equation. Therefore, the trajectory of the SAN model could be
represented in the 3D space (Fig. 2F). The coiled part and thestraight part of the trajectory correspond to the up state and thedown state, respectively (Fig. 2F). The Ca2+ concentration
gradually increased or decreased in the up or down state, respectively (Fig. 2F). To understand how this trajectory was
generated, we plotted the V nullcline, nK nullcline, and [Ca2+]
nullcline with the trajectory (SI Appendix, Fig. S3C). As thenullcline is a set of the points whose differential of a certain
direction is zero, its position describes which direction each pointin the space moves. Therefore, the nullcline enables us to understand the approximate dynamics of the differential equation.
To analyze this phase space in more detail, we then plotted thenullclines and the trajectory at four fixed Ca2+ concentrations(Fig. 2G). At the low Ca2+ concentration (7 μM in cluster 1), thetrajectory converged to the stable limit cycle (Fig. 2G), whichcorresponded to the up state of the SWS firing pattern. At the
middle Ca2+ concentrations (8.5 and 10.6 μM in cluster 1), there
were at least two stable states: a stable point (Fig. 2G, orange
point) and a stable limit cycle, which correspond to the downstate and the up state, respectively (Fig. 2G). At these Ca2+
concentrations, the stable limit cycle was not able to reach theorange point, and therefore still existed (Fig. 2G). At the highCa2+ concentration (12.5 μM in cluster 1), the stable limit cycle
disappeared, and the trajectory converged to the orange point(Fig. 2G). In summary, the phase planes of the SAN model
revealed that the transition between the up and the down state iscontrolled by a calcium-dependent bistable system, which is unYoshida et al.

derstood by the regulation of two nullclines, the V nullcline andthe nK nullcline. In other words, the bistability at the middle
Ca2+ concentrations enabled oscillation between the up anddown states.
To demonstrate the role of the leak K+ channels in regulating
the nullclines, and hence the transition from the up to the downstate, we down-regulated the conductance of leak K+ channels
and plotted the nullclines and the trajectory. This analysis
showed that down-regulating the conductance of leak K+ channels caused the V nullcline to shift in a positive direction on thenK axis, which prevented the transition from the up to the downstate and generated the awake firing patterns (Fig. 2H). Interestingly, down-regulating the conductance of Ca2+-dependent
K+ channels, the other major source of outward current during
the down state (Fig. 2C), also caused the V nullcline to shift in apositive direction on the nK axis, preventing the transition fromthe up to the down state, and generated the awake firing pattern(SI Appendix, Fig. S3G). These results indicated that the leak K+
channels and the Ca2+-dependent K+ channels cooperatively
control the shift of the V nullcline and hence the intersection
point of the V nullcline and the nK nullcline, thereby determining
the transition from the SWS to awake firing pattern by regulating
the transition from the up to the down state. Furthermore, theadditive control of the V nullcline by the leak K+ channels andCa2+-dependent K+ channels could explain the results of the 2D
bifurcation analysis against these channels (Fig. 2E).
The voltage-gated Ca2+ channels and the time constant of
Ca2+ efflux cooperatively regulate the Ca2+ concentration, which
also controls the transition from the SWS to awake firing pattern. This can also be elucidated by the control of the V nullcline
by the Ca2+-dependent K+ channels. Because the Ca2+-dependent
K+ channels are regulated by the Ca2+ concentration, the voltagegated Ca2+ channels and the time constant of Ca2+ efflux could
affect the transition from the up to the down state by regulating
the Ca2+-dependent K+ channels and thus the V nullcline (Fig.
2G). Therefore, the voltage-gated Ca2+ channels, the time constant of Ca2+ efflux, and the leak K+ channels could cooperatively
control the transition from the up to the down state and hencefrom the SWS to awake firing pattern.
To investigate the generality of these findings, we conducted abifurcation analysis of 480 parameter sets in cluster 1 of the SAN
model, in which the leak channels are divided into the leak K+
and Na+ channels. Down-regulating the conductance of leak K+
channels caused the transition from the SWS to awake firing
patterns in 96 parameter sets (20.0%) of cluster 1 (Fig. 2I), indicating that the leak K+ channels are involved in generating the
SWS firing pattern. Notably, phase plane analyses against the
representative parameter set and bifurcation analyses in cluster
2 also supported this finding (SI Appendix, Fig. S3 C–J).
Leak K+ Channels Play a Role in Generating the Down State of the
SWS Firing Pattern Based on the AN Model. To evaluate the role ofthe leak K+ channels in generating the SWS firing pattern undersynaptic connections (i.e., a model containing the AMPA receptors, the NMDA receptors, and the GABA receptors), we
conducted the current and bifurcation analyses in the AN model,
which contains 13 components including synaptic currents. In theprevious study, the leak K+ channels and the leak Na+ channels
were represented as one type of leak channels (13). To elucidate
the role of the leak K+ channels, we analyzed currents of the leakK+ and Na+ channels together with the other channels in theoriginal AN model, the KO of NMDA receptors (i.e., setting the
conductance of the NMDA receptors to zero in the original ANmodel) model, and the KO of voltage-gated Ca2+ channels (i.e.,setting the conductance of the voltage-gated Ca2+ channels tozero in the original AN model) model. The percentages of thecurrent of the leak K+ channels and the Ca2+-dependent K+
channels were high in the original AN model (Fig. 3 A and B and

PNAS | vol. 115 | no. 40 | E9463

NEUROSCIENCE

firing pattern. To investigate the role of the leak K+ channels ingenerating the SWS firing pattern, we conducted the bifurcation
analysis with the representative parameter sets. However, downregulating the conductance of the leak K+ channels alone did notcause a transition between the SWS and awake firing patterns(Fig. 2D and SI Appendix, Fig. S2 D and E). Considering that theleak K+ channels and the Ca2+-dependent K+ channels play
redundant roles in the down state (Fig. 2C and SI Appendix, Fig.
S2C), the leak K+ channels might regulate the threshold of thetransition derived from the Ca2+-dependent hyperpolarization
pathway. Therefore, we next conducted a 2D bifurcation analysis, in which two of four parameters (the conductance of the leakK+ channels, Ca2+-dependent K+ channels, and voltage-gated
Ca2+ channels; and the time constant of Ca2+ efflux) werechanged gradually and simultaneously. This analysis showed that
down-regulating the conductance of leak K+ channels caused atransition from SWS to awake firing patterns under the downregulated conductance of Ca2+-dependent K+ channels or voltagegated Ca2+ channels or the reduced time constant of Ca2+
efflux (Fig. 2E and SI Appendix, Fig. S3A). These results suggested that the leak K+ channels change the threshold for thetransition induced by the Ca2+-dependent hyperpolarization
pathway and are therefore important for generating the SWS firing pattern. Notably, the current and bifurcation analyses against
the representative parameter set in cluster 2 also supported this
scenario (Fig. 2 C and D and SI Appendix, Figs. S2 C and D
and S3B).

SI Appendix, Fig. S4 A–C), as with the SAN model (Fig. 2C andSI Appendix, Fig. S2C). Interestingly, the percentages of thecurrent of the leak K+ channels and Ca2+-dependent K+ channels in the down state were relatively high: 29.4% and 58.0%,
respectively, in the original AN model, 44.8% and 34.5% in the
KO of NMDA receptors model, and 30.2% and 58.2% in the KOof voltage-gated Ca2+ channels model (Fig. 3C and SI Appendix,
Fig. S4D). On the contrary, the percentages of the current of theleak K+ channels and Ca2+-dependent K+ channels in the up
state were relatively low: 12.7% and 18.6%, respectively, in theoriginal AN model, 8.7% and 8.2% in the KO of NMDA receptors model, and 13.8% and 10.2% in the KO of voltage-gated
Ca2+ channels model (Fig. 3D and SI Appendix, Fig. S4E). These
results suggested that the leak K+ channels and the Ca2+dependent K+ channels are involved in the down state rather
than the up state. To elucidate the role of the leak K+ channels
in the transition between the SWS and awake firing patterns, we
conducted the bifurcation analysis in the original AN model, theKO of NMDA receptors model, and the KO of voltage-gated
Ca2+ channels model. The result showed that down-regulating
the conductance of leak K+ channels and Ca2+-dependent K+
channels led to the transition from SWS to awake firing patterns
in 270 parameter sets (24.7%) and 828 parameter sets (75.9%) of
the original AN model, 395 parameter sets (33.8%) and
1,024 parameter sets (87.7%) of the KO of NMDA receptors
model, and 140 parameter sets (14.8%) and 851 parameter sets
(89.8%) of the KO of voltage-gated Ca2+ channels model (Fig.
3E). These results indicated that the Ca2+-dependent K+ channels are important (13) and that the leak K+ channels are also
involved in generating the SWS firing pattern. These results fromthe AN model also suggest that the currents of the leak K+
channels and Ca2+-dependent K+ channels cause the transition
from the up to the down state, and hence the transition between
the SWS and awake firing patterns. In support of these findings,
the KO of NMDA receptors model, which showed a higher
proportion of the leak K+ channels in the down state (44.8%;
Fig. 3C), exhibited a higher proportion of the transition from the
SWS to awake firing patterns (33.8%; Fig. 3E) in down-regulating
the conductance of the leak K+ channels than the other twomodels. These results indicated that the leak K+ channels areespecially involved in the down state and generated the SWS firing
pattern in the SAN and AN models.
Impairment of Leak K+ Channels Decreases Sleep Duration in Mice.

Analyses of the SAN and AN models indicated that downregulating the leak K+ channels led to the transition from the
SWS to awake firing patterns. This also suggested the possibility
that a mouse whose leak K+ channel is down-regulated might
have a stronger tendency to generate the awake firing patterns
than the SWS firing patterns and exhibit decreased sleep duration. Therefore, we tested the hypothesis that an impairment ofleak K+ channels would decrease the sleep duration. In themouse genomes, 14 genes [K2p1.1 (Kcnk1), K2p2.1 (Kcnk2),
K2p3.1 (Kcnk3), K2p4.1 (Kcnk4), K2p5.1 (Kcnk5), K2p6.1 (Kcnk6),
K2p7.1 (Kcnk7), K2p9.1 (Kcnk9), K2p10.1 (Kcnk10), K2p12.1(Kcnk12), K2p13.1 (Kcnk13), K2p15.1 (Kcnk15), K2p16.1(Kcnk16), and K2p18.1 (Kcnk18)] are categorized into the K2pchannel family, which is composed of the two-pore-domain potassium channels and is also known as the leak K+ channels
according to their sequence similarity.
To knock out all of the K2p channel family members, we used
the “triple-target CRISPR” method, which achieved almost
perfect efficiency in producing whole-body KO mice in a single
generation by using three guide RNAs targeting a single gene (SIAppendix, Fig. S5 A–N and Table S8) (20), and also used thehighly accurate sleep/wake recording system called the Snappy
Sleep Stager (SSS) (20). Sleep/wake phenotyping was performed
for animals with the KO genotype, which was confirmed by
E9464 | www.pnas.org/cgi/doi/10.1073/pnas.1806486115

quantitative PCR (qPCR; SI Appendix, Fig. S5 P–AC and Table
S9). The Kcnk9 and Kcnk16-KO mice exhibited significantly
decreased sleep duration (Fig. 4 A and C). Because the Kcnk9KO mice exhibited the greatest decrease in sleep duration
among the K2p channel family members, i.e., 606.2 ± 22.8 min(mean ± SEM; n = 8), which was 106.3 min (∼2.6 SD) shorter
than that of WT mice (P < 0.001), we focused on these mice. The
Kcnk9-KO mice exhibited shorter sleep durations than WT miceover most times of the day (Fig. 4B). The Kcnk9-KO mice also
exhibited a higher Psw (the transition probability from a sleep
state to an awake one) but a normal Pws (the transition probability from an awake state to a sleep one) compared with WT
mice (Fig. 4C), suggesting that a decrease in sleep-state stabilization, but not an increase in awake-state stabilization, underlies
their short-sleeper phenotype (i.e., shorter episode durations of asleep state). Although Kcnk9-KO mice were reported to have a(not-significant) tendency to have decreased sleep duration (21,22), our noninvasive sleep phenotyping clearly revealed a significant short-sleeper phenotype in the Kcnk9-KO mice. As
further validation, we generated another group of Kcnk9-KO
mice by using an independent CRISPR/Cas9 probe set (set 2;
SI Appendix, Fig. S5 O and AD and Tables S8 and S9), andconfirmed the observed short-sleeper phenotypes (Fig. 4 D–F).
These results indicated that the observed short-sleeper phenotype of the Kcnk9-KO mice was not caused by an off-target effectof the CRISPR, but rather by common genomic defects in theKcnk9 gene. We also note that the Kcnk9, Kcnk10, and Kcnk16KO mice exhibited increased relative amplitude, which suggested
that the variation of sleep duration within the daily sleep/wakecycle would increase in these mice and the KO of these genes
may have an effect on the regulation of daily sleep patterns (SIAppendix, Fig. S6A).
To exclude the possibility that the gene KO affected the respiration of the mice, causing the SSS system to miscalculate their
sleep duration, we performed EEG/electromyography (EMG)
recordings of the Kcnk9-KO mice. The observed short-sleep
phenotype of Kcnk9-KO mice was confirmed by this EEG/EMG
recording, which revealed that the Kcnk9-KO mice exhibited significant decreases in NREM sleep (i.e., SWS) duration and in total
sleep duration (Fig. 4 G and H). The differences in the transition
probabilities between Kcnk9-KO mice and WT mice were notsignificant (SI Appendix, Fig. S6B). Of particular note, the significant difference of Psw observed in the SSS recording was notconfirmed in the EEG/EMG recording (SI Appendix, Fig. S6B).
This might be because the ages of mice used in each measurement
were different (i.e., 8-wk-old mice were used in the SSS system and
10–14-wk-old mice were used in the EEG/EMG recording) or
because the surgery required for the EEG/EMG recording
might have a slight effect on the detailed structure of sleep/
wake transition. Further study will be needed to clarify whichtransition between NREM, REM, and awake state causes thedecreased NREM sleep duration in Kcnk9-KO mice. Notably,
the Kcnk9-KO mice exhibited a tendency of decreased delta
power, suggesting that Kcnk9 may play a role in the regulationof sleep homeostasis (SI Appendix, Fig. S6 C and D). However,as the difference of delta power between Kcnk9-KO mice andWT mice is not significant, further studies will be needed tovalidate the role of Kcnk9 in the regulation of sleep homeostasis. These results validated the hypothesis that leak K+
channels regulate sleep duration in mammals.
Discussion
Putative Mechanism for Generating the SWS Firing Pattern. In thisstudy, we developed the SAN model, which recapitulates the
SWS firing pattern using five channels and a pump. The SAN
model enabled us to conduct detailed mathematical analyses(i.e., current and phase plane analyses) against the SWS firing pattern, which revealed that the leak K+ channels andYoshida et al.

Parameter sets
20

G
AB
A

KS

g
A
g
AR
g

P

gA
K
g

AM

aL

g
N
a
g
N
a
gP

KL

N

g

ga

M

KC

g DA

N

g

0

NEUROSCIENCE

Voltage-gatedCa2+ -channel KO Model
100
80
60
40
20

G
AB
A

KS

g
A
g
AR
g

P

gA
K
g

AM

g
N
a
g
N
a
gP

aL
N

a

g
KL
g

M

g

0

Voltage-gatedCa2+ -channel KO Model
100
80
60
40
20

G
AB
A

g
K
g
KS
g
A
g
A
gR

PA
AM

N

aL

g
N
a
gN
a
gP

ga
KL
g

KC

A

0
D

Current proportion (%)

Current proportion (%)

40

NOutward current in the down state (%)

60

M

G
AB
A

g
A
g
AR
g

KS

g
K
g

PA
AM

N

aL

g
N
a
g
N
a
gP

0

Outward current in the up state (%)

G
AB
A

KS

g
A
g
A
gR

g
K
g

PA
AM

N

aL

g
N
a
g
N
aP
g

a

g
KL
g

a
C

KC

g
20

80

ΤC
a
g

G
AB
A

g
K
g
KS
g
A
g
A
gR

PA
AM

N

aL

g
N
ga
N
g aP

ga
KL
g

KC

ΤC A
a
g

D
M
N

g
C
ga

0

40

Voltage-gatedCa2+ -channel KO Model

N

20

60

6

100

g

40

80

Bifurcation from SWS to Awake (%)

60

100

G
AB
A

80

Bifurcation proportion
NMDAR KO Model

P

100

0

gA
K
g
KS
g
A
g
A
gR

Original Model

20

aL

E

40

KL

G
AB
A

KS

g
A
g
AR
g

P

gA
K
g

AM

aL

g
N
a
g
N
a
gP

N

g
KL
g

a
KC

N

M

D

gA

0

60

g

20

80

ga
KL
g

40

100

ga

60

Outward current in the up state
NMDAR KO Model

KC

80

0

KC

100

20

g

G
AB
A

KS

g
A
g
AR
g

P

gA
K
g

AM

aL
N

Original Modelg
Cga

Outward current in the up state (%)

D

g
N
a
g
N
a
gP

ga
KL
g

M

KC

g DA

N

g
C
ga

0

40

a

20

60

C

40

0
6
0
Bursting cycles

0

80

g

60

0

100

g

80

100

Outward current in the down state
NMDAR KO Model

g
C
a
ΤC
a
g

100

Outward current in the down state (%)

Original Model

6

Outward current in the up state (%)

C
Outward current in the down state (%)-3
60
Bursting cycles

0

Bifurcation from SWS to Awake (%)

-3

0

N

2+

Standardized Ca

0

100

3

g
N
a
gN
a
gP

Parameter sets
Standardized
MembranePotential

3

Bifurcation from SWS to Awake (%)

0

2

KC

-2

Original Model
Ca2+-dependent K+-current proportion
Leak K+-current proportionn = 1,091
n = 1,091
100
0
100

Parameter sets

2

Parameter sets

-2

B

Intracellular Calciumn = 1,091
-1
0
1

AM

Membrane Potential
n = 1,091
-1
0
1

g DA

Original Model

A

Fig. 3. Leak K+ channels play a role in generating the SWS firing pattern based on the AN model. (A) Standardized membrane potential and intracellular
Ca2+ concentration of 1,091 parameter sets obtained by the parameter searches in the original AN model. The line and shadow (Lower) indicate the mean andSD, respectively. (B) Normalized proportion of the current of leak K+ channels or Ca2+-dependent K+ channels of 1,091 parameter sets obtained by theparameter searches in the original AN model. (C and D) Proportion of each channel’s outward current in the down state (C) or the up state (D) under theoriginal AN model, the NMDAR KO model, or the voltage-gated Ca2+-channel KO model. The leak K+ channel data are in red. (E) Summary of the bifurcation
analysis under the original AN model, the NMDAR KO model, and the voltage-gated Ca2+ channel KO model. The leak K+ channel data are in red. Error bars
indicate SD.

Yoshida et al.

PNAS | vol. 115 | no. 40 | E9465

A

B

C

I

D

E

F

G

H

Fig. 4. Impairment of the leak K+ channels decreases sleep duration in mice. (A) Sleep duration (per hour) over 6 d for mice knocked out for the 14 leak K+
channels. Each row indicates the data from one KO mouse. (B) Sleep duration (per hour) over 24 h, averaged over 6 d in mice knocked out for leak K+
channels. Red lines indicate mean sleep duration at each time of day for each strain. Gray, WT (n = 101). Shaded area, SEM for each time point. (C) Distributions of sleep/wake parameters of the mice knocked out for leak K+ channels. The definition of each parameter (Psw and Pws) is described in SI Appendix,
SI Materials and Methods. Red dots and lines indicate the mean and SEM, respectively. Black dashed line and gray shading indicate the mean and 1 SD range
from the WT data. (D–F) Sleep-duration phenotype for another group of Kcnk9-KO mice (set 2). Sleep duration (per hour) over 6 d (D), sleep duration (perhour) over 24 h averaged over 6 d (E), and distributions of sleep/wake parameters (F). (G and H) Sleep-duration phenotype based on the EEG/EMG for WT (n =7) and Kcnk9-KO mice (set 1; n = 7). Sleep duration (per hour) over 6 d (G). NREM sleep duration, REM sleep duration, and sleep duration for 1 d are shown ina bar graph (H). (I) Schematic of the mechanisms generating the SWS firing pattern. Toward the end of the down state, the Ca2+ concentration is low and theCa2+-dependent K+ channels are inactivated, which initiates the up state. After initiation of the up state, Ca2+ enters through voltage-gated Ca2+ channels
and NMDA receptors, raising the intracellular Ca2+ concentration. Toward the end of the up state, the high Ca2+ concentration activates the Ca2+-dependent
K+ channels, which causes the transition from the up to the down state in cooperation with the leak K+ channels. After initiation of the down state, theintracellular Ca2+ concentration decreases as a result of low Ca2+ entry. When the conductance of leak K+ channels is down-regulated, the up state does nottransition to the down state because the leak K+ channels regulate the threshold for the transition from the up to the down state. WT, C57BL/6N male mice.
Error bars indicate SEM (*P < 0.05, **P < 0.01, and ***P < 0.001; SI Appendix, SI Materials and Methods).

E9466 | www.pnas.org/cgi/doi/10.1073/pnas.1806486115

Yoshida et al.

The SAN Model as a Tool for Elucidating the Mechanisms ofGenerating the SWS Firing Pattern and Regulating Sleep Duration.

In this study, the SAN model predicted that the leak K+ channels
are involved in generating the SWS firing pattern, which led tothe hypothesis that a mouse whose leak K+ channel is downregulated might have a stronger tendency to generate theawake firing patterns than the SWS firing patterns and exhibit
decreased sleep duration. This hypothesis was validated by acomprehensive sleep phenotyping of mice knocked out for leakK+ channel family members such as Kcnk9. Together with theresults of a previous study (13), five components—the leak K+
channel (Kcnk9), Ca2+-dependent K+ channels (Kcnn2 andKcnn3), voltage-gated Ca2+ channels (Cacna1g and Cacna1h),
NMDA receptors (probably Nr1/Nr2b), and Ca2+-exchanger/
pump (Atp2b3)—might be involved in generating the SWS firing
pattern and in the regulation of sleep duration. These results
suggested that generating the SWS firing patterns and regulating
sleep duration could share common molecular mechanisms.
Yoshida et al.

Therefore, elucidating the mechanism of generating the SWS
firing pattern with mathematical models might provide a clue
about the mechanism of regulating sleep duration.
Mathematical analyses using the phase plane method indicated that the intersection of the V nullcline and the nK nullcline at approximately −65 mV regulates the transition from theup to the down state (Fig. 2G and SI Appendix, Fig. S3E). Because the current of each channel around −65 mV determines
the position of this intersection, these currents might be essential
for the transition from the up to the down state. Thus, the balance between the outward currents through leak K+ channels
and Ca2+-dependent K+ channels and the inward currents
through leak Na+ channels could regulate the transition between
the SWS and awake firing patterns. Interestingly, these channels
are also thought to regulate the resting membrane potential,
suggesting that common mechanisms are involved in controlling
the SWS firing pattern and the resting membrane potential. Thisis consistent with the previous report that the resting membranepotential might affect the SWS firing pattern (10). In this sense,
other components that regulate the resting membrane potential(e.g., Na+-K+ ATPase) are candidates for molecular components
generating the SWS firing pattern and possibly regulating sleep
duration. Indeed, mice carrying an inactivating mutation in theneuron-specific Na+-K+ ATPase α3 subunit exhibit decreased
sleep duration (27). Thus, other components that regulate theresting membrane potential are possible molecular candidates
for regulating sleep duration. To assess the association of thesecomponents with the generation of SWS firing patterns, it might
be possible to add each of them to the SAN model and to analyze
its role mathematically in detail. Taken together, the SAN model
predicts that common mechanisms are involved in regulating the
SWS firing pattern and the resting membrane potential andcould be a versatile model for investigating the role of other ion
channels in generating the SWS firing pattern and regulating
sleep duration.
The Role of Leak K+ Channels in Generating the SWS Firing Pattern
and Regulating Sleep Duration. Although mathematical approaches

predicted that the leak K+ channels play a role in generating theup and down states of the SWS firing pattern, its underlying
mechanism is still a mystery. It would be a promising approach toclarify the mechanism of generating the up and down states incortical neurons electrophysiologically. In fact, the cortex has
numerous types of neurons, and it is still difficult to measure andperturb an activity of a specific neuron, such as the SWS firing
pattern, in vivo. If such technical limitations can be overcome, themechanism of generating the SWS firing pattern as well as its
relationship with the regulation of sleep duration will be uncovered. Given the current difficulty of in vivo validation, the next
step will be conducting an electrophysiological experiment toevaluate the role of leak K+ channels in a cortical slice. It would
be a good surrogate to measure the membrane potential of a
cortical slice under Kcnk9-KO or Kcnk9 knock-down conditions.
One interesting question rising from this study is why Kcnk9KO mice exhibited a severe decrease in sleep duration whereas
the other KO mice of the K2p channel family members did not.
One possible reason for this result is a difference in neuronal
expression between the K2p channel family members. Although
Kcnk9 has expression in various brain regions including cerebral
cortex, other members have different expression patterns in mice(28), which implies that some members might have weak expression in cerebral cortex. In addition, cell-type specificity of
Kcnk9 expression could be also related to the sleep phenotype of
its KO mice. Investigating where and in what type of neuron
Kcnk9 has strong expression could be an interesting future direction to validate the role of Kcnk9 in the regulation of sleep
duration. Another possible reason is the compensation of therole of a K2p channel family member by other members. The KO
PNAS | vol. 115 | no. 40 | E9467

NEUROSCIENCE

Ca2+-dependent K+ channels cooperatively induce the downstate, and hence the SWS firing pattern. According to the ANmodel (13, 14, 23, 24), the putative intracellular mechanism forgenerating the SWS firing pattern is as follows (Fig. 4I). (i) During
the up state, Ca2+ enters mainly through NMDA receptors andvoltage-gated Ca2+ channels to increase the intracellular Ca2+
concentration. (ii) The transition from the up to the down state
occurs when the intracellular Ca2+ concentration reaches a certain
threshold, which activates Ca2+-dependent K+ channels. (iii) In
the down state, Ca2+ exits through Ca2+-pump/exchangers to decrease the intracellular Ca2+ concentration. (iv) The transition
from the down to the up state occurs when the intracellular Ca2+
concentration decreases, which inactivates the Ca2+-dependent
K+ channels. The SAN model further indicated that the conductance of leak K+ channels alters the threshold for the transition
from the up to the down state mediated by the Ca2+-dependent
K+ channels. When the conductance of the leak K+ channels is
relatively high, a relatively lower conductance of the Ca2+dependent K+ channels is sufficient to induce the down state.
Notably, this putative mechanism for generating the SWS
firing pattern is consistent with the reported extracellular K+ bias
between sleep and wakefulness: the extracellular K+ concentration is relatively low during sleep and high during wakefulness invivo (25). Because an extracellular K+ bias will alter the tendencyof the K+ flux, the relatively low extracellular K+ concentration
during sleep could drive K+ efflux through the leak K+ channels
and Ca2+-dependent K+ channels, thereby generating the downstate. Consistent with this scenario, the extracellular Ca2+ concentration is relatively high during sleep (25), which could drive
Ca2+ influx through the voltage-dependent Ca2+ channels and
NMDA receptors, thereby generating the down state by activating the Ca2+-dependent K+ channels during sleep. Indeed, byimplementing ion concentration dependency in the AN model, a
recent computational study indicated that the extracellular ion
concentration and the conductance of the Ca2+-dependent K+
channels could coordinately regulate the transition between the
SWS and awake firing patterns (26). The SAN model also predicted that persistent Na+ channels are involved in initiating theup state (SI Appendix, Fig. S2C), which is a proposed mechanism
for rhythmic bursting patterns in the pre-Bötzinger complex or
hippocampal pyramidal neurons (15, 16, 18, 19). These results
show that the common mechanisms may underlie the generationof the SWS firing pattern and other rhythmic bursting patterns.
Collectively, we hypothesize that the Ca2+-dependent hyperpolarization pathway and leak K+ channels are involved in generating the SWS firing pattern. A remaining challenge will be to
combine this putative intracellular mechanism with the accumulated knowledge of intercellular mechanisms in sleep/wakecycle regulation.

effect of a K2p channel family member could be reduced if theother members have stronger expression than WT mice andcompensate the gene’s role. Further study such as investigating
the expression of other K2p channel family members in KO miceor measuring sleep duration in double-KO mice of K2p channelfamily members will be needed to examine this possibility.
Materials and Methods
Details of the study materials and methods are supplied in the SI Appendix, SI
Materials and Methods.
SAN Model. The SAN model was constructed by simplifying the AN model. The
full system of equations is presented in SI Appendix, SI Materials andMethods. Parameter searches for the SWS firing patterns were conducted,
and the SAN model was determined by comparing the hit rates of various
models, which are shown in SI Appendix, Table S7, and conducting thebifurcation analyses.

NMDA receptors model, and the KO of voltage-gated Ca2+ channels model.
The bifurcation analyses were conducted with all of the parameter sets inthe original AN model, the KO of NMDA receptors model, and the KO ofvoltage-gated Ca2+ channels model, in all of which the leak channels weredivided into the leak K+ channels and the leak Na+ channels.
Animals and Sleep Phenotyping. All experimental procedures and housing
conditions involving animals and their care were approved by the institutional animal care and use committee. Sleep phenotyping was performed noninvasively by SSS. KO mice were generated and subjected to sleep
phenotyping at the University of Tokyo, and C57BL/6N mice (n = 101) wereused as controls (Fig. 4 B, C, E, and F and SI Appendix, Fig. S6A). The EEG/
EMG recording for Kcnk9-KO mice was performed at the University ofTokyo. CRISPR-KO animals were generated by one-cell embryo microinjection of synthesized Cas9 mRNA and gRNAs (mix of three targeting sequences; template sequences are shown in SI Appendix, Table S8) into
C57BL/6N fertilized eggs. KO genotypes were confirmed by qPCR.

Further Analysis of the AN Model. The currents in the down and up states were
calculated with all of the parameter sets in the original AN model, the KO of

ACKNOWLEDGMENTS. We thank all the laboratory members at RIKEN
Center for Biosystems Dynamics Research and the University of Tokyo,
particularly S. Tomita, for their kind help in the genotyping of animals andsupporting experiments; K. Kuwana and M. Shiokawa, for their kind help insynthesizing gRNAs; and K. Ode, F. Tatsuki, and S. Sugai for helpful discussion.
This work was supported by Japan Agency for Medical Research and Development (AMED)-Core Research for Evolutional Science and Technology (CREST)[AMED/Ministry of Education, Culture, Sports, Science and Technology (MEXT)
Grant JP17GM0610006 to H.R.U.], Brain Mapping by Integrated Neurotechnologies for Disease Studies (Brain/MINDS) (AMED/MEXT Grant JP17DM0207049 to
H.R.U.), Basic Science and Platform Technology Program for Innovative Biological
Medicine (AMED/MEXT Grant JP17AM0301025 to H.R.U.), World Premier International Research Center Initiative (MEXT; H.R.U.), Grant-in-Aid for Research
Activity Start-up [Japan Society for the Promotion of Science (JSPS) KAKENHI
Grant 17H06623 to S.S.], Grant-in-Aid for Scientific Research (S) (JSPS KAKENHI
Grant 25221004 to H.R.U.), and Grant-in-Aid for Research from the Son
Masayoshi Foundation (to K.Y.).

1. Siegel JM (2005) Clues to the functions of mammalian sleep. Nature 437:1264–1271.
2. Krueger JM, et al. (2008) Sleep as a fundamental property of neuronal assemblies. Nat
Rev Neurosci 9:910–919.
3. Saper CB, Scammell TE, Lu J (2005) Hypothalamic regulation of sleep and circadianrhythms. Nature 437:1257–1263.
4. Caton R (1875) Electrical currents of the brain. Chicago J Nerv Mental Dis 2:610.
5. Berger H (1929) U€ ber das Elektrenkephalogramm des Menschen. Arch Psychiatr
Nervenkr 87:527–570.
6. Steriade M, Timofeev I, Grenier F (2001) Natural waking and sleep states: A view frominside neocortical neurons. J Neurophysiol 85:1969–1985.
7. Steriade M (2003) The corticothalamic system in sleep. Front Biosci J Virtual Library 8:
d878–d899.
8. Bazhenov M, Timofeev I, Steriade M, Sejnowski TJ (2002) Model of thalamocorticalslow-wave sleep oscillations and transitions to activated States. J Neurosci 22:
8691–8704.
9. Chen JY, Chauvette S, Skorheim S, Timofeev I, Bazhenov M (2012) Interneuronmediated inhibition synchronizes neuronal activity during slow oscillation. J Physiol
590:3987–4010.
10. Compte A, Sanchez-Vives MV, McCormick DA, Wang XJ (2003) Cellular and network
mechanisms of slow oscillatory activity (<1 Hz) and wave propagations in a cortical
network model. J Neurophysiol 89:2707–2725.
11. Hill S, Tononi G (2005) Modeling sleep and wakefulness in the thalamocortical system.
J Neurophysiol 93:1671–1698.
12. Sanchez-Vives MV, et al. (2010) Inhibitory modulation of cortical up states.
J Neurophysiol 104:1314–1324.
13. Tatsuki F, et al. (2016) Involvement of Ca(2+)-dependent hyperpolarization in sleep
duration in mammals. Neuron 90:70–85.
14. Ode KL, Katsumata T, Tone D, Ueda HR (2017) Fast and slow Ca2+-dependent hyperpolarization mechanisms connect membrane potential and sleep homeostasis.
Curr Opin Neurobiol 44:212–221.
15. Butera RJ, Jr, Rinzel J, Smith JC (1999) Models of respiratory rhythm generation in the
pre-Bötzinger complex. I. Bursting pacemaker neurons. J Neurophysiol 82:382–397.

16. Golomb D, Yue C, Yaari Y (2006) Contribution of persistent Na+ current and M-type
K+ current to somatic bursting in CA1 pyramidal cells: Combined experimental andmodeling study. J Neurophysiol 96:1912–1926.
17. Izhikevich EM (2007) Dynamical Systems in Neuroscience: The Geometry of Excitability
and Bursting (MIT Press, Cambridge, MA).
18. Del Negro CA, Koshiya N, Butera RJ, Jr, Smith JC (2002) Persistent sodium current,
membrane properties and bursting behavior of pre-Bötzinger complex inspiratoryneurons in vitro. J Neurophysiol 88:2242–2250.
19. Su H, Alroy G, Kirson ED, Yaari Y (2001) Extracellular calcium modulates persistentsodium current-dependent burst-firing in hippocampal pyramidal neurons. J Neurosci
21:4173–4182.
20. Sunagawa GA, et al. (2016) Mammalian reverse genetics without crossing reveals
Nr3a as a short-sleeper gene. Cell Rep 14:662–677.
21. Pang DS, et al. (2009) An unexpected role for TASK-3 potassium channels in network
oscillations with implications for sleep mechanisms and anesthetic action. Proc Natl
Acad Sci USA 106:17546–17551.
22. Gotter AL, et al. (2011) TASK-3 as a potential antidepressant target. Brain Res 1416:
69–79.
23. Shi S, Ueda HR (2018) Ca2+-dependent hyperpolarization pathways in sleep homeostasis and mental disorders. Bioessays 40.
24. Tatsuki F, Ode KL, Ueda HR (2017) Ca2+-dependent hyperpolarization hypothesis for
mammalian sleep. Neurosci Res 118:48–55.
25. Ding F, et al. (2016) Changes in the composition of brain interstitial ions control thesleep-wake cycle. Science 352:550–555.
26. Rasmussen R, Jensen MH, Heltberg ML (2017) Chaotic dynamics mediate brain statetransitions, driven by changes in extracellular ion concentrations. Cell Syst 5:591–603.e4.
27. Kirshenbaum GS, et al. (2011) Mania-like behavior induced by genetic dysfunction of
the neuron-specific Na+,K+-ATPase α3 sodium pump. Proc Natl Acad Sci USA 108:
18144–18149.
28. Aller MI, Wisden W (2008) Changes in expression of some two-pore domain potassium channel genes (KCNK) in selected brain regions of developing mice.
Neuroscience 151:1154–1172.

Analysis of the SAN Model. We divided all of the parameter sets with the SWS
firing patterns into two clusters according to the conductance of the voltagegated K+ channels and voltage-gated Ca2+ channels. The representative parameter set was determined by using kernel density estimation in each cluster.
First, the proportion of the current of each channel was analyzed. Second, innormal bifurcation and 2D analyses, each channel’s conductance was gradually changed from the representative parameter sets. Third, the trajectory ofthese differential equations was analyzed by the phase plane. The changes inthe trajectory and the nullclines caused by down-regulating the conductanceof leak K+ channels or Ca2+-dependent K+ channels were analyzed. Fourth,bifurcation analysis of all of the parameter sets was conducted, in which thetransitions from the SWS to awake firing patterns were counted.

E9468 | www.pnas.org/cgi/doi/10.1073/pnas.1806486115

Yoshida et al.

Supplementary Information forLeak potassium channels regulate sleep duration
Kensuke Yoshida1,4, Shoi Shi1,2,3,4, Maki Ukai-Tadenuma2,3, Hiroshi Fujishima2, Rei-ichiro
Ohno1, and Hiroki R. Ueda1,2,3*
1

Department of Systems Pharmacology, Graduate School of Medicine, The University ofTokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
2
Laboratory for Synthetic Biology, RIKEN Center for Biosystems Dynamics Research, 1-3
Yamadaoka, Suita, Osaka 565-0871, Japan
3
International Research Center for Neurointelligence (WPI-IRCN), UTIAS, The Universityof Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
4
These authors contributed equally to this work.
Correspondence: Hiroki R. Ueda
Email: uedah-tky@umin.ac.jp
This PDF file includes:
Supplementary text
Figs. S1 to S6Tables S1 to S9References for SI reference citations

www.pnas.org/cgi/doi/10.1073/pnas.1806486115

Sl Materials and Methods
Software for computer simulation
The probability density analysis for determining the representative parameter set was
conducted by using MATLAB R2016b. All other numeric calculations were conducted byusing Mathematica software version 10.0 (Wolfram Research).

Original model and KO models of channels
To analyze the SWS firing pattern, a computational model was constructed based on theprevious study (1). The differential equations for this model are described below.
𝐶𝐴𝑑𝑉= −𝐴(𝐼𝐿 (𝑉) + 𝐼𝑁𝑎 (𝑉, ℎ𝑛𝑎 ) + 𝐼𝐾 (𝑉, 𝑛𝑘 ) + 𝐼𝐴 (𝑉, ℎ𝐴 ) + 𝐼𝐾𝑆 (𝑉, 𝑚𝐾𝑆 ) + 𝐼𝐶𝑎 (𝑉)𝑑𝑡+ 𝐼𝐾𝐶𝑎 (𝑉, [𝐶𝑎2+ ]) + 𝐼𝑁𝑎𝑃 (𝑉) + 𝐼𝐴𝑅 (𝑉)) − 𝐼𝑁𝑀𝐷𝐴 (𝑉, 𝑠𝑁𝑀𝐷𝐴 , 𝑥𝑁𝑀𝐷𝐴 )− 𝐼𝐴𝑀𝑃𝐴 (𝑉, 𝑠𝐴𝑀𝑃𝐴 ) − 𝐼𝐺𝐴𝐵𝐴 (𝑉, 𝑠𝐺𝐴𝐵𝐴 )𝑑ℎ𝑛𝑎= 4(𝛼ℎ (𝑉)(1 − ℎ𝑁𝑎 ) − 𝛽ℎ (𝑉)ℎ𝑁𝑎 )𝑑𝑡𝑑𝑛𝑘= 4(𝛼𝑛 (𝑉)(1 − 𝑛𝐾 ) − 𝛽𝑛 (𝑉)𝑛𝐾 )𝑑𝑡𝑑ℎ𝐴= (ℎ𝐴∞ (𝑉) − ℎ𝐴 )/𝜏ℎ𝐴
𝑑𝑡
𝑑𝑚𝐾𝑆= (𝑚𝐾𝑆∞ (𝑉) − 𝑚𝐾𝑆 )/𝜏𝑚𝐾𝑆 (𝑉)𝑑𝑡𝑑𝑠𝐴𝑀𝑃𝐴
𝑠𝐴𝑀𝑃𝐴= 3.48𝑓(𝑉) −𝑑𝑡
𝜏𝐴𝑀𝑃𝐴
𝑑𝑠𝑁𝑀𝐷𝐴𝑠𝑁𝑀𝐷𝐴= 0.5𝑥𝑁𝑀𝐷𝐴 (1 − 𝑠𝑁𝑀𝐷𝐴 ) −𝑑𝑡
𝜏𝑠𝑁𝑀𝐷𝐴
𝑑𝑥𝑁𝑀𝐷𝐴
𝑥𝑁𝑀𝐷𝐴= 3.48𝑓(𝑉) −𝑑𝑡
𝜏𝑥𝑁𝑀𝐷𝐴

𝑑𝑠𝐺𝐴𝐵𝐴𝑠𝐺𝐴𝐵𝐴= 𝑓(𝑉) −𝑑𝑡
𝜏𝐺𝐴𝐵𝐴[𝐶𝑎2+ ]𝑑[𝐶𝑎2+ ]= −𝛼𝐶𝑎 (𝐴𝐼𝐶𝑎 (𝑉) + 𝐼𝑁𝑀𝐷𝐴 (𝑉, 𝑠𝑁𝑀𝐷𝐴 , 𝑥𝑁𝑀𝐷𝐴 )) −𝑑𝑡
𝜏𝐶𝑎
, where 𝐶 is the membrane capacitance, 𝐴 is the area of a single neuron, 𝑉 is themembrane potential, [𝐶𝑎2+ ] is the intracellular calcium concentration, and 𝐼𝑋 (𝑋: each ion
channel) denotes the current of each ion channel. Each function is listed below.
𝐼𝐿 (𝑉) = 𝑔𝐿 (𝑉 − 𝑉𝐿 )3
𝐼𝑁𝑎 (𝑉, ℎ𝑛𝑎 ) = 𝑔𝑁𝑎 𝑚𝑁𝑎∞
(𝑉)ℎ𝑁𝑎 (𝑉 − 𝑉𝑁𝑎 )𝑚𝑁𝑎∞ (𝑉) = 𝛼𝑚 (𝑉)/(𝛼𝑚 (𝑉) + 𝛽𝑚 (𝑉))𝛼𝑚 (𝑉) = 0.1(𝑉 + 33)/[1 − exp⁡(−(𝑉 + 33)/10)]
𝛽𝑚 (𝑉) = 4exp⁡(−(𝑉 + 53.7)/12)𝛼ℎ (𝑉) = 0.07exp⁡(−(𝑉 + 50)/10)𝛽ℎ (𝑉) = 1/[1 + exp⁡(−(𝑉 + 20)/10)]
𝐼𝐾 (𝑉, 𝑛𝑘 ) = 𝑔𝐾 𝑛𝐾4 (𝑉 − 𝑉𝐾 )𝛼𝑛 (𝑉) = 0.01(𝑉 + 34)/[1 − exp(−(𝑉 + 34)/10)]
𝛽𝑛 (𝑉) = 0.125exp⁡(−(𝑉 + 44)/25)
3
𝐼𝐴 (𝑉, ℎ𝐴 ) = 𝑔𝐴 𝑚𝐴∞(𝑉)ℎ𝐴 (𝑉 − 𝑉𝐾 )𝑚𝐴∞ (𝑉) = 1/[1 + exp⁡(−(𝑉 + 50)/20]ℎ𝐴∞ (𝑉) = 1/[1 + exp⁡((𝑉 + 80)/6)]
𝐼𝐾𝑆 (𝑉, 𝑚𝐾𝑆 ) = 𝑔𝐾𝑆 𝑚𝐾𝑆 (𝑉 − 𝑉𝐾 )𝑚𝐾𝑆∞ (𝑉) = 1/[1 + exp⁡(−(𝑉 + 34)/6.5)]
𝜏𝑚𝐾𝑆 (𝑉) = 8/[exp(−(𝑉 + 55)/30) + exp⁡((𝑉 + 55)/30)]
2(𝑉)(𝑉 − 𝑉𝐶𝑎 )𝐼𝐶𝑎 (𝑉) = 𝑔𝐶𝑎 𝑚𝐶𝑎∞

𝑚𝐶𝑎∞ (𝑉) = 1/[1 + exp(−(𝑉 + 20)/9)]

𝐼𝐾𝐶𝑎 (𝑉, [𝐶𝑎2+ ]) = 𝑔𝐾𝐶𝑎 𝑚𝐾𝐶𝑎∞ ([𝐶𝑎2+ ])(𝑉 − 𝑉𝐾 )𝑚𝐾𝐶𝑎∞ ([𝐶𝑎2+ ]) = 1/[1 + (𝐾𝐷 /[𝐶𝑎2+ ])3.5 ]3𝐼𝑁𝑎𝑃 (𝑉) = 𝑔𝑁𝑎𝑃 𝑚𝑁𝑎𝑃∞(𝑉)(𝑉 − 𝑉𝑁𝑎 )𝑚𝑁𝑎𝑃∞ (𝑉) = 1/[1 + exp⁡(−(𝑉 + 55.7)/7.7)]
𝐼𝐴𝑅 (𝑉) = 𝑔𝐴𝑅 ℎ𝐴𝑅∞ (𝑉)(𝑉 − 𝑉𝐾 )ℎ𝐴𝑅∞ (𝑉) = 1/[1 + exp⁡((𝑉 + 75)/4)]𝑓(𝑉) = 1/[1 + exp⁡(−(𝑉 − 20)/2)]
𝐼𝐴𝑀𝑃𝐴 (𝑉, 𝑠𝐴𝑀𝑃𝐴 ) = 𝑔𝐴𝑀𝑃𝐴 𝑠𝐴𝑀𝑃𝐴 (𝑉 − 𝑉𝐴𝑀𝑃𝐴 )𝐼𝑁𝑀𝐷𝐴 (𝑉, 𝑠𝑁𝑀𝐷𝐴 , 𝑥𝑁𝑀𝐷𝐴 ) = 𝑔𝑁𝑀𝐷𝐴 𝑠𝑁𝑀𝐷𝐴 (𝑉 − 𝑉𝑁𝑀𝐷𝐴 )𝐼𝐺𝐴𝐵𝐴 (𝑉, 𝑠𝐺𝐴𝐵𝐴 ) = 𝑔𝐺𝐴𝐵𝐴 𝑠𝐺𝐴𝐵𝐴 (𝑉 − 𝑉𝐺𝐴𝐵𝐴 )Please note that intrinsic (non-synaptic) ion currents (e.g. 𝐼𝐶𝑎 ) should be multipliedby 10 to adjust its unit to nanoampere (nA) if the numerical values listed in Tables S1 and S2
are directly used in the numerical simulation. The constant values are listed in Tables S1 andS2 and the initial values of these differential equations are listed in Table S3. All differential
equations were solved in t = [0, 10] except for t = [0, 30] and t = [0, 90] in Fig. 3 C
and D and Fig. S4 D and E. The plotted time ranges are described in Table S4.
In the KO models of channels, the conductance of the target channels was set to zero.
For example, 𝑔𝑁𝑎 , 𝑔𝐴 , 𝑔𝐾𝑆 , ⁡𝑔𝐴𝑅 , 𝑔𝐴𝑀𝑃𝐴 , 𝑔𝑁𝑀𝐷𝐴⁡ , and 𝑔𝐺𝐴𝐵𝐴 were set to zero in the
SAN model.
In some analyses, the conductance of leak channels 𝑔𝐿 was divided into the
conductance of leak K+ channels (𝑔𝐾𝐿 ) and the conductance of leak Na+ channels (𝑔𝑁𝑎𝐿 ) as
follows:
𝑔𝐾𝐿 = 𝑔𝐿 ×

𝑉𝐿 − 𝑉𝑁𝑎𝐿
𝑉𝐾 − 𝑉𝑁𝑎𝐿

𝑔𝑁𝑎𝐿 = 𝑔𝐿 ×

𝑉𝐿 − 𝑉𝐾𝑉𝑁𝑎𝐿 − 𝑉𝐾

These definitions satisfy the following conditions:
𝐼𝐿 (𝑉) = 𝑔𝐿 (𝑉 − 𝑉𝐿 ) = 𝑔𝐾𝐿 (𝑉 − 𝑉𝐾 ) + 𝑔𝑁𝑎𝐿 (𝑉 − 𝑉𝑁𝑎𝐿 )𝑔𝐿 = 𝑔𝐾𝐿 + 𝑔𝑁𝑎𝐿

SAN model
The SAN model is described on the condition that 𝑔𝑁𝑎 , 𝑔𝐴 , 𝑔𝐾𝑆 , ⁡𝑔𝐴𝑅 , 𝑔𝐴𝑀𝑃𝐴 , 𝑔𝑁𝑀𝐷𝐴⁡ ,and 𝑔𝐺𝐴𝐵𝐴 are set to zero in the original model. Thus, the differential equations of the SAN
model are as follows.
𝐶𝑑𝑉= −(𝐼𝐿 (𝑉) + 𝐼𝐾 (𝑉, 𝑛𝑘 ) + 𝐼𝐶𝑎 (𝑉) + 𝐼𝐾𝐶𝑎 (𝑉, [𝐶𝑎2+ ]) + 𝐼𝑁𝑎𝑃 (𝑉))𝑑𝑡𝑑𝑛𝑘
= 4(𝛼𝑛 (𝑉)(1 − 𝑛𝐾 ) − 𝛽𝑛 (𝑉)𝑛𝐾 )𝑑𝑡[𝐶𝑎2+ ]𝑑[𝐶𝑎2+ ]= −𝛼𝐶𝑎 𝐴𝐼𝐶𝑎 (𝑉) −𝑑𝑡
𝜏𝐶𝑎𝐼𝐿 (𝑉) = 𝑔𝐿 (𝑉 − 𝑉𝐿 )𝐼𝐾 (𝑉, 𝑛𝑘 ) = 𝑔𝐾 𝑛𝐾4 (𝑉 − 𝑉𝐾 )𝛼𝑛 (𝑉) = 0.01(𝑉 + 34)/[1 − exp(−(𝑉 + 34)/10)]
𝛽𝑛 (𝑉) = 0.125exp⁡(−(𝑉 + 44)/25)
2(𝑉)(𝑉 − 𝑉𝐶𝑎 )𝐼𝐶𝑎 (𝑉) = 𝑔𝐶𝑎 𝑚𝐶𝑎∞

𝑚𝐶𝑎∞ (𝑉) = 1/[1 + exp(−(𝑉 + 20)/9)]
𝐼𝐾𝐶𝑎 (𝑉, [𝐶𝑎2+ ]) = 𝑔𝐾𝐶𝑎 𝑚𝐾𝐶𝑎∞ ([𝐶𝑎2+ ])(𝑉 − 𝑉𝐾 )𝑚𝐾𝐶𝑎∞ ([𝐶𝑎2+ ]) = 1/[1 + (𝐾𝐷 /[𝐶𝑎2+ ])3.5 ]3𝐼𝑁𝑎𝑃 (𝑉) = 𝑔𝑁𝑎𝑃 𝑚𝑁𝑎𝑃∞(𝑉)(𝑉 − 𝑉𝑁𝑎 )𝑚𝑁𝑎𝑃∞ (𝑉) = 1/[1 + exp⁡(−(𝑉 + 55.7)/7.7)]

Parameter search
The method of the parameter searches was based on the previous study (1). A parameter set
includes 13 values (𝑔𝐿 , 𝑔𝑁𝑎 , 𝑔𝐾 , 𝑔𝐴 , 𝑔𝐾𝑆 , 𝑔𝐶𝑎 , 𝑔𝐾𝐶𝑎 , 𝑔𝑁𝑎𝑃⁡ ,⁡𝑔𝐴𝑅 , 𝑔𝐴𝑀𝑃𝐴 , 𝑔𝑁𝑀𝐷𝐴⁡ ,𝑔𝐺𝐴𝐵𝐴⁡ , 𝜏𝐶𝑎 ). All parameter sets were generated so that the logarithm of each parameter has auniform distribution. Let U(𝑎, 𝑏) be a distribution whose probability density function (𝑓(𝑥))is as follows.
1⁡⁡⁡⁡⁡⁡⁡for⁡𝑎 ≤ 𝑥 ≤ 𝑏⁡⁡⁡⁡𝑓(𝑥) = { (log 𝑏 − log 𝑎)𝑥⁡0⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡for⁡𝑥 < 𝑎, 𝑥 > 𝑏
The values of 𝑔𝐿 , 𝑔𝑁𝑎 , 𝑔𝐾 , 𝑔𝐴 , 𝑔𝐾𝑆 , 𝑔𝐶𝑎 , 𝑔𝐾𝐶𝑎 , 𝑔𝑁𝑎𝑃⁡ , and 𝑔𝐴𝑅 (mS/cm2) weregenerated according to U(0.01, 100), 𝑔𝐴𝑀𝑃𝐴 , 𝑔𝑁𝑀𝐷𝐴 , and 𝑔𝐺𝐴𝐵𝐴 (µS) were generated
according to U(0.002, 20), and⁡𝜏𝐶𝑎 (ms) was generated according to U(10, 1000). For eachparameter set, the differential equations were solved by using the NDSolve function in
Mathematica 10.0. The parameter sets that did not yield real solutions under the definedcomputer memory usage (1GB) and calculation time (50 seconds) were excluded. Let 𝑉(𝑡)be a membrane potential obtained by solving the differential equation, and 𝑣𝑖 was defined as𝑉(4.999 + 0.001𝑖) (𝑖 = 1, 2, … , 5001). First, the parameter sets that at least one of the value|𝑣𝑖 | (𝑖 = 1, 2, … , 5001) was bigger than 200 mV were excluded. Then the solutions wereclassified into four categories based on the peak frequency and the number of spikes (Table
S5). The lists of 𝑣𝑖 were analyzed by the discrete Fourier transform after normalization anddetrending. The frequency whose Fourier power is the largest between 0 Hz and 500 Hz was
defined as the peak frequency. The number of spikes was defined as half the number of theelements in the following set S.
S = {(𝑣𝑖 , 𝑣𝑖+1 )|1 ≤ 𝑖 ≤ 5000, (𝑣𝑖 + 20)(𝑣𝑖+1 + 20) < 0}
All the solutions classified into “SWS” were then checked manually and some of them wereclassified into “True SWS”. A hit rate was calculated as a proportion of the parameter sets
classified into “True SWS” in all the randomly-generated parameter sets. We interpreted that
two hit rates were different when one of them was more than 100 times as high as the other.

Normalizing the frequency of membrane potential
Normalizing the frequency of membrane potential was executed by the following procedure.
Let 𝑀, 𝑚 be maximum value and minimum value of 𝑣𝑖 respectively. The set 𝑆𝑖 (𝑖 =0, 1, 2, … , [𝑀 − 𝑚] ([] is a floor function)) was defined as follows.
𝑆𝑖 = {𝑗 ∊ ℤ⁡⁡|1 ≤ 𝑗 ≤ 5000, (𝑣𝑗 − (𝑀 − 𝑖)) ・ (𝑣𝑗+1 − (𝑀 − 𝑖)) < 0}

Let 𝑐𝑖,𝑗 be the 𝑗-th minimum integer of 𝑆𝑖 . Let 𝑑𝑖 be the maximum value of the set{𝑐𝑖,𝑗+1 − 𝑐𝑖,𝑗 |1 ≤ 𝑗 ≤ |𝑆𝑖 | − 1}. The set A was defined as a set of the integer 𝑘 that |𝑆𝑘 |was the largest in all⁡𝑆𝑖 . The integer ℎ was defined as the minimum number of the set{𝑘 ∈ 𝐴|𝑑𝑘 = min 𝑑𝑖 }. I set the 𝑀 − ℎ as the threshold voltage. The set 𝐸 was defined as𝑖∈𝐴follows.
𝐸 = {𝑗|𝑐ℎ,𝑗+1 − 𝑐ℎ,𝑗 ≥ 0.5𝑑ℎ }If |𝐸| ≤ 6, the time range in which the differential equation was solved was spread and thesame procedure was conducted. If |𝐸| ≥ 7, let 𝑠, 𝑓 as the the minimum, 7-th minimum
integer of E, and the standardized membrane potential, the standardized Ca2+ concentration,or the proportion of ion currents of 𝑡 = [5 + 0.001𝑠, 5 + 0.001𝑓] were shown (Figs. 1M
and 3 A and B and Figs. S1 D-H and S4 A and B). The membrane potential and Ca2+
concentration are standardized in 𝑡 = [5 + 0.001𝑠, 5 + 0.001𝑓] to set the mean andstandard deviation zero and one, respectively. In some parameter sets, the value 0.5 in thedefinition of E or the time range in which the differential equation was solved was changed inorder to judge the bursting cycle correctly.

Calculation of the representative parameter sets
The multivariate kernel density was estimated in log scale for each cluster by using the
mvksdensity function in MATLAB. Let 𝜎𝑖 and 𝑛 be the corrected sample standard
deviation of each parameter and the number of parameter sets respectively. In reference to theSilverman’s rule (2), the bandwidth 𝑏𝑖 of each parameter was set as follows:1

1 10
𝑏𝑖 = 𝜎𝑖 { } .
2𝑛The 216 grid points were set by the 21 geometric sequences whose minimum was equal to theminimum parameter and whose maximum was equal to the maximum parameter for eachparameter. The point whose multivariable kernel density was maximum was defined as therepresentative parameter in each cluster.

The proportion of currents
In the Fig. 3 C and D and Fig. S4 D and E, the differential equations were solved in t =[0,10], t = [0,30], or t = [0,90]. Some parameter sets were excluded from this analysis
because their frequency could not be normalized or they did not show the SWS firing patterns
when the differential equations were solved in⁡t = [0,30] or t = [0,90]. Let 𝑊 be the setof all 13 channels, receptors and a pump. The proportion of channel X’s inward current 𝑃𝐸𝑋
and outward currents 𝑃𝐼𝑋 (Figs. 2C and 3B and Figs. S2C and S4 A-C) was defined asfollows:
𝑃𝐸𝑋 = (𝐼𝑋 + |𝐼𝑋 |)/ ∑ (𝐼𝐶 + |𝐼𝐶 |)𝐶∊𝑊𝑃𝐼𝑋 = (𝐼𝑋 − |𝐼𝑋 |)/ ∑ (𝐼𝐶 − |𝐼𝐶 |)𝐶∊𝑊

Let 𝑚1 , 𝑚2 be the minimum and the second minimum element of the set 𝐸. The average
current 𝑆𝐼𝑋 in the down state (Fig. 3C and Fig. S4D) and the average current 𝐵𝐼𝑋 in the up
state (Fig. 3D and Fig. S4E) are respectively defined as follows:1
𝑆𝐼𝑋 =𝑐ℎ,𝑚2 +1 − 𝑐ℎ,𝑚2

1
𝐵𝐼𝑋 =𝑐ℎ,𝑚2 − 𝑐ℎ,𝑚1 +1

𝑐ℎ,𝑚2 +1

𝐼𝑋 (4.999 + 0.001𝑡)∑𝑡=𝑐ℎ,𝑚2 +1
𝑐ℎ,𝑚2

∑

𝐼𝑋 (4.999 + 0.001𝑡)𝑡=𝑐ℎ,𝑚1 +1 +1

The proportion of channel X’s outward current 𝑃𝑆𝐸𝑋 and inward currents 𝑃𝑆𝐼𝑋 in thedown state (Fig. 3C and Fig. S4D), and the proportion of channel X’s outward current 𝑃𝐵𝐸𝑋
and inward currents 𝑃𝐵𝐼𝑋 in the up state (Fig. 3D and Fig. S4E) were defined as follows:
𝑃𝑆𝐸𝑋 = (𝑆𝐼𝑋 + |𝑆𝐼𝑋 |)/ ∑ (𝑆𝐼𝐶 + |𝑆𝐼𝐶 |)𝐶∊𝑊

𝑃𝑆𝐼𝑋 = (𝑆𝐼𝑋 − |𝑆𝐼𝑋 |)/ ∑ (𝑆𝐼𝐶 − |𝑆𝐼𝐶 |)𝐶∊𝑊𝑃𝐵𝐸𝑋 = (𝐵𝐼𝑋 + |𝐵𝐼𝑋 |)/ ∑ (𝐵𝐼𝐶 + |𝐵𝐼𝐶 |)𝐶∊𝑊𝑃𝐵𝐼𝑋 = (𝐵𝐼𝑋 − |𝐵𝐼𝑋 |)/ ∑ (𝐵𝐼𝐶 − |𝐵𝐼𝐶 |)𝐶∊𝑊

Bifurcation analysis of the representative parameter set
In the normal bifurcation analysis, each parameter was gradually changed from 0.001 to 10
times its original value and the solutions were automatically classified into four categories by
the criteria (Table S5). In the two-dimensional bifurcation analysis, two parameters were
gradually changed in certain area, and the solutions were automatically classified into four
categories by the criteria (Table S5).

Bifurcation analysis of all parameter sets
When the channel conductance of excitatory synaptic currents (IAMPA, INMDA) decreased, thesolution sometimes showed a cyclic firing pattern but weak synaptic currents. These
parameter sets were not classified into the ’True awake’ based on the previous study (1). The
definition of the ’True awake’ is described below.
1. If a model contains the AMPA receptors and does not contain the NMDA receptors, asolution is defined as the ’True awake’ when it satisfies the following conditions:・It is classified into the ‘Awake’ described above
・Maximum current of the AMPA receptors in 5 seconds to 10 seconds is equal or bigger thanhalf of its original value
2. In the other cases, a solution is defined as ’True awake’ when it satisfied the following
conditions:・It is classified into the ‘Awake’ described above
・Maximum current of the AMPA or NMDA receptors in 5 seconds to 10 seconds is equal or
bigger than half of its original value
In the bifurcation analysis of all parameter sets, the number of the solutionsclassified into ‘True Awake’ was counted as follows. Let 𝑋 be a channel, and 𝑠 be aparameter set and 𝐵(𝑋, 𝑠) was defined as follows. First, the conductance of the channel 𝑋of a parameter set 𝑠 was changed to have 0.001 times its original value. When the real
solution could be solved under the defined computer memory usage (3GB) and calculation

time (400 seconds) and it satisfied the condition of ‘True awake’, we defined 𝐵(𝑋, 𝑠) = 1.
When the real solution could not be solved under the defined computer memory usage (3GB)and calculation time (400 seconds) or at least one of the value |𝑣𝑖 | (𝑖 = 1, 2, … , 5001) was
bigger than 200 mV, the conductance of the channel X was gradually decreased from its
original value until the solution failed to be solved. One of those solutions satisfied thecondition of ‘True awake’, we defined⁡𝐵(𝑋, 𝑠) = 1. Otherwise, we defined 𝐵(𝑋, 𝑠) = 0. Let
𝑃 be the set of all parameter sets. We defined the bifurcation proportion of the channel X as
follows:
𝑅𝑋 = ∑ 𝐵(𝑋, 𝑠)𝑠∊𝑃

The 𝑅𝑋 values of all channels were plotted (Figs. 1 H-J, 2I, and 3E and Figs. S1 A-C andS3 I and J).

Parameter sets of the original AN model, the KO of NMDA receptors model, and the
KO of voltage-gated Ca2+ channels model
In this study, the parameter sets with the SWS firing patterns used in the previous study (1)were re-analyzed in Mathematica 10.0. Some parameter sets were not classified into ‘SWS’or failed to be solved in the bifurcation analysis due to the difference between Mathematica
10.0 and 9.0. Excluding these parameter sets, 1091, 1167, or 948 parameter sets wereanalysed in the original AN model, the KO of NMDA receptors model, or the KO ofvoltage-gated Ca2+ channels model, respectively.

Phase space
We fixed the Ca2+ concentration in Fig. 2G and Fig. S3E as follows:𝑑[𝐶𝑎2+ ]=0𝑑𝑡
The other equations were the same as the original equations. The nullclines of 𝑉 and 𝑛𝐾(Fig. 2 G and H and Fig. S3 E-H) are described as follows:
𝐼𝐿 (𝑉) + 𝐼𝐾 (𝑉, 𝑛𝑘 ) + 𝐼𝐶𝑎 (𝑉) + 𝐼𝐾𝐶𝑎 (𝑉, 𝐶𝑎𝑓𝑖𝑥) + 𝐼𝑁𝑎𝑃 (𝑉) = 0𝛼𝑛 (𝑉)(1 − 𝑛𝐾 ) − 𝛽𝑛 (𝑉)𝑛𝐾 = 0

𝐶𝑎𝑓𝑖𝑥 value is as follows:
Figure𝐶𝑎𝑓𝑖𝑥

2G

7, 8.5, 10.6, 12.5

2H

19.4068

S3E

5, 25, 50, 70

S3F

96.5712S3G

12.6475

S3H

74.5802

The Ca2+ concentration was not fixed in Fig. S3C. The nullclines of 𝑉 and 𝑛𝑘 and [𝐶𝑎2+ ]are described as follows:
𝐼𝐿 (𝑉) + 𝐼𝐾 (𝑉, 𝑛𝑘 ) + 𝐼𝐶𝑎 (𝑉) + 𝐼𝐾𝐶𝑎 (𝑉, [𝐶𝑎2+ ]) + 𝐼𝑁𝑎𝑃 (𝑉) = 0𝛼𝑛 (𝑉)(1 − 𝑛𝐾 ) − 𝛽𝑛 (𝑉)𝑛𝐾 = 0−𝛼𝐶𝑎 𝐴𝐼𝐶𝑎 (𝑉) −[𝐶𝑎2+ ]=0𝜏𝐶𝑎

Linear stability analysis
To elucidate whether each fixed point is stable or not in the SAN model, we conducted thelinear stability analysis. The function 𝑓1 , 𝑔1 are defined as below.
1𝑓1 (𝑉, 𝑛𝑘 ) = − (𝐼𝐿 (𝑉) + 𝐼𝐾 (𝑉, 𝑛𝑘 ) + 𝐼𝐶𝑎 (𝑉) + 𝐼𝐾𝐶𝑎 (𝑉, 𝐶𝑎𝑓𝑖𝑥) + 𝐼𝑁𝑎𝑃 (𝑉))𝐶
𝑔1 (𝑉, 𝑛𝑘 ) = 4(𝛼𝑛 (𝑉)(1 − 𝑛𝐾 ) − 𝛽𝑛 (𝑉)𝑛𝐾 )In the condition that the Ca2+ concentration was fixed, the fixed point (𝑉 ∗ , 𝑛𝑘 ∗ ) was
calculated by solving 𝑓1 (𝑉, 𝑛𝑘 ) = 𝑔1 (𝑉, 𝑛𝑘 ) = 0. Then the Jacobian matrix at (𝑉 ∗ , 𝑛𝑘 ∗ ) was
described below.
𝜕𝑓1 ∗ ∗(𝑉 , 𝑛𝑘 )𝜕𝑉𝐽=𝜕𝑔1 ∗ ∗(𝑉 , 𝑛𝑘 )( 𝜕𝑉𝜕𝑓1 ∗ ∗(𝑉 , 𝑛𝑘 )𝜕𝑛𝑘
𝜕𝑔1 ∗ ∗(𝑉 , 𝑛𝑘 )𝜕𝑛𝑘
)The eigenvalues of the matrix 𝐽 was described in Table S6.

Animals and sleep phenotyping
All experimental procedures and housing conditions involving animals and their care wereapproved by the Animal Care and the Use Committee of Graduate School of Medicine, theUniversity of Tokyo, and all of the animals were cared for and treated humanely inaccordance with the Institutional Guidelines for Experiments using animals. All mice were
given food and water ad libitum, and kept at an ambient temperature and humidity. The light
was controlled under 12-hr light/12-hr dark cycle. Sleep phenotyping was conducted at
8-weeks old in Snappy Sleep Stager (SSS) (3) by using C57BL/6N mice (n = 101) as control
information (Fig. 4 B, C, E, and F). SSS is a respiration-based fully automated non-invasive
sleep phenotyping system, whose accuracy was reported to be 95.3±0.4% (3). In SSS, a
mouse was placed in a chamber connected to a sensor, which detected the pressure difference
between the outside and inside of the chamber. The respiration time domain data wasconverted by fast Fourier transformation and the power spectrum was then subjected to theprincipal-component analysis. The probability density clustering was conducted by using thefirst three principal components. Finally, each cluster was annotated according to the power
spectrum of each cluster. The detailed method was described in the previous study (3).

Design of target sequences for guide RNA (gRNA)
The target sequences for Kcnk1, Kcnk2, Kcnk3, Kcnk4, Kcnk5, Kcnk6, Kcnk7, Kcnk9, Kcnk10,
Kcnk12, Kcnk13, Kcnk15, Kcnk16, Kcnk18, and Kcnk9 (set 2) (Fig. S5 A-O) were selectedfrom the list resulting from the mm10 CRISPR/Cas9 database (3)(http://www.crispr.riken.jp/) or selected manually by using the on-line tools including the
CRISPR guide RNA Design tool (4) (http://cas9.cbi.pku.edu.cn/index.jsp), the CRISPR
Design Tool (5) (http://tools.genome-engineering.org), the CRISPR-ERA: a comprehensive
designer tool for CRISPR genome editing, (gene) repression, and activation(http://crispr-era.stanford.edu/), the CRISPRdirect (6) (http://crispr.dbcls.jp/), and the
UNAFold Web Server (http://unafold.rna.albany.edu/).

gRNA synthesis
The gRNA templates for Kcnk1, Kcnk2, Kcnk3, Kcnk4, Kcnk5, Kcnk6, Kcnk7, Kcnk9,
Kcnk10, Kcnk12, Kcnk13, Kcnk15, Kcnk16, Kcnk18, and Kcnk9 (set 2) (Fig. S5 A-O) weredirectly synthesized and fused to the T7 promoter by PCR. First, partial fragments of thegRNA templates including each target sequence were amplified from the pX330 plasmids(Addgene, #42230) by PCR with the Common Reverse primer

(5'-AAAAGCACCGACTCGGTGCC-3', Fasmac Co., Ltd.) (7) and Forward primer-1 (Table
S8, Fasmac Co., Ltd.) were used for each target sequence. The T7 promoter-fused gRNA
templates were amplified from the diluted PCR products by PCR with the Common Reverse
primer and Forward primer-2 for each target sequence (Table S8, Fasmac Co., Ltd.). Further
information on gRNA synthesis was published previously (3).

Cas9 mRNA synthesis
p3s-Cas9HC (8) (Addgene, #43945), which includes a T7 promoter-fused Cas9 coding
region, was digested with XbaI (TaKaRa), and used as the template for in vitro transcription
using the mMESSAGE mMACHINE T7 kit (Life Technologies). The Cas9 mRNA was
puriﬁed using the MEGAclear kit (Life Technologies).

One-cell embryo microinjection
C57BL/6N females (4-6 weeks old, CLEA Japan Inc) were superovulated and mated with
C57BL/6N males (CLEA Japan Inc). Cas9 mRNA (100 ng/µl) and gRNAs (150 ng/µl in
total) were co-injected into the cytoplasm of fertilized eggs in M2 medium (ARK Resource)
at room temperature. Further information on one-cell embryo microinjection was publishedpreviously (1, 3).

Genotyping of knockout mice by quantitative PCR (qPCR)
Genotyping was performed by using qPCR. Primers for qPCR (Table S9, Fasmac Co., Ltd.)
were annealed to the targeting sequences. The absolute target site abundance was calculatedusing a standard curve obtained from wild-type genomic DNA. The amount of Tbp (9) was
quantified and used as an internal control. When the amplified intact DNA by qPCR is lessthan 0.5% of wild-type genome, we judged that the target DNA is not detectable. When either
three targets of the gene were not detectable, we classified the animal as a KO. When we
could not confirm KO genotype by qPCR, we performed 2nd qPCR using the alternative
primer which was independent of 1st qPCR. Further information on genotyping was publishedpreviously (1, 3).

EEG/EMG recording
The recorded EEG/EMG data were analysed by the FASTER method. All results werechecked manually. Sleep phenotyping was conducted among 10-14 weeks old. Further
information on implantation of telemetry devices and data analysis of EEG/EMG recordings
has been published previously (10).

Statistical analyses
Statistical analyses were performed by R version 3.4.1. or 3.4.3.
For comparing more than two samples (i.e., CRISPR mutant mice) against identical
one sample (i.e., C57BL/6N control mice), the Steel test, a nonparametric multiple
comparison test, was applied after the confirmation of nonnormality and unequal variance by
the Kolmogorov-Smirnov test and the Bartlett's test, respectively.
For paired samples (i.e., mice in EEG/EMG recording), the normality was tested by
the Shapiro test at a significance level of 0.05. When normality was rejected in each group, a
two-sample Wilcoxon t-test was applied. When normality was not rejected in both groups,
the homogeneity of variance was tested by an F-test at a significance level of 0.05. When thenull hypothesis of the normal distribution with the equal variance for the two groups was not
rejected, a Student's paired t-test was used.
In this study, P < 0.05 was considered significant (*P < 0.05, **P < 0.01, ***P <0.001, and n.s. for not significant).

Calculation of the EEG spectrum and the delta power
The recorded EEG data was normalized by the standard deviation of the EEG voltage on thefirst day. This data was analyzed by fast Fourier transformation and the average power
spectrum during NREM sleep was plotted (Fig. S6C). The delta power is defined as the sumof delta frequency (0.5-4 Hz) power (Fig. S6D)

Definitions of sleep/wake parameters
Sleep duration is an averaged total sleep duration per day over six days.
Pws and Psw are the transition probability from an awake state to a sleep one and thetransition probability from a sleep state to an awake one, respectively. The mathematicaldefinition is as follows. s and w denote Sleep and Wake respectively in this section. Let 𝐸𝑖
be the state of the i-th epoch, which is s or w. Let 𝑁𝑋𝑌 (𝑋, 𝑌 ∊ {𝑠, 𝑤}) be the number of
elements of the set {(𝐸𝑖 , ⁡𝐸𝑖+1 )⁡|𝐸𝑖 = 𝑋, 𝐸𝑖+1 = 𝑌⁡}. Pws and Psw are defined as 𝑁𝑤𝑠 /(𝑁𝑤𝑠 +𝑁𝑤𝑤 ) and 𝑁𝑠𝑤 /(𝑁𝑠𝑤 + 𝑁𝑠𝑠 ), respectively.
The sleep/wake relative amplitude is defined as the coefficient of variation (CV, thestandard deviation divided by the mean) of sleep time for each 10-min bin for 24 hours.

Definitions of transition parameters in EEG/EMG
n, r, and w denote NREM, REM and Wake respectively in this section. Let 𝐸𝑖 be
the state of the i-th epoch, which is n, r, or w. Let 𝑁𝑋𝑌 (𝑋, 𝑌 ∊ {𝑛, 𝑟, 𝑤}) be the number of
elements of the set {(𝐸𝑖 , ⁡𝐸𝑖+1 )⁡|𝐸𝑖 = 𝑋, 𝐸𝑖+1 = 𝑌⁡}. The transition probability from 𝑋 to 𝑌
(𝑋, 𝑌 ∊ {𝑛, 𝑟, 𝑤}) is defined as 𝑁𝑋𝑌 /(𝑁𝑋𝑛 + 𝑁𝑋𝑟 + 𝑁𝑋𝑤 ).

Fig. S1. The membrane potential and the Ca2+ concentration of the various models,related to Fig. 1(A-C). Summary of the bifurcation analysis in each condition: Model 1 (A), Model 2 plus
gAMPA model (B), or Model 2 plus gA model (C). The conductance and the time constantrelated to the Ca2+-dependent hyperpolarization pathway are in red.
(D-H) The example of membrane potential and intracellular Ca2+ concentration (upper), thestandardized membrane potential, and the standardized intercellular Ca2+ concentration(lower) in Model 1 (D), Model 2 plus gNa model (E), Model 2 plus gGABA model (F), Model 2plus gAMPA model (G), or Model 2 plus gA model (H). The number of the parameter sets used
for calculating the standardized membrane potential and the intercellular Ca2+ concentration
was 279 for Model 1 (D), 86 for Model 2 plus gNa model (E), 137 for Model 2 plus gGABA
model (F), 19 for Model 2 plus gAMPA model (G), and 20 for Model 2 plus gA model (H),respectively. The line and shadow (lower) indicate the mean and standard deviation,
respectively.

Model 1 (NMDAR KO Model)
22

-80

60

5

40

-2

aP

g

A

Time (sec)Intracellular Calciumn = 86
0

1

2

5

-2

Time (sec)Membrane Potential
n = 137
-1

0

1

0

-3
6

0

5

Time (sec)Intracellular Calciumn = 19
-1

0

1

5

-2

2

Time (sec)Membrane Potential
n = 20
-1

0

1

10

2

0

-3

2

6

Standardized
MembranePotential

0

0

-3

-3
0

6

0

Fig. S1.

5-2Time (sec)Intracellular Calciumn = 20
-1

0

1

10

2

3

0

-3

6
Bursting cycles

Bursting cycles

0

Parameter sets
3

3

6

1

20

Standardized Ca2+

-3

0

0

0

-80

10

Parameter sets
0

-1

10

3

Intracellular
Calcium
(µM)

0

-2

Standardized Ca2+

3

-2Time (sec)Intracellular Calciumn = 137

Model 2
+ gA

25

2

Parameter sets

1

10

5

Bursting cycles

Parameter sets

Time (sec)Membrane Potential
n = 190

6H

Intracellular
Calcium(µM)

MembranePotential(mV)

20Standardized
MembranePotential

6

Parameter sets

-30Model 2
+ gAMPA

2

Standardized Ca2+Standardized
MembranePotential

0

MembranePotential
(mV)

G

10

3

Bursting cycles

0

-3

20

Parameter sets

-1

-100

10

3

6

0

0

Intracellular
Calcium(µM)

5Parameter sets
0

-1

2

Model 2
+ gGABA
20-2Standardized Ca2+Standardized
MembranePotential

-3

-2

1

0

MembranePotential(mV)

Intracellular
Calcium(µM)

MembranePotential(mV)

2Parameter sets

1

0

5

0

3

6

F

10

10

3

-100

-1

Bursting cycles

Time (sec)Membrane Potential
n = 860Standardized Ca2+
0

50

-1

0

-3

0

a

0

N

C

g

ΤC

AM

0

L

0

g

0
a

PA

0

50

-2

-2

10

Parameter sets

Parameter sets
Standardized
MembranePotential

20

Model 2
+ gNa

5

2

Time (sec)Intracellular Calciumn = 279

AB
A
G

g

40

g

0

E

-100

1

3

60

g

L

g

N

g

0
aP

0

a
KC

a

0

g

C

g

ΤC

a

0

0

5

80

a

20

100

KC

40

Model 2
+ gA

g

60

0

A

Bifurcation from SWS to Awake (%)

80

0

AR

C

Model 2
+ gAMPA

100

g

K

g

g

g

KS

N
aP

L

N
a

g

g

a

KC
a

C
a

g

g

ΤC

B

-1

10

10

0

g

0
PA

0

AM

20

Time (sec)Membrane Potential
n = 279

Intracellular
Calcium(µM)

MembranePotential
(mV)

40

80

0

Bifurcation from SWS to Awake (%)

D

Model 1 (NMDAR KO Model)
100g

Bifurcation from SWS to Awake (%)

A

0

6

Fig. S2. Detail analysis of the SAN model, related to Fig. 2(A) The histogram and the kernel density estimator about the conductance of each channel of
1,008 parameter sets obtained by the parameter searches. Let 𝜎′𝑖 and 𝑛′ be the uncorrected
sample standard deviation of each parameter and the number of parameter sets (i.e.
𝑛′=1,008) respectively. The bandwidth 𝑏′𝑖 of each parameter was set as follows:1

4 5
𝑏′𝑖 = 𝜎′𝑖 { } .
3𝑛′All 1,008 parameter sets were divided into cluster 1 with 480 parameter sets and cluster 2
with 528 parameter sets.
(B) All parameter sets projected on the principal component space. Cluster 1 and cluster 2 areindicated as blue and red points, respectively.
(C) Membrane potential, proportion of each channel’s inward current, and proportion of eachchannel’s outward current of the representative parameter set in each cluster over time.
(D) The transition from the SWS (red) to awake (blue) firing patterns by the downregulationof the conductance of Ca2+-dependent K+ channels (gKCa) and voltage-gated Ca2+ channels(gCa), and the time constant for Ca2+ efflux (τCa) in the representative parameter set in eachcluster.
(E) The examples of SWS, Awake, Resting, and Slow-wave with few spikes.

A

gL

gK

0.20

0.20

0.10

0.00

-2

-1

0

1

2

0.20

-2

-1

0

1

2

0.20

Cluster 1 0.10

0.15

0.00

-2

-1

0

1

2

0.00

0.20

Cluster 2 0.10

0.15

0

1

2

-2

-1

0

1

2

-2

-1

0

1

2

0.15

-1

0

1

log each value

2

Cluster 2

-2

-1

0

1

2

log each value

0.00

-2

-1

0

1

log each value

0.05

0.00

0.00

-2

-1

0

1

2

0.06

0.10

0.05

0.08

2

2

1.0

1.5

2.0

2.5

3.0

0.00

0.10

0.02

0.05

0.00

0.00

0

1

2

0.06

0.10

0.05

0.08

1.0

1.5

2.0

2.5

3.0

0.00

0.10

0.02

0.05

0.00

0.00

0

1

2

-2

-1

0

1

2

-2

-1

0

1

2

0.15

0.01

log each value

2

0.20

0.04

-1

1

0.25

0.06

-2

0

0.15

0.01
-1

-1

0.20

0.04

-2

-2

0.25

0.06

C

Cluster 1

Proportion of
inward current
(%)

-2

0.00

0.10

0.02

0.05

0.15

0.02

0.03

0.05

0.20

0.01

0.04

0.10

0.25

0.04

0.02

0.00

gNaP

0.06

0.04

0.05

0.10
0.05

0.08

0.03

0.25

0.15

B

-1

0.05

0.20

0.00

-2

0.10

0.10
0.05

0.00

0.15

0.25

0.15

0.05

0.02

0.05

0.05
0.00

0.10

0.03

0.10
0.05

0.06

0.04

0.10

0.15

ΤCa

gKCa

0.15

0.25

0.15

All

gCa

1.0

1.5

2.0

2.5

log each value

3.0

0.00

log each value

Cluster 1 (representative)

100

Proportion ofoutward current(%)

0

Proportion of
inward current
(%)

PC2

6

6.5

6

6.5

Cluster 2 (representative)

100

-2

3

PC1

Proportion ofoutward current
(%)

-2

6.2

6.4

6.2

6.4

Time (sec)
gNaP

gNaL

gCa

gK

gKCa

gKL

E

Cluster 1 (representative)

gKCa (Original = 1)

10-0.60-20

gCa (Original = 1)

10-0.150

5Time (sec)

10 5Time (sec)

-80

10

ΤCa (Original = 1)

10-0.3150

SWSTime (sec)

10 5Time (sec)(gCa = Original)
-20

0
-20

-80

Awake

(gCa = 10-0.3×Original)-80

5Cluster 1 (representative)

10-0.2

MembranePotential(mV)

MembranePotential
(mV)

-20

MembranePotential
(mV)

-80

10-0.25150

Peak Freq.
(Hz)

10-0.7200

Peak Freq.
(Hz)

Peak Freq.
(Hz)

D

0

10

5

Time (sec)

10 5Time (sec)

10

-20

-80

Resting(gCa = 10-1×Original)

gCa (Original = 1)

100

5Time (sec)

10 5Time (sec)

10

0
20

-100

10-0.65

ΤCa (Original = 1)

10-0.55200

MembranePotential
(mV)

0
20

-100

10-0.05
300

Peak Freq.
(Hz)

100

-20

-80

Slow-wave with few spikes(gCa = 100.5×Original)

0
20

-20

MembranePotential(mV)

gKCa (Original = 1)

Peak Freq.
(Hz)

10-0.1250

MembranePotential(mV)

Peak Freq.
(Hz)

Cluster 2 (representative)

5

Time (sec)

10 5Time (sec)

10

-100

Fig. S2.

5Time (sec)

10 5Time (sec)

10

-80

Fig. S3. Leak K+ channels play a role in generating the SWS firing pattern, based on the
SAN model, related to Fig. 2(A-B) Two-dimensional bifurcation analysis of the representative parameter set in cluster 1(A) or in cluster 2 (B). The definitions and examples of SWS, Awake, and Resting aredescribed in Table S5 and Fig. S2E, respectively.
(C) The three-dimensional phase space of the representative parameter set in cluster 1 (upper)
or cluster 2 (lower). The red line, the orange plane, the blue plane, and the green plane
indicate the trajectory of the differential equation, the V nullcline, the nK nullcline, and theCa2+ nullcline, respectively.
(D) Three-dimensional (upper) and two-dimensional (lower) plots of the trajectory of therepresentative parameter set in cluster 2.
(E) Phase planes of the representative parameter set in cluster 2 in the case where the Ca2+
concentration is fixed at 5 μM, 25 μM, 50 μM, or 70 μM. The red, green, and blue lines
indicate the trajectory of the differential equation, V nullcline, and nK nullcline, respectively.
(F) Membrane potential and phase plane of the parameter set whose gKCa was lower (left) or
whose gKCa and gKL were lower (right) than the representative parameter set in cluster 2.
(G-H) Membrane potential and phase plane of the parameter set whose gKCa was normal (left)
or gKCa was lower (right) than the representative parameter set in each cluster, cluster 1 (G)and cluster 2 (H).
(I) Summary of the bifurcation analysis in cluster 2. Conductance of the leak K+ channels isin red.
(J) Summary of the bifurcation analysis in each cluster, where the leak K+ channels and theleak Na+ channels are represented as one type of leak channels.

AC

Cluster 1 (representative)

-0.70
-0.30

log g Ca value (Original = 1)

-0.30

-0.30
-0.30

0

Cluster 2 (representative)

log g Ca value (Original = 1)

SWS

Awake

0
80

2+

[Ca ]
(µM)

-0.30
-0.70

0

V
(mV)
20

nK 1.0 -80

0log ΤCa value (Original = 1)0log ΤCa value (Original = 1)

log g KCa value (Original = 1)

0

D

Cluster 1 (representative)

log g KCa value (Original = 1)0

nK[Ca2+]

Cluster 2 (representative)

V

Resting
0

B

1.0

-0.50

log ΤCa value (Original = 1)

-0.10
-0.05

log g Ca value (Original = 1)

0

-0.70
-0.05

0

log g Ca value (Original = 1)

SWSF

log g KL value (Original = 1)

0-80
V nullcline

0

1.0

1.0

20
V
nK nullcline(mV)
2+[Ca ] =50 (µM)

log g KCa value (Original = 1)

0

0

-80

-0.035gKL = representativeG×representative)gKL = 10-1 ×representative

20

V
(mV)
2+
[Ca ] =70 (µM)

-80

20

V
(mV)

Cluster 1 (representative)
gKCa = 10-3 ×representative
gKCa = representative

20

-20

-20

V
(mV)

V
(mV)

20

-80

0

20

V
(mV)

Awake

Cluster 2 ( gKCa =10

0

-80

V
(mV)

20

V
(mV)

I

5

V
(mV)

nK
0

-50

-80

J

Cluster 2

10Time (sec)

0.6

nK

0
-60

20

V
(mV)

10

Time (sec)

0.6

1.0

-80

-80
5

0

-20

V
(mV)

-80

Cluster 1

-20

V
(mV)

Cluster 2

40
20

aP
N

g

g

L

g

C

g

K

0
a

0
a

0

ΤC

N

g

aP

g

a
L

a

KC

g

C

g

0
K

0

0
a

Fig. S3.

0

g

20

60

a

40

80

KC

60

100

g

80

Bifurcation from SWS to Awake (%)

100

ΤC

K

g

aL

aP

N

N

g

0

0

g

V
(mV)

20

0
a

-80

0
KL

0

20

g

V
(mV)

20

40

a

-80

60

KC

nK

nK

1.0

80

g

10

a

1.0

Time (sec)

C

-100
5
10

g

Time (sec)

100

ΤC

V
(mV)
5

Bifurcation from SWS to Awake (%)

20

20

0

0

-50

10

Time (sec)

nK

0
-60

Cluster 2 (representative)
gKCa = representative gKCa = 10-3 ×representative

-100

-80
5

1.0

nK

1.0

nK

1.0

-100

10

Time (sec)nK

5Bifurcation from SWS to Awake (%)

-100

H

nK

0

0

-0.10
-0.70

0

2+

[Ca ] =25 (µM)

1.0

nK

1.0

-0.65
-1.00

0

20

V
(mV)

nK

log ΤCa value (Original = 1)

0

log g KL value (Original = 1)

-80

Cluster 2 (representative)
2+[Ca ] =5 (µM)

nK

-0.04
-1.00

0

log g KL value (Original = 1)0.0Elog ΤCa value (Original = 1)

-0.07
-1.00

log g KCa value (Original = 1)

nK

-0.02

-0.04

log g Ca value (Original = 1)

log g KCa value (Original = 1)

Cluster 2 (representative)

Fig. S4. Current analyses in the AN model, related to Fig. 3(A) Normalized proportion of the other outward currents except for the current of leak K+
channels or Ca2+-dependent K+ channels of 1,091 parameter sets obtained by the parametersearches in the original AN model. The line and shadow (lower) indicate the mean andstandard deviation, respectively.
(B) Normalized proportion of the current of leak Na+ channels or NMDA receptors, or theother inward currents of 1,091 parameter sets obtained by the parameter searches in theoriginal AN model.
(C) Membrane potential, proportion of each channel’s inward current, and proportion of eachchannel’s outward current in the original AN model.
(D-E) Proportion of each channel’s inward current in the down state (D) or the up state (E)under the original AN model, the NMDAR KO model, or the voltage-gated Ca2+ channel KOmodel.

aA

Fig. S4.

A

AB

G

g

AR

g

A

g

KS

g

K

g

A

PA

A

G
AB
A

AR

g

g

PA

KS

g

g

K

g

AM

g

aP

N

g

a

aL

N

g

N

g

KL

g

A

a

0

AM

g

aP

N

g

a

0

N

20

0

g

40

Current proportion (%)

Parameter sets

0

aL

60

6

KL

80

D

100

N

Inward current in the up state

Voltage-gated Ca2+-channel KO Model

100

KC

Parameter sets

Original Model
NMDA receptor inward current proportionn = 1,091g

NMDAR KO Model

g

100

g

100

a

0

M

N

NMDAR KO Model

KC

g

20

A

40

D

60
100

M

Time (sec)g

80

Inward current in the down state (%)

100

Inward current in the up state (%)

AB
A

Bursting cycles

N

G

g

Current proportion (%)

0

g

A

AB

G

g

A

AR

g

0

AR

g

KS

g

0

g

0

KS

0

g

PA

5.7

K

5.7

g

5.7

g

PA

Proportion ofoutward current
(%)

0

K

Proportion ofoutward current
(%)

6

g

AM

Proportion ofoutward current
(%)

100

AM

100

g

a

aP

N

100

g

100

aP

-80

N

20

g

aL

N

Bursting cyclesaOriginal ModelN

Original Modelg+

Leak Na -current proportionn = 1,091

N

0

g

gGABA

aL

20

a

gCa

N

40

KL

C

g

60

g

5.0

g

80

KL

E

g

0

a

20

C

40

KC

60

a

D

KC

80

g

gNMDA

g

gNaL

a

0

C

0

g

100

Inward current in the down state (%)

Parameter sets
0

g

100

Inward current in the up state (%)

A

AB

G

Current proportion (%)

B

g

A

6

AB

g

Parameter sets

100

G

g

A

gAMPA

AR

g

Current proportion (%)
100

AR

gNaP

g

g

KS

g

gNa

KS

PA

Bursting cycles

g

K

g

AM

g

MembranePotential(mV)
0

Original Modelthe other outward current proportionn = 1,091

g

K

PA

a

aP

N

g

5.0

g

AM

g

aP

N

g

N

aL

KL

g

g

5.0

N

g

a

N

g

N

aL

KL

g

A

a

KC

g

a

D

M

C

g

N

g

Proportion of
inward current
(%)
0

g

A

a

KC

g

D

M

C

0

g

Proportion of
inward current
(%)
0

N

Inward current in the down state (%)
0g

Inward current in the up state (%)

A
the other inward current proportionn = 1,091
100

100

Bursting cycles
6

Original Model
100

100

100
5.0

gK
gKCa
gKL
5.7

5.0

gKS
gA
gAR
5.7

5.0

gAMPA
gNMDA
gGABA
5.7

Time (sec)Inward current in the down state

Voltage-gated Ca2+-channel KO Model

80

60

40

20

0

80

60

40

20

0

Fig. S5. Impairment of the leak K+ channels decreases sleep duration in mice, related toFig. 4(A-O) Target sequences of the gRNAs for KO mice of the leak K+ channel family [Kcnk1 (A),Kcnk2 (B), Kcnk3 (C), Kcnk4 (D), Kcnk5 (E), Kcnk6 (F), Kcnk7 (G), Kcnk9 (set 1, H),
Kcnk10 (I), Kcnk12 (J), Kcnk13 (K), Kcnk15 (L), Kcnk16 (M), Kcnk18 (N), and Kcnk9 (set 2,O)]; each gene had three target sequences. Mouse genomic sequence data were obtained from
GRCm38/mm10 via the UCSC Genome Browser (11) (http://genome.ucsc.edu/). Colored
letters (blue, orange, and green) show the 20-base target sequences. The target sequenceswere designed on the sense (+) or antisense (-) strand of genomic DNA.
(P-AD) The genotype of KO mice for the leak K+ channel family [Kcnk1 (P), Kcnk2 (Q),Kcnk3 (R), Kcnk4 (S), Kcnk5 (T), Kcnk6 (U), Kcnk7 (V), Kcnk9 (set 1, W), Kcnk10 (X),
Kcnk12 (Y), Kcnk13 (Z), Kcnk15 (AA), Kcnk16 (AB), Kcnk18 (AC), and Kcnk9 (set 2, AD)];The relative amount of intact DNA for each target sequence was measured by quantitative
PCR (qPCR). The genomic DNA was purified from the tail of each mouse. The relative
amount of intact DNA for each target sequence was scaled so that the level in a wild-type
mouse was defined as 100%. The KO mice, which are confirmed by 1st qPCR or 2nd qPCR,
are labeled in red.

B

Kcnk1
chr8(+)

126,000,000

126,030,000

C

Kcnk2
chr1(-)

189,300,000

189,250,000

gtcccaattccaatttcccg (+)
Kcnk1 target sequence-2

D

attgtccatgccgttctgct (-)
Kcnk1 target sequence-3

atctccccacgaactgaagg (-)
Kcnk2 target sequence-1

gttatcactctgacgaccat (-)
Kcnk2 target sequence-2

E

Kcnk4
chr19(-)

6,930,000

G

5,705,500

F

aaaaagccttccagaccccg (-)
Kcnk4 target sequence-3

accgatactttgtggagctt (-)
Kcnk5 target sequence-1

gctgaagatgttgacgtcct (+)
Kcnk5 target sequence-2

chr14(-)

H

5,707,000

J

acccagccatttaccacctt (+)
Kcnk7 target sequence-2

chr7(-)

gaagtcagctcgatcgtgtt (+)
Kcnk7 target sequence-3

72,540,000

ggtttggtaactatccttgg (+)
Kcnk5 target sequence-3

ggtgtccagtggaaattcgc (-)
Kcnk9 target sequence-2

40

40

20

20

20

0

0

0

60

60

40

40

20
0

wt 1 2 3 4 5 6 7

wt 5 6

Genome DNA/ Tbp(wild-type =100%)

60
40

0

Kcnk6

80

60

60

40

40

20

20

20

20

0

0

0

0

wt 1 2 3 4 5 6 7 8 9 10 11 12

Genome DNA/ Tbp(wild-type =100%)

Genome DNA/ Tbp(wild-type =100%)

Kcnk16 (1st qPCR)

80
60
40

60
40
20
wt 1 2 3 4 5 6 7 8 9 10 11 12 13 14

Kcnk9 (set1, 1st qPCR)
120

100
80
60
40

120

100
80
60
40
20

wt 1

2

3

4

0

5

wt 1 2 3 4 5 6 7 8

Kcnk15

80
60
40

80
60
40

0

wt 1 2 3 4 5 6 7

AD

100

80

80

60

60

40

40

20
0

wt 1 2 3 4 5 6 7 8 9 10 11 12

Kcnk9 (set2, 1st qPCR)

120

120

100

Fig. S5.

100

20

Kcnk18(2nd qPCR)

wt 1 2 3(1st qPCR)

120

100

0

Kcnk18(1st qPCR)

wt 2

AA

wt 1 2 3 4 5 6 7 8AC0

wt 1 2 3 4 5 6 7 8

20Genome DNA/ Tbp(wild-type =100%)

80

20

120

100

0

40

20

Kcnk13 (1st qPCR)

20
wt 1 2 3 4 5 6 7

60

40

120

Z
120

100

0

wt 3 5 9

Kcnk12 (1st qPCR)

20

80

60

W

100

40

80

Kcnk7 (1st qPCR)

nd

80

60

100

0

100

80

100

V

(2 qPCR)

Kcnk4(2nd qPCR)

120

wt 1 2 3 4 5 6

120

100

120

40

120

120

AB

60

wt 1 2 3 4

Y

0

80

gaataccgatgaggagctac (-)
Kcnk9 target sequence-6

120

20

Kcnk6 (1st qPCR)

Kcnk10 (1st qPCR)
Genome DNA/ Tbp(wild-type =100%)

80

100

ctgtatgttctacgctgtgc (-)
Kcnk9 target sequence-5

Kcnk4 (1st qPCR)

120

100

U

wt 4 7

S

Kcnk3 (1st qPCR)

Genome DNA/ Tbp(wild-type =100%)

80

72,520,000

aggaaccggcgaatttccac (+)
Kcnk9 target sequence-4

Genome DNA/ Tbp(wild-type =100%)

100

80

72,540,000

Kcnk9

aaagggaagtttcgctacct (+)
Kcnk18 target sequence-3

Genome DNA/ Tbp(wild-type =100%)100

actctcaccaccataggctt (+)
Kcnk15 target sequence-3

Kcnk9 (set2)
chr15(-)

Genome DNA/ Tbp(wild-type =100%)

60

120

ctccttctacttcgccatca (+)
Kcnk15 target sequence-2

O

R

120
Genome DNA/ Tbp(wild-type =100%)

80

60

120

cggagaagcggtactttctg (-)
Kcnk15 target sequence-1

59,235,000

cactgaccaccatcgggttt (+)
Kcnk18 target sequence-2

Kcnk2 (1st qPCR)

80

X

Genome DNA/ Tbp(wild-type =100%)

Q

100

Kcnk5
(2nd qPCR)

163,858,000
Kcnk15cacctgttccggatgttacc (-)
Kcnk13 target sequence-3

59,225,000

gatgtatgcaaagacgacca (-)
Kcnk18 target sequence-1

Kcnk1(2nd qPCR)

Genome DNA/ Tbp(wild-type =100%)

Genome DNA/ Tbp(wild-type =100%)
Genome DNA/ Tbp(wild-type =100%)

agcattacatcgctgtgtat (-)
Kcnk16 target sequence-3

120

Kcnk5 (1st qPCR)

163,855,000

Kcnk18100T

ctatccgaagtatggagcgc (-)
Kcnk10 target sequence-3

Kcnk15chr2(+)

Kcnk18chr9(+)

120

wt 1 2 3 4 5 6 7 8 9 10

gttcttctgactgctctcga (+)
Kcnk10 target sequence-2

L

aaaaatcttccctcccgttg (+)
Kcnk13 target sequence-2

N

Kcnk1 (1st qPCR)

agtggttgccatcttcgtgg (-)
Kcnk10 target sequence-1

97,550,000

Kcnk16P

Kcnk10ggaacataaccagcgtcagg (+)
Kcnk9 target sequence-3

97,500,000

ggtacagtggtttctaccat (-)
Kcnk13 target sequence-1

20,263,000

acgcgaagtagaagccctca (+)
Kcnk16 target sequence-2

98,450,000

Kcnk13

Kcnk16

gatggggatcccactcaatg (-)
Kcnk16 target sequence-1

ggtgctgatgtagatctgtc (+)
Kcnk6 target sequence-3

Kcnk1098,500,000

chr12(-)

Kcnk13chr12(+)

cgctgggcatcatgaacaac (-)
Kcnk12 target sequence-3

20,266,000

chr14(-)

actcatcttgctgcccgatc (-)
Kcnk6 target sequence-2

I

Kcnk12M

agtctgcaacacgagcacca (+)
Kcnk6 target sequence-1

72,520,000K

gctgttatctccgcaggttt (-)
Kcnk12 target sequence-2

29,225,000
Kcnk6

Kcnk9 (set1)
chr15(-)

ggtgaaggtacaggcgatca (+)
Kcnk9 target sequence-1

87,660,000

ctcacctatggtcgacacca (+)
Kcnk12 target sequence-1

29,230,000

Kcnk9

Kcnk12chr17(-) 87,700,000

attcatgaccatgaacgccg (+)
Kcnk3 target sequence-3

Kcnk6

20,150,000

Kcnk7

gaatggctacccagtcacct (-)
Kcnk7 target sequence-1

gaaaccgatgagcaccatgt (-)
Kcnk3 target sequence-2

Kcnk5

Kcnk7
chr19(+)

aacgcatcaacaccttcgtg (+)
Kcnk3 target sequence-1

Kcnk5
20,160,000

6,927,500

gtctgaggactcgtcgatga (+)
Kcnk4 target sequence-2

30,620,000
Kcnk3

gtgtggttctggatcctcgt (-)
Kcnk2 target sequence-3

Kcnk4

gaaatacctctgcccgagtt (+)
Kcnk4 target sequence-1

30,600,000

Kcnk2

Kcnk1

caccgacactccataattgc (+)
Kcnk1 target sequence-1

Kcnk3
chr5(+)

Genome DNA/ Tbp(wild-type =100%)

A

Kcnk9
(set2, 2nd qPCR)
120

100

100

80

80

60

60

40

40

20

20

20

0

0

wt 2

wt 1 2 3 4 5 6 70

wt 2 4

Fig. S6. Impairment of the leak K+ channels decreases sleep duration in mice, related toFig. 4(A) Distributions of sleep/wake parameter of the KO mice of leak K+ channels. The definitionof relative amplitude is described in SI Materials and Methods. Red dots and lines indicate
the mean and SEM, respectively. Black dashed line and gray shading: the mean and 1 SD
range from the WT data.
(B) The transition probability was shown in a bar graph. The definition of transition
probability is described in SI Materials and Methods. Error bars indicate SEM.
(C) The power spectrum of Kcnk9 KO mice. The mean power density at each frequency wasplotted for mutant mice (magenta line, n=7) and for control mice (black line, n=7). All thetime series of EEG data were normalized so that the standard deviation of the first 24 hours
becomes 1 before calculating the power spectrum. Shaded area represents SEM.
(D) The sum of power density between 0.5 Hz to 4.0 Hz was calculated and plotted formutant mice (magenta box) and for control mice (black box).
WT: C57BL/6N male mice.

*

0.02

Pws

Psw

0.02

0.01

0.01

**

0.00

**

0.2 0.4 0.6 0.8
Relative amplitude

C

0.100
0.075
0.050
0.025
0.000
WT

WT

Kcnk9

Kcnk9

Transition Probability from REM to Awake

Kcnk9(set2)

Transition Probability from REM to NREM**

0.2 0.4 0.6 0.8
Relative amplitude

0.00
WT

Kcnk9

0.06
0.04
0.02
0.00
WT

Kcnk9

0.020
0.015
0.010
0.005
0.000
WT

Kcnk9

WT

Kcnk9

Transition Probability from NREM to REM

Kcnk1
Kcnk2
Kcnk3
Kcnk4
Kcnk5
Kcnk6
Kcnk7
Kcnk9
Kcnk10
Kcnk12
Kcnk13
Kcnk15
Kcnk16
Kcnk18Transition Probability from Awake to NREM Transition Probability from NREM to Awake

B

A

0.015

0.010

0.005

0.000
WT

Kcnk9

0.02

0.01

0.00

D

Power Spectrum during NREM sleep

0.6
WT

0.2

Kcnk9

0.1

Delta power

power density (normalized mV^2/Hz)

0.3

WT
0.4

0.2

0.0

0.0
0

10

frequency (Hz)20

30

Fig. S6.

Kcnk9

Supplemental Tables
Table S1. Parameter values in the model
Values from C to τGABA except for VNaL were based on the previous study (1). The sodium
leak reversal potential was different from the sodium reverse potential because NALCN
channels’ reverse potential is near 0 mV (12). The values indicated by asterisks weredifferent in each simulation, which are listed in Table S2.

CMembrane capacitance

1 μF/cm2

A

Area of neuron

0.02 mm2

VL

Leak reversal potential-60.95 mV

VNaL

Sodium leak reversal potential

0 mV

VNa

Sodium reversal potential

55 mV

VK

Potassium reversal potential-100 mVτhA

Time-constant of hA

15 ms

2+

VCa

Ca reversal potential

120 mVKD

Dissociation constant of Ca2+-dependent K+
channels

30 μM

VAMPA

AMPA receptor reversal potential

0 mV

VNMDA

NMDA receptor reversal potential

0 mV

VGABA

GABAA receptor reversal potential-70 mVαCa

Coefficient of Ca2+-entry

0.5 μM/(nA ms)

τAMPA

Time-constant of sAMPA

2 ms

τsNMDA

Time-constant of sNMDA

100 msτxNMDA

Time-constant of xNMDA

2 ms

τGABA

Time-constant of sGABA

10 ms

gL

The conductance of leak K+ channels*

gNa

The conductance of voltage-gated Na+ channels

*

+

gK

The conductance of voltage-gated K channels*

gA

The conductance of fast A-type K+ channels*

gKS

The conductance of slowly inactivating K
channels

+

*

gNaP

The conductance of persistent Na+ channels*

gAR

The conductance of inwardly rectifying K
channels

+

*

gCa

The conductance of voltage-gated Ca2+ channels

*

gKCa

The conductance of Ca2+-dependent K+ channels

*

gAMPA

The conductance of AMPA receptors

*

gNMDA

The conductance of NMDA receptors

*

gGABA

The conductance of GABAA receptors

*

τCa

2+

The time constant of Ca efflux*

Table S2. Parameter sets in each simulation
The parameter set used in each simulation is as follows.

Figure 1L

Cluster1

Cluster2(representative)(representative)

Figure S4C

gL [mS/cm2]

10-1.7876

10-1.7982

10-1.60899

0.03573

gNa [mS/cm2]

0

0

0

12.24382

1.2834

1.269074

0.473475gK [mS/cm ]

10

10

10

2.61868

gA [mS/cm2]

0

0

0

1.79259

20

gKS [mS/cm ]20-0.19850-0.18345

0.0350135-0.13573

gNaP [mS/cm ]

10

10

10

0.0717984

gAR [mS/cm2]

0

0

0

0.0166454

2-0.7895-0.8362-0.0337

gCa [mS/cm ]

10

10

10

0.0256867

gKCa [mS/cm2]

10-0.1246

10-0.01853

10-0.6427

2.34906

gAMPA [µS]

0

0

0

0.513425

gNMDA [µS]

0

0

0

0.00434132

gGABA [µS]0

τCa [ms]

010

2.86870

102.87528

0.00252916

10

1.908778

121.403

Table S3. Initial values of the differential equations
The initial values of differential equations are described in the first table. 𝑉, 𝑛𝑘 ,and
[𝐶𝑎2+ ] were different in simulations and are described in the second, third, and fourth table,
respectively.

𝑉

*

ℎ𝑛𝑎

0.045

𝑛𝑘

*

ℎ𝐴

0.045

𝑚𝐾𝑆

0.34

𝑠𝑁𝑀𝐷𝐴

0.01

𝑥𝑁𝑀𝐷𝐴

0.01

𝑠𝐴𝑀𝑃𝐴

0.01

𝑠𝐺𝐴𝐵𝐴
[𝐶𝑎2+ ]

0.01*
𝑉

Figure1L-78.060990002692 mV

The other figures-45 mV𝑛𝑘

Figure1L

0.01099578591813208
2+2+

2G ([Ca ]=12.5 µM), S3E ([Ca ]=12.5 µM)

0.1

The other figures

0.54[𝐶𝑎2+ ]

Figure1L

9.050249227774513 µM

2G7, 8.5, 10.6, 12.5 µM

S3E

5, 25, 50, 70 µM

The other figures

1 µM

Table S4. Plotted area in each simulation
All differential equations were solved in t = [0, 10] except for t = [0, 30] and t = [0, 90]
in Fig. 3 C and D and Fig. S4 D and E. The plotted time ranges are described in thefollowing table.

FigureInitial time [sec]Final time [sec]

1L

0

10

2B (Cluster 2)

5.0

7.0

2C (Cluster 1 / upper)

6.0

7.1

2C (Cluster 1 / lower)

6.11

6.16

2C (Cluster 2 / upper)

6.0

6.4

2C (Cluster 2 / lower)

6.28

6.30

2F

5.05

5.60

5.0

5.1

2G ([Ca ]=12.5 µM)

00.1

2H (lower)

5.0

5.7

S2C (Cluster 1)

6.00

6.55

S2C (Cluster 2)

6.2

6.4

S3C (Cluster1)

5.1

5.8

S3C (Cluster2)

5.00

5.25

S3D

5.00

5.25

5.0

5.1

S3E ([Ca ]=70 µM)

0

0.1

S3F

5.0

5.3

S3G

5.1

5.7

S3H

5.0

5.3

S4C

5.0

5.7

The other figures5

10

2G ([Ca2+]=7, 8.5, 10.6 µM)
2+S3E ([Ca2+]=5, 25, 50 µM)
2+

Table S5. Classification criteria of Resting, SWS, Awake, and Slow-wave with few
spikes
The classification criteria of Resting, SWS, Awake, and Slow-wave with few spikes is
described in the table.

Resting

Spike number per second< 2 or peak frequency=0.2Hz

SWS

0.2<peak frequency<10.2 Hz
and spike number per second >= 5 × peak frequency -1

Awake

peak frequency>=10.2Hz

Slow-wave withfew spikes

0.2<peak frequency<10.2 Hz
and spike number per second < 5 × peak frequency -1

Table S6. Eigenvalues of fixed points
The eigenvalues of the fixed points are as follows. When more than one fixed points exist,
their eigenvalues are listed in an ascending order of the membrane potential.

FigureEigenvalues
2+{−0.11 + 2.52𝑖, −0.11 − 2.52𝑖}

2+

2G ([Ca ]=8.5 µM){−1.86, −0.02}, {−1.46, 0.04}, {−0.11 + 2.51𝑖, −0.11 − 2.51𝑖}

2G ([Ca2+]=10.6 µM){−2.55, −0.04}, {−1.25, 0.19}, {−0.10 + 2.48𝑖, −0.10 − 2.48𝑖}

2G ([Ca ]=7 µM)

2+

2G ([Ca ]=12.5 µM){−3.08, −0.06}, {−1.15, 0.36}, {−0.09 + 2.45𝑖, −0.09 − 2.45𝑖}

S3E ([Ca2+]=5 µM){0.30 + 1.97𝑖, 0.30 − 1.97𝑖}

S3E ([Ca2+]=25 µM){−3.25, −0.10}, {−1.07, 0.67}, {0.22 + 1.96𝑖, 0.22 − 1.96𝑖}

2+

S3E ([Ca ]=50 µM){−3.95, −0.22}, {−0.93, 1.47}, {0.09 + 1.91𝑖, 0.09 − 1.91𝑖}

S3E ([Ca2+]=70 µM){−4.01, −0.24}, {−0.91, 1.59}, {0.07 + 1.90𝑖, 0.07 − 1.90𝑖}

Table S7. The number of the randomly-generated parameter sets and the parameter
sets classified into ‘True SWS’
The number of the randomly-generated parameter sets and the parameter sets classified intothe ‘True SWS’ and the hit rates are listed in the table below. In order to avoid theunanticipated errors during simulations, the kernels were shut down and again run at hourlyintervals in parameter searches. Since it was not clear whether the calculation with the lastparameter of each hour in each kernel was solved or not, the number of all
randomly-generated parameters was not strictly accurate and had errors estimated within
100,000 parameters (5.3×10-2% of total parameters).
*The number of all randomly-generated parameters** The number of parameters classified into the ‘True SWS’

Channels in the model

All*

SWS**

Hit rate

gCa, τCa, gKCa, gL, gNa, gNaP, gAMPA, gK, gKS, gA, gAR, gGABA(Model 1)

5690995279

0.00490%

gCa, τCa, gKCa, gNa, gNaP, gAMPA, gK, gKS, gA, gAR, gGABA

517477310.00002%

gCa, τCa, gKCa, gL, gK, gKS, gA, gAR, gGABA

535878200.00000%

gCa, τCa, gKCa, gL, gNaP, gAMPA, gK, gKS, gA, gAR, gGABA

507332859

0.00116%

gCa, τCa, gKCa, gL, gNa, gAMPA, gK, gKS, gA, gAR, gGABA

741496834

0.00046%

gCa, τCa, gKCa, gL, gNa, gNaP, gK, gKS, gA, gAR, gGABA

581747026

0.00045%

gCa, τCa, gKCa, gL, gNa, gNaP, gAMPA, gGABA

5036704116

0.00230%

gCa, τCa, gKCa, gL, gNa, gNaP, gAMPA, gKS, gA, gAR, gGABA

6874041160.00023%

gCa, τCa, gKCa, gL, gNa, gNaP, gAMPA, gK, gA, gAR, gGABA

5044511527

0.01045%

gCa, τCa, gKCa, gL, gNa, gNaP, gAMPA, gK, gKS, gAR, gGABA

5184372287

0.00554%

gCa, τCa, gKCa, gL, gNa, gNaP, gAMPA, gK, gKS, gA, gGABA

5120092905

0.01768%

gCa, τCa, gKCa, gL, gNa, gNaP, gAMPA, gK, gKS, gA, gAR

5056561340

0.00672%

gCa, τCa, gKCa, gL, gNaP (Model 2)502042200.00000%

gCa, τCa, gKCa, gL, gNa, gNaP

546315286

0.00157%

gCa, τCa, gKCa, gL, gNaP, gAMPA

5359992190.00035%

gCa, τCa, gKCa, gL, gNaP, gK (SAN model)

81487909

1008

0.00124%

gCa, τCa, gKCa, gL, gNaP, gKS

601917100.00000%

gCa, τCa, gKCa, gL, gNaP, gA

5088113200.00039%

gCa, τCa, gKCa, gL, gNaP, gAR

585482300.00000%

gCa, τCa, gKCa, gL, gNaP, gGABA

5208280137

0.00263%

Table S8. Oligonucleotide sequences for the T7-gRNA templates
gRNAs were produced by the T7-gRNA templates of Kcnk1, Kcnk2, Kcnk3, Kcnk4, Kcnk5,
Kcnk6, Kcnk7, Kcnk9 (set 1), Kcnk9 (set 2), Kcnk10, Kcnk12, Kcnk13, Kcnk15, Kcnk16, andKcnk18 (Fig. S5 A-O).

Gene

Kcnk1

gRNA

ForwardNo.

primer No.

1

1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2

2
3
1

Kcnk2

2
3
1

Kcnk3

2
3
1

Kcnk4

2
3
1

Kcnk5

2
3
1

Kcnk6

2
3
1

Kcnk7

2
3
1

Kcnk9
2

Sequence (5' to 3')
CACTATAGGACCGACACTCCATAATTGCGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGACCGACACTCCATAATTGCGCACTATAGGTCCCAATTCCAATTTCCCGGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGTCCCAATTCCAATTTCCCGG
CACTATAGGTTGTCCATGCCGTTCTGCTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGTTGTCCATGCCGTTCTGCTG
CACTATAGGTCTCCCCACGAACTGAAGGGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGTCTCCCCACGAACTGAAGGG
CACTATAGGTTATCACTCTGACGACCATGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGTTATCACTCTGACGACCATGCACTATAGGTGTGGTTCTGGATCCTCGTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGTGTGGTTCTGGATCCTCGTG
CACTATAGGACGCATCAACACCTTCGTGGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGACGCATCAACACCTTCGTGG
CACTATAGGAAACCGATGAGCACCATGTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGAAACCGATGAGCACCATGTG
CACTATAGGTTCATGACCATGAACGCCGGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGTTCATGACCATGAACGCCGGCACTATAGGAAATACCTCTGCCCGAGTTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGAAATACCTCTGCCCGAGTTG
CACTATAGGTCTGAGGACTCGTCGATGAGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGTCTGAGGACTCGTCGATGAGCACTATAGGAAAAGCCTTCCAGACCCCGGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGAAAAGCCTTCCAGACCCCGG
CACTATAGGCCGATACTTTGTGGAGCTTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGCCGATACTTTGTGGAGCTTG
CACTATAGGCTGAAGATGTTGACGTCCTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGCTGAAGATGTTGACGTCCTG
CACTATAGGGTTTGGTAACTATCCTTGGGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGGTTTGGTAACTATCCTTGGG
CACTATAGGGTCTGCAACACGAGCACCAGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGGTCTGCAACACGAGCACCAGCACTATAGGCTCATCTTGCTGCCCGATCGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGCTCATCTTGCTGCCCGATCGCACTATAGGGTGCTGATGTAGATCTGTCGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGGTGCTGATGTAGATCTGTCGCACTATAGGAATGGCTACCCAGTCACCTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGAATGGCTACCCAGTCACCTG
CACTATAGGCCCAGCCATTTACCACCTTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGCCCAGCCATTTACCACCTTG
CACTATAGGAAGTCAGCTCGATCGTGTTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGAAGTCAGCTCGATCGTGTTG
CACTATAGGGTGAAGGTACAGGCGATCAGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGGTGAAGGTACAGGCGATCAG
CACTATAGGGTGTCCAGTGGAAATTCGCGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGGTGTCCAGTGGAAATTCGCG

3
1
Kcnk9
2
(set 2)
3
1

Kcnk102
31

Kcnk122
3
1

Kcnk132
31

Kcnk152
3
1

Kcnk162
3
1

Kcnk18

2
3

1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2

CACTATAGGGAACATAACCAGCGTCAGGGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGGAACATAACCAGCGTCAGGG
CACTATAGGGGAACCGGCGAATTTCCACGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGGGAACCGGCGAATTTCCACG
CACTATAGGTGTATGTTCTACGCTGTGCGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGTGTATGTTCTACGCTGTGCGCACTATAGGAATACCGATGAGGAGCTACGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGAATACCGATGAGGAGCTACGCACTATAGGGTGGTTGCCATCTTCGTGGGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGGTGGTTGCCATCTTCGTGGG
CACTATAGGTTCTTCTGACTGCTCTCGAGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGTTCTTCTGACTGCTCTCGAGCACTATAGGTATCCGAAGTATGGAGCGCGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGTATCCGAAGTATGGAGCGCGCACTATAGGTCACCTATGGTCGACACCAGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGTCACCTATGGTCGACACCAGCACTATAGGCTGTTATCTCCGCAGGTTTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGCTGTTATCTCCGCAGGTTTG
CACTATAGGGCTGGGCATCATGAACAACGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGGCTGGGCATCATGAACAACGCACTATAGGGTACAGTGGTTTCTACCATGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGGTACAGTGGTTTCTACCATGCACTATAGGAAAATCTTCCCTCCCGTTGGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGAAAATCTTCCCTCCCGTTGG
CACTATAGGACCTGTTCCGGATGTTACCGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGACCTGTTCCGGATGTTACCGCACTATAGGGGAGAAGCGGTACTTTCTGGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGGGAGAAGCGGTACTTTCTGG
CACTATAGGTCCTTCTACTTCGCCATCAGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGTCCTTCTACTTCGCCATCAG
CACTATAGGCTCTCACCACCATAGGCTTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGCTCTCACCACCATAGGCTTG
CACTATAGGATGGGGATCCCACTCAATGGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGATGGGGATCCCACTCAATGG
CACTATAGGCGCGAAGTAGAAGCCCTCAGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGCGCGAAGTAGAAGCCCTCAGCACTATAGGGCATTACATCGCTGTGTATGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGGCATTACATCGCTGTGTATGCACTATAGGATGTATGCAAAGACGACCAGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGATGTATGCAAAGACGACCAGCACTATAGGACTGACCACCATCGGGTTTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGACTGACCACCATCGGGTTTG
CACTATAGGAAGGGAAGTTTCGCTACCTGTTTTAGAGCTAGAAATAGC
GGGCCTAATACGACTCACTATAGGAAGGGAAGTTTCGCTACCTG

Table S9. Primer sequences used in quantitative PCR (qPCR)
1st qPCR, the No. from one to three, and 2nd qPCR, the No. from four to six, primers (Fasmac
Co., Ltd.) were used in genotyping of knockout mice Kcnk1, Kcnk2, Kcnk3, Kcnk4, Kcnk5,
Kcnk6, Kcnk7, Kcnk9 (set 1), Kcnk9 (set 2), Kcnk10, Kcnk12, Kcnk13, Kcnk15, Kcnk16, andKcnk18.

Gene

No.

primerSequence (5' to 3')

1

Forward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForwardTGCTGAGCACCGACACTCCATAATTGC
AGCTGCCTTATGAGGACCTGCTGCG
AGGTGAAGTCCCAATTCCAATTTCCCGCAGCAATTATGGAGTGTCGGTGCTCAG
GTGGCCATTGTCCATGCCGTTCTGC
GCTGGGATGAAGAAGAAGCAGGAAACG
CCGCGTGCTGGAGGCCAGCAATTA
CCGAGGTGAAGTCCCAATTCCAATTTCCC
TCGGTGCTCAGCAACGCCTCGGGAAGCGAAGAAGAGCGCCGAGGTGAAGTCCAGGAAACGGTGACAAATCCCAGCAGACTTCTCCAAGCAGGTGGTGGCCATTGTGGAAACATCTCCCCACGAACTGAAGGGGAATTCCCAGCAAGGCATAGATGATG
ATTTTGTGGTTATCACTCTGACGACCAGAAGCCTTACCTGCCACGTAGTCTC
CCTGTGGTGTGGTTCTGGATCCTCG
CCAATCATGCTCAGAACAGCTGCAAAGAGCCTGGGTGAACGCATCAACACCTTCGTG
TGCTGATGCACGACACGAAACCGATGAGCA
ATGCACGACACGAAACCGATGAGCACCATG
CCGCTCACACTAGTCATGTTCCAGAGCCTG
GGTGCTGCGATTCATGACCATGAACGCCGCGACAGCCTGGCCGTTGTGCGTGAG
CACGTAGAAATACCTCTGCCCGAGTTCCTAGCCTACTTCGCCTCAGTGCTCACTGCGTGTCTGAGGACTCGTCGATGACAGAGAAGGTTGAGACTCCGTCCCC
CATCCAAAAAGCCTTCCAGACCCCG
GGTCCAGAGATCCTGCCCCTACACCTGGTTGCGAGCAGTGTCCCGCCGAACTCG
CCGGGGATAAGGGCACAACACAAGAGACATGCA
GGATTACCCCAGTGAGAATCTGGCCTTCATCG
TCTGGAAGGCTTTTTGGATGGGTTGGGTCG
TCTCGGAGTCGCCCAGGACCCCGGGTGAGAATCTGGCCTTCATCGACGAGTCGCCCTATACCGATACTTTGTGGAGCTTTG
CTCACCTTCCAGTTGACAAAGAGGGAC
GGATAGGAAGCTGAAGATGTTGACGTCCGAAGGCCCTACAGATGGCTGGAAGTAC
ACACTTCAGAGGTTTGGTAACTATCCTTGG
CAGCTTGGAAGAGGTATCTCAGACTCCCAGCCCCAGGTAGATCCAAAGCTCCAGTGCCAACTACCACGCCCTATACGCTGGAAGTACAGCATCCAAGGACG
GCCTTCTTCCCAATCTGTTTGATGAGG
GCTGGGGCACAGGCCCCCAAGGA

2
3
Kcnk1
4
5
6
1
Kcnk2

2
3
1

Kcnk3

2
3
1
2
3

Kcnk4
4
5
6
1
2

Kcnk5

3
4
5
6

1
2
3
Kcnk6
4
5
6
1
Kcnk7

2
3
1

Kcnk9
2
(set 1)3
1
2
Kcnk9(set 2)3
4
5
6
1

Kcnk102
31

Kcnk12

2
3

ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseTGGCATTTTGGAAGATGAGCGGGTGGTAGCGGAAAGTCTGCAACACGAGCACCA
CCCATGTCTACAGCGTACCTCTTCCTCTCACAGAACTCATCTTGCTGCCCGATCTACCTGGGGATGGAGGCGTAGTCAGGCAGAGAGGTGCTGATGTAGATCTGTCCAGCCAGGATGAGGATGATCAGGTG
TCCTCGGCCTGGTTGCCATGGTGC
ATCGGGCAGCAAGATGAGTTCTGTGAGTGAGGCTGGCAGGGTCCGGATCGG
GTGTTGCAGACTTTCCGCCGTGTTTCTGGCCGTTCTGGATGCCCGGACAGACTACCTGGGGATGGAGGCGTAGTCAGC
CCAGCGAATGGCTACCCAGTCACC
CTTCTGTGTGGTCTATGCAGCCCTTG
CTGCACCCAGCCATTTACCACCTTG
CTCTCTACCCTCCTCTTTTCTGGACTCGACAGGAAGTCAGCTCGATCGTGTTG
CTTCCTGCTATGTTCTTCAAGCACTCC
GCAGGTAGGTGAAGGTACAGGCGATCAACTTCCCTTTAGCATCTCCTTCTTCGC
CGCTGGTGTCCAGTGGAAATTCGC
ACACCCATTCCTTACCGATAGTTGTGAGGCTCTGGAACATAACCAGCGTCAGG
AAGGCCTTCTGTATGTTCTACGCTGTGCGAAGTAGAAGGAACCGGCGAATTTCCACCTTAAAGCCGAAGAAGTCCGTCTCAGAGCAAGGCCTTCTGTATGTTCTACGCTGTGCTGGAACATAACCAGCGTCAGGGGGATACCCTGACCATGAATACCGATGAGGAGCTAC
GGATTTCCAGCAAGAATCTCCGCAACTTCCCACCGCGCTGGTGTCCAGTGGAAGATGACAGTGATGGCGAAGTAGAAGGAACCG
CCAGCGTCAGGGGGATACCCAGCACAG
ACCCCACAGAACCTTTCCTAACTATTGTTCCCAATCTCCGCAACTTCTCCCTCCAGTAGCTCTCAATCTTGTGGTCCTGCGATTCCTGACGGAAAACAGTGGTTGCCATCTTCGTGG
TTCTGACTGCTCTCGAAGGGCTGTTCCGATTGTGTTCTTCTGACTGCTCTCGACTACCTCGTCACTGGTGGCCTTGTCAGCCACTATCCGAAGTATGGAGCGC
CTCAGGGGATAGCATGTCCAGTGAGTG
GCCGCTCACCTATGGTCGACACCA
CCACCGTGTTCTCTGCGCTCGAAAG
CCTGGCTGTTATCTCCGCAGGTTTCAAAGAGGTTGAAGAACAAGATGGTCCC
GCGTGGGCGCGCTGGGCATCATGAACAAC
CCGCTCTACCTGGAGGCGCTGGTCTCTG

1
Kcnk132
31

Kcnk152
3
1

Kcnk162
31
2
3

Kcnk184
5
6

Forward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseForward
ReverseACTTCGTGGGTACAGTGGTTTCTACCATA
CCAGGAACTAGGGGCATGGTGTGACGAAAAATCTTCCCTCCCGTTGTGG
TGGTTTTGTCTTCTTCAGGGTTTGGGAGATGTTGCACCTGTTCCGGATGTTACCGTGCTTCCCACCATGCCAAAGAGGAC
AGTCGTCGGCGGAGAAGCGGTACTTTCTG
GCTGTTTTCGATGCGCTGGAATCCG
CTTCGCAGGCTCCTTCTACTTCGCCATCA
GCTGGCAAGCCCTTTCCAAAGTCAAAGTTCATCACTCTCACCACCATAGGCTTCCTTCCTCTGCAGCGCCTCATCTCTCATGCTCTGATGGGGATCCCACTCAATGGAGGACGATGCCTTACCTGGGAATGGTGATAAACGCGAAGTAGAAGCCCTCA
CTGGTCATTCTCATCTTCCCACCCATG
CCAGCAAGCATTACATCGCTGTGTATCTAAGGTCCAGGCCTCTGATGAGCTG
CGCAGGAGATGTATGCAAAGACGACCAATTCTCAGCAATCTGGATGAAGTGGGCGTGACACTGACCACCATCGGGTTTG
ACGATGATGTAGATGGAGAAGAAGAGGATTGCCAAAGGGAAGTTTCGCTACCTTG
CTCACATCTGTAGGTCACGGAAATGC
CCTCTCCCCGTCATCGCCCTGGTCTCGAAGCCTAGCTCGGTCTCCCAGAAGGCACCAGCACGATGTCCCCAAACCC
ACCGAGCTAGGCTTCGAGGATGCTTTCCGGAACCCTCTCCATTTACCAAGGTAG
TCCTGCACACCTACAAAACCCTCATG

References
1.
2.
3.
4.

5.
6.

7.
8.
9.

10.

11.
12.

Tatsuki F, et al. (2016) Involvement of Ca(2+)-Dependent Hyperpolarization in Sleep
Duration in Mammals. Neuron 90(1):70-85.
Silverman BW (1986) Density estimation for statistics and data analysis(Chapman&Hall, London).
Sunagawa GA, et al. (2016) Mammalian Reverse Genetics without Crossing Reveals
Nr3a as a Short-Sleeper Gene. Cell reports 14(3):662-677.
Ma M, Ye AY, Zheng W, & Kong L (2013) A guide RNA sequence design platform
for the CRISPR/Cas9 system for model organism genomes. BioMed researchinternational 2013:270805.
Ran FA, et al. (2013) Genome engineering using the CRISPR-Cas9 system. Nature
protocols 8(11):2281-2308.
Naito Y, Hino K, Bono H, & Ui-Tei K (2015) CRISPRdirect: software for designing
CRISPR/Cas guide RNA with reduced off-target sites. Bioinformatics
31(7):1120-1123.
Wang H, et al. (2013) One-step generation of mice carrying mutations in multiple
genes by CRISPR/Cas-mediated genome engineering. Cell 153(4):910-918.
Cho SW, Kim S, Kim JM, & Kim JS (2013) Targeted genome engineering in humancells with the Cas9 RNA-guided endonuclease. Nature biotechnology 31(3):230-232.
Tsujino K, et al. (2013) Establishment of TSH beta real-time monitoring system inmammalian photoperiodism. Genes to cells : devoted to molecular & cellularmechanisms 18(7):575-588.
Sunagawa GA, Sei H, Shimba S, Urade Y, & Ueda HR (2013) FASTER: an
unsupervised fully automated sleep staging method for mice. Genes to cells : devoted
to molecular & cellular mechanisms 18(6):502-518.
Rhead B, et al. (2010) The UCSC Genome Browser database: update 2010. Nucleic
Acids Res 38(Database issue):D613-619.
Lu B, et al. (2007) The Neuronal Channel NALCN Contributes Resting Sodium
Permeability and Is Required for Normal Respiratory Rhythm. Cell 129(2):371-383.