start metric for cardio respiratory synchronization

examples/example_06_cardio_respiratory_synchronization.py

@@ -0,0 +1,160 @@
+'''
+
+Cardio-respiratory synchronization
+==================================
+
+RSA is not the only form of cardio-respiratory coupling.
+Bartsch et al, 2012 (https://doi.org/10.1073/pnas.1204568109) described another form that they called the cardio-respi phase synchronisation (CRPS).
+CRPS leads to clustering of heartbeats at certain phases of the breathing cycle.
+We developed a way of studying such coupling that is presented in this example
+'''
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+from pathlib import Path
+from pprint import pprint
+
+import physio
+
+
+
+##############################################################################
+# 
+# Read data
+# ---------
+#  
+# 
+
+raw_resp = np.load('resp2.npy')
+raw_ecg = np.load('ecg2.npy')
+srate = 1000.
+times = np.arange(raw_resp.size) / srate
+
+
+##############################################################################
+# 
+# Compute respiration cycles and ecg R peaks
+# ------------------------------------------
+#  
+# 
+
+
+ecg, ecg_peaks = physio.compute_ecg(raw_ecg, srate)
+resp, resp_cycles = physio.compute_respiration(raw_resp, srate)
+
+
+##############################################################################
+# 
+# First plot
+# ------------------------------------------
+#  
+# Just a figure to explore the question : 
+#   * Are R peaks clusterized at a certain phase/time of the respiratory cycle ?
+
+
+
+fig, axs = plt.subplots(nrows=2, sharex=True, figsize=(15, 10))
+ax = axs[0]
+ax.plot(times, resp)
+ax.scatter(resp_cycles['inspi_time'], resp[resp_cycles['inspi_index']], color='g')
+ax.scatter(resp_cycles['expi_time'], resp[resp_cycles['expi_index']], color='r')
+for t in ecg_peaks['peak_time']:
+    ax.axvline(t, color='m')
+
+ax = axs[1]
+ax.plot(times, ecg)
+ax.scatter(ecg_peaks['peak_time'], ecg[ecg_peaks['peak_index']], color='m')
+
+ax.set_xlim(0,  40)
+
+
+##############################################################################
+# 
+# ECG peaks are transformed according to their relative position during their corresponding respiratory cycle phase
+# ---------------------------------------
+#  
+# 
+# physio.time_to_cycle() is the key function that transform times of ECG peaks into phases
+# It requires :
+#   * time of ECG R peaks (ecg_peaks['peak_time'].values)
+#   * respiratory cycle times (resp_cycles[['inspi_time', 'expi_time', 'next_inspi_time']].values)
+#   * segment ratios (the phase ratio of inhalation to exhalation transition. Could be computed via resp_cycles['cycle_ratio'].median())
+# It returns the R peaks phases as floats :
+#   * 4.32 means for example that the current R peak occured during the 4th respiratory cycle at 32% of the duration of the current respiratory cycle.
+# Pooled phases can be obtained by the 1 modulo and representend on a raster plot or a phase histogram.
+
+# sphinx_gallery_thumbnail_number = 2
+
+
+inspi_ratio = resp_cycles['cycle_ratio'].median()
+
+cycle_times = resp_cycles[['inspi_time', 'expi_time', 'next_inspi_time']].values
+
+rpeak_phase = physio.time_to_cycle(ecg_peaks['peak_time'].values, cycle_times, segment_ratios=[inspi_ratio])
+
+
+count, bins = np.histogram(rpeak_phase % 1, bins=np.linspace(0, 1, 101))
+
+fig, axs = plt.subplots(nrows=2, sharex=True, figsize=(8, 6))
+ax = axs[0]
+ax.scatter(rpeak_phase % 1, np.floor(rpeak_phase))
+ax.set_ylabel('resp cycle')
+
+ax = axs[1]
+ax.bar(bins[:-1], count, width=bins[1] - bins[0], align='edge')
+ax.set_xlabel('resp phase')
+ax.set_ylabel('count R peaks')
+
+for ax in axs:
+    ax.axvline(0, color='g', label='inspiration', alpha=.6)
+    ax.axvline(inspi_ratio, color='r', label='expiration', alpha=.6)
+    ax.axvline(1, color='g', label='next_inspiration', alpha=.6)
+ax.legend()
+ax.set_xlim(-0.01, 1.01)
+
+##############################################################################
+# 
+# Cross-correlogram between expiration/inspiration times and ECG peak times
+# -------------------------------------------------------------
+# 
+# Another way of exploring preferential clustering of R peaks is to compare their timing 
+# to a given respiratory cycle time like inspiration of expiration time.
+# 
+# The key function is physio.crosscorrelogram() that :
+#   * computes a combinatorial difference between all R peak times and all given respiratory times
+#   * binarizes the obtained time differences according to the provided bins.
+#
+# An histogram can be plotted to present the possible non-uniformity of the distribution,
+# and in such case, the "attraction" of the R peaks by a given respiratory time.
+
+bins = np.linspace(-3, 3, 100)
+
+
+count, bins = physio.crosscorrelogram(ecg_peaks['peak_time'].values, 
+                               resp_cycles['expi_time'].values,
+                               bins=bins)
+
+fig, axs = plt.subplots(nrows=2, sharex=True)
+ax = axs[0]
+ax.bar(bins[:-1], count, align='edge', width=bins[1] - bins[0])
+ax.set_xlabel('time lag')
+ax.set_ylabel('count')
+ax.axvline(0, color='r', label='expiration', alpha=.6)
+ax.legend()
+
+
+ax = axs[1]
+count, bins = physio.crosscorrelogram(ecg_peaks['peak_time'].values, 
+                               resp_cycles['inspi_time'].values,
+                               bins=bins)
+ax.bar(bins[:-1], count, align='edge', width=bins[1] - bins[0])
+ax.set_xlabel('time lag')
+ax.set_ylabel('count')
+ax.axvline(0, color='g', label='inspiration', alpha=.6)
+ax.legend()
+
+
+plt.show()
+
+ax.set_xlim(-0.01, 1.01)

