FASTER: an unsupervised fully automated sleep staging
method for mice
Genshiro A. Sunagawa1, Hiroyoshi S
ei2, Shigeki Shimba3, Yoshihiro Urade4 andHiroki R. Ueda1,5,6,7,8*
1

Laboratory for Systems Biology, RIKEN Center for Developmental Biology, Kobe 650-0047, Japan
Department of Integrative Physiology, Institute of Health Biosciences, The University of Tokushima Graduate School,
Tokushima 770-8503, Japan3
Department of Health Science, College of Pharmacy, Nihon University, Funabashi 274-8555, Japan
4
Department of Molecular Behavioral Biology, Osaka Bioscience Institute, Suita 565-0874, Japan
5
Laboratory for Synthetic Biology, RIKEN Quantitative Biology Center, Kobe 650-0047, Japan
6
Department of Biological Sciences, Graduate School of Science, Osaka University, Toyonaka 560-0043, Japan
7
Graduate School of Frontier Biosciences, Osaka University, Suita 565-0871, Japan
8
Department of Systems Pharmacology, Graduate School of Medicine, University of Tokyo, Tokyo 113-0033, Japan
2

Identifying the stages of sleep, or sleep staging, is an unavoidable step in sleep research andtypically requires visual inspection of electroencephalography (EEG) and electromyography(EMG) data. Currently, scoring is slow, biased and prone to error by humans and thus is themost important bottleneck for large-scale sleep research in animals. We have developed an
unsupervised, fully automated sleep staging method for mice that allows less subjective andhigh-throughput evaluation of sleep. Fully Automated Sleep sTaging method via EEG/EMG
Recordings (FASTER) is based on nonparametric density estimation clustering of comprehensive EEG/EMG power spectra. FASTER can accurately identify sleep patterns in mice thathave been perturbed by drugs or by genetic modification of a clock gene. The overall accuracy is over 90% in every group. 24-h data are staged by a laptop computer in 10 min, which
is faster than an experienced human rater. Dramatically improving the sleep staging process inboth quality and throughput FASTER will open the door to quantitative and comprehensive
animal sleep research.

Introduction
Ever since the discovery of sleep/wake status relationship to electroencephalography (EEG) in the early
20th century, sleep staging based on EEG has been
the standard method to evaluate sleep/wake status inanimals. In mice, states of consciousness are classified
into at least three stages: nonrapid eye movement
sleep (NREM sleep), rapid eye movement sleep(REM sleep) and wake. Determination of an animal’s
sleep/wake stage is based on visual inspection of EEG
and electromyography (EMG) of the animal by welltrained human raters with some or no computational
support. The classical sleep scoring criteria for miceCommunicated by: Ryoichiro Kageyama*Correspondence: uedah-tky@umin.ac.jp

502

Genes to Cells (2013) 18, 502–518

use the amplitude of selected frequency bands of
EEG and EMG (Fig. 1A). NREM sleep is characterized by large and slow EEG waves with low EMG
amplitude, REM sleep has lower and faster EEG withvery low EMG amplitude, and wake has an EEG pattern similar to REM sleep with very high amplitudeof EMG.
There are two problems scoring sleep/wake stages
by visual inspection: quality control and throughput.
Scoring sleep stages depends on human rater bias andwill differ between two human raters (inter-rater variance) and a human rater analyzing the same data multiple times (intrarater variance). The main reason why
visual inspection by a human rater is still used for sleep
staging is related to the fundamental difficulty of sleep
staging in drawing boundaries between different
stages. There are many factors causing this difficulty,

DOI: 10.1111/gtc.12053© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Fully automated sleep staging method
(A)

(B)
NREM
Frequency domain data
50

0
–0.5

0

2

4

6

8

0
–0.5
0.5
0

PC01
–3.5

0
50

EMG

0

8

16

30

0
–0.5

2

4

6

8

PC02

Voltage (mV)

0

0

2

4

Time (s)

0

–5

–10

–10

–10

–15

–15
–20

5

0

20

6

8

–20

5
4

–5
2

0

0

20

–40

–20

PC01

1

–20

60

120

Time (m)

–15
–40

0

–15
–40

0

0

–5

–5

–10

–0.5

–35
0 50 0 50

Frequency (Hz)

5

0

10
–0.5
0.5

–0.1

10

(3) Annotaion

0

–35
35

5

PC01

Wake

0.5

–0.1
0.1

10

–40

Time (s)

–35
35

5

3

20

PC01

Log10 (EEG Delta power (mV2))

0

–0.1
0.1

10

PC02

PC02

Voltage (mV)

–0.5
0.5

120

Principal
components
35

Nonparametric density estimation clustering(2) Clustering

0

90

0.1

Eigenvectors
EEG EMG

PC02

REM

60

Time (m)Time (s)
0.5

PC03

–3.5
0

24

PC02

–1.5

0

–0.5

Time (s)

–1.0

0

20

–40

NREM

4

10

0.5
0.0

2

5

1

REM

Wake

–2

–1

0

–0.5

0

20

PC01

3

–3

–20

PC01

PC02

–0.5
0.5

Frequency (Hz)0.5

Principal component analysis

EEGTime domain data

Normalized power (A.U.)(1) Character extraction0

Voltage (mV)

Voltage (mV)

0.5

REM
Wake

5
0
–5
–10

–1.0
1

Log10 (EMG power (mV2))

NREM

–15
–40

–20

0

20

PC01

Figure 1 (A) Representative EEG and EMG time-domain data during 8 s of NREM sleep, REM sleep and wake status. NREM
sleep is characterized by high-amplitude EEG delta waves (0.5–4 Hz) with low EMG power, REM sleep is characterized by lowamplitude, high-frequency EEG theta waves (6–10 Hz) with very low EMG power, and wake is characterized by high and varying
EMG power. (B) An overview of FASTER, the unsupervised fully automated sleep staging. (1) Each epoch’s power spectrum ofboth the EEG and EMG is calculated by FFT. EEG/EMG power spectra are transformed into values of principal components byprinciple component analysis (PCA), and the values of top four principle components are used for further analysis. (2) Epochs are
clustered by the nonparametric density estimation clustering using the principal components. This clustering method estimates theprobability density of the dataset and defines a cluster as a high-density area that is connected by Delaunay triangulation. The
number of the clusters is automatically calculated by the method. (3) The clusters are annotated according to the median logarithmof EMG power and EEG delta power of each cluster (middle panel). Histograms around the middle panel represent the relative
counts of each sleep/wake stage by the logarithm of EMG power (upper histogram in the middle panel) and the logarithm of
EEG delta power (left histogram in the middle panel). In this representative case, cluster 1 is first annotated as ‘wake’ by the median logarithm of EMG power. Then, the remaining clusters are annotated as ‘NREM sleep’ and ‘REM sleep’ by the median logarithm of EEG delta power. The final staging is shown in the right panel.

such as the differences between individual animals, thenoise of EEG/EMG recording device and variance insurgical techniques used in electrode implantation thatcan make the definition of rigid boundaries problematic. Visual inspection aims to compensate for thesesources of variance in the EEG/EMG data, sacrificing
objectivity in favor of the human brain’s pattern-finding abilities. The other problem with sleep staging
based on visual inspection is low throughput. Visualinspection of long-term EEG/EMG data is time-consuming and tedious. Even an experienced human rater
requires hours to score 24 h of mouse sleep data. This

makes it difficult to carry out quantitative and comprehensive sleep research in animals and is hindering
sleep biology from becoming a truly data-driven field.
To overcome these problems, many types of sleep
staging programs have been developed since the1960s (Drane et al. 1969; Larsen & Walter 1970;
Smith & Karacan 1971; Martin et al. 1972; Kohn
et al. 1974; Winson 1976). These methods can bedivided into two processes: character extraction andsleep staging.
Character extraction identifies features from EEG/
EMG to discriminate sleep stages. Along with the© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Genes to Cells (2013) 18, 502–518

503

GA Sunagawa et al.

development of automated sleep staging programs,
efforts have been made to identify the proper features for extraction. Many automated staging programs use specific bands of the EEG power
spectrum as features (Johns et al. 1977; Chouvet
et al. 1980; Van Luijtelaar & Coenen 1984; Stanus
et al. 1987; Mamelak et al. 1988; Benington et al.
1994; Berthomier et al. 2007; Kohtoh et al. 2008;
Brankack et al. 2010; Gilmour et al. 2010; Vural &Yildiz 2010; Pan et al. 2012). The most common
power spectra bands in rodents are referred to as
delta (0.5–4 Hz), theta (6–10 Hz) and sigma (10–15 Hz). Other characters have also been proposed,
including information that is harder for humans tointuitively interpret. These include coefficients ofwavelet analysis (Ebrahimi et al. 2008; Gabran et al.
2008; Sinha 2008; Fraiwan et al. 2010, 2011; Garces
Correa & Laciar Leber 2010; Agrawal et al. 2011),
bispectral density (Acharya et al. 2010; Swarnkar
et al. 2010), parameters of multichannel autoregressive modeling (Zhovna & Shallom 2008) and matching pursuit method with slow wave patterns (Picot
et al. 2011). Recently, two groups have proposedusing a broader range of the EEG spectrum, more
finely binned than the three classical bands, andreported higher-quality sleep stage classification(Vivaldi & Bassi 2006; Rytk€
onen et al. 2011). No
matter what kind of staging methods are used, thenumber of features passed to the sleep staging step isusually important for performance. Some groups
have achieved reduction of feature dimension without losing critical information using principal component analysis (Vivaldi & Bassi 2006; Gilmour et al.
2010) or support vector machine-based recursive feature elimination (Koley & Dey 2012).
The features abstracted from the raw dataset arethen used for scoring sleep stages. Automated staging
programs can be divided into two groups according
to their annotation methods: unsupervised andsupervised. Unsupervised staging programs use ‘hard’rules to annotate the data. Such rules are defined
before the staging process. Therefore, unsupervised
staging programs do not use rules based on specificdata derived from each subject. This makes theresults reproducible, but sensitive to data outside theanticipated input range. However, supervised staging
programs, which have recently become a popular
approach, use ‘soft’ rules for annotation. They learn
the proper decision rules based on a training datasetof each subject and then analyze the remaining dataof the subject. This makes the staging program resilient to irregular data because the training dataset
504

