OSCR

Spiking neurons as predictive controllers of linear systems.

Code ↔ Paper

8 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 8 matches
  1. [1] § Results › Reactive vs. predictive spiking control ↔ lib_plot.py, lines 520–557 · score 0.71 · green arrow, red arrows, reactive spiking, predictive spiking, rightward, upward
  2. [2] § Results › Arms reaching task simulation ↔ fig9.ipynb, lines 239–282 · score 0.67 · fiber stretches, SMD responses, system matrix, damping, stiffness, Figure 9
  3. [3] § Materials and methods › Filtered spiking controller ↔ fig7,11.ipynb, lines 24–114 · score 0.66 · optimal feedback gain, control actuation, state deviations, Spiking activity, Filtered, LQR
  4. [4] § Results › Controlling coupled oscillators ↔ fig9.ipynb, lines 1–21 · score 0.66 · arm movement, neuromuscular units, point mass, arm position, point arm, space
  5. [5] § Materials and methods › Simulation and network setup › Continuous and filtered-spiking parameters. ↔ lib_sim.py, lines 261–361 · score 0.61 · decay factor, filtered spiking control, spiking cost, control signal, LQR, maps
  6. [6] § Results › Reactive vs. predictive spiking control ↔ lib_plot.py, lines 520–557 · score 0.59 · green arrow, red arrows, reactive spiking, rightward, upward, trajectory
  7. [7] § Materials and methods › Simulation and network setup › Target settings. ↔ fig9.ipynb, lines 92–163 · score 0.52 · curved arch, angles, arm, Figure 9, zero, trajectories
  8. [8] § Results › Controlling coupled oscillators ↔ fig9.ipynb, lines 24–89 · score 0.51 · coupled oscillator system, point arm, fiber, mass, connected, space

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

