Heart rate n-variability (HRnV) measures for prediction of mortality in sepsis patients presenting at the emergency department

Sepsis is a potentially life-threatening condition that requires prompt recognition and treatment. Recently, heart rate variability (HRV), a measure of the cardiac autonomic regulation derived from short electrocardiogram tracings, has been found to correlate with sepsis mortality. This paper presents using novel heart rate n-variability (HRnV) measures for sepsis mortality risk prediction and comparing against current mortality prediction scores. This study was a retrospective cohort study on patients presenting to the emergency department of a tertiary hospital in Singapore between September 2014 to April 2017. Patients were included if they were above 21 years old and were suspected of having sepsis by their attending physician. The primary outcome was 30-day in-hospital mortality. Stepwise multivariable logistic regression model was built to predict the outcome, and the results based on 10-fold cross-validation were presented using receiver operating curve analysis. The final predictive model comprised 21 variables, including four vital signs, two HRV parameters, and 15 HRnV parameters. The area under the curve of the model was 0.77 (95% confidence interval 0.70–0.84), outperforming several established clinical scores. The HRnV measures may have the potential to allow for a rapid, objective, and accurate means of patient risk stratification for sepsis severity and mortality. Our exploration of the use of wealthy inherent information obtained from novel HRnV measures could also create a new perspective for data scientists to develop innovative approaches for ECG analysis and risk monitoring.


Introduction
Sepsis is a potentially life-threatening condition caused by the body's dysregulated response to infection [1]. Every year, over 50 million people are affected, resulting in over five million a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 deaths worldwide [2]. Prompt recognition and treatment of sepsis has been shown to impact patient outcomes, and guidelines have been developed for its management [3]. There is, however, a need for a rapid method to grade sepsis severity and prognosticate the risk for mortality in septic patients. A quick and accurate triage tool for risk stratification of septic patients presenting at the emergency department (ED) would be invaluable, allowing for greater confidence in clinical decisions, and in guiding management.
Several common disease severity scoring systems that have been utilised in the ED for the prediction of sepsis mortality including the Mortality in ED Sepsis (MEDS) score [4], quick SOFA (qSOFA) [5], and intensive care unit (ICU)-based scores such as the Sequential Organ Failure Assessment (SOFA) score [6], and the well-established Acute Physiology and Chronic Health Evaluation II (APACHE II) score [7]. Although these scoring systems have shown good predictive value, certain limitations have prevented their widespread adoption [8][9][10][11]. In recent years, heart rate variability (HRV) measurements derived from electrocardiogram (ECG) tracings have allowed for an alternative and complementary approach to predict sepsis mortality. HRV analysis measures the beat-to-beat variation between each R-R interval on an ECG tracing and reflects the autonomic regulation of the cardiovascular system [12]. Being a non-invasive tool that can be rapidly obtained even from patients who are unable to give a history, HRV has been shown to be dysregulated in sepsis [13] and correlates well with subsequent mortality [14,15]. Indeed, scoring systems that incorporate HRV parameters among its predictors have outperformed traditional clinical indicators and established disease severity scores in predicting sepsis mortality [16][17][18][19]. The use of HRV may thus further enhance our ability to stratify for risk of sepsis mortality.
In our previous work [20], we invented novel heart rate n-variability (HRnV) parameters to provide enhanced prognostic information to complement traditional HRV parameters. The proposed HRnV has two measures-HR n V and HR n V m . HR n V is derived from non-overlapping R-R intervals, while HR n V m is computed from overlapping R-R intervals. For each of the traditional HRV, HR n V, and HR n V m measures, time domain, frequency domain, and nonlinear analysis will yield its respective set of parameters. An application of the novel HRnV variables demonstrated improved predictive ability for major adverse cardiac events among patients with chest pain presenting at the ED [20].
This paper aims to study the prognostic ability of HRnV measures alongside traditional HRV parameters in predicting the outcomes in septic patients presenting at the ED and comparing the HRnV-based model with existing mortality prediction scores.

Study design and clinical setting
We conducted a retrospective cohort analysis on a convenience sample of patients presenting to Singapore General Hospital (SGH) between September 2014 to April 2017. SGH is the largest hospital in Singapore, with its ED seeing 300 to 500 patients daily. Patients are triaged on presentation at the ED according to a symptom-based Patient Acuity Category Scale (PACS).

