Stochastic modeling of central apnea events in preterm infants

The University of Virginia (UVA) Institutional Review Board approved this study under a waiver of consent. Computerized apnea detection results were not available to the patient care team, and decisions about apnea management and discharge readiness were based on nursing documentation of events and standard NICU monitor alarm data.

We analyzed cardiorespiratory waveforms from 298 VLBW infants consecutively admitted to the UVA NICU between January 2009 and June 2011. Demographic information, including birth weight, gestational age at birth, birth date, and types and times of respiratory support, were collected from patient medical records. Table 1 shows the demographic characteristics of the infants. Median gestational age at birth was 27 weeks and median birth weight 1010 g. The number of infants discharged home includes those that were discharged with respiratory support. Figure 1 shows a histogram of the number of patients with data available at each PMA (gestational age plus age from birth), with distinction made by disposition at discharge.

Zoom In Zoom Out Reset image size

We also analyzed patients within the 8 d preceding discharge to home: we found 196 VLBW infants that had data available for analysis immediately preceding discharge to home. Table 1 shows the demographic characteristics of this sample of infants.

All patients in the UVA NICU have continuous cardiorespiratory monitoring. We collected signals from three electrocardiograph leads digitized at 240 Hz, the chest impedance pneumograph digitized at 60 Hz, and pulse oximetry digitized at 120 Hz from the GE bedside monitors using a central network server (Bedmaster Ex, Excel Medical, Jupiter, FL) behind the clinical firewall. We also collected heart rate, respiratory rate, and oxygen saturation vital signs derived by the monitor. Data were transmitted to our centralized computing and storage cluster and analyzed in parallel.

Of the total NICU stay for all patients (55 patient years) (Clark et al 2012), 51% of data were deemed unsuitable for analysis due to epochs where recording was interrupted. Data were analyzed in 16 min epochs, and we required continuous time series of 3 ECG leads, chest impedance, and heart rate from the physiological monitor to define epochs as suitable for analysis. Missing chest impedance and dropped data packets from the data acquisition system were the most common causes for excluding an epoch. Of the remaining data (26.5 patient years), 10% were excluded to avoid the times during which the infants were being treated with mechanical ventilation.

Apnea events were identified in this data set as part of a previous study by our group—the apnea detection algorithm has been validated with accuracy  >90% by comparison with clinical inspection of waveform time series (Lee et al 2012, Vergales et al 2013). Briefly, apnea was detected as low variance epochs in the chest impedance pneumograph after notch filtering in heart-clock time to eliminate cardiac artifact and high-pass filtering to remove movement artifact. Heartbeats were detected using a threshold-based algorithm (Pan and Tompkins 1985) as implemented by Clifford and co-workers (Tarassenko et al 2001, Li et al 2008). For the current analysis, apnea was defined as breathing cessation of at least 10 s with both associated bradycardia and desaturations. We consider this definition for apnea clinically relevant because the breathing cessation caused both a decrease in heart rate and a decline in oxygen level in the blood—this is a stricter definition than that given by Finer et al (2006) in order to assure clinical relevance. The thresholds for bradycardia (HR  <  100, within 50 s from the beginning of the apnea or 25 s from the end of the apnea) and desaturation (SpO₂  <  80%, within 55 s from the beginning of the apnea or 38 s from the end of the apnea) are based on detailed inspection of 932 central apnea events by clinical personnel, see the appendix of Lee et al (2012). This apnea detection method is based on chest movement and therefore identifies only central apnea events (obstructive events are identified as breathing). While this algorithm does not explicitly exclude epochs of periodic breathing (Barrington and Finer 1990, Poets and Southall 1991), breathing pauses in periodic breathing epochs are shorter than 10 s and seldom induce sufficient bradycardia and desaturation to meet our thresholds (Khan et al 2005, Wilkinson et al 2007).

The results of Lee et al (2012) show an acceptable false positive rate (5%) throughout a NICU stay. Near discharge, however, the rate of false positives may be higher due to the reduced apnea rate and increased noise due to movement. We filtered out false positive apnea events within 8 d of NICU discharge as detailed in the appendix, section Filtering out false positive apnea events near discharge.

