OSCR

Divergent periodic and aperiodic EEG signatures of Propofol versus sevoflurane anesthesia: a comparative neurophysiological study.

Code ↔ Paper

1 match 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 1 match
  1. [1] § Methods › Multivariate statistical framework and interpretability analysis ↔ LocalResources/ConsciousnessClassifier.py, lines 323–429 · score 0.54 · cross validation, fold, probabilities, score, predicted, model

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 authors' code

Python · 824 lines · 32 KB · no license · 1 match

  1. """
  2. OOP code for consciousness classification.
  4. JHA
  6. sklearn needs to be version 0.21.2
  7. """
  8. import os
  9. import warnings
  10. import pickle
  11. from time import time
  12. from itertools import combinations
  13. import numpy as np
  14. from sklearn.metrics import roc_curve, auc, accuracy_score
  15. from sklearn import decomposition, discriminant_analysis, svm, linear_model
  16. import matplotlib.pyplot as plt
  17. import pandas as pd
  18. from hmmlearn import hmm as hmm_model
  20. # ignore FutureWarning errors which clutter the outputs
  21. warnings.simplefilter(action='ignore', category=FutureWarning)
  23. ##########################
  24. ###### GLOBAL VARS #######
  25. ##########################
  26. # CNN features are the only ones not being generated online so we need file paths to where its stored on disk
  27. FP_BTLNCKS_VOLUNTEER = 'Data/Volunteer_CNN/btlnck_df.csv'
  29. n_pcs_cnn = 10 #number of principal components for PCA of CNN bottlenecks
  31. if not os.path.exists(FP_BTLNCKS_VOLUNTEER):
  32. warnings.warn(f"filepath to CNN bottlenecks for volunteer data doesn't exist\nFP passed: {FP_BTLNCKS_VOLUNTEER}")
  34. class ConsciousnessClassifier(object):
  35. """Consciousness classifier object that can use a variety of algorithms
  36. to calculate the probability of consciousness from a spectogram of Fp1.
  37. Has been tested on Volunteer data from Purdon et al 2013 and on OR data.
  38. """
  40. def __init__(self, name=None):
  41. """
  42. An object for writing, fitting, and assessing models
  43. """
  44. self.name = name
  45. self.datasets = {}
  46. self.models = {} # option for multiple models
  48. def load_dataset(self, path, names, dsetname, subjects='OR'):
  49. """loads dataset as a dictionary and gets added to the self.datasets attribute
  50. NOTE: This object is built to have separate data sets for train and test (internal validation or external validation)
  52. Parameters
  53. ----------
  54. path : str
  55. directory path to folder containing by case preprocessed data
  56. e.g.~/Dropbox (Partners HealthCare)/HumanSignalsData/projects/consciousness_classifier/volunteer/by_case/
  57. names : list(str)
  58. list of case ids to be loaded from path
  59. dsetname : str
  60. key to add the resulting data_dict to the self.datasets atribute
  61. subjects : str, optional
  62. category of subjects in dataset, either 'OR' or 'volunteer' by default 'OR'
  64. Raises
  65. ------
  66. SystemExit
  67. if something other than 'OR; or 'volunteer' is passed as subjects
  68. """
  69. data_dicts = []
  70. if subjects == 'volunteer':
  71. for tdn in names:
  72. dd = load_volunteer_todict(path, tdn)
  73. if dd is not None:
  74. data_dicts.append(dd)
  76. elif subjects == 'OR':
  77. for tdn in names:
  78. dd = load_OR_todict(path, tdn)
  79. if dd is not None:
  80. data_dicts.append(dd)
  82. # assemble info from each data dict into a dataset
  83. caseid = []
  84. sdata = []
  85. times = []
  86. labels = []
  87. egq = []
  88. for dd in data_dicts:
  89. caseid = np.hstack([caseid,np.full(len(dd['t']),dd['name'])])
  90. dd_sdb = dd['Sdb']
  92. # check for quality
  93. if subjects is 'OR':
  94. dd_egq = dd['egq']
  95. elif subjects is 'volunteer':
  96. dd['egq'] = np.ones(len(dd['l']))
  97. dd_egq = dd['egq']
  98. else:
  99. raise SystemExit(f"Error: subjects parameter is incorrect \n subjects is {subjects} but it must be 'volunteer' or 'OR'")
  101. # apply norm
  102. sdata.append(dd_sdb)
  103. times += list(dd['t'])
  104. labels += list(dd['l'])
  105. egq += list(dd_egq)
  107. ddf = pd.DataFrame(data = np.hstack([np.array([times]).T,
  108. np.hstack(sdata).T,
  109. np.array([egq]).T,
  110. np.array([labels]).T]),
  111. columns=np.hstack(['times', dd['f'], 'egq', 'l']))
  113. ddf['caseid'] = caseid
  114. ddf = ddf.astype({"egq":'int64',"caseid":'category'})
  115. # assemble the dict
  116. dataset_dict = {}
  117. dataset_dict['fs'] = dd['f']
  118. dataset_dict['caseid'] = np.array(caseid)
  119. dataset_dict['labels'] = np.array(labels)
  120. dataset_dict['times'] = np.array(times)
  121. dataset_dict['Sdb'] = np.hstack(sdata)
  122. dataset_dict['subjects'] = subjects
  123. dataset_dict['egq'] = np.array(egq).astype(bool)
  124. dataset_dict['df'] = ddf
  125. self.datasets[dsetname] = dataset_dict
  127. def featurize(self, t_df, fs, v_df=None, which_features=['Sdb','bands','PCA','LDA','CNN']):
  128. """Function for generating the features that a classifier actually uses.
  129. This is called when a classifier is fit and is applied to only the training data.
  131. returns training_features, pca, lda, validation_features
  132. TODO
  133. Parameters
  134. ----------
  135. which_features : list (optional)
  136. which features to featurize
  137. default is all so ['Sdb','bands','PCA','LDA','CNN']
  140. Returns
  141. -------
  142. [type]
  143. [description]
  144. """
  146. # filter out signal dropout, unknown labels ONLY FOR FITTING
  147. good_tdf = t_df[(t_df['l'].isin([0,1])) & (t_df['egq']==1)]
  149. # do PCA and LDA
  150. pca_sdb = decomposition.pca.PCA(n_components=3)
  151. pca_sdb.fit(good_tdf[np.array(fs, dtype='str')].values)
  152. pca_10 = decomposition.pca.PCA(n_components=10)
  153. pca_10.fit(good_tdf[np.array(fs, dtype='str')].values)
  154. lda_sdb = discriminant_analysis.LinearDiscriminantAnalysis(n_components=1)
  155. lda_sdb.fit(good_tdf[np.array(fs, dtype='str')].values, good_tdf['l'].values)
  156. if 'CNN' in which_features:
  157. pca_cnn = fit_pca_cnn(np.unique(good_tdf.caseid))
  158. else:
  159. pca_cnn = []
  161. # temp storage of pca_sdb and lda_sbd
  162. self.pca_lda = [pca_10, lda_sdb]
  164. # apply to the training data (regardless of label)
  165. train_feats = self._featurize_df(t_df,which_features,fs,pca_sdb,lda_sdb,pca_cnn)
  167. # check whether there is a validation df to featurize
  168. if isinstance(v_df, pd.DataFrame):
  169. val_feats = self._featurize_df(v_df,which_features,fs,pca_sdb,lda_sdb,pca_cnn)
  171. return train_feats, val_feats
  172. else:
  173. return train_feats
  175. def _featurize_df(self, df, which_features, fs, pca_sdb, lda_sdb, pca_cnn):
  176. """
  177. Helper function to featurize a data frame
  178. which_features is a list that may contain any of ['Sdb','bands','PCA','LDA','CNN']
  179. assumes features have already been fit.
  180. """
  181. feats_dict = {}
  182. if 'Sdb' in which_features:
  183. # Sdb df unchanged
  184. feats_dict['Sdb'] = df
  186. if 'bands' in which_features:
  187. # tb df from bands
  188. bands = pd.DataFrame(data = Sdb_to_bands(df[np.array(fs, dtype='str')].values, fs))
  189. dfbands = pd.concat([df['times'], bands, df['egq'], df['l'], df['caseid']], axis=1)
  190. feats_dict['bands'] = dfbands
  192. if 'PCA' in which_features:
  193. # same for pca
  194. PCA = pd.DataFrame(data = pca_sdb.transform(df[np.array(fs, dtype='str')].values))
  195. dfPCA = pd.concat([df['times'], PCA, df['egq'], df['l'], df['caseid']], axis=1)
  196. feats_dict['PCA'] = dfPCA
  198. if 'LDA' in which_features:
  199. # same for LDA
  200. LDA = pd.DataFrame(data = lda_sdb.transform(df[np.array(fs, dtype='str')].values))
  201. dfLDA = pd.concat([df['times'], LDA, df['egq'], df['l'], df['caseid']], axis=1)
  202. feats_dict['LDA'] = dfLDA
  204. if 'CNN' in which_features:
  205. # helper fcn for CNN
  206. dfCNN = apply_pca_cnn(pca_cnn,np.unique(df.caseid))
  207. feats_dict['CNN'] = dfCNN
  209. return feats_dict
  211. def add_model(self, ctype='lr', ftype='Sdb', ttype='standard',
  212. mname='lr_standard_Sdb'):
  213. """Adds a model comprised of: a featureset, a classifier, and a method for handling the timeseries.
  214. Classifier types: 'lr', 'svr'
  215. Timeseries types: 'standard', 'hmm2', 'hmmfree'
  216. Feature types: 'Sdb', 'bands', 'PCs', 'LD', 'CNN'
  218. Parameters
  219. ----------
  220. ctype : str, optional
  221. classifier type, by default 'lr'
  222. ftype : str, optional
  223. features for the model, by default 'Sdb'
  224. ttype : str, optional
  225. time series treatment, by default 'standard'
  226. mname : str, optional
  227. model name, by default 'lr_standard_Sdb'
  228. """
  229. # the ultimate classification approach
  230. if ctype == 'svr':
  231. model = svm.SVR(kernel='linear', cache_size=8000)
  232. elif ctype == 'lr':
  233. model = linear_model.LogisticRegression()
  235. # the method for treatment of the timeseries
  236. model.timeseries = ttype
  237. if ttype == 'standard':
  238. model.hmm = None
  239. elif ttype == 'hmm2':
  240. hmm = hmm_model.GaussianHMM(2, algorithm='viterbi', n_iter=10)
  241. model.hmm = hmm
  242. elif ttype == 'hmmfree':
  243. hmm = hmm_model.GaussianHMM(6, algorithm='viterbi', n_iter=10)
  244. model.hmm = hmm
  245. else:
  246. raise Exception("The given timeseries treatment does not fall into an approved type (standard, hmm2, hmmfree).")
  248. # the features fed into the timeseries treatment
  249. model.ftype = ftype
  251. # and now append it
  252. self.models[mname] = model
  254. def train_model(self, mname, training_dname, timer=False, which_features=['Sdb','bands','PCA','LDA','CNN']):
  255. """Fits model. The classifier knows which features to take and whether
  256. or not to HMM it.
  258. Parameters
  259. ----------
  260. mname : str
  261. model name
  262. training_dname : str
  263. dataset to fit on
  264. timer : bool, optional
  265. whether we want it trained, by default False
  266. """
  268. if timer is True:
  269. lap = laptimer()
  271. # collect the data
  272. model = self.models[mname]
  273. dataset = self.datasets[training_dname]
  274. ddf = dataset['df']
  275. fs = dataset['fs']
  277. # do the featurization step
  278. train_feats = self.featurize(ddf, fs, which_features=which_features)
  280. # get the specific feature dict
  281. train_feat = train_feats[model.ftype]
  283. # drop cols
  284. nonvalue_cols = ['times', 'egq', 'l', 'caseid']
  286. # perform the timeseries analysis by taking only eeg quality spots
  287. if model.timeseries == 'standard':
  288. # no treatment of the timeseries as a timeseries
  289. training_series = train_feat[train_feat['egq']==1].drop(nonvalue_cols, axis=1).values
  290. training_labels = train_feat[train_feat['egq']==1]['l']
  292. else:
  293. # get the training values from the HMM timeseries
  294. hmm = model.hmm
  295. train_lengths = _continuous_lengths(train_feat)
  296. hmm.fit(train_feat[train_feat['egq']==1].drop(nonvalue_cols, axis=1).values,
  297. train_lengths)
  298. # calculate posterior probabilities for each state in order to train logistic regression
  299. posteriors = hmm.score_samples(
  300. train_feat[train_feat['egq']==1].drop(nonvalue_cols, axis=1).values,
  301. train_lengths)[1]
  303. # ## calcualte AIC for model parameterized in this way
  304. # logprob = hmm.decode(train_feat, algorithm='viterbi')[0]
  305. # n_params = 2*hmm.n_components*hmm.n_features +(hmm.n_components)**2 -1
  306. # aic = 2*(n_params) - 2*logprob
  307. # hmm.aic = aic
  309. training_series = posteriors
  310. training_labels = train_feat[train_feat['egq']==1]['l']
  312. # perform training, then get val py
  313. model.fit(training_series, training_labels)
  314. model.isfit = True
  316. # used to featurize the validation data
  317. model.training_info = [ddf, fs]
  319. # give the time of the fitting
  320. if timer is True:
  321. print(f"Processing time: {np.round(lap(),3)}")
  323. def crossvalidation(self, mname, dname, n_held=1, timer=False,
  324. features=['Sdb', 'bands', 'PCA', 'LDA', 'CNN']):
  325. """Performs cross-validation noting that each time, a different n_held must be held out. We want to do every combination available.
  327. Parameters
  328. ----------
  329. mname : str
  330. model name
  331. dname : str
  332. dataset used for crossval
  333. n_held : int, optional
  334. number held out in crossval, by default 1
  335. timer : bool, optional
  336. whether to print the time taken to run, by default False
  337. features : list, optional
  338. which features should be generated, by default ['Sdb', 'bands', 'PCA', 'LDA', 'CNN']
  340. Returns
  341. -------
  342. list
  343. crossvalidation AUCs
  344. """
  346. # collect the data
  347. model = self.models[mname]
  348. dataset = self.datasets[dname]
  349. ddf = dataset['df']
  350. fs = dataset['fs']
  352. # get combinations for training / val splits
  353. unique_trials = np.unique(dataset['caseid'])
  354. combs = list(combinations(unique_trials, len(unique_trials)-n_held))
  356. # what we collect from each crossval iteration
  357. model_performance = [] # AUC over single left out case
  358. for fold in combs:
  359. # split for featurization
  360. train_df = ddf[ddf['caseid'].isin(fold)].reset_index(drop=True)
  361. val_df = ddf[~ddf['caseid'].isin(fold)].reset_index(drop=True)
  363. # do the featurization step
  364. train_feats, val_feats = self.featurize(train_df, fs, v_df=val_df,
  365. which_features=features)
  367. # get the specific feature dict
  368. train_feat = train_feats[model.ftype]
  369. val_feat = val_feats[model.ftype]
  371. # drop cols
  372. nonvalue_cols = ['times', 'egq', 'l', 'caseid']
  374. # perform the timeseries analysis by taking only eeg quality spots
  375. if model.timeseries == 'standard':
  376. # no treatment of the timeseries as a timeseries
  377. training_series = train_feat[train_feat['egq']==1].drop(nonvalue_cols, axis=1).values
  378. training_labels = train_feat[train_feat['egq']==1]['l']
  379. validation_series = val_feat[val_feat['egq']==1].drop(nonvalue_cols, axis=1).values
  380. validation_labels = val_feat[val_feat['egq']==1]['l']
  382. else:
  383. # get the training values from the HMM timeseries
  384. hmm = model.hmm
  385. train_lengths = _continuous_lengths(train_feat)
  386. hmm.fit(train_feat[train_feat['egq']==1].drop(nonvalue_cols, axis=1).values,
  387. train_lengths)
  388. # calculate posterior probabilities for each state in order to train logistic regression
  389. posteriors = hmm.score_samples(
  390. train_feat[train_feat['egq']==1].drop(nonvalue_cols, axis=1).values,
  391. train_lengths)[1]
  393. # ## calcualte AIC for model parameterized in this way
  394. # logprob = hmm.decode(train_feat, algorithm='viterbi')[0]
  395. # n_params = 2*hmm.n_components*hmm.n_features +(hmm.n_components)**2 -1
  396. # aic = 2*(n_params) - 2*logprob
  397. # hmm.aic = aic
  399. training_series = posteriors
  400. training_labels = train_feat[train_feat['egq']==1]['l']
  402. val_lengths = _continuous_lengths(val_feat)
  403. try:
  404. val_posteriors = hmm.score_samples_fwd(
  405. val_feat[val_feat['egq']==1].drop(nonvalue_cols, axis=1).values,
  406. val_lengths)[1]
  408. except:
  409. print('WARNING: You are not using the modified version of HMM learn')
  410. print('Your classifier may be using the backward algorithm to predict consciousness')
  411. print('This does not affect the performance of the model. It only means this classifier could not be used in real time')
  412. print('For access to the forward-only hmmlearn see https://github.com/benyameister/hmmlearn/blob/master/README.rst')
  413. val_posteriors = hmm.score_samples(
  414. val_feat[val_feat['egq']==1].drop(nonvalue_cols, axis=1).values,
  415. val_lengths)[1]
  416. validation_series = val_posteriors
  417. validation_labels = val_feat[val_feat['egq']==1]['l']
  419. # perform training, then get val py
  420. model.fit(training_series, training_labels)
  421. model.isfit = True
  422. py = model.predict_proba(validation_series)[:,1]
  424. # save roc from each split
  425. fpr, tpr = roc_curve(validation_labels, py)[:2]
  426. auc_split = auc(fpr, tpr)
  428. model_performance.append(auc_split)
  429. return model_performance
  431. def save_model(self, path, mname):
  432. """
  433. Deprecated.
  434. """
  435. model = self.models[mname]
  436. pickle.dump(model, open(path, 'wb'))
  438. def load_model(self, path, mname):
  439. """Deprecated."""
  440. self.models[mname] = pickle.load(open(path, 'rb'))
  443. def validate_model(self, mname, validation_dname, timer=False):
  444. """applies a model to data within dsetname dataset.
  445. the resulting predictions are appended to the dataset with mname.
  446. ensures that norming is same.
  448. Note that this predicts *everything*, and nan labels must be stripped later (similarly hmm does predict nan labels but not egq=0 because data must be continuous).
  450. Also note that the features are generated using the PCA or LDA fit from the training data, which is now part of the model.
  451. """
  453. if timer is True:
  454. lap = laptimer()
  456. # collect the data
  457. model = self.models[mname]
  458. dataset = self.datasets[validation_dname]
  459. ddf = dataset['df']
  460. fs = dataset['fs']
  461. ftype=model.ftype
  463. # ensure model has been fit
  464. assert model.isfit, "Model is not yet fit!"
  466. # do the featurization step
  467. tddf, tfs = model.training_info
  468. assert all(tfs==fs), "Training/validation frequencies of MTSGs must be same!"
  470. val_feats = self.featurize(tddf, tfs, v_df=ddf, which_features=[ftype])[1]
  472. # get the specific feature dict
  473. val_feat = val_feats[model.ftype]
  475. # drop cols
  476. nonvalue_cols = ['times', 'egq', 'l', 'caseid']
  478. # perform the timeseries analysis by taking only eeg quality spots
  479. if model.timeseries == 'standard':
  480. # no treatment of the timeseries as a timeseries
  481. validation_series = val_feat[val_feat['egq']==1].drop(nonvalue_cols, axis=1).values
  483. else:
  484. # get the training values from the HMM timeseries
  485. hmm = model.hmm
  486. val_lengths = _continuous_lengths(val_feat[val_feat['egq']==1])
  487. try:
  488. val_posteriors = hmm.score_samples_fwd(
  489. val_feat[val_feat['egq']==1].drop(nonvalue_cols, axis=1).values,
  490. val_lengths)[1]
  492. except:
  493. print('WARNING: You are not using the modified version of HMM learn')
  494. print('Your classifier may be using the backward algorithm to predict consciousness')
  495. print('This does not affect the performance of the model. It only means this classifier could not be used in real time')
  496. print('For access to the forward-only hmmlearn see https://github.com/benyameister/hmmlearn/blob/master/README.rst')
  497. val_posteriors = hmm.score_samples(
  498. val_feat[val_feat['egq']==1].drop(nonvalue_cols, axis=1).values,
  499. val_lengths)[1]
  501. validation_series = val_posteriors
  503. py = model.predict_proba(validation_series)[:,1]
  505. # only return wehre we made predictions
  506. ddf_val = ddf[ddf['egq']==1].copy()
  507. ddf_val['py'] = py
  509. # append the validation result to the model itself
  510. if not hasattr(model, "val_result"):
  511. model.val_result = {}
  512. model.val_result[validation_dname] = ddf_val
  513. model.isval = True
  514. self.models[mname] = model
  516. # give the time of the fitting
  517. if timer is True:
  518. print(f"Processing time: {np.round(lap(),3)}")
  520. def roc_auc(self, mname, dname, cases='all'):
  521. """returns fpr, tpr, auc
  522. strips regions with bad eeg quality or NANs
  523. """
  525. ddf = self.models[mname].val_result[dname]
  526. if cases=='all':
  527. validation_l = ddf[ddf['egq']==1]['l'].values
  528. validation_py = ddf[ddf['egq']==1]['py'].values
  529. vl = validation_l[~np.isnan(validation_l)]
  530. py = validation_py[~np.isnan(validation_l)]
  531. elif type(cases) is list:
  532. validation_l = ddf[(ddf['egq']==1) & (ddf['caseid'].isin(cases))]['l'].values
  533. validation_py = ddf[(ddf['egq']==1) & (ddf['caseid'].isin(cases))]['py'].values
  534. vl = validation_l[~np.isnan(validation_l)]
  535. py = validation_py[~np.isnan(validation_l)]
  537. fpr, tpr, thr = roc_curve(vl, py)[:3]
  538. auct = auc(fpr, tpr)
  539. sens_plus_spec = (1-fpr)+tpr
  540. thr_opt = thr[np.argmax(sens_plus_spec)]
  541. return fpr, tpr, thr, auct, thr_opt
  543. def acc(self, mname, dname, thr=0.5, cases='all'):
  544. """returns fpr, tpr, auc
  545. strips regions with bad eeg quality or NANs
  546. """
  548. ddf = self.models[mname].val_result[dname]
  549. if cases=='all':
  550. validation_l = ddf[ddf['egq']==1]['l'].values
  551. validation_py = ddf[ddf['egq']==1]['py'].values
  552. vl = validation_l[~np.isnan(validation_l)]
  553. py = validation_py[~np.isnan(validation_l)]
  554. elif type(cases) is list:
  555. validation_l = ddf[(ddf['egq']==1) & (ddf['caseid'].isin(cases))]['l'].values
  556. validation_py = ddf[(ddf['egq']==1) & (ddf['caseid'].isin(cases))]['py'].values
  557. vl = validation_l[~np.isnan(validation_l)]
  558. py = validation_py[~np.isnan(validation_l)]
  560. ys = py>=thr
  561. acc = accuracy_score(ys, vl)
  562. return acc
  565. def plot_roc(self, mname, dsetname, ax=None, label='', c='b', ls='-'):
  566. """ plots the roc """
  567. if ax is None:
  568. ax = plt.subplot()
  570. fpr, tpr, thr, auct, _ = self.roc_auc(mname, dsetname)
  571. ax.plot(fpr, tpr, label=f'{label} AUC = {np.round(auct, 3)}', c=c,
  572. ls=ls)
  573. ax.plot([0, 1], [0, 1], 'k:')
  574. ax.set_xlabel('FPR')
  575. ax.set_ylabel('TPR')
  576. ax.legend()
  578. def save_inference_table(self, path, mname, dname):
  579. """ Deprecated. """
  580. inf_tab = self.models[mname].val_result[dname]
  581. inf_tab.to_csv(path+mname+'.csv')
  585. # utility functions
  587. class laptimer:
  588. """
  589. Whenever you call it, it times laps.
  590. """
  592. def __init__(self):
  593. self.time = time()
  595. def __call__(self):
  596. ret = time() - self.time
  597. self.time = time()
  598. return ret
  600. def __str__(self):
  601. return "%.3E" % self()
  603. def __repr__(self):
  604. return "%.3E" % self()
  607. def save_nparray(filename, nparray, colnames=None):
  608. """
  609. Uses pandas to save a numpy array with column headers.
  610. """
  611. assert(len(colnames)) == nparray.shape[1], "columns do not match table"
  612. output_df = pd.DataFrame(data=nparray, columns=colnames)
  613. output_df.to_csv(filename, index=False)
  616. def save_inference_table(filename, table):
  617. """helper function for saving inference tables"""
  618. save_nparray(filename, table, colnames=['case_id', 't', 'p_y', 'y'])
  621. def load_volunteer_todict(path, name):
  622. """
  623. Loads .csv files into a dict. I find this a useful utility.
  624. """
  625. f = np.genfromtxt(path + '/' + name + '_f.csv', delimiter=',')
  626. t = np.genfromtxt(path + '/' + name + '_t.csv', delimiter=',')
  627. l = np.genfromtxt(path + '/' + name + '_l.csv', delimiter=',')
  628. Sdb = np.genfromtxt(path + '/' + name + '_Sdb.csv', delimiter=',')
  629. egq = np.ones(len(t))
  630. return {'f': f, 't': t, 'Sdb': Sdb, 'l': l, 'name': name,
  631. 'egq': egq}
  634. def load_OR_todict(path, name):
  635. """
  636. Loads .csv files into a dict including events information from CL files
  637. """
  638. try:
  639. f = np.genfromtxt(path + '/' + name + '_f.csv', delimiter=',')
  640. t = np.genfromtxt(path + '/' + name + '_t.csv', delimiter=',')
  641. l = np.genfromtxt(path + '/' + name + '_l.csv', delimiter=',')
  642. Sdb = np.genfromtxt(path + '/' + name + '_Sdb.csv', delimiter=',')
  643. # events = np.genfromtxt(path +'/' + name + '_events.csv',
  644. # delimiter=',', usecols=[0, 1], dtype="str")
  645. egq = np.genfromtxt(path + '/' + name + '_EEGquality.csv',
  646. delimiter=',')
  647. return_dict = {'f': f, 't': t, 'Sdb': Sdb, 'l':l, #'events': events,
  648. 'name': name, 'egq':egq}
  650. # get bolus if it's real
  651. if os.path.exists(path + '/' + name + '_bolus.csv'):
  652. bolus = np.genfromtxt(path + '/' + name + '_bolus.csv',
  653. delimiter=',', dtype='str')
  654. return_dict['bolus'] = bolus
  656. # and infusion if it's real
  657. if os.path.exists(path + '/' + name + '_infusion.csv'):
  658. infusion = np.genfromtxt(path + '/' + name + '_infusion.csv',
  659. delimiter=',', dtype='str')
  660. return_dict['infusion'] = infusion
  662. # and gas if it's real
  663. if os.path.exists(path + '/' + name + '_gas.csv'):
  664. gas = np.genfromtxt(path + '/' + name + '_gas.csv',
  665. delimiter=',', dtype='str')
  666. return_dict['gas'] = gas
  667. return return_dict
  669. except OSError:
  670. print(f"Incomplete data for {name}")
  671. except ValueError:
  672. print(f"Incorrect CSV format for {name}")
  675. def _continuous_lengths(data_df, dt=2):
  676. """
  677. calcualtes lengths of continuous observations of eeg data for HMM processing
  678. This method is written for external validation data which has spontaneous
  679. drop out of signal. That is why a 'egq' array is expected from the data _dict
  680. This method still works for volunteer data, though it is analogous to stacking
  681. the lengths of each case
  683. Parameters
  684. ----------
  685. data_df : pandas dataframe
  686. should have entries 'egq' and 'caseid'
  687. dt: float
  688. spacing in t
  690. Returns
  691. -------
  692. lengths : numpy array
  693. array of lengths. Each length is an integer number that is the number of continuous samples.
  694. Easy check that it is being computed correctly is sum(lengths) = len(goodeeg_times)
  695. """
  697. lengths = []
  698. goodeeg_times = data_df[data_df['egq']==1]['times'].values
  699. length = 1
  700. for ti, t in enumerate(goodeeg_times[:-1]):
  701. if goodeeg_times[ti+1]-t==2:
  702. length+=1
  703. else:
  704. lengths.append(length)
  705. length=1
  706. lengths.append(length)
  708. return np.array(lengths)
  710. def hmm_aic(n_states_options,training_data,training_lengths,timer=True,plot=True):
  711. """
  712. calcualte AIC for a number of different states based on training data
  713. to see if there is an ideal number of states for HMMfree
  714. TODO - params and results
  715. Parameters
  716. ----------
  718. Results
  719. -------
  720. """
  721. AIC = []
  722. for k in n_states_options:
  723. if timer is True:
  724. lap = laptimer()
  725. model = hmm_model.GaussianHMM(
  726. k,
  727. algorithm='viterbi',
  728. n_iter=10)
  729. model.fit(training_data.transpose(), training_lengths)
  732. logprob = model.decode(training_data.transpose(),algorithm='viterbi')[0]
  733. n_params = 2*model.n_components*model.n_features +(model.n_components)**2 -1
  734. aic = aic = 2*(n_params) - 2*logprob
  735. AIC.append(aic)
  736. if timer is True and k>15:
  737. print(f"finished generating model with {k} states")
  738. print(f"Processing time: {np.round(lap(), 1)} seconds")
  739. plt.figure(figsize=(4,1.5))
  740. plt.plot(n_states_options,AIC)
  741. plt.title("AIC")
  742. return AIC
  745. def Sdb_to_bands(Sdb, fi):
  746. """Gets total power in slow-delta, theta, alpha, beta, gamma(?) range.
  747. """
  748. S = 10**(Sdb.T/10)
  750. sr = np.logical_and(fi>=0, fi<=1.)
  751. dr = np.logical_and(fi>1, fi<4)
  752. tr = np.logical_and(fi>=4, fi<8)
  753. ar = np.logical_and(fi>=8, fi<13)
  754. br = np.logical_and(fi>=13, fi<25)
  755. gr = np.logical_and(fi>=25, fi<50)
  757. slow_db = 10*np.log10(S[sr,:].sum(0)+1E-10)
  758. delta_db = 10*np.log10(S[dr,:].sum(0)+1E-10)
  759. theta_db = 10*np.log10(S[tr,:].sum(0)+1E-10)
  760. alpha_db = 10*np.log10(S[ar,:].sum(0)+1E-10)
  761. beta_db = 10*np.log10(S[br,:].sum(0)+1E-10)
  762. gamma_db = 10*np.log10(S[gr,:].sum(0)+1E-10)
  763. return np.vstack([slow_db, delta_db, theta_db, alpha_db, beta_db, gamma_db]).T
  765. def fit_pca_cnn(case_ids_train):
  766. """transform CNN botlenecks into its PC's
  767. NOTE: only uses volunteer bottlenecks right now because OR ones havent been generated
  769. Parameters
  770. ----------
  771. case_id_train : np.ndarray
  772. list of integers corresponding to unique case ids
  775. Returns
  776. -------
  777. cnn_pcs : pd.DataFrame
  778. first 3 principal component of CNN bottleneck values for each window
  779. also contains times and labels because they don't line up with other features
  780. pca_cnn : decomposition.pca.PCA
  781. fitted PCA to be used for validation data
  782. """
  783. btlncks = pd.read_csv(FP_BTLNCKS_VOLUNTEER)
  785. # filter btlncks to only use cases in training set
  786. train_case_inds = [case_id in case_ids_train for case_id in btlncks.case_id]
  787. btlncks = btlncks.loc[train_case_inds,:]
  789. pca_cnn = decomposition.pca.PCA(n_components=n_pcs_cnn)
  790. pca_cnn.fit(btlncks.loc[:,'btlnck_0':])
  791. return pca_cnn
  793. def apply_pca_cnn(pca_cnn, case_ids):
  794. """calculate principal components for validation data using prefitted pca
  795. NOTE: only uses volunteer bottlenecks right now because OR ones havent been generated
  797. Parameters
  798. ----------
  799. pca_cnn : decomposition.pca.PCA
  800. fitted PCA to be used for validation data
  801. case_id_train : np.ndarray
  802. list of strings corresponding to unique case ids
  805. Returns
  806. -------
  807. cnn_pcs : pd.DataFrame
  808. first 3 principal component of CNN bottleneck values for each window
  809. also contains times and labels because they don't line up with other features
  810. """
  811. btlncks = pd.read_csv(FP_BTLNCKS_VOLUNTEER)
  813. # filter btlncks to only use cases in training set
  814. train_case_inds = [case_id in case_ids for case_id in btlncks.case_id]
  815. btlncks = btlncks.loc[train_case_inds,:]
  816. pcs = pca_cnn.transform(btlncks.loc[:,'btlnck_0':])
  817. # output as dataframe so it keeps associated t and labels because has different alignment than other features
  818. column_names = [f'PC{n+1}' for n in range(n_pcs_cnn)]
  819. cnn_pcs = pd.DataFrame(pcs, columns=column_names)
  820. cnn_pcs['times'] = btlncks.t.values
  821. cnn_pcs['l'] = btlncks.is_conscious.values
  822. cnn_pcs['egq'] = np.array([1]*len(cnn_pcs['l']))
  823. cnn_pcs['caseid'] = btlncks.case_id
  824. return cnn_pcs

