# 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.")