We studied the duration and timing of 26 088 apnea events occurring from admission through discharge for the patient population in table 1. Figure 2(a) shows the cumulative distribution of the apnea durations (dots). For each patient we identified sequential apnea events and calculated the time between each—the inter-apnea epochs. Inter-apnea epochs that started or ended at missing data were excluded as these are periods where apnea may occur but not be detected. Figure 2(b) shows the cumulative distribution of observed inter-apnea epochs for all infants regardless of PMA (dots).

Zoom In Zoom Out Reset image size

We modeled apnea of prematurity as a Markov process based on the interpretation of apnea as the result of transitions through decreasingly stable breathing states. One type of Markov model is a Markov chain, which has a large or infinite number of states and transitions only between adjacent neighbors. These models have diffusion-like behavior and have been used to model the tail of inter-breath interval distributions in infants (Frey et al 1998). These models result in power law distributions that may be useful for fitting apnea durations (Frey et al 1998) but are not appropriate for fitting multimodal distributions (McManus et al 1988)—no single power law can fit the distribution of inter-apnea epochs in figure 2(b), thus no solvable diffusion or fractal model describes the statistics of apnea in VLBW infants.

In the second type of Markov model a small number of states are postulated, with transitions occurring randomly from each state to others. Generally transitions are allowed directly from every state to every other state. Models of this type are used in a wide range of fields, e.g. speech recognition (Lee and Jean 2013, Schafer and Jin 2014), disease tracking (Detilleux 2008, Robertson et al 2011), and ion channel recording (Colquhoun and Hawkes 1982). In ion channel models, states are defined for open, closed, and blocked channels and transition rates are determined empirically. Knowledge of these transition rates allows for determination of important behaviors, for example mean open lifetime, mean time between openings, and number of ion channel openings per cluster. This type of Markov modeling also has important applications in medicine (Sonnenberg and Beck 1993) and has been used to evaluate obstructive sleep apnea in adults (Kim et al 2009) and detect apnea of prematurity in short heart rate time series for preterm infants (Altuve et al 2011).

Parameters of Markov models were derived following methods from analysis of single ion channel recordings (Colquhoun and Hawkes 1977, Colquhoun and Hawkes 1982). A Markov model is specified by the number of distinct observables M, the number of states N, the observation symbol probability distribution B = b_j, the state transition probability distribution C or the related transition rate matrix K, and the initial state distribution P(0) (Rabiner 1989). We define M  =  2 distinct observables for neonatal breathing: breathing and apnea. Derivation of parameters N and b_j are given in the following section. Details on the remaining parameters are provided in the appendix.

The number of states N in a Markov model is chosen from study of the distribution of dwell times and transition rates (Sigworth and Sine 1987), in our case apnea durations and time between apnea events, the dots in figure 2. If a Markov model with discrete states and continuous time is correct, then the probability q that an infant will stay in any given state for time t with no transition to another state decreases with time exponentially, q(t)  =  exp(−t/τ) where τ is the expected duration in the state. The corresponding cumulative distribution, which in our case is the probability that the duration is less than t, is given by Q(t)  =  1  −  exp(–t/τ). Similarly, the probability that the infant will stay among m hidden breathing states is a superposition of exponentials, and the corresponding cumulative distribution for the time between apneas is

$$\begin{eqnarray}&&Q(t)=1-\underset{i=1}{\overset{m}{\sum}}\,{{\alpha}_{i}}\exp \left(-t/{{\tau}_{i}}\right).\end{eqnarray} \tag{ 1 }$$

Here τ_i and α_i are the mean of the distribution and relative magnitude for exponential mode i, respectively.

The cumulative distribution of apnea durations was modeled with a single exponential, equation (1) with α  =  m  =  1 (figure 2(a), solid curve), with parameter τ_A, the mean apnea duration, by maximizing the log likelihood of the parameter given the observed durations (Sigworth and Sine 1987). The mean apnea duration based on the model for apnea events greater than 10 s was τ  =  24.9 s, and a single exponential function provided an accurate fit (R²  =  0.99).

