A Persistent Homology Approach to Heart Rate Variability Analysis With an Application to Sleep-Wake Classification

Persistent homology is a recently developed theory in the field of algebraic topology to study shapes of datasets. It is an effective data analysis tool that is robust to noise and has been widely applied. We demonstrate a general pipeline to apply persistent homology to study time series, particularly the instantaneous heart rate time series for the heart rate variability (HRV) analysis. The first step is capturing the shapes of time series from two different aspects—the persistent homologies and hence persistence diagrams of its sub-level set and Taken's lag map. Second, we propose a systematic and computationally efficient approach to summarize persistence diagrams, which we coined persistence statistics. To demonstrate our proposed method, we apply these tools to the HRV analysis and the sleep-wake, REM-NREM (rapid eyeball movement and non rapid eyeball movement) and sleep-REM-NREM classification problems. The proposed algorithm is evaluated on three different datasets via the cross-database validation scheme. The performance of our approach is better than the state-of-the-art algorithms, and the result is consistent throughout different datasets.


INTRODUCTION
Heart rate variability (HRV) is the physiological phenomenon of variation in the lengths of consecutive cardiac cycles, or the rhythm of heart rate (Draghici and Taylor, 2016). Interest in HRV has a long history (Billman, 2011), and there have been several theories describing how the heart rate rhythm, including, for example, the polyvagal theory (Porges, 2009) and the model of neurovisceral integration (Thayer and Sternberg, 2006). In short, HRV results from an integration of complicated interactions between various physiological systems and external stimuli (Vanderlei et al., 2009;Shaffer et al., 2014;Draghici and Taylor, 2016) on various scales, and the autonomic nervous system (ANS) plays a critical role (Thayer and Sternberg, 2006;Porges, 2009). A correct quantification of HRV yields dynamical information of various physiological systems and has various clinical applications (Stys and Stys, 1998), including improving diagnostic accuracy and treatment quality (Vanderlei et al., 2009).
In practice, the heart rhythm is quantified by the time series called instantaneous heart rate (IHR) coming from intervals between consecutive pairs of heart beats, which is usually determined from the R peak to R peak interval (RRI) by reading the electrocardiogram (ECG). See Figure 1 for an illustration of the ECG, R peaks, and RRI. To quantify HRV, a common approach is studying various statistics of IHR. There have been a lot of efforts trying to quantify HRV, and proposed statistics could be briefly classified into four major categories-time domain approach, frequency domain approach (Task Force of the European Society of Cardiology and others, 1996), nonlinear geometric approach (Marwan et al., 2002;Voss et al., 2008), and information theory based approach (Costa et al., 2002). It is worth mentioning that while there have been a lot of researches in this direction with several proposed statistics, there is limited consensus and it is still an active research field due to the non-stationarity nature of the IHR time series (Pincus and Goldberger, 1994;Glass, 2009).
Topological data analysis (TDA) is a data analysis framework based on tools from algebraic topology (Carlsson, 2009;Epstein et al., 2011). In the past decades, its theoretical foundation has been actively established, and various algorithms have been proposed to study datasets from different fields. The basic idea underlying TDA is that the data organization can be well-captured by counting holes. Theoretically, the number of holes of different dimensions characterizes how the data is organized. Thus, researchers design useful statistics based on the information of holes. This simple yet powerful idea has been applied to different fields. Specifically, there have been several efforts applying TDA to analyze time series. For example, the Vietoris-Rips (VR) complex filtration and the bottleneck or Wasserstein distances among persistence diagrams are applied to study voices and body motions (Seversky et al., 2016;Venkataraman et al., 2016). A transformation of the persistence diagram, called persistence landscapes (Bubenik, 2015), has been applied to study trading records (Gidea and Katz, 2018), electroencephalogram (EEG) signals (Piangerelli et al., 2018;Wang et al., 2018;Wang et al., 2019), and cryptocurrency trend forecasting (Kim et al., 2018). Sliding Windows and 1-Persistence Scoring (Perea, 2019) offers both theoretical and practical TDA method to detect the periodicity of a time series. A brief overview of common techniques on the usages of TDA to time series is summarized in a preprint (Ravishanker and Chen, 2019).
In this article, motivated by the complicated interaction among different physiological systems over various scales and inter-individual variability, the need for a useful tool for the HRV analysis, and the numerical limitation of the recently developed TDA tools, we hypothesize that topological information could be useful to quantify the HRV, and propose a computationally efficient approach to analyze time series via TDA.

Our Contribution
Based on the flexibility of TDA tools, and due to the nonstationarity of complicated time series we commonly encounter in real life, like the IHR, we propose a systematic, principled, and computationally efficient approach to study complicated time series by the TDA tools.
Our main scheme for studying a complicated time series is shown in Figure 2, which is divided into three steps that we will detail later. First, consider two filtrations, the Vietoris-Rips (VR) complex filtration of the Takens' lag map (Takens, 1981) and the sub-level set filtration of the time series, and persistent homology. Second, compute corresponding persistence diagrams. Finally, calculate persistence statistics (PS) as a novel statistic of the time series of interest. We mention that compared with existing TDA approach for time series analysis, our proposed persistence statistics features based on both sub-level set and Vietoris-Rips complexes filtrations are intuitive, straightforward to implement, and also computationally efficient.

Application-Sleep Dynamics
To demonstrate the usefulness of the proposed persistence statistics, we apply it to study IHR time series recorded during sleep, and use obtained statistics to classify sleep stages. Sleep is a universal recurrent physiological phenomenon. Sleep impacts the whole body, so we can read sleep via reading different physiological signals. Taking ECG into account is specifically attractive, since the ECG sensor is easy to install, and it is now widely available in mobile health devices. HRV of a subject is usually quantified by analyzing ECG, and it has been shown to be related to sleep dynamics (Zemaityte and Varoneckas, 1984;Vaughn et al., 1995;Toscani et al., 1996;Bonnet and Arand, 1997;Elsenbruch et al., 1999;Chouchou and Desseilles, 2014;Penzel et al., 2016). In other words, the heart rate rhythm provides a non-invasive window for researchers to study sleep. While there have been several studies trying to classify sleep stages based on HRV (Lewicke et al., 2008;Mendez and Matteucci, 2010;Long et al., 2012;Xiao et al., 2013;Aktaruzzaman et al., 2015;Ye et al., 2016;Malik et al., 2018), it still remains a challenging problem in the field. The challenge and difficulty of this mission can be appreciated from the reported results. In this article, we apply the proposed persistence statistics to quantify HRV during sleep, and propose a new prediction algorithm for the sleep stage; for example, an automatic classification of wake and sleep, REM and NREM, and wake, REM and NREM. We remark that while we focus on the HRV and sleep stage classification, the result FIGURE 1 | An illustration of ECG (black curve), R peaks (red circles), and RRI (the numbers between two consecutive R peaks). It is clear that the RRI changes from time to time, which form a new time series.
FIGURE 2 | The scheme of our proposed time series analysis can be separated by three steps: constructing filtrations, computing persistence diagrams, and extracting persistence statistics as features. The features are applied to train a machine learning model for the classification purpose.
indicates the potential of applying TDA-based approaches to study other complicated time series.

