OSCR

NMDA receptor kinetics drive distinct routes to chaotic firing in pyramidal neurons.

Code ↔ Paper

26 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 26 matches
  1. [1] § Materials and methods › Data analysis pipeline › Frequency band classification ↔ v_StimW_Beta_5_6_1D.py, lines 399–414 · score 0.96 · 7–33 ms, 33–77 ms, 77–125 ms, 125–250 ms, 250–2000 ms, 0.5–4 Hz
  2. [2] § Materials and methods › Data analysis pipeline › Frequency band classification ↔ v_StimW_Beta_5_6_1D_CaMKII_6_5_GABA_4_Complete_6.py, lines 516–531 · score 0.96 · 7–33 ms, 33–77 ms, 77–125 ms, 125–250 ms, 250–2000 ms, 0.5–4 Hz
  3. [3] § Materials and methods › Initial conditions robustness ↔ Sensitivity.py, lines 1–19 · score 0.91 · phase plane trajectories, mAMPA, mGABA, mNMDA, Steady state, synaptic gating
  4. [4] § Results › Differential excitatory and inhibitory neuron responses ↔ reviewer_response_simulations_final.py, lines 635–778 · score 0.82 · Conductance tuning, balance heatmaps, Glutamatergic scaling, inhibitory neuron, firing rate, GABAergic
  5. [5] § Results › Oscillatory band evolution and spectral consolidation ↔ v4_6_8_CaMKII_3.py, lines 90–144 · score 0.79 · neural oscillations, 0.5–4 Hz, 13–30 Hz, 8–13 Hz, 30–100 Hz, Frequency band
  6. [6] § Results › Oscillatory band evolution and spectral consolidation ↔ v4_6.py, lines 636–666 · score 0.78 · 0.5–4 Hz, 13–30 Hz, 8–13 Hz, 30–100 Hz, Frequency band, 4–8 Hz
  7. [7] § Results › Oscillatory band evolution and spectral consolidation ↔ v4_6_8_CaMKII_3.py, lines 90–144 · score 0.76 · 0.5–4 Hz, 13–30 Hz, 8–13 Hz, 30–100 Hz, Frequency band, 4–8 Hz
  8. [8] § Results › Oscillatory band evolution and spectral consolidation ↔ v4_6.py, lines 636–666 · score 0.74 · 0.5–4 Hz, 13–30 Hz, 8–13 Hz, 30–100 Hz, Frequency band, 4–8 Hz
  9. [9] § Materials and methods › Data analysis pipeline › Dynamical analysis measures › Maximum Lyapunov exponents ↔ v_StimW_Beta_5_6_1D.py, lines 141–175 · score 0.72 · nearest neighbor, KD tree, Lyapunov exponent, efficient, computationally
  10. [10] § Results ↔ v4_6_8_CaMKII_sub_21.py, lines 569–640 · score 0.70 · GABA mediated modulation, state transitions, Parameter space, plasticity region, E1, GABA frequencies
  11. [11] § Materials and methods › Data analysis pipeline › Bifurcation diagram construction ↔ gaba_bifurcation.py, lines 185–232 · score 0.67 · inter spike intervals, chaotic regions, branches, background, scatter, Bifurcation
  12. [12] § Materials and methods › Data analysis pipeline › Bifurcation diagram construction ↔ v_StimW_Beta_5_6_1D.py, lines 141–175 · score 0.65 · nearest neighbor, KD tree, Lyapunov exponent, simulations
  13. [13] § Results › Bifurcation analysis reveals period-doubling routes to chaos ↔ v4_6.py, lines 745–852 · score 0.63 · bifurcation diagram, frequency bands, Lyapunov exponent, chaotic regions, chaos, alpha
  14. [14] § Materials and methods › Data analysis pipeline › Dynamical analysis measures › Maximum Lyapunov exponents ↔ v4_6.py, lines 745–852 · score 0.61 · chaotic dynamics, Lyapunov exponent, bifurcation diagrams, chaos, ISI
  15. [15] § Materials and methods › Neuronal model overview ↔ v_StimW_Beta_5_6_1D_CaMKII_6_5_GABA_4_Complete_6.py, lines 286–331 · score 0.60 · isAHP, GABAergic, IKdr, IL, INaP, ICa
  16. [16] § Results › Differential excitatory and inhibitory neuron responses ↔ reviewer_response_simulations_final.py, lines 1–36 · score 0.60 · excitatory neuron, inhibitory neuron, GABAergic, sensitivity, glutamate, stimulation
  17. [17] § Materials and methods › Neuronal model overview ↔ Code1.py, lines 540–580 · score 0.59 · isAHP, GABAergic, IKdr, IL, INaP, ICa
  18. [18] § Results › GABA modulates CaMKII-mediated plasticity states ↔ v4_6_8_CaMKII_3.py, lines 456–521 · score 0.58 · phosphorylation states, dependent plasticity, GABA frequency, classified, CaMKII, space
  19. [19] § Results › GABAergic inhibition provides frequency-selective stabilization ↔ v4_6_8_CaMKII_3.py, lines 315–342 · score 0.58 · moderate inhibitory, low GABA frequencies, 20 ms, modulated, NMDA
  20. [20] § Results › Differential excitatory and inhibitory neuron responses ↔ reviewer_response_simulations_final.py, lines 781–819 · score 0.57 · fast spiking, inhibitory neurons, neuron responses, glutamate, GABA, NMDA
  21. [21] § Results › Model validation and basic neuronal response ↔ Code1.py, lines 582–656 · score 0.57 · spike detection, post transient, membrane potential, CaMKII, models, AMPA
  22. [22] § Results › Bifurcation analysis reveals period-doubling routes to chaos ↔ gaba_bifurcation.py, lines 131–183 · score 0.55 · bifurcation diagram, Lyapunov exponent, chaotic regions, doubled, scatter, intervals
  23. [23] § Results › Model validation and basic neuronal response ↔ v4_6_8_CaMKII_3.py, lines 736–779 · score 0.54 · transition zone, LTP thresholds, LTD, modulation, plasticity, calcium
  24. [24] § Materials and methods › Synaptic currents › NMDA receptors ↔ reviewer_response_simulations_final.py, lines 124–168 · score 0.53 · decay rate, NMDA receptor, intrinsic, neurons, synaptic
  25. [25] § Materials and methods › Data analysis pipeline › Statistical analysis and validation ↔ reviewer_response_simulations_final.py, lines 1–36 · score 0.53 · response curves, parameter sweeps, heatmaps, simulated
  26. [26] § Materials and methods › Data analysis pipeline › Bifurcation diagram construction ↔ gaba_bifurcation.py, lines 185–232 · score 0.50 · inter spike interval, chaotic regions, bifurcation, chaos

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,127 lines · 48 KB · MIT · 5 matches

  1. # -*- coding: utf-8 -*-
  2. """
  3. Excitatory/Inhibitory Neuron Simulations with Parameter Sensitivity
  4. ====================================================================
  6. Supplementary analyses:
  7. 1) Excitatory vs inhibitory neuron responses under different synaptic
  8. conditions (glutamate-only, GABA-only, combined)
  9. 2) β_NMDA parameter sweep across neuron types
  10. 3) Conductance tuning of GABAergic and glutamatergic stimulation
  11. 4) β_AMPA (AMPA receptor closing rate) sensitivity analysis
  13. Figures produced:
  14. Fig 1 — Synaptic input pulse profiles (glutamate, GABA, combined)
  15. Fig 2 — Three stimulation conditions: excitatory neuron
  16. Fig 3 — Excitatory vs inhibitory neuron comparison under combined input
  17. Fig 4 — β_NMDA parameter sweep: firing rate, CaMKII, E/I balance
  18. Fig 5 — Conductance tuning: dose-response curves + E/I balance heatmaps
  19. Fig 6 — Three stimulation conditions: inhibitory neuron
  20. Fig 7 — β_NMDA × neuron type × input condition interaction matrix
  21. Fig 8 — β_AMPA sensitivity: EPSC waveforms, firing rate, CaMKII,
  22. E/I ratio, and β_AMPA × β_NMDA interaction heatmaps
  24. Core simulation model: UNCHANGED from original manuscript.
  26. @author: Mehdi Borjkhani
  27. """
  29. import numpy as np
  30. import matplotlib.pyplot as plt
  31. import matplotlib.gridspec as gridspec
  32. from matplotlib.lines import Line2D
  33. import os
  34. import warnings
  36. warnings.filterwarnings('ignore')
  38. # == Publication-quality settings =============================================
  39. plt.rcParams.update({
  40. 'font.size': 11,
  41. 'axes.titlesize': 12,
  42. 'axes.labelsize': 11,
  43. 'xtick.labelsize': 9,
  44. 'ytick.labelsize': 9,
  45. 'legend.fontsize': 9,
  46. 'font.family': 'serif',
  47. 'font.serif': ['Times New Roman', 'DejaVu Serif', 'serif'],
  48. 'mathtext.fontset': 'dejavuserif',
  49. 'axes.linewidth': 1.0,
  50. 'axes.spines.top': False,
  51. 'axes.spines.right': False,
  52. 'lines.linewidth': 1.5,
  53. 'legend.frameon': False,
  54. 'axes.grid': False,
  55. 'savefig.dpi': 300,
  56. 'savefig.bbox': 'tight',
  57. })
  59. save_directory = os.path.join(os.getcwd(), "supplementary_figures")
  60. os.makedirs(save_directory, exist_ok=True)
  62. # =============================================================================
  63. # MODEL PARAMETERS (Same as original manuscript)
  64. # =============================================================================
  65. # Stimulation parameters
  66. FreqS = 10; StimW = 1000 / FreqS; App = 1; Puls = 5
  67. FreqS_GABA = 10; StimW_GABA = 1000 / FreqS_GABA; App_GABA = 1; Puls_GABA = 5
  69. # NMDA/AMPA/GABA synapse parameters
  70. alpha_nmda = 0.38
  71. alpha_ampa = 1.1; beta_ampa = 0.67; k_ampa = 1; g_ampa = 0.35
  72. k_nmda = 1
  73. g_gaba = 0.35; alpha_gaba = 5; beta_gaba = 0.18; k_gaba = 1
  75. # Simulation time
  76. Tmax = 8.0e3 # 8 seconds (allows CaMKII to equilibrate)
  77. deltaT = 0.05
  78. tn = int(Tmax / deltaT)
  79. t = np.arange(0, Tmax + deltaT, deltaT)
  81. # Reversal potentials
  82. vna = 55; vk = -90; vl = -70; vca = 120; E_gaba = -80
  84. # Conductance densities (excitatory pyramidal neuron defaults)
  85. gna_exc = 50; gnap_exc = 0.2; gl_exc = 0.1; gkdr_exc = 10
  86. gA_exc = 1; gM_exc = 0.5; gca_exc = 0.1; gc_exc = 5; gsAHP_exc = 2
  88. # Gating parameters
  89. theta_m, tau_m = -35, 8
  90. theta_h, tau_h = -50, -6
  91. theta_n, tau_n = -40, 12
  92. theta_a, tau_a = -55, 15
  93. theta_z, tau_z = -45, 4
  94. theta_p, tau_p = -40, 2
  95. theta_ht, theta_nt = -45, -30
  96. theta_b, tau_b = -85, -5
  97. theta_r, tau_r = -25, 8
  98. theta_c, tau_c = -35, 5
  100. # CaMKII parameters
  101. k_I = 0.001; nh1 = 3; vCaN = 2.0e-3; kM = 20
  102. kh1 = 4.0; kh2 = 0.7; vPKA = 0.45e-3; kpk = 0.0059
  103. ck1 = 0.5e-3; ck2 = 10e-3; ck3 = 1e-3; ck4 = 1e-6
  105. # --- Inhibitory interneuron parameter rationale ---
  106. # Fast-spiking (FS) interneurons differ from pyramidal neurons in:
  107. # - Higher g_Kdr: faster repolarisation -> narrower action potentials
  108. # (cf. Erisir et al., 1999, J Neurophysiol)
  109. # - Reduced g_M, g_sAHP: minimal spike-frequency adaptation
  110. # (cf. Bhatt et al., 2019, Front Cell Neurosci)
  111. # - Lower g_NaP: reduced persistent Na+ current
  112. # - Lower g_Ca, g_C: weaker Ca2+-dependent K+ coupling
  113. # - Slightly higher g_leak: lower input resistance
  114. # - Higher resting potential (~-65 mV vs -70 mV)
  115. # These adjustments produce the characteristic FS phenotype: high-frequency
  116. # non-adapting firing with narrow spikes.
  117. INHIBITORY_PARAMS = {
  118. 'gna': 55, 'gnap': 0.05, 'gl': 0.15, 'gkdr': 15,
  119. 'gA': 0.3, 'gM': 0.1, 'gca': 0.02, 'gc': 1.0, 'gsAHP': 0.2,
  120. 'v_init': -65,
  121. }
  124. # =============================================================================
  125. # CORE SIMULATION FUNCTION (parameterised -- UNCHANGED from original)
  126. # =============================================================================
  127. def run_neuron_simulation(beta_nmda, neuron_type='excitatory',
  128. enable_glutamate=True, enable_gaba=True,
  129. g_ampa_scale=1.0, g_gaba_scale=1.0,
  130. g_nmda_scale=1.0, beta_ampa_val=None):
  131. """
  132. Run a single neuron simulation under specified conditions.
  134. Parameters
  135. ----------
  136. beta_nmda : float
  137. NMDA receptor decay rate parameter.
  138. neuron_type : str
  139. 'excitatory' or 'inhibitory'. Adjusts intrinsic conductances.
  140. enable_glutamate : bool
  141. Enable glutamatergic (AMPA + NMDA) input.
  142. enable_gaba : bool
  143. Enable GABAergic input.
  144. g_ampa_scale, g_gaba_scale, g_nmda_scale : float
  145. Scaling factors for tuning synaptic conductances.
  146. beta_ampa_val : float or None
  147. AMPA receptor closing (unbinding) rate. If None, uses global beta_ampa.
  149. Returns
  150. -------
  151. dict : All recorded traces and metadata.
  152. """
  153. # -- Adjust intrinsic properties by neuron type --
  154. if neuron_type == 'inhibitory':
  155. gna = INHIBITORY_PARAMS['gna']
  156. gnap = INHIBITORY_PARAMS['gnap']
  157. gl = INHIBITORY_PARAMS['gl']
  158. gkdr = INHIBITORY_PARAMS['gkdr']
  159. gA = INHIBITORY_PARAMS['gA']
  160. gM = INHIBITORY_PARAMS['gM']
  161. gca = INHIBITORY_PARAMS['gca']
  162. gc = INHIBITORY_PARAMS['gc']
  163. gsAHP = INHIBITORY_PARAMS['gsAHP']
  164. v_init = INHIBITORY_PARAMS['v_init']
  165. else:
  166. gna = gna_exc; gnap = gnap_exc; gl = gl_exc; gkdr = gkdr_exc
  167. gA = gA_exc; gM = gM_exc; gca = gca_exc; gc = gc_exc; gsAHP = gsAHP_exc
  168. v_init = -70
  170. # Scale synaptic conductances
  171. g_ampa_eff = g_ampa * g_ampa_scale
  172. g_gaba_eff = g_gaba * g_gaba_scale
  173. beta_ampa_eff = beta_ampa_val if beta_ampa_val is not None else beta_ampa
  175. # -- Allocate arrays --
  176. temp = np.zeros(tn + 1)
  177. v = temp.copy(); v[0] = v_init
  178. h = temp.copy(); h[0] = 0.98
  179. n = temp.copy(); n[0] = 0.01
  180. b = temp.copy(); b[0] = 0
  181. z = temp.copy(); z[0] = 0.050
  182. r = temp.copy()
  183. c = temp.copy()
  184. q = temp.copy()
  185. Ca = temp.copy()
  186. m_ampa = temp.copy()
  187. m_nmda = temp.copy()
  188. g_VDpost = temp.copy()
  189. I_ampa = temp.copy()
  190. I_nmda = temp.copy()
  191. m_gaba = temp.copy()
  192. I_gaba = temp.copy()
  193. Iapp = np.zeros(tn)
  194. Iapp_GABA = np.zeros(tn)
  196. # CaMKII
  197. CaMKII_p = np.zeros(tn + 1)
  198. cam = np.zeros(tn + 1)
  199. ep = np.zeros(tn + 1); ep[0] = 1.0
  200. I_camkii = np.zeros(tn + 1)
  201. P0 = 1.0
  202. P1 = P2 = P3 = P4 = P5 = P6 = P7 = P8 = P9 = P10 = 0.0
  204. # -- Integration loop --
  205. for i in range(tn):
  206. # Stimulus pulses
  207. if enable_glutamate and (0 <= (t[i] % StimW) <= Puls):
  208. Iapp[i] = App
  209. if enable_gaba and (0 <= (t[i] % StimW_GABA) <= Puls_GABA):
  210. Iapp_GABA[i] = App_GABA
  212. # Gating steady-states
  213. m_inf = 1 / (1 + np.exp(-(v[i] - theta_m) / tau_m))
  214. h_inf = 1 / (1 + np.exp(-(v[i] - theta_h) / tau_h))
  215. n_inf = 1 / (1 + np.exp(-(v[i] - theta_n) / tau_n))
  216. a_inf = 1 / (1 + np.exp(-(v[i] - theta_a) / tau_a))
  217. b_inf = 1 / (1 + np.exp(-(v[i] - theta_b) / tau_b))
  218. z_inf = 1 / (1 + np.exp(-(v[i] - theta_z) / tau_z))
  219. p_inf = 1 / (1 + np.exp(-(v[i] - theta_p) / tau_p))
  220. r_inf = 1 / (1 + np.exp(-(v[i] - theta_r) / tau_r))
  221. tau_ht_mod = 0.1 + 0.75 / (1 + np.exp(-(v[i] - theta_ht) / (-6)))
  222. tau_nt_mod = 0.1 + 0.5 / (1 + np.exp(-(v[i] - theta_nt) / (-15)))
  223. d_inf = 1 / (1 + 6 / (Ca[i] + 1e-6))
  224. c_inf = 1 / (1 + np.exp(-(v[i] - theta_c) / tau_c))
  225. q_inf = 1 / (1 + (2 / (Ca[i] + 1e-6))**4)
  227. # Gating variable updates
  228. h[i+1] = h[i] + deltaT * ((h_inf - h[i]) / tau_ht_mod)
  229. n[i+1] = n[i] + deltaT * ((n_inf - n[i]) / tau_nt_mod)
  230. b[i+1] = b[i] + deltaT * ((b_inf - b[i]) / 15)
  231. z[i+1] = z[i] + deltaT * ((z_inf - z[i]) / 75)
  232. r[i+1] = r[i] + deltaT * ((r_inf - r[i]) / 1)
  233. c[i+1] = c[i] + deltaT * ((c_inf - c[i]) / 2)
  234. q[i+1] = q[i] + deltaT * ((q_inf - q[i]) / 450)
  236. # Ionic currents
  237. ina = gna * m_inf**3 * h[i] * (v[i] - vna)
  238. inap = gnap * p_inf * (v[i] - vna)
  239. ikdr = gkdr * n[i]**4 * (v[i] - vk)
  240. iA = gA * a_inf**3 * b[i] * (v[i] - vk)
  241. iM = gM * z[i] * (v[i] - vk)
  242. il = gl * (v[i] - vl)
  243. ica = gca * r[i]**2 * (v[i] - vca)
  244. ic = gc * d_inf * c[i] * (v[i] - vk)
  245. isAHP = gsAHP * q[i] * (v[i] - vk)
  247. # Synaptic currents
  248. G_syn = Iapp[i]
  249. m_ampa[i+1] = m_ampa[i] + deltaT * (alpha_ampa * G_syn * (1 - m_ampa[i]) - beta_ampa_eff * m_ampa[i])
  250. I_ampa[i] = k_ampa * g_ampa_eff * m_ampa[i] * (v[i] - 55)
  252. m_nmda[i+1] = m_nmda[i] + deltaT * (alpha_nmda * G_syn * (1 - m_nmda[i]) - beta_nmda * m_nmda[i])
  253. g_conpost = g_VDpost[i] + 1
  254. g_VDinfpost = 0.0007 * (v[i] - (-100))
  255. g_VDpost[i+1] = g_VDpost[i] + deltaT * ((g_VDinfpost - g_VDpost[i]) / 0.05)
  256. MG = 1 / (1 + (1.4 / 3.75) * np.exp(-0.062 * v[i]))
  257. I_nmda[i] = k_nmda * g_nmda_scale * g_conpost * MG * m_nmda[i] * (v[i] - 55)
  259. G_syn_gaba = Iapp_GABA[i]
  260. m_gaba[i+1] = m_gaba[i] + deltaT * (alpha_gaba * G_syn_gaba * (1 - m_gaba[i]) - beta_gaba * m_gaba[i])
  261. I_gaba[i] = k_gaba * g_gaba_eff * m_gaba[i] * (v[i] - E_gaba)
  263. # Calcium update
  264. Ca[i+1] = Ca[i] + deltaT * (-0.13 * ica - 0.0012 * I_ampa[i] - 0.012 * I_nmda[i] - Ca[i] / 13)
  266. # CaMKII
  267. Ca_uM = Ca[i] * 1
  268. num_cv1 = 10 * ck1 * P0 * (Ca_uM / kh1) ** (2 * nh1)
  269. den_cv1 = (1 + (Ca_uM / kh1) ** nh1) ** 2
  270. cv1 = num_cv1 / den_cv1
  271. cv2 = ck1 * (Ca_uM / kh1) ** nh1 / (1 + (Ca_uM / kh1) ** nh1)
  272. sum_k_Pk = P1 + 2*P2 + 3*P3 + 4*P4 + 5*P5 + 6*P6 + 7*P7 + 8*P8 + 9*P9 + 10*P10
  273. cv3 = ck2 * ep[i] / (kM + sum_k_Pk)
  275. P0 += deltaT * (-cv1 + cv3 * P1)
  276. P1 += deltaT * (cv1 - cv3*P1 - cv2*P1 + 2*cv3*P2)
  277. P2 += deltaT * (cv2*P1 - 2*cv3*P2 - 1.8*cv2*P2 + 3*cv3*P3)
  278. P3 += deltaT * (1.8*cv2*P2 - 3*cv3*P3 - 2.3*cv2*P3 + 4*cv3*P4)
  279. P4 += deltaT * (2.3*cv2*P3 - 4*cv3*P4 - 2.7*cv2*P4 + 5*cv3*P5)
  280. P5 += deltaT * (2.7*cv2*P4 - 5*cv3*P5 - 2.8*cv2*P5 + 6*cv3*P6)
  281. P6 += deltaT * (2.8*cv2*P5 - 6*cv3*P6 - 2.7*cv2*P6 + 7*cv3*P7)
  282. P7 += deltaT * (2.7*cv2*P6 - 7*cv3*P7 - 2.3*cv2*P7 + 8*cv3*P8)
  283. P8 += deltaT * (2.3*cv2*P7 - 8*cv3*P8 - 1.8*cv2*P8 + 9*cv3*P9)
  284. P9 += deltaT * (1.8*cv2*P8 - 9*cv3*P9 - cv2*P9 + 10*cv3*P10)
  285. P10 += deltaT * (cv2*P9 - 10*cv3*P10)
  287. cam[i+1] = P1+P2+P3+P4+P5+P6+P7+P8+P9+P10
  288. total_P = P0+P1+P2+P3+P4+P5+P6+P7+P8+P9+P10
  289. CaMKII_p[i+1] = cam[i+1] / total_P
  291. ep[i+1] = ep[i] + deltaT * (-ck3*I_camkii[i]*ep[i] + ck4*(ep[0]-ep[i]) + k_I*1.0)
  292. I_camkii[i+1] = I_camkii[i] + deltaT * (
  293. -ck3*I_camkii[i]*ep[i] + ck4*(ep[0]-ep[i]) +
  294. vPKA*(1.0/(kpk+1.0)) -
  295. (vCaN*I_camkii[i]*(Ca_uM/kh2)**3 / (1+(Ca_uM/kh2)**3))
  296. )
  298. I_syn = I_ampa[i] + I_nmda[i] + I_gaba[i]
  299. v[i+1] = v[i] + deltaT * (-(ina + inap + ikdr + iA + iM + il + ica + ic + isAHP + I_syn))
  301. # Spike detection
  302. spike_indices = np.where((v[:-1] < 0) & (v[1:] >= 0))[0]
  303. spike_times = t[spike_indices]
  304. firing_rate = len(spike_times) / (Tmax / 1000) if Tmax > 0 else 0
  306. return {
  307. 't': t, 'v': v, 'Ca': Ca, 'CaMKII_p': CaMKII_p,
  308. 'Iapp': Iapp, 'Iapp_GABA': Iapp_GABA,
  309. 'I_ampa': I_ampa, 'I_nmda': I_nmda, 'I_gaba': I_gaba,
  310. 'm_ampa': m_ampa, 'm_nmda': m_nmda, 'm_gaba': m_gaba,
  311. 'spike_times': spike_times, 'firing_rate': firing_rate,
  312. 'neuron_type': neuron_type, 'beta_nmda': beta_nmda,
  313. 'enable_glutamate': enable_glutamate, 'enable_gaba': enable_gaba,
  314. }
  317. # =============================================================================
  318. # ANALYSIS HELPERS
  319. # =============================================================================
  320. def count_spikes(v_trace, t_array, t_start=2000):
  321. """Count spikes after t_start ms to exclude transients."""
  322. idx_start = int(t_start / deltaT)
  323. v_seg = v_trace[idx_start:]
  324. spikes = np.where((v_seg[:-1] < 0) & (v_seg[1:] >= 0))[0]
  325. duration_s = (t_array[-1] - t_start) / 1000
  326. return len(spikes), len(spikes) / duration_s if duration_s > 0 else 0
  329. def compute_ei_ratio(res, t_start=2000):
  330. """
  331. Compute excitation/inhibition (E/I) current ratio in steady state.
  333. E/I ratio = |mean(I_AMPA + I_NMDA)| / |mean(I_GABA)|
  335. Returns np.inf if GABA current is zero.
  336. """
  337. idx_ss = int(t_start / deltaT)
  338. I_exc = np.abs(np.mean(res['I_ampa'][idx_ss:] + res['I_nmda'][idx_ss:]))
  339. I_inh = np.abs(np.mean(res['I_gaba'][idx_ss:]))
  340. if I_inh < 1e-12:
  341. return np.inf if I_exc > 1e-12 else 1.0
  342. return I_exc / I_inh
  345. def compute_mean_currents(res, t_start=2000):
  346. """Compute mean absolute synaptic currents in steady state."""
  347. idx_ss = int(t_start / deltaT)
  348. return {
  349. 'I_ampa': np.mean(np.abs(res['I_ampa'][idx_ss:])),
  350. 'I_nmda': np.mean(np.abs(res['I_nmda'][idx_ss:])),
  351. 'I_gaba': np.mean(np.abs(res['I_gaba'][idx_ss:])),
  352. 'I_exc_total': np.mean(np.abs(res['I_ampa'][idx_ss:] + res['I_nmda'][idx_ss:])),
  353. }
  356. def save_fig(fig, name):
  357. """Save figure in both PNG and PDF format."""
  358. fig.savefig(os.path.join(save_directory, f'{name}.png'), dpi=300)
  359. fig.savefig(os.path.join(save_directory, f'{name}.pdf'), dpi=300)
  360. plt.close(fig)
  363. # =============================================================================
  364. # FIGURE 1: PULSE PROFILES OF GABA AND GLUTAMATE
  365. # =============================================================================
  366. def figure1_pulse_profiles():
  367. """Show the synaptic input pulse profiles clearly."""
  368. print(" Figure 1: Pulse profiles ...")
  369. t_short = np.arange(0, 300 + deltaT, deltaT)
  371. glu_pulse = np.array([App if 0 <= (ti % StimW) <= Puls else 0 for ti in t_short])
  372. gaba_pulse = np.array([App_GABA if 0 <= (ti % StimW_GABA) <= Puls_GABA else 0 for ti in t_short])
  374. fig, axes = plt.subplots(3, 1, figsize=(7, 5.5), sharex=True)
  376. # Glutamate
  377. axes[0].fill_between(t_short, 0, glu_pulse, color='#2E86C1', alpha=0.35, step='mid')
  378. axes[0].step(t_short, glu_pulse, color='#2E86C1', linewidth=1.8, where='mid')
  379. axes[0].set_ylabel('Glutamate\nInput (a.u.)')
  380. axes[0].set_title('(A) Glutamatergic Pulse Profile', loc='left', fontweight='bold')
  381. axes[0].set_ylim(-0.1, 1.5)
  382. axes[0].text(0.75, 0.8, f'Freq = {FreqS} Hz\nWidth = {Puls} ms',
  383. transform=axes[0].transAxes, fontsize=9,
  384. bbox=dict(boxstyle='round', fc='white', alpha=0.8))
  386. # GABA
  387. axes[1].fill_between(t_short, 0, gaba_pulse, color='#E74C3C', alpha=0.35, step='mid')
  388. axes[1].step(t_short, gaba_pulse, color='#E74C3C', linewidth=1.8, where='mid')
  389. axes[1].set_ylabel('GABA\nInput (a.u.)')
  390. axes[1].set_title('(B) GABAergic Pulse Profile', loc='left', fontweight='bold')
  391. axes[1].set_ylim(-0.1, 1.5)
  392. axes[1].text(0.75, 0.8, f'Freq = {FreqS_GABA} Hz\nWidth = {Puls_GABA} ms',
  393. transform=axes[1].transAxes, fontsize=9,
  394. bbox=dict(boxstyle='round', fc='white', alpha=0.8))
  396. # Combined
  397. axes[2].fill_between(t_short, 0, glu_pulse, color='#2E86C1', alpha=0.3, step='mid', label='Glutamate')
  398. axes[2].fill_between(t_short, 0, gaba_pulse, color='#E74C3C', alpha=0.3, step='mid', label='GABA')
  399. axes[2].step(t_short, glu_pulse, color='#2E86C1', linewidth=1.5, where='mid')
  400. axes[2].step(t_short, gaba_pulse, color='#E74C3C', linewidth=1.5, where='mid', linestyle='--')
  401. axes[2].set_ylabel('Input (a.u.)')
  402. axes[2].set_xlabel('Time (ms)')
  403. axes[2].set_title('(C) Combined Input Profile', loc='left', fontweight='bold')
  404. axes[2].set_ylim(-0.1, 1.5)
  405. axes[2].legend(loc='upper right', fontsize=9)
  407. fig.tight_layout()
  408. save_fig(fig, 'Fig1_pulse_profiles')
  411. # =============================================================================
  412. # FIGURE 2: GLUTAMATE-ONLY vs GABA-ONLY vs COMBINED (Excitatory)
  413. # =============================================================================
  414. def figure2_three_conditions(beta_nmda_val=0.03):
  415. """
  416. Separate panels showing excitatory neuron under three conditions:
  417. (i) Glutamate only, (ii) GABA only, (iii) Combined.
  418. """
  419. print(" Figure 2: Three stimulation conditions (excitatory neuron) ...")
  421. res_glu = run_neuron_simulation(beta_nmda_val, 'excitatory', enable_glutamate=True, enable_gaba=False)
  422. res_gaba = run_neuron_simulation(beta_nmda_val, 'excitatory', enable_glutamate=False, enable_gaba=True)
  423. res_both = run_neuron_simulation(beta_nmda_val, 'excitatory', enable_glutamate=True, enable_gaba=True)
  425. t0, t1 = 7000, 7500
  426. i0, i1 = int(t0/deltaT), int(t1/deltaT)
  428. conditions = [
  429. (res_glu, 'Glutamate Only (AMPA + NMDA)', '#2E86C1'),
  430. (res_gaba, 'GABA Only', '#E74C3C'),
  431. (res_both, 'Combined (Glutamate + GABA)', '#2C3E50'),
  432. ]
  434. fig, axes = plt.subplots(6, 1, figsize=(10, 14), sharex=True)
  435. labels = 'ABCDEF'
  437. for j, (res, title, color) in enumerate(conditions):
  438. ax_v = axes[2*j]
  439. ax_c = axes[2*j + 1]
  440. tt = res['t'][i0:i1]
  442. # Membrane potential
  443. ax_v.plot(tt, res['v'][i0:i1], color=color, linewidth=1.2)
  444. ax_v.axhline(0, color='gray', ls='--', lw=0.7, alpha=0.5)
  445. n_spk, rate = count_spikes(res['v'], res['t'], t_start=t0)
  446. ax_v.set_ylabel('V (mV)')
  447. ax_v.set_title(f'({labels[2*j]}) {title} | Firing rate: {rate:.1f} Hz',
  448. loc='left', fontweight='bold', fontsize=10)
  449. ax_v.set_ylim(-95, 55)
  451. # Synaptic currents
  452. ax_c.plot(tt, res['I_ampa'][i0:i1], color='#27AE60', lw=1.2, label='$I_{AMPA}$')
  453. ax_c.plot(tt, res['I_nmda'][i0:i1], color='#8E44AD', lw=1.2, label='$I_{NMDA}$')
  454. ax_c.plot(tt, res['I_gaba'][i0:i1], color='#F39C12', lw=1.2, label='$I_{GABA}$')
  455. ax_c.set_ylabel('Current (nA)')
  456. ax_c.set_title(f'({labels[2*j+1]}) Synaptic Currents | {title}',
  457. loc='left', fontsize=10)
  458. ax_c.legend(loc='upper right', ncol=3, fontsize=8)
  460. axes[-1].set_xlabel('Time (ms)')
  461. fig.suptitle(f'Excitatory Neuron Responses Under Different Synaptic Conditions '
  462. r'($\beta_{NMDA}$' + f' = {beta_nmda_val})', fontsize=13, fontweight='bold', y=1.01)
  463. fig.tight_layout()
  464. save_fig(fig, 'Fig2_three_conditions_excitatory')
  465. return res_glu, res_gaba, res_both
  468. # =============================================================================
  469. # FIGURE 3: EXCITATORY vs INHIBITORY NEURON COMPARISON
  470. # =============================================================================
  471. def figure3_exc_vs_inh(beta_nmda_val=0.03):
  472. """
  473. Side-by-side comparison of excitatory and inhibitory neurons
  474. under combined GABA + Glutamate input.
  475. """
  476. print(" Figure 3: Excitatory vs Inhibitory neuron comparison ...")
  478. res_exc = run_neuron_simulation(beta_nmda_val, 'excitatory', True, True)
  479. res_inh = run_neuron_simulation(beta_nmda_val, 'inhibitory', True, True)
  481. t0, t1 = 7000, 7500
  482. i0, i1 = int(t0/deltaT), int(t1/deltaT)
  484. fig, axes = plt.subplots(4, 2, figsize=(12, 10), sharex='col')
  486. for col, (res, ntype, clr) in enumerate([
  487. (res_exc, 'Excitatory (Pyramidal)', '#2E86C1'),
  488. (res_inh, 'Inhibitory (Fast-Spiking)', '#E74C3C')]):
  490. tt = res['t'][i0:i1]
  491. _, rate = count_spikes(res['v'], res['t'], t_start=t0)
  493. axes[0, col].plot(tt, res['v'][i0:i1], color=clr, lw=1.2)
  494. axes[0, col].axhline(0, color='gray', ls='--', lw=0.6, alpha=0.5)
  495. axes[0, col].set_ylabel('V (mV)')
  496. axes[0, col].set_title(f'{ntype}\nFiring rate: {rate:.1f} Hz', fontweight='bold')
  498. axes[1, col].plot(tt, res['I_ampa'][i0:i1], '#27AE60', lw=1, label='AMPA')
  499. axes[1, col].plot(tt, res['I_nmda'][i0:i1], '#8E44AD', lw=1, label='NMDA')
  500. axes[1, col].plot(tt, res['I_gaba'][i0:i1], '#F39C12', lw=1, label='GABA')
  501. axes[1, col].set_ylabel('I$_{syn}$ (nA)')
  502. axes[1, col].legend(fontsize=7, ncol=3)
  504. axes[2, col].plot(tt, res['Ca'][i0:i1], color='#16A085', lw=1.2)
  505. axes[2, col].set_ylabel('[Ca$^{2+}$]$_i$ (a.u.)')
  507. axes[3, col].plot(tt, res['CaMKII_p'][i0:i1], color='black', lw=1.2)
  508. axes[3, col].set_ylabel('CaMKII$_p$')
  509. axes[3, col].set_xlabel('Time (ms)')
  511. fig.suptitle(f'Excitatory vs Inhibitory Neuron | Combined Input '
  512. r'($\beta_{NMDA}$' + f' = {beta_nmda_val})',
  513. fontsize=13, fontweight='bold')
  514. fig.tight_layout()
  515. save_fig(fig, 'Fig3_exc_vs_inh')
  518. # =============================================================================
  519. # FIGURE 4: NMDA beta PARAMETER SWEEP (ENHANCED)
  520. # =============================================================================
  521. def figure4_nmda_parameter_sweep():
  522. """
  523. Sweep beta_NMDA over a wider range and show:
  524. (A) Firing rate vs beta_NMDA for both neuron types
  525. (B) CaMKII phosphorylation vs beta_NMDA
  526. (C) E/I current ratio vs beta_NMDA
  527. (D-F) Representative voltage traces at low, mid, high beta_NMDA
  528. """
  529. print(" Figure 4: NMDA parameter sweep (enhanced) ...")
  531. # Wider range with finer resolution
  532. beta_vals = np.arange(0.005, 0.121, 0.005)
  533. n_beta = len(beta_vals)
  535. rates_exc = np.zeros(n_beta)
  536. rates_inh = np.zeros(n_beta)
  537. camkii_exc = np.zeros(n_beta)
  538. camkii_inh = np.zeros(n_beta)
  539. ei_ratio_exc = np.zeros(n_beta)
  540. ei_ratio_inh = np.zeros(n_beta)
  542. # Store full results at selected beta values for traces
  543. beta_trace_vals = [0.01, 0.04, 0.10]
  544. trace_results = {bv: {} for bv in beta_trace_vals}
  546. idx_ss = int(5000 / deltaT)
  548. for k, bv in enumerate(beta_vals):
  549. print(f" beta_NMDA = {bv:.3f} ({k+1}/{n_beta})")
  550. re = run_neuron_simulation(bv, 'excitatory', True, True)
  551. ri = run_neuron_simulation(bv, 'inhibitory', True, True)
  553. _, rates_exc[k] = count_spikes(re['v'], re['t'])
  554. _, rates_inh[k] = count_spikes(ri['v'], ri['t'])
  555. camkii_exc[k] = np.mean(re['CaMKII_p'][idx_ss:])
  556. camkii_inh[k] = np.mean(ri['CaMKII_p'][idx_ss:])
  557. ei_ratio_exc[k] = compute_ei_ratio(re)
  558. ei_ratio_inh[k] = compute_ei_ratio(ri)
  560. # Save traces at selected values
  561. for btv in beta_trace_vals:
  562. if abs(bv - btv) < 1e-6:
  563. trace_results[btv] = {'exc': re, 'inh': ri}
  565. # -- Create figure: 2 rows x 3 columns --
  566. fig = plt.figure(figsize=(14, 9))
  567. gs = gridspec.GridSpec(2, 3, hspace=0.45, wspace=0.35)
  569. clr_exc = '#2E86C1'
  570. clr_inh = '#E74C3C'
  572. # (A) Firing rate
  573. ax_a = fig.add_subplot(gs[0, 0])
  574. ax_a.plot(beta_vals, rates_exc, 'o-', color=clr_exc, lw=2, ms=4, label='Excitatory')
  575. ax_a.plot(beta_vals, rates_inh, 's--', color=clr_inh, lw=2, ms=4, label='Inhibitory')
  576. ax_a.set_xlabel(r'$\beta_{NMDA}$')
  577. ax_a.set_ylabel('Firing Rate (Hz)')
  578. ax_a.set_title(r'(A) Firing Rate vs $\beta_{NMDA}$', loc='left', fontweight='bold')
  579. ax_a.legend(fontsize=8)
  580. # Mark trace locations
  581. for btv in beta_trace_vals:
  582. ax_a.axvline(btv, color='gray', ls=':', lw=0.8, alpha=0.5)
  584. # (B) CaMKII
  585. ax_b = fig.add_subplot(gs[0, 1])
  586. ax_b.semilogy(beta_vals, camkii_exc, 'o-', color=clr_exc, lw=2, ms=4, label='Excitatory')
  587. ax_b.semilogy(beta_vals, camkii_inh, 's--', color=clr_inh, lw=2, ms=4, label='Inhibitory')
  588. ax_b.set_xlabel(r'$\beta_{NMDA}$')
  589. ax_b.set_ylabel('Mean CaMKII$_p$ (log scale)')
  590. ax_b.set_title('(B) CaMKII Phosphorylation', loc='left', fontweight='bold')
  591. ax_b.legend(fontsize=8)
  593. # (C) E/I ratio
  594. ax_c = fig.add_subplot(gs[0, 2])
  595. # Clip inf values for plotting
  596. ei_exc_plot = np.clip(ei_ratio_exc, 0, 50)
  597. ei_inh_plot = np.clip(ei_ratio_inh, 0, 50)
  598. ax_c.plot(beta_vals, ei_exc_plot, 'o-', color=clr_exc, lw=2, ms=4, label='Excitatory')
  599. ax_c.plot(beta_vals, ei_inh_plot, 's--', color=clr_inh, lw=2, ms=4, label='Inhibitory')
  600. ax_c.axhline(1.0, color='gray', ls='--', lw=1, alpha=0.6, label='E/I = 1')
  601. ax_c.set_xlabel(r'$\beta_{NMDA}$')
  602. ax_c.set_ylabel('E/I Current Ratio')
  603. ax_c.set_title('(C) Excitation/Inhibition Balance', loc='left', fontweight='bold')
  604. ax_c.legend(fontsize=8)
  606. # (D-F) Representative traces
  607. trace_labels = [r'(D) Low $\beta_{NMDA}$', r'(E) Mid $\beta_{NMDA}$', r'(F) High $\beta_{NMDA}$']
  608. t0_tr, t1_tr = 7000, 7400 # 400 ms window
  609. i0_tr = int(t0_tr / deltaT)
  610. i1_tr = int(t1_tr / deltaT)
  612. for j, btv in enumerate(beta_trace_vals):
  613. ax_tr = fig.add_subplot(gs[1, j])
  614. if btv in trace_results and trace_results[btv]:
  615. tt = trace_results[btv]['exc']['t'][i0_tr:i1_tr]
  616. ax_tr.plot(tt, trace_results[btv]['exc']['v'][i0_tr:i1_tr],
  617. color=clr_exc, lw=1.2, label='Exc')
  618. ax_tr.plot(tt, trace_results[btv]['inh']['v'][i0_tr:i1_tr],
  619. color=clr_inh, lw=1.2, alpha=0.8, label='Inh')
  620. ax_tr.axhline(0, color='gray', ls='--', lw=0.5, alpha=0.4)
  622. ax_tr.set_ylabel('V (mV)')
  623. ax_tr.set_xlabel('Time (ms)')
  624. ax_tr.set_title(f'{trace_labels[j]} = {btv}', loc='left', fontweight='bold')
  625. ax_tr.set_ylim(-95, 55)
  626. if j == 0:
  627. ax_tr.legend(fontsize=8)
  629. fig.suptitle(r'Effect of $\beta_{NMDA}$ on Excitatory and Inhibitory Neuron Responses',
  630. fontsize=14, fontweight='bold')
  631. fig.tight_layout(rect=[0, 0, 1, 0.96])
  632. save_fig(fig, 'Fig4_nmda_sweep_enhanced')
  635. # =============================================================================
  636. # FIGURE 5: CONDUCTANCE TUNING -- DOSE-RESPONSE + DUAL HEATMAPS
  637. # =============================================================================
  638. def figure5_conductance_tuning():
  639. """
  640. (A) GABA dose-response (glutamate fixed at 1x) for both cell types
  641. (B) Glutamate dose-response (GABA fixed at 1x) for both cell types
  642. (C) E-I balance heatmap: excitatory neuron
  643. (D) E-I balance heatmap: inhibitory neuron
  644. """
  645. print(" Figure 5: Conductance tuning (enhanced) ...")
  647. beta_val = 0.03
  648. clr_exc = '#2E86C1'
  649. clr_inh = '#E74C3C'
  651. # -- (A) GABA dose-response --
  652. gaba_scales_1d = np.array([0, 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 4.0])
  653. rates_gaba_exc = np.zeros(len(gaba_scales_1d))
  654. rates_gaba_inh = np.zeros(len(gaba_scales_1d))
  656. print(" GABA dose-response ...")
  657. for k, gs in enumerate(gaba_scales_1d):
  658. en_gab = gs > 0
  659. re = run_neuron_simulation(beta_val, 'excitatory', True, en_gab,
  660. g_gaba_scale=max(gs, 1e-6))
  661. ri = run_neuron_simulation(beta_val, 'inhibitory', True, en_gab,
  662. g_gaba_scale=max(gs, 1e-6))
  663. _, rates_gaba_exc[k] = count_spikes(re['v'], re['t'])
  664. _, rates_gaba_inh[k] = count_spikes(ri['v'], ri['t'])
  666. # -- (B) Glutamate dose-response --
  667. glu_scales_1d = np.array([0, 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 4.0])
  668. rates_glu_exc = np.zeros(len(glu_scales_1d))
  669. rates_glu_inh = np.zeros(len(glu_scales_1d))
  671. print(" Glutamate dose-response ...")
  672. for k, gs in enumerate(glu_scales_1d):
  673. en_glu = gs > 0
  674. re = run_neuron_simulation(beta_val, 'excitatory', en_glu, True,
  675. g_ampa_scale=max(gs, 1e-6),
  676. g_nmda_scale=max(gs, 1e-6))
  677. ri = run_neuron_simulation(beta_val, 'inhibitory', en_glu, True,
  678. g_ampa_scale=max(gs, 1e-6),
  679. g_nmda_scale=max(gs, 1e-6))
  680. _, rates_glu_exc[k] = count_spikes(re['v'], re['t'])
  681. _, rates_glu_inh[k] = count_spikes(ri['v'], ri['t'])
  683. # -- (C, D) 2D heatmaps for both neuron types --
  684. gaba_scales_2d = np.array([0, 0.5, 1.0, 1.5, 2.0, 3.0])
  685. glu_scales_2d = np.array([0, 0.5, 1.0, 1.5, 2.0, 3.0])
  687. rate_map_exc = np.zeros((len(gaba_scales_2d), len(glu_scales_2d)))
  688. rate_map_inh = np.zeros((len(gaba_scales_2d), len(glu_scales_2d)))
  690. print(" 2D heatmaps ...")
  691. for gi, gs_gaba in enumerate(gaba_scales_2d):
  692. for gj, gs_glu in enumerate(glu_scales_2d):
  693. en_glu = gs_glu > 0
  694. en_gab = gs_gaba > 0
  695. re = run_neuron_simulation(
  696. beta_val, 'excitatory', en_glu, en_gab,
  697. g_ampa_scale=max(gs_glu, 1e-6),
  698. g_gaba_scale=max(gs_gaba, 1e-6),
  699. g_nmda_scale=max(gs_glu, 1e-6))
  700. ri = run_neuron_simulation(
  701. beta_val, 'inhibitory', en_glu, en_gab,
  702. g_ampa_scale=max(gs_glu, 1e-6),
  703. g_gaba_scale=max(gs_gaba, 1e-6),
  704. g_nmda_scale=max(gs_glu, 1e-6))
  705. _, rate_map_exc[gi, gj] = count_spikes(re['v'], re['t'])
  706. _, rate_map_inh[gi, gj] = count_spikes(ri['v'], ri['t'])
  708. # -- Plot --
  709. fig = plt.figure(figsize=(14, 10))
  710. gs_fig = gridspec.GridSpec(2, 2, hspace=0.38, wspace=0.3)
  712. # (A) GABA dose-response
  713. ax_a = fig.add_subplot(gs_fig[0, 0])
  714. ax_a.plot(gaba_scales_1d, rates_gaba_exc, 'o-', color=clr_exc, lw=2, ms=5, label='Excitatory')
  715. ax_a.plot(gaba_scales_1d, rates_gaba_inh, 's--', color=clr_inh, lw=2, ms=5, label='Inhibitory')
  716. ax_a.axvline(1.0, color='gray', ls=':', lw=0.8, alpha=0.6, label='Baseline')
  717. ax_a.set_xlabel('GABAergic Conductance Scale ($g_{GABA}$ / $g_{GABA,0}$)')
  718. ax_a.set_ylabel('Firing Rate (Hz)')
  719. ax_a.set_title('(A) GABA Dose-Response\n(Glutamate fixed at 1$\\times$)',
  720. loc='left', fontweight='bold')
  721. ax_a.legend(fontsize=8)
  723. # (B) Glutamate dose-response
  724. ax_b = fig.add_subplot(gs_fig[0, 1])
  725. ax_b.plot(glu_scales_1d, rates_glu_exc, 'o-', color=clr_exc, lw=2, ms=5, label='Excitatory')
  726. ax_b.plot(glu_scales_1d, rates_glu_inh, 's--', color=clr_inh, lw=2, ms=5, label='Inhibitory')
  727. ax_b.axvline(1.0, color='gray', ls=':', lw=0.8, alpha=0.6, label='Baseline')
  728. ax_b.set_xlabel('Glutamatergic Conductance Scale ($g_{AMPA/NMDA}$ / $g_0$)')
  729. ax_b.set_ylabel('Firing Rate (Hz)')
  730. ax_b.set_title('(B) Glutamate Dose-Response\n(GABA fixed at 1$\\times$)',
  731. loc='left', fontweight='bold')
  732. ax_b.legend(fontsize=8)
  734. # (C) Heatmap -- excitatory
  735. ax_c = fig.add_subplot(gs_fig[1, 0])
  736. vmax = max(rate_map_exc.max(), rate_map_inh.max())
  737. im_c = ax_c.imshow(rate_map_exc, origin='lower', aspect='auto',
  738. extent=[glu_scales_2d[0]-0.25, glu_scales_2d[-1]+0.25,
  739. gaba_scales_2d[0]-0.25, gaba_scales_2d[-1]+0.25],
  740. cmap='RdYlBu_r', interpolation='bilinear',
  741. vmin=0, vmax=vmax)
  742. for gi, gs_gaba in enumerate(gaba_scales_2d):
  743. for gj, gs_glu in enumerate(glu_scales_2d):
  744. val = rate_map_exc[gi, gj]
  745. ax_c.text(gs_glu, gs_gaba, f'{val:.0f}',
  746. ha='center', va='center', fontsize=7,
  747. color='white' if val > vmax*0.5 else 'black')
  748. fig.colorbar(im_c, ax=ax_c, label='Firing Rate (Hz)', shrink=0.85)
  749. ax_c.set_xlabel('Glutamatergic Scale')
  750. ax_c.set_ylabel('GABAergic Scale')
  751. ax_c.set_title('(C) E-I Balance | Excitatory Neuron', loc='left', fontweight='bold')
  752. ax_c.set_xticks(glu_scales_2d)
  753. ax_c.set_yticks(gaba_scales_2d)
  755. # (D) Heatmap -- inhibitory
  756. ax_d = fig.add_subplot(gs_fig[1, 1])
  757. im_d = ax_d.imshow(rate_map_inh, origin='lower', aspect='auto',
  758. extent=[glu_scales_2d[0]-0.25, glu_scales_2d[-1]+0.25,
  759. gaba_scales_2d[0]-0.25, gaba_scales_2d[-1]+0.25],
  760. cmap='RdYlBu_r', interpolation='bilinear',
  761. vmin=0, vmax=vmax)
  762. for gi, gs_gaba in enumerate(gaba_scales_2d):
  763. for gj, gs_glu in enumerate(glu_scales_2d):
  764. val = rate_map_inh[gi, gj]
  765. ax_d.text(gs_glu, gs_gaba, f'{val:.0f}',
  766. ha='center', va='center', fontsize=7,
  767. color='white' if val > vmax*0.5 else 'black')
  768. fig.colorbar(im_d, ax=ax_d, label='Firing Rate (Hz)', shrink=0.85)
  769. ax_d.set_xlabel('Glutamatergic Scale')
  770. ax_d.set_ylabel('GABAergic Scale')
  771. ax_d.set_title('(D) E-I Balance | Inhibitory Neuron', loc='left', fontweight='bold')
  772. ax_d.set_xticks(glu_scales_2d)
  773. ax_d.set_yticks(gaba_scales_2d)
  775. fig.suptitle(r'Tuning GABAergic and Glutamatergic Stimulation ($\beta_{NMDA}$' + f' = {beta_val})',
  776. fontsize=14, fontweight='bold')
  777. fig.tight_layout(rect=[0, 0, 1, 0.96])
  778. save_fig(fig, 'Fig5_conductance_tuning_enhanced')
  781. # =============================================================================
  782. # FIGURE 6: THREE CONDITIONS FOR INHIBITORY NEURON
  783. # =============================================================================
  784. def figure6_three_conditions_inhibitory(beta_nmda_val=0.03):
  785. """Same three-condition comparison but for the inhibitory neuron."""
  786. print(" Figure 6: Three conditions (inhibitory neuron) ...")
  788. res_glu = run_neuron_simulation(beta_nmda_val, 'inhibitory', True, False)
  789. res_gaba = run_neuron_simulation(beta_nmda_val, 'inhibitory', False, True)
  790. res_both = run_neuron_simulation(beta_nmda_val, 'inhibitory', True, True)
  792. t0, t1 = 7000, 7500
  793. i0, i1 = int(t0/deltaT), int(t1/deltaT)
  795. conditions = [
  796. (res_glu, 'Glutamate Only', '#2E86C1'),
  797. (res_gaba, 'GABA Only', '#E74C3C'),
  798. (res_both, 'Combined', '#2C3E50'),
  799. ]
  801. fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)
  803. for j, (res, title, color) in enumerate(conditions):
  804. tt = res['t'][i0:i1]
  805. ax = axes[j]
  806. ax.plot(tt, res['v'][i0:i1], color=color, lw=1.2)
  807. ax.axhline(0, color='gray', ls='--', lw=0.6, alpha=0.5)
  808. _, rate = count_spikes(res['v'], res['t'], t_start=t0)
  809. ax.set_ylabel('V (mV)')
  810. ax.set_title(f'({chr(65+j)}) Inhibitory Neuron | {title} | '
  811. f'Rate: {rate:.1f} Hz', loc='left', fontweight='bold', fontsize=10)
  812. ax.set_ylim(-95, 55)
  814. axes[-1].set_xlabel('Time (ms)')
  815. fig.suptitle(f'Inhibitory (Fast-Spiking) Neuron Responses '
  816. r'($\beta_{NMDA}$' + f' = {beta_nmda_val})',
  817. fontsize=13, fontweight='bold')
  818. fig.tight_layout()
  819. save_fig(fig, 'Fig6_three_conditions_inhibitory')
  822. # =============================================================================
  823. # FIGURE 7: beta_NMDA x NEURON TYPE x INPUT CONDITION INTERACTION
  824. # =============================================================================
  825. def figure7_interaction_summary():
  826. """
  827. Comprehensive summary: how beta_NMDA modulates firing rate under each
  828. input condition (Glu only, GABA only, Combined) for both neuron types.
  829. """
  830. print(" Figure 7: Interaction summary ...")
  832. beta_vals = np.arange(0.005, 0.101, 0.01)
  833. n_beta = len(beta_vals)
  835. conditions = [
  836. ('Glutamate Only', True, False),
  837. ('GABA Only', False, True),
  838. ('Combined', True, True),
  839. ]
  840. neuron_types = [
  841. ('Excitatory', 'excitatory', '#2E86C1'),
  842. ('Inhibitory', 'inhibitory', '#E74C3C'),
  843. ]
  845. # results[condition_name][neuron_label] = array of firing rates
  846. results = {}
  847. camkii_results = {}
  848. idx_ss = int(5000 / deltaT)
  850. for cond_name, en_glu, en_gab in conditions:
  851. results[cond_name] = {}
  852. camkii_results[cond_name] = {}
  853. for n_label, n_type, _ in neuron_types:
  854. rates = np.zeros(n_beta)
  855. camkii = np.zeros(n_beta)
  856. for k, bv in enumerate(beta_vals):
  857. res = run_neuron_simulation(bv, n_type, en_glu, en_gab)
  858. _, rates[k] = count_spikes(res['v'], res['t'])
  859. camkii[k] = np.mean(res['CaMKII_p'][idx_ss:])
  860. results[cond_name][n_label] = rates
  861. camkii_results[cond_name][n_label] = camkii
  863. # -- Plot: 2 rows x 3 columns --
  864. fig, axes = plt.subplots(2, 3, figsize=(14, 8), sharey='row')
  865. panel_labels = iter('ABCDEF')
  867. for j, (cond_name, _, _) in enumerate(conditions):
  868. ax_rate = axes[0, j]
  869. ax_cam = axes[1, j]
  870. lbl = next(panel_labels)
  871. lbl2 = next(panel_labels)
  873. for n_label, n_type, clr in neuron_types:
  874. ax_rate.plot(beta_vals, results[cond_name][n_label],
  875. 'o-' if n_label == 'Excitatory' else 's--',
  876. color=clr, lw=2, ms=4, label=n_label)
  877. # Guard against zero/negative values for log scale
  878. # Clamp floor at 1e-10 (values below are numerical zero)
  879. cam_data = camkii_results[cond_name][n_label]
  880. cam_data_safe = np.clip(cam_data, 1e-10, None)
  881. ax_cam.semilogy(beta_vals, cam_data_safe,
  882. 'o-' if n_label == 'Excitatory' else 's--',
  883. color=clr, lw=2, ms=4, label=n_label)
  885. ax_rate.set_title(f'({lbl}) {cond_name}', loc='left', fontweight='bold')
  886. ax_rate.set_xlabel(r'$\beta_{NMDA}$')
  887. if j == 0:
  888. ax_rate.set_ylabel('Firing Rate (Hz)')
  889. ax_rate.legend(fontsize=8)
  891. ax_cam.set_title(f'({lbl2}) {cond_name}', loc='left', fontweight='bold')
  892. ax_cam.set_xlabel(r'$\beta_{NMDA}$')
  893. if j == 0:
  894. ax_cam.set_ylabel('Mean CaMKII$_p$ (log)')
  895. ax_cam.set_ylim(bottom=1e-10) # floor: values below are numerical zero
  897. fig.suptitle(r'$\beta_{NMDA}$ $\times$ Neuron Type $\times$ Input Condition Interaction',
  898. fontsize=14, fontweight='bold')
  899. fig.tight_layout(rect=[0, 0, 1, 0.96])
  900. save_fig(fig, 'Fig7_interaction_summary')
  903. # =============================================================================
  904. # FIGURE 8: AMPA RECEPTOR CLOSING RATE (β_AMPA) SENSITIVITY ANALYSIS
  905. # =============================================================================
  906. def figure8_ampa_closing_rate_sensitivity():
  907. """
  908. AMPA receptor closing rate (β_AMPA) sensitivity analysis.
  910. β_AMPA (AMPA unbinding rate) controls AMPA EPSC decay kinetics.
  911. Physiological range: ~0.2–2.0 ms⁻¹ (decay τ ≈ 0.5–5 ms).
  912. Default value: 0.67 ms⁻¹ (τ ≈ 1.5 ms).
  914. Panels:
  915. (A) AMPA EPSC waveforms at selected β_AMPA values
  916. (B) Firing rate vs β_AMPA — excitatory & inhibitory neurons
  917. (C) Mean CaMKII phosphorylation vs β_AMPA — both neuron types
  918. (D) E/I current ratio vs β_AMPA — both neuron types
  919. (E) 2D heatmap: β_AMPA × β_NMDA interaction (firing rate, excitatory)
  920. (F) 2D heatmap: β_AMPA × β_NMDA interaction (CaMKII, excitatory)
  921. """
  922. print(" Figure 8: β_AMPA (AMPA closing rate) sensitivity ...")
  924. # --- β_AMPA values (finer resolution around 0.3–0.7 transition) ---
  925. beta_ampa_1d = np.array([0.1, 0.2, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6,
  926. 0.67, 0.8, 1.0, 1.5, 2.0, 3.0])
  927. beta_nmda_fixed = 0.03 # fixed for 1D sweeps
  929. clr_exc = '#2E86C1'
  930. clr_inh = '#E74C3C'
  931. idx_ss = int(5000 / deltaT) # steady-state index
  933. # ── 1D sweeps ──────────────────────────────────────────────────────
  934. fr_exc = np.zeros(len(beta_ampa_1d))
  935. fr_inh = np.zeros(len(beta_ampa_1d))
  936. camk_exc = np.zeros(len(beta_ampa_1d))
  937. camk_inh = np.zeros(len(beta_ampa_1d))
  938. ei_exc = np.zeros(len(beta_ampa_1d))
  939. ei_inh = np.zeros(len(beta_ampa_1d))
  941. print(" 1D β_AMPA sweep ...")
  942. for k, ba in enumerate(beta_ampa_1d):
  943. re = run_neuron_simulation(beta_nmda_fixed, 'excitatory',
  944. beta_ampa_val=ba)
  945. ri = run_neuron_simulation(beta_nmda_fixed, 'inhibitory',
  946. beta_ampa_val=ba)
  947. _, fr_exc[k] = count_spikes(re['v'], re['t'])
  948. _, fr_inh[k] = count_spikes(ri['v'], ri['t'])
  949. camk_exc[k] = np.mean(re['CaMKII_p'][idx_ss:])
  950. camk_inh[k] = np.mean(ri['CaMKII_p'][idx_ss:])
  951. ei_exc[k] = compute_ei_ratio(re)
  952. ei_inh[k] = compute_ei_ratio(ri)
  953. if k % 3 == 0:
  954. print(f" β_AMPA = {ba:.2f}: FR(exc) = {fr_exc[k]:.1f} Hz, "
  955. f"CaMKII = {camk_exc[k]:.2e}")
  957. # ── 2D heatmap: β_AMPA × β_NMDA ──────────────────────────────────
  958. beta_ampa_2d = np.array([0.2, 0.4, 0.67, 1.0, 1.5, 2.0])
  959. beta_nmda_2d = np.array([0.005, 0.01, 0.02, 0.03, 0.05, 0.08, 0.12])
  961. rate_map = np.zeros((len(beta_ampa_2d), len(beta_nmda_2d)))
  962. camk_map = np.zeros((len(beta_ampa_2d), len(beta_nmda_2d)))
  964. print(" 2D β_AMPA × β_NMDA heatmap ...")
  965. for ai, ba in enumerate(beta_ampa_2d):
  966. for nj, bn in enumerate(beta_nmda_2d):
  967. re = run_neuron_simulation(bn, 'excitatory', beta_ampa_val=ba)
  968. _, rate_map[ai, nj] = count_spikes(re['v'], re['t'])
  969. camk_map[ai, nj] = np.mean(re['CaMKII_p'][idx_ss:])
  971. # ── AMPA EPSC waveforms (analytical, for illustration) ────────────
  972. # Single-pulse AMPA gating: m_AMPA(t) = (α/(α+β)) * exp(-β*t) for t > pulse
  973. t_epsc = np.linspace(0, 30, 600) # 30 ms window
  974. pulse_dur = 5.0 # ms
  975. epsc_betas = [0.2, 0.67, 1.5, 3.0]
  976. epsc_colors = ['#2C3E50', '#2E86C1', '#E67E22', '#E74C3C']
  978. # ── Plot ───────────────────────────────────────────────────────────
  979. fig = plt.figure(figsize=(18, 11))
  980. gs_fig = gridspec.GridSpec(2, 3, hspace=0.38, wspace=0.32)
  982. # (A) AMPA EPSC waveforms
  983. ax_a = fig.add_subplot(gs_fig[0, 0])
  984. for ba_val, col in zip(epsc_betas, epsc_colors):
  985. # Simulate single-pulse AMPA gating
  986. dt_ep = 0.05
  987. tn_ep = int(30 / dt_ep)
  988. t_ep = np.arange(tn_ep + 1) * dt_ep
  989. m_ep = np.zeros(tn_ep + 1)
  990. for i_ep in range(tn_ep):
  991. G_ep = 1.0 if t_ep[i_ep] <= pulse_dur else 0.0
  992. m_ep[i_ep + 1] = m_ep[i_ep] + dt_ep * (
  993. alpha_ampa * G_ep * (1 - m_ep[i_ep]) - ba_val * m_ep[i_ep])
  994. # Normalise to peak for shape comparison
  995. if m_ep.max() > 0:
  996. m_norm = m_ep / m_ep.max()
  997. else:
  998. m_norm = m_ep
  999. tau_decay = 1.0 / ba_val
  1000. ax_a.plot(t_ep, m_norm, color=col, lw=2,
  1001. label=f'β = {ba_val:.2f} (τ ≈ {tau_decay:.1f} ms)')
  1002. ax_a.set_xlabel('Time (ms)')
  1003. ax_a.set_ylabel('Normalised AMPA conductance')
  1004. ax_a.set_title('(A) AMPA EPSC Waveforms\n(single pulse, normalised)',
  1005. loc='left', fontweight='bold')
  1006. ax_a.legend(fontsize=8, title='$β_{AMPA}$', title_fontsize=9)
  1008. # (B) Firing rate vs β_AMPA
  1009. ax_b = fig.add_subplot(gs_fig[0, 1])
  1010. ax_b.plot(beta_ampa_1d, fr_exc, 'o-', color=clr_exc, lw=2, ms=6,
  1011. label='Excitatory')
  1012. ax_b.plot(beta_ampa_1d, fr_inh, 's--', color=clr_inh, lw=2, ms=6,
  1013. label='Inhibitory')
  1014. ax_b.axvline(0.67, color='gray', ls=':', lw=1, alpha=0.6, label='Default (0.67)')
  1015. ax_b.set_xlabel(r'$\beta_{AMPA}$ (ms$^{-1}$)')
  1016. ax_b.set_ylabel('Firing Rate (Hz)')
  1017. ax_b.set_title(f'(B) Firing Rate vs $β_{{AMPA}}$\n'
  1018. f'($β_{{NMDA}}$ = {beta_nmda_fixed})',
  1019. loc='left', fontweight='bold')
  1020. ax_b.legend(fontsize=8)
  1022. # (C) CaMKII vs β_AMPA
  1023. ax_c = fig.add_subplot(gs_fig[0, 2])
  1024. ax_c.semilogy(beta_ampa_1d, camk_exc, 'o-', color=clr_exc, lw=2, ms=6,
  1025. label='Excitatory')
  1026. ax_c.semilogy(beta_ampa_1d, camk_inh, 's--', color=clr_inh, lw=2, ms=6,
  1027. label='Inhibitory')
  1028. ax_c.axvline(0.67, color='gray', ls=':', lw=1, alpha=0.6, label='Default')
  1029. ax_c.set_xlabel(r'$\beta_{AMPA}$ (ms$^{-1}$)')
  1030. ax_c.set_ylabel('Mean CaMKII$_p$ (log scale)')
  1031. ax_c.set_title(f'(C) CaMKII Phosphorylation vs $β_{{AMPA}}$',
  1032. loc='left', fontweight='bold')
  1033. ax_c.legend(fontsize=8)
  1035. # (D) E/I ratio vs β_AMPA
  1036. ax_d = fig.add_subplot(gs_fig[1, 0])
  1037. ax_d.plot(beta_ampa_1d, ei_exc, 'o-', color=clr_exc, lw=2, ms=6,
  1038. label='Excitatory')
  1039. ax_d.plot(beta_ampa_1d, ei_inh, 's--', color=clr_inh, lw=2, ms=6,
  1040. label='Inhibitory')
  1041. ax_d.axvline(0.67, color='gray', ls=':', lw=1, alpha=0.6, label='Default')
  1042. ax_d.axhline(1.0, color='black', ls=':', lw=0.8, alpha=0.4)
  1043. ax_d.set_xlabel(r'$\beta_{AMPA}$ (ms$^{-1}$)')
  1044. ax_d.set_ylabel('E/I Current Ratio')
  1045. ax_d.set_title('(D) Excitation–Inhibition Balance vs $β_{AMPA}$',
  1046. loc='left', fontweight='bold')
  1047. ax_d.legend(fontsize=8)
  1049. # (E) 2D heatmap: β_AMPA × β_NMDA → Firing Rate
  1050. ax_e = fig.add_subplot(gs_fig[1, 1])
  1051. im_e = ax_e.imshow(rate_map, origin='lower', aspect='auto',
  1052. extent=[beta_nmda_2d[0], beta_nmda_2d[-1],
  1053. beta_ampa_2d[0], beta_ampa_2d[-1]],
  1054. cmap='viridis', interpolation='bilinear')
  1055. for ai, ba in enumerate(beta_ampa_2d):
  1056. for nj, bn in enumerate(beta_nmda_2d):
  1057. val = rate_map[ai, nj]
  1058. ax_e.text(bn, ba, f'{val:.0f}', ha='center', va='center',
  1059. fontsize=6, color='white' if val > rate_map.max()*0.5 else 'black')
  1060. fig.colorbar(im_e, ax=ax_e, label='Firing Rate (Hz)', shrink=0.85)
  1061. ax_e.set_xlabel(r'$\beta_{NMDA}$')
  1062. ax_e.set_ylabel(r'$\beta_{AMPA}$ (ms$^{-1}$)')
  1063. ax_e.set_title('(E) $β_{AMPA}$ × $β_{NMDA}$ → FR\n(Excitatory)',
  1064. loc='left', fontweight='bold')
  1065. # Mark default
  1066. ax_e.plot(0.03, 0.67, 'r*', markersize=12, zorder=5)
  1068. # (F) 2D heatmap: β_AMPA × β_NMDA → CaMKII
  1069. ax_f = fig.add_subplot(gs_fig[1, 2])
  1070. # Use log scale for CaMKII
  1071. camk_log = np.log10(camk_map + 1e-30)
  1072. im_f = ax_f.imshow(camk_log, origin='lower', aspect='auto',
  1073. extent=[beta_nmda_2d[0], beta_nmda_2d[-1],
  1074. beta_ampa_2d[0], beta_ampa_2d[-1]],
  1075. cmap='magma', interpolation='bilinear')
  1076. for ai, ba in enumerate(beta_ampa_2d):
  1077. for nj, bn in enumerate(beta_nmda_2d):
  1078. val = camk_map[ai, nj]
  1079. ax_f.text(bn, ba, f'{val:.0e}', ha='center', va='center',
  1080. fontsize=5, color='white')
  1081. fig.colorbar(im_f, ax=ax_f, label='log₁₀(CaMKII$_p$)', shrink=0.85)
  1082. ax_f.set_xlabel(r'$\beta_{NMDA}$')
  1083. ax_f.set_ylabel(r'$\beta_{AMPA}$ (ms$^{-1}$)')
  1084. ax_f.set_title('(F) $β_{AMPA}$ × $β_{NMDA}$ → CaMKII\n(Excitatory)',
  1085. loc='left', fontweight='bold')
  1086. ax_f.plot(0.03, 0.67, 'r*', markersize=12, zorder=5)
  1088. fig.suptitle(
  1089. r'AMPA Receptor Closing Rate ($\beta_{AMPA}$) Sensitivity Analysis',
  1090. fontsize=14, fontweight='bold')
  1091. fig.tight_layout(rect=[0, 0, 1, 0.96])
  1092. save_fig(fig, 'Fig8_beta_AMPA_sensitivity')
  1094. # ── Print summary ──────────────────────────────────────────────────
  1095. print("\n ── β_AMPA Sensitivity Summary ──")
  1096. print(f" Range tested: {beta_ampa_1d[0]:.2f} – {beta_ampa_1d[-1]:.2f} ms⁻¹")
  1097. print(f" Default value: 0.67 ms⁻¹ (τ_decay ≈ 1.5 ms)")
  1098. fr_range = fr_exc.max() - fr_exc.min()
  1099. fr_pct = 100 * fr_range / fr_exc[beta_ampa_1d == 0.67][0] if 0.67 in beta_ampa_1d else 0
  1100. print(f" Excitatory FR range: {fr_exc.min():.1f} – {fr_exc.max():.1f} Hz "
  1101. f"(Δ = {fr_pct:.1f}%)")
  1102. print(f" Inhibitory FR range: {fr_inh.min():.1f} – {fr_inh.max():.1f} Hz")
  1103. camk_ratio = camk_exc.max() / (camk_exc.min() + 1e-30)
  1104. print(f" CaMKII range (exc): {camk_exc.min():.2e} – {camk_exc.max():.2e} "
  1105. f"({camk_ratio:.1f}× variation)")
  1108. # =============================================================================
  1109. # RUN ALL FIGURES
  1110. # =============================================================================
  1111. if __name__ == '__main__':
  1112. print("=" * 60)
  1113. print("Generating Supplementary Figures")
  1114. print("=" * 60)
  1116. figure1_pulse_profiles()
  1117. figure2_three_conditions()
  1118. figure3_exc_vs_inh()
  1119. figure4_nmda_parameter_sweep()
  1120. figure5_conductance_tuning()
  1121. figure6_three_conditions_inhibitory()
  1122. figure7_interaction_summary()
  1123. figure8_ampa_closing_rate_sensitivity()
  1125. print("\n" + "=" * 60)
  1126. print(f"All figures saved to: {save_directory}")
  1127. print("=" * 60)