In contrast, multiple distributions are evident in the cumulative distribution of inter-apnea epochs, figure 2(b) dots. We fit the data with increasing number of exponentials until the goodness-of-fit R²  ⩾  99.9%. The goal was to minimize the complexity of the model (i.e. as few states as possible) while accurately fitting the data. We optimized the parameters using maximum likelihood and found that the sum of four exponentials, equation (1) with m  =  4 (figure 2(b), solid curve), adequately describes the data. The parameters are provided in table 2. Based on these results we defined the number of states in the Markov model for neonatal breathing to be N  =  5, including four breathing states and one apnea state, S  =  [B₁, B₂, B₃, B₄, A].

As stated above, the observables are that the infant is breathing or that the infant is not breathing, and the observation matrix B gives the probability of observing O_j given that the model is in state S_i. For the apnea model this is trivial; the probability of observing breathing is 1 in each breathing state and 0 in the apnea state,

$$\begin{eqnarray}B=\left[\begin{array}{c} 1 & 0 \\ 1 & 0 \\ 1 & 0 \\ 1 & 0 \\ 0 & 1 \end{array}\right].\end{eqnarray} \tag{ 2 }$$

In order to characterize changes in breathing control with post-natal development, we created Markov models at weekly PMA—for each, we used a 3 week epoch centered on the one of interest. Specifically, we optimized the transition rate matrix K of a Markov model for each PMA from 25 to 40 weeks, 16 models in all. At each PMA a time series was created for each infant with data at that PMA, where each time series was a number each second: '0' if the infant is breathing and '1' if the infant is apneic. It is possible to have multiple time series per patient due to missing data or excluded epochs of mechanical ventilation. Each model was developed using B from equation (2) and the initial state distribution P(0) was taken to be,

$$\begin{eqnarray}&&{{P}_{i}}(0)=\frac{{{\tau}_{i}}}{\underset{i}{\sum}\,{{\tau}_{i}}}\end{eqnarray} \tag{ 3 }$$

where τ_i were determined using the exponential distribution fits to the data in figure 2, and are listed in table 2. Note that the same observation matrix B and initial state distribution P(0) were used for each PMA and these were not allowed to vary during optimization. The transition rate matrix K for each PMA was determined using the optimization technique described in the appendix using the initial estimate from equation (A.4). The transition probability density C was calculated as the matrix exponential of the transition rate matrix K.

At each PMA 40 bootstrap runs were completed (Feng et al 1996, Harrell 2001). For each run the n patients at that PMA (see figure 1) were randomly resampled with replacement—for each run some patients were used multiple times and some patients not used. For each bootstrap run we calculated the transition rates using the same procedure described above and in the appendix. At each PMA, and for each transition rate, we calculated the 95% confidence interval by identifying the 2.5 and 97.5%-tiles of the 40 samples—for 40 samples, the 95% confidence interval is defined by the minimum and maximum values.

Inspection of the waveform data from bedside monitors revealed clusters of apnea separated by longer periods without apnea. Figure 3(a) shows a cluster of three clinically significant apneas within three minutes. Only the events labeled with black bars were identified as apnea by the algorithm (because, as stated above, the algorithm selects breathing cessations of at least 10 s accompanied by bradycardia and desaturation). Throughout, we use ABDX to indicate apnea events of at least X seconds accompanied by bradycardia and oxygen desaturation. The first apnea is an ABD20. This is followed by two subclinical apneas, or breathing cessations that do not cause critical slowing of heart rate and declining oxygen saturation in the blood. Subsequently there are two clinical apneas; the first is 30 s and the second is 25 s, an ABD30 and an ABD25. Figure 3(b) depicts the number of ABD10 events for a single preterm infant from 33 to 35 weeks PMA. Each half hour is represented by a vertical band, with the intensity of the band indicating the number of ABD events longer than 10 s during that time period (zero events per half hour is white, six is black). In this example apnea occurs at highly variable intervals; there are instances of multiple apneas within a half-hour as well as days without apnea. From day 21–24 there are no apnea events (the infant was not on mechanical ventilation during this time, and data were available more than 90% of the time). Subsequently, there is another cluster of apnea. These results suggested breathing states of varying stability, with apnea accessible from any state.

Zoom In Zoom Out Reset image size

