714 lines
33 KiB
Python
714 lines
33 KiB
Python
# Tripartite Synapse - Multi-Scale Computational Model
|
|
# =====================================================
|
|
# Presynaptic perspective, with all missing behaviors integrated.
|
|
#
|
|
# Change log:
|
|
# ORIG - present from the original document
|
|
# NEW - added in the missing-behavior integration pass
|
|
# DET - this pass: deterministic Ca2+-driven vesicle release
|
|
# NKA - this pass: explicit Na/K-ATPase V_pre decay and ATP cost
|
|
#
|
|
# Clock structure:
|
|
# Loop 1 - dt = 1 ms (Ca2+, vesicle release, short-term traces)
|
|
# Loop 2 - dt = 1000 ms (astrocyte clearance, eCB, mGluR feedback)
|
|
# Loop 3 - dt = 60000 ms (glutamine shuttle, metabolic health)
|
|
#
|
|
# =======================================================================
|
|
# METABOLIC SILENCING CASCADE - variable map
|
|
# =======================================================================
|
|
# Each step of the cascade is tagged inline with [CASCADE n].
|
|
#
|
|
# [CASCADE 1] HIGH FIRING RATE -> VESICLE DEPLETION (fast, ~seconds)
|
|
# Driver : spike_train density -> released_NT per ms
|
|
# Victim : N_RRP, N_RP
|
|
# Bottleneck : RP->RRP recruitment cannot keep up at high rates.
|
|
# Outcome : N_RRP -> 0, released_NT -> 0, NT_cleft collapses.
|
|
#
|
|
# [CASCADE 2] HIGH FIRING RATE -> ATP DEPLETION (slow, ~minutes)
|
|
# Driver : Na/K-ATPase recharge cost per AP (NKA_cost_per_AP)
|
|
# + PMCA/SERCA pump load + vesicle re-docking
|
|
# Victim : ATP_demand accumulator -> ATP_level
|
|
# Bottleneck : Glucose_level sets replenishment ceiling.
|
|
# Outcome : ATP_level < 1, pump_scale < 1.
|
|
#
|
|
# [CASCADE 3] LOW ATP -> PUMP FAILURE (PMCA / SERCA slow)
|
|
# Driver : pump_scale = Hill(ATP_level)
|
|
# Victim : cleared_PMCA, cleared_SERCA
|
|
# Outcome : total Ca2+ clearance rate drops.
|
|
#
|
|
# [CASCADE 4] PUMP FAILURE -> RESIDUAL Ca2+ STAYS HIGH
|
|
# Driver : reduced clearance + saturated buffer
|
|
# Victim : Ca_micro
|
|
# Outcome : Ca_micro persists between spikes.
|
|
#
|
|
# [CASCADE 5] RESIDUAL Ca2+ -> CDI STAYS ACTIVE (VGCCs lock shut)
|
|
# Driver : Ca_micro > 0 between spikes
|
|
# Victim : CDI_factor
|
|
# Outcome : CDI_factor -> 1, effective_conductance -> 0.
|
|
#
|
|
# [CASCADE 6] RESULT - SYNAPSE SILENCES (excitotoxicity protection)
|
|
# Driver : CDI_factor ~= 1
|
|
# Mechanism : effective_conductance = N_VGCC
|
|
# * (1 - eCB_level)
|
|
# * (1 - CDI_factor)
|
|
# * (1 - mGluR*alpha_mGluR)
|
|
# -> raw_influx ~= 0 -> no release.
|
|
# Auto-reset : NCX keeps clearing; CDI recovers once drive stops.
|
|
# =======================================================================
|
|
|
|
import numpy as np
|
|
|
|
# -----------------------------------------------------------------------
|
|
# PARAMETERS
|
|
# -----------------------------------------------------------------------
|
|
|
|
dt = 1.0 # ms - high-freq timestep
|
|
dt_slow = 1000.0 # ms - astrocyte / slow loop timestep
|
|
dt_meta = 60_000.0 # ms - metabolic loop timestep
|
|
|
|
High_Freq_Multiplier = int(dt_slow / dt) # 1000
|
|
Metabolic_Multiplier = int(dt_meta / dt) # 60000
|
|
|
|
# Unit-conversion scalar: biological rate constants in /s, timestep in ms.
|
|
# increment = k [/s] * signal * dt_s [s/step] -> dimensionless per step
|
|
dt_s = dt / 1000.0 # 0.001 s per 1 ms step
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- Voltage / membrane (NKA) --
|
|
# -----------------------------------------------------------------------
|
|
tau_V_pre = 2.0 # NKA: ms - membrane repolarisation time constant after AP
|
|
# The AP waveform decays with this time constant.
|
|
# Controls how long Ca2+ channels see a depolarised membrane.
|
|
V_pre_peak = 1.0 # NKA: normalised peak depolarisation on each AP (dimensionless)
|
|
V_rest = 0.0 # NKA: resting membrane potential (normalised)
|
|
V_pre_voltage = -10.0 # ORIG: driving force used in compute_flux (mV, kept for continuity)
|
|
|
|
# Na/K-ATPase ATP cost
|
|
NKA_cost_per_AP = 0.002 # NKA: ATP units consumed per AP for Na/K-ATPase recharge.
|
|
# This is the largest single ATP cost per spike.
|
|
# [CASCADE 2] Accumulates in ATP_demand each time V_pre = 1;
|
|
# at 20 Hz that is 0.002 * 20 = 0.04 ATP units/s,
|
|
# dominating pump and docking costs at high rates.
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- Calcium influx & buffering --
|
|
# -----------------------------------------------------------------------
|
|
N_VGCC = 100 # ORIG: number of presynaptic VGCCs
|
|
# [CASCADE 6] Absolute ceiling of Ca2+ influx.
|
|
k_flux = 0.05 # ORIG: Ca2+ influx per open channel per unit driving force (a.u.)
|
|
|
|
B_total = 1.0 # ORIG: total buffer capacity (normalised)
|
|
# [CASCADE 4] Saturates during bursting; B_free -> 0.
|
|
tau_buffer_rebind = 200.0 # NEW: ms - buffer recharge time constant
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- Ca2+ clearance rate constants (/ms, already per-millisecond) --
|
|
# -----------------------------------------------------------------------
|
|
# [CASCADE 3] All three define maximum clearance; PMCA and SERCA are ATP-gated.
|
|
k_PMCA = 0.03 # NEW: plasma membrane Ca-ATPase - primary, ATP-dependent
|
|
k_NCX = 0.10 # NEW: sodium-calcium exchanger - fast, NOT ATP-dependent
|
|
k_SERCA = 0.01 # NEW: ER pump - slowest, ATP-dependent
|
|
ATP_half = 0.3 # NEW: Hill half-saturation for ATP-dependent pumps
|
|
|
|
# ATP cost of pumping (per unit Ca2+ cleared per ms)
|
|
ATP_cost_PMCA = 0.0005 # NKA: ATP units per unit Ca2+ extruded by PMCA per ms
|
|
# [CASCADE 2] Second-largest ongoing ATP drain after NKA.
|
|
ATP_cost_SERCA = 0.0002 # NKA: ATP units per unit Ca2+ pumped into ER per ms
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- Deterministic vesicle release (replaces stochastic_release) --
|
|
# -----------------------------------------------------------------------
|
|
# The Ca2+-sensor model (synaptotagmin cooperativity):
|
|
#
|
|
# release_fraction = k_rel * Ca_micro^n / (Ca_micro^n + KD_rel^n)
|
|
# * (1 - NT_suppression_weight * NT_cleft_norm)
|
|
#
|
|
# Term 1 - Hill equation:
|
|
# k_rel : maximum fraction of RRP releasable per spike (saturation ceiling)
|
|
# KD_rel : [Ca2+] at half-maximum release (half-saturation constant)
|
|
# n_rel : Hill cooperativity exponent (~4 for synaptotagmin-1)
|
|
#
|
|
# Term 2 - NT suppression modulation:
|
|
# NT_suppression_weight : how strongly accumulated cleft NT reduces the
|
|
# releasable fraction (autoreceptor-independent brake,
|
|
# represents physical occlusion of release sites and
|
|
# depletion-sensing).
|
|
# NT_cleft_norm : NT_cleft normalised to NT_suppression_sat (0 -> 1).
|
|
#
|
|
# The product is clamped to [0, 1] before multiplying by N_RRP,
|
|
# so released_NT is always a real non-negative number <= N_RRP.
|
|
|
|
k_rel = 0.5 # DET: max releasable fraction of RRP per spike (0->1)
|
|
# Lower = more reluctant synapse at any [Ca2+].
|
|
KD_rel = 1.0 # DET: half-saturation [Ca2+] (same a.u. as Ca_micro)
|
|
# Higher KD = release only at high [Ca2+] peaks.
|
|
n_rel = 4 # DET: Hill exponent - cooperativity of Ca2+ sensor
|
|
# n=4 matches synaptotagmin-1 (four C2 domain sites).
|
|
# Higher n = sharper threshold, more digital release.
|
|
NT_suppression_weight = 0.3 # DET: strength of NT_cleft feedback on release fraction
|
|
# 0 = no suppression; 1 = full block at saturation.
|
|
NT_suppression_sat = 50.0 # DET: NT_cleft level that saturates the suppression term
|
|
# (same units as NT_cleft vesicle count)
|
|
ATP_cost_docking = 0.001 # NKA: ATP units per vesicle docked (RP -> RRP per step)
|
|
# [CASCADE 2] Docking is ATP-dependent (NSF/SNAPs);
|
|
# adds to ATP_demand proportionally to refill_amount.
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- CDI (calcium-dependent inactivation) --
|
|
# -----------------------------------------------------------------------
|
|
k_CDI_rise = 0.8 # ORIG: CDI build rate (/s per unit Ca_micro)
|
|
Ca_micro_saturation = 2.0 # NEW: normalisation ceiling for CDI recovery
|
|
k_CDI_rec = 0.015 # NEW: CDI de-inactivation rate (/s)
|
|
# Both expressed in /s; applied with * dt_s in the loop.
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- Vesicle pools --
|
|
# -----------------------------------------------------------------------
|
|
Max_RRP = 20 # ORIG: [CASCADE 1] ceiling of the firing-ready pool
|
|
Max_RP = 200 # ORIG: [CASCADE 1] ceiling of the reserve pool
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- Trace integrator (RP->RRP recruitment speed) --
|
|
# -----------------------------------------------------------------------
|
|
tau_Tr_Ca = 1000.0 # ORIG: ms - calcium trace decay
|
|
T_high = 0.6 # ORIG: trace threshold -> fast recruitment
|
|
T_low = 0.2 # ORIG: trace threshold -> slow recruitment
|
|
k_rec_fast = 5.0 # /s — fast RP->RRP recruitment rate
|
|
# At high Tr_Ca: refills ~5% of headroom per second.
|
|
# Full RRP recovery from empty takes ~4-5 s of sustained activity.
|
|
k_rec_slow = 0.5 # /s — slow RP->RRP recruitment rate
|
|
# At low Tr_Ca: refills ~0.5% of headroom per second.
|
|
# Matches the ~30-60 s recovery seen after deep depletion.
|
|
|
|
dt_slow_s = dt_slow / 1000.0 # 1.0 s — the slow loop timestep in seconds
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- Postsynaptic --
|
|
# -----------------------------------------------------------------------
|
|
tau_membrane = 20.0 # ORIG: ms
|
|
tau_desensitization = 500.0 # ORIG: ms
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- eCB retrograde brake --
|
|
# -----------------------------------------------------------------------
|
|
# [CASCADE 6] Slow retrograde suppressor of effective_conductance.
|
|
tau_eCB_rise = 2000.0 # ORIG: ms
|
|
tau_eCB_decay = 10_000.0 # ORIG: ms
|
|
eCB_threshold = 0.7 # ORIG: V_post level that triggers eCB synthesis
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- mGluR presynaptic autoreceptor --
|
|
# -----------------------------------------------------------------------
|
|
# [CASCADE 6] Fastest conductance brake; reads NT_cleft directly.
|
|
Km_mGluR = 0.5 # NEW: Michaelis-Menten half-saturation for NT_cleft
|
|
tau_mGluR = 2000.0 # NEW: ms
|
|
alpha_mGluR = 0.4 # NEW: max fractional VGCC suppression
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- Astrocyte / IP3 --
|
|
# -----------------------------------------------------------------------
|
|
tau_IP3 = 3000.0 # ORIG: ms
|
|
IP3_threshold = 0.8 # ORIG
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- Glutamine shuttle --
|
|
# -----------------------------------------------------------------------
|
|
conversion_efficiency_base = 0.8 # ORIG: fraction of Gln pool converted per cycle
|
|
|
|
# -----------------------------------------------------------------------
|
|
# -- NT cleft --
|
|
# -----------------------------------------------------------------------
|
|
tau_NT_decay = 5.0 # ms - passive NT diffusion / dilution out of cleft
|
|
|
|
|
|
# -----------------------------------------------------------------------
|
|
# HELPER FUNCTIONS
|
|
# -----------------------------------------------------------------------
|
|
|
|
def compute_flux(conductance, voltage):
|
|
# ORIG: Ca2+ influx into microdomain.
|
|
# [CASCADE 4] Collapses to near zero once CASCADE 6 locks conductance.
|
|
return k_flux * conductance * abs(voltage)
|
|
|
|
|
|
def deterministic_release(N_RRP, Ca_micro, NT_cleft):
|
|
# DET: Deterministic Ca2+-driven vesicle release.
|
|
#
|
|
# Replaces stochastic_release (binomial with p = p_base * Ca_micro).
|
|
#
|
|
# Biology: synaptotagmin-1 is the fast Ca2+ sensor. Its four C2-domain
|
|
# Ca2+-binding sites give steep, cooperative Ca2+ sensitivity (Hill n~4).
|
|
# Release is not random per vesicle but is driven by the microdomain [Ca2+]
|
|
# that the sensor actually sees. At low [Ca2+] essentially no vesicles fuse;
|
|
# at high [Ca2+] a saturating fraction fuses within the AP window.
|
|
#
|
|
# Step 1: Hill equation gives Ca2+-dependent release fraction (0 -> k_rel).
|
|
Ca_n = Ca_micro ** n_rel
|
|
release_frac = k_rel * Ca_n / (Ca_n + KD_rel ** n_rel)
|
|
#
|
|
# Step 2: NT suppression modulation (0 -> 1, then inverted).
|
|
# Accumulated NT_cleft represents: (a) depletion sensing via presynaptic
|
|
# autoreceptors that are faster than the mGluR loop, and (b) physical
|
|
# competition for release site access at the active zone.
|
|
# This is distinct from the mGluR brake which reduces VGCC conductance;
|
|
# this term reduces the fraction of already-docked vesicles that fuse.
|
|
NT_norm = min(1.0, NT_cleft / NT_suppression_sat)
|
|
suppression = NT_suppression_weight * NT_norm
|
|
release_frac = release_frac * (1.0 - suppression)
|
|
release_frac = np.clip(release_frac, 0.0, 1.0)
|
|
#
|
|
# Step 3: Apply to pool size and floor at zero.
|
|
released = release_frac * N_RRP
|
|
return max(0.0, released)
|
|
|
|
|
|
def map_trace_to_speed(Tr_Ca):
|
|
# ORIG: Maps calcium trace to RP->RRP recruitment rate.
|
|
# [CASCADE 1] k_rec_fast lags release at high firing; collapse is
|
|
# self-accelerating as N_RP and headroom both shrink.
|
|
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 map_calcium_to_inactivation(Ca_micro):
|
|
# ORIG: Ca2+ drives CDI increment each ms.
|
|
# k_CDI_rise is /s; * dt_s gives dimensionless per-step increment.
|
|
# [CASCADE 5] Accumulates between spikes under pump failure.
|
|
return k_CDI_rise * Ca_micro * dt_s
|
|
|
|
|
|
def compute_pump_atp_factor(ATP_level):
|
|
# NEW: Hill function - ATP gates PMCA and SERCA speed.
|
|
# [CASCADE 3] Bridge from ATP_level (Loop 3) to clearance rate (Loop 1).
|
|
# ATP=1.0 -> ~0.92 (near full) | ATP=0.3 -> 0.50 | ATP=0.1 -> ~0.10
|
|
return (ATP_level ** 2) / (ATP_level ** 2 + ATP_half ** 2)
|
|
|
|
|
|
def compute_EPSP(receptor_conductance):
|
|
# ORIG: postsynaptic potential increment.
|
|
return receptor_conductance * 0.1
|
|
|
|
|
|
def compute_postsynaptic_eCB_signal(V_post_history):
|
|
# ORIG: eCB synthesis from sustained postsynaptic activity.
|
|
# [CASCADE 6] Slow retrograde brake; persists 10 s after burst ends.
|
|
recent_mean = (np.mean(V_post_history[-2000:])
|
|
if len(V_post_history) >= 2000
|
|
else np.mean(V_post_history))
|
|
if recent_mean > eCB_threshold:
|
|
return recent_mean - eCB_threshold
|
|
return 0.0
|
|
|
|
|
|
def compute_astrocyte_metabolic_health(Glucose_level, ATP_demand_accumulated):
|
|
# ORIG + NKA: Converts Glucose_level and accumulated ATP demand into:
|
|
# ATP_level -> [CASCADE 2->3 bridge] read every ms in Loop 1
|
|
# conversion_efficiency -> [CASCADE 1] gates glutamine shuttle throughput
|
|
#
|
|
# NKA change: ATP_demand_accumulated (summed in Loop 1 each ms) is now
|
|
# subtracted from health before clamping, so high firing rates visibly
|
|
# reduce ATP_level within the same metabolic window.
|
|
# The demand term is normalised by a scale factor so that a physiological
|
|
# 20 Hz burst causes a ~20% ATP drop over one minute.
|
|
ATP_demand_scale = 50.0 # normalisation: demand at 20 Hz over 60 s ~ 1.0
|
|
health = np.clip(Glucose_level - ATP_demand_accumulated / ATP_demand_scale,
|
|
0.0, 1.0)
|
|
return health, health # (conversion_efficiency, ATP_level)
|
|
|
|
|
|
def trigger_slow_astrocyte_calcium_wave():
|
|
# ORIG: placeholder - would release gliotransmitters over ~10 s.
|
|
pass
|
|
|
|
|
|
# -----------------------------------------------------------------------
|
|
# STATE VARIABLES (initial values)
|
|
# -----------------------------------------------------------------------
|
|
|
|
# -- Voltage / membrane --
|
|
V_pre_state = 0.0 # NKA: continuous membrane potential (0=rest, 1=peak AP)
|
|
# Decays with tau_V_pre after each spike.
|
|
# Controls the effective driving window for Ca2+ influx.
|
|
|
|
# -- Presynaptic Ca2+ --
|
|
Ca_micro = 0.0 # ORIG: free cytosolic [Ca2+] in microdomain [CASCADE 4]
|
|
Ca_ER = 0.5 # NEW: Ca2+ stored in ER
|
|
Ca_buffer_bound = 0.0 # NEW: Ca2+ bound to buffer proteins
|
|
B_free = 1.0 # NEW: free buffer sites [CASCADE 4]
|
|
|
|
# -- CDI --
|
|
CDI_factor = 0.0 # ORIG [CASCADE 5,6]
|
|
|
|
# -- Vesicle pools --
|
|
N_RRP = 15.0 # ORIG: readily-releasable pool [CASCADE 1] (float for deterministic)
|
|
N_RP = 150.0 # ORIG: reserve pool [CASCADE 1]
|
|
|
|
# -- Calcium trace --
|
|
Tr_Ca = 0.0 # ORIG: integrative Ca2+ memory
|
|
|
|
# -- NT in cleft --
|
|
NT_cleft = 0.0 # ORIG [CASCADE 6]
|
|
|
|
# -- Postsynaptic --
|
|
V_post = 0.0
|
|
receptor_conductance = 0.0
|
|
Desensitization_level = 0.0
|
|
V_post_history = []
|
|
|
|
# -- Retrograde / autoreceptor --
|
|
eCB_level = 0.0 # ORIG [CASCADE 6]
|
|
mGluR_activation = 0.0 # NEW [CASCADE 6]
|
|
|
|
# -- Astrocyte --
|
|
IP3 = 0.0
|
|
Glutamine_pool = 50.0
|
|
|
|
# -- Metabolic --
|
|
ATP_level = 1.0 # NEW [CASCADE 2->3]
|
|
ATP_demand = 0.0 # NKA: accumulated ATP demand within current metabolic window
|
|
conversion_efficiency = 0.8 # ORIG
|
|
Glucose_level = 1.0 # ORIG: set < 1.0 to engage metabolic silencing cascade
|
|
|
|
|
|
# -----------------------------------------------------------------------
|
|
# MAIN SIMULATION LOOP
|
|
# -----------------------------------------------------------------------
|
|
|
|
def run_simulation(spike_train, total_steps):
|
|
"""
|
|
spike_train : list of int timestep indices at which an AP arrives
|
|
total_steps : int number of 1 ms steps to simulate
|
|
"""
|
|
|
|
global V_pre_state
|
|
global Ca_micro, Ca_ER, Ca_buffer_bound, B_free
|
|
global CDI_factor
|
|
global N_RRP, N_RP
|
|
global Tr_Ca, NT_cleft
|
|
global V_post, receptor_conductance, Desensitization_level, V_post_history
|
|
global eCB_level, mGluR_activation
|
|
global IP3, Glutamine_pool
|
|
global ATP_level, ATP_demand, conversion_efficiency, Glucose_level
|
|
|
|
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", "eCB_level", "mGluR_activation",
|
|
"released_NT", "ATP_level", "ATP_demand",
|
|
]}
|
|
|
|
spike_set = set(spike_train)
|
|
|
|
for step in range(total_steps):
|
|
|
|
# ==============================================================
|
|
# LOOP 1 - HIGH-FREQUENCY (dt = 1 ms)
|
|
# ==============================================================
|
|
|
|
V_pre = 1 if step in spike_set else 0
|
|
released_NT = 0.0
|
|
|
|
# -- 1A. MEMBRANE VOLTAGE / Na-K-ATPase ------------------------
|
|
#
|
|
# NKA: V_pre_state is now a continuous variable.
|
|
# On each AP it jumps to V_pre_peak, then decays exponentially
|
|
# toward V_rest with time constant tau_V_pre (~2 ms).
|
|
# This models the width of the depolarisation window that keeps
|
|
# VGCCs open after the peak, giving Ca2+ influx a temporal profile
|
|
# rather than a single instantaneous pulse.
|
|
if V_pre == 1:
|
|
V_pre_state = V_pre_peak # AP fires: membrane jumps to peak
|
|
|
|
# Exponential decay toward rest (Na/K-ATPase restores the gradient)
|
|
V_pre_state += (V_rest - V_pre_state) * dt / tau_V_pre
|
|
|
|
# NKA: ATP cost of Na/K-ATPase recharge on each AP.
|
|
# Fired once per spike, not per ms — the cost is per action potential.
|
|
# [CASCADE 2] This is the dominant ATP drain at high firing rates.
|
|
if V_pre == 1:
|
|
ATP_demand += NKA_cost_per_AP
|
|
|
|
# -- 1B. PRESYNAPTIC Ca2+ INFLUX & VESICLE RELEASE -------------
|
|
#
|
|
# Ca2+ influx now uses V_pre_state (the continuous voltage) instead
|
|
# of the binary V_pre flag, so influx has the same temporal profile
|
|
# as the depolarisation window and tapers as the membrane repolarises.
|
|
# [CASCADE 6] OUTCOME: effective_conductance collapses when any of
|
|
# the three suppression terms approaches 1.
|
|
effective_conductance = (
|
|
N_VGCC
|
|
* (1.0 - eCB_level) # [CASCADE 6] retrograde brake
|
|
* (1.0 - CDI_factor) # [CASCADE 5->6] CDI lock-out
|
|
* (1.0 - mGluR_activation * alpha_mGluR) # [CASCADE 6] autoreceptor brake
|
|
)
|
|
|
|
# Ca2+ influx is gated by V_pre_state: significant only while the
|
|
# membrane is depolarised; tapers to zero as V_pre_state -> V_rest.
|
|
raw_influx = compute_flux(effective_conductance, V_pre_state)
|
|
|
|
# NEW: Buffer proteins capture a fraction of raw_influx immediately.
|
|
# [CASCADE 4] B_free -> 0 during sustained bursting (saturation).
|
|
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)
|
|
|
|
# DET: Deterministic Ca2+-driven vesicle release.
|
|
# Released on every ms that Ca_micro > 0 (not only on the spike flag),
|
|
# giving a smooth release profile that follows the Ca2+ transient.
|
|
# This is more physically accurate than gating release on V_pre == 1:
|
|
# in biology, vesicles fuse throughout the Ca2+ microdomain lifetime
|
|
# (~1-2 ms), not only at the exact moment of depolarisation.
|
|
# [CASCADE 1] released_NT draws from N_RRP; the Ca2+ and NT_cleft
|
|
# dependence means release self-limits as both Ca_micro
|
|
# falls (clearance) and NT_cleft rises (suppression).
|
|
if N_RRP > 0 and Ca_micro > 0:
|
|
released_NT = deterministic_release(N_RRP, Ca_micro, NT_cleft)
|
|
released_NT = min(released_NT, N_RP)
|
|
N_RRP -= released_NT
|
|
NT_cleft += released_NT
|
|
NT_released_this_window += released_NT # NEW: accumulator for Loop 2
|
|
|
|
# Passive diffusion — this is fine as-is, represents lateral escape
|
|
NT_cleft *= (1.0 - dt / tau_NT_decay)
|
|
NT_cleft = max(0.0, NT_cleft)
|
|
|
|
# -- 1C. Ca2+ CLEARANCE (every ms) ----------------------------
|
|
#
|
|
# [CASCADE 3] pump_scale: ATP bridge from Loop 3 to clearance.
|
|
pump_scale = compute_pump_atp_factor(ATP_level)
|
|
|
|
# [CASCADE 3->4] PMCA: primary, ATP-dependent.
|
|
cleared_PMCA = k_PMCA * Ca_micro * pump_scale
|
|
# [CASCADE 3 note] NCX: fast, NOT ATP-dependent. Floor, not rescue.
|
|
cleared_NCX = k_NCX * Ca_micro
|
|
# [CASCADE 3->4] SERCA: slowest, ATP-dependent.
|
|
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 # ER store loaded while ATP is available
|
|
|
|
# NKA: ATP cost of PMCA and SERCA clearing Ca2+ this step.
|
|
# [CASCADE 2] Ongoing drain proportional to Ca2+ load; highest
|
|
# during the early burst when Ca_micro peaks.
|
|
ATP_demand += ATP_cost_PMCA * cleared_PMCA
|
|
ATP_demand += ATP_cost_SERCA * cleared_SERCA
|
|
|
|
# NEW: Buffer recharge - bound Ca2+ slowly re-releases to cytosol.
|
|
# [CASCADE 4] Sustains Ca_micro elevation during pump failure.
|
|
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
|
|
|
|
# -- 1D. CDI - RISE AND RECOVERY --------------------------------
|
|
#
|
|
# [CASCADE 5] Rise: proportional to Ca_micro, fires every ms.
|
|
CDI_factor += map_calcium_to_inactivation(Ca_micro)
|
|
# [CASCADE 5] Recovery: suppressed when Ca_micro is high.
|
|
# Self-locking: pump failure -> Ca_micro high ->
|
|
# recovery_rate -> 0 -> CDI_factor -> 1 -> silence.
|
|
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)
|
|
|
|
# -- 1E. TRACE INTEGRATOR -------------------------------------
|
|
# ORIG: integrates Ca_micro; drives RP->RRP recruitment speed.
|
|
Tr_Ca = Tr_Ca + (Ca_micro - Tr_Ca / tau_Tr_Ca) * dt
|
|
|
|
# -- 1G. POSTSYNAPTIC FAST RESPONSE ---------------------------
|
|
# ORIG: receptor activation and desensitization.
|
|
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)
|
|
|
|
# ORIG: NT diffuses / dilutes out of cleft each ms.
|
|
NT_cleft *= (1.0 - dt / tau_NT_decay)
|
|
NT_cleft = max(0.0, NT_cleft)
|
|
|
|
V_post_history.append(V_post)
|
|
if len(V_post_history) > 5000:
|
|
V_post_history.pop(0)
|
|
|
|
# -- 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["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)
|
|
|
|
# ==============================================================
|
|
# LOOP 2 - SLOW / ASTROCYTE (dt_slow = 1 s)
|
|
# ==============================================================
|
|
|
|
if (step % High_Freq_Multiplier) == 0:
|
|
|
|
# EAAT clearance: a fraction of what remains in the cleft
|
|
# (this is now meaningful because diffusion hasn't zeroed it yet
|
|
# at realistic tau_NT_decay values — but we also fix IP3 sourcing)
|
|
cleared_NT = NT_cleft * 0.3
|
|
NT_cleft = max(0.0, NT_cleft - cleared_NT)
|
|
|
|
# IP3 integrates total release load, not post-diffusion residual.
|
|
# This correctly represents the astrocyte sensing cumulative activity.
|
|
IP3 += NT_released_this_window - (IP3 / tau_IP3) * dt_slow
|
|
IP3 = max(0.0, IP3)
|
|
NT_released_this_window = 0.0 # reset accumulator for next window
|
|
|
|
if IP3 > IP3_threshold:
|
|
trigger_slow_astrocyte_calcium_wave()
|
|
|
|
# ORIG: eCB retrograde signal.
|
|
# [CASCADE 6] Second conductance brake (~2 s onset, 10 s decay).
|
|
eCB_signal = compute_postsynaptic_eCB_signal(V_post_history)
|
|
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)
|
|
|
|
|
|
|
|
|
|
# ── LOOP 2 — SLOW / ASTROCYTE (dt_slow = 1 s) ──────────────────────
|
|
|
|
# -- RP -> RRP RECRUITMENT (with pool guards) --------------------
|
|
# Runs once per second. k_rec_fast and k_rec_slow are in /s,
|
|
# so multiplying by dt_slow_s = 1.0 s gives a dimensionless
|
|
# per-step fraction — no hidden unit scaling needed.
|
|
#
|
|
# [CASCADE 1] Recruitment is the only counter-force to depletion.
|
|
# Even k_rec_fast fills only ~5% of headroom per second,
|
|
# lagging well behind release at high firing rates.
|
|
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 cost of docking moves here too — it is per recruitment event,
|
|
# not per millisecond.
|
|
ATP_demand += ATP_cost_docking * refill_amount
|
|
|
|
# ==============================================================
|
|
# LOOP 3 - METABOLIC (dt_meta = 1 min)
|
|
# ==============================================================
|
|
|
|
if (step % Metabolic_Multiplier) == 0:
|
|
|
|
# NKA + ORIG: ATP_level is now driven by both Glucose_level
|
|
# (supply) and ATP_demand (demand accumulated over this window).
|
|
# [CASCADE 2] ATP_demand encodes all three cost sources:
|
|
# NKA_cost_per_AP - dominant at high firing rates
|
|
# ATP_cost_PMCA/SERCA - proportional to Ca2+ load
|
|
# ATP_cost_docking - proportional to recruitment rate
|
|
conversion_efficiency, ATP_level = compute_astrocyte_metabolic_health(
|
|
Glucose_level, ATP_demand
|
|
)
|
|
|
|
# Reset demand accumulator for the next metabolic window.
|
|
ATP_demand = 0.0
|
|
|
|
# [CASCADE 1 - slow refill] Glutamine shuttle rebuilds N_RP.
|
|
refill_RP = Glutamine_pool * conversion_efficiency
|
|
N_RP = min(Max_RP, N_RP + refill_RP)
|
|
Glutamine_pool = max(0.0, Glutamine_pool - refill_RP)
|
|
|
|
return log
|
|
|
|
|
|
# -----------------------------------------------------------------------
|
|
# EXAMPLE USAGE
|
|
# -----------------------------------------------------------------------
|
|
|
|
if __name__ == "__main__":
|
|
import matplotlib.pyplot as plt
|
|
|
|
total_steps = 10_000 # 10 seconds of simulated time
|
|
|
|
# 20 Hz burst for 2 s, then silence.
|
|
# Set Glucose_level = 0.2 to engage the full metabolic cascade.
|
|
spike_train = list(range(0, 2000, 50))
|
|
|
|
results = run_simulation(spike_train, total_steps)
|
|
|
|
t = np.arange(total_steps) * dt
|
|
|
|
fig, axes = plt.subplots(7, 1, figsize=(12, 16), sharex=True)
|
|
fig.suptitle("Tripartite Synapse - Presynaptic Model", fontsize=13)
|
|
|
|
axes[0].plot(t, results["V_pre_state"], color="slateblue", lw=0.8)
|
|
axes[0].set_ylabel("V_pre (a.u.)")
|
|
axes[0].set_title("NKA - membrane repolarises with tau_V_pre after each AP",
|
|
fontsize=9, loc="left")
|
|
|
|
axes[1].plot(t, results["Ca_micro"], color="darkorange", lw=0.8)
|
|
axes[1].set_ylabel("[Ca2+] free (a.u.)")
|
|
axes[1].set_title("CASCADE 4 - residual Ca2+ profile follows V_pre decay",
|
|
fontsize=9, loc="left")
|
|
|
|
axes[2].plot(t, results["CDI_factor"], color="firebrick", lw=0.8, label="CDI factor")
|
|
axes[2].plot(t, results["B_free"], color="steelblue", lw=0.8, label="Buffer free")
|
|
axes[2].set_ylabel("CDI / Buffer (0-1)")
|
|
axes[2].set_title("CASCADE 5 - CDI lock-out | CASCADE 4 - buffer saturation",
|
|
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 - deterministic depletion (Hill Ca2+ sensor + NT suppression)",
|
|
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 multiplicative brakes on effective_conductance",
|
|
fontsize=9, loc="left")
|
|
axes[4].legend(fontsize=8)
|
|
|
|
axes[5].plot(t, results["V_post"], color="navy", lw=0.8)
|
|
axes[5].set_ylabel("V_post (a.u.)")
|
|
axes[5].set_title("CASCADE 6 result - postsynaptic silence", fontsize=9, loc="left")
|
|
|
|
axes[6].plot(t, results["ATP_level"], color="goldenrod", lw=0.8, label="ATP level")
|
|
axes[6].plot(t, results["ATP_demand"], color="tomato", lw=0.8, label="ATP demand (cumul.)")
|
|
axes[6].set_ylabel("ATP (a.u.)")
|
|
axes[6].set_title("CASCADE 2 - NKA + pump + docking demand drives ATP depletion",
|
|
fontsize=9, loc="left")
|
|
axes[6].set_xlabel("Time (ms)")
|
|
axes[6].legend(fontsize=8)
|
|
|
|
plt.tight_layout()
|
|
plt.savefig("./synapse_simulation.png", dpi=150)
|
|
plt.close()
|
|
print("Done.") |