Genes to Cells (2013) 18, 502–518

includes uncommon characters from the originaldata. However, a human rater should annotate thetraining dataset for each subject, which introduces
subjectivity.
Classically, automated sleep staging programs werebased on unsupervised algorithms (Larsen & Walter
1970; Smith & Karacan 1971; Martin et al. 1972;
Kohn et al. 1974; Winson 1976; Johns et al. 1977;
Mendelson et al. 1980; Van Luijtelaar & Coenen
1984; Kuwahara et al. 1988; Neckelmann et al. 1994;
Doman et al. 1995; Veasey et al. 2000; Kohtoh et al.
2008). Since the late 1980s, supervised sleep staging
programs incorporating sophisticated machine learning algorithms have become more popular, due totheir ability to improve the accuracy of staging by
including the variances between subject animals. The
machine learning algorithms include neural networks(Mamelak et al. 1991; Robert et al. 1996, 1997;
Schaltenbrand et al. 1996; Baumgart-Schmitt et al.
1997; Gr€ozinger & R€
oschke 2002; Ebrahimi et al.
2008; Gabran et al. 2008; Hassaan & Morsy 2008;
Sinha 2008; Garces Correa & Laciar Leber 2010;
Fraiwan et al. 2011), support vector machines (Crisler
et al. 2008; Koley & Dey 2012), training hidden
Markov models (Grube et al. 2002; Pan et al. 2012)and parameter estimation of Gaussian mixture models(Acharya et al. 2010).
Although many of the supervised sleep staging
programs have improved the throughput of sleep
staging because only small amount of manual staging
is required, subjectivity still remains as a problem.
These methods let the machine ‘learn’ the desired
sleep stage recognition rules from the dataset of
EEG/EMG tagged with sleep/wake status scored by‘human’ experts. Sleep staging by a human rater orby supervised programs can absorb many kinds ofvariance in the EEG/EMG data although it sacrifices
the objectivity because of inter- or inner-rater differences. Therefore, an automated sleep staging programthat can fully substitute manual scoring requires minimization of subjectiveness due to inter- and innerrater differences and maximization of robustness
against many sources of variability inherited in EEG/
EMG data. The past automated sleep staging programs have either of four major problems in this
aspect.
The first problem is related to subjectiveness.
Modern automated sleep staging programs are oftentied with machine learning algorithms using ‘soft’rules, which requires supervision by a human rater.
They are very robust against variance, but subjective
in nature. To solve this subjectiveness problem, we© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Fully automated sleep staging method

adopted the classical unsupervised approach using‘hard’ rules, which is objective in principle.
The second problem is related to robustness, oneof the weaknesses in unsupervised algorithms. Mostof the automated sleep staging programs, either supervised or unsupervised, attempt to stage each ‘epoch’of EEG/EMG data, which has the minimal time
length, typically 4–10 s, for staging sleep. It is notvery difficult to stage an epoch that has typical waveform for a specific stage. The toughest cases are whenthe epoch has borderline characteristics of multiple
stages because the borderline between multiple stages
is very variable and therefore hard to determine. To
solve this ‘borderline’ problem, we adopted the strategy introduced by Gilmour et al. (Gilmour et al.
2010). In this strategy, epochs are first clustered intoa limited number of groups and then annotated basedon the statistical value of each group. This ‘clusterfirst’ strategy divides the borderline problem into twoparts, clustering and annotation, and makes the annotation much easier and more robust because the borderline cases are already clustered into a group thathas similar characters if the clustering algorithm isappropriate. Therefore, this ‘cluster-first’ strategy canconvert the difficult classification problem in staging
of borderline cases into the clustering problem. This
conversion is beneficial because we can use recently
developed sophisticated algorithms for clustering.
The third problem is also related to robustness,
especially involved in clustering. Many automatedsleep staging programs adopt model-based algorithms
in clustering or classification. Model-based algorithms
work well when the model is known a priori. Whenthe system steps out from the model, characters thatare extracted through the false model are useless. To
solve this problem, we adopted a model-free clustering algorithm which is based on nonparametric density estimation (Azzalini & Torelli 2007).
The remaining forth problem is involved in character extraction and again related to robustness of
sleep staging. Many automated sleep staging programs use limited power bands of EEG/EMG for
sleep staging. They may cause problems becausemouse inbred strains are known to have different
distributions of EEG power bands (Franken et al.
1998). As mentioned above, some groups reported
that the broader range of the EEG spectrum withfiner bin than the three classical band contains useful
information for sleep staging (Vivaldi & Bassi 2006;
Rytk€
onen et al. 2011). Thus, we used comprehensive EEG/EMG power spectra with the finest bin
the recorder may take. We expected that a wider

range and finer bin of EEG/EMG power spectra
might cover the individual variance of the subject.
The use of broader and finer EEG/EMG power
spectra increases the dimension of the data, which
consumes the computation power and takes longer
time for sleep staging. To reduce the dimension ofdata, we therefore performed principle component
analysis of EEG/EMG power spectra in the character
extraction step.
In this study, we combined solutions to all of thefour problems and developed an unsupervised fully
automated sleep staging program, FASTER (Fully
Automated Sleep sTaging method via EEG/EMG
Recordings). All efforts were made to minimize subjectivity and increase robustness against many sourcesof variance accompanying to the EEG/EMG data.
FASTER was first developed by using the EEG/
EMG data from wild-type mice and then tested formice with drug-induced alterations in sleep/wakepatterns or arrhythmic genetic modifications in aclock gene. The source code of FASTER is open
source (GNU General Public License) and available
for future improvement in the sleep research field.

Results
Basic structure of the FASTER algorithm

FASTER is composed of three major steps, character
extraction, clustering and annotation (Fig. 1B). In thecharacter extraction step (Fig. 1B, top panels), both
EEG and EMG time-domain data are first split into
epochs with constant time length of 8 s. This EEG/
EMG time-domain data are next converted by fast
Fourier transform (FFT) into frequency-domain data,
the EEG/EMG power spectra. We used the power
up to the maximum frequency, called Nyquist frequency (i.e., half of the sampling frequency), whichcorresponds to 50 Hz in this study because EEG/
EMG data were recorded at 100 Hz. Characters of
EEG/EMG power spectra are then extracted by principal component analysis (PCA). In the following
clustering step (Fig. 1B, middle panels), all epochs aregrouped together into a limited number of clusters
according to the characters of EEG/EMG power
spectra obtained in the character extraction step. We
used the nonparametric density estimation clustering
method (Azzalini & Torelli 2007) to cluster theepochs. This method uses two pieces of information
for clustering. First, it estimates a probability density
for each epoch by a Gaussian kernel method. Second,
all epochs are connected by Delaunay triangulation,© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Genes to Cells (2013) 18, 502–518

505

GA Sunagawa et al.

which is a method to triangulate a set of points (i.e.,
epochs in this study) under the criteria that the circumcircle does not include any other points. Oncethe information is calculated, the clusters are definedby finding subsets of high-density regions that areconnected by Delaunay triangulation. Because thenonparametric density estimation clustering method
itself selects the number of clusters, it is powerful tocluster certain datasets without a priori model or without information on the number of clusters. In thisstudy, the model-free clustering was proper because
we wanted to keep the clustering as objective as possible. In the final annotation step (Fig. 1B, bottom
panels), each cluster is annotated without supervision
by a human rater according to the statistical value ofcharacteristics for each cluster. In this study, we use
relative differences in the median logarithm of EMG
power and EEG delta power of each cluster. Overall,
all steps were fully automated and unsupervised.
Optimization of FASTER

To optimize the parameters for hard rules in the
FASTER algorithm, we first prepared an EEG/EMG
dataset with manual annotation of sleep/wake stages.
We recorded continuous EEG/EMG data for 6 days
from four C57BL/6J mice and manually scored into
three stages (NREM sleep, REM sleep and wake)according to the criteria described in the Experimental Procedures section. These annotated data werethen used to determine the optimal parameter values
for the FASTER algorithm by assessing its performance and computation time for different values ofparameters.
Parameter values for the FASTER algorithm were
optimized in three steps: character extraction, clustering and annotation. First, we optimized a parametervalue for character extraction. Because the optimal
degree of information extracted from EEG/EMG
power spectrum is not immediately known, we optimized the number of principal components in character extraction to achieve higher performance withless computation time for sleep staging. Second, we
optimized parameter values for clustering. The developers of nonparametric density estimation clustering
supply a recommended shrinkage factor for thesmoothing bandwidth and the optimal number ofgrids for the given dataset size (Azzalini & Torelli2007). However, the clustering results are heavily
dependent on the smoothing factor (bandwidth ofthe density estimation kernel) and the grid width ofcluster core search. Therefore, the parameter values
506