We also identified breathing cessations of substantial duration, many times longer than the 30 s duration characterized as ALTE (Ramanathan et al 2001). Figure 3(c) shows an example of an epoch of abnormal breathing control lasting nearly two minutes with bradycardia and oxygen desaturation. The detected event is labeled with a horizontal black bar. This epoch comes from an infant born at 27 weeks with a birth weight of 850 g who was eventually discharged to home from the NICU. This event was detected as an ABD110 by the automated detection algorithm, starting at 20 s and ending at 131 s. It is important to note that, although this episode represents a period of apnea detected as an ABD, isolated tiny fluctuations in the impedance signal may occur in long ABDs. This may permit these episodes to reach substantial lengths (Mohr et al 2014).

We studied Markov models with the state diagram shown in figure 4(a), with one apneic state (A) and four breathing states (B_i, i  =  1, 2, 3, 4) numbered from most to least stable. For each PMA we optimized the transition rate matrix K of a Markov model. Figures 4(b) and (c) show the dominant transitions (arrows) and states (bold) for 27 and 40 weeks PMA, respectively. The important findings are that the earlier PMA has higher probabilities of transitions towards and residence in less stable breathing and apnea states. The model at higher PMA is characterized by larger probability of the second-most stable breathing state, and lower probability of apnea.

Zoom In Zoom Out Reset image size

The Markov model transition rates are shown in figure 5 as a function of PMA. We found the average durations of these breathing states to be on the order of 10 s (B₄), 10 min (B₃), 2 h (B₂), and 12 h (B₁), see figure 5 and appendix section Calculating apnea burden and mean apnea duration from a Markov model. Confidence intervals are shown for the transition rates and percent residence for each state.

Zoom In Zoom Out Reset image size

The suitability of the model can be assessed by comparing its predictions to the observed data. We tested the goodness-of-fit of the model by comparing predicted and observed apnea burdens and the average duration of individual apneas. For each PMA, we calculated the apnea burden based on observed data as the number of seconds of apnea divided by the total monitoring time in days. We repeated the process for each Markov model being careful to censor apneas shorter than 10 s first, as apnea shorter than 10 s is not clinically significant (Finer et al 2006). Details are described in the appendix sections Calculating apnea burden and mean apnea duration from a Markov model and Censoring short apnea from the model. Figure 6 shows very good agreement between the predictions of the model and the observed data for apnea duration (figure 6(a)) and apnea burden (figure 6(b)). Observed data are shown as mean (circles) and 95% confidence interval (error bars) for the sample of patients at each PMA. This agreement is evidence that a Markov model with five (four breathing and one apnea) states and 2 observables (breathing and not breathing) can accurately characterize the kinetics of apnea in VLBW infants.

Zoom In Zoom Out Reset image size

Rare but clinically significant apnea events occur after discharge from the NICU. These ALTEs often lead to rehospitalization. We wish to estimate the rate of ALTE after NICU discharge. Thus we constructed a Markov model based on data from the eight days before discharge home, see table 1. We optimized the transition rate matrix K of a Markov model using the same methods described above: the mean residence time in each Markov state in this model are 18.1 h (B₁), 2.2 h (B₂), 13.9 min (B₃), and 12.3 s (B₄).

We used the optimal transition rate constants from the discharge Markov model to estimate the rate of prolonged apnea after discharge. An analytical expression to estimate the rate of extreme apnea after NICU discharge is derived in the appendix section Estimating the rate of extreme apneas after discharge. This estimate is shown in figure 7—the figure addresses the question, 'given that an infant is discharged from the NICU to home, what is the probability that that infant will go X days and have at least one apnea of duration at least Y '. The rate of an apnea event is shown as a function of days after NICU discharge on the abscissa and duration on the ordinate. For example, our model estimates that the rate at which infants discharged to home have at least one apnea  >128 s within the first 6 months after NICU discharge is 1 in 100. This extrapolation predicts that the rate of severe apnea (~180 s) occurring after discharge is on the order of 1 per 1000. This is on the order of the occurrence rate of SIDS in the US population (Moon and the Task Force on Sudden Infant Death Syndrome 2011).

Zoom In Zoom Out Reset image size

