Cortical motor activity modulates respiration and reduces apnoea in neonates.

Code ↔ Paper

15 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 15 matches · 7 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted

[1] § Materials and methods › Breath and apnoea identification ↔ ip_to_interbreath_intervals.m, lines 2–140 · score 0.79 · shallow breathing, identify breaths, standard deviation, IP signals, adaptive, machine
[2] § Materials and methods › Data recordings ↔ functions/align_using_xcorr.m, the whole file · a weak match · score 0.79 · PiNe, ECG epochs, cross correlations, vital signs, temporal, box
[3] § Materials and methods › Statistical analysis of relationships between cortico-respiratory coupling and apnoea ↔ functions/plot_coherence.m, lines 79–121 · score 0.73 · linear models, confounding factors, linear mixed, apnoea rate, coefficient, predictors
[4] § Materials and methods › Phase-amplitude coupling ↔ functions/data_to_wavelet_time_coh.m, the whole file · a weak match · score 0.72 · Hanning window, cross spectrum, standard deviation, wavelet, epochs, median
[5] § Materials and methods › Phase-amplitude coupling ↔ functions/data_to_wavelet_coh.m, the whole file · a weak match · score 0.71 · Hanning window, cross spectrum, standard deviation, wavelet, epochs, median
[6] § Materials and methods › Statistical analysis of relationships between cortico-respiratory coupling and apnoea ↔ run_coherence_age_and_apnoea_rate.m, lines 13–95 · score 0.65 · linear mixed, apnoea rate, Rstudio, analyse, ID, ventilation
[7] § Materials and methods › Surrogate analysis ↔ functions/permutation_testing_cfc.m, the whole file · a weak match · score 0.63 · surrogate spectra, FieldTrip, permutation, FDR, alpha, cross
[8] § Materials and methods › Directionality of PAC ↔ functions/data_to_pac_and_psi.m, the whole file · a weak match · score 0.61 · frequency resolution, phase slope, cross spectrum, coherency, PAC
[9] § Materials and methods › Changes in PAC throughout the respiratory cycle ↔ compute_time_pac.m, lines 23–147 · score 0.59 · surrogate PAC, EEG channels, cycle, theta, delta, modulation
[10] § Materials and methods › Data recordings ↔ ip_to_interbreath_intervals.m, lines 2–140 · score 0.57 · oxygen saturation, heart rate, Breathing, infants, respiration
[11] § Materials and methods › Data recordings ↔ subfunctions/remove_artifactual_longpauses.m, the whole file · a weak match · score 0.56 · oxygen saturation, heart rate, Breathing, infants, respiration
[12] § Materials and methods › Breath and apnoea identification ↔ subfunctions/remove_artifactual_longpauses.m, the whole file · a weak match · score 0.56 · shallow breathing, IP signals, machine, vector, intervals, model
[13] § Materials and methods › Surrogate analysis ↔ compute_time_pac.m, lines 23–147 · score 0.56 · surrogate PAC, EEG channels, permutation, alpha, coherence, amplitude
[14] § Results and discussion › Respiratory and cortical activity exhibit PAC in infants ↔ compute_cortical_amplitude_modulation.m, lines 154–220 · score 0.52 · amplitude modulation, EEG amplitude, locked, spectral, 13 Hz, alpha
[15] § Results and discussion › Cortico-respiratory coupling is strongest during inspiration and may be driven by cortical activity ↔ compute_cortical_amplitude_modulation.m, lines 154–220 · score 0.52 · amplitude modulation, cortical activity, spectral, cycle, EEG, respiratory

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

MATLAB · 1,678 lines · 54 KB · no license · 2 matches

