OSCR

The Virtual Brain links transcranial magnetic stimulation evoked potentials and inhibitory neurotransmitter changes in major depressive disorder.

Code ↔ Paper

6 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 6 matches
  1. [1] § Methods and Materials › Whole-brain simulations with the Virtual Brain › Neural mass model ↔ functions.py, lines 88–96 · score 0.84 · firing threshold, neuronal population, rate operator, firing rate, steepness, NMM
  2. [2] § Methods and Materials › Data ↔ functions.py, lines 155–256 · score 0.67 · Rit neural mass, SC matrices, structural connectivity, transformed, Jansen, excitatory
  3. [3] § Methods and Materials › Whole-brain simulations with the Virtual Brain › Neural mass model ↔ functions.py, lines 88–96 · score 0.64 · neuronal populations, sigmoidal function, operator, NMM, Jansen, Rit
  4. [4] § Methods and Materials › Whole-brain simulations with the Virtual Brain › Neural mass model ↔ functions.py, lines 155–256 · score 0.63 · coupled activity, inhibitory population, transformed, noise, sigmoid, excitatory
  5. [5] § Methods and Materials › Whole-brain simulations with the Virtual Brain › Neural mass model ↔ 4_bids_conversion.py, lines 212–283 · score 0.60 · neural populations, sigmoid function, external, coupled, excitatory, synapses
  6. [6] § Methods and Materials › Data ↔ 4_bids_conversion.py, lines 212–283 · score 0.57 · connection weights, Rit neural mass, ms, excitatory, global, synapses

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 · 600 lines · 30 KB · EUPL-1.2 · 4 matches

  1. # ------------------------------------------------------------------------------
  2. # functions.py
  3. # Author: Dr. Timo Hofsähs
  4. #
  5. # Description:
  6. # This Python script is part of the code accompanying the scientific publication:
  7. # The Virtual Brain links transcranial magnetic stimulation evoked potentials and
  8. # neurotransmitter changes in major depressive disorder
  9. # Dr. Timo Hofsähs, Marius Pille, Dr. Jil Meier, Prof. Petra Ritter
  10. # (in prep)
  11. #
  12. # This code provides all functions necessary to perform the fitting in 'fitting.py'.
  13. # The fitting method and all functions in this script except from 'gmfa' and
  14. # 'gmfa_timepoint' are based on the code provided with the following publication:
  15. # Momi D, Wang Z, Griffiths JD. 2023. TMS-EEG evoked responses are driven by recurrent
  16. # large-scale network dynamics. eLife2023;12:e83232 DOI: https://doi.org/10.7554/eLife.83232
  17. # Licensed under a Creative Commons Attribution license (CC-BY)
  18. # The original code can be found at:
  19. # https://github.com/GriffithsLab/PyTepFit/blob/main/tepfit/fit.py
  20. #
  21. # Copyright (c) 2024 Dr. Timo Hofsähs. All rights reserved.
  22. #
  23. # License: This code is licensed under the Creative Commons Attribution 4.0 International
  24. # License (CC-BY 4.0), which allows for redistribution, adaptation, and use in source
  25. # and binary forms, with or without modification, provided proper credit is given to
  26. # the original authors. You can view the full terms of this license at:
  27. # https://creativecommons.org/licenses/by/4.0/
  28. # ------------------------------------------------------------------------------
  30. import numpy as np
  31. import scipy.stats
  32. import torch
  33. import torch.optim as optim
  34. from torch.nn.parameter import Parameter
  35. from sklearn.metrics.pairwise import cosine_similarity
  36. import pickle
  38. class OutputNM():
  39. mode_all = ['train', 'test']
  40. stat_vars_all = ['m', 'v']
  42. def __init__(self, model_name, node_size, param, fit_weights=False):
  43. self.loss = np.array([])
  44. if model_name == "JR":
  45. state_names = ['E', 'Ev', 'I', 'Iv', 'P', 'Pv']
  46. self.output_name = "eeg"
  47. for name in state_names + [self.output_name]:
  48. for m in self.mode_all:
  49. setattr(self, name + '_' + m, [])
  51. vars = [a for a in dir(param) if not a.startswith('__') and not callable(getattr(param, a))]
  52. for var in vars:
  53. if np.any(getattr(param, var)[1] > 0):
  54. if var != 'std_in':
  55. setattr(self, var, np.array([]))
  56. for stat_var in self.stat_vars_all:
  57. setattr(self, var + '_' + stat_var, [])
  58. else:
  59. setattr(self, var, [])
  61. self.weights = []
  63. def save(self, filename):
  64. with open(filename, 'wb') as f:
  65. pickle.dump(self, f)
  68. class ParamsJR():
  69. '''Jansen & Rit neural mass model'''
  70. def __init__(self, model_name, **kwargs):
  71. if model_name == "JR":
  72. param = {"A ": [3.25, 0], "a": [100, 0.], "B": [22, 0], "b": [50, 0], "g": [100, 0],"c1": [135, 0.], "c2": [135 * 0.8, 0.],"c3": [135 * 0.25, 0.],"c4": [135 * 0.25, 0.],"std_in": [100, 0], "vmax": [5, 0], "v0": [6, 0], "r": [0.56, 0], "mu": [.5, 0],"speed": [2.5, 0],"k": [5, 0],"ki": [1, 0]}
  73. for var in param:
  74. setattr(self, var, param[var])
  75. for var in kwargs:
  76. setattr(self, var, kwargs[var])
  78. def sys2nd(A, a, u, x, v):
  79. '''
  80. Jansen & Rit neural mass model equation
  81. A = maximum amplitude of PSP (A=EPSP, B=IPSP)
  82. a = maximum firing rate of populations (a = excitatory = PC, EIN; b = inhibitory = IIN)
  83. u = input activity
  84. x,v = state variable activity
  85. '''
  86. return A*a*u -2*a*v-a**2*x
  88. def sigmoid(x, vmax, v0, r):
  89. '''
  90. Jansen & Rit neural mass model sigmoid function (potential to rate operator)
  91. x = input membrane potential
  92. vmax = maximum firing rate of all neuronal populations (default = 5 s^-1)
  93. v0 = PSP for which half of the maximum firing rate of the neuronal population is achieved; can be interpreted as excitability of all neuronal populations (default = 6mV)
  94. r = steepness at the firing threshold; represents the variance of firing thresholds within the NMM (default = 0.56 mV^-1)
  95. '''
  96. return vmax/(1+torch.exp(r*(v0-x)))
  99. class RNNJANSEN(torch.nn.Module):
  100. '''
  101. Embedding of whole-brain simulation parameters in fitting algorithm.
  102. Defines parameters of the whole-brain simulation and if and to what extent they are included in the fitting.
  103. '''
  104. state_names = ['E', 'Ev', 'I', 'Iv', 'P', 'Pv']
  105. model_name = "JR"
  107. def __init__(self,
  108. input_size: int,
  109. node_size: int,
  110. batch_size: int,
  111. step_size: float,
  112. output_size: int,
  113. tr: float,
  114. sc: float,
  115. lm: float,
  116. dist: float,
  117. param: ParamsJR,
  118. seed_int: int) -> None:
  120. super(RNNJANSEN, self).__init__()
  121. self.state_size = 6 # number of state variables
  122. self.input_size = input_size # 1 or 2 or 3
  123. self.tr = tr # tr ms (integration step 0.1 ms)
  124. self.step_size = torch.tensor(step_size, dtype=torch.float32) # integration step 0.1 ms
  125. self.hidden_size = int(tr / step_size)
  126. self.batch_size = batch_size # size of batch in ms applied at each fittin iteration
  127. self.node_size = node_size # number of regions of the whole-brain simulation
  128. self.output_size = output_size # number of EEG channels
  129. self.sc = sc # structural connectivity weights matrix (shaped node_size^2)
  130. self.dist = torch.tensor(dist, dtype=torch.float32) # structural connectivity tract length matrix (shaped node_size^2)
  131. self.param = param # parameters of the whole-brain simulation
  132. self.seed_int = seed_int # number of the seed
  134. np.random.seed(self.seed_int) # define seed
  135. self.lm = torch.tensor(lm, dtype=torch.float32) # leadfield matrix from sourced data to eeg
  137. # define if a parameter is fitted or not and add variance of default parameter if defined
  138. vars = [a for a in dir(param) if not a.startswith('__') and not callable(getattr(param, a))]
  139. for var in vars:
  140. if np.any(getattr(param, var)[1] > 0):
  141. setattr(self, var, Parameter(torch.tensor(getattr(param, var)[0] + getattr(param, var)[2] * np.random.randn(1,)[0], dtype=torch.float32)))
  142. dict_nv = {}
  143. dict_nv['m'] = getattr(param, var)[0]
  144. dict_nv['v'] = getattr(param, var)[1]
  146. dict_np = {}
  147. dict_np['m'] = var + '_m'
  148. dict_np['v'] = var + '_v'
  150. for key in dict_nv:
  151. setattr(self, dict_np[key], Parameter(torch.tensor(dict_nv[key], dtype=torch.float32)))
  152. else:
  153. setattr(self, var, torch.tensor(getattr(param, var)[0], dtype=torch.float32))
  155. def forward(self, input, noise_in, noise_out, hx, hE):
  156. '''
  157. Forward function of the concurrent simulation and fitting.
  158. Calculates the state variables of the neural mass model for every time step.
  159. '''
  160. bound_coupling = 500 # upper and lower limit of the coupling function
  161. dt = self.step_size
  162. self.delays = (self.dist / self.speed).type(torch.int64) # delays of the coupling per region
  164. w = torch.exp(self.w_bb) * torch.tensor(self.sc, dtype=torch.float32) # multiplication of SC weights matrix with fitting transformation factor w_bb
  165. w_n = torch.log1p(0.5 * (w + w.T)) / torch.linalg.norm(torch.log1p(0.5 * (w + w.T))) # normalization, logarithmization, symmetrication of altered SC weights matrix
  166. self.sc_m = w_n
  167. dg = -torch.diag(torch.sum(w_n, axis=1)) # negative of the sum of all SC weights per region
  169. M = hx[:, 0:1] # current of main population
  170. E = hx[:, 1:2] # current of excitory population
  171. I = hx[:, 2:3] # current of inhibitory population
  172. Mv = hx[:, 3:4] # voltage of main population
  173. Ev = hx[:, 4:5] # voltage of exictory population
  174. Iv = hx[:, 5:6] # voltage of inhibitory population
  176. current_state = torch.zeros_like(hx)
  177. next_state = {}
  178. eeg_batch = []
  179. E_batch = []
  180. I_batch = []
  181. M_batch = []
  182. Ev_batch = []
  183. Iv_batch = []
  184. Mv_batch = []
  187. for i_batch in range(self.batch_size):
  188. # noiseEEG = noise_out[:, i_batch:i_batch + 1]
  189. for i_hidden in range(self.hidden_size):
  190. # Define SC matrix & coupling
  191. hE_new = hE.clone() # hE = history object, history of excitatory state variable
  192. Ed = torch.tensor(np.zeros((self.node_size, self.node_size)), dtype=torch.float32) # delays of E
  193. Ed = hE_new.gather(1, self.delays) # delay of E as input through structural connectivity
  194. LEd = torch.reshape(torch.sum(w_n * torch.transpose(Ed, 0, 1), 1), (self.node_size, 1)) # weighted delayed input from E
  196. # Define noise
  197. noiseE = noise_in[:, i_hidden, i_batch, 0:1] * self.std_in
  198. # noiseI = noise_in[:, i_hidden, i_batch, 1:2] * self.std_in
  199. # noiseM = noise_in[:, i_hidden, i_batch, 2:3] * self.std_in
  201. # Define stimulus
  202. u = input[:, i_hidden:i_hidden + 1, i_batch]
  203. stimulus = self.k * self.ki * u
  205. # Define activity input for Jansen & Rit equations
  206. rM = sigmoid(E - I, self.vmax, self.v0, self.r)
  207. rE = self.c2 * sigmoid(self.c1 * M , self.vmax, self.v0, self.r)
  208. rI = self.c4 * sigmoid(self.c3 * M, self.vmax, self.v0, self.r)
  210. # Define coupling
  211. coupled_activity = self.g * (LEd + torch.matmul(dg, E-I))
  212. coupling = bound_coupling * torch.tanh(coupled_activity / bound_coupling)
  214. # Jansen & Rit neural mass model equations
  215. ddM = M + dt * Mv
  216. ddE = E + dt * Ev
  217. ddI = I + dt * Iv
  218. ddMv = Mv + dt * sys2nd(self.A, self.a, rM, M, Mv)
  219. ddEv = Ev + dt * sys2nd(self.A, self.a, (rE + coupling + stimulus + noiseE + self.mu), E, Ev)
  220. ddIv = Iv + dt * sys2nd(self.B, self.b, rI, I, Iv)
  222. E = ddE
  223. I = ddI
  224. M = ddM
  225. Ev = ddEv
  226. Iv = ddIv
  227. Mv = ddMv
  229. hE[:, 0] = E[:, 0] - I[:, 0] # create new time step for history of E object
  231. # Put M E I Mv Ev and Iv at every tr to the placeholders for checking them visually.
  232. M_batch.append(M)
  233. I_batch.append(I)
  234. E_batch.append(E)
  235. Mv_batch.append(Mv)
  236. Iv_batch.append(Iv)
  237. Ev_batch.append(Ev)
  239. hE = torch.cat([(E - I), hE[:, :-1]], axis=1) # update placeholders for E buffer
  241. # Put the EEG signal each tr to the placeholder being used in the cost calculation.
  242. eeg_currently = 0.0005 * torch.matmul(self.lm, E - I)
  243. eeg_batch.append(eeg_currently)
  245. # Update the current state.
  246. current_state = torch.cat([M, E, I, Mv, Ev, Iv], axis=1)
  247. next_state['current_state'] = current_state
  248. next_state['eeg_batch'] = torch.cat(eeg_batch, axis=1)
  249. next_state['E_batch'] = torch.cat(E_batch, axis=1)
  250. next_state['I_batch'] = torch.cat(I_batch, axis=1)
  251. next_state['P_batch'] = torch.cat(M_batch, axis=1)
  252. next_state['Ev_batch'] = torch.cat(Ev_batch, axis=1)
  253. next_state['Iv_batch'] = torch.cat(Iv_batch, axis=1)
  254. next_state['Pv_batch'] = torch.cat(Mv_batch, axis=1)
  256. return next_state, hE
  258. class Costs:
  259. '''
  260. Class that defines equations for the cost function of the fittin algorithm
  261. '''
  262. def __init__(self, method):
  263. self.method = method
  265. def cost_dist(self, sim, emp):
  266. losses = torch.sqrt(torch.mean((sim - emp) ** 2))
  267. return losses
  269. def cost_beamform(self, model, emp):
  270. corr = torch.matmul(emp, emp.T)
  271. corr_inv = torch.inverse(corr)
  272. corr_inv_s = torch.inverse(torch.matmul(model.lm.T, torch.matmul(corr_inv, model.lm)))
  273. W = torch.matmul(corr_inv_s, torch.matmul(model.lm.T, corr_inv))
  274. return torch.trace(torch.matmul(W, torch.matmul(corr, W.T)))
  276. def cost_r(self, logits_series_tf, labels_series_tf):
  277. node_size = logits_series_tf.shape[0]
  278. labels_series_tf_n = labels_series_tf - torch.reshape(torch.mean(labels_series_tf, 1), [node_size, 1]) # remove mean across time
  279. logits_series_tf_n = logits_series_tf - torch.reshape(torch.mean(logits_series_tf, 1), [node_size, 1])
  280. cov_sim = torch.matmul(logits_series_tf_n, torch.transpose(logits_series_tf_n, 0, 1)) # correlation
  281. cov_def = torch.matmul(labels_series_tf_n, torch.transpose(labels_series_tf_n, 0, 1))
  282. FC_sim_T = torch.matmul(torch.matmul(torch.diag(torch.reciprocal(torch.sqrt(torch.diag(cov_sim)))), cov_sim), torch.diag(torch.reciprocal(torch.sqrt(torch.diag(cov_sim))))) # fc for sim and empirical BOLDs
  283. FC_T = torch.matmul(torch.matmul(torch.diag(torch.reciprocal(torch.sqrt(torch.diag(cov_def)))), cov_def), torch.diag(torch.reciprocal(torch.sqrt(torch.diag(cov_def)))))
  284. ones_tri = torch.tril(torch.ones_like(FC_T), -1) # mask for lower triangle without diagonal
  285. zeros = torch.zeros_like(FC_T) # create a tensor all ones
  286. mask = torch.greater(ones_tri, zeros) # boolean tensor, mask[i] = True iff x[i] > 1
  287. FC_tri_v = torch.masked_select(FC_T, mask) # mask out fc to vector with elements of the lower triangle
  288. FC_sim_tri_v = torch.masked_select(FC_sim_T, mask)
  289. FC_v = FC_tri_v - torch.mean(FC_tri_v) # remove the mean across the elements
  290. FC_sim_v = FC_sim_tri_v - torch.mean(FC_sim_tri_v)
  291. corr_FC = torch.sum(torch.multiply(FC_v, FC_sim_v)) * torch.reciprocal(torch.sqrt(torch.sum(torch.multiply(FC_v, FC_v)))) * torch.reciprocal(torch.sqrt(torch.sum(torch.multiply(FC_sim_v, FC_sim_v)))) # corr_coef
  292. losses_corr = -torch.log(0.5000 + 0.5 * corr_FC) # use surprise: corr to calculate probability and -log
  293. return losses_corr
  295. def cost_eff(self, model, sim, emp):
  296. if self.method == 0: # in the current version this methd is applied
  297. return self.cost_dist(sim, emp)
  298. elif self.method == 1:
  299. return self.cost_beamform(model, emp) + self.cost_dist(sim, emp)
  300. else:
  301. return self.cost_r(sim, emp)
  304. class Model_fitting:
  305. '''Fitting algorithm equations'''
  306. def __init__(self, model, ts, num_epoches, cost):
  307. self.model = model
  308. self.num_epoches = num_epoches
  309. self.ts = ts
  310. self.cost = Costs(cost)
  312. def save(self, filename):
  313. with open(filename, 'wb') as f:
  314. pickle.dump(self, f)
  316. def train(self, u=0):
  317. '''Train function of the fitting algorithm'''
  318. # define some constants
  319. lb = 0.001
  320. delays_max = 500
  321. state_ub = 2 # lower bound for initial conditions
  322. state_lb = 0.5 # upper bound for initial conditions
  323. w_cost = 1.0 # factor to scale cost of similarity against cost of parameter change
  324. epoch_min = 200 # minimum amount of epochs to run, stop criterium
  325. r_lb = 0.85 # minimum pearson correlation value, stop criterium
  326. self.u = u # stimulus
  328. self.output_sim = OutputNM(self.model.model_name, # placeholder for output(EEG and histoty of model parameters and loss)
  329. self.model.node_size,
  330. self.model.param)
  332. # define an optimizor(ADAM)
  333. optimizer = optim.Adam(self.model.parameters(), lr=0.05, eps=1e-7)
  335. X = torch.tensor(np.random.uniform(state_lb, state_ub, (self.model.node_size, self.model.state_size)), dtype=torch.float32) # initial conditions
  336. hE = torch.tensor(np.random.uniform(state_lb, state_ub, (self.model.node_size, delays_max)), dtype=torch.float32) # history object
  338. # define masks for geting lower triangle matrix
  339. mask = np.tril_indices(self.model.node_size, -1)
  340. mask_e = np.tril_indices(self.model.output_size, -1)
  342. # placeholders for the history of model parameters
  343. fit_param = {}
  344. fit_sc = [self.model.sc[mask].copy()] # sc weights history
  345. for key, value in self.model.state_dict().items():
  346. fit_param[key] = [value.detach().numpy().ravel().copy()]
  348. loss_his = []
  349. num_batches = int(self.ts.shape[2] / self.model.batch_size) # define num_batches
  351. # empty dicts and arrays for course of parameter values, gradients, similarity metrics and loss
  352. parameters_fitted_initial = {}
  353. for key, value in self.model.state_dict().items():
  354. if '_' not in key:
  355. parameters_fitted_initial[key] = value.detach().numpy().astype(np.float32)
  356. parameters_fitted_export = {}
  357. parameters_fitted_export[-1] = parameters_fitted_initial
  358. # parameters_gradients_export = {}
  359. course_loss = np.zeros((3, self.num_epoches))
  360. course_cos_sim = np.zeros((self.num_epoches,))
  361. course_pcc = np.zeros((self.num_epoches,))
  362. course_sc_values = np.zeros((self.num_epoches+1, self.model.node_size, self.model.node_size))
  363. course_sc_values[0] = self.model.sc
  365. for i_epoch in range(self.num_epoches):
  366. eeg = self.ts[i_epoch % self.ts.shape[0]]
  367. # Create placeholders for the simulated EEG E I M Ev Iv and Mv of entire time series.
  368. for name in self.model.state_names + [self.output_sim.output_name]:
  369. setattr(self.output_sim, name + '_train', [])
  371. external = torch.tensor(np.zeros([self.model.node_size, self.model.hidden_size, self.model.batch_size]), dtype=torch.float32)
  373. # Perform the training in batches
  374. for i_batch in range(num_batches):
  376. optimizer.zero_grad() # Reset the gradient to zeros after update model parameters.
  378. # Get the input and output noises for the module.
  379. noise_in = torch.tensor(np.random.randn(self.model.node_size, self.model.hidden_size, self.model.batch_size, self.model.input_size), dtype=torch.float32)
  380. noise_out = torch.tensor(np.random.randn(self.model.node_size, self.model.batch_size), dtype=torch.float32)
  382. if not isinstance(self.u, int):
  383. external = torch.tensor((self.u[:, :, i_batch * self.model.batch_size:(i_batch + 1) * self.model.batch_size]), dtype=torch.float32)
  385. # Use the model.forward() function to update next state and get simulated EEG in this batch.
  386. next_batch, hE_new = self.model(external, noise_in, noise_out, X, hE)
  388. # Get the batch of emprical EEG signal.
  389. ts_batch = torch.tensor((eeg.T[i_batch * self.model.batch_size:(i_batch + 1) * self.model.batch_size, :]).T, dtype=torch.float32)
  391. loss_prior = []
  392. m = torch.nn.ReLU() # define the relu function
  393. variables_p = [a for a in dir(self.model.param) if
  394. not a.startswith('__') and not callable(getattr(self.model.param, a))]
  396. # get penalty on each fitted model parameter based on the derivation from default
  397. for var in variables_p:
  398. if np.any(getattr(self.model.param, var)[1] > 0):
  399. dict_np = {}
  400. dict_np['m'] = var + '_m'
  401. dict_np['v'] = var + '_v'
  402. loss_prior.append(torch.sum((lb + (1 / m(self.model.get_parameter(dict_np['v'])))) * (m(self.model.get_parameter(var)) - m(self.model.get_parameter(dict_np['m']))) ** 2))
  403. # calculate total loss
  404. loss_similarity = w_cost * self.cost.cost_eff(self.model, next_batch['eeg_batch'], ts_batch) # loss from difference between simulated and empirical timeseries
  405. loss_complexity = sum(loss_prior) # loss from derivation of fitted parameters from default
  406. loss = loss_similarity + loss_complexity
  408. # Put the batch of the simulated EEG, E I M Ev Iv Mv in to placeholders for entire time-series.
  409. for name in self.model.state_names + [self.output_sim.output_name]:
  410. name_next = name + '_batch'
  411. tmp_ls = getattr(self.output_sim, name + '_train')
  412. tmp_ls.append(next_batch[name_next].detach().numpy())
  413. setattr(self.output_sim, name + '_train', tmp_ls)
  415. loss_his.append(loss.detach().numpy())
  416. loss.backward(retain_graph=True) # Calculate gradient using backward (backpropagation) method of the loss function.
  417. optimizer.step() # Optimize the model based on the gradient method in updating the model parameters.
  418. for key, value in self.model.state_dict().items(): # Put the updated model parameters into the history placeholders.
  419. fit_param[key].append(value.detach().numpy().ravel().copy())
  420. fit_sc.append(self.model.sc_m.detach().numpy()[mask].copy()) # add newly fitted SC to list
  421. X = torch.tensor(next_batch['current_state'].detach().numpy(), dtype=torch.float32)
  422. hE = torch.tensor(hE_new.detach().numpy(), dtype=torch.float32) # update history object
  424. fc = np.corrcoef(self.ts.mean(0))
  425. tmp_ls = getattr(self.output_sim, self.output_sim.output_name + '_train')
  426. ts_sim = np.concatenate(tmp_ls, axis=1)
  427. fc_sim = np.corrcoef(ts_sim[:, 10:])
  429. for name in self.model.state_names + [self.output_sim.output_name]:
  430. tmp_ls = getattr(self.output_sim, name + '_train')
  431. setattr(self.output_sim, name + '_train', np.concatenate(tmp_ls, axis=1))
  432. self.output_sim.loss = np.array(loss_his)
  434. # store current SC and current SC gradients
  435. current_sc = np.zeros((200,200))
  436. current_sc_mask = np.tril_indices(200,-1)
  437. current_sc[current_sc_mask] = np.array(fit_sc)[-10:,:].mean(0)
  438. current_sc = current_sc+current_sc.T
  439. course_sc_values[i_epoch+1] = current_sc
  440. # course_sc_gradients[i_epoch] = self.model.w_bb.grad.detach().numpy()
  442. # store current parameter values and gradients for export
  443. parameters_fitted_epoch = {}
  444. # parameters_gradients_epoch = {}
  445. for key, value in self.model.state_dict().items():
  446. if '_' not in key:
  447. parameters_fitted_epoch[key] = value.detach().numpy().astype(np.float32)
  448. # parameters_gradients_epoch[key] = np.array(value.grad, dtype=np.float32)
  449. parameters_fitted_export[i_epoch] = parameters_fitted_epoch
  450. # parameters_gradients_export[i_epoch] = parameters_gradients_epoch
  452. # store current loss values
  453. course_loss[0] = loss_similarity.detach().numpy()
  454. course_loss[1] = loss_complexity.detach().numpy()
  455. course_loss[2] = loss.detach().numpy()
  457. # store and print current cosine similarity value
  458. course_cos_sim[i_epoch] = np.round(np.diag(cosine_similarity(ts_sim, self.ts.mean(0))).mean(), 4)
  459. print(f'epoch: {i_epoch} - cosine similarity: {np.round(course_cos_sim[i_epoch], 4)}')
  461. pcc_value = np.zeros((62,))
  462. for j in range(62):
  463. channel_emp = self.ts.mean(0).T[:,j]
  464. channel_sim = ts_sim.T[:,j]
  465. r, p = scipy.stats.pearsonr(channel_emp, channel_sim)
  466. pcc_value[j] = r
  467. course_pcc[i_epoch] = np.round(np.mean(pcc_value), 4)
  468. # print(f'epoch: {i_epoch} - pcc: {course_pcc[i_epoch]}')
  470. # FITTING CHECK
  471. # if i_epoch > 0 and i_epoch % 10 == 0:
  472. # # PLOT FITTED HETEROGENEOUS PARAMETERS
  473. # fig, axs = plt.subplots(len(parameters_fitted_epoch), 1, figsize=(5, 7), dpi=150)
  474. # for i, key in enumerate(parameters_fitted_epoch):
  475. # axs[i].plot(np.array(parameters_fitted_epoch[key]))
  476. # axs[i].set_ylabel(key)
  477. # yticks = axs[i].get_yticks()
  478. # axs[i].set_yticklabels(['{:,.2f}'.format(y) for y in yticks])
  479. # plt.tight_layout()
  480. # plt.show()
  482. # # PLOT FITTED SC
  483. # plt.figure(figsize=(5, 5), dpi=200)
  484. # plt.imshow(current_sc)
  485. # plt.colorbar()
  486. # plt.title('SC')
  487. # plt.show()
  489. # # PLOT TIMESERIES
  490. # fig, axs = plt.subplots(1, 2, figsize=(8, 2), dpi=70)
  491. # for j in range(200):
  492. # axs[0].plot(np.arange(-100, 300), (E_train[j] - I_train[j]).T, alpha=0.2, c='red' if self.model.ki[j] != 0 else 'black')
  493. # axs[0].set_title('fitted RAW')
  494. # axs[0].set_xlim(-100, 300)
  495. # axs[1].plot(np.arange(-100, 300), ts_sim.T)
  496. # axs[1].set_title('fitted EEG')
  497. # axs[1].set_xlim(-100, 300)
  498. # plt.tight_layout()
  499. # plt.show()
  501. if i_epoch > epoch_min and np.corrcoef(fc_sim[mask_e], fc[mask_e])[0, 1] > r_lb:
  502. break
  504. # store outputs
  505. self.output_sim.weights = np.array(fit_sc)
  506. self.output_sim.course_parameter_values = parameters_fitted_export
  507. # self.output_sim.course_parameter_gradients = parameters_gradients_export
  508. self.output_sim.course_sc_values = course_sc_values
  509. # self.output_sim.course_sc_gradients = course_sc_gradients
  510. self.output_sim.course_loss = course_loss
  511. self.output_sim.course_cos_sim = course_cos_sim
  512. self.output_sim.course_pcc = course_pcc
  514. for key, value in fit_param.items():
  515. setattr(self.output_sim, key, np.array(value))
  517. def test(self, x0, he0, base_batch_num, u=0):
  518. '''Test function of the fitting algorithm'''
  519. transient_num = 10
  520. self.u = u
  522. # initial state
  523. X = torch.tensor(x0, dtype=torch.float32)
  524. hE = torch.tensor(he0, dtype=torch.float32)
  526. num_batches = int(self.ts.shape[2] / self.model.batch_size) + base_batch_num
  527. # Create placeholders for the simulated BOLD E I x f and q of entire time series.
  528. for name in self.model.state_names + [self.output_sim.output_name]:
  529. setattr(self.output_sim, name + '_test', [])
  531. u_hat = np.zeros((self.model.node_size, self.model.hidden_size, base_batch_num * self.model.batch_size + self.ts.shape[2]))
  532. u_hat[:, :, base_batch_num * self.model.batch_size:] = self.u
  534. # Perform the training in batches.
  535. for i_batch in range(num_batches):
  536. # Get the input and output noises for the module.
  537. noise_in = torch.tensor(np.random.randn(self.model.node_size, self.model.hidden_size, self.model.batch_size, self.model.input_size), dtype=torch.float32)
  538. noise_out = torch.tensor(np.random.randn(self.model.node_size, self.model.batch_size), dtype=torch.float32)
  539. external = torch.tensor((u_hat[:, :, i_batch * self.model.batch_size:(i_batch + 1) * self.model.batch_size]), dtype=torch.float32)
  541. # Use the model.forward() function to update next state and get simulated EEG in this batch.
  542. next_batch, hE_new = self.model(external, noise_in, noise_out, X, hE)
  544. if i_batch > base_batch_num - 1:
  545. for name in self.model.state_names + [self.output_sim.output_name]:
  546. name_next = name + '_batch'
  547. tmp_ls = getattr(self.output_sim, name + '_test')
  548. tmp_ls.append(next_batch[name_next].detach().numpy())
  549. setattr(self.output_sim, name + '_test', tmp_ls)
  551. X = torch.tensor(next_batch['current_state'].detach().numpy(), dtype=torch.float32)
  552. hE = torch.tensor(hE_new.detach().numpy(), dtype=torch.float32)
  554. fc = np.corrcoef(self.ts.mean(0))
  555. tmp_ls = getattr(self.output_sim, self.output_sim.output_name + '_test')
  556. ts_sim = np.concatenate(tmp_ls, axis=1)
  558. fc_sim = np.corrcoef(ts_sim[:, transient_num:])
  560. for name in self.model.state_names + [self.output_sim.output_name]:
  561. tmp_ls = getattr(self.output_sim, name + '_test')
  562. setattr(self.output_sim, name + '_test', np.concatenate(tmp_ls, axis=1))
  565. # GMFA
  566. def gmfa(data, start_time, end_time):
  567. '''
  568. Calculates the Global Mean Field Amplitude of an EEG time series for a defined time interval
  569. data: EEG data (n_timepoints, n_channels)
  570. start_time: Start index in the time series
  571. end_time: End index in the time series
  572. all: array shaped (n_timepoints, ), contains the gmfp for every time point
  573. '''
  574. interval_data = data[start_time:end_time, :]
  575. mean_values = np.mean(interval_data, axis=1)
  576. deviations = interval_data - mean_values[:, np.newaxis]
  577. squared_deviations = deviations ** 2
  578. sum_squared_deviations = np.sum(squared_deviations, axis=1)
  579. all = np.sqrt(sum_squared_deviations / interval_data[0].shape[0])
  580. return all
  582. def gmfa_timepoint(data, timepoint):
  583. '''
  584. Calculates the Global Mean Field Amplitude of an EEG time series for a defined time interval
  585. data: EEG data (n_timepoints, n_channels)
  586. start_time: Start index in the time series
  587. end_time: End index in the time series
  589. Returns 2 arrays
  590. all: array shaped (n_timepoints, ), contains the gmfp for every time point
  591. avg: average value of all time points
  592. '''
  593. data_timepoint = data[timepoint, :]
  594. mean_value = np.mean(data_timepoint)
  595. deviations = data_timepoint - mean_value
  596. squared_deviations = deviations ** 2
  597. sum_squared_deviations = np.sum(squared_deviations)
  598. division = np.sqrt(sum_squared_deviations / data.shape[1])
  600. return division

