Spindles and slow oscillations (SO)
Command  Description 

SPINDLES 
Waveletbased spindle detection 
SO 
Slow oscillation detection 
SPINDLES
Waveletbased sleep spindle detection
This command detects spindles using a waveletbased approach.
Optionally, it also detects slow oscillations (SO) and the temporal
coupling between spindles and SO. A target frequency is specified,
which corresponds to the center frequency of the Morlet wavelet. A
single SPINDLES
command can detect spindles at different
frequencies, and on different channels; there are also options for
collating putative spindles that are coincident in time but
observed on different channels and/or at neighbouring frequencies.
The SPINDLES
command has quite a large number of options (although
most of them are not necessary for its basic operation). For clarity,
this section is split into a number of subsections that describe
different aspects of this command:
 Basic usage
 Determining thresholds
 Collating spindles
 Additional spindle morphology metrics
 Generating spindle annotations and viewing spindles
 Spindle/SO coupling
After running the wavelet transform, spindles are detected based on the following rules:

the magnitude of the wavelet transform is smoothed (default sliding window of 0.1,
win
) and scaled by the wholesignal mean (or median, if themedian
flag is used) 
here, wholesignal means all unmasked epochs available to this command (i.e. may be only all N2 epochs rather than all epochs)

a spindle is called when a core interval of at least 0.3 seconds (
min0
) exhibits a signal at least 4.5 times (th
) the mean (or median) 
furthermore, all spindles must also exhibit an extended interval of at least 0.5 seconds (
min
) with a signal at least 2 times (th2
) the mean (or median) 
putative spindle intervals cannot be longer than 3 seconds (
max
); further, spindles within 0.5 seconds (merge
) are merged (unless the resulting spindle is greater thanmax
seconds) 
putative spindles are also filtered based on a quality metric q, the threshold of which is given by the
q
parameter
Following this heuristic procedure, in addition to the count and density (count per minute) of spindles, the following quantities are estimated to characterize the morphology of each spindle, after bandpassfiltering the signal for F_{C} ± 2 Hz:

spindle duration (of core plus flanking interval) in seconds, and number of oscillations

spindle amplitude, based on the maximum peaktopeak amplitude

mean and maximum of the transformed wavelet statistic

spindle frequency, based on counting zerocrossings, and secondarily based on the modal frequency from a FFT

spindle symmetry, based on the relative location of the spindle's central "peak" (the point of maximum peaktopeak amplitude), and a second "folded" symmetry index, calculated as 2S0.5

spindle chirp, based on the log ratio of mean peaktopeak time intervals in the first versus the second half of the spindle

a measure of integrated spindle activity (ISA), based on the sum of normalized wavelet coefficients in the spindle interval

a quality measure q based on the ratio of 1) the relative enrichment of power in the sigma band (for the spindle compared to baseline) versus 2) the comparable relative enrichment in power across all bands (i.e. again for the spindle compared to baseline)
Averaging over all spindles in a given run, Luna reports the mean spindle density, amplitude, duration, frequency, number of oscillations, chirp, symmetry indices, etc. Luna also reports perspindle statistics, as well as perepoch counts of spindles.
Basic usage
Parameters
The most basic parameter is fc
, which specifies the target frequency
(or frequencies) for the wavelet(s). Combined with the cycles
parameter, this defines the distribution of spindle frequencies that
are targeted by the wavelet. See also the
CWTDESIGN
command, which gives further information on
wavelet properties.
The primary output variables for an individual are spindle density
(DENS
), mean duration (DUR
) and frequency (FRQ
or FFT
, both
are typically very similar). An additional metric of interest is
integrated spindle activity (ISA), which, for a given spindle, is
the sum of the normalized wavelet coefficients in the spindle
interval. Higher values indicate more evidence in favour of a
spindle, and longer and higher amplitude spindles. Averaging across
spindles detected, for an individual, ISA_{S} is the
mean ISA, which reflects the average amplitude and duration of that
individual's spindles; ISA_{M} is the mean ISA per
minute, which reflects the combination of spindle amplitude, duration
and density; ISA_{T} is the total ISA, which
reflects the combination of spindle amplitude, duration and count.
Because it is normalized by the individual's baseline,
ISA_{S} may be a better metric on which to base
interindividual comparisons, compared to AMP
which is more likely
to reflect technical differences in recordings.
Primary options
Parameter  Example  Description 

sig 
sig=C3,F3 
Restrict analysis to these channels (otherwise, all channels are included) 
fc 
fc=11,15 
F_{C} value(s), e.g. to find 11 and 15 Hz spindles (default 13.5 Hz) 
cycles 
cycles=12 
Number of cycles, where more cycles gives better frequency but poorer temporal resolution (default 7) 
th 
th=6 
Multiplicative threshold for core spindle detection (default 4.5) 
th2 
th2=3 
Multiplicative threshold for noncore spindle detection (default=2) 
median 
median 
Flag to indicate that the median, not mean, is used for thresholding 
q 
q=0.3 
Quality metric criterion for individual spindles (default 0) 
Secondary parameters
Most users will not need to alter these.
Parameter  Example  Description 

