Second round of revisions

Author: Max Czapanskiy

analysis/paper/paper.Rmd

@@ -64,7 +64,7 @@ Recent advances in physio-logging (recording physiological variables using anima
 
 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 latter phases of dives, consistent with the progressive increase in heartrate routinely observed prior to and during ascent [@goldbogenExtremeBradycardiaTachycardia2019; @mcdonaldDeepdivingSeaLions2014].
+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
 
@@ -76,7 +76,7 @@ A 3868 kg adult female killer whale in managed care at SeaWorld of California, S
 
 **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 [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].
+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
 
@@ -85,7 +85,7 @@ The BCG waveform is three dimensional, but strongest in the cranio-caudal axis [
 **Procedure**
 
 1.  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.
-2.  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, the signal should be de-noised prior to differentiation. A moving average smoother can 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).
+2.  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).
 3.  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).
 4.  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.
 5.  The BCG contained major peaks (corresponding to heartbeats) and minor peaks (residual noise) (Fig S1**A**). 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 S1**B**) and estimated the density distribution of these distances (Fig S1**C**). 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 S1**D**).
@@ -175,7 +175,7 @@ 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 (Fig. \@ref(fig:bw-bcg)**A-C**).
+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**).
 
 ```{r bw-beats}
 bw_400hz <- bw180905_53_400hz %>%
@@ -219,7 +219,7 @@ walk(unique(na.omit(bw_400hz$region_id)), function(rid) {
 })
 ```
 
-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. The signal-to-noise ratio was 2.00 for the 3-dimensional BCG, compared to 0.17 for the 1-dimensional BCG (Fig. \@ref(fig:validation-plots)**A**). 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.
+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.
 
 ```{r bw-psd}
 bcg1d_psd <- with(bw_400hz, pspectrum(surge_smooth[!is.na(region_id)], 400))
@@ -275,11 +275,11 @@ bpm_pred <- tibble(dive_norm = seq(0, 1, length.out = 100)) %>%
   mutate(bcg_bpm = predict(bpm_ts, newdata = .))
 ```
 
-3-dimensional BCG-derived heart rates exhibited a relaxation of bradycardia 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**).
+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 [@inanBallistocardiographySeismocardiographyReview2015], 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 unlikely to quantify neither the magnitude of surface tachycardia [@goldbogenExtremeBradycardiaTachycardia2019] nor 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.
+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 [@inanBallistocardiographySeismocardiographyReview2015], 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).