Anisotropic conductivity modeling for tDCS in Parkinson’s disease using multidimensional diffusion MRI

Code ↔ Paper

21 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 21 matches · 6 of them tie a paragraph to a whole file without given lines: weak matches, whose lines are not tinted

[1] § Methods › Registration and tensor reorientation ↔ analysis/qc_harness.py, lines 167 to 234 · score 0.90 · superior inferior, cerebral peduncles, edge alignment, corpus callosum, boundary, frame · Report a mismatch
[2] § Methods › Association with mechanical properties ↔ analysis/09_mre_cohort_corr.py, lines 133 to 274 · score 0.90 · CSF morphology, partial Pearson, CSF fraction, CSF volume, added age, coefficients · Report a mismatch
[3] § Methods › Registration and tensor reorientation ↔ pipeline/02_register_dmri_to_T1.sh, lines 37 to 119 · score 0.87 · rigid pre, DOF affine, charm T2, mesh elements, conductivity tensors, FSL · Report a mismatch
[4] § Methods › E-field metrics and statistics ↔ analysis/aggregate_matched_control.py, lines 117 to 157 · score 0.85 · subject clustered bootstrap, absolute percentage, V1 angle, orientation divergence, replacement, interval · Report a mismatch
[5] § Methods › E-field metrics and statistics ↔ analysis/06_cohort_stats.py, lines 61 to 200 · score 0.83 · Mann Whitney, rank biserial, linear age, Wilcoxon, residuals, Cohen · Report a mismatch
[6] § Methods › Regions of interest ↔ analysis/build_rois.py, lines 1 to 25 · score 0.82 · Desikan Killiany, FreeSurfer, WM lobes, parietal, wmparc, frontal · Report a mismatch
[7] § Methods › Registration and tensor reorientation ↔ pipeline/08b_matched_reg_control.sh, lines 50 to 111 · score 0.80 · rigid pre, mutual information, DOF affine, FSL, FLIRT, alignment · Report a mismatch
[8] § Methods › Conductivity models ↔ config/config.example.sh, the whole file · a weak match · score 0.80 · spherical encoding series, brain mask, covariance fit, single shell, oriented, MD dMRI · Report a mismatch
[9] § Methods › Conductivity models ↔ pipeline/prepare_dmri_tensor.py, lines 65 to 88 · score 0.76 · brain voxels, positive definite, cumulant, plausible, physiological, ms · Report a mismatch
[10] § Results › Association between predicted E-field and MRE stiffness ↔ analysis/09_mre_cohort_corr.py, lines 133 to 274 · score 0.74 · partial Pearson, CSF fraction, CSF volume, attenuated, collinear, regression · Report a mismatch
[11] § Methods › Conductivity models ↔ pipeline/run_qti_cov_cohort.m, the whole file · a weak match · score 0.71 · covariance fit, motion, elastix, eddy, topup, Synb0 · Report a mismatch
[12] § Results › Participants ↔ analysis/qc_harness.py, lines 167 to 234 · score 0.71 · superior inferior, cerebral peduncle, corpus callosum, reoriented, segmentation, dMRI · Report a mismatch
[13] § Methods › Conductivity models ↔ pipeline/01_dwi2cond.sh, the whole file · a weak match · score 0.71 · T1 coregistration, dwi2cond, single shell, DTI model, FSL, dtifit · Report a mismatch
[14] § Results › Orientation divergence between the two anisotropic models ↔ analysis/aggregate_matched_control.py, lines 1 to 25 · score 0.69 · subject clustered bootstrap, absolute percentage, V1 angle, CI, orientation, dMRI · Report a mismatch
[15] § Methods › Electrode montages ↔ pipeline/04_run_simulations.py, lines 1 to 35 · score 0.63 · tDCS, cathode, anode, pads, simulated, DLPFC · Report a mismatch
[16] § Methods › Head model ↔ analysis/qc_harness.py, lines 101 to 140 · score 0.63 · SimNIBS, compact, spongy, bone, tetrahedral, scalp · Report a mismatch
[17] § Methods › Conductivity models ↔ config/config.example.sh, the whole file · a weak match · score 0.61 · covariance fit, single shell, eddy, rotated, topup, FSL · Report a mismatch
[18] § Methods › Registration and tensor reorientation ↔ pipeline/run_qti_cov_cohort.m, the whole file · a weak match · score 0.59 · covariance fit, MD dMRI tensors, spherical, smoothing, encoded, linear · Report a mismatch
[19] § Results › Group comparison (PD vs. HC) ↔ analysis/06_cohort_stats.py, lines 61 to 200 · score 0.58 · Mann Whitney, age adjusted, residuals, Cohen, HC, dMRI · Report a mismatch
[20] § Methods › Head model ↔ pipeline/00_charm.sh, the whole file · a weak match · score 0.56 · head model, SimNIBS, tetrahedral, pipeline, charm, mesh · Report a mismatch
[21] § Methods › Study design ↔ pipeline/04_run_simulations.py, lines 1 to 35 · score 0.54 · tDCS, SimNIBS, cathodes, anode, pads, DLPFC · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

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

The authors' code

Python · 587 lines · 29 KB · MIT · 3 matches

