Random walk model simulates the increased drowsiness of children with obstructive sleep apnea

Human sleep is characterized by rapid eye movement (REM) sleep as well as by light and deep sleep (non-REM). In addition, sleep is often interrupted by brief awakenings (arousals) that last from few seconds to minutes [1,2]. Most of the modeling approaches to sleep focus on the diurnal and circadian sleep/wake cycle [3–8], and models based on flip-flop switches were developed to account for the transitions between sleep and wake state as well as between REM and non-REM sleep [8–10]. There have been a few studies analyzing the probability distributions of arousal durations during sleep revealing a power-law behavior [1,11,12], which is consistently observed for several mammalian species [13] and seems to be related to different stages of maturation [14,15]. First attempts to simulate the sleep microarchitecture and the underlying power-law in arousal durations utilize stochastic and random walk based models [1,16], however, the complex dynamics of sleep-stage transitions and their alteration with sleep disorders such as obstructive sleep apnea (OSA) is not yet understood.

In this letter, we focus on non-REM sleep transitions and changes in their dynamics that possibly occur in children with OSA. Here, we propose a new parameter L that is related to the depth of light sleep, and we show its utility in characterizing the severity of OSA. In the current guidelines of the American Academy of Sleep Medicine (AASM), non-REM sleep is divided into three stages according to sleep depth: stages N1 and N2 (light sleep) as well as stage N3 (deep sleep) [17]. The sleep stage N1 represents sleep in its lightest form and is therefore considered a transitional state between wakefulness and sleep. N1 sleep is characterized by drowsiness [18] during which individuals usually respond to auditory or visual stimuli [19]. The N1 stage accounts for only a small portion of the total sleep time (TST), and therefore has often been ignored or grouped together with stage N2 into a single light sleep stage [20,21]. However, since N2 is more dominant (>50% of TST) than N1 (<10% of TST) [22], combining N1 with N2 might abolish any N1 effects that are related, for example, to cognitive performance. Indeed, very recent studies have shown a relationship between N1 (but not N2) percentage and cognitive impairment. Specifically, high percentage of N1 was found to be a good predictor of declining cognitive function in elderly subjects [23,24]. Additionally, it has been shown that N1 percentage increases for OSA patients (for both adults [25–27] and children [28,29]), whereas no difference has been found for N2 [25,30]. Furthermore, for children, OSA and sleep-disordered breathing in general, have been linked to hyperactivity [31] and poor school performance [32,33], and timely diagnosis and treatment of OSA can lead to immediate improvements [34,35]. In this work we focus on N1 sleep bout durations $(\Delta t_{N1})$, and based on a statistical physics model, we derive N1 sleep characteristics that can aid in better understanding sleep in OSA patients.

Recently, a model has been developed to simulate brief arousals from sleep [12,16]. The model is based on a physiological mechanism that takes into account the intrinsic noise ("sub-threshold voltage fluctuations") of wake-promoting neurons (WPN) in the brain stem. Although WPN are suppressed during sleep [10,36], the integrated uncorrelated neuronal noise from WPN can from time to time cross the excitability threshold and through projections to the cortex provoke an arousal. Because of the excitation along the ascending arousal pathway, sleep-promoting neurons (SPN) in the ventrolateral (VLPO) and median (MnPO) preoptic nuclei are stimulated [37,38], and can in turn actively inhibit the WPN (negative feedback) through an inhibitory current [12]. The integrated voltage V(t) of WPN can be modeled as a random walk with drift (during arousals/short awakenings) or without drift (during sleep) as shown in fig. 1. During sleep, the random walker transitions through the different non-REM sleep stages, and in order to quantify the probability of V(t) being in a state corresponding to N1, here we introduce a threshold L (fig. 1) that distinguishes between shallow and deeper light sleep, N1 $(0\le V <-L)$ and N2 $(-L\le V <-D)$ light sleep, respectively. In the following, we propose an analytical method to evaluate L from standard sleep recordings (hypnograms), and we determine L for healthy children and children with OSA. We find that OSA yields higher values of L that can be reduced through proper treatment.

Zoom In Zoom Out Reset image size

We analyze hypnograms of two groups of children that are age- and BMI-matched: i) healthy children (n = 38; age: $6.7 \pm 1.9$ years; BMI: $17.7 \pm 4.3\ \text{kg/m}^{2}$), and ii) OSA children (n = 58; age: $6.5 \pm 2.3$ years; BMI: $18.6 \pm 5.0\ \text{kg/m}^{2}$). In addition, we analyze the hypnograms of n = 9 children that received successful OSA treatment of either CPAP or surgery (age: $7.5 \pm 2.3$, BMI: $18.1 \pm 5.8\ \text{kg/m}^{2}$). Hypnograms were scored in 30 s epochs [17] and contain the following stages: wakefulness, REM sleep, and non-REM sleep stages N1, N2 and N3.