function [ibi,ibi_time,filtered_ip,filtered_ip_time,filtered_ip_fs] = ip_to_interbreath_intervals(ecg,ecg_time,ecg_fs,ip,ip_time,ip_fs,hr,sats,hr_time,sats_time,model_input)
%% This function takes an impedance pneumograph signal as input, and outputs a series of interbreath-intervals, and the cleaned impedance pneumograph signal, with ECG artefact removed.
%
%*** ALL INPUTS SHOULD BE LONGER THAN 5 MINUTES, AND BE IN SECONDS***
%The following toolboxes are required to run this function:
% *Signal Processing Toolbox
% *Statistics and Machine Learning Toolbox
% *DSP System Toolbox
%
%Inputs:
% ecg________________ecg signal. Use [], if there is no ecg input.
% ecg_time___________times corresponding to the ecg signal. Use [], if there is no ecg_time input.
% ecg_fs_____________sampling frequency of the ecg signal. Use [], if there is no ecg_fs input.
% ip_________________impedance pneumograph signal
% ip_time____________times corresponding to the impedance pneumograph signal
% ip_fs______________sampling frequency of the impedance pneumograph signal
% hr_________________heart rate signal
% sats_______________oxygen saturation values
% sats_time__________times corresponding to the sats values
% hr_time____________times corresponding to the hr values
% model_input________enter 1 if you want to the run the machine
% learning model to remove ibi values greater than 20 which are
% artifactual due to shallow breathing or poor electrode
% placement. If this is not wanted, enter 0.
%
%Outputs: ibi__________________________interbreath-interval series consisting of the times between each breath
% ibi_time_____________________times corresponding to the end breath of each interval
% filtered_ip__________________cleaned impedance pneumograph signal
% filtered_ip_time_____________time vector corresponding to the cleaned impedance pneumograph signal
% filtered_ip_fs_______________sampling frequency of the cleaned impedance pneumograph signal
%
%
% The function is made up of the following sections;
% 1a) Check if the ECG signal is inverted (ie upside-down): r-peaks should be the highest peaks
% b) NaNs in the ECG must be found, such that the ip at these times can be kept as the original ip signal, without filtering out the r-peaks
% c) Find any segments of ECG where the ECG value does not change
% d) Find any segments of ECG where samples were missed
% e) Using c) and d), the times at which the ECG values did not change or were missed are concatenated such that the original ip can be inserted into the ip at these times (see 9a)
% f) Insert NaNs into the ECG where there are gaps in time
%
% 2) Detrend the ECG and identify the ECG R-peaks to derive RR intervals
%
% 3) Clean RR intervals
%
% 4) The rr interval (rri) series is be stretched or compressed so that there are an equal number of steps between each rri
%
% 5) Start of ip signal preprocessing: Remove NaNs from the ip signal
%
% 6) Resample the ip signal in accordance to the rri
%
% 7a) Filter ip signals
% b)filter out frequencies between 0.9-1.1 Hz to take out 1 Hz data, as this will represent the r-peak artefact
% c)filter out frequencies between 1.9-2.1 Hz to take out 2 Hz data, as this will represent the r-peak artefact
%
% 8) Resample filtered_ip2 to be back in time
%
% 9a) Reinsert the original ip signal where the ECG was NaN
% b) Filter out low frequencies below 0.5 Hz to tidy up signal, and remove outliers
%
% 10a) Make all data points in filtered_ip which correspond to missing ip data NaNs
% b) NaNs in the original ip must be found, such that the filtered_ip at these times can be made NaNs
%
% 11a) This part of the script analyses the ip signal to identify breaths using an adaptive threshold
% b) Identify breath peaks and compute the interbreath-interval (IBI)
% c) The interbreath interval series must be altered to reflect NaNs, where the ip signal contains NaNs
% d) This part of the code joins pauses which occur within 2 seconds of each other
%
%12) If requested, run the machine learning model using the remove_artifactual_longpauses function
%
% Tricia Adjei March 2021
% If you use this code please cite it as
% Adjei T, Purdy R, Jorge J, et al
% New method to measure interbreath intervals in infants for the assessment of apnoea and respiration
% BMJ Open Respiratory Research 2021;8:e001042. doi: 10.1136/bmjresp-2021-001042
%
% For further help please email Associate Prof. Caroline Hartley ([email hidden]).
%If the inputs do not include an ecg signal, an adaptive threshold of 0.5
%multiplied by the standard deviation of the cleaned ip signal is used. If
%the inputs include an ecg signal, an adaptive threshold of 0.4
%multiplied by the standard deviation of the cleaned ip signal is used.
%First reshape all input vectors
if isempty(ip)
disp('No IP signal')
ibi=[]; ibi_time=[]; filtered_ip=[]; filtered_ip_time=[]; filtered_ip_fs=[];
return
end
if exist('ecg')
if ~isempty(ecg)
ecg=reshape(ecg,length(ecg),1);
ecg_time=reshape(ecg_time,length(ecg_time),1);
[ecg_time,iecg,~]=unique(ecg_time);
ecg=ecg(iecg);
end
end
if exist('ip')
if ~isempty(ip)
ip=reshape(ip,length(ip),1);
ip_time=reshape(ip_time,length(ip_time),1);
[ip_time,iip,~]=unique(ip_time);
ip=ip(iip);
end
end
if exist('hr')
if ~isempty(hr)
hr=reshape(hr,length(hr),1);
hr_time=reshape(hr_time,length(hr_time),1);
[hr_time,ihr,~]=unique(hr_time);
hr=hr(ihr);
end
end
if exist('sats')
if ~isempty(sats)
sats=reshape(sats,length(sats),1);
sats_time=reshape(sats_time,length(sats_time),1);
[sats_time,isats,~]=unique(sats_time);
sats=sats(isats);
end
end
if isempty(ecg)
thresh=0.5;
final_ip_fs=50;
final_ip=spline(ip_time,ip,[ip_time(1):1/final_ip_fs:ip_time(end)]);
final_ip_time=[ip_time(1):1/final_ip_fs:ip_time(end)];
original_ip=ip;
original_ip_time=ip_time;
else
%% 1a) Check if the ECG signal is inverted (ie upside-down): r-peaks should be the highest peaks, and append NaNs to the end of the ECG if it is shorter than the ip signal
% Take a sample of the ecg signal to check if the signal could be inverted:
if (ecg_time(end)-ecg_time(1))<=(15*60)
test_window=length(ecg);
else
test_window_sec=300;
test_window=test_window_sec*ecg_fs;
end
ecg_sample=ecg(1:test_window);
ecg_sample(isnan(ecg_sample))=[];
if length(ecg_sample)>(12*ecg_fs)
%find the envelopes of the sample signal. The magnitude of the upper and
%lower envelopes will indicate which way 'up' the signal should be.
[upper,lower]=envelope(ecg_sample,(4*ecg_fs),'peak');
upper_avg=abs(mode(upper));
lower_avg=abs(mode(lower));
if upper_avg>lower_avg
ecg=ecg;
elseif lower_avg>upper_avg
ecg=-1*ecg;
end
else ecg=ecg;
end
%Append NaNs to the ECG if it is more than 1 minute shorter than the ip signal
if (ip_time(end)-ecg_time(end))>=60
shorter_by=ip_time(end)-ecg_time(end);
shorter_by_time=[ecg_time(end)+(1/ecg_fs):(1/ecg_fs):ecg_time(end)+shorter_by]';
shorter_by_length=length(shorter_by_time);
shorter_by_nans=nan(shorter_by_length,1);
ecg=[ecg;shorter_by_nans];
ecg_time=[ecg_time;shorter_by_time];
end
%% b) NaNs in the ECG must be found, such that the ip at these times can be kept as the original ip signal, without filtering out the r-peaks
ecgNaNs=find(isnan(ecg));
if ~isempty(ecgNaNs)
ecgNaNtime=ecg_time(find(isnan(ecg)));
ecgNaN_loc_st=zeros(1,length(ecgNaNs));
ecgNaN_loc_en=zeros(1,length(ecgNaNs));
for i=1:length(ecgNaNs)
if i==1
ecgNaN_loc_st(i)=ecgNaNs(i);
elseif ecgNaNs(i-1)~=ecgNaNs(i)-1
ecgNaN_loc_st(i)=ecgNaNs(i);
end
if i==length(ecgNaNs)
ecgNaN_loc_en(i)=ecgNaNs(i);
elseif ecgNaNs(i+1)~=ecgNaNs(i)+1
ecgNaN_loc_en(i)=ecgNaNs(i);
end
end
clear i
ecgNaN_loc_st(ecgNaN_loc_st==0)=[];
ecgNaN_loc_en(ecgNaN_loc_en==0)=[];
ecgNaN_time_st=ecg_time(ecgNaN_loc_st);
ecgNaN_time_en=ecg_time(ecgNaN_loc_en);
end
%% c) Find any segments of ECG where the ECG value does not change
nonconstantECG_loc=find(diff(ecg)~=0);
if ~isempty(nonconstantECG_loc)
nonconstantECG_loc=nonconstantECG_loc+1;
nonconstantECG_loc(find(diff(nonconstantECG_loc)~=1))=[];
nonconstantECG_loc(end)=[];
if nonconstantECG_loc(1)==2
nonconstantECG_loc=[1;nonconstantECG_loc];
end
constantECG_loc=1:length(ecg);
constantECG_loc(nonconstantECG_loc)=[];
ecg_repeats=ecg(constantECG_loc);
ecg_repeats_time=ecg_time(constantECG_loc);
constantECG_st=nan(length(ecg_repeats),1);
constantECG_en=nan(length(ecg_repeats),1);
for i=1:length(ecg_repeats)
if i==1
constantECG_st(i)=ecg_repeats_time(i);
elseif i~=1 & ecg_repeats(i-1)~=ecg_repeats(i)
constantECG_st(i)=ecg_repeats_time(i);
end
if i==length(constantECG_loc)
constantECG_en(i)=ecg_repeats_time(i);
elseif i~=length(constantECG_loc) & ecg_repeats(i+1)~=ecg_repeats(i)
constantECG_en(i)=ecg_repeats_time(i);
end
end
clear i
constantECG_st(isnan(constantECG_st))=[];
constantECG_en(isnan(constantECG_en))=[];
end
%% d) Find any segments of ECG where samples were missed
miss_data_loc=find(diff(ecg_time)>=(2/ecg_fs));
if ~isempty(miss_data_loc)
new_miss_data_loc=miss_data_loc;
miss_start_time=ecg_time(miss_data_loc); %Find start time of missing segments
miss_end_time=ecg_time(miss_data_loc+1);
miss_start_time=miss_start_time;
miss_end_time=miss_end_time;
end
%% e) constantECG_st and ecgNaN_time_st and miss_start_time need to be concatenated so that they
% can be used to reinsert the original ip signal at step 9a.
if exist('ecgNaN_time_st')==1 & exist('constantECG_st')==1 & exist('miss_start_time')==1
[insert_time_st,original_order,~]=unique([ecgNaN_time_st;constantECG_st;miss_start_time]);
insert_time_en=[ecgNaN_time_en;constantECG_en;miss_end_time];
insert_time_en=insert_time_en(original_order);
elseif exist('ecgNaN_time_st')==1 & exist('constantECG_st')==1 & exist('miss_start_time')==0
[insert_time_st,original_order,~]=unique([ecgNaN_time_st;constantECG_st]);
insert_time_en=[ecgNaN_time_en;constantECG_en];
insert_time_en=insert_time_en(original_order);
elseif exist('ecgNaN_time_st')==1 & exist('constantECG_st')==0 & exist('miss_start_time')==0
insert_time_st=ecgNaN_time_st;
insert_time_en=ecgNaN_time_en;
elseif exist('ecgNaN_time_st')==1 & exist('constantECG_st')==0 & exist('miss_start_time')==1
[insert_time_st,original_order,~]=unique([ecgNaN_time_st;miss_start_time]);
insert_time_en=[ecgNaN_time_en;miss_end_time];
insert_time_en=insert_time_en(original_order);
elseif exist('ecgNaN_time_st')==0 & exist('constantECG_st')==0 & exist('miss_start_time')==1
insert_time_st=miss_start_time;
insert_time_en=miss_end_time;
elseif exist('ecgNaN_time_st')==0 & exist('constantECG_st')==1 & exist('miss_start_time')==1
[insert_time_st,original_order,~]=unique([constantECG_st;miss_start_time]);
insert_time_en=[constantECG_en;miss_end_time];
insert_time_en=insert_time_en(original_order);
elseif exist('ecgNaN_time_st')==0 & exist('constantECG_st')==1 & exist('miss_start_time')==0
insert_time_st=constantECG_st;
insert_time_en=constantECG_en;
end
%% f) Insert NaNs into the ECG where there are gaps in time
%define gaps in the ECG as 2 missing data values, and find the gaps
if ~isempty(miss_data_loc)
if length(miss_data_loc)==1
miss_start_time=ecg_time(miss_data_loc); %Find start time of missing segments
miss_end_time=ecg_time(miss_data_loc+1);
else
% Only keep the last index of a missing segment
for md=1:length(miss_data_loc)-1
if md~=1 && miss_data_loc(md)+1==(miss_data_loc(md+1))
new_miss_data_loc(md)=NaN;
else
new_miss_data_loc(md)=miss_data_loc(md);
end
end
new_miss_data_loc(find(isnan(new_miss_data_loc)))=[];
miss_data_loc=new_miss_data_loc;
miss_start_time=ecg_time(miss_data_loc);
miss_end_time=ecg_time(miss_data_loc+1);
end
insert_times1=cell(1,length(miss_start_time));
%Find the missing times which correspond to the missing data
for m=1:length(miss_start_time)
insert_time=[(miss_start_time(m)+(1/ecg_fs)):(1/ecg_fs):(miss_end_time(m))];
if insert_time(end)==miss_end_time(m)
insert_times1{m}=insert_time(1:(length(insert_time)-1));
else
insert_times1{m}=insert_time;
end
end
%Insert NaNs into ecg and insert the missing times into ecg_time
insert_times=cat(2,insert_times1{:});
new_ecg_time=[ecg_time;insert_times'];
[sorted_ecg_time,order]=sort(new_ecg_time);
new_ecg=[ecg;NaN(length(insert_times),1)];
ecg_time=sorted_ecg_time;
ecg=new_ecg(order);
end
clear miss_data_loc new_miss_data_loc miss_start_time miss_end_time
%% 2) Detrend the ECG and identify the ECG R-peaks to derive RR intervals
%The ecg signal should be a row vector.
dim=size(ecg);
[M,I]=max(dim);
if I==1
ecg=reshape(ecg,1,[]);
ecg_time=reshape(ecg_time,1,[]);
end
window_sec=5; %detrend the ecg in a window, 5 seconds in length.
window=window_sec*ecg_fs;
lower=zeros(1,floor((length(ecg))/window)-1);
upper=zeros(1,floor((length(ecg))/window)-1);
nan_start=zeros(1,floor((length(ecg))/window)-1);
nan_end=zeros(1,floor((length(ecg))/window)-1);
detrend_ecg1=cell(1,floor((length(ecg))/window)-1);
detrend_ecg_time1=cell(1,floor((length(ecg))/window)-1);
time_of_r_peaks1=cell(1,floor((length(ecg))/window)-1);
time_of_r_peaks2=cell(1,floor((length(ecg))/window)-1);
r_peaks1=cell(1,floor((length(ecg))/window)-1);
r_peak_loc1=cell(1,floor((length(ecg))/window)-1);
time_of_volt_peaks1=cell(1,floor((length(ecg))/window)-1);
for j=1:floor((length(ecg))/window)-1
lower(j)=floor(1+((j-1)*window)); %lower window limit
upper(j)=floor(j*window); %upper window limit
windowed_signal1=ecg(lower(j):upper(j));
%tell matlab to move to the next loop if there are NaNs in the
%windowed segment, otherwise, carry on. The NaNs will be re-inserted
%into the rr series later.
if sum(isnan(windowed_signal1))~=0
nan_start(j)=ecg_time(lower(j));
nan_end(j)=ecg_time(upper(j));
detrend_ecg1{j}=windowed_signal1;
detrend_ecg_time1{j}=ecg_time(lower(j):upper(j));
time_of_r_peaks1{j}=NaN;
time_of_r_peaks2{j}=ecg_time(lower(j));
r_peak_loc1{j}=NaN;
r_peaks1{j}=NaN;
continue
end
windowed_signal=[0 0 0 windowed_signal1 0 0 0]; %add zeros before and after the windowed signal segment to help the detrending process
%find the polynomial coefficients which describe the wander in the
%windowed signal.
[p,s,mu]=polyfit((1:numel(windowed_signal)),windowed_signal,20);
trend=polyval(p,(1:numel(windowed_signal)),[],mu); %compute the polynomial which describes the wander
%subtract the wander polynomial from the windowed signal
detrend_ecg_windowed=windowed_signal-trend;
%fit an envelope to the detrended signal segment and find the
%indices where the envelope and signal intersect.
[upper_env,lower_env]=envelope(detrend_ecg_windowed,round(ecg_fs/(250/65)),'peak');
detrend_ecg_windowed=detrend_ecg_windowed(4:(length(detrend_ecg_windowed)-3));
detrend_ecg_time_windowed=ecg_time(lower(j):upper(j));
[touch,touchloc]=intersect(detrend_ecg_windowed,upper_env,'stable');
touchloc=touchloc';
%delete touches which are definitely not r_peaks
touchloc(touch<(0.9*std(touch)))=[];
touch(touch<(0.9*std(touch)))=[];
%move to the next segment if there are too few touches, or if the
%touches are too noisy.
if length(touch)<3 | std(touch)>0.7
detrend_ecg1{j}=windowed_signal1;
detrend_ecg_time1{j}=ecg_time(lower(j):upper(j));
time_of_r_peaks1{j}=NaN;
time_of_r_peaks2{j}=ecg_time(lower(j));
r_peak_loc1{j}=NaN;
r_peaks1{j}=NaN;
continue
end
%if the last touched peak is short, delete that touch
if touch(end)<0.5*touch(length(touch)-1)
touch(end)=[];
touchloc(end)=[];
end
touch_loc=touchloc;
%search to see if there is a potential r-peak which follows the last touch
remainder=detrend_ecg_windowed(max(touch_loc)+3:end);
remainder_pk=max(remainder);
if remainder_pk>0.7*detrend_ecg_windowed(max(touch_loc))
touch_loc=[touch_loc find(detrend_ecg_windowed==remainder_pk)];
end
%search the gaps between the touches, and define good touches as
%those which are of a consistent size.
gaps=diff(touch_loc);
gap_length=1:length(gaps);
good_touches=[touch_loc(1) touch_loc(gap_length(gaps>2*std(gaps))+1)];
check_gap=gap_length;
touch_locs=good_touches;
%if a gap is large, search for missing peaks in the gap
for bb=1:length(check_gap)-1
num_miss_val=round(gaps(check_gap(bb)+1)/gaps(check_gap(bb)))-1; %This denotes the number of missing peaks
if num_miss_val+1<2
continue
end
approx_ind=round((touch_loc(check_gap(bb)+2)-touch_loc(check_gap(bb)+1))/(num_miss_val+1)); %approx_ind approximates the jump to the next r peak
%for the number of num_miss_vals, the location and
%value of r-peaks missed by the envelope are found.
for mi=1:length(num_miss_val)
miss_ind(mi)=[touch_loc(check_gap(bb)+1)+(mi*approx_ind)];
approx_range=[miss_ind(mi)-5:miss_ind(mi)+5]; %this is the approximate range of indices in which to search for a missed r-peak
miss_dat=detrend_ecg_windowed(approx_range);
[highest,highest_loc]=sort(miss_dat,'descend');
miss_val(mi)=highest(mi);
new_touch_loc1(mi)=(approx_range(mi)-1)+highest_loc(mi);
end
touch_locs=[touch_locs new_touch_loc1];
end
r_peaks_locs_windowed=sort(touch_locs);
additional=lower(j)-1; %this is the signal index immediately before the start of the data window.
%If the number of data points between additional and the first
%value of r_peaks_locs_windowed is too small to be an rr interval, delete the first value of r_peaks_locs_windowed.
if j~=1 & ~isempty(r_peak_loc1{j-1}) & ((r_peaks_locs_windowed(1)+additional)-(r_peak_loc1{j-1}(end)))<(mean(diff(r_peak_loc1{j-1}))*0.5)
r_peaks_locs_windowed(1)=[];
end
%% Identify r-peaks
detrend_ecg_time1{j}=detrend_ecg_time_windowed;
detrend_ecg1{j}=detrend_ecg_windowed;
r_peak_loc1{j}=r_peaks_locs_windowed+additional; %location of r-peaks
r_peaks1{j}=detrend_ecg_windowed(r_peaks_locs_windowed); %r-peaks
time_of_r_peaks1{j}=detrend_ecg_time_windowed(r_peaks_locs_windowed); %Time of r-peaks
time_of_r_peaks2{j}=detrend_ecg_time_windowed(r_peaks_locs_windowed); %Time of r-peaks
%If the std of the r-peaks is too high, the window is too noisy
if std(r_peaks1{j})>=1
time_of_r_peaks1{j}=NaN(1,length(time_of_r_peaks2{j}));
r_peak_loc1{j}=NaN(1,length(r_peak_loc1{j}));
r_peaks1{j}=NaN(1,length(r_peaks1{j}));
continue
end
r_peak_loc1{j}=r_peaks_locs_windowed+additional; %location of r-peaks
r_peaks1{j}=detrend_ecg_windowed(r_peaks_locs_windowed); %r-peaks
time_of_volt_peaks1{j}=detrend_ecg_time_windowed(r_peaks_locs_windowed); %Time of all peaks
end
clear windowed_signal
detrended_ecg=cat(2,detrend_ecg1{:}); %concatenate the windowed segments
detrended_ecg_time=cat(2,detrend_ecg_time1{:});
detrended_ecg_fs=ecg_fs;
r_peak_loc=cat(2,r_peak_loc1{:});
r_peaks=cat(2,r_peaks1{:});
time_of_r_peaks_nan=cat(2,time_of_r_peaks1{:});
time_of_r_peaks=cat(2,time_of_r_peaks2{:});
rri=abs(diff(time_of_r_peaks_nan)); %time intervals between r-peaks
rri_time=(time_of_r_peaks(2:end)); %timestamps of rr intervals
%% 3) Clean rr intervals
rri(rri<=0.25| rri>=1.5)=NaN;
%% 4) The rr interval (rri) series is to be stretched or compressed so that there are an equal number of steps between each rri
stretched_rri_time=cell(1,length(rri_time));
division=50; %this defines the number of steps between each rri
division1=division-1;
for i=1:(length(rri_time)-1)
increment_steps=(rri_time(i+1)-rri_time(i))/division1;
for ii=1:(division1)
stretched_rri_time{i}(ii)=rri_time(i)+((ii-1)*increment_steps);
end
end
clear i ii
stretched_rri_time=cat(2,stretched_rri_time{:});
%% 5) Start of ip signal preprocessing: Remove NaNs from the ip signal
% First, get rid of NaNs in the ip signal and remove any square edges in the ip signal
original_ip=ip;
original_ip_time=ip_time;
ip_nan_indices=find(isnan(ip));
ip(ip_nan_indices)=[];
ip_time(ip_nan_indices)=[];
ip_copy=ip;
ip_copy_time=ip_time;
gapps=find(diff(ip)==0);
gap_st_loc=zeros(1,length(gapps));
gap_en_loc=zeros(1,length(gapps));
for g=1:length(gapps)
if g==1
gap_st_loc(g)=gapps(g);
elseif gapps(g-1)~=gapps(g)-1
gap_st_loc(g)=gapps(g);
end
if g==length(gapps)
gap_en_loc(g)=gapps(g);
elseif gapps(g+1)~=gapps(g)+1
gap_en_loc(g)=gapps(g);
end
end
gap_st_loc(gap_st_loc==0)=[];
gap_en_loc(gap_en_loc==0)=[];
gap_en_loc=gap_en_loc+1;
for h=1:length(gap_st_loc)
ip_copy([gap_st_loc(h):gap_en_loc(h)])=NaN;
ip_copy_time([gap_st_loc(h):gap_en_loc(h)])=NaN;
end
new_ip=ip_copy;
new_ip_time=ip_copy_time;
new_ip(isnan(ip_copy))=[];
new_ip_time(isnan(ip_copy))=[];
%% 6) Resample the ip signal in accordance to the rri
% (Lee, H., Rusin, C.G., Lake, D.E., Clark, M.T., Guin, L., Smoot, T.J., Paget-Brown, A.O., Vergales, B.D., Kattwinkel, J., Moorman, J.R. and Delos, J.B., 2011.
% A new algorithm for detecting central apnea in neonates. Physiological measurement, 33(1), p.1).
rri_sampled_ip=spline(new_ip_time,new_ip,stretched_rri_time); %Resample rri according to stretched_rri_time
%% 7a) Filter ip signals
% filter out frequencies less than 0.1 Hz to take out noise
fltr=designfilt('highpassfir','FilterOrder',floor(length(rri_sampled_ip)/1000),'CutoffFrequency',0.1,'DesignMethod','window','Window','Hamming','ScalePassband',false,'SampleRate',division);
filt_ip=filtfilt(fltr,rri_sampled_ip);
delay=mean(grpdelay(fltr,length(filt_ip),division));
filt_ip(1:delay)=[];
filt_ip_time=stretched_rri_time;
filt_ip_time(1:delay)=[];
%% b)filter out frequencies between 0.9-1.1 Hz to take out 1 Hz data, as this will represent the r-peak artefact
nyquist=division/2;
wo=1/nyquist;
bw=wo/4;
[SOS1,G1]=iirnotch(wo,bw,10);
filtered_ip1=filtfilt(SOS1,G1,filt_ip);
delay=mean(grpdelay(SOS1,length(filt_ip),division));
filtered_ip1(1:delay)=[];
filtered_ip_time1=filt_ip_time;
filtered_ip_time1(1:delay)=[];
%% c)filter out frequencies between 1.9-2.1 Hz to take out 2 Hz data, as this will represent the r-peak artefact
wo=2/nyquist;
bw=wo/4;
[SOS2,G2]=iirnotch(wo,bw,10);
filtered_ip2=filtfilt(SOS2,G2,filtered_ip1);
delay=mean(grpdelay(SOS2,length(filtered_ip1),division));
filtered_ip2(1:delay)=[];
filtered_ip_time2=filtered_ip_time1;
filtered_ip_time2(1:delay)=[];
%% 8) Resample filtered_ip2 to be back in time
final_ip_fs=50; %Resample the ip signal at a frequency of 50 Hz
X=filtered_ip_time2;
Y=filtered_ip2;
final_ip_time=filtered_ip_time2(1):1/final_ip_fs:filtered_ip_time2(end);
final_ip=spline(X,Y,final_ip_time);
%% 9a) Reinsert the original ip signal where the ECG was NaN or did not change, or missed samples
if exist('insert_time_st')==1
insertip=spline(original_ip_time,original_ip,final_ip_time);
diff_insert=insert_time_en-insert_time_st;
f=find(diff_insert>=60); %The original ip signal should only be inserted if the segment to be inserted is longer than 60s
for i=1:length(f)
insertloc1=find(final_ip_time<=insert_time_st(f(i)),1,'last');
insertloc2=find(final_ip_time>=insert_time_en(f(i)),1,'first');
if ~isempty(insertloc1) && ~isempty(insertloc2)
final_ip(insertloc1:insertloc2)=insertip(insertloc1:insertloc2);
end
end
end
end
clear i
%% b) Filter out low frequencies below 0.5 Hz to tidy up signal, and remove outliers
high_pass_filt=designfilt('highpassfir', 'FilterOrder', 300, 'CutoffFrequency', 0.5, 'SampleRate', final_ip_fs);
filtered_final_ip=filter(high_pass_filt,final_ip);
delay=mean(grpdelay(high_pass_filt));
delay_time=delay*(1/final_ip_fs);
filtered_final_ip(1:delay)=[];
filtered_final_ip_time=final_ip_time-delay_time;
filtered_final_ip_time(1:delay)=[];
filtered_final_ip_fs=final_ip_fs;
%Outliers will be identified and removed from every 12 hours of the ip signal.
window_pts=12*3600*filtered_final_ip_fs;
windowed_ip1=cell(1,ceil(length(filtered_final_ip)/window_pts));
windowed_ip_time1=cell(1,ceil(length(filtered_final_ip_time)/window_pts));
for j=1:ceil(length(filtered_final_ip)/window_pts)
if j<ceil(length(filtered_final_ip)/window_pts)
lower=1+((j-1)*window_pts); %lower window limit
upper=j*window_pts; %upper window limit
elseif (length(filtered_final_ip)-((j-1)*window_pts))>=(60*filtered_final_ip_fs)
lower=1+((j-1)*window_pts);
upper=length(filtered_final_ip);
else
break
end
windowed_ip=filtered_final_ip(lower:upper);
windowed_ip_time=filtered_final_ip_time(lower:upper);
temp_windowed_ip1=windowed_ip(windowed_ip>0);
high_outlier=6*prctile(temp_windowed_ip1,90);
outliers1=find(windowed_ip>high_outlier);
windowed_ip(outliers1)=high_outlier;
temp_windowed_ip2=windowed_ip(windowed_ip<0);
low_outlier=6*prctile(temp_windowed_ip2,10);
outliers2=find(windowed_ip<low_outlier);
windowed_ip(outliers2)=low_outlier;
windowed_ip1{j}=windowed_ip;
windowed_ip_time1{j}=windowed_ip_time;
end
filtered_ip=cat(2,windowed_ip1{:});
filtered_ip_time=cat(2,windowed_ip_time1{:});
filtered_ip_fs=filtered_final_ip_fs;
%% 10a) Make all data points in filtered_ip which correspond to missing ip data NaNs
%define gaps in the original ip signal as 2 missing data values,
%and find the gaps
miss_data_loc=find(abs(diff(original_ip_time))>=(2/ip_fs));
if ~isempty(miss_data_loc)
new_miss_data_loc=miss_data_loc;
for md=1:length(miss_data_loc)-1
if md~=1 & miss_data_loc(md)+1==(miss_data_loc(md+1))
new_miss_data_loc(md)=NaN;
end
end
new_miss_data_loc(find(isnan(new_miss_data_loc)))=[];
miss_data_loc=new_miss_data_loc;
miss_start_time=original_ip_time(miss_data_loc);
miss_end_time=original_ip_time(miss_data_loc+1);
nearest_start_loc=zeros(1,length(miss_start_time));
nearest_end_loc=zeros(1,length(miss_start_time));
check_start_loc=zeros(1,length(miss_start_time));
check_end_loc=zeros(1,length(miss_start_time));
missing_periods=cell(1,length(miss_start_time));
check_start_periods=cell(1,length(miss_start_time));
check_end_periods=cell(1,length(miss_start_time));
%find nearest times in filtered_ip_time which correspond to the
%missing times.
for m=1:length(miss_start_time)
[~,nearest_start_loc(m)]=min(abs(miss_start_time(m)-filtered_ip_time));
[~,nearest_end_loc(m)]=min(abs(miss_end_time(m)-filtered_ip_time));
if (nearest_start_loc(m)-(60*filtered_ip_fs))>=1
check_start_loc(m)=nearest_start_loc(m)-(60*filtered_ip_fs);
else check_start_loc(m)=1;
end
if (nearest_end_loc(m)+(60*filtered_ip_fs))<=length(filtered_ip)
check_end_loc(m)=nearest_end_loc(m)+(60*filtered_ip_fs);
else check_end_loc(m)=length(filtered_ip);
end
%Define indices of periods of missing data
missing_periods{m}=[nearest_start_loc(m):nearest_end_loc(m)];
check_start_periods{m}=[check_start_loc(m):nearest_start_loc(m)-1];
check_end_periods{m}=[nearest_end_loc(m)+1:check_end_loc(m)];
%remove outliers occurring immediatley before a segment of missing
%data.
to_check=filtered_ip(check_start_periods{m});
take_out1=find(to_check>=high_outlier | to_check<=low_outlier,1,'first');
take_out2=find(to_check(1:take_out1)<0,1,'last');
take_out3=check_start_periods{m}(take_out2:end);
filtered_ip(take_out3)=NaN;
clear take_out1 take_out2 take_out3
%remove outliers occurring immediatley after a segment of missing
%data.
to_check1=filtered_ip(check_end_periods{m});
take_out1=find(to_check1>=high_outlier | to_check1<=low_outlier,1,'last');
take_out2=find(to_check1(take_out1:end)<0,1,'first');
if (take_out1+take_out2) <= length(to_check1)
take_out3=check_end_periods{m}(1:(take_out1+take_out2));
else
take_out3=check_end_periods{m}(1:take_out1);
end
filtered_ip(take_out3)=NaN;
%filtered_ip currently contains interpolated values where data
%should be missing, so replace data during missing periods with NaNs
filtered_ip(missing_periods{m})=NaN;
end
data_restart_loc=nearest_end_loc+1; %Time at which data restarts after a period of missing data
end
%% b) NaNs in the original ip must be found, such that the filtered_ip at these times can be made NaNs
ipNaNs=find(isnan(original_ip));
if ~isempty(ipNaNs)
ipNaNtime=original_ip_time(find(isnan(ip)));
ipNaN_loc_st=zeros(1,length(ipNaNs));
ipNaN_loc_en=zeros(1,length(ipNaNs));
for i=1:length(ipNaNs)
if i==1
ipNaN_loc_st(i)=ipNaNs(i);
elseif ipNaNs(i-1)~=ipNaNs(i)-1
ipNaN_loc_st(i)=ipNaNs(i);
end
if i==length(ipNaNs)
ipNaN_loc_en(i)=ipNaNs(i);
elseif ipNaNs(i+1)~=ipNaNs(i)+1
ipNaN_loc_en(i)=ipNaNs(i);
end
end
clear i
ipNaN_loc_st(ipNaN_loc_st==0)=[];
ipNaN_loc_en(ipNaN_loc_en==0)=[];
ipNaN_time_st=original_ip_time(ipNaN_loc_st);
ipNaN_time_en=original_ip_time(ipNaN_loc_en);
ip_nan_times2=cell(1,length(ipNaN_time_st));
for i=1:length(ipNaN_time_st)
ip_nan_times2{i}=find(filtered_ip_time>=ipNaN_time_st(i) & filtered_ip_time<=ipNaN_time_en(i));
end
clear i
nan_time_loc=cat(2,ip_nan_times2{:});
filtered_ip(nan_time_loc)=NaN;
end
%% 11a) This part of the script analyses the ip signal
clearvars -except ecgNaN_time_st ecgNaN_time_en thresh ip ip_time ip_fs original_ip original_ip_time filtered_ip_time filtered_ip filtered_ip_fs data_restart_times hr sats sats_time hr_time model_input
signal=filtered_ip;
time=filtered_ip_time;
%a) Set signal to NaN if the hr signal has NaNs.
checkHR=find(isnan(hr));
%but not if there are large chunks of missing HR data
diffcheckHR=diff(checkHR);
fdhr=find(diffcheckHR>1);
if isempty(fdhr)
checkHR=[];
else
todelete=[];
for ind=1:length(fdhr)
if ind==1
len=checkHR(fdhr(ind))-checkHR(1);
if len>60
todelete=[1:fdhr(ind)];
end
else
len=checkHR(fdhr(ind))-checkHR(fdhr(ind-1)+1);
if len>60
todelete=[todelete,fdhr(ind-1)+1:fdhr(ind)];
end
end
end
checkHR(todelete)=[];
end
if ~isempty(checkHR)
time_nanHR=hr_time(checkHR);
for i=1:length(checkHR) %for each NaN change the ip signal to NaN for 5 seconds around it
st_nan=time_nanHR(i)-2.5;
en_nan=time_nanHR(i)+2.5;
ind=find(time>=st_nan & time<=en_nan);
signal(ind)=NaN;
end
end
clear i;
diff_ip=diff(original_ip);
ind=find(diff_ip==0);
diff_ind=diff(ind);
i=1;
while i<length(diff_ind)
if diff_ind(i)==1
f=find(diff_ind(i:end)>1,1,'first');
if f>1*ip_fs %second is fs
%Find other segments of data which reach the cap within 5
%seconds of the first segment.
jump=1;
while jump==1
g=find(diff_ind(i+f:end)>1,1,'first');
if ~isempty(g)
if g>1*ip_fs && diff_ind(i+f-1)<5*ip_fs
f=f+g;
else
jump=2;
end
else
jump=2;
end
end
ind_nan=ind(i:i+f-1);
st_time=original_ip_time(ind_nan(1))-2.5; %take 2.5 seconds either side
en_time=original_ip_time(ind_nan(end))+2.5;
t_ind=find(time>=st_time & time<=en_time);
signal(t_ind)=NaN;
end
else
f=1;
end
i=i+f;
end
clear i f;
%Find short segments of signal left over and set them to NaN
f=find(isnan(signal));
d=diff(f);
ff=find(d>1 & d<2*ip_fs);
for i=1:length(ff)
signal(f(ff(i)):f(ff(i)+1))=NaN;
end
clear i f d ff;
f=find(isnan(signal));
if ~isempty(f)
d=diff(f);
ff=find(d>1);
data_restart_times=time(f(ff)+1);
data_end=time(f(ff+1)-1);
if f(1)==1
data_end=[time(1),data_end];
else
data_end=[time(f(1)-1),data_end];
end
data_end(end)=[];
duration_stop=data_restart_times-data_end;
clear f d ff;
end
%% b) Identify breath peaks and compute the interbreath-interval (IBI)
if (time(end)-time(1))<=15*60
initial_thresh_window=30; % window length in which to compute ip threshold
else
initial_thresh_window=600;
end
window_pts=initial_thresh_window*filtered_ip_fs;
sig_start=find(~isnan(signal),1,'first');
temp_sig=signal(sig_start:end);
temp_tim=time(sig_start:end);
all_thresholds=cell(1,ceil(length(signal)));
crossings_windowed=zeros(1,ceil(length(signal)));
x=cell(1,ceil(length(signal)/window_pts));
%Detect breaths and interbreath intervals in nonoverlapping windows
window_pts_sum=window_pts;
j=1;
extra=0;
n_breaths=15;
while window_pts_sum<=length(temp_sig)
lower=window_pts_sum-window_pts+1;
upper=window_pts_sum;
windowed_ip=temp_sig(lower:upper);
windowed_ip_time=temp_tim(lower:upper);
if exist('ecgNaN_time_st')==1
for i=1:length(ecgNaN_time_st)
if windowed_ip_time(1)>=ecgNaN_time_st(i) & windowed_ip_time(end)<=ecgNaN_time_en(i)
thresh=0.5;
break
else
thresh=0.4;
end
end
else thresh=0.4;
end
%If the windowed signal only contains NaNs,
%skip to the next window
if sum(isnan(windowed_ip))==length(windowed_ip) || sum(diff(windowed_ip))==0
threh=all_thresholds{1:j};
if isempty(threh)
extra=1;
end
window_pts_sum=window_pts_sum+window_pts;
continue
end
if lower==1 | extra==1
threshold=thresh*std(windowed_ip,'omitnan');
rel_thres=windowed_ip-threshold;
pos=find(rel_thres>=0);
diffposs=diff(pos);
f=find(diffposs~=1);
% length(pos) denotes the last threshold crossing.
% If the signal in the window ends at a point above
% the threshold, f is left as f. But if it ends at
% another point, f must also include length(pos).
if ~ismember(window_pts,pos)
f=[f length(pos)];
end
breath_time_windowed=windowed_ip_time(pos(f));
%Delete if only 1 crossing within 5 seconds either side
f=find(diff(breath_time_windowed)>5);
if ~isempty(f)
g=find(diff(f)==1);
if ~isempty(g)
breath_time_windowed(f(g)+1)=[];
end
end
x{j}=find(abs(diff(breath_time_windowed))<0.3); %Ignore breaths which occur less than 0.3 seconds after the previous peak.
breath_time_windowed(x{j}+1)=[];
crossings_windowed(1:length(breath_time_windowed))=breath_time_windowed; %Store times at which the ip signal crossed the threshold
all_thresholds{j}=repmat(threshold,1,length(breath_time_windowed)); %Store the threshold used in each window
num_breath=find(crossings_windowed~=0);
num_breath=length(num_breath); %Number of breaths identified
min_num_breaths=num_breath;
%If the number of identified breaths is less than n_breath (ie the required number of previous breaths),
%a fixed threshold, calculated over the same length as initial_thresh_window, is used.
if min_num_breaths<n_breaths
clear all_thresholds
crossings_windowed1=cell(1,ceil(length(temp_sig)/window_pts));
all_thresholds=cell(1,ceil(length(temp_sig)/window_pts));
%Detect breaths in nonoverlapping windows if the number of
%breaths computed in the first window is less than that
%required to update the threshold
for j=1:ceil(length(temp_sig)/window_pts)
if j<ceil(length(temp_sig)/window_pts)
lower=1+((j-1)*window_pts); %lower window limit
upper=j*window_pts; %upper window limit
elseif (length(temp_sig)-((j-1)*window_pts))>=(5*filtered_ip_fs)
lower=1+((j-1)*window_pts);
upper=length(temp_sig);
else
break
end
windowed_ip=temp_sig(lower:upper);
windowed_ip_time=temp_tim(lower:upper);
if sum(isnan(windowed_ip))==length(windowed_ip) || sum(diff(windowed_ip))==0
continue
end
threshold=thresh*std(windowed_ip,'omitnan');
rel_thres=windowed_ip-threshold;
pos=find(rel_thres>=0);
diffposs=diff(pos);
f=find(diffposs~=1);
if ~ismember(window_pts,pos)
f=[f length(pos)];
end
breath_time_windowed=windowed_ip_time(pos(f));
%Delete if only 1 crossing within 5 seconds either side
f=find(diff(breath_time_windowed)>5);
if ~isempty(f)
g=find(diff(f)==1);
if ~isempty(g)
breath_time_windowed(f(g)+1)=[];
end
end
crossings_windowed1{j}=breath_time_windowed; %Store times at which the ip signal crossed the threshold
all_thresholds{j}=repmat(threshold,1,length(crossings_windowed1{j})); %Store the threshold used in each window
end
crossings_windowed=cat(2,crossings_windowed1{:});
x=find(abs(diff(crossings_windowed))<0.3); %Ignore breaths which occur less than 0.3 seconds after the previous peak.
crossings_windowed(x+1)=[];
break
end
%This part of the loop finds the first breath, using an adaptive
%threshold.
extra=0;
nbreaths=n_breaths-1; %n_breaths is the number of previous breaths over which to compute the adaptive threshold
j=j+1;
lastbreath=crossings_windowed(num_breath-nbreaths);
lastbreath_loc1=find(temp_tim==lastbreath);
lastbreath_loc2=find(temp_tim==crossings_windowed(num_breath));
temp_thresh=thresh*std(temp_sig(lastbreath_loc1:lastbreath_loc2),'omitnan');
sig_len=(lastbreath_loc2+1):length(temp_sig);
pot_cross=sig_len(temp_sig((lastbreath_loc2+1):end)>=temp_thresh); %pot_cross contains the indices of all signal values which are above the threshold
gd_point=find(diff(pot_cross)~=1); %find pot_cross values which are not preceded by another pot_cross value
gd_point(gd_point==(pot_cross(gd_point)-lastbreath_loc2))=[];
pot_cross=pot_cross(gd_point);
pot_cross(pot_cross==length(temp_sig))=[];
pot_cross=pot_cross((temp_tim(pot_cross)-temp_tim(lastbreath_loc2))>0.3);
cross_loc=find(temp_sig(pot_cross+1)<temp_thresh,1,'first');
if isempty(cross_loc)
break
end
next_breath_loc=pot_cross(cross_loc);
space=find(crossings_windowed==0,1,'first');
crossings_windowed(space)=temp_tim(next_breath_loc);
all_thresholds{j}=temp_thresh;
previous_breath_loc=lastbreath_loc2;
window_pts_sum=previous_breath_loc; %window_pts_sum must be defined so that the next window can be found
window_pts=next_breath_loc-previous_breath_loc;
window_pts_sum=window_pts_sum+window_pts;
if window_pts_sum>length(temp_sig)
break
end
else
j=j+1;
num_breath=find(crossings_windowed~=0);
num_breath=length(num_breath);
lastbreath=crossings_windowed(num_breath-nbreaths);
lastbreath_loc1=find(temp_tim==lastbreath);
lastbreath_loc2=find(temp_tim==crossings_windowed(num_breath));
temp_thresh=thresh*std(temp_sig(lastbreath_loc1:lastbreath_loc2),'omitnan');
sig_len=(lastbreath_loc2+1):length(temp_sig);
pot_cross=sig_len(temp_sig((lastbreath_loc2+1):end)>=temp_thresh);
gd_point=find(diff(pot_cross)~=1);
gd_point(gd_point==(pot_cross(gd_point)-lastbreath_loc2))=[];
pot_cross=pot_cross(gd_point);
pot_cross(pot_cross==length(temp_sig))=[];
pot_cross=pot_cross((temp_tim(pot_cross)-temp_tim(lastbreath_loc2))>0.3);
cross_loc=find(temp_sig(pot_cross+1)<temp_thresh,1,'first');
if isempty(cross_loc)
break
end
next_breath_loc=pot_cross(cross_loc);
space=find(crossings_windowed==0,1,'first');
crossings_windowed(space)=temp_tim(next_breath_loc);
all_thresholds{j}=temp_thresh;
previous_breath_loc=lastbreath_loc2;
window_pts=next_breath_loc-previous_breath_loc;
window_pts_sum=window_pts_sum+window_pts;
if window_pts_sum>length(temp_sig)
break
end
end
end
all_thresholds=cat(2,all_thresholds{:});
good=find(crossings_windowed~=0);
crossings=crossings_windowed(good);
threshold_crossings=crossings;
%Compute the interbreath intervals
ibi=abs(diff(crossings));
ibi_time=crossings(2:end);
%% c) The interbreath interval series must be altered to reflect NaNs, where the ip signal contains NaNs to identify breaths using an adaptive threshold
if sum(isnan(signal))~=0
time_aft_nan=zeros(1,length(signal));
time_bef_nan=zeros(1,length(signal));
getrid_ibi_loc=zeros(1,length(signal));
%The times at which sequences of NaNs start or end are found as time_bef_nan and time_aft_nan
for i=1:length(time)-1
if i==1 && isnan(signal(i)) && ~isnan(signal(i+1))
time_bef_nan(i)=time(i);
time_aft_nan(i)=time(i);
elseif i==1 && isnan(signal(i)) && isnan(signal(i+1))
time_bef_nan(i)=time(i);
elseif i~=1 && ~isnan(signal(i-1)) && isnan(signal(i))
time_bef_nan(i)=time(i);
elseif i~=1 && isnan(signal(i)) && ~isnan(signal(i+1))
time_aft_nan(i)=time(i);
if find(crossings>time_aft_nan(i),1,'first')~=1 & crossings(end)>time_aft_nan(i)
getrid_ibi_loc(i)=find(ibi_time==(crossings(find(crossings>time_aft_nan(i),1,'first'))));
end
elseif i==length(time)-1 && isnan(signal(i+1))
time_aft_nan(i)=time(i);
end
end
clear i;
getrid_ibi_loc(getrid_ibi_loc==0)=[];
if ~isempty(getrid_ibi_loc)
getrid_ibi_loc=unique(getrid_ibi_loc);
ibi(getrid_ibi_loc)=[];
ibi_time(getrid_ibi_loc)=[];
end
time_bef_nan(time_bef_nan==0)=[];
time_aft_nan(time_aft_nan==0)=[];
ibi_time=[ibi_time time_bef_nan];
[nanibi_time,ibi_ord,new_ord]=unique(ibi_time);
ibi=[ibi NaN(1,length(time_bef_nan))];
nanibi=ibi(ibi_ord);
%To account for breath pauses which are followed by
%segments of NaNs, 'additional' ibi values must be found. In other words,
%pauses followed by NaNs and not by threshold crossings must be found.
additional_ibi_time=zeros(1,length(time_bef_nan));
additional_ibi=zeros(1,length(time_bef_nan));
for i=1:length(time_bef_nan)
if ~isempty(find(crossings<time_bef_nan(i),1,'last'))
last_cross_loc=find(time==crossings(find(crossings<time_bef_nan(i),1,'last')));
if sum(isnan(signal(last_cross_loc:find(time<time_bef_nan(i),1,'last'))))==0
if time(find(time<time_bef_nan(i),1,'last'))-crossings(find(crossings<time_bef_nan(i),1,'last')) >=3
additional_ibi_time(i)=time(find(time<time_bef_nan(i),1,'last'));
additional_ibi(i)=time(find(time<time_bef_nan(i),1,'last'))-crossings(find(crossings<time_bef_nan(i),1,'last'));
end
end
end
end
additional_ibi_time(additional_ibi_time==0)=[];
additional_ibi(additional_ibi==0)=[];
ibi_time1=[nanibi_time additional_ibi_time];
ibi1=[nanibi additional_ibi];
clear i
%To account for breath pauses which are preceded and followed by
%segments of NaNs, 'additional' ibi values must be found; pauses preceded
%and followed by NaNs and not by threshold crossings must be found.
additional_ibi_time=zeros(1,length(time_aft_nan));
additional_ibi=zeros(1,length(time_aft_nan));
if length(time_bef_nan)>=2
for i=2:length(time_bef_nan)
if isempty(find(crossings>time_aft_nan(i-1) & crossings<time_bef_nan(i)))
if time(find(time<time_bef_nan(i),1,'last'))-time(find(time>time_aft_nan(i-1),1,'first')) >=3
additional_ibi_time(i)=time(find(time<time_bef_nan(i),1,'last'));
additional_ibi(i)=time(find(time<time_bef_nan(i),1,'last'))-time(find(time>time_aft_nan(i-1),1,'first'));
end
end
end
additional_ibi_time(additional_ibi_time==0)=[];
additional_ibi(additional_ibi==0)=[];
ibi_time2=[ibi_time1 additional_ibi_time];
ibi2=[ibi1 additional_ibi];
clear i
else
additional_ibi_time(additional_ibi_time==0)=[];
additional_ibi(additional_ibi==0)=[];
ibi_time2=[ibi_time1 additional_ibi_time];
ibi2=[ibi1 additional_ibi];
end
%Pauses preceded by NaNs must be accounted for.
additional_ibi_time=zeros(1,length(time_aft_nan));
additional_ibi=zeros(1,length(time_aft_nan));
for i=1:length(time_aft_nan)
if crossings(end)>=time_aft_nan(i)
nxt_cross_loc=find(time==crossings(find(crossings>time_aft_nan(i),1,'first')));
if sum(isnan(signal(find(time>time_aft_nan(i),1,'first'):nxt_cross_loc)))==0
if crossings(find(crossings>time_aft_nan(i),1,'first'))-time(find(time>time_aft_nan(i),1,'first')) >=3
additional_ibi_time(i)=crossings(find(crossings>time_aft_nan(i),1,'first'));
additional_ibi(i)=crossings(find(crossings>time_aft_nan(i),1,'first'))-time(find(time>time_aft_nan(i),1,'first'));
end
end
end
end
additional_ibi_time(additional_ibi_time==0)=[];
additional_ibi(additional_ibi==0)=[];
ibi_time3=[ibi_time2 additional_ibi_time];
ibi3=[ibi2 additional_ibi];
[final_ibi_time,prev_ord]=unique(ibi_time3);
final_ibi=ibi3(prev_ord);
ibi=final_ibi;
ibi_time=final_ibi_time;
clear i
end
filtered_ip=signal;
filtered_ip_time=time;
%% d) This part of the code joins pauses which occur within 2 seconds of each other
pause_ibi_loc=find(ibi>=5); %Identify pauses by finding ibi greater than a specified value
if ~isempty(pause_ibi_loc)
pause_start=ibi_time(pause_ibi_loc)-ibi(pause_ibi_loc); %start of pause
pause_end=ibi_time(pause_ibi_loc); %end of pause
c=find((pause_start(2:end)-pause_end(1:end-1))<=2); %find pauses less than 2 seconds apart
if ~isempty(c)
%Find chains of pauses which are close; cindex1 contains the start indices of the chains, and index2 contains the end indices
for i=1:length(c)
if i==1 || c(i)-1~=c(i-1)
cindex1(i)=c(i);
end
if i==length(c) || c(i)+1~=c(i+1)
cindex2(i)=c(i)+1;
end
end
clear i
cindex1(cindex1==0)=[];
cindex2(cindex2==0)=[];
%alter the ibi and ibi_time values to join the close pauses
to_delete=cell(1,size(cindex1,2));
for i=1:size(cindex1,2)
index1=pause_ibi_loc(cindex1(i));
index2=pause_ibi_loc(cindex2(i));
ibi(index1)=sum(ibi(index1:index2));
ibi_time(index1)=ibi_time(index2);
to_delete{i}=[index1+1:index2];
end
to_delete=cat(2,to_delete{:});
ibi(to_delete)=[];
ibi_time(to_delete)=[];
clear pause_start pause_end
pause_ibi_loc2=find(ibi>=5);
if ~isempty(pause_ibi_loc2)
pause_start=ibi_time(pause_ibi_loc2)-ibi(pause_ibi_loc2);
pause_end=ibi_time(pause_ibi_loc2);
else
pause_start=[];
pause_end=[];
end
end
end
%% 12) If requested, run the machine learning model using the remove_artifactual_longpauses function
if model_input==1
current_folder=pwd;
addpath([current_folder '/subfunctions/']);
[ibi_corrected,response_predictions]=remove_artifactual_longpauses(ibi,ibi_time,...
sats_time,hr_time,sats,hr,filtered_ip_time,filtered_ip);
ibi=ibi_corrected;
end
end

