Patent Publication Number: US-11653871-B2

Title: Atrial fibrillation prediction using heart rate variability

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
     This application is a National Phase of PCT Patent Application No. PCT/IL2019/050256 having International filing date of Mar. 7, 2019, which claims the benefit of priority of U.S. Provisional Patent Application No. 62/639,543 filed Mar. 7, 2018. The content of the above applications is incorporated herein by reference in its entirety. 
    
    
     BACKGROUND 
     The heart rate never reaches a steady state, as it is controlled by complex, dynamic, chaotic processes and oscillates at different periods over different, continuously shifting time scales (Goldberger A L et al. Sci Am. 1990; 262:42-9). Therefore, it is not surprising that ECG recordings in mammals, even under resting conditions, reveal complex beat-to-beat variations in heartbeat intervals (Papaioannou V et al. Frontiers in physiology. 2013; 4:174). While a decrease in this complexity in humans with cardiovascular disease correlates with increased morbidity and mortality (Hillebrand S et al. Europace. 2013; 15:742-9), an increase in heart rate variability (HRV) above a certain threshold, or in mathematical terms, in system entropy, leads to the abnormal electrical impulse conductance that is defined as fibrillation. 
     Cardiac fibrillation is one of the leading causes of morbidity and mortality in the developed world, where atrial fibrillation (AF) is the most common sustained arrhythmia (Wolf P A et al. Stroke. 1991; 22:983-8), affecting approximately 2.2 million patients in the United States, and its prevalence increases with age (Feinberg W M et al. Archives of internal medicine. 1995; 155:469-73). Although AF is not considered a life-threating arrhythmia, it affects the circulatory system of patients and their general health and quality of life; further, it poses a significant burden on the healthcare system (Feinberg W M et al. Archives of internal medicine. 1995; 155:469-73): e.g., AF is a significant risk factor for stroke, with about 15% of strokes occurring in people with AF (Go et al. JAMA. 2001; 285:2370-5). Unfortunately, currently available drug treatments for AF are less than satisfactory. An expert panel convened by the Heart, Lung and Blood Institute to resolve this issue determined that “AF detection by using emerging methods and technologies” constitutes an important future direction in preventing AF (Benjamin E J et al. Circulation. 2009; 119:606-18). To date, there is no effective clinical method for early prediction and precise identification of AF events. 
     Interestingly, while persistent AF is associated with an increased HRV on average (Bergfeldt L et al. J Appl Physiol. 2003; 94:2217-24), reduced HRV indices are observed just prior to arrhythmogenic events (Akyurek O et al. Clinical cardiology. 2003; 26:196-200). Although the correlation between changes in heart beat complexity and prevalence of AF has been known for more than three decades (Akyurek O et al. Clinical cardiology. 2003; 26:196-200; Vaziri S M et al. Circulation. 1994; 89:724-30; Seaborn G E et al. Annals of noninvasive electrocardiology: the official journal of the International Society for Holter and Noninvasive Electrocardiology, Inc. 2014; 19:23-33), a method to detect HRV changes on line and in short period before AF event occurs does not exist. 
     SUMMARY 
     The following embodiments and aspects thereof are described and illustrated in conjunction with systems, tools and methods which are meant to be exemplary and illustrative, not limiting in scope. 
     There is provided, in an embodiment, a method comprising receiving a plurality of data segments each representing human cardiac activity; calculating, from each of said data segments, a first value representing a distribution of a heart rate variability (HRV) parameter; acquiring a signal representing cardiac activity in a patient; deriving a second value representing distribution of an HRV parameter in said signal; and indicating an atrial fibrillation (AF) event in said patient, when a difference between said first and second values is above a predetermined threshold. 
     There is also provided in an embodiment, a system comprising at least one hardware processor; and a non-transitory computer-readable storage medium having stored thereon program instructions, the program instructions executable by the at least one hardware processor to: receive a plurality of data segments each representing human cardiac activity, calculate, from each of said data segments, a first value representing a distribution of a heart rate variability (HRV) parameter, acquire a signal representing cardiac activity in a patient, derive a second value representing distribution of an HRV parameter in said signal, and indicate an atrial fibrillation (AF) event in said patient, when a difference between said first and second values is above a predetermined threshold. 
     There is further provided, in an embodiment, a computer program product comprising a non-transitory computer-readable storage medium having program instructions embodied therewith, the program instructions executable by at least one hardware processor to: receive a plurality of data segments each representing human cardiac activity; calculate, from each of said data segments, a first value representing a distribution of a heart rate variability (HRV) parameter; acquire a signal representing cardiac activity in a patient; derive a second value representing distribution of an HRV parameter in said signal; and indicate an atrial fibrillation (AF) event in said patient, when a difference between said first and second values is above a predetermined threshold. 
     In some embodiments, said plurality of data segments represent cardiac activity in a plurality of subjects. 
     In some embodiments, at least some of said data segments are acquired from said patient. 
     In some embodiments, said HRV parameter is an entropy parameter. 
     In some embodiments, said entropy parameter is modified Shannon entropy. 
     In some embodiments, said calculating comprises calculating a distribution of said HRV parameter in all said data segments. 
     In some embodiments, each of said data segments is acquired over a period of between 1 and 10 minutes. 
     In some embodiments, said signal is acquired over a period of between 0.5 and 5 minutes. 
     In some embodiments, the method further comprises training, and the program instructions are further executable to train, at a training stage, a machine learning classifier on a training set comprising: (i) said plurality of data segment, wherein at least some of said segments comprise cardiac activity reflecting an arrythmia event, and (ii) labels indicating the occurrence of an arrythmia event in said cardiac activity. 
     In some embodiments, said data segments are labelled with said labels. 
     In some embodiments, said indicating is performed by applying, at an inference stage, said machine learning classifier to said signal, to indicate an AF event in said patient. 
     In some embodiments, said indicating comprises at least one of detecting an ongoing AF event, and predicting an upcoming AF event. 
     In some embodiments, said upcoming AF event occurs between 1 and 5 minutes after said indicating. 
     In some embodiments, said detecting is based, at least in part, on calculating statistical estimators of said distribution. 
     In some embodiments, said statistical estimators comprise at least standard deviation and deviation from normal distribution. 
     In some embodiments, said deviation from normal distribution is based on the Kolmogorov-Smirnov statistic. 
     In some embodiments, the method further comprises applying, and the program instructions are further executable to apply, a clustering algorithm to said plurality of data segments, to cluster at least some of said data segments based, at least in part, on said HRV parameter. In some embodiments, said clustering algorithm is a spectral clustering algorithm. 
     In addition to the exemplary aspects and embodiments described above, further aspects and embodiments will become apparent by reference to the figures and by study of the following detailed description. 
    
    
     
       BRIEF DESCRIPTION OF THE FIGURES 
       Exemplary embodiments are illustrated in referenced figures. Dimensions of components and features shown in the figures are generally chosen for convenience and clarity of presentation and are not necessarily shown to scale. The figures are listed below. 
         FIG.  1    is a block diagram of a system for predicting and/or detecting atrial fibrillation or other irregular regularity, according to an embodiment; 
         FIG.  2    is a flowchart of a method for predicting and/or detecting atrial fibrillation or other irregular regularity, according to an embodiment; and 
         FIGS.  3 ,  4 A -B,  5 ,  6 ,  7 A-B,  8 ,  9 A-B and  10  illustrate experimental results of a method for predicting and/or detecting atrial fibrillation or other irregular regularity, according an embodiment. 
     
    
    
     DETAILED DESCRIPTION 
     A system and method are disclosed herein for detecting and/or predicting atrial fibrillation (AF) or other irregular regularity. In some embodiments, the present invention provides for real-time prediction or detection of AF or other irregular regularity by employing an algorithm based on pattern classification of the distribution of a modified Shannon entropy of beating intervals of heart beats intervals. The present method and system produce a real-time alarm when certain changes occur in the modified Shannon entropy of a patient. 
     The present invention is based on the finding that persistent AF is associated with an increased HRV on average (L. Bergfeldt and Y. Haga, “Power spectral and Poincare plot characteristics in sinus node dysfunction,”  J. Appl. Physiol ., vol. 94, no. 6, pp. 2217-2224, 2003). However, reduced HRV indices are observed just prior to events (O. Akyurek, E. Diker, M. Guldal, and D. Oral, “Predictive value of heart rate variability for the recurrence of chronic atrial fibrillation after electrical cardioversion,”  Clin Cardiol, vol.  26, no. 4, pp. 196-200, 2003). Therefore, an atypical decrease in HRV might be predictive for an upcoming episode of AF, and thus a useful marker to be used in the detection, prevention, research and treatment of AF. 
     The present invention provides clinical and diagnostic tools for early prediction and detection of AF events by examining changes in HRV. Real-time analysis of HRV of a patient, based on short segments of heartbeat intervals, according to the present method may help to detect and prevent AF events, thereby reducing the probability of associated events, such as strokes. In some embodiments, the disclosed invention allows for prediction of AF events shortly (e.g., 5 min.) before they occur. Early prediction of AF events thus will allow to determine which hemodynamic, mechanical or electrical parameters change just before each episode, and help during ablation procedure. 
     The present invention can be implemented, for example, as an early detection and prediction system for atrial fibrillation based on continuous heart rate monitoring, as well as for early prediction of arrhythmia episodes in general. In some embodiments, heart beat segments can be acquired from a mobile device, a wearable device, a device that is used to monitor or diagnose medical conditions by medical experts, and/or based on a smart implantable cardiac defibrillator. Because the disclosed method provides for quantifying the HRV accurately based on relatively short segments of heart beats, it can be implemented based on, e.g., a remote ECG acquiring system. Thus, the present method can provide for efficient and cost-effective local and/or remote monitoring of multiple AF patients. 
     In some embodiments, the present disclosure can be implemented, for example, as a diagnostic system used for assessing the underlying etiology and/or the arrythmia profile of each patient&#39;s AF events, allowing personalized treatment. 
     In some embodiments, methods for predicting and/or detecting AF according to the present disclosure comprise a first step of acquiring and/or receiving one or more prerecorded baseline segments of cardiac activity comprising multiple heart beat intervals of a subject. The baseline segments comprise a plurality of heart beat intervals taken, e.g., at rest, when no arrythmia event was in progress. In some embodiments, the recording of the baseline signal is performed over a period of, e.g., 1 to 10 minutes. A modified Shannon entropy distribution and/or other HRV indices are may then be calculated for the baseline signal. In some embodiments, the present method further comprises calculating a difference in Shannon entropies between every two or more consecutive heart beat intervals of the baseline signal, wherein the calculation of the baseline distribution is based on the calculated differences in Shannon entropies. 
     In some embodiments, the modified Shannon Entropy is used to separate patients into groups by, e.g., Spectral Clustering. 
     In real time, a second, shorter, segment of the patient cardiac activity measurement is acquired (the “real-time segment”), wherein a modified Shannon entropy distribution and/or other HRV indices are calculated for the real-time segment as well. The distributions of the baseline and real-time segments are compared to one another. In some embodiments, when a difference in distribution width between the baseline and second segments is above a predetermined threshold, an AF event is predicted and/or detected. In some embodiments, the recording of the real-time segment is performed over a period of 30 seconds to 5 minutes. 
     In some embodiments, the present disclosure may be implemented as a machine leanring classifier trained on a training set comprising a plurality of cardiac activity segments acquired from a population of subjects, wherein at an inference stage the classifier may be configured to predict and/or detect AF event in a real-time segment acquired from a monitored subject. 
     Accordingly, in some embodiments, a machine leanring classifier may be trained, e.g., on a training set comprising a plurality of recorded baseline segments measurements from, e.g., a group of subjects. In some embodiments, the baseline segments may comprise a plurality of heart beat intervals of each subject taken, e.g., at rest, when no arrythmia event was in progress. In some embodiments, each segment may comprise, e.g., 1 to 10 minutes of recorded cardiac activity. In some embodiments, an existing database of cardiac activity recordings form a plurality of subjects may comprise at least a portion of the training set. In some embodiments, a modified Shannon entropy distribution and/or other HRV indices are may be derived for each segment within the training set. In some embodiments, modified Shannon entropy may be derived, e.g., by calculating a difference in Shannon entropies between every two or more consecutive heart beat intervals of the baseline signal. 
     In some embodiments, the training set may be labelled with appropriate labels indicating the existence of an AF event within each segment. In some embodiments, the machine learning classifier may be trained on the training set to classify a segment of cardiac activity as indicative of a concurrent and/or imminent (e.g., within a time horizon of 1-5 mins.) AF event. In some embodiments, a trained classifier of the present disclosure may be used to detect and or predict an AF event, based, at least in part, on a segment of cardiac activity measurement acquired in real time from a subject, wherein the segment may comprise a recording of cardiac activity over a period of 30 seconds to 5 minutes. 
     With reference to  FIG.  1   , in some embodiments, an exemplary system  100  in accordance with the present disclosure comprises processing unit  102 , comprising at least one hardware processor. Processing unit  102  is communicatively connected to a non-transient computer readable storage device  104 . Processing unit  102  is configured to automatically control the operation of system  100 , based on one or more applications, sets of software instructions, and algorithms stored storage device  104 . Storage device  104  may be configured to have stored thereon digital representations of cardiac signals and/or measurements. In some embodiments, processing unit  102  is configured to receive and process one or more cardiac signals and/or measurements collected by one or more cardiac signal sensors  106 . System  100  may further comprise a communication module  108 , which may be configured to connect system  100  through a communication network, such as the Internet, a local area network, a wide area network, and/or a wireless network. For example, system  100  may be configured to receive cardiac signals from a remote acquisition unit over a network via communication module  108 . System  100  described herein is only an exemplary embodiment, and may have more or fewer components than shown, may combine two or more components, or a may have a different configuration or arrangement of the components. The various components of system  100  may be implemented in hardware, software, or a combination of both hardware and software. According to various other embodiments, processing unit  102  or processing tasks performed thereby may be implemented by a handheld or worn computing device such as, but not limited to, a smart phone, a tablet computer, a notepad computer, smart watch, and the like. In addition, aspects of the present system which can be implemented by computer program instructions, may be executed on a general-purpose computer, a special-purpose computer, or other programmable data processing apparatus. 
       FIG.  2    is a flowchart illustrating the functional steps of an exemplary prediction and/or detection machine leanring classifier  200 , according to some embodiments of the present disclosure. Classifier  200  is configured to calculate HRV over relatively short real-time segments of heart beat, based on pattern classification of distribution of the modified Shannon entropy of beating intervals. Classifier  200  works in two stages: a learning stage (steps  202 - 206 ) and an inference stage (steps  208 - 214 ). 
     At a learning stage, starting at  202 , a training set of baseline cardiac activity segments of a single or a plurality of subjects may be acquired or received. The baseline segments can be acquired using an ECG and/or similar acquiring methodology from the subjects at rest. In some embodiments, the baseline segments comprise a segment of ECG and/or other bio-signal measurement which can be used to calculate the heart beating intervals (e.g., saturated oxygen, pulse detector, blood pressure, etc.). In some embodiments, each segment is recorded over a time period that is sufficient to create a reliable distribution curve for statistics, for example, 5 minutes. In some embodiments, the recording of the baseline segments is performed over a period of 1 to 10 minutes where no arrhythmia events are present. 
     At  204 , modified Shannon entropy is calculated with respect to each segment in the training set. For example, the difference in Shannon entropy of a pair in each three consecutive beats (or more beats with N−1 combination) is calculated (to reduce noise in the signal). The absolute value between each pair of beats is calculated, to provide a modified Shannon entropy value. 
     At  206 , the distribution of modified Shannon entropy values is calculated for the entire segment. In some embodiments, the distribution is normalized for the average modified Shannon entropy, and the shape of this probability density function is calculated. In some embodiments, the distribution of the modified Shannon entropy of the segment is represented by estimation of statistical parameters. 
     At  208 , a machine learning classifier is trained on the training set acquired and processed in steps  202 - 206 . 
     At an inference stage  210 , a real-time sample segment is taken over a relatively short period of time, e.g., 1 min. In some embodiments, the real-time recording of is performed over a period of 30 seconds to 5 minutes. 
     At  212 , classifier  200  is applied to the real-time segment, to predict and/or detect an AF event, based, at least in part, on a distribution parameter of Shannon entropy and/or other HRV indices in the real-time segment. 
     At  214 , in case of significant differences from the baseline distribution, the real-time measurements and processing may be repeated, e.g., between 2 and 6 times. 
     In some embodiments, the present disclosure employs, at least in part, the following definitions and assumptions for predicting and/or detecting AF events. For example, in the case of a finite set of consecutive beat intervals from a certain subject (time intervals between consecutive heart beats):
 
