I feel like the rabbit in Alice in Wonderland

And you are all Alice…..

Basic idea of residual HR

High frequency power (HFP) of heart rate variability (HRV) is associated with parasympathetic cardiac regulation (PNS). HF has a stronger relationship to parasympathetic cardiac regulation than any other HRV index does to sympathetic regulation (SNS) (empirically demonstrated). In addition to the empirical demonstration this makes theoretical sense for a few reasons. First, because the PNS neurotransmitter acetylcholine (ACTH) can rapidly (no second messenger needed) impact the SA node membrane potential. Second, because the PNS connection to the SA node is neural (not hormonal) and thus even when “tonic” tends to be “tonic” via “phasic” activity. Third, because the SNS neurotransmitters epinephrine and norepinephrine require a second messenger to impact the membrane potential. Fourth, because the SNS neurotransmitters are also hormones, so there are two routes by which SNS impacts the membrane potential and one of them is truly “tonic”, not “tonic” via “phasic”.

HR is one of the end organ responses to whatever the balance may be between PNS and ANS tonic and phasic activity to the heart.

In essence, \(HR \propto SNS \times \frac{1}{PNS}\) (note this is being used conceptually and does not imply that this is the exact equation that would be used)

Solving for SNS we get: \(SNS \propto HR \times PNS\)

The idea is that if we have three related variables and we know two of them we can figure out the one unknown.

Empirical approach to residual HR

The empirical approach replaces \(PNS\) with \(LnHFP\) and \(HR\) with \(RR\) interval due to the linear relationship between \(RR\) and \(LnHFP\).

\(RR \propto \frac{1}{SNS} \times PNS\)

An empirical approach to this is required due to the fact that the actual relationship depicted with \(RR \propto \frac{1}{SNS} \times PNS\) very likely varies between people, and there is no “equation” with this relationship because we do not have units for “PNS” and “SNS” in this relationship (and I’m confident many additional arguments could be made).

One possible empirical approach is to use a regression to relate an individuals HF (as a measure of PNS) to their HR. This individualized regression model then provides a residual value for each datum that is how much higher, or lower, the RR is compared to the regression model’s prediction based on the HF. If the RR is higher than the HF would predict it is because there is less SNS activity, if the RR is lower than the HF would predict it is because there is more SNS activity.

\(RR = \beta_0 + \beta_{LnHFP} \times LnHFP + \epsilon\)

In that regression formula \(\epsilon\) is the error term, which is the residual, which we then consider to be proportional to SNS, \(SNS \propto \epsilon\)

\(RR \propto \frac{1}{\epsilon}\times LnHFP\)

Data Processing (draft)

Raw actiheart data prepared for Kubios

James (Karin’s colleague) wrote a script to transform the old Actiheart text files into single columns of RR interval data

HRV analysis in Kubios

Karin reads those files into Kubios, Time-varying analysis parameters, etc, exports a file for each subject with their estimated HF power and HR for each minute that they were wearing the monitor that is centered within a 5 minute window

Another way of saying this - the window of analysis is 5 minutes, and that window slides down the raw data by minute, so we have a 5 minute HR and HRV analysis every minute.

Exploration into the approach - “Consideration of concept”

A challenge with this analysis will be coming up with what we’re actually looking for…. but perhaps this is something you’re used to with your line of research in MECFS

Subject FL101

Day 1

# This section will need to be changed drastically with additional files

d1<-read.csv("~/Dropbox/MECFS/HRV-Data-Analysis/tv_hrv_full/FL101d1_hrv.csv")
# d1$Time <- lubridate::hms(d1$Time)

d1 <- mutate(d1, lnHFP = log(HF.power))
d1 <- mutate(d1, Beats.corrected.percent = (Beats.corrected / Beats.total)*100)
d1 <- select(d1, ID, CPET, Beats.corrected.percent, Mean.RR, lnHFP)
# This section will need to be changed drastically with additional files

d2<-read.csv("~/Dropbox/MECFS/HRV-Data-Analysis/tv_hrv_full/FL101d2_hrv.csv")
# d2$Time <- lubridate::hms(d2$Time)

d2 <- mutate(d2, lnHFP = log(HF.power))
d2 <- mutate(d2, Beats.corrected.percent = (Beats.corrected / Beats.total)*100)
d2 <- select(d2, ID, CPET, Beats.corrected.percent, Mean.RR, lnHFP)
# Combine d1 and d2 before the individual models are built

combined <- rbind(d1, d2)
d1_plot <- ggplot(d1, aes(lnHFP, Mean.RR)) +
  geom_point() +
  geom_smooth(method = lm)
d1_plot

Before we proceed we should consider what we’re looking at with these “Residual HR” models.