ip_to_interbreath_intervals.m at commit 72a9ae4, no license · at the source

Overview

Authors: Coen S Zandvoort¹, Fatima Usman¹, Shellie Robinson¹, Odunayo Fatunla¹, Eleri Adams², Kyle TS Pattinson^3,4, Simon F Farmer^5,6, Caroline Hartley¹

ORCID iDs: Coen S Zandvoort, Simon F Farmer, Caroline Hartley

Institutions: University of Oxford (United Kingdom); Oxford University Hospitals NHS Trust (United Kingdom); National Hospital for Neurology and Neurosurgery (United Kingdom); UCL Queen Square Institute of Neurology (United Kingdom)

Journal: eLife, volume 14, article RP107081

Dates: published online 12 May 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.7554/elife.107081 · PMID 42118681 · PMCID PMC13167109 · OpenAlex W4413342978

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: EEG (modality), human (organism), developmental (subfield)

Methods: Statistics, Spectral & time-frequency, Connectivity, Machine learning, Preprocessing, Physiology & signal measures

Keywords: EEG, cortical-respiratory coupling, neonate, brain-body interaction, Human

MeSH: Apnea*, Motor Cortex*, Respiration*, Electroencephalography, Female, Humans, Infant, Newborn, Infant, Premature, Male (* major topic)