D={BI i } i=1   N , BI i ∈ ,i∈ 
 
The present modified Shannon entropy of order 2 may be defined as
 
 Ŝ   i   2 ( D )=| BI   i−2   −BI   i−1   |−|BI   i−1   −BI   i |,
 
     And for a general order M where M∈ , M&gt;2:
 
Ŝ i   M =|Ŝ i−1   M−1 |−|Ŝ i   M=1 |
 
     It is further assumed that there exists a probability density function parametric family
 
 Y={ƒ   x ( x ,θ)|θ∈Θ},
 
     It should be noted that θ might be a multi-dimensional parameter set, for which distribution there will be used the notion Y(θ). Satisfying that for every possible set D of consecutive beat intervals from the same subject exists a θ∈Θ, such that
 
 Ŝ   i   M   ˜Y (θ)∀ BI∈D.  
 
     For the probability density function ƒ x (x, θ), there is defined a “distribution width” function:
 
 W:Y→     + 
 
where W(f x , (x, θ)) might be thought of as a generalized variance.
 
     At the learning stage, the prediction algorithm works to estimate the characteristic θ for a certain subject in an arrhythmia-free time interval. For example, in the case of a learning set
 
 D   l ={Bl i } i=1   N ,Bl i ∈ ∀i ∈ 
 