Study population and eligibility
Patients were included in the study if they were aged 21 years and above, triaged to either PACS 1 or 2 at the ED, suspected to have sepsis as determined by their attending physician, and if they met two or more out of four Systemic Inflammatory Response Syndrome (SIRS) criteria [21,22]. The SIRS criteria (temperature >38˚C or <36˚C, heart rate >90 beats per minute, respiratory rate >20 breaths per minute, and total white blood cell count >12,000/mm 3 or <4000/mm 3 ) were used despite recent revisions under the Sepsis-3 consensus that recommend for sepsis screening with qSOFA score [1]. This decision was made primarily to allow for comparability with the existing literature. Additionally, subsequent validation studies have disputed the utility of qSOFA over SIRS for sepsis screening in the ED due to its poor sensitivity for septic patients [23][24][25][26]. Patients were excluded if their ECGs had non-sinus rhythm, a high noise level (>30% of the entire recording), or if they had a pacemaker or were on mechanical ventilator support.

Data collection
Five-minute one-lead ECGs were performed on patients who met the inclusion criteria using the ZOLL X Series monitor/defibrillator (ZOLL Medical Corporation, Chelmsford, MA). In addition, patient demographics, vital signs taken at triage, medical history, and laboratory investigations performed in the ED were retrieved from the electronic medical records. We defined the primary outcome as 30-day in-hospital mortality (IHM).

