OSCR

Fast reconstruction of degenerate populations of conductance-based neuron models from spike times.

Code ↔ Paper

14 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 14 matches
  1. [1] § Materials and methods › Conductance-based models ↔ FIGURES_GENERATION/figures_bg_bio/stg.py, lines 1141–1249 · score 0.85 · delayed rectifier potassium, stomatogastric ganglion, intracellular calcium, STG model, neuron model, sodium
  2. [2] § Materials and methods › Conductance-based models ↔ FIGURES_GENERATION/figures_methods/stg.py, lines 1141–1249 · score 0.85 · delayed rectifier potassium, stomatogastric ganglion, intracellular calcium, STG model, neuron model, sodium
  3. [3] § Materials and methods › Training procedure ↔ FIGURES_GENERATION/figures_methods/Alternative_Figure_APpendix.ipynb, lines 28–59 · score 0.70 · cosine annealing learning, rate schedule, restarts
  4. [4] § Results › DICs provide a structured and learnable representation of neuronal activity ↔ FIGURES_GENERATION/figures_bg_bio/stg.py, lines 1740–1848 · score 0.67 · calcium concentration, compensated conductances, sensitivity matrix, target DIC, linear, STG
  5. [5] § Results › DICs provide a structured and learnable representation of neuronal activity ↔ FIGURES_GENERATION/figures_methods/stg.py, lines 1740–1850 · score 0.67 · calcium concentration, compensated conductances, sensitivity matrix, target DIC, linear, STG
  6. [6] § Results › General problem statement ↔ FIGURES_GENERATION/figures_bg_bio/stg.py, lines 1141–1249 · score 0.63 · gating variables, membrane potential, reversal potential, capacitance, external, intracellular
  7. [7] § Results › General problem statement ↔ FIGURES_GENERATION/figures_methods/stg.py, lines 1141–1249 · score 0.63 · gating variables, membrane potential, reversal potential, capacitance, external, intracellular
  8. [8] § Materials and methods › Generating degenerate populations from DICs ↔ FIGURES_GENERATION/figures_bg_bio/stg.py, lines 1740–1848 · score 0.62 · calcium conductances, calcium dynamics, sensitivity matrix, linear, STG, iteratively
  9. [9] § Materials and methods › Generating degenerate populations from DICs ↔ FIGURES_GENERATION/figures_methods/stg.py, lines 1740–1850 · score 0.62 · calcium conductances, calcium dynamics, sensitivity matrix, linear, STG, iteratively
  10. [10] § Materials and methods › Transfer to the DA model ↔ FIGURES_GENERATION/figures_results/architecture_transfer_final_evaluate.py, lines 190–287 · score 0.60 · attention layers, LoRA, linear layers, adapters, network, Transfer
  11. [11] § Materials and methods › Transfer to the DA model ↔ FIGURES_GENERATION/figures_results/results_inference_da.ipynb, lines 203–300 · score 0.60 · attention layers, LoRA, linear layers, adapters, network
  12. [12] § Materials and methods › Architecture overview ↔ FIGURES_GENERATION/figures_results/poisson_process.ipynb, lines 122–210 · score 0.52 · positional encoding, ISIs, stack, embedding, layer, encoder
  13. [13] § Materials and methods › Architecture overview ↔ FIGURES_GENERATION/figures_results/architecture_final_evaluate.py, lines 18–106 · score 0.52 · positional encoding, ISIs, stack, embedding, layer, encoder
  14. [14] § Materials and methods › Evaluation metrics ↔ FIGURES_GENERATION/figures_results/results.ipynb, lines 326–413 · score 0.51 · primary metric, hyperparameter, validation, regression, head, activity

Paper

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

The paper is loaded when this pane is shown.

The authors' code