Genes to Cells (2013) 18, 502–518

for smoothing factor and the grid width of cluster
search were optimized by assessing both performanceand computation time. Finally, we optimized theparameter values for annotation. As discussed above,
many unsupervised staging programs use a certain
threshold to classify (or annotate) each epoch into
different sleep stages. However, FASTER annotateseach cluster, not an epoch, by using two statisticalvalues of each cluster: the median logarithm of EMG
power and EEG delta power. We thus calculated
optimal values for these parameters in annotation. In
the following sections, we describe the detailed procedures to determine the optimal parameter values
for FASTER.
Optimization of character extraction: the numberof principal components

The parameter optimized in the character extraction
step is the number of principle components passed tothe subsequent steps. In this study, the top four principal components are sufficient to express over 38.4%of the original EEG/EMG power spectrum’s variance(Fig. S1 in Supporting Information). However, theoptimal number of components to adopt in the subsequent steps is not immediately obvious. Therefore,
we compared the performance of staging results andcomputation time with different principal component
numbers (PCN). When PCN is in the range of 2–6,
the accuracy, sensitivity of NREM sleep and wakeand specificity of all three stages reach to nearplateau, especially when PCN is greater than 4.
Meanwhile, computation time exponentially increases
along the increase in PCN. However, the sensitivityof REM sleep showed a peak at PCN = 4. According to these results (Fig. 2A), PCN = 4 was chosen
as an optimal number of principal components.
Optimization of clustering: the smoothing factor
of density estimation

The clustering algorithm first estimates the probability
density from the EEG/EMG dataset. Therefore, wefirst optimized the smoothing factor, which is a critical parameter in probability density estimation. A
Gaussian kernel method is used for the nonparametricdensity estimation in the clustering step. Estimating
the probability density of a dataset is like rescaling
and smoothing the dataset. When the dataset x has N
points of d dimensional data and the Gaussian kernel
which has a bandwidth h, the estimated probability
density ^f ðyÞ is given by:© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Fully automated sleep staging method(A)

Number of principal components
NREM

1.25

1.00

REM

Wake

1.00

REM

Wake

4

5

6

0.9

1.0

1.00

1000

Specificity

Time (s)0.50

Sensitivity

2000

0.75

Accuracy

NREM

1.25

0.75
0.50

0.75
0.50

0.25
0.25
0.00

0
2

3

4

5

6

0.00
2

3

PCN

(B)

4

5

6

0.00
2

3

PCN

4

5

2

6

PCN

Smoothing factor of density estimation

0.25

40

1.25

1.25

1.00

1.00

Specificity

Time (s)0.50

Sensitivity

60

0.75

0.75
0.50

0.00

0
0.6

0.7

0.8

0.9

0.50
0.25
0.00

0.00
0.6

1.0

0.7

hmult

(C)

0.75

20
0.25

0.8

0.9

1.0

0.6

0.7

0.8

0.9

0.6

1.0

0.7

hmult

hmult

0.8

hmult

Grid number for cluster seek
1.00
300

1.25

1.25

1.00

1.00

200

Specificity

0.50

Sensitivity

Time (s)0.75

Accuracy

3

PCN

1.00

Accuracy

0.25

0.75
0.50

0.75
0.50

100

0.25

0.25
0.00

0
0

300

600

Ngrid

900

0.25

0.00
0

300

600

Ngrid

900

0.00
0

300

600

Ngrid

900

0

300

600

900

Ngrid

Figure 2 Optimization of parameter values in the character extraction and clustering steps. The panels show computation time,
accuracy, sensitivity and specificity for each stage (from left). In the sensitivity and specificity panel, red, green and blue dots
denote NREM sleep, REM sleep and wake, respectively. The points are mean and the shaded area denotes standard error of themean. Every optimization is carried out using 5400 epochs randomly sampled from the 6-day-length dataset of four C57BL/6J
mice. (A) Optimization of the number of principal components PCN in the character extraction step. While fixing hmult = 0.85
and Ngrid = 135 as arbitrary initial values, FASTER was tested in PCN from 2 to 6. Each PCN was evaluated four times for four
different mice, resulting in 16 tests. PCN = 4 was selected for the optimal value because nearly all of the performances (e.g., accuracy, sensitivity of REM sleep) peaked at PCN = 4. (B) Optimization of the smoothing factor hmult of probability density estimation in the clustering step. Fixing the Ngrid = 135 and PCN = 4 as arbitrary initial values, the performance of FASTER was tested
in hmult from 0.6 to 1.0 by 0.05. Each hmult was evaluated twelve times for four different mice, resulting in 48 tests. The computation time, accuracy and the specificity did not differ among different hmult. Because the sensitivity of REM sleep peaked at
hmult = 0.7 and the other two stages did not have change within this range, hmult = 0.7 was selected as an optimal value. (B) Optimization of the grids numbers Ngrid for cluster detection in the clustering step. While fixing hmult = 0.85 and PCN = 4 as arbitrary
initial values, FASTER was tested in Ngrid from 33 to 1080. Each Ngrid was evaluated twelve times for four different mice, resulting in 48 tests. Accuracy, sensitivity and specificity all reached near-plateau and keep increasing slightly when Ngrid is 540 orgreater. For an optimal number of grids, 540 was chosen for this data length. If there is enough computation power, picking Ngrid
for the same number of the dataset will be the best choice because theoretically it should provide the best results.

© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Genes to Cells (2013) 18, 502–518

507

GA Sunagawa et al.

^f ðyÞ ¼



N Y
d
X
yj  xi;j1;K
Nh1    hd i¼1 j¼1
hj

h ¼ ðh1 ;    ; hd ÞT ;
1

dþ4
4
hj ¼ hmult
rj ;ðd þ 2Þn
1
1 2
KðxÞ ¼ pﬃﬃﬃﬃﬃ e2x
2p
Azzalini and Torelli observed that shrinking hslightly toward zero is often advantageous and recommend the shrinkage factor hmult of 0.75 (Azzalini &Torelli 2007). Because the number of cluster that canbe detected is heavily dependent on smoothness ofthe estimated probability density, we optimized hmult
to improve the accuracy of the staging results(Fig. 2B). The smaller the bandwidth is, the shaggier
the estimated probability density would be. Especiallywhen hmult is smaller than 0.4, the number of theclusters, which are the number of peaks in the estimated probability density, did not match with thenumber of the detected saddle points (usually,
the number of peaks and saddles matches). Therefore,
the clustering program was not able to complete(Fig. S2 in Supporting Information). If the bandwidth
is too large, the cluster number will decrease anddetected clusters might not be sufficient to adequately
describe the dataset. We have chosen the optimal
hmult as 0.70 because the accuracy was constantly highthrough the range of 0.60 < hmult < 0.85 and thesensitivity of REM sleep was highest at hmult = 0.70,
whereas sensitivity of NREM sleep and wake did notchange much in this range (Fig. 2B).
Optimization of clustering: the grid numbers forcluster detection

After the calculation of estimated probability density
as described above, the clustering algorithm thensearches for cluster cores in the estimated probability
density. Therefore, the second optimized parameterin the clustering step is the number of grids Ngrid,
which is a critical parameter in the detection of cluster cores. Supposed the dataset has N points in total,
the clustering algorithm scans every N/Ngrid points inthe dataset from the point of maximum probability
density (Fig. S3 in Supporting Information). Theoretically, scanning every point in the dataset (i.e.,
Ngrid = N) will maximize the capability in detecting
cluster cores although it also costs the maximum
508

Genes to Cells (2013) 18, 502–518

computation power. Practically, sufficiently large
number of Ngrid grids will give acceptable results.
Therefore, we optimized Ngrid as smaller as possible
by maximizing performance of sleep staging alongwith minimizing the computation time (Fig. 2C).
The accuracy gets nearly constant when Ngrid is larger
than 135, as well as the sensitivities of NREM sleep
and wake. Meanwhile, REM sleep does not reach toplateau unless Ngrid is greater than 540. Considering
its performance and computation time, we chose 540
as the optimal Ngrid. Because the total points in thedataset were 5400 in this optimization, this grid number denotes that the probability density is scanned
through every 10 points. This is finer than the original clustering method, which uses default Ngrid of 50
for any size of the data which has more than 50 data
points (Azzalini et al. 2012).
Optimization of annotation