reviewer_response_simulations_final.py at commit 006879e, under MIT · at the source

Overview

Authors: Mehdi Borjkhani1,2, Hadi Borjkhani3, Morteza A Sharif4, Fariba Bahrami5, Mahyar Janahmadi6
ORCID iDs: Mehdi Borjkhani, Hadi Borjkhani
  1. International Centre for Translational Eye Research (ICTER), Institute of Physical Chemistry, Polish Academy of Sciences, Warsaw, Poland
  2. Institute of Physical Chemistry, Polish Academy of Sciences, Warsaw, Poland
  3. School of Engineering – Energy and Information, HTW Berlin – University of Applied Sciences, Berlin, Germany
  4. Optics and Laser Engineering Group, Faculty of Electrical Engineering, Urmia University of Technology, Urmia, Iran
  5. CIPCE, Motor Control and Computational Neuroscience Laboratory, School of ECE, College of Engineering, University of Tehran, Tehran, Iran
  6. Neuroscience Research Center and Department of Physiology, Medical School, Shahid Beheshti University of Medical Sciences, Tehran, Iran
Institutions: Institute of Physical Chemistry (Poland); Polish Academy of Sciences (Poland); HTW Berlin - University of Applied Sciences (Germany); Urmia University of Technology (Iran); University of Tehran (Iran); Shahid Beheshti University (Iran); Shahid Beheshti University of Medical Sciences (Iran)
Journal: Frontiers in computational neuroscience, volume 20, article 1753444
Dates: received 24 November 2025; accepted 21 April 2026; published online 3 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.3389/fncom.2026.1753444 · PMID 42318008 · PMCID PMC13272316 · OpenAlex W4414141469
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.09.06.674628
Status: code verified
Categories: none (in silico) (organism), other condition (population), cellular / molecular (subfield)
Methods: Smoothing, state filtering, decompositions, Connectivity, Single-unit activity, calcium imaging
Keywords: addiction, chaos theory, computational neuroscience, information theory, neuronal dynamics, NMDA receptors, synaptic plasticity, visual processing
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: European Regional Development Fund; Fundacja na rzecz Nauki Polskiej
Citations: not cited yet (Europe PMC); 83 references in the paper

