OSCR

Multimodal Image Guidance in Subthalamic Deep Brain Stimulation for Parkinson's Disease.

Code ↔ Paper

7 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 7 matches · 2 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Materials and Methods › DBS Electrode Localization and Electric Field Calculation ↔ ea_normalize_ants.m, the whole file · a weak match · score 0.75 · ANTs SyN, Advanced Normalization Tools, template space, subcortical, MNI, reconstructed
  2. [2] § Materials and Methods › DBS Electrode Localization and Electric Field Calculation ↔ ea_normalize_schoenecker.m, the whole file · a weak match · score 0.72 · ANTs SyN, co registered, template space, CT, subcortical, post
  3. [3] § Materials and Methods › DBS Electrode Localization and Electric Field Calculation ↔ explorers/unifiedmapping_explorer/ea_unifiedmapping.m, lines 125–236 · score 0.64 · electric field, native space, OSS DBS, voxels, Stimulation, modeling
  4. [4] § Results › Optimal DBS Target Sites Across Multiple Modalities ↔ explorers/unifiedmapping_explorer/ea_unifiedmapping.m, lines 1–116 · score 0.63 · fold CV, structural connectivity, CVs, mirroring, peaked, thresholding
  5. [5] § Results › Optimal DBS Target Sites Across Multiple Modalities ↔ explorers/fiberfiltering_explorer/ea_disctract.m, lines 1–113 · score 0.62 · fold CV, structural connectivity, CVs, mirroring, overlapped, peaked
  6. [6] § Results › Optimal DBS Target Sites Across Multiple Modalities ↔ explorers/networkmapping_explorer/ea_networkmapping.m, lines 1–74 · score 0.59 · cross validation, DBS Network, network maps, fold, peaked, CV
  7. [7] § Materials and Methods › Statistical Analysis › Validation of the Models and Calculation of Surrogate Values ↔ explorers/networkmapping_explorer/ea_networkmapping.m, lines 1–74 · score 0.55 · fold CV, cross validation, networks, correlated, variables, patient

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