During sleep, the random walker transitions freely through the non-REM sleep stages N1, N2 and N3 (fig. 1), and changes in the integrated neuronal voltage V(t) of WPN can be calculated by the following differential equation:

$$\begin{equation} \text{d}V = \sigma\cdot \text{d}w,\quad \text{for } -\Delta\le V\le 0, \end{equation} \tag{ 1 }$$

where σ is the standard deviation of the neuronal voltage fluctuations, w is a standard Wiener process, and $-\Delta$ is the relative voltage value between threshold potential and potassium Nernst potential [12]. To model the probability distribution of $\Delta t_{N1}$, we introduce a N1 light sleep threshold, L, so that, for voltages between 0 and $-L$, N1 sleep occurs (fig. 1). In the following, we show how to estimate L from standard hypnogram data.

During N1 light sleep, the voltage is bounded within $-L\le V\le 0$ (fig. 1). Therefore, integrating eq. (1) in this range yields

$$\begin{equation} 0-(-L)=\sigma\cdot \int \text{d}w \end{equation} \tag{ 2 }$$

or

$$\begin{equation} {L}^2=\sigma^2 \cdot \left(\int \text{d}w \right)^2. \end{equation} \tag{ 3 }$$

Since the differentials $\text{d}w$ of a Wiener process are uncorrelated and have zero mean, averaging eq. (3) over all N1 sleep bouts yields

$$\begin{eqnarray} L^2&= \sigma^2 \cdot \int \langle (\text{d}w)^2 \rangle, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} L^2&= \sigma^2 \cdot \int \text{d}t, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} L^2&= \sigma^2 \cdot \tau_{N1}. \end{eqnarray} \tag{ 6 }$$

In the step from eq. (4) to eq. (5), we use a property of w (standard Wiener process) that the variance of $\text{d}w$ is the time interval $\text{d}t$ of the process, and therefore $\tau_{N1}$ (in eq. (6)) is the average duration of an N1 sleep bout. Thus, for the N1 light sleep threshold we obtain

$$\begin{equation} L=\sqrt{\tau_{N1}} \cdot \sigma. \end{equation} \tag{ 7 }$$

As seen in fig. 2, we find that the probability distribution of $\Delta t_{N1}$ decays exponentially as

$$\begin{equation} P(\Delta t_{N1})\propto \exp [{-(\Delta t_{N1})/\tau_{N1}}], \end{equation} \tag{ 8 }$$

and the characteristic time constant $\tau_{N1}$ equals the inverse slope of the probability distribution of $\Delta t_{N1}$ on a semi-log plot (fig. 2). The standard deviation of the neuronal voltage fluctuations σ can be thus calculated by using the characteristic time constant τ from all sleep bout durations: $\sigma=\Delta/\sqrt{\tau}$ with $\Delta=25\ \text{mV}$ for children and $\Delta=23\ \text{mV}$ for adults (for details on how to derive σ and Δ, see [16]).

Zoom In Zoom Out Reset image size

The probability distribution of $\Delta t_{N1}$ decays exponentially with characteristic time constant $\tau_{N1}$ (fig. 2 and eq. (8)) that can be calculated using the maximum likelihood estimation (MLE) [12,39]. An advantage of MLE analysis is its independence from a specific bin selection of the probability distribution.

The exponential probability distribution of N1 sleep bout durations can be approximated as $f(\Delta t_{N1})={1}/{\tau_{N1}} \cdot\exp [{-(\Delta t_{N1}-M)/\tau_{N1}}]$, where $\Delta t_{N1}$ is the variable of the N1 sleep bout durations, and M is the shortest value of $\Delta t_{N1}$. The likelihood ℓ for the data (with N variables) is

$$\begin{eqnarray*} \ell&= \prod_{i=1}^{N} \left(\frac{1}{\tau_{N1}} \cdot \exp [{-(\Delta t_{N1}^{i}-M)/\tau_{N1}}]\right),\\ \ell&= \left(\frac{1}{\tau_{N1}}\right)^N \cdot \exp \left[{\sum_{i} -(\Delta t_{N1}^{i}-M)/\tau_{N1}}\right],\\ \ln{\ell}&= -N \cdot \ln{\tau_{N1}} -\frac{ \sum_{i} (\Delta t_{N1}^{i}-M)}{\tau_{N1}}, \end{eqnarray*}$$