fclower 
fclower=9 
Lower limit if iterating over multiple F_{C} values 
fcupper 
fcupper=15 
Upper limit if iterating over multiple F_{C} values 
fcstep 
fcstep=2 
Increment step if iterating over multiple F_{C} values 
thmax 
thmax=10 
Maximum threshold for spindle core (default: none) 
min 
min=1 
Minimum duration for an entire spindle (default 0.5 seconds) 
min0 
min0=0.3 
Minimum duration for a spindle core (default 0.3 seconds) 
max 
max=2 
Maximum duration for an entire spindle (default 3 seconds) 
win 
win=0.2 
Smoothing window for wavelet coefficients (default 0.1 seconds) 
local 
local=120 
Use local window (in seconds) to define baseline for spindle detection 
Outputs
Individuallevel output (strata: F
x CH
)
Variable  Description 

DENS 
Spindle density (count per minute) 
N 
Total number of spindles 
AMP 
Mean spindle amplitude (uV or mV units) 
DUR 
Mean spindle duration (core+flanking region) 
NOSC 
Mean number of oscillations per spindle 
FWHM 
Mean spindle FWHM (full width at half maximum) 
ISA_S 
Mean integrated spindle activity (ISA) per spindle 
ISA_M 
Mean integrated spindle activity (ISA) per minute 
ISA_T 
Total integrated spindle activity (ISA) 
FRQ 
Mean spindle frequency (from counting zerocrossings) 
FFT 
Mean spindle frequency (from FFT) 
CHIRP 
Mean chirp metric per spindle 
SYMM 
Mean spindle symmetry metric 
SYMM2 
Mean spindle foldedsymmetry metric 
Q 
Mean spindle quality metric 
DISPERSION 
Mean dispersion index of epoch spindle count 
DISPERSION_P 
Pvalue for test of overdispersion 
MINS 
Total duration of signal entered into the analysis (minutes) 
NE 
Number of epochs 
N01 
Number of spindles prior to merging 
N02 
Number of spindles post merging, prior to QC 
Epochlevel output (strata: E
x F
x CH
)
Variable  Description 

N 
Number of spindles observed in that epoch (for that target frequency/channel) 
Spindlelevel output (strata: SPINDLE
x F
x CH
)
Variable  Description 