Jupyter notebook · 521 lines · 16 KB · CC-BY-4.0 · 4 matches

  1. # %% [markdown]
  2. # # Figure 9: Neuromuscular units control
  3. # %% [markdown]
  4. # In this Notebook, we want to showcase the activity of a bigger network of our spiking neurons that controls a a number of neuromuscular units (modeled as linear SMDs) which in turn affect the position of a 2D point arm in space. Interestingly, the target here is only set on the 2D arm position, and not on any state of the single neuromuscular units. The direct spiking control is actuated from our network onto the units, the state of which are in turn are coupled with the position of the arm. NOTE: the 2D arm is simply modeled as two point masses moving in space without any nonlinearity. This is of course not a realistic model of arm movement, but expanding the paradigm to nonlinear system is still in progress.
  5. # %% [markdown]
  6. # ### Imports
  7. #
  8. # We import numpy, scipy.linalg, pyplot, an rng from np.random, and our two libraries for simulation and plotting.
  9. # %%
  10. import numpy as np
  11. import scipy.linalg as la
  12. import matplotlib.pyplot as plt
  13. rng = np.random.default_rng()
  14. import lib_sim
  15. import lib_plot
  16. from importlib import reload
  17. reload(lib_sim)
  18. reload(lib_plot)
  19. # %% [markdown]
  20. # ### Initializations
  21. # We initialize all the elements we need to run the simulations.
  22. #
  23. #
  24. # 1. Number of Neurons *N*
  25. # 2. Number of coupled neuromuscular units (8) + the coordinates of the point arm they control (2) *M=10*
  26. # 2. **A** matrix. For a coupled oscillator system our state martix **A** will be slightly more complicated. We first fill the elements corresponding to the dynamics of each single unit, in a way that scales with our *M-1* units. Then, we fill in the cross terms that represent the effect of the state of a unit onto the states of the ones connected to it. Finallt, we use the value *gamma* to represent the continuous cross terms which allow the state of the units to be associated with a point mass (arm) which does not have the same oscillatory dynamics as the SMDs, but rather can be linearly controlled to move freely in space.
  27. # 3. **B** matrix. In this case, we have to scale **B** according to the neurons *N* and the units *M*. We initialize it with random values and then we normalize those values by columns, so that each neuron has a comparable effect onto the system. No connection is present between our neurons and the point arm, so the first columns will be zeros.
  28. # 4. **C** matrix. Cost martrix for state deviations. Also scaled based on *M*
  29. # %%
  30. N = 500
  31. M = 10 #Two arm coordinates and 8 MUs
  32. gamma_scale = 1
  33. A = np.zeros((2*M, 2*M))
  34. # Position equations
  35. for i in range(M):
  36. A[2*i, 2*i+1] = 1
  37. # Remove self dynamics for first two masses (arm)
  38. # (skip spring/damping for i=0,1)
  39. # Normal SMDs for fibers (i >= 2)
  40. for i in range(2, M):
  41. A[2*i+1, 2*i] = -12/6
  42. A[2*i+1, 2*i+1] = -3.8/6
  43. # --- All fibers influence arm ---
  44. for j in range(2, M):
  45. pos_index = 2*j
  46. # alternating sign
  47. sign = (-1)**(j)
  48. A[1, pos_index] = -gamma_scale * sign
  49. A[3, pos_index] = gamma_scale * sign
  50. B = np.zeros((2*M, N))
  51. # Fill only even rows with random values
  52. for i in range(1, 2 * M, 2):
  53. B[i] = np.random.randn(N)
  54. # Normalize each column
  55. B /= np.linalg.norm(B, axis=0)
  56. for i in range(1, 2 * M, 2):
  57. B[i] *= 6
  58. B[:3, :] = 0
  59. C = np.zeros((2*M, 2*M))
  60. # C[0, 0] = 1
  61. # C[2, 2] = 1
  62. for i in range(0, 2*M, 2):
  63. C[i, i] = 1
  64. C[0:, 0:] *= 0.001
  65. C[0, 0] = 1
  66. C[2, 2] = 1
  67. # %% [markdown]
  68. # We also need to initialize all the other elements of the simulation:
  69. #
  70. # 6. Time settings:
  71. # - Number of seconds (T) timestep lenght (dt), time array, number of timesteps (nT)
  72. # - Future time window $f$, here called *t_fut*
  73. # 7. Target settings:
  74. # - Here we want to set a different target trajectory for each coordinate of the arm. *angles* and *targets* store the eight angles and endpoints of the trajectory, while *z_base* is compiled to trace the trajectory. In the simulation, the target will iterate on the values of *z_base* over time.
  75. # 8. Threshold parameters:
  76. # - *mu*, spiking cost on the threshold (see math derivation)
  77. # - If input delays from the system to the network are not considered, we simulate with asyncronous firing (one neuron is selected to spikes at a time). Optionally (not done in the paper), we can take delays into consideration, which means allowing multiple neurons to fire in each timestep. To prevent overspiking in this case, we adjust the threshold based on previous spiking activity. The values of *a* and *lambda* scale this adjustments
  78. # - refractory period *ref_period_lenght* is set to zero. The algorithm selects one neuron to fire in each timestep but does not impose any non-firing period. Since neurons fire to improve a future loss, they already tend to not overspike. We kept the option, so that the refractory period can be introduced if needed (for instance if delays are present)
  79. # 9. Cell silencing parameters (optional):
  80. # - *cell_death_timings* contains the timesteps in which we want to silence neurons
  81. # - *k_per_step* manages how many neuron to silence in each of the silencing timetsep
  82. # %%
  83. # Time settings
  84. T = 100
  85. dt = 1e-2
  86. times = np.arange(0, T, dt)
  87. nT = len(times)
  88. t_fut = 2
  89. # Target parameters
  90. leak_z = 1
  91. z = np.zeros((2*M, nT))
  92. start_idx = 1000
  93. end_idx = 10000
  94. Tmove = end_idx - start_idx
  95. tau = np.linspace(0, 1, Tmove)
  96. w = 10*tau**3 - 15*tau**4 + 6*tau**5
  97. R = 20
  98. angles = np.linspace(0, 2*np.pi, 8, endpoint=False)
  99. targets = [(R*np.cos(a), R*np.sin(a)) for a in angles]
  100. z_bases = []
  101. for (wx1, wy1) in targets:
  102. z_base = np.zeros((2*M, nT))
  103. wx0, wy0 = 0, 0
  104. # midpoint for curved arch
  105. wxm = (wx0 + wx1)/2
  106. wym = (wy0 + wy1)/2 + 10 # vertical lift (constant arch height)
  107. for k in range(Tmove):
  108. wk = w[k]
  109. idx = start_idx + k
  110. z_base[0, idx] = (1-wk)**2 * wx0 + 2*(1-wk)*wk * wxm + wk**2 * wx1
  111. z_base[2, idx] = (1-wk)**2 * wy0 + 2*(1-wk)*wk * wym + wk**2 * wy1
  112. z_bases.append(z_base.copy())
  113. # theta = 0.2
  114. x = np.zeros((2*M, nT))
  115. x0 = np.zeros((2*M))
  116. #Threshold/firing parameters
  117. lam = 0.3 #0.25
  118. a = 0
  119. mu = 0.3
  120. ref_period_lenght = 0
  121. #Cell silencing parameters
  122. cell_death_timings = [3000, 7000] # Timesteps to kill neurons
  123. k_per_step = 180 # Kill 180 neurons at each step
  124. # %% [markdown]
  125. # ### Simulation
  126. #
  127. # We make use of the spiking control simulation function from lib_sim to simulate control of the SMDs and consequently the arm.
  128. # We repeat the experiment 8 times for the eight different reaching trajectories
  129. #
  130. # We extract:
  131. # - target z,
  132. # - spikes s,
  133. # - state x,
  134. # - voltages V,
  135. # - thresholds Th,
  136. # - predicted state x_pred,
  137. # - error E,
  138. # - predicted error pE.
  139. # - killed neurons indexes
  140. # %%
  141. x_list = []
  142. z_list = []
  143. for k in range(8):
  144. if k == 0:
  145. zk, sk, xk, Vk, Thk, x_predk, Ek, pEk, killed = lib_sim.simulate_coupled_spiking_control(nT, dt, t_fut, M, N, A, B, C, z.copy(), x.copy(), x0.copy(), lam, a, mu, ref_period_lenght, k_per_step, cell_death_timings, delay=False, Zstep=True, kill=False, perturb=True, musc=True, z_base=z_bases[k], leak_z=leak_z)
  146. else:
  147. zk, sk, xk, Vk, Thk, x_predk, Ek, pEk, killed = lib_sim.simulate_coupled_spiking_control(nT, dt, t_fut, M, N, A, B, C, z.copy(), x.copy(), x0.copy(), lam, a, mu, ref_period_lenght, k_per_step, cell_death_timings, delay=False, Zstep=True, kill=False, perturb=False, musc=True, z_base=z_bases[k], leak_z=leak_z)
  148. x_list.append(xk.copy())
  149. z_list.append(zk.copy())
  150. x_list.reverse()
  151. z_list.reverse()
  152. # %% [markdown]
  153. # ### Plotting Functions
  154. #
  155. # Since the figure is quite complex and requires external illustrations, we will plot all the axs separately and combined them in a separate time.
  156. # First, we need to define two ad-hoc function to
  157. # - Plot the trajectories of the point arm in 2D space
  158. # - Showcase the example response of one SMD to one spike input (for illustration)
  159. # %%
  160. def plot_arm_reach(x, z, hide=False, ax=None):
  161. # --- Target trajectory ---
  162. ax.plot(z[0, :], z[2, :],
  163. color='gray',
  164. linewidth=1.5,
  165. label='Target trajectory')
  166. # --- Actual arm trajectory ---
  167. ax.plot(x[0, :], x[2, :],
  168. color='indianred',
  169. linewidth=1.5,
  170. label='Arm trajectory')
  171. ax.scatter(x[0,0], x[2,0], color='green', s=80, label='Start')
  172. ax.scatter(x[0,-1], x[2,-1], color='red', s=80, label='End')
  173. # --- Labels ---
  174. ax.set_xlabel(r'$x_1$')
  175. ax.set_ylabel(r'$x_2$')
  176. #ax.legend()
  177. # ax.set_yticks(np.linspace(0, 20, 3))
  178. # ax.set_xticks(np.linspace(0, 20, 3))
  179. if hide:
  180. ax.spines['top'].set_visible(False)
  181. ax.spines['right'].set_visible(False)
  182. ax.set_yticks(np.linspace(-20, 20, 3))
  183. ax.set_xticks(np.linspace(-20, 20, 3))
  184. #ax.spines['bottom'].set_visible(False)
  185. #ax.tick_params(bottom=False, labelbottom=False)
  186. # %%
  187. dte = 0.01
  188. Te = 20
  189. nTe = int(Te/dte)
  190. se = np.zeros(nTe)
  191. def plot_single_smd_response(s, nT, T, dte, hide=True, ax=None):
  192. # --- Parameters ---
  193. k = 12/6 # stiffness
  194. d = 3.8/6 # damping
  195. # --- System matrices ---
  196. A = np.array([[0, 1],
  197. [-k, -d]])
  198. B = np.array([[0],
  199. [80]])
  200. # --- State ---
  201. x = np.zeros((2, nT))
  202. # --- Input (spike at specific time) ---
  203. spike_time = 2.0
  204. spike_index = int(spike_time / dte)
  205. s[spike_index] = 1 # spike amplitude
  206. # --- Simulate ---
  207. for t in range(nT-1):
  208. x[:, t+1] = x[:, t] + dt * (A @ x[:, t] + B.flatten() * s[t])
  209. # --- Plot ---
  210. time = np.linspace(0, T, nT)
  211. ax.plot(time, x[0, :], label='Fiber stretch', color='purple', alpha=0.7, linewidth=2)
  212. ax.plot(time, s[:], label='Spike Input', color='orange', linestyle='-.')
  213. if hide:
  214. ax.axis('off')
  215. ax.spines['top'].set_visible(False)
  216. ax.spines['right'].set_visible(False)
  217. ax.spines['bottom'].set_visible(False)
  218. ax.tick_params(bottom=False, labelbottom=False)
  219. # Create a floating legend
  220. ax.legend(loc='upper right', frameon=False, markerscale=1.5, fontsize='large', fancybox=False, framealpha=0.7, bbox_to_anchor=(1.01, 0.8))
  221. # %% [markdown]
  222. # ### Plotting
  223. #
  224. # Now we use the functions we created along wit the ones in our lib_plot library to plot all the individual axes of the figure.
  225. # %%
  226. plot_mosaic = [['B']
  227. ]
  228. #Generate figure and axes
  229. fig, axs = plt.subplot_mosaic(
  230. plot_mosaic,
  231. #layout='constrained',
  232. #empty_sentinel=None,
  233. figsize=(5, 2),
  234. dpi=300,
  235. gridspec_kw={'height_ratios': [1], 'width_ratios': [1]}
  236. )
  237. plot_single_smd_response(se, nTe, Te, dte, hide=True, ax=axs['B'])
  238. #fig.set_constrained_layout_pads(hspace=0.01)
  239. # for label, ax in axs.items():
  240. # if label != '.': # Skip spacers
  241. # ax.plot()
  242. # ax.set_xticks([])
  243. # ax.set_yticks([])
  244. # lib_plot.label_axes(axs, fig)
  245. plt.plot()
  246. plt.savefig('figs/fig9/fibersketch.svg', dpi=300)
  247. # %%
  248. plot_mosaic = [['B']
  249. ]
  250. #Generate figure and axes
  251. fig, axs = plt.subplot_mosaic(
  252. plot_mosaic,
  253. #layout='constrained',
  254. #empty_sentinel=None,
  255. figsize=(3.5, 3.5),
  256. dpi=300,
  257. gridspec_kw={'height_ratios': [1], 'width_ratios': [1]}
  258. )
  259. fig.set_constrained_layout_pads(hspace=0.01)
  260. for label, ax in axs.items():
  261. if label != '.': # Skip spacers
  262. ax.plot()
  263. ax.set_xticks([])
  264. ax.set_yticks([])
  265. lib_plot.label_axes(axs, fig)
  266. #axsB
  267. for k in range(0, 8):
  268. plot_arm_reach(x_list[k], z_list[k],
  269. hide=True,
  270. ax= axs['B'])
  271. handles, labels = axs['B'].get_legend_handles_labels()
  272. by_label = dict(zip(labels, handles))
  273. axs['B'].legend(by_label.values(), by_label.keys(),
  274. loc='upper right', frameon=False, markerscale=0.2,
  275. fontsize='x-small', fancybox=False, framealpha=0.7,
  276. bbox_to_anchor=(1.4, 0.3))
  277. axs['B'].set_aspect('equal')
  278. # plot_arm_reach(x, z, hide=True, ax=axs['B'])
  279. plt.plot()
  280. plt.savefig('figs/fig9/axBfig9.svg', dpi=300)
  281. plt.show()
  282. # %%
  283. #AxsC
  284. plot_mosaic = [['C'],
  285. ['Ci'],
  286. ['Cii'],
  287. ]
  288. #Generate figure and axes
  289. fig, axs = plt.subplot_mosaic(
  290. plot_mosaic,
  291. #layout='constrained',
  292. #empty_sentinel=None,
  293. figsize=(4, 5),
  294. dpi=300,
  295. gridspec_kw={'height_ratios': [1, 1, 1], 'width_ratios': [1]}
  296. )
  297. fig.set_constrained_layout_pads(hspace=0.01)
  298. for label, ax in axs.items():
  299. if label != '.': # Skip spacers
  300. ax.plot()
  301. ax.set_xticks([])
  302. ax.set_yticks([])
  303. lib_plot.label_axes(axs, fig)
  304. x_es_arm = x_list[6][:3, :]
  305. z_es_arm = z_list[6][:3, :]
  306. x_es_mus = x_list[6][4:, :]
  307. z_es_mus = z_list[6][3:, :]
  308. lib_plot.plot_coupled_activity(x_es_arm, z_es_arm, Ek, times, xpred=None, incxpred=False, musc=False, timelabels=False, hide=True, ax=axs['C'])
  309. lib_plot.plot_coupled_activity(x_es_mus, z_es_mus, Ek*0, times, xpred=None, incxpred=False, musc=True, timelabels=False, hide=True, ax=axs['Ci'])
  310. lib_plot.coupled_raster_plot(sk, times, slence=True, hide=True, ax=axs['Cii'], cell_death_timings=cell_death_timings, k_per_step=k_per_step, N=N)
  311. plt.plot()
  312. plt.savefig('figs/fig9/axCfig9.svg', dpi=300)
  313. plt.show()
  314. # %%
  315. # plot_mosaic = [['A']
  316. # ]
  317. # #Generate figure and axes
  318. # fig, axs = plt.subplot_mosaic(
  319. # plot_mosaic,
  320. # #layout='constrained',
  321. # empty_sentinel=None
  322. # figsize=(5, 5),
  323. # dpi=300,
  324. # gridspec_kw={'height_ratios': [1], 'width_ratios': [1]}
  325. # )
  326. # fig.set_constrained_layout_pads(hspace=0.01)
  327. # for label, ax in axs.items():
  328. # if label != '.': # Skip spacers
  329. # ax.plot()
  330. # ax.set_xticks([])
  331. # ax.set_yticks([])
  332. # lib_plot.label_axes(axs, fig)
  333. # plt.plot()
  334. # plt.savefig('figs/fig11/axAfig11.svg', dpi=300)
  335. # %% [markdown]
  336. # ## Extra: Perturbation on the arms trajectory
  337. #
  338. # We now run only the simulation for the trajctory for which perturb the velocity of the arm in the x and y coordinates as the tagret is being traced
  339. # %%
  340. zk, sk, xk, Vk, Thk, x_predk, Ek, pEk, killed = lib_sim.simulate_coupled_spiking_control(nT, dt, t_fut, M, N, A, B, C, z.copy(), x.copy(), x0.copy(), lam, a, mu, ref_period_lenght, k_per_step, cell_death_timings, delay=False, Zstep=True, kill=False, perturb=True, musc=True, z_base=z_bases[0], leak_z=leak_z)
  341. # %%
  342. plot_mosaic = [['B']
  343. ]
  344. #Generate figure and axes
  345. fig, axs = plt.subplot_mosaic(
  346. plot_mosaic,
  347. #layout='constrained',
  348. #empty_sentinel=None,
  349. figsize=(4, 3),
  350. dpi=300,
  351. gridspec_kw={'height_ratios': [1], 'width_ratios': [1]}
  352. )
  353. fig.set_constrained_layout_pads(hspace=0.01)
  354. for label, ax in axs.items():
  355. if label != '.': # Skip spacers
  356. ax.plot()
  357. ax.set_xticks([])
  358. ax.set_yticks([])
  359. lib_plot.label_axes(axs, fig)
  360. #axsB
  361. plot_arm_reach(xk, zk,
  362. hide=True,
  363. ax= axs['B'])
  364. handles, labels = axs['B'].get_legend_handles_labels()
  365. by_label = dict(zip(labels, handles))
  366. axs['B'].legend(by_label.values(), by_label.keys(),
  367. loc='upper right', frameon=False, markerscale=0.2,
  368. fontsize='x-small', fancybox=False, framealpha=0.7,
  369. bbox_to_anchor=(1.4, 0.3))
  370. axs['B'].set_aspect('equal')
  371. # plot_arm_reach(x, z, hide=True, ax=axs['B'])
  372. plt.plot()
  373. plt.savefig('figs/fig9/perturbed_ex_fig9.svg', dpi=300)
  374. plt.show()
  375. # %%
  376. #AxsC
  377. plot_mosaic = [['C'],
  378. ['Ci'],
  379. ['Cii'],
  380. ]
  381. #Generate figure and axes
  382. fig, axs = plt.subplot_mosaic(
  383. plot_mosaic,
  384. #layout='constrained',
  385. #empty_sentinel=None,
  386. figsize=(4, 5),
  387. dpi=300,
  388. gridspec_kw={'height_ratios': [1, 1, 1], 'width_ratios': [1]}
  389. )
  390. fig.set_constrained_layout_pads(hspace=0.01)
  391. for label, ax in axs.items():
  392. if label != '.': # Skip spacers
  393. ax.plot()
  394. ax.set_xticks([])
  395. ax.set_yticks([])
  396. lib_plot.label_axes(axs, fig)
  397. x_es_arm = xk[:3, :]
  398. z_es_arm = zk[:3, :]
  399. x_es_mus = xk[4:, :]
  400. z_es_mus = zk[3:, :]
  401. lib_plot.plot_coupled_activity(x_es_arm, z_es_arm, Ek, times, xpred=None, incxpred=False, musc=False, timelabels=False, hide=True, ax=axs['C'])
  402. lib_plot.plot_coupled_activity(x_es_mus, z_es_mus, Ek*0, times, xpred=None, incxpred=False, musc=True, timelabels=False, hide=True, ax=axs['Ci'])
  403. lib_plot.coupled_raster_plot(sk, times, slence=True, hide=True, ax=axs['Cii'], cell_death_timings=cell_death_timings, k_per_step=k_per_step, N=N)
  404. plt.plot()
  405. plt.savefig('figs/fig9/perturbed_ex_act_fig9.svg', dpi=300)
  406. plt.show()
  407. # %%

