An Electroencephalogram Metric of Temporal Complexity Tracks Psychometric Impairment Caused by Low-dose Nitrous Oxide

This multidose (in randomized order), single-blind, crossover trial took place at the laboratory at Waikato Clinical School, University of Auckland, in July and August 2018. The study protocol was approved by the Health and Disability Ethics Committee, Auckland, New Zealand (reference 16/NTA/93), and was registered with the Australian New Zealand Clinical Trial Registry (registry No. U1111-1181-9722) on December 3, 2018, by S.J.M. The sample size was based on similar studies previously published.^16,22–24 

This study was a prelude to further work investigating EEG effects of gas narcosis in divers, so participants were recruited from that community. Eligible subjects were certified, healthy divers (checked with the Recreational Scuba Training Council [Jacksonville, Florida] screening questionnaire for fitness), aged between 18 and 60 yr, with normal visual acuity, either corrected or uncorrected. Exclusion criteria were the use of recreational drugs, tobacco, psychoactive medication, excessive alcohol (more than 21 standard drinks per week), or over five caffeine-containing beverages a day. All participants provided written informed consent.

Participants abstained from any caffeinated drink on the measurement day and from alcohol at least 24 h before. Participants had at least 6 h of sleep and fasted for 4 h before the measurement.

Participants were seated and breathed from a closed-circuit anesthesia loop (Vital Signs, Mexico) attached to an anesthesia machine (S/5 Aespire, Datex-Ohmeda, USA). A mouthpiece and disposable anesthetic antibacterial filter (Ultipor 25, Pall, USA) were replaced for each participant. The nose was occluded with a nose clip. The inspired fraction of oxygen, end-tidal pressure of carbon dioxide (Petco₂), and end-tidal percentage of nitrous oxide were continuously sampled from the mouthpiece filter and were recorded every minute. Oxygen saturation and breathing frequency were monitored for participant safety.

Every participant started with a baseline measurement while breathing 50% oxygen (balance nitrogen) on the circuit. The subjects then breathed a titrated nitrous oxide mixture (balance oxygen) to achieve an end-tidal level of 20%, 30%, or 40% N₂O, followed by 20 min of breathing air between each nitrous oxide exposure. Participants were blinded to the dose of nitrous oxide, which was administered in random order of doses dictated by a ticket drawn by the researcher when the participant arrived (order: 20-30-40%; 20-40-30%; 30-20-40%; 30-40-20%; 40-20-30%; or 40-30-20%).

Every exposure started with a 3- to 5-min wash-in period, during which the inspired fraction of nitrous oxide was manually adjusted to establish and maintain the desired end-tidal level. Then a set of measurements was undertaken consisting of EEG recording over 1 min with eyes open and 1 min with eyes closed, completion of psychometric tests, and a pupillometry measurement (pupillometry results were published elsewhere²⁵ ). Finally, the 1-min eyes-open and -closed EEG recordings were repeated. At the end of each exposure, nitrous oxide was washed out using oxygen with a flow of 6 l · min⁻¹ and was followed by a 20-min air breathing rest period. A final measurement set during 50% oxygen (balance nitrogen) breathing was recorded 20 min after the last nitrous oxide exposure (fig. 1).

The EEG was recorded using a portable active electrode 32-channel system (ActiveTwo, BioSemi, The Netherlands). The electrodes were placed in a sized cap divided over the scalp based on the international 10-10 system.²⁶  Two additional electrodes were placed under the eyes to record the electrooculogram to filter ocular artefacts. The offset (impedance equivalent for active systems) was checked for all electrodes, and electrode placement and gelling (SignaGel, Parker Laboratories, USA) were adjusted if the offset was above 25 μV. All signals were recorded at a sample rate of 1,024 Hz in the BioSemi Data Format (BDF) format on a laptop (Macbook Pro, Apple Inc., USA) using ActiView software (BioSemi, The Netherlands) for offline analysis. EEG was recorded continuously from the wash-in period to the last measurement for each exposure.

Three conditions were identified for analysis: tasked (while doing psychometric tests), resting-state eyes open, and resting-state eyes closed. EEG recordings taken after the psychometric tests and pupillometry (fig. 1) were used as more stable nitrous oxide levels were achieved, and fewer artefacts were present in these recordings.

Two psychometric tests were administered on a 9.7-inch tablet computer (Galaxy Tab Active2, Samsung, South Korea). The test administration program (PenScreenSix version 2.1, Mobile Cognition Ltd., United Kingdom) stored average reaction time and number of errors made (accuracy) in each test. Two tests were selected that have shown sensitivity to narcotic effects because they test for higher-order cognitive functions.²⁷ 