Abstract

Introduction: Neuronal firing patterns emerge from complex interactions between intrinsic membrane properties and synaptic receptor dynamics. N-methyl-D-aspartate (NMDA) receptors critically shape calcium influx and synaptic plasticity through their voltage-dependent Mg2+ block and prolonged activation kinetics, yet how their closing kinetics interact with glutamatergic drive and GABAergic modulation to control neuronal dynamics and information processing remains incompletely understood.

Methods: We developed a Hodgkin–Huxley-type computational model incorporating NMDA, AMPA, and GABA receptor kinetics to investigate how the NMDA receptor closing rate βNMDA and glutamatergic stimulation frequency control neuronal dynamics. We performed a systematic analysis of over 2.9 million inter-spike intervals across a large multi-parameter sweep of NMDA kinetics, glutamatergic stimulation frequency, and GABAergic modulation. Dynamical behavior was characterized using entropy–Lyapunov correlation analysis and frequency-dependent bifurcation analysis, and CaMKII phosphorylation was quantified to link kinetic regimes to downstream plasticity signaling.

Results: The analysis revealed two mechanistically distinct pathways to firing irregularity. Pathway 1 (rapid-deactivation irregularity) emerged under relatively fast NMDA deactivation combined with specific input-frequency conditions, producing deterministic chaos with compromised information encoding. Pathway 2 (prolonged-activation irregularity) resulted from slow NMDA deactivation under weak drive, creating irregularity through sustained receptor activation and calcium influx. An optimal kinetic window emerged at βNMDA = 0.042 ms−1, maximizing information transfer (0.275 bits) while maintaining stable dynamics. Entropy–Lyapunov correlation analysis confirmed deterministic chaos, and frequency-dependent bifurcation analysis demonstrated progressive narrowing and displacement of chaotic windows across the analyzed βNMDA range as stimulation frequency increased. GABAergic inhibition provided frequency-selective stabilization, expanding the stable parameter space by 34.2% while preserving gamma oscillations. CaMKII phosphorylation analysis revealed that prolonged NMDA activation maintained elevated phosphorylation levels, creating conditions for pathological long-term potentiation.