physio/__init__.py

@@ -1,12 +1,13 @@
-from .tools import compute_median_mad, detect_peak, get_empirical_mode
+from .tools import compute_median_mad, detect_peak, get_empirical_mode, crosscorrelogram
 from .preprocess import preprocess, smooth_signal
 from .ecg import (compute_ecg, clean_ecg_peak, compute_ecg_metrics, 
                   compute_instantaneous_rate, compute_hrv_psd)
 from .respiration import (compute_respiration, get_respiration_baseline, detect_respiration_cycles, 
     clean_respiration_cycles, compute_respiration_cycle_features)
-from .cyclic_deformation import deform_traces_to_cycle_template
+from .cyclic_deformation import deform_traces_to_cycle_template, time_to_cycle
 from .rsa import compute_rsa
 
 from .reader import read_one_channel
 from .plotting import plot_cyclic_deformation
-from .parameters import get_respiration_parameters, get_ecg_parameters
+from .parameters import get_respiration_parameters, get_ecg_parameters
+# from .cardio_respiratory_synchronization import *

physio/cyclic_deformation.py

@@ -129,58 +129,75 @@ def deform_traces_to_cycle_template(data, times, cycle_times, points_per_cycle=4
 
 
 
-# def time_to_cycle(times, cycle_times,  inspi_ratio = 0.4):
-#     """
-#     Map absoulut event time to cycle position.
-#     Util for spike, events, trigs...
-
-#     Input:
-#     times: a times vector
-#     cycle_times: N*2 array columns are inspi and expi times. If expi is "nan", corresponding cycle is skipped
+def time_to_cycle(times, cycle_times,  segment_ratios = 0.4):
+    """
+    Map absolut event time to cycle position.
+    Useful for event to respiration cycle histogram
 
-#     Output:
-#     cycles: cycle position for times (same size than times) nan if outside.
+    Parameters
+    ----------
 
-#     Parameters
-#     ----------
+    segment_ratios: None or float or list of float
+        If multi segment deformation then a list of segmetn ratio must provived.
 
-#     Returns
-#     -------
+    Returns
+    -------
 
-#     """
-
-#     n = cycle_times.shape[0]
-#     num_seg_phase = cycle_times.shape[1]
-#     assert num_seg_phase in (1, 2)
+    """
+    if cycle_times.ndim == 1:
+        cycle_times_1d = cycle_times
+        cycle_times = np.zeros((cycle_times_1d.size - 1, 2), dtype=cycle_times_1d.dtype)
+        cycle_times[:, 0] = cycle_times_1d[:-1]
+        cycle_times[:, 1] = cycle_times_1d[1:]
+
+    # check that the end of a cycle is the same the start of the following cycle
+    assert (cycle_times[1:, 0] == cycle_times[:-1, -1]).all(), 'Start and end cycle times do not match'
 
+    num_seg_phase = cycle_times.shape[1] - 1
+
+    if num_seg_phase == 1:
+        assert segment_ratios is None
+        ratios = [0., 1.]
+    else:
+        assert segment_ratios is not None
+        if num_seg_phase == 2 and np.isscalar(segment_ratios):
+            segment_ratios = [segment_ratios]
+        assert len(segment_ratios) == num_seg_phase - 1
+        ratios = [0.] + list(segment_ratios) + [1.]
+
+
+
+
+
+    n = cycle_times.shape[0]
+
+    # num_seg_phase = cycle_times.shape[1]
+    # assert num_seg_phase in (1, 2)
 
-#     cycle_point = np.zeros_like(cycle_times)
-#     if num_seg_phase ==1:
-#         cycle_point[:, 0] = np .arange(n)
-#     elif num_seg_phase ==2:
-#         cycle_point = np.zeros_like(cycle_times)
-#         cycle_point[:, 0] = np .arange(n)
-#         cycle_point[:, 1] = np .arange(n) + inspi_ratio
 
-#     flat_cycle_times = cycle_times.flatten()
-#     flat_cycle_point = cycle_point.flatten()
-#     keep = ~np.isnan(flat_cycle_times)
-#     flat_cycle_times = flat_cycle_times[keep]
-#     flat_cycle_point = flat_cycle_point[keep]
+    cycle_point = np.zeros((cycle_times.shape[0], len(ratios) - 1))
+    for i in range(len(ratios) - 1):
+        cycle_point[:, i] = np .arange(n) + ratios[i]
+
 
-#     interp = scipy.interpolate.interp1d(flat_cycle_times, flat_cycle_point, kind='linear', axis=0, bounds_error=False, fill_value='extrapolate')
+    flat_cycle_times = cycle_times[:, :-1].flatten()
+    flat_cycle_point = cycle_point.flatten()
+    keep = ~np.isnan(flat_cycle_times)
+    flat_cycle_times = flat_cycle_times[keep]
+    flat_cycle_point = flat_cycle_point[keep]
+    interp = scipy.interpolate.interp1d(flat_cycle_times, flat_cycle_point, kind='linear', axis=0, bounds_error=False, fill_value='extrapolate')
 
 
-#     inside = (times>=cycle_times[0,0]) & (times<cycle_times[-1,0])
-#     cycles = np.zeros_like(times) * np.nan
-#     cycles[inside] = interp(times[inside])
+    inside = (times>=cycle_times[0,0]) & (times<cycle_times[-1,0])
+    cycles = np.zeros_like(times) * np.nan
+    cycles[inside] = interp(times[inside])
 
-#     # put nan when some times are in missing cycles
-#     if num_seg_phase == 2:
-#         ind_missing, = np.nonzero(np.isnan(cycle_times[:, 1]))
-#         in_missing = np.in1d(np.floor(cycles), ind_missing.astype(cycles.dtype))
-#         cycles[in_missing] = np.nan
+    # put nan when some times are in missing cycles
+    if num_seg_phase == 2:
+        ind_missing, = np.nonzero(np.isnan(cycle_times[:, 1]))
+        in_missing = np.in1d(np.floor(cycles), ind_missing.astype(cycles.dtype))
+        cycles[in_missing] = np.nan
 
 
-#     return cycles
+    return cycles
 

physio/tools.py

@@ -96,7 +96,16 @@ def get_empirical_mode(traces, nbins=200):
     mode = bins[ind_max]
 
     return mode
-
+
+
+def crosscorrelogram(a, b, bins):
+    """
+    Lazy implementation of crosscorrelogram.
+    """
+    diff = a[:, np.newaxis] - b[np.newaxis, :]
+    count, bins = np.histogram(diff, bins)
+    return count, bins
+