satisfying
 
                   ∑     i   =   1     N       BI   i       ≈     5   ⁢         minutes       ,         
there is calculated an estimation for θ, for which the notion {circumflex over (θ)} will be used. For {circumflex over (θ)} there is calculated the estimated distribution width
 
 Ŵ=W (ƒ x ,( x ,{circumflex over (θ)}))
 
     In the real-time calculation stage, the prediction algorithm defines a “narrow” distribution as:
 
 Y (θ)| W (ƒ x ( x ,θ))&lt; T ( Ŵ ),
 
Where, T(w) is a function T: →  used in threshold determination. One might think of T(w) as a function determining how low a value of W should be to satisfy P(w≤W)&lt;α. The real-time calculation set is then:
 
 D   RT ={Bl i } i=1   N ,Bl i ∈ ∀i ∈ ,
 
satisfying
 
     
       
         
           
             
               
                 
                   ∑ 
                   
                     i 
                     = 
                     1 
                   
                 
                 N 
               
               
                 BI 
                 i 
               
             
             ≈ 
             
               1 
               ⁢ 
                   
               
                 minute 
                 . 
               
             
           
         
       
     
     Accordingly, with every new beat, the prediction algorithm estimates the current Ŵ for n iterations. In some embodiments, n is a predetermined constant number, which can be chosen based on either heuristic methods, or through statistical analysis. If n or more consecutive real-time calculation sets are considered “narrow,” it is considered as an “alarm.” The implementation of the present prediction algorithm requires determination of an estimator for the distribution width W and a threshold estimating function T. 
     In some embodiments, the present disclosure can be implemented, for example, as a diagnostic system used for assessing the underlying etiology and/or the arrythmia profile of each patient&#39;s AF events, allowing personalized treatment. In some embodiments, the modified Shannon Entropy is used to separate patients into groups by Spectral Clustering (see, e.g., Andrew Y. Ng, Michael I. Jordan, Yair Weiss, “On spectral clustering: analysis and an algorithm”, Proceedings of the 14th International Conference on Neural Information Processing Systems: Natural and Synthetic, Dec. 3-8, 2001). In some embodiments, the spectral clustering is based on a Laplace Kernel distance metric: 
               LaplaceDistance   =     e       1   N     ⁢     ∑       ❘   &#34;\[LeftBracketingBar]&#34;           Patient   ′     ⁢   sEntropy     -   GroupEntropy       ❘   &#34;\[RightBracketingBar]&#34;               ,         