Python · 1,974 lines · 78 KB · no license · 4 matches

  1. import numpy as np
  2. from utils import gsigmoid
  3. from utils import get_w_factors, get_w_factors_constant_tau
  4. from utils import d_gsigmoid
  5. from utils import find_first_decreasing_zero_bisection
  6. from scipy.integrate import solve_ivp
  7. from utils import gamma_uniform_mean_std_matching
  9. # == simulation functions == #
  11. def simulate_individual(args):
  12. """
  13. Simulates the dynamics of a single STG neuron based on its conductance parameters.
  15. Parameters
  16. ----------
  17. args : tuple
  18. A tuple containing:
  19. - u0 : array-like
  20. Initial conditions for the state variables.
  21. - individual : array-like
  22. Conductance parameters for the neuron (e.g., maximal conductances for ion channels).
  23. - T_final : float
  24. Final time for the simulation.
  25. - dt : float
  26. Time step for evaluating the solution.
  27. - params : dict
  28. Dictionary of fixed neuron parameters (e.g., reversal potentials and calcium dynamics).
  30. Returns
  31. -------
  32. np.ndarray
  33. A 2D array where the first row contains the time points and the second row contains
  34. the voltage trace of the neuron during the simulation.
  35. """
  36. u0, individual, T_final, dt, params = args
  37. t_eval = np.arange(0, T_final, dt)
  38. return simulate_individual_t_eval((u0, individual, t_eval, params))
  41. def simulate_individual_t_eval(args):
  42. """
  43. Simulates the dynamics of a single STG neuron, evaluating the solution at specific time points.
  45. Parameters
  46. ----------
  47. args : tuple
  48. A tuple containing:
  49. - u0 : array-like
  50. Initial conditions for the state variables.
  51. - individual : array-like
  52. Conductance parameters for the neuron (e.g., maximal conductances for ion channels).
  53. - t_eval : array-like
  54. Time points at which to evaluate the solution.
  55. - params : dict
  56. Dictionary of fixed neuron parameters (e.g., reversal potentials and calcium dynamics).
  58. Returns
  59. -------
  60. np.ndarray
  61. A 2D array where the first row contains the specified time points and the second row contains
  62. the voltage trace of the neuron during the simulation.
  63. """
  64. u0, individual, t_eval, params = args
  65. sol = solve_ivp(
  66. ODEs,
  67. [0, t_eval[-1]],
  68. u0,
  69. t_eval=t_eval,
  70. args=(
  71. individual[0], individual[1], individual[2], individual[3],
  72. individual[4], individual[5], individual[6], individual[7],
  73. params['E_Na'], params['E_K'], params['E_H'], params['E_leak'],
  74. params['E_Ca'], params['alpha_Ca'], params['beta_Ca'], params['tau_Ca']
  75. ),
  76. method='BDF',
  77. dense_output=False,
  78. jac=jacobian
  79. )
  80. return np.array((sol.t, sol.y[0]))
  83. def get_u0(V0, Ca0):
  84. """
  85. Generates the initial conditions for the state variables of the STG neuron.
  87. Parameters
  88. ----------
  89. V0 : float
  90. Initial membrane voltage.
  91. Ca0 : float
  92. Initial intracellular calcium concentration.
  94. Returns
  95. -------
  96. np.ndarray
  97. An array of initial values for the state variables, including the gating variables and calcium concentration.
  98. """
  99. u0 = np.zeros(13)
  101. u0[0] = V0
  102. u0[1] = m_inf_Na(V0)
  103. u0[2] = h_inf_Na(V0)
  104. u0[3] = m_inf_Kd(V0)
  105. u0[4] = m_inf_CaT(V0)
  106. u0[5] = h_inf_CaT(V0)
  107. u0[6] = m_inf_CaS(V0)
  108. u0[7] = h_inf_CaS(V0)
  109. u0[8] = m_inf_KCa(V0, Ca0)
  110. u0[9] = m_inf_A(V0)
  111. u0[10] = h_inf_A(V0)
  112. u0[11] = m_inf_H(V0)
  113. u0[12] = Ca0
  115. return u0
  118. def get_default_parameters():
  119. """
  120. Provides the default neuron parameters for the STG model.
  122. Returns
  123. -------
  124. dict
  125. A dictionary containing the reversal potentials and calcium dynamics parameters.
  126. """
  127. params = {}
  129. params['E_leak'] = -50 # Leak reversal potential
  130. params['E_Na'] = 50 # Sodium reversal potential
  131. params['E_K'] = -80 # Potassium reversal potential
  132. params['E_H'] = -20 # H-current reversal potential
  133. params['E_Ca'] = 80 # Calcium reversal potential
  135. params['tau_Ca'] = 20 # Calcium decay time constant
  136. params['alpha_Ca'] = 0.94
  137. params['beta_Ca'] = 0.05
  139. return params
  142. def get_default_u0():
  143. """
  144. Provides default initial conditions for the STG neuron state variables.
  145. The default values are for a resting neuron with a membrane potential of -70 mV and a calcium concentration of 0.5 µM.
  146. Returns
  147. -------
  148. np.ndarray
  149. Default initial state variables, including resting membrane potential and calcium concentration.
  150. """
  151. V0 = -70 # Resting membrane potential
  152. Ca0 = 0.5 # Initial calcium concentration
  153. u0 = get_u0(V0, Ca0)
  154. return u0
  157. def get_best_set(g_s, g_u):
  158. """
  159. Determines the best set of conductances to neuromulate based on the slow and ultra-slow DIC conductance parameters.
  160. Derived from the reachability analysis of the STG model.
  162. Parameters
  163. ----------
  164. g_s : float
  165. Slow DIC conductance parameter.
  166. g_u : float
  167. Ultra-slow DIC conductance parameter.
  169. Returns
  170. -------
  171. list
  172. A list of strings identifying the best conductances ('A', 'H', or 'CaS').
  174. Notes
  175. -----
  176. - If `g_u` is negative, a warning is printed and the function returns ['CaS', 'A'].
  177. - If `g_s` is positive, the function returns ['A', 'H'].
  178. - If `g_s` is negative, the function returns ['CaS', 'H'].
  179. """
  180. if g_u < 0:
  181. print('Cautious, g_u is negative!')
  182. return ['CaS', 'A']
  184. if g_s >= 0:
  185. return ['A', 'H']
  186. if g_s < 0:
  187. return ['CaS', 'H']
  190. # == Gating variables functions == #
  192. def m_inf_Na(V):
  193. """
  194. Computes the steady-state activation variable (m) for sodium (Na) channels as a function of membrane potential.
  196. Parameters
  197. ----------
  198. V : float
  199. Membrane potential.
  201. Returns
  202. -------
  203. float
  204. The steady-state value of the m variable for Na channels.
  206. Reference
  207. ---------
  208. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  209. """
  210. return gsigmoid(V, A=0, B=1, C=-5.29, D=25.5)
  213. def h_inf_Na(V):
  214. """
  215. Computes the steady-state inactivation variable (h) for sodium (Na) channels as a function of membrane potential.
  217. Parameters
  218. ----------
  219. V : float
  220. Membrane potential.
  222. Returns
  223. -------
  224. float
  225. The steady-state value of the h variable for Na channels.
  227. Reference
  228. ---------
  229. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  230. """
  231. return gsigmoid(V, A=0, B=1, C=5.18, D=48.9)
  234. def tau_m_Na(V):
  235. """
  236. Computes the time constant for the activation variable (m) for sodium (Na) channels as a function of membrane potential.
  238. Parameters
  239. ----------
  240. V : float
  241. Membrane potential.
  243. Returns
  244. -------
  245. float
  246. The time constant of the m variable for Na channels.
  248. Reference
  249. ---------
  250. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  251. """
  252. return gsigmoid(V, A=1.32, B=-1.26, C=-25, D=120)
  255. def tau_h_Na(V):
  256. """
  257. Computes the time constant for the inactivation variable (h) for sodium (Na) channels as a function of membrane potential.
  259. Parameters
  260. ----------
  261. V : float
  262. Membrane potential.
  264. Returns
  265. -------
  266. float
  267. The time constant of the h variable for Na channels.
  269. Reference
  270. ---------
  271. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  272. """
  273. return gsigmoid(V, A=0, B=0.67, C=-10, D=62.9) * gsigmoid(V, A=1.5, B=1, C=3.6, D=34.9)
  276. def m_inf_Kd(V):
  277. """
  278. Computes the steady-state activation variable (m) for delayed rectifier potassium (Kd) channels as a function of membrane potential.
  280. Parameters
  281. ----------
  282. V : float
  283. Membrane potential.
  285. Returns
  286. -------
  287. float
  288. The steady-state value of the m variable for Kd channels.
  290. Reference
  291. ---------
  292. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  293. """
  294. return gsigmoid(V, A=0, B=1, C=-11.8, D=12.3)
  297. def tau_m_Kd(V):
  298. """
  299. Computes the time constant for the activation variable (m) for delayed rectifier potassium (Kd) channels as a function of membrane potential.
  301. Parameters
  302. ----------
  303. V : float
  304. Membrane potential.
  306. Returns
  307. -------
  308. float
  309. The time constant of the m variable for Kd channels.
  311. Reference
  312. ---------
  313. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  314. """
  315. return gsigmoid(V, A=7.2, B=-6.4, C=-19.2, D=28.3)
  318. def m_inf_CaT(V):
  319. """
  320. Computes the steady-state activation variable (m) for T-type calcium (CaT) channels as a function of membrane potential.
  322. Parameters
  323. ----------
  324. V : float
  325. Membrane potential.
  327. Returns
  328. -------
  329. float
  330. The steady-state value of the m variable for CaT channels.
  332. Reference
  333. ---------
  334. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  335. """
  336. return gsigmoid(V, A=0, B=1, C=-7.2, D=27.1)
  339. def h_inf_CaT(V):
  340. """
  341. Computes the steady-state inactivation variable (h) for T-type calcium (CaT) channels as a function of membrane potential.
  343. Parameters
  344. ----------
  345. V : float
  346. Membrane potential.
  348. Returns
  349. -------
  350. float
  351. The steady-state value of the h variable for CaT channels.
  353. Reference
  354. ---------
  355. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  356. """
  357. return gsigmoid(V, A=0, B=1, C=5.5, D=32.1)
  360. def tau_m_CaT(V):
  361. """
  362. Computes the time constant for the activation variable (m) for T-type calcium (CaT) channels as a function of membrane potential.
  364. Parameters
  365. ----------
  366. V : float
  367. Membrane potential.
  369. Returns
  370. -------
  371. float
  372. The time constant of the m variable for CaT channels.
  374. Reference
  375. ---------
  376. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  377. """
  378. return gsigmoid(V, A=21.7, B=-21.3, C=-20.5, D=68.1)
  381. def tau_h_CaT(V):
  382. """
  383. Computes the time constant for the inactivation variable (h) for T-type calcium (CaT) channels as a function of membrane potential.
  385. Parameters
  386. ----------
  387. V : float
  388. Membrane potential.
  390. Returns
  391. -------
  392. float
  393. The time constant of the h variable for CaT channels.
  395. Reference
  396. ---------
  397. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  398. """
  399. return gsigmoid(V, A=105, B=-89.8, C=-16.9, D=55)
  401. def m_inf_CaS(V):
  402. """
  403. Computes the steady-state activation variable (m) for S-type calcium (CaS) channels as a function of membrane potential.
  405. Parameters
  406. ----------
  407. V : float
  408. Membrane potential.
  410. Returns
  411. -------
  412. float
  413. The steady-state value of the m variable for CaS channels.
  415. Reference
  416. ---------
  417. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  418. """
  419. return gsigmoid(V, A=0, B=1, C=-8.1, D=33)
  422. def h_inf_CaS(V):
  423. """
  424. Computes the steady-state inactivation variable (h) for S-type calcium (CaS) channels as a function of membrane potential.
  426. Parameters
  427. ----------
  428. V : float
  429. Membrane potential.
  431. Returns
  432. -------
  433. float
  434. The steady-state value of the h variable for CaS channels.
  436. Reference
  437. ---------
  438. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  439. """
  440. return gsigmoid(V, A=0, B=1, C=6.2, D=60)
  443. def tau_m_CaS(V):
  444. """
  445. Computes the time constant for the activation variable (m) for S-type calcium (CaS) channels as a function of membrane potential.
  447. Parameters
  448. ----------
  449. V : float
  450. Membrane potential.
  452. Returns
  453. -------
  454. float
  455. The time constant of the m variable for CaS channels.
  457. Reference
  458. ---------
  459. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  460. """
  461. return 1.4 + 7/(np.exp((V+27)/10) + np.exp((V+70)/-13))
  464. def tau_h_CaS(V):
  465. """
  466. Computes the time constant for the inactivation variable (h) for S-type calcium (CaS) channels as a function of membrane potential.
  468. Parameters
  469. ----------
  470. V : float
  471. Membrane potential.
  473. Returns
  474. -------
  475. float
  476. The time constant of the h variable for CaS channels.
  478. Reference
  479. ---------
  480. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  481. """
  482. return 60 + 150/(np.exp((V+55)/9) + np.exp((V+65)/-16))
  485. def m_inf_KCa(V, Ca):
  486. """
  487. Computes the steady-state activation variable (m) for calcium-dependent potassium (KCa) channels as a function of membrane potential and calcium concentration.
  489. Parameters
  490. ----------
  491. V : float
  492. Membrane potential.
  493. Ca : float
  494. Calcium concentration.
  496. Returns
  497. -------
  498. float
  499. The steady-state value of the m variable for KCa channels.
  501. Reference
  502. ---------
  503. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  504. """
  505. return Ca/(Ca + 3) * gsigmoid(V, A=0, B=1, C=-12.6, D=28.3)
  508. def tau_m_KCa(V):
  509. """
  510. Computes the time constant for the activation variable (m) for calcium-dependent potassium (KCa) channels as a function of membrane potential.
  512. Parameters
  513. ----------
  514. V : float
  515. Membrane potential.
  517. Returns
  518. -------
  519. float
  520. The time constant of the m variable for KCa channels.
  522. Reference
  523. ---------
  524. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  525. """
  526. return gsigmoid(V, A=90.3, B=-75.1, C=-22.7, D=46)
  529. def m_inf_A(V):
  530. """
  531. Computes the steady-state activation variable (m) for A-type potassium (A) channels as a function of membrane potential.
  533. Parameters
  534. ----------
  535. V : float
  536. Membrane potential.
  538. Returns
  539. -------
  540. float
  541. The steady-state value of the m variable for A-type K channels.
  543. Reference
  544. ---------
  545. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  546. """
  547. return gsigmoid(V, A=0, B=1, C=-8.7, D=27.2)
  550. def h_inf_A(V):
  551. """
  552. Computes the steady-state inactivation variable (h) for A-type potassium (A) channels as a function of membrane potential.
  554. Parameters
  555. ----------
  556. V : float
  557. Membrane potential.
  559. Returns
  560. -------
  561. float
  562. The steady-state value of the h variable for A-type K channels.
  564. Reference
  565. ---------
  566. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  567. """
  568. return gsigmoid(V, A=0, B=1, C=4.9, D=56.9)
  571. def tau_m_A(V):
  572. """
  573. Computes the time constant for the activation variable (m) for A-type potassium (A) channels as a function of membrane potential.
  575. Parameters
  576. ----------
  577. V : float
  578. Membrane potential.
  580. Returns
  581. -------
  582. float
  583. The time constant of the m variable for A-type K channels.
  585. Reference
  586. ---------
  587. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  588. """
  589. return gsigmoid(V, A=11.6, B=-10.4, C=-15.2, D=32.9)
  592. def tau_h_A(V):
  593. """
  594. Computes the time constant for the inactivation variable (h) for A-type potassium (A) channels as a function of membrane potential.
  596. Parameters
  597. ----------
  598. V : float
  599. Membrane potential.
  601. Returns
  602. -------
  603. float
  604. The time constant of the h variable for A-type K channels.
  606. Reference
  607. ---------
  608. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  609. """
  610. return gsigmoid(V, A=38.6, B=-29.2, C=-26.5, D=38.9)
  613. def m_inf_H(V):
  614. """
  615. Computes the steady-state activation variable (m) for H-type potassium (H) channels as a function of membrane potential.
  617. Parameters
  618. ----------
  619. V : float
  620. Membrane potential.
  622. Returns
  623. -------
  624. float
  625. The steady-state value of the m variable for H-type K channels.
  627. Reference
  628. ---------
  629. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  630. """
  631. return gsigmoid(V, A=0, B=1, C=6, D=70)
  634. def tau_m_H(V):
  635. """
  636. Computes the time constant for the activation variable (m) for H-type potassium (H) channels as a function of membrane potential.
  638. Parameters
  639. ----------
  640. V : float
  641. Membrane potential.
  643. Returns
  644. -------
  645. float
  646. The time constant of the m variable for H-type K channels.
  648. Reference
  649. ---------
  650. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  651. """
  652. return gsigmoid(V, A=272, B=1499, C=-8.73, D=42.2)
  655. def tau_Ca_constant_function(V, tau=20):
  656. """
  657. Returns a constant time constant for calcium decay.
  659. Parameters
  660. ----------
  661. V : float
  662. Membrane potential (not used in this function, but kept for consistency with other functions).
  663. tau : float
  664. The constant time constant value.
  666. Returns
  667. -------
  668. float
  669. The constant time constant.
  671. Reference
  672. ---------
  673. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  674. """
  675. if np.isscalar(V):
  676. return tau
  677. return np.full_like(V, tau)
  679. # == Derivatives == #
  681. def d_m_inf_Na(V):
  682. """
  683. Compute the derivative of the Na activation gating variable, m_inf, with respect to membrane potential (V).
  685. Parameters
  686. ----------
  687. V : float or array
  688. Membrane potential in mV.
  690. Returns
  691. -------
  692. float or array
  693. Derivative of m_inf with respect to V.
  694. """
  695. return d_gsigmoid(V, A = 0, B = 1, C = -5.29, D = 25.5)
  697. def d_h_inf_Na(V):
  698. """
  699. Compute the derivative of the Na inactivation gating variable, h_inf, with respect to membrane potential (V).
  701. Parameters
  702. ----------
  703. V : float or array
  704. Membrane potential in mV.
  706. Returns
  707. -------
  708. float or array
  709. Derivative of h_inf with respect to V.
  710. """
  711. return d_gsigmoid(V, A = 0, B = 1, C = 5.18, D = 48.9)
  713. def d_m_inf_Kd(V):
  714. """
  715. Compute the derivative of the K delayed rectifier activation gating variable, m_inf, with respect to membrane potential (V).
  717. Parameters
  718. ----------
  719. V : float or array
  720. Membrane potential in mV.
  722. Returns
  723. -------
  724. float or array
  725. Derivative of m_inf with respect to V.
  726. """
  727. return d_gsigmoid(V, A = 0, B = 1, C = -11.8, D = 12.3)
  729. def d_m_inf_CaT(V):
  730. """
  731. Compute the derivative of the Ca T-type activation gating variable, m_inf, with respect to membrane potential (V).
  733. Parameters
  734. ----------
  735. V : float or array
  736. Membrane potential in mV.
  738. Returns
  739. -------
  740. float or array
  741. Derivative of m_inf with respect to V.
  742. """
  743. return d_gsigmoid(V, A = 0, B = 1, C = -7.2, D = 27.1)
  745. def d_h_inf_CaT(V):
  746. """
  747. Compute the derivative of the Ca T-type inactivation gating variable, h_inf, with respect to membrane potential (V).
  749. Parameters
  750. ----------
  751. V : float or array
  752. Membrane potential in mV.
  754. Returns
  755. -------
  756. float or array
  757. Derivative of h_inf with respect to V.
  758. """
  759. return d_gsigmoid(V, A = 0, B = 1, C = 5.5, D = 32.1)
  761. def d_m_inf_CaS(V):
  762. """
  763. Compute the derivative of the Ca S-type activation gating variable, m_inf, with respect to membrane potential (V).
  765. Parameters
  766. ----------
  767. V : float or array
  768. Membrane potential in mV.
  770. Returns
  771. -------
  772. float or array
  773. Derivative of m_inf with respect to V.
  774. """
  775. return d_gsigmoid(V, A = 0, B = 1, C = -8.1, D = 33)
  777. def d_h_inf_CaS(V):
  778. """
  779. Compute the derivative of the Ca S-type inactivation gating variable, h_inf, with respect to membrane potential (V).
  781. Parameters
  782. ----------
  783. V : float or array
  784. Membrane potential in mV.
  786. Returns
  787. -------
  788. float or array
  789. Derivative of h_inf with respect to V.
  790. """
  791. return d_gsigmoid(V, A = 0, B = 1, C = 6.2, D = 60)
  793. def d_m_inf_KCa_dV(V, Ca):
  794. """
  795. Compute the derivative of the K Ca-dependent activation gating variable, m_inf, with respect to membrane potential (V).
  797. Parameters
  798. ----------
  799. V : float or array
  800. Membrane potential in mV.
  801. Ca : float
  802. Intracellular calcium concentration in µM.
  804. Returns
  805. -------
  806. float or array
  807. Derivative of m_inf with respect to V.
  808. """
  809. return Ca/(Ca + 3) * d_gsigmoid(V, A = 0, B = 1, C = -12.6, D = 28.3)
  811. def d_m_inf_KCa_dCa(V, Ca):
  812. """
  813. Compute the derivative of the K Ca-dependent activation gating variable, m_inf, with respect to calcium concentration (Ca).
  815. Parameters
  816. ----------
  817. V : float or array
  818. Membrane potential in mV.
  819. Ca : float
  820. Intracellular calcium concentration in µM.
  822. Returns
  823. -------
  824. float or array
  825. Derivative of m_inf with respect to Ca.
  826. """
  827. return 3 * gsigmoid(V, A = 0, B = 1, C = -12.6, D = 28.3) / (Ca + 3)**2
  829. def d_Ca_inf_dV(V, alpha, E_Ca, g_CaT, g_CaS, m_inf_CaT_values, m_inf_CaS_values, h_inf_CaT_values, h_inf_CaS_values, d_m_inf_CaT_values, d_m_inf_CaS_values, d_h_inf_CaT_values, d_h_inf_CaS_values):
  830. """
  831. Compute the derivative of the intracellular calcium concentration with respect to membrane potential (V).
  833. Parameters
  834. ----------
  835. V : float or array
  836. Membrane potential in mV.
  837. alpha : float
  838. Constant for scaling the calcium concentration change.
  839. E_Ca : float
  840. Calcium reversal potential in mV.
  841. g_CaT : float
  842. Maximal conductance for Ca T-type channels.
  843. g_CaS : float
  844. Maximal conductance for Ca S-type channels.
  845. m_inf_CaT_values : array
  846. m_inf values for the Ca T-type channel.
  847. m_inf_CaS_values : array
  848. m_inf values for the Ca S-type channel.
  849. h_inf_CaT_values : array
  850. h_inf values for the Ca T-type channel.
  851. h_inf_CaS_values : array
  852. h_inf values for the Ca S-type channel.
  853. d_m_inf_CaT_values : array
  854. Derivative of m_inf values for Ca T-type.
  855. d_m_inf_CaS_values : array
  856. Derivative of m_inf values for Ca S-type.
  857. d_h_inf_CaT_values : array
  858. Derivative of h_inf values for Ca T-type.
  859. d_h_inf_CaS_values : array
  860. Derivative of h_inf values for Ca S-type.
  862. Returns
  863. -------
  864. float or array
  865. Derivative of calcium concentration with respect to V.
  867. Notes
  868. -----
  869. The equation is derived from the chain rule of differentiation. The original equation can be found in the STG model ODEs.
  870. """
  871. d = np.zeros_like(V)
  872. d = g_CaT * 3 * m_inf_CaT_values**2 * h_inf_CaT_values * d_m_inf_CaT_values * (V - E_Ca) +\
  873. g_CaT * m_inf_CaT_values**3 * d_h_inf_CaT_values * (V - E_Ca) +\
  874. g_CaT * m_inf_CaT_values**3 * h_inf_CaT_values +\
  875. g_CaS * 3 * m_inf_CaS_values**2 * h_inf_CaS_values * d_m_inf_CaS_values * (V - E_Ca) +\
  876. g_CaS * m_inf_CaS_values**3 * d_h_inf_CaS_values * (V - E_Ca) +\
  877. g_CaS * m_inf_CaS_values**3 * h_inf_CaS_values
  878. return - alpha * d
  880. def d_m_inf_A(V):
  881. """
  882. Compute the derivative of the A-type K activation gating variable, m_inf, with respect to membrane potential (V).
  884. Parameters
  885. ----------
  886. V : float or array
  887. Membrane potential in mV.
  889. Returns
  890. -------
  891. float or array
  892. Derivative of m_inf with respect to V.
  893. """
  894. return d_gsigmoid(V, A = 0, B = 1, C = -8.7, D = 27.2)
  896. def d_h_inf_A(V):
  897. """
  898. Compute the derivative of the A-type K inactivation gating variable, h_inf, with respect to membrane potential (V).
  900. Parameters
  901. ----------
  902. V : float or array
  903. Membrane potential in mV.
  905. Returns
  906. -------
  907. float or array
  908. Derivative of h_inf with respect to V.
  909. """
  910. return d_gsigmoid(V, A = 0, B = 1, C = 4.9, D = 56.9)
  912. def d_m_inf_H(V):
  913. """
  914. Compute the derivative of the H-type K activation gating variable, m_inf, with respect to membrane potential (V).
  916. Parameters
  917. ----------
  918. V : float or array
  919. Membrane potential in mV.
  921. Returns
  922. -------
  923. float or array
  924. Derivative of m_inf with respect to V.
  925. """
  926. return d_gsigmoid(V, A = 0, B = 1, C = 6, D = 70)
  928. # == UTILS == #
  930. def compute_equilibrium_Ca(alpha, I_Ca, beta):
  931. """
  932. Compute the equilibrium calcium concentration based on the calcium current (I_Ca),
  933. the scaling constant (alpha), and the calcium influx rate (beta).
  935. The equilibrium condition is given by:
  936. dCa/dt = 0 = -alpha * I_Ca - Ca + beta
  938. Parameters
  939. ----------
  940. alpha : float
  941. A constant scaling factor representing the effect of the calcium current (I_Ca) on calcium concentration.
  942. I_Ca : float
  943. The calcium current (typically in µA or similar units).
  944. beta : float
  945. The calcium influx rate (or leakage rate) in the system.
  947. Returns
  948. -------
  949. float
  950. The equilibrium calcium concentration at the steady-state (Ca).
  951. """
  952. return -alpha * I_Ca + beta
  955. def find_V_th_DICs(V, g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak,
  956. E_Na, E_K, E_H, E_leak, E_Ca, alpha_Ca, beta_Ca, tau_Ca,
  957. tau_f_stg = tau_m_Na, tau_s_stg = tau_m_Kd, tau_u_stg = tau_m_H, get_I_static = False, normalize = True, y_tol = 1e-6, x_tol=1e-6, max_iter = 1000, verbose=True):
  958. """
  959. Find the threshold voltage (V_th) for dynamic input conductances (DICs).
  961. This function uses a bisection method to find the first voltage where the total conductance
  962. (g_t) decreases to zero. It also returns the values of the DICs at this threshold voltage.
  964. Parameters
  965. ----------
  966. V : array-like
  967. Array of membrane potentials (in mV) to search for the threshold voltage.
  968. g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak : float
  969. Maximum conductances (in mS/cm²) for the respective ion channels.
  970. E_Na, E_K, E_H, E_leak, E_Ca : float
  971. Reversal potentials (in mV) for the respective ion channels.
  972. alpha_Ca : float
  973. Scaling factor for calcium influx.
  974. beta_Ca : float
  975. Rate constant for calcium extrusion.
  976. tau_Ca : float
  977. Time constant for calcium concentration dynamics.
  978. tau_f_stg, tau_s_stg, tau_u_stg : callable, optional
  979. Functions to compute the time constants for fast, slow, and ultra-slow dynamics.
  980. get_I_static : bool, optional
  981. If True, also compute the static current.
  982. normalize : bool, optional
  983. If True, normalize the sensitivity matrix by the leak conductance.
  984. y_tol : float, optional
  985. Tolerance for the y-axis (conductance) in the bisection method.
  986. x_tol : float, optional
  987. Tolerance for the x-axis (voltage) in the bisection method.
  988. max_iter : int, optional
  989. Maximum number of iterations for the bisection method.
  990. verbose : bool, optional
  991. If True, print additional information during the bisection process.
  993. Returns
  994. -------
  995. V_th : float
  996. The threshold voltage where the total conductance decreases to zero.
  997. values : tuple
  998. A tuple containing the values of the DICs (g_f, g_s, g_u, g_t) at the threshold voltage.
  1000. References
  1001. ----------
  1002. Fyon, A., Franci, A., Sacré, P., & Drion, G. (2024). Dimensionality reduction of neuronal degeneracy reveals two interfering physiological mechanisms. PNAS Nexus, 3(10), pgae415. https://doi.org/10.1093/pnasnexus/pgae415
  1003. """
  1004. g_t = lambda V_scalar : DICs(np.asarray([V_scalar,]), g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak, E_Na, E_K, E_H, E_leak, E_Ca, alpha_Ca, beta_Ca, tau_Ca, tau_f_stg, tau_s_stg, tau_u_stg, False, normalize)[3]
  1006. V_th = find_first_decreasing_zero_bisection(V, g_t, y_tol = y_tol, x_tol=x_tol, max_iter = max_iter, verbose=verbose)
  1007. V_th = np.asarray([V_th,], dtype=np.float64)
  1009. if V_th is None or np.isnan(V_th):
  1010. return V_th, (np.atleast_1d(np.nan), np.atleast_1d(np.nan), np.atleast_1d(np.nan), np.atleast_1d(np.nan))
  1012. values = DICs(V_th, g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak, E_Na, E_K, E_H, E_leak, E_Ca, alpha_Ca, beta_Ca, tau_Ca, tau_f_stg, tau_s_stg, tau_u_stg, get_I_static, normalize)
  1014. return V_th, values
  1017. # == ODEs == #
  1019. def jacobian(t, u, g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak, E_Na, E_K, E_H, E_leak, E_Ca, alpha_Ca, beta_Ca, tau_Ca):
  1020. """
  1021. Compute the Jacobian matrix for the STG model.
  1023. Parameters
  1024. ----------
  1025. t : float
  1026. Time variable (not used in this function, but kept for consistency with other functions).
  1027. u : array-like, shape (13,)
  1028. State vector containing the following variables:
  1029. - u[0]: V (membrane potential, mV)
  1030. - u[1]: m_Na (activation of sodium channel)
  1031. - u[2]: h_Na (inactivation of sodium channel)
  1032. - u[3]: m_Kd (activation of delayed rectifier potassium channel)
  1033. - u[4]: m_CaT (activation of T-type calcium channel)
  1034. - u[5]: h_CaT (inactivation of T-type calcium channel)
  1035. - u[6]: m_CaS (activation of S-type calcium channel)
  1036. - u[7]: h_CaS (inactivation of S-type calcium channel)
  1037. - u[8]: m_KCa (activation of calcium-activated potassium channel)
  1038. - u[9]: m_A (activation of A-type potassium channel)
  1039. - u[10]: h_A (inactivation of A-type potassium channel)
  1040. - u[11]: m_H (activation of H-current channel)
  1041. - u[12]: Ca (intracellular calcium concentration, μM)
  1042. g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak : float
  1043. Maximum conductances (mS/cm²) of the respective ion channels.
  1044. E_Na, E_K, E_H, E_leak, E_Ca : float
  1045. Reversal potentials (mV) of the respective ion channels.
  1046. alpha_Ca : float
  1047. Scaling factor for the calcium current.
  1048. beta_Ca : float
  1049. Calcium influx rate (or leakage rate).
  1050. tau_Ca : float
  1051. Time constant for calcium decay.
  1053. Returns
  1054. -------
  1055. J : ndarray, shape (13, 13)
  1056. The Jacobian matrix of the system, where J[i, j] is the partial derivative
  1057. of the i-th state variable's time derivative with respect to the j-th state variable.
  1058. The rows and columns of J correspond to the variables in `u`.
  1059. """
  1061. # Initialize Jacobian matrix (13x13 for your system)
  1062. J = np.zeros((13, 13))
  1064. # Unpack variables (for clarity)
  1065. V = u[0]
  1066. m_Na, h_Na = u[1], u[2]
  1067. m_Kd = u[3]
  1068. m_CaT, h_CaT = u[4], u[5]
  1069. m_CaS, h_CaS = u[6], u[7]
  1070. m_KCa = u[8]
  1071. m_A, h_A = u[9], u[10]
  1072. m_H = u[11]
  1073. Ca = u[12]
  1075. # 1. Compute partial derivatives for V equation (voltage dynamics)
  1076. J[0, 0] = -(
  1077. g_Na * m_Na**3 * h_Na + g_Kd * m_Kd**4 + g_CaT * m_CaT**3 * h_CaT +
  1078. g_CaS * m_CaS**3 * h_CaS + g_KCa * m_KCa**4 + g_A * m_A**3 * h_A +
  1079. g_H * m_H + g_leak
  1080. ) # ∂V_dot/∂V
  1082. J[0, 1] = -3 * g_Na * m_Na**2 * h_Na * (V - E_Na) # ∂V_dot/∂m_Na
  1083. J[0, 2] = -g_Na * m_Na**3 * (V - E_Na) # ∂V_dot/∂h_Na
  1084. J[0, 3] = -4 * g_Kd * m_Kd**3 * (V - E_K) # ∂V_dot/∂m_Kd
  1085. J[0, 4] = -3 * g_CaT * m_CaT**2 * h_CaT * (V - E_Ca) # ∂V_dot/∂m_CaT
  1086. J[0, 5] = -g_CaT * m_CaT**3 * (V - E_Ca) # ∂V_dot/∂h_CaT
  1087. J[0, 6] = -3 * g_CaS * m_CaS**2 * h_CaS * (V - E_Ca) # ∂V_dot/∂m_CaS
  1088. J[0, 7] = -g_CaS * m_CaS**3 * (V - E_Ca) # ∂V_dot/∂h_CaS
  1089. J[0, 8] = -4 * g_KCa * m_KCa**3 * (V - E_K) # ∂V_dot/∂m_KCa
  1090. J[0, 9] = -3 * g_A * m_A**2 * h_A * (V - E_K) # ∂V_dot/∂m_A
  1091. J[0, 10] = -g_A * m_A**3 * (V - E_K) # ∂V_dot/∂h_A
  1092. J[0, 11] = -g_H * (V - E_H) # ∂V_dot/∂m_H
  1094. # 2. Compute partial derivatives for calcium concentration dynamics
  1095. J[12, 0] = -(alpha_Ca * (g_CaT * m_CaT**3 * h_CaT + g_CaS * m_CaS**3 * h_CaS)) / tau_Ca # ∂Ca_dot/∂V
  1096. J[12, 4] = -(alpha_Ca * 3 * g_CaT * m_CaT**2 * h_CaT * (V - E_Ca)) / tau_Ca # ∂Ca_dot/∂m_CaT
  1097. J[12, 5] = -(alpha_Ca * g_CaT * m_CaT**3 * (V - E_Ca)) / tau_Ca # ∂Ca_dot/∂h_CaT
  1098. J[12, 6] = -(alpha_Ca * 3 * g_CaS * m_CaS**2 * h_CaS * (V - E_Ca)) / tau_Ca # ∂Ca_dot/∂m_CaS
  1099. J[12, 7] = -(alpha_Ca * g_CaS * m_CaS**3 * (V - E_Ca)) / tau_Ca # ∂Ca_dot/∂h_CaS
  1100. J[12, 12] = -1 / tau_Ca # ∂Ca_dot/∂Ca
  1102. # 3. Partial derivatives for gating variables
  1103. # m_Na and h_Na dynamics
  1104. J[1, 0] = d_m_inf_Na(V) / tau_m_Na(V) # ∂m_Na_dot/∂V
  1105. J[1, 1] = -1 / tau_m_Na(V) # ∂m_Na_dot/∂m_Na
  1106. J[2, 0] = d_h_inf_Na(V) / tau_h_Na(V) # ∂h_Na_dot/∂V
  1107. J[2, 2] = -1 / tau_h_Na(V) # ∂h_Na_dot/∂h_Na
  1109. # m_Kd dynamics
  1110. J[3, 0] = d_m_inf_Kd(V) / tau_m_Kd(V) # ∂m_Kd_dot/∂V
  1111. J[3, 3] = -1 / tau_m_Kd(V) # ∂m_Kd_dot/∂m_Kd
  1113. # m_CaT and h_CaT dynamics
  1114. J[4, 0] = d_m_inf_CaT(V) / tau_m_CaT(V) # ∂m_CaT_dot/∂V
  1115. J[4, 4] = -1 / tau_m_CaT(V) # ∂m_CaT_dot/∂m_CaT
  1116. J[5, 0] = d_h_inf_CaT(V) / tau_h_CaT(V) # ∂h_CaT_dot/∂V
  1117. J[5, 5] = -1 / tau_h_CaT(V) # ∂h_CaT_dot/∂h_CaT
  1119. # m_CaS and h_CaS dynamics
  1120. J[6, 0] = d_m_inf_CaS(V) / tau_m_CaS(V) # ∂m_CaS_dot/∂V
  1121. J[6, 6] = -1 / tau_m_CaS(V) # ∂m_CaS_dot/∂m_CaS
  1122. J[7, 0] = d_h_inf_CaS(V) / tau_h_CaS(V) # ∂h_CaS_dot/∂V
  1123. J[7, 7] = -1 / tau_h_CaS(V) # ∂h_CaS_dot/∂h_CaS
  1125. # m_KCa dynamics
  1126. J[8, 0] = d_m_inf_KCa_dV(V, Ca) / tau_m_KCa(V) # ∂m_KCa_dot/∂V
  1127. J[8, 8] = -1 / tau_m_KCa(V) # ∂m_KCa_dot/∂m_KCa
  1129. # m_A and h_A dynamics
  1130. J[9, 0] = d_m_inf_A(V) / tau_m_A(V) # ∂m_A_dot/∂V
  1131. J[9, 9] = -1 / tau_m_A(V) # ∂m_A_dot/∂m_A
  1132. J[10, 0] = d_h_inf_A(V) / tau_h_A(V) # ∂h_A_dot/∂V
  1133. J[10, 10] = -1 / tau_h_A(V) # ∂h_A_dot/∂h_A
  1135. # m_H dynamics
  1136. J[11, 0] = d_m_inf_H(V) / tau_m_H(V) # ∂m_H_dot/∂V
  1137. J[11, 11] = -1 / tau_m_H(V) # ∂m_H_dot/∂m_H
  1139. return J
  1141. def ODEs(t, u, g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak, E_Na, E_K, E_H, E_leak, E_Ca, alpha_Ca, beta_Ca, tau_Ca):
  1142. """
  1143. Compute the time derivatives of the state variables in the STG model.
  1145. This function implements the system of ordinary differential equations (ODEs) governing
  1146. the dynamics of the membrane potential, gating variables, and intracellular calcium concentration
  1147. in a neuron modeled after the stomatogastric ganglion (STG).
  1149. Parameters
  1150. ----------
  1151. t : float
  1152. Time variable (included for compatibility with ODE solvers, but not explicitly used in the equations).
  1153. u : array-like, shape (13,)
  1154. State vector containing the following variables in order:
  1155. - u[0] : V (membrane potential, mV)
  1156. - u[1] : m_Na (activation variable for sodium channel)
  1157. - u[2] : h_Na (inactivation variable for sodium channel)
  1158. - u[3] : m_Kd (activation variable for delayed rectifier potassium channel)
  1159. - u[4] : m_CaT (activation variable for T-type calcium channel)
  1160. - u[5] : h_CaT (inactivation variable for T-type calcium channel)
  1161. - u[6] : m_CaS (activation variable for S-type calcium channel)
  1162. - u[7] : h_CaS (inactivation variable for S-type calcium channel)
  1163. - u[8] : m_KCa (activation variable for calcium-activated potassium channel)
  1164. - u[9] : m_A (activation variable for A-type potassium channel)
  1165. - u[10]: h_A (inactivation variable for A-type potassium channel)
  1166. - u[11]: m_H (activation variable for H-current channel)
  1167. - u[12]: Ca (intracellular calcium concentration, μM)
  1168. g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak : float
  1169. Maximum conductances (μS/nF) of the respective ion channels.
  1170. E_Na, E_K, E_H, E_leak, E_Ca : float
  1171. Reversal potentials (mV) for sodium, potassium, leak, and calcium currents.
  1172. alpha_Ca : float
  1173. Proportionality constant for calcium influx (μM·cm²·ms⁻¹).
  1174. beta_Ca : float
  1175. Rate constant for calcium extrusion (ms⁻¹).
  1176. tau_Ca : float
  1177. Time constant for calcium concentration dynamics (ms).
  1179. Returns
  1180. -------
  1181. du : ndarray, shape (13,)
  1182. Time derivatives of the state variables. The output contains:
  1183. - du[0] : dV/dt (rate of change of membrane potential)
  1184. - du[1] : dm_Na/dt (rate of change of sodium channel activation)
  1185. - du[2] : dh_Na/dt (rate of change of sodium channel inactivation)
  1186. - du[3] : dm_Kd/dt (rate of change of delayed rectifier potassium activation)
  1187. - du[4] : dm_CaT/dt (rate of change of T-type calcium activation)
  1188. - du[5] : dh_CaT/dt (rate of change of T-type calcium inactivation)
  1189. - du[6] : dm_CaS/dt (rate of change of S-type calcium activation)
  1190. - du[7] : dh_CaS/dt (rate of change of S-type calcium inactivation)
  1191. - du[8] : dm_KCa/dt (rate of change of calcium-activated potassium activation)
  1192. - du[9] : dm_A/dt (rate of change of A-type potassium activation)
  1193. - du[10]: dh_A/dt (rate of change of A-type potassium inactivation)
  1194. - du[11]: dm_H/dt (rate of change of H-current activation)
  1195. - du[12]: dCa/dt (rate of change of intracellular calcium concentration)
  1197. References
  1198. ----------
  1199. The STG model is based on the following paper:
  1200. Liu, Z., Golowasch, J., Marder, E., & Abbott, L. F. (1998). A model neuron with activity-dependent conductances regulated by multiple calcium sensors. The Journal of neuroscience: the official journal of the Society for Neuroscience, 18(7), 2309–2320. https://doi.org/10.1523/JNEUROSCI.18-07-02309.1998
  1201. - But the equations are adapted and the g_bar dynamics is not considered in this implementation.
  1202. - The capacitance should be considered in the conductances directly.
  1203. - No external input current is considered in this implementation.
  1204. """
  1206. # Preallocate du (the vector of derivatives)
  1207. du = np.zeros_like(u)
  1208. # Extract variables (u is treated as a vector)
  1209. V = u[0]
  1210. m_Na, h_Na = u[1], u[2]
  1211. m_Kd = u[3]
  1212. m_CaT, h_CaT = u[4], u[5]
  1213. m_CaS, h_CaS = u[6], u[7]
  1214. m_KCa = u[8]
  1215. m_A, h_A = u[9], u[10]
  1216. m_H = u[11]
  1217. Ca = u[12]
  1219. # Compute currents using vectorized operations
  1220. I_Na = g_Na * m_Na**3 * h_Na * (V - E_Na)
  1221. I_Kd = g_Kd * m_Kd**4 * (V - E_K)
  1222. I_CaT = g_CaT * m_CaT**3 * h_CaT * (V - E_Ca)
  1223. I_CaS = g_CaS * m_CaS**3 * h_CaS * (V - E_Ca)
  1224. I_KCa = g_KCa * m_KCa**4 * (V - E_K)
  1225. I_A = g_A * m_A**3 * h_A * (V - E_K)
  1226. I_H = g_H * m_H * (V - E_H)
  1227. I_leak = g_leak * (V - E_leak)
  1229. # Compute the voltage derivative
  1230. du[0] = -(I_Na + I_Kd + I_CaT + I_CaS + I_KCa + I_A + I_H + I_leak)
  1232. # Calcium concentration derivative
  1233. I_Ca = I_CaT + I_CaS
  1234. du[12] = (-alpha_Ca * I_Ca - Ca + beta_Ca) / tau_Ca
  1236. # Gating variables (vectorized)
  1237. du[1] = (m_inf_Na(V) - m_Na) / tau_m_Na(V)
  1238. du[2] = (h_inf_Na(V) - h_Na) / tau_h_Na(V)
  1239. du[3] = (m_inf_Kd(V) - m_Kd) / tau_m_Kd(V)
  1240. du[4] = (m_inf_CaT(V) - m_CaT) / tau_m_CaT(V)
  1241. du[5] = (h_inf_CaT(V) - h_CaT) / tau_h_CaT(V)
  1242. du[6] = (m_inf_CaS(V) - m_CaS) / tau_m_CaS(V)
  1243. du[7] = (h_inf_CaS(V) - h_CaS) / tau_h_CaS(V)
  1244. du[8] = (m_inf_KCa(V, Ca) - m_KCa) / tau_m_KCa(V)
  1245. du[9] = (m_inf_A(V) - m_A) / tau_m_A(V)
  1246. du[10] = (h_inf_A(V) - h_A) / tau_h_A(V)
  1247. du[11] = (m_inf_H(V) - m_H) / tau_m_H(V)
  1249. return du
  1251. # == DICs related functions == #
  1253. def DICs(V, g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak,
  1254. E_Na, E_K, E_H, E_leak, E_Ca, alpha_Ca, beta_Ca, tau_Ca,
  1255. tau_f_stg=tau_m_Na, tau_s_stg=tau_m_Kd, tau_u_stg=tau_m_H, get_I_static=False, normalize=True):
  1256. """
  1257. Computes the dynamic input conductances (DICs) for a given set of membrane potentials and conductances.
  1259. Parameters
  1260. ----------
  1261. V : array-like
  1262. Membrane potentials at which to compute the DICs. Can be a scalar or a 1D array.
  1263. g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak : array-like
  1264. Maximum conductances for the respective ion channels. Each can be a scalar or a 1D array.
  1265. E_Na, E_K, E_H, E_leak, E_Ca : float
  1266. Reversal potentials for the respective ion channels.
  1267. alpha_Ca : float
  1268. Proportionality constant for calcium influx.
  1269. beta_Ca : float
  1270. Rate constant for calcium extrusion.
  1271. tau_Ca : float
  1272. Time constant for calcium concentration dynamics.
  1273. tau_f_stg, tau_s_stg, tau_u_stg : callable, optional
  1274. Functions to compute the time constants for fast, slow, and ultra-slow dynamics.
  1275. get_I_static : bool, optional
  1276. If True, also compute the static current.
  1277. normalize : bool, optional
  1278. If True, normalize the sensitivity matrix by the leak conductance.
  1280. Returns
  1281. -------
  1282. g_f, g_s, g_u, g_t : array-like
  1283. The fast, slow, ultra-slow, and total conductances. Each can be a scalar, a 1D array, or a 2D array.
  1284. I_static : array-like, optional
  1285. The static current, returned only if `get_I_static` is True.
  1287. Notes
  1288. -----
  1289. - The dimensions of the inputs can vary:
  1290. - If all inputs are scalars, the outputs will be scalars.
  1291. - If `V` is a 1D array and the conductances are scalars, the outputs will be 1D arrays of length N (N,).
  1292. - If `V` is a scalar and the conductances are 1D arrays of length M, the outputs will be 1D arrays of length M (M,).
  1293. - If both `V` and the conductances are 1D arrays of length N and M respectively, the outputs will be 2D arrays of shape (M, N).
  1295. References
  1296. ----------
  1297. Drion, G., Franci, A., Dethier, J., & Sepulchre, R. (2015). Dynamic Input Conductances Shape Neuronal Spiking. eNeuro, 2(1), ENEURO.0031-14.2015. https://doi.org/10.1523/ENEURO.0031-14.2015
  1298. """
  1300. # get the S matrix
  1301. if not get_I_static:
  1302. S = sensitivity_matrix(V, g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak,
  1303. E_Na, E_K, E_H, E_leak, E_Ca, alpha_Ca, beta_Ca, tau_Ca,
  1304. tau_f_stg, tau_s_stg, tau_u_stg, normalize, get_I_static)
  1305. else:
  1306. S, S_static = sensitivity_matrix(V, g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak,
  1307. E_Na, E_K, E_H, E_leak, E_Ca, alpha_Ca, beta_Ca, tau_Ca,
  1308. tau_f_stg, tau_s_stg, tau_u_stg, normalize, get_I_static)
  1310. if S.ndim == 3:
  1311. S = S[np.newaxis, :, :, :]
  1313. m = S.shape[0]
  1315. g_vec = np.array([g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak]).T
  1316. g_vec = np.atleast_2d(g_vec)
  1317. g_vec = g_vec[:, :, np.newaxis, np.newaxis].transpose(0, 2, 1, 3)
  1319. S_mul = np.sum(S * g_vec, axis=2)
  1321. g_f = S_mul[:, 0, :]
  1322. g_s = S_mul[:, 1, :]
  1323. g_u = S_mul[:, 2, :]
  1324. g_t = g_f + g_s + g_u
  1326. if get_I_static:
  1327. V_Na = V - E_Na
  1328. V_K = V - E_K
  1329. V_Ca = V - E_Ca
  1330. V_H = V - E_H
  1331. V_leak = V - E_leak
  1333. V_vec = np.array([V_Na, V_K, V_Ca, V_Ca, V_K, V_K, V_H, V_leak])
  1334. g_vec = g_vec[:, 0, :, 0][:, :, np.newaxis]
  1336. I_static = np.sum(S_static * g_vec * V_vec, axis=1)
  1338. if m == 1:
  1339. I_static = I_static[0]
  1340. g_f = g_f[0]
  1341. g_s = g_s[0]
  1342. g_u = g_u[0]
  1343. g_t = g_t[0]
  1345. return g_f, g_s, g_u, g_t, I_static
  1347. if m == 1:
  1348. g_f = g_f[0]
  1349. g_s = g_s[0]
  1350. g_u = g_u[0]
  1351. g_t = g_t[0]
  1353. return g_f, g_s, g_u, g_t
  1355. def sensitivity_matrix(V, g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak,
  1356. E_Na, E_K, E_H, E_leak, E_Ca, alpha_Ca, beta_Ca, tau_Ca,
  1357. tau_f_stg = tau_m_Na, tau_s_stg = tau_m_Kd, tau_u_stg = tau_m_H, normalize = True, get_I_static = False):
  1358. """
  1359. Computes the sensitivity matrix for the dynamic input conductances (DICs) of a neuron model.
  1361. Parameters
  1362. ----------
  1363. V : array-like
  1364. Membrane potentials at which to compute the sensitivity matrix. Can be a scalar or a 1D array.
  1365. g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, g_leak : array-like
  1366. Maximum conductances for the respective ion channels. Each can be a scalar or a 1D array.
  1367. E_Na, E_K, E_H, E_leak, E_Ca : float
  1368. Reversal potentials for the respective ion channels.
  1369. alpha_Ca : float
  1370. Proportionality constant for calcium influx.
  1371. beta_Ca : float
  1372. Rate constant for calcium extrusion.
  1373. tau_Ca : float
  1374. Time constant for calcium concentration dynamics.
  1375. tau_f_stg, tau_s_stg, tau_u_stg : callable, optional
  1376. Functions to compute the time constants for fast, slow, and ultra-slow dynamics.
  1377. normalize : bool, optional
  1378. If True, normalize the sensitivity matrix by the leak conductance.
  1379. get_I_static : bool, optional
  1380. If True, also compute the static current.
  1382. Returns
  1383. -------
  1384. S : ndarray
  1385. The sensitivity matrix of shape (m, 3, 8, n) where n is the length of V and m is the number of conductances.
  1386. S_static : ndarray, optional
  1387. The static current sensitivity matrix, returned only if `get_I_static` is True.
  1389. Notes
  1390. -----
  1391. - The dimensions of the inputs can vary:
  1392. - If all inputs are scalars, the outputs will be 3D arrays of shape (3, 8, 1).
  1393. - If `V` is a 1D array and the conductances are scalars, the outputs will be 3D arrays of shape (3, 8, N).
  1394. - If `V` is a scalar and the conductances are 1D arrays of length M, the outputs will be 3D arrays of shape (M, 3, 8, 1).
  1395. - If both `V` and the conductances are 1D arrays of length N and M respectively, the outputs will be 4D arrays of shape (M, 3, 8, N).
  1397. References
  1398. ----------
  1399. Drion, G., Franci, A., Dethier, J., & Sepulchre, R. (2015). Dynamic Input Conductances Shape Neuronal Spiking. eNeuro, 2(1), ENEURO.0031-14.2015. https://doi.org/10.1523/ENEURO.0031-14.2015
  1400. """
  1402. V = np.atleast_1d(V)
  1403. g_Na = np.atleast_1d(g_Na)
  1404. g_Kd = np.atleast_1d(g_Kd)
  1405. g_CaT = np.atleast_1d(g_CaT)
  1406. g_CaS = np.atleast_1d(g_CaS)
  1407. g_KCa = np.atleast_1d(g_KCa)
  1408. g_A = np.atleast_1d(g_A)
  1409. g_H = np.atleast_1d(g_H)
  1410. g_leak = np.atleast_1d(g_leak)
  1412. # m represents the number of conductances, it will be 1 for scalars, or match the size of arrays
  1413. m = g_Na.size # we assume all conductances have the same size ... should be the case with a proper call to the function !
  1414. n = V.size # n is the size of V
  1416. # S will be a 3xN matrix with N = 7+1. Each row will correspond to a different variable (f, s, u) and each column to a different channel.
  1417. S_Na = np.zeros((m, 3, n))
  1418. S_Kd = np.zeros((m, 3, n))
  1419. S_CaT = np.zeros((m, 3, n))
  1420. S_CaS = np.zeros((m, 3, n))
  1421. S_KCa = np.zeros((m, 3, n))
  1422. S_A = np.zeros((m, 3, n))
  1423. S_H = np.zeros((m, 3, n))
  1424. S_leak = np.zeros((m, 3, n))
  1426. V = np.atleast_2d(V).T
  1428. m_inf_Na_values = m_inf_Na(V)
  1429. h_inf_Na_values = h_inf_Na(V)
  1430. m_inf_Kd_values = m_inf_Kd(V)
  1431. m_inf_CaT_values = m_inf_CaT(V)
  1432. h_inf_CaT_values = h_inf_CaT(V)
  1433. m_inf_CaS_values = m_inf_CaS(V)
  1434. h_inf_CaS_values = h_inf_CaS(V)
  1436. I_CaT = g_CaT * m_inf_CaT_values**3 * h_inf_CaT_values * (V - E_Ca)
  1437. I_CaS = g_CaS * m_inf_CaS_values**3 * h_inf_CaS_values * (V - E_Ca)
  1439. I_Ca = I_CaT + I_CaS
  1440. Ca = compute_equilibrium_Ca(alpha_Ca, I_Ca, beta_Ca)
  1442. m_inf_KCa_values = m_inf_KCa(V, Ca)
  1443. m_inf_A_values = m_inf_A(V)
  1444. h_inf_A_values = h_inf_A(V)
  1445. m_inf_H_values = m_inf_H(V)
  1447. d_m_inf_Na_values = d_m_inf_Na(V)
  1448. d_h_inf_Na_values = d_h_inf_Na(V)
  1449. d_m_inf_Kd_values = d_m_inf_Kd(V)
  1450. d_m_inf_CaT_values = d_m_inf_CaT(V)
  1451. d_h_inf_CaT_values = d_h_inf_CaT(V)
  1452. d_m_inf_CaS_values = d_m_inf_CaS(V)
  1453. d_h_inf_CaS_values = d_h_inf_CaS(V)
  1454. d_m_inf_A_values = d_m_inf_A(V)
  1455. d_h_inf_A_values = d_h_inf_A(V)
  1456. d_m_inf_H_values = d_m_inf_H(V)
  1457. d_m_inf_KCa_dCa_values = d_m_inf_KCa_dCa(V, Ca)
  1458. d_m_inf_KCa_dV_values = d_m_inf_KCa_dV(V, Ca)
  1459. d_Ca_inf_dV_values = d_Ca_inf_dV(V, alpha_Ca, E_Ca, g_CaT, g_CaS, m_inf_CaT_values, m_inf_CaS_values, h_inf_CaT_values, h_inf_CaS_values, d_m_inf_CaT_values, d_m_inf_CaS_values, d_h_inf_CaT_values, d_h_inf_CaS_values)
  1461. S_Na[:, 0, :] += (m_inf_Na_values**3 * h_inf_Na_values).T
  1462. S_Kd[:,0,:] += (m_inf_Kd_values**4).T
  1463. S_CaT[:,0,:] += (m_inf_CaT_values**3 * h_inf_CaT_values).T
  1464. S_CaS[:,0,:] += (m_inf_CaS_values**3 * h_inf_CaS_values).T
  1465. S_KCa[:,0,:] += (m_inf_KCa_values**4).T
  1466. S_A[:,0,:] += (m_inf_A_values**3 * h_inf_A_values).T
  1467. S_H[:,0,:] += (m_inf_H_values).T
  1468. S_leak[:,0,:] += 1.0
  1470. if get_I_static:
  1471. S_Na_static = S_Na[:, 0, :].copy()
  1472. S_Kd_static = S_Kd[:, 0, :].copy()
  1473. S_CaT_static = S_CaT[:, 0, :].copy()
  1474. S_CaS_static = S_CaS[:, 0, :].copy()
  1475. S_KCa_static = S_KCa[:, 0, :].copy()
  1476. S_A_static = S_A[:, 0, :].copy()
  1477. S_H_static = S_H[:, 0, :].copy()
  1478. S_leak_static = S_leak[:, 0, :].copy()
  1480. dV_dot_dm_Na = - 3 * m_inf_Na_values**2 * h_inf_Na_values * d_m_inf_Na_values * (V - E_Na)
  1481. dV_dot_dh_Na = - m_inf_Na_values**3 * d_h_inf_Na_values * (V - E_Na)
  1482. dV_dot_dm_Kd = - 4 * m_inf_Kd_values**3 * d_m_inf_Kd_values * (V - E_K)
  1483. dV_dot_dm_CaT = - 3 * m_inf_CaT_values**2 * h_inf_CaT_values * d_m_inf_CaT_values * (V - E_Ca)
  1484. dV_dot_dh_CaT = - m_inf_CaT_values**3 * d_h_inf_CaT_values * (V - E_Ca)
  1485. dV_dot_dm_CaS = - 3 * m_inf_CaS_values**2 * h_inf_CaS_values * d_m_inf_CaS_values * (V - E_Ca)
  1486. dV_dot_dh_CaS = - m_inf_CaS_values**3 * d_h_inf_CaS_values * (V - E_Ca)
  1487. dV_dot_dm_KCa = - 4 * m_inf_KCa_values**3 * d_m_inf_KCa_dV_values * (V - E_K)
  1488. dV_dot_dCa_KCa = - 4 * m_inf_KCa_values**3 * d_m_inf_KCa_dCa_values * (V - E_K) * d_Ca_inf_dV_values
  1489. dV_dot_dm_A = - 3 * m_inf_A_values**2 * h_inf_A_values * d_m_inf_A_values * (V - E_K)
  1490. dV_dot_dh_A = - m_inf_A_values**3 * d_h_inf_A_values * (V - E_K)
  1491. dV_dot_dm_H = - d_m_inf_H_values * (V - E_H)
  1493. w_fs_m_Na, w_su_m_Na = get_w_factors(V, tau_m_Na, tau_f_stg, tau_s_stg, tau_u_stg)
  1494. w_fs_h_Na, w_su_h_Na = get_w_factors(V, tau_h_Na, tau_f_stg, tau_s_stg, tau_u_stg)
  1495. w_fs_m_Kd, w_su_m_Kd = get_w_factors(V, tau_m_Kd, tau_f_stg, tau_s_stg, tau_u_stg)
  1496. w_fs_m_CaT, w_su_m_CaT = get_w_factors(V, tau_m_CaT, tau_f_stg, tau_s_stg, tau_u_stg)
  1497. w_fs_h_CaT, w_su_h_CaT = get_w_factors(V, tau_h_CaT, tau_f_stg, tau_s_stg, tau_u_stg)
  1498. w_fs_m_CaS, w_su_m_CaS = get_w_factors(V, tau_m_CaS, tau_f_stg, tau_s_stg, tau_u_stg)
  1499. w_fs_h_CaS, w_su_h_CaS = get_w_factors(V, tau_h_CaS, tau_f_stg, tau_s_stg, tau_u_stg)
  1500. w_fs_m_KCa, w_su_m_KCa = get_w_factors(V, tau_m_KCa, tau_f_stg, tau_s_stg, tau_u_stg)
  1501. w_fs_m_KCa2, w_su_m_KCa2 = get_w_factors_constant_tau(V, tau_Ca, tau_f_stg, tau_s_stg, tau_u_stg) # TO BE CHANGED
  1502. w_fs_m_A, w_su_m_A = get_w_factors(V, tau_m_A, tau_f_stg, tau_s_stg, tau_u_stg)
  1503. w_fs_h_A, w_su_h_A = get_w_factors(V, tau_h_A, tau_f_stg, tau_s_stg, tau_u_stg)
  1504. w_fs_m_H, w_su_m_H = get_w_factors(V, tau_m_H, tau_f_stg, tau_s_stg, tau_u_stg)
  1506. S_Na[:,0,:] += (- w_fs_m_Na * dV_dot_dm_Na - w_fs_h_Na * dV_dot_dh_Na).T
  1507. S_Kd[:,0,:] += (- w_fs_m_Kd * dV_dot_dm_Kd).T
  1508. S_CaT[:,0,:] += (- w_fs_m_CaT * dV_dot_dm_CaT - w_fs_h_CaT * dV_dot_dh_CaT).T
  1509. S_CaS[:,0,:] += (- w_fs_m_CaS * dV_dot_dm_CaS - w_fs_h_CaS * dV_dot_dh_CaS).T
  1510. S_KCa[:,0,:] += (- w_fs_m_KCa * dV_dot_dm_KCa - w_fs_m_KCa2 * dV_dot_dCa_KCa).T
  1511. S_A[:,0,:] += (- w_fs_m_A * dV_dot_dm_A - w_fs_h_A * dV_dot_dh_A).T
  1512. S_H[:,0,:] += (- w_fs_m_H * dV_dot_dm_H).T
  1514. S_Na[:,1,:] += (- (w_su_m_Na - w_fs_m_Na) * dV_dot_dm_Na - (w_su_h_Na - w_fs_h_Na) * dV_dot_dh_Na).T
  1515. S_Kd[:,1,:] += (- (w_su_m_Kd - w_fs_m_Kd) * dV_dot_dm_Kd).T
  1516. S_CaT[:,1,:] += (- (w_su_m_CaT - w_fs_m_CaT) * dV_dot_dm_CaT - (w_su_h_CaT - w_fs_h_CaT) * dV_dot_dh_CaT).T
  1517. S_CaS[:,1,:] += (- (w_su_m_CaS - w_fs_m_CaS) * dV_dot_dm_CaS - (w_su_h_CaS - w_fs_h_CaS) * dV_dot_dh_CaS).T
  1518. S_KCa[:,1,:] += (- (w_su_m_KCa - w_fs_m_KCa) * dV_dot_dm_KCa - (w_su_m_KCa2 - w_fs_m_KCa2) * dV_dot_dCa_KCa).T
  1519. S_A[:,1,:] += (- (w_su_m_A - w_fs_m_A) * dV_dot_dm_A - (w_su_h_A - w_fs_h_A) * dV_dot_dh_A).T
  1520. S_H[:,1,:] += (- (w_su_m_H - w_fs_m_H) * dV_dot_dm_H).T
  1522. S_Na[:,2,:] += (- (1 - w_su_m_Na) * dV_dot_dm_Na - (1 - w_su_h_Na) * dV_dot_dh_Na).T
  1523. S_Kd[:,2,:] += (- (1 - w_su_m_Kd) * dV_dot_dm_Kd).T
  1524. S_CaT[:,2,:] += (- (1 - w_su_m_CaT) * dV_dot_dm_CaT - (1 - w_su_h_CaT) * dV_dot_dh_CaT).T
  1525. S_CaS[:,2,:] += (- (1 - w_su_m_CaS) * dV_dot_dm_CaS - (1 - w_su_h_CaS) * dV_dot_dh_CaS).T
  1526. S_KCa[:,2,:] += (- (1 - w_su_m_KCa) * dV_dot_dm_KCa - (1 - w_su_m_KCa2) * dV_dot_dCa_KCa).T
  1527. S_A[:,2,:] += (- (1 - w_su_m_A) * dV_dot_dm_A - (1 - w_su_h_A) * dV_dot_dh_A).T
  1528. S_H[:,2,:] += (- (1 - w_su_m_H) * dV_dot_dm_H).T
  1530. if normalize:
  1531. S_Na /= g_leak[:, np.newaxis, np.newaxis]
  1532. S_Kd /= g_leak[:, np.newaxis, np.newaxis]
  1533. S_CaT /= g_leak[:, np.newaxis, np.newaxis]
  1534. S_CaS /= g_leak[:, np.newaxis, np.newaxis]
  1535. S_KCa /= g_leak[:, np.newaxis, np.newaxis]
  1536. S_A /= g_leak[:, np.newaxis, np.newaxis]
  1537. S_H /= g_leak[:, np.newaxis, np.newaxis]
  1538. S_leak /= g_leak[:, np.newaxis, np.newaxis]
  1540. S_Na = S_Na[:, :, np.newaxis, :]
  1541. S_Kd = S_Kd[:, :, np.newaxis, :]
  1542. S_CaT = S_CaT[:, :, np.newaxis, :]
  1543. S_CaS = S_CaS[:, :, np.newaxis, :]
  1544. S_KCa = S_KCa[:, :, np.newaxis, :]
  1545. S_A = S_A[:, :, np.newaxis, :]
  1546. S_H = S_H[:, :, np.newaxis, :]
  1547. S_leak = S_leak[:, :, np.newaxis, :]
  1549. if not get_I_static:
  1550. if m == 1:
  1551. return np.concatenate((S_Na[0], S_Kd[0], S_CaT[0], S_CaS[0], S_KCa[0], S_A[0], S_H[0], S_leak[0]), axis=1)
  1552. else:
  1553. return np.concatenate((S_Na, S_Kd, S_CaT, S_CaS, S_KCa, S_A, S_H, S_leak), axis=2)
  1554. else:
  1555. if m == 1:
  1556. return np.concatenate((S_Na[0], S_Kd[0], S_CaT[0], S_CaS[0], S_KCa[0], S_A[0], S_H[0], S_leak[0]), axis=1), np.stack((S_Na_static, S_Kd_static, S_CaT_static, S_CaS_static, S_KCa_static, S_A_static, S_H_static, S_leak_static), axis=1)
  1557. else:
  1558. return np.concatenate((S_Na, S_Kd, S_CaT, S_CaS, S_KCa, S_A, S_H, S_leak), axis=2), np.stack((S_Na_static, S_Kd_static, S_CaT_static, S_CaS_static, S_KCa_static, S_A_static, S_H_static, S_leak_static), axis=1)
  1560. # == Compensation algorithms and generation functions == #
  1562. def generate_population(n_cells, V_th, g_f_target, g_s_target, g_u_target, g_bar_range_leak, g_bar_range_Na, g_bar_range_Kd, g_bar_range_CaT, g_bar_range_CaS, g_bar_range_KCa, g_bar_range_A, g_bar_range_H, params, default_g_CaS_for_Ca = 10., default_g_CaT_for_Ca = 6.0, distribution='uniform', normalize_by_leak=True):
  1563. """
  1564. Generates a population of neurons with specified dynamic input conductances (DICs).
  1566. Parameters
  1567. ----------
  1568. n_cells : int
  1569. Number of neurons to generate.
  1570. V_th : float
  1571. Threshold voltage for dynamic input conductances (DICs).
  1572. g_f_target : float
  1573. Target fast DIC.
  1574. g_s_target : float
  1575. Target slow DIC.
  1576. g_u_target : float
  1577. Target ultra-slow DIC.
  1578. g_bar_range_leak : list
  1579. Range of leak conductances.
  1580. g_bar_range_Na : list
  1581. Range of sodium conductances.
  1582. g_bar_range_Kd : list
  1583. Range of delayed rectifier potassium conductances.
  1584. g_bar_range_CaT : list
  1585. Range of T-type calcium conductances.
  1586. g_bar_range_CaS : list
  1587. Range of S-type calcium conductances.
  1588. g_bar_range_KCa : list
  1589. Range of calcium-activated potassium conductances.
  1590. g_bar_range_A : list
  1591. Range of A-type potassium conductances.
  1592. g_bar_range_H : list
  1593. Range of H-current conductances.
  1594. params : dict
  1595. Dictionary of fixed neuron parameters (e.g., reversal potentials and calcium dynamics).
  1596. default_g_CaS_for_Ca : float, optional
  1597. Default S-type calcium conductance for calcium dynamics.
  1598. default_g_CaT_for_Ca : float, optional
  1599. Default T-type calcium conductance for calcium dynamics.
  1600. distribution : str, optional
  1601. Distribution type for generating conductances ('uniform' or 'gamma').
  1602. normalize_by_leak : bool, optional
  1603. If True, normalize the conductances by the leak conductance.
  1605. Returns
  1606. -------
  1607. np.ndarray
  1608. Array of generated neuron conductances.
  1610. Notes
  1611. -----
  1612. - The population generation from this method can fail. In practice, the method using neuromodulation introduced by A. Fyon et al. (2015) and improved by this work can be used to generate a population of neurons with specified DICs.
  1613. - The method using neuromodulation is implemented in the `generate_neuromodulated_population` function and using the `get_best_set` function ensure reachability from a spiking population.
  1614. - The method is based on the work of Drion et al. (2015) and the compensation algorithm for DICs.
  1616. - If g_CaS and g_CaT are among the conductances to be compensated, the system is non-linear and providing a default value for the conductances is necessary. The compensation is, in this case, approximate and the reached DICs may not be exactly the target DICs.
  1617. References
  1618. ----------
  1619. For the compensation algorithm:
  1620. - Drion, G., Franci, A., Dethier, J., & Sepulchre, R. (2015). Dynamic Input Conductances Shape Neuronal Spiking. eNeuro, 2(1), ENEURO.0031-14.2015. https://doi.org/10.1523/ENEURO.0031-14.2015
  1621. For the population generation:
  1622. - Fyon, A., Franci, A., Sacré, P., & Drion, G. (2024). Dimensionality reduction of neuronal degeneracy reveals two interfering physiological mechanisms. PNAS Nexus, 3(10), pgae415. https://doi.org/10.1093/pnasnexus/pgae415
  1623. """
  1624. g_Na = np.random.uniform(g_bar_range_Na[0], g_bar_range_Na[1], n_cells) if g_bar_range_Na is not None else np.full(n_cells, np.nan)
  1625. g_Kd = np.random.uniform(g_bar_range_Kd[0], g_bar_range_Kd[1], n_cells) if g_bar_range_Kd is not None else np.full(n_cells, np.nan)
  1626. g_CaT = np.random.uniform(g_bar_range_CaT[0], g_bar_range_CaT[1], n_cells) if g_bar_range_CaT is not None else np.full(n_cells, np.nan)
  1627. g_CaS = np.random.uniform(g_bar_range_CaS[0], g_bar_range_CaS[1], n_cells) if g_bar_range_CaS is not None else np.full(n_cells, np.nan)
  1628. g_KCa = np.random.uniform(g_bar_range_KCa[0], g_bar_range_KCa[1], n_cells) if g_bar_range_KCa is not None else np.full(n_cells, np.nan)
  1629. g_A = np.random.uniform(g_bar_range_A[0], g_bar_range_A[1], n_cells) if g_bar_range_A is not None else np.full(n_cells, np.nan)
  1630. g_H = np.random.uniform(g_bar_range_H[0], g_bar_range_H[1], n_cells) if g_bar_range_H is not None else np.full(n_cells, np.nan)
  1632. if distribution == 'uniform':
  1633. # here we assume that the ranges are the min and max values for the uniform distribution
  1634. mean_leak = (g_bar_range_leak[0] + g_bar_range_leak[1]) / 2
  1635. g_leak = np.random.uniform(g_bar_range_leak[0], g_bar_range_leak[1], n_cells)
  1637. elif distribution == 'gamma':
  1638. g_bar_range_leak = gamma_uniform_mean_std_matching(*g_bar_range_leak)
  1639. mean_leak = g_bar_range_leak[0] * g_bar_range_leak[1]
  1640. g_leak = np.random.gamma(g_bar_range_leak[0], g_bar_range_leak[1], n_cells)
  1641. else:
  1642. raise ValueError('Invalid distribution type ! Please use either "uniform" or "gamma".')
  1644. if normalize_by_leak:
  1645. f = g_leak/mean_leak
  1646. g_Na *= f
  1647. g_Kd *= f
  1648. g_CaT *= f
  1649. g_CaS *= f
  1650. g_KCa *= f
  1651. g_A *= f
  1652. g_H *= f
  1654. x = general_compensation_algorithm(V_th, [g_f_target, g_s_target, g_u_target], g_leak, g_Na, g_Kd, g_CaT, g_CaS, g_KCa, g_A, g_H, params['E_Na'], params['E_K'], params['E_H'], params['E_leak'], params['E_Ca'], params['alpha_Ca'], params['beta_Ca'], params['tau_Ca'], default_g_CaS_for_Ca=default_g_CaS_for_Ca, default_g_CaT_for_Ca=default_g_CaT_for_Ca)
  1656. return x
  1658. def modulate_population(population, V_th, g_f_target, g_s_target, g_u_target, params, set_to_compensate, default_g_CaS_for_Ca = 10., default_g_CaT_for_Ca = 6.0, iterations=0):
  1659. """
  1660. Modulates a population of neurons to achieve specified dynamic input conductances (DICs).
  1662. Parameters
  1663. ----------
  1664. population : array-like
  1665. Array of neuron conductances to be modulated.
  1666. V_th : float
  1667. Threshold voltage for dynamic input conductances (DICs).
  1668. g_f_target : float or None
  1669. Target fast DIC. If None, it will be compensated.
  1670. g_s_target : float or None
  1671. Target slow DIC. If None, it will be compensated.
  1672. g_u_target : float or None
  1673. Target ultra-slow DIC. If None, it will be compensated.
  1674. params : dict
  1675. Dictionary of fixed neuron parameters (e.g., reversal potentials and calcium dynamics).
  1676. set_to_compensate : list
  1677. List of conductances to be compensated (e.g., ['Na', 'Kd', 'CaT']).
  1678. default_g_CaS_for_Ca : float or array-like, optional
  1679. Default S-type calcium conductance for calcium dynamics.
  1680. default_g_CaT_for_Ca : float or array-like, optional
  1681. Default T-type calcium conductance for calcium dynamics.
  1682. iterations : int, optional
  1683. Number of iterations for the compensation algorithm. Useless if g_CaS and g_CaT are not among the conductances to compensate.
  1685. Returns
  1686. -------
  1687. np.ndarray
  1688. Array of modulated neuron conductances.
  1690. Notes
  1691. -----
  1692. - The function uses a compensation algorithm to adjust the conductances of the neurons to achieve the specified DICs.
  1693. - The number of conductances to compensate should be equal to the number of target DICs that are None or NaN.
  1695. - If g_CaS and g_CaT are among the conductances to be compensated, the system is non-linear and providing a default value for the conductances is necessary. The compensation is, in this case, approximate and the reached DICs may not be exactly the target DICs. Iterations can be used to refine the compensation.
  1697. References
  1698. ----------
  1699. Drion, G., Franci, A., Dethier, J., & Sepulchre, R. (2015). Dynamic Input Conductances Shape Neuronal Spiking. eNeuro, 2(1), ENEURO.0031-14.2015. https://doi.org/10.1523/ENEURO.0031-14.2015
  1700. """
  1701. population = np.asarray(population)
  1703. number_none_dics = 0
  1704. if g_f_target is None or np.isnan(g_f_target):
  1705. number_none_dics += 1
  1706. if g_s_target is None or np.isnan(g_s_target):
  1707. number_none_dics += 1
  1708. if g_u_target is None or np.isnan(g_u_target):
  1709. number_none_dics += 1
  1711. if 3 - number_none_dics != len(set_to_compensate):
  1712. raise ValueError('Number of conductances to compensate should be equal to the number of target DICS')
  1714. while iterations >= 0:
  1715. if 'Na' in set_to_compensate:
  1716. population[:, 0] = np.nan
  1717. if 'Kd' in set_to_compensate:
  1718. population[:, 1] = np.nan
  1719. if 'CaT' in set_to_compensate:
  1720. population[:, 2] = np.nan
  1721. if 'CaS' in set_to_compensate:
  1722. population[:, 3] = np.nan
  1723. if 'KCa' in set_to_compensate:
  1724. population[:, 4] = np.nan
  1725. if 'A' in set_to_compensate:
  1726. population[:, 5] = np.nan
  1727. if 'H' in set_to_compensate:
  1728. population[:, 6] = np.nan
  1730. if 'CaS' not in set_to_compensate and 'CaT' not in set_to_compensate:
  1731. iterations = 0
  1733. population = general_compensation_algorithm(V_th, [g_f_target, g_s_target, g_u_target], population[:, 7], population[:, 0], population[:, 1], population[:, 2], population[:, 3], population[:, 4], population[:, 5], population[:, 6], params['E_Na'], params['E_K'], params['E_H'], params['E_leak'], params['E_Ca'], params['alpha_Ca'], params['beta_Ca'], params['tau_Ca'], default_g_CaS_for_Ca=default_g_CaS_for_Ca, default_g_CaT_for_Ca=default_g_CaT_for_Ca)
  1734. iterations -= 1
  1735. default_g_CaT_for_Ca = population[:, 2].copy()
  1736. default_g_CaS_for_Ca = population[:, 3].copy()
  1738. return population
  1740. def general_compensation_algorithm(V_th, target_DICs, new_g_leak, new_g_Na, new_g_Kd, new_g_CaT, new_g_CaS, new_g_KCa, new_g_A, new_g_H, E_Na, E_K, E_H, E_leak, E_Ca, alpha_Ca, beta_Ca, tau_Ca, tau_f_stg = tau_m_Na, tau_s_stg = tau_m_Kd, tau_u_stg = tau_m_H, default_g_CaS_for_Ca = 10., default_g_CaT_for_Ca = 6.0):
  1741. """
  1742. General compensation algorithm to adjust conductances to achieve target dynamic input conductances (DICs).
  1744. Parameters
  1745. ----------
  1746. V_th : float or array-like
  1747. Threshold voltage for dynamic input conductances (DICs).
  1748. target_DICs : list of float
  1749. List containing the target fast, slow, and ultra-slow DICs. Use None or NaN for the DICs to be compensated.
  1750. new_g_leak, new_g_Na, new_g_Kd, new_g_CaT, new_g_CaS, new_g_KCa, new_g_A, new_g_H : array-like or float
  1751. Initial conductances for the respective ion channels. Use None or NaN for the conductances to be compensated.
  1752. E_Na, E_K, E_H, E_leak, E_Ca : float
  1753. Reversal potentials for the respective ion channels.
  1754. alpha_Ca : float
  1755. Proportionality constant for calcium influx.
  1756. beta_Ca : float
  1757. Rate constant for calcium extrusion.
  1758. tau_Ca : float
  1759. Time constant for calcium concentration dynamics.
  1760. tau_f_stg, tau_s_stg, tau_u_stg : callable, optional
  1761. Functions to compute the time constants for fast, slow, and ultra-slow dynamics.
  1762. default_g_CaS_for_Ca : float or array-like, optional
  1763. Default S-type calcium conductance for calcium dynamics.
  1764. default_g_CaT_for_Ca : float or array-like, optional
  1765. Default T-type calcium conductance for calcium dynamics.
  1767. Returns
  1768. -------
  1769. np.ndarray
  1770. Array of compensated conductances.
  1772. Notes
  1773. -----
  1774. - The number of conductances to compensate should be equal to the number of target DICs that are None or NaN.
  1775. - If g_CaS and g_CaT are among the conductances to be compensated, the system is non-linear and providing a default value for the conductances is necessary. The compensation is, in this case, approximate and the reached DICs may not be exactly the target DICs. Iterative compensation can be used to refine the compensation.
  1777. References
  1778. ----------
  1779. Drion, G., Franci, A., Dethier, J., & Sepulchre, R. (2015). Dynamic Input Conductances Shape Neuronal Spiking. eNeuro, 2(1), ENEURO.0031-14.2015. https://doi.org/10.1523/ENEURO.0031-14.2015
  1780. """
  1781. none_index = []
  1783. new_g_dict = {
  1784. 'Na': new_g_Na,
  1785. 'Kd': new_g_Kd,
  1786. 'CaT': new_g_CaT,
  1787. 'CaS': new_g_CaS,
  1788. 'KCa': new_g_KCa,
  1789. 'A': new_g_A,
  1790. 'H': new_g_H,
  1791. 'leak': new_g_leak
  1792. }
  1794. # Handle None and NaN values
  1795. for idx, (key, g_value) in enumerate(new_g_dict.items()):
  1796. g_value = np.atleast_1d(g_value)
  1797. if g_value is None or np.isnan(g_value).all():
  1798. none_index.append(idx)
  1799. new_g_dict[key] = np.zeros_like(g_value) if key not in ['CaT', 'CaS'] else np.full_like(g_value, default_g_CaT_for_Ca if key == 'CaT' else default_g_CaS_for_Ca)
  1801. if len(none_index) == 0 or len(none_index) > 3:
  1802. raise ValueError('Number of conductances to compensate should be between 1 and 3')
  1804. target_DICs = np.array(target_DICs, dtype=np.float64)
  1805. not_none_index_target = np.where(~np.isnan(target_DICs))[0]
  1806. target_DICs = target_DICs[not_none_index_target]
  1808. # verify if the number of none index is equal to the number of none in target_DICs
  1809. if len(none_index) != len(not_none_index_target):
  1810. raise ValueError('Number of None in target_DICs should be equal to the number of None in the conductances')
  1812. new_gs = np.asarray([new_g_dict[key] for key in new_g_dict.keys()])
  1813. copy_new_gs = new_gs.copy().T
  1815. if not isinstance(V_th, np.ndarray):
  1816. V_th = np.array([V_th,])
  1818. S_full = sensitivity_matrix(V_th, *new_gs, E_Na, E_K, E_H, E_leak, E_Ca, alpha_Ca, beta_Ca, tau_Ca, tau_f_stg, tau_s_stg, tau_u_stg)
  1820. S_full = S_full.squeeze()
  1821. if S_full.ndim == 2:
  1822. S_full = S_full[np.newaxis, :, :]
  1824. not_none_index = [i for i in range(len(new_g_dict)) if i not in none_index]
  1826. S_random = S_full[:, not_none_index_target, :][:, :, not_none_index]
  1827. S_compensated = S_full[:, not_none_index_target, :][:, :, none_index]
  1829. new_gs = new_gs[not_none_index]
  1831. new_gs = new_gs.T[:, np.newaxis, :]
  1832. A = S_compensated
  1833. result_dot_product = np.sum(S_random * new_gs, axis=2)
  1834. b = target_DICs[np.newaxis, :] - result_dot_product
  1836. # add a dimension to b - I have to introduce this because of the new version of numpy ...
  1837. b = b[:, :, np.newaxis]
  1839. try:
  1840. x = np.linalg.solve(A, b)
  1841. except np.linalg.LinAlgError:
  1842. x = np.full((1, len(none_index)), -np.inf)
  1844. x = x.squeeze()
  1846. # refill new_gs with the new values
  1847. copy_new_gs[:, none_index] = x
  1848. return copy_new_gs
  1850. def generate_spiking_population(n_cells, V_th=-51.0):
  1851. """
  1852. Generates a population of spiking neurons (corresponding to g_f = -6.2, g_s = 4.0, g_u = 5.0).
  1854. Parameters
  1855. ----------
  1856. n_cells : int
  1857. Number of neurons to generate.
  1858. V_th : float, optional
  1859. Threshold voltage for dynamic input conductances (DICs). Default is -51.0 mV.
  1861. Returns
  1862. -------
  1863. np.ndarray
  1864. Array of generated spiking neuron conductances, shape (n_cells, 8).
  1866. Notes
  1867. -----
  1868. - The population is generated with the following conductance ranges:
  1869. - Leak: [0.007, 0.014]
  1870. - (Na: [2500, 4500]) [compensated]
  1871. - Kd: [70, 140]
  1872. - CaT: [3, 7]
  1873. - CaS: [6, 22]
  1874. - KCa: [140, 180]
  1875. - (A: [200, 400]) [compensated]
  1876. - (H: [0.25, 0.5]) [compensated]
  1877. - Fyon et al. (2024) used a uniform distribution to generate the conductances. We use a gamma distribution to generate the conductances in this work. Results are similar.
  1879. References
  1880. ----------
  1881. Fyon, A., Franci, A., Sacré, P., & Drion, G. (2024). Dimensionality reduction of neuronal degeneracy reveals two interfering physiological mechanisms. PNAS Nexus, 3(10), pgae415. https://doi.org/10.1093/pnasnexus/pgae415
  1882. """
  1884. # FROM Fyon et al. 2024
  1885. g_bar_range_Kd2 = [70, 140]
  1886. g_bar_range_CaT2 = [3, 7]
  1887. g_bar_range_CaS2 = [6, 22]
  1888. g_bar_range_KCa2 = [140, 180]
  1889. g_bar_range_leak2 = [0.007, 0.014]
  1891. # g_Na is compensated
  1892. g_bar_range_Kd = [g_bar_range_Kd2[0], g_bar_range_Kd2[1]]
  1893. g_bar_range_CaT = [g_bar_range_CaT2[0], g_bar_range_CaT2[1]]
  1894. g_bar_range_CaS = [g_bar_range_CaS2[0], g_bar_range_CaS2[1]]
  1895. g_bar_range_KCa = [g_bar_range_KCa2[0], g_bar_range_KCa2[1]]
  1896. # g_A is compensated
  1897. # g_H is compensated
  1898. g_bar_range_leak = [g_bar_range_leak2[0], g_bar_range_leak2[1]]
  1900. g_s_spiking = 4.
  1901. g_u_spiking = 5.
  1902. g_f_spiking = -g_s_spiking - 2.2
  1904. PARAMS = get_default_parameters()
  1905. spiking_population = generate_population(n_cells, V_th, g_f_spiking, g_s_spiking, g_u_spiking, g_bar_range_leak, None, g_bar_range_Kd, g_bar_range_CaT, g_bar_range_CaS, g_bar_range_KCa, None, None, params=PARAMS, distribution="gamma")
  1907. return spiking_population
  1909. def generate_neuromodulated_population(n_cells, V_th_target, g_s_target, g_u_target, set_to_compensate=None, clean=True, use_fitted_gCaS = lambda g_s, g_u : 34.12021074772369 -2.3296612301271464*g_s, use_fitted_gCaT = lambda g_s, g_u : 24.6 - 5.14 * g_s, iterations=0):
  1910. """
  1911. Generates a population of neurons with specified dynamic input conductances (DICs) using neuromodulation from a spiking population.
  1913. Parameters
  1914. ----------
  1915. n_cells : int
  1916. Number of neurons to generate.
  1917. V_th_target : float
  1918. Target threshold voltage for dynamic input conductances (DICs).
  1919. g_s_target : float
  1920. Target slow DIC.
  1921. g_u_target : float
  1922. Target ultra-slow DIC.
  1923. set_to_compensate : list, optional
  1924. List of conductances to be compensated (e.g., ['A', 'H']). If None, the best set will be determined automatically from the target DICs and the reachability from a spiking population.
  1925. clean : bool, optional
  1926. If True, remove any neuron with negative conductances.
  1927. use_fitted_gCaS : callable, optional
  1928. Function to determine the default S-type calcium conductance for calcium dynamics based on the target DICs.
  1929. The default function is a linear fit performed in this work and is associated with the best set of conductances to compensate.
  1930. use_fitted_gCaT : callable, optional
  1931. Function to determine the default T-type calcium conductance for calcium dynamics based on the target DICs.
  1932. The default function is a linear fit performed in this work and is associated with the best set of conductances to compensate.
  1933. iterations : int, optional
  1934. Number of iterations for the compensation algorithm. Useless if g_CaS and g_CaT are not among the conductances to compensate.
  1936. Returns
  1937. -------
  1938. np.ndarray
  1939. Array of generated neuromodulated neuron conductances.
  1941. Notes
  1942. -----
  1943. - The function first generates a population of spiking neurons and then modulates them to achieve the specified DICs.
  1944. - The number of conductances to compensate should be equal to the number of target DICs that are None or NaN.
  1945. - If g_CaS and g_CaT are among the conductances to be compensated, the system is non-linear and providing a default value for the conductances is necessary. The compensation is, in this case, approximate and the reached DICs may not be exactly the target DICs. Iterative compensation can be used to refine the compensation.
  1946. - In practice, the method using neuromodulation introduced by A. Fyon et al. (2015) and improved by this work can be used to generate a population of neurons with specified DICs without the failure of the direct population generation method.
  1948. References
  1949. ----------
  1950. - Fyon, A., Franci, A., Sacré, P., & Drion, G. (2024). Dimensionality reduction of neuronal degeneracy reveals two interfering physiological mechanisms. PNAS Nexus, 3(10), pgae415. https://doi.org/10.1093/pnasnexus/pgae415
  1951. """
  1952. g_f_target = None
  1953. spiking_population = generate_spiking_population(n_cells, V_th_target)
  1955. if set_to_compensate is None:
  1956. set_to_compensate = get_best_set(g_s_target, g_u_target)
  1958. if use_fitted_gCaS:
  1959. d_gCaS = use_fitted_gCaS(g_s_target, g_u_target)
  1960. else:
  1961. d_gCaS = 10.
  1963. if use_fitted_gCaT:
  1964. d_gCaT = use_fitted_gCaT(g_s_target, g_u_target)
  1965. else:
  1966. d_gCaT = 6.0
  1968. neuromodulated_population = modulate_population(spiking_population, V_th_target, g_f_target, g_s_target, g_u_target, get_default_parameters(), set_to_compensate, default_g_CaS_for_Ca=d_gCaS, default_g_CaT_for_Ca=d_gCaT, iterations=iterations)
  1970. if clean:
  1971. # remove any neuron with < 0 conductances
  1972. neuromodulated_population = neuromodulated_population[np.all(neuromodulated_population >= 0, axis=1)]
  1974. return neuromodulated_population