fig9.ipynb at commit ba2f4fc, under CC-BY-4.0 · at the source

Overview

Authors: Paolo Agliati1, André Urbano2, Pablo Lanillos2, Nasir Ahmad1, Marcel van Gerven1, Sander Keemink1
  1. Department of Machine Learning and Neural Computing, Donders Institute for Brain, Cognition and Behaviour, Radboud University, Nijmegen, The Netherlands
  2. Neuro AI and Robotics group, Cajal International Neuroscience Center, Spanish National Research Council, Madrid, Spain
Journal: PLoS computational biology, volume 22, issue 7, article e1014432
Dates: received 20 October 2025; accepted 11 June 2026; published online 9 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pcbi.1014432 · PMID 42424433 · PMCID PMC13384403 · OpenAlex W4414422982
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: computational modeling (no new data) (modality), human (organism), computational (subfield)
Methods: Single-unit activity, calcium imaging, Machine learning
MeSH: Action Potentials*, Models, Neurological*, Neurons*, Animals, Computational Biology, Computer Simulation, Humans, Linear Models, Nerve Net, Neural Networks, Computer (* major topic)
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Ministerie van Onderwijs, Cultuur en Wetenschap (024005022)
Citations: not cited yet (Europe PMC); 51 references in the paper

Abstract