Journal subjects: Neuroscience

Topic: Neuroscience of respiration and sleep (Endocrine and Autonomic Systems, Neuroscience), according to OpenAlex

Funding: Wellcome Trust (Sir Henry Dale Fellowship 213486/Z/18/Z); Commonwealth Scholarship Commission; University of Oxford; Royal Society (Sir Henry Dale Fellowship 213486/Z/18/Z)

Citations: cited by 1 paper (Europe PMC); 39 references in the paper

Abstract

Respiration is governed by a widespread network of cortical and subcortical structures. This complex communication between the brain and lungs is altered in pathological conditions. Apnoea – the cessation of respiration – is a common condition in infants, particularly those born prematurely. Apnoea in infants is believed to relate to immaturity of brainstem respiratory centres; involvement of the cortex in respiration in infants has yet to be explored. We investigated if there was any evidence for cortical coupling with respiration in newborn humans and whether it relates to apnoea. Using simultaneous electroencephalography (EEG) and impedance pneumography, we investigated interactions between cortical and respiratory activity (known as cortico-respiratory coupling) using phase-amplitude coupling. We show that cortico-respiratory coupling is present in premature and term newborns (104 recordings from 68 infants; 34.5±2.6 weeks postmenstrual age), identifying an interplay between breathing phase and EEG amplitude. We further shed light on the biological meaning by revealing that the strongest coupling occurs during inspiration and that cortical activity precedes respiration, with coupling strongest over frontocentral regions. Whilst our study was limited in spatial resolution, and determining causality is challenging, we believe these findings support the notion that the cortico-respiratory coupling observed here constitutes communication between cortical motor areas and lung effectors. Moreover, we show that cortico-respiratory coupling is negatively correlated with the rate of apnoea, revealing novel insight into this common and potentially life-threatening neonatal pathology.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 15 matches between paragraphs and lines of code.