A limitation of our study is the inconsistent quality of neonatal breathing measurements—a large proportion of patient data were missing or not of suitable quality for analysis. More accurate and robust monitoring would improve the estimate for severe apnea occurring after NICU discharge, and would allow the stability of neonatal breathing control to be more accurately assessed. Another limitation is that additional filtering was performed only for apnea events within 8 d of discharge. Filtering all 26 088 would likely increase the transition rates out of apnea in figure 5 and reduce the residual apnea burden and duration at higher PMA in figure 6. Another limitation is that, were data excluded due to movement, this may have led to a non-random distribution of missing data. If the more active infant is less likely to experience apnea, this would tend to skew our results towards shorter time periods between apneas. Whatever the imperfections of our algorithm, we emphasize that this continuous-monitoring system is far more accurate than nursing records.

Another limitation of our approach is that central apnea events are identified without distinction by cause. Some apneas near discharge may reflect poor sucking-swallowing-breathing coordination during feeds, or changing sleep architecture. Other normal physiologic processes contribute to apnea, notably anemia (Zagol et al 2012) and gastro-esophageal reflux (Poets 2004). Both of these normal physiologic phenomena may be progressing at the time of NICU discharge and may play a role in the stochastic distribution of apnea, with the possibility of very long time between apnea events as the chemoreceptors and other elements of the respiratory control system mature. Another limitation of our study is that we have considered only one of many manifestations of pathologic cardiorespiratory control. Future models for transitions between stable and abnormal cardiorespiratory control should include not only apnea (Martin et al 2002, Darnall 2009, Di Fiore et al 2013), but also, for example, periodic breathing (Khan et al 2005, Wilkinson et al 2007, Mohr et al 2015), cardiorespiratory uncoupling (Godin and Buchman 1996, Clark et al 2012), and transient heart rate decelerations (Dorostkar et al 2005, Moorman et al 2011). A final limitation is that we model the rate of prolonged apnea after discharge using data from the last 8 d in the NICU, thus extrapolating the data as if neurological development were complete at discharge.

A major strength of this study is the use of computer algorithm-detected apnea based on chest impedance waveforms rather than nursing documentation or a tally of bedside monitor alarms, thus allowing us to accurately quantify occurrence and duration of events (Southall et al 1983, Muttitt et al 1988, Vergales et al 2013). Automated detection algorithms could be continuously implemented at the bedside of preterm infants in intensive care, where all such infants already have routine cardiorespiratory monitoring. Such implementation would allow online estimation of the risk for apparently random transitions to isolated and severe apnea events following extended periods without apnea.

Hidden Markov models represent time series data as random transitions between states that are not visible to the observer. At each instant the probability that an infant is in any one of N states is described by a vector with N components, P(t)  =  (P₁, P₂, ..., P_N). In the context of neonatal apnea, there are states that represent apnea and others that represent breathing states with varying levels of stability that are hidden from observers. The change of that probability vector is given by,

$$\begin{eqnarray}&&\mathbf{P}(t+1)\,=\mathbf{P}(t)\mathbf{C}\end{eqnarray} \tag{ A.1 }$$

where C is a matrix that represents transition probabilities: the probability P that the system is in state j at time $t+1$ given that it was in state i at time $t$ is

$$\begin{eqnarray}&&{{\mathbf{C}}_{ij}}=P\left(S(t+1)=j|S(t)=i\right)\end{eqnarray} \tag{ A.2 }$$

where S(t) represents the state of the system at time $t$. In a stationary Markov model, this transition probability matrix does not depend on time. After each step, there is a finite probability that the infant remains in the same state (i  =  j), and for each initial state ${{S}_{i}}$

$$\begin{eqnarray}&&\underset{j=1}{\overset{N}{\sum}}\,{{C}_{ij}}=1.\end{eqnarray} \tag{ A.3 }$$

The general theory of hidden Markov models includes also a collection of observables (O_k, k  =  1, 2, ...), and there is another matrix $\mathbf{B}$ that converts the state probability vector P(t) into the probability of observing each particular outcome (O_k). In our case this matrix is trivial: the two outcomes are breathing or not-breathing. If an infant is in any of the breathing states the probability of observing breathing is 1, while if the infant is in an apnea state the probability of observing breathing is 0. Thus, 'hidden' refers to the fact that we cannot identify the stability of an infant's respiratory control system solely by observing that the infant is breathing.