The shape-recognition test is supposed to measure the effect on short-term memory.²⁸  However, we found erratic results unsuitable for consistent detection of impairment, and this test was not included in our analysis.

The serial sevens test measures information processing, in particular mathematics, memory, and decision-making.²⁹  The participant had to decide whether the current three-digit number was the previously shown number minus seven (yes or no). The test consisted of a series of 16 descending numbers (15 questions). The maximum allowed response time per question was 10 s.

Participants attended a training session for the psychometric tests to mitigate a learning effect during the actual measurements. The tests were practiced until the results were stable.

The EEG data were cleaned using the FieldTrip toolbox (version c6d58e9).³⁰  The data for each condition and exposure were cut out of the continuous recording, rereferenced to the average, demeaned, detrended, and resampled to 256 Hz. Line (including higher harmonic) and low-frequency noise (less than 1 Hz) were filtered out. Next, independent component analysis was used to select out noise components from eye blinks, high- frequency noise, nonphysiologic noise, and bad channels. An algorithm was used to advise on the manual selection of components. The data were cut in 2-s epochs and manually inspected for remaining artefacts, with an algorithm indicating bad segments for remaining eye blinks (correlation with the electrooculogram channels) and muscle artefacts (based on high-frequency content of 105 to 120 Hz). Whole epochs were discarded if they were marked as containing artefacts. The remaining epochs for that condition and exposure were stored for further analysis. See appendix 1 for the script.

Frequency power was estimated using a multitaper Fourier transformation implemented in the FieldTrip toolbox using a Hanning window from 1 to 30 Hz in 1-Hz increments. Average per frequency band was calculated for delta (1 to 4 Hz), theta (4 to 8 Hz), alpha (8 to 14 Hz), and beta (14 to 30 Hz) bands during each exposure.

The sequence to obtain our novel complexity metric is shown in figures 2 and 3. In brief, it quantifies the variability in repeated short motifs of EEG signal over a time scale of seconds to tens of seconds (fig. 2),¹⁸ i.e., the recurrence properties of the EEG signal. Recurrence plots are widely used in physiology,¹⁹  but when applied to EEG analysis, they provide a graphical indication of periods when a section of the EEG signal is similar to subsequent sections. The presence of a long diagonal line in the recurrence plot indicates a slowly evolving similarity in EEG pattern (i.e., to produce a diagonal line, the motif at time = 1 is similar to that at time = 2, and the motif at time = 3 is similar to that of time = 2). Blocks of high correlation (yellow in fig. 3B) are stable periods. The standard recurrence plot summary statistic (Shannon entropy of the diagonal line length) can be seen as a measure of the variation in the duration of metastable states, i.e., the temporal complexity of the signal for each channel. A high entropy suggests the presence of a mixture of slowly evolving metastable brain states, coexistent with short-lived brain states. As seen in the histograms of figure 3C, a low entropy is indicative of predominantly short-lived recurrences. To condense this into a single spatial summary statistic for the whole scalp, we found the median complexity of a subset of high-complexity channels. This subset of channels was found to roughly overlie the default mode network. This is in agreement with previous work showing that the default mode network has high complexity based on its critical functional role in resting-state networks.³¹  The following steps give the metric:

1. To obtain the motifs, the EEG signals were split into samples of 0.25-s duration. For each channel, the absolute values of the Pearson correlations between all of these 0.25-s samples were calculated (fig. 3A).¹⁸ 

2. The median value of Pearson correlation over all channels and cross-correlations of the baseline exposure was taken as the threshold to produce a dichotomized matrix for the recurrence analysis (fig. 3B).

3. The Shannon entropy of the probability distribution of the diagonal line lengths³²  for each channel was calculated to capture the variability in temporal changes (fig. 3C).

4. A region of interest was defined as the electrodes with a complexity value larger than 1 at baseline. As described above, based on the regional distribution of our results, we named the spatial distribution of the temporal complexity metric, “default-mode-network complexity.” The metric was calculated for each exposure and condition (tasked, eyes open, eyes closed) as the median of the complexity of the electrodes in the set region of interest (red-encircled area in fig. 3D).

The results were robust to a range of differing trialed values for the sample length, dichotomization threshold, and region-of-interest threshold (results not shown). See appendix 2 for the script. Regional complexity was calculated by taking the mean of the Shannon entropy of the channels of that region.