The annotation step in FASTER involves assigning
the proper sleep/wake stages to the limited numberof clusters. Sleep stages can be characterized by relative difference in EMG power and EEG delta power.
If each cluster is mainly made of the same stage, themedian logarithm of EMG power and EEG delta
power are sufficient to discriminate sleep/wake stages.
Therefore, we adopted the method to annotate theclusters by their statistical values (i.e., median logarithm of EMG power and EEG delta power). The
advantage of this approach is that we do not need tocare about the borderline epochs; at least the cluster
is guaranteed to have mostly epochs with the identical stage. If the median logarithm of EMG power ofa cluster is larger than the PEMG quartile of the total
EMG power, the cluster is annotated as ‘wake’. The
remaining clusters are annotated as ‘NREM sleep’ if
the median logarithm of EEG delta power of thecluster is greater than Pdelta quartile of the nonwake
epochs and as ‘REM sleep’ when it is not. We compared the accuracy of staging results for different
PEMG and Pdelta (Fig. 3A). For annotating ‘wake’ and‘nonwake’ status, the accuracy was maximized when
0.45 < PEMG < 0.70. For annotating ‘NREM sleep’and ‘REM sleep’, the accuracy was maximized when
0.05 < Pdelta < 0.45. In this study, we chose 0.5 for
PEMG and 0.1 for Pdelta because these values give themaximum accuracy when we used these values for
nonclustered epochs (Fig. 3B).
Overall, we optimized FASTER for the numberof principal components to pass to subsequent procedures, the smoothing factor for probability density© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Fully automated sleep staging method

(A)

1.25

1.00

NREM

REM

Wake

1.25

1.00

NREM

REM

Wake

0.50

0.75

1.00

0.50

0.90

0.75
0.50

0.75
0.50

0.25
0.5

0.25

0.90 0.50

0.0

0.00

0.25

0.00
0.25

0.10

Pdelta

0.50

0.75

0.00
0.25

PEMG

–0.5

–1.0

0.50

0.75

0.25

PEMG

1.00

PEMG

1.25

1.25

1.00

1.00

–2

–1

0

1

Log10 of EMG power (mV2)

0.50

Specificity–3

Sensitivity

0.75

Accuracy

Log10 of EEG Delta power (mV2)

0.10

0.50

Specificity

PEMG

Sensitivity

Accuracy

0.75

0.75
0.50

0.75
0.50

0.25
0.25
0.00

0.25

0.00
0.25

0.50

0.75

0.00
0.25

Pdelta

(B)

0.50

0.75

0.25

Pdelta

1.25

1.00

NREM

0.50

0.75

Pdelta

REM

Wake

1.25

1.00

NREM

REM

Wake

0.50

0.75

1.00

0.90

0.10

–0.5

0.25
0.00
0.50

0.75

0

1

0.50

0.75

0.25

PEMG

PEMG

1.25

1.25

1.00

1.00

0.75

Sensitivity–1

Accuracy

–2

0.00
0.25

PEMG
1.00

Log10 of EMG power (mV2)

0.50
0.25

0.00
0.25

–1.0
–3

0.50

0.75

0.25

Pdelta

0.0

0.75

0.50

Specificity

0.5

0.50

Specificity

Accuracy

0.50

0.90 0.50

Log10 of EEG Delta power (mV2)

0.10

Sensitivity

0.75

PEMG

0.75
0.50

0.75
0.50

0.25
0.25
0.00

0.25

0.00
0.25

0.50

Pdelta

0.75

0.00
0.25

0.50

0.75

Pdelta

0.25

0.50

0.75

PdeltaFigure 3 Optimization of parameter values in the annotation step. In the sensitivity and specificity panel, red, green and blue dots
denote NREM sleep, REM sleep and wake, respectively. The points and the shaded area in the three right columns denote meanand standard error of the mean, respectively. Every optimization is carried out using 5400 epochs randomly sampled from the
6-day-length dataset of four C57BL/6J mice. (A) The left panel shows a schematic view of the relationship between median logarithm of EMG power and EEG delta power and quartile of logarithm of EMG power (PEMG) and nonwake epochs’ EEG delta
power (Pdelta). Each dot represents different clusters. In this case, there are five clusters. The upper row of the right panels shows
results of the annotation when PEMG is tested from 0.05 to 0.95 by 0.05 and the Pdelta is fixed to 0.1. The lower row of the rightpanels shows results of annotation when Pdelta is tested from 0.05 to 0.95 by 0.05 and the PEMG is fixed to 0.5. Each PEMG and
Pdelta was evaluated twelve times for four different mice, resulting in 48 tests. The accuracy was maximized when Pdelta is from
0.45 to 0.70 and Pdelta is from 0.10 to 0.45, respectively. (B) Each dot in the left panel represents different epochs before clustering. The upper row of the right panels shows results of annotation when PEMG is tested from 0.05 to 0.95 by 0.05 and the Pdelta
is fixed to 0.1. The lower row of the right panels shows results of annotation when Pdelta is tested from 0.05 to 0.95 by 0.05 andthe PEMG is fixed to 0.5. Each PEMG and Pdelta had evaluated 12 times for four different mice, resulting in 48 tests. The accuracy
was maximized when PEMG is 0.5 and Pdelta is 0.10, respectively. We chose this number because this value gave optimal results inthe annotation of the clustered dataset (see (A)).

© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Genes to Cells (2013) 18, 502–518

509

GA Sunagawa et al.

(Table 1 and Table S1 in Supporting Information).
The average time spent to analyze a 24 h of mouse
data was 605.7  13.9 s using a laptop computerwith a single processor. One of the features in wildtype animals is that they exhibit their circadian sleep/
wake rhythm, even in constant darkness. They have arobust internal time, which is called circadian time inthis study. Figure 4E shows that FASTER was able
to detect the NREM sleep time difference in thesubjective day and night under constant darkness condition. The high accuracy suggests combining in-parallel optimization worked well in FASTER.

estimation, the number of grids to seek the cluster
peaks and the thresholds for annotation. Because
these optimizations were carried out in parallel, every
optimized parameter was combined to evaluate themaximal functionality of the FASTER algorithm, andthe analysis was redone in C57BL/6J male mice thatwere used for optimization. These animals wererecorded in a 12-h light/dark condition for 3 days,
followed by a constant darkness for 3 days (Fig. 4A).
The first PC has larger vector elements in the EMG
than in the EEG (upper row of Fig. 4B and Fig. S4A
in Supporting Information). The second, third andforth PCs have relatively smaller elements in EMG
and unique distribution of elements in EEG (thelower three rows of Fig. 4B and Fig. S4A in Supporting Information). Using these principal components, FASTER divided the epochs into three sleep/
wake stages (Fig. 4C, D and Fig. S4A in Supporting
Information) and the accuracy was 94.6  1.4%, thesensitivity and specificity for each stage were
96.2  1.7% and 95.6  3.1%, 87.1  3.3% and97.8  1.1%, 94.2  3.3% and 99.2  0.3% foreach NREM, REM and wake status, respectively

2

3

4

5

6

7

8

9

Observation

Start EEG/EMG
recordings

End EEG/EMG
recordings

EEG

(C)

EMG

0.1
0
–0.1

0.1
0
–0.1
0.1
0
–0.1

0
–10

Wake
NREM

–20
50 0

50

Frequency (Hz)Manual staging

FASTER method

–20

0

REM

10
0
–10

Wake
NREM

–20

20

–20

PC01

0

20

PC01

(E)

ManualFASTER

0-12 12-24

0-12 12-24

1500

NREM time (m)60Hourly mean (m)FASTER

REM

10

0.1
0
–0.1

0

(D)

ManualPC02

1

Because our final goal is to evaluate animals withunnatural sleep/wake status by FASTER, we tested
our method on animals with drug-induced sleep/
wake status. In other words, this examines the robustness of FASTER when the sleep/wake ratio haschanged temporally in a certain period of the experiment. First, we examined prolonged wakefulness by

PC02

Day

PC04 PC03 PC02 PC01

(B)
Change of light condition

Normalized power (A.U.)(A)

Staging drug-induced overwaking andoversleeping mice with FASTER

0 NREM
20
REM
0
60

Wake

1000

500

0

0

3

4

5

6

Day

7

8

Subjective time

Figure 4 Results of staging C57BL/6J mice with FASTER. (B), (C) and (D) show a representative mouse. Results for the otherthree animals are found in Supporting Information (Fig. S3A). In (A) and (D), orange and black boxes denote light on and off,
respectively. (A) Experimental scheme. (B) Eigenvectors of the top four principal components. (C) Scatter plot of the first two
principal components with staging results. The left and right panels show the results from manual staging and FASTER, respectively. (D) Time spent in each stage per hour during experiments. The dashed and solid lines show the manual staging and thestaging by FASTER, respectively. (E) The total NREM time during subjective day and night during 3 days of constant darkness.
Both in manual and in FASTER staging, mean NREM sleep time in subjective day and night was significantly different(P = 0.0012 and P = 0.0012 for manual staging and FASTER by the unpaired Student’s t-test).

510

Genes to Cells (2013) 18, 502–518© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

94.6  1.4%
91.9  2.5%
93.2  1.0%
92.6  1.1%
99.2  0.3%
98.1  0.8%
98.4  1.5%
98.5  0.7%
94.2  3.3%
92.2  2.2%
92.8  0.8%
94.3  0.8%
97.8  1.1%
98.7  0.2%
98.2  0.2%
98.4  2.2%
95.6  3.1%
90.9  4.0%
93.7  0.6%
92.4  2.6%
Basal
MOD-IP
DIP-IP
Basal
C57BL/6J
C57BL/6J
C57BL/6J
Bmal1/

