Incorporating Todd's feedback

Author: Max Czapanskiy

analysis/paper/ballistocardiogram_v3.docx

@@ -20,7 +20,7 @@ author:
 institute:
   - Hopkins: Hopkins Marine Station, Department of Biology, Stanford University, USA
   - Scripps: Scripps Institution of Oceanography, University of California San Diego, USA
-  - SeaWorld: Veterinary Services, SeaWorld of California, USA
+  - SeaWorld: Animal Health Department, SeaWorld of California, USA
 date: "`r format(Sys.time(), '%d %B, %Y')`"
 output: 
     bookdown::word_document2:
@@ -59,9 +59,9 @@ 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]. 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. 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.
+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. 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. 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 the medical community has recently returned to ballistocardiography as a potential means of passive monitoring of heart function in at-risk populations [@giovangrandi2011ballistocardiography], which has led to substantial progress in signal processing methodology for generating interpretable BCGs [@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 are the IJK complex, which progressively occurs during systole [@pinheiroTheoryDevelopmentsUnobtrusive2010]. The BCG J wave is the most robust feature in the waveform and typically used for detecting heart beats.
+The ballistocardiogram (BCG) has potential applications to using accelerometers as heartrate monitors in both the wild and in managed care. 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 are the IJK complex, which progressively occurs during systole [@pinheiroTheoryDevelopmentsUnobtrusive2010]. The BCG J wave is the most robust feature in the waveform and typically used for detecting heart beats.
 
 Here we present a method for generating a BCG from bio-logger accelerometry. We validated our method with a simultaneously recorded electrocardiogram (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 the tag on the body is uncertain in cetacean bio-logging in the wild [@johnsonDigitalAcousticRecording2003], so in addition to a one-dimensional BCG based solely on cranio-caudal acceleration, we also generated a three-dimensional BCG, which we expected would be more robust in a field setting. Specifically, we tested three hypotheses to validate our method. First, a one-dimensional BCG would, in a controlled setting, produce statistically equivalent instantaneous heartrates as an ECG. Second, a three-dimensional BCG would, in a field setting, produce a more robust signal than a one-dimensional 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].
 
@@ -71,19 +71,19 @@ Here we present a method for generating a BCG from bio-logger accelerometry. We
 
 **Killer whale**
 
-A 3807 kg (**needs confirmation**) killer whale in managed care at the SeaWorld, San Diego, CA marine facilities was double-tagged with Customized Animal Tracking Solutions IMU (CATS, www.cats.is) and electrocardiogram (ECG) tags on August 16, 2021 as part of clinical animal cardiac evaluations under the SeaWorld display permit. We attached the CATS tag on the mid-lateral left chest posterior to the pectoral fin (Movie S1). The CATS tag recorded acceleration at 400 Hz, magnetometer and gyroscope at 50 Hz, pressure at 10 Hz, and video at 30 fps. 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. The ECG tag hardware and data processing followed the methods in [@bickettHeartRatesHeart2019]. Briefly, the tag was attached approximately midline on the ventral chest just caudal (posterior) to the axilla and we recorded the ECG at 100 Hz. Individual heart beats 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.
+A 3868 kg adult female killer whale in managed care at SeaWorld of California, San Diego, CA was double-tagged with Customized Animal Tracking Solutions IMU (CATS, www.cats.is) and electrocardiogram (ECG) tags on August 16, 2021 as part of clinical animal cardiac evaluations under the SeaWorld USDA APHIS display permit. We attached the CATS tag on the mid-lateral left chest posterior to the pectoral fin (Movie S1). The CATS tag recorded acceleration at 400 Hz, magnetometer and gyroscope at 50 Hz, pressure at 10 Hz, and video at 30 fps. 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. The ECG tag hardware and data processing followed the methods in [@bickettHeartRatesHeart2019]. Briefly, the tag was attached approximately midline on the ventral chest just caudal (posterior) to the axilla and we recorded the ECG at 100 Hz. Individual heart beats 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 a 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 in [@goughScalingSwimmingPerformance2019]). 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. 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.
+A 24.5 m blue whale was tagged with a 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]. 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. 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.
 
 ## 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, we used 0.5 s for killer whale data and 2.0 s for blue whale data.
 
 **Procedure**
 
