Skip to content

commonality

Description

This module contains the functions needed to estimate common drive.


pooled_intramuscular_coherence(emgfile, number_of_mus_cst=-1, number_iterations=100, window_type='hanning', window_duration_seconds=1, overlap_seconds=0.5, nfft=None, random_state=42)

Estimate common synaptic input oscillations using coherence analysis between two cumulative spike trains.

The procedure follows these steps:

  1. Randomly split MUs into two groups.
  2. Sum the binary spike trains within each group to obtain two cumulative spike trains (CSTs).
  3. Estimate auto-spectra of CST 1 and CST 2.
  4. Estimate cross-spectrum between CST 1 and CST 2.
  5. Repeat over random permutations (default: 100).
  6. Pool spectra across iterations.
  7. Compute coherence from the pooled spectra.

References:

PARAMETER DESCRIPTION
emgfile

The dictionary containing the emgfile.

TYPE: dict

number_of_mus_cst

Number of MUs included in each group for the cumulative spike train. If -1, defaults to floor(number_of_mus / 2).

TYPE: int DEFAULT: -1

number_iterations

Number of random permutations.

TYPE: int DEFAULT: 100

window_type

Type of window used for spectral estimation. Can be:

hanning A Hann window.

hamming A Hamming window

TYPE: str {"hanning", "hamming"} DEFAULT: "hanning"

window_duration_seconds

Size of the window for coherence estimation in seconds.

TYPE: float DEFAULT: 1

overlap_seconds

Number of seconds of overlap between adjoining segments. Considering a window_duration_seconds of 1 second, an overlap of 0.5 seconds means a 50% overlap.

TYPE: float DEFAULT: 0.5

nfft

Number of FFT points. If None, defaults to 10 * fsamp.

TYPE: int or None DEFAULT: None

random_state

Random seed for reproducible random permutations.

TYPE: int DEFAULT: 42

RETURNS DESCRIPTION
coherence_results

Dictionary containing the following keys:

  • "coherence_all_pairs": list of numpy arrays containing the coherence signals for all CST pairs (list of numpy arrays).
  • "coherence": the coherence curve (numpy array).
  • "frequency": the frequency values corresponding to the coherence curve (numpy array).
  • "window_duration_seconds": the window duration in seconds (float).
  • "window": the window used for spectral estimation (numpy array).
  • "overlap_seconds": the overlap between segments in seconds (float).
  • "cst_duration_seconds": the duration of the cumulative spike trains in seconds (float).
  • "fsamp": the sampling frequency in Hz (float).

TYPE: dict

Examples:

Estimate pooled intramuscular coherence with overlapping windows.

>>> import openhdemg.library as emg
>>> emgfile = emg.emg_from_samplefile()
>>> coherence_results = emg.pooled_intramuscular_coherence(
...     emgfile=emgfile,
...     number_iterations=100,
...     window_duration_seconds=1,
...     overlap_seconds=0.5,
... )
>>> coherence_results["coherence"]

Estimate coherence without overlapping windows.

>>> coherence_results = emg.pooled_intramuscular_coherence(
...     emgfile=emgfile,
...     number_iterations=100,
...     window_duration_seconds=1,
...     overlap_seconds=0,
... )
>>> coherence_results["coherence"]


z_score_coherence(coherence_results)

Z-score coherence following two possible cases:

  1. No overlap:

    z = sqrt(2 * L) * atanh(sqrt(coherence))

    where L is the number of non-overlapping time segments used in the Welch method.

  2. With overlap:

    z = sqrt(2 * K_tilde) * atanh(sqrt(coherence))

    where K_tilde is the corrected number of independent segments accounting for overlap between Welch segments.

PARAMETER DESCRIPTION
coherence_results

Dictionary returned by pooled_intramuscular_coherence.

TYPE: dict

RETURNS DESCRIPTION
z_score_results

Dictionary containing the same keys of pooled_intramuscular_coherence plus the following keys:

  • "z_score": the z-scored coherence curve (numpy array)
  • "factor": the factor used for z-scoring (float)