4
3
3
3

96.2  1.7%
96.6  1.5%
96.1  1.8%
95.1  4.3%

87.1  3.3%
60.5  13.3%
71.7  6.1%
63.4  17.6%

AccuracySpecificity
Sensitivity
Specificity
Experiment
StrainNumberof mice

Sensitivity

Sensitivity

Specificity

Wake
REM
NREM

Table 1 Summary of the FASTER performance tests. Basal: spontaneous sleep/wake cycle under 12:12-h light/dark cycle for the first 3 days followed by constantdarkness for 3 days; MOD-IP: modafinil IP experiment to induce excessive wakefulness; DIP-IP: diphenhydramine IP experiment to induce excessive sleepiness. Eachvalue represents the mean  SD.

Fully automated sleep staging method

injecting modafinil into the intraperitoneal cavity to
induce wakefulness through dopaminergic transporter
inhibition (Qu et al. 2008). During EEG/EMG
recording, we injected modafinil intraperitoneally to
mice at the beginning of the light phase (Fig. 5A)
when mice have tendency to sleep. After manuallystaging the data (8-day-length, n = 3), we compared
the staging results derived from FASTER (Fig. S4B
in Supporting Information). The increase in wakeinduced by modafinil injection was detected both inmanual staging and in FASTER (Fig. 5B). The total
accuracy was 91.9  2.5% and the sensitivity andspecificity for each stage were 96.6  1.5% and90.9  4.0%, 60.5  13.3% and 98.7  0.2%,
92.2  2.2% and 98.1  0.8% for each NREM,
REM and wake status, respectively (mean  SD, 8day-length, n = 3, Table 1 and Table S1 in Supporting Information). Next, we examined prolonged
sleepiness by using diphenhydramine, which is a histamine H1 receptor antagonist. Diphenhydraminecauses sleep by blocking the wakefulness maintenancepathway of histamine (Saitou et al. 1999). Under
EEG/EMG recording, we injected diphenhydramine
intraperitoneally to mice at the beginning of the darkphase, which is the active phase for nocturnal animals(Fig. 5C). After manual staging (8-day-length, n = 3),
the results were compared with FASTER-derived
stages (Fig. S4C in Supporting Information). The
increase in NREM sleep time after diphenhydramine
injection was clearly detected by both manual staging
and FASTER (Fig. 5D). The accuracy for these micewas 93.2  1.0% and the sensitivity and specificity
for each stage were 96.1  1.8% and 93.7  0.6%,
71.7  6.1% and 98.2  0.2%, 92.8  0.8% and98.4  1.5% for each NREM, REM and wakestatus, respectively (mean  SD, 8-day-length, n = 3,
Table 1 and Table S1 in Supporting Information).
These experiments show that FASTER is capable ofdetecting not only normally distributed sleep/
wake status but also temporarily altered sleep/wakestatus.
The EEG/EMG recording data of six C57BL/6J
mice obtained in these overwaking or oversleeping
drug administration experiments can be also used forthe validation of FASTER in wild-type mice at normal condition. One-day-length of EEG/EMG
recordings during the previous day of drug administration that was recorded for control were evaluatedfor this purpose. Because each mouse had two
chances of administration (vehicle vs. overwaking
drug, or vehicle vs. oversleeping drug), total lengthof 2 days were used for the test per animal. In these© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Genes to Cells (2013) 18, 502–518

511

GA Sunagawa et al.

(A)

(B)

ManualFASTER

150

IP-1
DAY

1

2

IP-2

3

4

5

6

Observation-1

7

8

9

10

Observation-2

Start EEG/EMG
recordings

End EEG/EMG
recordings

Wake time (m)Light conditions
100

50

0

VEH

MOD

VEH

MOD

Drug types(C)

Manual(D) 300

FASTER

IP-1
DAY

1

2

3

NREM time (m)Light conditions
IP-2
4

5

6

Observation-1

7

8

9

10

Observation-2

Start EEG/EMG
recordings

200

100

End EEG/EMG
recordings
0

VEH

DIP

VEH

DIP

Drug types

(E)

(F)

Day

1

2

3

4

5

6

7

8

NREM time (m)Change of light condition

9

Observation
Start EEG/EMG
recordings

End EEG/EMG
recordings

ManualFASTER

0-12 12-24

0-12 12-24

1500

1000

500

0

Subjective time

Figure 5 Results of performance test of FASTER in various conditions. (A) Experimental scheme of staging modafinil induced
prolonged wakefulness in C57BL/6J mice with FASTER. (B) The total wake time for three hours after the IP. Both in manualand in FASTER staging, it was possible to detect significantly longer mean wake time in the modafinil-administered group(P = 0.016 and P = 0.018 for manual staging and FASTER by the unpaired Student’s t-test.). (C) Experimental scheme of staging
diphenhydramine-induced prolonged sleepiness in C57BL/6J mice with FASTER. (D) The total NREM time for six hours after
the IP. Both in manual and in FASTER staging, it was possible to detect longer NREM time in the diphenhydramine-administered group (P = 0.016 and P = 0.0064 for manual staging and FASTER by the unpaired Student’s t-test.). (E) Experimental
scheme of staging genetically modified circadian mutant Bmal1/ mice with FASTER. (F) The total NREM time during subjective day and night during 3 days under constant darkness. Both in manual and in FASTER staging, it was unable to detect significant difference in NREM time between subjective day and night (P = 0.77 and P = 0.42 for manual staging and FASTER by the
unpaired Student’s t-test.).

mice, the total accuracy was 92.5  1.3% and thesensitivity and specificity for each stage were95.4  1.2% and 93.6  2.6%, 66.5  10.1% and98.5  0.5%, 93.2  1.2% and 98.6  1.0% foreach NREM, REM and wake status, respectively(mean  SD, 2-day-length, n = 6). Taken together,
512

Genes to Cells (2013) 18, 502–518

these results confirm that FASTER is able to detect
sleep/wake stages accurately in wild-type animals, ineither normal or drug-administrated conditions, without human assistance. These results are much objective and the time spent for staging is faster incomparison with manual staging by human raters.

© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Fully automated sleep staging method

Staging genetically modified animals with
FASTEROne of the motivations behind FASTER is toincrease the throughput of sleep staging, for example,
to screen a large number of genetically modified micefor defects in sleep. As a proof of principle, we examined a genetically modified strain with severe circadian abnormalities – mice that lack circadian rhythmin their sleep/wake behavior. Testing FASTER in
these arrhythmic mice will evaluate the robustness of
this method against genetically modified mice withsleep/wake phenotypes. We used Bmal1/ mice,
which are known to exhibit an arrhythmic sleep/
wake pattern (Laposky et al. 2005). EEG/EMG fromthree male Bmal1/ mice was recorded while they
were kept in 12-h light/dark condition for the first
3 days followed by another 3 days under constantdarkness (Fig. 5E). Manual staging and results
obtained by FASTER were compared (Fig. S4D in
Supporting Information). FASTER showed no significant difference in NREM sleep time between subjective day and night in Bmal1/ mice during theconstant darkness condition (Fig. 5F), which is oneof the major characteristics found in mice lacking circadian clock. This arrhythmic phenotype is markedly
different from the circadian sleep/wake pattern inwild-type mice under the constant darkness condition(Fig. 4E). The total accuracy in Bmal1/ mice was
92.6  1.1%. Taken together, our results show that
FASTER can reliably detect unnatural sleep/wakedurations and distributions in both drug- and genetically perturbed animals, which are key factors forscreening mice with sleep/wake abnormalities.

Discussion
We have developed FASTER, an unsupervised fully
automated sleep staging method for mice. Thismethod reports sleep/wake stages by analyzing thecomprehensive power spectrum of EEG and EMG.
Although there are other semi-automated staging
programs, no system has achieved a fully automatedsleep staging functionality. FASTER is not only fully
automated, but the accuracy is comparable to thesemi-automated staging programs and it is much faster than manual staging by visual inspection.
The basic strategy of FASTER is to use the classical‘hard’ rule, which is objective but prone to variancewithin subjects. To adopt the classical hard rule, wehave deliberately divided the ‘classification’ into ‘clustering’ and ‘annotation’. The FASTER algorithm

