slow waves and sleep spindles are the essential elements in the definition of sleep stages .
they recently received additional attention because of their hypothesized involvement in processes related to neuronal plasticity , learning , and memory [ 24 ] . furthermore , slow oscillations ( <  1  hz ) are considered as a central feature of sleep .
these intracortically generated fluctuations consist at the cellular level of rhythmic depolarizing components ( up states ) separated by prolonged hyperpolarizations ( down states ) .
there is evidence that they are involved in the temporal organization of other sleep rhythms such as spindles and delta waves ( for a review , see ) .
also , it is hypothesized that the slow oscillation is involved in processes of learning and memory during sleep .
tononi and cirelli proposed that slow oscillations were an indication of synaptic downscaling in the cortex occurring during nonrapid eye movement sleep in order to counteract synaptic potentiation resulting from learning processes occurring during wakefulness [ 6 , 7 ] .
evidence for this hypothesis comes from data of a motor learning task , arm immobilization , animal studies , and computational modeling .
we shall address in this paper only one particular aspect in this context , which is , however , related to the interplay of sleep spindles and slow oscillations : i.e. , the temporal organization of spindles .
it was postulated that these thalamically generated sleep spindles are grouped by cortically generated slow oscillations and more generally that the slow oscillations also organize the occurrence of other sleep oscillations .
this hypothesis is based on experimental and theoretical evidence : sleep spindles lost their coherence and their temporal 
however , there are also inconsistencies in relation to this hypothesis : first , the frequency of the slow oscillations was originally considered as being only slightly <  1  hz . achermann and borbly , for instance , reported a peak in the power spectrum at 0.7  hz .
this would correspond to a shorter typical interevent interval between sleep spindles than 4  s. thus , only approximately every third slow oscillation would trigger a spindle .
second , the rhythmic occurrence of spindles is most pronounced in sleep stage 2 ( see fig .  4 ) , but slow oscillations are most prominent in deep sleep .
thus , one should also take alternative explanations into account . in the present paper , we investigated whether the  periodicity  in the occurrence of sleep spindles and other sleep oscillations might be a pure stochastic effect that occurs without assuming an additional slow process . to test this hypothesis , we compared interevent interval distributions of the events detected in the human sleep eeg with events detected in artificial data generated with stationary stochastic processes , which resemble to some extent the short - term properties of the eeg signal .
the analysis is based on a recently proposed method to detect and analyze oscillatory events in the human sleep eeg .
eeg data were modeled on short time scales (   1s ) by linear autoregressive ( ar ) models in a sliding window fashion .
systems described by such linear models can be interpreted as a set of coupled stochastically driven harmonic oscillators with time - dependent frequencies and damping coefficients . in the case of the human sleep eeg , this frequency might represent resonances of the underlying thalamocortical networks .
in , it was shown that the modulation of the frequency and the damping constant of the resonance of a simple excitatory 
inhibitory network can be detected from the mean activity of the network by applying this method .
we analyzed 8-h sleep eeg data ( single channel , derivation c3a2 ) of eight healthy young male subjects ( mean age 24.1 0.6  years ) , each contributing four baseline nights [ 17 , 19 ] .
the sampling rate was 128  hz and sleep stages were visually scored according to standard criteria .
oscillatory events were detected by fitting ar models of order 8 to subsegments of 1  s in length in a sliding window fashion resulting in a time dependence of the parameters of the ar model .
the signal generated by an ar(p)-model can be interpreted as the superposition of the output of n oscillators and p    2n relaxators .
the p parameters of the ar(p)-model are transformed into n frequencies and damping constants of the oscillators , and p    2n damping constants of the purely decaying nonoscillating modes .
thus , time - dependent frequencies and damping constants result from the sliding window approach .
figure  1 shows a 20-s segment of eeg ( stage 2 ) and the time - dependent modules r and the corresponding frequencies \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{mathrsfs } 
 \usepackage{upgreek } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document}$f\:$\end{document } of the poles .