Organization
In section 2, we review the mathematical background of the persistent homology and persistence diagram. In section 3, we demonstrate two ways to use the persistent homology to study time series, and propose a new approach to summarize the persistence diagram, called the persistence statistics. The classification model based on the persistence statistics for the sleep stage classification will be discussed in detail in section 4. The discussion of our classification performance and a comparison with the state-of-arts results will be included in section 5. More technical details and numerical results are postponed to the Supplementary Material.

MATHEMATICAL BACKGROUND
In this section, we describe the mathematical background, including simplicial complex, homology, filtration of sets, and the persistent homology. Although these topics can be studied in an abstract and general way (see e.g., Munkres, 1993), to enhance the readability, we present them in a relatively concrete way without losing critical information.

Simplicial Complexes
To investigate the complicated structure of an object, an intuitive way is to use simple objects as building blocks to approximate the original object. In TDA, the main building block is the simplicial complex, which we briefly recall now. See Supplementary Material (section 1.1) for more detailed mathematical background and illustrative examples.
We start with the simplex. Intuitively, a simplex is a "triangle" of different dimension. Let x 0 , x 1 , . . . , x q be affinely independent points in R d , where d, q ∈ N and d ≥ q. The q-simplex, denoted by σ := x 0 , x 1 , . . . , x q , is defined to be the convex hull of x 0 , x 1 , . . . , x q . Denote Vert(σ ) := {x 0 , x 1 , . . . , x q }. Any q-simplex is a q-dimensional object consisting of lower degree simplexes. We are interested in the relation among simplexes of different dimensions. Since any V ⊂ Vert(σ ) is also affinely independent, the convex hull of V, called a face of σ , forms a simplex of dimension |V| ≤ q, where |V| is the cardinality of V. If |V| = k, the face τ = V is called a k-face of σ . A simplicial complex K in R d is a collection of finite simplexes σ in R d so that any intersection of two arbitrary simplexes is a face to each of them; that is, • If σ ∈ K and τ is a face of σ , then τ ∈ K; • If σ 1 , σ 2 ∈ K, then σ 1 ∩ σ 2 is a face of σ 1 and σ 2 . In particular, σ 1 ∩ σ 2 ∈ K.

Homology and Betti Numbers
In order to study the topological information of a given simplicial complex, we study relations among simplexes of different dimensions, and hence the "holes." Homology and Betti numbers are classic subjects in the algebraic topology (Munkres, 1993), which capture "holes" of geometric objects of different dimensions. While we can discuss these topics in a more general setup, in this work, we mainly consider simplicial complexes as our target object. For example, the shape of the notation "∞" contains two 1-dimensional holes (Supplementary Figure 2C) and the empty void surrounded by the unit sphere S 2 = {(x, y, z) : x 2 + y 2 + z 2 = 1} in R 2 is a 2-dimensional hole (Supplementary Figure 2E). Moreover, the 0-dimensional holes of an object are defined to be its disjoint connected components (Supplementary Figure 2E). See Supplementary Material (section 1.2) for more information and illustrative examples. We need an algebraic structure of simplexes. Given qsimplexes σ 1 , σ 2 , . . . , σ n in a simplicial complex K, define the sum over Z 2 as c = n i=1 ν i σ i , where ν i ∈ Z 2 . This formal sum is commonly known as a q-chain. One could also define an addition operator as We consider the collection of all q-chains, denoted as One could prove that C q (K) is actually a vector space over Z 2 with the above addition. There is a natural relation between C q (K) and C q−1 (K), called the boundary map (Munkres, 1993, section 1.5, p. 30). The qth boundary map ∂ q : C q (K) → C q−1 (K) over Z 2 is defined by x 0 , · · · , x i , · · · x q , where x 0 , x 1 , · · · , x q ∈ K and the • denotes the drop-out operation. With the boundary maps, there is a nested relation among chains This nested relation among chains is known as the chain complex, which is denoted as C = {C q , ∂ q } q∈Z .
A fundamental result in the homology theory (Munkres, 1993 Lemma 5.3 section 1.5, p. 30) is that the composition of any two consecutive boundary maps is trivial, i.e., ∂ q−1 • ∂ q = 0. This result allows one to define the following quotient space. Denote cycles and boundaries by Z q (K) and B q (K), respectively, which are defined as Note that each B q (K) is a subspace of Z q (K) since ∂ q−1 • ∂ q = 0. Therefore, we can define the qth homology to be the quotient space which is again a vector space. For instance, if K = K 3 in Figure 3, then H 0 (K 3 ) ≃ Z 2 and H 1 (K 3 ) ≃ Z 2 because it contains one connected component and one 1-dimensional hole. More precisely, the 1-dimensional hole in K 3 is represented by the 1-cycle Actually, a q-hole may be represented by different q-cycles. For example, if K = K 3 in Figure 3, then the 1-cycles is the boundary of the 2simplex v 1 , v 2 , v 4 . This gives us an intuition about the algebraic structure (1). The qth Betti number is defined to be the dimension of the qth homology; that is, which measures the number of q-dimensional holes. As a result, given a simplicial complex K, the homology of K is a collection of vector spaces {H q (K)} ∞ q=0 , and its Betti numbers is denoted as

Persistent Homology
We now introduce a natural generalization of homology, the persistent homology, that is suitable for data analysis. persistent homology is more suitable for data analysis than homology due to this capability of dealing with inevitable noise in real world dataset. It depends on the notion of filtration to handle noise. In general, filtration is a sequence of simplicial complexes (see Figure 3 for an example). We are interested in how the "holes" vary in the filtration. Intuitively, if certain holes are "robust, " they will survive in the filtration. For an index set I, a filtration is a sequence of simplicial complexes, {K t } t∈I , satisfying From the previous discussion, for each K t in a filtration, one could compute its homology group and Betti number. Because of the nested subset relation in a filtration, there exist relations among simplicial complexes. This allows one to track and record the changes of the homology group and the Betti numbers, which we detail now. Given a fixed q ≥ 0, each K i induces homology H q (K i ). Denote ι i : K i → K i+1 to be the inclusion map. Then ι i (Z q (K i )) ⊆ Z q (K i+1 ) and ι i (B q (K i )) ⊆ B q (K i+1 ) (Edelsbrunner and Harer, 2010 section 4.4, p. 95). Therefore, the mapping induced by ι i is a well-defined linear transformation over Z 2 . We also define a linear transformation which maps H q (K i ) to H q (K i+k ). The following definition is crucial for defining lifespans of connected components or holes in homology theory.