AMP 
Spindle amplitude (uV or mV units) 
CHIRP 
Spindle chirp (1 to +1) 
DUR 
Spindle duration (seconds) 
FWHM 
Spindle FWHM (seconds) 
NOSC 
Number of oscillations 
FRQ 
Spindle frequency based on counting zerocrossings in bandpass filtered signal 
FFT 
Spindle frequency based on FFT 
ISA 
Integrated spindle activity 
MAXTSTAT 
Maximum wavelet statistic 
MEANSTAT 
Mean wavelet statistic 
Q 
Quality metric 
PASS 
Flag (0/1) for whether this spindle passes the quality metric criterion 
START 
Start position of the spindle (seconds elapsed since start of EDF) 
STOP 
Stop position of the spindle (seconds elapsed since start of EDF) 
START_SP 
Start position of the spindle (in sampleunits relative to current inmemory EDF) 
STOP_SP 
Stop position of the spindle (in sampleunits relative to the current inmemory EDF) 
SYMM 
Symmetry index (relative position of peak) 
SYMM2 
Folded symmetry index (0=symmetrical, 1=asymmetrical) 
Example
Here we estimate spindles for all NREM2 sleep for the
tutorial individual nsrr02
:
luna s.lst nsrr02 o out.db s "MASK ifnot=NREM2 & RE & SPINDLES fc=11,15 sig=EEG"
We see some output is sent to the console describing the process. For slow (i.e. target frequency of 11 Hz) spindles:
CMD #3: SPINDLES
detecting spindles around F_C 11Hz
wavelet with 7 cycles
smoothing window = 0.1s
detection thresholds (core, flank, max) = 4.5, 2x
core duration threshold (core, min, max) = 0.3, 0.5, 3s
pruned spindle count from 296 to 275
filtering at 9 to 13
filtering channel EEG_BP_11_4, bandpass FIR order 71
estimated spindle density is 1.37845
The relevant strata in out.db
are given by destrat
, namely
perindividual, perepoch and perspindle respectively:
destrat out.db
[SPINDLES] : F CH : 2 level(s) : AMP CHIRP DENS DISPERSION DISPERSION_P
: : : DUR FFT FRQ FWHM ISA_M ISA_S ISA_T
: : : MINS N N01 N02 NE NOSC Q SYMM SYMM2
: : :
: : :
[SPINDLES] : E F CH : (...) : N
: : :
[SPINDLES] : F CH SPINDLE : 692 level(s) : AMP CHIRP DUR FFT FRQ FWHM ISA
: : : MAXSTAT MEANSTAT NOSC PASS Q START
: : : START_SP STOP STOP_SP SYMM SYMM2
: : :
To get the spindle density and mean spindle duration and frequency for this individual:
destrat out.db +SPINDLES r F CH v DENS DUR FFT p 3
ID CH F DENS DUR FFT
nsrr02 EEG 11 1.378 0.847 11.280
nsrr02 EEG 15 2.090 0.847 13.475
We will further query the output using lunaR and
ldb()
. In R:
library(luna)
k < ldb( "out.db" )
We can use lx()
to see the strata (factor combinations) that
correspond to the destrat
output above:
lx(k)
MASK : EPOCH_MASK
RE : BL
SPINDLES : CH_F CH_E_F CH_F_SPINDLE
Looking at the perepoch data:
d < lx( k , "SPINDLES" , "CH_E_F" )
head(d)
> head(d)
ID CH E F N
1 nsrr02 EEG 92 11 2
2 nsrr02 EEG 92 15 0
3 nsrr02 EEG 93 11 4
4 nsrr02 EEG 93 15 0
5 nsrr02 EEG 98 11 2
6 nsrr02 EEG 98 15 1
Not that it is particularly informative in itself, but we can plot the number of spindles per epoch for fast and slow spindles:
par(mfcol=c(2,1), mar=c(4,4,1,1))
plot( d$E[ d$F == 11 ] , d$N[ d$F == 11 ] , pch= 20 , col="orange" ,
xlab="Epoch" , ylab="N spindles" , main="11 Hz spindles" )
plot( d$E[ d$F == 15 ] , d$N[ d$F == 15 ] , pch= 20 , col="skyblue" ,
xlab="Epoch" , ylab="N spindles" , main="15 Hz spindles" )
Luna provides a rough sanitycheck to test for unusual distributions
of spindles per epoch, as might occur if a movement artifact, for
example, causes multiple spindles to be scored in one epoch. To a
first approximation, we can treat the distribution of spindles per
epoch as Poissondistributed, and test for overdispersion (i.e. if the variance
is significantly greater than the mean). This is provided by the DISPERSION
and DISPERSION_P
variables in the individuallevel output:
destrat out.db +SPINDLES r F CH v DISPERSION DISPERSION_P p 3
ID CH F DISPERSION DISPERSION_P
nsrr02 EEG 11 1.143 0.606
nsrr02 EEG 15 1.243 0.038
In practice, we take values of DISPERSION
above 2.0 to reflect
potentially problematic recordings/patterns of spindle detection, so these
don't look too bad. Potentially the epochs with more than five spindles may
questionable (especially as we did not apply any type of artifact correction prior
to detecting spindles):
table( d$F , d$N )
0 1 2 3 4 5 6 7
11 206 131 50 7 3 1 1 0
15 154 139 62 27 14 2 0 1
Finally, here is an example spindle detected in this individual, shown alongside the bandpass filtered EEG for delta, theta, alpha and beta power as well as sigma power. As expected (i.e. by definition), the spindle corresponds to a brief increase in sigma power:
Note
Primarily in the context of lunaR, in the near future we'll be adding tutorials and helper functions to facilitate automatically making these types of plots quickly for all spindles in a dataset.
Thresholds
The core of the spindle detection algorithm involves imposing discrete thresholds on the quantitative metrics derived from the continuouslyvarying EEG. As such, there are a lot of magic numbers inherent in the spindle detection process. One of the key parameters is the threshold which the normalized wavelet coefficient must exceed for a putative spindle to be scored (the default of which is 4.5 times the mean). But is 4.5 a good magic number? What about 2 instead? Or what about 20? And so on...
To get some insight into whether this may be an appropriate choice for the particular individuals and analyses under consideration, Luna has a set of functions that evaluate different thresholds based on Otsu's method, a technique from computer vision that attempts to set a threshold so as to minimize the withinclass variance (and so correspondingly, maximize the between class variance) in the two classes that are formed by placing a dichotomizing threshold on a continuous measure. In our scenario, the measure is the wavelet coefficient, and we are trying to place a threshold that roughly corresponds to "spindle" and "nonspindle".
If the continuous metric exhibits clear bimodality, this approach should easily find a threshold that effectively splits the two peaks. In the context of spindles, where things are likely to be a bit messier, we consider the distribution of wavelet coefficients to be a mixture of two classes: the null distribution of noise, i.e. where there is no true spindle, and a second, less frequency class representing points where a spindle is present (and where we will expect higher values of the wavelet coefficient).
Parameters
Parameter  Example  Description 

empirical 
empirical 
Empirically determine thresholds 
setempirical 
setempirical 
Use empirically determined thresholds for spindle detection 
verboseempirical 
verboseempirical 
Output extensive information on threshold estimation 
Output
Individuallevel output (option: empirical
, strata: F
x CH
)
Variable  Description 

EMPTH 
Empiricallydetermined threshold 
EMPF 
Relative frequency of abovethresholds points based on EMPTH 
MEAN_OVER_MEDIAN 
Ratio of mean to median, to index skewness of the wavelet coefficients 
Betweenclass variance over range of thresholds (option: empirical
, strata: TH
x F
x CH
)
Variable  Description 