ConsciousnessClassifier.py at commit fe3bbdc, no license · at the source

Overview

Authors: Ping Shen1, Miao Da1, Zhongxia Shen1
  1. Huzhou Third Municipal Hospital, The Affiliated Hospital of Wenzhou Medical University, Huzhou, China
Institutions: Wenzhou Medical University (China)
Journal: Frontiers in medicine, volume 13, article 1794626
Dates: received 23 January 2026; accepted 23 March 2026; published online 28 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3389/fmed.2026.1794626 · PMID 42131603 · PMCID PMC13161752 · OpenAlex W7156607278
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: EEG (modality), clinical / translational (subfield)
Methods: Spectral & time-frequency, Statistics, Machine learning, Preprocessing, Complexity
Keywords: anesthesia, aperiodic, components, electroencephalography, general, Propofol, sevoflurane
Topic: Anesthesia and Sedative Agents (Anesthesiology and Pain Medicine, Medicine), according to OpenAlex
Citations: not cited yet (Europe PMC); 48 references in the paper

Abstract

Background: Current clinical anesthesia monitors often utilize drug-invariant indices that simplify cortical dynamics, potentially overlooking pharmacological nuances. While Propofol and Sevoflurane are both GABAergic, they may induce distinct neural states. This study aimed to identify the divergent periodic and aperiodic EEG signatures that distinguish these two regimens during the steady-state maintenance phase.

