Nonlinear heart rate variability biomarkers for gastric cancer severity: A pilot study

Identifying prognostic factors by affordable tools is crucial for guiding gastric cancer (GC) treatments especially at earlier stages for timing interventions. The autonomic function that is clinically assessed by heart rate variability (HRV) is involved in tumorigenesis. This pilot study was aimed to examine whether nonlinear indices of HRV can be biomarkers of GC severity. Sixty-one newly-diagnosed GC patients were enrolled. Presurgical serum fibrinogen (FIB), carcinoembryonic antigen (CEA), and carbohydrate antigen 19-9 (CA199) were examined. Resting electrocardiogram (ECG) of 5-min was collected prior to surgical treatments to enable the HRV analysis. Twelve nonlinear HRV indices covering the irregularity, complexity, asymmetry, and temporal correlation of heartbeat fluctuations were obtained. Increased short-range temporal correlations, decreased asymmetry, and increased irregularity of heartbeat fluctuations were associated with higher FIB level. Increased irregularity and decreased complexity were also associated with higher CEA level. These associations were independent of age, sex, BMI, alcohol consumption, history of diabetes, left ventricular ejection fraction, and anemia. The results support the hypothesis that perturbations in nonlinear dynamical patterns of HRV predict increased GC severity. Replication in larger samples as well as the examination of longitudinal associations of HRV nonlinear features with cancer prognosis/survival are warranted.

The fact of late-stage presentation and inaccessible treatment is urging an early diagnosis of cancer, the second leading cause of death worldwide. Such a malignancy spreads equally without preference on human beings all over the world but gastric cancer (GC) acts rather eccentrically that has been pushed out as an exception. It is so common in China and other East Asia countries as well, ranking the second in cancer death as opposed to the fifth globally 1 . To improve GC prognosis and facilitate a better treatment planning, early and sensitive diagnosis with feasible and affordable clinical measurements is essential.
Converging evidence has suggested a pivotal role of the autonomic control in tumor progression, in particular the contribution of the vagal nerve activity through many tumor-inhibiting mechanisms 2,3 . As a clinical routine, the measurement of electrocardiogram (ECG) or more specifically the analysis of the beat-to-beat ECG RR interval variations -the heart rate variability (HRV) -is an optimal noninvasive biomarker for the autonomic regulation 4,5 . In previous studies, reduced HRV was found in cancer patients compared to healthy peers 6 . Lower HRV at baseline was also reported to predict increased carcinoembryonic antigen (CEA) months later in a historical-prospective study 7 . In addition, in a most recent systematic review on HRV and cancer prognosis 8 , 19 studies that involved various kinds of cancer patients were included and appraised. Regardless of the cancer types, this review concluded an adverse effect of lower HRV towards shorter survival, higher tumor burden, or more advanced metastasis stage. Consistently, in another recent clinical study of GC patients 9 , lower HRV was found to be associated with advanced clinical stage, increased tumor size, tumor infiltration, lymph node metastasis, and involvement of distant organs.
Surprisingly, all studies reviewed above used only traditional linear HRV measures, albeit a commonly accepted nonlinear nature of HRV 10,11 . It is considered highly complex owning to the competition between spontaneity and adaptability of the heart beat regulation. Across a variety of studies in the field of cardiovascular diseases, nonlinear dynamical HRV analysis has shown a tremendous advantage over these linear time-and Figure 1 shows examples of RR interval time-series that illustrate the construction of RR interval time-series from ECG without ectopic beats and with ectopic beats, respectively. Table 1 summarizes the demographics and the clinical and HRV measures of patients. Pearson correlation analyses resulted in six nonlinear HRV features, i.e., FuzzyEn, PermEn, DistEn, PI, α 1 , that had significant (p < 0.1) correlations with at least one of the three clinical GC parameters (Table 2). Based on these results, linear regression models of FIB with separately FuzzyEn, PermEn, PI, or α 1 were performed; linear regression models of CEA with PermEn, DistEn, or PI separately were performed; linear regression models of CA199 with PI or GI were separately examined. Results from linear regressions were summarized in Table 3.

