Heart Rate Variability

Our heart rate varies from moment to moment. That variability is driven by many different causes and mainly aims to ensure homeostasis – that the body’s stability is preserved in the face of fluctuating environmental demands.

When we think of variability, we think naturally of measuring it with variance or more usually its square root – standard deviation. That is often sufficient for heart rate variability. We can record someone’s heart rate for 5 minutes and calculate the sd. The higher the value, the more variable the heart rate has been over the 5 minutes.

However, these two 5 minutes recordings of heart rate have the same sd, but are obviously quite different.

The left one has slowly varying heart rate, going through a complete cycle of speed up and slow down in 120 seconds; the right one has much faster variation with a cycle period of 30 seconds. Technically these are sine waves and we talk about their period (120s and 30s). We also talk about their frequency: how many cycles per second. The left one has a frequency of 0.5 cycles per second; the right one is 2 cycles per second or 0.5Hz and 2Hz, to use the proper unit.

Now look at this third one, which has both the slow component and the fast component:

In practice, heart rate recordings have many different components, varying over different time scales, from vey rapid moment to moment changes, to a diurnal pattern of sleeping, waking, eating etc. An example is shown next: it is easy to get a sense by eye that there are smaller fast variations and larger, slower variations. These variations with different time scales might well have very different causes and therefore different meanings.

That poses two issues for an analysis of heart rate variations:
1. How do we disentangle the different time-scales?
2. Which time-scales are we interested in?

Disentangling Time-Scales: Fourier Transform

Now we have the most technical part. Easy steps, but skip to the end if you prefer.

1. We start with an observation about any infinitely long and periodic signal – just a signal that keeps on repeating exactly the same for ever.
   The observation is that we can make an exact copy of that signal by adding together a set of sine waves with different frequencies (like the two at the top) in appropriate amounts.
2. And, happily for us, we can also reverse that process: from the signal itself we can make the large set of individual sine waves that we would require.
3. We never have an infinite periodic signal. But we can pretend that we do. We just imagine that the signal we have got repeats itself for ever. There are a few things we need to do to ensure that the process doesn’t then introduce any components that are not really there.
4. The process to obtain the set of sine waves from a signal is called the Fourier Transform.

The next figure shows the Fourier Transform of the heart rate signal from above. This graph is called a spectrogram.
Along the x-axis is the frequency of the sine waves that are put together. Frequency is how many repeats per second (measured in Hertz). Since we are looking at heart rate and that never goes much above 1 or 2 repeats per second, the x-axis only goes up to 1Hz.
The y-axis shows the power for each sine wave: the more power, the more of that sine wave we need.

Looking at this graph, we can see that there is a bit of a raised bump around 0.04-0.1Hz. This tells us that there is more variation in the heart rate signal that has that time-scale than at the time-scales either side. This is quite a long time-scale – it is changes in heart rate that are taking place over 10-25 secs.

So we can see what we have done, we can show what the heart rate signal looks like when we just consider those time-scales and no others. That is show in this next figure. The red line is what the heart rate signal looks like after we filter out all the very short time-scale changes and focus on changes that are sustained for 10 seconds or more.

Please note that this last graph is really a deviation from where we are headed, but it gives some sense of what we are doing.

Which time-scales are we interested in?

There has been a lot of research into this. The broad consensus is that there are two bands of heart-rate frequency that are of interest.

1. The low frequency band which ranges from 0.04 to 0.15Hz. This band probably (almost certainly) corresponds to activity in the sympathetic nervous system (SNS) as manifest in vagus nerve activity and can be thought, loosely speaking, as a response to uncertainty and other stressors. The higher the power in this band of frequencies is, the more tense the person is.
2. A higher frequency band which ranges from 0.15 to 0.4Hz. This band is less easily interpreted, but has some of the characteristics one might expect if it were related to the parasympathetic nervous system (PNS).

The commonly used mnemonic for this is:
SNS: fight or flight
PNS: rest and digest
Although we think of these as being opposites, they are not so clearly opposites in reality.

The Next (and Last) Step

There is one more step in the analysis that we need. Look at the low frequency heart rate graph again, focussing on two segments of time in it.

1. Orange: the time between 25 and 35secs. Here the heart rate isn’t varying much, either in the original signal or in the low frequency signal.
2. Green: the time between 75 and 85secs. Here, despite the heart rate being lower than in the orange region, the heart rate is changing much more. Something more interesting is happening here.

That difference between the green and orange regions is the information we want – that says that the heart rate variability (in the low frequency) is higher in the green region than the orange region. We have done this informally by eye and now we need to do it by a calculation. This calculation is the final step.

Let’s start with just those two regions and the heart-rate signal within each. It’s really just as easy as this. Think of just taking all the heart rate values in each region and calculating the variance (or standard deviation) of those values. The orange region would give us a small variance and the green region would give us a much larger variance. We can make a graph of this, with just those two points on:

and the result confirms our observation that there is more variability in the 10 secs before time=85secs than in the 10 seconds before time=35secs. Note that the points we have drawn here give the heart rate variability over the preceding 10 seconds – so we have drawn them at time=35secs and 85secs.