paediatric_neuroimaging/cortico_respiratory_coupling

License: GPL-3.0

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: d4370e2c53c38e523f1149d12e0b95f7b4572002, 11 May 2026

Languages: MATLAB (27)

Size: 41 files, 27 scripts

Software Heritage: archived

Found in: “Data availability”

Holds: README, license file

Not found: CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Signal Processing Toolbox (5 files), FieldTrip (4 files), EEGLAB (2 files), Image Processing Toolbox (2 files), Statistics and Machine Learning Toolbox (2 files), Brainstorm (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

29 files

compute_cortical_amplitude_modulation.m, MATLAB, 220 lines, 2 matches
compute_pac.m, MATLAB, 104 lines
compute_respiration_metrics.m, MATLAB, 134 lines
compute_time_pac.m, MATLAB, 152 lines, 2 matches
functions/align_using_events.m, MATLAB, 221 lines
functions/align_using_xcorr.m, MATLAB, 132 lines, 1 match
functions/combine_vectors.m, MATLAB, 33 lines
functions/data_to_pac_and_psi.m, MATLAB, 104 lines, 1 match
functions/data_to_time_pac.m, MATLAB, 92 lines
functions/data_to_wavelet_coh.m, MATLAB, 124 lines, 1 match
functions/data_to_wavelet_time_coh.m, MATLAB, 124 lines, 1 match
functions/ft_multiplotTFR.m, MATLAB, 782 lines
functions/ft_prepare_layout.m, MATLAB, 1,621 lines
functions/get_filtered_data.m, MATLAB, 54 lines
functions/load_filter_data.m, MATLAB, 122 lines
functions/load_preprocess_ecg.m, MATLAB, 89 lines
functions/make_ibi_into_spiketrain.m, MATLAB, 41 lines
functions/permutation_testing_cfc.m, MATLAB, 95 lines, 1 match
functions/plot_coherence.m, MATLAB, 185 lines, 1 match
functions/preprocess_ibi.m, MATLAB, 46 lines
functions/read_channels.m, MATLAB, 42 lines
functions/read_events_fromeventfile.m, MATLAB, 55 lines
functions/time_normalise_spectrum.m, MATLAB, 55 lines
prepare_data_for_pac.m, MATLAB, 193 lines
prepare_eeg_data_continuous.m, MATLAB, 39 lines
run_coherence_age_and_apnoea_rate.m, MATLAB, 95 lines, 1 match
stats_and_plot_pac.m, MATLAB, 243 lines
LICENSE, License, 674 lines
README.md, Text, 156 lines

At the source: gitlab.com/paediatric_neuroimaging/cortico_respiratory_coupling

paediatric_neuroimaging/identify_ibi_from_ip

License: none: the authors keep all their rights

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: 72a9ae4bc288257158bdbce4f382b36f1142a387, 1 June 2022

Languages: MATLAB (2)

Size: 6 files, 2 scripts

Software Heritage: archived

Found in: the references

Holds: README

Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: Signal Processing Toolbox (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

3 files

ip_to_interbreath_intervals.m, MATLAB, 1,678 lines, 2 matches
subfunctions/remove_artifactual_longpauses.m, MATLAB, 116 lines, 2 matches
README.md, Text, 14 lines

At the source: gitlab.com/paediatric_neuroimaging/identify_ibi_from_ip

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
29 scripts, each with its path and the digest of its content;
15 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

Ethical restrictions stipulate that data can be shared with other researchers in adherence with ethical constraints (e.g. time limits on data usage). Due to these restrictions, the raw data (raw EEG and vital signs recordings) that support the findings of this study are available on reasonable request. Requests can be made directly to the corresponding author or directed to the Paediatric Neuroimaging Group (University of Oxford) through the institutional will be required to agree to destroy the data within the time limits stipulated in the ethically approved study protocols. Upon agreement, data will be provided within 1 month. Access can only be provided to other researchers, owing to ethical constraints. Summary-level data used to produce figures, and channel- and recording-specific PAC spectra are provided at https://gitlab.com/paediatric_neuroimaging/cortico_respiratory_coupling. The codes for all analyses are available at: https://gitlab.com/paediatric_neuroimaging/cortico_respiratory_coupling (copy archived at Hartley, 2026). The algorithm for inter-breath interval detection is online at: https://gitlab.com/paediatric_neuroimaging/identify_ibi_from_ip (Hartley, 2022).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 8 authors, 5 keywords, 9 MeSH terms, 4 funders, 37 references.

Cite

This paper

Zandvoort, C. S., Usman, F., Robinson, S., Fatunla, O., Adams, E., Pattinson, K. T., Farmer, S. F., & Hartley, C. (2026). Cortical motor activity modulates respiration and reduces apnoea in neonates. eLife, 14, RP107081. https://doi.org/10.7554/elife.107081

BibTeX

@article{zandvoort2026cortical, author = {Zandvoort, Coen S and Usman, Fatima and Robinson, Shellie and Fatunla, Odunayo and Adams, Eleri and Pattinson, Kyle TS and Farmer, Simon F and Hartley, Caroline}, title = {{Cortical motor activity modulates respiration and reduces apnoea in neonates}}, journal = {eLife}, year = {2026}, month = may, volume = {14}, pages = {RP107081}, publisher = {eLife Sciences Publications, Ltd}, issn = {2050-084X}, doi = {10.7554/elife.107081}, url = {https://doi.org/10.7554/elife.107081}, pmid = {42118681}, pmcid = {PMC13167109} }

RIS

TY - JOUR AU - Zandvoort, Coen S AU - Usman, Fatima AU - Robinson, Shellie AU - Fatunla, Odunayo AU - Adams, Eleri AU - Pattinson, Kyle TS AU - Farmer, Simon F AU - Hartley, Caroline TI - Cortical motor activity modulates respiration and reduces apnoea in neonates T2 - eLife J2 - Elife PY - 2026 DA - 2026/05/12 VL - 14 SP - RP107081 SN - 2050-084X PB - eLife Sciences Publications, Ltd DO - 10.7554/elife.107081 UR - https://doi.org/10.7554/elife.107081 LA - en ER -

CSL-JSON

{ "id": "10.7554/elife.107081", "type": "article-journal", "title": "Cortical motor activity modulates respiration and reduces apnoea in neonates", "container-title": "eLife", "author": [ { "family": "Zandvoort", "given": "Coen S" }, { "family": "Usman", "given": "Fatima" }, { "family": "Robinson", "given": "Shellie" }, { "family": "Fatunla", "given": "Odunayo" }, { "family": "Adams", "given": "Eleri" }, { "family": "Pattinson", "given": "Kyle TS" }, { "family": "Farmer", "given": "Simon F" }, { "family": "Hartley", "given": "Caroline" } ], "container-title-short": "Elife", "volume": "14", "page": "RP107081", "DOI": "10.7554/elife.107081", "PMID": "42118681", "PMCID": "PMC13167109", "ISSN": "2050-084X", "publisher": "eLife Sciences Publications, Ltd", "URL": "https://doi.org/10.7554/elife.107081", "language": "en", "issued": { "date-parts": [ [ 2026, 5, 12 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1371/journal.pbio.3003982 [code]: No evidence for modulation of the readiness potential by respiratory phase during natural breathing.
Journal: PLoS biology
In common: FieldTrip, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, EEG, 7 references
[2] doi:10.1038/s41467-026-73828-0 [code]: Human forebrain neural synchronization and entrainment to breathing during wakefulness, sleep, and external mechanical ventilation.
Journal: Nature communications
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 7 references
[3] doi:10.1038/s41467-026-71604-8 [code]: Respiration as a dynamic modulator of sensory sampling.
Journal: Nature communications
In common: FieldTrip, Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 5 references
[4] doi:10.1111/psyp.70271 [code]: Disentangling Respiratory Phase-Dependent and Phase-Independent Components of Anticipatory Cardiac Deceleration.
Journal: Psychophysiology
In common: EEGLAB, FieldTrip, Image Processing Toolbox, 2 other tools, EEG, 3 references
[5] doi:10.1038/s42003-026-10164-5 [code]: The influence of nonlinear resonance on human cortical oscillations.
Journal: Communications biology
In common: Brainstorm, EEGLAB, FieldTrip, 2 other tools, EEG, 1 reference
[6] doi:10.1038/s41467-026-73106-z [code]: Respiratory pauses highlight sleep architecture in mice.
Journal: Nature communications
In common: EEGLAB, FieldTrip, Image Processing Toolbox, 2 other tools, EEG, 2 references
[7] doi:10.1038/s41598-026-49900-6 [code]: Global neural oscillations underlie performance variability and attentional state fluctuations in humans.
Journal: Scientific reports
In common: Brainstorm, EEGLAB, FieldTrip, 3 other tools
[8] doi:10.1162/imag.a.1201 [code]: All spectral frequencies of neural activity reveal semantic representation in the human anterior ventral temporal cortex.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Brainstorm, EEGLAB, FieldTrip, 3 other tools
[9] doi:10.1371/journal.pcbi.1014672 [code]: Robust circular cluster-based statistics for respiration-brain coupling.
Journal: PLoS computational biology
In common: Signal Processing Toolbox, Statistics and Machine Learning Toolbox, 5 references
[10] doi:10.1371/journal.pbio.3003938 [code]: Theta oscillations tag episodic memories for sleep-dependent consolidation.
Journal: PLoS biology
In common: EEGLAB, FieldTrip, Image Processing Toolbox, 2 other tools, EEG, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.