A combined psychometric-impairment metric (scaled 0 to 1) was calculated using both the mean reaction time and error rate to counteract the speed-accuracy trade-off (equation 1).³³ 

All data were analyzed with Matlab version 2018a (Mathworks, USA) except the carbon dioxide data, which were analyzed with SPSS version 25 (IBM, USA). The idea of EEG data analysis using temporal complexity was preplanned, based on previous published work.^20,34,35  However, the specific values of various parameters in the algorithm and the default mode network region of interest were data-driven, as part of the exploratory study. The comparison between the EEG metric and psychometric test score was predefined. All outcome measures were tested for normality and subsequently characterized by their mean and SD. Comparisons of regional complexity were done using a two-tailed paired t test and reported as mean difference and 95% CI. P values were regarded as significant at P < 0.05, with Bonferroni correction for multiple comparisons.

For all exposures and individuals, a linear mixed-effects model was used to calculate the relationship between the EEG default-mode-network complexity metric (DMNcomplexity) and the psychometric test metric (S77) using equation 2.

The between-participant variation was included as a random effect, with a random intercept and slope. This design can capture variation in baseline complexity and dose response for each subject. The models for the three conditions were compared using a likelihood ratio test with 1,000 simulations.

To research the possible influence of hypercapnia, carbon dioxide levels were analyzed. The Petco₂ value at the start of the psychometric test was added as a fixed-effect parameter to the mixed-effects model for the tasked condition (EEG recording during the psychometric test). This model was compared to the model without carbon dioxide level (same method as above).

The influence of the power in each frequency band (delta [1 to 4 Hz], theta [4 to 8 Hz], alpha [8 to 14 Hz], beta [14 to 30 Hz], and gamma [30 to 45 Hz]) on the final model was also investigated. For this subanalysis, the EEG data of the eyes-open condition were bandpass-filtered with a Butterworth filter (order 4) in both directions using a Hamming window for each frequency band.

A receiver operating characteristic analysis was used to compare the default-mode-network complexity metric against the serial sevens test results. The serial sevens data were dichotomized to “impaired” (greater than 0.16) and “not impaired” (less than or equal to 0.16), based on the upper value of the interquartile range (fig. 4) of the baseline exposure.