We can do this calculation for every point along the time axis – find the variance of the heart rate over the preceding 10 seconds. That gives us a continuous line for the graph:

This is the graph we want. It shows us how heart rate variability in the low frequency band changes over time. The value at each point along the time axis is how much variability in heart rate there has been over the preceding 10 seconds.

Footnote: the variance calculation can be done in the longhand way, or as a shortcut it is also given by the area under the curve in the spectrogram. This speeds the calculation up dramatically without changing the result.

R Code

The code for the analysis has 3 main steps:

1. Read the data from files that are created during data acquisition.
2. Process the data to produce the parameters/quantities of interest.
3. Display the results

For plots, we will use the library ggplot2. To load this you must start a session with the command:
library(ggplot2)
and then this command:
source("https://github.com/rjwatt42/hrv/raw/main/hr_start.R")

Some of this code is highly specific for the application we have – reading files that are created by the Biopac and by the eprime experiment protocol. We will deal with that first.

Read the data

All of the data that we require is saved in two text files for each participant. A text file is one where everything is stored as a printable character and each character has a unique code. So, for example, the number 3.14 is stored as a sequence of codes: 51 46 49 52.

The first text file is the one produced by the Biopac. The code to read this is:
hdata<-hr_load_acq_txt("filename1")
where “filename1” refers to which file is to be read. This uses a function named hr_load_acq_txt which takes the filename and returns an object hdata that holds the contents of the file.

For our purposes, an object is like a container that has slots for various different things. The slots are called fields. In this case, the object hdata has fields for 2 things:

1. an object containing the heart rate recording data hdata$hr as
   1. times: hdata$hr$hr_times
   2.  heart rates: hdata$hr$hr_signal
2. the timings of events that were sent to the Biopac hdata$stimuli

If we want to see the variation in heart rate with time, then we can just plot:
ggplot(hdata$hr)+geom_line(aes(x=hr_time,y=hr_signal))

The second text file is the one produced by the eprime system. This contains a lot of information about the session that was run. The part of it that is most relevant for us is the part that lists the sequence of stimuli (the pieces of music). We already know exactly when they occurred, but not what order they came in. So we know that a stimulus started at 10.3 seconds, but not which one it was. That extra information is stored in this eprime file.

The code to read this is:
stimulusdata<-hr_load_stimulus_txt("filename2")
and it also returns an object, which in this case has just one field (slot) for the sequence of stimulus names.

The final step in reading the data is to combine the information from the two functions, which we can do with this function:
hrobject<-hr_eventProcess(hdata,stimulusdata)
The object that this returns holds everything we need and will be what we use for here onwards. This function also holds information about how long after the start of each stimulus, the unexpected event(s) occurred. It uses this information to calculate where in the heart rate timeline those events occur.

The hrobject has these fields (slots):
hrobject$hr contains the same fields and data as are in hdata$hr, so the heart rate timeline and heart rate values
hrobject$stimuli contains a list of the start times of each stimulus
hrobject$eventTimes contains a list of unexpected event times
hrobject$eventList contains the corresponding list of unexpected event names

Note: the details of how these two functions work are probably not worth spending time on. They are, in the best sense of the word, hacks. They work because of very specific details of how those files are set up. For example, the actual time & hr data in the Biopac file are all on lines that use tab characters (with code 9) to separate the columns. No other lines in the file have that character.

Displaying the Heart Rate & Events

We showed above how the heart rate could be plotted. There is a convenience function that does more and make the graph more useful:
hr_plot_hr(hrobject)
This produces a plot that is suitable for further use (it was used for several of the graphs above). At its heart, this function uses the same code that we saw above:
ggplot(hdata$hr)+geom_line(aes(x=hr_time,y=hr_signal))
However, it has a few extra lines. It uses a theme that makes the graph look clearer; it makes the graph axes more suitable and labels them; and it has the option to also draw the stimulus starts and the unexpected events.

A look at the HRV Data

The first graph here shows the response in the low frequency HRV to the onset of a stimulus (a piece of music).

The stimulus starts at time zero. The heavy line is the average response, across participants and stimuli. The paler lines are standard errors of measurement (so they are equivalent to error bars).

There are three things to notice

1. The HRV before the onset of the stimulus is flat – that is useful.
2. The HRV for the first 3 seconds or so of the stimulus reduces. This corresponds to the participants becoming more relaxed. This is what we would have expected/hoped for.
3. And then the HRV stops reducing and after about 7 or 8 seconds rises again. This might be due to two things: either participants begin to expect a change; or one of the unexpected events occurs. A further analysis might help with this.

The second graph shows the response to the unexpected event:

1. This time the HRV before the unexpected event is not flat: but instead shows the falling pattern we would expect, given the graph above.
2. At time zero (when the unexpected event occurs), that fall in HRV comes to an end.
3. By 3 seconds after the unexpected event, HRV is rising and continues to do so for at least a further 5 seconds.

This is what we would expect if the unexpected event introduces uncertainty and tension.