TYPE: dict

Examples:

Z-score the coherence curve returned by pooled_intramuscular_coherence().

>>> import openhdemg.library as emg
>>> emgfile = emg.emg_from_samplefile()
>>> coherence_results = emg.pooled_intramuscular_coherence(
...     emgfile=emgfile,
...     overlap_seconds=0.5,
... )
>>> z_score_results = emg.z_score_coherence(coherence_results)
>>> z_score_results["z_score"]


calculate_coherence_bands(coherence_results, confidence_level_type='high_frequency_bias', bias_frequency_range=[250, 500], delta_frequency_range=[1, 5], alpha_frequency_range=[5, 15], beta_frequency_range=[15, 35], gamma_frequency_range=[35, 60])

Calculate band-specific parameters from either a z-scored coherence curve or a raw coherence curve.

If coherence_results contains "z_score", the function uses it. Otherwise, the function uses raw coherence curve.

For each frequency band, this function calculates:

  1. Average coherence/z-score above the confidence level.
  2. Area under the coherence/z-score curve above the confidence level.
PARAMETER DESCRIPTION
coherence_results

Dictionary containing either keys {"z_score", "frequency"} or {"coherence", "frequency"}.

TYPE: dict

confidence_level_type

Type of confidence level used to determine which parts of the coherence/z-score curve are considered significant.

high_frequency_bias The confidence level is defined as the average coherence/z-score in the high-frequency bias band defined by bias_frequency_range.

theoretical The confidence level is defined as 95% confidence level based on the theoretical distribution.

- For z-score: 1.65 (95th percentile of the standard normal
distribution).

- For coherence: 1 - 0.05 ** (1 / (L - 1)), where L is the number
of non-overlapping segments used in the Welch method.

TYPE: str {"high_frequency_bias", "theoretical"} DEFAULT: "high_frequency_bias"

bias_frequency_range

Frequency range for the bias band in Hz when confidence_level_type is "high_frequency_bias".

TYPE: list DEFAULT: [250, 500]

delta_frequency_range

Frequency range for the delta band in Hz.

TYPE: list DEFAULT: [1, 5]

alpha_frequency_range

Frequency range for the alpha band in Hz.

TYPE: list DEFAULT: [5, 15]

beta_frequency_range

Frequency range for the beta band in Hz.

TYPE: list DEFAULT: [15, 35]

gamma_frequency_range

Frequency range for the gamma band in Hz.

TYPE: list DEFAULT: [35, 60]

RETURNS DESCRIPTION
coherence_band_results

Dataframe with the following columns:

  • "coherence_type": "z_score" or "coherence".
  • "confidence_level": confidence level used (float).
  • "average_delta": average coherence/z-score in the delta band above the confidence level.
  • "auc_delta": area under the coherence/z-score curve in the delta band above the confidence level.
  • "average_alpha": average coherence/z-score in the alpha band above the confidence level.
  • "auc_alpha": area under the coherence/z-score curve in the alpha band above the confidence level.
  • "average_beta": average coherence/z-score in the beta band above the confidence level.
  • "auc_beta": area under the coherence/z-score curve in the beta band above the confidence level.
  • "average_gamma": average coherence/z-score in the gamma band above the confidence level.
  • "auc_gamma": area under the coherence/z-score curve in the gamma band above the confidence level.

TYPE: DataFrame

Examples:

Calculate frequency-band summaries from a z-scored coherence curve.

>>> import openhdemg.library as emg
>>> emgfile = emg.emg_from_samplefile()
>>> coherence_results = emg.pooled_intramuscular_coherence(
...     emgfile=emgfile,
...     overlap_seconds=0.5,
... )
>>> z_score_results = emg.z_score_coherence(coherence_results)
>>> coherence_band_results = emg.calculate_coherence_bands(
...     coherence_results=z_score_results,
...     confidence_level_type="theoretical",
... )
>>> coherence_band_results