where N is the length of the patient&#39;s entropy signal.
 
     In some embodiments, the threshold estimating function, T, may be implemented using classic statistical estimators. In some embodiments, threshold estimating function, T, may be implemented using a machine learning algorithm. 
     In some embodiments, an exemplary detection algorithm according to some embodiments is based on the assumptions and definitions above, with the following additional definitions. For a set of N observations D={X i }N i=1   N , the empirical distribution function is defined as: 
     
       
         
           
             
               
                 F 
                 N 
               
               ( 
               x 
               ) 
             
             = 
             
               
                 1 
                 N 
               
               ⁢ 
               
                 
                   
                     ∑ 
                     
                       i 
                       = 
                       1 
                     
                   
                   N 
                 
                 
                   
                     1 
                     
                       { 
                       
                         x 
                         &gt; 
                         
                           X 
                           i 
                         
                       
                       } 
                     
                   
                   ⁢ 
                   
                     ( 
                     x 
                     ) 
                   
                 
               
             
           
         
       
     
     The estimators for the expectancy standard deviation are defined as the maximum likelihood estimators: 
                 μ   ^     =       1   N     ⁢       ∑     i   =   1     N       X   i           ,       σ   ^     =         1     N   -   1       ⁢         ∑     i   =   1       N       (       X   i   2     -     μ   ^       )                   