functions.py at commit 78d5511, under EUPL-1.2 · at the source

Overview

Authors: Timo Hofsähs1,2, Marius Pille1,2, Lucas Kern1,2,3, Anuja Negi1,2,4,5, Jil Mona Meier1,2, Petra Ritter1,2,4,6,7
ORCID iDs: Timo Hofsähs, Marius Pille, Lucas Kern, Anuja Negi, Jil Mona Meier, Petra Ritter
  1. Berlin Institute of Health at Charité - Universitätsmedizin Berlin, Berlin, Germany
  2. Charité - Universitätsmedizin Berlin, Corporate Member of Freie Universität Berlin and Humboldt-Universität zu Berlin, Department of Neurology with Experimental Neurology, Brain Simulation Section, Berlin, Germany
  3. FIL Methods Group - Department for Imaging Neuroscience (formerly The Wellcome Centre for Human Neuroimaging), UCL Queen Square Institute of Neurology, University College London, London, United Kingdom
  4. Bernstein Focus State Dependencies of Learning and Bernstein Center for Computational Neuroscience, Berlin, Germany
  5. Technical University Berlin, Berlin, Germany
  6. Einstein Center for Neuroscience Berlin, Berlin, Germany
  7. Einstein Center Digital Future, Berlin, Germany
Institutions: Humboldt-Universität zu Berlin (Germany); Berlin Institute of Health at Charité - Universitätsmedizin Berlin (Germany); Freie Universität Berlin (Germany); Charité - Universitätsmedizin Berlin (Germany); Wellcome Centre for Human Neuroimaging (United Kingdom); UCL Queen Square Institute of Neurology (United Kingdom); University College London (United Kingdom); Bernstein Center for Computational Neuroscience Berlin (Germany); Technische Universität Berlin (Germany); Einstein Center Digital Future (Germany); Einstein Center for Neurosciences Berlin (Germany)
Journal: Imaging neuroscience (Cambridge, Mass.), volume 4, article IMAG.a.1147
Dates: received 16 April 2025; accepted 29 January 2026; published online 3 March 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1162/imag.a.1147 · PMID 41799679 · PMCID PMC12961309 · OpenAlex W7128439752
Open access: diamond, a free copy (OpenAlex)
Status: code verified
Categories: other (modality), human (organism), depression (population), computational (subfield)
Methods: Single-unit activity, calcium imaging
Keywords: GABA, major depressive disorder, neural mass model, TEP, TMS, whole-brain simulations
Topic: Transcranial Magnetic Stimulation Studies (Neurology, Neuroscience), according to OpenAlex
Funding: Deutsche Forschungsgemeinschaft (SFB-1436, SFB 936, SFB 1315, RI 2073/6-1, 504745852, 785907, 327654276, TRR-295, RI 2073/10-2, 424778381, RI 2073/9-1, 178316478, 425899996); Berlin Institute of Health; European Regional Development Fund (101147319, #945539, EU‐H2020, 826421, 785907); HORIZON EUROPE Framework Programme (101057655, 945539, 785907, SGA3 945539, 101058516, 101147319, 101137289); Horizon 2020 (EU H2020 Virtual Brain Cloud 826421, AISN – 101057655, EBRAIN-Health 101058516, Digital Europe TEF-Health 101100700, Virtual Brain Twin (101137289), EBRAINS2.0 (101147319), EBRAINS-PREP 101079717)
Citations: cited by 1 paper (Europe PMC); 108 references in the paper

Abstract

Transcranial magnetic stimulation evoked potentials (TEPs) show promise as a biomarker in major depressive disorder (MDD), but the origin of the increased TEP amplitude in these patients remains unclear. Gamma aminobutyric acid (GABA) may be involved, as TEP peak amplitude is known to increase with GABAergic activity in healthy controls. We employed a computational modeling approach to investigate this phenomenon. Whole-brain simulations in ‘The Virtual Brain’ (thevirtualbrain.org), employing the Jansen and Rit neural mass model, were optimized to simulate TEPs of healthy individuals (Nsubs = 20, 14 females, 24.5 ± 4.9 years). To mimic MDD-like impaired inhibition, a GABAergic deficit was introduced to the simulations by altering one of two selected inhibitory parameters, the inhibitory synaptic decay rate b or the number of inhibitory synapses C4. The TEP amplitude was quantified and compared for all simulations. The inhibitory synaptic decay rate showed a quadratic correlation (r = 0.99, p < 0.001) and the number of inhibitory synapses a negative exponential correlation (r = 0.99, p < 0.001) with the TEP amplitude. Moreover, significant correlations between these simulation-derived values and all TEP peaks and troughs were detected (p < 0.001). Thus, under local parameter changes, we were able to alter the TEP amplitude toward pathological levels, that is, creating an MDD-like increase of the global mean field amplitude in line with empirical results. Our model suggests specific GABAergic deficits as the cause of increased TEP amplitude in MDD patients, which may serve as therapeutic targets. This work highlights the potential of whole-brain simulations in the investigation of neuropsychiatric diseases.

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

Repository

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

virtual-twin/TMS_MDD

License: EUPL-1.2
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 78d55110979aa6621c2e25ce96a029b26c76f3b3, 8 January 2026
Languages: Python (6)
Size: 14 files, 6 scripts
Software Heritage: not archived
Found in: “Data and Code Availability”
Holds: README, license file, environment (environment.yml, requirements.txt)
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (6 files), pandas (4 files), SciPy (4 files), Numba (2 files), The Virtual Brain (2 files), Matplotlib (1 file), MNE-Python (1 file), PyTorch (1 file), scikit-learn (1 file), seaborn (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
8 files
At the source: github.com/virtual-twin/TMS_MDD

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

Tracing map

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

What the map holds:

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

The empirical TEP dataset (Biabani & Rogasch, 2019) is publicly available for download (https://bridges.monash.edu/articles/dataset/TEPs-SEPs/7440713?file=13772894). The entire code to generate the results presented in this study is openly available (https://github.com/virtual-twin/TMS_MDD).

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

Versions

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

Version 1, 30 September 2026: the first record

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

Cite

This paper

Hofsähs, T., Pille, M., Kern, L., Negi, A., Meier, J. M., & Ritter, P. (2026). The Virtual Brain links transcranial magnetic stimulation evoked potentials and inhibitory neurotransmitter changes in major depressive disorder. Imaging neuroscience (Cambridge, Mass.), 4, IMAG.a.1147. https://doi.org/10.1162/imag.a.1147

BibTeX

@article{hofsahs2026virtual,
author = {Hofsähs, Timo and Pille, Marius and Kern, Lucas and Negi, Anuja and Meier, Jil Mona and Ritter, Petra},
title = {{The Virtual Brain links transcranial magnetic stimulation evoked potentials and inhibitory neurotransmitter changes in major depressive disorder}},
journal = {Imaging neuroscience (Cambridge, Mass.)},
year = {2026},
month = mar,
volume = {4},
pages = {IMAG.a.1147},
publisher = {MIT Press},
issn = {2837-6056},
doi = {10.1162/imag.a.1147},
url = {https://doi.org/10.1162/imag.a.1147},
pmid = {41799679},
pmcid = {PMC12961309}
}

RIS

TY - JOUR
AU - Hofsähs, Timo
AU - Pille, Marius
AU - Kern, Lucas
AU - Negi, Anuja
AU - Meier, Jil Mona
AU - Ritter, Petra
TI - The Virtual Brain links transcranial magnetic stimulation evoked potentials and inhibitory neurotransmitter changes in major depressive disorder
T2 - Imaging neuroscience (Cambridge, Mass.)
J2 - Imaging Neurosci (Camb)
PY - 2026
DA - 2026/03/03
VL - 4
SP - IMAG.a.1147
SN - 2837-6056
PB - MIT Press
DO - 10.1162/imag.a.1147
UR - https://doi.org/10.1162/imag.a.1147
LA - en
ER -

CSL-JSON

{
"id": "10.1162/imag.a.1147",
"type": "article-journal",
"title": "The Virtual Brain links transcranial magnetic stimulation evoked potentials and inhibitory neurotransmitter changes in major depressive disorder",
"container-title": "Imaging neuroscience (Cambridge, Mass.)",
"author": [
{
"family": "Hofsähs",
"given": "Timo"
},
{
"family": "Pille",
"given": "Marius"
},
{
"family": "Kern",
"given": "Lucas"
},
{
"family": "Negi",
"given": "Anuja"
},
{
"family": "Meier",
"given": "Jil Mona"
},
{
"family": "Ritter",
"given": "Petra"
}
],
"container-title-short": "Imaging Neurosci (Camb)",
"volume": "4",
"page": "IMAG.a.1147",
"DOI": "10.1162/imag.a.1147",
"PMID": "41799679",
"PMCID": "PMC12961309",
"ISSN": "2837-6056",
"publisher": "MIT Press",
"URL": "https://doi.org/10.1162/imag.a.1147",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
3
]
]
}
}

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.1371/journal.pcbi.1013463 [code]
A multi-frequency whole-brain neural mass model with homeostatic feedback inhibition.
Journal: PLoS computational biology
In common: Numba, scikit-learn, SciPy, 2 other tools, computational, 9 references
[2] doi:10.1038/s41467-026-71918-7 [code]
Developmental disinhibition gates language lateralization in childhood.
Journal: Nature communications
In common: MNE-Python, PyTorch, seaborn, 5 other tools, 5 references
[3] doi:10.1073/pnas.2532072123 [code]
Spatially structured heterogeneity shapes large-scale cortical dynamics in a model of the human cortex.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: The Virtual Brain, SciPy, Matplotlib, 1 other tool, 5 references
[4] doi:10.1371/journal.pdig.0001445 [code]
A digital twin approach for simultaneous reconstruction of brain anatomy and dynamics from neural data.
Journal: PLOS digital health
In common: 8 references
[5] doi:10.1371/journal.pcbi.1014222 [code]
Neural population models for EEG: From Canonical models to alternative model structures.
Journal: PLoS computational biology
In common: computational, 7 references
[6] doi:10.1038/s41398-026-04107-1
Abnormal left prefrontal N100 and its relationship with fronto-limbic metabolism in major depressive disorder.
Journal: Translational psychiatry
In common: depression, other, 6 references
[7] doi:10.1162/netn.a.543 [code]
High-resolution Bayesian Virtual Epileptic Patient using neural field models.
Journal: Network neuroscience (Cambridge, Mass.)
In common: The Virtual Brain, MNE-Python, SciPy, 2 other tools, computational, 2 references
[8] doi:10.1038/s41593-026-02205-3 [code]
Competitive interactions shape mammalian brain network dynamics and computation.
Journal: Nature neuroscience
In common: PyTorch, seaborn, scikit-learn, 4 other tools, 3 references
[9] doi:10.1162/imag.a.1356
Sources of the N15 TMS-evoked potential following motor cortex stimulation localize rostrals to TMS-induced electric fields and depend on dose.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: other, 6 references
[10] doi:10.1162/imag.a.1222 [code]
Network-based near-scalp personalized brain stimulation targets.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: PyTorch, scikit-learn, SciPy, 2 other tools, other, 3 references

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.