thereby absorbs the individual variance as much as possible in the character extraction and clustering steps tosafely adopt the classical ‘hard’ rule in annotation step.
For example, in the character extraction step, FASTER uses comprehensive EEG/EMG power insteadof band-specific powers (e.g., EEG delta power),which is often used in many manual and automatedstagers. The use of comprehensive EEG/EMG power
could cover and hence absorb the individual varianceamong subjects. This ‘full power spectrum’ strategy isbased on the report on the interstrain variance in distribution of classical power bands (Franken et al.
1998), as well as other reports on that expanded bands
beyond the classical ones are informative for staging
accuracy (Vivaldi & Bassi 2006; Rytk€
onen et al. 2011).
In this aspect, two improvements might be possible in
the future to further cover and absorb the varianceamong subjects. First, using faster sampling to gather
higher powers as characters might detect intersubject
variance. In the current study, we have used 100-Hz
sampling for EEG/EMG recordings that limits theavailable power bands of both signals up to the
Nyquist frequency, 50 Hz. Because some eigenvectorsof the principal components include large elements in
the high end of the power spectrum (e.g., see the third
eigenvector in Fig. 4B), higher-frequency powers
might have useful information for sleep staging. Second, information on the transition of states might
absorb variance among subjects. In FASTER, theorder of the epochs does not affect the staging results.
When human raters stage the time series of EEG/
EMG, it is natural to take information on the previous
epochs (e.g., REM sleep tends to occur after NREM
sleep, but not after wake stage). Therefore, it will be
an important future work to take the order information into account in the character extraction step.
In the clustering step, we have adopted the clustering algorithm based on nonparametric estimation ofprobability density (Azzalini & Torelli 2007) to avoid
subjectiveness from holding up a model before clustering. This clustering method makes no assumptions
about the distribution of data and chooses the numberof clusters automatically. One of the few disadvantagesof this algorithm is that it is fairly computationallyexpensive. Related to this higher computational cost,
we need to point out two limitations in this study,
both of which can be solved in near future by the
improvement of computers. First, because we aimed ahandy system, we implemented FASTER in a free
software environment R (R Core Team 2012) andran on a laptop computer with a single processor. This
resulted in the maximum clustering size at once to be© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Genes to Cells (2013) 18, 502–518

513

GA Sunagawa et al.

5400 epochs practically. For example, if the total data
length is 86 400 epochs (8-day-length data when one
epoch is 8 s), the data are divided into 16 groups andclustered independently. In this way, 24 h of EEG/
EMG recorded from single mouse can be staged in
10 min. However, if 8-day-length data are clustereddirectly, approximately 6 h is expected for staging
1-day-length of data even ignoring other hardware
concerns (e.g., memory size). Ideally, all data should
be processed at once which will be practical in the
future with the increase in speed of computers. Second, the number of grid for scanning clusters was optimized maximizing the performance while minimizing
the computation time. Theoretically, the capability indetection of cluster cores is maximized when all of thedata points are evaluated through the cluster scan.
Therefore, the recommended number of grids of thisclustering might be the number of the data points inthe future.
In the annotation step, we adopted a classical ‘hard’rule-based selection, which is objective but sensitive
to variance among subjects. This employment waspossible in FASTER algorithm as the results of theabsorbance of the individual variance in the character
extraction and clustering steps as discussed above. This
allowed us to use very simple ‘hard’ rules to annotate
clusters with satisfactory performance.
The robustness of FASTER was tested in druginduced unusual sleep/wake conditions and genetically perturbed animals. The results of satisfactory
detection of drug-induced prolonged wakefulness andsleepiness show the ability to stage temporarily sleep/
wake differences. Because we have used circadianmutant Bmal1/ mice, which have arrhythmic
sleep/wake activity, the high performance of FASTER against these genetically modified animals not
only shows potential to evaluate mutants but also
implies the ability to detect unnatural sleep/wake distributions. In each experiment, FASTER was comparable to manual staging, which suggests the robustnessof the method and the potential of this method toanalyze animals with unknown sleep/wake status.
In comparison with the latest supervised sleep programs for rodents, FASTER has comparable accuracy
but lower sensitivity in REM sleep detection. For
example, using approximately 5% of human rater
scored data as a training dataset, Rytk€
onen et al. (2011)
achieved accuracy of 93% in rats and in mice, which iscomparable to our results. Their method, however,
was able to detect REM sleep in sensitivity of 89%,
which is higher than our optimized results (Table. 1
and Table S1 in Supporting Information). In this
514

Genes to Cells (2013) 18, 502–518

study, we used full EEG/EMG power spectra in thecharacter extraction step and then maximized not only
total accuracy but also sensitivity of REM sleep by
optimizing possible critical parameters in this scheme.
Therefore, we expected that, to further increase thesensitivity of REM sleep, it might be important to use
additional characters, which can efficiently distinguish
REM sleep from other stages. As mentioned above,the candidates for additional characters are higher-frequency power spectrum or order information of stages.
The improvement of REM sleep sensitivity is at the
top of the list for future works.
We believe at least two points need to be tested in
the future to further confirm the robustness of FASTER. First, we used diphenhydramine to induce prolonged sleep. However, this drug is not a major
choice for insomnia patients. To confirm the robustness of FASTER, other drugs for insomnia, that is,benzodiazepines, should be tested in the future. Second, we used Bmal1/ mice as a sleep disordermodel. In this strain, the sleep/wake distribution isvery different from wild-type mice because of theclear disorientation in their circadian rhythm. Therefore, other models characterized by disorganized or
fragmented sleep architecture, which has intact circadian rhythms, such as narcoleptic mice, should betested in the future to challenge the robustness of
FASTER in genetically modified mice.
In this study, we have developed FASTER, which
is an unsupervised fully automated sleep staging
method for mice based on comprehensive EEG/EMG
recordings. Full automation was achieved combining
the classical ‘hard’ rule-based annotation with themodern comprehensive character extraction alongwith model-free clustering. FASTER has comparable
accuracy to conventional visual inspection-based manual sleep staging method, and it is quicker than manual
inspection. All of the source codes of FASTER are
fully available because we believe this automated staging program has potential to free many sleep researchers from manual sleep staging labor and making thesource code freely available is a simple, but effective
way to dedicate to the field of animal sleep research.
Therefore, FASTER has the potential to enable datadriven quantitative and comprehensive sleep research.

Experimental procedures
Animal preparation and data collection
Four C57BL/6J mice (14 weeks old at recording) were used
for optimization of FASTER. To test the method, a total of 9

© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Fully automated sleep staging method
mice were used to represent different animal groups. These
mice include two groups of C57BL/6J mice that were administrated drugs (modafinil, diphenhydramine, 3 mice each, 11–13.4 weeks old at recording) by intraperitoneal injection andthree Bmal1/ mice (Shimba et al. 2011) (12–14 weeks old atrecording), which are circadian mutants. All C57BL/6J mice
were purchased from Oriental Yeast Co., Ltd. (Itabashi-ku,
Tokyo, Japan). See Table 1 for a summary of experiments. Procedures involving animals and their care were performed
according to the RIKEN Regulations for Animal Experiments(approval ID: AH18-01-19). We anesthetized animals andimplanted telemetry transmitters (Model F20-EET, DSI, St.
Paul, MN, USA) for simultaneous recording of two biopotentials (EEG and EMG). Two stainless-steel screws (1 mm diameter) were soldered to the wires of telemetry transmitters andinserted through the skull of the cortex (anteroposterior,+1.0 mm; right, +1.5 mm from bregma or lambda) and served
as EEG electrodes. Two other wires from the transmitter wereplaced into trapezius muscles serving as EMG electrodes. Animals were allowed at least 10 days to recover from surgery.
EEG/EMG data were recorded wirelessly with food and water
available ad libitum. The EEG/EMG data collecting system consisted of transmitters, an analog-digital converter and a recording
computer. Sampling rate was 100 Hz for both EEG and EMG.
Gold Acquisition (version 4.00, DSI) was used for the recording. The EEG/EMG data were converted to ASCII format, andboth the manual and the automated sleep staging were carriedout in the originally developed software on Ruby on Rails, ver.
3.1 (Rails Core Team 2011) and R (R Core Team 2012).

Optimization of FASTER
The animals were housed in an insulated soundproof recording
chamber maintained at an ambient temperature of 21 °C witha relative humidity of 50%. The chamber was light controlled
under 12-h light/12-h dark cycle (light on at 6:00 A.M.) for
the first 5 days followed by 4 days of constant darkness. The
recording was started on the second day. The data analysis was
carried out by the data recorded 6 days from the 6 A.M. ofthe third day (Fig. 4A).