Neurons communicate with downstream systems via sparse and incredibly brief electrical pulses, or spikes. Using these events, they control various targets such as neuromuscular units, neurosecretory systems, and other neurons in connected circuits. This gave rise to the idea of spiking neurons as controllers, in which spikes are the control signal. Using instantaneous events directly as the control inputs, also called ‘impulsive control’, is challenging as it does not scale well to larger networks and has low analytical tractability. Therefore, current spiking control usually relies on filtering the spike signal to approximate analog control. This ultimately means spiking neural networks (SNNs) have to output a continuous control signal, necessitating continuous energy input into downstream systems. Here, we circumvent the need for rate-based representations, providing a scalable method for task-specific spiking control with sparse neural activity. In doing so, we take inspiration from both control theory and neuroscience, and define a spiking rule where spikes are only emitted if they bring a dynamical system closer to a target. From this principle, we derive the required connectivity for an SNN, and show that it can successfully control linear systems. We show that for physically constrained systems, predictive control is required, and the control signal ends up exploiting the passive dynamics of the downstream system to reach a target. Finally, we show that the paradigm scales to both high-dimensional systems and bio-inspired motor control tasks. Importantly, in all cases, we maintain a closed-form mathematical derivation of the network connectivity, the network dynamics and the control objective. This work advances the understanding of SNNs as biologically-inspired controllers, providing insight into how real neurons could exert control, and enabling applications in neuromorphic hardware design.

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