stg.py at commit 42c572c, no license · at the source

Overview

Authors: Julien Brandoit1, Damien Ernst1,2, Guillaume Drion1, Arthur Fyon1
ORCID iDs: Julien Brandoit, Arthur Fyon
  1. Department of Electrical Engineering and Computer Science, University of Liège, Liège, Belgium
  2. LTCI, Telecom Paris, Institut Polytechnique de Paris, Palaiseau, France
Institutions: University of Liège (Belgium); Télécom Paris (France); Institut Polytechnique de Paris (France); Laboratoire Traitement et Communication de l’Information (France)
Journal: PLoS computational biology, volume 22, issue 5, article e1014337
Dates: received 17 September 2025; accepted 15 May 2026; published online 21 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pcbi.1014337 · PMID 42166500 · PMCID PMC13241015 · OpenAlex W4415316636
Open access: gold, a free copy (OpenAlex)
Preprint: arxiv:2509.12783
Status: code verified
Categories: computational modeling (no new data) (modality), human (organism), computational (subfield)
Methods: Source localization, Machine learning, Single-unit activity, calcium imaging
MeSH: Action Potentials*, Models, Neurological*, Neurons*, Algorithms, Animals, Computational Biology, Computer Simulation, Deep Learning, Humans, Ion Channels, Neural Networks, Computer (* major topic)
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: Belgian Federal Science Policy Office (NEMODEI2); Fonds De La Recherche Scientifique - FNRS (ASP-REN40024838); Fonds pour la Formation à la Recherche dans l’Industrie et dans l’Agriculture (FRIA40038025)
Citations: cited by 2 papers (Europe PMC); 58 references in the paper

