Shulgach JA, Beam DW, Nanivadekar AC, Miller DM, Fulton S, Sciullo M, Ogren J, Wong L, McLaughlin BL, Yates BJ, Horn CC, Fisher LE. Selective stimulation of the ferret abdominal vagus nerve with multi-contact nerve cuff electrodes. Sci Rep. 2021 Jun 21;11(1):12925. doi: 10.1038/s41598-021-91900-1. PMID: 34155231; PMCID: PMC8217223.
This is the repository for analysis code with selective stimulation. Loading and analyzing the data can be performed together or separately. Code folder uses a file organizing structure inspired by BIDS (Brain Imaging Data Structure).
- At least Matlab 2019b
- Box Drive (optional to access raw data)
- Stastistics and Machine Learning Toolbox
- Signal Processing Toolbox
Download code as .ZIP file and extract to desired directory, or open a terminal and use git as illustrated below:
git clone https://github.com/cchorn/SPARC_vagus_stim_optimization
Open Matlab and add folder to matlab path (adds support functions to matlab internal file searchpath) by inserting the following into the Matlab command window:
cd SPARC_vagus_stim_optimization
addpath(pwd)
savepath
- Online - Scripts used during experiments
- Processing Pipeline - Primary files for data analysis using Matlab
- Supporting Functions - Helper functions for figures, data indexing, and wrappers
- analysis_sctipt.m - Data analysis pipeline with analysis steps similar to readme file.
- expmt_list.mat - Header file containing animal metadata
The "expmt_list" header contains all of the metadata associated with each animal including the experiment day, number of stim trials, stim amplitudes, pulse widths, etc. It keeps track of the last stim analysis run with channnel stim response detections across all trials contained in the field last_run, which can be used to compare header file versions if running the analysis yourself.
The following header fields will be useful to the user:
cohorttrial_listcuff_listpulseWidthstim_histselectivityGridminthreshstd_gaincv
The name of the animal (as a note, normally cohort refers to a group...)
An array of up to 6 values representing the pulse width used while performing a binary search for threshold stim. Each pulse width can be referred to as a session.
A 6x2 double with each row as a pulse width session, column 1 as starting trial number and 2nd column as ending trial number
A 6x3 char array with the cuff electrode contact pairs used in each pulse width session
An 6xN array of numbers representing the stimulation amplitude used for a particular trial and session. N is the length of the longest session with the most trials. Each row covers the stimulation amplitudes used for one session, with the first column as the starting stim value. To read animal stim history, as an example using F21-19, the first pulse width tested was 100uA using trial 2 (in trial_list, first row == first session).
A 6xNx32 array of binary values containing response information across all trials, with 6 sessions per animal, and 32 channels from the MEA. N is the length of the longest session with the most trials.
A 32x6 array of values where each column is a session and each row represents one of the 32 MEA channels, numbered from 1-32 read top-down (easy referencing). If a NaN is present in one of the columns then no response was detected throughout that session.
The gain value multiplied with the standard deviation of the baseline data. The baseline data is collected from all trials, however ROC curves were generated by running the analysis with varying gain values adjusting the response detection threshold line to determine the gain value that resulted in a 90% True Positive rate.
A struct containing the response "conduction velocity" (cv) metadata with each channel and trial, subpoppulated with vectors of binary data indicating detected responses within 0.5m/s intervals and time indices when those cv intervals begin/end. Within a 500ms detection period, from right to left, the cv will increase, and the space between time indices for cv steps exponentially decreases.
Analyzing the raw data for stim responses, performing classification, or parameter optimization for selectivity pipelines can be done separately. Running the stim analysis pipeline will overwrite the channel response detection history in the expmt_list header file, so make sure to save a backup before running the analyze_stim_data_finalfunction or if you want to see the previous history.
Animals are analyzed in the chronological order of the experiments - 1 : F20-19, 2 : F21-19, 3 : F19-19, 4: F25-19, 6 : F26-19, 7 : F34-19
Load "expmt_list.mat" header variable into matlab workspace, set "expmt" to number from 2 to 7 (referring to animal you want to process), and set data directory (source folder) where the data is located.
Note: Because animal 1 (F20-19) has only 2PW tested with one cuff pair, analysis on this animal is ignored, and performed on animals 2-7 Ex: expmt = 2 : (F21-19)
load(expmt_list.mat);
expmt = 2;
Having box installed on your desktop is preffered. Adjust the directory as necessary
source_path = 'C:\Users\{PC_username}\box\path\to\raw\data';
If you want to view the response data, skip to step 6.
To establish baseline, each trial and channel is looked at to inspect noise levels.
baseline = visualizeNoise(source_path, expmt_list, expmt);
create_baseline_figure(expmt_list, expmt, baseline); % to visualize baseline across trials
Inspect figure for any trials with large noise. Check if any noise levels cross a set threshold. This is only needed if changes to the analysis pipeline occur, otherwise all files have been manually checked beforehand.
thresh = 1; %stdev
[good_trials, bad_trials, trial_results] = analyze_baseline_values(baseline, thresh, expmt);
(Optional) remove any noisy trials or channels that are NOT CURRENTLY in the header file.
manual_trial_skip = []; % add in any abnormal trials here
manual_chan_skip = []; % add in any abormal channels here
expmt_list{expmt,1}.exclude_trials = sort([excluded_trials, bad_trials, manual_trial_skip]);
expmt_list{expmt,1}.exclude_channels = manual_chan_skip
Load data into workspace for reading files from F21-19, F22-19,..., or F34-19. Must read from only one animal at a time. Repeat for new animal if desired
data_table = segment_raw_data(source_path, expmt_list, expmt);
[new_signal_data, avg_signal_data] = convert_seg_to_table(expmt_list, expmt, data_table);
Responses are detected using the threshold-crossing RMS of the averaged segmented channel snippets, taken from time indices of stimulation events. Threshold is obtained by baseline mean + chan_gain x baseline std for each channel.
expmt_list = analyze_stim_data_final(expmt_list, baseline, avg_signal_data, expmt);
Count all responses and compare MEA performance across all PW combinations between cuff pairs. The Selective Index function takes the stimulation response data across trials stored in the expmt_list header and, for each animal, collects metadata from shared or unique channel responses as well as calculates the selectivity index with the following function: SI = a*((Total_Response - Total_Shared)/32) + b*(abs(Total_A_Response - Total_B_Response)/Total_Response). "a" and "b" are coefficients corresponding to providing weight distribution between efficiency and balance. (a + b = 1)
[SI_combos_list, overlap_combos_list] = SI_search(expmt_list, 2:7)
Generate plots to show the optimal stimulation parameter settings to get the most "selective" responses, using the "Selectivity Index" function output. e.g, PW####.max_SI.max_SI_selec_chan contains 90% selectivity results (Note: May need to set patrameters inside function to plot desired figures)
[avg_line_data, expmt_list] = get_selec_stim_bound(SI_combos_list, expmt_list);
Generate plots to show the optimal stimulation parameter settings to get the most "selective" responses, using the "Selectivity Index" function output.
opt_params_list = get_opt_params(SI_combos_list, expmt_list, 2:7);
Generate cuff threshold RC curves with x-axis as PW duration and y-axis as stim amplitude
min_thresh_plot(expmt_list)
Generate cv figures showing response behavior at max stim and cuff threshold
[cv_min_data, cv_max_data] = get_conduction_velocity(expmt_list, 2:7, true);
The nearest active MEA channels are collected running the script nearest_neighbor.m. The centroid is calculated for each animal from the MEA response running the script Centroids.m. For a specific stim and pulse width, manually set the index variables for either max stim or threshold like below. (Note: Threshold settings are referencing the header file and position of values within stim_hist and pulseWidth fields)
expmt = 2; % choose animal (2 to 7)
PW_idx_1 = 1; % PW session for cuff pair 1:2 (choose from 1 to 3)
stim_idx_1 = 4; % stimulation amplitudes (select from PW session)
PW_idx_2 = 4; % PW session for cuff pair 3:4 (choose from 4 to 6)
stim_idx_2 = 4;
To generate violin plots, run the script Violin_Plotting.m after collecting the centroid data. The variables PW100. PW500, and PW1000 contain the centroid locations and are included in the script to save time.
Generate plots to show selective boundaries with x-axis and y-axis as stim amplitude for each cuff pair. (Note: change parameters inside function to choose different classifier methods)
[avg_line_data, expmt_list] = get_selec_stim_bound(SI_combos_list, expmt_list);
This work was supported by funding from the NIH Common Fund SPARC Program (award#U18TR002205; Closed-loop neuroelectric control of emesis and gastric motility)
SPARC, Stimulating Peripheral activity to relieve conditions
