Clinical potential of pupillary light reflex and heart rate variability parameters as biomarkers for assessing pain relief effects on autonomic function: a prospective longitudinal study

Patients with acute inflammatory disorders in the otolaryngology—head and neck region, treated at the Hyogo Prefectural Amagasaki General Medical Center between August 2018 and March 2019, were enrolled in a prospective longitudinal study.

The acute inflammatory disorders were defined as acute tonsillitis, peritonsillar abscess, acute epiglottitis, acute sinusitis, and deep neck space abscess. Inclusion criteria were patients aged older than 17 years, requiring hospitalization for medical treatment of acute inflammatory disorders in the otolaryngology—head and neck region. Experienced otolaryngologists (more than 10 years' experience) judged the necessity of the hospitalization if the following symptoms and findings were present: dysphagia due to pain; abscess requiring surgical drainage; potential of sudden respiratory obstruction. We collected the following sociodemographic and medical data: sex; age; body mass index (BMI); alcohol consumption; smoking intensity; asthma status; diabetes status, by means of questionnaires. BMI was categorized as follows: (1) Underweight and healthy, BMI < 25; (2) Overweight and obese, 25.0 $\leqq$ BMI. Based on alcohol consumption, individuals were classified as follows: (1) Non-drinker; (2) Light/Moderate drinker; (3) Heavy drinker. Light/Moderate drinker was defined as drinking one (included) to five drinks per day, while heavy drinker as drinking more than five drinks per day. Smoking intensity based on Brinkman index (number of cigarettes smoked per day multiplied by the number of years of smoking),

All patients gave informed consent prior to initiating treatment. After performing pain assessments using subjective and objective tools, patients received treatments, mainly by way of intravenous and oral antibiotic agents. When diagnostic imaging showed an abscess in the tonsils or neck, the patients underwent abscess drainage prior to antibiotic treatment. All patients were hospitalized for at least 3 days. As they obtained sufficient symptomatic relief, patients were discharged from hospital, with subsequent ambulant follow-up. All patients were followed-up from day 7 to day 28 post-treatment.

This prospective study was approved by the Research Ethics Committee of Hyogo Prefectural Amagasaki General Medical Center.

We used portable infrared pupillometer (PLR-3000, NeurOptics, Irvine, California, USA) and HRV devices (CheckMyHeart, Daily Care BioMedical, Taiwan) for the PLR and HRV tests, respectively.

In PLR tests, we determined the following eight parameters: (1) initial pupil size (INIT, mm); (2) minimum pupil size after the light stimulus (END, mm); (3) averaged constriction velocity (ACV, mm s⁻¹); (4) maximum constriction velocity (MCV, mm s⁻¹); (5) constriction ratio, defined as the difference between INIT and END divided by INIT (DELTA, %); (6) latency for constriction (LAT, ms); (7) averaged dilation velocity (ADV, mm s⁻¹); (8) time taken for 75% recovery of INIT (T75, s). PLR parameters are markedly affected by the stimulus intensity range (Ellis 1981, Gavriysky 1991, Bitsios et al 1996, Ishikawa 2020). To assess how a range of stimulus intensity can affect PLR parameters, we prepared four different stimulus intensities (10, 50, 121, and 180 μW) by a transient flashlight with 800 ms pulse duration. We set the stimulation duration of three seconds to 10 μW, while that of five seconds to the other stimulus intensities. To measure parameters in response to each stimulus intensity, we adopted a step-up method to both eyes, as reported previously (Ishikawa 2020).

HRV data were recorded from both forearms for five minutes, and then analyzed automatically by the device, using the built-in software provided by the manufacturer. For time domain analysis, we obtained three detrend parameters of standard deviation of normal R-R intervals (SDNN), root mean square of successive R-R interval differences (rMSSD), pNN50 (percentage of adjacent R-R intervals that differed from each other by more than 50 ms). Normal values of HRV parameters obtained from analysis of the time domain are known to have a universal exponential decay-like relationship with heart rate, and can result in misleading conclusions about scientific findings (Monfredi et al 2014). Therefore, it is recommended to correct the values for heart rate level (Monfredi et al 2014, van Roon et al 2016, van den Berg et al 2018). We obtained corrected SDNN (SDNNc) and rMSSD (rMSSDc) using a correction procedure reported previously (van den Berg et al 2018). For frequency domain analysis, we obtained 5 detrend parameters, quantified with an autoregressive method: (1) the power of low frequency HRV (LF power; 0.04 to 0.15 Hz); (2) the power of high frequency HRV (HF power; 0.15 to 0.4 Hz); (3) normalized values of low frequency HRV (LF norm; LF power × 100/Total power—very low frequency); (4) normalized values of high frequency HRV (HF norm; HF power × 100/Total power—very low frequency); and (5) LF/HF values. Previous studies pointed out the disadvantage of raw values (LF power and HF power) because of the large individual variability, which leads to a long-tailed right-skewed exponential distribution (Burr 2007, Ishikawa 2020). Finally, we adopted SDNN, SDNNc, rMSSD, rMSSDc, pNN50, HFnorm, LFnorm, and LF/HF for HRV analysis.

Four time-points were set for assessment during the time course of treatments: (1) pre-treatment; (2) day 1 post-treatment; (3) day 2 post-treatment; (4) last follow-up during the days 7–28 post-treatment period. The last follow-up was defined as the final consultation during the follow-up term.

For the assessments, all patients completed questionnaires on pain intensity, followed by PLR and HRV tests. In the questionnaire, the NRS, ranging from 0 to 10, was used as a general index of pain intensity, with 0 representing no pain and 10 representing the worst imaginable pain.

We performed PLR and HRV tests in a quiet consulting room. Patients underwent the HRV test in the dim room, followed by the PLR test. The data collection protocol was time-sensitive, particularly at the pre-treatment and at the last follow-up. Pre-treatment measurements were performed between 8 am and 11 pm, while those at the last follow-up were performed between 9 am and 3 pm. At days 1 and 2 post-treatment, measurements were performed between 8 am and 11 am.

We expressed variables as means and 95% confidence intervals (CI) or as medians and interquartile ranges (IQR). We compared NRS data over the time course using a Friedman test with Dunn's multiple comparison test.

To explore the association between NRS and PLR/HRV parameters, we used four LMMs. LMMs consist of fixed or random effects. In fixed effects of HRV data, we specified predictor variables in the following: NRS and sociodemographic factors such as age, sex, categorized BMI, alcohol consumption, and smoking intensity. Random effects consisted of random intercept and slope. In fitting only a random intercept, we assumed values of PLR/HRV parameters can differ among subjects; however, the alteration of the parameters for a one unit change in NRS can be common: different intercepts with common slopes among subjects. In fitting random intercepts and slopes, we assumed that the alteration of PLR/HRV parameters to the NRS alteration can differ among subjects: different intercepts and slopes among subjects. Model 1 consisted of NRS. In Model 1 of PLR parameters, the stimulus intensity was added to fixed effects. Model 2 contained all the terms of Model 1, with an additional term for sociodemographic factors such as age, sex, categorized BMI, alcohol consumption, and smoking intensity. Model 3 (defined as random intercept model) added random intercept to all the terms of Model 2. Model 4 (defined as random slope model) added random slope to all the terms of Model 3. The NRS was specified as the random slope by adding it to the structure. We fit subject ID (n = 35 subjects for HRV data) or subject's eye ID (n = 70 eyes of 35 subjects for PLR data) as a random intercept. We compared values of Akaike's information criterion (AIC) among the four models. The better-fit models were defined as the lower AIC model. A two-sided p-value < 0.05 was considered statistically significant. All statistical analyses were performed with R in the statistical analyses (version 3.6.1, R Foundation for Statistical Computing, Vienna, Austria, lme4 and lmerTest package). The formula of the four models are shown using notation equivalent to fitting models in the lme4 package code.

Model 1: ANS parameter values = Intercept + NRS (+ stimulus intensity in PLR)

Model 2: ANS parameter values = Model 1 +sociodemographic factors

Model 3: ANS parameter values = Model 2 +(1 ∣ subject ID in HRV or subject's eye ID in PLR)

Model 4: ANS parameter values = Model 3 + (1 + NRS ∣ subject ID in HRV or subject's eye ID in PLR)

Thirty-five Japanese patients were enrolled in the current study. Table 1 shows the sociodemographic and medical characteristics of the patients. There were no patients suffering from asthma and diabetes, and lost to follow-up.

We compared the NRS scores along the 4 time-points. The median values (IQR) were 8 (5, 9), 5 (3, 6), 2 (5, 1), and 0 (0, 0) at pre-treatment, Day 1, Day 2, and last follow-up, respectively. Statistical analysis revealed a significant decrease in the NRS over time (p < 0.001). A post-hoc test revealed a significant decrease of NRS values at last follow-up in comparison to values at pre-treatment (p < 0.001).

In PLR data, more than 32% of the T75 data were unanalyzable because of an insufficient test duration to recover 75% of INIT, and therefore we excluded the parameter from the analysis. In Model 2 of LMMs, we observed collinearity between age and smoking intensity through the model selection process. Therefore, we excluded the latter from sociodemographic factors.

We compared AIC among the four models in PLR/HRV parameters (table 2). All parameters except LF/HF showed the lower AIC due to adding random effects. Based on the lowest AIC, we selected, as better-fit models, Model 3 in all HRV parameters except LF/HF, Model2 in LF/HF, and Model 4 in all PLR parameters.

Table 3 shows how regression estimate and 95% CI can differ among the four models when we expected change in values of PLR parameters for one microwatt change in stimulus intensity. In all PLR parameters, we observed same regression estimates among four models, and narrower 95% CI due to adding predictor variables such as sociodemographic factors and random effects. Model 4, the better-fit model, showed all significant PLR parameter shifts due to an increased stimulus intensity in the following: (1) the decrease of INIT, END, and ADV, (2) the shortening of LAT, and (3) the increase of ACV, MCV, and DELTA.

In PLR analysis, we observed significant associations between NRS and INIT, END, MCV, LAT, and ADV (table 4). One unit decrease in NRS was significantly associated with a mean increase of 2.4 × 10⁻² mm in INIT (figures 1(a)) and 2.1 × 10⁻² mm in END (figure 1(b)), mean decrease of 2.8 × 10⁻² mm s⁻¹ in MCV (figure 1(c)) and of 1.3 × 10⁻² mm s⁻¹ in ADV (figure 1(d)), and mean prolongation of 78.6 ×10⁻² ms in LAT (figure 1(e)).

Zoom In Zoom Out Reset image size

In HRV analysis, we observed significant associations between NRS and pNN50, HF norm, LF norm, and LF/HF (table 5). Decreased one unit in NRS caused mean increase of 0.56% in pNN50 (figure 1(f)) and of 1 in HF norm (figure 1(g)), while mean decrease of 1 in LF norm (figure 1(h)) and of 0.1 in LF/HF (figure 1(i)).