Abstract

Inferring the biophysical parameters of conductance-based models (CBMs) from experimentally accessible recordings remains a central challenge in computational neuroscience. Spike times are the most widely available data, yet they reveal little about which combinations of ion channel conductances generate the observed activity. This inverse problem is further complicated by neuronal degeneracy, where multiple distinct conductance sets yield similar spiking patterns. We introduce a method that addresses this challenge by combining deep learning with Dynamic Input Conductances (DICs), a theoretical framework that reduces complex CBMs to three interpretable feedback components governing excitability and firing patterns. Our approach first maps spike times directly to DIC densities at threshold using a lightweight neural network that learns a low-dimensional representation of neuronal activity. The predicted DIC values are then used to generate degenerate CBM populations via an iterative compensation algorithm, ensuring compatibility with the intermediate target DICs, and thereby reproducing the corresponding firing patterns, even in high-dimensional models. Applied to two neuronal models, this algorithmic pipeline reconstructs spiking, bursting, and irregular regimes with high accuracy and robustness to variability, including spike trains generated under noisy current injection mimicking physiological stochasticity. It produces diverse degenerate populations within milliseconds on standard hardware, enabling scalable and efficient inference from spike recordings alone. Beyond methodological advances, we provide an open-source software package with a graphical interface that allows experimentalists to generate and explore CBM populations directly from spike trains without requiring programming expertise. Together, this work positions DICs as a practical and interpretable link between experimentally observed activity and mechanistic models. By enabling fast and scalable reconstruction of degenerate populations directly from spike times, our approach provides a powerful way to investigate how neurons exploit conductance variability to achieve reliable computation and provides the foundation for experimental applications that span from neuromodulation studies to real-time model-guided interventions.

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