"""
qc_harness.py - post-pipeline per-subject QC: one qc_summary.csv row per subject with sanity
metrics, a PASS/FLAG per stage, and an overall verdict. Beyond the absolute thresholds, any
metric > 3 MAD from the cohort median is also flagged. An outlier is a signal to inspect/fix.
Usage:
simnibs_python analysis/qc_harness.py # the subject in config/config.sh
simnibs_python analysis/qc_harness.py --cohort config/cohort.json # {"subjects":[{"id",...}]} or [{"id","work","m2m"}...]
simnibs_python analysis/qc_harness.py --calibrate # suggest cohort-percentile thresholds, then exit
"""
import os, sys, csv, json, glob, argparse
import numpy as np
import nibabel as nib
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "pipeline"))
from _config import cfg # noqa: E402
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from _sims import sim_mesh, MODELS # noqa: E402 (shared montage-aware mesh lookup)
from _rois import eigvals_6comp, fa_from_evals # noqa: E402 (shared eigendecomp + FA; one source of truth)
# Absolute thresholds seeded from a single-subject pilot (n=1). [PROV] items are provisional
# single-subject cutoffs to recalibrate as a cohort percentile (--calibrate) once 5-8 subjects
# exist. The MAD>3 outlier pass in main() (active only at n>=4) is the real per-subject guard;
# on one subject only these absolute thresholds are in force.
THR = dict(
charm_tet_q_med_min=0.10, charm_tet_q_badfrac_max=0.05, # broadly-poor mesh / too many slivers
tissue_plausible={"WM": (3e5, 7.5e5), "GM": (4e5, 8e5), "CSF": (5e4, 6e5),
"bone": (2e5, 8.5e5), "scalp": (5e5, 2.5e6)}, # mm3; upper bounds recalibrated to the n=29 cohort
cc_fa_min=0.25, csf_fa_max=0.25, # dwi2cond; cc_fa_min recalibrated to the n=29 cohort min
reg_containment_min=0.95, # dMRI FOV covers brain (gross misreg / cutoff)
cc_v1x_min=0.6, # [PROV] CC V1 left-right (tensor reorientation)
peduncle_siz_min=0.55, # [PROV] cerebral-peduncle V1 superior-inferior in the SimNIBS frame
reg_grad_corr_min=0.08, # [PROV] b0<->T1 edge alignment (mid-range shift)
efield_p95_range=(0.10, 0.60), efield_spike=8.0, # 03: p95 |E| range; max/p95 spike [PROV] default mesh ~6 (tail artifact, not the median readout); recalibrate at cohort scale
electrode_tol_mm=15.0,
)
SIGMA0 = {"WM": 0.126, "GM": 0.275} # tissue baseline conductivities (S/m), vn anchors
SAMSEG = dict(wm=(2, 41), gm=(3, 42), csf=(4, 43, 14, 15, 24), brainstem=(16,))
NUCLEI = ["Pu", "Ca", "NAC", "GPe", "GPi", "SNc", "SNr", "VTA", "RN", "STN"]
def _load(p):
return nib.load(p) if os.path.exists(p) else None
def _data(p):
img = _load(p)
return np.asarray(img.dataobj) if img is not None else None
def _roi_dir(reg):
"""Group 1 ROI dir: the recon-all parcellation (freesurfer_rois) built by analysis/build_rois.py.
Mirrors analysis/_rois.load_labeled so QC reads the same masks."""
return os.path.join(reg, "freesurfer_rois")
def _eigs_masked(t, sel):
"""Ascending eigenvalues (N,3) of the 6-comp tensor `t` over the selected voxels only. Delegates to
_rois.eigvals_6comp (one source of truth for the FSL pack + the eigvalsh values-only fast path)."""
return eigvals_6comp(t, sel)
def _tensor_eigs(path, mask=None):
"""Return sorted eigenvalues (N,3) of a 6-comp tensor NIfTI, optionally within mask."""
t = _data(path)
if t is None:
return None
sel = (np.abs(t[..., 0]) > 1e-9) if mask is None else mask
return _eigs_masked(t, sel)
def _fa(eigs):
"""Mean FA over voxels, from ascending eigenvalues (N,3); uses _rois.fa_from_evals (per-voxel)."""
return np.nanmean(fa_from_evals(eigs))
def _vn_check(e, sigma0):
"""Volume-normalization sanity checks on diffusion-tensor eigenvalues `e` (N,3, ascending).
SimNIBS's vn step scales each voxel's D to sigma = D * sigma0 / geomean(eig D). Reconstructing
sigma from the SAME eigenvalues makes geomean(eig sigma) == sigma0 hold by construction, so that
identity is NOT a check (it can only be violated by reading SimNIBS's own assigned conductivity,
which this function does not have). What this DOES verify:
Returns (degfrac, sigma_hi, sigma_lo):
degfrac - fraction NOT positive-definite (vn undefined there). MUST be read from the raw
smallest eigenvalue e[:,0]: a clamped geom-mean is always > 0 and would silently
report 0 (the bug this function exists to prevent).
sigma_hi/lo - 99th/1st percentile of the reconstructed conductivity eigenvalues, i.e. whether
the tensor shape scaled to sigma0 lands in a physiological conductivity range.
"""
good = np.isfinite(e).all(axis=1) & (e[:, 0] > 0)
degfrac = float(np.mean(~good))
if not good.any():
return degfrac, float("nan"), float("nan")
eg = e[good] # positive-definite -> no clamp needed
g = np.exp(np.mean(np.log(eg), axis=1))
sig = eg * (sigma0 / g[:, None]) # reconstructed conductivity eigenvalues
return degfrac, float(np.percentile(sig, 99)), float(np.percentile(sig, 1))
# stage checks: each returns (metrics dict, list of flag strings)
def qc_charm(P):
m, f = {}, []
m2m = P["m2m"]
mesh = next(iter(glob.glob(os.path.join(m2m, "*.msh"))), None)
ft = os.path.join(m2m, "final_tissues.nii.gz")
m["charm_m2m"] = int(os.path.isdir(m2m)); m["charm_mesh"] = int(mesh is not None)
m["charm_tissues"] = int(os.path.exists(ft))
if not (m["charm_m2m"] and mesh and m["charm_tissues"]):
f.append("charm:missing_output"); return m, f
try:
from simnibs import mesh_io
msh = mesh_io.read_msh(mesh)
m["charm_n_tet"] = int((msh.elm.elm_type == 4).sum())
try:
# tetrahedra_quality() -> (radius_edge_ratio, inscribed/circumscribed ratio); [1] is the
# standard tet quality (1.0 = ideal, low = sliver), sized to ALL elements so triangles are
# NaN -> keep finite (the tets). Gate on median + sliver fraction, not min (a few slivers
# are normal and the FEM handles them).
q = msh.tetrahedra_quality()[1].value
qt = q[np.isfinite(q)]
m["charm_tet_q_med"] = round(float(np.median(qt)), 3)
m["charm_tet_q_badfrac"] = round(float(np.mean(qt < 0.1)), 4)
if np.median(qt) < THR["charm_tet_q_med_min"] or np.mean(qt < 0.1) > THR["charm_tet_q_badfrac_max"]:
f.append("charm:low_tet_quality")
except Exception:
m["charm_tet_q_med"] = np.nan; m["charm_tet_q_badfrac"] = np.nan
except Exception as e:
m["charm_n_tet"] = np.nan; m["charm_tet_q_med"] = np.nan; m["charm_tet_q_badfrac"] = np.nan
f.append(f"charm:mesh_read({type(e).__name__})")
seg = _data(ft); vx = abs(np.linalg.det(_load(ft).affine[:3, :3]))
# SimNIBS charm final_tissues labels: WM=1, GM=2, CSF=3, scalp=5; bone is SPLIT into 7=compact +
# 8=spongy (generic bone label 4 is EMPTY in charm output), so total bone = 7+8.
lut = {"WM": [1], "GM": [2], "CSF": [3], "bone": [7, 8], "scalp": [5]}
for name, labs in lut.items():
v = float(np.isin(seg, labs).sum() * vx); m[f"vol_{name}_mm3"] = round(v)
lo, hi = THR["tissue_plausible"].get(name, (0, 1e12))
if not (lo <= v <= hi):
f.append(f"charm:vol_{name}_implausible")
return m, f
def qc_dwi2cond(P):
m, f = {}, []
dti = os.path.join(P["m2m"], "DTI_coregT1_tensor.nii.gz")
if not os.path.exists(dti):
f.append("dwi2cond:missing"); return m, f
seg = _data(P["samseg"])
cc = _data(P.get("cc_mask")) if P.get("cc_mask") else None
cc_mask = (cc > 0) if cc is not None else None
csf_mask = np.isin(seg, SAMSEG["csf"]) if seg is not None else None
if cc_mask is not None:
ec = _tensor_eigs(dti, cc_mask); m["dwi2cond_cc_fa"] = round(_fa(ec), 3)
if m["dwi2cond_cc_fa"] < THR["cc_fa_min"]:
f.append("dwi2cond:cc_fa_low")
if csf_mask is not None:
es = _tensor_eigs(dti, csf_mask); m["dwi2cond_csf_fa"] = round(_fa(es), 3)
if m["dwi2cond_csf_fa"] > THR["csf_fa_max"]:
f.append("dwi2cond:csf_fa_high")
eall = _tensor_eigs(dti)
m["dwi2cond_neg_eig_frac"] = round(float(np.mean(eall[:, 0] <= 0)), 4) if eall is not None else np.nan
if eall is not None and (eall[:, 0] <= 0).mean() > 0.001:
f.append("dwi2cond:neg_eigenvalues")
return m, f
def qc_register(P):
m, f = {}, []
seg = _data(P["samseg"])
dm = _data(os.path.join(P["reg"], "dMRI_mask_T1.nii.gz"))
if seg is None or dm is None:
f.append("register:missing"); return m, f
# The registered dMRI brain mask is ~2x the SAMSEG brain (includes neck/inferior FOV), so Dice or
# centroid distance is size-confounded (caps ~0.62 / ~8 mm even when aligned). Use size-INDEPENDENT
# signals instead:
# (1) containment - does the dMRI FOV cover the brain? Catches gross misreg / FOV cutoff.
# (2) CC reorientation - do corpus-callosum V1 point left-right? Catches tensor-reorient errors
# (highest-risk step: a flip/rotation collapses |V1_x| in the CC).
brain = np.isin(seg, [2, 41, 3, 42, 10, 49, 11, 50, 12, 51, 13, 52, 16, 17, 18, 53, 54, 26, 58, 28, 60])
dmm = dm > 0
m["reg_containment"] = round(float((dmm & brain).sum() / max(int(brain.sum()), 1)), 4)
if m["reg_containment"] < THR["reg_containment_min"]:
f.append("register:brain_not_covered")
v1 = _data(os.path.join(P["reg"], "v1_T1.nii.gz"))
cc = _data(P.get("cc_mask")) if P.get("cc_mask") else None
if v1 is not None and cc is not None and (cc > 0).any():
ccm = cc > 0
nrm = np.linalg.norm(v1[ccm], axis=1)
ok = nrm > 1e-6
if ok.any():
v1x = np.abs(v1[ccm][ok, 0]) / nrm[ok]
m["reg_cc_v1x"] = round(float(np.median(v1x)), 3) # median: robust to CC-edge partial volume
if np.median(v1x) < THR["cc_v1x_min"]:
f.append("register:cc_not_LR") # tensor reorientation failed
# (2b) cerebral-peduncle reorientation in the SimNIBS frame. CC V1x above is voxel-frame and the
# CC L-R signal is resolution-limited at 2.5 mm; the peduncle is a thick coherent S-I tract, the
# reliable orientation gate. SimNIBS rotates the tensor by M = colnorm(affine).diag(-1,1,1 if det>0)
# (cond_utils.cond2elmdata) before the FEM, so we score M@v1 against world-z (superior-inferior).
ped = _data(os.path.join(_roi_dir(P["reg"]), "roi_Mesencephalon.nii.gz"))
if v1 is not None and ped is not None and (ped > 0).any():
aff = nib.load(os.path.join(P["reg"], "v1_T1.nii.gz")).affine[:3, :3]
M = aff / np.linalg.norm(aff, axis=0)[:, None]
if np.linalg.det(M) > 0:
M = M.dot(np.diag([-1.0, 1.0, 1.0]))
vw = v1[ped > 0] @ M.T
nrm = np.linalg.norm(vw, axis=1); ok = nrm > 1e-6
if ok.any():
siz = np.abs(vw[ok, 2]) / nrm[ok]
m["reg_peduncle_siz"] = round(float(np.median(siz)), 3)
if np.median(siz) < THR["peduncle_siz_min"]:
f.append("register:peduncle_not_SI") # reorientation / registration off
# (3) edge alignment - gradient-magnitude correlation of registered b0 vs T1 inside the brain.
# Containment and CC-V1 miss a few-mm rigid shift that still covers the brain; a shift misaligns
# tissue boundaries, so |grad| correlation drops sharply (BBR principle): ~0.19 aligned, ~-40% at
# a 4 mm shift. (NMI on this b0/T1 pair sits near the information floor and barely moves, so unused.)
try:
from scipy import ndimage
b0 = _data(os.path.join(P["reg"], "s0_T1.nii.gz")) # registered dMRI S0 in T1 space
t1 = _data(os.path.join(P["m2m"], "T1.nii.gz"))
if b0 is not None and t1 is not None:
inner = ndimage.binary_erosion(brain, iterations=2) # avoid the mask boundary itself
def _gmag(v):
gx, gy, gz = np.gradient(ndimage.gaussian_filter(v.astype(float), 1.0))
return np.sqrt(gx * gx + gy * gy + gz * gz)
a = _gmag(b0)[inner]; b = _gmag(t1)[inner]
a = a - a.mean(); b = b - b.mean()
gc = float((a * b).sum() / (np.sqrt((a * a).sum() * (b * b).sum()) + 1e-30))
m["reg_grad_corr"] = round(gc, 3)
if gc < THR["reg_grad_corr_min"]:
f.append("register:edge_misaligned") # [PROV] likely a few-mm shift
except Exception:
pass # scipy missing or b0 absent -> skip; containment+CC still gate
return m, f
def qc_conductivity(P):
m, f = {}, []
seg = _data(P["samseg"])
# The MD-dMRI conductivity input is the plain QTI mean tensor <D> (tensor_MD_dMRI.nii.gz).
tensors = [("MDdMRI", "tensor_MD_dMRI.nii.gz")]
for tag, fname in tensors:
path = os.path.join(P["work"], fname)
if not os.path.exists(path):
f.append(f"cond:{tag}_missing"); continue
t = _data(path) # load once; reuse for global + per-tissue
naninf = bool(np.isnan(t).any() or np.isinf(t).any())
m[f"cond_{tag}_nan"] = int(naninf)
if naninf:
f.append(f"cond:{tag}_naninf")
valid = np.abs(t[..., 0]) > 1e-9
if not valid.any():
f.append(f"cond:{tag}_empty"); continue
eigs = _eigs_masked(t, valid)
negfrac = float(np.mean(eigs[:, 0] <= 0))
m[f"cond_{tag}_neg_eig_frac"] = round(negfrac, 5)
if negfrac > 0:
f.append(f"cond:{tag}_neg_eig")
ratio = eigs[:, 2] / np.maximum(eigs[:, 0], 1e-9)
m[f"cond_{tag}_cap8_frac"] = round(float(np.mean(ratio > 8)), 4)
# Per-tissue mean diffusivity = geom-mean of the D eigenvalues. tensor_MD_dMRI is the DIFFUSION
# tensor; SimNIBS applies vn (geom-mean sigma -> sigma0) at sim time, so sigma0 is verified
# INDIRECTLY via the E-field range in qc_sims. Here we sanity-check the tensor: positive, real
# MD-scale, not corrupt/zeroed.
if seg is not None:
for tis, labs in [("WM", SAMSEG["wm"]), ("GM", SAMSEG["gm"])]:
tm = np.isin(seg, labs) & valid
if not tm.any():
continue
e = _eigs_masked(t, tm)
md_med = float(np.median(np.exp(np.mean(np.log(np.maximum(e, 1e-12)), axis=1))))
m[f"cond_{tag}_{tis}_md"] = round(md_med, 3)
if not (0.2 <= md_med <= 3.0): # plausible mean diffusivity band (um^2/ms)
f.append(f"cond:{tag}_{tis}_md_implausible")
# vn check (see _vn_check): count voxels that cannot be vn-normalized (not
# positive-definite) and confirm the reconstructed sigma is physiological.
s0 = SIGMA0.get(tis)
if s0:
degfrac, sig_hi, sig_lo = _vn_check(e, s0)
m[f"cond_{tag}_{tis}_vn_degfrac"] = round(degfrac, 5)
if degfrac > 0.001:
f.append(f"cond:{tag}_{tis}_vn_degenerate") # D cannot be vn-normalized here
if np.isfinite(sig_hi):
m[f"cond_{tag}_{tis}_sigma_hi"] = round(sig_hi, 3)
if sig_hi > 2.0 or sig_lo < 1e-3:
f.append(f"cond:{tag}_{tis}_sigma_nonphysio") # runaway / insulating sigma
return m, f
def qc_sims(P):
m, f = {}, []
try:
from simnibs import mesh_io
except Exception:
f.append("sims:no_simnibs"); return m, f
montage = P.get("montage", "M1")
for model in MODELS: # ISO, DTI, MD-dMRI
mod = model.replace("-", "") # metric label: ISO / DTI / MDdMRI
mp = sim_mesh(P["work"], montage, model, P["id"])
if not mp:
f.append(f"sims:{mod}_missing"); continue
msh = mesh_io.read_msh(mp)
gm = (msh.elm.tag1 == 2)
E = msh.field["magnE"].value
p95 = float(np.percentile(E[gm], 95)); mx = float(np.max(E[gm]))
m[f"E_{mod}_p95"] = round(p95, 3); m[f"E_{mod}_spike"] = round(mx / max(p95, 1e-9), 1)
if not (THR["efield_p95_range"][0] <= p95 <= THR["efield_p95_range"][1]):
f.append(f"sims:{mod}_p95_out_of_range")
if mx / max(p95, 1e-9) > THR["efield_spike"]:
f.append(f"sims:{mod}_spike")
return m, f
def qc_rois(P):
m, f = {}, []
n_empty = 0; n_total = 0
for p in glob.glob(os.path.join(_roi_dir(P["reg"]), "roi_*.nii.gz")):
n_total += 1
if (_data(p) > 0).sum() == 0:
n_empty += 1; f.append(f"rois:empty({os.path.basename(p)})")
m["roi_n_total"] = n_total; m["roi_n_empty"] = n_empty
if n_total == 0:
f.append("rois:none_built")
return m, f
def qc_nuclei(P):
# Group 2 deep nuclei (CIT168/Pauli) are OVERLAP-ALLOWED separate binary masks (07_build_nuclei.sh).
# Glob the per-nucleus masks; do NOT read nuclei_labeled.nii.gz (winner-take-all int labels destroy
# the intended overlap). Pairwise overlap is expected (midbrain/STN), logged not flagged.
m, f = {}, []
nd = os.path.join(P["reg"], "atlas_rois", "nuclei")
paths = sorted(glob.glob(os.path.join(nd, "roi_*.nii.gz")))
paths = [p for p in paths if os.path.basename(p) != "nuclei_labeled.nii.gz"]
masks = {}
for p in paths:
d = _data(p)
if d is None:
continue
b = d > 0
masks[os.path.basename(p)] = b
if b.sum() == 0:
f.append(f"nuclei:empty({os.path.basename(p)})")
m["nuclei_n_total"] = len(masks)
m["nuclei_n_empty"] = sum(1 for b in masks.values() if b.sum() == 0)
names = sorted(masks)
n_ov = 0
for i in range(len(names)):
for j in range(i + 1, len(names)):
if (masks[names[i]] & masks[names[j]]).any():
n_ov += 1
m["nuclei_n_overlap_pairs"] = n_ov # expected (overlap-allowed), reported not flagged
if paths and not masks:
f.append("nuclei:unreadable")
# presence: each nucleus should have an L and R mask (07 emits roi_{L,R}_{nucleus}.nii.gz)
if masks:
for nuc in NUCLEI:
for side in ("L", "R"):
if f"roi_{side}_{nuc}.nii.gz" not in masks:
f.append(f"nuclei:missing(roi_{side}_{nuc})")
return m, f
def qc_roi_efield_csv(P):
# Mask presence (qc_rois/qc_nuclei) is necessary but not sufficient: an ROI whose mask exists but
# overlaps no GM/WM tetrahedra reads back as all-NaN in the per-subject roi_efield CSV and would
# otherwise pass QC silently. Open the actual CSV that 04_extract_roi_efield.py wrote and flag any
# ROI that is all-NaN across a model's columns.
m, f = {}, []
montage = P.get("montage", "M1")
csv_path = os.path.join(P["results_root"], P["id"], f"roi_efield_{montage}.csv")
if not os.path.exists(csv_path):
f.append(f"roi_efield:missing(roi_efield_{montage}.csv)"); return m, f
rows = list(csv.DictReader(open(csv_path)))
if not rows:
f.append("roi_efield:empty"); return m, f
# Per-model columns are "<model>_p95" / "<model>_median" (model names from _sims.MODELS).
cols = {mod: [f"{mod}_p95", f"{mod}_median"] for mod in MODELS}
def _is_nan(v):
if v is None or v == "":
return True
try:
return not np.isfinite(float(v))
except ValueError:
return True
n_nan_rows = 0
for r in rows:
roi = r.get("ROI", "?")
for mod, mcols in cols.items():
present = [c for c in mcols if c in r]
if present and all(_is_nan(r[c]) for c in present):
n_nan_rows += 1
f.append(f"roi_efield:all_nan({roi},{mod})") # mask exists but no GM/WM overlap
m["roi_efield_n_rows"] = len(rows); m["roi_efield_n_nan"] = n_nan_rows
return m, f
STAGES = [("00_charm", qc_charm), ("01_dwi2cond", qc_dwi2cond), ("02_register", qc_register),
("03_conductivity", qc_conductivity), ("04_sims", qc_sims), ("05_rois", qc_rois),
("06_nuclei", qc_nuclei), ("07_roi_efield", qc_roi_efield_csv)]
def _cohort_local_root():
"""Root holding the per-subject dirs the cohort.json ids index into: data/cohort_local/<id>/
(each with work/m2m_<id> + registration). Resolved relative to the repo root, NOT cfg["WORK_DIR"]
(which tracks the active single-subject config and may point at the pilot, not the cohort layout).
Override with COHORT_LOCAL_ROOT if the cohort lives elsewhere."""
env = os.environ.get("COHORT_LOCAL_ROOT")
if env:
return env
repo = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
return os.path.join(repo, "data", "cohort_local")
def _subject_paths(sid, root):
"""Per-subject path dict from the data/cohort_local/<id>/ layout (work/m2m_<id>, registration)."""
base = os.path.join(root, sid)
work = os.path.join(base, "work")
m2m = os.path.join(work, f"m2m_{sid}")
reg = os.path.join(base, "registration")
return {"id": sid, "work": work, "m2m": m2m, "reg": reg,
"samseg": os.path.join(m2m, "segmentation", "labeling.nii.gz"),
"cc_mask": os.path.join(_roi_dir(reg), "roi_CC.nii.gz")}
def resolve_subjects(args):
if args.cohort and os.path.exists(args.cohort):
data = json.load(open(args.cohort))
# Two accepted shapes:
# (a) config/cohort.json: {"montages":[...], "subjects":[{"id","group",...}]} - derive
# per-subject paths from the data/cohort_local/<id>/ layout.
# (b) a flat list of explicit path dicts ([{"id","work","m2m",...}]) - used as-is.
if isinstance(data, dict) and "subjects" in data:
root = _cohort_local_root()
return [_subject_paths(s["id"], root) for s in data["subjects"]]
return data
return [{"id": cfg["SUBJECT"], "work": cfg["WORK_DIR"], "m2m": cfg["M2M_DIR"],
"reg": cfg["REG_DIR"],
"samseg": os.path.join(cfg["M2M_DIR"], "segmentation", "labeling.nii.gz"),
"cc_mask": os.path.join(_roi_dir(cfg["REG_DIR"]), "roi_CC.nii.gz")}]
# Directional metrics -> which tail to guard once a cohort exists ("low" = flag below the lower tail,
# "high" = above the upper). Used by --calibrate to turn single-subject cutoffs into cohort percentiles.
CALIB = {"reg_grad_corr": "low", "reg_cc_v1x": "low", "reg_containment": "low",
"dwi2cond_cc_fa": "low", "charm_tet_q_med": "low",
"dwi2cond_csf_fa": "high", "charm_tet_q_badfrac": "high"}
def _calibrate(out):
p = os.path.join(out, "qc_summary.csv")
if not os.path.exists(p):
print(f"no cohort CSV at {p} - run the harness over the cohort first"); return
rows = list(csv.DictReader(open(p)))
print(f"Cohort threshold calibration from {len(rows)} subject(s) in {p}")
if len(rows) < 5:
print(" WARNING: < 5 subjects - percentile tails are unstable; treat as indicative only.")
for k, side in CALIB.items():
vals = np.array([float(r[k]) for r in rows
if r.get(k) not in (None, "", "nan")], float)
vals = vals[np.isfinite(vals)]
if len(vals) < 3:
continue
med = float(np.median(vals))
if side == "low":
print(f" {k:22s} flag-below: median={med:.3f} p10={np.percentile(vals,10):.3f} "
f"p5={np.percentile(vals,5):.3f} <- set THR near p5")
else:
print(f" {k:22s} flag-above: median={med:.3f} p90={np.percentile(vals,90):.3f} "
f"p95={np.percentile(vals,95):.3f} <- set THR near p95")
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--cohort"); ap.add_argument("--out", default="analysis/qc")
ap.add_argument("--results-root",
default=os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))),
"analysis", "results"),
help="dir holding <id>/roi_efield_<montage>.csv from 04_extract_roi_efield.py")
ap.add_argument("--montage", default="M1", help="which montage's sims to QC (default M1)")
ap.add_argument("--calibrate", action="store_true",
help="read out/qc_summary.csv and suggest cohort-percentile thresholds, then exit")
ap.add_argument("--emit-metrics", help="write the per-subject numeric metrics to this JSON and exit")
args = ap.parse_args()
if args.calibrate:
_calibrate(args.out); return
subjects = resolve_subjects(args)
os.makedirs(args.out, exist_ok=True)
rows, stage_flags = [], []
for P in subjects:
P.setdefault("reg", os.path.join(P["work"], "registration"))
P.setdefault("samseg", os.path.join(P["m2m"], "segmentation", "labeling.nii.gz"))
P.setdefault("cc_mask", os.path.join(_roi_dir(P["reg"]), "roi_CC.nii.gz"))
P.setdefault("montage", args.montage)
P.setdefault("results_root", args.results_root)
row = {"subject": P["id"]}; flags_by_stage = {}
for stage_name, fn in STAGES:
try:
metrics, flags = fn(P)
except Exception as e:
metrics, flags = {}, [f"{stage_name}:exception({type(e).__name__})"]
row.update(metrics)
flags_by_stage[stage_name] = flags
row[f"{stage_name}_PASS"] = "PASS" if not flags else "FLAG"
row["overall"] = "PASS" if all(row[f"{s}_PASS"] == "PASS" for s, _ in STAGES) else "FLAG"
rows.append(row); stage_flags.append(flags_by_stage)
# cohort MAD outlier flagging, gated at n>=4 (MAD is meaningless below that); a lone subject is
# judged ONLY by the absolute THR. The `numeric` list is still built for the CSV.
numeric = sorted({k for r in rows for k, v in r.items() if isinstance(v, (int, float))})
if len(rows) >= 4:
for k in numeric:
vals = np.array([r.get(k, np.nan) for r in rows], float)
med = np.nanmedian(vals); mad = np.nanmedian(np.abs(vals - med)) * 1.4826
if mad > 0:
for i, r in enumerate(rows):
if np.isfinite(vals[i]) and abs(vals[i] - med) / mad > 3:
stage_flags[i].setdefault("MAD", []).append(f"{k}_outlier")
r["overall"] = "FLAG"
# write CSV
cols = ["subject"] + numeric + [f"{s}_PASS" for s, _ in STAGES] + ["overall"]
csv_path = os.path.join(args.out, "qc_summary.csv")
with open(csv_path, "w", newline="") as fh:
w = csv.DictWriter(fh, fieldnames=cols, extrasaction="ignore"); w.writeheader()
for r in rows:
w.writerow(r)
if args.emit_metrics:
snap = {r["subject"]: {k: r[k] for k in numeric if k in r} for r in rows}
with open(args.emit_metrics, "w") as fh:
json.dump(snap, fh, indent=2, sort_keys=True)
print(f"metrics -> {args.emit_metrics}")
return
# PNGs: flagged subjects + random 10% sample
rng = np.random.default_rng(0)
flagged = [i for i, r in enumerate(rows) if r["overall"] == "FLAG"]
sample = set(flagged) | set(rng.choice(len(rows), max(1, len(rows) // 10), replace=False).tolist())
png_dir = os.path.join(args.out, "qc_report"); os.makedirs(png_dir, exist_ok=True)
for i in sorted(sample):
try:
_overlay_png(subjects[i], png_dir)
except Exception as e:
print(f" PNG failed for {subjects[i]['id']}: {type(e).__name__}")
# console summary
print(f"\nQC summary ({len(rows)} subject(s)) -> {csv_path}")
for i, r in enumerate(rows):
if r["overall"] == "FLAG":
fl = [f for s in stage_flags[i].values() for f in s]
print(f" FLAG {r['subject']}: {', '.join(fl)}")
else:
print(f" PASS {r['subject']}")
print(f"PNGs ({len(sample)}: flagged + 10% sample) -> {png_dir}")
def _overlay_png(P, png_dir):
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
t1 = _data(os.path.join(P["m2m"], "T1.nii.gz"))
if t1 is None:
return
lab = _data(os.path.join(_roi_dir(P["reg"]), "roi_labels_meshspace.nii.gz")) # recon-all ROI labels
# center slices on the ROI-label centroid so small ROIs (CC, mesencephalon) are in-plane
if lab is not None and (lab > 0).any():
ci, cj, ck = (int(round(c)) for c in np.array(np.where(lab > 0)).mean(axis=1))
else:
ci, cj, ck = (s // 2 for s in t1.shape[:3])
fig, ax = plt.subplots(1, 3, figsize=(12, 4))
planes = [(t1[ci, :, :], None if lab is None else lab[ci, :, :]),
(t1[:, cj, :], None if lab is None else lab[:, cj, :]),
(t1[:, :, ck], None if lab is None else lab[:, :, ck])]
for a, (bg, ov) in zip(ax, planes):
a.imshow(np.rot90(bg), cmap="gray")
if ov is not None:
a.imshow(np.rot90(np.ma.masked_equal(ov, 0)), cmap="tab10", alpha=0.5, vmin=1, vmax=10)
a.axis("off")
fig.suptitle(f"{P['id']} - T1 + Group 1 volume ROIs"); fig.tight_layout()
fig.savefig(os.path.join(png_dir, f"{P['id']}_overlay.png"), dpi=80); plt.close(fig)
if __name__ == "__main__":
main()

qc_harness.py at commit f4ad4bf, under MIT · at the source

Overview

Authors: Santiago Osorio Jurado^1,2, Mikael Skorpil³, Per Svenningsson⁴, Rodrigo Moreno^2,5, Christoffer Olsson²

Department of Biomedical Engineering, The City College of New York, City University of New York, New York, NY, United States
Department of Biomedical Engineering and Health Systems, KTH Royal Institute of Technology, Stockholm, Sweden
Department of Molecular Medicine and Surgery, Karolinska Institutet, Stockholm, Sweden
Department of Clinical Neuroscience, Karolinska Institutet, Stockholm, Sweden
Department of Neurobiology, Care Sciences and Society, Karolinska Institutet, Stockholm, Sweden

Journal: Frontiers in human neuroscience, volume 20, article 1972676

Dates: received 19 August 2026; accepted 6 October 2026; published online 21 October 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI · PMCID PMC13730951

Status: code verified (how the code is found and verified)

Categories: structural MRI / diffusion (modality), other (modality), human (organism), Parkinson's (population)

Methods: Connectivity, Statistics, Preprocessing, Machine learning, fMRI & imaging, Physiology & signal measures

Keywords: diffusion MRI, electric field modeling, finite element method, Parkinson disease, q-space trajectory imaging, transcranial direct current stimulation

Citations: not cited yet (Europe PMC); 68 references in the paper

Abstract

Introduction: Predictions of the electric field (E-field) induced by transcranial direct current stimulation (tDCS) depend on how brain conductivity is modeled. White matter anisotropy is conventionally estimated from single-shell diffusion tensor imaging (DTI). Multidimensional diffusion MRI (MD-dMRI), specifically q-space trajectory imaging (QTI), instead gives a mean tensor expected to carry less kurtosis bias. Our primary question was whether replacing the conventional single-shell tensor with this mean tensor would change the predicted E-field.

Methods: We built, to our knowledge, the first MD-dMRI tDCS conductivity model and compared it against DTI and isotropic models in 29 participants (12 with Parkinson’s disease and 17 controls) across four montages, with the same mesh, electrodes, and solver.

Results: The two anisotropic models differed mainly in tensor orientation (18–25 degrees across the white matter), with E-field magnitudes differing by 6.6% or less. E-fields were not significantly different between patients and controls in any region or montage (which was an exploratory, underpowered comparison). Whole-brain E-field correlated with MR elastography stiffness (partial r = 0.58) but attenuated to non-significance once cerebrospinal fluid (CSF) morphology was accounted for (r = 0.06 to 0.09).

Discussion: With no ground-truth E-field or conductivity available, the study establishes the feasibility of the MD-dMRI model and characterizes E-field sensitivity rather than improved dosimetry accuracy. The choice of diffusion tensor is second order for the predicted E-field, which is primarily influenced by individual anatomy. For Parkinson’s disease, modeling efforts should focus on CSF- and atrophy-aware head models and dose normalization, rather than a more complex diffusion tensor.

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

santioso-hue/MRE-tDCS-PD

License: MIT

State: the link answers, verified on 7 November 2026

Evidence: files inventoried

Commit: f4ad4bfe2893d5f92fc7ced3a401ff5e44b1ee02, 16 July 2026

Languages: Python (27), Shell (10), MATLAB (1)

Size: 46 files, 38 scripts

Software Heritage: not archived

Found in: “Data availability statement”

Holds: README, license file, CITATION.cff, environment (requirements.txt), tests

Not found: continuous integration, documentation

Tools: NumPy (22 files), NiBabel (12 files), FSL (9 files), SciPy (7 files), FreeSurfer (1 file), Matplotlib (1 file), MRtrix3 (1 file)

Availability: 1 check, the latest on 7 November 2026: the link answers

7 November 2026: the link answers

40 files

analysis/04_extract_roi_efield.py, Python, 128 lines
analysis/05_mre_efield_comparison.py, Python, 182 lines
analysis/06_cohort_stats.py, Python, 204 lines, 2 matches
analysis/07_build_nuclei.sh, Shell, 54 lines
analysis/08_tensor_divergence.py, Python, 131 lines
analysis/09_mre_cohort_corr.py, Python, 278 lines, 2 matches
analysis/_build_nuclei_labels.py, Python, 82 lines
analysis/_rois.py, Python, 164 lines
analysis/_sims.py, Python, 21 lines
analysis/_stats.py, Python, 243 lines
analysis/aggregate_matched_control.py, Python, 240 lines, 2 matches
analysis/build_cohort_manifest.py, Python, 96 lines
analysis/build_rois.py, Python, 173 lines, 1 match
analysis/charm_qc_index.py, Python, 68 lines
analysis/qc_harness.py, Python, 587 lines, 3 matches
analysis/registration_crosscheck.py, Python, 75 lines
config/config.example.sh, Shell, 50 lines, 2 matches
pipeline/00_charm.sh, Shell, 73 lines, 1 match
pipeline/01_dwi2cond.sh, Shell, 67 lines, 1 match
pipeline/02_register_dmri_to_T1.sh, Shell, 119 lines, 1 match
pipeline/03_build_conductivity_tensor.py, Python, 122 lines
pipeline/04_run_simulations.py, Python, 134 lines, 2 matches
pipeline/05_register_mre_to_T1.sh, Shell, 80 lines
pipeline/08b_matched_reg_control.sh, Shell, 111 lines, 1 match
pipeline/_config.py, Python, 35 lines
pipeline/_tensormath.py, Python, 66 lines
pipeline/prepare_dmri_tensor.py, Python, 113 lines, 1 match
pipeline/run_cohort.sh, Shell, 149 lines
pipeline/run_qti_cov_cohort.m, MATLAB, 74 lines, 2 matches
pipeline/run_qti_cov_cohort.sh, Shell, 26 lines
pipeline/run_subject.sh, Shell, 140 lines
tests/test_bootstrap_ci.py, Python, 16 lines
tests/test_cohort_stats.py, Python, 92 lines
tests/test_lobe_grouping.py, Python, 47 lines
tests/test_matched_reg_control.py, Python, 85 lines
tests/test_qc_harness.py, Python, 68 lines
tests/test_tensor_divergence.py, Python, 77 lines
tests/validate_mean_tensor.py, Python, 121 lines
LICENSE, License, 21 lines
README.md, Text, 127 lines

At the source: github.com/santioso-hue/MRE-tDCS-PD

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;
38 scripts, each with its path and the digest of its content;
21 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 statement

The MR imaging data analyzed here are subject to ethical and legal restrictions and cannot be shared publicly. Anonymized data can be made available to researchers who meet the applicable legal and ethical requirements, on reasonable request to the corresponding author. The analysis and modeling code is openly available at https://github.com/santioso-hue/MRE-tDCS-PD.

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

Recorded: type, language, journal, volume, pages, dates, 5 authors, 6 keywords, 67 references.

Cite

This paper

Osorio Jurado, S., Skorpil, M., Svenningsson, P., Moreno, R., & Olsson, C. (2026). Anisotropic conductivity modeling for tDCS in Parkinson’s disease using multidimensional diffusion MRI. Frontiers in human neuroscience, 20, 1972676.

BibTeX

@article{osoriojurado2026anisotropic, author = {Osorio Jurado, Santiago and Skorpil, Mikael and Svenningsson, Per and Moreno, Rodrigo and Olsson, Christoffer}, title = {{Anisotropic conductivity modeling for tDCS in Parkinson’s disease using multidimensional diffusion MRI}}, journal = {Frontiers in human neuroscience}, year = {2026}, month = nov, volume = {20}, pages = {1972676}, publisher = {Frontiers Media SA}, issn = {1662-5161}, pmcid = {PMC13730951} }

RIS

TY - JOUR AU - Osorio Jurado, Santiago AU - Skorpil, Mikael AU - Svenningsson, Per AU - Moreno, Rodrigo AU - Olsson, Christoffer TI - Anisotropic conductivity modeling for tDCS in Parkinson’s disease using multidimensional diffusion MRI T2 - Frontiers in human neuroscience J2 - Front Hum Neurosci PY - 2026 DA - 2026/11/05 VL - 20 SP - 1972676 SN - 1662-5161 PB - Frontiers Media SA LA - en ER -

CSL-JSON

{ "id": "pmcid:PMC13730951", "type": "article-journal", "title": "Anisotropic conductivity modeling for tDCS in Parkinson’s disease using multidimensional diffusion MRI", "container-title": "Frontiers in human neuroscience", "author": [ { "family": "Osorio Jurado", "given": "Santiago" }, { "family": "Skorpil", "given": "Mikael" }, { "family": "Svenningsson", "given": "Per" }, { "family": "Moreno", "given": "Rodrigo" }, { "family": "Olsson", "given": "Christoffer" } ], "container-title-short": "Front Hum Neurosci", "volume": "20", "page": "1972676", "PMCID": "PMC13730951", "ISSN": "1662-5161", "publisher": "Frontiers Media SA", "language": "en", "issued": { "date-parts": [ [ 2026, 11, 5 ] ] } }

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

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.

Add the badge to the code's README

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown

<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/repo/github.com/santioso-hue/mre-tdcs-pd.svg)](https://openscicode.org/badge/repo/github.com/santioso-hue/mre-tdcs-pd)
<!-- oscr-badge:end -->

reStructuredText

.. oscr-badge:start

.. image:: https://openscicode.org/badge/repo/github.com/santioso-hue/mre-tdcs-pd.svg
   :target: https://openscicode.org/badge/repo/github.com/santioso-hue/mre-tdcs-pd
   :alt: OSCR

.. oscr-badge:end

HTML

<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/repo/github.com/santioso-hue/mre-tdcs-pd"><img src="https://openscicode.org/badge/repo/github.com/santioso-hue/mre-tdcs-pd.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

santioso-hue/MRE-tDCS-PD: README.md on GitHub

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge 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: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

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