691 lines
30 KiB
Python
691 lines
30 KiB
Python
|
# Tripartite Synapse - Multi-Scale Computational Model
|
||
|
|
# =====================================================
|
||
|
|
# Presynaptic + Postsynaptic perspectives, fully integrated.
|
||
|
|
#
|
||
|
|
# Change log:
|
||
|
|
# ORIG - present from the original document
|
||
|
|
# NEW - added in the missing-behavior integration pass
|
||
|
|
# DET - deterministic Ca2+-driven vesicle release
|
||
|
|
# NKA - explicit Na/K-ATPase V_pre decay and ATP cost
|
||
|
|
# POST-ATP - postsynaptic Ca2+ dynamics and ATP loop
|
||
|
|
# FIX - corrections applied in this pass:
|
||
|
|
# * NT_released_this_window accumulator (was missing entirely)
|
||
|
|
# * k_rec_fast / k_rec_slow converted to /s, recruitment moved to Loop 2
|
||
|
|
# * dt_slow_s added
|
||
|
|
# * mGluR now reads NT_released_this_window (not NT_cleft snapshot)
|
||
|
|
# * IP3 now reads NT_released_this_window (not cleared_NT residual)
|
||
|
|
# * wave_active flag + conversion_efficiency boost on astrocyte wave
|
||
|
|
# * CDI rise gated to spike window only
|
||
|
|
#
|
||
|
|
# Clock structure:
|
||
|
|
# Loop 1 - dt = 1 ms (Ca2+, vesicle release, traces, postsynaptic fast)
|
||
|
|
# Loop 2 - dt = 1000 ms (astrocyte clearance, eCB, mGluR, recruitment)
|
||
|
|
# Loop 3 - dt = 60000 ms (glutamine shuttle, metabolic health)
|
||
|
|
#
|
||
|
|
# =======================================================================
|
||
|
|
# THREE CLOSED LOOPS
|
||
|
|
# =======================================================================
|
||
|
|
#
|
||
|
|
# PRESYNAPTIC:
|
||
|
|
# NT loop : release (ms) -> cleft -> astrocyte clearance (s) ->
|
||
|
|
# glutamine shuttle (min) -> RP refill -> RRP -> release
|
||
|
|
# Ca2+ loop : VGCC influx (ms) -> Tr_Ca -> recruitment speed (s) ->
|
||
|
|
# eCB retrograde from post (s) -> VGCC suppression
|
||
|
|
# ATP loop : NKA + pump costs (ms) -> ATP_demand (min) -> ATP_level ->
|
||
|
|
# pump_scale -> Ca2+ clearance rate -> CDI recovery
|
||
|
|
#
|
||
|
|
# POSTSYNAPTIC:
|
||
|
|
# NT detection loop : NT_cleft -> AMPA -> V_post -> desensitization ->
|
||
|
|
# reduces next response
|
||
|
|
# Ca2+ coincidence : NMDA (NT + V_post) -> Ca_post -> eCB -> pre brake
|
||
|
|
# ATP loop : NKA + PMCA costs (ms) -> ATP_demand_post (min) ->
|
||
|
|
# ATP_level_post -> pump_scale_post -> Ca_post clearance
|
||
|
|
#
|
||
|
|
# SHARED:
|
||
|
|
# eCB_level : post synthesises -> pre reads (retrograde brake)
|
||
|
|
# NT_cleft : pre releases -> post detects -> astrocyte clears
|
||
|
|
# Glucose_level : astrocyte supplies both sides from same budget
|
||
|
|
#
|
||
|
|
# =======================================================================
|
||
|
|
# METABOLIC SILENCING CASCADE (presynaptic)
|
||
|
|
# =======================================================================
|
||
|
|
# [CASCADE 1] HIGH FIRING -> VESICLE DEPLETION (~seconds)
|
||
|
|
# release rate >> recruitment rate -> N_RRP -> 0
|
||
|
|
# [CASCADE 2] HIGH FIRING -> ATP DEPLETION (~minutes)
|
||
|
|
# NKA + PMCA + docking demand > glucose-driven supply
|
||
|
|
# [CASCADE 3] LOW ATP -> PUMP FAILURE
|
||
|
|
# pump_scale = Hill(ATP_level) -> cleared_PMCA/SERCA fall
|
||
|
|
# [CASCADE 4] PUMP FAILURE -> RESIDUAL Ca2+ STAYS HIGH
|
||
|
|
# Ca_micro persists between spikes
|
||
|
|
# [CASCADE 5] RESIDUAL Ca2+ -> CDI LOCKS VGCCs SHUT
|
||
|
|
# CDI rise (spike only) + recovery blocked by Ca2+ -> CDI -> 1
|
||
|
|
# [CASCADE 6] SYNAPSE SILENCES (excitotoxicity protection)
|
||
|
|
# effective_conductance = N_VGCC*(1-eCB)*(1-CDI)*(1-mGluR*alpha)
|
||
|
|
# -> 0; NCX auto-reset when drive stops
|
||
|
|
#
|
||
|
|
# POSTSYNAPTIC ATP CASCADE (no CDI equivalent -> dangerous):
|
||
|
|
# [POST-ATP 1] HIGH V_post + NMDA -> ATP_demand_post rises
|
||
|
|
# [POST-ATP 2] ATP_level_post falls -> pump_scale_post falls
|
||
|
|
# [POST-ATP 3] Ca_post clearance slows -> Ca_post stays elevated
|
||
|
|
# [POST-ATP 4] Ca_post > eCB_threshold without real coincidence
|
||
|
|
# -> false retrograde signal suppresses presynapse
|
||
|
|
# [POST-ATP 5] Critically low ATP_post -> runaway Ca_post -> excitotoxicity
|
||
|
|
# =======================================================================
|
||
|
|
|
||
|
|
import numpy as np
|
||
|
|
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
# CLOCK
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
dt = 1.0 # ms
|
||
|
|
dt_slow = 1000.0 # ms
|
||
|
|
dt_meta = 60_000.0 # ms
|
||
|
|
|
||
|
|
High_Freq_Multiplier = int(dt_slow / dt) # 1000
|
||
|
|
Metabolic_Multiplier = int(dt_meta / dt) # 60000
|
||
|
|
|
||
|
|
dt_s = dt / 1000.0 # 0.001 s/step - for /s rate constants in Loop 1
|
||
|
|
dt_slow_s = dt_slow / 1000.0 # 1.0 s/step - for /s rate constants in Loop 2
|
||
|
|
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
# PRESYNAPTIC PARAMETERS
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
|
||
|
|
# -- Voltage / membrane --
|
||
|
|
tau_V_pre = 2.0 # ms - AP waveform decay (Na/K-ATPase recharge)
|
||
|
|
V_pre_peak = 1.0 # a.u. - normalised AP peak
|
||
|
|
V_rest = 0.0 # a.u. - resting potential
|
||
|
|
V_pre_voltage = -10.0 # mV - driving force for compute_flux
|
||
|
|
|
||
|
|
NKA_cost_per_AP = 0.002 # ATP units per AP (dominant drain at high rates)
|
||
|
|
|
||
|
|
# -- Ca2+ influx & buffering --
|
||
|
|
N_VGCC = 100 # number of VGCCs (ceiling of effective_conductance)
|
||
|
|
k_flux = 0.05 # Ca2+ influx per open channel per unit driving force
|
||
|
|
B_total = 1.0 # total buffer capacity (normalised)
|
||
|
|
tau_buffer_rebind = 200.0 # ms - buffer recharge time constant
|
||
|
|
|
||
|
|
# -- Ca2+ clearance (/ms constants) --
|
||
|
|
k_PMCA = 0.03 # ATP-dependent primary pump
|
||
|
|
k_NCX = 0.10 # ATP-independent floor
|
||
|
|
k_SERCA = 0.01 # ATP-dependent ER pump
|
||
|
|
ATP_half = 0.3 # Hill half-saturation for presynaptic pumps
|
||
|
|
|
||
|
|
ATP_cost_PMCA = 0.0005 # ATP per unit Ca2+ extruded by PMCA
|
||
|
|
ATP_cost_SERCA = 0.0002 # ATP per unit Ca2+ pumped into ER
|
||
|
|
ATP_cost_docking = 0.001 # ATP per vesicle docked (RP->RRP)
|
||
|
|
|
||
|
|
# -- Deterministic release (Hill + NT suppression) --
|
||
|
|
k_rel = 0.5 # max releasable fraction of RRP per spike
|
||
|
|
KD_rel = 1.0 # half-saturation [Ca2+]
|
||
|
|
n_rel = 4 # Hill cooperativity (synaptotagmin-1)
|
||
|
|
NT_suppression_weight = 0.3 # max NT_cleft brake on release fraction
|
||
|
|
NT_suppression_sat = 50.0 # NT_cleft level that saturates suppression
|
||
|
|
|
||
|
|
# -- CDI --
|
||
|
|
k_CDI_rise = 0.8 # /s - CDI build rate (applied * dt_s, spike only)
|
||
|
|
Ca_micro_saturation = 2.0 # normalisation ceiling for CDI recovery
|
||
|
|
k_CDI_rec = 0.015 # /s - CDI de-inactivation rate (applied * dt_s)
|
||
|
|
|
||
|
|
# -- Vesicle pools --
|
||
|
|
Max_RRP = 20
|
||
|
|
Max_RP = 200
|
||
|
|
|
||
|
|
# -- Calcium trace --
|
||
|
|
tau_Tr_Ca = 1000.0 # ms
|
||
|
|
T_high = 0.6 # Tr_Ca threshold -> fast recruitment
|
||
|
|
T_low = 0.2 # Tr_Ca threshold -> slow recruitment
|
||
|
|
|
||
|
|
# -- RP->RRP recruitment (/s, runs in Loop 2) --
|
||
|
|
k_rec_fast = 5.0 # /s - fast recruitment (at Tr_Ca > T_high)
|
||
|
|
k_rec_slow = 0.5 # /s - slow recruitment (at Tr_Ca < T_low)
|
||
|
|
|
||
|
|
# -- NT accumulator for Loop 2 signals --
|
||
|
|
NT_window_sat = 40.0 # vesicles/s that saturates mGluR and IP3
|
||
|
|
# at 20 Hz releasing ~2/spike = 40/s
|
||
|
|
|
||
|
|
# -- eCB retrograde brake --
|
||
|
|
tau_eCB_rise = 2000.0
|
||
|
|
tau_eCB_decay = 10_000.0
|
||
|
|
eCB_threshold = 0.7 # Ca_post level that triggers eCB synthesis
|
||
|
|
|
||
|
|
# -- mGluR presynaptic autoreceptor --
|
||
|
|
Km_mGluR = 0.5
|
||
|
|
tau_mGluR = 2000.0 # ms
|
||
|
|
alpha_mGluR = 0.4 # max fractional VGCC suppression
|
||
|
|
|
||
|
|
# -- Astrocyte / IP3 --
|
||
|
|
tau_IP3 = 3000.0 # ms
|
||
|
|
IP3_threshold = 0.8
|
||
|
|
wave_boost = 0.2 # conversion_efficiency boost when wave fires
|
||
|
|
tau_wave_decay = 2 # metabolic cycles before boost decays back
|
||
|
|
|
||
|
|
# -- Glutamine shuttle --
|
||
|
|
conversion_efficiency_base = 0.8
|
||
|
|
|
||
|
|
# -- NT cleft --
|
||
|
|
tau_NT_decay = 5.0 # ms
|
||
|
|
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
# POSTSYNAPTIC PARAMETERS
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
|
||
|
|
# -- NMDA coincidence detection --
|
||
|
|
k_NMDA = 0.08 # Ca_post influx per unit NT * (1 - Mg_block) per ms
|
||
|
|
V_NMDA_half = 0.3 # V_post at which Mg block is 50% lifted
|
||
|
|
|
||
|
|
# -- Ca_post clearance --
|
||
|
|
k_Ca_post_clear = 0.05 # /ms - ATP-dependent PMCA in spine
|
||
|
|
k_Ca_post_NCX = 0.02 # /ms - ATP-independent NCX floor
|
||
|
|
ATP_half_post = 0.3 # Hill half-saturation for postsynaptic pumps
|
||
|
|
|
||
|
|
# -- Postsynaptic ATP costs --
|
||
|
|
NKA_cost_per_bAP_post = 0.002 # ATP per unit V_post per s (continuous)
|
||
|
|
ATP_cost_Ca_post_pump = 0.0005 # ATP per unit Ca_post cleared
|
||
|
|
ATP_demand_scale_post = 50.0 # normalisation (same as presynaptic)
|
||
|
|
|
||
|
|
# -- Receptor desensitization --
|
||
|
|
tau_membrane = 20.0 # ms
|
||
|
|
tau_desensitization = 500.0 # ms
|
||
|
|
|
||
|
|
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
# HELPER FUNCTIONS
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
|
||
|
|
def compute_flux(conductance, voltage):
|
||
|
|
return k_flux * conductance * abs(voltage)
|
||
|
|
|
||
|
|
|
||
|
|
def deterministic_release(N_RRP, Ca_micro, NT_cleft):
|
||
|
|
# Hill equation: Ca2+ sensor cooperativity (synaptotagmin-1, n=4)
|
||
|
|
Ca_n = Ca_micro ** n_rel
|
||
|
|
release_frac = k_rel * Ca_n / (Ca_n + KD_rel ** n_rel)
|
||
|
|
# NT suppression: physical crowding + fast local autoreceptors
|
||
|
|
NT_norm = min(1.0, NT_cleft / NT_suppression_sat)
|
||
|
|
release_frac = release_frac * (1.0 - NT_suppression_weight * NT_norm)
|
||
|
|
release_frac = np.clip(release_frac, 0.0, 1.0)
|
||
|
|
return max(0.0, release_frac * N_RRP)
|
||
|
|
|
||
|
|
|
||
|
|
def map_trace_to_speed(Tr_Ca):
|
||
|
|
# Returns /s recruitment rate based on Tr_Ca level
|
||
|
|
if Tr_Ca > T_high:
|
||
|
|
return k_rec_fast
|
||
|
|
elif Tr_Ca < T_low:
|
||
|
|
return k_rec_slow
|
||
|
|
else:
|
||
|
|
t = (Tr_Ca - T_low) / (T_high - T_low)
|
||
|
|
return k_rec_slow + t * (k_rec_fast - k_rec_slow)
|
||
|
|
|
||
|
|
|
||
|
|
def compute_pump_atp_factor(atp, atp_half):
|
||
|
|
# Hill function: ATP gates pump speed (shared by pre and post)
|
||
|
|
return (atp ** 2) / (atp ** 2 + atp_half ** 2)
|
||
|
|
|
||
|
|
|
||
|
|
def compute_EPSP(receptor_conductance):
|
||
|
|
return receptor_conductance * 0.1
|
||
|
|
|
||
|
|
|
||
|
|
def compute_astrocyte_metabolic_health(Glucose_level, ATP_demand_accumulated,
|
||
|
|
demand_scale=50.0):
|
||
|
|
# Converts glucose supply and accumulated demand into ATP_level (0->1)
|
||
|
|
# and conversion_efficiency (0->1). Both sides use this function with
|
||
|
|
# their own demand accumulators but the same Glucose_level — shared
|
||
|
|
# metabolic vulnerability.
|
||
|
|
health = np.clip(Glucose_level - ATP_demand_accumulated / demand_scale,
|
||
|
|
0.0, 1.0)
|
||
|
|
return health, health # (conversion_efficiency, ATP_level)
|
||
|
|
|
||
|
|
|
||
|
|
def trigger_slow_astrocyte_calcium_wave():
|
||
|
|
# Placeholder - gliotransmitter release over ~10 s
|
||
|
|
pass
|
||
|
|
|
||
|
|
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
# STATE VARIABLES
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
|
||
|
|
# -- Presynaptic membrane --
|
||
|
|
V_pre_state = 0.0
|
||
|
|
|
||
|
|
# -- Presynaptic Ca2+ --
|
||
|
|
Ca_micro = 0.0
|
||
|
|
Ca_ER = 0.5
|
||
|
|
Ca_buffer_bound = 0.0
|
||
|
|
B_free = B_total
|
||
|
|
|
||
|
|
# -- CDI --
|
||
|
|
CDI_factor = 0.0
|
||
|
|
|
||
|
|
# -- Vesicle pools --
|
||
|
|
N_RRP = 15.0
|
||
|
|
N_RP = 150.0
|
||
|
|
|
||
|
|
# -- Calcium trace --
|
||
|
|
Tr_Ca = 0.0
|
||
|
|
|
||
|
|
# -- NT cleft --
|
||
|
|
NT_cleft = 0.0
|
||
|
|
|
||
|
|
# -- NT accumulator for slow signals --
|
||
|
|
# FIX: this was missing. Accumulates every ms in Loop 1,
|
||
|
|
# consumed by mGluR and IP3 in Loop 2, reset each second.
|
||
|
|
NT_released_this_window = 0.0
|
||
|
|
|
||
|
|
# -- Postsynaptic membrane + receptors --
|
||
|
|
V_post = 0.0
|
||
|
|
receptor_conductance = 0.0
|
||
|
|
Desensitization_level = 0.0
|
||
|
|
V_post_history = []
|
||
|
|
|
||
|
|
# -- Postsynaptic Ca2+ (spine compartment) --
|
||
|
|
Ca_post = 0.0
|
||
|
|
# Driven by NMDA coincidence (NT + V_post). Cleared by PMCA (ATP-gated)
|
||
|
|
# and NCX (always). Drives eCB synthesis. No CDI equivalent ->
|
||
|
|
# elevated Ca_post under ATP failure has no self-limiting mechanism.
|
||
|
|
|
||
|
|
# -- Retrograde / autoreceptor --
|
||
|
|
eCB_level = 0.0
|
||
|
|
mGluR_activation = 0.0
|
||
|
|
|
||
|
|
# -- Astrocyte --
|
||
|
|
IP3 = 0.0
|
||
|
|
wave_active = 0 # countdown: cycles remaining of wave boost
|
||
|
|
Glutamine_pool = 50.0
|
||
|
|
|
||
|
|
# -- Presynaptic ATP --
|
||
|
|
ATP_level = 1.0
|
||
|
|
ATP_demand = 0.0
|
||
|
|
conversion_efficiency = conversion_efficiency_base
|
||
|
|
Glucose_level = 1.0 # set < 1.0 to engage metabolic silencing
|
||
|
|
|
||
|
|
# -- Postsynaptic ATP --
|
||
|
|
ATP_level_post = 1.0 # separate pool; same glucose budget as presynaptic
|
||
|
|
ATP_demand_post = 0.0 # accumulates from NKA (V_post) and PMCA (Ca_post)
|
||
|
|
|
||
|
|
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
# MAIN SIMULATION LOOP
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
|
||
|
|
def run_simulation(spike_train, total_steps, bAP_train=None):
|
||
|
|
"""
|
||
|
|
spike_train : list of int - presynaptic AP timestep indices
|
||
|
|
total_steps : int
|
||
|
|
bAP_train : list of int - postsynaptic bAP timestep indices (optional)
|
||
|
|
if None, no bAPs are delivered
|
||
|
|
"""
|
||
|
|
|
||
|
|
global V_pre_state
|
||
|
|
global Ca_micro, Ca_ER, Ca_buffer_bound, B_free
|
||
|
|
global CDI_factor
|
||
|
|
global N_RRP, N_RP, Tr_Ca, NT_cleft, NT_released_this_window
|
||
|
|
global V_post, receptor_conductance, Desensitization_level, V_post_history
|
||
|
|
global Ca_post
|
||
|
|
global eCB_level, mGluR_activation
|
||
|
|
global IP3, wave_active, Glutamine_pool
|
||
|
|
global ATP_level, ATP_demand, conversion_efficiency, Glucose_level
|
||
|
|
global ATP_level_post, ATP_demand_post
|
||
|
|
|
||
|
|
log = {k: [] for k in [
|
||
|
|
"V_pre_state", "Ca_micro", "Ca_ER", "CDI_factor", "B_free",
|
||
|
|
"N_RRP", "N_RP", "Tr_Ca", "NT_cleft",
|
||
|
|
"V_post", "Ca_post", "eCB_level", "mGluR_activation",
|
||
|
|
"released_NT", "ATP_level", "ATP_demand",
|
||
|
|
"ATP_level_post", "ATP_demand_post",
|
||
|
|
]}
|
||
|
|
|
||
|
|
spike_set = set(spike_train)
|
||
|
|
bAP_set = set(bAP_train) if bAP_train else set()
|
||
|
|
|
||
|
|
for step in range(total_steps):
|
||
|
|
|
||
|
|
# ==============================================================
|
||
|
|
# LOOP 1 — HIGH-FREQUENCY (dt = 1 ms)
|
||
|
|
# ==============================================================
|
||
|
|
|
||
|
|
V_pre = 1 if step in spike_set else 0
|
||
|
|
bAP = 1 if step in bAP_set else 0
|
||
|
|
released_NT = 0.0
|
||
|
|
|
||
|
|
# -- 1A. PRESYNAPTIC MEMBRANE / Na-K-ATPase -------------------
|
||
|
|
# AP fires: membrane jumps to peak, then decays with tau_V_pre.
|
||
|
|
# Ca2+ influx uses V_pre_state (continuous) not binary V_pre,
|
||
|
|
# giving a temporal influx profile that tapers as membrane repolarises.
|
||
|
|
if V_pre == 1:
|
||
|
|
V_pre_state = V_pre_peak
|
||
|
|
ATP_demand += NKA_cost_per_AP # dominant presynaptic ATP cost
|
||
|
|
|
||
|
|
V_pre_state += (V_rest - V_pre_state) * dt / tau_V_pre
|
||
|
|
|
||
|
|
# -- 1B. PRESYNAPTIC Ca2+ INFLUX ------------------------------
|
||
|
|
# Three multiplicative brakes on effective_conductance:
|
||
|
|
# eCB_level : retrograde brake from postsynapse (Loop 2)
|
||
|
|
# CDI_factor : Ca2+-dependent inactivation (below)
|
||
|
|
# mGluR_activation : autoreceptor brake (Loop 2)
|
||
|
|
effective_conductance = (
|
||
|
|
N_VGCC
|
||
|
|
* (1.0 - eCB_level)
|
||
|
|
* (1.0 - CDI_factor)
|
||
|
|
* (1.0 - mGluR_activation * alpha_mGluR)
|
||
|
|
)
|
||
|
|
raw_influx = compute_flux(effective_conductance, V_pre_state)
|
||
|
|
|
||
|
|
# Buffer proteins capture a fraction immediately (fast sponge).
|
||
|
|
# B_free -> 0 during sustained bursting -> capture_fraction -> 0
|
||
|
|
# -> full raw_influx enters Ca_micro (CASCADE 4 acceleration).
|
||
|
|
capture_fraction = B_free / B_total
|
||
|
|
captured = raw_influx * capture_fraction
|
||
|
|
B_free = max(0.0, B_free - captured)
|
||
|
|
Ca_buffer_bound += captured
|
||
|
|
Ca_micro += (raw_influx - captured)
|
||
|
|
|
||
|
|
# -- 1C. VESICLE RELEASE --------------------------------------
|
||
|
|
# Deterministic: Hill Ca2+ sensor * NT suppression * N_RRP.
|
||
|
|
# Runs every ms that Ca_micro > 0 (release profile follows Ca2+
|
||
|
|
# transient, not locked to spike flag).
|
||
|
|
if N_RRP > 0 and Ca_micro > 0:
|
||
|
|
released_NT = deterministic_release(N_RRP, Ca_micro, NT_cleft)
|
||
|
|
released_NT = min(released_NT, N_RRP)
|
||
|
|
N_RRP -= released_NT
|
||
|
|
NT_cleft += released_NT
|
||
|
|
# FIX: accumulate for Loop 2 mGluR and IP3 signals.
|
||
|
|
# This is the only correct way to feed slow signals from fast
|
||
|
|
# events — snapshot of NT_cleft at Loop 2 time would be ~0
|
||
|
|
# because passive diffusion has already cleared it.
|
||
|
|
NT_released_this_window += released_NT
|
||
|
|
|
||
|
|
# Passive NT diffusion out of cleft each ms.
|
||
|
|
NT_cleft *= (1.0 - dt / tau_NT_decay)
|
||
|
|
NT_cleft = max(0.0, NT_cleft)
|
||
|
|
|
||
|
|
# -- 1D. PRESYNAPTIC Ca2+ CLEARANCE ---------------------------
|
||
|
|
# pump_scale: Hill(ATP_level) — bridges Loop 3 ATP to Loop 1 clearance.
|
||
|
|
# NCX is ATP-independent (floor); PMCA and SERCA are ATP-gated.
|
||
|
|
pump_scale = compute_pump_atp_factor(ATP_level, ATP_half)
|
||
|
|
cleared_PMCA = k_PMCA * Ca_micro * pump_scale
|
||
|
|
cleared_NCX = k_NCX * Ca_micro
|
||
|
|
cleared_SERCA = k_SERCA * Ca_micro * pump_scale
|
||
|
|
|
||
|
|
Ca_micro -= (cleared_PMCA + cleared_NCX + cleared_SERCA)
|
||
|
|
Ca_micro = max(0.0, Ca_micro)
|
||
|
|
Ca_ER += cleared_SERCA
|
||
|
|
|
||
|
|
ATP_demand += ATP_cost_PMCA * cleared_PMCA
|
||
|
|
ATP_demand += ATP_cost_SERCA * cleared_SERCA
|
||
|
|
|
||
|
|
# Buffer recharge: bound Ca2+ slowly re-releases back to cytosol.
|
||
|
|
# During pump failure this sustains Ca_micro elevation (CASCADE 4).
|
||
|
|
rebind_flux = Ca_buffer_bound * dt / tau_buffer_rebind
|
||
|
|
Ca_micro += rebind_flux
|
||
|
|
Ca_buffer_bound = max(0.0, Ca_buffer_bound - rebind_flux)
|
||
|
|
B_free = B_total - Ca_buffer_bound
|
||
|
|
|
||
|
|
# -- 1E. CDI — RISE (spike only) AND RECOVERY (every ms) ------
|
||
|
|
# RISE: Ca2+ entering through open channels inactivates them locally.
|
||
|
|
# Gated to spike window — requires channels to be open.
|
||
|
|
# (Running every ms was wrong: CDI needs Ca2+ flowing through
|
||
|
|
# the channel, not ambient cytosolic Ca2+.)
|
||
|
|
if V_pre == 1:
|
||
|
|
CDI_factor += k_CDI_rise * Ca_micro * dt_s
|
||
|
|
|
||
|
|
# RECOVERY: continuous, suppressed when Ca_micro is high.
|
||
|
|
# Self-locking: pump failure -> Ca_micro high -> recovery ~0
|
||
|
|
# -> CDI_factor -> 1 -> effective_conductance -> 0 (CASCADE 5-6).
|
||
|
|
CDI_recovery_rate = k_CDI_rec * (1.0 - Ca_micro / Ca_micro_saturation)
|
||
|
|
CDI_factor = np.clip(CDI_factor - CDI_recovery_rate * dt_s, 0.0, 1.0)
|
||
|
|
|
||
|
|
# -- 1F. CALCIUM TRACE ----------------------------------------
|
||
|
|
# Leaky integrator — integrates full Ca2+ waveform every ms
|
||
|
|
# including inter-spike clearance. Drives Loop 2 recruitment speed.
|
||
|
|
Tr_Ca = Tr_Ca + (Ca_micro - Tr_Ca / tau_Tr_Ca) * dt
|
||
|
|
|
||
|
|
# -- 1G. POSTSYNAPTIC: NT DETECTION & AMPA --------------------
|
||
|
|
# Desensitization reduces effective NT — sustained NT exposure
|
||
|
|
# progressively silences receptors (postsynaptic equivalent of CDI).
|
||
|
|
effective_NT = released_NT * (1.0 - Desensitization_level)
|
||
|
|
receptor_conductance += effective_NT * 0.05
|
||
|
|
receptor_conductance *= (1.0 - dt / tau_membrane)
|
||
|
|
|
||
|
|
V_post += compute_EPSP(receptor_conductance) - (V_post / tau_membrane) * dt
|
||
|
|
V_post = max(0.0, V_post)
|
||
|
|
|
||
|
|
Desensitization_level += NT_cleft * 0.001 * dt
|
||
|
|
Desensitization_level -= (Desensitization_level / tau_desensitization) * dt
|
||
|
|
Desensitization_level = np.clip(Desensitization_level, 0.0, 1.0)
|
||
|
|
|
||
|
|
V_post_history.append(V_post)
|
||
|
|
if len(V_post_history) > 5000:
|
||
|
|
V_post_history.pop(0)
|
||
|
|
|
||
|
|
# -- 1H. POSTSYNAPTIC: NMDA COINCIDENCE DETECTION -------------
|
||
|
|
# Ca_post enters only when BOTH conditions hold simultaneously:
|
||
|
|
# (1) NT_cleft > 0 — ligand gate (glutamate present)
|
||
|
|
# (2) V_post elevated — voltage gate (Mg2+ block lifted)
|
||
|
|
# Mg block removal is a sigmoid of V_post.
|
||
|
|
# bAP (backpropagating AP) boosts V_post further, enabling
|
||
|
|
# full NMDA opening and larger Ca_post surge.
|
||
|
|
V_post_effective = V_post + (bAP * 0.5) # bAP adds depolarisation
|
||
|
|
Mg_block_removal = V_post_effective / (V_post_effective + V_NMDA_half)
|
||
|
|
NMDA_Ca_influx = k_NMDA * NT_cleft * Mg_block_removal
|
||
|
|
Ca_post += NMDA_Ca_influx
|
||
|
|
|
||
|
|
# Postsynaptic NKA: membrane recharge cost proportional to V_post.
|
||
|
|
# [POST-ATP 1] Dominant postsynaptic ATP drain at high activity.
|
||
|
|
ATP_demand_post += NKA_cost_per_bAP_post * V_post * dt_s
|
||
|
|
|
||
|
|
# -- 1I. POSTSYNAPTIC: Ca_post CLEARANCE ----------------------
|
||
|
|
# pump_scale_post: Hill(ATP_level_post) — same structure as presynaptic.
|
||
|
|
# NCX is ATP-independent floor (enables auto-reset after ATP recovery).
|
||
|
|
# [POST-ATP 3] When pump_scale_post falls, Ca_post stays elevated ->
|
||
|
|
# eCB threshold crossed without genuine coincidence -> false retrograde.
|
||
|
|
pump_scale_post = compute_pump_atp_factor(ATP_level_post, ATP_half_post)
|
||
|
|
cleared_Ca_post_pump = k_Ca_post_clear * Ca_post * pump_scale_post
|
||
|
|
cleared_Ca_post_NCX = k_Ca_post_NCX * Ca_post
|
||
|
|
Ca_post -= (cleared_Ca_post_pump + cleared_Ca_post_NCX)
|
||
|
|
Ca_post = max(0.0, Ca_post)
|
||
|
|
|
||
|
|
# [POST-ATP 2] ATP cost of postsynaptic PMCA.
|
||
|
|
ATP_demand_post += ATP_cost_Ca_post_pump * cleared_Ca_post_pump
|
||
|
|
|
||
|
|
# -- RECORD ---------------------------------------------------
|
||
|
|
log["V_pre_state"].append(V_pre_state)
|
||
|
|
log["Ca_micro"].append(Ca_micro)
|
||
|
|
log["Ca_ER"].append(Ca_ER)
|
||
|
|
log["CDI_factor"].append(CDI_factor)
|
||
|
|
log["B_free"].append(B_free)
|
||
|
|
log["N_RRP"].append(N_RRP)
|
||
|
|
log["N_RP"].append(N_RP)
|
||
|
|
log["Tr_Ca"].append(Tr_Ca)
|
||
|
|
log["NT_cleft"].append(NT_cleft)
|
||
|
|
log["V_post"].append(V_post)
|
||
|
|
log["Ca_post"].append(Ca_post)
|
||
|
|
log["eCB_level"].append(eCB_level)
|
||
|
|
log["mGluR_activation"].append(mGluR_activation)
|
||
|
|
log["released_NT"].append(released_NT)
|
||
|
|
log["ATP_level"].append(ATP_level)
|
||
|
|
log["ATP_demand"].append(ATP_demand)
|
||
|
|
log["ATP_level_post"].append(ATP_level_post)
|
||
|
|
log["ATP_demand_post"].append(ATP_demand_post)
|
||
|
|
|
||
|
|
# ==============================================================
|
||
|
|
# LOOP 2 — SLOW / ASTROCYTE (dt_slow = 1 s)
|
||
|
|
# ==============================================================
|
||
|
|
|
||
|
|
if (step % High_Freq_Multiplier) == 0:
|
||
|
|
|
||
|
|
# Astrocyte EAAT clearance — active NT removal from cleft.
|
||
|
|
cleared_NT = NT_cleft * 0.3
|
||
|
|
NT_cleft = max(0.0, NT_cleft - cleared_NT)
|
||
|
|
|
||
|
|
# FIX: IP3 integrates NT_released_this_window (total release
|
||
|
|
# since last Loop 2), not the post-diffusion NT_cleft residual
|
||
|
|
# which is ~0 by the time Loop 2 runs.
|
||
|
|
IP3 += NT_released_this_window - (IP3 / tau_IP3) * dt_slow
|
||
|
|
IP3 = max(0.0, IP3)
|
||
|
|
|
||
|
|
if IP3 > IP3_threshold:
|
||
|
|
trigger_slow_astrocyte_calcium_wave()
|
||
|
|
# FIX: wave boosts conversion_efficiency in the next mins cycle.
|
||
|
|
# The astrocyte responds to heavy load by upregulating its
|
||
|
|
# recycling machinery — shipping more glutamine back to the
|
||
|
|
# presynapse. Boost decays over tau_wave_decay metabolic cycles.
|
||
|
|
wave_active = tau_wave_decay
|
||
|
|
|
||
|
|
# FIX: mGluR reads NT_released_this_window (accumulated release
|
||
|
|
# load), not NT_cleft snapshot. NT_cleft is ~0 at Loop 2 time
|
||
|
|
# due to diffusion; the accumulator correctly represents the
|
||
|
|
# burst load the autoreceptor has sensed during this window.
|
||
|
|
NT_window_norm = min(1.0, NT_released_this_window / NT_window_sat)
|
||
|
|
mGluR_target = NT_window_norm
|
||
|
|
mGluR_activation += (mGluR_target - mGluR_activation) * (dt_slow / tau_mGluR)
|
||
|
|
mGluR_activation = np.clip(mGluR_activation, 0.0, 1.0)
|
||
|
|
|
||
|
|
# FIX: reset accumulator for next window.
|
||
|
|
NT_released_this_window = 0.0
|
||
|
|
|
||
|
|
# eCB retrograde synthesis: now driven by Ca_post (spine Ca2+),
|
||
|
|
# not V_post_history. The actual eCB synthesis in the spine is
|
||
|
|
# triggered by Ca2+-dependent enzymes (DAGL, PLC), not voltage.
|
||
|
|
# Under normal conditions Ca_post only rises with coincidence.
|
||
|
|
# Under POST-ATP failure Ca_post stays elevated without genuine
|
||
|
|
# coincidence -> false retrograde signal (POST-ATP 4).
|
||
|
|
recent_Ca_post = (np.mean(log["Ca_post"][-2000:])
|
||
|
|
if len(log["Ca_post"]) >= 2000
|
||
|
|
else (np.mean(log["Ca_post"]) if log["Ca_post"] else 0.0))
|
||
|
|
eCB_signal = max(0.0, recent_Ca_post - eCB_threshold)
|
||
|
|
if eCB_signal > 0:
|
||
|
|
eCB_level += eCB_signal * (dt_slow / tau_eCB_rise)
|
||
|
|
else:
|
||
|
|
eCB_level -= eCB_level * (dt_slow / tau_eCB_decay)
|
||
|
|
eCB_level = np.clip(eCB_level, 0.0, 1.0)
|
||
|
|
|
||
|
|
# FIX: RP->RRP recruitment moved here from Loop 1.
|
||
|
|
# Biological timescale: vesicle docking and priming take seconds,
|
||
|
|
# not milliseconds. k_rec_fast/slow are /s; * dt_slow_s = 1.0 s
|
||
|
|
# gives dimensionless per-step fraction — no hidden unit scaling.
|
||
|
|
current_recruitment_rate = map_trace_to_speed(Tr_Ca) # /s
|
||
|
|
refill_amount = (current_recruitment_rate * dt_slow_s
|
||
|
|
* N_RP * (Max_RRP - N_RRP) / Max_RRP)
|
||
|
|
refill_amount = max(0.0, refill_amount)
|
||
|
|
refill_amount = min(refill_amount, N_RP)
|
||
|
|
|
||
|
|
N_RRP = min(N_RRP + refill_amount, Max_RRP)
|
||
|
|
N_RP = max(0.0, N_RP - refill_amount)
|
||
|
|
ATP_demand += ATP_cost_docking * refill_amount
|
||
|
|
|
||
|
|
# ==============================================================
|
||
|
|
# LOOP 3 — METABOLIC (dt_meta = 1 min)
|
||
|
|
# ==============================================================
|
||
|
|
|
||
|
|
if (step % Metabolic_Multiplier) == 0:
|
||
|
|
|
||
|
|
# Presynaptic ATP: glucose supply minus accumulated demand.
|
||
|
|
conversion_efficiency, ATP_level = compute_astrocyte_metabolic_health(
|
||
|
|
Glucose_level, ATP_demand
|
||
|
|
)
|
||
|
|
ATP_demand = 0.0
|
||
|
|
|
||
|
|
# FIX: wave boost applied to conversion_efficiency.
|
||
|
|
# Astrocyte calcium wave (triggered by high IP3) upregulates
|
||
|
|
# glutamine synthetase -> faster NT recycling -> more RP refill.
|
||
|
|
# Boost decays over tau_wave_decay cycles.
|
||
|
|
if wave_active > 0:
|
||
|
|
conversion_efficiency = min(1.0, conversion_efficiency + wave_boost)
|
||
|
|
wave_active -= 1
|
||
|
|
|
||
|
|
# Glutamine shuttle: astrocyte converts cleared NT to glutamine,
|
||
|
|
# presynapse repackages it into vesicles -> N_RP replenished.
|
||
|
|
refill_RP = Glutamine_pool * conversion_efficiency
|
||
|
|
N_RP = min(Max_RP, N_RP + refill_RP)
|
||
|
|
Glutamine_pool = max(0.0, Glutamine_pool - refill_RP)
|
||
|
|
|
||
|
|
# Postsynaptic ATP: same glucose budget, own demand accumulator.
|
||
|
|
# Both sides draw from Glucose_level -> shared metabolic vulnerability.
|
||
|
|
# Presynaptic silence reduces NT -> less NMDA -> less Ca_post ->
|
||
|
|
# less ATP_demand_post: presynaptic protection indirectly
|
||
|
|
# protects the postsynapse.
|
||
|
|
_, ATP_level_post = compute_astrocyte_metabolic_health(
|
||
|
|
Glucose_level, ATP_demand_post, ATP_demand_scale_post
|
||
|
|
)
|
||
|
|
ATP_demand_post = 0.0
|
||
|
|
|
||
|
|
return log
|
||
|
|
|
||
|
|
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
# EXAMPLE USAGE
|
||
|
|
# -----------------------------------------------------------------------
|
||
|
|
|
||
|
|
if __name__ == "__main__":
|
||
|
|
import matplotlib.pyplot as plt
|
||
|
|
|
||
|
|
total_steps = 10_000 # 10 seconds
|
||
|
|
|
||
|
|
# Presynaptic 20 Hz burst for 2 s.
|
||
|
|
spike_train = list(range(0, 2000, 50))
|
||
|
|
|
||
|
|
# Postsynaptic bAPs coincident with every 5th presynaptic spike
|
||
|
|
# (simulates partial coincidence for NMDA activation).
|
||
|
|
bAP_train = list(range(0, 2000, 250))
|
||
|
|
|
||
|
|
results = run_simulation(spike_train, total_steps, bAP_train)
|
||
|
|
t = np.arange(total_steps) * dt
|
||
|
|
|
||
|
|
fig, axes = plt.subplots(8, 1, figsize=(12, 18), sharex=True)
|
||
|
|
fig.suptitle("Tripartite Synapse — Presynaptic + Postsynaptic", fontsize=13)
|
||
|
|
|
||
|
|
axes[0].plot(t, results["V_pre_state"], color="slateblue", lw=0.8)
|
||
|
|
axes[0].set_ylabel("V_pre")
|
||
|
|
axes[0].set_title("Presynaptic membrane (AP waveform)", fontsize=9, loc="left")
|
||
|
|
|
||
|
|
axes[1].plot(t, results["Ca_micro"], color="darkorange", lw=0.8)
|
||
|
|
axes[1].set_ylabel("[Ca2+] pre")
|
||
|
|
axes[1].set_title("CASCADE 4 — presynaptic Ca2+", fontsize=9, loc="left")
|
||
|
|
|
||
|
|
axes[2].plot(t, results["CDI_factor"], color="firebrick", lw=0.8, label="CDI")
|
||
|
|
axes[2].plot(t, results["B_free"], color="steelblue", lw=0.8, label="Buffer free")
|
||
|
|
axes[2].set_ylabel("CDI / Buffer")
|
||
|
|
axes[2].set_title("CASCADE 5 — CDI lock-out", fontsize=9, loc="left")
|
||
|
|
axes[2].legend(fontsize=8)
|
||
|
|
|
||
|
|
axes[3].plot(t, results["N_RRP"], color="teal", lw=0.8, label="RRP")
|
||
|
|
axes[3].plot(t, results["N_RP"], color="purple", lw=0.8, label="RP")
|
||
|
|
axes[3].set_ylabel("Vesicles")
|
||
|
|
axes[3].set_title("CASCADE 1 — vesicle depletion", fontsize=9, loc="left")
|
||
|
|
axes[3].legend(fontsize=8)
|
||
|
|
|
||
|
|
axes[4].plot(t, results["NT_cleft"], color="darkgreen", lw=0.8, label="NT cleft")
|
||
|
|
axes[4].plot(t, results["mGluR_activation"], color="saddlebrown", lw=0.8, label="mGluR")
|
||
|
|
axes[4].plot(t, results["eCB_level"], color="crimson", lw=0.8, label="eCB")
|
||
|
|
axes[4].set_ylabel("Cleft / Feedback")
|
||
|
|
axes[4].set_title("CASCADE 6 — three brakes on conductance", fontsize=9, loc="left")
|
||
|
|
axes[4].legend(fontsize=8)
|
||
|
|
|
||
|
|
axes[5].plot(t, results["V_post"], color="navy", lw=0.8, label="V_post")
|
||
|
|
axes[5].plot(t, results["Ca_post"], color="coral", lw=0.8, label="Ca_post (spine)")
|
||
|
|
axes[5].set_ylabel("Postsynaptic")
|
||
|
|
axes[5].set_title("Postsynaptic potential + NMDA spine Ca2+", fontsize=9, loc="left")
|
||
|
|
axes[5].legend(fontsize=8)
|
||
|
|
|
||
|
|
axes[6].plot(t, results["ATP_level"], color="goldenrod", lw=0.8, label="ATP pre")
|
||
|
|
axes[6].plot(t, results["ATP_level_post"], color="darkorange", lw=0.8, label="ATP post")
|
||
|
|
axes[6].set_ylabel("ATP level")
|
||
|
|
axes[6].set_title("CASCADE 2 / POST-ATP — presynaptic and postsynaptic ATP", fontsize=9, loc="left")
|
||
|
|
axes[6].legend(fontsize=8)
|
||
|
|
|
||
|
|
axes[7].plot(t, results["ATP_demand"], color="tomato", lw=0.8, label="demand pre")
|
||
|
|
axes[7].plot(t, results["ATP_demand_post"], color="orangered", lw=0.8, label="demand post")
|
||
|
|
axes[7].set_ylabel("ATP demand")
|
||
|
|
axes[7].set_title("Accumulated ATP demand (resets each min cycle)", fontsize=9, loc="left")
|
||
|
|
axes[7].set_xlabel("Time (ms)")
|
||
|
|
axes[7].legend(fontsize=8)
|
||
|
|
|
||
|
|
plt.tight_layout()
|
||
|
|
plt.savefig("./synapse_simulation.png", dpi=150)
|
||
|
|
plt.close()
|
||
|
|
print("Done.")
|