Discussion: These findings establish NMDA receptor kinetics as fundamental controllers of cortical excitability and information processing. The dual-pathway framework provides mechanistic insights into addiction-related memory formation, where prolonged NMDA activation enables pathological plasticity, and into visual processing disorders, where altered kinetics disrupt retinal function and cortical oscillatory balance. The identification of optimal kinetic windows and frequency-selective GABA modulation suggests therapeutic strategies based on kinetically specific interventions for neuropsychiatric disorders involving NMDA dysfunction.

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

Repository

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

borjkhani/Bifurcation_NMDA_FCN

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 006879efec965df60df84b7d103bc00b714a8983, 27 February 2026
Languages: Python (9)
Size: 11 files, 9 scripts
Software Heritage: not archived
Found in: “Reproducibility statement”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (9 files), NumPy (9 files), SciPy (5 files), seaborn (5 files), pandas (2 files), scikit-learn (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
11 files
At the source: github.com/borjkhani/Bifurcation_NMDA_FCN

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

Tracing map

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

What the map holds:

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

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

Data

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

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/borjkhani/Bifurcation_NMDA_FCN.

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 2, 28 September 2026

  • Authors: added Mehdi Borjkhani (0000-0003-0469-0902); Hadi Borjkhani (0000-0001-5495-3964); removed Mehdi Borjkhani; Hadi Borjkhani
  • Funding: added European Commission; Fundacja na rzecz Nauki Polskiej

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 5 authors, 8 keywords, 82 references.

Cite

This paper

Borjkhani, M., Borjkhani, H., Sharif, M. A., Bahrami, F., & Janahmadi, M. (2026). NMDA receptor kinetics drive distinct routes to chaotic firing in pyramidal neurons. Frontiers in computational neuroscience, 20, 1753444. https://doi.org/10.3389/fncom.2026.1753444

BibTeX

@article{borjkhani2026nmda,
author = {Borjkhani, Mehdi and Borjkhani, Hadi and Sharif, Morteza A and Bahrami, Fariba and Janahmadi, Mahyar},
title = {{NMDA receptor kinetics drive distinct routes to chaotic firing in pyramidal neurons}},
journal = {Frontiers in computational neuroscience},
year = {2026},
month = jun,
volume = {20},
pages = {1753444},
publisher = {Frontiers Media SA},
issn = {1662-5188},
doi = {10.3389/fncom.2026.1753444},
url = {https://doi.org/10.3389/fncom.2026.1753444},
pmid = {42318008},
pmcid = {PMC13272316}
}

RIS

TY - JOUR
AU - Borjkhani, Mehdi
AU - Borjkhani, Hadi
AU - Sharif, Morteza A
AU - Bahrami, Fariba
AU - Janahmadi, Mahyar
TI - NMDA receptor kinetics drive distinct routes to chaotic firing in pyramidal neurons
T2 - Frontiers in computational neuroscience
J2 - Front Comput Neurosci
PY - 2026
DA - 2026/06/03
VL - 20
SP - 1753444
SN - 1662-5188
PB - Frontiers Media SA
DO - 10.3389/fncom.2026.1753444
UR - https://doi.org/10.3389/fncom.2026.1753444
LA - en
ER -

CSL-JSON

{
"id": "10.3389/fncom.2026.1753444",
"type": "article-journal",
"title": "NMDA receptor kinetics drive distinct routes to chaotic firing in pyramidal neurons",
"container-title": "Frontiers in computational neuroscience",
"author": [
{
"family": "Borjkhani",
"given": "Mehdi"
},
{
"family": "Borjkhani",
"given": "Hadi"
},
{
"family": "Sharif",
"given": "Morteza A"
},
{
"family": "Bahrami",
"given": "Fariba"
},
{
"family": "Janahmadi",
"given": "Mahyar"
}
],
"container-title-short": "Front Comput Neurosci",
"volume": "20",
"page": "1753444",
"DOI": "10.3389/fncom.2026.1753444",
"PMID": "42318008",
"PMCID": "PMC13272316",
"ISSN": "1662-5188",
"publisher": "Frontiers Media SA",
"URL": "https://doi.org/10.3389/fncom.2026.1753444",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
3
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41598-026-42120-y [code]
Qualitative EEG abnormalities in ASD reflect inhibition-dominated brain dynamics.
Journal: Scientific reports
In common: seaborn, pandas, SciPy, 2 other tools, 3 references
[2] doi:10.1016/j.patter.2026.101619 [code]
Sampling bias corrections for discrete and Gaussian partial information decompositions.
Journal: Patterns (New York, N.Y.)
In common: seaborn, scikit-learn, pandas, 3 other tools, 2 references
[3] doi:10.7554/elife.105482 [code]
Principles of gamma synchrony predict figure-ground perception in texture stimuli.
Journal: eLife
In common: seaborn, scikit-learn, pandas, 3 other tools, 2 references
[4] doi:10.1162/imag.a.1252 [code]
Does the brain's E:I balance really shape long-range temporal correlations? Lessons learned from 3T MRI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: seaborn, scikit-learn, pandas, 3 other tools, 2 references
[5] doi:10.1371/journal.pcbi.1014304 [code]
Linking reduced prefrontal microcircuit inhibition in schizophrenia to EEG biomarkers in silico.
Journal: PLoS computational biology
In common: scikit-learn, pandas, SciPy, 2 other tools, 2 references
[6] doi:10.1371/journal.pcbi.1014391 [code]
Multi-stable oscillations in cortical networks with two classes of inhibition.
Journal: PLoS computational biology
In common: seaborn, scikit-learn, pandas, 3 other tools, none (in silico), 1 reference
[7] doi:10.1371/journal.pcbi.1014571 [code]
SynAPSeg: A novel dataset and image analysis framework for deep learning-based synapse detection and quantification.
Journal: PLoS computational biology
In common: seaborn, scikit-learn, pandas, 3 other tools, 1 reference
[8] doi:10.1038/s41467-026-74227-1 [code]
Age-related changes in behavioural and neural variability in a decision-making task.
Journal: Nature communications
In common: seaborn, scikit-learn, pandas, 3 other tools, 1 reference
[9] doi:10.1038/s41586-026-10331-y [code]
Active dissociation of intracortical spiking and high gamma activity.
Journal: Nature
In common: scikit-learn, pandas, SciPy, 2 other tools, 2 references
[10] 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: seaborn, scikit-learn, pandas, 3 other tools, cellular / molecular, 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.