Several algorithms are available for obtaining maximum likelihood estimates of the one-step transition probability matrix C and from that its matrix logarithm K. Here K is a time-invariant N  ×  N matrix of rate coefficients. The diagonal elements of K are negative, all others are positive, and each row of the matrix sums to zero. We used the Baum–Welch maximization algorithm to determine the optimal transition rate matrix K (Baum et al 1970, Juang and Rabiner 1991) as implemented in MatLab^® (The MathWorks, Natick, MA), though the Viterbi algorithm is also an option (Shinghal and Toussaint 1979). These optimization algorithms allow for non-zero transition rates between any pair of states. This method requires an initial estimate of C and of the initial vector of probabilities P(0). The initial estimate C₀ was derived from an initial estimate K₀ that was chosen to be,

$$\begin{eqnarray}{{K}_{0}}=\left[\begin{array}{c} -{}^{1}/{}_{{{\tau}_{1}}} & {}^{1}/{}_{{{\tau}_{1}}} & 0 & 0 & 0 \\ {}^{1}/{}_{2{{\tau}_{2}}} & -{}^{1}/{}_{{{\tau}_{2}}} & {}^{1}/{}_{2{{\tau}_{2}}} & 0 & 0 \\ 0 & {}^{1}/{}_{2{{\tau}_{3}}} & -{}^{1}/{}_{{{\tau}_{3}}} & {}^{1}/{}_{2{{\tau}_{3}}} & 0 \\ 0 & 0 & {}^{1}/{}_{2{{\tau}_{4}}} & -{}^{1}/{}_{{{\tau}_{4}}} & {}^{1}/{}_{2{{\tau}_{4}}} \\ 0 & 0 & 0 & {}^{1}/{}_{{{\tau}_{A}}} & -{}^{1}/{}_{{{\tau}_{A}}} \end{array}\right].\end{eqnarray} \tag{ A.4 }$$

The transition rates were determined as a function of PMA using these parameters. The resulting coefficients of the transition rate matrix K are shown in figures 5(a)–(e). Note that the diagonal elements of K are negative to maintain a zero row sum. The coefficients are shown for the longest lifetime breathing state B₁ (a), the next longest lifetime breathing state B₂ (b), the third longest lifetime breathing B₃ (c), the shortest lifetime breathing state B₄ (d), and the apneic state A (e). The percentage of time spent in each state is shown in figure 5(f).

To check the consistency between the model and the data, we estimate the apnea burden and the mean apnea duration. At steady state, the probability of apnea is the value of the corresponding element of the normalized left-eigenvector of the transition rate matrix K having eigenvalue equal to zero. If this probability is multiplied by the number of seconds in a day, it gives the apnea burden (seconds of apnea per days of monitoring) predicted by the Markov model. The average apnea duration predicted by the model is  −1/K_AA, where K_AA is the transition rate from the apnea state into the apnea state. For any state S_i the average residence time in that state is  −1/K_ii.

A Markov model assumes the dwell time in each state is exponentially distributed. Apnea durations are censored below 10 s. Censoring an exponential probability distribution below a value x has the effect of increasing the mean of that distribution by x. The average residence time in the apnea state (A) must be decreased by 10 s to account for this fact. This is achieved by modifying transition rates K_Ai, i  =  1...N representing transitions out of apnea as,

$$\begin{eqnarray}&&K_{\text{A}i}^{\prime}=\frac{{{K}_{\text{A}i}}}{1-10|{{K}_{\text{A}i}}|}.\end{eqnarray} \tag{ A.5 }$$

When modeling apnea during the 8 d prior to discharge home, the apnea detection algorithm produced a larger false positive rate than that during the rest of the stay. This is a result of the lower apnea rate and increased infant activity at this point in neonatal care. In order to reduce the rate of false positives and more accurately estimate the rate of apnea, we filtered events using an automated algorithm. For each apnea event we calculated two parameters: the area under the apnea probability (Lee et al 2012) divided by the apnea duration, and the power in the chest impedance between 0.5 Hz and 2.0 Hz. These two parameters were calculated for 500 automatically-identified apnea events randomly selected throughout the entire stay of all infants in the study. Of these, 471 events were adjudicated as true apnea or false positive by agreement of two neonatologists. A logistic regression model was created using maximum likelihood estimation. The receiver operating characteristic (ROC) area of this model is 0.93.