The estimated normal distribution is defined as:
 
                 F   ~     (   x   )     =     Φ   ⁡   (       x   -     μ   ˆ         σ   ^       )           
where Φ is the cumulative distribution function of a normal standard random variable and the deviation from normality is defined as the Kolmogorov-Smirnov statistic:
 
     
       
         
           
             
               KS 
               ⁡ 
               ( 
               D 
               ) 
             
             = 
             
               
                 sup 
                 x 
               
               ⁢ 
               
                 
                   ❘ 
                   &#34;\[LeftBracketingBar]&#34; 
                 
                 
                   
                     
                       F 
                       N 
                     
                     ( 
                     x 
                     ) 
                   
                   - 
                   
                     
                       F 
                       ~ 
                     
                     ( 
                     x 
                     ) 
                   
                 
                 
                   ❘ 
                   &#34;\[RightBracketingBar]&#34; 
                 
               
             
           
         
       
     
     Experimental Results 
     Atrial Fibrillation Prediction 
     Although under the present algorithm assumptions described above, the probability density function of the modified entropy exists, even if it is assumed that beat intervals are normally distributed, the use of a sum of absolute values in the modified entropy calculation might result in a different parametric family for the modified entropy. Moreover, the calculation of consecutive modified entropies used overlapping beat intervals, thus, one cannot assume Ŝ i  and Ŝ i+1  to be statistically independent. Thus, under different assumptions, different parametric families of distribution functions may be used. 
     This theoretical result was clearly observed in numerical attempts to fit a normal distribution density function to sets of modified entropy samples, even when calculated from signals taken from the healthy subjects.  FIG.  3    illustrates a density function fitting calculated and plotted by a Matlab Distribution fitting toolbox, in which a discrepancy between the histogram and the best fitting is evident. Accordingly, instead of finding a parametric family and truly calculating the distribution width from the estimated probability density function, in some embodiments, the algorithm is configured to try to estimate the distribution width directly from the samples, some examples for such an estimator are described below. In some embodiments the algorithm directly assesses the need for an alarm from the beat intervals sequence. Nevertheless, in some embodiments, the full process of distribution parametric family determination, parameter estimation and threshold estimation is utilized. 
     The “W” Estimator 
     For a set of N modified entropy samples D={Ŝ i } i=1   N ,  first calculate the empirical cumulative distribution function: 
                     F   ^     n     (     S   ^     )     =       1   N     ⁢         ∑     i   =   1       N       1     {         S   ^     i     &lt;     S   ^       }             ,         