Drug administration
Mice were housed in the same environment as in the recordings for FASTER optimization for 10 days, and EEG/EMG
was recorded for 8 days from 6:00 A.M. of the second day.
The light was controlled under 12-h light/12-h dark cycle(light on at 6:00 A.M.) through the experiment. In the prolonged wakefulness experiment (Fig. 5A), modafinil (M6940-50MG, Lot#029K4618, Sigma-Aldrich, St Louis, MO, USA)was dissolved to be 1% in sterile natural saline containing 10%
DMSO (06593-54, Lot#LZR7072, NACALAI TESQUE,
INC., Nakagyo-ku, Kyoto, Japan) and 2% cremophor EL(C5135-500G, 1439553-13509161, Sigma-Aldrich) immediately before use and administered intraperitoneally at 8:00
A.M. (2:00 in circadian time) on the experimental day at doseof 0.3 mL per mouse. The control group was administeredvehicle at a dose of 0.3 mL per mouse. Every mouse had two
opportunities for intraperitoneal injection during the experiment (2:00 in circadian time on the third or seventh day), andif the animal was injected modafinil on the first chance, vehicle was injected on the second and vice versa. In the prolonged sleepiness experiment (Fig. 5C), diphenhydramine(D3630-5G, Lot#040M0205V, Sigma-Aldrich) was dissolved
to be 0.2% to sterile water immediately before use and administered intraperitoneally at 7:00 P.M. (13:00 in circadian time)on the experimental day at a dose of 0.3 mL per mouse. The
control group was administered vehicle at a dose of 0.3 mL
per mice. As in the modafinil injection, all mice had two
opportunities for injection (13:00 in circadian time on thethird or seventh day), and if the animal was injected diphenhydramine on the first chance, vehicle was injected on thesecond and vice versa.

Genetically modified animals
Mice were housed in the same environment as in the recordings for FASTER optimization for nine days. The light control of the chamber was identical to the FASTER
optimization recordings as well (Fig. 5E).

Data analysis
We have used R (R Core Team 2012) on MacBook Air(1.8 GHz Intel Core i7, 4 GB 1333 MHz DDR3, Mac OS X
Lion 10.7.5) for every analysis. For both manual and automatedstaging, data were scored in 8-s epochs. In manual scoring, eachepoch was staged by visual inspection as NREM, REM andwake by the following criteria: NREM sleep was characterizedby high-amplitude EEG delta waves (0.5–4 Hz) with low
EMG power, REM sleep was characterized by low-amplitude,high-frequency EEG theta waves (6–10 Hz) with very low
EMG power, and wake was characterized by high and varying
EMG power. In the automated staging, EEG/EMG data weredivided into 8-s epochs and the power spectrum was computed
after detrended by subtracting estimated simple linear regressions followed by the production of a Hann window. BothEEG and EMG powers from the same epoch were connected
in tandem, resulting in combined EEG/EMG data. Because thesampling rate was 100 Hz, these 8 s of EEG/EMG combined
power spectra had 798 columns after removal of direct current
amplitudes (the first column in FFT results). If the logarithm tobase 10 of the sum of EEG or EMG power is either above 1 or
2, the epoch was annotated as ‘dirty’ data and ignored for automated staging. The logarithm to base 10 of every power for
every epoch was computed. By subtracting the mean and dividing by the standard deviation of these values, every power was
normalized within the epoch. Then, the principal component
analysis was performed on the normalized EEG/EMG combined power dataset to reduce the dimension. The top four
principal components were used for clustering giving two
parameters, the smoothing factor and the grid numbers, to the© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Genes to Cells (2013) 18, 502–518

515

GA Sunagawa et al.
pdfCluster library (Azzalini & Torelli 2007). The number ofprincipal components handed to the clustering library was set to
4 based on optimization. Furthermore, the parameters given to
the clustering library were optimized by simulation study. The
smoothing factor hmult was set to 0.7 and the number of grids
Ngrid was set to 540 when the dataset size is 5400. Because theclustering time increases in an exponential manner against thedata points, we decided to resample the data and divide thedataset into subsets which has less size than 5400 points. Whenthe original data are X = {x1,…, xM} and dividing X into ksubsets, the subset Xi that has N data points (N = 5400 in thisstudy) will be defined as Xi = {xi, xi+k,…, xi+(N-1)k} wherek = M/N and i = {1,…, k}. The clusters were annotated by thefollowing rules which were also optimized based on simulation.
First, the median logarithm of EMG power of each clusters wascomputed and clusters that have higher median logarithm of
EMG power than the optimized threshold 0.5 quartile of thatof all data points were annotated as ‘wake’. Within the remaining clusters, the median logarithm of EEG delta power (0.5–
4 Hz) was calculated, and if it was greater than the optimized
threshold 0.1 quartile of EEG delta power among all nonwake
data points, it was annotated as ‘NREM sleep’ and the rest of
the clusters were annotated as ‘REM sleep’.

Performance tests
The accuracy is defined by the ratio of epochs that haveagreement with manual and automated staging within the total
number of epochs. The sensitivity of stage X is defined by theratio of correctly staged X by automated staging within the
total number of X of manual staging. The specificity of stage
X is defined by the ratio of correctly staged non-X of automated staging within the total number of non-X of manual
staging.

Acknowledgements
We thank Yohei Koyama, Akinori Baba and Tatsuto Shibata
for helpful discussion. We also thank Masahiro Hayashi for
introducing the EEG/EMG recording technique to our laboratory. We thank Kenichi Inoue, Hiroshi Kiyonari and Kazuki
Nakao for providing us Bmal1/ mice with IVF technique.
We are grateful to the Laboratory for Animal Resources andGenetic Engineering for housing the mice. This work was
supported by a Grant-in-Aid for Scientific Research on Innovative Areas ‘Spying minority in biological phenomena(No.3306)’ (23115006) of MEXT, Japan, by Research Program of Innovative Cell Biology by Innovative Technology of
MEXT, Japan (to H. R. Ueda), by RIKEN Junior Research
Associate program for graduate students (to G. A. Sunagawa),
by Sky Orthopedic Clinic (to G. A. Sunagawa), by an intramural Grant-in-Aid from the RIKEN Quantitative Biology
Center (to H. R. Ueda), by an intramural Grant-in-Aid fromthe RIKEN Center for Developmental Biology (to H. R.
Ueda), by a Grant-in-Aid for Scientific Research (23590091)
JSPS, MEXT, Japan (to S. Shimba), and by The ‘‘High-Tech

516

Genes to Cells (2013) 18, 502–518

Research Center’’ Project for Private Universities, a matching
fund subsidy from the MEXT 2007 (to S. Shimba).

References
Acharya, U.R., Chua, E.C.-P., Chua, K.C., Min, L.C. &Tamura, T. (2010) Analysis and automatic identification of
sleep stages using higher order spectra. Int. J. Neural Syst.
20, 509–521.
Agrawal, G., Modarres, M., Zikov, T. & Bibian, S. (2011)
NREM sleep staging using WAV(CNS) index. J. Clin.
Monit. Comput. 25, 137–142.
Azzalini, A., Menardi, G. & Rosolin, T. (2012) R Package“pdfCluster”: Cluster Analysis via Nonparametric Density
Estimation (Version 1.0-0). http://cran.r-project.org/web/
packages/pdfCluster/index.html.
Azzalini, A. & Torelli, N. (2007) Clustering via nonparametricdensity estimation. Stat. Comput. 17, 71–80.
Baumgart-Schmitt, R., Herrmann, W.M., Eilers, R. & Bes, F.
(1997) On the use of neural network techniques to analyse
sleep EEG data. First communication: application of evolutionary and genetic algorithms to reduce the feature space andto develop classification rules. Neuropsychobiology 36, 194–210.
Benington, J.H., Kodali, S.K. & Heller, H.C. (1994) Scoringtransitions to REM sleep in rats based on the EEG phenomena of pre-REM sleep: an improved analysis of sleep
structure. Sleep 17, 28–36.
Berthomier, C., Drouot, X., Herman-Sto€ıca, M., Berthomier,P., Prado, J., Bokar-Thire, D., Benoit, O., Mattout, J. &d’Ortho, M.-P. (2007) Automatic analysis of single-channel
sleep EEG: validation in healthy individuals. Sleep 30,
1587–1595.
Brankack, J., Kukushka, V.I., Vyssotski, A.L. & Draguhn, A.
(2010) EEG gamma frequency and sleep-wake scoring in
mice: comparing two types of supervised classifiers. Brain
Res. 1322, 59–71.
Chouvet, G., Odet, P., Valatx, J.L. & Pujol, J.F. (1980) An
automatic sleep classifier for laboratory rodents. Waking
Sleeping 4, 9–31.
Crisler, S., Morrissey, M.J., Anch, A.M. & Barnett, D.W.
(2008) Sleep-stage scoring in the rat using a support vector
machine. J. Neurosci. Methods 168, 524–534.
Doman, J., Detka, C., Hoffman, T., Kesicki, D., Monahan,
J.P., Buysse, D.J., Reynolds, C.F., Coble, P.A., Matzzie, J.
& Kupfer, D.J. (1995) Automating the sleep laboratory:
implementation and validation of digital recording and analysis. Int. J. Biomed. Comput. 38, 277–290.
Drane, D.B., Martin, W.B. & Viglione, S.S. (1969) Patternrecognition applied to sleep state classification. Electroencephalogr. Clin. Neurophysiol. 26, 238.
Ebrahimi, F., Mikaeili, M., Estrada, E. & Nazeran, H. (2008)
Automatic sleep stage classification based on EEG signals byusing neural networks and wavelet packet coefficients. Conf.
Proc. IEEE Eng. Med. Biol. Soc. 2008, 1151–1154.
Fraiwan, L., Lweesy, K., Khasawneh, N., Fraiwan, M., Wenz,
H. & Dickhaus, H. (2010) Classification of sleep stages© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Fully automated sleep staging method
using multi-wavelet time frequency entropy and LDA.
Methods Inf. Med. 49, 230–237.
Fraiwan, L., Lweesy, K., Khasawneh, N., Fraiwan, M., Wenz,
H. & Dickhaus, H. (2011) Time frequency analysis forautomated sleep stage identification in fullterm and preterm
neonates. J. Med. Syst. 35, 693–702.
Franken, P., Malafosse, A. & Tafti, M. (1998) Genetic variation in EEG activity during sleep in inbred mice. Am. J.
Physiol. 275, R1127–R1137.
Gabran, S.R.I., Zhang, S., Salama, M.M.A., Mansour, R.R.
& George, C. (2008) Real-time automated neural-network
sleep classifier using single channel EEG recording fordetection of narcolepsy episodes. Conf. Proc. IEEE Eng.
Med. Biol. Soc. 2008, 1136–1139.
Garces Correa, A. & Laciar Leber, E. (2010) An automatic
detector of drowsiness based on spectral analysis and wavelet
decomposition of EEG records. Conf. Proc. IEEE Eng. Med.
Biol. Soc. 2010, 1405–1408.
Gilmour, T.P., Fang, J., Guan, Z. & Subramanian, T. (2010)
Manual rat sleep classification in principal component space.
Neurosci. Lett. 469, 97–101.
Gr€ozinger, M. & R€
oschke, J. (2002) The automatic recognition of REM sleep: a challenge and some answers. Methods
Find. Exp. Clin. Pharmacol. 24(Suppl D), 33–35.
Grube, G., Flexer, A. & Dorffner, G. (2002) Unsupervised
continuous sleep analysis. Methods Find. Exp. Clin. Pharmacol. 24(Suppl D), 51–56.
Hassaan, A.A. & Morsy, A.A. (2008) Adaptive hybrid system
for automatic sleep staging. Conf. Proc. IEEE Eng. Med. Biol.
Soc. 2008, 1631–1634.
Johns, T.G., Piper, D.C., James, G.W.W., Birtley, R.D.D. &Fischer, M. (1977) Automated analysis of sleep in the rat.
Electroencephalogr. Clin. Neurophysiol. 43, 103–105.
Kohn, M., Litchfield, D., Branchey, M. & Brebbia, D. (1974)
An automatic hybrid analyzer of sleep stages in the rat. Electroencephalogr. Clin. Neurophysiol. 37, 518–520.
Kohtoh, S., Taguchi, Y., Matsumoto, N., Wada, M., Huang,
Z.-L. & Urade, Y. (2008) Algorithm for sleep scoring inexperimental animals based on fast Fourier transform power
spectrum analysis of the electroencephalogram. Sleep Biol.
Rhythms 6, 163–171.
Koley, B. & Dey, D. (2012) An ensemble system for automatic sleep stage classification using single channel EEG signal. Comput. Biol. Med. 42, 1186–1195.
Kuwahara, H., Higashi, H., Mizuki, Y., Matsunari, S.,
Tanaka, M. & Inanaga, K. (1988) Automatic real-time analysis of human sleep stages by an interval histogram method.
Electroencephalogr. Clin. Neurophysiol. 70, 220–229.
Laposky, A., Easton, A., Dugovic, C., Walisser, J., Bradfield,
C. & Turek, F. (2005) Deletion of the mammalian circadian clock gene BMAL1/Mop3 alters baseline sleep architecture and the response to sleep deprivation. Sleep 28,
395–409.
Larsen, L.E. & Walter, D.O. (1970) On automatic methods ofsleep staging by EEG spectra. Electroencephalogr. Clin. Neurophysiol. 28, 459–467.

Mamelak, A., Quattrochi, J.J. & Hobson, J.A. (1988) A
microcomputer-based system for automated EEG collection
and scoring of behavioral state in cats. Brain Res. Bull. 21,
843–849.
Mamelak, A.N., Quattrochi, J.J. & Hobson, J.A. (1991) Automated staging of sleep in cats using neural networks. Electroencephalogr. Clin. Neurophysiol. 79, 52–61.
Martin, W.B., Johnson, L.C., Viglione, S.S., Naitoh, P.,
Joseph, R.D. & Moses, J.D. (1972) Pattern recognition of
EEG-EOG as a technique for all-night sleep stage scoring.
Electroencephalogr. Clin. Neurophysiol. 32, 417–427.
Mendelson, W.B., Vaughn, W.J., Walsh, M.J. & Wyatt, R.J.
(1980) A signal analysis approach to rat sleep scoring instrumentation. Waking Sleeping 4, 1–8.
Neckelmann, D., Olsen, O.E., Fagerland, S. & Ursin, R.
(1994) The reliability and functional validity of visual andsemiautomatic sleep/wake scoring in the Møll-Wistar rat.
Sleep 17, 120–131.
Pan, S.-T., Kuo, C.-E., Zeng, J.-H. & Liang, S.-F. (2012) A
transition-constrained discrete hidden Markov model for
automatic sleep staging. Biomed. Eng. Online 11, 52.
Picot, A., Whitmore, H. & Chapotot, F. (2011) Automateddetection of sleep EEG slow waves based on matching pursuit using a restricted dictionary. Conf. Proc. IEEE Eng.
Med. Biol. Soc. 2011, 4824–4827.
Qu, W.-M., Huang, Z.-L., Xu, X.-H., Matsumoto, N. & Urade,Y. (2008) Dopaminergic D1 and D2 receptors are essential forthe arousal effect of modafinil. J. Neurosci. 28, 8462–8469.
R Core Team (2012) R: A Language and Environment for Statistical Computing., Vienna, Austria: R Foundation for Statistical Computing.
Rails Core Team (2011) Ruby on Rails (ver. 3.1). http://rubyonrails.org.
Robert, C., Guilpin, C. & Limoge, A. (1997) Comparison
between conventional and neural network classifiers for rat
sleep-wake stage discrimination. Neuropsychobiology 35, 221–
225.
Robert, C., Karasinski, P., Natowicz, R. & Limoge, A.
(1996) Adult rat vigilance states discrimination by artificial
neural networks using a single EEG channel. Physiol. Behav.
59, 1051–1060.
Rytk€
onen, K.-M., Zitting, J. & Porkka-Heiskanen, T. (2011)
Automated sleep scoring in rats and mice using the naive
Bayes classifier. J. Neurosci. Methods 202, 60–64.
Saitou, K., Kaneko, Y., Sugimoto, Y., Chen, Z. & Kamei, C.
(1999) Slow wave sleep-inducing effects of first generation
H1-antagonists. Biol. Pharm. Bull. 22, 1079–1082.
Schaltenbrand, N., Lengelle, R., Toussaint, M., Luthringer,
R., Carelli, G., Jacqmin, A., Lainey, E., Muzet, A. &Macher, J.P. (1996) Sleep stage scoring using the neuralnetwork model: comparison between visual and automatic
analysis in normal subjects and patients. Sleep 19, 26–35.
Shimba, S., Ogawa, T., Hitosugi, S., Ichihashi, Y., Nakadaira,Y., Kobayashi, M., Tezuka, M., Kosuge, Y., Ishige, K., Ito,Y., Komiyama, K., Okamatsu-Ogura, Y., Kimura, K. &Saito, M. (2011) Deficient of a clock gene, brain and© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd

Genes to Cells (2013) 18, 502–518

517

GA Sunagawa et al.
muscle Arnt-like protein-1 (BMAL1), induces dyslipidemia
and ectopic fat formation. PLoS ONE 6, e25231.
Sinha, R.K. (2008) Artificial neural network and wavelet
based automated detection of sleep spindles, REM sleep andwake states. J. Med. Syst. 32, 291–299.
Smith, J.R. & Karacan, I. (1971) EEG sleep stage scoring by
an automatic hybrid system. Electroencephalogr. Clin. Neurophysiol. 31, 231–237.
Stanus, E., Lacroix, B., Kerkhofs, M. & Mendlewicz, J. (1987)
Automated sleep scoring: a comparative reliability study of twoalgorithms. Electroencephalogr. Clin. Neurophysiol. 66, 448–456.
Swarnkar, V., Abeyratne, U. & Hukins, C. (2010) Objective
measure of sleepiness and sleep latency via bispectrum analysis of EEG. Med. Biol. Eng. Comput. 48, 1203–1213.
Van Luijtelaar, E.L. & Coenen, A.M. (1984) An EEG averaging technique for automated sleep-wake stage identificationin the rat. Physiol. Behav. 33, 837–841.
Veasey, S.C., Valladares, O., Fenik, P., Kapfhamer, D., Sanford, L., Benington, J. & Bucan, M. (2000) An automatedsystem for recording and analysis of sleep in mice. Sleep 23,
1025–1040.
Vivaldi, E.A. & Bassi, A. (2006) Frequency domain analysis ofsleep EEG for visualization and automated state detection.
Conf. Proc. IEEE Eng. Med. Biol. Soc. 1, 3740–3743.
Vural, C. & Yildiz, M. (2010) Determination of sleep stage
separation ability of features extracted from EEG signals
using principle component analysis. J. Med. Syst. 34, 83–89.
Winson, J. (1976) A simple sleep stage detector for the rat.
Electroencephalogr. Clin. Neurophysiol. 41, 179–182.

518

Genes to Cells (2013) 18, 502–518

Zhovna, I. & Shallom, I.D. (2008) Automatic detection andclassification of sleep stages by multichannel EEG signal
modeling. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2008,
2665–2668.
Received: 7 February 2013
Accepted: 3 March 2013

Supporting Information
Additional Supporting Information may be found in the
online version of this article at the publisher’s web site:
Figure S1 Ratio of eigenvalues of each principal component
within the variance of the original signal.
Figure S2 Optimization results when smoothing factor of
density estimation hmult is selected from 0.2 to 1.4 by 0.2.
Figure S3 Schematic view of nonparametric density estimation clustering (Azzalini & Torelli 2007).
Figure S4 Individual results of FASTER, which were notshown in main figures.
Table S1 Results from all recordings in this study
Data S1 Source code of FASTER written in R (R CoreTeam 2012).

© 2013 The Authors
Genes to Cells © 2013 by the Molecular Biology Society of Japan and Wiley Publishing Asia Pty Ltd