The 12 participants (8 male), aged between 23 and 55 yr (mean. 36 yr), had an average body mass index of 26.3 kg/m² (22.7 to 31.3). On average, 1.7 of 30, 1.3 of 30, and 19.5 of 108 (eyes open, eyes closed, and tasked, respectively) of the 2-s EEG samples were removed due to artefacts. The latter were mostly muscle artefacts. EEG frequency analysis shows an increase in power in the beta band, with a relative decrease in the delta band, most prominent in the frontal region and most visible in the 20% and 40% exposures (Supplemental Digital Content 1, https://links.lww.com/ALN/C512, shows the frequency band analysis). This is similar to previous studies.^10–14  In general, the serial sevens test results worsened (fig. 4A) and the complexity metric decreased (fig. 4B) at increased end-tidal concentrations of nitrous oxide, particularly once the 40% concentration was reached. End-tidal nitrous oxide concentrations exhibited a moderate correlation with both the serial sevens scores and the EEG default-mode-network complexity (r = 0.50, P < 0.001 and r = −0.55, P < 0.001, respectively).

At baseline (fig. 5), regional temporal complexity was greater in the midline (1.33 ± 0.14 bits), which approximately overlies the default mode network, than in channels overlying lateral temporal brain regions (0.81 ± 0.12 bits; mean difference, 0.52 [0.39 to 0.65], P < 0.001). Increasing nitrous oxide exposure caused a decrease in brain complexity predominantly in the midline (baseline vs. 40% N₂O: mean difference, 0.20 bits [0.09 to 0.31], P = 0.002) and prefrontal (baseline vs. 40% N₂O: mean difference, 0.17 bits [0.08 to 0.27], P = 0.002) regions, whereas complexity in the lateral temporal region did not change significantly (baseline vs. 40% N₂O: mean difference, 0.14 bits [−0.03 to 0.30], P = 0.100). These changes were reversed after cessation of nitrous oxide.

At baseline, the EEG default-mode-network complexity metric ranged from 1.13 to 1.38 bits between subjects. The effects of the different doses of nitrous oxide were variable between subjects. For some participants, the ranking order of the nitrous oxide concentrations did not correlate with the ranking order of both the psychometric tests (y axis) and EEG default-mode-network complexity (x axis; fig. 6). There was a significant difference between the models for only the eyes-closed versus tasked conditions (P = 0.003), but the log likelihoods of the models for the three conditions were very similar (52.3, 47.4, and 42.7). Therefore, figure 7 shows the results of the model based on the eyes-open data.

At an individual level, for each subject there was a significant linear trend between the EEG default-mode-network complexity metric and cognitive decline (fig. 6 blue lines and table 1). The goodness of fit of the model was confirmed by a good agreement with fitted versus measured values (fig. 7A, r² = 0.67) and lack of autocorrelation in the residuals (fig. 7B). On average, psychometric impairment increased 10% with every 38% decrease in EEG default-mode-network complexity. The receiver operating characteristic curve of the complexity metric versus the serial sevens test showed an area under the curve of 0.72 (95% CI, 0.59 to 0.85; P < 0.001; fig. 8).

Petco₂ was within normal ranges at group level (5.5 ± 0.6 kPa [41 ± 4.5 mmHg]), and there was no significant difference between exposures. However, some participants were above normal (greater than 5.7 kPa [43 mmHg]) at the start of the psychometric test. The addition of Petco₂ values as a fixed effect in the model did not change the results (P = 0.460). The frequency subanalysis revealed that there was no specific frequency band driving the default-mode-network complexity (results not shown).

In this exploratory study, our novel complexity metric tracked psychometric impairment during low-dose nitrous oxide exposure. The EEG default-mode-network complexity metric is a simple and fast algorithm that incorporates two main concepts: cross-correlation–based recurrence quantification analysis³²  and specific regional changes.³⁶ 

For each channel, we calculated the complexity of that signal by determining the correlation between short successive EEG samples. The concept behind this is that if there is a high correlation over time, the EEG signal is repeating itself, indicating the presence of prolonged or similar metastable brain states. Variation in runs of repetition of EEG patterns is seen as an increase in the variability of diagonal line lengths in the recurrence plots and consequently an increase in the Shannon entropy of these plots (also seen as an increase in width of the histograms; fig. 3C). The increased entropy of the recurrences is thus primarily driven by many periods of prolonged (more than approximately 1,000 ms) EEG similarity, which widen the recurrence histograms. The neurophysiologic correlates of this pattern are not well defined, but it is well established that conscious perceptions are associated with ignition of sustained (more than 250 ms) neural activity³⁷  and prolonged (more than approximately 1,000 ms) complex responses to transcranial magnetic stimulation.³⁸  Although our study does not look at the transition to unconsciousness, we can speculate that inability of the brain to sustain recurrence patterns is a signature of loss of diversity of metastable states, which seems to be an indication of impairment of higher-level function.

This first part of our analysis method utilizes methodologies that originate from recurrence quantification analysis, but in our study, these methodologies were applied to cross-correlation plots.¹⁸  These methods quantify nonlinear dynamic phenomena.¹⁹  Based on the global neuronal workspace hypothesis, normal cognitive performance could be associated with recurrence loops for two reasons.³⁷  First, recurrence loops can help cortical processors sustain a signal, and hence the information can be maintained in the working memory. Second, by recurrent excitation, a signal can be amplified, helping information sharing between cortical processors. Hence, a certain amount of recurrence is needed to sustain normal cognitive functioning. However, it is unclear whether the EEG recurrence is driven by local or global reverberant networks. In the neuroscience field, recurrence analysis has been utilized in epilepsy detection,³⁹  sleep-stage recognition,^40,41  and depth of anesthesia detection.^35,42,43  However, the depth-of-anesthesia articles all report results with GABAergic agents, whereas our study used an N-methyl-d-aspartate antagonist.

The other observation from this work is the marked difference in regional baseline complexity and the regional effect of the nitrous oxide. Prefrontal and midline regional complexity reduced under nitrous oxide exposure, whereas complexity of the lateral fronto/temporal/parietal regions showed relatively minor change. The serial sevens test requires coordination of several cognitive tasks, e.g., mathematics, memory, and decision-making,²⁹  which were impaired by nitrous oxide. The regional differences correlate well with the concept of the fragmentation of selfhood,³⁶  where higher-order thinking (affected in our study) is primarily located at the (pre)frontal region, whereas the sentience and salience (located in the temporal/insular region) are less affected by low-dose nitrous oxide. Numerous studies have supported this division of functions in the brain, all concluding that higher cognitive functions are predominantly located in the (pre)frontal region and with their connections to posterior medial regions.^44,45  Similar reductions in frontal connectivity have been found during subanesthetic ketamine exposure.⁴⁶ 

A recent study showed a significant decrease in EEG network connectivity (lagged phase coherence) induced by 50% inhaled nitrous oxide.⁴⁷  This disruptive effect of relatively high-dose nitrous oxide on the flow of information between brain regions is consistent with the default mode network regional decrease in complexity found in our study. The nitrous oxide–impaired cognitive performance was strongly associated with loss of complexity in the medial default mode network regions (the so-called “hot zone”).³⁸  The exact role of these brain regions (precuneus, posterior, and anterior cingulate cortices) is still intensely debated,⁴⁸  but seems to be more closely aligned with overall higher-order brain integration and coordination⁴⁹  rather than specific tasks. In particular, activity is associated with semantic processing and memory retrieval, which are necessary for the serial sevens test.⁵⁰ 

Our study had some strengths. It was performed in a laboratory where the environment was well controlled to optimize the EEG recording, with full-scalp 32 channels. By controlling the end-tidal concentrations of nitrous oxide, we were able to control the exposure to a narcotic agent precisely at three distinct levels, without the interference of other drugs. However, the psychometric tests showed a high level of personal variability in participants’ responses to nitrous oxide. The variability could be explained by a carryover or learning effect between the exposures, although we did incorporate a 20-min pause between exposures and randomized the order of exposures to minimize this effect. We measured end-tidal carbon dioxide and found that an increased Petco₂ above a normal value of 5.7 kPa (43 mmHg), which was found in some participants/exposures, did not contribute to the narcotic effects found, similar to the results of Foster and Liley.¹¹ 

Our study also had some weaknesses. First, the participant group was relatively small. In mitigation, it was reassuring that the EEG default-mode-network complexity metric consistently decreased in every subject as psychometric impairment increased. Second, our exploratory analysis was directed toward a time-domain analysis. We did not use source localization to refine the locations of temporal complexity further; hence, we chose to use general descriptors of locations. This increased the speed of analysis, kept the analysis method simple, and avoided the assumptions underlying all source localization methods. We also acknowledge that our limited spatial resolution means that we have used the term “default-mode-network complexity” as a descriptive shorthand rather than being able to establish its full anatomical detail.

In conclusion, our novel quantitative EEG default-mode-network complexity metric based on temporal complexity is sensitive to psychometric impairment caused by low-dose nitrous oxide. However, further research is needed to evaluate its functionality, because this was an exploratory study. If robust, this algorithm may be useful in quantifying and studying mild narcosis in situations where patients are less self-aware due to the dissociative effects of nitrous oxide or other non-GABAergic sedatives.

The authors are grateful to all participants in this study. Furthermore, the authors acknowledge Jonathan Termaat, Dip.R.N., Gay Mans, N.Z.R.G.O.N., Amy Gaskell, M.B.Ch.B., F.A.N.Z.C.A. (all Department of Anaesthesia, Waikato Hospital, Hamilton, New Zealand), Rebecca Pullon, Ph.D., and Marta Seretny, M.D., M.P.H., Ph.D., F.R.C.A. (both Department of Anaesthesiology, University of Auckland, Auckland, New Zealand), for their support during the data collection.

This study was supported by funding from the Office for Naval Research Global (Tokyo, Japan), U.S. Navy grant No. N62909-18-1-2007.

The authors declare no competing interests.

Full protocol available at: x.vrijdag@auckland.ac.nz. Raw data available at: x.vrijdag@auckland.ac.nz.

close all

clear

clc

epochtype= ‘closed end’;%’math’;% ‘open end’;

[datasetselection,datapath] = data_selector;

for i=1:length(datasetselection) %all selected data folders

filepath=[datapath datasetselection{i} filesep];

filelist=dir([filepath ‘*.bdf’]);%find all.bdf files

%creating folder for cleaned files to be stored

newfilepath=[filepath filesep ‘cleaned’ filesep];

mkdir(newfilepath);

for j=2:length({filelist.name}) % all BDF files in the folder

%define epoch trial data

cfg = [ ];

cfg.dataset = [filepath filelist(j).name];

cfg.trialfun = ‘EEG_trialfunc_narcosis’; % my function to read epochtimes from excel file

cfg.trialdef.epochtype = epochtype;

cfg.trialdef.pretrial = 1; % add sec of data before trial

cfg.trialdef.posttrial = 1; % add sec of data after trial

cfg = ft_definetrial(cfg);

%read dataset based on trial definition

cfg.channel = 1:34;

cfg.reref = ‘yes’;

cfg.refchannel = ‘EEG’; % average of all EEG channels

cfg.demean = ‘yes’;

cfg.detrend = ‘yes’;

data = ft_preprocessing(cfg);

% downsample data

cfg = [ ];

cfg.resamplefs = 256; %Hz

data = ft_resampledata(cfg, data);

% removing line and low freq noise

cfg = [ ];

cfg.channel = ‘all’;

cfg.bsfilter = ‘yes’; % band-stop method

cfg.bsfreq = [49 51];

cfg.hpfilter = ‘yes’;

cfg.hpfreq = 1; % removing all activity below 1 Hz

data = ft_preprocessing(cfg,data);

% clean high harmonic of line noise

cfg=[ ];

cfg.channel = channels;

cfg.bsfilter = ‘yes’; % band-stop method

cfg.bsfreq = [99 101];

cfg.demean = ‘yes’;

cfg.detrend = ‘yes’;

data = ft_preprocessing(cfg,cleandata);

% calculate ICA components to project out eye and cardiac artefacts

cfg = [ ];

cfg.method = ‘runica’;

comp = ft_componentanalysis(cfg, data);

cfg = [ ];

cfg.layout = ‘biosemi32.lay’;

cfg.viewmode = ‘component’;

cfg.ylim = [-0.0002 0.0002];

cfg.zlim = ‘maxmin’;

cfg.compscale = ‘local’; % scale each component separately

cfg.blocksize = 62;

cfg.artifactalpha = 0.8;

cfg.position = [0 0 3000 2000];

ft_databrowser(cfg, comp);

% calculating properties to suggest, select and reject bad components

EOG_chans = [33 34];

[bad_comps,reasons] = ICAbadICdetection(data,comp, EOG_chans);

for k=1:length(bad_comps)

ft_info(‘Suggested bad component: %d because of %s’,bad_comps(k),reasons{k})

end

ic.artifact = input(‘ICs to reject (i.e. [8] or 0 for all suggested bad components): ‘);

if ic.artifact==0

ic.artifact=bad_comps;

end

cfg = [ ];

cfg.component = ic.artifact;

data = ft_rejectcomponent(cfg, comp, data);

% remove padding of 1 sec before and after and correct time axes

cfg = [ ];

cfg.toilim = [1 data.time{1,1}(end)-1];

data = ft_redefinetrial(cfg, data);

cfg = [ ];

cfg.offset = -data.fsample;

data = ft_redefinetrial(cfg, data);

%cut in 2 second “trials”

cfg = [ ];

cfg.length = 2;

cfg.overlap = 0;

data_segmented = ft_redefinetrial(cfg, data);

% detect artefacts

cfg=[ ];

cfg.artfctdef.eog.channel = {‘EXG1’, ‘EXG2’};

cfg.artfctdef.eog.cutoff = 7;

cfg.artfctdef.muscle.channel = ‘EEG’;

cfg.artfctdef.muscle.bpfreq = [105 120];

cfg.artfctdef.muscle.cutoff = 10;

[cfg, ~] = ft_artifact_eog(cfg, data_segmented);

[cfg, ~] = ft_artifact_muscle(cfg, data_segmented);

%visually inspect data and defined artefacts

cfg.channel = ‘all’;

cfg.viewmode = ‘vertical’;

cfg.ylim = [-40 40];

cfg.position = [0 0 3000 2000];

cfg.blocksize = 10; % time window to browse

cfg.artifactalpha = 0.8; % this make the colors less transparent and thus more vibrant

cfg = ft_databrowser(cfg, data);

%rejecting trials with artefacts

cfg.artftdef.reject = ‘complete’;

cleandata=ft_rejectartifact(cfg,data_segmented);

%store cleaned data

[~,newfilename,~]=fileparts(filelist(j).name);

newfilename=[‘CLEAN’ newfilename ‘ ‘ epochtype ‘.mat’];

save(strcat(newfilepath,newfilename),’cleandata’,’-v7.3’)

end

end

function [bad_comps,reasons] = ICAbadICdetection(data,comp, EOG_chans)

% ICA bad components detection function

% detection scripts are from the FASTER toolbox

% utilising functions from the EEGlab toolbox

list_properties = component_properties(data,comp,EOG_chans);

rejection_options.measure=ones(1,size(list_properties,2)); %use all properties

rejection_options.z=3*ones(1,size(list_properties,2));% z-value of 3 for each property

[lengths,all_l] = min_z(list_properties,rejection_options);

bad_comps=find(lengths);

properties={‘High freq noise’ ‘Kurtosis’ ‘Hurst’ ‘Eyeblink’};

for i=1:length(bad_comps)

reasons{i}=cell2str(properties(all_l(bad_comps(i),:)));

end

end

function string1 = cell2str(cellarray)

string1=[ ];

for i=1:length(cellarray)

string1=[string1 cellarray{i}];

if length(cellarray)>i

string1=[string1 ‘ & ‘];

end

end

end

function list_properties = component_properties (data,comp, blink_chans)

list_properties = zeros(size(comp.trial,1),4); %This 4 corresponds to number of properties.

for u=1:size(comp.trial{1},1)

measure = 1;

% TEMPORAL PROPERTIES

% Median gradient value, for high frequency artefacts

list_properties(u,measure) = median(diff(comp.trial{1}(u,:)));

measure = measure + 1;

% SPATIAL PROPERTIES

% Kurtosis of spatial map (if v peaky, i.e. one or two points high

% and everywhere else low, then it’s probably noise on a single

% channel)

list_properties(u,measure) = kurt(comp.topo(:,u));

measure = measure + 1;

% OTHER PROPERTIES

% Hurst exponent - detects nonphysiological signals

list_properties(u,measure) = hurst_exponent(comp.trial{1}(u,:));

measure = measure + 1;

% Eyeblink correlations

if (exist(‘blink_chans’,’var’) && ~isempty(blink_chans))

for v = 1:length(blink_chans)

if ~(max(data.trial{1}(blink_chans(v),:))==0 && min(data.trial{1}(blink_chans(v),:))==0)

f = corrcoef(comp.trial{1}(u,:),data.trial{1}(blink_chans(v),:));

x(v) = abs(f(1,2));

else

x(v) = v;

end

end

list_properties(u,measure) = max(x);

measure = measure + 1;

end

end

for u = 1:size(list_properties,2)

list_properties(isnan(list_properties(:,u)),u)=nanmean(list_properties(:,u));

list_properties(:,u) = list_properties(:,u) - median(list_properties(:,u));

end

end

function [lengths,all_l] = min_z(list_properties,rejection_options)

rejection_options.measure=logical(rejection_options.measure);

zs=list_properties-repmat(mean(list_properties,1),size(list_properties,1),1);

zs=zs./repmat(std(zs,[],1),size(list_properties,1),1);

zs(isnan(zs))=0;

all_l = abs(zs) > repmat(rejection_options.z,size(list_properties,1),1);

lengths = any(all_l(:,rejection_options.measure),2);

end

close all; clear; clc

%settings

Ts=0.25; %duration of sample in sec

RQAselect=3; % entropy

Rnpercentile=50; % percentile for threshold

channels= 1:32;

%select dataset

%for model figure select participant 1

[datasetselection,datapath] = data_selector;

epochtype= ‘open end’; % ‘closed end’;% ‘math’; % ‘open end’; %

exposures={‘baseline’ ‘20%’ ‘30%’ ‘40%’ ‘final’};

for k=1:length(datasetselection) %all selected data folders

filepath=[datapath datasetselection{k} filesep ‘cleaned’ filesep];

RnThreshold=[ ];

for l=1:length(exposures) % loop over all exposures

% load datafile

filename=[‘CLEAN’ datasetselection{k} ‘ ‘ exposures{l} ‘ ‘ epochtype ‘.mat’];

load([filepath filename])

% split previously defined trials (2 second) into smaller samples

a=1;

samples=[ ];

N=Ts*data.fsample;

for i=1:length(data.trial)

for j=1:(floor(length(data.trial{1})/N))

samples(a,:,:)=data.trial{i}(:,1+(j-1)*N:N*j);

a=a+1;

end

end

%DIMORD samples = samples * channels * time

% calculate the Pierson correlation between the samples for all channels

Rn=[ ];

for i=1:size(samples,2)

Rn(i,:,:)=corr(squeeze(samples(:,i,:))’);

end

% convert the correlation values

Rn=abs(Rn);

% set threshold at baseline exposure for all exposures for each defined percentile

if l==1

RnT=Rn;

RnT(RnT==0)=NaN;

RnThreshold=prctile(reshape(RnT,1,[ ]),Rnpercentile);

end

% Convert to logical values based on threshold

RnLogic = (Rn > RnThreshold);

%calculate the defined RQA metrics for each channel

for i=1:size(RnLogic,1)

RQresults(k,l,i) = Recu_RQA(squeeze(RnLogic (i,:,:)),1,RQAselect);

end

% DIMORD RQresults participants, exposures, channels

% calculate DMN complexity: median of the region of

% interest which are the electrodes that have entropy >1 at baseline

if l==1

DMN=find(RQresults(k,l,:)>1);

end

DMNcomplexity(k,l)=median(RQresults(k,l,DMN));

% DIMORD participants, exposures

end

end

function [RQA] = Recu_RQA(RP,I,metrics)

% Recurrence quantification analysis of recurrence plots

% RP: the Recurrence Plot

% I: the indication marks (I=0 RP is the symmetry matrix

%    I=1 RP is the asymmetry matrix)

%outputs in RQA

% 1 RR: Recurrence rate RR, The percentage of recurrence points in an RP

% Corresponds to the correlation sum;

% 2 DET: Determinism DET, The percentage of recurrence points which form

% diagonal lines

% 3 ENTR: Entropy ENTR, The Shannon entropy of the probability distribution of the diagonal

% line lengths p(l)

% 4 L: Averaged diagonal line length L, The average length of the diagonal lines

% 5 LAM: Laminnarity, The percentage of recurrence points which form vertical lines

% 6 TT: Trapping time, The average length of the vertical lines

% 7 Clust:global clustering coefficient, Network analysis

% 8 Trans:Transitivity, Network analysis

% If you need these codes that implement critical functions with (fast) C code, please visit my website:

% http://www.escience.cn/people/gxouyang/Tools.html

% revise time: May 5 2014, Ouyang,Gaoxiang

% Email: ouyang@bnu.edu.cn

%% minimal length of diagonal and vertical line structures

Lmin=2; %diagonal

Vmin=2; %vertical

if nargin < 2

I=0;

end

%% calculate diagonals

N1=size(RP,1);

Yout=zeros(1,N1);

for k=2:N1

On=1;

while On<=N1 + 1-k

if RP(On,k+On-1)==1

A=1;off=0;

while off==0 & On~=N1 + 1-k

if RP(On+1,k+On)==1

A=A+1;On=On+1;

else

off=1;

end

end

Yout(A)=Yout(A)+1;

end

On=On+1;

end

end

if I==0

S=2*Yout;

end

if I==1

RP=RP’;

for k=2:N1

On=1;

while On<=N1 + 1-k

if RP(On,k+On-1)==1

A=1;off=0;

while off==0 & On~=N1 + 1-k

if RP(On+1,k+On)==1

A=A+1;On=On+1;

else

off=1;

end

end

Yout(A)=Yout(A)+1;

end

On=On+1;

end

end

S=Yout;

end

%% calculate the recurrence rate (RR)

SR=0;

for i=1:N1

SR=SR+i*S(i);

end

RR=SR/(N1*(N1-1));

%% calculate the determinism (%DET)

if SR==0

DET=0;

else

DET=(SR-sum(S(1:Lmin-1)))/SR;

end

%% calculate the ENTR = entropy (ENTR)

pp=S/sum(S);

entropy=0;

F=find(S(Lmin:end));

l=length(F);

if l==0

ENTR=0;

else

F=F+Lmin-1;

ENTR=-sum(pp(F).*log(pp(F)));

end

%% calculate Averaged diagonal line length (L)

L=(SR-sum([1:Lmin-1].*S(1:Lmin-1)))/sum(S(Lmin:end));

%% calculate Laminarity (LAM) and Trapping time (TT) (Marwan et al. 2002)

[~, d, ~]=tt(RP);

d(d<Vmin)=[];

if sum(RP(:))>0

LAM=sum(d)/sum(sum(RP));

else

LAM=NaN;

end

TT=mean(d);

%% network measures

% global clustering coefficient (Marwan et al., 2009)

kv = sum(RP,1); % degree of nodes

Clust = nanmean(diag(double(RP)*double(RP)*double(RP))’ ./ (kv.* (kv-1)));

% transitivity

denom = sum(sum(double(RP) * double(RP)));

Trans = trace(double(RP)*double(RP)*double(RP))/denom;

aa=1;

if ismember(1,metrics);RQA(aa)=RR;aa=aa+1;end

if ismember(2,metrics);RQA(aa)=DET;aa=aa+1;end

if ismember(3,metrics);RQA(aa)=ENTR;aa=aa+1;end

if ismember(4,metrics);RQA(aa)=L;aa=aa+1;end

if ismember(5,metrics);RQA(aa)=LAM;aa=aa+1;end

if ismember(6,metrics);RQA(aa)=TT;aa=aa+1;end

if ismember(7,metrics);RQA(aa)=Clust;aa=aa+1;end

if ismember(8,metrics);RQA(aa)=Trans;end