Definition 2 (Edelsbrunner and Harer, 2010 section 7.1, p. 151).
Let {K i } n i=0 be a filtration of simplicial complexes. For q ∈ Z ≥0 and i, j ∈ Z ≥0 with i ≤ j, we define the persistent homology as Since K 0 ⊆ K 1 ⊆ · · · ⊆ K n , we have inclusions of q-chains: C q (K 0 ) ⊆ C q (K 1 ) ⊆ · · · ⊆ C q (K n ) for all q ≥ 0. Hence, the intersection B q (K j ) ∩ Z q (K i ) is a well-defined subspace of Z q (K i ). Moreover, for i ≤ j, the kernel of the linear transformation By the first isomorphism theorem, we obtain an injective linear transformation Via the one-to-one linear mapping φ , which means that the persistent homology is a generalization of the homology. With the inclusion H i,j q ֒− → H q (K j ), we define the moments of birth and death of a "hole" in the filtration.
The death d = ∞ means that the q-hole never dies in the filtration.
We use Figure 3 to explain the relation between these two abstract definitions. For instance, the non-trivial element c represented by 1-chain 1 (c)}, thus c dies at K 5 . We refer readers with interest to Edelsbrunner and Harer (2010) for more details in persistent homology.

Persistence Diagram
Persistence diagram proposed in Edelsbrunner et al. (2000) or equivalently persistence barcodes proposed in Carlsson et al. (2005) is a tool to visualize the complicated lifespans of holes in a given filtration for data analysis. We use persistence diagram in this paper.
The persistence diagram possesses the desired stability property (Cohen-Steiner et al., 2007)-a bounded perturbation of a given filtration leads to a bounded perturbation of the corresponding persistence diagram. Due to the inevitable noise in real data, this stability property renders persistence diagram based approaches suitable for data analysis. The bottleneck and Wasserstein distances (Cohen-Steiner et al., 2007) are typical ways to measure differences among persistence diagrams. The formal statements of the stability property based on these two distances are provided in sections 3.1, 3.2. We refer readers with interest to Edelsbrunner and Harer (2010) for details in persistence diagram.

