OSCR

Effects of Transcranial Direct Current Stimulation over the Left Sensorimotor Cortex on Bimanual Force Control: A Computational and Experimental Investigation.

Code ↔ Paper

10 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 10 matches
  1. [1] § 2. Materials and Methods › 2.8. Computational Model › 2.8.1. State Equations ↔ tDCS_motorcontrol_03012026.py, lines 348–493 · score 0.89 · delayed proprioceptive, K_base, target drift, proprioceptive feedback gain, internal target, G_proprio
  2. [2] § 2. Materials and Methods › 2.8. Computational Model › 2.8.1. State Equations ↔ tDCS_motorcontrol_03012026.py, lines 348–493 · score 0.84 · K_vis, common drive weight, visual feedback gain, 10.2 Hz, noise, tremor
  3. [3] § 2. Materials and Methods › 2.1. Study Design ↔ tDCS_motorcontrol_03012026.py, lines 521–647 · score 0.77 · electrode montage, force matching task, visual feedback, anode, cathode, window
  4. [4] § 2. Materials and Methods › 2.7. Statistical Analysis ↔ tDCS_motorcontrol_03012026.py, lines 930–1007 · score 0.63 · post hoc, retest reliability, Epoch interaction, change scores, Cohen, Holm
  5. [5] § 2. Materials and Methods › 2.5. Outcome Measures › 2.5.3. Inter-Hand Coherence ↔ tDCS_motorcontrol_03012026.py, lines 650–732 · score 0.58 · 7–12 Hz, 3–7 Hz, 0–1 Hz, bands, squared, coherence
  6. [6] § 3. Results › 3.2. Computational Model Results ↔ tDCS_motorcontrol_03012026.py, lines 848–928 · score 0.57 · Dose response, G_proprio, corrective power, computationally, N2, model
  7. [7] § 2. Materials and Methods › 2.6. Primary and Secondary Outcomes ↔ tDCS_motorcontrol_03012026.py, lines 650–732 · score 0.52 · inter hand coherence, force undershoot, spectral power, band, 1–3 Hz, RMSE
  8. [8] § 2. Materials and Methods › 2.7. Statistical Analysis ↔ tDCS_motorcontrol_03012026.py, lines 244–268 · score 0.52 · Pearson correlations, change scores, RMSE, PRE, power, undershoot
  9. [9] § 2. Materials and Methods › 2.3. Transcranial Direct Current Stimulation Protocol ↔ tDCS_motorcontrol_03012026.py, lines 1–83 · score 0.52 · left sensorimotor cortex, stimulator, Transcranial, tDCS
  10. [10] § 3. Results › 3.1. Experimental Results › 3.1.7. Sensitivity and Reliability Analyses ↔ tDCS_motorcontrol_03012026.py, lines 311–340 · score 0.51 · retest reliability, PRE POST, ICC, Coherence, 1–3 Hz, RMSE

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 · 1,011 lines · 42 KB · MIT · 10 matches

  1. #!/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3. """
  4. ================================================================================
  5. Reproducible Analysis Script
  7. Effects of Transcranial Direct Current Stimulation over the Left Sensorimotor
  8. Cortex on Bimanual Force Control: A Computational and Experimental Investigation
  10. Lima, V.M.S., Arthur, E.F., Gonzaga, R.R.D., Diniz, L.F.,
  11. Pedreiro, R.C.M., Pinto Neto, O.
  13. Correspondence: [email hidden]
  14. ================================================================================
  16. INSTRUCTIONS:
  17. 1. Place the data file 'data_tdcs_12212025.csv' and the raw force files
  18. (e.g., DAS_CONSTANTE_21_10_PRE.txt) in a folder of your choice.
  19. 2. Update DATA_PATH below to point to that folder.
  20. 3. Run this script: python reproducible_analysis.py
  21. 4. All figures and statistical output will be saved to DATA_PATH/results/
  23. Requirements:
  24. pip install numpy pandas scipy matplotlib
  26. ================================================================================
  27. """
  29. import numpy as np
  30. import pandas as pd
  31. import matplotlib.pyplot as plt
  32. from matplotlib.patches import FancyBboxPatch
  33. from matplotlib.gridspec import GridSpec
  34. from scipy.signal import butter, filtfilt, welch, coherence
  35. from scipy import stats
  36. from scipy.stats import f as f_dist
  37. import warnings
  38. import os
  39. import sys
  41. warnings.filterwarnings('ignore')
  43. # =============================================================================
  44. # USER CONFIGURATION — CHANGE THIS PATH
  45. # =============================================================================
  47. DATA_PATH = r"C:\Users\osmar\OneDrive\Documents\PESQUISAS\Vinicius"
  49. # Raw force files for Figure 1 (optional — schematic used if not found)
  50. # NOTE: On the authors' system these are in a different directory tree
  51. RAW_DATA_PATH = DATA_PATH
  52. EXAMPLE_PARTICIPANT = "DAS" # participant ID for Figure 1
  53. EXAMPLE_DATE = "21_10" # date string in filename (underscore, not dot)
  55. # =============================================================================
  56. # SETUP
  57. # =============================================================================
  59. DATA_FILE = os.path.join(DATA_PATH, "data_tdcs_12212025.csv")
  60. RESULTS_DIR = os.path.join(DATA_PATH, "results")
  61. os.makedirs(RESULTS_DIR, exist_ok=True)
  63. # Publication-quality figure defaults
  64. plt.rcParams.update({
  65. 'font.family': 'DejaVu Sans',
  66. 'font.size': 12,
  67. 'axes.titlesize': 14,
  68. 'axes.labelsize': 13,
  69. 'xtick.labelsize': 11,
  70. 'ytick.labelsize': 11,
  71. 'legend.fontsize': 10,
  72. 'legend.frameon': False,
  73. 'figure.dpi': 300,
  74. 'savefig.dpi': 300,
  75. 'savefig.bbox': 'tight',
  76. 'axes.spines.top': False,
  77. 'axes.spines.right': False,
  78. })
  81. # =============================================================================
  82. # 1. DATA LOADING AND PREPARATION
  83. # =============================================================================
  85. def load_data(filepath):
  86. """Load experimental data and compute derived metrics."""
  87. if not os.path.exists(filepath):
  88. print(f"\nERROR: Data file not found at:\n {filepath}")
  89. print(f"\nPlease update DATA_PATH at the top of this script.")
  90. sys.exit(1)
  92. df = pd.read_csv(filepath, encoding='utf-8-sig')
  94. # Standardize group label
  95. df['Group'] = df['Group'].replace('Placebo', 'Sham')
  97. # Compute force undershoot (%)
  98. df['Undershoot_pct'] = (
  99. 100 * (df['OF_F_alvo_mean'] - df['OF_Total_mean_raw'])
  100. / df['OF_F_alvo_mean']
  101. )
  103. print(f"Loaded {len(df)} observations from {df['participant'].nunique()} "
  104. f"participants ({', '.join(df['Group'].unique())})")
  105. return df
  108. def load_raw_force_file(filepath):
  109. """Load raw force data from Vernier .txt file."""
  110. with open(filepath, 'r', encoding='utf-8-sig') as f:
  111. lines = f.readlines()
  112. time_v, f1_v, f2_v, ft_v = [], [], [], []
  113. target = None
  114. for line in lines[7:]:
  115. line = line.strip().replace('\r', '')
  116. if not line:
  117. continue
  118. parts = line.split('\t')
  119. if len(parts) >= 4:
  120. try:
  121. time_v.append(float(parts[0].replace(',', '.')))
  122. f1_v.append(float(parts[1].replace(',', '.')))
  123. f2_v.append(float(parts[2].replace(',', '.')))
  124. ft_v.append(float(parts[3].replace(',', '.')))
  125. if len(parts) >= 5 and parts[4].strip():
  126. try:
  127. target = float(parts[4].replace(',', '.'))
  128. except ValueError:
  129. pass
  130. except ValueError:
  131. continue
  132. t_arr = np.array(time_v)
  133. fs = 1.0 / np.mean(np.diff(t_arr)) if len(t_arr) > 1 else 100.0
  134. return {'time': t_arr, 'force1': np.array(f1_v), 'force2': np.array(f2_v),
  135. 'total': np.array(ft_v), 'target': target, 'fs': fs}
  138. def load_participant_raw(participant, date_str):
  139. """Load PRE and POST raw force files for a participant."""
  140. result = {}
  141. for epoch, suffix in [('PRE', 'PRE'), ('POST', 'POS')]:
  142. fpath = os.path.join(RAW_DATA_PATH,
  143. f"{participant}_CONSTANTE_{date_str}_{suffix}.txt")
  144. if os.path.exists(fpath):
  145. result[epoch] = load_raw_force_file(fpath)
  146. return result if result else None
  149. # =============================================================================
  150. # 2. STATISTICAL ANALYSES
  151. # =============================================================================
  153. def compute_descriptives(df):
  154. """Compute cell means and SDs for all primary metrics."""
  155. metrics = {
  156. 'Undershoot_pct': 'Undershoot (%)',
  157. 'OF_Total_RMSE_raw': 'RMSE (N)',
  158. 'OF_Total_P_1_3Hz': 'Power 1-3 Hz (N²/Hz)',
  159. }
  160. rows = []
  161. for col, label in metrics.items():
  162. for grp in ['Sham', 'tDCS']:
  163. for ep in ['PRE', 'POS']:
  164. vals = df[(df['Group'] == grp) & (df['Epoch'] == ep)][col]
  165. rows.append({
  166. 'Metric': label, 'Group': grp, 'Epoch': ep,
  167. 'Mean': vals.mean(), 'SD': vals.std(),
  168. 'SEM': vals.std() / np.sqrt(len(vals)), 'N': len(vals)
  169. })
  170. return pd.DataFrame(rows)
  173. def run_interaction_tests(df):
  174. """
  175. Test Group × Epoch interaction using independent t-test on change scores.
  176. F(1,22) = t² for the between-group comparison of (POST − PRE) differences.
  177. """
  178. metrics = {
  179. 'Undershoot_pct': 'Undershoot (%)',
  180. 'OF_Total_RMSE_raw': 'RMSE (N)',
  181. 'OF_Total_P_1_3Hz': 'Power 1-3 Hz',
  182. }
  183. results = []
  184. for col, label in metrics.items():
  185. wide = df.pivot_table(index='participant', columns='Epoch',
  186. values=col, aggfunc='first')
  187. wide['Group'] = (df.drop_duplicates('participant')
  188. .set_index('participant')['Group'])
  189. wide['change'] = wide['POS'] - wide['PRE']
  191. tdcs = wide[wide['Group'] == 'tDCS']['change']
  192. sham = wide[wide['Group'] == 'Sham']['change']
  193. t_stat, p_val = stats.ttest_ind(tdcs, sham)
  195. results.append({
  196. 'Metric': label,
  197. 'tDCS_change': f"{tdcs.mean():.3f} ± {tdcs.std():.3f}",
  198. 'Sham_change': f"{sham.mean():.3f} ± {sham.std():.3f}",
  199. 'F(1,22)': t_stat ** 2,
  200. 'p_interaction': p_val,
  201. })
  202. return pd.DataFrame(results)
  205. def run_posthoc_paired(df):
  206. """Within-group paired t-tests (PRE vs POST) with Holm correction."""
  207. metrics = {
  208. 'Undershoot_pct': 'Undershoot (%)',
  209. 'OF_Total_RMSE_raw': 'RMSE (N)',
  210. 'OF_Total_P_1_3Hz': 'Power 1-3 Hz',
  211. }
  212. results = []
  213. for col, label in metrics.items():
  214. wide = df.pivot_table(index='participant', columns='Epoch',
  215. values=col, aggfunc='first')
  216. wide['Group'] = (df.drop_duplicates('participant')
  217. .set_index('participant')['Group'])
  218. for grp in ['tDCS', 'Sham']:
  219. sub = wide[wide['Group'] == grp]
  220. t_stat, p_val = stats.ttest_rel(sub['PRE'], sub['POS'])
  221. diff = sub['POS'] - sub['PRE']
  222. d = diff.mean() / diff.std()
  223. results.append({
  224. 'Metric': label, 'Group': grp,
  225. 'PRE_mean': sub['PRE'].mean(), 'PRE_sd': sub['PRE'].std(),
  226. 'POST_mean': sub['POS'].mean(), 'POST_sd': sub['POS'].std(),
  227. 't': t_stat, 'df': len(sub) - 1,
  228. 'p_uncorrected': p_val, 'Cohen_d': d,
  229. })
  231. res_df = pd.DataFrame(results)
  232. # Holm–Bonferroni correction across all 6 tests
  233. sorted_idx = res_df['p_uncorrected'].argsort().values
  234. n = len(res_df)
  235. p_holm = np.ones(n)
  236. for rank, idx in enumerate(sorted_idx):
  237. p_holm[idx] = min(res_df['p_uncorrected'].iloc[idx] * (n - rank), 1.0)
  238. for i in range(1, n):
  239. p_holm[sorted_idx[i]] = max(p_holm[sorted_idx[i]], p_holm[sorted_idx[i - 1]])
  240. res_df['p_holm'] = p_holm
  241. return res_df
  244. def run_correlations(df):
  245. """Pearson correlations between change scores, by group."""
  246. wide_u = df.pivot_table(index='participant', columns='Epoch',
  247. values='Undershoot_pct', aggfunc='first')
  248. wide_p = df.pivot_table(index='participant', columns='Epoch',
  249. values='OF_Total_P_1_3Hz', aggfunc='first')
  250. wide_r = df.pivot_table(index='participant', columns='Epoch',
  251. values='OF_Total_RMSE_raw', aggfunc='first')
  252. groups = (df.drop_duplicates('participant')
  253. .set_index('participant')['Group'])
  255. delta_u = wide_u['POS'] - wide_u['PRE']
  256. delta_p = wide_p['POS'] - wide_p['PRE']
  257. delta_r = wide_r['POS'] - wide_r['PRE']
  259. results = []
  260. for grp in ['tDCS', 'Sham']:
  261. mask = groups == grp
  262. r1, p1 = stats.pearsonr(delta_p[mask], delta_u[mask])
  263. r2, p2 = stats.pearsonr(delta_p[mask], delta_r[mask])
  264. results.append({'Group': grp, 'Comparison': 'ΔPower vs ΔUndershoot',
  265. 'r': r1, 'p': p1})
  266. results.append({'Group': grp, 'Comparison': 'ΔPower vs ΔRMSE',
  267. 'r': r2, 'p': p2})
  268. return pd.DataFrame(results)
  271. def run_ancova(df):
  272. """ANCOVA: POST ~ Group + PRE (baseline covariate)."""
  273. metrics = {
  274. 'Undershoot_pct': 'Undershoot (%)',
  275. 'OF_Total_RMSE_raw': 'RMSE (N)',
  276. 'OF_Total_P_1_3Hz': 'Power 1-3 Hz',
  277. }
  278. results = []
  279. for col, label in metrics.items():
  280. wide = df.pivot_table(index='participant', columns='Epoch',
  281. values=col, aggfunc='first')
  282. wide['Group'] = (df.drop_duplicates('participant')
  283. .set_index('participant')['Group'])
  284. wide['Group_code'] = (wide['Group'] == 'tDCS').astype(float)
  286. X = np.column_stack([
  287. np.ones(len(wide)),
  288. wide['Group_code'].values,
  289. wide['PRE'].values,
  290. ])
  291. y = wide['POS'].values
  292. # Manual OLS for minimal dependencies
  293. beta = np.linalg.lstsq(X, y, rcond=None)[0]
  294. y_hat = X @ beta
  295. residuals = y - y_hat
  296. n, k = X.shape
  297. mse = np.sum(residuals ** 2) / (n - k)
  298. se = np.sqrt(mse * np.diag(np.linalg.inv(X.T @ X)))
  299. t_vals = beta / se
  300. p_vals = 2 * stats.t.sf(np.abs(t_vals), df=n - k)
  302. results.append({
  303. 'Metric': label,
  304. 'Group_beta': beta[1], 'Group_SE': se[1],
  305. 'Group_t': t_vals[1], 'Group_p': p_vals[1],
  306. 'Baseline_beta': beta[2], 'Baseline_p': p_vals[2],
  307. })
  308. return pd.DataFrame(results)
  311. def compute_reliability(df):
  312. """Test-retest reliability (ICC) from Sham group PRE-POST."""
  313. metrics = {
  314. 'Undershoot_pct': 'Undershoot (%)',
  315. 'OF_Total_RMSE_raw': 'RMSE (N)',
  316. 'OF_Total_P_1_3Hz': 'Power 1-3 Hz',
  317. 'OF_Coh_1_3Hz': 'Coherence 1-3 Hz',
  318. }
  319. results = []
  320. sham = df[df['Group'] == 'Sham']
  321. for col, label in metrics.items():
  322. wide = sham.pivot_table(index='participant', columns='Epoch',
  323. values=col, aggfunc='first')
  324. pre, post = wide['PRE'].values, wide['POS'].values
  325. n = len(pre)
  326. # ICC(3,1) — two-way mixed, single measures, consistency
  327. grand_mean = np.mean(np.concatenate([pre, post]))
  328. ss_between = 2 * np.sum((np.mean([pre, post], axis=0) - grand_mean) ** 2)
  329. ss_within = np.sum((pre - np.mean([pre, post], axis=0)) ** 2 +
  330. (post - np.mean([pre, post], axis=0)) ** 2)
  331. ms_between = ss_between / (n - 1)
  332. ms_within = ss_within / n
  333. icc = (ms_between - ms_within) / (ms_between + ms_within)
  334. results.append({'Metric': label, 'ICC': icc, 'N': n})
  335. return pd.DataFrame(results)
  338. # =============================================================================
  339. # 3. COMPUTATIONAL MODEL
  340. # =============================================================================
  342. def butter_lowpass(x, fs, fc, order=4):
  343. """Apply a low-pass Butterworth filter."""
  344. b, a = butter(order, fc / (fs / 2), btype='low')
  345. return filtfilt(b, a, x)
  348. def simulate_bimanual_force(G_proprio=0.5, target_N=45.0, duration=40.0,
  349. fs=100.0, seed=None):
  350. """
  351. Closed-loop computational model of bimanual isometric force control.
  353. Matches the model described in Section 2.6 of the manuscript:
  354. - Vision ON (t < 20 s): visual feedback control with gain K_vis
  355. - Vision OFF (t >= 20 s): proprioceptive feedback with gain G_proprio,
  356. internal target drift (lambda_drift_base = 0.0045), and
  357. proprioceptive delay (delta)
  358. - Between-trial variability in lambda_drift and K_correction
  360. Parameters
  361. ----------
  362. G_proprio : float
  363. Proprioceptive feedback gain [0.2, 0.8]. Sham ≈ 0.25, tDCS ≈ 0.70.
  364. target_N : float
  365. Per-hand target force in Newtons (total target = 2 × target_N).
  366. duration : float
  367. Trial duration in seconds.
  368. fs : float
  369. Sampling rate in Hz.
  370. seed : int or None
  371. Random seed for reproducibility.
  373. Returns
  374. -------
  375. dict with keys: time, force_L, force_R, total_force, target,
  376. undershoot_pct, rmse, power_1_3Hz, freqs_psd, psd,
  377. freqs_coh, coh_spectrum, G_proprio
  378. """
  379. if seed is not None:
  380. np.random.seed(seed)
  382. dt = 1.0 / fs
  383. n_samples = int(duration * fs)
  384. t = np.arange(n_samples) / fs
  385. target_total = target_N * 2
  387. force_L = np.zeros(n_samples)
  388. force_R = np.zeros(n_samples)
  390. # ---- Fixed parameters (Table in Section 2.6) ----
  391. K_vis = 0.12 # Visual feedback gain
  392. K_base = 0.15 # Base correction gain
  393. w_common = 0.15 # Common drive weight
  394. delay_samples = int(0.1 * fs) # Proprioceptive delay δ = 100 ms
  395. sigma_proprio = 0.15 # Proprioceptive noise SD
  396. sigma_motor_von = 0.22 # Motor noise during Vision ON
  397. sigma_motor_voff = 0.25 # Motor noise during Vision OFF
  399. # ---- Between-trial variability (Section 2.6) ----
  400. lambda_drift_base = 0.0045 * (1 - 0.5 * G_proprio)
  401. lambda_drift = max(0.0005,
  402. lambda_drift_base + np.random.randn() * 0.0008)
  403. K_correction = max(0.02,
  404. K_base * G_proprio + np.random.randn() * 0.008)
  406. # ---- Shared signals ----
  407. common_drive = butter_lowpass(np.random.randn(n_samples) * 0.25, fs, 2)
  408. tremor_L = 0.06 * np.sin(2 * np.pi * 10.0 * t +
  409. np.random.rand() * 2 * np.pi)
  410. tremor_R = 0.06 * np.sin(2 * np.pi * 10.2 * t +
  411. np.random.rand() * 2 * np.pi)
  413. # ---- Simulation ----
  414. ramp_end = int(3 * fs)
  415. vision_off_start = int(20 * fs)
  416. internal_target = target_N
  418. for i in range(1, n_samples):
  419. if i < ramp_end:
  420. # Ramp-up phase
  421. frac = i / ramp_end
  422. force_L[i] = target_N * frac + np.random.randn() * 0.15
  423. force_R[i] = target_N * frac + np.random.randn() * 0.15
  425. elif i < vision_off_start:
  426. # Vision ON — Eq. (1) in manuscript
  427. err_L = target_N - force_L[i - 1]
  428. err_R = target_N - force_R[i - 1]
  429. force_L[i] = (force_L[i - 1]
  430. + K_vis * err_L
  431. + common_drive[i] * w_common
  432. + np.random.randn() * sigma_motor_von
  433. + tremor_L[i])
  434. force_R[i] = (force_R[i - 1]
  435. + K_vis * err_R
  436. + common_drive[i] * w_common
  437. + np.random.randn() * sigma_motor_von
  438. + tremor_R[i])
  440. else:
  441. # Vision OFF — Eqs. (2-4) in manuscript
  442. internal_target *= (1 - lambda_drift * dt)
  444. # Delayed proprioceptive estimate
  445. idx_delay = max(0, i - delay_samples)
  446. proprio_L = force_L[idx_delay] + np.random.randn() * sigma_proprio
  447. proprio_R = force_R[idx_delay] + np.random.randn() * sigma_proprio
  449. err_L = internal_target - proprio_L
  450. err_R = internal_target - proprio_R
  452. force_L[i] = (force_L[i - 1]
  453. + K_correction * err_L
  454. + common_drive[i] * w_common
  455. + np.random.randn() * sigma_motor_voff
  456. + tremor_L[i])
  457. force_R[i] = (force_R[i - 1]
  458. + K_correction * err_R
  459. + common_drive[i] * w_common
  460. + np.random.randn() * sigma_motor_voff
  461. + tremor_R[i])
  463. force_L[i] = np.clip(force_L[i], target_N * 0.7, target_N * 1.3)
  464. force_R[i] = np.clip(force_R[i], target_N * 0.7, target_N * 1.3)
  466. # ---- Compute metrics from Vision OFF analysis window (23–40 s) ----
  467. total_force = force_L + force_R
  468. of_start = int(23 * fs)
  469. of_end = int(40 * fs)
  470. of_force = total_force[of_start:of_end]
  472. undershoot = 100 * (target_total - np.mean(of_force)) / target_total
  473. rmse = np.sqrt(np.mean((of_force - target_total) ** 2))
  475. freqs_psd, psd = welch(of_force - np.mean(of_force), fs=fs,
  476. nperseg=min(1024, len(of_force)), noverlap=512)
  477. mask_1_3 = (freqs_psd >= 1) & (freqs_psd <= 3)
  478. _trapz = np.trapezoid if hasattr(np, 'trapezoid') else np.trapz
  479. power_1_3 = _trapz(psd[mask_1_3], freqs_psd[mask_1_3]) if mask_1_3.sum() > 1 else 0
  481. freqs_coh, coh = coherence(force_L[of_start:of_end],
  482. force_R[of_start:of_end],
  483. fs=fs, nperseg=min(512, len(of_force)))
  485. return {
  486. 'time': t, 'force_L': force_L, 'force_R': force_R,
  487. 'total_force': total_force, 'target': target_total,
  488. 'undershoot_pct': undershoot, 'rmse': rmse,
  489. 'power_1_3Hz': power_1_3,
  490. 'freqs_psd': freqs_psd, 'psd': psd,
  491. 'freqs_coh': freqs_coh, 'coh_spectrum': coh,
  492. 'G_proprio': G_proprio,
  493. }
  496. def run_dose_response(G_values, n_seeds=20, target_N=45.0):
  497. """Run simulations across a range of G_proprio values."""
  498. rows = []
  499. for G in G_values:
  500. for seed in range(n_seeds):
  501. sim = simulate_bimanual_force(G_proprio=G, target_N=target_N,
  502. seed=seed * 1000 + int(G * 10000))
  503. rows.append({
  504. 'G_proprio': G, 'seed': seed,
  505. 'undershoot_pct': sim['undershoot_pct'],
  506. 'rmse': sim['rmse'],
  507. 'power_1_3Hz': sim['power_1_3Hz'],
  508. })
  509. return pd.DataFrame(rows)
  512. # =============================================================================
  513. # 4. FIGURES
  514. # =============================================================================
  516. # Color palette
  517. C_SHAM = '#3182bd'
  518. C_TDCS = '#e31a1c'
  521. def fig1_experimental_design(raw_data, save_path):
  522. """Figure 1: Experimental design (montage + timeline + real force trace)."""
  523. print(" Creating Figure 1: Experimental Design...")
  525. fig = plt.figure(figsize=(12, 8))
  526. gs = GridSpec(2, 2, height_ratios=[1, 1.2], width_ratios=[1, 1.2])
  528. # Panel A: tDCS Montage schematic
  529. ax = fig.add_subplot(gs[0, 0])
  530. ax.set_xlim(-1.5, 1.5); ax.set_ylim(-1.5, 1.5)
  531. ax.set_aspect('equal'); ax.axis('off')
  532. ax.set_title('A. tDCS Electrode Montage', fontweight='bold')
  533. head = plt.Circle((0, 0), 1, fill=False, color='black', linewidth=2)
  534. ax.add_patch(head)
  535. ax.plot([0, 0], [1, 1.2], 'k-', lw=2)
  536. ax.text(0, 1.3, 'Nose', ha='center', fontsize=9)
  537. anode = plt.Circle((-0.25, 0.2), 0.22, color='red', alpha=0.7)
  538. ax.add_patch(anode)
  539. ax.text(-0.25, 0.2, '+', ha='center', va='center', fontsize=20,
  540. color='white', fontweight='bold')
  541. ax.text(-0.25, -0.12, 'Left of Cz\n(Anode)', ha='center', fontsize=9)
  542. cathode = plt.Circle((0.5, 0.7), 0.15, color='blue', alpha=0.7)
  543. ax.add_patch(cathode)
  544. ax.text(0.5, 0.7, '\u2013', ha='center', va='center', fontsize=18,
  545. color='white', fontweight='bold')
  546. ax.text(0.72, 0.85, 'Fp2\n(Cathode)', ha='center', fontsize=9)
  547. ax.plot(-1, 0, 'ko', ms=4); ax.text(-1.15, 0, 'T3', ha='right', fontsize=8)
  548. ax.plot(1, 0, 'ko', ms=4); ax.text(1.15, 0, 'T4', ha='left', fontsize=8)
  549. ax.text(-1.3, -0.8, "2 mA, 20 min\n35 cm\u00B2 electrodes\nSaline sponges",
  550. fontsize=9, bbox=dict(boxstyle='round', fc='lightgray', alpha=0.8),
  551. va='top')
  553. # Panel B: Task Timeline
  554. ax = fig.add_subplot(gs[0, 1])
  555. ax.set_xlim(-2, 42); ax.set_ylim(-0.5, 2); ax.axis('off')
  556. ax.set_title('B. Task Timeline', fontweight='bold')
  557. ax.add_patch(FancyBboxPatch((0, 0.5), 20, 1, boxstyle="round,pad=0.02",
  558. fc='lightblue', ec='blue', lw=2))
  559. ax.text(10, 1, 'Vision ON\n(0\u201320 s)', ha='center', va='center',
  560. fontsize=11, fontweight='bold')
  561. ax.add_patch(FancyBboxPatch((20, 0.5), 20, 1, boxstyle="round,pad=0.02",
  562. fc='lightyellow', ec='orange', lw=2))
  563. ax.text(30, 1, 'Vision OFF\n(20\u201340 s)', ha='center', va='center',
  564. fontsize=11, fontweight='bold')
  565. ax.add_patch(FancyBboxPatch((23, 0.1), 17, 0.35, boxstyle="round,pad=0.02",
  566. fc='lightgreen', ec='green', lw=1.5))
  567. ax.text(31.5, 0.27, 'Analysis Window (23\u201340 s)', ha='center',
  568. va='center', fontsize=9)
  569. ax.annotate('', xy=(42, 0), xytext=(-2, 0),
  570. arrowprops=dict(arrowstyle='->', color='black', lw=1.5))
  571. ax.text(42, -0.15, 'Time (s)', ha='right', fontsize=10)
  572. for t_mark in [0, 3, 20, 23, 40]:
  573. ax.plot([t_mark, t_mark], [-0.1, 0.1], 'k-', lw=1)
  574. ax.text(t_mark, -0.25, str(t_mark), ha='center', fontsize=9)
  576. # Panel C: Real (or schematic) force trace
  577. ax = fig.add_subplot(gs[1, :])
  578. if raw_data is not None and 'POST' in raw_data:
  579. data = raw_data['POST']
  580. time_arr = data['time']
  581. total = data['total']
  582. target = data['target']
  583. fs_raw = data['fs']
  584. total_sm = butter_lowpass(total, fs_raw, 5)
  585. mask40 = time_arr <= 40
  586. time_arr, total_sm = time_arr[mask40], total_sm[mask40]
  587. ax.set_title('C. Bimanual Force-Matching Task (Representative Trial)',
  588. fontweight='bold')
  589. ax.axvspan(0, 20, alpha=0.15, color='blue', label='Vision ON')
  590. ax.axvspan(20, 40, alpha=0.15, color='orange', label='Vision OFF')
  591. ax.plot(time_arr, total_sm, 'b-', lw=1, label='Total Force')
  592. ax.axhline(target, color='green', ls='--', lw=2,
  593. label=f'Target ({target:.1f} N)')
  594. of_mask = (time_arr >= 23) & (time_arr <= 40)
  595. if of_mask.sum() > 0:
  596. mean_of = np.mean(total_sm[of_mask])
  597. u_pct = 100 * (target - mean_of) / target
  598. ax.annotate('', xy=(32, mean_of-10), xytext=(32, target-10),
  599. arrowprops=dict(arrowstyle='<->', color='red', lw=2.5))
  600. ax.text(33.5, (mean_of + target) / 2-10,
  601. f'Undershoot\n{u_pct:.1f}%', fontsize=14,
  602. color='red', va='center', fontweight='bold')
  603. ax.axhline(mean_of, xmin=0.575, xmax=1.0, color='red', ls=':',
  604. lw=1.5, alpha=0.7)
  605. ax.text(10, target + 4, 'Visual feedback\navailable',
  606. ha='center', fontsize=10)
  607. ax.text(30, target + 4, 'Proprioceptive\nfeedback only',
  608. ha='center', fontsize=10)
  609. ax.set_ylim(np.min(total_sm) - 5, target + 40)
  610. else:
  611. ax.set_title('C. Bimanual Force-Matching Task (Schematic)',
  612. fontweight='bold')
  613. fs_s = 100; t_s = np.arange(0, 40, 1 / fs_s); tgt = 45.0
  614. force = np.zeros_like(t_s)
  615. force[t_s < 3] = tgt * (t_s[t_s < 3] / 3)
  616. m_on = (t_s >= 3) & (t_s < 20)
  617. force[m_on] = tgt + np.random.randn(m_on.sum()) * 0.8
  618. force[m_on] = butter_lowpass(force[m_on], fs_s, 3)
  619. force[m_on] += tgt - np.mean(force[m_on])
  620. m_off = t_s >= 20
  621. force[m_off] = tgt - 0.15 * (t_s[m_off] - 20) + np.random.randn(m_off.sum()) * 1
  622. force[m_off] = butter_lowpass(force[m_off], fs_s, 3)
  623. force = butter_lowpass(force, fs_s, 5)
  624. ax.axvspan(0, 20, alpha=0.15, color='blue', label='Vision ON')
  625. ax.axvspan(20, 40, alpha=0.15, color='orange', label='Vision OFF')
  626. ax.plot(t_s, force, 'b-', lw=1.5, label='Total Force')
  627. ax.axhline(tgt, color='green', ls='--', lw=2, label='Target (30% MVC)')
  628. m_of = np.mean(force[t_s >= 23])
  629. ax.annotate('', xy=(35, m_of), xytext=(35, tgt),
  630. arrowprops=dict(arrowstyle='<->', color='red', lw=2))
  631. ax.text(36.5, (m_of + tgt) / 2, 'Undershoot', fontsize=14,
  632. color='red', va='center')
  633. ax.text(10, tgt + 3, 'Visual feedback\navailable', ha='center',
  634. fontsize=10)
  635. ax.text(30, tgt + 3, 'Proprioceptive\nfeedback only', ha='center',
  636. fontsize=10)
  637. ax.set_ylim(38, 50)
  638. ax.set_xlabel('Time (s)')
  639. ax.set_ylabel('Force (N)')
  640. ax.set_xlim(0, 42)
  641. ax.legend(loc='lower right')
  642. plt.tight_layout()
  643. fpath = os.path.join(save_path, 'Figure1_Experimental_Design.png')
  644. plt.savefig(fpath, facecolor='white')
  645. plt.savefig(fpath.replace('.png', '.svg'), format='svg')
  646. print(f" Saved: {fpath}")
  647. plt.close()
  650. def fig2_experimental_results(df, save_path):
  651. """Figure 2: Experimental results (Group × Epoch bar plots)."""
  652. print(" Creating Figure 2: Experimental Results...")
  654. fig, axes = plt.subplots(2, 2, figsize=(10, 8))
  655. epochs = ['PRE', 'POS']
  656. x = np.array([0, 1])
  657. w = 0.35
  659. metric_specs = [
  660. ('Undershoot_pct', 'Undershoot (%)', 'A. Force Undershoot', axes[0, 0]),
  661. ('OF_Total_RMSE_raw', 'RMSE (N)', 'B. Root Mean Square Error', axes[0, 1]),
  662. ('OF_Total_P_1_3Hz', 'Power (N²/Hz)', 'C. Spectral Power (1–3 Hz)', axes[1, 0]),
  663. ]
  665. for col, ylabel, title, ax in metric_specs:
  666. for gi, (grp, color) in enumerate(
  667. [('Sham', C_SHAM), ('tDCS', C_TDCS)]):
  668. means, sems = [], []
  669. for ep in epochs:
  670. vals = df[(df['Group'] == grp) & (df['Epoch'] == ep)][col]
  671. means.append(vals.mean())
  672. sems.append(vals.std() / np.sqrt(len(vals)))
  673. offset = -w / 2 if gi == 0 else w / 2
  674. ax.bar(x + offset, means, w, yerr=sems, capsize=5,
  675. color=color, alpha=0.6, label=grp, edgecolor='black',
  676. error_kw={'elinewidth': 1.5, 'capthick': 1.5})
  678. # Add within-tDCS significance bracket
  679. posthoc = run_posthoc_paired(df)
  680. row = posthoc[(posthoc['Metric'] == ylabel.split(' (')[0] +
  681. (' (%)' if '%' in ylabel else ' (N)' if 'N' in ylabel else ''))
  682. & (posthoc['Group'] == 'tDCS')]
  683. # Simpler approach: just mark tDCS POST bar
  684. tdcs_post = df[(df['Group'] == 'tDCS') & (df['Epoch'] == 'POS')][col]
  685. tdcs_pre = df[(df['Group'] == 'tDCS') & (df['Epoch'] == 'PRE')][col]
  686. t_val, p_val = stats.ttest_rel(tdcs_pre, tdcs_post)
  687. if p_val < 0.05:
  688. ymax = max(tdcs_pre.mean(), tdcs_post.mean()) + \
  689. max(tdcs_pre.std(), tdcs_post.std()) / np.sqrt(12)
  690. ax.plot([0 + w/2, 1 + w/2], [ymax * 1.08, ymax * 1.08],
  691. 'k-', linewidth=1.5)
  692. stars = '***' if p_val < 0.001 else '**' if p_val < 0.01 else '*'
  693. ax.text(0.5 + w/2, ymax * 1.10, stars, ha='center', fontsize=14)
  695. ax.set_ylabel(ylabel)
  696. ax.set_title(title, fontweight='bold')
  697. ax.set_xticks(x)
  698. ax.set_xticklabels(['PRE', 'POST'])
  699. ax.legend()
  701. # Panel D: Coherence
  702. ax = axes[1, 1]
  703. bands = ['0–1 Hz', '1–3 Hz', '3–7 Hz', '7–12 Hz']
  704. coh_cols = ['OF_Coh_0_1Hz', 'OF_Coh_1_3Hz', 'OF_Coh_3_7Hz', 'OF_Coh_7_12Hz']
  705. x_coh = np.arange(len(bands))
  706. bw = 0.2
  708. for ei, (ep, alpha) in enumerate(zip(['PRE', 'POS'], [0.4, 0.8])):
  709. for gi, (grp, color) in enumerate([('Sham', C_SHAM), ('tDCS', C_TDCS)]):
  710. means = [df[(df['Group'] == grp) & (df['Epoch'] == ep)][c].mean()
  711. for c in coh_cols]
  712. sems = [df[(df['Group'] == grp) & (df['Epoch'] == ep)][c].std()
  713. / np.sqrt(12) for c in coh_cols]
  714. offset = (-1.5 + ei + gi * 2) * bw
  715. label = f'{grp} {ep.replace("POS", "POST")}'
  716. ax.bar(x_coh + offset, means, bw, yerr=sems, capsize=3,
  717. color=color, alpha=alpha, label=label, edgecolor='black',
  718. error_kw={'elinewidth': 1, 'capthick': 1})
  720. ax.set_ylabel('Coherence')
  721. ax.set_title('D. Inter-hand Coherence', fontweight='bold')
  722. ax.set_xticks(x_coh)
  723. ax.set_xticklabels(bands)
  724. ax.legend(ncol=2, fontsize=9)
  725. ax.set_ylim(0, 0.7)
  727. plt.tight_layout()
  728. fpath = os.path.join(save_path, 'Figure2_Experimental_Results.png')
  729. plt.savefig(fpath, facecolor='white')
  730. plt.savefig(fpath.replace('.png', '.svg'), format='svg')
  731. print(f" Saved: {fpath}")
  732. plt.close()
  735. def fig3_correlations(df, save_path):
  736. """Figure 3: Individual-level correlations (ΔPower vs ΔUndershoot/ΔRMSE)."""
  737. print(" Creating Figure 3: Correlations...")
  739. wide_u = df.pivot_table(index='participant', columns='Epoch',
  740. values='Undershoot_pct', aggfunc='first')
  741. wide_p = df.pivot_table(index='participant', columns='Epoch',
  742. values='OF_Total_P_1_3Hz', aggfunc='first')
  743. wide_r = df.pivot_table(index='participant', columns='Epoch',
  744. values='OF_Total_RMSE_raw', aggfunc='first')
  745. groups = df.drop_duplicates('participant').set_index('participant')['Group']
  747. delta_u = wide_u['POS'] - wide_u['PRE']
  748. delta_p = wide_p['POS'] - wide_p['PRE']
  749. delta_r = wide_r['POS'] - wide_r['PRE']
  751. fig, axes = plt.subplots(1, 2, figsize=(12, 5))
  752. colors = {'Sham': C_SHAM, 'tDCS': C_TDCS}
  754. for ax, (delta_y, ylabel, title) in zip(axes, [
  755. (delta_u, 'Δ Undershoot (%)', 'A. Δ Power vs Δ Undershoot'),
  756. (delta_r, 'Δ RMSE (N)', 'B. Δ Power vs Δ RMSE'),
  757. ]):
  758. for grp in ['Sham', 'tDCS']:
  759. mask = groups == grp
  760. xv, yv = delta_p[mask], delta_y[mask]
  761. ax.scatter(xv, yv, c=colors[grp], s=100, edgecolors='k',
  762. label=grp, zorder=5)
  763. if len(xv) > 2:
  764. slope, intercept, r, p, _ = stats.linregress(xv, yv)
  765. xline = np.array([xv.min(), xv.max()])
  766. ls = '-' if p < 0.05 else '--'
  767. ax.plot(xline, intercept + slope * xline,
  768. color=colors[grp], linestyle=ls, alpha=0.7, lw=2)
  769. sig = '*' if p < 0.05 else ''
  770. ax.text(0.98, 0.95 if grp == 'tDCS' else 0.85,
  771. f'{grp}: r = {r:.2f}, p = {p:.3f}{sig}',
  772. transform=ax.transAxes, ha='right', fontsize=10,
  773. color=colors[grp])
  775. ax.axhline(0, color='gray', ls=':', alpha=0.5)
  776. ax.axvline(0, color='gray', ls=':', alpha=0.5)
  777. ax.set_xlabel('Δ Power 1–3 Hz (N²/Hz)')
  778. ax.set_ylabel(ylabel)
  779. ax.set_title(title, fontweight='bold')
  780. ax.legend()
  782. plt.tight_layout()
  783. fpath = os.path.join(save_path, 'Figure3_Correlations.png')
  784. plt.savefig(fpath, facecolor='white')
  785. plt.savefig(fpath.replace('.png', '.svg'), format='svg')
  786. print(f" Saved: {fpath}")
  787. plt.close()
  790. def fig4_model_results(df, save_path):
  791. """Figure 4: Computational model simulated time series (2 panels)."""
  792. print(" Creating Figure 4: Model Results...")
  794. sim_sham = simulate_bimanual_force(G_proprio=0.25, seed=42)
  795. sim_tdcs = simulate_bimanual_force(G_proprio=0.70, seed=42)
  796. t = sim_sham['time']
  797. target = sim_sham['target']
  799. fig, axes = plt.subplots(1, 2, figsize=(13, 5))
  801. # Panel A: Full time series
  802. ax = axes[0]
  803. ax.axvspan(0, 20, alpha=0.12, color='blue', label='Vision ON')
  804. ax.axvspan(20, 40, alpha=0.12, color='orange', label='Vision OFF')
  805. ax.plot(t, sim_sham['total_force'], color=C_SHAM, lw=0.8, alpha=0.8,
  806. label='Sham (G = 0.25)')
  807. ax.plot(t, sim_tdcs['total_force'], color=C_TDCS, lw=0.8, alpha=0.8,
  808. label='tDCS (G = 0.70)')
  809. ax.axhline(target, color='green', ls='--', lw=2, label='Target')
  810. ax.set_xlabel('Time (s)')
  811. ax.set_ylabel('Total Force (N)')
  812. ax.set_title('A. Simulated Force Time Series', fontweight='bold')
  813. ax.legend(loc='lower right', fontsize=9)
  814. ax.set_xlim(0, 40)
  815. ax.set_ylim(target * 0.87, target * 1.07)
  817. # Panel B: Vision OFF detail
  818. ax = axes[1]
  819. of_mask = (t >= 23) & (t <= 40)
  820. ax.plot(t[of_mask], sim_sham['total_force'][of_mask], color=C_SHAM,
  821. lw=1.2, label='Sham')
  822. ax.plot(t[of_mask], sim_tdcs['total_force'][of_mask], color=C_TDCS,
  823. lw=1.2, label='tDCS')
  824. ax.axhline(target, color='green', ls='--', lw=2)
  826. m_s = np.mean(sim_sham['total_force'][of_mask])
  827. m_t = np.mean(sim_tdcs['total_force'][of_mask])
  828. u_s = 100 * (target - m_s) / target
  829. u_t = 100 * (target - m_t) / target
  830. ax.axhline(m_s, color=C_SHAM, ls=':', alpha=0.7,
  831. label=f'Sham mean: {m_s:.1f} N ({u_s:.1f}%)')
  832. ax.axhline(m_t, color=C_TDCS, ls=':', alpha=0.7,
  833. label=f'tDCS mean: {m_t:.1f} N ({u_t:.1f}%)')
  834. ax.set_xlabel('Time (s)')
  835. ax.set_ylabel('Total Force (N)')
  836. ax.set_title('B. Vision OFF Epoch Detail', fontweight='bold')
  837. ax.legend(fontsize=9)
  838. ax.set_ylim(target * 0.90, target * 1.05)
  840. plt.tight_layout()
  841. fpath = os.path.join(save_path, 'Figure4_Model_Results.png')
  842. plt.savefig(fpath, facecolor='white')
  843. plt.savefig(fpath.replace('.png', '.svg'), format='svg')
  844. print(f" Saved: {fpath}")
  845. plt.close()
  848. def fig5_model_dose_response(save_path):
  849. """Figure 5: Computational model dose-response."""
  850. print(" Creating Figure 5: Model Dose-Response...")
  852. G_values = [0.20, 0.25, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80]
  853. sweep = run_dose_response(G_values, n_seeds=20)
  854. grouped = sweep.groupby('G_proprio')
  855. G_vals = sorted(sweep['G_proprio'].unique())
  857. fig, axes = plt.subplots(1, 3, figsize=(15, 5))
  859. for ax, (metric, ylabel, title, scale) in zip(axes, [
  860. ('undershoot_pct', 'Undershoot (%)', 'A. Undershoot vs G_proprio', 1),
  861. ('power_1_3Hz', 'Power 1–3 Hz (×10 N²/Hz)',
  862. 'B. Corrective Power vs G_proprio', 10),
  863. ]):
  864. means = [grouped.get_group(g)[metric].mean() * scale for g in G_vals]
  865. stds = [grouped.get_group(g)[metric].std() * scale for g in G_vals]
  866. ax.errorbar(G_vals, means, yerr=stds, fmt='o-', capsize=5,
  867. color='gray', lw=2, ms=8, alpha=0.8, capthick=1.5,
  868. elinewidth=1.5, label='Model sweep')
  870. # Highlight Sham and tDCS
  871. for G, color, marker, label in [
  872. (0.25, C_SHAM, 'o', 'Sham (G = 0.25)'),
  873. (0.70, C_TDCS, 's', 'tDCS (G = 0.70)'),
  874. ]:
  875. g = grouped.get_group(G)[metric]
  876. ax.errorbar([G], [g.mean() * scale], yerr=[g.std() * scale],
  877. fmt=marker, capsize=8, color=color, ms=14,
  878. markeredgewidth=2, markeredgecolor='black',
  879. capthick=2.5, elinewidth=2.5, label=label)
  881. ax.set_xlabel('G_proprio (Proprioceptive Gain)')
  882. ax.set_ylabel(ylabel)
  883. ax.set_title(title, fontweight='bold')
  884. ax.legend()
  885. ax.grid(True, alpha=0.3)
  887. # Panel C: Combined dose-response
  888. ax = axes[2]
  889. ax2 = ax.twinx()
  890. means_u = [grouped.get_group(g)['undershoot_pct'].mean() for g in G_vals]
  891. means_p = [grouped.get_group(g)['power_1_3Hz'].mean() * 10 for g in G_vals]
  893. l1, = ax.plot(G_vals, means_u, 'o-', color=C_SHAM, lw=2, ms=7,
  894. label='Undershoot (%)')
  895. l2, = ax2.plot(G_vals, means_p, 's-', color=C_TDCS, lw=2, ms=7,
  896. label='Power 1–3 Hz (×10)')
  898. for G, label_txt in [(0.25, 'Sham'), (0.70, 'tDCS')]:
  899. g = grouped.get_group(G)
  900. ax.plot(G, g['undershoot_pct'].mean(), 'o', color=C_SHAM, ms=14,
  901. markeredgecolor='black', markeredgewidth=2, zorder=10)
  902. ax2.plot(G, g['power_1_3Hz'].mean() * 10, 's', color=C_TDCS, ms=14,
  903. markeredgecolor='black', markeredgewidth=2, zorder=10)
  904. ax.annotate(label_txt, xy=(G, g['undershoot_pct'].mean()),
  905. xytext=(0, 12), textcoords='offset points',
  906. ha='center', fontweight='bold', fontsize=11)
  908. ax.set_xlabel('G_proprio (Proprioceptive Gain)')
  909. ax.set_ylabel('Undershoot (%)', color=C_SHAM)
  910. ax2.set_ylabel('Power 1–3 Hz (×10)', color=C_TDCS)
  911. ax.set_title('C. Dose-Response Relationship', fontweight='bold')
  912. ax.tick_params(axis='y', labelcolor=C_SHAM)
  913. ax2.tick_params(axis='y', labelcolor=C_TDCS)
  914. ax.legend([l1, l2], [l1.get_label(), l2.get_label()], loc='center right')
  915. ax.grid(True, alpha=0.3)
  917. plt.tight_layout()
  918. fpath = os.path.join(save_path, 'Figure5_Model_DoseResponse.png')
  919. plt.savefig(fpath, facecolor='white')
  920. plt.savefig(fpath.replace('.png', '.svg'), format='svg')
  921. print(f" Saved: {fpath}")
  922. plt.close()
  926. # =============================================================================
  927. # 5. MAIN EXECUTION
  928. # =============================================================================
  930. def main():
  931. print("=" * 75)
  932. print("REPRODUCIBLE ANALYSIS — tDCS Bimanual Force Control")
  933. print("=" * 75)
  935. # ---- Load data ----
  936. df = load_data(DATA_FILE)
  937. print()
  939. # ---- Statistical analyses ----
  940. print("Running statistical analyses...")
  942. desc = compute_descriptives(df)
  943. print("\n--- Descriptive Statistics ---")
  944. print(desc.to_string(index=False))
  946. interactions = run_interaction_tests(df)
  947. print("\n--- Group × Epoch Interactions ---")
  948. print(interactions.to_string(index=False))
  950. posthoc = run_posthoc_paired(df)
  951. print("\n--- Post-hoc Paired t-tests (Holm-corrected) ---")
  952. cols = ['Metric', 'Group', 'PRE_mean', 'POST_mean', 't', 'df',
  953. 'p_uncorrected', 'p_holm', 'Cohen_d']
  954. print(posthoc[cols].to_string(index=False))
  956. corr = run_correlations(df)
  957. print("\n--- Correlations (change scores) ---")
  958. print(corr.to_string(index=False))
  960. ancova = run_ancova(df)
  961. print("\n--- ANCOVA (POST ~ Group + Baseline) ---")
  962. print(ancova.to_string(index=False))
  964. reliability = compute_reliability(df)
  965. print("\n--- Test-Retest Reliability (Sham ICC) ---")
  966. print(reliability.to_string(index=False))
  968. # ---- Save statistical tables ----
  969. desc.to_csv(os.path.join(RESULTS_DIR, 'descriptives.csv'), index=False)
  970. interactions.to_csv(os.path.join(RESULTS_DIR, 'interactions.csv'), index=False)
  971. posthoc.to_csv(os.path.join(RESULTS_DIR, 'posthoc.csv'), index=False)
  972. corr.to_csv(os.path.join(RESULTS_DIR, 'correlations.csv'), index=False)
  973. ancova.to_csv(os.path.join(RESULTS_DIR, 'ancova.csv'), index=False)
  974. reliability.to_csv(os.path.join(RESULTS_DIR, 'reliability.csv'), index=False)
  975. print(f"\nStatistical tables saved to: {RESULTS_DIR}")
  977. # ---- Generate figures ----
  978. print("\nGenerating figures...")
  980. # Figure 1: Experimental Design (uses raw force files if available)
  981. raw_data = None
  982. try:
  983. raw_data = load_participant_raw(EXAMPLE_PARTICIPANT,
  984. EXAMPLE_DATE)
  985. if raw_data:
  986. print(f" Loaded raw data for {EXAMPLE_PARTICIPANT} "
  987. f"({len(raw_data)} epochs)")
  988. except Exception as e:
  989. print(f" Raw data not found ({e}); using schematic for Figure 1")
  990. fig1_experimental_design(raw_data, RESULTS_DIR)
  992. # Figure 3: Experimental Results
  993. fig2_experimental_results(df, RESULTS_DIR)
  995. # Figure 4: Correlations
  996. fig3_correlations(df, RESULTS_DIR)
  998. # Figure 5: Computational Model Results
  999. fig4_model_results(df, RESULTS_DIR)
  1001. # Figure 6: Model Dose-Response
  1002. fig5_model_dose_response(RESULTS_DIR)
  1004. print("\n" + "=" * 75)
  1005. print("ANALYSIS COMPLETE")
  1006. print(f"All outputs saved to: {RESULTS_DIR}")
  1007. print("=" * 75)
  1010. if __name__ == "__main__":
  1011. main()

tDCS_motorcontrol_03012026.py at commit 5053987, under MIT · at the source

Overview

Authors: Vinicius de Moura Silva Lima1, Eduarda Faria Arthur1,2, Rafaela Rodrigues Dousseau Gonzaga1,2, Luan Faria Diniz1, Rodrigo Cunha de Mello Pedreiro3,4, Osmar Pinto Neto1,2,5
ORCID iDs: Rodrigo Cunha de Mello Pedreiro, Osmar Pinto Neto
  1. Biomedical Engineering Department, Anhembi Morumbi University, São José dos Campos 12247-016, SP, Brazil; (V.d.M.S.L.); (E.F.A.); (R.R.D.G.); (L.F.D.)
  2. Arena235Research Lab., São José dos Campos 12246-876, SP, Brazil
  3. School of Physical Education and Sports, Federal University of Rio de Janeiro, Rio de Janeiro 21941-901, RJ, Brazil
  4. Physical Education Department, Estácio de Sá University, Teresópolis 25963-150, RJ, Brazil
  5. Kinesiology Department, California State University San Marcos, San Marcos, CA 92096, USA
Institutions: Anhembi Morumbi University (Brazil); Universidade Federal do Rio de Janeiro (Brazil); Estácio (Brazil) (Brazil); California State University, San Marcos (United States)
Journal: Bioengineering (Basel, Switzerland), volume 13, issue 5, article 502
Dates: received 2 March 2026; accepted 23 April 2026; published online 26 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3390/bioengineering13050502 · PMID 42194259 · PMCID PMC13203247 · OpenAlex W7159775081
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: other (modality), human (organism), systems (subfield)
Methods: Spectral & time-frequency, Statistics, Connectivity, Preprocessing
Keywords: transcranial direct current stimulation, primary motor cortex, lateral premotor cortex, supplementary motor area, bimanual coordination, proprioception, force control, computational modeling
Topic: Transcranial Magnetic Stimulation Studies (Neurology, Neuroscience), according to OpenAlex
Funding: Anima Institute (#12/2025)
Citations: not cited yet (Europe PMC); 30 references in the paper

Abstract

Transcranial direct current stimulation (tDCS) over motor–premotor regions may modulate motor performance, though underlying mechanisms remain unclear. Twenty-four athletes (9 females, 15 males) were randomly assigned to receive anodal tDCS (2 mA, 20 min) over the left sensorimotor cortex (n = 12) or sham stimulation (n = 12). Participants performed a bimanual isometric force-matching task at 30% maximal voluntary contraction, with visual feedback initially provided and then removed. Force undershoot, root mean square error (RMSE), spectral power (1–3 Hz), and inter-hand coherence were analyzed. A computational model was developed to test whether enhanced proprioceptive feedback processing could account for observed effects. Following tDCS, force undershoot decreased significantly (p = 0.002, d = −1.15) and RMSE improved (p = 0.010, d = −0.91). Spectral power in the 1–3 Hz band increased (p = 0.012, d = 0.87), suggesting enhanced corrective oscillations. These within-group changes were absent in the sham group (all p > 0.20), although Group × Epoch interactions did not reach significance (all p > 0.05), likely due to limited statistical power. Inter-hand coherence remained unchanged. The computational model demonstrated that enhanced proprioceptive feedback gain qualitatively reproduces the observed behavioral pattern. Anodal tDCS over the left sensorimotor/premotor region may enhance bimanual force control under conditions requiring proprioceptive feedback. Replication with larger samples is needed to confirm between-group specificity.

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

Zenodo 19745973

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

osmar235/tdcs_left_sensorimotor_cortex_bimanual_force_control

License: MIT
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 5053987cc6098c994c615ca0555a0dd364601fdc, 1 March 2026
Languages: Python (1)
Size: 8 files, 1 script
Software Heritage: not archived
Found in: the Zenodo archive record
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (1 file), NumPy (1 file), pandas (1 file), SciPy (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
3 files
At the source: github.com/osmar235/tDCS_Left_Sensorimotor_Cortex_Bimanual_Force_Control/tree/v.01

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;
  • 10 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

Anonymized experimental data and computational model code are available at Zenodo (DOI: 10.5281/zenodo.19745973) and can also be obtained upon request from the corresponding author.

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

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 8 keywords, 1 funder, 30 references.

Cite

This paper

Lima, V. d. M. S., Arthur, E. F., Gonzaga, R. R. D., Diniz, L. F., Pedreiro, R. C. d. M., & Pinto Neto, O. (2026). Effects of Transcranial Direct Current Stimulation over the Left Sensorimotor Cortex on Bimanual Force Control: A Computational and Experimental Investigation. Bioengineering (Basel, Switzerland), 13(5), 502. https://doi.org/10.3390/bioengineering13050502

BibTeX

@article{lima2026effects,
author = {Lima, Vinicius de Moura Silva and Arthur, Eduarda Faria and Gonzaga, Rafaela Rodrigues Dousseau and Diniz, Luan Faria and Pedreiro, Rodrigo Cunha de Mello and Pinto Neto, Osmar},
title = {{Effects of Transcranial Direct Current Stimulation over the Left Sensorimotor Cortex on Bimanual Force Control: A Computational and Experimental Investigation}},
journal = {Bioengineering (Basel, Switzerland)},
year = {2026},
month = apr,
volume = {13},
number = {5},
pages = {502},
publisher = {Multidisciplinary Digital Publishing Institute (MDPI)},
issn = {2306-5354},
doi = {10.3390/bioengineering13050502},
url = {https://doi.org/10.3390/bioengineering13050502},
pmid = {42194259},
pmcid = {PMC13203247}
}

RIS

TY - JOUR
AU - Lima, Vinicius de Moura Silva
AU - Arthur, Eduarda Faria
AU - Gonzaga, Rafaela Rodrigues Dousseau
AU - Diniz, Luan Faria
AU - Pedreiro, Rodrigo Cunha de Mello
AU - Pinto Neto, Osmar
TI - Effects of Transcranial Direct Current Stimulation over the Left Sensorimotor Cortex on Bimanual Force Control: A Computational and Experimental Investigation
T2 - Bioengineering (Basel, Switzerland)
J2 - Bioengineering (Basel)
PY - 2026
DA - 2026/04/26
VL - 13
IS - 5
SP - 502
SN - 2306-5354
PB - Multidisciplinary Digital Publishing Institute (MDPI)
DO - 10.3390/bioengineering13050502
UR - https://doi.org/10.3390/bioengineering13050502
LA - en
ER -

CSL-JSON

{
"id": "10.3390/bioengineering13050502",
"type": "article-journal",
"title": "Effects of Transcranial Direct Current Stimulation over the Left Sensorimotor Cortex on Bimanual Force Control: A Computational and Experimental Investigation",
"container-title": "Bioengineering (Basel, Switzerland)",
"author": [
{
"family": "Lima",
"given": "Vinicius de Moura Silva"
},
{
"family": "Arthur",
"given": "Eduarda Faria"
},
{
"family": "Gonzaga",
"given": "Rafaela Rodrigues Dousseau"
},
{
"family": "Diniz",
"given": "Luan Faria"
},
{
"family": "Pedreiro",
"given": "Rodrigo Cunha de Mello"
},
{
"family": "Pinto Neto",
"given": "Osmar"
}
],
"container-title-short": "Bioengineering (Basel)",
"volume": "13",
"issue": "5",
"page": "502",
"DOI": "10.3390/bioengineering13050502",
"PMID": "42194259",
"PMCID": "PMC13203247",
"ISSN": "2306-5354",
"publisher": "Multidisciplinary Digital Publishing Institute (MDPI)",
"URL": "https://doi.org/10.3390/bioengineering13050502",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
26
]
]
}
}

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.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: Matplotlib, NumPy, other, 3 references
[2] 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, Matplotlib, 1 other tool, other, systems, 1 reference
[3] doi:10.1038/s41598-026-45549-3
Effects of single-session anodal transcranial direct current stimulation (tDCS) on cognitive and motor performance in athletes and healthy adults: a systematic review and meta-analysis.
Journal: Scientific reports
In common: other, 3 references
[4] doi:10.1038/s41467-026-75799-8 [code]
Behaviourally driven closed-loop beta-tACS enhances beta activity and motor behaviour.
Journal: Nature communications
In common: other, systems, 3 references
[5] doi:10.1038/s41467-026-72917-4 [code]
OFC-induced network modularity improves positive symptoms and attentional alertness in schizophrenia: a combined rTMS-fMRI study.
Journal: Nature communications
In common: pandas, SciPy, Matplotlib, 1 other tool, other, 1 reference
[6] doi:10.1162/imag.a.1222 [code]
Network-based near-scalp personalized brain stimulation targets.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: SciPy, Matplotlib, NumPy, other, systems, 1 reference
[7] doi:10.3390/s26103131 [code]
Neuromagnetism "On the Cheap": Evaluating a Combined Cylindrical Shield and Partial-Coverage OPM-MEG System for Detecting Sensorimotor Responses in Humans.
Journal: Sensors (Basel, Switzerland)
In common: pandas, SciPy, Matplotlib, 1 other tool, systems, 1 reference
[8] doi:10.1162/imag.a.1256 [code]
Gamer in the scanner: Event-related analysis of fMRI activity during retro videogame play guided by automated annotations of game content.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: pandas, SciPy, Matplotlib, 1 other tool, 1 reference
[9] doi:10.3389/fnins.2026.1803897 [code]
Multimodal imaging-based targeting approach for network-level brain stimulation.
Journal: Frontiers in neuroscience
In common: SciPy, Matplotlib, NumPy, other, systems, 1 reference
[10] doi:10.1093/psyrad/kkag015 [code]
Neurocomputational mechanisms of reward-based online mood regulation in adolescents with bipolar disorder and major depressive disorder.
Journal: Psychoradiology
In common: pandas, SciPy, Matplotlib, 1 other tool, 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.