smooth_spiketrains_convolution(binary_mus_firing, fsamp, window_type='hanning', window_duration_seconds=0.4)

Smooth a binary spike train by convolving it with a window.

PARAMETER DESCRIPTION
binary_mus_firing

The binary MU firing data.

TYPE: ndarray

window_type

Type of window used for spectral estimation. Can be:

hanning A Hann window.

hamming A Hamming window

TYPE: str {"hanning", "hamming"} DEFAULT: "hanning"

window_duration_seconds

Smoothing window duration in seconds.

TYPE: float DEFAULT: 0.4

RETURNS DESCRIPTION
smoothed_mus_firing

Smoothed MU discharge rates obtained by convolving the binary spike train with the specified window. Arrangement of the output array is the same as binary_mus_firing (samples x MUs).

TYPE: ndarray

Examples:

Smooth the binary spike trains stored in an EMG file.

>>> import openhdemg.library as emg
>>> emgfile = emg.emg_from_samplefile()
>>> binary_mus_firing = emgfile["BINARY_MUS_FIRING"].to_numpy()
>>> smoothed_mus_firing = emg.smooth_spiketrains_convolution(
...     binary_mus_firing=binary_mus_firing,
...     fsamp=emgfile["FSAMP"],
... )
>>> smoothed_mus_firing


smoothed_dr_pca(emgfile, window_type='hanning', window_duration_seconds=0.4, filter_highcut=0.75, filter_order=3, remove_edge_seconds=1, method_n_components='parallel_analysis', number_iterations_pa=1000, percentile_pa=95, variance_threshold=80.0, random_state=42)

Estimate low-dimensional components from MU spike trains using PCA applied to smoothed discharge-rate profiles.

The procedure follows this pipeline:

  1. Convert binary MU spike trains to smoothed discharge-rate profiles.
  2. High-pass filter the smoothed discharge rates to remove offsets and slow trends.
  3. Standardize each MU discharge profile.
  4. Calculate the covariance/correlation matrix.
  5. Extract eigenvalues and eigenvectors.
  6. Use parallel analysis to estimate the number of components to retain.
  7. Project the standardized discharge profiles onto the retained eigenvectors.

References:

PARAMETER DESCRIPTION
emgfile

The dictionary containing the emgfile.

TYPE: dict

window_duration_seconds

Duration of the Hann window used to smooth binary spike trains.

TYPE: float DEFAULT: 0.4

filter_highcut

High-pass cutoff frequency used to remove offsets and slow trends.

TYPE: float DEFAULT: 0.75

filter_order

Butterworth filter order.

TYPE: int DEFAULT: 3

remove_edge_seconds

Number of seconds to remove from the beginning and end of the smoothed discharge profiles to avoid edge artifacts from convolution and filtering.

TYPE: float DEFAULT: 1

method_n_components

Method used to estimate the number of components to retain.

parallel_analysis Retain the number of components that present eigenvalues greater than the percentile (e.g., 95th) of eigenvalues generated from randomly shuffling the data.

variance_greater_than_threshold Retain the number of components that explain more than a certain percentage (e.g., 80%) of the variance.