SIGMAB 
Betweenclass variance for given threshold 
Example
Here we use this approach on the three tutorial individuals: we run a
basic command to estimate spindles (for all NREM2 sleep, with no other
artifact detection in place), adding the empirical
threshold:
luna s.lst sig=EEG o out.db s "MASK ifnot=NREM2 & RE & SPINDLES fc=11,15 empirical"
In lunaR, we use ldb()
to load the resulting out.db
file:
k < ldb("out.db")
TH
):
d < lx( k , "SPINDLES" , "CH_F_TH" )
lid()
and basic R functions to create three plots separately, of the betweenclass variance as a function of threshold, split by fast and slow spindles. First, we get a list of the individual IDs:
ids < unique( d$ID )
par(mfcol=c(1,3))
for (id in ids) {
t < lid( d , id )
plot( t$TH[ t$F==11 ] , t$SIGMAB[ t$F == 11 ] ,
type="l" , lwd=2 , col="orange" ,
main=id, xlab="Threshold" , ylab="Betweenclass variance" , ylim=c(0,1))
lines( t$TH[ t$F==15 ] , t$SIGMAB[ t$F == 15 ] , type="l" , lwd=2 , col="skyblue" )
legend( 5,0.2,c("11 Hz","15 Hz"),fill=c("orange","skyblue"))
}
For the first two individuals, based on Otsu's criterion at least,
this suggests that thresholds around 4 or 5 are optimal, in the sense
of maximizing the between class variance (i.e. in wavelet coefficients
between "spindle" and "nonspindle" classes). For the third
individual, the results look strange, suggesting that extreme threshold
such as 20 times the mean is optimal. This is mirrored in the
individuallevel output for EMPTH
and EMPF
: the optimal thresholds
and the corresponding area under the curve:
destrat out.db +SPINDLES r F CH v EMPTH EMPF p 3
ID CH F EMPF EMPTH
nsrr01 EEG 11 0.034 4.000
nsrr01 EEG 15 0.067 3.250
nsrr02 EEG 11 0.041 4.000
nsrr02 EEG 15 0.024 6.250
nsrr03 EEG 11 0.001 20.000
nsrr03 EEG 15 0.002 20.000
This suggests there may be something off with the EEG channel for
nsrr03
. As well as looking at the raw data, one approach may be to
impose some automated artifact detection prior to spindle detection
and threshold estimation (i.e. as we would typically recommend).
Rerunning with aberrant epochs removed prior to spindle detection
(i.e. by including ARTIFACTS
and
SIGSTATS
commands:
luna s.lst sig=EEG o out2.db s "MASK ifnot=NREM2 & RE \
& ARTIFACTS mask & SIGSTATS mask th=3,3,3 & RE \
& SPINDLES fc=11,15 empirical"
This appears to normalize the "optimal" thresholds somewhat, in that the area under the curve (i.e. roughly corresponding to the portion of time in "spindle" versus "nonspindle" states) are ballpark what we may expect to observe in a typical individual who has around 2 or so spindles per minute detected at a central scalp EEG channel, with each spindle lasting approximately 1 second (i.e. 2/60 = 0.033):
destrat out2.db +SPINDLES r F CH v EMPTH EMPF p 3
ID CH F EMPF EMPTH
nsrr01 EEG 11 0.046 3.500
nsrr01 EEG 15 0.092 2.750
nsrr02 EEG 11 0.041 4.000
nsrr02 EEG 15 0.029 5.500
nsrr03 EEG 11 0.033 4.000
nsrr03 EEG 15 0.066 3.000
Returning to the plots, things look a bit better:
What to make of all this? We do not advise selecting different "optimal" thresholds for each individual based on these analyses. Rather, we suggest that looking at these (and other) metrics can be useful in basic sanitychecking for whether a given analysis for a given individual is likely to yield meaningful results.
Quality metrics
Luna calculates a simple QC metric for each spindle, by considering
the relative increase of nonspindle activity within the spindle
interval, relative to spindleactivity. We use five fixed bands,
which unlike the PSD
command, are not altered by changing special
variables. The three
nonspindle bands are: delta (0.54 Hz), theta (48 Hz), and
beta (defined here as 2030 Hz). The two spindle bands are slow
sigma (1013.5 Hz) and fast sigma (13.516 Hz).
For each channel, Luna calculates the baseline power for each band, P^{band}. For the interval spanning the i^{th} spindle, we calculate band power P_{i}^{band}. We then calculate relative enrichment for each band (on a log10 scale), as:
E_{i}^{band} = log10(P_{i}^{band})  log10(P^{band}).
The Q score for each spindle is based on the difference between the largest spindleband relative enrichment versus the largest nonspindle band relative enrichment:
Q_{i} = max( E_{i}^{slowsigma} , E_{i}^{fastsigma} )  max( E_{i}^{delta} , E_{i}^{theta} , E_{i}^{beta} )
That is, numbers less than zero indicate that a nonspindle band is showing greater relative enrichment (relative to the baseline logscaled power for that band) than either spindle band. Such spindles are more likely to reflect artifact which had a nonspecific effect on the power spectrum.
To illustrate this QC metric, we'll detect spindles for the second tutorial individual:
luna s.lst 2 sig=EEG o out.db \
s "MASK ifnot=NREM2 & RE & \
SPINDLES fc=11,15 enrich q=9 ftr=sp1 cycles=12"
Here, we detect both fast and slow spindles, targeting 11 and 15 Hz
respectively; we've set cycles to 12 instead of the default 7, in
order to increase frequencyspecificity. We've added the enrich
option to output
E_{i}^{slowsigma}
and
E_{i}^{fastsigma}
along with
E_{i}^{delta}
E_{i}^{theta}
and
E_{i}^{beta}
for each spindle/channel/targetfrequency/band (and so stratified by
CH
xF
xSPINDLE
xB
). We also set q
equal to 9, which effectively turns off any
QC filtering, so that we can see what some spindles that would not otherwise have been called
with the default threshold (q=0
) look like. Finally, the ftr
option writes out two files
that detail where spindles occur.
In lunaR:
library(luna)
lattach( lsl( "s.lst" ) , 2 )
ladd.annot.file( "id_nsrr02_feature_SPINDLESspindlesEEGwavelet11sp1.ftr" )
ladd.annot.file( "id_nsrr02_feature_SPINDLESspindlesEEGwavelet15sp1.ftr" )
out.db
:
k < ldb( "out.db" )
d < lx( k , "SPINDLES" , "F" , "CH" , "SPINDLE" )
par(mfcol=c(1,2))
hist( d$Q[ d$F == 11 ] , breaks=30 , main="11 Hz spindles" , xlab="Q" , xlim=c(2,2) )
abline(v=0,col="red",lwd=2)
hist( d$Q[ d$F == 15 ] , breaks=30 , main="15 Hz spindles" , xlab="Q" , xlim=c(2,2) )
abline(v=0,col="red",lwd=2)
Of the full set (i.e. running without the default q=0
threshold in place), we see that approximately 10% of spindles have q
less than 0:
table( d$Q > 0 , d$F )
11 15
FALSE 68 36
TRUE 405 370
Let's visualize some q<0
and q>0
spindles (we'll actually select
qc>1
spindles, to ensure that we don't randomly select spindles with
q
very near 0). For simplicity, we'll do this just for fast spindles.
First, we'll save the labels of the two new annotation tracks we
appended (i.e. you can see the names from running lannots()
):
a15 = "SPINDLESspindlesEEGwavelet15sp1"
s15 < lannots( a15 )
length(s15)
is indeed 406.
For visualization, here we'll just select five low Q and five high Q spindles:
lo15 < sample(d$SPINDLE[ d$Q < 0 & d$F == 15 ])[1:5]
hi15 < sample(d$SPINDLE[ d$Q > 1 & d$F == 15 ])[1:5]
Using lunaR functions, we can then write a function to plot the raw EEG signal for each spindle, with a 4second flanking window on either side:
f1 < function(i,a) {
x < ldata.intervals( i , chs = "EEG" , annots = a , w = 4 )
y < x$EEG
plot( x$SEC , y , type="l" , axes=F,xlab="",ylab="" , ylim=c(125,125) )
t < x$SEC[ x[,lsanitize(a)] == 1 ]
points( t , rep( 125,length(t)) , col=rgb(0,0,255,150,max=255) , pch="" )
}
Finally, we'll generate a matrix of plots, with the five low Q spindles in the left column, and the five high Q spindles in the right column:
par(mfrow=c(5,2),mar=c(0,0,0,0) )
for (i in 1:5) {
f1( i = s15[ lo15[ i ] ] , a = a15 )
f1( i = s15[ hi15[ i ] ] , a = a15 )
}
As you can see, the low Q spindles (that would not have been
detected using Luna's default q=0
setting) appear to reflect some
sources of artifact, whereas, as a whole, the high Q spindles look
more like canonical spindles.
Be aware that setting q
too high may unnecessarily miss many true
spindles though. A value of q=0.3
(i.e. log10(2)) requires that the
spindleband enrichment is at least twice nonspindle band enrichment;
a value of q=1
requires a 10fold enrichment. Considering the
distribution of mean Q over all individuals in the study may provide a good starting point,
if you want to select nondefault values.
Warning
Naturally, neither this particular automated spindle detection algorithm, nor this particular quality metric, are guaranteed to always get it right... We fully expect that this algorithm will flag as a "spindle" events that most scorers would not, and vice versa. When first approaching this tool, we'd advise experimenting with the different thresholds and parameters, and spotcheck studies in the way we have above. (The lunaR tutorial shows how to automatically iterate through all spindles detected.)
In general, if visual review seems to indicate likely artifact, the spindle is likely to be a false positive. In contrast, however, otherwise high Q spindles may appear relatively modest on visual review, and may not have been detected by some human scorers. Human scoring does not constitute a real goldstandard, however, and both human and automatic approaches impose inherently arbitrary thresholds on what is likely a continuouslyvarying and heterogeneous set of physiological processes. As such, we'd advise against automatically assuming the automatic detector is wrong, at least for this latter class of events ("low amplitude spindles").
Overall, the usual principles of a) trying to ensure that highquality data go into the analysis, and b) working with large samples and appropriate statistical techniques downstream, are likely to be just as (if not more) important than spending too much tweaking the detection procedure, or agonizing over what is or is not "a true spindle".
Collating spindles
When detecting spindles on multiple channels or at multiple, nearby frequencies, what presumably is the same underlying physiological event may often be detected multiple times. These options assist in collating the results of spindle detection across channels and/or frequencies:

you can control whether two temporally close spindles detected on the same channel/targetfrequency should be merged and considered as one event or not (
merge
parameter, the default is 0.5 seconds) 
second, within channel, you can collate temporally overlapping spindles detected at nearby frequencies (i.e. we would expect a great deal of overlap in the results of an analysis with F_{C}=10 Hz and F_{C}=11 Hz).

third, across channels, you can collate temporally overlapping spindles detected at nearby target frequencies; these may or may not represent physiologically distinct events depending on the montage and channels being collated.
These commands essentially take a list of spindles and generate another list of socalled mspindles: groups of spindles that, on some definition, "overlap", i.e. represent the same underlying event.
The primary output from this analysis is an estimate of mspindle
density (MSP_DENS
), i.e. which is taken to reflect the total number
of unique spindles (i.e. allowing for potentially overlapping
definitions). For each mspindle (set a spindles), Luna also
calculates an mean frequency (a weighted mean of the constituent
spindle frequencies, with weights given by the ISA of each spindle).
This is output for each spindle and also used to generate a
frequencyconditioned estimate of mspindle density.
Parameters
Parameter  Example  Description 

merge 
merge=0.2 
Merge two putative spindles if within this interval (default 0.5 seconds) 
collate 
collate 
Within each channel, collate overlapping spindles with similar frequencies (to MSPINDLES ) 
collatechannels 
collatechannels 
As above, except merge across channels also 
thfrq 
thfrq=1 
Criterion for merging two spindles if difference in frequency is within this range (default 2 Hz) 
listallspindles 
listallspindles 
List all spindles that comprise each mspindle 
Secondary parameters:
Parameter  Example  Description 