Methods: A retrospective analysis was conducted using data from an open-access clinical database comprising 44 surgical patients (Propofol group, n = 27; Sevoflurane group, n = 17). EEG data were extracted during the pharmacological steady-state (20 min post-loss of consciousness to 10 min pre-end of surgery). Seventeen features, including relative band power, alpha peak frequency, and aperiodic components, were then derived. A multivariate statistical framework utilizing subject-independent cross-validation and SHapley Additive exPlanations (SHAP) analysis was implemented to identify and rank the most discriminatory biological markers.

Results: The multivariate model achieved high discriminatory performance with a rigorous subject-level accuracy of 91.43%. Relative theta power, theta-to-alpha ratio, and alpha peak frequency were identified as the primary differentiators, occupying the top tiers of the SHAP importance ranking. Specifically, the Sevoflurane group exhibited a distinct elevation in theta-band prominence and a significant downward shift in alpha peak frequency (8.78 Hz vs. 10.88 Hz for Propofol). Furthermore, the aperiodic exponent emerged as a critical discriminatory feature, demonstrating a significantly steeper background spectral slope under Sevoflurane (2.37 vs. 2.07 for Propofol, p = 0.039). Conversely, alpha bandwidth (p = 0.263) and signal complexity measures (e.g., spectral entropy, p = 0.721) provided negligible discriminatory value.