the damping constant  is a function of r ,   =    fs logr with fs denoting the sampling frequency .
the case r  =  1 corresponds to the limiting case of an undamped oscillation .
the threshold for the event detection was set to r1  =  0.95 ( fig . 
=  0.9 was applied : in order to account for statistical fluctuations , the event was not necessarily terminated whenever r became smaller than r1 .
it was considered as a single event if the corresponding r became larger than r1 again without falling below r2 .
the value r1 was chosen such that clearly visible sleep spindles were reliably detected by the algorithm , while r2 had to be low enough so that a single spindle was detected as a single event , but high enough so that subsequent spindles were recognized as distinct events .
the time between the first and the final crossing of the threshold r1 plus 1  s ( the length of the segment used to fit the ar model ) was considered as the duration of the oscillatory event . the time when r was maximal was defined as the time of occurrence of the event and the corresponding frequency
the interevent intervals were estimated as the difference of the time of occurrence of successive oscillatory events . 

1detection of oscillatory events in a 20-s eeg segment of stage 2 non - rem sleep .
the upper panel shows the eeg signal , the middle panel shows r , and the lower panel the corresponding frequencies f , with the colors indicating corresponding oscillatory modes .
the detection threshold r1  =  0.95 is indicated in the middle panel ( see section  2.2 ) detection of oscillatory events in a 20-s eeg segment of stage 2 non - rem sleep .
the upper panel shows the eeg signal , the middle panel shows r , and the lower panel the corresponding frequencies f , with the colors indicating corresponding oscillatory modes .
the detection threshold r1  =  0.95 is indicated in the middle panel ( see section  2.2 ) for the analysis of the interevent intervals and the event durations , we assigned the events into different frequency bands .
we are aware that this is not an optimal approach because of interindividual variation in mean event frequencies and their dependence on sleep stage and cycle .
however , at the moment we have no method available to analyze the events on the basis of individually determined time - dependent frequency bands .
therefore , we used the frequency bands employed previously , i.e. , delta (  : 0    f  <  4.5  hz ) , alpha (  : 8    f  <  11.5  hz ) , and sigma (  : 11.5    f  <  16  hz ; sleep spindles ) bands .
figure  2 illustrates the detected events of a single night and the corresponding histogram of the event frequencies .
2detected oscillatory events in a single night ( black dots represent single events ) , hypnogram ( mt , movement time ; w , waking , rem , rapid eye movement sleep ; 1 to 4 , nonrapid eye movement sleep stages 1 to 4 ) , and corresponding histogram of the event frequencies .
colors indicate the different frequency bands detected oscillatory events in a single night ( black dots represent single events ) , hypnogram ( mt , movement time ; w , waking , rem , rapid eye movement sleep ; 1 to 4 , nonrapid eye movement sleep stages 1 to 4 ) , and corresponding histogram of the event frequencies .
colors indicate the different frequency bands artificial eeg data were generated using an ar(4 ) model with random inputs and parameters corresponding to a frequency in the spindle frequency band ( f  =  14  hz ) and in the delta band ( f  =  3  hz ) , assuming a sampling frequency of 128  hz .
we restricted the analysis to data from an ar(4 ) model because the two higher frequencies that can be seen in fig . 
we analyzed data sets with different values of r for the spindle frequency because the spindle oscillations are our main focus .
for each set of parameter values , a time series of 6 10 data points ( corresponding to approximately 13  h of eeg data ) was generated .
3twenty - second simulated ar(4 ) data with time - dependent modules and frequencies derived by fitting an ar(4 ) model similar to fig .  1 .
the upper panel shows the signal , the middle panel shows r , and the lower panel the corresponding frequencies f , with the colors indicating corresponding oscillatory modes .
the ar(4 ) model contained two oscillatory modes , one in the spindle frequency band , f  =  14 hz , r  =  0.94 , and one in the delta frequency band , f  =  3hz , r  =  0.91 twenty - second simulated ar(4 ) data with time - dependent modules and frequencies derived by fitting an ar(4 ) model similar to fig .  1 .
the upper panel shows the signal , the middle panel shows r , and the lower panel the corresponding frequencies f , with the colors indicating corresponding oscillatory modes .
the ar(4 ) model contained two oscillatory modes , one in the spindle frequency band , f  =  14 hz , r  =  0.94 , and one in the delta frequency band , f  =  3hz , r  =  0.91
we analyzed 8-h sleep eeg data ( single channel , derivation c3a2 ) of eight healthy young male subjects ( mean age 24.1 0.6  years ) , each contributing four baseline nights [ 17 , 19 ] .
the sampling rate was 128  hz and sleep stages were visually scored according to standard criteria .
oscillatory events were detected by fitting ar models of order 8 to subsegments of 1  s in length in a sliding window fashion resulting in a time dependence of the parameters of the ar model .
the signal generated by an ar(p)-model can be interpreted as the superposition of the output of n oscillators and p    2n relaxators .
the p parameters of the ar(p)-model are transformed into n frequencies and damping constants of the oscillators , and p    2n damping constants of the purely decaying nonoscillating modes .
thus , time - dependent frequencies and damping constants result from the sliding window approach .
figure  1 shows a 20-s segment of eeg ( stage 2 ) and the time - dependent modules r and the corresponding frequencies \documentclass[12pt]{minimal } 
 \usepackage{amsmath } 
 \usepackage{wasysym } 
 \usepackage{amsfonts } 
 \usepackage{amssymb } 
 \usepackage{amsbsy } 
 \usepackage{mathrsfs } 
 \usepackage{upgreek } 
 \setlength{\oddsidemargin}{-69pt } 
 \begin{document}$f\:$\end{document } of the poles .
the damping constant  is a function of r ,   =    fs logr with fs denoting the sampling frequency .
the case r  =  1 corresponds to the limiting case of an undamped oscillation .
the threshold for the event detection was set to r1  =  0.95 ( fig .  1 ) .
to determine the event duration , a second threshold r2  =  0.9 was applied : in order to account for statistical fluctuations , the event was not necessarily terminated whenever r became smaller than r1 .
it was considered as a single event if the corresponding r became larger than r1 again without falling below r2 .
the value r1 was chosen such that clearly visible sleep spindles were reliably detected by the algorithm , while r2 had to be low enough so that a single spindle was detected as a single event , but high enough so that subsequent spindles were recognized as distinct events .
the time between the first and the final crossing of the threshold r1 plus 1  s ( the length of the segment used to fit the ar model ) was considered as the duration of the oscillatory event . the time when r was maximal was defined as the time of occurrence of the event and the corresponding frequency
the interevent intervals were estimated as the difference of the time of occurrence of successive oscillatory events . 

1detection of oscillatory events in a 20-s eeg segment of stage 2 non - rem sleep .
the upper panel shows the eeg signal , the middle panel shows r , and the lower panel the corresponding frequencies f , with the colors indicating corresponding oscillatory modes .
the detection threshold r1  =  0.95 is indicated in the middle panel ( see section  2.2 ) detection of oscillatory events in a 20-s eeg segment of stage 2 non - rem sleep .
the upper panel shows the eeg signal , the middle panel shows r , and the lower panel the corresponding frequencies f , with the colors indicating corresponding oscillatory modes
. the detection threshold r1  =  0.95 is indicated in the middle panel ( see section  2.2 ) for the analysis of the interevent intervals and the event durations , we assigned the events into different frequency bands .
we are aware that this is not an optimal approach because of interindividual variation in mean event frequencies and their dependence on sleep stage and cycle .
however , at the moment we have no method available to analyze the events on the basis of individually determined time - dependent frequency bands
. therefore , we used the frequency bands employed previously , i.e. , delta (  : 0    f  <  4.5  hz ) , alpha (  : 8    f  <  11.5  hz ) , and sigma (  : 11.5    f  <  16  hz ; sleep spindles ) bands .
figure  2 illustrates the detected events of a single night and the corresponding histogram of the event frequencies .
2detected oscillatory events in a single night ( black dots represent single events ) , hypnogram ( mt , movement time ; w , waking , rem , rapid eye movement sleep ; 1 to 4 , nonrapid eye movement sleep stages 1 to 4 ) , and corresponding histogram of the event frequencies .
colors indicate the different frequency bands detected oscillatory events in a single night ( black dots represent single events ) , hypnogram ( mt , movement time ; w , waking , rem , rapid eye movement sleep ; 1 to 4 , nonrapid eye movement sleep stages 1 to 4 ) , and corresponding histogram of the event frequencies .
artificial eeg data were generated using an ar(4 ) model with random inputs and parameters corresponding to a frequency in the spindle frequency band ( f  =  14  hz ) and in the delta band ( f  =  3  hz ) , assuming a sampling frequency of 128  hz .
we restricted the analysis to data from an ar(4 ) model because the two higher frequencies that can be seen in fig . 
we analyzed data sets with different values of r for the spindle frequency because the spindle oscillations are our main focus .
for each set of parameter values , a time series of 6 10 data points ( corresponding to approximately 13  h of eeg data ) was generated .
3twenty - second simulated ar(4 ) data with time - dependent modules and frequencies derived by fitting an ar(4 ) model similar to fig .  1 .
the upper panel shows the signal , the middle panel shows r , and the lower panel the corresponding frequencies f , with the colors indicating corresponding oscillatory modes .
the ar(4 ) model contained two oscillatory modes , one in the spindle frequency band , f  =  14 hz , r  =  0.94 , and one in the delta frequency band , f  =  3hz , r  =  0.91 twenty - second simulated ar(4 ) data with time - dependent modules and frequencies derived by fitting an ar(4 ) model similar to fig .  1 .
the upper panel shows the signal , the middle panel shows r , and the lower panel the corresponding frequencies f , with the colors indicating corresponding oscillatory modes .
the ar(4 ) model contained two oscillatory modes , one in the spindle frequency band , f  =  14 hz , r  =  0.94 , and one in the delta frequency band , f  =  3hz , r  =  0.91
figure  4 shows the interevent interval distributions of sleep spindles ( 11.516  hz ) .
the histograms of the empirical data were estimated based on the four recordings of all subjects ( pooled data ) . in sleep stage 2 ,
the maximum is at t    4 s and is much more pronounced than in deep sleep [ slow wave sleep ( sws ) ] .
the distribution of the interevent intervals detected by our algorithm was similar to the distribution reported in   for spindles identified by visual detection .
figure  4 also contains the distributions of events detected from simulated data of stationary ar(4 ) models for two different values of the module r for the spindle frequency .
events were detected with the same algorithm and parameters as for the eeg data . also , for the simulated data , the interevent interval distributions showed a peak around 4  s ( 2.5  s for r  =  0.94 and between 2.5 and 4 for r  =  0.89 ) .
the results of the simulated data indicate that stationary stochastic processes also exhibit a periodicity in the occurrence of events with a maximum similar to empirical data without the assumption of an additional slow process grouping the events .
thus , in order to understand the character of the maxima in the empirical distributions , it is necessary to investigate these distributions in more detail . 

4normalized histograms ( pooled data ) of interevent intervals ( t ) in the spindle frequency range ( 11.516  hz ) in nonrapid eye movement sleep stage 2 ( st2 ; top left ) and slow wave sleep ( sws , stages 3 and 4 ; top right ) and of events in the spindle frequency range detected from data generated by stationary ar(4 ) models with r  =  0.94 ( bottom left ) and r  =  0.89 ( bottom right ) for the spindle frequency normalized histograms ( pooled data ) of interevent intervals ( t ) in the spindle frequency range ( 11.516  hz ) in nonrapid eye movement sleep stage 2 ( st2 ; top left ) and slow wave sleep ( sws , stages 3 and 4 ; top right ) and of events in the spindle frequency range detected from data generated by stationary ar(4 ) models with r  =  0.94 ( bottom left ) and r  =  0.89 ( bottom right ) for the spindle frequency in fig .  5 ,
an exponential decay for large interevent intervals was present in both the synthetic and the empirical data . for events occurring randomly with a certain constant rate ( poisson process )
, one would expect that the interevent intervals are exponentially distributed , with the exponent being equal to the event rate ( number of events per unit of time ) .
in our situation , this can not be exactly the case because of the finite duration of the oscillatory events .
however , for large interevent intervals compared to the event duration , this might be a reasonable assumption .
the empirical distributions differ substantially from an exponential distribution not only for times below the maximum of the distribution but also for larger time intervals .
thus , the region around the maximum can not be explained in the same way as the distributions for the ar processes .
this is the first indication of an additional process contributing to the peak in the empirical interevent distributions . 

5normalized histograms of interevent intervals as in fig .  4 but plotted on a semilogarithmic scale including additional values of r. the straight lines show fits to the exponential tails of these distributions normalized histograms of interevent intervals as in fig .  4 but plotted on a semilogarithmic scale including additional values of r. the straight lines show fits to the exponential tails of these distributions the systematic dependence of the exponential decay on r of the ar model and , therefore , on the damping of the spindle frequency pole , is illustrated in fig . 
6a . for r    0.93 , the exponent becomes smaller with smaller r. a similar relationship holds for the event rate ( fig . 
the inverse relationship is observed only for large r , i.e. , lower damping leads to a decrease in the event rate , which is overcompensated by an increased event duration .
the relationship between the exponent of the exponential decay and the event rate is shown for r    0.93 in fig . 
systematic deviations for larger event rates might be explained by more events of longer duration as illustrated in fig  7
6a exponents from the fit of the exponential tails of the interevent interval distributions ( fig .  5 ) of the simulated ar(4 ) data as a function of r of the spindle frequency pole .
b event rate as a function of r. c exponent of the exponential decay ( mean  sd ) as a function of the event rate shown only for 0.85    r    0.93 .
values for the empirical distributions ( red ) from sleep spindles shown separately for nonrapid eye movement sleep stages 2 , 3 , and 4 .
7normalized histograms ( pooled data ) of the event duration ( td ) in the spindle frequency range ( 11.516  hz ) in nonrapid eye movement sleep stage 2 ( st2 ; top left ) and slow wave sleep ( sws , stages 3 and 4 ; top right ) and events in the spindle frequency range from a stationary ar(4 ) process with different values of r for the spindle frequency a exponents from the fit of the exponential tails of the interevent interval distributions ( fig .  5 ) of the simulated ar(4 ) data as a function of r of the spindle frequency pole .
b event rate as a function of r. c exponent of the exponential decay ( mean  sd ) as a function of the event rate shown only for 0.85    r    0.93 .
values for the empirical distributions ( red ) from sleep spindles shown separately for nonrapid eye movement sleep stages 2 , 3 , and 4 .
the straight line indicates the linear relationship as expected for a poisson process normalized histograms ( pooled data ) of the event duration ( td ) in the spindle frequency range ( 11.516  hz ) in nonrapid eye movement sleep stage 2 ( st2 ; top left ) and slow wave sleep ( sws , stages 3 and 4 ; top right ) and events in the spindle frequency range from a stationary ar(4 ) process with different values of r for the spindle frequency the exponent of the empirical data of sleep spindles ( fig . 
6c ) suggests an increased damping ( corresponding to a smaller r ) with increasing sleep depth from stage 2 through stage 3 to stage 4 , which meets our expectations .
the spindles exhibited , however , a larger exponent than expected from the resulting event rate . or , to put it the other way around , the event rate is smaller than one would expect from the exponent of the exponential decay . from the 4-s  periodicity ,  one would have expected the contrary , i.e. , a larger event rate , because an interevent interval of 4  s corresponds to an event rate of 15 events per minute . a possible explanation for this discrepancy might be the estimation of the event density of spindles : the number of detected spindles during a particular sleep stage was divided by the total time spent in this stage .
this event rate is not equal to the inverse of the mean of the interevent interval distribution because only interevent intervals with no sleep stage transition between the events were included .
the normalized distribution of the event duration ( fig .  7 ) plotted on a logarithmic scale seems to be less specific for sleep spindles .
first , the maximum around t  =  4  s occurs in the simulated data for a wide range of values of r. second , there is no large difference between the distributions of the occurrence of spindles in sleep stage 2 and sws .
there is a slightly larger probability to have longer spindles in sleep stage 2 compared to sws ( fig . 
7 ) , which would correspond to a higher value of r if we assume that spindles were generated by a stationary ar process .
note , however , that the probability of long spindles increased strongly in data generated by ar models of large r. this increase might be the reason for the slower exponential decrease in the interevent interval distribution for r  =  0.95 compared to r  =  0.93 .
we have demonstrated that important properties of the temporal organization of the occurrence of sleep spindles are already present for oscillatory events detected in data generated with stationary stochastic processes .
such properties are the position of the maximum of the interevent interval distribution , i.e. , the  periodicity  in the occurrence of events between 2 and 5  s and the maximum around 2  s in the distribution of event durations .
a more detailed analysis , however , revealed that the exponential decay in the interevent distributions , which can be interpreted as reflecting the random occurrence of oscillatory events , starts in the case of the sleep spindles only at much longer interevent intervals compared to data from stationary ar processes .
for interevent intervals around the maximum of 4  s , the empirical distributions of sleep spindles clearly differ from the distributions of the simulated events and show , in particular , no exponential decay .
one could assume in a first approximation a gaussian distribution around 4  s , which would indicate that an additional process is responsible for the observed  periodicity  in the occurrence of spindles .
including such a process in modeling the time series , i.e. , modeling explicitly also the dynamics of the parameters of the ar model , would result in a nonlinear model .
the extent to which such a process is related to cortical slow oscillations has to be the subject of further research .
in particular , the present analysis could be extended by ( 1 ) also varying the delta frequency and/or its damping in the simulated data and ( 2 ) taking into account that single sleep stages are not dynamically homogeneous , i.e. , not assuming that the spindles in a particular sleep stage were generated by a single stationary process but rather by a suitable ensemble of such processes .
moreover , it is known that there are strong interindividual differences in the event rates of sleep spindles ( e.g. , ) .
thus , one should also investigate to which extent interindividual differences in the interevent interval distributions are present .
our analysis demonstrated , nonetheless , how phenomenological modeling of the short - term dynamics of the eeg using linear models can be used to detect and to analyze oscillatory events and to test hypotheses about the origin of some of their properties .
 the sleep electroencephalogram ( eeg ) is characterized by typical oscillatory patterns such as sleep spindles and slow waves . 
 recently , we proposed a method to detect and analyze these patterns using linear autoregressive models for short (   1  s ) data segments . 
 we analyzed the temporal organization of sleep spindles and discuss to what extent the observed interevent intervals correspond to properties of stationary stochastic processes and whether additional slow processes , such as slow oscillations , have to be assumed . 
 we have found evidence for such an additional slow process , most pronounced in sleep stage 2 . 
Introduction
Data and Methods
Data
Oscillatory Events
Simulated Data
Results
Discussion