MATLAB · 1,712 lines · 85 KB · GPL-3.0 · 2 matches

  1. classdef ea_unifiedmapping < handle
  2. % Discriminative fiber class to handle visualizations of discriminative fibers in lead dbs resultfig / 3D Matlab figures
  3. % A. Horn
  4. properties (SetObservable)
  5. tool = 1; %option of 1 = sweetspotmapping, 2 = fiberfiltering, 3 = networkmapping
  6. fileformatversion; % 1.2 is current format
  7. M % content of lead group project
  8. resultfig % figure handle to plot results
  9. ID % name / ID of unified analysis object
  10. posvisible = 1 % pos tract/vox visible
  11. negvisible = 0 % neg tract/vox visible
  12. roivisible = 0 % show ROI (usually VTAs)
  13. connfibvisible = 0 % show all connected tracts in white
  14. showposamount = [25 25] % two entries for right and left
  15. shownegamount = [25 25] % two entries for right and left
  16. statmetric = nan;
  17. statsettings = struct;
  18. calcsettings = struct; %calcthreshold, connectivity_type
  19. threshstrategy = 'Percentage Relative to Peak'; % can be 'Relative to Amount' or 'Fixed Amount'
  20. multi_pathways = 0 % if structural connectome is devided into pathways (multiple .mat in dMRI_MultiTract)
  21. map_list % list that contains global indices of the first fibers in each pathway (relevant when multi_pathways = 1)
  22. pathway_list % list that contains names of pathways (relevant when multi_pathways = 1)
  23. connFiberInd % legacy
  24. switch_connectivity = 0 % flag if connectivity type was changed in the GUI
  25. nestedLOO = false % if true, will conducted LOO in the training set
  26. corrtype = 'Spearman' % correlation strategy in case of statmetric == Correlations / E-fields (Irmen 2020). In case of one-sample tests used for 'T-Tests' vs 'Wicoxon Signed Rank Tests'.
  27. SigmoidTransform = 0;
  28. multitractmode = 'Single Tract Analysis' % multi mode now coded by this value %should we use abreviation?
  29. numpcs = 4; % standard value of how many PCs to compute in case of PCA mode
  30. doactualprediction = 0; % set up nested CVs to carry out actual predictions of response variables
  31. predictionmodel = 'Linear'; % type of glm used to fit fiber values to actual scores
  32. showsignificantonly = 0
  33. alphalevel = 0.05
  34. multcompstrategy = 'FDR'; % could be 'Bonferroni'
  35. subscore
  36. explorerdrawn
  37. results = struct
  38. customRoi = struct % struct used only for pseudoM case (customRoi.isbinary and customRoi.minmax)
  39. activateby={}; % entry to use to show fiber activations
  40. cvlivevisualize = 0; % if set to 1 shows crossvalidation results during processing.
  41. basepredictionon = 'Mean of Scores';
  42. fiberdrawn % struct contains fibercell and vals drawn in the resultfig
  43. spotdrawn % struct contains nifti of the sweetspot
  44. networkdrawn %struct contains nifti of the networkmap
  45. vizmode
  46. smooth_fp = 0; %for networkmapping, smooth fingerprints
  47. normalize_fp = 0; %for networkmapping, normalize fingerprints
  48. cvmask = 'Gray Matter';
  49. model='Smoothed'; %for networkmapping
  50. drawobject = struct % actual streamtube handle
  51. drawvals % weights of the fibers drawn
  52. connfiberdrawn % struct contains white connected fibers
  53. conndrawobject % actial streamtube handle for the latter
  54. roidrawobject % actual patch handle for ROI/VTAs
  55. roidata % data used for ROIs (nifti file cell usually w/ Efields
  56. roiprotocol % protocol for drawing rois used to check if need to be redrawn.
  57. patientselection % selected patients to include. Note that connected fibers are always sampled from all (& mirrored) VTAs of the lead group file
  58. setlabels={};
  59. setselections={};
  60. customselection % selected patients in the custom test list
  61. allpatients % list of all patients (as from M.patient.list)
  62. mirrorsides = 0 % flag to mirror VTAs / Efields to contralateral sides using ea_flip_lr_nonlinear()
  63. responsevar % response variable
  64. responsevarlabel % label of response variable
  65. covars = {} % covariates
  66. covarlabels = {} % covariate labels
  67. analysispath % where to store results
  68. leadgroup % redundancy protocol only, path to original lead group project
  69. useExternalModel = false
  70. ExternalModelFile = 'None'
  71. % NM visualization
  72. NMviz=struct;
  73. % NMviz.vizmode='Regions'; % way to plot results
  74. % NMviz.model='Smoothed'; % in case of surface above, on which surface to plot.
  75. % NMviz.modelLH=1; % show left hemisphere
  76. % NMviz.modelRH=1; % show right hemisphere
  77. % stats: (how many fibers available and shown etc for GUI)
  78. modelNormalization = 'None';
  79. numBins=15;
  80. stats
  81. % additional settings:
  82. rngseed = 'default';
  83. Nperm = 1000 % how many permutations in leave-nothing-out permtest strategy
  84. kfold = 5 % divide into k sets when doing k-fold CV
  85. Nsets = 5 % divide into N sets when doing Custom (random) set test
  86. adjustforgroups = 1 % adjust correlations for group effects
  87. kIter = 1;
  88. roiintersectdata = {}; %roi, usually efield with which you can calculate fiber intersection
  89. roithresh = 200; % threshold above which efield metrics are considered
  90. drawTool = 'sweetspotmapping'; % active draw tool: sweetspotmapping, fiberfiltering, networkmapping
  91. activated = struct; % activated.fiberfiltering, activated.sweetspotmapping, activated.networkmapping % for activation status mapping.
  92. % misc
  93. runwhite = 0; % flag to calculate connected tracts instead of stat tracts
  94. e_field_metric = 'Magnitude'; % 'Magnitude' or 'Projection'
  95. %color
  96. colorbar % colorbar information
  97. posBaseColor = [1,1,1] % positive main color
  98. poscolor = [0.9176,0.2000,0.1373] % positive peak color
  99. negBaseColor = [1,1,1] % negative main color
  100. negcolor = [0.2824,0.6157,0.9725] % negative peak color
  101. hasResults = false; % results check
  102. end
  103. properties (Access = private)
  104. switchedFromSpace=3 % if switching space, this will protocol where from
  105. end
  106. methods
  107. function obj=ea_unifiedmapping(analysispath) % class constructor
  108. if exist('analysispath', 'var') && ~isempty(analysispath)
  109. obj.analysispath = analysispath;
  110. [~, ID] = fileparts(obj.analysispath);
  111. obj.ID = ID;
  112. end
  113. end
  114. function initialize(obj,datapath,resultfig)
  115. % statsettings
  116. % initial hard threshold to impose on (absolute) nifti files only when calculating the data
  117. obj.statsettings.doVoxels = 1;
  118. obj.statsettings.doFibers = 1;
  119. obj.statsettings.outcometype = 'gradual';
  120. obj.statsettings.stimulationmodel = 'Electric Field';
  121. obj.statsettings.efieldmetric = 'Peak'; % if statmetric == ;Correlations / E-fields (Irmen 2020)’, efieldmetric can calculate sum, mean or peak along tracts
  122. obj.statsettings.efieldthreshold = 200;
  123. obj.statsettings.nanthreshold = 0; % set values below this number to nan on the fly when calculating sweetspot statistics
  124. obj.statsettings.sweetspotresolution = 0.5; % resolution of sweetspot in mm
  125. obj.statsettings.connthreshold = 20;
  126. obj.statsettings.statfamily = 'Correlations'; % the
  127. obj.statsettings.stattest = 'Spearman';
  128. obj.statsettings.H0 = 'Average';
  129. obj.calcsettings.selectedTool = 1; %1 = sweetspotmapping, 2 = fiberfiltering, 3 = networkmapping;
  130. obj.calcsettings.calcthreshold = 100;
  131. obj.calcsettings.switch_connectivity = 1;
  132. obj.calcsettings.connectivity_type = 1; %1 = vta, 2 = PAM
  133. obj.calcsettings.functionalresolution = '2 mm';
  134. obj.calcsettings.structuralresolution = '2 mm';
  135. obj.calcsettings.calcmethod = 1; %1 = e-field based method, 2 = fiber based method
  136. obj.calcsettings.calcspace = 1; %0 = native space, 1 = MNI space
  137. obj.calcsettings.netmap_connectome = '';
  138. obj.calcsettings.fibfilt_connectome = '';
  139. if ~isfield(obj.NMviz,'modelRH')
  140. obj.NMviz.modelRH=1;
  141. end
  142. if ~isfield(obj.NMviz,'modelLH')
  143. obj.NMviz.modelLH=1;
  144. end
  145. obj.subscore.vars = {};
  146. obj.subscore.labels = {};
  147. obj.subscore.pcavars = {};
  148. obj.subscore.weights = [];
  149. obj.subscore.colors{1,1} = ea_color_wes('all');
  150. obj.subscore.colors{1,2} = flip(ea_color_wes('all'));
  151. obj.subscore.vis.showposamount = repmat([25,25],10,1); %total of 10 subscores - will delete when we know the total number of subscores
  152. obj.subscore.vis.shownegamount = repmat([25,25],10,1);
  153. obj.subscore.vis.pos_shown = repmat([25,25],10,1);
  154. obj.subscore.vis.neg_shown = repmat([25,25],10,1);
  155. obj.subscore.negvisible = zeros(10,1);
  156. obj.subscore.posvisible = ones(10,1);
  157. obj.subscore.splitbysubscore = 0;
  158. obj.subscore.special_case = 0;
  159. obj.covarlabels={};
  160. obj.fiberdrawn = struct();
  161. obj.spotdrawn = struct();
  162. obj.networkdrawn = struct();
  163. obj.vizmode = 'Regions';
  164. obj.model = 'Smoothed';
  165. obj.smooth_fp = 0;
  166. obj.normalize_fp = 0;
  167. obj.cvmask = 'Gray Matter';
  168. obj.fileformatversion=1.2; % new current version with settings to harmonize stats.
  169. datapath = GetFullPath(datapath);
  170. U = load(datapath, '-mat');
  171. if isfield(U, 'M') % Lead Group analysis path loaded
  172. obj.M = U.M;
  173. obj.leadgroup = datapath;
  174. testID = obj.M.guid;
  175. ea_mkdir([fileparts(obj.leadgroup),filesep,'UnifiedMappingExplorer',filesep]);
  176. id = 1;
  177. while exist([fileparts(obj.leadgroup),filesep,'UnifiedMappingExplorer',filesep,testID,'.explorer'],'file')
  178. testID = [obj.M.guid, '_', num2str(id)];
  179. id = id + 1;
  180. end
  181. obj.ID = testID;
  182. obj.resultfig = resultfig;
  183. if isfield(obj.M,'pseudoM')
  184. obj.allpatients = obj.M.ROI.list;
  185. obj.patientselection = 1:length(obj.M.ROI.list);
  186. obj.M.root = [fileparts(datapath),filesep];
  187. obj.M.patient.list = cell(size(obj.M.ROI.list,1), 1);
  188. for i = 1:size(obj.M.ROI.list,1)
  189. obj.M.patient.list{i,1} = obj.M.ROI.list{i,1};
  190. end
  191. obj.M.patient.group=obj.M.ROI.group; % copies
  192. else
  193. obj.allpatients = obj.M.patient.list;
  194. obj.patientselection = obj.M.ui.listselect;
  195. end
  196. obj.responsevar = obj.M.clinical.vars{1};
  197. obj.responsevarlabel = obj.M.clinical.labels{1};
  198. elseif isfield(U, explorer) % Saved explorer class loaded
  199. props = properties(U.explorer);
  200. for p = 1:length(props) %copy all public properties
  201. if ~(strcmp(props{p}, 'analysispath') && ~isempty(obj.analysispath) ...
  202. || strcmp(props{p}, 'ID') && ~isempty(obj.ID))
  203. obj.(props{p}) = U.explorer.(props{p});
  204. end
  205. end
  206. clear U
  207. else
  208. ea_error('You have opened a file of unknown type.')
  209. return
  210. end
  211. obj.compat_statmetric; % check and resolve for old statmetric code (which used to be integers)
  212. addlistener(obj,'activateby','PostSet',@activatebychange);
  213. % added a check here otherwise errors out for files w/o vatmodels
  214. if ~isfield(obj.M,'pseudoM')
  215. if ~isempty(obj.M.vatmodel) && contains(obj.M.vatmodel, 'OSS-DBS (Butenko 2020)')
  216. obj.statmetric = 3;
  217. end
  218. end
  219. end
  220. function compat_statmetric(obj)
  221. if ~ischar(obj.statmetric) % old language used:
  222. switch obj.statmetric % 3 was never used
  223. case 1
  224. obj.statmetric='Two-Sample T-Tests / VTAs (Baldermann 2019) / PAM (OSS-DBS)';
  225. case 2
  226. obj.statmetric='Correlations / E-fields (Irmen 2020)';
  227. case 4
  228. obj.statmetric='Proportion Test (Chi-Square) / VTAs (binary vars)';
  229. case 5
  230. obj.statmetric='Binomial Tests / VTAs (binary vars)';
  231. case 6
  232. obj.statmetric='Reverse T-Tests / E-Fields (binary vars)';
  233. case 7
  234. obj.statmetric='Plain Connections';
  235. case 8
  236. obj.statmetric='Odds Ratios / EF-Sigmoid (Jergas 2023)';
  237. case 9
  238. obj.statmetric='Weighted Linear Regression / EF-Sigmoid (Dembek 2023)';
  239. end
  240. end
  241. ea_unifiedmapping_compat_statmetrics2statsettings(obj);
  242. end
  243. function calculate(obj)
  244. % check that this has not been calculated before:
  245. %first store the rois for automatic calculations
  246. if ~isfield(obj.results,'roi')
  247. if isfield(obj.M,'pseudoM')
  248. vatlist = obj.M.ROI.list;
  249. else
  250. vatlist = ea_sweetspotmapping_getvats(obj);
  251. end
  252. for vat=1:length(vatlist)
  253. for side = 1:2
  254. vta_nii = ea_load_nii(vatlist{vat,side});
  255. obj.results.roi{vat,side} = vta_nii;
  256. end
  257. end
  258. end
  259. if obj.calcsettings.selectedTool == 1 % sweetspotmapping
  260. % in case of the sweetspot explorer, calculate rather means to
  261. % gather all E-Fields. To keep consistency of the logic with
  262. % discfiberexplorer and networkmappingexplorer, we will keep
  263. % the same name (calculate) for the function, nonetheless.
  264. % check that results aren't already there
  265. if isfield(obj.M,'pseudoM')
  266. vatlist = obj.M.ROI.list;
  267. else
  268. vatlist = ea_sweetspotmapping_getvats(obj);
  269. end
  270. [AllX,space] = ea_unifiedmapping_exportefieldmap(vatlist,obj);
  271. % Apply threshold: set all values below nanthreshold to
  272. % NaN, like in fiberfiltering
  273. for i = 1:numel(AllX)
  274. if ~isempty(AllX{i})
  275. AllX{i}(AllX{i} < obj.calcsettings.calcthreshold) = nan;
  276. end
  277. end
  278. obj.results.sweetspotmapping.efield = AllX;
  279. obj.results.sweetspotmapping.space = space;
  280. if ~isfield(obj.M,'pseudoM')
  281. try
  282. % get active coordinates, as well
  283. for pt=1:length(obj.M.patient.list)
  284. for side=1:2
  285. obj.results.sweetspotmapping.activecnt{side}(pt,:)=...
  286. mean(obj.M.elstruct(pt).coords_mm{side}(find(obj.M.S(pt).activecontacts{side}),:),1); %#ok<FNDSB> % find is necessary here
  287. end
  288. end
  289. for side=1:2
  290. obj.results.sweetspotmapping.activecnt{side}=[obj.results.sweetspotmapping.activecnt{side};ea_flip_lr_nonlinear(obj.results.sweetspotmapping.activecnt{side})];
  291. end
  292. end
  293. end
  294. elseif obj.calcsettings.selectedTool == 2 % Fiber Filtering
  295. % if multi_pathways = 1, assemble cfile from multiple
  296. % pathway.dat files in dMRI_MultiTract/Connectome_name/
  297. % stores the result in the LeadGroup folder
  298. % also merges fiberActivation_.._.mat and stores them in
  299. % stimulation folders
  300. if ~isempty(obj.results) % something has been calculated
  301. if isfield(obj.results,'fiberfiltering')
  302. fname = ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome);
  303. if isfield(obj.results.fiberfiltering, fname)
  304. connField = obj.results.fiberfiltering.(fname);
  305. % now check the specific subfields safely
  306. if (isfield(connField,'PAM_Ttest') && obj.calcsettings.connectivity_type==2) || ...
  307. (isfield(connField,'efield_mean') && obj.calcsettings.connectivity_type==1)
  308. answ = questdlg('This has already been calculated. Are you sure you want to re-calculate everything?', ...
  309. 'Recalculate Results','No','Yes','No');
  310. if ~strcmp(answ,'Yes')
  311. return
  312. end
  313. end
  314. end
  315. end
  316. end
  317. if obj.calcsettings.multi_pathways == 1
  318. [cfile, obj.map_list, obj.pathway_list] = ea_unifiedmapping_mergePathways(obj);
  319. else
  320. cfile = [ea_getconnectomebase('dMRI'), obj.calcsettings.fibfilt_connectome, filesep, 'data.mat'];
  321. end
  322. %check if files exist
  323. FilesExist = check_stimvols(obj);
  324. if isfield(obj.M,'pseudoM') % failsave - this should not be necessary but still making sure things are set correctly for the pseudoM case.
  325. obj.calcsettings.connectivity_type=1;
  326. obj.calcsettings.calcmethod=1;
  327. end
  328. switch obj.calcsettings.connectivity_type
  329. case 2 % if PAM, then just extracts activation states from fiberActivation.mat
  330. fprintf("Calculating using the PAM method. Using dMRI connectome: %s",obj.calcsettings.fibfilt_connectome);
  331. if all(FilesExist)
  332. calculate_on_pam(obj,cfile)
  333. end
  334. otherwise % check fiber recruitment via intersection with VTA
  335. if obj.calcsettings.calcmethod == 1 %'E-field/Voxel Based Method'
  336. fprintf("Calculating using the traditional E-field based method. Using dMRI connectome: %s",obj.calcsettings.fibfilt_connectome);
  337. if all(FilesExist)
  338. calculate_on_efield(obj,cfile)
  339. end
  340. elseif obj.calcsettings.calcmethod == 2 %'Fiber Based Method'
  341. % check whether to use new (calc_on_fibers) or old method:
  342. fprintf("Calculating using the Fiber based method. Using dMRI connectome: %s",obj.calcsettings.fibfilt_connectome);
  343. if all(FilesExist)%why can this not be psuedo M?
  344. calculate_on_fibers(obj,cfile)
  345. end
  346. end
  347. end
  348. elseif obj.calcsettings.selectedTool == 3 % network mapping
  349. if ~isempty(obj.results) % something has been calculated
  350. if isfield(obj.results,'networkmapping')
  351. if isfield(obj.results.networkmapping,ea_unifiedmapping_conn2connid(obj.calcsettings.netmap_connectome))
  352. answ=questdlg('This has already been calculated. Are you sure you want to re-calculate everything?','Recalculate Results','No','Yes','No');
  353. if ~strcmp(answ,'Yes')
  354. return
  355. end
  356. end
  357. end
  358. end
  359. if isfield(obj.M,'pseudoM')
  360. vatlist = obj.M.ROI.list;
  361. else
  362. %TODO:I have removed this from the networkmapping explorer folder and added it to the unified mapping explorer. Please adjust based on the future of the tool. I refrained from making a copy since the name of this script makes sense and would be redundant to change the name
  363. vatlist = ea_unified_nm_getvats(obj);
  364. end
  365. %TODO:I have removed this from the networkmapping explorer folder and added it to the unified mapping explorer. Please adjust based on the future of the tool. I refrained from making a copy since the name of this script makes sense and would be redundant to change the name
  366. [AllX] = ea_unified_nm_calcvals(vatlist, obj.calcsettings.netmap_connectome);
  367. obj.results.networkmapping.(ea_unifiedmapping_conn2connid(obj.calcsettings.netmap_connectome)).connval = AllX;
  368. % Functional connectome, add spacedef to results
  369. if contains(obj.calcsettings.netmap_connectome, ' > ')
  370. connid = (ea_unifiedmapping_conn2connid(obj.calcsettings.netmap_connectome));
  371. connName = regexprep(obj.calcsettings.netmap_connectome, ' > .*$', '');
  372. load([ea_getconnectomebase('fmri'), connName, filesep, 'dataset_volsurf.mat'], 'vol');
  373. obj.results.networkmapping.(connid).space = vol.space;
  374. if isfield(vol,'cifti')
  375. obj.results.networkmapping.(connid).space.cifti=vol.cifti; % add cifti space as well (hidden in nii space)
  376. obj.results.networkmapping.(connid).space.outidx=vol.outidx;
  377. obj.results.networkmapping.(connid).space.inidx=vol.inidx;
  378. end
  379. end
  380. end
  381. obj.hasResults = true;
  382. end
  383. function FilesExist = check_stimvols(obj)
  384. switch obj.calcsettings.connectivity_type
  385. case 2
  386. [~,FilesExist] = ea_unifiedmapping_getpams(obj);
  387. if ~all(FilesExist)
  388. answ=questdlg('It seems like PAM has not been (completely) run. We can initiate the process now, but this will take some time. Proceed?','PAM not run','yes','no','yes');
  389. switch answ
  390. case 'yes'
  391. options=ea_defaultoptions;
  392. options.prefs.machine.vatsettings.butenko_calcPAM=1;
  393. options.prefs.machine.vatsettings.butenko_calcVAT=0;
  394. options.prefs.machine.vatsettings.butenko_connectome=obj.calcsettings.fibfilt_connectome;
  395. options.groupdir=fileparts(obj.leadgroup);
  396. obj.M.vatmodel='OSS-DBS (Butenko 2020)';
  397. if isfield(obj.M.ui, 'stimSetMode') && obj.M.ui.stimSetMode
  398. options.stimSetMode = 1;
  399. else
  400. options.stimSetMode = 0;
  401. end
  402. filesToCalc = find(sum(FilesExist(1:length(obj.M.patient.list),:),2)<2)';
  403. calc_biophysical(obj,options,filesToCalc);
  404. case 'no'
  405. return
  406. end
  407. end
  408. %recheck
  409. [~,FilesExist] = ea_unifiedmapping_getpams(obj);
  410. otherwise
  411. if obj.calcsettings.calcmethod == 2 %Fiber based method
  412. [~,FilesExist] = ea_unifiedmapping_getlattice(obj);
  413. else
  414. if isfield(obj.M,'pseudoM')
  415. for entry=1:length(obj.M.ROI.list)
  416. FilesExist(entry)=exist(obj.M.ROI.list{entry},'file');
  417. end
  418. else
  419. [~,FilesExist] = ea_unifiedmapping_getvats(obj);
  420. end
  421. end
  422. while ~all(FilesExist(:))
  423. answ=questdlg('It seems like not all stimulation volumes have been calculated. We can initiate the process now, but this will take some time. Proceed?','Stimvolumes not calculated','yes','no','yes');
  424. switch answ
  425. case 'yes'
  426. if obj.calcsettings.calcmethod == 2 %Fiber based method
  427. obj.M.vatmodel='OSS-DBS (Butenko 2020)';
  428. switch obj.calcsettings.calcspace
  429. case 0
  430. space = 'native';
  431. case 1
  432. space = 'MNI';
  433. end
  434. end
  435. options=ea_defaultoptions;
  436. options.prefs.machine.vatsettings.butenko_calcPAM=0;
  437. options.prefs.machine.vatsettings.butenko_calcVAT=1;
  438. options.groupdir=fileparts(obj.leadgroup);
  439. if isfield(obj.M.ui, 'stimSetMode') && obj.M.ui.stimSetMode
  440. options.stimSetMode = 1;
  441. else
  442. options.stimSetMode = 0;
  443. end
  444. filesToCalc = find(sum(FilesExist(1:length(obj.M.patient.list),:),2)<2)';
  445. calc_biophysical(obj,options,filesToCalc);
  446. [~,FilesExist] = ea_unifiedmapping_getvats(obj);
  447. case 'no'
  448. return
  449. end
  450. end
  451. %recheck
  452. end
  453. return
  454. end
  455. function calculate_on_pam(obj,cfile)
  456. connid = (ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome));
  457. [pamlist,~] = ea_unifiedmapping_getpams(obj);
  458. [fibsvalBin, fibsvalprob,~, ~, ~, fibcell_pam, connFiberInd, totalFibers] = ea_unifiedmapping_calcvals_pam_prob(pamlist, obj, cfile);
  459. obj.results.fiberfiltering.(connid).('PAM_probA').fibsval = fibsvalprob;
  460. obj.results.fiberfiltering.(connid).connFiberInd_PAM = connFiberInd;
  461. obj.results.fiberfiltering.(connid).totalFibers = totalFibers; % total number of fibers in the connectome to work with global indices
  462. obj.results.fiberfiltering.(connid).('pam_fibers').fibcell= fibcell_pam;
  463. % temp. duplicate fibcell, will be fixed in the new explorer
  464. obj.results.fiberfiltering.(ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome)).fibcell = obj.results.fiberfiltering.(ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome)).('pam_fibers').fibcell;
  465. %add a provision for the results
  466. obj.results.fiberfiltering.(ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome)).calculationMethod = obj.calcsettings.calcmethod;
  467. end
  468. function calculate_on_efield(obj,cfile)
  469. connid = ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome);
  470. if isfield(obj.M,'pseudoM')
  471. vatlist=obj.M.ROI.list;
  472. [obj.customRoi.isbinary,obj.customRoi.minmax]=ea_unifiedmapping_checkcustomNii(vatlist);
  473. if obj.customRoi.isbinary
  474. obj.statsettings.stimulationmodel='VTA';
  475. end
  476. else
  477. [vatlist,~] = ea_unifiedmapping_getvats(obj);
  478. end
  479. [fibsvalBin, fibsvalSum, fibsvalMean, fibsvalPeak, fibsval5Peak, fibcell_efield, connFiberInd, totalFibers] = ea_fiberfiltering_calcvals(vatlist, cfile, obj.calcsettings.calcthreshold);
  480. obj.results.fiberfiltering.(connid).('VAT_Ttest').fibsval = fibsvalBin;
  481. obj.results.fiberfiltering.(connid).connFiberInd_VAT = connFiberInd; % old fiberfiltering files do not have these data and will fail when using pathway atlases
  482. obj.results.fiberfiltering.(connid).totalFibers = totalFibers; % total number of fibers in the connectome to work with global indices
  483. % only for e-fields
  484. obj.results.fiberfiltering.(connid).('efield_sum').fibsval = fibsvalSum;
  485. obj.results.fiberfiltering.(connid).('efield_mean').fibsval = fibsvalMean;
  486. obj.results.fiberfiltering.(connid).('efield_peak').fibsval = fibsvalPeak;
  487. obj.results.fiberfiltering.(connid).('efield_5peak').fibsval = fibsval5Peak;
  488. obj.results.fiberfiltering.(connid).('plainconn').fibsval = fibsvalBin;
  489. obj.results.fiberfiltering.(connid).('efield_fibers').fibcell= fibcell_efield;
  490. % temp. duplicate fibcell, will be fixed in the new explorer
  491. obj.results.fiberfiltering.(connid).fibcell = obj.results.fiberfiltering.(ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome)).('efield_fibers').fibcell;
  492. %add a provision for results
  493. obj.results.fiberfiltering.(connid).calculationMethod = obj.calcsettings.calcmethod;
  494. end
  495. function calculate_on_fibers(obj,cfile)
  496. connid = (ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome));
  497. % OSS-DBS E-field should be computed (not just warped!) in this space
  498. % get VAT list
  499. if isfield(obj.M,'pseudoM')
  500. vatlist = obj.M.ROI.list;
  501. else
  502. [vatlist,~] = ea_unifiedmapping_getlattice(obj);
  503. end
  504. % warp connectome to native space and compute E-field metrics
  505. ea_unified_get_Eproj(obj,vatlist)
  506. % define space again
  507. switch obj.calcsettings.calcspace
  508. case 0
  509. space = 'native';
  510. case 1
  511. space = 'MNI';
  512. end
  513. % load e-field projection metrics
  514. [fibsvalBin_proj, fibsvalSum_proj, fibsvalMean_proj, fibsvalPeak_proj, fibsval5Peak_proj, fibcell_proj, connFiberInd_proj,fibsvalBin_magn, fibsvalSum_magn, fibsvalMean_magn, fibsvalPeak_magn, fibsval5Peak_magn, fibcell_magn, connFiberInd_magn, totalFibers] = ea_unifiedmapping_native_calcvals(vatlist, cfile, space, obj);
  515. obj.results.fiberfiltering.(connid).totalFibers = totalFibers; % total number of fibers in the connectome to work with global indices
  516. obj.results.fiberfiltering.(connid).('VAT_Ttest').fibsval = fibsvalBin_magn;
  517. obj.results.fiberfiltering.(connid).('efield_sum').fibsval = fibsvalSum_magn;
  518. obj.results.fiberfiltering.(connid).('efield_mean').fibsval = fibsvalMean_magn;
  519. obj.results.fiberfiltering.(connid).('efield_peak').fibsval = fibsvalPeak_magn;
  520. obj.results.fiberfiltering.(connid).('efield_5peak').fibsval = fibsval5Peak_magn;
  521. obj.results.fiberfiltering.(connid).('plainconn').fibsval = fibsvalBin_magn;
  522. obj.results.fiberfiltering.(connid).('efield_fibers').fibcell = fibcell_magn;
  523. obj.results.fiberfiltering.(connid).('efield_fibers').connFiberInd_VAT = connFiberInd_magn; % old fiberfiltering files do not have these data and will fail when using pathway atlases
  524. obj.results.fiberfiltering.(connid).('VAT_Ttest_proj').fibsval = fibsvalBin_proj;
  525. obj.results.fiberfiltering.(connid).('efield_proj_sum').fibsval = fibsvalSum_proj;
  526. obj.results.fiberfiltering.(connid).('efield_proj_mean').fibsval = fibsvalMean_proj;
  527. obj.results.fiberfiltering.(connid).('efield_proj_peak').fibsval = fibsvalPeak_proj;
  528. obj.results.fiberfiltering.(connid).('efield_proj_5peak').fibsval = fibsval5Peak_proj;
  529. obj.results.fiberfiltering.(connid).('plainconn_proj').fibsval = fibsvalBin_proj;
  530. obj.results.fiberfiltering.(connid).('efield_proj').fibcell = fibcell_proj;
  531. obj.results.fiberfiltering.(connid).('efield_proj').connFiberInd_VAT = connFiberInd_proj; % old fiberfiltering files do not have these data and will fail when using pathway atlases
  532. obj.results.fiberfiltering.(connid).calculationMethod = 'Fiber Based Method';
  533. if strcmp(obj.e_field_metric,'Magnitude')
  534. obj.results.fiberfiltering.(connid).fibcell = obj.results.fiberfiltering.(ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome)).('efield_fibers').fibcell;
  535. obj.results.fiberfiltering.(connid).connFiberInd_VAT = obj.results.fiberfiltering.(ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome)).('efield_fibers').connFiberInd_VAT;
  536. else
  537. obj.results.fiberfiltering.(connid).fibcell = obj.results.fiberfiltering.(ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome)).('efield_proj').fibcell;
  538. obj.results.fiberfiltering.(connid).connFiberInd_VAT = obj.results.fiberfiltering.(ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome)).('efield_proj').connFiberInd_VAT;
  539. end
  540. end
  541. function recalculate_fiberfiltering_threshold(obj)
  542. % Recomputes fiber connectivity with current obj.calcsettings.calcthreshold.
  543. % Call this after changing the E-field threshold so results and drawing use
  544. % the new threshold. Applies when Fiber Filtering is selected and E-field
  545. % (voxel-based) or Fiber-based method is used.
  546. % After calling, the GUI should refresh stats and redraw (e.g. obj.draw()).
  547. if obj.calcsettings.selectedTool ~= 2
  548. ea_cprintf('CmdWinWarnings', 'recalculate_fiberfiltering_threshold: Fiber Filtering is not selected. No action.\n');
  549. return;
  550. end
  551. connid = ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome);
  552. if ~isfield(obj.results,'fiberfiltering') || ~isfield(obj.results.fiberfiltering, connid)
  553. ea_cprintf('CmdWinWarnings', 'No prior fiberfiltering results for current connectome. Run Calculate first.\n');
  554. return;
  555. end
  556. if obj.calcsettings.multi_pathways == 1
  557. [cfile, obj.map_list, obj.pathway_list] = ea_unifiedmapping_mergePathways(obj);
  558. else
  559. cfile = [ea_getconnectomebase('dMRI'), obj.calcsettings.fibfilt_connectome, filesep, 'data.mat'];
  560. end
  561. FilesExist = check_stimvols(obj);
  562. if ~all(FilesExist(:))
  563. ea_cprintf('CmdWinWarnings', 'Not all stimulation volumes exist. Recalculation aborted.\n');
  564. return;
  565. end
  566. switch obj.calcsettings.calcmethod
  567. case 1
  568. fprintf('Recalculating fiber connectivity with E-field threshold %g ...\n', obj.calcsettings.calcthreshold);
  569. calculate_on_efield(obj, cfile);
  570. case 2
  571. fprintf('Recalculating fiber connectivity (fiber-based) with E-field threshold %g ...\n', obj.calcsettings.calcthreshold);
  572. calculate_on_fibers(obj, cfile);
  573. otherwise
  574. ea_cprintf('CmdWinWarnings', 'recalculate_fiberfiltering_threshold: unsupported calcmethod. No action.\n');
  575. end
  576. end
  577. function results = calc_biophysical(obj,options,filesToCalc)
  578. for pt = filesToCalc
  579. [options.root, options.patientname] = fileparts(obj.M.patient.list{pt});
  580. options.root = [options.root, filesep];
  581. options = ea_getptopts(fullfile(options.root, options.patientname), options);
  582. fprintf('\nProcessing %s...\n\n', options.patientname);
  583. if ~isfield(obj.M,'S')
  584. ea_error(['Stimulation parameters for ', options.subj.subjId, ' are not set.']);
  585. end
  586. vfs = ea_regexpdir(ea_getearoot, 'ea_genvat_.*\.m$', 0);
  587. vfs = regexp(vfs, '(ea_genvat_.*)(?=\.m)', 'match', 'once');
  588. [vfnames,~,~] = cellfun(@(x) eval([x, '(''prompt'');']), vfs, 'Uni', 0);
  589. [~,ix]=ismember(obj.M.vatmodel,vfnames);
  590. try
  591. ea_genvat=eval(['@',vfs{ix}]);
  592. catch
  593. keyboard
  594. end
  595. if ~isfield(options.subj, 'norm')
  596. ea_cprintf('CmdWinWarnings', 'Running in Miniset mode: %s...\n', options.subj.subjId);
  597. volumespresent=0;
  598. elseif isempty(dir([options.subj.norm.transform.inverseBaseName, '*']))
  599. ea_cprintf('CmdWinWarnings', 'Tranformation not found for %s...\n', options.subj.subjId);
  600. volumespresent=0;
  601. else
  602. volumespresent=1;
  603. end
  604. options.orignative=options.native; % backup
  605. options.native=~ea_getprefs('vatsettings.estimateInTemplate'); % see whether VTAs should be directly estimated in template space or not
  606. if options.native && ~volumespresent
  607. ea_cprintf('CmdWinWarnings', 'Calculating VTA in template space since patient folder %s is incomplete.\n', options.subj.subjId);
  608. options.native=0;
  609. end
  610. if options.native % Reload native space coordinates
  611. coords = ea_load_reconstruction(options);
  612. else
  613. coords = obj.M.elstruct(pt).coords_mm;
  614. end
  615. if strcmp(obj.M.vatmodel, 'OSS-DBS (Butenko 2020)')
  616. if options.prefs.machine.vatsettings.butenko_calcAxonActivation
  617. feval(ea_genvat,obj.M.S(pt),options);
  618. ea_cprintf('CmdWinWarnings', 'OSS-DBS axon activation mode detect, skipping calc stats for %s!\n', options.patientname);
  619. continue;
  620. else
  621. [vatCalcPassed, ~] = feval(ea_genvat,obj.M.S(pt),options);
  622. end
  623. else
  624. for side=1:2
  625. try
  626. [vtafv,vtavolume] = feval(ea_genvat,coords,obj.M.S(pt),side,options,['gs_',obj.M.guid]);
  627. vatCalcPassed(side) = 1;
  628. catch
  629. vatCalcPassed(side) = 0;
  630. end
  631. if ~vatCalcPassed(side)
  632. ea_cprintf('CmdWinWarnings', 'VTA calculation failed for %s!\n', options.patientname);
  633. end
  634. end
  635. end
  636. options.native=options.orignative; % restore
  637. end
  638. end
  639. function Amps = getstimamp(obj)
  640. Amps=zeros(length(obj.M.patient.list),2);
  641. for pt=1:length(obj.M.patient.list)
  642. for side=1:2
  643. thisamp=obj.M.stats(pt).ea_stats.stimulation.vat(side).amp;
  644. thisamp(thisamp==0)=nan;
  645. Amps(pt,side)=ea_nanmean(thisamp');
  646. end
  647. end
  648. end
  649. function VTAvolumes = getvtavolumes(obj)
  650. if ~isfield(obj.M.stats(1).ea_stats.stimulation.vat(1),'volume')
  651. VTAvolumes = obj.getstimamp;
  652. warning('No VTA volumes found. Using stimulation amplitudes instead. Re-run stats in Lead-group to obtain volumes.');
  653. return
  654. end
  655. VTAvolumes=zeros(length(obj.M.patient.list),2);
  656. for pt=1:length(obj.M.patient.list)
  657. for side=1:2
  658. VTAvolumes(pt,side)=obj.M.stats(pt).ea_stats.stimulation.vat(side).volume;
  659. end
  660. end
  661. end
  662. function Efieldmags = getefieldmagnitudes(obj)
  663. if ~isfield(obj.M.stats(1).ea_stats.stimulation.efield(1),'volume')
  664. Efieldmags = obj.getstimamp;
  665. warning('No Efield magnitude sums found. Using stimulation amplitudes instead. Re-run stats in Lead-group to obtain values.');
  666. return
  667. end
  668. Efieldmags=zeros(length(obj.M.patient.list),2);
  669. for pt=1:length(obj.M.patient.list)
  670. for side=1:2
  671. try
  672. if isempty(obj.M.stats(pt).ea_stats.stimulation.efield(side).volume)
  673. val=0;
  674. else
  675. val=obj.M.stats(pt).ea_stats.stimulation.efield(side).volume;
  676. end
  677. Efieldmags(pt,side)=val;
  678. catch % could be efield(side) is not defined.
  679. Efieldmags(pt,side)=0;
  680. end
  681. end
  682. end
  683. end
  684. function refreshlg(obj)
  685. if ~exist(obj.leadgroup,'file')
  686. msgbox('Groupan alysis file has vanished. Please select file.');
  687. [fn,pth]=uigetfile();
  688. obj.leadgroup=fullfile(pth,fn);
  689. end
  690. U = load(obj.leadgroup);
  691. obj.M = U.M;
  692. obj.allpatients=obj.M.patient.list;
  693. end
  694. function coh = getcohortregressor(obj)
  695. coh=ea_cohortregressor(obj.M.patient.group(obj.patientselection));
  696. end
  697. function [I, Ihat] = loocv(obj,silent)
  698. if ~exist('silent','var')
  699. silent=0;
  700. end
  701. rng(obj.rngseed);
  702. cvp = cvpartition(length(obj.patientselection), 'LeaveOut');
  703. [I, Ihat] = crossval(obj, cvp,[],0,silent);
  704. end
  705. function [I, Ihat] = lococv(obj,silent)
  706. if length(unique(obj.M.patient.group(obj.patientselection))) == 1
  707. ea_error(sprintf(['Only one cohort in the analysis.\n', ...
  708. 'Leave-One-Cohort-Out-validation not possible.']));
  709. end
  710. [I, Ihat] = crossval(obj, obj.M.patient.group(obj.patientselection),[],0,silent);
  711. end
  712. function [I, Ihat, val_struct] = kfoldcv(obj,silent)
  713. if ~exist('silent','var')
  714. silent=0;
  715. end
  716. I_iter = {};
  717. Ihat_iter = {};
  718. rng(obj.rngseed);
  719. iter = obj.kIter;
  720. if iter == 1
  721. cvp = cvpartition(length(obj.patientselection),'KFold',obj.kfold);
  722. [I,Ihat, val_struct] = crossval(obj,cvp,[],0,silent);
  723. else
  724. % plot some statistics over shuffles
  725. r_over_iter = zeros(iter,1);
  726. p_over_iter = zeros(iter,1);
  727. for i=1:iter
  728. cvp = cvpartition(length(obj.patientselection), 'KFold', obj.kfold);
  729. if ~silent
  730. fprintf("Iterating fold set: %d",i)
  731. end
  732. [I_iter{i}, Ihat_iter{i},val_struct] = crossval(obj, cvp, [], 1,silent);
  733. if ~silent
  734. switch obj.multitractmode
  735. case 'Split & Color By PCA'
  736. disp("Fold Agreement is not evaluated for PCA")
  737. otherwise
  738. inx_nnan = find(isnan(I_iter{i}) ~= 1);
  739. [r_over_iter(i),p_over_iter(i)]=ea_permcorr(I_iter{i}(inx_nnan),Ihat_iter{i}(inx_nnan),'spearman');
  740. end
  741. end
  742. end
  743. % check model agreement over shuffles using Sequential Rank Agreement
  744. % disabled for PCA
  745. switch obj.multitractmode
  746. case 'Split & Color By PCA'
  747. if ~silent
  748. disp("Fold Agreement is not evaluated for PCA")
  749. end
  750. otherwise
  751. if ~silent
  752. r_Ihat = zeros(size(Ihat_iter,2));
  753. for i = 1:size(r_Ihat,1)
  754. for j = 1:size(r_Ihat,1)
  755. [r_Ihat(i,j),~]=ea_permcorr(Ihat_iter{i},Ihat_iter{j},'spearman');
  756. end
  757. end
  758. % plot correlation matrix
  759. figure('Name','Patient scores'' correlations','Color','w','NumberTitle','off')
  760. imagesc(triu(r_Ihat));
  761. title('Patient scores'' correlations over K-fold shuffles', 'FontSize', 16); % set title
  762. colormap('bone');
  763. cb = colorbar;
  764. set(cb)
  765. % plot r-vals over shuffles
  766. p_above_05 = p_over_iter(find(p_over_iter>0.05),:);
  767. p_above_01 = p_over_iter(find(p_over_iter>0.01),:);
  768. h = figure('Name','Over-fold analysis','Color','w','NumberTitle','off');
  769. g = ea_raincloud_plot(r_over_iter,'box_on',1);
  770. a1=gca;
  771. set(a1,'ytick',[])
  772. a1.XLabel.String='Spearman''s R of model and clinical scores';
  773. if min(r_over_iter) >= -0.9
  774. r_lower_lim = min(r_over_iter) - 0.1;
  775. else
  776. r_lower_lim = -1.0;
  777. end
  778. if max(r_over_iter) <= 0.9
  779. r_upper_lim = max(r_over_iter) + 0.1;
  780. else
  781. r_upper_lim = 1.0;
  782. end
  783. a1.XLim=([r_lower_lim r_upper_lim]);
  784. text(0.25,0.9,['N(p>0.05) = ',sprintf('%d',length(p_above_05))],'FontWeight','bold','FontSize',14,'HorizontalAlignment','right','Units','normalized');
  785. text(0.25,0.83,['N(p>0.01) = ',sprintf('%d',length(p_above_01))],'FontWeight','bold','FontSize',14,'HorizontalAlignment','right','Units','normalized');
  786. end
  787. end
  788. % we should think about this part
  789. I_iter = cell2mat(I_iter);
  790. Ihat_iter = cell2mat(Ihat_iter);
  791. I = mean(I_iter,2,'omitnan');
  792. Ihat = mean(Ihat_iter,2,'omitnan');
  793. end
  794. end
  795. function [I, Ihat, val_struct] = lno(obj, Iperm, silent)
  796. if ~exist('silent','var')
  797. silent=0;
  798. end
  799. rng(obj.rngseed);
  800. cvp = cvpartition(length(obj.patientselection), 'resubstitution');
  801. if ~exist('Iperm', 'var') || isempty(Iperm)
  802. [I, Ihat, val_struct] = crossval(obj, cvp, [], [], silent);
  803. else
  804. [I, Ihat, val_struct] = crossval(obj, cvp, Iperm, [], silent);
  805. end
  806. end
  807. function [Improvement, Ihat, val_struct] = crossval(obj, cvp, Iperm, shuffle, silent)
  808. if ~exist('silent','var')
  809. silent=0;
  810. end
  811. if ~exist('shuffle','var') || isempty(shuffle)
  812. shuffle=0;
  813. end
  814. if isnumeric(cvp) % cvp is crossvalind
  815. cvIndices = cvp;
  816. cvID = unique(cvIndices);
  817. cvp = struct;
  818. cvp.NumTestSets = length(cvID);
  819. for i=1:cvp.NumTestSets
  820. cvp.training{i} = cvIndices~=cvID(i);
  821. cvp.test{i} = cvIndices==cvID(i);
  822. end
  823. end
  824. % Check if patients are selected in the custom training/test list
  825. if isempty(obj.customselection)
  826. patientsel = obj.patientselection;
  827. else
  828. patientsel = obj.customselection;
  829. end
  830. switch obj.multitractmode
  831. case 'Split & Color By PCA'
  832. if ~exist('Iperm', 'var') || isempty(Iperm)
  833. %Improvement = obj.subscore.vars;
  834. for i=1:length(obj.subscore.vars)
  835. Improvement{i} = obj.subscore.vars{i}(patientsel);
  836. end
  837. else
  838. for i=1:length(obj.subscore.vars)
  839. Improvement{i} = Iperm(patientsel,i);
  840. end
  841. end
  842. otherwise
  843. if ~exist('Iperm', 'var') || isempty(Iperm)
  844. Improvement = obj.responsevar(patientsel,:);
  845. else
  846. Improvement = Iperm(patientsel,:);
  847. end
  848. end
  849. % Ihat is the estimate of improvements (not scaled to real improvements)
  850. if strcmp(obj.multitractmode,'Single Tract Analysis')
  851. Ihat = nan(length(patientsel),2);
  852. Ihat_train_global = nan(cvp.NumTestSets,length(patientsel),2);
  853. else
  854. Ihat = nan(length(patientsel),2,length(obj.subscore.vars));
  855. Ihat_train_global = nan(cvp.NumTestSets,length(patientsel),2,length(obj.subscore.vars));
  856. end
  857. if strcmp(obj.drawTool,'fiberfiltering')
  858. if obj.useExternalModel == true && ~strcmp(obj.ExternalModelFile, 'None')
  859. S = load(obj.ExternalModelFile);
  860. if ~strcmp(ea_unifiedmapping_method2methodid(obj),S.fibsvalType)
  861. waitfor(msgbox('Change the Model Setup! See terminal'));
  862. disp('The loaded model uses: ')
  863. disp(S.fibsvalType)
  864. end
  865. fibsval = full(obj.results.fiberfiltering.(ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome)).(S.fibsvalType).fibsval);
  866. else
  867. fibsval = full(obj.results.fiberfiltering.(ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome)).(ea_unifiedmapping_method2methodid(obj)).fibsval);
  868. end
  869. else
  870. fibsval = {};
  871. end
  872. % for nested LOO, store some statistics
  873. if obj.nestedLOO
  874. Abs_pred_error = zeros(cvp.NumTestSets, 1);
  875. Predicted_scores = zeros(length(patientsel), 1);
  876. Slope = zeros(cvp.NumTestSets, 1);
  877. Intercept = zeros(cvp.NumTestSets, 1);
  878. end
  879. for c=1:cvp.NumTestSets
  880. if cvp.NumTestSets ~= 1
  881. if ~silent
  882. fprintf(['\nIterating set: %0',num2str(numel(num2str(cvp.NumTestSets))),'d/%d\n'], c, cvp.NumTestSets);
  883. end
  884. end
  885. if isobject(cvp)
  886. training = cvp.training(c);
  887. test = cvp.test(c);
  888. elseif isstruct(cvp)
  889. training = cvp.training{c};
  890. test = cvp.test{c};
  891. end
  892. % now do LOO within the training group
  893. if obj.nestedLOO
  894. % use all patients, but outer loop left-out is always 0
  895. if strcmp(obj.multitractmode,'Single Tract Analysis')
  896. Ihat_inner = nan(length(patientsel),2);
  897. Ihat_train_global_inner = nan(cvp.NumTestSets,length(patientsel),2);
  898. else
  899. Ihat_inner = nan(length(patientsel),2,length(obj.subscore.vars));
  900. Ihat_train_global_inner = nan(cvp.NumTestSets,length(patientsel),2,length(obj.subscore.vars));
  901. end
  902. for test_i = 1:length(training)
  903. training_inner = training;
  904. training_inner(test_i) = 0;
  905. % check if inner and outer left-out match
  906. if all(training_inner == training)
  907. continue
  908. end
  909. test_inner = logical(zeros(length(training), 1));
  910. test_inner(test_i) = logical(training(test_i));
  911. % updates Ihat_inner(test_inner)
  912. if ~exist('Iperm', 'var') || isempty(Iperm)
  913. [Ihat_inner, ~, ~] = ea_compute_unified_model(c,obj, fibsval, Ihat_inner, Ihat_train_global_inner, patientsel, training_inner, test_inner);
  914. else
  915. [Ihat_inner, ~, ~] = ea_compute_unified_model(c, obj, fibsval, Ihat_inner, Ihat_train_global_inner, patientsel, training_inner, test_inner,Iperm);
  916. end
  917. end
  918. % fit the linear model based on inner loop fibscores
  919. predictor=squeeze(ea_nanmean(Ihat_inner,2));
  920. % iterating over all test_inner gives us training
  921. mdl=fitglm(predictor(training),Improvement(training),lower(obj.predictionmodel));
  922. Intercept(c) = mdl.Coefficients.Estimate(1);
  923. Slope(c) = mdl.Coefficients.Estimate(2);
  924. end
  925. % now compute Ihat for the true 'test' left out
  926. % updates Ihat(test)
  927. if ~exist('Iperm', 'var') || isempty(Iperm)
  928. [Ihat, Ihat_train_global, val_struct{c}] = ea_compute_unified_model(c,obj, fibsval, Ihat, Ihat_train_global, patientsel, training, test);
  929. else
  930. [Ihat, Ihat_train_global, val_struct{c}] = ea_compute_unified_model(c,obj, fibsval, Ihat, Ihat_train_global, patientsel, training, test, Iperm);
  931. end
  932. % predict the improvement in the left-out patient (fold) of
  933. % the outer loop
  934. if obj.nestedLOO
  935. predictor=squeeze(ea_nanmean(Ihat,2));
  936. Ihat_voters_prediction = repmat(predict(mdl,predictor(test)),1,2);
  937. %Abs_pred_error(c) = abs(Improvement(test) - Ihat_voters_prediction(test));
  938. Predicted_scores(test) = Ihat_voters_prediction(1:end,1); % only one value here atm
  939. end
  940. end
  941. % check if binary variable and not permutation test
  942. if (~exist('Iperm', 'var') || isempty(Iperm)) && all(ismember(Improvement(:,1), [0,1])) && size(val_struct{c}.vals,1) == 1
  943. % average across sides. This might be wrong for capsular response.
  944. Ihat_av_sides = ea_nanmean(Ihat,2);
  945. if isobject(cvp)
  946. % In-sample
  947. AUC = ea_logit_regression(0 ,Ihat_av_sides, Improvement, 1:size(Improvement,1), 1:size(Improvement,1));
  948. elseif isstruct(cvp)
  949. % actual training and test
  950. Ihat_train_global_av_sides = ea_nanmean(Ihat_train_global,3); % in this case, dimens is (1, N, sides)
  951. AUC = ea_logit_regression(Ihat_train_global_av_sides(training)', Ihat_av_sides, Improvement, training, test);
  952. end
  953. end
  954. if ~silent
  955. % plot patient score correlation matrix over folds
  956. if (~exist('shuffle', 'var')) || shuffle == 0 || isempty(shuffle)
  957. if cvp.NumTestSets ~= 1 && (strcmp(obj.multitractmode,'Single Tract Analysis') || strcmp(obj.multitractmode,'Single Tract Analysis Button'))
  958. % put training and test scores together
  959. Ihat_combined = cell(1,cvp.NumTestSets);
  960. %Ihat_combined = Ihat_train_global;
  961. for c=1:cvp.NumTestSets
  962. if isobject(cvp)
  963. training = cvp.training(c);
  964. test = cvp.test(c);
  965. elseif isstruct(cvp)
  966. training = cvp.training{c};
  967. test = cvp.test{c};
  968. end
  969. Ihat_combined{c}(training,1) = Ihat_train_global(c,training,1)';
  970. Ihat_combined{c}(test,1) = Ihat(test,1);
  971. end
  972. r_Ihat = zeros(size(Ihat_combined,2));
  973. for i = 1:size(r_Ihat,1)
  974. for j = 1:size(r_Ihat,1)
  975. [r_Ihat(i,j),~]=ea_permcorr(Ihat_combined{i},Ihat_combined{j},'spearman');
  976. end
  977. end
  978. figure('Name','Patient scores'' correlations','Color','w','NumberTitle','off')
  979. imagesc(triu(r_Ihat)); % Display correlation matrix as an image
  980. title('Patient scores'' correlations over folds', 'FontSize', 16); % set title
  981. colormap('bone');
  982. cb = colorbar;
  983. % set(cb)
  984. end
  985. end
  986. end
  987. if obj.nestedLOO
  988. % cvs = 'L-O-O-O';
  989. % h = ea_corrbox(Improvement,Predicted_dif_models,'permutation',{['Disc. Fiber prediction ',upper(cvs)],empiricallabel,fibscorelabel});
  990. LM_values_slope = [num2str(mean(Slope)) ' ' char(177) ' ' num2str(std(Slope))];
  991. LM_values_intercept = [num2str(mean(Intercept)) ' ' char(177) ' ' num2str(std(Intercept))];
  992. disp('Mean and STD for slopes and intercepts of LMs')
  993. disp(LM_values_slope)
  994. disp(LM_values_intercept)
  995. % visualize lms and CIs for 5-fold or less
  996. if cvp.NumTestSets < 6
  997. groups_nested = zeros(length(Predicted_scores),1);
  998. for group_idx = 1:cvp.NumTestSets
  999. groups_nested(cvp.test(group_idx)) = group_idx;
  1000. end
  1001. side = 1;
  1002. plotName = 'Fitting of linear models for K-folds using nested LOO';
  1003. empiricallabel = 'Empirical score';
  1004. pred_label = 'Predicted score';
  1005. h=ea_corrbox(Improvement,Predicted_scores,'permutation',{['Disc. Fiber prediction ',plotName],empiricallabel,pred_label, plotName, LM_values_slope, LM_values_intercept},groups_nested);
  1006. end
  1007. end
  1008. if obj.doactualprediction % repeat loops partly to fit to actual response variables:
  1009. Ihat_voters_prediction=nan(size(Ihat));
  1010. %add some warnings
  1011. switch obj.multitractmode
  1012. case 'Single Tract Analysis'
  1013. if obj.useExternalModel && size(val_struct{c}.vals,1) > 1
  1014. ea_error("You can only use the Fit-to-Score feature with a Single Tract Analysis analysis model");
  1015. end
  1016. otherwise
  1017. if obj.useExternalModel
  1018. ea_error("You can only use the Fit-to-Score feature with Single Tract Analysis");
  1019. end
  1020. end
  1021. numVoters = size(val_struct{c}.vals,1);
  1022. for c=1:cvp.NumTestSets
  1023. if isobject(cvp)
  1024. training = cvp.training(c);
  1025. test = cvp.test(c);
  1026. elseif isstruct(cvp)
  1027. training = cvp.training{c};
  1028. test = cvp.test{c};
  1029. end
  1030. for voter=1:numVoters
  1031. switch obj.multitractmode
  1032. case 'Split & Color By Subscore'
  1033. if ~exist('Iperm', 'var') || isempty(Iperm)
  1034. useI=obj.subscore.vars{voter}(patientsel);
  1035. else
  1036. % to be added by Nanditha
  1037. end
  1038. case 'Split & Color By PCA'
  1039. if ~exist('Iperm', 'var') || isempty(Iperm)
  1040. useI=obj.subscore.pcavars{voter}(patientsel);
  1041. else
  1042. PCscores = ea_nanzscore(Iperm(patientsel, : ))*obj.subscore.pcacoeff;
  1043. useI = PCscores(:, voter);
  1044. end
  1045. otherwise
  1046. if ~exist('Iperm', 'var') || isempty(Iperm)
  1047. useI=obj.responsevar(patientsel);
  1048. else
  1049. useI=Iperm(patientsel);
  1050. end
  1051. end
  1052. if size(useI,2)>1
  1053. ea_error('This has not been implemented for hemiscores.');
  1054. end
  1055. % these predictors are defined within the same fiberfiltering model
  1056. % of iteration 'c'
  1057. % do not get rid of the first dimension when it has size of 1
  1058. Ihat_train_global_av_sides = ea_nanmean(Ihat_train_global,3);
  1059. predictor_training = reshape(Ihat_train_global_av_sides, ...
  1060. size(Ihat_train_global_av_sides,1),...
  1061. size(Ihat_train_global_av_sides,2),...
  1062. size(Ihat_train_global_av_sides,4));
  1063. predictor_test = squeeze(ea_nanmean(Ihat,2));
  1064. %predictor=squeeze(ea_nanmean(Ihat_voters,2));
  1065. covariates=[];
  1066. for cv = 1:length(obj.covars)
  1067. covariates = [covariates,obj.covars{cv}(patientsel)];
  1068. end
  1069. if obj.useExternalModel == true %only use for single tract analysis
  1070. if ~strcmp(obj.multitractmode,'Single Tract Analysis')
  1071. ea_error("Sorry, you cannot use exported model and fit-to-scores for multi-tract model");
  1072. else
  1073. mdl = S.mdl;
  1074. end
  1075. else
  1076. if ~isempty(covariates)
  1077. mdl=fitglm([predictor_training(c,training,voter)',covariates(training,:)],useI(training),lower(obj.predictionmodel));
  1078. else
  1079. mdl=fitglm([predictor_training(c,training,voter)],useI(training),lower(obj.predictionmodel));
  1080. end
  1081. end
  1082. if size(useI,2) == 1 % global scores
  1083. if ~isempty(covariates)
  1084. Ihat_voters_prediction(test,:,voter)=repmat(predict(mdl,[predictor_test(test,voter),covariates(test,:)]),1,2); % fill both sides equally
  1085. else
  1086. Ihat_voters_prediction(test,:,voter)=repmat(predict(mdl,[predictor_test(test,voter)]),1,2); % fill both sides equally
  1087. end
  1088. elseif size(useI,2)==2 % bihemispheric scores
  1089. ea_error('Fitting to scores has not been implemented for bihemispheric scores.');
  1090. end
  1091. end
  1092. end
  1093. % quantify the prediction accuracy (if Train-Test)
  1094. if cvp.NumTestSets == 1 && voter == 1 && size(obj.responsevar,2) == 1 && (~exist('Iperm', 'var') || isempty(Iperm))
  1095. side = 1;
  1096. SS_tot = var(useI(test)) * (length(useI(test)) - 1); % just a trick to use one line
  1097. SS_res = sum((Ihat_voters_prediction(test,side,1) - useI(test)).^2);
  1098. R2 = 1 - SS_res/SS_tot;
  1099. RMS = sqrt(mean((Ihat_voters_prediction(test,side,1) - useI(test)).^2));
  1100. MAD = median(abs(Ihat_voters_prediction(test,side,1) - useI(test)));
  1101. MAE = mean(abs(Ihat_voters_prediction(test,side,1) - useI(test)));
  1102. plotName = 'TRAIN-TEST';
  1103. R2_label = ['R2 = ', sprintf('%.3f',R2)];
  1104. RMS_label = ['RMS = ', sprintf('%.3f',RMS)];
  1105. MAD_label = ['MAD = ', sprintf('%.3f',MAD)];
  1106. empiricallabel = 'Empirical score';
  1107. pred_label = 'Predicted score';
  1108. h = ea_corrbox(useI(test),Ihat_voters_prediction(test,side,1),'permutation',{['Disc. Fiber prediction ',plotName],empiricallabel,pred_label, plotName, R2_label, RMS_label, MAD_label});
  1109. % h2 = ea_corrbox(-1*useI(test),Ihat_voters_prediction(test,side,1),'permutation',{['Disc. Fiber prediction ',plotName],empiricallabel,pred_label, plotName, R2_label, RMS_label, MAD_label});
  1110. end
  1111. Ihat=Ihat_voters_prediction; % replace with actual response variables.
  1112. end
  1113. switch obj.multitractmode
  1114. case 'Split & Color By Subscore'
  1115. if ~obj.CleartuneOptim
  1116. % here we map back to the single response variable using a
  1117. % weightmatrix
  1118. if isempty(obj.customselection)
  1119. selected_pts = obj.patientselection;
  1120. else
  1121. selected_pts = obj.customselection;
  1122. end
  1123. weightmatrix=zeros(size(Ihat));
  1124. for voter=1:size(Ihat,3)
  1125. if ~isnan(obj.subscore.weights(voter)) % same weight for all subjects in that voter (slider was used)
  1126. weightmatrix(:,:,voter)=obj.subscore.weights(voter);
  1127. else % if the weight value is nan, this means we will need to derive a weight from the variable of choice
  1128. weightmatrix(:,:,voter)=repmat(ea_minmax(obj.subscore.weightvars{voter}(selected_pts)),1,size(weightmatrix,2)/size(obj.subscore.weightvars{voter}(selected_pts),2));
  1129. weightmatrix(:,:,voter)=weightmatrix(:,:,voter)./max(obj.subscore.weightvars{voter}(selected_pts)); % weight for unnormalized data across voters *
  1130. % *) e.g. in case one symptom - bradykinesia -
  1131. % has a max of 20, while a second - tremor -
  1132. % will have a max of 5, we want to equilize
  1133. % those. We use minmax() in the line above to
  1134. % get rid of negative values and use
  1135. % ./ea_nansum below to take the average.
  1136. end
  1137. end
  1138. for xx=1:size(Ihat,1) % make sure voter weights sum up to 1
  1139. for yy=1:size(Ihat,2)
  1140. % for xx=1:size(Ihat_voters,1) % make sure voter weights sum up to 1
  1141. % for yy=1:size(Ihat_voters,2)
  1142. weightmatrix(xx,yy,:)=weightmatrix(xx,yy,:)./ea_nansum(weightmatrix(xx,yy,:));
  1143. end
  1144. end
  1145. Ihat=ea_nansum(Ihat.*weightmatrix,3);
  1146. else
  1147. Ihat = Ihat(test,:,:);
  1148. Ihat = reshape(Ihat,2,length(obj.subscore.vars))';
  1149. Improvement = Improvement(test);
  1150. return;
  1151. end
  1152. case 'Split & Color By PCA'
  1153. Ihat=squeeze(ea_nanmean(Ihat,2));
  1154. %Ihat_voters=squeeze(ea_nanmean(Ihat_voters,2)); % need to assume global scores here for now.
  1155. % map back to PCA:
  1156. for i=1:length(obj.subscore.vars)
  1157. selected_subscores{i} = obj.subscore.vars{i}(patientsel);
  1158. end
  1159. subvars=ea_nanzscore(cell2mat(selected_subscores));
  1160. if size(subvars,2) <= 2
  1161. ea_warndlg("You may not have enough subscores & this might result in errors. Please consider selecting more subscores.")
  1162. end
  1163. % [coeff,score,latent,tsquared,explained,mu]=pca(subvars,'Rows','complete');
  1164. % use saved weights to ensure consistency
  1165. coeff = obj.subscore.pcacoeff;
  1166. if ~silent
  1167. % show predictions for PC scores
  1168. if ~exist('Iperm', 'var') || isempty(Iperm) % avoid plotting for each permutation if using permutations!
  1169. for pcc=1:obj.numpcs
  1170. if obj.subscore.posvisible(pcc)==1 || obj.subscore.negvisible(pcc)==1 % don't try to plot if not showing any fibers for this PC
  1171. ea_corrplot(obj.subscore.pcavars{pcc}(patientsel),Ihat(:,pcc), 'noperm', ...
  1172. {['Disc. Fiber prediction for PC ',num2str(pcc)],'PC score (Empirical)','PC score (Predicted)'},...
  1173. [], [], obj.subscore.pcacolors(pcc, :));
  1174. % sum(obj.subscore.pcavars{pcc}(obj.patientselection) - score(:,pcc)) % quick check
  1175. end
  1176. end
  1177. end
  1178. end
  1179. % data is zscored, such as mu is 0 (+ some computer rounding error)
  1180. % then adding mean is not required
  1181. % also, we want to take scores of the chosen PCs ONLY,
  1182. % and multiply by coeff of these PCs (= how they map to
  1183. % the variables) to get estimated clinical scores
  1184. Ihatout = Ihat(:,1:obj.numpcs)*coeff(:,1:obj.numpcs)';
  1185. %Ihatout = Ihat*coeff(:,1:obj.numpcs)' + repmat(mu,size(score,1),1);
  1186. %Ihatout = Ihat_voters*coeff(:,1:obj.numpcs)' + repmat(mu,size(score,1),1);
  1187. Ihat = mat2cell(Ihatout, size(Ihatout,1), ones(1,length(obj.subscore.vars)));
  1188. otherwise
  1189. Ihat=squeeze(Ihat);
  1190. %Ihat=squeeze(Ihat_voters);
  1191. end
  1192. if ~iscell(Ihat)
  1193. if cvp.NumTestSets == 1
  1194. Ihat = Ihat(test,:);
  1195. Improvement = Improvement(test);
  1196. end
  1197. if size(obj.responsevar,2)==2 % hemiscores
  1198. Ihat = Ihat(:); % compare hemiscores (electrode wise)
  1199. Improvement = Improvement(:);
  1200. else
  1201. Ihat = ea_nanmean(Ihat,2); % compare bodyscores (patient wise)
  1202. end
  1203. end
  1204. % restore original view in case of live drawing
  1205. if obj.cvlivevisualize
  1206. obj.draw;
  1207. end
  1208. end
  1209. function [Iperm, Ihat, R0, R1, pperm, Rp95, val_struct] = lnopb(obj, corrType, silent)
  1210. if ~exist('corrType', 'var')
  1211. corrType = 'Spearman';
  1212. end
  1213. if ~exist('silent','var')
  1214. silent=0;
  1215. end
  1216. numPerm = obj.Nperm;
  1217. if strcmp(obj.multitractmode,'Split & Color By PCA')
  1218. Iperm = ea_shuffle(cell2mat(obj.subscore.vars'), numPerm, obj.patientselection, obj.rngseed);
  1219. Iperm(2:numPerm+1,:,:) = Iperm;
  1220. Iperm(1,:,:) = cell2mat(obj.subscore.vars');
  1221. Ihat = cell(numPerm+1,1);
  1222. R = zeros(numPerm+1, length(obj.subscore.vars));
  1223. for perm=1:numPerm+1
  1224. if perm==1
  1225. if ~silent; fprintf('Calculating without permutation\n\n'); end
  1226. [~, Ihat{perm},thisval_struct] = lno(obj, [], silent);
  1227. else
  1228. if ~silent; fprintf('Calculating permutation: %d/%d\n\n', perm-1, numPerm); end
  1229. [~, Ihat{perm},thisval_struct] = lno(obj, squeeze(Iperm(perm,:,:)), silent);
  1230. end
  1231. val_struct{perm}=thisval_struct{1};
  1232. for subvar = 1:length(obj.subscore.vars)
  1233. R(perm,subvar) = corr(Iperm(perm, obj.patientselection, subvar)',...
  1234. Ihat{perm}{subvar},'type',corrType,'rows','pairwise');
  1235. end
  1236. end
  1237. R(isnan(R)) = 1e-5; % do not get rid of Nans
  1238. % generate null distribution
  1239. R1 = R(1,:);
  1240. for subvar = 1:length(obj.subscore.vars)
  1241. R0(:,subvar) = sort(R(2:end,subvar), 'descend');
  1242. Rp95(subvar) = R0(round(0.05*numPerm),subvar);
  1243. pperm(subvar) = mean(abs(R0(:,subvar))>=abs(R1(subvar)));
  1244. if ~silent; fprintf(['Permuted p for ' obj.subscore.labels{subvar} ' = ' num2str(pperm(subvar)) '.\n']); end
  1245. end
  1246. % Return only selected I
  1247. Iperm = Iperm(:,obj.patientselection,:);
  1248. else % any mode except PCA
  1249. Iperm = ea_shuffle(obj.responsevar, numPerm, obj.patientselection, obj.rngseed)';
  1250. Iperm = [obj.responsevar, Iperm];
  1251. Ihat = cell(numPerm+1, 1);
  1252. R = zeros(numPerm+1, 1);
  1253. for perm=1:numPerm+1
  1254. if perm==1
  1255. if ~silent; fprintf('Calculating without permutation\n\n'); end
  1256. [~, Ihat{perm},val_struct{perm}] = lno(obj, [], silent);
  1257. else
  1258. if ~silent; fprintf('Calculating permutation: %d/%d\n\n', perm-1, numPerm); end
  1259. [~, Ihat{perm},val_struct{perm}] = lno(obj, Iperm(:, perm), silent);
  1260. end
  1261. R(perm) = corr(Iperm(obj.patientselection,perm),Ihat{perm},'type',corrType,'rows','pairwise');
  1262. end
  1263. R(isnan(R)) = 1e-5;
  1264. % generate null distribution
  1265. R1 = R(1);
  1266. R0 = sort((R(2:end)),'descend');
  1267. Rp95 = R0(round(0.05*numPerm));
  1268. pperm = mean(abs(R0)>=abs(R1));
  1269. if ~silent; disp(['Permuted p = ',sprintf('%0.2f',pperm),'.']); end
  1270. % Return only selected I
  1271. Iperm = Iperm(obj.patientselection,:);
  1272. end
  1273. end
  1274. function save(obj)
  1275. % Create a temporary object with only the required fields
  1276. explorer = ea_unifiedmapping;
  1277. Incprops = {'results','calcsettings','leadgroup','ID','M'};
  1278. for i = 1:length(Incprops)
  1279. explorer.(Incprops{i}) = obj.(Incprops{i});
  1280. end
  1281. % Get all properties of the object
  1282. %This is necessary to match the settings file
  1283. %only save results in this
  1284. if isempty(obj.analysispath)
  1285. [pth,~,~] = fileparts(obj.leadgroup);
  1286. obj.analysispath=[pth,filesep,'UnifiedMappingExplorer',filesep,obj.ID,'.explorer'];
  1287. ea_mkdir([pth,filesep,'UnifiedMappingExplorer']);
  1288. end
  1289. rf=obj.resultfig; % need to stash fig handle for saving.
  1290. rd=obj.drawobject; % need to stash handle of drawing before saving.
  1291. try % could be figure is already closed.
  1292. setappdata(rf,['dt_',explorer.ID],rd); % store handle of tract to figure.
  1293. end
  1294. save(obj.analysispath,'explorer','-v7.3');
  1295. saveObjectToJson(obj);
  1296. obj.resultfig=rf;
  1297. obj.drawobject=rd;
  1298. end
  1299. function saveObjectToJson(obj)
  1300. % Convert object to a struct (including nested objects)
  1301. voxtractsettings = objectToStruct(obj);
  1302. % Convert struct to JSON
  1303. jsonStr = jsonencode(voxtractsettings, 'PrettyPrint', true);
  1304. %define filepaths
  1305. if isempty(obj.analysispath)
  1306. [DBSMappingfolder,~,~] = fileparts(obj.leadgroup);
  1307. else
  1308. [DBSMappingfolder,~,~] = fileparts(obj.analysispath);
  1309. end
  1310. conn_val = 'default';
  1311. switch obj.drawTool
  1312. case 'sweetspotmapping'
  1313. conn_val = 'default';
  1314. case 'fiberfiltering'
  1315. conn_val = ea_unifiedmapping_conn2connid(obj.calcsettings.fibfilt_connectome);
  1316. case 'networkmapping'
  1317. conn_val = ea_unifiedmapping_conn2connid(obj.calcsettings.netmap_connectome);
  1318. end
  1319. if ~isfolder(DBSMappingfolder)
  1320. ea_mkdir(DBSMappingfolder)
  1321. end
  1322. jsonPath=[DBSMappingfolder,filesep,'Settings-',obj.ID,'_conn-',conn_val,'.json'];
  1323. % Write JSON to a file
  1324. fileID = fopen(jsonPath, 'w');
  1325. if fileID == -1
  1326. error('Cannot open file for writing.');
  1327. end
  1328. fprintf(fileID, '%s', jsonStr);
  1329. fclose(fileID);
  1330. end
  1331. function s = objectToStruct(obj)
  1332. % Convert an object to a struct, handling nested objects
  1333. if nargin < 2
  1334. ignoreList = {'results','resultfig','drawobject','M'}; %M should be present in the explorer mat file. This is because there are some complicated structures in M files that are not well translated in struct (for json encoding). % Default: Do not ignore any properties unless specified
  1335. end
  1336. if isobject(obj)
  1337. props = properties(obj);
  1338. s = struct();
  1339. for i = 1:length(props)
  1340. propName = props{i};
  1341. propValue = obj.(props{i});
  1342. if ismember(propName, ignoreList) %skip some of the properties.
  1343. continue;
  1344. end
  1345. if isa(propValue, 'matlab.ui.Figure') || isa(propValue, 'handle') %also skip handles
  1346. continue;
  1347. end
  1348. if isobject(propValue) % Recursively convert nested objects
  1349. s.(props{i}) = objectToStruct(propValue);
  1350. else
  1351. s.(props{i}) = propValue;
  1352. end
  1353. end
  1354. else
  1355. s = obj; % If it's not an object, return as is (handles arrays, numbers, strings)
  1356. end
  1357. end
  1358. function draw(obj)
  1359. if ~isfield(obj.activated,'sweetspotmapping')
  1360. obj.activated.sweetspotmapping='Off';
  1361. end
  1362. if ~isfield(obj.activated,'fiberfiltering')
  1363. obj.activated.fiberfiltering='Off';
  1364. end
  1365. if ~isfield(obj.activated,'networkmapping')
  1366. obj.activated.networkmapping='Off';
  1367. end
  1368. % delete prior spots:
  1369. if isfield(obj.drawobject,'sweetspotmapping')
  1370. if ~isempty(obj.drawobject.sweetspotmapping)
  1371. for s=1:numel(obj.drawobject.sweetspotmapping)
  1372. for ins=1:numel(obj.drawobject.sweetspotmapping{s})
  1373. try delete(obj.drawobject.sweetspotmapping{s}{ins}.toggleH); end
  1374. try delete(obj.drawobject.sweetspotmapping{s}{ins}.patchH); end
  1375. try delete(obj.drawobject.sweetspotmapping{s}{ins}); end
  1376. end
  1377. end
  1378. end
  1379. end
  1380. % plot new spots
  1381. switch lower(obj.activated.sweetspotmapping)
  1382. case 'on'
  1383. obj.drawTool='sweetspotmapping';
  1384. ea_unified_draw(obj);
  1385. end
  1386. % delete prior tracts
  1387. if isfield(obj.drawobject,'fiberfiltering')
  1388. if ~isempty(obj.drawobject.fiberfiltering)
  1389. for s=1:numel(obj.drawobject.fiberfiltering)
  1390. surfArray=obj.drawobject.fiberfiltering{s};
  1391. for ins=1:numel(surfArray)
  1392. try delete(surfArray(ins).toggleH); end
  1393. try delete(surfArray(ins).patchH); end
  1394. try delete(surfArray(ins)); end
  1395. end
  1396. end
  1397. obj.drawobject.fiberfiltering{s} = [];
  1398. end
  1399. end
  1400. % plot new tracts
  1401. switch lower(obj.activated.fiberfiltering)
  1402. case 'on'
  1403. obj.drawTool='fiberfiltering';
  1404. ea_unified_draw(obj);
  1405. end
  1406. % delete prior nets
  1407. if isfield(obj.drawobject,'networkmapping')
  1408. if ~isempty(obj.drawobject.networkmapping)
  1409. for s=1:numel(obj.drawobject.networkmapping)
  1410. for ins=1:numel(obj.drawobject.networkmapping{s})
  1411. try delete(obj.drawobject.networkmapping{s}{ins}.toggleH); end
  1412. try delete(obj.drawobject.networkmapping{s}{ins}.patchH); end
  1413. try delete(obj.drawobject.networkmapping{s}{ins}); end
  1414. end
  1415. end
  1416. end
  1417. end
  1418. % plot new nets
  1419. switch lower(obj.activated.networkmapping)
  1420. case 'on'
  1421. obj.drawTool='networkmapping';
  1422. ea_unified_draw(obj);
  1423. end
  1424. end
  1425. end
  1426. methods (Static)
  1427. function changeevent(~,event)
  1428. update_trajectory(event.AffectedObject,event.Source.Name);
  1429. end
  1430. end
  1431. end
  1432. function activatebychange(~,event)
  1433. % activate_tractset();
  1434. end
  1435. function calculateIntersection(obj)
  1436. for nroi = 1:length(obj.roiintersectdata)
  1437. vat = ea_load_nii(obj.roiintersectdata{nroi}); %use only one, otherwise drawing doesn't make sense
  1438. thresh = obj.roithresh;
  1439. vatInd = find(abs(vat.img(:))>thresh);
  1440. [xvox, yvox, zvox] = ind2sub(size(vat.img), vatInd);
  1441. vatmm = ea_vox2mm([xvox, yvox, zvox], vat.mat);
  1442. for side = 1:2
  1443. for i=1:size(obj.drawobject,1)
  1444. vals = {};
  1445. valsPeak = {};
  1446. connected = [];
  1447. trimmedFiberInd = [];
  1448. resultFibers = obj.fiberdrawn.fibcell{i,side};
  1449. if isempty(resultFibers)
  1450. continue
  1451. end
  1452. fibers=ea_fibcell2fibmat(resultFibers);
  1453. filter = all(fibers(:,1:3)>=min(vatmm),2) & all(fibers(:,1:3)<=max(vatmm), 2);
  1454. if ~any(filter)
  1455. zeros_arr = zeros(size(obj.drawobject{i,side},1),1);
  1456. normwts = mat2cell(zeros_arr,ones(size(obj.drawobject{i,side},1),1));
  1457. [obj.drawobject{i,side}.FaceAlpha]=normwts{:};
  1458. continue
  1459. end
  1460. trimmedFiber = fibers(filter,:);
  1461. % Map mm connectome fibers into VAT voxel space
  1462. [trimmedFiberInd, ~, trimmedFiberID] = unique(trimmedFiber(:,4), 'stable');
  1463. fibVoxInd = splitapply(@(fib) {ea_mm2uniqueVoxInd(fib, vat)}, trimmedFiber(:,1:3), trimmedFiberID);
  1464. % Remove outliers
  1465. fibVoxInd(cellfun(@(x) any(isnan(x)), fibVoxInd)) = [];
  1466. trimmedFiberInd(cellfun(@(x) any(isnan(x)), fibVoxInd)) = [];
  1467. connected = cellfun(@(fib) any(ismember(fib, vatInd)), fibVoxInd);
  1468. vals = cellfun(@(fib) vat.img(intersect(fib, vatInd)), fibVoxInd(connected), 'Uni', 0);
  1469. valsPeak{1}(trimmedFiberInd(connected)) = cellfun(@mean, vals);
  1470. wts = cell2mat(valsPeak)';
  1471. if ~isempty(wts)
  1472. if length(wts) ~= size(obj.drawobject{i,side},1)
  1473. diff = length(wts) - size(obj.drawobject{i,side},1);
  1474. if diff < 0
  1475. wts = [wts;zeros(abs(diff),1)];
  1476. end
  1477. end
  1478. normwts = normalize(ea_contrast(wts,10,0),'range');
  1479. normwts = mat2cell(normwts,ones(size(normwts,1),1));
  1480. if ~isempty(normwts) && ~isempty(obj.drawobject{i,side})
  1481. try
  1482. [obj.drawobject{i,side}.FaceAlpha]=normwts{:};
  1483. disp(['Changed alpha of tract',num2str(i)]);
  1484. normwts = {};
  1485. end
  1486. end
  1487. else %if it is not connected then they should have zero alpha!!
  1488. zeros_arr = zeros(size(obj.drawobject{i,side},1),1);
  1489. normwts = mat2cell(zeros_arr,ones(size(obj.drawobject{i,side},1),1));
  1490. [obj.drawobject{i,side}.FaceAlpha]=normwts{:};
  1491. end
  1492. end
  1493. end
  1494. end
  1495. end
  1496. function check_and_update_visibility(obj, field, ~, condition, group)
  1497. if nargin < 5 % If group is not provided, operate on obj directly
  1498. group = [];
  1499. end
  1500. if eval(sprintf('obj.%s%s && all(values %s)', field, group_access(group), condition))
  1501. eval(sprintf('obj.%s%s = 0;', field, group_access(group)));
  1502. fprintf('\n')
  1503. warning('off', 'backtrace');
  1504. warning('No %s values found, %s is set to 0 now!', condition_description(condition), field);
  1505. warning('on', 'backtrace');
  1506. fprintf('\n')
  1507. end
  1508. end
  1509. function access = group_access(group)
  1510. if isempty(group)
  1511. access = '';
  1512. else
  1513. access = sprintf('(group)');
  1514. end
  1515. end
  1516. function desc = condition_description(condition)
  1517. if strcmp(condition, '<0')
  1518. desc = 'positive';
  1519. else
  1520. desc = 'negative';
  1521. end
  1522. end
  1523. function fibers=ea_fibcell2fibmat(fibers)
  1524. [idx,~]=cellfun(@size,fibers);
  1525. fibers=cell2mat(fibers);
  1526. idxv=zeros(size(fibers,1),1);
  1527. lid=1; cnt=1;
  1528. for id=idx'
  1529. idxv(lid:lid+id-1)=cnt;
  1530. lid=lid+id;
  1531. cnt=cnt+1;
  1532. end
  1533. fibers=[fibers,idxv];
  1534. end

ea_unifiedmapping.m at commit 5b1008d, under GPL-3.0 · at the source

Overview

Authors: Patricia Zvarova1,2,3,4, Christina van der Linden5, Ningfei Li1,4, Konstantin Butenko1, Thea Berger5, Garance M Meyer3, Ilkem Aysu Sahin1,2,4, Lukas L Goede1,3, Bahne H Bahners3,6,7, Barbara Hollunder1,2,8, Till A Dembek5, Andrew R Pines3,9, Martin Reich10, Jens Volkmann10, Vincent J J Odekerken11, Rob M A de Bie11, Xin Xu12, Zhipei Ling13, Chen Yao14, Andrea A Kühn1,2,8,15,16
and 8 other authorsSurjo R Soekadar2,17, Kerstin Ritter2,18,19, Michael T Barbe5, Veerle Visser‐Vandewalle20, Michael D Fox3, Jan Niklas Petry‐Schmelzer5, Nanditha Rajamani1,3, Andreas Horn3,4,21
21 affiliations
  1. Movement Disorder and Neuromodulation Unit, Department of Neurology, Charité‐Universitätsmedizin Berlin, Corporate Member of Freie Universität Berlin and Humboldt‐Universität zu Berlin, Berlin, Germany
  2. Einstein Center for Neurosciences Berlin, Charité – Universitätsmedizin Berlin, Berlin, Germany
  3. Center for Brain Circuit Therapeutics, Department of Neurology, Brigham and Women's Hospital, Harvard Medical School, Boston, MA, USA
  4. Network Stimulation Institute, Department of Stereotactic and Functional Neurosurgery, University Hospital Cologne, Cologne, Germany
  5. Department of Neurology, Faculty of Medicine and University Hospital Cologne, University of Cologne, Cologne, Germany
  6. Department of Neurology, Center for Movement Disorders and Neuromodulation, Medical Faculty and University Hospital Düsseldorf, Heinrich Heine University Düsseldorf, Düsseldorf, Germany
  7. Institute of Clinical Neuroscience and Medical Psychology, Medical Faculty and University Hospital Düsseldorf, Heinrich Heine University Düsseldorf, Düsseldorf, Germany
  8. Berlin School of Mind and Brain, Humboldt‐Universität zu Berlin, Berlin, Germany
  9. Department of Psychiatry, Brigham & Women's Hospital, Harvard Medical School, Boston, MA, USA
  10. Department of Neurology, University Clinic of Würzburg, Würzburg, Germany
  11. Department of Neurology, Amsterdam University Medical Center, University of Amsterdam, Amsterdam, The Netherlands
  12. Department of Neurosurgery, Chinese PLA General Hospital, Beijing, China
  13. Department of Neurosurgery, Hainan Hospital of Chinese PLA General Hospital, Sanya, China
  14. Department of Neurosurgery, The National Key Clinic Specialty, Shenzhen Key Laboratory of Neurosurgery, the First Affiliated Hospital of Shenzhen University, Shenzhen Second People's Hospital, Shenzhen, China
  15. Bernstein Center for Computational Neuroscience Berlin, Berlin, Germany
  16. NeuroCure Clinical Research Centre, Charité – Universitätsmedizin Berlin, corporate member of Freie Universität Berlin and Humboldt‐Universität zu Berlin, Berlin, Germany
  17. Clinical Neurotechnology Laboratory, Department of Psychiatry and Neurosciences (CCM), Charité ‐ Universitätsmedizin Berlin, Berlin, Germany
  18. Berlin Center for Advanced Neuroimaging (BCAN), Charité — Universitätsmedizin Berlin, Berlin, Germany
  19. Hertie Institute for AI in Brain Health, University of Tübingen, Tübingen, Germany
  20. Department of Stereotactic and Functional Neurosurgery, Faculty of Medicine and University Hospital Cologne, University of Cologne, Cologne, Germany
  21. MGH Neurosurgery and Center for Neurotechnology and Neurorecovery (CNTR) at MGH Neurology Massachusetts General Hospital, Harvard Medical School, Boston, MA, USA
Journal: Annals of neurology, volume 100, issue 1, pages 22-35
Dates: received 1 September 2025; accepted 25 February 2026; published online 17 April 2026; in print July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1002/ana.78206 · PMID 41992934 · PMCID PMC13327559 · OpenAlex W4414141573
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: human (organism), Parkinson's (population), clinical / translational (subfield)
Methods: Machine learning, Statistics, Connectivity
MeSH: Deep Brain Stimulation*, Multimodal Imaging*, Parkinson Disease*, Subthalamic Nucleus*, Aged, Female, Humans, Magnetic Resonance Imaging, Male, Middle Aged, Neuroimaging, Retrospective Studies, Treatment Outcome (* major topic)
Topic: Neurological disorders and treatments (Neurology, Medicine), according to OpenAlex
Funding: NIH HHS (1R01NS127892-01, R01MH130666, UM1NS132358, 2R01 MH113929); National Institutes of Health (R01MH130666, UM1NS132358, 1R01NS127892‐01, 2R01 MH113929); Koeln Fortune Scholarship; Einstein Stiftung Berlin (424778381); Deutsche Forschungsgemeinschaft (502436811); Einstein Center for Neurosciences Berlin
Citations: not cited yet (Europe PMC); 43 references in the paper

Abstract

Objective: Accurate electrode placement and individual stimulation parameters influence the outcomes of subthalamic deep brain stimulation in Parkinson's disease. Neuroimaging‐based models can help evaluate how electrode placement impacts improvement, aiming to reduce the burden of programming. However, most existing models have been developed to explain differences between patients rather than differences between contacts within the same patient, leaving the clinical relevance of image‐guided programming unclear.

Methods: We analyzed data from patients with Parkinson's disease treated with subthalamic deep brain stimulation to develop and validate a neuroimaging‐informed model of motor improvement measured by the Unified Parkinson's Disease Scale. Five approaches were tested: active contact coordinates, electric fields, tract activations, as well as structural and functional networks. All approaches were integrated into a combined ridge regression model and validated using 2 hold‐out datasets.

Results: The sample included 236 patients (604 stimulation sites), divided into a training cohort (N = 129), a retrospective validation cohort (N = 89), and a prospectively acquired validation cohort (N = 21 electrodes). Consistent with expectations, our model explained approximately 12% of the variance in unseen group‐level data (R 2 = 0.12, p = 0.001). At the individual level, the model identified the optimal clinical contact or its neighboring contact in all but one case (mixed‐effects R 2 = 0.31, p = 3.67 × 10−10).

Interpretation: An imaging‐informed model explained the expected variance at the group level and demonstrated potential for guiding stimulation programming, suggesting that image‐guided approaches may improve clinical decision making while reducing the need for lengthy postoperative testing. ANN NEUROL 2026;100:22–35

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 7 matches between paragraphs and lines of code.

leaddbs/leaddbs

License: GPL-3.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 5b1008d705e97fe0c8f693dece14455607f4afac, 2 March 2026
Languages: MATLAB (2362), Python (101), C++ (30), C (29), Shell (10), C/C++ (4), Java (3)
Size: 6,260 files, 2,539 scripts
Software Heritage: archived
Found in: “Data availability”
Holds: README, license file, CITATION.cff, environment (ext_libs/PaCER/docs/requirements.txt), documentation
Not found: tests, continuous integration
Tools: SPM (121 files), Statistics and Machine Learning Toolbox (84 files), Tools for NIfTI and ANALYZE image (MATLAB) (38 files), Image Processing Toolbox (36 files), NumPy (35 files), FieldTrip (24 files), h5py (12 files), SciPy (8 files), FreeSurfer (7 files), cifti-matlab (6 files), Signal Processing Toolbox (6 files), Matplotlib (6 files), Parallel Computing Toolbox (5 files), pandas (5 files), TensorFlow (4 files), Keras (3 files), Optimization Toolbox (3 files), ANTs (2 files), export_fig (2 files), GIfTI library for MATLAB (2 files), NiBabel (2 files), Psychtoolbox (2 files), CAT12 (1 file), Curve Fitting Toolbox (1 file), pydicom (1 file), seaborn (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
2,000 files

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;
  • 1,998 scripts, each with its path and the digest of its content;
  • 7 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

All code used to analyze the dataset is openly available within Lead‐DBS (including sweet spot mapping, fiber filtering, and network mapping software: https://github.com/leaddbs/leaddbs). Cohort‐wise demographic and clinical outcomes from the training and the test cohorts are available in Supplementary Tables S2 and S3. We cannot openly share patient imaging data due to data sharing and privacy regulations, but they can be made available upon request to the corresponding primary investigator. The corresponding author and the principal investigator (P.Z. and A.H.) commit to returning data requests within a time frame of 30 days.

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

Recorded: type, language, journal, volume, issue, pages, dates, 28 authors, 13 MeSH terms, 6 funders, 43 references.

Cite

This paper

Zvarova, P., van der Linden, C., Li, N., Butenko, K., Berger, T., Meyer, G. M., Sahin, I. A., Goede, L. L., Bahners, B. H., Hollunder, B., Dembek, T. A., Pines, A. R., Reich, M., Volkmann, J., Odekerken, V. J. J., de Bie, R. M. A., Xu, X., Ling, Z., Yao, C., . . . Horn, A. (2026). Multimodal Image Guidance in Subthalamic Deep Brain Stimulation for Parkinson's Disease. Annals of neurology, 100(1), 22-35. https://doi.org/10.1002/ana.78206

BibTeX

@article{zvarova2026multimodal,
author = {Zvarova, Patricia and van der Linden, Christina and Li, Ningfei and Butenko, Konstantin and Berger, Thea and Meyer, Garance M and Sahin, Ilkem Aysu and Goede, Lukas L and Bahners, Bahne H and Hollunder, Barbara and Dembek, Till A and Pines, Andrew R and Reich, Martin and Volkmann, Jens and Odekerken, Vincent J J and de Bie, Rob M A and Xu, Xin and Ling, Zhipei and Yao, Chen and Kühn, Andrea A and Soekadar, Surjo R and Ritter, Kerstin and Barbe, Michael T and Visser‐Vandewalle, Veerle and Fox, Michael D and Petry‐Schmelzer, Jan Niklas and Rajamani, Nanditha and Horn, Andreas},
title = {{Multimodal Image Guidance in Subthalamic Deep Brain Stimulation for Parkinson's Disease}},
journal = {Annals of neurology},
year = {2026},
month = apr,
volume = {100},
number = {1},
pages = {22--35},
publisher = {Wiley},
issn = {0364-5134},
doi = {10.1002/ana.78206},
url = {https://doi.org/10.1002/ana.78206},
pmid = {41992934},
pmcid = {PMC13327559}
}

RIS

TY - JOUR
AU - Zvarova, Patricia
AU - van der Linden, Christina
AU - Li, Ningfei
AU - Butenko, Konstantin
AU - Berger, Thea
AU - Meyer, Garance M
AU - Sahin, Ilkem Aysu
AU - Goede, Lukas L
AU - Bahners, Bahne H
AU - Hollunder, Barbara
AU - Dembek, Till A
AU - Pines, Andrew R
AU - Reich, Martin
AU - Volkmann, Jens
AU - Odekerken, Vincent J J
AU - de Bie, Rob M A
AU - Xu, Xin
AU - Ling, Zhipei
AU - Yao, Chen
AU - Kühn, Andrea A
AU - Soekadar, Surjo R
AU - Ritter, Kerstin
AU - Barbe, Michael T
AU - Visser‐Vandewalle, Veerle
AU - Fox, Michael D
AU - Petry‐Schmelzer, Jan Niklas
AU - Rajamani, Nanditha
AU - Horn, Andreas
TI - Multimodal Image Guidance in Subthalamic Deep Brain Stimulation for Parkinson's Disease
T2 - Annals of neurology
J2 - Ann Neurol
PY - 2026
DA - 2026/04/17
VL - 100
IS - 1
SP - 22
EP - 35
SN - 0364-5134
PB - Wiley
DO - 10.1002/ana.78206
UR - https://doi.org/10.1002/ana.78206
LA - en
ER -

CSL-JSON

{
"id": "10.1002/ana.78206",
"type": "article-journal",
"title": "Multimodal Image Guidance in Subthalamic Deep Brain Stimulation for Parkinson's Disease",
"container-title": "Annals of neurology",
"author": [
{
"family": "Zvarova",
"given": "Patricia"
},
{
"family": "van der Linden",
"given": "Christina"
},
{
"family": "Li",
"given": "Ningfei"
},
{
"family": "Butenko",
"given": "Konstantin"
},
{
"family": "Berger",
"given": "Thea"
},
{
"family": "Meyer",
"given": "Garance M"
},
{
"family": "Sahin",
"given": "Ilkem Aysu"
},
{
"family": "Goede",
"given": "Lukas L"
},
{
"family": "Bahners",
"given": "Bahne H"
},
{
"family": "Hollunder",
"given": "Barbara"
},
{
"family": "Dembek",
"given": "Till A"
},
{
"family": "Pines",
"given": "Andrew R"
},
{
"family": "Reich",
"given": "Martin"
},
{
"family": "Volkmann",
"given": "Jens"
},
{
"family": "Odekerken",
"given": "Vincent J J"
},
{
"family": "de Bie",
"given": "Rob M A"
},
{
"family": "Xu",
"given": "Xin"
},
{
"family": "Ling",
"given": "Zhipei"
},
{
"family": "Yao",
"given": "Chen"
},
{
"family": "Kühn",
"given": "Andrea A"
},
{
"family": "Soekadar",
"given": "Surjo R"
},
{
"family": "Ritter",
"given": "Kerstin"
},
{
"family": "Barbe",
"given": "Michael T"
},
{
"family": "Visser‐Vandewalle",
"given": "Veerle"
},
{
"family": "Fox",
"given": "Michael D"
},
{
"family": "Petry‐Schmelzer",
"given": "Jan Niklas"
},
{
"family": "Rajamani",
"given": "Nanditha"
},
{
"family": "Horn",
"given": "Andreas"
}
],
"container-title-short": "Ann Neurol",
"volume": "100",
"issue": "1",
"page": "22-35",
"DOI": "10.1002/ana.78206",
"PMID": "41992934",
"PMCID": "PMC13327559",
"ISSN": "0364-5134",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/ana.78206",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
17
]
]
}
}

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.1016/j.celrep.2026.117404 [code]
Action and rest tremor map to distinct networks within the primary motor cortex.
Journal: Cell reports
In common: CAT12, cifti-matlab, export_fig, 23 other tools, 7 references, 3 authors
[2] doi:10.1038/s41586-026-10631-3 [code]
A prognostic human brain network for diffuse midline glioma.
Journal: Nature
In common: CAT12, cifti-matlab, export_fig, 23 other tools, clinical / translational, 4 references, author Andreas Horn
[3] doi:10.1111/ene.70678 [code]
Who Falls After a Stroke? Evidence From a Prospective Stroke Cohort.
Journal: European journal of neurology
In common: CAT12, cifti-matlab, export_fig, 23 other tools, clinical / translational, 2 references, author Andrea A Kühn
[4] doi:10.1002/hbm.70602 [code]
Neuroimaging Correlates of Post-Stroke Pain After Ischemic Stroke: Secondary Analysis of the INSPiRE-TMS Trial.
Journal: Human brain mapping
In common: CAT12, cifti-matlab, export_fig, 23 other tools, clinical / translational, 1 reference
[5] doi:10.1038/s41467-026-71151-2 [code]
Common and distinct neural correlates of social interaction processing and theory of mind in narratives.
Journal: Nature communications
In common: cifti-matlab, Psychtoolbox, GIfTI library for MATLAB, 15 other tools
[6] doi:10.1162/imag.a.1262 [code]
Frame-wise multi-echo distortion correction for superior functional MRI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: pydicom, GIfTI library for MATLAB, Tools for NIfTI and ANALYZE image (MATLAB), 13 other tools, 1 reference
[7] doi:10.1162/imag.a.1222 [code]
Network-based near-scalp personalized brain stimulation targets.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: GIfTI library for MATLAB, Tools for NIfTI and ANALYZE image (MATLAB), ANTs, 9 other tools, 1 reference, author Michael D Fox
[8] doi:10.1038/s42003-026-10276-y [code]
The cellular correlates and adolescent reorganisation of cortical myelination networks in the common marmoset.
Journal: Communications biology
In common: CAT12, GIfTI library for MATLAB, Optimization Toolbox, 12 other tools, 1 reference
[9] doi:10.1080/07853890.2026.2685416 [code]
Pulmonary and cerebral damage in COVID-19 survivors: is there any association?
Journal: Annals of medicine
In common: pydicom, Curve Fitting Toolbox, Tools for NIfTI and ANALYZE image (MATLAB), 12 other tools, clinical / translational
[10] doi:10.1038/s41593-026-02345-6 [code]
Human hippocampal ripples tune cortical responses based on predicted uncertainty.
Journal: Nature neuroscience
In common: Curve Fitting Toolbox, FieldTrip, Image Processing Toolbox, 2 other tools, 2 references, 2 authors

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.