We used this model as an additional filter for apnea events, and determined a threshold below which apnea events were censored from our analysis. This threshold was determined by randomly selecting 100 events from the 8 d prior to discharge in 196 infants, of which 96 events were adjudicated. We selected a threshold for the apnea filter such that the false positive rate was equal to the false negative rate: we used an apnea filter threshold of 0.81.

We also computed a K matrix based on time series from the last 8 d prior to discharge home. We extrapolated the rate of apnea occurring after discharge from the NICU. We wish to estimate the probability that an infant will have at least one apnea event within N days after discharge, and that at least one of those events will have duration greater than d. Figure A1 shows the distribution of apnea durations in the 8 d before discharge (circles) and the modeled distribution (dashed). The rate of apneas longer than about 1½ min is extrapolated.

Zoom In Zoom Out Reset image size

We now model the probability of having n apnea events in N days using a Poisson distribution. We consider each day as a trial, and estimate the probability of n 'successes' (apnea events) based on the expected number of successes v. In this context, the probability of having exactly n apneas in N days is the Poisson probability density distribution,

$$\begin{eqnarray}&&{{P}_{v}}(n)=\frac{{{v}^{n}}}{n!}{{\text{e}}^{-v}}.\end{eqnarray} \tag{ A.6 }$$

Here v, the expected number of apneas in N days, is

$$\begin{eqnarray}&&v=-\alpha N{{K}_{\text{AA}}}{{p}_{\text{A}}}.\end{eqnarray} \tag{ A.7 }$$

where α  =  86 400 s per day, K_AA is the diagonal element of the rate matrix corresponding to the apneic state (i.e. the negative inverse of the mean apnea duration), and p_A is the probability of apnea. At steady state, this probability of apnea p_A is the value of the corresponding element of the normalized left-eigenvector of the transition rate matrix K having eigenvalue equal to zero. The probability that an apnea is shorter than d seconds is the cumulative of the exponential distribution,

$$\begin{eqnarray}&&P(D < d)=1-{{\text{e}}^{d{{K}_{\text{AA}}}}}.\end{eqnarray} \tag{ A.8 }$$

The proposition that there are no apneas within N days of discharge, or that any apnea events occurring within N days are shorter than d, is,

$$\begin{eqnarray}&&\left(n=0\vee \left[n>0\wedge \forall D:D<d\right]\right)\end{eqnarray} \tag{ A.9 }$$

We use the distributions in equations (A.6) and (A.8) to find the probability of this proposition,

$$\begin{eqnarray}&&P\left(n=0\vee \left[n>0\wedge \forall D:D<d\right]\right)=\underset{n=0}{\overset{\infty}{\sum}}\,P(n)P{{(D < d)}^{n}}.\end{eqnarray} \tag{ A.10 }$$

This simplifies to,

$$\begin{eqnarray}&&P\left(n=0\vee \left[n>0\wedge \forall D:D<d\right]\right)={{\text{e}}^{-\beta}}.\end{eqnarray} \tag{ A.11 }$$

$$\begin{eqnarray}&&\beta =v{{\text{e}}^{d{{K}_{\text{AA}}}}}.\end{eqnarray} \tag{ A.12 }$$

The proposition of having at least one apnea with duration greater than d is the negation of the proposition in equation (A.9),

$$\begin{eqnarray}&&n>0\wedge \exists D:D\geqslant d\equiv \neg \left(n=0\vee \left[n>0\wedge \forall D:D<d\right]\right).\end{eqnarray} \tag{ A.13 }$$

The probability of prolonged apnea after discharge is therefore one minus the probability in equation (A.11),

$$\begin{eqnarray}&&P\left(n>0\wedge \exists D:D\geqslant d\right)=1-{{\text{e}}^{-\beta}}\end{eqnarray} \tag{ A.14 }$$

This result is shown in figure 7.