Conclusion: Propofol and Sevoflurane maintain unconsciousness via distinct neurophysiological regimes. The differentiation between these two agents is primarily driven by structural oscillatory shifts, specifically theta-band prominence and alpha peak deceleration, along with steepened aperiodic background dynamics, rather than periodic bandwidth or overall signal complexity. These findings underscore the distinct cortical modulation patterns of different GABAergic anesthetics and support the development of agent-specific, multidimensional monitoring protocols to enhance precision in individualized brain state assessment.

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

Repository

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

johnabel/GABAergic_unconsciousness

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: fe3bbdc5dc63c2aad076e7a1bc099c8b6baa866b, 20 February 2021
Languages: Python (12)
Size: 22 files, 12 scripts
Software Heritage: not archived
Found in: “Data availability statement”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (12 files), NumPy (12 files), pandas (8 files), SciPy (8 files), scikit-learn (7 files), seaborn (6 files)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
13 files
At the source: github.com/johnabel/GABAergic_unconsciousness

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:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 12 scripts, each with its path and the digest of its content;
  • 1 match 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 statement

Publicly available datasets were analyzed in this study. This data can be found at: “GABAergic Anesthetic-Induced Unconsciousness and Recovery” https://github.com/johnabel/GABAergic_unconsciousness.

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, 30 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 3 authors, 7 keywords, 48 references.

