OSCR

Cerebellar microcircuits enable robust evidence-based decisions through cortico-cerebellar coupling.

Code ↔ Paper

The paper beside its authors' code: matches between them have not been computed for this paper yet.

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 · 542 lines · 23 KB · no license

  1. # coding=gbk
  2. import multiprocessing
  3. import os
  4. import matplotlib
  5. matplotlib.use('Agg')
  6. import matplotlib.pyplot as plt
  7. import numpy as np
  8. import matplotlib.gridspec as gridspec
  9. from brian2 import *
  10. from collections import OrderedDict
  11. import cupy as cp
  13. continuous_model_params = dict(
  14. V_L = -70*mV, Vth = -50*mV, Vreset = -55*mV,
  15. gE = 25*nS, tau_m_E = 20*ms, tau_ref_E = 2*ms,
  16. gI = 20*nS, tau_m_I = 10*ms, tau_ref_I = 1*ms,
  17. V_E = 0*mV, V_I = -70*mV,
  18. a = 0.062*mV**-1, b = 3.57,
  19. tauAMPA = 2*ms, tau_x = 2*ms, tauNMDA = 100*ms, alpha = 0.5*kHz, tauGABA = 5*ms, delay = 0*ms,
  20. gAMPA_ext_E = 1*nS, gAMPA_ext_E_input = 150*nS, gAMPA_ext_I = 1.62*nS,
  21. gAMPA_E = 200*nS, gNMDA_E = 400*nS, gGABA_E = 10*nS,
  22. gAMPA_I = 300*nS, gNMDA_I = 0*nS, gGABA_I = 400*nS,
  23. nu_ext = 100*Hz,
  24. g_SFA = 100*nS, tau_w = 80*ms, b_SFA = 0.1, V_K = -80*mV,
  25. N_E = 1600, N_I = 400,
  26. fsel = 0.15,
  27. wp = 35, sigma_w=2.4
  28. )
  30. gc_par = dict(
  31. C = 100 * pF, vr = -60 * mV, vt = -40 * mV, k = 0.7 * pA/(mV**2),
  32. a = 0.03 / ms , b = -2 * nS , vpeak = 0* mV , c_reset = -50 * mV , d_reset = 100 * pA ,
  33. V_E=0*mV, gAMPA=50*nS, f = 0.062*mV**-1, g = 3.57, tauAMPA = 2*ms,
  34. possion_noise=0.5*kHz, gAMPA_noise=1*nS, tau_c=2*second,
  35. N_gc = 1000, mossy_fiber_input_num=4
  36. )
  38. PC_par=dict(
  39. N=400, c = 268 * pF, gl = 8.47 * nS, el = -51.31 * mV, vt = -53.23 * mV, delta = 0.85 * mV, vreset = -60.35 * mV,
  40. a = 37.79 * nS, tauw = 20.76 * ms, b = 441.12 * pA,
  41. Ihold = -140 * pA, VSpike=60*mV, V_E=0*mV, V_I=-80*mV,
  42. sigma=40*pA, corr=2*ms,
  43. tauAMPA = 2*ms, tauGABA = 5*ms,
  44. gAMPA=3*nS, gGABA =17*nS
  45. )
  47. pulse_params = dict(
  48. start_time=1000*ms, end_time=4800*ms, stim_num_sum=19, left_right_diff=1, start_index=800, end_index=900
  49. )
  51. sim_params=dict(
  52. sim_purations=7000*ms
  53. )
  55. E_equations = '''
  56. dV/dt = (-(V - V_L) - Isyn/gE) / tau_m_E : volt (unless refractory)
  57. Isyn = I_AMPA_ext + I_AMPA + I_NMDA + I_GABA+I_SFA+I_AMPA_input: amp
  58. I_AMPA_ext = gAMPA_ext_E*sAMPA_ext*(V - V_E) : amp
  59. I_AMPA_input = gAMPA_ext_E_input*sAMPA_input*(V - V_E) : amp
  60. I_AMPA = gAMPA_E*S_AMPA*(V - V_E) : amp
  61. I_NMDA = gNMDA_E*S_NMDA*(V - V_E)/(1 + exp(-a*V)/b) : amp
  62. I_GABA = gGABA_E*S_GABA*(V - V_I) : amp
  63. I_SFA = g_SFA * w_SFA * (V - V_K) : amp
  64. dsAMPA_ext/dt = -sAMPA_ext/tauAMPA : 1
  65. dsAMPA_input/dt = -sAMPA_input/tauAMPA : 1
  66. dsAMPA/dt = -sAMPA/tauAMPA : 1
  67. dx/dt = -x/tau_x : 1
  68. dw_SFA/dt = -w_SFA / tau_w : 1
  69. dsNMDA/dt = -sNMDA/tauNMDA + alpha*x*(1 - sNMDA) : 1
  70. S_AMPA : 1
  71. S_NMDA : 1
  72. S_GABA : 1
  73. '''
  75. I_equations = '''
  76. dV/dt = (-(V - V_L) - Isyn/gI) / tau_m_I : volt (unless refractory)
  77. Isyn = I_AMPA_ext + I_AMPA + I_NMDA + I_GABA : amp
  78. I_AMPA_ext = gAMPA_ext_I*sAMPA_ext*(V - V_E) : amp
  79. I_AMPA = gAMPA_I*S_AMPA*(V - V_E) : amp
  80. I_NMDA = gNMDA_I*S_NMDA*(V - V_E)/(1 + exp(-a*V)/b) : amp
  81. I_GABA = gGABA_I*S_GABA*(V - V_I) : amp
  82. dsAMPA_ext/dt = -sAMPA_ext/tauAMPA : 1
  83. dsGABA/dt = -sGABA/tauGABA : 1
  84. S_AMPA: 1
  85. S_NMDA: 1
  86. S_GABA: 1
  87. '''
  89. eqs_gc = '''
  90. dv/dt = (k*(v-vr)*(v-vt)- u + I_syn)/C : volt
  91. du/dt = a*(b*(v-vr)-u) : amp
  92. I_syn=I_AMPA+I_AMPA_noise: amp
  93. I_AMPA= gAMPA*sAMPA*(V_E-v) : amp
  94. I_AMPA_noise= gAMPA_noise*sAMPA_noise*(V_E-v) : amp
  95. dsAMPA/dt= -sAMPA/tauAMPA : 1
  96. dsAMPA_noise/dt= -sAMPA_noise/tauAMPA : 1
  97. '''
  99. eqs_pc='''
  100. dv/dt = (gl * delta * exp((v - vt) / delta) - gl * (v - el) - w + Ihold+I_ou+I_ext) / c : volt
  101. dw/dt = (a * (v - el) - w) / tauw : amp
  102. dtemp/dt = -temp/corr + (sqrt(2/(corr*dt)))*randn() : 1
  103. I_ou =temp * sigma : amp
  104. I_ext=I_GABA+I_AMPA:amp
  105. dsAMPA/dt = -sAMPA/tauAMPA : 1
  106. dsGABA/dt= -sGABA/tauGABA : 1
  107. I_AMPA=gAMPA*sAMPA*(V_E-v):amp
  108. I_GABA=gGABA*S_GABA*(V_I-v) : amp
  109. S_GABA:1
  110. '''
  112. current_script_dir = os.path.dirname(__file__)
  113. base_data_dir = os.path.join(current_script_dir, '..', 'evidence_accumulation', 'Data','Figure1')
  115. def calculate_sliding_window_rate(spike_times, num_neurons_in_group, total_duration, window_size, step_size):
  116. spike_times_ms = spike_times / ms
  117. total_duration_ms = total_duration / ms
  118. window_size_ms = window_size / ms
  119. step_size_ms = step_size / ms
  121. time_points_ms = np.arange(0, total_duration_ms, step_size_ms)
  122. firing_rates_Hz = np.zeros(len(time_points_ms))
  124. if num_neurons_in_group == 0:
  125. return time_points_ms, firing_rates_Hz
  127. for i, t_start_ms in enumerate(time_points_ms):
  128. t_end_ms = t_start_ms + window_size_ms
  129. spikes_in_window = np.sum((spike_times_ms >= t_start_ms) & (spike_times_ms < t_end_ms))
  130. if window_size_ms > 0:
  131. firing_rates_Hz[i] = spikes_in_window / (num_neurons_in_group * (window_size_ms / 1000.0))
  132. else:
  133. firing_rates_Hz[i] = 0 # If window size is 0, rate is 0
  135. return time_points_ms, firing_rates_Hz
  137. def calculate_weight(theta_i, theta_j, J_plus, J_minus, sigma_w_deg):
  138. delta_theta_deg = abs(theta_i - theta_j)
  139. delta_theta_deg = np.minimum(delta_theta_deg, 360 - delta_theta_deg)
  140. return J_minus + (J_plus - J_minus) * np.exp(-(delta_theta_deg**2) / (2 * sigma_w_deg**2))
  143. def generate_random_pulse_times_strict_constrained(start_t_ms, end_t_ms, num_pulses, proportional_jitter_factor=0.5):
  144. if num_pulses <= 0:
  145. return np.array([])
  146. if num_pulses == 1:
  147. return np.array([start_t_ms])
  149. min_interval = 200
  150. total_duration_ms = end_t_ms - start_t_ms
  152. if total_duration_ms < (num_pulses - 1) * min_interval:
  153. raise ValueError(
  154. f"²ÎÊý³åÍ»: ÎÞ·¨ÔÚ {total_duration_ms:.2f}ms ÄÚ·ÅÖà {num_pulses} ¸öÂö³å "
  155. f"ÇÒ±£³Ö×îС¼ä¸ô {min_interval}ms¡£"
  156. )
  158. free_time = total_duration_ms - (num_pulses - 1) * min_interval
  160. random_shifts = np.random.uniform(0, free_time, num_pulses - 1)
  161. random_shifts.sort()
  163. pulse_times = np.zeros(num_pulses)
  164. pulse_times[0] = start_t_ms
  166. for i in range(1, num_pulses):
  167. pulse_times[i] = start_t_ms + i * min_interval + random_shifts[i - 1]
  169. return pulse_times
  171. def run_continuous_layer(args):
  172. left_right_diff_val, iteration_num, gpu_id = args
  174. output_dir = os.path.join(base_data_dir, f'{left_right_diff_val}')
  175. os.makedirs(output_dir, exist_ok=True)
  177. txt_save_path = os.path.join(output_dir, f'E_spikes_run_{iteration_num}.txt')
  179. if os.path.exists(txt_save_path):
  180. print(f"Skipping Phase 1 for diff={left_right_diff_val}, run={iteration_num} (File already exists: {txt_save_path})")
  181. return
  183. start_scope()
  185. np.random.seed(iteration_num)
  187. try:
  188. cp.cuda.Device(gpu_id).use()
  189. except Exception as e:
  190. print(f"Phase 1 Warning: Could not set GPU {gpu_id}. Error: {e}")
  191. pass
  193. params = continuous_model_params.copy()
  194. pulse_params_local = pulse_params.copy()
  195. sim_params_local = sim_params.copy()
  196. pulse_params_local['left_right_diff'] = left_right_diff_val
  198. for par in ['gAMPA_E', 'gAMPA_I', 'gNMDA_E', 'gNMDA_I']:
  199. params[par] /= params['N_E']
  200. for par in ['gGABA_E', 'gGABA_I']:
  201. params[par] /= params['N_I']
  203. N_E = params['N_E']
  204. fsel = params['fsel']
  205. wp = params['wp']
  206. delay = params['delay']
  207. N1 = int(fsel*N_E)
  208. N2 = N1
  209. N0 = N_E - (N1 + N2)
  210. params['N0'] = N0
  211. params['N1'] = N1
  212. params['N2'] = N2
  213. wm = 0
  214. params['wm'] = wm
  216. weight_matrix = np.zeros((N_E, N_E))
  217. excitatory_thetas = np.linspace(0, 360, N_E, endpoint=False)
  218. for l in range(N_E):
  219. for m in range(N_E):
  220. weight_matrix[l, m] = calculate_weight(excitatory_thetas[l], excitatory_thetas[m], params['wp'], params['wm'], params['sigma_w'])
  221. weight_matrix_gpu = cp.asarray(weight_matrix)
  223. start_time_ms = pulse_params_local['start_time']/ms
  224. end_time_ms = pulse_params_local['end_time']/ms
  225. stim_num_left = (pulse_params_local['stim_num_sum'] - pulse_params_local['left_right_diff']) // 2
  226. stim_num_right = (pulse_params_local['stim_num_sum'] + pulse_params_local['left_right_diff']) // 2
  227. stim_num_left = max(0, stim_num_left)
  228. stim_num_right = max(0, stim_num_right)
  230. input_neuron_indices_left = np.array([], dtype=int)
  231. spike_times_left = np.array([]) * ms
  232. if stim_num_left > 0:
  233. try:
  234. pulse_times_ms_left = generate_random_pulse_times_strict_constrained(start_time_ms, end_time_ms, stim_num_left)
  235. input_neuron_indices_left = np.repeat(np.arange(N1), stim_num_left)
  236. spike_times_left = np.tile(pulse_times_ms_left, N1)*ms
  237. except ValueError as e:
  238. print(f"Simulation (diff={left_right_diff_val}, run={iteration_num}) for LEFT side failed: {e}")
  239. return
  241. input_neuron_indices_right = np.array([], dtype=int)
  242. spike_times_right = np.array([]) * ms
  243. if stim_num_right > 0:
  244. try:
  245. pulse_times_ms_right = generate_random_pulse_times_strict_constrained(start_time_ms, end_time_ms, stim_num_right)
  246. input_neuron_indices_right = np.repeat(np.arange(N1), stim_num_right)
  247. spike_times_right = np.tile(pulse_times_ms_right, N1)*ms
  248. except ValueError as e:
  249. print(f"Simulation (diff={left_right_diff_val}, run={iteration_num}) for RIGHT side failed: {e}")
  250. return
  252. net = OrderedDict()
  253. net['input_neurons_left'] = SpikeGeneratorGroup(N1, input_neuron_indices_left, spike_times_left)
  254. net['input_neurons_right'] = SpikeGeneratorGroup(N1, input_neuron_indices_right, spike_times_right)
  256. net['E']=NeuronGroup(2*params['N_E'],E_equations,method='rk2',threshold='V > Vth',reset='V = Vreset;w_SFA+=b_SFA',refractory=params['tau_ref_E'],namespace=params)
  257. net['E_left']= net['E'][0:params['N_E']]
  258. net['E_right']= net['E'][params['N_E']:]
  260. net['I_left']=NeuronGroup(params['N_I'],I_equations,method='rk2',threshold='V > Vth',reset='V = Vreset',refractory=params['tau_ref_I'],namespace=params)
  261. net['I_right']=NeuronGroup(params['N_I'],I_equations,method='rk2',threshold='V > Vth',reset='V = Vreset',refractory=params['tau_ref_I'],namespace=params)
  263. net['pgE_left'] = PoissonGroup(params['N_E'], params['nu_ext'])
  264. net['pgI_left'] = PoissonGroup(params['N_I'], params['nu_ext'])
  265. net['pgE_right'] = PoissonGroup(params['N_E'], params['nu_ext'])
  266. net['pgI_right'] = PoissonGroup(params['N_I'], params['nu_ext'])
  268. if len(input_neuron_indices_left) > 0:
  269. net['inE_left']=Synapses(net['input_neurons_left'], net['E_left'][pulse_params_local['start_index']:pulse_params_local['end_index']],on_pre='sAMPA_input += 1',delay=delay)
  270. net['inE_left'].connect(condition='i == j')
  271. if len(input_neuron_indices_right) > 0:
  272. net['inE_right']=Synapses(net['input_neurons_right'], net['E_right'][pulse_params_local['start_index']:pulse_params_local['end_index']],on_pre='sAMPA_input += 1',delay=delay)
  273. net['inE_right'].connect(condition='i == j')
  275. net['icE_left'] = Synapses(net['pgE_left'], net['E_left'],on_pre='sAMPA_ext += 1', delay=delay)
  276. net['icE_left'].connect(condition='i == j')
  277. net['icE_right'] = Synapses(net['pgE_right'], net['E_right'],on_pre='sAMPA_ext += 1', delay=delay)
  278. net['icE_right'].connect(condition='i == j')
  280. net['icI_left'] = Synapses(net['pgI_left'], net['I_left'],on_pre='sAMPA_ext += 1', delay=delay)
  281. net['icI_left'].connect(condition='i == j')
  282. net['icI_right'] = Synapses(net['pgI_right'], net['I_right'],on_pre='sAMPA_ext += 1', delay=delay)
  283. net['icI_right'].connect(condition='i == j')
  285. net['icAMPA_left'] = Synapses(net['E_left'], net['E_left'], on_pre='sAMPA += 1', delay=delay)
  286. net['icAMPA_left'].connect(condition='i == j')
  287. net['icNMDA_left'] = Synapses(net['E_left'], net['E_left'], on_pre='x += 1', delay=delay)
  288. net['icNMDA_left'].connect(condition='i == j')
  289. net['icGABA_left'] = Synapses(net['I_left'], net['I_left'], on_pre='sGABA += 1', delay=delay)
  290. net['icGABA_left'].connect(condition='i == j')
  292. net['icAMPA_right'] = Synapses(net['E_right'], net['E_right'], on_pre='sAMPA += 1', delay=delay)
  293. net['icAMPA_right'].connect(condition='i == j')
  294. net['icNMDA_right'] = Synapses(net['E_right'], net['E_right'], on_pre='x += 1', delay=delay)
  295. net['icNMDA_right'].connect(condition='i == j')
  296. net['icGABA_right'] = Synapses(net['I_right'], net['I_right'], on_pre='sGABA += 1', delay=delay)
  297. net['icGABA_right'].connect(condition='i == j')
  299. @network_operation(when='start')
  300. def recurrent_input():
  301. sAMPA_gpu_left = cp.asarray(net['E_left'].sAMPA)
  302. sNMDA_gpu_left = cp.asarray(net['E_left'].sNMDA)
  303. sAMPA_gpu_right = cp.asarray(net['E_right'].sAMPA)
  304. sNMDA_gpu_right = cp.asarray(net['E_right'].sNMDA)
  306. net['E_left'].S_AMPA = cp.dot(sAMPA_gpu_left, weight_matrix_gpu).get()
  307. net['E_left'].S_NMDA = cp.dot(sNMDA_gpu_left, weight_matrix_gpu).get()
  308. net['I_left'].S_AMPA=sum(net['E_left'].sAMPA)
  309. net['I_left'].S_NMDA=sum(net['E_left'].sNMDA)
  310. net['E_right'].S_AMPA = cp.dot(sAMPA_gpu_right, weight_matrix_gpu).get()
  311. net['E_right'].S_NMDA = cp.dot(sNMDA_gpu_right, weight_matrix_gpu).get()
  312. net['I_right'].S_AMPA=sum(net['E_right'].sAMPA)
  313. net['I_right'].S_NMDA=sum(net['E_right'].sNMDA)
  315. net['E_left'].S_GABA=sum(net['I_left'].sGABA)
  316. net['E_right'].S_GABA=sum(net['I_right'].sGABA)
  318. net['I_left'].S_GABA=sum(net['I_left'].sGABA)
  319. net['I_right'].S_GABA=sum(net['I_right'].sGABA)
  321. net['E_left'].V = np.random.uniform(params['Vreset'],params['Vth'],size=params['N_E']) * volt
  322. net['I_left'].V = np.random.uniform(params['Vreset'],params['Vth'],size=params['N_I']) * volt
  323. net['E_right'].V = np.random.uniform(params['Vreset'],params['Vth'],size=params['N_E']) * volt
  324. net['I_right'].V = np.random.uniform(params['Vreset'],params['Vth'],size=params['N_I']) * volt
  326. spikemon_E = SpikeMonitor(net['E'])
  328. defaultclock.dt=0.05*ms
  329. total_duration = sim_params_local['sim_purations']
  331. network = Network(list(net.values()) + [recurrent_input, spikemon_E])
  332. print(f"Continuous Model (diff={left_right_diff_val}, run={iteration_num}): Starting...")
  333. network.run(total_duration, report='text')
  335. spike_data = np.column_stack((spikemon_E.i, spikemon_E.t/ms))
  336. np.savetxt(txt_save_path, spike_data, fmt=['%d', '%.4f'])
  338. print(f"Continuous Model (diff={left_right_diff_val}, run={iteration_num}): Saved to {txt_save_path}")
  340. del weight_matrix_gpu
  341. try:
  342. cp.get_default_memory_pool().free_all_blocks()
  343. except:
  344. pass
  347. def run_purkinje_layer(args):
  348. left_right_diff_val, iteration_num, gGABA_val = args
  350. start_scope()
  352. output_dir = os.path.join(base_data_dir, f'{gGABA_val}nS', f'{left_right_diff_val}')
  353. os.makedirs(output_dir, exist_ok=True)
  355. txt_load_path = os.path.join(base_data_dir, f'{left_right_diff_val}', f'E_spikes_run_{iteration_num}.txt')
  357. if not os.path.exists(txt_load_path):
  358. print(f"Error: Continuous model data not found at {txt_load_path}")
  359. return
  361. np.random.seed(iteration_num)
  362. params = continuous_model_params.copy()
  363. gc_par_local = gc_par.copy()
  364. PC_par_local = PC_par.copy()
  365. sim_params_local = sim_params.copy()
  367. PC_par_local['gGABA'] = gGABA_val * nS
  369. PC_par_local['gAMPA'] = PC_par_local['gAMPA'] / gc_par_local['N_gc']
  371. try:
  372. loaded_data = np.loadtxt(txt_load_path)
  373. if loaded_data.size == 0:
  374. E_replay_indices = np.array([], dtype=int)
  375. E_replay_times = np.array([]) * ms
  376. elif loaded_data.ndim == 1:
  377. E_replay_indices = np.array([int(loaded_data[0])])
  378. E_replay_times = np.array([loaded_data[1]]) * ms
  379. else:
  380. E_replay_indices = loaded_data[:, 0].astype(int)
  381. E_replay_times = loaded_data[:, 1] * ms
  382. except Exception as e:
  383. print(f"Error reading file {txt_load_path}: {e}")
  384. return
  386. net = OrderedDict()
  388. net['E_replay'] = SpikeGeneratorGroup(2 * params['N_E'], E_replay_indices, E_replay_times)
  389. net['E'] = net['E_replay']
  391. net['pgr'] = PoissonGroup(gc_par_local['N_gc'], rates=gc_par_local['possion_noise'])
  393. net['gc_layer'] = NeuronGroup(gc_par_local['N_gc'], model=eqs_gc,threshold='v >= vpeak',reset='v = c_reset; u += d_reset',method='euler',namespace=gc_par_local)
  394. net['gc_layer'].v = gc_par_local['vr']
  395. net['gc_layer'].u = gc_par_local['b'] * (net['gc_layer'].v - gc_par_local['vr'])
  397. net['PC_neuron_left'] = NeuronGroup(PC_par_local['N'], eqs_pc, threshold='v >= VSpike', reset='v = vreset; w += b', method='euler',namespace=PC_par_local)
  398. net['PC_neuron_right'] = NeuronGroup(PC_par_local['N'], eqs_pc, threshold='v >= VSpike', reset='v = vreset; w += b', method='euler',namespace=PC_par_local)
  399. net['PC_neuron_left'].v = -55 * mV
  400. net['PC_neuron_left'].w = 0* pA
  401. net['PC_neuron_right'].v = -55 * mV
  402. net['PC_neuron_right'].w = 0* pA
  403. net['gc_left']=net['gc_layer'][0:gc_par_local['N_gc'] // 2]
  404. net['gc_right']=net['gc_layer'][gc_par_local['N_gc'] // 2:]
  406. net['icgr']=Synapses(net['pgr'],net['gc_layer'],on_pre="sAMPA_noise+=1",namespace=gc_par_local)
  407. net['icgr'].connect(condition='i == j')
  409. pre_indices = []
  410. post_indices = []
  411. E_left_global_indices = np.arange(params['N_E'])
  412. E_right_global_indices = np.arange(params['N_E'], 2 * params['N_E'])
  413. for l in range(gc_par_local['N_gc'] // 2):
  414. chosen_pre_indices = np.random.choice(E_left_global_indices, size=gc_par_local['mossy_fiber_input_num'], replace=False)
  415. pre_indices.extend(chosen_pre_indices)
  416. post_indices.extend([l] * gc_par_local['mossy_fiber_input_num'])
  417. for l in range(gc_par_local['N_gc'] // 2, gc_par_local['N_gc']):
  418. chosen_pre_indices = np.random.choice(E_right_global_indices, size=gc_par_local['mossy_fiber_input_num'], replace=False)
  419. pre_indices.extend(chosen_pre_indices)
  420. post_indices.extend([l] * gc_par_local['mossy_fiber_input_num'])
  422. net['icEgr']=Synapses(net['E'],net['gc_layer'],on_pre="sAMPA+=1",namespace=gc_par_local)
  423. net['icEgr'].connect(i=pre_indices, j=post_indices)
  425. net['icgrPc_left']=Synapses(net['gc_left'], net['PC_neuron_left'], on_pre='sAMPA_post += 1')
  426. net['icgrPc_left'].connect()
  427. net['icgrPc_right']=Synapses(net['gc_right'], net['PC_neuron_right'], on_pre='sAMPA_post += 1')
  428. net['icgrPc_right'].connect()
  430. net['icPcPc_left']=Synapses(net['PC_neuron_left'],net['PC_neuron_left'],on_pre='sGABA+=1')
  431. net['icPcPc_left'].connect('i==j')
  432. net['icPcPc_right']=Synapses(net['PC_neuron_right'],net['PC_neuron_right'],on_pre='sGABA+=1')
  433. net['icPcPc_right'].connect('i==j')
  435. @network_operation(when='start')
  436. def recurrent_input_pc():
  437. net['PC_neuron_left'].S_GABA=sum(net['PC_neuron_right'].sGABA)/PC_par_local['N']
  438. net['PC_neuron_right'].S_GABA=sum(net['PC_neuron_left'].sGABA)/PC_par_local['N']
  440. spikemon_PC_left = SpikeMonitor(net['PC_neuron_left'])
  441. spikemon_PC_right = SpikeMonitor(net['PC_neuron_right'])
  443. defaultclock.dt=0.05*ms
  444. total_duration = sim_params_local['sim_purations']
  446. network = Network(list(net.values()) + [recurrent_input_pc, spikemon_PC_left, spikemon_PC_right])
  448. print(f"Purkinje Simulation (gGABA={gGABA_val}nS, diff={left_right_diff_val}, run={iteration_num}): Starting run...")
  449. try:
  450. network.run(total_duration, report='text')
  451. print(f"Purkinje Simulation (gGABA={gGABA_val}nS, diff={left_right_diff_val}, run={iteration_num}): Finished run.")
  452. except Exception as e:
  453. print(f"Purkinje Simulation (gGABA={gGABA_val}nS, diff={left_right_diff_val}, run={iteration_num}): Failed with error: {e}")
  454. return
  456. PC_spike_times_left = spikemon_PC_left.t
  457. PC_spike_indices_left = spikemon_PC_left.i
  458. PC_spike_times_right = spikemon_PC_right.t
  459. PC_spike_indices_right = spikemon_PC_right.i
  461. np.save(os.path.join(output_dir, f'spikes_PC_left_t_run_{iteration_num}.npy'), PC_spike_times_left/ms)
  462. np.save(os.path.join(output_dir, f'spikes_PC_left_i_run_{iteration_num}.npy'), PC_spike_indices_left)
  463. np.save(os.path.join(output_dir, f'spikes_PC_right_t_run_{iteration_num}.npy'), PC_spike_times_right/ms)
  464. np.save(os.path.join(output_dir, f'spikes_PC_right_i_run_{iteration_num}.npy'), PC_spike_indices_right)
  466. window_size = 100*ms
  467. step_size = 2*ms
  468. time_points, rate_PC_left = calculate_sliding_window_rate(PC_spike_times_left, PC_par_local['N'], total_duration, window_size, step_size)
  469. _, rate_PC_right = calculate_sliding_window_rate(PC_spike_times_right, PC_par_local['N'], total_duration, window_size, step_size)
  471. np.save(os.path.join(output_dir, f'rates_time_points_run_{iteration_num}.npy'), time_points)
  472. np.save(os.path.join(output_dir, f'rates_PC_left_run_{iteration_num}.npy'), rate_PC_left)
  473. np.save(os.path.join(output_dir, f'rates_PC_right_run_{iteration_num}.npy'), rate_PC_right)
  475. if __name__ == '__main__':
  476. try:
  477. multiprocessing.set_start_method('spawn', force=True)
  478. except RuntimeError:
  479. pass
  481. # phase1_gpu_config = {
  482. # 1: 3,
  483. # 2: 3,
  484. # 3: 3,
  485. # 4: 3,
  486. # 5: 3
  487. # }
  488. phase1_gpu_config = {
  489. 0:3
  490. }
  491. phase1_num_processes_to_use = sum(list(phase1_gpu_config.values()))
  493. phase1_gpu_pool = []
  494. for gid, threads in phase1_gpu_config.items():
  495. phase1_gpu_pool.extend([gid] * threads)
  496. # -------------------------------------------------------------
  498. left_right_diff_values = [1,3,5,9,13]
  499. num_runs_per_diff = 255
  500. gGABA_values = range(2, 31)
  501. num_processes_to_use = 255
  502. # -------------------------------------------------------------
  504. print(f"Using {phase1_num_processes_to_use} processes for Phase 1 multiprocessing.")
  505. print(f"Using {num_processes_to_use} processes for Phase 2 multiprocessing.")
  507. for diff_val in left_right_diff_values:
  508. print(f"\n=== Phase 1: Running Continuous Model for diff = {diff_val} ===")
  509. phase1_tasks = []
  510. for run_idx in range(1, num_runs_per_diff + 1):
  512. expected_txt_path = os.path.join(base_data_dir, f'{diff_val}', f'E_spikes_run_{run_idx}.txt')
  513. if os.path.exists(expected_txt_path):
  514. print(f"File exists, skipping Phase 1 task submission for diff={diff_val}, run={run_idx}")
  515. continue
  517. assigned_gpu = phase1_gpu_pool[run_idx % len(phase1_gpu_pool)]
  518. phase1_tasks.append((diff_val, run_idx, assigned_gpu))
  520. if len(phase1_tasks) > 0:
  521. print(f"Submitting {len(phase1_tasks)} Phase 1 tasks...")
  522. with multiprocessing.Pool(processes=phase1_num_processes_to_use) as pool:
  523. for _ in pool.imap_unordered(run_continuous_layer, phase1_tasks):
  524. pass
  525. else:
  526. print("All Phase 1 tasks already completed. Skipping simulation.")
  528. for g_val in gGABA_values:
  529. print(f"\n=== Processing gGABA = {g_val}nS, left_right_diff = {diff_val} ===")
  531. print(f"--- Phase 2: Running Purkinje Model (Multi Threaded) ---")
  532. current_tasks = []
  533. for run_idx in range(1, num_runs_per_diff + 1):
  534. current_tasks.append((diff_val, run_idx, g_val))
  536. print(f"Submitting {len(current_tasks)} tasks for gGABA = {g_val}nS...")
  537. with multiprocessing.Pool(processes=num_processes_to_use) as pool:
  538. for _ in pool.imap_unordered(run_purkinje_layer, current_tasks):
  539. pass
  540. print(f"--- Finished gGABA = {g_val}nS ---")
  542. print(f"\nAll simulations completed across all gGABA and diff_values.")

cortico-cerebellar.py at commit e03ac74, no license · at the source

Overview

Authors: Yeyao Bao1,2, Liao Yu3, Liangfu Lu1,2, Zhuoqin Yang3, Yunliang Zang1,2
ORCID iDs: Yeyao Bao, Yunliang Zang
  1. Academy of Medical Engineering and Translational Medicine, Medical Faculty, Tianjin University, Tianjin 300072, China
  2. State Key Laboratory of Advanced Medical Materials and Medical Devices, Tianjin University, Tianjin 300072, China
  3. School of Mathematical Sciences, Beihang University, Beijing 100191, China
Institutions: Tianjin University (China); Beihang University (China)
Journal: Proceedings of the National Academy of Sciences of the United States of America, volume 123, issue 35, article e2616911123
Dates: received 13 May 2026; accepted 22 July 2026; published online 27 August 2026; in print 1 September 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1073/pnas.2616911123 · PMID 42658769 · PMCID PMC13534631 · OpenAlex W7204485119
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: computational modeling (no new data) (modality), human (organism), cognitive (subfield)
Methods: Smoothing, state filtering, decompositions, Single-unit activity, calcium imaging, Connectivity
Keywords: cerebellar learning, decision-making, feedforward network, primacy bias, pattern separation
MeSH: Cerebellar Cortex*, Cerebellum*, Decision Making*, Models, Neurological*, Animals, Computer Simulation, Humans, Nerve Net, Purkinje Cells (* major topic)
Topic: Vestibular and auditory disorders (Neurology, Neuroscience), according to OpenAlex
Funding: The National Key Research and Development Program of China (2023YFF1204200); the National Natural Science Foundation of China (12372060, 62476197)
Citations: cited by 1 paper (Europe PMC); 88 references in the paper

Abstract

The cerebellum is increasingly implicated in perceptual decision-making, yet it remains unclear how the cerebellar cortex can support sparse evidence accumulation over behavioral time scales without assuming cortical-style dense local excitatory recurrence. We present a biologically constrained modeling framework showing that cerebellar microcircuits can implement graded accumulation and competition in the absence of such recurrence. In our model, type-II Purkinje-neuron excitability generates firing-rate hysteresis that prolongs the impact of brief inputs far beyond intrinsic membrane and synaptic time constants, enabling accumulation across long interevent intervals. Purkinje neuron collateral inhibition produces competitive divergence and tunes temporal evidence weighting, revealing a trade-off between commitment and primacy bias. In a bidirectionally coupled cortico-cerebello-cortical model, cerebellar processing reduces primacy while cortical processing reduces indecision, improving robustness. Finally, granule-layer sparsification improves the separability of correlated inputs, enhancing discrimination under biologically realistic stimulus statistics. Together, these simulation results propose a mechanistic division of labor that positions the cerebellum as an active computational partner in perceptual decisions.

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

Repository

Its files are read in the Code ↔ Paper reader above.

yyb177/cerebellar-decision-making

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: e03ac74ba0c66ccb4bb275d1b81b8defd1122c7d, 13 May 2026
Languages: Python (3)
Size: 4 files, 3 scripts
Software Heritage: not archived
Found in: the references
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Brian 2 (3 files), NumPy (3 files), CuPy (2 files), Matplotlib (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
4 files
At the source: github.com/yyb177/cerebellar-decision-making

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;
  • 3 scripts, each with its path and the digest of its content;
  • no match between paragraphs and code yet;
  • 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, Materials, and Software Availability

The simulation code used in this study is publicly available at GitHub (88).

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

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 5 keywords, 9 MeSH terms, 2 funders, 87 references.

Cite

This paper

Bao, Y., Yu, L., Lu, L., Yang, Z., & Zang, Y. (2026). Cerebellar microcircuits enable robust evidence-based decisions through cortico-cerebellar coupling. Proceedings of the National Academy of Sciences of the United States of America, 123(35), e2616911123. https://doi.org/10.1073/pnas.2616911123

BibTeX

@article{bao2026cerebellar,
author = {Bao, Yeyao and Yu, Liao and Lu, Liangfu and Yang, Zhuoqin and Zang, Yunliang},
title = {{Cerebellar microcircuits enable robust evidence-based decisions through cortico-cerebellar coupling}},
journal = {Proceedings of the National Academy of Sciences of the United States of America},
year = {2026},
month = aug,
volume = {123},
number = {35},
pages = {e2616911123},
publisher = {National Academy of Sciences},
issn = {0027-8424},
doi = {10.1073/pnas.2616911123},
url = {https://doi.org/10.1073/pnas.2616911123},
pmid = {42658769},
pmcid = {PMC13534631}
}

RIS

TY - JOUR
AU - Bao, Yeyao
AU - Yu, Liao
AU - Lu, Liangfu
AU - Yang, Zhuoqin
AU - Zang, Yunliang
TI - Cerebellar microcircuits enable robust evidence-based decisions through cortico-cerebellar coupling
T2 - Proceedings of the National Academy of Sciences of the United States of America
J2 - Proc Natl Acad Sci U S A
PY - 2026
DA - 2026/08/27
VL - 123
IS - 35
SP - e2616911123
SN - 0027-8424
PB - National Academy of Sciences
DO - 10.1073/pnas.2616911123
UR - https://doi.org/10.1073/pnas.2616911123
LA - en
ER -

CSL-JSON

{
"id": "10.1073/pnas.2616911123",
"type": "article-journal",
"title": "Cerebellar microcircuits enable robust evidence-based decisions through cortico-cerebellar coupling",
"container-title": "Proceedings of the National Academy of Sciences of the United States of America",
"author": [
{
"family": "Bao",
"given": "Yeyao"
},
{
"family": "Yu",
"given": "Liao"
},
{
"family": "Lu",
"given": "Liangfu"
},
{
"family": "Yang",
"given": "Zhuoqin"
},
{
"family": "Zang",
"given": "Yunliang"
}
],
"container-title-short": "Proc Natl Acad Sci U S A",
"volume": "123",
"issue": "35",
"page": "e2616911123",
"DOI": "10.1073/pnas.2616911123",
"PMID": "42658769",
"PMCID": "PMC13534631",
"ISSN": "0027-8424",
"publisher": "National Academy of Sciences",
"URL": "https://doi.org/10.1073/pnas.2616911123",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
27
]
]
}
}

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.1014595 [code]
Cerebellum-inspired neural network of supervised learning with tensor-based sparse coding for multi-class classification.
Journal: PLoS computational biology
In common: Matplotlib, NumPy, 6 references
[2] doi:10.1038/s41593-026-02255-7 [code]
Neural circuits encode prior knowledge of temporal statistics.
Journal: Nature neuroscience
In common: Matplotlib, NumPy, 6 references
[3] doi:10.1016/j.celrep.2026.117262 [code]
Cerebellar acceleration of learning in an evidence-accumulation task.
Journal: Cell reports
In common: Matplotlib, NumPy, cognitive, 5 references
[4] doi:10.1126/sciadv.aed6417 [code]
Intrinsic timing, not temporal prediction, underlies ramping dynamics in visual and parietal cortex during passive behavior.
Journal: Science advances
In common: Matplotlib, NumPy, 5 references
[5] doi:10.1038/s41467-026-75924-7 [code]
Data-driven reduced modeling of neural dynamics.
Journal: Nature communications
In common: Matplotlib, NumPy, computational modeling (no new data), 4 references
[6] doi:10.1038/s41467-026-75371-4 [code]
Cerebellar neural populations orchestrate dopamine reward signaling with single-trial precision.
Journal: Nature communications
In common: Matplotlib, NumPy, 4 references
[7] doi:10.1038/s41467-026-75704-3 [code]
A minimal model of working memory in neural systems and neuromorphic circuits.
Journal: Nature communications
In common: Brian 2, Matplotlib, NumPy, computational modeling (no new data), 1 reference
[8] doi:10.7554/elife.102813
Multiple functions of cerebello-thalamic neurons in learning and offline consolidation of a motor skill in mice.
Journal: eLife
In common: 4 references
[9] doi:10.1038/s42003-026-10034-0 [code]
Region- and cell type-specific changes in gene expression in the cerebellum after classical fear conditioning.
Journal: Communications biology
In common: Matplotlib, NumPy, 3 references
[10] doi:10.1038/s41593-026-02388-9 [code]
Hippocampal CA3 connectomics reveals a gradient of mossy fiber inputs and selective feedforward inhibition onto pyramidal cells.
Journal: Nature neuroscience
In common: Matplotlib, NumPy, 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.