gitlab.socsci.ru.nl/2023-phd-paolo-agliati/spiking-neurons-as-predictive-controllers-of-linear-systems

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: ba2f4fca389ff483ac0e1c2ca3afe5c6acc2c951, 3 June 2026
Languages: Jupyter (8), Python (2)
Size: 33 files, 10 scripts
Software Heritage: not archived
Found in: “Code availability.”
Holds: README, 8 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (10 files), Matplotlib (9 files), SciPy (5 files), seaborn (2 files), Pillow (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
11 files

Code availability

All simulations were run with Python 3.11.5. The source code and all the necessary dependencies are available at https://gitlab.socsci.ru.nl/2023-phd-paolo-agliati/spiking-neurons-as-predictive-controllers-of-linear-systems.

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

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;
  • 10 scripts, each with its path and the digest of its content;
  • 8 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

There are no primary data in the paper; All relevant data are within the manuscript and its Supporting Information files. All simulations were run with Python 3.11.5. The source code and all the necessary dependencies are available at https://gitlab.socsci.ru.nl/2023-phd-paolo-agliati/spiking-neurons-as-predictive-controllers-of-linear-systems.

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, 6 authors, 10 MeSH terms, 1 funder, 45 references.

Cite

This paper

Agliati, P., Urbano, A., Lanillos, P., Ahmad, N., van Gerven, M., & Keemink, S. (2026). Spiking neurons as predictive controllers of linear systems. PLoS computational biology, 22(7), e1014432. https://doi.org/10.1371/journal.pcbi.1014432

BibTeX

@article{agliati2026spiking,
author = {Agliati, Paolo and Urbano, André and Lanillos, Pablo and Ahmad, Nasir and van Gerven, Marcel and Keemink, Sander},
title = {{Spiking neurons as predictive controllers of linear systems}},
journal = {PLoS computational biology},
year = {2026},
month = jul,
volume = {22},
number = {7},
pages = {e1014432},
publisher = {PLOS},
issn = {1553-734X},
doi = {10.1371/journal.pcbi.1014432},
url = {https://doi.org/10.1371/journal.pcbi.1014432},
pmid = {42424433},
pmcid = {PMC13384403}
}

RIS

TY - JOUR
AU - Agliati, Paolo
AU - Urbano, André
AU - Lanillos, Pablo
AU - Ahmad, Nasir
AU - van Gerven, Marcel
AU - Keemink, Sander
TI - Spiking neurons as predictive controllers of linear systems
T2 - PLoS computational biology
J2 - PLoS Comput Biol
PY - 2026
DA - 2026/07/09
VL - 22
IS - 7
SP - e1014432
SN - 1553-734X
PB - PLOS
DO - 10.1371/journal.pcbi.1014432
UR - https://doi.org/10.1371/journal.pcbi.1014432
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pcbi.1014432",
"type": "article-journal",
"title": "Spiking neurons as predictive controllers of linear systems",
"container-title": "PLoS computational biology",
"author": [
{
"family": "Agliati",
"given": "Paolo"
},
{
"family": "Urbano",
"given": "André"
},
{
"family": "Lanillos",
"given": "Pablo"
},
{
"family": "Ahmad",
"given": "Nasir"
},
{
"family": "van Gerven",
"given": "Marcel"
},
{
"family": "Keemink",
"given": "Sander"
}
],
"container-title-short": "PLoS Comput Biol",
"volume": "22",
"issue": "7",
"page": "e1014432",
"DOI": "10.1371/journal.pcbi.1014432",
"PMID": "42424433",
"PMCID": "PMC13384403",
"ISSN": "1553-734X",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pcbi.1014432",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
9
]
]
}
}

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.1014577 [code]
Combining sampling and attractor dynamics in spiking models of head direction systems.
Journal: PLoS computational biology
In common: SciPy, Matplotlib, NumPy, computational, computational modeling (no new data), 4 references
[2] doi:10.1371/journal.pcbi.1014337 [code]
Fast reconstruction of degenerate populations of conductance-based neuron models from spike times.
Journal: PLoS computational biology
In common: Pillow, seaborn, SciPy, 2 other tools, computational, computational modeling (no new data), 1 reference
[3] doi:10.1016/j.celrep.2026.117420 [code]
Neural population dynamics of direct electrical stimulation of neocortex.
Journal: Cell reports
In common: Pillow, seaborn, SciPy, 2 other tools, 1 reference
[4] doi:10.1038/s41467-026-75924-7 [code]
Data-driven reduced modeling of neural dynamics.
Journal: Nature communications
In common: SciPy, Matplotlib, NumPy, computational modeling (no new data), 2 references
[5] doi:10.1038/s41467-026-71961-4 [code]
Spatiotemporal asymmetries on brain energy landscape uncover system entrapment related to depression severity.
Journal: Nature communications
In common: Pillow, seaborn, SciPy, 2 other tools, computational
[6] doi:10.1038/s41467-026-74357-6 [code]
Hippocampo-neocortical interaction as compressive retrieval-augmented generation.
Journal: Nature communications
In common: Pillow, seaborn, SciPy, 2 other tools, computational modeling (no new data)
[7] doi:10.1038/s42003-026-10094-2 [code]
E/I imbalance and internal noise cause weak neural representations and face recognition challenges in ASD.
Journal: Communications biology
In common: Pillow, seaborn, SciPy, 2 other tools, computational modeling (no new data)
[8] doi:10.1007/s00422-026-01048-2 [code]
A quality measure for repeating multiple-unit spike patterns.
Journal: Biological cybernetics
In common: Matplotlib, NumPy, computational, 2 references
[9] doi:10.1016/j.isci.2026.116705 [code]
Spatiotemporal characteristics of visual cortical responses to transpalpebral electrical stimulation.
Journal: iScience
In common: Pillow, SciPy, Matplotlib, 1 other tool, computational, computational modeling (no new data)
[10] doi:10.3758/s13428-026-03122-w [code]
Modeling intraindividual variability in affect (MIVA): Formalized theoretical approach, computational model, and parameter recovery study.
Journal: Behavior research methods
In common: seaborn, SciPy, Matplotlib, 1 other tool, computational, computational modeling (no new data)

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.