Repositories

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

Zenodo 16912160

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

julienbrandoit/automatic-degenerate-cbm

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 42c572c9f0ca4abbdef15b6201468c9fd0323ed6, 16 September 2025
Languages: Jupyter (23), Python (18)
Size: 247 files, 41 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README, 23 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (40 files), pandas (26 files), Matplotlib (23 files), seaborn (20 files), PyTorch (13 files), SciPy (12 files), scikit-learn (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
42 files
At the source: github.com/julienbrandoit/automatic-degenerate-cbm

julienbrandoit/Spike2Pop---Bridging-Experimental-Neuroscience-and-Computational-Modeling

License: MIT
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: a2400d09a1b9ff1528ebdd6b2d8c47d2380fce0e, 22 August 2025
Languages: Python (9)
Size: 27 files, 9 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (7 files), pandas (5 files), PyTorch (3 files), SciPy (2 files), Matplotlib (1 file), Pillow (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
11 files
At the source: github.com/julienbrandoit/Spike2Pop---Bridging-Experimental-Neuroscience-and-Computational-Modeling

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:

  • 3 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 50 scripts, each with its path and the digest of its content;
  • 14 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability

All training and testing datasets are available on Zenodo at link https://doi.org/10.5281/zenodo.16912160 All code used for running experiments, model fitting, and plotting is available on a GitHub repository at https://github.com/julienbrandoit/automatic-degenerate-cbm The GUI software package is available on a GitHub repository at https://github.com/julienbrandoit/Spike2Pop---Bridging-Experimental-Neuroscience-and-Computational-Modeling.

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

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 11 MeSH terms, 3 funders, 53 references.

Cite

This paper

Brandoit, J., Ernst, D., Drion, G., & Fyon, A. (2026). Fast reconstruction of degenerate populations of conductance-based neuron models from spike times. PLoS computational biology, 22(5), e1014337. https://doi.org/10.1371/journal.pcbi.1014337

BibTeX

@article{brandoit2026fast,
author = {Brandoit, Julien and Ernst, Damien and Drion, Guillaume and Fyon, Arthur},
title = {{Fast reconstruction of degenerate populations of conductance-based neuron models from spike times}},
journal = {PLoS computational biology},
year = {2026},
month = may,
volume = {22},
number = {5},
pages = {e1014337},
publisher = {PLOS},
issn = {1553-734X},
doi = {10.1371/journal.pcbi.1014337},
url = {https://doi.org/10.1371/journal.pcbi.1014337},
pmid = {42166500},
pmcid = {PMC13241015}
}

RIS

TY - JOUR
AU - Brandoit, Julien
AU - Ernst, Damien
AU - Drion, Guillaume
AU - Fyon, Arthur
TI - Fast reconstruction of degenerate populations of conductance-based neuron models from spike times
T2 - PLoS computational biology
J2 - PLoS Comput Biol
PY - 2026
DA - 2026/05/21
VL - 22
IS - 5
SP - e1014337
SN - 1553-734X
PB - PLOS
DO - 10.1371/journal.pcbi.1014337
UR - https://doi.org/10.1371/journal.pcbi.1014337
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pcbi.1014337",
"type": "article-journal",
"title": "Fast reconstruction of degenerate populations of conductance-based neuron models from spike times",
"container-title": "PLoS computational biology",
"author": [
{
"family": "Brandoit",
"given": "Julien"
},
{
"family": "Ernst",
"given": "Damien"
},
{
"family": "Drion",
"given": "Guillaume"
},
{
"family": "Fyon",
"given": "Arthur"
}
],
"container-title-short": "PLoS Comput Biol",
"volume": "22",
"issue": "5",
"page": "e1014337",
"DOI": "10.1371/journal.pcbi.1014337",
"PMID": "42166500",
"PMCID": "PMC13241015",
"ISSN": "1553-734X",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pcbi.1014337",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
21
]
]
}
}

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.1014177 [code]
Activity-dependent neuromodulation and calcium homeostasis cooperate to produce robust and modulable neuronal function.
Journal: PLoS computational biology
In common: computational modeling (no new data), 13 references, author Arthur Fyon
[2] doi:10.1371/journal.pcbi.1014458 [code]
Neuronal excitability and parameter variability in the Hodgkin-Huxley model.
Journal: PLoS computational biology
In common: pandas, SciPy, Matplotlib, 1 other tool, computational, computational modeling (no new data), 7 references
[3] doi:10.1007/s10827-026-00936-7 [code]
When can neuronal activity-dependent homeostatic plasticity maintain circuit-level properties?
Journal: Journal of computational neuroscience
In common: Matplotlib, NumPy, computational modeling (no new data), 6 references
[4] doi:10.1016/j.isci.2026.115488 [code]
An integrated &lt;i&gt;i&lt;/i&gt; &lt;i&gt;n vitro&lt;/i&gt; platform and biophysical modeling approach for studying synaptic transmission in isolated neuronal pairs.
Journal: iScience
In common: Pillow, PyTorch, seaborn, 5 other tools, 2 references
[5] doi:10.1126/sciadv.aee9425 [code]
Probabilistic inference of homonymous and heteronymous recurrent inhibition in human muscles from large-scale motor neuron recordings.
Journal: Science advances
In common: PyTorch, seaborn, scikit-learn, 4 other tools, 3 references
[6] doi:10.3389/frai.2026.1771088 [code]
Few-shot deployment of pretrained MRI transformers in brain imaging tasks.
Journal: Frontiers in artificial intelligence
In common: Pillow, PyTorch, seaborn, 5 other tools, 2 references
[7] doi:10.1016/j.patter.2026.101538 [code]
A multi-modal foundation model for brain disease diagnosis and medical imaging.
Journal: Patterns (New York, N.Y.)
In common: Pillow, PyTorch, scikit-learn, 4 other tools, 2 references
[8] doi:10.1093/bioinformatics/btag652 [code]
mmVelo: a deep generative model for estimating cell state-dependent dynamics across multiple modalities.
Journal: Bioinformatics (Oxford, England)
In common: Pillow, PyTorch, seaborn, 5 other tools, 1 reference
[9] doi:10.1162/imag.a.1164 [code]
Bias and generalizability of brain age prediction models: A multi-cohort evaluation with anatomical and interpretability insights.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Pillow, PyTorch, seaborn, 5 other tools, 1 reference
[10] doi:10.1098/rstb.2024.0461 [code]
Shallow recurrent decoders for neural and behavioural dynamics.
Journal: Philosophical transactions of the Royal Society of London. Series B, Biological sciences
In common: Pillow, PyTorch, scikit-learn, 4 other tools, computational, 1 reference

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

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

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

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