eigenvalue_greater_than_one (Kaiser's criterion). Retain the number of components with eigenvalues greater than 1.

TYPE: str DEFAULT: "parallel_analysis"

number_iterations_pa

Number of iterations for parallel analysis.

TYPE: int DEFAULT: 1000

percentile_pa

Percentile used in parallel analysis.

TYPE: float DEFAULT: 95

variance_threshold

Threshold if "variance_greater_than_threshold" is used.

TYPE: float DEFAULT: 80.0

random_state

Seed for reproducibility.

TYPE: int DEFAULT: 42

RETURNS DESCRIPTION
pca_results

Dictionary containing the following keys:

  • "smoothed_mus_firing": the smoothed binary spike trains calculated by window convolution (numpy array) smoothed_mus_firing has the same shape as the input binary_mus_firing (samples x MUs).
  • "kmo": the Kaiser-Meyer-Olkin measure of sampling adequacy (float).
  • "method_n_components": the selected method to estimate the number of components to retain (str).
  • "eigenvalues_threshold_pa": the eigenvalue threshold based on parallel analysis (float).
  • "variance_threshold": the variance threshold used if method_n_components is "variance_greater_than_threshold" (float).
  • "number_of_components_retained": the number of components retained based on parallel analysis (int).
  • "eigenvalues": the eigenvalues of the covariance/correlation matrix (numpy array).
  • "eigenvectors": the eigenvectors of the covariance/correlation matrix (numpy array).
  • "explained_variance": the explained variance in percentage of each component (numpy array).
  • "cumulative_explained_variance": the cumulative explained variance in percentage (numpy array).
  • "low_dimensional_components": the low-dimensional components (numpy array).

TYPE: dict

Notes
  1. Kaiser-Meyer-Olkin (KMO) measure is calculated to assess the factorability of the dataset before performing PCA. Typically, a KMO value above 0.7 is considered acceptable (Hoelzle & Meyer, 2012; https://doi.org/10.1002/9781118133880.hop202006).
  2. Altough other methods to define the number of components are available, we reccomend using parallel analysis as it is more robust and data-driven than arbitrary thresholds (e.g., eigenvalue > 1 or variance explained > 80%).

Examples:

Estimate PCA components from smoothed MU discharge-rate profiles.

>>> import openhdemg.library as emg
>>> emgfile = emg.emg_from_samplefile()
>>> pca_results = emg.smoothed_dr_pca(
...     emgfile=emgfile,
...     method_n_components="parallel_analysis",
... )
>>> pca_results["number_of_components_retained"]
>>> pca_results["explained_variance"]


common_drive_index(emgfile, window_type='hanning', window_duration_seconds=0.4, filter_highcut=0.75, filter_order=3, remove_edge_seconds=1.0, window_corr_seconds=5.0, overlap_corr_seconds=0.0, max_lag_seconds=0.1, number_of_surrogates=2, confidence_percentile=95, random_state=42)

Calculate the common drive index as the average pairwise cross-correlation between smoothed MU discharge profiles.

Cross-correlation is calculated in windows of window_corr_seconds. For each MU pair, the maximum absolute correlation within ± max_lag_seconds is extracted for each segment and then averaged across segments.

References:

PARAMETER DESCRIPTION
emgfile

The dictionary containing the emgfile.

TYPE: dict

window_duration_seconds

Duration of the Hann window used to smooth binary spike trains.

TYPE: float DEFAULT: 0.4

filter_highcut

High-pass cutoff frequency used to remove offsets and slow trends.

TYPE: float DEFAULT: 0.75

filter_order

Butterworth filter order.

TYPE: int DEFAULT: 3

remove_edge_seconds

Number of seconds to remove from the beginning and end of the smoothed discharge profiles to avoid edge artifacts from convolution and filtering.

TYPE: float DEFAULT: 1

window_corr_seconds

Segment duration in seconds for pairwise cross-correlation. If -1, the entire signal is used.

TYPE: float DEFAULT: 5.0

overlap_corr_seconds

Overlap duration in seconds for pairwise cross-correlation. Segments will overlap by overlap_corr_seconds. 0.0 means no overlap.

TYPE: float DEFAULT: 0.0

max_lag_seconds

Maximum lag for cross-correlation, in seconds.

TYPE: float DEFAULT: 0.1

number_of_surrogates

Number of surrogate spike-train sets used to estimate the confidence level. Surrogates are generated by independently shuffling the interspike intervals of each MU.

TYPE: int DEFAULT: 2

confidence_percentile

Percentile of surrogate peak correlations used as the confidence level to determine whether each pairwise correlation is significant.

TYPE: float DEFAULT: 95

random_state

Random seed for reproducible surrogate generation.

TYPE: int DEFAULT: 42

RETURNS DESCRIPTION
cdi_results

Dictionary containing pairwise correlations and common drive index with keys:

  • "smoothed_mus_firing": the smoothed binary spike trains calculated by window convolution (numpy array). smoothed_mus_firing has the same shape as the input binary_mus_firing (samples x motor units).
  • "pairwise_correlation": dataframe of average pairwise correlation and if it is significant for each MU pair (dataframe).
  • "pairwise_correlation_matrix": square matrix of average pairwise correlations between all MU pairs (numpy array). pairwise_correlation_matrix[i,j] contains the average pairwise correlation between MU i and j.
  • "pairwise_correlation_matrix_thresholded": thresholded pairwise_correlation_matrix (numpy array). correlation values that are not significant are set to zero.
  • "common_drive_index_mean": mean of cross-correlation values across all MU pairs (float).
  • "common_drive_index_std": standard deviation of cross-correlation values across all MU pairs (float).
  • "common_drive_index_thresholded_mean": mean of cross-correlation values across all MU pairs that are above the significance threshold (float).
  • "common_drive_index_thresholded_std": standard deviation of cross-correlation values across all MU pairs that are above the significance threshold (float).
  • "confidence_level": the confidence level used to determine if a pairwise correlation is significant (float).
  • "percentage_significant_pairs": percentage of MU pairs with significant correlation (float).
  • "corr_signals_all_pairs": list of numpy arrays containing the cross-correlation signals for all pairs (list of numpy arrays).
  • "grand_average_corr_values": grand average of the corr_signals_all_pairs (numpy array).
  • "lags_seconds": lags used to calculate the cross-correlation signals (numpy array).

TYPE: dict

Examples:

Calculate the common drive index using non-overlapping 5-second correlation windows.

>>> import openhdemg.library as emg
>>> emgfile = emg.emg_from_samplefile()
>>> cdi_results = emg.common_drive_index(
...     emgfile=emgfile,
...     window_corr_seconds=5,
...     overlap_corr_seconds=0,
... )
>>> cdi_results["pairwise_correlation"]
>>> cdi_results["common_drive_index_mean"]


pci_index(emgfile, number_iterations=100, window_type='hanning', window_duration_seconds=1.0, overlap_seconds=0.0, nfft=None, delta_frequency_range=[1, 5], alpha_frequency_range=[5, 15], beta_frequency_range=[15, 35], random_state=42)

Estimate the proportion of common input index (PCI) from MU spike trains.

This implementation follows the logic of Negro et al. (2016), where coherence is calculated between two CSTs containing n MUs each. The number of MUs per CST is varied from 1 to floor (total number of MUs / 2).

The relationship between average coherence and number of units per CST is fitted using Eq. 4:

C(n) = ((n^2 * A)^2) / ((n * B + n^2 * A)^2)

PCI is then calculated as:

PCI = sqrt(A / B)

References:

PARAMETER DESCRIPTION
emgfile

The dictionary containing the emgfile.

TYPE: dict

number_iterations

Number of random permutations. If the number of unique splits is smaller than number_iterations, all unique splits are used.

TYPE: int DEFAULT: 100

window_type

Type of window used for spectral estimation. Can be:

hanning A Hann window.

hamming A Hamming window

TYPE: str {"hanning", "hamming"} DEFAULT: "hanning"

window_duration_seconds

Size of the window for coherence estimation in seconds.

TYPE: float DEFAULT: 1

overlap_seconds

Number of seconds of overlap between adjoining segments.

TYPE: float DEFAULT: 0

nfft

Number of FFT points. If None, defaults to 10 * fsamp.

TYPE: int or None DEFAULT: None

delta_frequency_range

Delta frequency range in Hz.

TYPE: list DEFAULT: [1, 5]

alpha_frequency_range

Alpha frequency range in Hz.

TYPE: list DEFAULT: [5, 15]

beta_frequency_range

Beta frequency range in Hz.

TYPE: list DEFAULT: [15, 35]

random_state

Random seed for reproducible random permutations.

TYPE: int DEFAULT: 42

RETURNS DESCRIPTION
pci_results

Dictionary containing coherence-vs-CST-size values and fitted PCI values. Keys are:

  • "average_coherence_df": dataframe containing the number of MUs per CST and the corresponding average coherence values for each frequency band (dataframe).
  • "fitted_average_coherence_df": dataframe containing the number of MUs per CST and the corresponding fitted coherence values using the fitted A and B parameters for each frequency band (dataframe).
  • "pci_fit_df": dataframe containing the fitted A and B parameters, and the PCI index for each frequency band (dataframe).

TYPE: dict

Notes

The PCI index is typically implemented only for the delta band (see Negro et al., 2016), but this function also fits the model and estimates the PCI for alpha and beta bands.

Examples:

Estimate the proportion of common input index.

>>> import openhdemg.library as emg
>>> emgfile = emg.emg_from_samplefile()
>>> pci_results = emg.pci_index(
...     emgfile=emgfile,
...     number_iterations=100,
... )
>>> pci_results["average_coherence_df"]


smoothed_dr_mutualinformation(emgfile, window_type='hanning', window_duration_seconds=0.4, filter_highcut=0.75, filter_order=3, remove_edge_seconds=1.0, type_clustering='single')

Estimate non-linear dependencies between MUs using pairwise mutual information.

This method:

  1. Smooths binary MU spike trains.
  2. High-pass filters the smoothed discharge rates.
  3. Copula-normalizes each MU discharge-rate profile.
  4. Estimates pairwise Gaussian mutual information between MUs.
  5. Constructs a symmetric mutual-information network.
  6. Applies modified percolation analysis to identify a network threshold.
  7. Applies link-community detection to identify overlapping modules.

References:

PARAMETER DESCRIPTION
emgfile

The dictionary containing the emgfile.

TYPE: dict

window_duration_seconds

Duration of the Hann window used to smooth binary spike trains.

TYPE: float DEFAULT: 0.4

filter_highcut

High-pass cutoff frequency used to remove offsets and slow trends.

TYPE: float DEFAULT: 0.75

filter_order

Butterworth filter order.

TYPE: int DEFAULT: 3

remove_edge_seconds

Number of seconds to remove from the beginning and end of the smoothed discharge profiles to avoid edge artifacts from convolution and filtering.

TYPE: float DEFAULT: 1

type_clustering

Link-community clustering type.

single

complete

TYPE: str {"single", "complete"} DEFAULT: "single"

RETURNS DESCRIPTION
network_results

Dictionary containing keys:

  • smoothed_mus_firing: the smoothed binary spike trains calculated by window convolution (numpy array). smoothed_mus_firing has the same shape as the input binary_mus_firing (samples x MUs).
  • pairwise_mi: dataframe of average pairwise mutual information and if it is significant for each MU pair (dataframe).
  • pairwise_mi_matrix: square matrix of average pairwise mutual information between all MU pairs (numpy array). pairwise_mi_matrix[i,j] contains the average pairwise mutual information between MU i and j.
  • pairwise_mi_matrix_thresholded : thresholded mutual-information matrix (numpy array). Mutual-information values that are not significant are set to zero.
  • network_density : network density after thresholding (float)
  • confidence_level : the confidence level calculated using the percolation analysis to determine if a mutual-information value is significant (float).
  • percentage_significant_pairs: percentage of MU pairs with significant mutual information (float).
  • module_affiliation_matrix: binary matrix indicating the module affiliation of each MU (numpy array). module_affiliation_matrix[i,j] is 1 if MU j belongs to component (or community) i, and 0 otherwise.
  • number_of_components: number of low-dimensional components (communities) of the mutual-information network (int).
  • percentage_mus_first_component: percentage of MUs that belong to the first component (float).

TYPE: dict

Examples:

Estimate a mutual-information network from smoothed discharge rates.

>>> import openhdemg.library as emg
>>> emgfile = emg.emg_from_samplefile()
>>> network_results = emg.smoothed_dr_mutualinformation(
...     emgfile=emgfile,
... )
>>> network_results["pairwise_mi"]
>>> network_results["module_affiliation_matrix"]