d1_model <- lm(Mean.RR ~ lnHFP, d1)
summary(d1_model)
## 
## Call:
## lm(formula = Mean.RR ~ lnHFP, data = d1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -101.864  -16.062    8.306   29.469   67.709 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  369.691     19.674   18.79  < 2e-16 ***
## lnHFP         60.658      4.382   13.84 2.62e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 45.7 on 29 degrees of freedom
## Multiple R-squared:  0.8685, Adjusted R-squared:  0.864 
## F-statistic: 191.6 on 1 and 29 DF,  p-value: 2.617e-14
d1$ResHR <- d1_model$residuals
plot(d1$Mean.RR)

plot(d1$lnHFP)

plot(d1$ResHR)

resHF1_plot <- ggplot(d1, aes(lnHFP, ResHR)) +
  geom_point() +
  geom_smooth(method = lm)
resHF1_plot

resHR_plot <- ggplot(d1, aes(Mean.RR, ResHR)) +
  geom_point() +
  geom_smooth(method = lm)
resHR_plot

Day 2

d2_plot <- ggplot(d2, aes(Mean.RR, lnHFP)) +
  geom_point() +
  geom_smooth(method = lm)
d2_plot

d2_model <- lm(Mean.RR ~ lnHFP, d2)
summary(d2_model)
## 
## Call:
## lm(formula = Mean.RR ~ lnHFP, data = d2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -80.42 -20.40  12.52  29.27  59.68 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  409.019     15.751   25.97  < 2e-16 ***
## lnHFP         53.872      3.394   15.87 7.74e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 41.64 on 29 degrees of freedom
## Multiple R-squared:  0.8968, Adjusted R-squared:  0.8932 
## F-statistic:   252 on 1 and 29 DF,  p-value: 7.739e-16
d2$ResHR <- d2_model$residuals
plot(d2$Mean.RR)

plot(d2$lnHFP)

plot(d2$ResHR)

resHF2_plot <- ggplot(d2, aes(lnHFP, ResHR)) +
  geom_point() +
  geom_smooth(method = lm)
resHF2_plot

resHR2_plot <- ggplot(d2, aes(Mean.RR, ResHR)) +
  geom_point() +
  geom_smooth(method = lm)
resHR2_plot

Compare Days

FL101_df <- rbind(d1, d2)
FL101_df$CPET <- as.factor(FL101_df$CPET)
PNSplot <- ggplot(FL101_df, aes(lnHFP, Mean.RR, color = CPET)) +
  geom_point() +
  geom_smooth(method = lm)
PNSplot

SNSplot <- ggplot(FL101_df, aes(ResHR, Mean.RR, color = CPET)) +
  geom_point() +
  geom_smooth(method = lm)
SNSplot

ANSplot <- ggplot(FL101_df, aes(lnHFP, ResHR, color = CPET)) +
  geom_point() +
  geom_smooth(method = lm)
ANSplot

Need a sufficiently high level metric to capture the possible different response dynamics between individuals. Which is a challenge without knowing what that possible different response dynamics may be….

So - what could they be and in turn what would we expect?

For example, if there was a difference between CPET 1 and 2 for the above plot?

A sobering thought - there is no difference between days because the models have been constructed for each day. What if the model is created using data from both tests?

Start over

A combined regression model (meaning, regression and residuals of all data)

combined$CPET <- as.factor(combined$CPET)
combined_plot <- ggplot(combined, aes(lnHFP, Mean.RR, color = CPET)) +
  geom_point() +
  geom_smooth(method = lm)
combined_plot
## `geom_smooth()` using formula 'y ~ x'

combined_model <- lm(Mean.RR ~ lnHFP, combined)
summary(combined_model)
## 
## Call:
## lm(formula = Mean.RR ~ lnHFP, data = combined)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -108.303  -24.196    7.412   27.676   65.305 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  391.569     12.460   31.43   <2e-16 ***
## lnHFP         56.721      2.729   20.79   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 43.94 on 60 degrees of freedom
## Multiple R-squared:  0.8781, Adjusted R-squared:  0.876 
## F-statistic:   432 on 1 and 60 DF,  p-value: < 2.2e-16
combined$ResHR <- combined_model$residuals
ANS_Combinedplot <- ggplot(combined, aes(lnHFP, ResHR, color = CPET)) +
  geom_point() +
  geom_smooth(method = lm)
ANS_Combinedplot
## `geom_smooth()` using formula 'y ~ x'

ANS_CombinedModel <- lm(ResHR ~ lnHFP + CPET + lnHFP*CPET, combined)
summary(ANS_CombinedModel)
## 
## Call:
## lm(formula = ResHR ~ lnHFP + CPET + lnHFP * CPET, data = combined)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -101.86  -17.87   10.41   30.05   67.71 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -21.877     18.821  -1.162    0.250
## lnHFP          3.937      4.192   0.939    0.352
## CPET2         39.327     25.054   1.570    0.122
## lnHFP:CPET2   -6.786      5.502  -1.233    0.222
## 
## Residual standard error: 43.71 on 58 degrees of freedom
## Multiple R-squared:  0.04318,    Adjusted R-squared:  -0.006306 
## F-statistic: 0.8726 on 3 and 58 DF,  p-value: 0.4606