thinterval 
thinterval=0.5 
Merge two overlapping spindles if the ratio of intersection to union is at least this value (default 0, i.e. any overlap) 
thintervalcrosschannel 

thintervalwithinchannel 

window 
window=0.5 
Set window around each spindle when defining temporal overlap 
hms 
hms 
Show clocktime of each mspindle 
Output
Individuallevel summaries of mspindles (option: collate
, strata: none)
Variable  Description 

MSP_DENS 
mspindle density 
MSP_N 
mspindle count 
MSP_MINS 
Denominator for density, i.e. minutes of signal analyzed 
mspindle density stratified by mspindle frequency (option: collate
, strata: F
)
Variable  Description 

MSP_DENS 
mspindle density conditional on mspindle frequency 
Mergedspindle output (option: collate
, strata: MSPINDLE
)
or Mergedspindle output (option: collatewithinchannel
, strata: CH
x MSPINDLE
)
Variable  Description 

MSP_DUR 
Duration of this mspindle 
MSP_F 
Estimated frequency of this mspindle 
MSP_FL 
Lower frequency of this mspindle 
MSP_FU 
Upper frequency of this mspindle 
MSP_SIZE 
Number of spindles in this mspindle 
MSP_STAT 
Statistic for mspindle 
MSP_START 
Start time (seconds elapsed from EDF start) of mspindle 
MSP_STOP 
Stop time (seconds elapsed from EDF start) of mspindle 
Spindle to mspindle mappings (option: listallspindles
, strata: SPINDLE
x MSPINDLE
)
Variable  Description 

SCH 
Spindle label (channel:target frequency) 
FFT 
Spindle estimated frequency (via FFT) 
START 
Spindle start time (elapsed seconds from EDF start) 
STOP 
Spindle stop time (elapsed seconds from EDF start) 
Additional output (option: hms
, strata: MSPINDLE
or CH
x MSPINDLE
))
Variable  Description 

MSP_START_HMS 
Merged spindle start clocktime 
MSP_STOP_HMS 
Merged spindle stop clocktime 
Example
luna s.lst 2 sig=EEG o out.db s "MASK ifnot=NREM2 & RE & \
SPINDLES fclower=10 fcupper=16 fcstep=0.25 \
collate listallspindles"
The total mspindle density is 4.89 spindles per second, i.e. this estimates the total unique number of spindles over this broad frequency range ( with target frequencies of 1016 Hz and each wavelet spanning approximately plus/minus 2 Hz):
destrat out.db +SPINDLES
ID MSP_DENS MSP_MINS MSP_N
nsrr02 4.8972 199.5 977
We can also view the mspindle density as a function of the
estimated mspindle frequency: i.e. these numbers will sum to the
total MSP_DENS
estimate (or close to it, as we only consider the
fixed range of 8 to 16 Hz in 0.5 Hz bins, by default, when outputing
this):
destrat out.db +SPINDLES r F
Plotting the output of this command, we potentially see evidence of a bimodality in spindle frequencies for this individual:
Relationship to individualized spindle detectors
Rather than using fixed target frequencies (e.g. 13.5 Hz, 11 Hz, or 15 Hz, etc) to define "spindles", some approaches to spindle detection seek to find, empirically, the "peak" or "peaks" (which may or may not be obviously present in the power spectra) for each particular individual, and then target spindles of those particular frequencies for that individual. In contrast, this approach is effectively agnostic to an individual's "true" characteristic spindle frequency (or frequencies), or to having to define spindle frequencies up front, by virtue of simply considering all possible frequencies in the broad sigma range, and then combining results. (Arguably this approach might be preferred a priori in any case, on the grounds that even within a single individual we can see marked shifts in spindle frequency over the course of the night.)
Other metrics
The following options report some additional spindle morphology
characteristics (these are unlikely to be desired as part of any
primary spindle analysis, and so can be ignored by most users). The
if
option estimates spindle instantaneous frequency, averaged over
spindles, as a function of the relative location within the spindle.
The tlock
option produces an averaged EEG signature of detected
spindles, synced to the spindle peak (point of max peaktopeak
amplitude).
Parameters
Parameter  Example  Description 

if 
if 
Estimate instantaneous frequency of spindles 
iffrq 
iffrq=1 
Window around target frequency (default 2 hz) 
tlock 
tlock 
Flag to request (verbose) average, peaklocked waveforms 
Instantaneous frequency is estimated via the filterHilbert method,
where the filter is F_{C} +/ H Hz (where H is
set by iffrq
, default 2).
Output
Instantaneous frequency (IF) per spindle (option: if
, strata: CH
x F
x SPINDLE
)
Variable  Description 

IF 
Mean frequency per spindle over duration 
Mean IF stratified by relative location in spindle (option: if
, strata CH
x F
x RELLOC
)
Variable  Description 