Results
After correcting for multiple comparisons, α 1 , PI, and PermEn were significantly associated with FIB, specifically, FIB increased by 0.41 ± 0.10, −0.35 ± 0.10, and 0.30 ± 0.11, respectively, for each 1-SD increase in α 1 , PI, and PermEn (all false discovery rate [FDR]-corrected p < 0.05). PermEn and DistEn were significantly associated with CEA, specifically, CEA increased by 0.36 ± 0.15 and −0.32 ± 0.14, respectively, for each 1-SD increase in PermEn and DistEn (both FDR-corrected p < 0.05). These five significant associations are further explained by the partial correlation plots as shown in Fig. 2 (the corresponding correlation plots without adjust for demographics were shown in Fig. S1 documented in the online Supplemental Materials). No parameters were significantly associated with CA199. After further adjusting for BMI, alcohol consumption, history of diabetes, Hb, and LVEF, all these associations still held with slightly changes in the estimated coefficients as shown in Table 4.   www.nature.com/scientificreports www.nature.com/scientificreports/ The associations of FIB with α 1 , PI, and PermEn were also significant in secondary analysis using logistic regression with dichotomized outcomes, specifically, the odds of having higher FIB increased by 168% ( Table 5). The associations of CEA with PermEn and DistEn became not significant in the logistic regression models. However, PI showed significant associations with CA199 with an odds ratio of 0.53 (95% CI: [0.27, 0.95]) for 1-SD increase in PI (p = 0.03; Table 5).