-1.  Remove noise and de-trend the acceleration signal with a 5th order Butterworth band-pass filter (killer whale: [1-25Hz], blue whale: [1-10Hz]) [R package `signal` @R-signal].
+1.  Remove noise and de-trend the acceleration signal with a 5th order Butterworth band-pass filter (killer whale: [1-25Hz], blue whale: [1-10Hz]) (R package `signal`) [@R-signal].
 2.  Enhance the IJK complex by differentiating acceleration using a 4th order Savitzky-Golay filter (R package `signal`). Differentiation exaggerates impulses like the J wave.
 3.  Further enhance the peaks by calculating the Shannon entropy ($H_i=-\sum_{k} |a_{ik}| \times ln(|a_{ik}|)$, where $k$ is the acceleration axis). Additionally, the Shannon entropy is strictly positive, which facilitates peak detection. In the 1-dimensional case, $k$ is surge (cranio-caudal acceleration) only.
 4.  Remove noise by applying a triangular moving average smoother.
@@ -99,13 +99,13 @@ We fit ordinary least squares regression to BCG-derived instantaneous heart rate
 
 Dynamic body movements produce an acceleration signal that masks the ballistocardiogram, so we limited our analyses to motionless periods (Fig. S2). These periods occured 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 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.
+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.
 
-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.
+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 [R package `rrtools`; @marwickPackagingDataAnalytical2018; @rrtools2019]. The research compendium was written as an R package so other researchers can read, run, and modify the methods described here.
+The data and code used in this analysis were packaged as a research compendium (R package `rrtools`) [@marwickPackagingDataAnalytical2018; @rrtools2019]. The research compendium was written as an R package so other researchers can read, run, and modify the methods described here.
 
 # Results and discussion
 
@@ -286,9 +286,7 @@ Here we presented a ballistocardiogram method for detecting resting apneic heart
 
 # Acknowledgements
 
--   Everyone who helped collect and process the blue whale data.
--   The Sea World trainers.
--   Anna Krystalli, Ben Marwick, Karthik Ram, Nicholas Tierney, and other members of the open science R community for developing tools and educational resources that facilitate open science practices.
+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.
 
 # Footnotes
 

analysis/paper/ballistocardiogram_v4.docx

@@ -551,3 +551,27 @@ @article{wright2014
 	url = {http://dx.doi.org/10.3354/meps10592},
 	langid = {en}
 }
+
+@article{williams2021,
+	title = {Field Physiology: Studying Organismal Function in the Natural Environment},
+	author = {{Williams}, {Cassondra L.} and {Hindle}, {Allyson G.}},
+	year = {2021},
+	month = {06},
+	date = {2021-06-30},
+	journal = {Comprehensive Physiology},
+	pages = {1979--2015},
+	doi = {10.1002/cphy.c200005},
+	url = {http://dx.doi.org/10.1002/cphy.c200005}
+}
+
+@article{andreozzi2021,
+	title = {A Novel Broadband Forcecardiography Sensor for Simultaneous Monitoring of Respiration, Infrasonic Cardiac Vibrations and Heart Sounds},
+	author = {{Andreozzi}, {Emilio} and {Gargiulo}, {Gaetano D.} and {Esposito}, {Daniele} and {Bifulco}, {Paolo}},
+	year = {2021},
+	month = {11},
+	date = {2021-11-18},
+	journal = {Frontiers in Physiology},
+	volume = {12},
+	doi = {10.3389/fphys.2021.725716},
+	url = {http://dx.doi.org/10.3389/fphys.2021.725716}
+}