and the maximum of ℓ is attained when $\text{d}(\ln{\ell})/\text{d}\tau_{N1}=0$:

$$\begin{equation*} \frac{\text{d}(\ln{\ell})}{\text{d}\tau_{N1}}=0=-\frac{N}{\tau_{N1}} + \frac{ \sum_{i} (\Delta t_{N1}^{i}-M)}{(\tau_{N1})^2}. \end{equation*}$$

Therefore, the estimated value for $\tau_{N1}$ is (as expected for exponential distributions)

$$\begin{equation} \widehat{\tau}_{N1}=\frac{ \sum_{i} (\Delta t_{N1}^{i}-M)}{N}=\langle \Delta t_{N1}-M \rangle. \end{equation} \tag{ 9 }$$

Here, $\widehat{\tau}_{N1}$ was calculated from the $\Delta t_{N1}$ data in the interval [1 min, 10 min], therefore, $M=1\ \text{min}$ in eq. (9).

Using the Cramér-Rao bound [40], the estimated error of $\tau_{N1}$ is

$$\begin{eqnarray*} J_{{\tau}_{N1}}&= -\left\langle \frac{\ \text{d}^2 (\ln \ell)}{\ \text{d} (\tau_{N1})^2} \right\rangle\\ &= -\left\langle\frac{N}{(\tau_{N1})^2} - 2 \cdot \frac{ \sum_{i} (\Delta t_{N1}^{i}-M)}{(\tau_{N1})^3} \right\rangle\\ &= -\frac{N}{(\tau_{N1})^2}+ 2 \cdot \frac{N \cdot \tau_{N1}}{(\tau_{N1})^3} = \frac{N}{(\tau_{N1})^2}\approx\frac{N}{(\widehat{\tau}_{N1})^2}, \end{eqnarray*}$$

where $J_{{\tau}_{N1}}$ is the Fisher information of $\tau_{N1}$. Therefore, the lower bound of the standard deviation of $\tau_{N1}$ is

$$\begin{equation} \sigma_{\tau_{N1}}=\sqrt{\frac{1}{J_{\tau_{N1}}}} = \frac{\widehat{\tau}_{N1}}{\sqrt{N}}. \end{equation} \tag{ 10 }$$

Figure 2 shows the cumulative probability distributions of $\Delta t_{N1}$ for the healthy children and for the OSA group on semi-log plots. The cumulative distribution function is defined as CDF $\equiv \int_{\Delta t_{N1}}^\infty f(r) \ \text{d}r$, where $f(\Delta t_{N1})$ is the probability density function for the occurrence of sleep stage N1 with duration $\Delta t_{N1}$. Both distributions decay approximately exponentially with Akaike weights [41] for an exponential distribution > 0.9 and for a power-law distribution < 0.1. We find that the characteristic time constant $\tau_{N1}$ tends to be larger for OSA children, however, the difference is not statistically significant when comparing the values of each individual (see inset table in fig. 2(a)). Our model simulations result in probability distributions for $\Delta t_{N1}$ similar to the empirically obtained distributions (see fig. 2(b)).

We next calculate the N1 light sleep threshold L from sleep bout distributions for each group according to eq. (7). Interestingly, although $\tau_{N1}$ values are not significantly different when comparing healthy to OSA children, the obtained values of L show significant differences as can be seen in fig. 3. We find that L is significantly larger for OSA children, particularly for severe cases of OSA (Obstructive Apnea Index (OAI)${}\ge10$ events/hour [34]). We relate the higher values of L for OSA to the significantly longer N1 durations seen in those children [28,29].

Zoom In Zoom Out Reset image size

Previous studies have found that N1 sleep duration increases under OSA [25,30]. For children, this N1 increase correlates with OSA severity, i.e., children with higher OAI tend to have more N1 sleep (fig. 4(a)) [27,29]. The most important consequences of OSA, however, are oxygen desaturation episodes and hypoxemia that cause a multitude of neuropathological changes and neuropsychological impairments [42]. These episodes are detected by a fingertip pulse oximeter usually in combination with a polysomnogram (PSG), and the average number of desaturations per hour is quantified by the oxygen desaturation index (ODI). Since ODI is sensitive to apnea depth and length, it can be used as a prognostic marker for OSA complications and mortality [43]. In our data of OSA children, ODI shows good correlation with OAI (fig. 4(b)), similar to what was found recently for OSA adults [44]. Interestingly, as we show in fig. 4(c), the N1 light sleep threshold L, a physiological parameter proposed here, can be used to evaluate OSA severity. Plotting L as a function of the patients OAI shows a strong correlation of 0.71 (fig. 4(c)) for moderate OSA $(5\le \text{OAI} <10)$ to severe cases of OSA ($\text{OAI} \ge 10$) [34]). This is an even higher correlation than what we find for N1 and ODI, indicating that L may be useful as a diagnostic marker of OSA. Moreover, calculating L only requires hypnograms and no additional device to measure oxygen desaturation (which is necessary to determine ODI).