Discussion
With 61 pathologically-diagnosed GC patients in this pilot study, for the first time we demonstrated significant associations between clinical cancer markers and several nonlinear HRV measures after accounting for multiple comparisons. Specifically, the increase of short-range temporal correlations in heartbeat fluctuations (i.e., increase in α 1 which was calculated within time scales 4-16 beats), the decrease of the asymmetry in heartbeat acceleration/deceleration patterns (i.e., PI), and the increase of the irregularity of heartbeat fluctuations (i.e., PermEn) were associated with higher serum FIB level. The increase in PermEn as well as the decrease of the complexity of heartbeat dynamics (i.e., DistEn) were also associated with higher serum CEA level. Importantly, these associations were independent of several potential confounding factors including age, sex, BMI, alcohol consumption, history of diabetes, Hb, and LVEF.  Table 3. Results from linear regression models (adjusted for age and sex). * Effects for 1-standard deviation increase in the predictor adjusted for covariates. Abbreviations: CA199 = carbohydrate antigen 19-9; CEA = carcinoembryonic antigen; DistEn = distribution entropy; FDR: false discovery rate; FIB = fibrinogen; FuzzyEn = fuzzy entropy; GI = Guzik's index; PermEn = permutation entropy; PI = Porta's index; SE: standard error.
Re{PermEn ~ Age + Sex}  22 . There is mounting evidence supporting the vagal neuromodulation hypothesis in tumorigenesis through its effects of anti-inflammation, antioxidative stress, and sympathetic activity 2,3,23 including animal studies 24-26 that demonstrated a causal pathway.
Through the analysis of vagal and sympathetic modulation on heartbeat, HRV is a well-known and widely-applied noninvasive tool for assessing the autonomic nervous function. Increased HRV has consistently been associated to better prognosis in cancer patients 7,8,22,[27][28][29][30] . There is also an increasing preference in the biomedical science/engineering communities of using nonlinear analysis approaches as complements to the traditional linear methods. Coming from the fields of statistical physics and nonlinear dynamics, these nonlinear approaches could uniquely capture the information content (i.e., entropy-based measures), asymmetry, or scaling invariant properties (i.e., DFA), all of which has been shown to offer additional, valuable knowledge to the underlying control mechanism, i.e., the autonomic regulation 12,31-35 . Nonlinear HRV measures for the organisms' plasticity and adaptability coping with stress.
Although being not immediately interpretable with regards to the vagal or sympathetic regulation, there are published pilot studies that have already explored the effects of vagal and sympathetic outflows on several nonlinear measures of heartbeat dynamics including entropy measures and DFA scaling exponents 36,37 . Using both human and animal models, they have offered the direct evidence of autonomic control influencing the complex behavior of the heart. A more traditional or systemic level viewpoint is that the nonlinear behavior of HRV is attributed to the competing regulation on the heart coming from the two branches of the autonomic nervous system and the spontaneity of the organism itself. Such competition renders healthy organisms high complexity, enabling daunting plasticity and adaptability to the stresses/perturbations of everyday life 38,39 . In parallel with this complex  Table 5. Results from Logistic regression models (adjusted for age and sex). Results presented in the same order as in Table 3 physiology hypothesis, aging and disease progressions are usually accompanied by a progressive reduction of the complexity 12,34,38,40 ; and the other way around the degradation of the complexity also predicts future incidence of disease 41 .
In keeping with the complexity loss theory, our results suggest that GC patients with worse prognosis showed lower autonomic control complexity even though each different nonlinear metric showed different changing directions. Theoretically, the most complex system should be neither too random nor too regular; it should correspond to a critical point in-between 42 . The departure of α 1 from the value 1, whichever direction, both imply a reduced complexity 43 . In our case, it was a reduction towards the regular side (i.e., increase of α 1 towards 1.5). The decrease of PI indicates a loss of time irreversibility, an important property of complex system which suggests an evolution of the system to equilibrium or a loss of hysteresis 44 . The increase of PermEn indicates higher irregularity which suggests a seemingly controversial behavior as compared with α 1 . However, the calculation of PermEn focused on the fluctuation motifs composed of 3 heart beats that were not included in the calculation of α 1 for sake of a robust fitting. The decrease of DistEn, although it is in nature an entropy metric, directly suggests a decrease in complexity as evidenced by the simulation analysis in the original DistEn study 45 . potential clinical relevance and usefulness. With the rising of global GC epidemic, the assessment of HRV may potentially meet the clinical urgency in three ways: (1) HRV analysis may offer a sensitive and noninvasive tool targeting an early GC diagnosis. On one hand, the nonlinear approaches used in this study are necessary in the way that they cope well with the nonlinear and nonstationary nature, resulting thus in a more robust assessment as compared with the existing linear methods. On the other hand, it might be possible to leverage these nonlinear indices with the existing linear measures to construct an integrated biomarker for early GC diagnosis. Further studies with larger samples are required to test this hypothesis. (2) HRV could help with the evaluation of prognosis and treatment planning for GC patients. The link between HRV and survival time suggests a role of HRV in helping screen the general health status and prognosis of cancer patients 8,29,30 . Previous studies also discovered a link between serum FIB and adjacent organ involvement 18 . Given the strong associations of serum FIB with HRV nonlinear indices reported in the current study, HRV analysis may thrive as a sensitive tool for the surgical planning of GC patients. Further clinicopathological studies are warranted to formally examine their associations. (3) HRV might be a target for interventions to prevent the disease or to slow down the progression. Although only cross-sectional associations were reported in this current study, previous animal studies have established a causal link between vagal activity and tumor genesis [24][25][26] . A case-control study that focuses on exposure to interventions improving HRV especially in terms of nonlinear properties is required to validate this causal pathway in humans.