where 1 A  is the indicator function of event A. Then, define a tail
 
               0   &lt;   tail   &lt;     1   2       ,         
and find the samples Ŝ low , Ŝ high  satisfying:
 
     
       
         
           
             
               
                 
                   
                     
                       
                         s 
                         ˆ 
                       
                       low 
                     
                     = 
                     
                       
                         max 
                         
                           
                             i 
                             = 
                             1 
                           
                           , 
                           2 
                           , 
                           
                             … 
                             ⁢ 
                                  
                             N 
                           
                         
                       
                       
                         { 
                         
                           
                             S 
                             ˆ 
                           
                           i 
                         
                         } 
                       
                     
                   
                   , 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                         : 
                       
                       ⁢ 
                       
                         
                           S 
                           ˆ 
                         
                         i 
                       
                     
                     &lt; 
                     tail 
                   
                 
               
             
             
               
                 
                   
                     
                       
                         S 
                         ˆ 
                       
                       
                         h 
                         ⁢ 
                         i 
                         ⁢ 
                         g 
                         ⁢ 
                         h 
                       
                     
                     = 
                     
                       
                         min 
                         
                           
                             i 
                             = 
                             1 
                           
                           , 
                           2 
                           , 
                           
                             … 
                             ⁢ 
                                  
                             N 
                           
                         
                       
                       
                         { 
                         
                           
                             S 
                             ˆ 
                           
                           i 
                         
                         } 
                       
                     
                   
                   , 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                         : 
                       
                       ⁢ 
                       
                         
                           S 
                           ˆ 
                         
                         i 
                       
                     
                     &gt; 
                     
                       1 
                       - 
                       tail 
                     
                   
                 
               
             
           
         
       
     
     The estimator to the distribution width is then:
 
 W ( D )= Ŝ   high   −Ŝ   low  
 
     The use of empirical distribution functions yielded an interesting observation.  FIG.  4    shows multiple empirical distribution functions of 1-minute sets of beat intervals from the same subject on the same figure. Difference can be noted between panel A, plotted using a recording from a subject with occasional AF episodes, and panel B, plotted using a recording from a healthy subject. In the healthy subject, it can qualitatively be seen that although some variability between the distribution of different sets exists, the variability is low and there are no distinct groups of curves. Conversely, in panel A, representing the AF patient&#39;s records, larger variability can be seen with  2  distinguishable groups of distribution width. Another, more subtle difference is the presence of some very narrow curves in the AF patient&#39;s results, which are absent in the healthy patient&#39;s results. These differences suggest that the empirical distribution width might hold the information needed to differentiate patients with and without AF. 
     The {circumflex over (σ)} 2  Estimator 
     A natural choice for a parameter to estimate distribution width is variance. The estimation of variance was made by the unbiased version of the maximum likelihood estimator: 
     
       
         
           
             
               
                 σ 
                 ^ 
               
               2 
             
             = 
             
               
                 1 
                 
                   n 
                   - 
                   1 
                 
               
               ⁢ 
               
                 
                   
                     ∑ 
                     
                       i 
                       = 
                       1 
                         
                     
                   
                   n 
                 
                 
                   
                     
                       ( 
                       
                         
                           x 
                           i 
                         
                         - 
                         
                           x 
                           _ 
                         
                       
                       ) 
                     
                     2 
                   
                   . 
                 
               
             
           
         
       
     
     To reduce the noise, a moving average filter may be used, which produce a smoother signal.  FIG.  5    shows a time plot with the {circumflex over (σ)} 2  estimator, while  FIG.  6    shows the same time interval with the filtered estimator, wherein annotations of AF onset are marked by Xs. Higher values during arrhythmia episodes can be seen, with and the low values often preceding their onset. 
     An example of the threshold determination from a learning set {BI i } i=1   N , comprised of beat intervals from a patient, may be:
         Calculate the modified entropy for each possible time point to yield a set of modified entropies {Ŝ i   2 } i=2   N−2 .  Calculate the average heart rate from the first minute of the sample, HR, and then, for each beat from the HR+1 beat and onward, calculate the distribution width estimation:       

     
       
         
           
             
               
                 σ 
                 ^ 
               
               i 
               2 
             
             = 
             
               
                 1 
                 HR 
               
               ⁢ 
               
                 
                   
                     ∑ 
                     
                       k 
                       = 
                       
                         i 
                         - 
                         HR 
                       
                     
                   
                   i 
                 
                 
                   
                     ( 
                     
                       
                         
                           S 
                           ^ 
                         
                         k 
                       
                       - 
                       
                         
                           1 
                           HR 
                         
                         ⁢ 
                         
                           
                             
                               ∑ 
                               
                                 l 
                                 = 
                                 
                                   i 
                                   - 
                                   HR 
                                 
                               
                             
                             i 
                           
                           
                             
                               S 
                               ^ 
                             
                             k 
                           
                         
                       
                     
                     ) 
                   
                   2 
                 
               
             
           
         
       
         
         
           
             to yield the set:
 
 W={{circumflex over (σ)}   i   2 } HR+1   N ,
 
             Then calculate the mean of the variance estimator: 
           
         
       
    
     
       
         
           
             
               
                 
                   
                     σ 
                     ^ 
                   
                   2 
                 
                 _ 
               
               = 
               
                 
                   1 
                   
                     N 
                     - 
                     HR 
                     - 
                     1 
                   
                 
                 ⁢ 
                 
                   
                     
                       ∑ 
                       
                         i 
                         = 
                         
                           HR 
                           - 
                           1 
                         
                       
                     
                     N 
                   
                   
                     
                       σ 
                       ^ 
                     
                     i 
                     2 
                   
                 
               
             
             , 
           
         
       
         
         
           
             where the threshold is defined as:
 
 T =α {circumflex over (σ)} 2   ,0&lt;α&lt;1.
 
           
         
       
    
     Using this threshold, the preliminary definition for an alarm is:
 