Definition 4 (Edelsbrunner and Harer, 2010 section 7.1, p. 152).
Let {K i } n i=0 be a filtration of simplicial complexes and q ∈ Z ≥0 . The qth persistence diagram, denoted as P q ({K i } n i=0 ), of the filtration is the multiset of q-dimension holes in the filtration.
In other words, a q-dimensional hole in a filtration is recorded by a pair (b, d) of integers where b and d are called the birth and death of the hole, respectively . Although the above definition of persistence diagram seems technical, its interpretation is intuitive. For instance, consider the filtration shown in Figure 3. We look for the "changes" of topological structure (holes). Note that since a connected component is born at K 1 (specifically, v 1 , v 2 ), its birth value is b = 1; since it lives throughout the filtration, its death value is ∞. We now turn our focus to the 1-dimensional hole. Note that a 1-dimensional hole (specifically, v 1 , v 2 + v 2 , v 3 + v 1 , v 3 ) is formed at K 3 , so its birth value is 3; note also that this hole is filled at K 5 , so its death value is 5. Since there is no more change of holes, we have the persistence diagrams P 0 ( Before closing this subsection, we illustrate how persistent homology and persistence diagram work by taking a noisy point cloud sampled from a circle contaminated by Gaussian noise shown in Figure 4A. If there is no noise, the 1st Betti number of the circle is β 1 = 1. In the noisy case, the Betti number information is contained in the form of the persistence diagram as shown in Figure 4B, where each point represents one 1dimensional hole associated with its birth and death value. In Figure 4B, we observe that there is an outstanding point with long lifespan (located around birth value 0.05 and death value 0.25), while lifespans for other points are very small. This suggests that the noisy point cloud has a strong/robust 1-dimensional hole. This captures the main topology information, β 1 = 1, about this data.

Data Analysis With Persistence Diagram and Commonly Considered TDA Statistics
Usually, researchers design statistics on the persistence diagram of a given dataset via the chosen filtration. One basic result supporting this approach is Mileyko et al. (2011), where authors showed that the space of persistence diagrams with certain metric is complete and separable. This result forms a theoretic foundation for any statistical methods. In Fasy et al. (2014) and Blumberg et al. (2014), authors derived confidence sets of persistence diagrams in order to separate the long lifespan holes from noisy ones, and also proposed four ways to estimated them. While these theoretical results shed light on applying TDA to analyze complex data, however, any operation in the space of persistence diagrams is complicated and difficult to compute. For example, computing bottleneck or Wasserstein distances among persistence diagrams is a difficult task and can be time consuming, even for the state-of-art algorithm (Kerber et al., 2017). Another result indicates that the mean in the space of persistence diagrams may not be unique (Turner et al., 2014a). This computational burden renders it less applicable to data analysis.
To get around the computational issue when working with those distances, one major approach is to "vectorize" persistence diagrams; that is, researchers map the space of persistence diagrams into another space. For example, persistence landscapes (Bubenik, 2015) map persistence diagrams into a Banach space, specifically L p space. More examples include persistence image (Adams et al., 2017), generalized persistence landscapes (Berry et al., 2020), persistence path (Chevyrev et al., 2018), persistence codebook (Zelinski et al., 2020), persistence curves (Chung and Lawson, 2019), kernel based methods (Reininghaus et al., 2015;Kusano et al., 2016), and persistent entropy (Chintakunta et al., 2015;Atienza et al., 2019b). These methods have been studied and applied to different applications. In Figure 5, we provide a chart depicting the relationship among existing TDA tools. We mention that the proposed persistence statistics in section 3 could be viewed as a computationally efficient vectorization of persistence diagrams.

TDA FOR TIME SERIES ANALYSIS AND FEATURES EXTRACTION
Armed with the theoretical background in section 2, we are ready to describe how to apply TDA for time series analysis. To apply the persistent homology to analyze complicated time series, we introduce two useful filtrations, the sub-level set filtration and the Vietoris-Rips complexes filtration. With these two filtrations, we introduce a novel features extraction methods, coined persistence statistics, based on the persistence diagrams of the sub-level set filtration and the Vietoris-Rips complexes filtration.

First Useful Filtration-Sub-Level Set Filtration
To simplify the discussion and illustrate the idea, we identify a time series as a discretization of a continuous function f :[0, T] → R, where T is some fixed constant. For each h ∈ R, the sub-level set of f is defined as Clearly, f h 1 ⊆ f h 2 whenever h 1 ≤ h 2 . Therefore, for any increasing sequence {h i } i , the collection of sub-level sets, {f h i }, forms a filtration. Intuitively, the sub-level set filtration reveals  the oscillating information of the functions. Since each f h is a subset of [0, T] ⊆ R, it only contains 0-dimensional structures, i.e., connected components. Hence, the only nontrivial persistence diagram in this case is P 0 . For simplicity, when there is no danger of confusion, for a given function f , we use P 0 (f ) to denote P 0 ({f h i }), the persistence diagram associated with the sub-level sets filtration of f . As discussed in Edelsbrunner and Harer (2010), each element in P 0 is a min-max pair in the original function f (t). The concept of this filtration is closely related to the size function theory (see Biasotti et al., 2008 and references therein) and is commonly used as a shape descriptor (Biasotti et al., 2008). In practice, persistence diagram is robust to noise under the bottleneck distance. This fact renders persistence diagram an useful data analysis quantity. A precise statement of this robustness is below.
Theorem 3.1.1. Let X be an n-dimensional rectangle in R n . Take two continuous functions f , g : X → R with finitely many local extremums (minimums or maximums). Then, we have for q ∈ N, where γ ranges over all bijections from P q (f ) to P q (g) considering the infinite points on the diagonal.

Second Useful Filtration-Vietoris-Rips Complexes Filtration
To introduce Vietoris-Rips (VR) complexes filtration for a given time series, we first embed the time series into a high dimension point cloud via Taken's lag map (Takens, 1981), which is constructed in the following way. Take p ∈ N to be the dimension of the embedding, and τ ∈ N to be the lag step. For a given time series x : Z → R, the lag map with lag τ and dimension p is defined as which is a subset of R p . We postpone details of Taken's lag map to Supplementary Material (section 1.2). With the point cloud R p,τ (x) ⊆ R p , we are ready to introduce the Vietoris-Rips complex.
In general, given a point cloud X = {x 1 , . . . , x N } ⊂ R p , the main idea of Vietoris-Rips complex is to build simplicial complexes from X if points in X are closed enough. A formal definition is given below.

Definition 5 (Edelsbrunner and Harer, 2010 section 3.2, p. 61).
Clearly, for an increasing sequence ǫ 1 < ǫ 2 < · · · < ǫ N , the corresponding sequence of Vietoris-Rips complexes forms a filtration: After determining the representation rules of connected components, the lifespan of holes of different dimensions can be computed easily. See Supplementary Material (section 1.3) for an illustrative example of the Vietoris-Rips filtration. For simplicity, we denote the q-th persistence diagram associated with the Vietoris-Rips filtration as P q (R p,τ (x)) := P q ({VR(R p,τ (x); ǫ)} ǫ ). In parallel with Theorem 3.1.1, the stability of persistence diagrams extracted from a Vietoris-Rips filtration has been discussed in Chazal et al. (2014). Theorem 3.2.1 (Chazal et al., 2014, Theorem 5.2). For finite metric spaces (X, d X ) and (Y, d Y ), then for q ∈ N, where d B is the bottleneck distance and d GH is the Gromov-Hausdorff distance.
The formal definition of Gromov-Hausdorff distance can be found in Chazal et al. (2014) and Burago et al. (2001), and conceptually, it measures the similarity between two metric spaces under distance-preserving transformations.

Persistence Statistics
We now introduce a set of new features to summarize persistence diagrams. It is computationally efficient and straightforward to implement. We propose to explore distributions of the birth b and the death d of all possible holes, and calculate their statistic measurements. This idea is considered one of the most straightforward way to extract features from persistence diagrams (Pun et al., 2018). Despite its simplicity, it has been used in several studies, such as skin lesions classification , bifurcations analysis in dynamical systems (Mittal and Gupta, 2017), and protein classification (Cang et al., 2015).
To be more specific, given a persistence diagram P, we transform it into two multi-sets of numbers, M and L, defined as (11) Note that for the Vietoris-Rips complex filtration, d+b 2 captures the "size" of the associated hole, and d−b captures the robustness of the associated hole. On the other hand, for the sublevel set filtration, d+b 2 reveals the locations of holes, and d − b captures the differences between low and high peaks in a time series. Note that since the hole (0, ∞) always exists in the persistence diagram as is shown in the previous section, it is omitted.
In this paper, for each persistence diagram, we consider eight summary statistics to represent the multi-set M, including mean, standard deviation, skewness, kurtosis, 25th, 50th, 75th percentile, and the persistent entropy (Chintakunta et al., 2015). We number them from 1 to 8. We consider the same summary statistics for the multi-set L, and number them from 9 to 16.
Definition 6 (Persistence Statistics). Given a persistence diagram, the persistence statistics (PS) is defined as a map, (PS) , that transforms the persistence diagram to a point in R 16 .
As shown in Algorithm 1, for the Vietoris-Rips complex filtration, we consider 0-th and 1-th persistence diagrams; for the sub-level set filtration, we consider 0-th persistence diagram.
The persistent entropy of M and L, denoted as E(M) (Chung and Lawson, 2019) and E(L) (Atienza et al., 2020) respectively, describes the complexity of M and L. They are formally defined by E(L) has been used to study the cell arrangements (Atienza et al., 2019a), emotion recognition (Gonzalez-Diaz et al., 2019), and epileptic seizures detection in EEG signals (Piangerelli et al., 2018). From the theoretical perspective, E(L) is a stable measurement (Theorem 3.12 in Atienza et al., 2020). E(M) was first appeared in Chung et al. (2018) and the discussion about the stability of E(M) can be found in Chung and Lawson (2019).
Note that while intuitively, holes with long lifespans are considered important features and those with short lifespans are considered noises, in our proposed features, we do not discriminate holes with long or short lifespans. In other words, we take all holes into consideration. This approach is supported by a recent discovery that those considered as noisy holes might actually contain important information. For example, in the drivers' behavior classification (Bendich et al., 2016), authors transformed the space of persistence diagrams into "binned" diagrams, and found that the main differences occurred in those short lifespan holes. Another work on the leave classification (Patrangenaru et al., 2018) also suggested that holes with short lifespans could better distinguish different types of leaves.

APPLICATION TO SLEEP STAGE CLASSIFICATION
In recent decades, a growing body of evidence shows that sleep is not only intimately related to personal health (Karni et al., 1994;Kang et al., 2009) but also has a direct impact on public health (Colten and Altevogt, 2006). In clinics, sleep experts score sleep stage by reading the electroencephalogram (EEG), electrooculogram (EOG), and electromyogram (EMG) based on the American Academy of Sleep Medicine (AASM) criteria (Iber et al., 2007;Berry et al., 2012). Sleep, however, impacts the whole body, and we can read sleep via reading physiological signals other than EEG, particularly ECG and HRV mentioned in section 1. The relationship between HRV and sleep dynamics has been widely studied in the physiology society (Zemaityte and Varoneckas, 1984;Vaughn et al., 1995;Toscani et al., 1996;Bonnet and Arand, 1997;Elsenbruch et al., 1999;Chouchou and Desseilles, 2014;Penzel et al., 2016). Specifically, when a subject is awake, since the sympathetic tone of the ANS is dominant, he/she has a higher heart rate and a less stable heart rhythm due to external stimuli (Somers et al., 1993). When a subject is asleep, the heart rate is lower, and it reaches its lowest value during deep (slow wave) sleep (Snyder et al., 1964). During NREM (non-rapid eye movement) sleep, the parasympathetic nervous system dominates the sympathetic tone and the energy restoration and metabolic rates reach their lowest levels, so the heart rate decreases and the rhythm of the heart stabilizes (Somers et al., 1993).
The above physiological facts indicate that the heart rate rhythm provides a non-invasive window for researchers to study sleep. There have been several studies trying to classify sleep stages based solely on HRV. Most of them focus on classifying wake and sleep (Lewicke et al., 2008;Long et al., 2012;Aktaruzzaman et al., 2015;Ye et al., 2016;Malik et al., 2018), some focus on detecting drowsiness (Vicente et al., 2016), and some focus on classifying rapid eye movement (REM) and NREM (Mendez and Matteucci, 2010), or wake, REM, and NREM (Xiao et al., 2013). The challenge and difficulty of this mission can be appreciated from the reported results. In this section, we apply the TDA tool and the proposed persistence statistics to study this problem.

Datasets
The databases we use here are the same as those used in Malik et al. (2018). Here we summarize them and refer readers with interest in the database details to Malik et al. (2018). The CGMH-training database consists of standard polysomnogram (PSG) signals on patients suspicious of sleep apnea syndrome at the sleep center in Chang Gung Memorial Hospital (CGMH), Linkou, Taoyuan, Taiwan. The Institutional Review Board of CGMH approved the study protocol (No. 101-4968A3). All recordings were acquired on the Alice 5 data acquisition system (Philips Respironics, Murrysville, PA). Each recording lasts for at least 5 h. The sleep stages, including wake, REM, and NREM (REM and NREM constitute the sleep stage), were annotated by two experienced sleep specialists according to the AASM 2007 guidelines (Iber et al., 2007), and a consensus was reached. According to the protocol, the sleep specialists provide annotation for non-overlapping 30 s long epochs. In this study, we focus on the second lead of the ECG recording, which was sampled at 200 Hz. There are 90 participants without sleep apnea [each with apnea-hypopnea index (AHI) <5] in this database, among which we consider only 56 participants who have at least 10% epochs labeled as wake to avoid the imbalanced data issue.
We consider three validation databases. The first one is the CGMH-validation database. This database consists of 27 participants acquired independently of CGMH-training from the same sleep laboratory in CGMH under the same Institutional Review Board. The other two validation databases are publicly available. The DREAMS Subjects Database 1 (DREAMS), consists of 20 recordings from healthy participants, where the ECG recordings were acquired by the Brainnet TM system (Medatec, Brussels, Belgium). The sampling rate is 200 Hz, and the minimum recording duration is 7 h. Although the race information is not provided, we may assume that its population constitution is different from that of the CGMH databases since it is collected from Belgium. This database is chosen to assess the model's performance on participants of a different race recorded from different recording machine. The third database is the St. Vincent's University Hospital/University College Dublin Sleep Apnea Database (UCDSADB) from Physionet (Goldberger et al., 2000) 2 . It consists of 25 participants with sleep apnea of various severities. The ECG signal was recorded by Holter monitor at the sampling rate of 128 Hz. The minimum recording is 6 h long. We focus on the first ECG lead in this study. The UCDSADB is chosen to assess the model's performance on recordings which come from participants with sleep disorders. We remark the these validation databases are not used to tune the model's parameters, and no subject is rejected.

Time Series to Analyze-Instantaneous Heart Rate
The data preprocessing steps are the same as those shown in Malik et al. (2018). Here we summarize those steps and refer readers to Malik et al. (2018) for more details. First, apply a standard automatic R peak detection algorithm (Elgendi, 2013). Suppose there are n k R peaks in the k-th subject's ECG recording. Denote {r k,i } n k i=1 the location in time (sec) of the detected R peaks of the k-th subject. We apply the 5-beat median filter to remove artifacts in the detected R peaks; that is, if a detected beat is too close or too far from their preceding beats, it is removed or interpolated. Then, the IHR of the k-th recording, denoted as x k , is determined by the shape-preserving piecewise cubic interpolation (Task Force of the European Society of Cardiology and others, 1996) over the nonuniform sampling x k (r k,i ) = 60(r k,i − r k,i−1 ) −1 . (12) x k describes the IHR at each time in beats-per-minute. The IHR is sampled at a sampling rate of 4 Hz. We break the IHR signal into 30-s epochs following the same epoch segmentation in the experts' annotations. We discard all epochs with fewer than 5 detected R peaks. This step is adjusted by physiological knowledge. For each labeled epoch, we build a time series of 90 s in length by concatenating the epoch with the preceding 2 epochs. For the sake of handling the inter-individual variance, each 90 s time series is normalized by subtracting its median value. Thus, for the j-th epoch of the k-th recording, the associated time series we consider is

IHR Time Series and Their Persistence Diagrams
Following the discussion in section 3, we apply TDA to IHR time series defined in (13), x (k,j) . More precisely, we consider P 0 (x (k,j) ) via the sub-level set filtration, and P i (R 120,1 (x (k,j) )) for i = 0, 1, via the Vietoris-Rips complex filtration. We extract persistence statistics from both P 0 (x (k,j) ) and P i (R 120,1 (x (k,j) )), where i = 0, 1. We summarize section 3 and highlight our approach in the following pseudocode. See also Figure 2 for an illustration.
Algorithm 1: Feature Extraction Scheme Input: A time series, x(t). Output: Topological features used in this article.
We illustrate the IHR time series and their persistence diagrams with different filtrations in Figures 6, 7. From a IHR time series during a wake (resp. sleep) epoch shown in Figure 6A (resp. Figure 7), we observe that these IHR's seem to be different: wake epoch seems to have more variability than sleep one does. Sub-level set filtration captures such variability in the form of the persistence diagram. As shown in Figures 6B, 7B, their persistence diagram's of sub-level set filtration are different. Points in Figure 6B spread widely while most points in Figure 7B are clustered around lower left portion of the diagram. Moreover, Figure 6B seems to have more long-lived points than Figure 7B does. Next, we examine the persistence diagrams of Vietoris-Rips complex filtration. In this work, we take (p, τ ) = (120, 1), where p = 120 is equivalent to a 30 s long time series (since the sampling rate is 4 Hz). This set of parameters is motivated/guided by the AASM criteria where the sleep stage is assigned based on the 30-s readings. The parameters (120, 1) can be thought as sliding a window of a 30-s long time series, and P(R 120,1 (x)) stores information about changes of this 30-s sliding window over time. Figures 6C, 7C show examples of R 120,1 (x (k,j) ) projected onto their first three principal components. Visually, the point clouds of Figures 6C, 7C have different shapes (former seems to have a "lamp" shape while the latter does not), and their persistence diagrams shown in Figures 6D, 7D are also different. For instance, the red points in Figure 6D cluster around birth values 10 ∼ 20, the red points in Figure 7D have three clusters around the birth values 15 ∼ 25, 30 ∼ 40, and 55. It is important to note that the computations on P 0 (R 120,1 (x (k,j) )) are done on the R 120 space, and projection onto their first three principal components is merely for the visualization purpose.
As discussed in section 2.4, while it is possible to analyze the data via persistence diagrams, it is usually computationally challenging. The proposed persistence statistics allows us to further summarize the persistence diagrams and quantify the above observations. To examine the persistence statistics features, FIGURE 6 | An illustration of the IHR during the wake stage. (A) The IHR signal x (k,j) ; (B) the persistence diagram of the sub-level set filtration, P 0 (x (k,j) ); (C) The first three principal components of the point cloud R 120,1 (x (k,j) ); (D) the persistence diagrams of the Vietoris-Rips filtration, P 0 (R 120,1 (x (k,j) )) and P 1 (R 120,1 (x (k,j) )) are superimposed, where blue and red points represent q = 0 and q = 1, respectively. take k { (PS) (P 0 (x (k,j) ))} n k j=1 as an example. In order to compare them on the same scale, we perform the standard z-score normalization for each subject. We abuse the notation and use (PS) (P 0 (x (k,j) )) to denote the normalized parameters. In Figure 8A, we show the boxplot of each normalized persistence statistics parameter, where blue (red) bars represent the persistence statistics associated with an IHR time series associated with the sleep (wake) stage. We performed a rank sum test with the null hypothesis that two samples have equal medians with a significance level of 0.05 with the Bonferroni correction. We found that there are significant differences between waking and sleeping features for all persistence statistics parameters, except for the kurtosis of M (labeled as 4 in Figure 8A), and the median of L (labeled as 14 in Figure 8A). The boxplot as shown in Figure 8A shows that the mean and standard deviation of M are the most distinguishable persistence statistics parameters between sleep and wake epochs. To further visualize these features, we apply the principle component analysis (PCA) to k { (PS) (P 0 (x (k,j) ))} n k j=1 , and plot the first three principal components in R 3 as shown Figure 8B. We observe a separation between sleep and wake features. The visualization of (PS) (P i (R 120,1 (x (k,j) ))), where i = 0, 1, is shown in Supplementary Figure 5.

Support Vector Machine as the Learning Model
We consider the widely applied classifier with a solid theoretical foundation, the support vector machine (SVM), to establish the heartbeat classification model. This is Step 4 (machine learning step) shown in Figure 2. The linear function kernel is considered in this work and we use the Matlab built-in function fitcsvm with default parameters. The input data are features shown in (14), which are calculated by the publicly available libraries DIPHA (https://github.com/DIPHA/dipha) and Ripser (https://github.com/Ripser/ripser). When there are more than 2 classes, we apply the Error-Correcting Output Codes (ECOC) Dietterich and Bakiri (1994) with one-versus-one design. Specifically, we use the Matlab built-in function fitcecoc with default parameters. The Matlab version is 2019b.

Statistics
We carry out the cross-database validation. Specifically, we train the SVM model on one database, and evaluate the performance The persistence diagram of the sub-level set filtration, P 0 (x (k,j) ); (C) The first three principal components of the point cloud R 120,1 (x (k,j) ); (D) the persistence diagrams of the Vietoris-Rips filtration P 0 (R 120,1 (x (k,j) )) and P 1 (R 120,1 (x (k,j) )) are superimposed, where blue and red points represent q = 0 and q = 1, respectively.
FIGURE 8 | Distribution of normalized persistence statistics features, (PS) (P 0 (x (k,j) )). (A) Boxplot of the (PS) (P 0 (x (k,j) )). The numbers listed on the horizontal axis indicates the number of persistence statistics. *Indicates that the feature fails to reject the null hypothesis (that two samples have equal medians) of the significance level of 0.05 on the rank sum test with Bonferroni correction. (B) Visualization of k { (PS) (P 0 (x (k,j) ))} nk j=1 by the first three principal components.
on the other databases. One of the main challenges in this automatic annotation problem is that the datasets are usually imbalanced; for example, the number of wake epochs is usually much smaller than that of sleep epochs (e.g., in the CGMHtraining, the total number of wake epochs is 9, 150, while the total number of sleep epochs is 54, 547). Learning on imbalanced datasets is one of challenging topics in machine learning (He and Ma, 2013;Kuhn and Johnson, 2013;Fernández et al., 2018). Typically, the accuracy would heavily skew toward the majority case. Taking the above CGMH-training as an example, if a model predicts all epochs as sleep, then its accuracy is 54547/(9150 + 54547) ≈ 85%, which may seem high at the first glance. However, this model is clearly useless because it has no predictability for the wake, which can be seen through the sensitivity, 0/9150 = 0. Therefore, in the case of imbalanced datasets, both sensitivity and specificity would be important indicators to evaluate the performance of a model. They both should be as high as possible, and they should be on the similar level. In order to account for the imbalanced dataset, we adopt a down-sampling process. Let E s and E w be the collection of all sleep and wake epochs, respectively, across all subjects in the training set, and denote their cardinality by |E s | and |E w |, respectively. We take all epochs in E w , and randomly select |E w | epochs from E s . The SVM model will then be built on these balanced epochs. Once the model is built on the training dataset, we test it on the entire testing dataset.
We report the following performance measurement indices. When there are m labels, denote M ∈ R m×m to be the confusion matrix of the automatic classification model, where M kl represents the count of epochs whose known group labels are k and whose predicted group labels are l. The sensitivity (SE), positive predictivity (+P) and F1 for the k-th class, the Cohen kappa, and the overall accuracy (Acc) are defined as (16) When m = 3, k = 1 means wake, k = 2 means REM, and k = 3 means NREM. When classifying wake and sleep stages, k = 1 means wake, and k = 2 means sleep; when classifying REM and NREM stages, k = 1 means REM, and k = 2 means NREM. When m = 2, SE 1 is reduced to the usual sensitivity (SE), SE 2 is reduced to the usual specificity (SP), and +P 1 is reduced to the precision (PR). For each database and each performance measurement, we report the mean ± standard deviation of all subjects. All experiments in this and next sections were done using Windows 7 operating system environment equipped with i5-4570 CPU and 32 GB RAM. Under this computational environment, given a random seed, the whole training process of an SVM model takes 5-7 min on average. For the reproducibility purpose, the Matlab code is available in the GitHub repository website 3 .

Automatic Sleep Stage Classification Result
We performed three classification tasks-sleep v.s. wake, REM v.s. NREM, and finally wake v.s. REM v.s. NREM. The random seed is fixed to 1 in all cases when we ran the subsampling scheme. The results are shown in Tables 1-3, where the SVM model was 3 https://github.com/peterbillhu/TDA_for_SleepWake_Classifications  trained on the CGMH-training dataset and tested on CGMHvalidation, DREAMS, and UCDSADB, respectively. For the interested readers, we also include extensive experimental results with different settings in Supplementary Material (section 3), such as results of training on different datasets, and different random seeds. All results are similar to those reported in the main article. Table 1 lists the result of classifying wake and sleep stages with different testing sets. For each testing database, we show the mean±standard deviation of each prediction outcome measurement of all subjects in that database. Table 1 shows the performances of training the model on CGMH-training and testing it on CGMH-validation, DREAMS, and UCDSADB. When considering the CGMH database, the (SE, SP) pair for CGMH-training and CGMH-validation are (78.3 ± 14.7%, 76.0 ± 6.1%) and (70.9 ± 16.0%, 78.9 ± 5.4%), respectively. When testing  Table 2 shows the performance for the REM and NREM classification. In this task, since the number of NREM epochs is much more than that of REM epochs, we apply the same downsampling process to NREM as discussed in section 4.4.2. Tables 1, 2 have several similarities.
Finally, Table 3 shows the performance for the wake, REM, and NREM classification. In this experiment, since the number of NREM is much more than those of wake and REM, the down-sampling scheme is applied to NREM. The Acc's in all cases are about 60%, except the UCDSADB. The SE's of wake, REM, and NREM are balanced and consistent across databases, except UCDSADB. Again, this result might be due to the fact that UCDSADB contains subjects with sleep apnea. On the other hand, note that the +P of NREM is higher than other classes, which is expected due to the dependence of +P on the database prevalence. In Supplementary Material, we provide more crossdatabase validation results.

DISCUSSION AND CONCLUSION
In this work, the TDA tools are considered to analyze time series. Specifically, we propose a set of novel persistence statistics features to quantify HRV by analyzing IHR time series by TDA tools. The proposed HRV features are applied to predict sleep stages, ranging from wake, REM, and NREM. In addition to being computationally efficient, the algorithm is theoretically sound supported by mathematical and statistical results. Note that while we focus on the HRV analysis for the sleep stage annotation, the proposed algorithm has a potential to be applied to analyze other time series and study the HRV for other clinical problems.

Theoretical Supports and Open Problems
We find that empirically, M and L are simple yet effective representations of the persistence diagram and reveal signatures about the underlying object. It would be interesting to investigate, in theory, the probability distribution of M and L for a given simplicial complex. However, to the best of our knowledge, while there have been several works in this direction, it is still a relatively open problem. Recently, there has been some theoretical work toward this direction, namely the theory of random complexes (see e.g., survey papers Kahle, 2014;Bobrowski and Kahle, 2018 and references therein). In order to understand the role of noises in the persistence diagram, there have been studies on the topology of the noise. In the theory of random fields, authors in Mischaikow et al. (2010) used sub-level sets as filtration to study the number of components (β 0 ) with various random processes; in Adler et al. (2010), authors studied the relation between random fields and the persistent homology in general. In particular, as mentioned in Adler et al. (2010): "It would be interesting to know more about the real distributions lying behind the persistence diagram, but at this point we know very little." There is also a result in random cubical complexes (Hiraoka and Tsunoda, 2018), and a few work on the limiting theorem of total sum persistence (Owada, 2018) and persistence diagrams (Hiraoka and Tsunoda, 2018). It would also be interesting to study the stability of each persistence statistics. As of now, only the sum of L, the max of L, and the entropy of L have been shown to be stable (Cohen-Steiner et al., 2010;Atienza et al., 2020). However, the rest of persistence statistics is still unknown. Another interesting direction, instead of focusing on each persistence statistics, is to study the probability distributions of M and L. For instance, let ρ M and ρ L be the empirical probability density function of the sets M and L, respectively. For S = M or S = L, could one establish ρ S (f ) − ρ S (g) D ≤ d B (P q (f ), P q (g)), where · D is some suitable statistical distance? We leave those interesting theoretical problems to future work.

Comparison With Existing Automatic Sleep Stage Annotation Results
There have been several results in automatic sleep stage annotation by taking solely the HRV into account. A common conclusion is that classifying sleep-wake by quantifying HRV is a challenging job. In general, due to the heterogeneity of the data sets, various evaluation criteria and different features and models used in these publications, it is difficult to have a direct comparison. But to be fair, below we summarize some related existing literatures for a discussion. To the best of our knowledge, except (Malik et al., 2018), there is no result reporting a cross-database validation. For those running validation on a single database, we shall distinguish two common cross validation (CV) schemesleave-one-subject-out CV (LOSOCV) and non-LOSOCV. When the validation set and the training set come from different subjects, we call it the LOSOCV scheme; otherwise we call it the non-LOSOCV scheme. The LOSOCV scheme is in general challenging due to the uncontrollable interindividual variability, while the non-LOSOCV scheme tends to over-estimate. Therefore, for a fair comparison, below we only summarize papers considering only the IHR features and carrying out the LOSOCV scheme.
In Xiao et al. (2013), the database was composed of healthy participants aged 16 − 61 years. A random forest model was established to differentiate between the wake, REM, and NREM stages for those epochs labeled as "stationary." Based on the confusion matrix provided in Xiao et al. (2013), the SE, SP, Acc, and F1 for detecting the wake stage are 51.2, 90.2, 84.0, and 0.50%. The authors also provided the SE of wake, REM, and NREM, which are 53.72±20.15, 59.01±19.72, and 79.50±7.82% respectively. In Table 3, our validation on CGMH-validation is 61.1 ± 19.0, 67.1 ± 20.9, and 72.6 ± 6.4%, respectively. Observe that our SE of wake and REM are better, and SE of NREM is on the similar level. In addition to the balance of all classes due to the sub-sampling scheme in our result, note that we focus on all epochs but not "stationary epochs, " and the subjects in CGMHvalidation are not healthy but simply without sleep apnea.
In Mendez and Matteucci (2010), the database was composed of 24 participants aged 40 − 50 years with 0 AHI. The authors took the temporal information and the phase and magnitude of the "sleepy pole" as features to train a hidden Markov model to differentiate REM and NREM stages. The reported SE, SP, and Acc were 70.2, 85.1, and 79.3%, respectively. Our results outperform theirs. Our SE, SP and Acc of the REM and NREM classification in CGMH-validation shown in Table 2 are 78.1 ± 17.4, 79.6 ± 6.5, and 77.8 ± 8.3%, respectively. Observe that both Accs are similar which means portion of correct predictions are similar. Not only our SE is better, but SE and SP are also balanced.
In Lewicke et al. (2008), the database is composed of 190 infants. A variety of features and classification algorithms were considered and the wake and sleep classes were balanced for the analysis. The SE and SP of their multi-layer perceptron model without rejection was 79.0 and 77.5%, respectively. In Table 1, the SE and SP of our result on CGMH-validation is 70.9 ± 16.0 and 78.9 ± 5.4%. Our performance is comparable to theirs. However, there is a fundamental difference between their experiments and ours-the sleep dynamics of infants and adults are different.
In Aktaruzzaman et al. (2015), the database is composed of 20 participants aged 49-68 years with varying degrees of sleep apnea. Detrended fluctuation analysis and a feed-forward neural network were applied to differentiate the wake and sleep stages. Various epoch lengths were considered, and the highest performance was recorded on an epoch length of 5 min. The Acc, SE, SP, and Cohen's kappa were 71.9 ± 18.2, 43.7 ± 27.3, 89.0 ± 7.8, and 0.29 ± 0.24, respectively. We consider UCDSADB for a comparison. In Table 1, the Acc, SE, SP, and Cohen's kappa of our testing result on UCDSADB is 70.6 ± 5.4, 57.6 ± 15.5, 75.3 ± 5.5, and 0.238 ± 0.133%. Their Acc and ours are on the same level, our SE is better than theirs, while their SP is better than ours. However, our SE and SP are balanced compared with theirs. A major difference is that our standard deviations for Acc, SE, SP are much smaller. Thus, our performance is comparable to theirs.
To make a conclusion, we emphasize that all those results under comparison are not carried out in the cross-database scheme. Also, usually the SE and SP are not balanced with high SP, which leads to the high accuracy. Therefore, the results suggest that the proposed persistence statistics features and chosen learning model lead a better, or at least similar, performance compared with the state-of-the-art results. The cross-database validation further suggests the usability of the persistence statistics features and the proposed learning scheme in clinical setups. Last but not the least, due to the numerical efficiency of the proposed persistence statistics features, it is potential to apply it to analyze large scale time series.

Technical Issues
We remark that although it is possible to include P i (R 120,1 (x (k,j) )) for i ≥ 2 in (14), in practice, it is a challenging task due to its computational complexity. Its computation is known to be poorly scalable in dimension and memory-intensive. We refer readers to Otter et al. (2017) for more details and comparisons among state-of-arts TDA packages and extensive benchmark. To get an idea of the computational cost, for any epoch, the computational time by the state-of-art package Ripser for P 1 (R 120,1 (x (k,j) )), P 2 (R 120,1 (x (k,j) )), and P 3 (R 120,1 (x (k,j) )) are 0.06, 1.7, and 106 s in a standard laptop, respectively. This echos the fact that the computation of P i (R 120,1 (x (k,j) )) does not scale well in i Otter et al. (2017). We demonstrated on adding features P 2 (R 120,1 (x (k,j) )) and tested classification performance on DREAMS and UCDSADB datasets. The results are listed in Supplementary Tables 18, 19. Comparing  Supplementary Table 18 with Supplementary Table 1 and  Supplementary Table 19 with Supplementary Table 2, we find that the improvements are too marginal to justify the additional computational time.
Thus, it would be inefficient to obtain the higher dimensional persistence features. A possible approach for tackling the computational inefficiency is to reduce the Taken's embedding dimension. In Myers et al. (2019), the authors discussed how to find a proper embedding dimension p by considering the false nearest neighbors (Fraser and Swinney, 1986). It would be interesting to combining this embedding technique to our future works. Also, finding another reduction criterion for the Taken's embedding is also an important future direction.

Limitations and Future Directions
In addition to the theoretical development discussed above, there are several interesting practical problems left untouched. While we systematically consider the inter-individual variance, the race, the machine, and the sleep disorder by taking three different databases into account, we acknowledge the fact that the data is collected from the sleep lab. When the data is collected from the real-world mobile device, it is not clear if the algorithm could perform as well and run in real-time. Moreover, its performance for the home-based screening needs to be further evaluated. Yet, in the current mobile health market, the photoplethysmography (PPG) sensor has been widely applied, and its applicability for the sleep-wake classification has been reported in Malik et al. (2018). It is interesting to see how the TDA approach could be applied to analyze the HRV from the PPG for the sleep stage classification mission. From the data analysis perspective, it would be interesting to perform a more sophisticated analysis and take other features from the persistence diagram. For instance, the persistent homology transformation (Turner et al., 2014b) was recently developed and proven to be a sufficient statistic, and had been successfully applied to the shape analysis. It would be interesting to combine the persistent homology transformation and persistence statistics. IHR is a well-known non-stationary time series. Based on the encouraging results of applying the TDA, we suspect that the persistence statistics features could be applied to study other clinical problems related to HRV, and furthermore, analyze other physiological time series, including the multivariate ones. There has been some work using TDA tools to analyzing the multivariate time series (Merelli et al., 2015;Gidea and Katz, 2018;Wu and Hargreaves, 2021) where the critical step is to transform multivariate time series into a point cloud so that Vietoris-Rips complex persistent homology can be computed. It would also be interesting to investigate ways to use the sublevel set filtration in this context. We will explore those limitations/directions in our future work.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.