HRnV measure and analysis
We processed the ECGs and detected QRS complex to convert the original ECG signals into R-R interval (RRI) sequences (i.e., intervals of consecutive R peaks in ECGs). Fig 1 illustrates the definitions of RRI and the derived RR n I and RR n I m sequences. Conventional HRV analysis evaluates consecutive single RRIs in ECGs. Novel HRnV measures (HR n V and HR n V m ) analyse consecutive combined RRIs (RR n I and RR n I m ) [20].
To define the HR n V measure, a new type of RRI called RR n I is obtained, where n is an integer between 1 and N and N (the number of conventional RRIs combined to form a new RR n I; for example, RR 2 I is a combination of 2 consecutive RRIs) is much smaller thanN (total number of RRIs). With newly generated RR n I sequences, traditional time and frequency domains, and nonlinear analyses [27,28] are applied to calculate HR n V parameters. In addition to conventional HRV parameters, HR n V also evaluates two newly created parameters: NN50n and pNN50n. These two parameters differ from the traditional NN50 and pNN50 parameters in that the threshold is changed from 50 ms to 50×n ms in describing the absolute difference between successive RR n I sequences.
Similarly, HR n V m is a measure derived from RR n I m , where m is the number determining non-overlapping RRIs for each RR n I. When m = n, RR n I m becomes RR n I as there are no overlapping RRIs, resulting in an upper limit of N-1 for m. Fig 1 depicts a scenario when n = 3 and m can be 1 or 2. Utilising all permissible combinations of n and m, N(N+1)/2 sets of traditional HRV, novel HR n V and HR n V m parameters can be generated from a single RRI sequence. Our analysis set the upper limit of N as three due to the relatively short duration of collected ECG samples. As a result, one set of HRV parameters, two sets of HR n V (HR 2 V and HR 3 V) parameters, and three sets of HR n V m (HR 2 V 1 , HR 3 V 1 , and HR 3 V 2 ) parameters were calculated. The HRnV-Calc software suite (https://github.com/nliulab/HRnV) was used for calculating the HRV and HRnV parameters, in which the functions from PhysioNet Cardiovascular Signal Toolbox [29] were performed for ECG signal processing.

Statistical analysis
Categorical variables were compared between patients who did and did not meet the primary outcome (30-day IHM) using χ 2 test or Fisher's exact test where appropriate. Continuous variables were checked for normality with the Kolmogorov-Smirnov test. Subsequently, normally distributed variables were presented as mean and standard deviation (SD) and were compared with independent two-tailed t test between groups, while non-normally distributed variables were presented as median and interquartile range (IQR; 25 th to 75 th percentiles) and compared using the Mann-Whitney U test.
Univariable regression analysis was conducted on traditional HRV parameters, novel HRnV parameters and demographic and clinical variables. Each variable was evaluated as an individual predictor of the primary outcome (30-day IHM) using binary logistic regression with odds ratio (OR), 95% confidence interval (CI), and p-value reported. For multivariable regression analysis, we adjusted for age, temperature, systolic blood pressure, heart rate, and Glasgow Coma Scale (GCS) as these variables were either shown to be significant predictors of sepsis mortality in previous literature [15,[30][31][32], or are included in well-established sepsis scoring systems such as the National Early Warning Score (NEWS) [33], Modified Early Warning Score (MEWS) [34], qSOFA, or APACHE II. HRV and HRnV parameters were included in the multivariable analysis if they achieved p<0.2 in the univariable analysis. Included variables were then checked for collinearity using Pearson's R correlation. For each collinear pair, the variable with the higher p-value on univariate analysis was eliminated until no collinear pairs remained.
The remaining variables were then fed into a backward stepwise multivariable logistic regression model, which used p<0.1 as an endpoint. We took statistical significance at p<0.05. Backward elimination was chosen for our stepwise variable selection because it has the advantage to assess the joint predictive ability of variables, and it removes the least essential variables. However, the eliminated variables cannot re-enter the model [35]. In comparison, all possible subset selection examines every combination of variables, requiring tremendous computing resources yet likely overfitting the model when the number of variables is large [35].
In predictive modelling with the selected variables, we conducted 10-fold cross-validation to avoid overfitting in evaluating models. We split the entire dataset into 10 non-overlapping subsets of equivalent size and then used nine subsets to build a model and validated the model with the remaining one subset. We repeated the above process ten times to ensure that each of the ten subsets could be validated. Subsequently, a receiver operating characteristic (ROC) curve was plotted to assess the predictive ability of the multivariable regression model and compared against other established disease scoring systems on their area under the curve (AUC).
Missing data were addressed by median imputation, in consideration of the low proportion of missing data (<0.3%) for each variable, the nature of variables, and recommendations for missing data in clinical trials [36]. There were three missing observations for which the median value was imputed; one patient had an unknown medical history of cancer, and another patient was missing both initial and worst qSOFA scores.
All statistical analyses were carried out using Python version 3.8.0 (Python Software Foundation, Delaware, USA) using the SciPy library (version 1.3.1). Regression models were built using the StatsModels library (version 0.10.2) and scikit-learn library (version 0.22). All methods were implemented in accordance with relevant guidelines and regulations.

Patient recruitment
Fig 2 presents the patient recruitment flowchart. Of the 659 patients that were initially recruited, 190 patients did not meet the SIRS criteria, and 127 patients had inapplicable ECG readings. Three hundred forty-two patients were included for analysis and classified depending on whether they met the primary outcome of 30-day IHM (n = 66, 19%) or did not meet the primary outcome (n = 276, 81%).

PLOS ONE
Heart rate n-variability (HRnV) measures for prediction of mortality in sepsis patients Table 1 illustrates baseline characteristics and clinical parameters of patients who met and did not meet with 30-day IHM. Patients who met with 30-day IHM were older and presented with higher respiratory rates but lower temperatures, systolic blood pressures (SBP) and GCS scores, compared to patients who did not meet with 30-day IHM. The worst recorded values of respiratory rate, GCS, and SBP during each patient's ED stay were also significantly more abnormal in patients that met with 30-day IHM. The difference in disposition from the ED was significant, with a larger proportion of patients who eventually met with 30-day IHM requiring admission to the ICU as compared to patients who did not meet with 30-day IHM (16.7% vs 4.3%, p = 0.001). Additionally, a larger proportion of patients who met with 30-day IHM had a respiratory source of infection (45.5% vs 27.2%, p = 0.006) while a smaller proportion had a source of infection originating from the urinary tract (7.6% vs 25.7%, p = 0.003) when compared to patients who did not meet with 30-day IHM. No significant differences were detected in gender, PACS status, ethnicity, or medical history between both groups. Table 2 presents the descriptive analysis of HRV and HRnV parameters. In this study, N was set as 3 and HR 2 V, HR 2 V 1 , HR 3 V, HR 3 V 1 and HR 3 V 2 parameters were calculated. Among time domain parameters such as mean NN and SDNN, HR n V and HR n V m values are generally directly proportional to n and increase when n increases. HR 2 V SampEn and HR 3 V SampEn were considerably larger than SampEn parameters of HRV, HR 2 V 1 , HR 3 V 1 , and HR 3 V 2 . This was because of insufficient data points since our ECG recordings were only five minutes long. HR 2 V 1 , HR 3 V 1 and HR 3 V 2 did not encounter this limitation as more data points were available from a calculation using overlapping RR n I m sequences [20]. Table 3 shows the results of univariable analysis of HRV and HRnV parameters. Of 142 HRV and HRnV parameters, 85 were significantly different between the two outcome groups. Specifically, 14 HRV, 14 HR 2 V, 16 HR 2 V 1 , 11 HR 3 V, 16 HR 3 V 1 , and 14 HR 3 V 2 parameters were statistically significant. In at least four out of six HRnV measures, RMSSD, kurtosis, NN50, pNN50, NN50n, pNN50n, HF power, HF power norm, Poincare SD1, and Poincare SD1/SD2 were significantly higher, while LF power norm and DFA α2 were significantly lower in patients who met the primary outcome compared to those who did not. Additionally, VLF power and DFA α1 were not significant in HRV analysis but were statistically significant in several HRnV measures.

HRV and HRnV parameter description and univariable analysis
Overall, six baseline characteristics (age and vital signs at triage including temperature, respiratory rate, SpO 2 , SBP and GCS), 17 HRV parameters, and 96 HRnV parameters had p<0.2 on univariable analysis. After collinearity assessment, the remaining 87 variables were entered into a stepwise-selection regression model. Table 4 presents the multivariable analysis of variables found to be significantly different on univariable analysis. A total of 21 out of 87 variables were selected through stepwise selection.

Discussion
In recent years, there has been a surge in research interest in HRV and its ability to prognosticate for adverse patient outcomes across various disease processes [18,19,28]. To improve the predictive power of HRV, several studies have sought to utilise advanced nonlinear techniques   to derive novel HRV parameters [29,30]. Indeed, we previously employed the novel HRnV measures to assess the risk of 30-day major adverse cardiac events in patients presenting to the ED with chest pain [20]. This study evaluated the predictive value of novel HRnV measures (HR n V and HR n V m ) in estimating the risk of 30-day IHM in patients presenting to the ED with sepsis. In addition to the 22 traditional HRV parameters, we derived an additional 120 HRnV parameters, 71 of which were found statistically significant in their association with the primary outcome. The newly generated HRnV parameters greatly augment the number of candidate predictors and have demonstrated improved predictive ability for sepsis mortality. Although the physiological meaning and clinical interpretation of some HRnV parameters are yet to discover, the rich inherent information obtained from novel HRnV measures could create a new perspective for data scientists and machine learning researchers to investigate innovative approaches for ECG analysis and risk monitoring.
In HRnV measures, the newly added parameters, NN50n and pNN50n, were significantly associated with mortality in the univariate analysis. They characterise the number of times that the absolute difference between two successive RR n I sequences exceeds 50×n ms, by assuming the absolute difference may be magnified when the corresponding RR n I is n times longer than RRI [20]. The composite HRnV model derived from multivariable logistic regression achieved the highest AUC on ROC analysis and outperformed other established disease scoring systems such as NEWS, MEWS, SOFA, and APACHE II for the prediction of 30-day IHM in patients presenting to the ED with sepsis. In addition to demonstrating the superior predictive ability for sepsis mortality, the HRnV model is made even more relevant in its capacity for rapid and objective prognostication where only vital signs and parameters calculated from five-minute ECG tracings are needed. Many established disease severity scores require invasive tests, which need long turnaround time and resources to obtain or include subjective parameters that involve interrater variability while scoring. Among disease severity scores, the MEDS score, explicitly developed for risk stratification of septic patients in the ED, suffers from some of these limitations and its adoption has thus not been widespread. Consequently, MEDS was not included in our comparison. APACHE II and SOFA scores initially designed for use in the intensive care unit (ICU) setting similarly require invasive investigations to calculate its score. In these aspects, the HRnV model which is derived from vital signs taken on ED presentation, and HRV and HRnV parameters calculated from five-minute ECG tracings, can overcome these limitations and provide a rapid, objective, and accurate risk assessment of the septic patient. A triage tool with these characteristics would be invaluable to the physician and can aid in risk stratification, clinical management, patient disposition, and accurate patient classification for administrative or research purposes. Furthermore, our HRnV analysis and modelling modules can be readily integrated into a monitoring device, making the real-time prediction of sepsis severity a feasible task.

Limitations
First, a majority of the HRnV variables were removed from predictive modelling with traditional logistic regression, which hindered the release of the power of the novel HRnV representation. Moving forward, we endeavour to explore the use of black box or interpretable machine learning algorithms [37][38][39] for full utilisation of the HRnV parameters. Second, the difficulty in interpreting the HRnV parameters may pose a challenge in their clinical implementation and adoption. However, data scientists and machine learning practitioners may find these parameters valuable in data mining tasks. Third, we were only able to recruit a convenience sample of all suspected sepsis patients at the ED due to resource constraints and the difficulty of continuous ECG measuring in an emergency setting. Moreover, there was a selection bias in patient recruitment as we only included patients from PACS 1 or 2 triage category. To address this issue, we are planning a prospective study to include all ED sepsis patients. Last, sepsis is a seasonal illness that varies throughout the year, but we were unable to examine the impact of seasonality because Singapore's weather is warm and humid all year round.

Conclusions
The use of novel HRV measures (HR n V and HR n V m ) can provide additional power to predictive models in the risk stratification of patients who present to the ED with sepsis. When included in a model with other clinical variables, the HRnV model outperforms traditional risk stratification scoring systems as shown in our preliminary results. Prospective multi-centre cohort studies would be valuable in validating the effectiveness of the HRnV parameters. The use of HRnV may allow for a rapid, objective, and accurate means of patient risk stratification for sepsis severity and mortality.