Alarms={{circumflex over (σ)} 2   ∈W|{circumflex over (σ)}   2   ≤T} 
 
       FIGS.  7 A- 7 B  shows a plot of the {circumflex over (σ)} 2  estimator after a running average filter in a representative time interval from an AF subject. Annotations of AF onset are marked by Xs, and alarms are marked by Xs of lighter shade. 
       FIG.  7 A  shows the same interval from the previous figures, now with the alarms (using 
               α   =     1   200       )         
marked on the plot.  FIG.  7 B  shows a similar plot of another time interval from the same subject.
 
     As a point of comparison,  FIG.  8    shows a similar plot from a healthy subject. Even though some fluctuating does occur in healthy subjects, extremely low values are rare, and therefore in most of healthy subjects, no alarms are produced. 
     Atrial Fibrillation Patient Clustering 
     Using the Laplace Kernel distance metric for the Spectral Clustering algorithm, we clustered a patients cohort of both AF and healthy patients into 3 groups. Than we calculated the mean Multi-Scale Entropy (see, Costa M. at. Al., Physical review letters 89, no. 6 (2002): 068102) measure for each group. In  FIGS.  9 A- 9 B , one can see the correlation between the 3 groups from a controlled rabbits&#39; experiment ( FIG.  9 A , see, e.g., Rosenberg, A. A. (2018), “Non-Invasive In-Vivo Analysis of Intrinsic Clock-Like Pacemaker Mechanisms: Decoupling Neural Input Using Heart Rate Variability Measurements” (master&#39;s thesis)) to that of the clustering ( FIG.  9 B ). On further investigation, the clustering was found to group patients with similar pattern of arrhythmia. 
     In some embodiments, the clustering method herein is used to classify patients to groups to assist the personalization of medical treatment. 
     Atrial Fibrillation Detection 
     Atrial fibrillation is classically defined as “ . . . a supraventricular arrhythmia characterized electrocardiographically by low-amplitude baseline oscillations (fibrillatory or f waves) and an irregularly irregular ventricular rhythm . . . ” (see, Mann, Douglas L. et al., “Braunwald&#39;s Heart Disease: A Textbook of Cardiovascular Medicine.” 10th edition. Philadelphia, Pa.: Elsevier/Saunders, 2015). Thus, in some embodiments, an algorithm identifying irregular irregularity of heart rate is used to detect atrial fibrillation episodes. 
     Heart rate irregularity show a strong correlation with the standard deviation of the modified Shannon entropy, but does not tell regular irregularity and irregular irregularity apart. The normality measure described above tell regular irregularity and irregular irregularity apart but is not specific for irregular irregularity and may give similar results in irregular irregularity and in normal sinus rhythm with complex dynamics. Thus, in some embodiments, both measures are used together to detect AF episodes. In  FIG.  10    one may see that although both the standard deviation (denoted “var” in the legend) and the normality measure (denoted “ks” in the legend) can rise separately on non-AF times, they rise together only during the AF episodes. 
     As will be appreciated by one skilled in the art, aspects of the present invention may be embodied as a system, method or computer program product. Accordingly, aspects of the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment (including firmware, resident software, micro-code, etc.) or an embodiment combining software and hardware aspects that may all generally be referred to herein as a “circuit,” “module” or “system.” Furthermore, aspects of the present invention may take the form of a computer program product embodied in one or more computer readable medium(s) having computer readable program code embodied thereon. 
     Any combination of one or more computer readable medium(s) may be utilized. The computer readable medium may be a computer readable signal medium or a computer readable storage medium. A computer readable storage medium may be, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the computer readable storage medium would include the following: an electrical connection having one or more wires, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the context of this document, a computer readable storage medium may be any tangible medium that can contain, or store a program for use by or in connection with an instruction execution system, apparatus, or device. 
     A computer readable signal medium may include a propagated data signal with computer readable program code embodied therein, for example, in baseband or as part of a carrier wave. Such a propagated signal may take any of a variety of forms, including, but not limited to, electromagnetic, optical, or any suitable combination thereof. A computer readable signal medium may be any computer readable medium that is not a computer readable storage medium and that can communicate, propagate, or transport a program for use by or in connection with an instruction execution system, apparatus, or device. 
     Program code embodied on a computer readable medium may be transmitted using any appropriate medium, including but not limited to wireless, wireline, optical fiber cable, RF, etc., or any suitable combination of the foregoing. 
     Computer program code for carrying out operations for aspects of the present invention may be written in any combination of one or more programming languages, including an object-oriented programming language such as Java, Smalltalk, C++ or the like and conventional procedural programming languages, such as the “C” programming language or similar programming languages. The program code may execute entirely on the user&#39;s computer, partly on the user&#39;s computer, as a stand-alone software package, partly on the user&#39;s computer and partly on a remote computer or entirely on the remote computer or server. In the latter scenario, the remote computer may be connected to the user&#39;s computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet Service Provider). 
     Aspects of the present invention are described herein with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems) and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and/or block diagrams, and combinations of blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a hardware processor of a general-purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks. 
     These computer program instructions may also be stored in a computer readable medium that can direct a computer, other programmable data processing apparatus, or other devices to function in a particular manner, such that the instructions stored in the computer readable medium produce an article of manufacture including instructions which implement the function/act specified in the flowchart and/or block diagram block or blocks. 
     The computer program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other devices to cause a series of operational steps to be performed on the computer, other programmable apparatus or other devices to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide processes for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks. 
     The flowcharts and block diagrams in the Figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods and computer program products according to various embodiments of the present invention. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of code, which comprises one or more executable instructions for implementing the specified logical function(s). It should also be noted that, in some alternative implementations, the functions noted in the block may occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowchart illustration, and combinations of blocks in the block diagrams and/or flowchart illustration, can be implemented by special purpose hardware-based systems that perform the specified functions or acts, or combinations of special purpose hardware and computer instructions. 
     The descriptions of the various embodiments of the present invention have been presented for purposes of illustration, but are not intended to be exhaustive or limited to the embodiments disclosed. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the described embodiments. The terminology used herein was chosen to best explain the principles of the embodiments, the practical application or technical improvement over technologies found in the marketplace, or to enable others of ordinary skill in the art to understand the embodiments disclosed herein. 
     In the description and claims of the application, each of the words “comprise” “include” and “have”, and forms thereof, are not necessarily limited to members in a list with which the words may be associated. In addition, where there are inconsistencies between this application and any document incorporated by reference, it is hereby intended that the present application controls.