In FlukeAndFeather/cetaceanbcg: Cetacean Heart Beats Detected by Ballistocardiogram
Keywords: r rmarkdown::metadata$keywords
knitr::opts_chunk$set(
collapse = TRUE,
warning = FALSE,
message = FALSE,
echo = FALSE,
cache = TRUE,
comment = "#>",
fig.path = "../figures/",
fig.height = 6,
fig.width = 6,
dev = "pdf",
dpi = 300
)
library(cetaceanbcg)
library(tidyverse)
Introduction
Recent advances in physio-logging (recording physiological variables using animal-borne devices) have largely been driven by new developments in sensor technology [@hawkesIntroductionThemeIssue2021; @williams2021]. For example, new physio-logging tags can detect regional changes in blood flow by incorporating functional near-infrared spectroscopy sensors [@mcknightShiningNewLight2021]. However, traditional inertial measurement unit (IMU) tags equipped with accelerometers and other inertial sensors can also measure important physiological and related variables, such as wing beat frequency [@patterson2019] and feeding rate [@divirgilio2018]. Through careful inspection and analysis of high-resolution acceleration, scientists have measured elevated respiration rates following record-breaking dives [@satoStrokeRatesDiving2011], near-continuous feeding by small cetaceans [@wisniewskaUltraHighForagingRates2016], social interactions between large cetaceans [@goldbogenUsingAccelerometersDetermine2014], and important biomechanical variables including movement speed [@cadeDeterminingForwardSpeed2018]. While physio-logging tags with cutting-edge biomedical technologies push the boundaries of physiological field research, simpler IMU tags have fewer logistical constraints and provide access to more species and larger sample sizes. This is particularly important for species that cannot be restrained or studied in managed care. For example, of the sixteen species of baleen whales (Mysticeti), heart rate has only been recorded with an electrocardiogram tag in the wild for one blue whale (Balaenoptera musculus) [@goldbogenExtremeBradycardiaTachycardia2019; but see @ponganisHeartRateElectrocardiogram1999]. Conversely, IMU tags have been deployed on hundreds of individuals of nearly every species in the clade for the last twenty years [@nowacekBuoyantBalaenidsUps2001]. These existing datasets (and future IMU tag deployments) could hold additional valuable physiological information, awaiting proper computational methods for mining them.
The ballistocardiogram (BCG) has potential applications to using accelerometers as heartrate monitors in both the wild and in managed care [@giovangrandi2011ballistocardiography; @sadekBallistocardiogramSignalProcessing2019; @inanBallistocardiographySeismocardiographyReview2015]. Ballistocardiography is a noninvasive method for measuring cardiac function based on the ballistic forces involved in the heart ejecting blood into the major vessels. The BCG originated as a clinical tool in the first half of the 20th century [@starrStudiesEstimationCardiac1939], but was largely superseded by electro- and echocardiography. However, potential novel applications like passive monitoring of heart function in at-risk populations [@giovangrandi2011ballistocardiography] has led to a recent resurgence of ballistocardiography research, with advances in hardware [@andreozzi2021] and signal processing methodology [@sadekBallistocardiogramSignalProcessing2019]. While the BCG is a three-dimensional phenomenon, it is strongest in the cranio-caudal axis [@inanBallistocardiographySeismocardiographyReview2015]. Along this axis, the waveform is composed of multiple peaks and valleys; most prominent of these is the so-called IJK complex [@pinheiroTheoryDevelopmentsUnobtrusive2010]. The precise physiological mechanism underlying the BCG waveform has not been fully resolved [@kim2016], but it has been established that the IJK complex occurs during systole and, in humans, occurs at approximately the same time as the T-wave in an electrocardiogram (ECG) [@inanBallistocardiographySeismocardiographyReview2015]. The BCG J wave is the most robust feature in the waveform and is typically used for detecting heart beats [@inanBallistocardiographySeismocardiographyReview2015].
Here we present a method for generating a BCG from bio-logger cranio-caudal acceleration. We validated our method with a simultaneously recorded ECG on an adult killer whale in managed care (Orcinus orca) and applied it to detect heartrate in a blue whale. The relative orientation of a tag on a cetacean's body is often uncertain when bio-loggers are deployed in the wild [@johnsonDigitalAcousticRecording2003], so isolating acceleration along the cranio-caudal axis is subject to error. Therefore, we also compared a tri-axial BCG to the cranio-caudal BCG. Specifically, we tested three hypotheses to validate our method. First, a cranio-caudal (1D) BCG would, in a controlled setting, produce instantaneous heartrates that are statistically equivalent to ECG instantaneous heartrates. Second, a tri-axial (3D) BCG would, in a field setting, produce a more robust signal than a 1D BCG. Third, BCG-derived heartrates would increase during the later phases of dives, consistent with the progressive increase in heartrate routinely observed prior to and during ascent [@goldbogenExtremeBradycardiaTachycardia2019; @mcdonaldDeepdivingSeaLions2014].
Materials and methods
Animal tagging
Killer whale
A 3868 kg adult female killer whale in managed care at SeaWorld of California, San Diego, CA was double-tagged with an archival Customized Animal Tracking Solutions IMU (CATS, www.cats.is) tag and a custom-built, archival ECG tag on August 16, 2021 as part of clinical animal cardiac evaluations under the SeaWorld USDA APHIS display permit. The ECG tag hardware and data processing procedures were previously described by @bickettHeartRatesHeart2019. Both tags were deployed by hand and attached with suction cups. We attached the CATS tag on the mid-lateral left chest posterior to the pectoral fin (Movie S1). The CATS tag recorded tri-axial acceleration at 400 Hz, tri-axial magnetometer and tri-axial gyroscope at 50 Hz, pressure at 10 Hz, and video at 30 fps. The IMU in the CATS tag was a MPU-9250 (InvenSense, San Jose, CA; www.invensense.com). The accelerometer had dynamic range of ±4 g, sensitivity of 8,192 LSB g^-1^, and accuracy of 6.1 × 10^-6^ g. All sensors were rotated from the tag's frame of reference to that of the whale using MATLAB (MathWorks, Inc., v2020b) tools for processing CATS data [@cadeToolsIntegratingInertial2021]. This rotation aligned the tag's x-, y-, and z- axes with the cranio-caudal, lateral, and dorso-ventral axes of the whale, respectively. We attached the ECG tag approximately midline on the ventral chest just caudal (posterior) to the axilla and we recorded the ECG at 100 Hz. Individual heart beats in the ECG record were identified from visually verified R-waves using a customized peak detection program (K. Ponganis; Origin 2017, OriginLab Co., Northampton, MA). ECG and IMU were recorded during a spontaneous breath hold while the whale rested at the surface.
Blue whale
A 24.5 m blue whale was tagged with an archival, suction-cup CATS IMU tag on September 5, 2018 in Monterey Bay, CA under permits MBNMS-MULTI-2017-007, NMFS 21678, and Stanford University IACUC 30123 [previously published by @goughScalingSwimmingPerformance2019]. We deployed the tag using a 4 m fiberglass pole from a 6.3 m rigid-hulled inflatable boat and recovered it via radio VHF tracking [as described by @goldbogen2006]. The tag slid behind the left pectoral flipper, similar to the placement of the CATS tag on the killer whale. Tag configuration and data processing followed the same procedure as the killer whale, including accelerometer specification and sampling rates for inertial sensors and video. The 400 Hz acceleration data was used for ballistocardiography (see section Signal processing). We downsampled the multi-sensor data to 10 Hz for movement analysis using the MATLAB CATS tools [@cadeToolsIntegratingInertial2021].
Signal processing
The BCG waveform is three dimensional, but strongest in the cranio-caudal axis [@inanBallistocardiographySeismocardiographyReview2015]. We tested both 1-dimensional (cranio-caudal only) and 3-dimensional metrics for identifying heartbeats in acceleration data based on the methods of @leePhysiologicalSignalMonitoring2016. For windowed operations (such as moving averages and signal filters), we used 0.5 s windows for killer whale data and 2.0 s windows for blue whale data, corresponding to 200 and 800 data points, respectively. The different window sizes were determined through trial and error to remove noise while retaining signal shape. Generally, longer windows will be necessary for larger animals due to their slower heartrates [@Stahl1967b].
Procedure
- First we removed noise and de-trended the acceleration signal with a 5th order Butterworth band-pass filter (killer whale: [1-25 Hz], blue whale: [1-10 Hz]) (R package
signal) [@R-signal]. The lower cut-off frequency de-trended the data. 1 Hz should be appropriate for most marine mammal species. The upper cut-off frequency removed noise. A suitable upper cut-off frequency will depend on body size; larger species' bodies produce lower magnitude accelerations [@martínlópez2021], so more conservative upper cut-off frequencies may be applied to remove more noise without sacrificing signal shape clarity.
- Then we enhanced the IJK complex by differentiating acceleration using a 4th order Savitzky-Golay filter (R package
signal). Differentiation (i.e., $a_{t+1}-a_t$, where $a_t$ is the observed acceleration at time step $t$) exaggerates peaks, like the J wave, but it is sensitive to noisy signals. Therefore, additional noise reduction is necessary prior to differentiation. A moving average smoother could remove noise, but it would also reduce the amplitude of peaks. Hence, differentiating Savitzky-Golay filters are preferred in peak-detection algorithms because they remove noise while retaining the general shape of peaks [@samann2019]. We described the resulting signal as "differenced acceleration", rather than jerk, because we did not take the derivative of acceleration with respect to time. The purpose of this signal was to exaggerate a phenomenon in the signal (i.e., the J wave), not to describe a physical quantity (i.e., jerk).
- We further enhanced the peaks in the differenced acceleration signal by calculating the Shannon entropy ($H_i=-\sum_{k} |a_{ik}| \times ln(|a_{ik}|)$, where $k$ is the acceleration axis). Additionally, Shannon entropy is strictly positive, which facilitated peak detection. In the 1D BCG, $k$ was surge (cranio-caudal acceleration). In the 3D BCG, $k$ included surge, sway (lateral acceleration), and heave (dorso-ventral acceleration).
- After enhancing the peaks through differentiation and entropy calculation, we removed residual noise by applying a triangular moving average (TMA) smoother. TMAs are equivalent to applying a simple moving average in two passes, which applies greater weight to the middle part of the window and retains peaks and valleys more clearly. After steps 2 and 3, the signal was clear enough that TMAs provided satisfactory results, obviating the need for a more complex algorithm like a Savitzky-Golay filter at this stage. We described the resulting signal as the BCG.
- The BCG contained major peaks (corresponding to heartbeats) and minor peaks (residual noise) (Fig S1A). We extracted all peaks from the BCG and applied a clustering algorithm to retain major peaks and reject minor peaks. First, we extracted all peaks in the BCG signal using
findpeaks() (R package pracma) [@R-pracma] with a minimum peak distance equivalent to the window size (0.5 s for the killer whale, 2.0 s for the blue whale). For each peak, we calculated its absolute height and its prominence (i.e., height relative to the lowest valley between a peak and its higher neighbors). Then, we calculated each peak's Euclidean distance in height-prominence space from the highest peak (Fig S1B) and estimated the density distribution of these distances (Fig S1C). The density distribution was bimodal, with a low-distance peak corresponding to major peaks and a high-distance peak corresponding to minor peaks. We used the distance corresponding to the valley between the two peaks as a threshold for rejecting minor peaks (Fig S1D).
This procedure may be applied to either 1-dimensional (i.e., cranio-caudal only) or 3-dimensional acceleration. In the case of 3-dimensional acceleration, the band-pass and Savitzky-Golay filters were applied to each axis independently.
BCG validation with killer whale ECG
We fit ordinary least squares regression to BCG-derived instantaneous heartrates with respect to ECG-derived heartrate and tested 1) if the intercept was significantly different than 0 and 2) if the slope was significantly different than 1. We calculated the mean and standard deviation of absolute error as an equivalence measure (1-dimensional BCG only).
BCG application to blue whale
Dynamic body movements produce an acceleration signal that masks the ballistocardiogram, so we limited our analyses to motionless periods (Fig. \@ref(fig:bw-bcg)B-C). These periods occurred during or near the bottom phase of dives between fluke strokes. Strokes were detected from visual examination of the rotational velocity around the lateral axis recorded by gyroscope [sensu @goughScalingSwimmingPerformance2019]. We used gyroscopes for stroke detection because they were 1) separate sensors from the accelerometers and 2) strokes are clearly visible in gyroscope signals and are robust to tag placement.
We tested whether the 3-dimensional BCG was more robust than 1-dimensional BCG in field data by comparing the signal-to-noise ratios. For both BCGs, we calculated the power spectral density (R package psd) [@Barbour2014]. Previously recorded blue whale apneic heart rate was 4-8 beats per minute (bpm) [@goldbogenExtremeBradycardiaTachycardia2019], so we quantified signal as the integration of the power spectral density curve from 4-8 bpm and noise as the integrated remainder, up to 60 bpm. The sample size recorded by @goldbogenExtremeBradycardiaTachycardia2019 was one individual, so we could not account for potential inter-individual variation. Nonetheless, 4-8 bpm was the best available estimate for typical blue whale apneic heart rate.
We also tested whether BCG-derived instantaneous heart rates were consistent with the range and pattern of heart rates previously observed in the blue whale and other marine mammals; namely a gradual increase in heart rate later in the dive, especially during the final ascent [@goldbogenExtremeBradycardiaTachycardia2019; @mcdonaldDeepdivingSeaLions2014]. We assigned dive start and end times when the whale swam deeper than 2 m, retaining dives that exceeded 10 m depth and 5 minutes duration. Dive times were normalized from 0 (start of dive) to 1 (end of dive). We regressed instantaneous heart rate against normalized dive time using robust Theil-Sen regression (to account for heteroscedascity) (R package RobustLinearReg) [@R-RobustLinearReg; @Sen-1968; @theilRankInvariantMethodLinear1992] and tested whether the slope was greater than 0.
Reproducibility
The data and code used in this analysis were packaged as a research compendium, containing the data, code, and an executable version of this manuscript. We used the R package rrtools [@marwickPackagingDataAnalytical2018; @rrtools2019] to initialize the compendium, which was written as an R package. This approach promotes reproducibility and facilitates adoption by other researchers [@alstonBeginnerGuideConducting2021; @powersOpenScienceReproducibility2019; @stoddenEmpiricalAnalysisJournal2018]. The steps described in section Signal processing were implemented as functions in the R package, and the executable manuscript demonstrates how to use those functions to perform the analyses presented in this study.
Results and discussion
BCG validation with killer whale ECG
ca_tz <- "Etc/GMT+7"
bcg_begin <- lubridate::ymd_hms("2021-08-16 10:09:24.5", tz = ca_tz)
bcg_end <- lubridate::ymd_hms("2021-08-16 10:09:38", tz = ca_tz)
oo_400hz <- corky_400hz %>%
filter(between(dt, bcg_begin - 10, bcg_end + 10)) %>%
mutate(surge_filt = filter_acc(surge, 400, 25.0),
surge_diff = lead(surge_filt) - surge_filt,
surge_se = shannon_entropy(surge_diff),
surge_smooth = tma(surge_se, 0.5 * 400)) %>%
filter(between(dt, bcg_begin, bcg_end)) %>%
mutate(bcg_beat = find_beats(surge_smooth, 400, 0.5),
bcg_bpm = bpm(bcg_beat, dt))
# Align ECG-detected heart beats
oo_400hz$ecg_beat <- FALSE
beat_idx <- approx(oo_400hz$dt,
seq_along(oo_400hz$dt),
xout = ecg_beats$dt,
method = "constant")$y %>%
na.omit()
oo_400hz$ecg_beat[beat_idx] = TRUE
oo_400hz$ecg_bpm <- with(oo_400hz, bpm(ecg_beat, dt))
bcg_bpm <- filter(oo_400hz, bcg_beat) %>%
select(dt, bcg_bpm)
nearest_bcg <- function(t) {
bcg_bpm$bcg_bpm[which.min(abs(t - bcg_bpm$dt))]
}
ecg_bpm <- filter(oo_400hz, ecg_beat) %>%
mutate(bcg_bpm = map_dbl(dt, nearest_bcg),
bpm_diff = abs(bcg_bpm - ecg_bpm) / ecg_bpm) %>%
drop_na(bpm_diff)
mean_bpm_diff <- mean(ecg_bpm$bpm_diff)
sd_bpm_diff <- sd(ecg_bpm$bpm_diff)
bpm_lm <- lm(bcg_bpm ~ ecg_bpm, ecg_bpm)
bpm_stderr <- sqrt(diag(vcov(bpm_lm)))
The ECG and BCG yielded nearly identical heart rate estimations (Fig. \@ref(fig:oo-bcg-ecg)). We collected 14 s of simultaneous ECG and BCG data during a motionless breath hold at the surface. Logistical constraints prevented us from gathering a longer sample, as these data were collected secondary to other projects. BCG-derived instantaneous heart rates were within r sprintf("%0.1f%% \u00B1 %0.1f%%", mean_bpm_diff * 100, sd_bpm_diff * 100) of the ECG-derived rates (mean ± standard deviation). Ordinary least squares regression of BCG heartrates on ECG heartrates yielded a slope of r sprintf("%0.2f \u00B1 %0.2f", coef(bpm_lm)[2], bpm_stderr[2]) and intercept of r sprintf("%0.2f \u00B1 %0.2f", coef(bpm_lm)[1], bpm_stderr[1]) (mean ± standard error), which were not significantly different from the hypothesized 1 and 0, respectively. (Fig. \@ref(fig:validation-plots)C).
BCG application to blue whale
bw_elg <- bw180905_53_elg %>%
mutate(
depth_min = map2_dbl(start, stop, function(t1, t2) {
bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>% min()
}),
depth_max = map2_dbl(start, stop, function(t1, t2) {
bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>% max()
})
)
n_motionless <- nrow(bw_elg)
motionless_minutes <- as.numeric(sum(bw_elg$stop - bw_elg$start), unit = "mins")
We generated 1-dimensional and 3-dimensional BCGs for 2 hours of data, including 10 rest dives and r n_motionless motionless periods totaling r format(motionless_minutes, digits = 3) minutes (r sprintf("%0.1f%%", motionless_minutes / 120 * 100) of the 2-hour record) (Fig. \@ref(fig:bw-bcg)A-C).
bw_400hz <- bw180905_53_400hz %>%
mutate(
across(surge:heave,
filter_acc, fs = 400, upper = 10.0,
.names = "{.col}_filt"),
surge_diff = lead(surge_filt) - surge_filt,
surge_se = shannon_entropy(surge_diff),
surge_smooth = tma(surge_se, 0.5 * 400),
jerk = jerk(cbind(surge_filt, sway_filt, heave_filt),
fs = 400, p = 4, n = 2 * 400 + 1),
jerk_se = shannon_entropy(jerk),
jerk_smooth = tma(jerk_se, 2 * 400),
# Annotate regions
rid_left = approx(bw_elg$start,
bw_elg$region_id,
dt,
"constant")$y,
rid_right = approx(bw_elg$stop,
bw_elg$region_id,
dt,
"constant",
yleft = 0)$y + 1,
region_id = ifelse(rid_left == rid_right, rid_left, NA),
# Zero-out signal in non-valid regions (i.e. remove movement artifacts)
jerk_smooth = ifelse(is.na(region_id), 0, jerk_smooth),
bcg_beat = find_beats(jerk_smooth, 400, 2)
) %>%
# Calculate heart rate within each region
group_by(region_id) %>%
mutate(bcg_bpm = bpm(bcg_beat, dt)) %>%
ungroup() %>%
select(-c(rid_left, rid_right))
walk(unique(na.omit(bw_400hz$region_id)), function(rid) {
p <- bw_bcg_plot(rid)
ggsave(here::here(sprintf("analysis/figures/bw_bcg/bw_bcg_%02i.png", rid)), p)
})
The 3-dimensional BCG (Fig. \@ref(fig:bw-bcg)) produced a more robust signal (i.e., higher signal-to-noise ratio) than the 1-dimensional BCG, which used only cranio-caudal acceleration (Fig. \@ref(fig:validation-plots)A). The signal-to-noise ratio was 2.00 for the 3-dimensional BCG, compared to 0.17 for the 1-dimensional BCG. Although the power spectral density curve for the 1-dimensional BCG had a peak in the 4-8 bpm frequency range, most of the signal's power was concentrated in lower frequencies. Conversely, the 3-dimensional BCG's power was concentrated precisely in the 4-8 bpm frequency range, with only a smaller peak in the lower frequencies.
bcg1d_psd <- with(bw_400hz, pspectrum(surge_smooth[!is.na(region_id)], 400))
bcg3d_psd <- with(bw_400hz, pspectrum(jerk_smooth[!is.na(region_id)], 400))
bw_psd <- tibble(
metric = "Surge only",
freq_hz = bcg1d_psd$freq,
spec = bcg1d_psd$spec
) %>%
rbind(
tibble(
metric = "Triaxial",
freq_hz = bcg3d_psd$freq,
spec = bcg3d_psd$spec
)
) %>%
mutate(freq_bpm = freq_hz * 60) %>%
filter(freq_bpm <= 60)
psd_s2n <- bw_psd %>%
group_by(metric) %>%
summarize(signal = pracma::trapz(freq_bpm[freq_bpm >= 4 & freq_bpm <= 8],
spec[freq_bpm >= 4 & freq_bpm <= 8]),
noise = pracma::trapz(freq_bpm, spec) - signal) %>%
mutate(s2n = signal / noise)
bw10hz <- bw180905_53_10hz %>%
mutate(
dive_id = split_dives(dt, depth,
surface = 2, min_depth = 10, min_dur = 5 * 60),
dive_norm = normalize_dives(dive_id),
rid_left = approx(bw_elg$start,
bw_elg$region_id,
dt,
"constant")$y,
rid_right = approx(bw_elg$stop,
bw_elg$region_id,
dt,
"constant",
yleft = 0)$y + 1
)
bw_400hz$dive_norm <- approx(bw10hz$dt, bw10hz$dive_norm, bw_400hz$dt)$y
bpm_data <- bw_400hz %>%
drop_na(bcg_bpm)
bpm_ts <- RobustLinearReg::theil_sen_regression(bcg_bpm ~ dive_norm, bpm_data)
bpm_pred <- tibble(dive_norm = seq(0, 1, length.out = 100)) %>%
mutate(bcg_bpm = predict(bpm_ts, newdata = .))
The 3-dimensional BCG exhibited increasing heart rates over the course of dives. Average heart rate increased from r format(bpm_pred$bcg_bpm[1], digits = 2) bpm at the start of dives to r format(bpm_pred$bcg_bpm[nrow(bpm_pred)], digits = 2) bpm at the end of dives (Theil-Sen regression, $p < 10^{-10}$) (Fig. \@ref(fig:validation-plots)B).
Limitations and considerations for future applications
While the BCG method presented here holds the potential to mine existing and future marine mammal bio-logging datasets for information about cardiovascular function, it has several limitations compared to ECG methods. Most importantly, BCGs are highly sensitive to movement artifacts [@inanBallistocardiographySeismocardiographyRemovview2015], so only motionless periods are valid for analysis. This limits the behavioral and physiological contexts in which heartrate may be measured. For example, the BCG is probably an inappropriate method for quantifying the magnitude of surface tachycardia [@goldbogenExtremeBradycardiaTachycardia2019] and exercise modulation of bradycardia [@noren2012], due to movement artifacts during those activities. Secondarily, we did not test whether the BCG is robust to tag placement location. The blue whale data presented in this study was collected when a dorsally-deployed tag slipped to the lateral chest cavity behind a flipper, where it is reasonable to expect greater accelerations caused by heart beats than from a tag farther from the animal's center of mass. It is possible that the ballistic forces generated by heart beats are strong enough to produce an interpretable BCG for a variety of potential tag deployment locations, but this likely varies with animal body size, as well as accelerometer sampling rate and sensitivity.
When auditing existing bio-logging data and planning future tag deployments for BCG analysis, careful consideration should be paid to sampling rate. As a rule of thumb in signal processing, the sampling rate should be at least twice the frequency of the phenomenon of interest. In the case of the BCG, the relevant frequency is that of the BCG waveform, not the heartrate. In humans, the power of the IJK-complex (the part of the BCG waveform used for heart beat detection) is concentrated between 4-7 Hz [@moukadem2018]. It is unlikely that marine mammal BCG waveforms have a higher frequency than humans, owing to their generally larger body sizes. Therefore, it is possible that BCGs may be generated for accelerometer sampling rates as low as 10-15 Hz. Conservatively, the authors recommend a sampling rate of no less than 50 Hz (i.e., twice the upper cut-off frequency of the widest bandpass filter used in this study).
Future bio-logging BCG methodology research should address the limitations imposed by tag-placement and movement artifacts. We used accelerometers in this study because of their prevalence in bio-logging research, but it is possible that other widely used bio-logger sensors, such as gyroscopes and/or magnetometers, could produce a clearer signal in a greater variety of contexts. Alternative bio-logger housing designs, such as limpet-style tags [@andrews2008] or "marine skin" [@nassar2018], could reduce noise, boost the signal-to-noise ratio, and make the method more widely applicable.
Conclusions
Here we presented a ballistocardiogram method for detecting resting apneic heartrate in cetaceans using accelerometers. We validated the method in a controlled setting with simultaneous ECG and in a field setting by confirming expected physiological patterns. As accelerometer tags have been deployed on many cetacean species for multiple decades, this method may be applied to mine existing datasets and better understand how heartrate scales with body size and other biological factors. It may also provide additional data for conservation physiology applications. For example, BCGs extracted from gliding phases before and after controlled sonar exposure experiments could quantify the physiological response to anthropogenic disturbance [@southall2019]. Even as the field of physio-logging progresses with new hardware innovations [@williams2021; @fahlman2021; @hawkesIntroductionThemeIssue2021], this method demonstrates that computational advances can derive new insights from traditional sensors.
Acknowledgements
The authors are grateful to the SeaWorld of California Killer Whale training staff for their efforts and support. We also thank Anna Krystalli, Ben Marwick, Karthik Ram, Nicholas Tierney, and other members of the R community for developing tools and educational resources to facilitate open science practices. This is a SeaWorld Parks and Entertainment Technical Contribution number 2021-12. We thank Lucy Hawkes and two anonymous reviewers for comments on the manuscript.
Footnotes
Author contributions
Conceptualization: M.F.C.,J.A.F.,P.J.P.,J.A.G.; Methodology: M.F.C.,J.A.F.,P.J.P.; Software: M.F.C.; Formal analysis: M.F.C.,P.J.P.; Investigation: M.F.C.,J.A.F.,P.J.P.; Resources: P.J.P.,J.A.G.,T.L.S.; Writing - original draft: M.F.C; Writing - review & editing: M.F.C.,J.A.F.,P.J.P.,J.A.G.,T.L.S.; Supervision: P.J.P.,J.A.G.; Project administration: P.J.P.,J.A.G.; Funding acquisition: J.A.G.
Funding
This work was supported by grant N000141912455 from the Office of Naval Research. M.F.C. was supported by the Stanford University William R. and Sara Hart Kimball Fellowship and a Stanford Data Science Scholar Fellowship.
Data availability
All data and code used in this analysis are available on Zenodo [@cetaceanbcgzenodo].
Competing interests
The authors declare no competing interests.
\newpage
Figures
t <- theme_classic(base_size = 8) +
theme(axis.title.x = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_blank())
t2 <- t + theme(plot.margin = margin(l = 0.5, unit = "cm"))
t_yvoid <- t2 +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
strip.text = element_blank())
# A: ECG
ecg_sample <- corky_ecg %>%
filter(between(dt, first(oo_400hz$dt), last(oo_400hz$dt)))
pA <- ggplot(ecg_sample, aes(dt, ecg_uV)) +
geom_line(size = 0.2) +
geom_point(data = filter(ecg_sample, ecg_beat), color = "red") +
scale_x_datetime(date_labels = "%H:%M:%S") +
labs(y = bquote(ECG~(mu * V))) +
expand_limits(x = max(oo_400hz$dt) + 1) +
t2
# B: Filtered surge
ijk_range <- as.POSIXct("2021-08-16 10:09:31.4", tz = "US/Pacific") + c(0, 0.4)
ijk_complex <- oo_400hz %>%
filter(between(dt, ijk_range[1], ijk_range[2])) %>%
mutate(t = as.numeric(dt - min(dt), unit = "secs"))
ijk_labels <- slice(ijk_complex, c(41, 53, 67, 85, 110)) %>%
mutate(wave = c("H", "I", "J", "K", "L"))
pB_inset <- ggplot(ijk_complex, aes(dt, surge_filt)) +
geom_line() +
geom_text(aes(label = wave),
ijk_labels,
vjust = c(-0.2, 1.2, 0.75, -0.6, -0.2),
hjust = c(0.5, 0.5, -0.6, 0.5, 0.5),
fontface = "bold",
size = 3) +
theme_void()
inset_range <- as.POSIXct("2021-08-16 10:09:29.9", tz = "US/Pacific") + c(0, 4)
pB <- ggplot(oo_400hz, aes(dt, surge_filt)) +
geom_line(size = 0.2) +
annotation_custom(grob = ggplotGrob(pB_inset),
xmin = inset_range[1], xmax = inset_range[2],
ymin = 0.25, ymax = 0.5) +
annotate("rect",
xmin = ijk_range[1],
xmax = ijk_range[2],
ymin = min(ijk_complex$surge_filt) - 0.02,
ymax = max(ijk_complex$surge_filt) + 0.02,
fill = NA,
color = "grey50",
lty = 3) +
scale_x_datetime(date_labels = "%H:%M:%S") +
expand_limits(y = 0.5) +
labs(y = "Filtered\nacceleration") +
t_yvoid
# C: Surge diff
pC <- ggplot(oo_400hz, aes(dt, surge_diff)) +
geom_line(size = 0.2) +
scale_x_datetime(date_labels = "%H:%M:%S") +
labs(y = "Differenced\nacceleration") +
t_yvoid
# D: Shannon entropy
pD <- ggplot(oo_400hz, aes(dt, surge_se)) +
geom_line(size = 0.2) +
scale_x_datetime(date_labels = "%H:%M:%S") +
labs(y = "Shannon\nentropy") +
t_yvoid
# E: Smoothed
pE <- ggplot(oo_400hz, aes(dt, surge_smooth)) +
geom_line(size = 0.2) +
geom_point(data = filter(oo_400hz, bcg_beat), color = "blue") +
scale_x_datetime(date_labels = "%H:%M:%S") +
labs(y = "BCG") +
t_yvoid
cowplot::plot_grid(pA, pB, pC, pD, pE,
align = "v",
axis = "lr",
ncol = 1,
rel_heights = c(1, 1.5, 0.75, 0.75, 1),
labels = "AUTO")
# Plot the BCG for a motionless period in blue whale data
bw_bcg_plot <- function(rid) {
bcg_sample <- filter(bw_400hz, region_id == rid)
palette <- rgb(c(230, 86, 0), c(159, 180, 158), c(0, 233, 115), maxColorValue = 255)
t_lmargin <- t_yvoid +
theme(plot.margin = margin(l = 0.5, unit = "cm"))
# A: Filtered acceleration
acceleration_data <- bcg_sample %>%
rename(`Cranio-caudal` = surge_filt,
Lateral = sway_filt,
`Dorso-ventral` = heave_filt) %>%
pivot_longer(`Cranio-caudal`:`Dorso-ventral`,
names_to = "axis",
values_to = "acc") %>%
mutate(axis = factor(axis, levels = c("Cranio-caudal",
"Lateral",
"Dorso-ventral")))
pA <- ggplot(acceleration_data, aes(dt, acc, color = axis)) +
geom_line(size = 0.2) +
scale_x_datetime(date_labels = "%H:%M:%S") +
scale_y_continuous(n.breaks = 3) +
scale_color_manual(values = palette) +
facet_grid(rows = vars(axis)) +
labs(y = "Filtered\nacceleration") +
t_lmargin +
theme(legend.position = "none",
panel.spacing = unit(0, "lines"))
# B: Jerk vectors
jerk_data <- bcg_sample %>%
mutate(jerk = set_names(as_tibble(jerk), c("Cranio-caudal",
"Lateral",
"Dorso-ventral"))) %>%
unpack(jerk) %>%
pivot_longer(`Cranio-caudal`:`Dorso-ventral`,
names_to = "axis",
values_to = "jerk") %>%
mutate(axis = factor(axis, levels = c("Cranio-caudal",
"Lateral",
"Dorso-ventral")))
pB <- ggplot(jerk_data, aes(dt, jerk, color = axis)) +
geom_line(size = 0.2) +
scale_x_datetime(date_labels = "%H:%M:%S") +
scale_y_continuous(n.breaks = 3) +
scale_color_manual(values = palette) +
facet_grid(rows = vars(axis)) +
labs(y = "Differenced\nacceleration") +
t_lmargin +
theme(legend.position = "none",
panel.spacing = unit(0, "lines"))
# C: Shannon entropy of jerk (NOT jerk norm)
pC <- ggplot(bcg_sample, aes(dt, jerk_se)) +
geom_line(size = 0.2) +
scale_x_datetime(date_labels = "%H:%M:%S") +
labs(y = "Shannon\nentropy") +
t_lmargin
# D: Smoothed
pD <- ggplot(bcg_sample, aes(dt, jerk_smooth)) +
geom_line(size = 0.2) +
geom_point(data = filter(bcg_sample, bcg_beat), color = "blue") +
scale_x_datetime(date_labels = "%H:%M:%S") +
labs(y = "BCG") +
t_lmargin
cowplot::plot_grid(pA, pB, pC, pD,
align = "v",
axis = "lr",
ncol = 1,
rel_heights = c(1.5, 1.5, 1, 1),
labels = c("D", "E", "F", "G"),
label_size = 10)
}
# Plot the depth profile and motionless periods for the blue whale
bw_10hz <- bw180905_53_10hz %>%
mutate(rid_left = approx(bw_elg$start,
bw_elg$region_id,
dt,
"constant")$y,
rid_right = approx(bw_elg$stop,
bw_elg$region_id,
dt,
"constant",
yleft = 0)$y + 1,
region_id = ifelse(rid_left == rid_right, rid_left, NA),
dive_id = split_dives(dt, depth,
surface = 2, min_depth = 10, min_dur = 5 * 60))
divesample <- filter(bw_10hz, dive_id == 5)
p_depth <- ggplot(bw_10hz, aes(dt, depth)) +
geom_line(size = 0.2) +
annotate("rect",
xmin = min(divesample$dt) - 60,
xmax = max(divesample$dt) + 60,
ymin = -1,
ymax = max(divesample$depth) + 1,
fill = NA,
color = "red",
linetype = "dotted") +
scale_x_datetime(date_labels = "%H:%M") +
scale_y_reverse("Depth (m)") +
theme_classic(base_size = 10) +
theme(axis.title.x = element_blank())
# Plot the depth and y-axis gyroscope for a motionless period
norm_y <- function(ymin, ymax, sign) {
(ymin + ymax) / 2 + sign * max(ymax - ymin) / 2
}
bw_elg <- bw180905_53_elg %>%
mutate(
depth_min = map2_dbl(start, stop, function(t1, t2) {
bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>% min()
}),
depth_max = map2_dbl(start, stop, function(t1, t2) {
bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>% max()
}),
depth_med = map2_dbl(start, stop, function(t1, t2) {
bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>%
first() %>%
median()
}),
ygyr_min = map2_dbl(start, stop, function(t1, t2) {
bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), ]$gw[,2] %>% min()
}),
ygyr_max = map2_dbl(start, stop, function(t1, t2) {
bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), ]$gw[,2] %>% max()
}),
ygyr_med = map2_dbl(start, stop, function(t1, t2) {
bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), ]$gw[,2] %>%
first() %>%
median()
}),
)
diveelg <- bw_elg %>%
filter(start >= first(divesample$dt), stop <= last(divesample$dt)) %>%
mutate(depth_min2 = norm_y(depth_min, depth_max, -1),
depth_max2 = norm_y(depth_min, depth_max, 1),
ygyr_min2 = norm_y(ygyr_min, ygyr_max, -1),
ygyr_max2 = norm_y(ygyr_min, ygyr_max, 1))
depth_profile <- ggplot(divesample, aes(dt, depth)) +
geom_line(size = 0.5) +
geom_rect(aes(xmin = start,
xmax = stop,
ymin = depth_min2,
ymax = depth_max2),
diveelg,
inherit.aes = FALSE,
fill = "orchid1",
color = NA,
alpha = 0.5) +
geom_rect(aes(xmin = start,
xmax = stop,
ymin = depth_min2,
ymax = depth_max2),
filter(diveelg, region_id == 23),
inherit.aes = FALSE,
fill = NA,
color = "red",
linetype = "dotted",
alpha = 0.5) +
scale_x_datetime(date_labels = "%H:%M") +
scale_y_reverse("Depth (m)") +
theme_classic(base_size = 10) +
theme(axis.title.x = element_blank(),
plot.margin = margin(l = 0.5, unit = "cm"))
ygyr_profile <- ggplot(divesample, aes(dt, gw[, 2])) +
geom_line(size = 0.2) +
geom_rect(aes(xmin = start,
xmax = stop,
ymin = ygyr_min2,
ymax = ygyr_max2),
diveelg,
inherit.aes = FALSE,
fill = "orchid1",
color = NA,
alpha = 0.5) +
geom_rect(aes(xmin = start,
xmax = stop,
ymin = ygyr_min2,
ymax = ygyr_max2),
filter(diveelg, region_id == 23),
inherit.aes = FALSE,
fill = NA,
color = "red",
linetype = "dotted",
alpha = 0.5) +
scale_x_datetime(date_labels = "%H:%M") +
scale_y_continuous(bquote(omega~(rad~s^{-1}))) +
theme_classic(base_size = 10) +
theme(axis.title.x = element_blank(),
plot.margin = margin(l = 0.5, unit = "cm"))
p_motionless <- cowplot::plot_grid(depth_profile, ygyr_profile,
ncol = 1,
align = "v",
labels = c("B", "C"),
label_size = 10)
cowplot::plot_grid(
cowplot::plot_grid(p_depth, p_motionless,
nrow = 1,
labels = c("A", ""),
label_size = 10),
bw_bcg_plot(23),
ncol = 1,
rel_heights = c(2, 3)
)
labels <- bw_psd %>%
filter(between(freq_bpm, 7.99, 8.01)) %>%
group_by(metric) %>%
slice(1) %>%
ungroup() %>%
left_join(psd_s2n, by = "metric") %>%
mutate(s2n_text = sprintf("Signal:noise = %0.2f", s2n))
bcg_labels <- bw_psd %>%
filter(between(freq_bpm, 7.99, 8.01)) %>%
group_by(metric) %>%
slice(1) %>%
ungroup() %>%
mutate(bcg_lbl = c("1D BCG", "3D BCG"),
spec = c(1.75e-7, 3.2e-9))
pA <- ggplot(bw_psd, aes(x = freq_bpm, y = spec)) +
geom_area(data = filter(bw_psd, between(freq_bpm, 4, 8)), fill = "grey80") +
geom_path() +
geom_text(aes(label = s2n_text), labels, hjust = -0.1, vjust = -0.05) +
geom_text(aes(label = bcg_lbl), bcg_labels, hjust = 0, vjust = 1) +
labs(x = "Frequency (bpm)",
y = "Power spectral density") +
facet_grid(rows = vars(metric),
scales = "free_y") +
theme_classic() +
theme(strip.background = element_blank(),
strip.text = element_blank())
pB <- ggplot(bpm_data, aes(dive_norm, bcg_bpm)) +
geom_point(size = 0.5) +
geom_line(data = bpm_pred, color = "blue", size = 1) +
scale_x_continuous() +
labs(x = "Dive time (normalized)",
y = "Heartrate (bpm)") +
expand_limits(y = 0) +
theme_classic()
pC <- ggplot(ecg_bpm, aes(ecg_bpm, bcg_bpm)) +
geom_abline(slope = 1, intercept = 0, size = 1, linetype = "dashed") +
geom_smooth(method = "lm", se = FALSE, size = 1) +
geom_point(shape = 21) +
scale_x_continuous("ECG (bpm)", breaks = seq(55, 70, by = 5)) +
scale_y_continuous("BCG (bpm)", breaks = seq(55, 70, by = 5)) +
coord_fixed() +
theme_classic() +
theme(axis.title.x = element_text())
cowplot::plot_grid(
pA,
cowplot::plot_grid(pB, pC, ncol = 1, labels = c("B", "C")),
nrow = 1,
labels = c("A", "")
)
References
::: {#refs}
:::
\newpage
Colophon
This report was generated on r Sys.time() using the following computational environment and dependencies:
# which R packages and versions?
if ("devtools" %in% installed.packages()) devtools::session_info()
The current Git commit details are:
# what commit is this file at?
if ("git2r" %in% installed.packages() & git2r::in_repository(path = ".")) git2r::repository(here::here())
FlukeAndFeather/cetaceanbcg documentation built on July 7, 2022, 12:36 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Keywords: r rmarkdown::metadata$keywords
knitr::opts_chunk$set( collapse = TRUE, warning = FALSE, message = FALSE, echo = FALSE, cache = TRUE, comment = "#>", fig.path = "../figures/", fig.height = 6, fig.width = 6, dev = "pdf", dpi = 300 ) library(cetaceanbcg) library(tidyverse)
Introduction
Recent advances in physio-logging (recording physiological variables using animal-borne devices) have largely been driven by new developments in sensor technology [@hawkesIntroductionThemeIssue2021; @williams2021]. For example, new physio-logging tags can detect regional changes in blood flow by incorporating functional near-infrared spectroscopy sensors [@mcknightShiningNewLight2021]. However, traditional inertial measurement unit (IMU) tags equipped with accelerometers and other inertial sensors can also measure important physiological and related variables, such as wing beat frequency [@patterson2019] and feeding rate [@divirgilio2018]. Through careful inspection and analysis of high-resolution acceleration, scientists have measured elevated respiration rates following record-breaking dives [@satoStrokeRatesDiving2011], near-continuous feeding by small cetaceans [@wisniewskaUltraHighForagingRates2016], social interactions between large cetaceans [@goldbogenUsingAccelerometersDetermine2014], and important biomechanical variables including movement speed [@cadeDeterminingForwardSpeed2018]. While physio-logging tags with cutting-edge biomedical technologies push the boundaries of physiological field research, simpler IMU tags have fewer logistical constraints and provide access to more species and larger sample sizes. This is particularly important for species that cannot be restrained or studied in managed care. For example, of the sixteen species of baleen whales (Mysticeti), heart rate has only been recorded with an electrocardiogram tag in the wild for one blue whale (Balaenoptera musculus) [@goldbogenExtremeBradycardiaTachycardia2019; but see @ponganisHeartRateElectrocardiogram1999]. Conversely, IMU tags have been deployed on hundreds of individuals of nearly every species in the clade for the last twenty years [@nowacekBuoyantBalaenidsUps2001]. These existing datasets (and future IMU tag deployments) could hold additional valuable physiological information, awaiting proper computational methods for mining them.
The ballistocardiogram (BCG) has potential applications to using accelerometers as heartrate monitors in both the wild and in managed care [@giovangrandi2011ballistocardiography; @sadekBallistocardiogramSignalProcessing2019; @inanBallistocardiographySeismocardiographyReview2015]. Ballistocardiography is a noninvasive method for measuring cardiac function based on the ballistic forces involved in the heart ejecting blood into the major vessels. The BCG originated as a clinical tool in the first half of the 20th century [@starrStudiesEstimationCardiac1939], but was largely superseded by electro- and echocardiography. However, potential novel applications like passive monitoring of heart function in at-risk populations [@giovangrandi2011ballistocardiography] has led to a recent resurgence of ballistocardiography research, with advances in hardware [@andreozzi2021] and signal processing methodology [@sadekBallistocardiogramSignalProcessing2019]. While the BCG is a three-dimensional phenomenon, it is strongest in the cranio-caudal axis [@inanBallistocardiographySeismocardiographyReview2015]. Along this axis, the waveform is composed of multiple peaks and valleys; most prominent of these is the so-called IJK complex [@pinheiroTheoryDevelopmentsUnobtrusive2010]. The precise physiological mechanism underlying the BCG waveform has not been fully resolved [@kim2016], but it has been established that the IJK complex occurs during systole and, in humans, occurs at approximately the same time as the T-wave in an electrocardiogram (ECG) [@inanBallistocardiographySeismocardiographyReview2015]. The BCG J wave is the most robust feature in the waveform and is typically used for detecting heart beats [@inanBallistocardiographySeismocardiographyReview2015].
Here we present a method for generating a BCG from bio-logger cranio-caudal acceleration. We validated our method with a simultaneously recorded ECG on an adult killer whale in managed care (Orcinus orca) and applied it to detect heartrate in a blue whale. The relative orientation of a tag on a cetacean's body is often uncertain when bio-loggers are deployed in the wild [@johnsonDigitalAcousticRecording2003], so isolating acceleration along the cranio-caudal axis is subject to error. Therefore, we also compared a tri-axial BCG to the cranio-caudal BCG. Specifically, we tested three hypotheses to validate our method. First, a cranio-caudal (1D) BCG would, in a controlled setting, produce instantaneous heartrates that are statistically equivalent to ECG instantaneous heartrates. Second, a tri-axial (3D) BCG would, in a field setting, produce a more robust signal than a 1D BCG. Third, BCG-derived heartrates would increase during the later phases of dives, consistent with the progressive increase in heartrate routinely observed prior to and during ascent [@goldbogenExtremeBradycardiaTachycardia2019; @mcdonaldDeepdivingSeaLions2014].
Materials and methods
Animal tagging
Killer whale
A 3868 kg adult female killer whale in managed care at SeaWorld of California, San Diego, CA was double-tagged with an archival Customized Animal Tracking Solutions IMU (CATS, www.cats.is) tag and a custom-built, archival ECG tag on August 16, 2021 as part of clinical animal cardiac evaluations under the SeaWorld USDA APHIS display permit. The ECG tag hardware and data processing procedures were previously described by @bickettHeartRatesHeart2019. Both tags were deployed by hand and attached with suction cups. We attached the CATS tag on the mid-lateral left chest posterior to the pectoral fin (Movie S1). The CATS tag recorded tri-axial acceleration at 400 Hz, tri-axial magnetometer and tri-axial gyroscope at 50 Hz, pressure at 10 Hz, and video at 30 fps. The IMU in the CATS tag was a MPU-9250 (InvenSense, San Jose, CA; www.invensense.com). The accelerometer had dynamic range of ±4 g, sensitivity of 8,192 LSB g^-1^, and accuracy of 6.1 × 10^-6^ g. All sensors were rotated from the tag's frame of reference to that of the whale using MATLAB (MathWorks, Inc., v2020b) tools for processing CATS data [@cadeToolsIntegratingInertial2021]. This rotation aligned the tag's x-, y-, and z- axes with the cranio-caudal, lateral, and dorso-ventral axes of the whale, respectively. We attached the ECG tag approximately midline on the ventral chest just caudal (posterior) to the axilla and we recorded the ECG at 100 Hz. Individual heart beats in the ECG record were identified from visually verified R-waves using a customized peak detection program (K. Ponganis; Origin 2017, OriginLab Co., Northampton, MA). ECG and IMU were recorded during a spontaneous breath hold while the whale rested at the surface.
Blue whale
A 24.5 m blue whale was tagged with an archival, suction-cup CATS IMU tag on September 5, 2018 in Monterey Bay, CA under permits MBNMS-MULTI-2017-007, NMFS 21678, and Stanford University IACUC 30123 [previously published by @goughScalingSwimmingPerformance2019]. We deployed the tag using a 4 m fiberglass pole from a 6.3 m rigid-hulled inflatable boat and recovered it via radio VHF tracking [as described by @goldbogen2006]. The tag slid behind the left pectoral flipper, similar to the placement of the CATS tag on the killer whale. Tag configuration and data processing followed the same procedure as the killer whale, including accelerometer specification and sampling rates for inertial sensors and video. The 400 Hz acceleration data was used for ballistocardiography (see section Signal processing). We downsampled the multi-sensor data to 10 Hz for movement analysis using the MATLAB CATS tools [@cadeToolsIntegratingInertial2021].
Signal processing
The BCG waveform is three dimensional, but strongest in the cranio-caudal axis [@inanBallistocardiographySeismocardiographyReview2015]. We tested both 1-dimensional (cranio-caudal only) and 3-dimensional metrics for identifying heartbeats in acceleration data based on the methods of @leePhysiologicalSignalMonitoring2016. For windowed operations (such as moving averages and signal filters), we used 0.5 s windows for killer whale data and 2.0 s windows for blue whale data, corresponding to 200 and 800 data points, respectively. The different window sizes were determined through trial and error to remove noise while retaining signal shape. Generally, longer windows will be necessary for larger animals due to their slower heartrates [@Stahl1967b].
Procedure
- First we removed noise and de-trended the acceleration signal with a 5th order Butterworth band-pass filter (killer whale: [1-25 Hz], blue whale: [1-10 Hz]) (R package
signal) [@R-signal]. The lower cut-off frequency de-trended the data. 1 Hz should be appropriate for most marine mammal species. The upper cut-off frequency removed noise. A suitable upper cut-off frequency will depend on body size; larger species' bodies produce lower magnitude accelerations [@martínlópez2021], so more conservative upper cut-off frequencies may be applied to remove more noise without sacrificing signal shape clarity. - Then we enhanced the IJK complex by differentiating acceleration using a 4th order Savitzky-Golay filter (R package
signal). Differentiation (i.e., $a_{t+1}-a_t$, where $a_t$ is the observed acceleration at time step $t$) exaggerates peaks, like the J wave, but it is sensitive to noisy signals. Therefore, additional noise reduction is necessary prior to differentiation. A moving average smoother could remove noise, but it would also reduce the amplitude of peaks. Hence, differentiating Savitzky-Golay filters are preferred in peak-detection algorithms because they remove noise while retaining the general shape of peaks [@samann2019]. We described the resulting signal as "differenced acceleration", rather than jerk, because we did not take the derivative of acceleration with respect to time. The purpose of this signal was to exaggerate a phenomenon in the signal (i.e., the J wave), not to describe a physical quantity (i.e., jerk). - We further enhanced the peaks in the differenced acceleration signal by calculating the Shannon entropy ($H_i=-\sum_{k} |a_{ik}| \times ln(|a_{ik}|)$, where $k$ is the acceleration axis). Additionally, Shannon entropy is strictly positive, which facilitated peak detection. In the 1D BCG, $k$ was surge (cranio-caudal acceleration). In the 3D BCG, $k$ included surge, sway (lateral acceleration), and heave (dorso-ventral acceleration).
- After enhancing the peaks through differentiation and entropy calculation, we removed residual noise by applying a triangular moving average (TMA) smoother. TMAs are equivalent to applying a simple moving average in two passes, which applies greater weight to the middle part of the window and retains peaks and valleys more clearly. After steps 2 and 3, the signal was clear enough that TMAs provided satisfactory results, obviating the need for a more complex algorithm like a Savitzky-Golay filter at this stage. We described the resulting signal as the BCG.
- The BCG contained major peaks (corresponding to heartbeats) and minor peaks (residual noise) (Fig S1A). We extracted all peaks from the BCG and applied a clustering algorithm to retain major peaks and reject minor peaks. First, we extracted all peaks in the BCG signal using
findpeaks()(R packagepracma) [@R-pracma] with a minimum peak distance equivalent to the window size (0.5 s for the killer whale, 2.0 s for the blue whale). For each peak, we calculated its absolute height and its prominence (i.e., height relative to the lowest valley between a peak and its higher neighbors). Then, we calculated each peak's Euclidean distance in height-prominence space from the highest peak (Fig S1B) and estimated the density distribution of these distances (Fig S1C). The density distribution was bimodal, with a low-distance peak corresponding to major peaks and a high-distance peak corresponding to minor peaks. We used the distance corresponding to the valley between the two peaks as a threshold for rejecting minor peaks (Fig S1D).
This procedure may be applied to either 1-dimensional (i.e., cranio-caudal only) or 3-dimensional acceleration. In the case of 3-dimensional acceleration, the band-pass and Savitzky-Golay filters were applied to each axis independently.
BCG validation with killer whale ECG
We fit ordinary least squares regression to BCG-derived instantaneous heartrates with respect to ECG-derived heartrate and tested 1) if the intercept was significantly different than 0 and 2) if the slope was significantly different than 1. We calculated the mean and standard deviation of absolute error as an equivalence measure (1-dimensional BCG only).
BCG application to blue whale
Dynamic body movements produce an acceleration signal that masks the ballistocardiogram, so we limited our analyses to motionless periods (Fig. \@ref(fig:bw-bcg)B-C). These periods occurred during or near the bottom phase of dives between fluke strokes. Strokes were detected from visual examination of the rotational velocity around the lateral axis recorded by gyroscope [sensu @goughScalingSwimmingPerformance2019]. We used gyroscopes for stroke detection because they were 1) separate sensors from the accelerometers and 2) strokes are clearly visible in gyroscope signals and are robust to tag placement.
We tested whether the 3-dimensional BCG was more robust than 1-dimensional BCG in field data by comparing the signal-to-noise ratios. For both BCGs, we calculated the power spectral density (R package psd) [@Barbour2014]. Previously recorded blue whale apneic heart rate was 4-8 beats per minute (bpm) [@goldbogenExtremeBradycardiaTachycardia2019], so we quantified signal as the integration of the power spectral density curve from 4-8 bpm and noise as the integrated remainder, up to 60 bpm. The sample size recorded by @goldbogenExtremeBradycardiaTachycardia2019 was one individual, so we could not account for potential inter-individual variation. Nonetheless, 4-8 bpm was the best available estimate for typical blue whale apneic heart rate.
We also tested whether BCG-derived instantaneous heart rates were consistent with the range and pattern of heart rates previously observed in the blue whale and other marine mammals; namely a gradual increase in heart rate later in the dive, especially during the final ascent [@goldbogenExtremeBradycardiaTachycardia2019; @mcdonaldDeepdivingSeaLions2014]. We assigned dive start and end times when the whale swam deeper than 2 m, retaining dives that exceeded 10 m depth and 5 minutes duration. Dive times were normalized from 0 (start of dive) to 1 (end of dive). We regressed instantaneous heart rate against normalized dive time using robust Theil-Sen regression (to account for heteroscedascity) (R package RobustLinearReg) [@R-RobustLinearReg; @Sen-1968; @theilRankInvariantMethodLinear1992] and tested whether the slope was greater than 0.
Reproducibility
The data and code used in this analysis were packaged as a research compendium, containing the data, code, and an executable version of this manuscript. We used the R package rrtools [@marwickPackagingDataAnalytical2018; @rrtools2019] to initialize the compendium, which was written as an R package. This approach promotes reproducibility and facilitates adoption by other researchers [@alstonBeginnerGuideConducting2021; @powersOpenScienceReproducibility2019; @stoddenEmpiricalAnalysisJournal2018]. The steps described in section Signal processing were implemented as functions in the R package, and the executable manuscript demonstrates how to use those functions to perform the analyses presented in this study.
Results and discussion
BCG validation with killer whale ECG
ca_tz <- "Etc/GMT+7" bcg_begin <- lubridate::ymd_hms("2021-08-16 10:09:24.5", tz = ca_tz) bcg_end <- lubridate::ymd_hms("2021-08-16 10:09:38", tz = ca_tz) oo_400hz <- corky_400hz %>% filter(between(dt, bcg_begin - 10, bcg_end + 10)) %>% mutate(surge_filt = filter_acc(surge, 400, 25.0), surge_diff = lead(surge_filt) - surge_filt, surge_se = shannon_entropy(surge_diff), surge_smooth = tma(surge_se, 0.5 * 400)) %>% filter(between(dt, bcg_begin, bcg_end)) %>% mutate(bcg_beat = find_beats(surge_smooth, 400, 0.5), bcg_bpm = bpm(bcg_beat, dt))
# Align ECG-detected heart beats oo_400hz$ecg_beat <- FALSE beat_idx <- approx(oo_400hz$dt, seq_along(oo_400hz$dt), xout = ecg_beats$dt, method = "constant")$y %>% na.omit() oo_400hz$ecg_beat[beat_idx] = TRUE oo_400hz$ecg_bpm <- with(oo_400hz, bpm(ecg_beat, dt))
bcg_bpm <- filter(oo_400hz, bcg_beat) %>% select(dt, bcg_bpm) nearest_bcg <- function(t) { bcg_bpm$bcg_bpm[which.min(abs(t - bcg_bpm$dt))] } ecg_bpm <- filter(oo_400hz, ecg_beat) %>% mutate(bcg_bpm = map_dbl(dt, nearest_bcg), bpm_diff = abs(bcg_bpm - ecg_bpm) / ecg_bpm) %>% drop_na(bpm_diff) mean_bpm_diff <- mean(ecg_bpm$bpm_diff) sd_bpm_diff <- sd(ecg_bpm$bpm_diff) bpm_lm <- lm(bcg_bpm ~ ecg_bpm, ecg_bpm) bpm_stderr <- sqrt(diag(vcov(bpm_lm)))
The ECG and BCG yielded nearly identical heart rate estimations (Fig. \@ref(fig:oo-bcg-ecg)). We collected 14 s of simultaneous ECG and BCG data during a motionless breath hold at the surface. Logistical constraints prevented us from gathering a longer sample, as these data were collected secondary to other projects. BCG-derived instantaneous heart rates were within r sprintf("%0.1f%% \u00B1 %0.1f%%", mean_bpm_diff * 100, sd_bpm_diff * 100) of the ECG-derived rates (mean ± standard deviation). Ordinary least squares regression of BCG heartrates on ECG heartrates yielded a slope of r sprintf("%0.2f \u00B1 %0.2f", coef(bpm_lm)[2], bpm_stderr[2]) and intercept of r sprintf("%0.2f \u00B1 %0.2f", coef(bpm_lm)[1], bpm_stderr[1]) (mean ± standard error), which were not significantly different from the hypothesized 1 and 0, respectively. (Fig. \@ref(fig:validation-plots)C).
BCG application to blue whale
bw_elg <- bw180905_53_elg %>% mutate( depth_min = map2_dbl(start, stop, function(t1, t2) { bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>% min() }), depth_max = map2_dbl(start, stop, function(t1, t2) { bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>% max() }) ) n_motionless <- nrow(bw_elg) motionless_minutes <- as.numeric(sum(bw_elg$stop - bw_elg$start), unit = "mins")
We generated 1-dimensional and 3-dimensional BCGs for 2 hours of data, including 10 rest dives and r n_motionless motionless periods totaling r format(motionless_minutes, digits = 3) minutes (r sprintf("%0.1f%%", motionless_minutes / 120 * 100) of the 2-hour record) (Fig. \@ref(fig:bw-bcg)A-C).
bw_400hz <- bw180905_53_400hz %>% mutate( across(surge:heave, filter_acc, fs = 400, upper = 10.0, .names = "{.col}_filt"), surge_diff = lead(surge_filt) - surge_filt, surge_se = shannon_entropy(surge_diff), surge_smooth = tma(surge_se, 0.5 * 400), jerk = jerk(cbind(surge_filt, sway_filt, heave_filt), fs = 400, p = 4, n = 2 * 400 + 1), jerk_se = shannon_entropy(jerk), jerk_smooth = tma(jerk_se, 2 * 400), # Annotate regions rid_left = approx(bw_elg$start, bw_elg$region_id, dt, "constant")$y, rid_right = approx(bw_elg$stop, bw_elg$region_id, dt, "constant", yleft = 0)$y + 1, region_id = ifelse(rid_left == rid_right, rid_left, NA), # Zero-out signal in non-valid regions (i.e. remove movement artifacts) jerk_smooth = ifelse(is.na(region_id), 0, jerk_smooth), bcg_beat = find_beats(jerk_smooth, 400, 2) ) %>% # Calculate heart rate within each region group_by(region_id) %>% mutate(bcg_bpm = bpm(bcg_beat, dt)) %>% ungroup() %>% select(-c(rid_left, rid_right))
walk(unique(na.omit(bw_400hz$region_id)), function(rid) { p <- bw_bcg_plot(rid) ggsave(here::here(sprintf("analysis/figures/bw_bcg/bw_bcg_%02i.png", rid)), p) })
The 3-dimensional BCG (Fig. \@ref(fig:bw-bcg)) produced a more robust signal (i.e., higher signal-to-noise ratio) than the 1-dimensional BCG, which used only cranio-caudal acceleration (Fig. \@ref(fig:validation-plots)A). The signal-to-noise ratio was 2.00 for the 3-dimensional BCG, compared to 0.17 for the 1-dimensional BCG. Although the power spectral density curve for the 1-dimensional BCG had a peak in the 4-8 bpm frequency range, most of the signal's power was concentrated in lower frequencies. Conversely, the 3-dimensional BCG's power was concentrated precisely in the 4-8 bpm frequency range, with only a smaller peak in the lower frequencies.
bcg1d_psd <- with(bw_400hz, pspectrum(surge_smooth[!is.na(region_id)], 400)) bcg3d_psd <- with(bw_400hz, pspectrum(jerk_smooth[!is.na(region_id)], 400)) bw_psd <- tibble( metric = "Surge only", freq_hz = bcg1d_psd$freq, spec = bcg1d_psd$spec ) %>% rbind( tibble( metric = "Triaxial", freq_hz = bcg3d_psd$freq, spec = bcg3d_psd$spec ) ) %>% mutate(freq_bpm = freq_hz * 60) %>% filter(freq_bpm <= 60) psd_s2n <- bw_psd %>% group_by(metric) %>% summarize(signal = pracma::trapz(freq_bpm[freq_bpm >= 4 & freq_bpm <= 8], spec[freq_bpm >= 4 & freq_bpm <= 8]), noise = pracma::trapz(freq_bpm, spec) - signal) %>% mutate(s2n = signal / noise)
bw10hz <- bw180905_53_10hz %>% mutate( dive_id = split_dives(dt, depth, surface = 2, min_depth = 10, min_dur = 5 * 60), dive_norm = normalize_dives(dive_id), rid_left = approx(bw_elg$start, bw_elg$region_id, dt, "constant")$y, rid_right = approx(bw_elg$stop, bw_elg$region_id, dt, "constant", yleft = 0)$y + 1 ) bw_400hz$dive_norm <- approx(bw10hz$dt, bw10hz$dive_norm, bw_400hz$dt)$y bpm_data <- bw_400hz %>% drop_na(bcg_bpm) bpm_ts <- RobustLinearReg::theil_sen_regression(bcg_bpm ~ dive_norm, bpm_data) bpm_pred <- tibble(dive_norm = seq(0, 1, length.out = 100)) %>% mutate(bcg_bpm = predict(bpm_ts, newdata = .))
The 3-dimensional BCG exhibited increasing heart rates over the course of dives. Average heart rate increased from r format(bpm_pred$bcg_bpm[1], digits = 2) bpm at the start of dives to r format(bpm_pred$bcg_bpm[nrow(bpm_pred)], digits = 2) bpm at the end of dives (Theil-Sen regression, $p < 10^{-10}$) (Fig. \@ref(fig:validation-plots)B).
Limitations and considerations for future applications
While the BCG method presented here holds the potential to mine existing and future marine mammal bio-logging datasets for information about cardiovascular function, it has several limitations compared to ECG methods. Most importantly, BCGs are highly sensitive to movement artifacts [@inanBallistocardiographySeismocardiographyRemovview2015], so only motionless periods are valid for analysis. This limits the behavioral and physiological contexts in which heartrate may be measured. For example, the BCG is probably an inappropriate method for quantifying the magnitude of surface tachycardia [@goldbogenExtremeBradycardiaTachycardia2019] and exercise modulation of bradycardia [@noren2012], due to movement artifacts during those activities. Secondarily, we did not test whether the BCG is robust to tag placement location. The blue whale data presented in this study was collected when a dorsally-deployed tag slipped to the lateral chest cavity behind a flipper, where it is reasonable to expect greater accelerations caused by heart beats than from a tag farther from the animal's center of mass. It is possible that the ballistic forces generated by heart beats are strong enough to produce an interpretable BCG for a variety of potential tag deployment locations, but this likely varies with animal body size, as well as accelerometer sampling rate and sensitivity.
When auditing existing bio-logging data and planning future tag deployments for BCG analysis, careful consideration should be paid to sampling rate. As a rule of thumb in signal processing, the sampling rate should be at least twice the frequency of the phenomenon of interest. In the case of the BCG, the relevant frequency is that of the BCG waveform, not the heartrate. In humans, the power of the IJK-complex (the part of the BCG waveform used for heart beat detection) is concentrated between 4-7 Hz [@moukadem2018]. It is unlikely that marine mammal BCG waveforms have a higher frequency than humans, owing to their generally larger body sizes. Therefore, it is possible that BCGs may be generated for accelerometer sampling rates as low as 10-15 Hz. Conservatively, the authors recommend a sampling rate of no less than 50 Hz (i.e., twice the upper cut-off frequency of the widest bandpass filter used in this study).
Future bio-logging BCG methodology research should address the limitations imposed by tag-placement and movement artifacts. We used accelerometers in this study because of their prevalence in bio-logging research, but it is possible that other widely used bio-logger sensors, such as gyroscopes and/or magnetometers, could produce a clearer signal in a greater variety of contexts. Alternative bio-logger housing designs, such as limpet-style tags [@andrews2008] or "marine skin" [@nassar2018], could reduce noise, boost the signal-to-noise ratio, and make the method more widely applicable.
Conclusions
Here we presented a ballistocardiogram method for detecting resting apneic heartrate in cetaceans using accelerometers. We validated the method in a controlled setting with simultaneous ECG and in a field setting by confirming expected physiological patterns. As accelerometer tags have been deployed on many cetacean species for multiple decades, this method may be applied to mine existing datasets and better understand how heartrate scales with body size and other biological factors. It may also provide additional data for conservation physiology applications. For example, BCGs extracted from gliding phases before and after controlled sonar exposure experiments could quantify the physiological response to anthropogenic disturbance [@southall2019]. Even as the field of physio-logging progresses with new hardware innovations [@williams2021; @fahlman2021; @hawkesIntroductionThemeIssue2021], this method demonstrates that computational advances can derive new insights from traditional sensors.
Acknowledgements
The authors are grateful to the SeaWorld of California Killer Whale training staff for their efforts and support. We also thank Anna Krystalli, Ben Marwick, Karthik Ram, Nicholas Tierney, and other members of the R community for developing tools and educational resources to facilitate open science practices. This is a SeaWorld Parks and Entertainment Technical Contribution number 2021-12. We thank Lucy Hawkes and two anonymous reviewers for comments on the manuscript.
Footnotes
Author contributions
Conceptualization: M.F.C.,J.A.F.,P.J.P.,J.A.G.; Methodology: M.F.C.,J.A.F.,P.J.P.; Software: M.F.C.; Formal analysis: M.F.C.,P.J.P.; Investigation: M.F.C.,J.A.F.,P.J.P.; Resources: P.J.P.,J.A.G.,T.L.S.; Writing - original draft: M.F.C; Writing - review & editing: M.F.C.,J.A.F.,P.J.P.,J.A.G.,T.L.S.; Supervision: P.J.P.,J.A.G.; Project administration: P.J.P.,J.A.G.; Funding acquisition: J.A.G.
Funding
This work was supported by grant N000141912455 from the Office of Naval Research. M.F.C. was supported by the Stanford University William R. and Sara Hart Kimball Fellowship and a Stanford Data Science Scholar Fellowship.
Data availability
All data and code used in this analysis are available on Zenodo [@cetaceanbcgzenodo].
Competing interests
The authors declare no competing interests.
\newpage
Figures
t <- theme_classic(base_size = 8) + theme(axis.title.x = element_blank(), panel.grid.minor = element_blank(), strip.background = element_blank()) t2 <- t + theme(plot.margin = margin(l = 0.5, unit = "cm")) t_yvoid <- t2 + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), strip.text = element_blank()) # A: ECG ecg_sample <- corky_ecg %>% filter(between(dt, first(oo_400hz$dt), last(oo_400hz$dt))) pA <- ggplot(ecg_sample, aes(dt, ecg_uV)) + geom_line(size = 0.2) + geom_point(data = filter(ecg_sample, ecg_beat), color = "red") + scale_x_datetime(date_labels = "%H:%M:%S") + labs(y = bquote(ECG~(mu * V))) + expand_limits(x = max(oo_400hz$dt) + 1) + t2 # B: Filtered surge ijk_range <- as.POSIXct("2021-08-16 10:09:31.4", tz = "US/Pacific") + c(0, 0.4) ijk_complex <- oo_400hz %>% filter(between(dt, ijk_range[1], ijk_range[2])) %>% mutate(t = as.numeric(dt - min(dt), unit = "secs")) ijk_labels <- slice(ijk_complex, c(41, 53, 67, 85, 110)) %>% mutate(wave = c("H", "I", "J", "K", "L")) pB_inset <- ggplot(ijk_complex, aes(dt, surge_filt)) + geom_line() + geom_text(aes(label = wave), ijk_labels, vjust = c(-0.2, 1.2, 0.75, -0.6, -0.2), hjust = c(0.5, 0.5, -0.6, 0.5, 0.5), fontface = "bold", size = 3) + theme_void() inset_range <- as.POSIXct("2021-08-16 10:09:29.9", tz = "US/Pacific") + c(0, 4) pB <- ggplot(oo_400hz, aes(dt, surge_filt)) + geom_line(size = 0.2) + annotation_custom(grob = ggplotGrob(pB_inset), xmin = inset_range[1], xmax = inset_range[2], ymin = 0.25, ymax = 0.5) + annotate("rect", xmin = ijk_range[1], xmax = ijk_range[2], ymin = min(ijk_complex$surge_filt) - 0.02, ymax = max(ijk_complex$surge_filt) + 0.02, fill = NA, color = "grey50", lty = 3) + scale_x_datetime(date_labels = "%H:%M:%S") + expand_limits(y = 0.5) + labs(y = "Filtered\nacceleration") + t_yvoid # C: Surge diff pC <- ggplot(oo_400hz, aes(dt, surge_diff)) + geom_line(size = 0.2) + scale_x_datetime(date_labels = "%H:%M:%S") + labs(y = "Differenced\nacceleration") + t_yvoid # D: Shannon entropy pD <- ggplot(oo_400hz, aes(dt, surge_se)) + geom_line(size = 0.2) + scale_x_datetime(date_labels = "%H:%M:%S") + labs(y = "Shannon\nentropy") + t_yvoid # E: Smoothed pE <- ggplot(oo_400hz, aes(dt, surge_smooth)) + geom_line(size = 0.2) + geom_point(data = filter(oo_400hz, bcg_beat), color = "blue") + scale_x_datetime(date_labels = "%H:%M:%S") + labs(y = "BCG") + t_yvoid cowplot::plot_grid(pA, pB, pC, pD, pE, align = "v", axis = "lr", ncol = 1, rel_heights = c(1, 1.5, 0.75, 0.75, 1), labels = "AUTO")
# Plot the BCG for a motionless period in blue whale data bw_bcg_plot <- function(rid) { bcg_sample <- filter(bw_400hz, region_id == rid) palette <- rgb(c(230, 86, 0), c(159, 180, 158), c(0, 233, 115), maxColorValue = 255) t_lmargin <- t_yvoid + theme(plot.margin = margin(l = 0.5, unit = "cm")) # A: Filtered acceleration acceleration_data <- bcg_sample %>% rename(`Cranio-caudal` = surge_filt, Lateral = sway_filt, `Dorso-ventral` = heave_filt) %>% pivot_longer(`Cranio-caudal`:`Dorso-ventral`, names_to = "axis", values_to = "acc") %>% mutate(axis = factor(axis, levels = c("Cranio-caudal", "Lateral", "Dorso-ventral"))) pA <- ggplot(acceleration_data, aes(dt, acc, color = axis)) + geom_line(size = 0.2) + scale_x_datetime(date_labels = "%H:%M:%S") + scale_y_continuous(n.breaks = 3) + scale_color_manual(values = palette) + facet_grid(rows = vars(axis)) + labs(y = "Filtered\nacceleration") + t_lmargin + theme(legend.position = "none", panel.spacing = unit(0, "lines")) # B: Jerk vectors jerk_data <- bcg_sample %>% mutate(jerk = set_names(as_tibble(jerk), c("Cranio-caudal", "Lateral", "Dorso-ventral"))) %>% unpack(jerk) %>% pivot_longer(`Cranio-caudal`:`Dorso-ventral`, names_to = "axis", values_to = "jerk") %>% mutate(axis = factor(axis, levels = c("Cranio-caudal", "Lateral", "Dorso-ventral"))) pB <- ggplot(jerk_data, aes(dt, jerk, color = axis)) + geom_line(size = 0.2) + scale_x_datetime(date_labels = "%H:%M:%S") + scale_y_continuous(n.breaks = 3) + scale_color_manual(values = palette) + facet_grid(rows = vars(axis)) + labs(y = "Differenced\nacceleration") + t_lmargin + theme(legend.position = "none", panel.spacing = unit(0, "lines")) # C: Shannon entropy of jerk (NOT jerk norm) pC <- ggplot(bcg_sample, aes(dt, jerk_se)) + geom_line(size = 0.2) + scale_x_datetime(date_labels = "%H:%M:%S") + labs(y = "Shannon\nentropy") + t_lmargin # D: Smoothed pD <- ggplot(bcg_sample, aes(dt, jerk_smooth)) + geom_line(size = 0.2) + geom_point(data = filter(bcg_sample, bcg_beat), color = "blue") + scale_x_datetime(date_labels = "%H:%M:%S") + labs(y = "BCG") + t_lmargin cowplot::plot_grid(pA, pB, pC, pD, align = "v", axis = "lr", ncol = 1, rel_heights = c(1.5, 1.5, 1, 1), labels = c("D", "E", "F", "G"), label_size = 10) } # Plot the depth profile and motionless periods for the blue whale bw_10hz <- bw180905_53_10hz %>% mutate(rid_left = approx(bw_elg$start, bw_elg$region_id, dt, "constant")$y, rid_right = approx(bw_elg$stop, bw_elg$region_id, dt, "constant", yleft = 0)$y + 1, region_id = ifelse(rid_left == rid_right, rid_left, NA), dive_id = split_dives(dt, depth, surface = 2, min_depth = 10, min_dur = 5 * 60)) divesample <- filter(bw_10hz, dive_id == 5) p_depth <- ggplot(bw_10hz, aes(dt, depth)) + geom_line(size = 0.2) + annotate("rect", xmin = min(divesample$dt) - 60, xmax = max(divesample$dt) + 60, ymin = -1, ymax = max(divesample$depth) + 1, fill = NA, color = "red", linetype = "dotted") + scale_x_datetime(date_labels = "%H:%M") + scale_y_reverse("Depth (m)") + theme_classic(base_size = 10) + theme(axis.title.x = element_blank()) # Plot the depth and y-axis gyroscope for a motionless period norm_y <- function(ymin, ymax, sign) { (ymin + ymax) / 2 + sign * max(ymax - ymin) / 2 } bw_elg <- bw180905_53_elg %>% mutate( depth_min = map2_dbl(start, stop, function(t1, t2) { bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>% min() }), depth_max = map2_dbl(start, stop, function(t1, t2) { bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>% max() }), depth_med = map2_dbl(start, stop, function(t1, t2) { bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), "depth"] %>% first() %>% median() }), ygyr_min = map2_dbl(start, stop, function(t1, t2) { bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), ]$gw[,2] %>% min() }), ygyr_max = map2_dbl(start, stop, function(t1, t2) { bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), ]$gw[,2] %>% max() }), ygyr_med = map2_dbl(start, stop, function(t1, t2) { bw180905_53_10hz[between(bw180905_53_10hz$dt, t1, t2), ]$gw[,2] %>% first() %>% median() }), ) diveelg <- bw_elg %>% filter(start >= first(divesample$dt), stop <= last(divesample$dt)) %>% mutate(depth_min2 = norm_y(depth_min, depth_max, -1), depth_max2 = norm_y(depth_min, depth_max, 1), ygyr_min2 = norm_y(ygyr_min, ygyr_max, -1), ygyr_max2 = norm_y(ygyr_min, ygyr_max, 1)) depth_profile <- ggplot(divesample, aes(dt, depth)) + geom_line(size = 0.5) + geom_rect(aes(xmin = start, xmax = stop, ymin = depth_min2, ymax = depth_max2), diveelg, inherit.aes = FALSE, fill = "orchid1", color = NA, alpha = 0.5) + geom_rect(aes(xmin = start, xmax = stop, ymin = depth_min2, ymax = depth_max2), filter(diveelg, region_id == 23), inherit.aes = FALSE, fill = NA, color = "red", linetype = "dotted", alpha = 0.5) + scale_x_datetime(date_labels = "%H:%M") + scale_y_reverse("Depth (m)") + theme_classic(base_size = 10) + theme(axis.title.x = element_blank(), plot.margin = margin(l = 0.5, unit = "cm")) ygyr_profile <- ggplot(divesample, aes(dt, gw[, 2])) + geom_line(size = 0.2) + geom_rect(aes(xmin = start, xmax = stop, ymin = ygyr_min2, ymax = ygyr_max2), diveelg, inherit.aes = FALSE, fill = "orchid1", color = NA, alpha = 0.5) + geom_rect(aes(xmin = start, xmax = stop, ymin = ygyr_min2, ymax = ygyr_max2), filter(diveelg, region_id == 23), inherit.aes = FALSE, fill = NA, color = "red", linetype = "dotted", alpha = 0.5) + scale_x_datetime(date_labels = "%H:%M") + scale_y_continuous(bquote(omega~(rad~s^{-1}))) + theme_classic(base_size = 10) + theme(axis.title.x = element_blank(), plot.margin = margin(l = 0.5, unit = "cm")) p_motionless <- cowplot::plot_grid(depth_profile, ygyr_profile, ncol = 1, align = "v", labels = c("B", "C"), label_size = 10) cowplot::plot_grid( cowplot::plot_grid(p_depth, p_motionless, nrow = 1, labels = c("A", ""), label_size = 10), bw_bcg_plot(23), ncol = 1, rel_heights = c(2, 3) )
labels <- bw_psd %>% filter(between(freq_bpm, 7.99, 8.01)) %>% group_by(metric) %>% slice(1) %>% ungroup() %>% left_join(psd_s2n, by = "metric") %>% mutate(s2n_text = sprintf("Signal:noise = %0.2f", s2n)) bcg_labels <- bw_psd %>% filter(between(freq_bpm, 7.99, 8.01)) %>% group_by(metric) %>% slice(1) %>% ungroup() %>% mutate(bcg_lbl = c("1D BCG", "3D BCG"), spec = c(1.75e-7, 3.2e-9)) pA <- ggplot(bw_psd, aes(x = freq_bpm, y = spec)) + geom_area(data = filter(bw_psd, between(freq_bpm, 4, 8)), fill = "grey80") + geom_path() + geom_text(aes(label = s2n_text), labels, hjust = -0.1, vjust = -0.05) + geom_text(aes(label = bcg_lbl), bcg_labels, hjust = 0, vjust = 1) + labs(x = "Frequency (bpm)", y = "Power spectral density") + facet_grid(rows = vars(metric), scales = "free_y") + theme_classic() + theme(strip.background = element_blank(), strip.text = element_blank()) pB <- ggplot(bpm_data, aes(dive_norm, bcg_bpm)) + geom_point(size = 0.5) + geom_line(data = bpm_pred, color = "blue", size = 1) + scale_x_continuous() + labs(x = "Dive time (normalized)", y = "Heartrate (bpm)") + expand_limits(y = 0) + theme_classic() pC <- ggplot(ecg_bpm, aes(ecg_bpm, bcg_bpm)) + geom_abline(slope = 1, intercept = 0, size = 1, linetype = "dashed") + geom_smooth(method = "lm", se = FALSE, size = 1) + geom_point(shape = 21) + scale_x_continuous("ECG (bpm)", breaks = seq(55, 70, by = 5)) + scale_y_continuous("BCG (bpm)", breaks = seq(55, 70, by = 5)) + coord_fixed() + theme_classic() + theme(axis.title.x = element_text()) cowplot::plot_grid( pA, cowplot::plot_grid(pB, pC, ncol = 1, labels = c("B", "C")), nrow = 1, labels = c("A", "") )
References
::: {#refs} :::
\newpage
Colophon
This report was generated on r Sys.time() using the following computational environment and dependencies:
# which R packages and versions? if ("devtools" %in% installed.packages()) devtools::session_info()
The current Git commit details are:
# what commit is this file at? if ("git2r" %in% installed.packages() & git2r::in_repository(path = ".")) git2r::repository(here::here())
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.