Bridging between ByCycle and NeuroDSP functionalities

[davidglitvin]

Dear ByCycle team,

I am writing to inquire about how to bridge/transition between ByCycle and NeuroDSP. More specifically, I have an array with multiple sub-arrays corresponding to different experimental subjects that I pass through the cycle-by-cycle algorithm, which I store in the "BycycleGroup" object. As a quick background -

• Here, a "BycycleGroup" object is initialized with specific thresholds and parameters. The fit method is then called on the object, passing the signal data (sigs), the sampling frequency (fs), and a frequency range (f_alpha).

bg_dict[key] = BycycleGroup(thresholds=thresholds, center_extrema='peak')
bg_dict[key].fit(sigs, fs, f_alpha, axis=0)

• I can then plot a specific sub-array and visualize cycle detection using the code-block below (see attached image for example of output).

# Specify the key for the sub-array you want to plot
sub_array_key = 'KO_m13_week_12_ch1_filtered_2-25Hz'

# Ensure the sub_array_key exists in bg_dict
if sub_array_key not in bg_dict:
    raise ValueError(f"Key {sub_array_key} not found in bg_dict")

# Get the BycycleGroup object for the specified sub-array
bg1 = bg_dict[sub_array_key]

# Number of channels in the selected BycycleGroup object
num_channels = len(bg1.df_features)
ch_names = [f"Channel {i+1}" for i in range(num_channels)]  # Generate channel names for plotting

# Plot the cycle analysis for each channel
for idx in tqdm(range(num_channels), desc='Processing channels'):
    bg1[idx].plot(xlim=(0, 5), figsize=(16, 3))
    plt.title(ch_names[idx])
    plt.show()
    print(f"Processed channel {ch_names[idx]}")

• My questions are as follows and relate to the properties of the "ByCycleGroup" object :

1. How do I switch to a burst visualization that includes only the black and red (overlaid) traces used in NeuroDSP (see below). Can you pass the ByCycleGroup object into something like the "plot_bursts" function in NeuroDSP (see below) and get just a red/black plot and then extract burst statistics? Or is there a way to extract the black and red traces directly from the "ByCycleGroup" object. Once I check the quality of the peak detection, I don't need to see the amp fraction, amp consistency, period consistency, and monotonicity, nor their corresponding rectangular colored highlighting.

2. How do I access segments that are, or are-not bursts from the "ByCycleGroup" object and compute lagged coherence on these segments. Do I need to pass the "ByCycleGroup" object through "NeuroDSP" - using the example provided in your tutorial (see below) ? Alternatively, if I have a data array containing so much information about every cycle, (i.e. time_peak, time_trough, volt_peak, volt_trough, time_decay, time_rise, volt_decay, volt_rise, volt_amp, time_rdsym, time_ptsym), can lagged coherence (or even window matching) be applied to the array instead?

• To be clear, I already can access the "ByCycleGroup" array, which has the "is_burst" column, and from which I can take all consecutive indices that are "TRUE" - and their corresponding time values. This allows me to highlight and extract the beginning and end of bursts from the filtered timeseries (used as input to ByCycle). Afterwards, I can apply lagged coherence analysis on just the epoch containing the burst, but this seems like a sloppy work-around.

[ryanhammonds]

Hi @davidglitvin, sorry for the delayed response. Lets start with some simulated data to answer your questions:

import numpy as np
from bycycle import BycycleGroup
from neurodsp.sim import sim_powerlaw
from neurodsp.plts import plot_bursts

# Simulate
n_seconds = 10
fs = 1000

osc = np.zeros(int(n_seconds * fs))
osc[4000:6000] = np.sin(2*np.pi*10*np.arange(0, 2000)/fs + np.pi)

n_rows = 2
sigs = np.zeros((n_rows, int(n_seconds * fs)))            
for i in range(n_rows):
    sigs[i] = sim_powerlaw(n_seconds, fs) + osc

# Fit Bycyle
bg = BycycleGroup(thresholds={})
bg.fit(sigs, fs, (5, 15))
bg.recompute_edges()

1) Below plots bursts similar to neurodsp. Unfortunately, the is_burst array below isn't easy accessible and has to be created. This could added more conveniently added with a property decorator on a is_burst method on the base or 1d class. Or even just a utility function that computes the is_osc variable below. PRs are welcome!

# Plot
sig = sigs[0]
df_features = bg[0].df_features

times = np.arange(0, len(sig)/fs, 1/fs)

# Determine which samples are defined as bursting
is_osc = np.zeros(len(sig), dtype=bool)
df_osc = df_features.loc[df_features['is_burst']]
start = 0

for cyc in df_osc.to_dict('records'):

    samp_start_burst = int(cyc['sample_last_trough']) - int(fs * start)
    samp_end_burst = int(cyc['sample_next_trough'] + 1) - int(fs * start)

    is_osc[samp_start_burst:samp_end_burst] = True

plot_bursts(times, sig, is_osc, labels=['Signal', 'Bursts'])

Alternatively, you could use this since it's built in, but it doesn't produce the same plots as plot_bursts:

bg[0].plot(plot_only_results=True)

2) For lagged coherence of non-bursts, I would create a chunk variable, this variable is defined in lagged_coherence_1freq, see neurodsp's source. I would copy these functions since there isn't an interface for user-defined chunk. Putative cycles are defined from bg[sig_ind].df_features.iloc[i_cyc_start].sample_last_trough to bg[sig_ind].df_features.iloc[i_cyc_start].sample_next_trough. These are associated with the is_burst column, which should be False for this.

# Determine these using the `is_burst` column 
i_cyc_start = ... # the first row is non-burst in the segment
i_cyc_end = ... # last row that is non-burst in the same segment

# Where non-burst segment starts and ends
sig_ind = 0
i_start = bg[sig_ind].df_features.iloc[i_cyc_start].sample_last_trough
i_end = bg[sig_ind].df_features.iloc[i_cyc_end].sample_next_trough

# Slice
sig_non_burst = sigs[sig_ind, i_start:i_end]

I would then chunk, e.g. reshape, sig_non_burst above based on the length of desired lagged coherence window. You'll probably need to slice off a few samples so it's evenly divisible by the desired window length. The code above is for one non-burst segment. The length of each chunk should be constant, so the end result will be a loop of the above and a single 2d array of non-burst chunks. Then continue from this line.

[davidglitvin]

Ok, thanks for the help! I had another question that is somewhat related, and probably doesn't warrant opening a new issue for this github repository. For the FOOOFGroup array containing information about the cycles that are (or are not) parts of bursts, there are a number of features that describe each cycles' shape (i.e. amp_fraction, amp_consistency, period_consistency, monotonicity, period, time_peak, time_trough, time_decay, time_rise, time_rdsym, time_ptsym, etc). I am noticing that in my data the bursts have several classes of cycle groups with common values for the attributes. Have you tried using something like ICA, or machine learning to classify discrete clusters of cycles within bursts?

[ryanhammonds]

You're welcome!

I've used kmeans in the past to cluster cycles. I did this by resampling cycles to have the same number of samples per cycle. Then passed the resampled cycles into kmeans. Resampling is probably not ideal unless each cycle has a sufficient number of samples and there is little variation in the number of samples per cycle.

Using bycycle features + a clustering method sounds like a good idea and would be a nice use case! Let me know how it works out. If the data is open source or simulate-able, then it may be a nice addition to the examples page.