Cite

This paper

Shen, P., Da, M., & Shen, Z. (2026). Divergent periodic and aperiodic EEG signatures of Propofol versus sevoflurane anesthesia: a comparative neurophysiological study. Frontiers in medicine, 13, 1794626. https://doi.org/10.3389/fmed.2026.1794626

BibTeX

@article{shen2026divergent,
author = {Shen, Ping and Da, Miao and Shen, Zhongxia},
title = {{Divergent periodic and aperiodic EEG signatures of Propofol versus sevoflurane anesthesia: a comparative neurophysiological study}},
journal = {Frontiers in medicine},
year = {2026},
month = apr,
volume = {13},
pages = {1794626},
publisher = {Frontiers Media SA},
issn = {2296-858X},
doi = {10.3389/fmed.2026.1794626},
url = {https://doi.org/10.3389/fmed.2026.1794626},
pmid = {42131603},
pmcid = {PMC13161752}
}

RIS

TY - JOUR
AU - Shen, Ping
AU - Da, Miao
AU - Shen, Zhongxia
TI - Divergent periodic and aperiodic EEG signatures of Propofol versus sevoflurane anesthesia: a comparative neurophysiological study
T2 - Frontiers in medicine
J2 - Front Med (Lausanne)
PY - 2026
DA - 2026/04/28
VL - 13
SP - 1794626
SN - 2296-858X
PB - Frontiers Media SA
DO - 10.3389/fmed.2026.1794626
UR - https://doi.org/10.3389/fmed.2026.1794626
LA - en
ER -

CSL-JSON

{
"id": "10.3389/fmed.2026.1794626",
"type": "article-journal",
"title": "Divergent periodic and aperiodic EEG signatures of Propofol versus sevoflurane anesthesia: a comparative neurophysiological study",
"container-title": "Frontiers in medicine",
"author": [
{
"family": "Shen",
"given": "Ping"
},
{
"family": "Da",
"given": "Miao"
},
{
"family": "Shen",
"given": "Zhongxia"
}
],
"container-title-short": "Front Med (Lausanne)",
"volume": "13",
"page": "1794626",
"DOI": "10.3389/fmed.2026.1794626",
"PMID": "42131603",
"PMCID": "PMC13161752",
"ISSN": "2296-858X",
"publisher": "Frontiers Media SA",
"URL": "https://doi.org/10.3389/fmed.2026.1794626",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
28
]
]
}
}

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.1002/mco2.70980 [code]
An Intraoperative EEG Biomarker for Postoperative Delirium Predicting Based on Interpretable Deep Learning Framework.
Journal: MedComm
In common: seaborn, scikit-learn, pandas, 3 other tools, EEG, clinical / translational, 2 references
[2] doi:10.1093/cercor/bhag113 [code]
Long-term reliability and stability of parameterized resting state EEG: evidence from a five-year follow-up.
Journal: Cerebral cortex (New York, N.Y. : 1991)
In common: seaborn, scikit-learn, pandas, 3 other tools, EEG, 2 references
[3] doi:10.1038/s42003-026-09769-7
Distinct origins of human low and high alpha rhythms revealed by simultaneous EEG-SEEG.
Journal: Communications biology
In common: EEG, 4 references
[4] doi:10.7554/elife.100605 [code]
Age-related changes in ‘cortical’ 1/f dynamics are linked to cardiac activity
Journal: n/a
In common: seaborn, scikit-learn, pandas, 3 other tools, 2 references
[5] doi:10.1111/ejn.70255 [code]
A Systematic Review of Aperiodic Neural Activity in Clinical Investigations
Journal: n/a
In common: seaborn, pandas, Matplotlib, 1 other tool, EEG, clinical / translational, 2 references
[6] doi:10.1038/s42003-026-10205-z [code]
Source-space EEG alpha activity reveals brain age gaps due to neurodegeneration and disparity.
Journal: Communications biology
In common: seaborn, scikit-learn, pandas, 3 other tools, EEG, 1 reference
[7] doi:10.1126/sciadv.adz6517 [code]
Corticosterone-linked microglial activity underpins sexually dimorphic neuroplasticity after ketamine anesthesia.
Journal: Science advances
In common: seaborn, scikit-learn, pandas, 3 other tools, 1 reference
[8] doi:10.1038/s41467-026-72454-0 [code]
Dynamic neuronal ensembles encode burst-suppression revealed by cortex-wide optical-electrical interfaces.
Journal: Nature communications
In common: seaborn, scikit-learn, pandas, 3 other tools, 1 reference
[9] doi:10.1038/s41593-026-02205-3 [code]
Competitive interactions shape mammalian brain network dynamics and computation.
Journal: Nature neuroscience
In common: seaborn, scikit-learn, pandas, 3 other tools, 1 reference
[10] doi:10.1038/s41467-026-74227-1 [code]
Age-related changes in behavioural and neural variability in a decision-making task.
Journal: Nature communications
In common: seaborn, scikit-learn, pandas, 3 other tools, 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.