Zoom In Zoom Out Reset image size

The treatment of OSA in children includes a combination of medications (such as topical nasal steroids), a surgical removal of tonsils and adenoids ("adenotonsillectomy"), therapy with a CPAP mask during sleep and/or oral appliances (such as dental devices or mouthpieces to keep the upper airway open) [45]. In our study, treatment strategies were chosen based on the children's clinical situation and each child's family decision, and mostly included CPAP and adenotonsillectomy.

In fig. 5 we show the effect of treatment on the N1 light sleep threshold L. Compared to children with untreated OSA, those who underwent treatment by CPAP or surgery have significantly lower values of L. In fact, their L values become similar to healthy children (fig. 5(a)). Moreover, we can reinforce this observation of the decrease of L with treatment on individual basis. For n = 3 children we have obtained longitudinal data before and after treatment (with either CPAP or surgery) and we calculated the corresponding L values. It is clearly seen that treatment of these children lowered L significantly by 30–40% (fig. 5(c)) and ODI by more than 50% (fig. 5(d)). However, these preliminary results need to be confirmed by follow-up studies on a larger study cohort.

Zoom In Zoom Out Reset image size

Traditionally, obstructive sleep apnea is diagnosed using polysomnography in combination with respiratory and/or cardiovascular data (heart rate [46] and blood pressure [47]), and conventional hypnograms alone are often considered ineffective in determining OSA severity [46]. In this letter, we introduce a method to estimate the OSA severity through a physiological parameter L that is derived solely from standard hypnograms (fig. 4). More specifically, we have analyzed sleep data from healthy children and children with OSA, and we show that N1 sleep bout durations $(\Delta t_{N1})$ can be described by an exponential distribution with higher characteristic time constant for OSA children. We suggest a random-walk model to simulate the exponential distribution of $\Delta t_{N1}$ by introducing a light sleep threshold L that distinguishes between more shallow and deeper light sleep. Using this model framework and the characteristics of a Wiener process, we can calculate L from the sleep data and find significantly higher L values for OSA compared to healthy children, which might explain the longer time OSA patients spend in light sleep N1. Moreover, we demonstrate that L correlates well with OSA severity, showing highest correlation values for children with the most severe Obstructive Apnea Index (OAI).

OSA is characterized by repeated episodes of upper airway obstruction during sleep, which can result in intermittent brain hypoxia and sleep disruptions. Epochs of intermittent hypoxia can cause cellular injury in various parts of the central nervous system [48]. Because remote brain regions (at the edge of major vascular supply) are more susceptible to hypoxia, the cerebral cortex is particularly vulnerable [49]. Therefore, in OSA patients the hypoxic damage may lead to cortical dysfunctions causing some degree of cognitive impairments in attention, learning, and motor performance, which may lead to a growth deficit and abnormal behavior in children [50]. Since most of these impairments are partially reversible with OSA treatment [51], early diagnosis is crucial, especially for young children. Our results suggest that the parameter L, which we propose and develop in this letter, could be used as a new diagnostic marker for OSA since it captures the sleep microarchitecture, and therefore complements ODI which measures breathing disturbances. The estimation of L can easily be incorporated into clinical studies since it is derived from standard hypnograms and does not require additional equipment.

Further studies are needed to evaluate the diagnostic significance of L on a larger cohort of OSA patients that also include different age groups of the adult population as well as different types of OSA treatment. Another important problem is to study the applicability of L for other sleep disorders, for example, central sleep apnea (CSA) [52] that, in contrast to OSA, is caused by absent respiratory drive from the brain stem [19]. Although more research is needed to find the exact physiological explanation of L, this parameter could already be used in sleep apnea diagnostics.

This research was supported by the Bar-Ilan University and the Planning and Budgeting Committee fellowship program of the council for higher education of Israel; Shulamit Aloni Fellowship for Advancing Women in Exact Sciences and Engineering, Ministry of Science and Technology, Israel (Grant No. 3-13276); the Colman-Soref Grant foundation Fellowship of the Council for Higher Education, Israel (Grant No. kra/colman/194); the Israel Science Foundation (Grant No. 1657/16); the German Israeli Foundation (Grant No. I-1372-303.7/2016).