OSCR

An automated derivative-based method for detection of motor evoked potential onset latencies in multi-muscle transcranial magnetic stimulation studies.

Code ↔ Paper

3 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 3 matches
  1. [1] § Methods › Computational runtime ↔ MEP_latency_derivative_ratio.py, lines 527–604 · score 0.59 · perf counter, Elapsed, export, runtime, epoched, EMG
  2. [2] § Methods › Derivative-ratio algorithm › Algorithm process in detail › Amplitude gate ↔ MEP_latency_derivative_ratio.py, lines 56–168 · score 0.55 · 45 ms, latency detection, amplitude, baseline, peak, gate
  3. [3] § Methods › Derivative-ratio algorithm › Algorithm process in detail › True MEP filter ↔ MEP_latency_derivative_ratio.py, lines 389–521 · score 0.54 · Latency cap, onset latencies, bound, MEP

Paper

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

The paper is loaded when this pane is shown.

The authors' code

Python · 748 lines · 28 KB · no license · 3 matches

  1. #!/usr/bin/env python3
  2. """
  3. advanced_latency_analysis_v11a_light_refactor_v5.py
  4. ----------------------------------------------------
  5. Publication-ready MEP-latency pipeline with optional resampling and
  6. sampling-rate-aware timing parameters.
  8. Detect onsets in EMG epochs stored as NumPy *.npy files and save one wide CSV
  9. per input file (frames × channels, latencies in ms).
  11. Key flags
  12. ---------
  13. --fs 2000 # native input sampling rate
  14. --resample-to-hz 5000 # optional target sampling rate
  15. --task-mode {rest,active,auto} # prestim outlier rule
  16. --parallel 4 # workers; 0 = serial
  18. Sensitivity-analysis flags
  19. --------------------------
  20. --ptp-factor 1.1
  21. --derivative-block-ms 2.5
  22. --search-back-factor 1.75
  23. --derivative-ratio-thresh 0.85
  24. --peak2trough-min-ms 5.0
  25. --peak2trough-max-ms 7.5
  27. Important implementation detail
  28. -------------------------------
  29. Parameters that were previously hard-coded as sample counts are now stored in
  30. milliseconds and converted to samples from cfg.fs. Because cfg.fs is replaced
  31. with the effective sampling rate after optional resampling, these parameters
  32. preserve their approximate duration at either native or resampled rates.
  33. """
  34. from __future__ import annotations
  36. import argparse
  37. import logging
  38. import math
  39. import re
  40. import sys
  41. import time
  42. from dataclasses import dataclass, replace
  43. from pathlib import Path
  44. from typing import Sequence
  46. import numpy as np
  47. import pandas as pd
  48. from joblib import Parallel, delayed
  49. from scipy.ndimage import gaussian_filter1d, uniform_filter1d
  50. from scipy.signal import filtfilt, iirnotch, resample_poly
  52. ###############################################################################
  53. # 1 | config #
  54. ###############################################################################
  56. @dataclass(frozen=True)
  57. class Cfg:
  58. """Configuration for the derivative-ratio latency detector.
  60. Notes
  61. -----
  62. ``fs`` should always represent the sampling rate currently used by the
  63. algorithm. At file-loading time this is the native sampling rate. If
  64. ``resample_to_hz`` is set, ``process_file()`` resamples the data and then
  65. replaces ``cfg.fs`` with the effective sampling rate before filtering,
  66. template construction, and latency detection.
  67. """
  69. # Native/effective sampling and event windows -----------------------------
  70. fs: int = 2000 # Hz before resampling; becomes effective Hz afterwards
  71. stim_on: float = 0.25 # seconds from epoch start to TMS pulse
  72. mep_on: float = 0.005 # seconds relative to pulse
  73. mep_off: float = 0.045 # seconds relative to pulse
  74. prestim_start: float = -0.101 # seconds relative to pulse
  75. prestim_end: float = -0.001 # seconds relative to pulse
  76. resample_to_hz: int | None = None # None = native fs; e.g. 5000 = resample to 5 kHz
  78. # Amplitude / derivative-ratio parameters ---------------------------------
  79. ptp_factor: float = 1.1
  80. derivative_block_ms: float = 2.5
  81. derivative_ratio_thresh: float = 0.85
  82. search_back_factor: float = 1.75
  83. latency_cap: float = 0.035 # seconds relative to pulse
  85. # Peak-to-trough constraints, now in time units ---------------------------
  86. peak2trough_min_ms: float = 5.0
  87. peak2trough_max_ms: float = 7.5
  89. # Filtering / smoothing ---------------------------------------------------
  90. mains_filter: bool = True
  91. notch_q: float = 30.0
  92. smoothing: str | None = "rolling" # {'gaussian', 'rolling', None}
  93. rolling_smooth_ms: float = 2.5
  94. gaussian_sigma_ms: float = 1.0
  96. # Refinement / template-anchor parameters, now in time units --------------
  97. refine_chunk_ms: float = 2.0
  98. tpl_anchor_tol_ms: float = 7.5
  100. # Session/gating knobs ----------------------------------------------------
  101. rms_multiplier: float = 1.5 # window RMS > X × baseline RMS
  102. task_mode: str = "rest" # {'rest', 'active', 'auto'}
  103. active_token: str = "act" # token that marks active files
  105. # Cached indices ----------------------------------------------------------
  106. @property
  107. def idx(self) -> dict[str, int]:
  108. """Return key analysis-window indices for the current effective fs."""
  109. s2i = lambda t: int(round((t + self.stim_on) * self.fs))
  110. return {
  111. "stim": int(round(self.stim_on * self.fs)),
  112. "mep_on": s2i(self.mep_on),
  113. "mep_off": s2i(self.mep_off),
  114. "pre_start": s2i(self.prestim_start),
  115. "pre_end": s2i(self.prestim_end),
  116. }
  118. def ms_to_samples(self, ms: float, *, minimum: int = 1) -> int:
  119. """Convert milliseconds to samples using the current effective fs."""
  120. return max(minimum, int(math.floor((ms / 1000.0) * self.fs + 0.5)))
  122. @property
  123. def derivative_block(self) -> int:
  124. """Derivative-ratio block length in samples."""
  125. return self.ms_to_samples(self.derivative_block_ms)
  127. @property
  128. def peak2trough_min(self) -> int:
  129. """Minimum peak-to-trough interval in samples."""
  130. return self.ms_to_samples(self.peak2trough_min_ms)
  132. @property
  133. def peak2trough_max(self) -> int:
  134. """Maximum peak-to-trough interval in samples."""
  135. return max(self.peak2trough_min, self.ms_to_samples(self.peak2trough_max_ms))
  137. @property
  138. def rolling_smooth_n(self) -> int:
  139. """Rolling smoothing window in samples."""
  140. return self.ms_to_samples(self.rolling_smooth_ms)
  142. @property
  143. def gaussian_sigma_samples(self) -> float:
  144. """Gaussian smoothing sigma in samples."""
  145. return max(1.0, (self.gaussian_sigma_ms / 1000.0) * self.fs)
  147. @property
  148. def refine_chunk(self) -> int:
  149. """Initial refinement chunk length in samples."""
  150. return self.ms_to_samples(self.refine_chunk_ms)
  152. @property
  153. def tpl_anchor_tol(self) -> int:
  154. """Allowed template-anchor tolerance in samples."""
  155. return self.ms_to_samples(self.tpl_anchor_tol_ms)
  157. def derived_sample_summary(self) -> dict[str, int | float]:
  158. """Return derived sample-count parameters for logging/reproducibility."""
  159. return {
  160. "fs": self.fs,
  161. "derivative_block": self.derivative_block,
  162. "peak2trough_min": self.peak2trough_min,
  163. "peak2trough_max": self.peak2trough_max,
  164. "rolling_smooth_n": self.rolling_smooth_n,
  165. "gaussian_sigma_samples": round(self.gaussian_sigma_samples, 3),
  166. "refine_chunk": self.refine_chunk,
  167. "tpl_anchor_tol": self.tpl_anchor_tol,
  168. }
  170. ###############################################################################
  171. # 2 | helpers #
  172. ###############################################################################
  174. def notch_50hz(x: np.ndarray, cfg: Cfg) -> np.ndarray:
  175. """Apply a 50 Hz notch filter along the sample axis if requested."""
  176. if not cfg.mains_filter:
  177. return x
  179. nyquist = cfg.fs / 2.0
  180. if 50.0 >= nyquist:
  181. raise ValueError(f"Cannot apply 50 Hz notch when Nyquist is only {nyquist:.1f} Hz")
  183. w0 = 50.0 / nyquist
  184. b, a = iirnotch(w0, cfg.notch_q)
  185. return filtfilt(b, a, x, axis=0)
  188. def resample_emg_block(emg: np.ndarray, cfg: Cfg) -> tuple[np.ndarray, int]:
  189. """Optionally resample an epoched EMG block along the sample axis.
  191. Parameters
  192. ----------
  193. emg
  194. Epoched EMG block shaped ``samples × frames × channels``.
  195. cfg
  196. Configuration object. ``cfg.fs`` is interpreted as the native sampling
  197. rate before resampling, and ``cfg.resample_to_hz`` is the optional target.
  199. Returns
  200. -------
  201. tuple[np.ndarray, int]
  202. The possibly resampled EMG block and the effective sampling rate to use
  203. for all downstream indexing, filtering, and latency conversion.
  204. """
  205. target_fs = cfg.resample_to_hz
  206. original_fs = int(cfg.fs)
  208. if target_fs is None or int(target_fs) == original_fs:
  209. return emg, original_fs
  211. target_fs = int(target_fs)
  212. if target_fs <= 0:
  213. raise ValueError("resample_to_hz must be a positive integer or None")
  215. common = math.gcd(original_fs, target_fs)
  216. up = target_fs // common
  217. down = original_fs // common
  219. emg_rs = resample_poly(emg, up=up, down=down, axis=0)
  220. return emg_rs, target_fs
  223. def smooth(sig: np.ndarray, cfg: Cfg) -> np.ndarray:
  224. """Smooth an EMG trace using windows defined in milliseconds.
  226. The actual sample width is derived from ``cfg.fs``, so smoothing is
  227. preserved in time units if the signal is resampled.
  228. """
  229. if cfg.smoothing == "gaussian":
  230. return gaussian_filter1d(sig, cfg.gaussian_sigma_samples, axis=0)
  232. if cfg.smoothing == "rolling":
  233. n = cfg.rolling_smooth_n
  234. return uniform_filter1d(uniform_filter1d(sig, n, axis=0), n, axis=0)
  236. return sig
  239. def prestim_mask(baseline: np.ndarray, mode: str) -> np.ndarray:
  240. """Return a Boolean frame mask based on pre-stimulus RMS outliers.
  242. Parameters
  243. ----------
  244. baseline
  245. Array shaped ``samples_in_prestim × n_frames``.
  246. mode
  247. ``'rest'`` or ``'active'``. Active sessions use a looser z-score band.
  248. """
  249. if baseline.shape[0] == 0:
  250. return np.ones(baseline.shape[1], dtype=bool)
  252. rms = np.sqrt((baseline ** 2).mean(axis=0))
  253. mu, sd = rms.mean(), rms.std(ddof=0)
  254. if sd == 0 or not np.isfinite(sd):
  255. return np.ones_like(rms, dtype=bool)
  257. thr = 3.5 if mode == "active" else 2.0
  258. z = (rms - mu) / sd
  259. return np.abs(z) <= thr
  262. def is_active_file(name: str, token: str) -> bool:
  263. """Return True if ``token`` appears as a standalone chunk in ``name``."""
  264. s = name.lower()
  265. pat = rf"(?<![a-z]){re.escape(token.lower())}(?![a-z])"
  266. return re.search(pat, s) is not None
  268. ###############################################################################
  269. # 3 | template builder #
  270. ###############################################################################
  272. def build_templates(emg: np.ndarray, chans: Sequence[str], cfg: Cfg):
  273. """Build one normalised average MEP template per channel."""
  274. idx, tpl = cfg.idx, {}
  275. for i, ch in enumerate(chans):
  276. sig = emg[:, :, i]
  277. base = sig[idx["pre_start"]:idx["pre_end"]]
  278. mep = sig[idx["mep_on"]:idx["mep_off"]]
  280. if base.shape[0] == 0 or mep.shape[0] == 0:
  281. tpl[ch] = None
  282. continue
  284. good = prestim_mask(base, cfg.task_mode)
  285. good &= np.ptp(mep, axis=0) > cfg.ptp_factor * np.ptp(base, axis=0)
  287. if not good.any():
  288. tpl[ch] = None
  289. continue
  291. waves = mep[:, good]
  292. waves = (waves - waves.mean(0)) / (waves.std(0) + 1e-9)
  293. tpl[ch] = waves.mean(1)
  295. return tpl
  297. ###############################################################################
  298. # 4 | candidate refinement #
  299. ###############################################################################
  301. def _refine(sig: np.ndarray, diff: np.ndarray, cand: int, start: int,
  302. mean_d: float, std_d: float, p2t: int, base_rms: float,
  303. cfg: Cfg) -> int | None:
  304. """Refine a candidate onset using local derivative/RMS criteria.
  306. This is a NumPy equivalent of the earlier pandas/DataFrame implementation.
  307. It keeps the same decision logic but avoids repeated DataFrame construction
  308. and ``.iloc`` slicing inside the per-frame loop.
  309. """
  310. chunk = cfg.refine_chunk
  311. win = p2t * 2
  312. span = max(1, p2t // 4)
  313. min_neg_samples = max(1, math.ceil(0.75 * chunk))
  314. mep_on = cfg.idx["mep_on"]
  316. def ok(j: int) -> bool:
  317. sig_w = sig[j:j + win]
  318. diff_w = diff[j:j + win]
  319. if sig_w.size == 0 or diff_w.size == 0:
  320. return False
  322. # Previous code used: d = mean_d - diff; cond = d < 0.
  323. # This is equivalent to diff > mean_d.
  324. cond = diff_w > mean_d
  325. rms = math.sqrt(float(np.mean(sig_w ** 2)))
  326. first_chunk = cond[:chunk]
  328. return (
  329. float(np.mean(cond)) > 0.5
  330. and int(np.sum(first_chunk)) >= min_neg_samples
  331. and (
  332. rms > cfg.rms_multiplier * base_rms
  333. or float(np.mean(diff_w)) > mean_d + 1.5 * std_d
  334. )
  335. )
  337. if ok(cand):
  338. return cand
  340. lower = max(cand - span, mep_on)
  341. for j in range(cand - 1, lower - 1, -1):
  342. if ok(j):
  343. return j
  345. for j in range(cand + 1, cand + span + 1):
  346. if j > start:
  347. break
  348. if ok(j):
  349. return j
  351. return None
  353. ###############################################################################
  354. # 5 | per-channel pipeline #
  355. ###############################################################################
  357. def _window_means_from_cumsum(values: np.ndarray, starts: np.ndarray,
  358. width: int) -> np.ndarray:
  359. """Return fixed-width window means for many start indices.
  361. Parameters
  362. ----------
  363. values
  364. One-dimensional numeric array.
  365. starts
  366. Start indices for each window.
  367. width
  368. Number of samples in each window.
  370. Returns
  371. -------
  372. np.ndarray
  373. Mean value for ``values[start:start + width]`` for each start.
  375. Notes
  376. -----
  377. This replaces repeated small pandas/Python slices in the derivative-ratio
  378. scan. ``diff`` can contain a NaN at index 0 from ``np.diff(..., prepend)``;
  379. that value is converted to zero before cumulative sums. The derivative-ratio
  380. search windows are post-stimulus and should not include index 0 in normal
  381. use, so this conversion preserves practical behaviour while keeping the
  382. vectorised implementation robust.
  383. """
  384. clean = np.nan_to_num(values, nan=0.0)
  385. cs = np.concatenate(([0.0], np.cumsum(clean, dtype=float)))
  386. return (cs[starts + width] - cs[starts]) / width
  389. def process_channel(sig: np.ndarray, tpl: np.ndarray | None, cfg: Cfg) -> np.ndarray:
  390. """Detect onset latencies for all frames in one channel.
  392. Optimisation notes
  393. ------------------
  394. The original implementation created a pandas ``DataFrame`` for every kept
  395. frame and then used repeated ``.iloc`` slices during the derivative-ratio
  396. search and refinement checks. This version keeps the same algorithmic steps
  397. but performs the inner-loop operations with NumPy arrays. It also smooths all
  398. kept frames for a channel in one vectorised call, rather than smoothing each
  399. frame separately.
  400. """
  401. idx, nF = cfg.idx, sig.shape[1]
  402. lat = np.full(nF, np.nan, dtype=object)
  404. if tpl is None:
  405. lat[:] = "NaN"
  406. return lat
  408. base = sig[idx["pre_start"]:idx["pre_end"]]
  409. mep = sig[idx["mep_on"]:idx["mep_off"]]
  411. if base.shape[0] == 0 or mep.shape[0] == 0:
  412. lat[:] = "NaN"
  413. return lat
  415. good_frames = prestim_mask(base, cfg.task_mode)
  416. lat[~good_frames] = "NaN"
  418. keep = good_frames & (np.ptp(mep, axis=0) > cfg.ptp_factor * np.ptp(base, axis=0))
  419. lat[~keep & good_frames] = "NaN"
  421. keep_idx = np.flatnonzero(keep)
  422. if keep_idx.size == 0:
  423. return lat
  425. tpl_anchor = min(np.argmax(tpl), np.argmin(tpl)) + idx["mep_on"]
  426. db = cfg.derivative_block
  427. eps = 1e-6
  429. # Smooth all kept frames at once. This is mathematically equivalent to
  430. # smoothing each frame independently because smoothing is along axis 0 only.
  431. sig_keep = sig[:, keep_idx]
  432. smoothed_keep = smooth(sig_keep, cfg)
  434. for local_col, f in enumerate(keep_idx):
  435. frame = sig_keep[:, local_col]
  436. smoothed = smoothed_keep[:, local_col]
  437. diff = np.abs(np.diff(smoothed, prepend=np.nan))
  439. pre_d = diff[idx["pre_start"]:idx["pre_end"]]
  440. mean_d = float(np.nanmean(pre_d))
  441. std_d = float(np.nanstd(pre_d))
  443. mwin = smoothed[idx["mep_on"]:idx["mep_off"]]
  444. if mwin.size == 0 or np.all(np.isnan(mwin)):
  445. lat[f] = "null_onset"
  446. continue
  448. peak = int(np.nanargmax(mwin) + idx["mep_on"])
  449. trough = int(np.nanargmin(mwin) + idx["mep_on"])
  450. p2t = int(np.clip(abs(peak - trough), cfg.peak2trough_min, cfg.peak2trough_max))
  451. start_idx = min(peak, trough)
  453. if not (tpl_anchor - cfg.tpl_anchor_tol <= start_idx <= tpl_anchor + cfg.tpl_anchor_tol):
  454. lat[f] = "null_onset"
  455. continue
  457. lower = max(start_idx - int(p2t * cfg.search_back_factor), idx["mep_on"])
  458. first_i = start_idx - db
  459. if first_i < lower:
  460. lat[f] = "null_onset"
  461. continue
  463. # Preserve the original scan order: start_idx - db, then step backward.
  464. i_vals = np.arange(first_i, lower - 1, -1, dtype=int)
  465. prev_starts = i_vals - db
  466. nxt_starts = i_vals
  468. # Defensive bounds check. In normal use the lower bound prevents this.
  469. valid = (prev_starts >= 0) & (nxt_starts + db <= diff.size)
  470. if not np.any(valid):
  471. lat[f] = "null_onset"
  472. continue
  474. i_vals = i_vals[valid]
  475. prev_starts = prev_starts[valid]
  476. nxt_starts = nxt_starts[valid]
  478. prev = _window_means_from_cumsum(diff, prev_starts, db)
  479. nxt = _window_means_from_cumsum(diff, nxt_starts, db)
  480. ratios = np.abs(nxt) / (np.abs(prev) + eps)
  482. if ratios.size == 0:
  483. lat[f] = "null_onset"
  484. continue
  486. max_pos = int(np.argmax(ratios))
  487. mx = int(i_vals[max_pos])
  488. mxr = float(ratios[max_pos])
  489. thr = cfg.derivative_ratio_thresh * mxr
  491. # Recreate the previous neighbour expansion around mx using a dictionary
  492. # for exact key lookup. This avoids subtle changes in candidate ordering.
  493. ratio_lookup = dict(zip(i_vals.tolist(), ratios.tolist()))
  494. cands = [mx]
  496. j = mx - 1
  497. while ratio_lookup.get(j, 0) >= thr:
  498. cands.append(j)
  499. j -= 1
  501. j = mx + 1
  502. while ratio_lookup.get(j, 0) >= thr:
  503. cands.append(j)
  504. j += 1
  506. base_rms = math.sqrt(float(np.mean(frame[idx["pre_start"]:idx["pre_end"]] ** 2)))
  507. refined = None
  509. for c in cands:
  510. if (c - idx["stim"]) / cfg.fs > cfg.latency_cap:
  511. continue
  512. refined = _refine(frame, diff, c, start_idx, mean_d, std_d, p2t, base_rms, cfg)
  513. if refined is not None:
  514. break
  516. lat[f] = (
  517. "null_onset" if refined is None
  518. else round((refined - idx["stim"]) / cfg.fs * 1000, 3)
  519. )
  521. return lat
  523. ###############################################################################
  524. # 6 | file pipeline #
  525. ###############################################################################
  527. def process_file(path: Path, out_dir: Path, chans: Sequence[str], cfg: Cfg,
  528. log: logging.Logger):
  529. """Process one ``.npy`` EMG block, write a wide latency CSV, and report runtime.
  531. Runtime is measured per input file and includes loading, optional resampling,
  532. filtering, template construction, latency detection, DataFrame construction,
  533. and CSV export. The printed per-frame value treats each stimulation frame as
  534. one trial containing all recorded channels. The per-channel-frame value gives
  535. a rough per-epoch speed, where one epoch is one frame from one channel.
  536. """
  537. file_timer_start = time.perf_counter()
  538. log.info("Processing %s", path.name)
  540. emg = np.load(path)
  541. if emg.ndim == 2:
  542. emg = emg[:, :, None]
  544. emg, effective_fs = resample_emg_block(emg, cfg)
  545. if effective_fs != cfg.fs:
  546. log.info("↳ resampled %s from %d Hz to %d Hz", path.name, cfg.fs, effective_fs)
  548. if cfg.task_mode == "auto":
  549. file_mode = "active" if is_active_file(path.stem, cfg.active_token) else "rest"
  550. else:
  551. file_mode = cfg.task_mode
  553. eff_cfg = replace(
  554. cfg,
  555. fs=effective_fs,
  556. task_mode=file_mode,
  557. # Preserve the existing behaviour: stricter RMS gate at rest, looser during contraction.
  558. rms_multiplier=(2.0 if file_mode == "rest" else 1.5),
  559. )
  561. log.debug("Task mode for %s → %s", path.name, file_mode)
  562. log.info("↳ effective config: %s", eff_cfg.derived_sample_summary())
  564. emg = notch_50hz(emg, eff_cfg)
  565. tpls = build_templates(emg, chans, eff_cfg)
  567. out = {
  568. ch: process_channel(emg[:, :, i], tpls[ch], eff_cfg)
  569. for i, ch in enumerate(chans)
  570. }
  572. df = pd.DataFrame(out)
  573. df.index = np.arange(1, len(df) + 1)
  574. out_csv = out_dir / f"{path.stem}_latencies.csv"
  575. df.to_csv(out_csv, index_label="frame")
  576. log.info("↳ saved %s", out_csv.name)
  578. elapsed_s = time.perf_counter() - file_timer_start
  579. n_frames = int(df.shape[0])
  580. n_channels = int(df.shape[1])
  581. n_channel_frames = n_frames * n_channels
  583. sec_per_frame = elapsed_s / n_frames if n_frames else float("nan")
  584. sec_per_channel_frame = (
  585. elapsed_s / n_channel_frames if n_channel_frames else float("nan")
  586. )
  588. runtime_msg = (
  589. f"↳ runtime {path.name}: {elapsed_s:.3f} s total | "
  590. f"{sec_per_frame * 1000:.2f} ms/frame | "
  591. f"{sec_per_channel_frame * 1000:.2f} ms/channel-frame "
  592. f"({n_frames} frames × {n_channels} channels)"
  593. )
  594. log.info(runtime_msg)
  596. return {
  597. "file": path.name,
  598. "seconds": elapsed_s,
  599. "frames": n_frames,
  600. "channels": n_channels,
  601. "channel_frames": n_channel_frames,
  602. "seconds_per_frame": sec_per_frame,
  603. "seconds_per_channel_frame": sec_per_channel_frame,
  604. }
  606. ###############################################################################
  607. # 7 | CLI #
  608. ###############################################################################
  610. def cli_args() -> argparse.Namespace:
  611. """Parse command-line arguments."""
  612. p = argparse.ArgumentParser(description="MEP latency detector (wide CSV)")
  614. p.add_argument("--in-dir", required=True, type=Path)
  615. p.add_argument("--out-dir", required=True, type=Path)
  616. p.add_argument("--channels", required=True, type=Path)
  618. p.add_argument("--fs", type=int, default=Cfg.fs,
  619. help=f"Native input sampling rate in Hz before optional resampling (default {Cfg.fs}).")
  620. p.add_argument("--resample-to-hz", type=int, default=Cfg.resample_to_hz,
  621. help="Optional target sampling rate in Hz before latency detection, e.g. 5000. "
  622. "Omit to analyse at native --fs.")
  624. p.add_argument("--parallel", type=int, default=0)
  625. p.add_argument("--log", default="INFO", choices=["DEBUG", "INFO", "WARNING", "ERROR"])
  627. p.add_argument("--task-mode", default=Cfg.task_mode, choices=["rest", "active", "auto"],
  628. help="Prestim outlier rule: 'rest', 'active', or 'auto' per file")
  629. p.add_argument("--active-token", default=Cfg.active_token,
  630. help="Token that marks active files when --task-mode=auto")
  632. p.add_argument("--ptp-factor", type=float, default=Cfg.ptp_factor,
  633. help="Amplitude gate: MEP ptp must exceed this × baseline ptp.")
  634. p.add_argument("--derivative-block-ms", type=float, default=Cfg.derivative_block_ms,
  635. help="Derivative-ratio block length in milliseconds.")
  636. p.add_argument("--derivative-ratio-thresh", type=float, default=Cfg.derivative_ratio_thresh,
  637. help="Candidate plateau threshold as a fraction of Rmax.")
  638. p.add_argument("--search-back-factor", type=float, default=Cfg.search_back_factor,
  639. help="Search-back limit as a multiple of peak-to-trough distance.")
  640. p.add_argument("--peak2trough-min-ms", type=float, default=Cfg.peak2trough_min_ms,
  641. help="Minimum peak-to-trough interval in milliseconds.")
  642. p.add_argument("--peak2trough-max-ms", type=float, default=Cfg.peak2trough_max_ms,
  643. help="Maximum peak-to-trough interval in milliseconds.")
  645. p.add_argument("--smoothing", default=Cfg.smoothing, choices=["rolling", "gaussian", "none"],
  646. help="Smoothing mode. Use 'none' to disable smoothing.")
  647. p.add_argument("--rolling-smooth-ms", type=float, default=Cfg.rolling_smooth_ms,
  648. help="Rolling smoothing window in milliseconds.")
  649. p.add_argument("--gaussian-sigma-ms", type=float, default=Cfg.gaussian_sigma_ms,
  650. help="Gaussian smoothing sigma in milliseconds.")
  651. p.add_argument("--refine-chunk-ms", type=float, default=Cfg.refine_chunk_ms,
  652. help="Initial refinement chunk length in milliseconds.")
  653. p.add_argument("--tpl-anchor-tol-ms", type=float, default=Cfg.tpl_anchor_tol_ms,
  654. help="Allowed template-anchor tolerance in milliseconds.")
  655. p.add_argument("--rms-multiplier", type=float, default=Cfg.rms_multiplier,
  656. help="Base RMS multiplier. Existing rest/active replacement is preserved in process_file().")
  658. return p.parse_args()
  661. def main(ns: argparse.Namespace | None = None):
  662. """Run batch latency detection."""
  663. ns = ns or cli_args()
  664. logging.basicConfig(level=getattr(logging, ns.log),
  665. format="%(levelname)s:%(name)s:%(message)s")
  666. log = logging.getLogger("latency")
  668. smoothing = None if ns.smoothing == "none" else ns.smoothing
  670. cfg = Cfg(
  671. fs=ns.fs,
  672. resample_to_hz=ns.resample_to_hz,
  673. task_mode=ns.task_mode,
  674. active_token=ns.active_token,
  675. rms_multiplier=ns.rms_multiplier,
  676. ptp_factor=ns.ptp_factor,
  677. derivative_block_ms=ns.derivative_block_ms,
  678. derivative_ratio_thresh=ns.derivative_ratio_thresh,
  679. search_back_factor=ns.search_back_factor,
  680. peak2trough_min_ms=ns.peak2trough_min_ms,
  681. peak2trough_max_ms=ns.peak2trough_max_ms,
  682. smoothing=smoothing,
  683. rolling_smooth_ms=ns.rolling_smooth_ms,
  684. gaussian_sigma_ms=ns.gaussian_sigma_ms,
  685. refine_chunk_ms=ns.refine_chunk_ms,
  686. tpl_anchor_tol_ms=ns.tpl_anchor_tol_ms,
  687. )
  689. chans = np.load(ns.channels, allow_pickle=True).tolist()
  690. ns.out_dir.mkdir(parents=True, exist_ok=True)
  692. files = sorted(ns.in_dir.glob("*.npy"))
  693. if not files:
  694. log.error("No .npy found in %s", ns.in_dir)
  695. sys.exit(1)
  697. def _safe_process(f: Path):
  698. try:
  699. return process_file(f, ns.out_dir, chans, cfg, log)
  700. except Exception as e:
  701. log.error("FAILED %s: %s", f.name, e, exc_info=True)
  702. (ns.out_dir / f"{f.stem}_latencies.csv").write_text("frame\n")
  703. return None
  705. if ns.parallel:
  706. runtime_results = Parallel(n_jobs=ns.parallel)(
  707. delayed(_safe_process)(f) for f in files
  708. )
  709. else:
  710. runtime_results = []
  711. for f in files:
  712. runtime_results.append(_safe_process(f))
  714. runtime_results = [r for r in runtime_results if isinstance(r, dict)]
  715. if runtime_results:
  716. total_seconds = sum(r["seconds"] for r in runtime_results)
  717. total_frames = sum(r["frames"] for r in runtime_results)
  718. total_channel_frames = sum(r["channel_frames"] for r in runtime_results)
  720. log.info(
  721. "Runtime summary: %.3f s total | %.2f ms/frame | %.2f ms/channel-frame "
  722. "(%d files, %d frames, %d channel-frames)",
  723. total_seconds,
  724. (total_seconds / total_frames * 1000) if total_frames else float("nan"),
  725. (total_seconds / total_channel_frames * 1000) if total_channel_frames else float("nan"),
  726. len(runtime_results),
  727. total_frames,
  728. total_channel_frames,
  729. )
  731. ###############################################################################
  732. # 8 | Spyder fallback #
  733. ###############################################################################
  735. if __name__ == "__main__":
  736. if len(sys.argv) == 1: # launched via F5 in Spyder
  737. sys.argv += [
  738. "--in-dir", r"C:\Users\rowbi\OneDrive - Imperial College London\BRC PhD Fellowship\Code\Neuromap\data\control_data\Practice",
  739. "--out-dir", r"C:\Users\rowbi\OneDrive - Imperial College London\BRC PhD Fellowship\Code\Neuromap\data\control_data\Practice",
  740. "--channels", "data/channels.npy",
  741. "--fs", "2000", # native input sampling rate
  742. # "--resample-to-hz", "5000", # remove these two entries to analyse at native --fs
  743. "--parallel", "4",
  744. "--task-mode", "auto",
  745. "--active-token", "act", # string to look for in filename denoting active trials
  746. "--rms-multiplier", "1.5",
  747. ]
  748. main()

MEP_latency_derivative_ratio.py at commit 5f95bbc, no license · at the source

Overview

Authors: Rowan Boyles1,2, Mikal Vicente1, Sofia Zibordi1, Jason Mallabone3, Paul H Strutton1
ORCID iDs: Rowan Boyles
  1. The Nick Davey Laboratory, Department of Surgery and Cancer, Imperial College London, London, UK
  2. Centre for Vestibular Neurology, Department of Brain Sciences, Imperial College London, London, UK
  3. Imperial College Healthcare NHS Trust, London, UK
Institutions: Imperial College London (United Kingdom); Imperial College Healthcare NHS Trust (United Kingdom)
Journal: Scientific reports, volume 16, issue 1, article 26643
Dates: received 15 April 2026; accepted 6 July 2026; published online 23 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41598-026-61560-0 · PMID 42642547 · PMCID PMC13507071 · OpenAlex W7170190401
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: other (modality), human (organism)
Methods: Spectral & time-frequency, Connectivity, Preprocessing, Evoked potentials, Physiology & signal measures, Statistics
Keywords: Transcranial magnetic stimulation, Motor cortex mapping, Motor evoked potentials, Upper limb, MEP latency, Automated detection, Neurology, Neuroscience
MeSH: Evoked Potentials, Motor*, Motor Cortex*, Muscle, Skeletal*, Transcranial Magnetic Stimulation*, Adult, Algorithms, Electromyography, Female, Humans, Male, Pilot Projects, Reaction Time, Reproducibility of Results, Young Adult (* major topic)
Topic: Transcranial Magnetic Stimulation Studies (Neurology, Neuroscience), according to OpenAlex
Funding: National Institute for Health and Care Research; Imperial College London; NIHR Imperial Biomedical Research Centre
Citations: not cited yet (Europe PMC); 35 references in the paper

Abstract

Motor evoked potential (MEP) onset latency is a useful neurophysiological measure, but manual measurement is time-consuming in transcranial magnetic stimulation (TMS) studies with large numbers of trials. This is particularly relevant in cortical mapping studies recording from multiple muscles simultaneously, where automated methods could support more scalable analyses. Existing onset-detection methods have shown promise in more restricted datasets, but their performance in heterogeneous multi-muscle mapping data remains uncertain. In this pilot validation study, three healthy adults underwent TMS cortical mapping with simultaneous EMG recording from eight upper-limb muscles during resting and active conditions, yielding 3,840 EMG epochs. Three independent raters classified MEP presence and marked onset latency for all trials. Inter-rater agreement was assessed using Fleiss’ κ for detection and ICC(2,1) for latency. Human majority vote for MEP presence and mean latency ratings were used as the reference standard to benchmark a novel derivative-ratio algorithm against an existing method. Human raters showed moderate to strong agreement for MEP detection (Fleiss’ κ = 0.69) and high reliability for latency ratings (ICC(2,1) = 0.95), with a pooled mean absolute pairwise difference of 0.88 ms. The derivative-ratio algorithm showed strong detection performance and human-like latency estimates, outperforming a previously published algorithm. These pilot validation data suggest that the derivative-ratio method provides promising automatic MEP onset latency detection in complex multi-muscle cortical mapping data and may provide a basis for scalable latency analysis. By enabling reproducible extraction of conduction-related MEP features, this approach may support future biomarker studies in neurological disorders characterised by altered corticospinal excitability or corticospinal conduction, pending further validation in larger and clinically diverse datasets.

Supplementary Information: The online version contains supplementary material available at 10.1038/s41598-026-61560-0.

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

Repositories

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

Zenodo 16966802

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Data availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (1 file), pandas (1 file), SciPy (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
2 files
At the source: doi.org/10.5281/zenodo.16966802

rowbe1/mep_latency

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 5f95bbca8f62ef92b7a53aef8ac6e20d55a0f61c, 27 May 2026
Languages: Python (1)
Size: 10 files, 1 script
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README, CITATION.cff
Not found: license file, environment file, tests, continuous integration, documentation
Tools: NumPy (1 file), pandas (1 file), SciPy (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
2 files
At the source: github.com/Rowbe1/MEP_latency/tree/v1.0.1

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

Tracing map

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

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 2 scripts, each with its path and the digest of its content;
  • 3 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 Python code implementing the derivative-ratio detector, together with anonymised example EMG, is openly available at https://doi.org/10.5281/zenodo.16966802.

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 2, 28 September 2026

  • Funding: added National Institute for Health and Care Research; Imperial College London; NIHR Imperial Biomedical Research Centre

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 8 keywords, 14 MeSH terms, 35 references.

Cite

This paper

Boyles, R., Vicente, M., Zibordi, S., Mallabone, J., & Strutton, P. H. (2026). An automated derivative-based method for detection of motor evoked potential onset latencies in multi-muscle transcranial magnetic stimulation studies. Scientific reports, 16(1), 26643. https://doi.org/10.1038/s41598-026-61560-0

BibTeX

@article{boyles2026automated,
author = {Boyles, Rowan and Vicente, Mikal and Zibordi, Sofia and Mallabone, Jason and Strutton, Paul H},
title = {{An automated derivative-based method for detection of motor evoked potential onset latencies in multi-muscle transcranial magnetic stimulation studies}},
journal = {Scientific reports},
year = {2026},
month = jul,
volume = {16},
number = {1},
pages = {26643},
publisher = {Nature Publishing Group},
issn = {2045-2322},
doi = {10.1038/s41598-026-61560-0},
url = {https://doi.org/10.1038/s41598-026-61560-0},
pmid = {42642547},
pmcid = {PMC13507071}
}

RIS

TY - JOUR
AU - Boyles, Rowan
AU - Vicente, Mikal
AU - Zibordi, Sofia
AU - Mallabone, Jason
AU - Strutton, Paul H
TI - An automated derivative-based method for detection of motor evoked potential onset latencies in multi-muscle transcranial magnetic stimulation studies
T2 - Scientific reports
J2 - Sci Rep
PY - 2026
DA - 2026/07/23
VL - 16
IS - 1
SP - 26643
SN - 2045-2322
PB - Nature Publishing Group
DO - 10.1038/s41598-026-61560-0
UR - https://doi.org/10.1038/s41598-026-61560-0
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41598-026-61560-0",
"type": "article-journal",
"title": "An automated derivative-based method for detection of motor evoked potential onset latencies in multi-muscle transcranial magnetic stimulation studies",
"container-title": "Scientific reports",
"author": [
{
"family": "Boyles",
"given": "Rowan"
},
{
"family": "Vicente",
"given": "Mikal"
},
{
"family": "Zibordi",
"given": "Sofia"
},
{
"family": "Mallabone",
"given": "Jason"
},
{
"family": "Strutton",
"given": "Paul H"
}
],
"container-title-short": "Sci Rep",
"volume": "16",
"issue": "1",
"page": "26643",
"DOI": "10.1038/s41598-026-61560-0",
"PMID": "42642547",
"PMCID": "PMC13507071",
"ISSN": "2045-2322",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41598-026-61560-0",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
23
]
]
}
}

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.1162/imag.a.1211 [code]
Effects of electric field direction on TMS-based motor cortex mapping.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: pandas, SciPy, NumPy, other, 2 references
[2] doi:10.1002/advs.202523009 [code]
Personalized Network-Guided Neuromodulation Enhances Human Working Memory.
Journal: Advanced science (Weinheim, Baden-Wurttemberg, Germany)
In common: pandas, SciPy, NumPy, other, 1 reference
[3] doi:10.1371/journal.pcbi.1014154 [code]
Complexity of resting cortical activity predicts neurophysiological responses to theta-burst stimulation but fails to generalize: A rigorous machine-learning approach.
Journal: PLoS computational biology
In common: pandas, SciPy, NumPy, other, 1 reference
[4] doi:10.1016/j.nicl.2026.104012 [code]
Structural-functional multilayer brain network properties and outcome of combined repetitive transcranial magnetic stimulation and psychotherapy for obsessive-compulsive disorder.
Journal: NeuroImage. Clinical
In common: pandas, SciPy, NumPy, other, 1 reference
[5] doi:10.1016/j.ebiom.2026.106293 [code]
Dynamic neural states underpin motor symptom severity in Parkinson's disease: a longitudinal analysis of chronic cortico-subthalamic nucleus recordings.
Journal: EBioMedicine
In common: pandas, SciPy, NumPy, 1 reference
[6] doi:10.1162/nol.a.244 [code]
A Novel Approach to Map the Causal Impact of Brain Stimulation on Semantic Processing With Language Models.
Journal: Neurobiology of language (Cambridge, Mass.)
In common: SciPy, NumPy, other, 1 reference
[7] doi:10.1113/jp291217 [code]
Transcranial direct current stimulation enhances delayed retention after 5 days of lower-limb motor skill learning.
Journal: The Journal of physiology
In common: NumPy, other, 1 reference
[8] doi:10.1113/jp290164 [code]
Improved subjective sleep quality in older adults by enhancing the GABAergic system in the sensorimotor cortex.
Journal: The Journal of physiology
In common: other, 1 reference
[9] doi:10.1111/ejn.70597
Superficial Ventral Premotor Pathways to Primary Motor Cortex Shape the Temporal Coordination of Precision Graspinnog.
Journal: The European journal of neuroscience
In common: other, 1 reference
[10] doi:10.1073/pnas.2604933123
Neuromotor modules revealed by direct electrical stimulation of the human primary motor cortex.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: other, 1 reference

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

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

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

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