Study limitations.
There are several notable limitations. First, the sample size is relatively small. Aside from the four indices (i.e., α 1 , PI, PermEn, DistEn) the other nonlinear HRV measures may also be correlated with the serum markers while their negative observations may simply be due to the power issue. Second, the cross-sectional nature of the study design limits our inference about longitudinal prediction ability of HRV nonlinear measures. It is of great clinical value to examine whether those presurgical HRV indices can predict longer term outcomes including treatment response and survival. Besides, it is also meaningful to check whether degradations of these nonlinear properties in otherwise normal people predict higher risk of developing GC later. Third, although the FIB, CEA, and CA199 are well developed serum markers of cancer severity or prognosis, the gold standard is pathological examination. We are still working to retrieve and pool detailed pathology data together in addition to the final diagnosis and are expecting to scrutinize their interrelationships with HRV nonlinear features in follow-up studies. Serum FIB, CEA, and CA199 levels were examined before breakfast 1-week before surgical treatments. FIB levels were determined using the Clauss method (Sysmex CS51000, Sysmex Corporation, Kobe, Japan). Chemiluminescent assays were used to determine the CEA and CA199 levels (Architect i2000sr, Abbott Diagnostics, Abbott Park, IL, USA). ECG data were recorded continuously for 5 min one day before treatments with patients lying down for at least 20 min before the formal collection (HeaLink-R211B, HeaLink Ltd., Bengbu, China). The sampling frequency of ECG collection was 400 Hz. The precordial V5 lead was configured and the Ag/AgCl disposable electrodes were used (Junkang Ltd., Shanghai, China).

Materials and
We further excluded participants with the following conditions: (1) recurrent GC (N = 1), (2) poor ECG quality (N = 2), (3) presence of ectopic beats (>10% of all beats; N = 8), or (4) administering of blood transfusion (N = 15) or chemotherapy (N = 3) prior to ECG collection. Therefore, data of 61 participants were presented and analyzed in this current study. Their demographics and clinical characteristics were shown in Table 1.
www.nature.com/scientificreports www.nature.com/scientificreports/ Nonlinear HRV analysis. R-wave peaks were located automatically using a template-matching approach 46 followed by visual inspections. Ectopic beats were identified manually during the visual inspection (1-3 ectopic beats, either premature atrial contraction or premature ventricular contraction, presented in 14 out of the 61 patients). The final HRV time-series were constructed for each patient by consecutive normal sinus R-R intervals. The normal to ectopic or ectopic to normal intervals were discarded and the corresponding two segments were stitched together to assure a reasonable length of the RR interval time-series. Figure 1 shows examples of the construction of RR interval time-series without and with ectopic beats.
The extraction of R-peaks, visual inspections of ectopic beats, and asymmetry analysis were done in MATLAB (Ver. R2018a, The MathWorks Inc., Natick, MA, US). Entropy analysis was performed using the EZ Entropy software 47 . The Kubios HRV software was used to perform the DFA 48 .
Statistical analysis. Histograms of CEA and CA199 both showed an obvious right skewness; they were thus natural-log-transformed prior to further analysis (unless otherwise indicated). Bivariate Pearson correlations of FIB, CEA, and CA199 with each of the 12 nonlinear HRV measures were performed to screen potential predictors. Conservatively, features with p level of <0.1 were considered significant. Linear regressions were then performed for these pairs that passed this screening process with either FIB, CEA, or CA199 as an outcome. These models were adjusted for age and sex. To avoid collinearity, the HRV features were each included in a separate regression model as a predictor. To determine the statistical significance, the Benjamini-Hochberg procedure was used to control for false discovery rate (FDR) within multiple comparisons 49 . FDR-corrected p < 0.05 is considered statistically significant. These models were then augmented by further adjusting for BMI 50 , alcohol consumption 51 , history of diabetes 52 , LVEF 53 , and anemia as assessed by Hb 54 followed by the same multiple comparison correction procedure. As an exploratory analysis, FIB, CEA, and CA199 were dichotomized each by a separate threshold value 19 : (1) FIB was considered high if >3.5 mg/mL and low if otherwise; (2) CEA was considered high if the original CEA level ≥ 5 and low if otherwise; and (3) CA199 was considered high if the original CA199 level ≥ 37 and low if otherwise. Logistic regression models were performed with each of the three dichotomized clinical parameters as an outcome and with each feature in the corresponding significant HRV feature set as a predictor while adjusted for age and sex. All statistical analyses were performed using the JMP Pro (ver. 14, SAS Institute, Cary, NC).

Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.