IF 
Mean frequency of all spindles, per relative position within the spindle (five bins) 
Example
With the tutorial data:
luna s.lst 2 o out.db s "MASK ifnot=NREM2 & RE & \
SPINDLES sig=EEG fc=11,15 if tlock q=0.3 cycles=12"
FRQ
and FFT
,
as they are all estimating essentially the same thing):
destrat out.db +SPINDLES r F CH SPINDLE v FFT FRQ IF
To get the mean IF stratified by relative position with the spindle (i.e. dividing each spindle into five, equallysized time bins):
destrat out.db +SPINDLES r F CH RELLOC p 3
ID CH F RELLOC IF
nsrr02 EEG 11 1 11.321
nsrr02 EEG 11 2 11.218
nsrr02 EEG 11 3 11.041
nsrr02 EEG 11 4 10.912
nsrr02 EEG 11 5 10.818
nsrr02 EEG 15 1 14.287
nsrr02 EEG 15 2 14.119
nsrr02 EEG 15 3 13.959
nsrr02 EEG 15 4 14.015
nsrr02 EEG 15 5 14.181
To plot the mean signal timelocked to each spindle peak (in lunaR):
k < ldb("out.db")
d < lx( k , "SPINDLES" , "CH" , "F" , "MSEC" )
par( mfcol=c(1,2) )
for (f in c(11,15)) {
plot( d$MSEC[ d$F==f ] , d$TLOCK[ d$F==f ] , ylim=c(20,20) , type="l"
xlab="Time (sec)",ylab="Mean EEG",main=paste(f,"Hz spindles"), lwd=2)
abline(h=0,col="red",lty=2)
}
Annotations
Luna can generate FTR files representing the spindles detected in a given run. These can subsequently be attached to an EDF via lunaR, for example, in order to visualize spindles, as in the examples above.
Parameters
Parameter  Example  Description 

ftr 
ftr=tag 
Produce FTR files for all spindles, with the tag in the filename 
ftrdir 
ftrdir=/path/to/folder 
Folder for FTR files 
showcoef 
showcoef 
Flag to request (very verbose) coefficient output (to stdout ) 
Output
The ftr
option generates one file per channel/target frequency
combination, with a filename in the form:
id_{ID}_feature_spindles{CH}wavelet{Fc}{tag}.ftr
where ID is the EDF ID, CH is the channel label,
F_{C} is the target frequency, and tag is the
usersupplied tag (from the ftr
option).
Spindle/SO coupling
The sw
option will use a heuristic method to detect slow
oscillations (SO) in the sleep EEG, along with spindles, and will
quantify the coupling between spindles and SO, via eventrelated phase
locking measures. There are a number of additional parameters that
modify the behavior of sw
.
Warning
The documentation for spindle/SO coupling, and for SO detection in general is still skeletal. We're working on it...
Parameters
Bandpass filtering for SO detection:
Parameter  Example  Description 

flwr 
flwr=0.2 
Bandpass filter lower transition frequency 
fupr 
fupr=4.5 
Bandpass filter upper transition frequency 
Time criteria for SO detection:
Parameter  Example  Description 

tlwr 
tlwr= 
Minimum duration of a SO 
tupr 
tupr= 
Maximum duration of a SO 
tneglwr 
tneglwr=0.5 
Lower threshold for negative peak (seconds) 
tnegupr 
tneglwr=2 
Upper threshold for negative peak (seconds) 
Amplitude criteria for SO detection:
Parameter  Example  Description 

uVneg 
uVneg=75 
Negative peak amplitude threshold ( default 0, i.e. no threshold ) 
uVp2p 
uVp2p=140 
Peaktopeak amplitude (default 0, i.e. no threshold) 
mag 
mag=1 
Use a relative amplitude threshold instead (see below) 
swmean 
swmean 
Use mean (instead of median) 
Additional SO detection options (currently not fully supported, do not use):
Parameter  Example  Description 

neg2pos 
neg2pos 
Define full SO based on consecutive negativetopositive zerocrossings 
halfwave 
halfwave 
Detect all half waves instead of full waves 
negativehalfwave 
negativehalfwave 
Detect negative halfwaves instead of full waves 
positivehalfwave 
positivehalfwave 
Detect positive halfwaves instead of full waves 
Output
See the output tables below for the SO
command for a
description of the additional output generated by the sw
option of
SPINDLES
. In addition, there are a number of additional tables
detailing spindle/SO coupling generated by the sw
option.
Primary individuallevel spindle/SO coupling output (option: sw
, strata: CH
x F
)
Variable  Description 

ITPC_ANGLE_PEAK 
Circular mean of SO phase at spindle peak (SW_SPHASE_PEAK below) 
ITPC_PEAK 
Intertrial (phase) coherence (ITPC) statistic 
ITPC_PV_PEAK 
ITPC pvalue 
PROP_ON_SW 
Proportion of spindles temporally overlapping a detected SO 
PROP_PEAK_ON_SW 
Proportion of spindles whose peak temporally overlaps a detected SO 
Spindlelevel spindle/SO coupling output (option: sw
, strata: CH
x F
x SPINDLE
)
Variable  Description 

SW_SPHASE_PEAK 
SO phase at spindle peak 
NEAREST_SW 
Time to nearest SO 
NEAREST_SW_NUM 
Number of nearest SO 
SO morphology (option: sw
, strata: CH
x PHASE
)
Variable  Description 

SWPL_EEG 
Phaselocked mean EEG value (as a function of SO phase, in 20deg bins) 
SO morphology (option: sw
, strata: CH
x SP
)
Variable  Description 

SWTL_EEG 
Timelocked mean EEG value (as a function of samplepoints before/after SO negative peak) 
Spindle/SO phase coupling (option: sw
, strata: CH
x F
x PHASE
Variable  Description 

SWPL_CWT 
SOphased locked mean spindle wavelet coefficient (18 20deg bins) 
Spindle/SO phase coupling (option: sw
, strata: CH
x F
x SP
Variable  Description 

SWTL_CWT 
Timelocked mean spindle wavelet coefficient 
Example
In the literature, different investigators have used different
definitions for SOs in different contexts. To detect SO based on
relative amplitude of waveforms (relative to the whole night/stage for
that individual), you could use something like (note, these are
options added to a SPINDLES
command, although we are not showing the
full command here):
sw flwr=0.3 fupr=4 tlwr=0.5 tupr=2 mag=1
Here, the EEG channel is first bandpass filtered with transition
frequencies 0.3 and 4 Hz. After counting all subsequent
positivetonegative zerocrossings, only intervals between 0.5 and 2
seconds in duration are retained as putative SO. Further, only
intervals where both the negative peak amplitude and peaktopeak
amplitude are at least 1.0 (mag
) times the median for that
individual are designated as the final set of SO.
Instead of using a relative amplitude threshold (mag
), you can set
absolute thresholds. For example, delta waves with a negative peak
of at least 40 uV, and at least 75 uV peaktopeak:
sw flwr=0.3 fupr=4 uVneg=40 uVp2p=75
sw flwr=0.3 fupr=4 uVneg=75 uVp2p=140
You can also specify separate thresholds for the duration of the negative and positive half waves:
tneglwr=0.3 tnegupr=1.5 tposlwr=0.0 tposupr=1.0
Warning
This section is still under development. We'll add a specific vignette on using Luna for spindle/SO coupling analyses in the near future.
Misc
showcoef
The showcoef
is primarily for development purposes and should not
be part of a standard spindleanalysis workflow. It produces a lot
of output: five values for every timepoint for every channel for
every spindle frequency: this will be hard work for the output
database, which is not really designed for this type of data: you are
better off using this without the o
database option; rather,
extract the desired information directly from standard output:
luna s.lst 1 < cmd.txt \
 awk ' $2 == "SPINDLES" && $5 == "CWT" { print $3,$4,$6 } ' OFS="\t"
Variable  Description 

RAWCWT 
Absolute, raw wavelet coefficient value 
CWT 
Normalized wavelet coefficient value 
CWT_TH 
Wavelet threshold for core spindle 
CWT_TH2 
Wavelet threshold for flanking spindle 
CWT_THMAX 
Wavelet max threshold for spindle 
SO
Detects slow oscillations via a simple heuristic approach
The SO
command uses a a heuristic approach to detect slow
oscillations in the EEG, combined with a filterHilbert approach to
estimate SO phase. Currently can only be applied to a single channel
at a time.
Parameters
see spindle/SO coupling section above
Additional options;
Parameter  Example  Description 

tl 
tl=sig1,sig2 
Other signals to be summarized w.r.t. timelocked SO 
onset 
onset 

pos 
pos 

window 
window=1 
Window around SO peak (default +/ 3 seconds) 
Outputs
Note: the output uses SW
(slow wave) as equivalent to SO (slow
oscillation). Also, additional variables starting SW_MD
reflect the
median rather than the mean across events.
SO information per channel (option: sw
, strata: CH
)
Variable  Description 

SW 
Number of SO detected 
SW_RATE 
SO per minute 
SW_TH_NEG 
Threshold used for negative peak 
SW_TH_P2P 
Threshold used for peaktopeak 
SW_DUR 
Median SO duration 
SW_NEG_DUR 
Median negative peak duration 
SW_POS_DUR 
Median positive peak duration 
SW_AMP 
Median amplitude (of negative peak) 
SW_P2P 
Median peaktopeak amplitude 
SW_SLOPE_NEG1 
Median slope (positivetonegative zerocrossing to negative peak) 
SW_SLOPE_NEG2 
Median slope (negative peak to negativetopositive zerocrossing) 
SW_SLOPE_POS1 
Median slope (negativetopositive zero crossing to positive peak) 
SW_SLOPE_POS2 
Median slope (positive peak to positivetonegative zerocrossing) 
Per SO output (strata: CH
x N
)
Variable  Description 

START_IDX 
Start of SO (in samplepoint units) 
STOP_IDX 
Stop of SO (in samplepoint units) 
START 
Start of SO (in seconds elapsed from start of EDF) 
STOP 
Stop of SO (in seconds elapsed from start of EDF) 
DUR 
SO duration 
UP_AMP 
SO positive peak amplitude 
DOWN_AMP 
SO negative peak amplitude 
P2P_AMP 
SO peaktopeak amplitude 
UP_IDX 
Sample index of positive peak 
DOWN_IDX 
Sample index of negative peak 
SLOPE_POS1 
Slope (as above) (unless halfwaves) 
SLOPE_POS2 
Slope (as above) (unless halfwaves) 
SLOPE_NEG1 
Slope (as above) (unless halfwaves) 
SLOPE_NEG2 
Slope (as above) (unless halfwaves) 
Verbose SOlevel out (strata: CH
x N
x SP
)
Variable  Description 

FLT 
Filtered waveform for each SO 
PH 
Instantaneous phase for each SO 
Epochlevel counts of SO (strata: CH
x E
)
Variable  Description 

N 
Number of SOs detected in that epoch 
From tl
( strata: CH
x CH2
x SP
)
Variable  Description 

SWTL_SIG 

SWTL 
Example
to be added