diff --git a/neuron/appunti/traditional_python_models/synapse_simulation.png b/neuron/appunti/traditional_python_models/synapse_simulation.png new file mode 100644 index 0000000..75385c4 Binary files /dev/null and b/neuron/appunti/traditional_python_models/synapse_simulation.png differ diff --git a/neuron/appunti/traditional_python_models/tripartite-new-1.py b/neuron/appunti/traditional_python_models/tripartite-new-1.py new file mode 100644 index 0000000..254fe55 --- /dev/null +++ b/neuron/appunti/traditional_python_models/tripartite-new-1.py @@ -0,0 +1,714 @@ +# 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.") \ No newline at end of file diff --git a/neuron/appunti/traditional_python_models/tripartite.py b/neuron/appunti/traditional_python_models/tripartite.py new file mode 100644 index 0000000..a618d33 --- /dev/null +++ b/neuron/appunti/traditional_python_models/tripartite.py @@ -0,0 +1,493 @@ +import numpy as np +import matplotlib as plt +import matplotlib.pyplot as plt # <-- this is essential +plt.rcdefaults() # Restore default rc parameters + +# ============================================================================ +# TRIPARTITE SYNAPSE MODEL +# ============================================================================ +# This code simulates a single synapse with presynapse, postsynapse, and astrocyte. +# Timescales range from milliseconds (AP, Ca2+ influx) to minutes (vesicle replenishment, +# astrocyte Ca2+ waves, synaptic weight changes). +# +# Units: +# - Time: seconds (s) or milliseconds (ms) as noted. Internal time step is in seconds. +# - Concentrations: micromolar (µM) +# - Vesicle pools: number of vesicles (can be fractional) +# - Conductances: nanosiemens (nS) +# - Membrane potential: millivolts (mV) +# ============================================================================ + +# ---------------------------------------------------------------------------- +# Presynapse Class +# ---------------------------------------------------------------------------- +class Presynapse: + """ + Models the presynaptic terminal. + State variables: + Cf (float): fast calcium concentration (µM) + Cs (float): slow calcium concentration (µM) + R (float): vesicles in Readily Releasable Pool (RRP) + P (float): vesicles in Reserve Pool (RP) + N_VGCC (int): number of functional VGCCs (can be modulated by astrocyte) + release_prob (float): probability of vesicle fusion per AP (computed) + NT_released (float): amount of neurotransmitter ejected (arbitrary units) + """ + def __init__(self, params): + self.p = params # store parameters + # Initial states + self.Cf = 0.0 + self.Cs = 0.0 + self.R = self.p['R_max'] # start with full RRP + self.P = self.p['P_max'] # start with full RP + self.N_VGCC = self.p['N_VGCC0'] # baseline VGCC number + self.release_prob = 0.0 + self.NT_released = 0.0 + + def ap_event(self): + """ + Called when an action potential arrives. + Increases calcium (fast and slow components) based on available VGCCs. + Triggers vesicle release and updates pools. + """ + # 1. Calcium influx (instantaneous jump) + self.Cf += self.p['alpha'] * self.N_VGCC + self.Cs += self.p['beta'] * self.N_VGCC # small slow component + + # 2. Compute release probability (Hill function of total calcium) + C_tot = self.Cf + self.Cs + self.release_prob = self.p['p_max'] * (C_tot**self.p['n_Hill']) / \ + (C_tot**self.p['n_Hill'] + self.p['Kd']**self.p['n_Hill']) + + # 3. Vesicles released from RRP + released = self.R * self.release_prob + self.NT_released = self.p['gamma'] * released + self.R -= released + + # Ensure pools don't go negative (shouldn't happen if dt small) + self.R = max(self.R, 0) + + def step(self, dt): + """ + Update presynaptic states between APs. + - Calcium decay (fast and slow) + - Vesicle mobilization (RP -> RRP) depending on calcium trace + - RP replenishment from an infinite source + """ + # Total calcium for trace classification + C_tot = self.Cf + self.Cs + + # 1. Calcium decay (first order) + self.Cf -= dt * (self.Cf / self.p['tau_f']) + self.Cs -= dt * (self.Cs / self.p['tau_s']) + + # 2. Vesicle mobilization rate based on calcium trace + if C_tot < self.p['theta_low']: + k_mob = self.p['k_slow'] + elif C_tot < self.p['theta_high']: + k_mob = self.p['k_med'] + else: + k_mob = self.p['k_fast'] + + # Move vesicles from RP to RRP + move = k_mob * self.P * dt + self.R += move + self.P -= move + + # 3. RP replenishment (slow refilling) + replenish = self.p['k_rep'] * (self.p['P_max'] - self.P) * dt + self.P += replenish + + # Optional: keep pools within bounds + self.R = min(self.R, self.p['R_max']) + self.P = min(self.P, self.p['P_max']) + self.P = max(self.P, 0) + + def set_vgcc_modulation(self, factor): + """ + Allow astrocyte to modulate VGCC availability. + factor: multiplier (0..1) applied to baseline N_VGCC. + """ + self.N_VGCC = self.p['N_VGCC0'] * factor + + +# ---------------------------------------------------------------------------- +# Postsynapse Class +# ---------------------------------------------------------------------------- +class Postsynapse: + """ + Models the postsynaptic density and spine. + State variables: + G (float): cleft glutamate concentration (a.u.) + gA (float): AMPA conductance (nS) + gN (float): NMDA conductance (nS) + Vm (float): membrane potential (mV) + C_post (float): postsynaptic calcium (µM) + w (float): synaptic weight (dimensionless, modulates AMPAR conductance) + """ + def __init__(self, params): + self.p = params + self.G = 0.0 + self.gA = 0.0 + self.gN = 0.0 + self.Vm = self.p['E_rest'] + self.C_post = self.p['C_post_rest'] + self.w = 1.0 # initial weight + + def receive_nt(self, NT_released): + """ + Neurotransmitter released into cleft. Instantaneous rise, + then decay handled in step(). + """ + self.G += NT_released # instantaneous jump + + def step(self, dt, astrocyte_dserine_factor=1.0): + """ + Update postsynaptic states. + - Glutamate clearance + - AMPA/NMDA conductance kinetics + - Membrane potential (Euler) + - Postsynaptic calcium dynamics + - Synaptic weight plasticity (slow) + """ + # 1. Glutamate clearance (first order) + self.G -= dt * (self.G / self.p['tau_glu']) + + # 2. AMPA receptor kinetics + # Target conductance = max conductance * (G/(G+K)) * weight + target_A = self.p['gA_max'] * (self.G / (self.G + self.p['K_glu'])) * self.w + self.gA += dt * ((target_A - self.gA) / self.p['tau_AMPA']) + + # 3. NMDA receptor kinetics (with voltage-dependent Mg block) + Mg_block = 1.0 / (1.0 + (self.p['Mg'] / 3.57) * np.exp(-0.062 * self.Vm)) + # astrocyte D-serine enhances NMDA conductance (multiplying max) + target_N = self.p['gN_max'] * astrocyte_dserine_factor * \ + (self.G / (self.G + self.p['K_glu'])) * Mg_block + self.gN += dt * ((target_N - self.gN) / self.p['tau_NMDA']) + + # 4. Membrane potential (current balance) + I_syn = self.gA * (self.Vm - self.p['E_AMPA']) + self.gN * (self.Vm - self.p['E_NMDA']) + I_leak = self.p['gL'] * (self.Vm - self.p['E_rest']) + dVm = (-I_leak - I_syn) / self.p['Cm'] + self.Vm += dt * dVm + + # 5. Postsynaptic calcium influx (via NMDA and VGCCs) + # Simplified: calcium current proportional to NMDA conductance and voltage driving force + I_Ca_NMDA = self.p['eta_N'] * self.gN * (self.Vm - self.p['E_Ca']) + # VGCC contribution (activated by depolarisation) + vgcc_act = 1.0 / (1.0 + np.exp(-(self.Vm - self.p['V_half'])/self.p['k_vgcc'])) + I_Ca_VGCC = self.p['eta_V'] * vgcc_act * (self.Vm - self.p['E_Ca']) + self.C_post += dt * (I_Ca_NMDA + I_Ca_VGCC - (self.C_post - self.p['C_post_rest'])/self.p['tau_Ca_post']) + + # 6. Synaptic weight plasticity (slow, calcium-driven) + # LTP and LTD thresholds based on calcium control hypothesis + if self.C_post > self.p['theta_LTP']: + dw = self.p['gamma_LTP'] * (1.0 - self.w) * dt + elif self.C_post > self.p['theta_LTD']: + dw = -self.p['gamma_LTD'] * self.w * dt + else: + dw = 0.0 + self.w += dw + self.w = np.clip(self.w, 0.0, 2.0) # keep within bounds + + +# ---------------------------------------------------------------------------- +# Astrocyte Class +# ---------------------------------------------------------------------------- +class Astrocyte: + """ + Models an astrocytic process enwrapping the synapse. + State variables: + I (float): IP3 concentration (µM) + A_Ca (float): astrocyte calcium (µM) + S (float): gliotransmitter (ATP) concentration (a.u.) + uptake_efficiency (float): modulates glutamate clearance (not used directly) + """ + def __init__(self, params): + self.p = params + self.I = 0.0 + self.A_Ca = self.p['A_Ca_rest'] + self.S = 0.0 + self.uptake = 1.0 # can be modulated + + def sense_glutamate(self, G): + """ + Glutamate spillover activates mGluRs, producing IP3. + IP3 production is proportional to low-pass filtered glutamate. + """ + # Low-pass filter of G (simple exponential moving average) + # We'll use a separate state for filtered G for simplicity. + # For now, we directly update IP3 based on G with a delay. + # Here we use a simplified approach: IP3 production rate depends on G. + pass + + def step(self, dt, G): + """ + Update astrocyte states. + - IP3 dynamics (production from mGluR activation, decay) + - Calcium release from internal stores (IP3-dependent) + - Gliotransmitter release when calcium high + """ + # 1. IP3 production (driven by glutamate spillover) and decay + # Use a Hill function for mGluR activation + mGluR_act = (G**self.p['n_mGluR']) / (G**self.p['n_mGluR'] + self.p['K_mGluR']**self.p['n_mGluR']) + prod = self.p['beta_IP3'] * mGluR_act + self.I += dt * (prod - self.I / self.p['tau_IP3']) + + # 2. Astrocyte calcium (simplified Li-Rinzel type) + # IP3 opens ER channels; calcium-induced calcium release not explicitly modeled + # We use a simple equation: increase when IP3 high, decay otherwise. + # More detailed: dA_Ca/dt = r_IP3 * (I^n/(I^n+K_I^n)) * (A_store - A_Ca) - A_Ca/tau_A_Ca + # For simplicity, use a direct dependence: + target = self.p['A_max'] * (self.I**self.p['n_IP3']) / (self.I**self.p['n_IP3'] + self.p['K_IP3']**self.p['n_IP3']) + self.A_Ca += dt * ((target - self.A_Ca) / self.p['tau_A_rise'] - (self.A_Ca - self.p['A_Ca_rest'])/self.p['tau_A_decay']) + + # 3. Gliotransmitter release (ATP) when calcium above threshold + if self.A_Ca > self.p['theta_S']: + self.S = self.p['S_max'] # instantaneous release, can be made graded + else: + self.S = 0.0 + # Gliotransmitter decays + self.S *= np.exp(-dt / self.p['tau_S']) + + def get_vgcc_modulation(self): + """ + Compute factor to modulate presynaptic VGCCs (e.g., via adenosine). + """ + # ATP released by astrocyte is converted to adenosine, which activates A1 receptors. + # We model it as a simple inhibition: factor = 1 - delta * S/(S+K_S) + factor = 1.0 - self.p['delta_A1'] * (self.S / (self.S + self.p['K_A1'])) + return max(factor, 0.0) + + def get_dserine_modulation(self): + """ + D-serine enhances NMDA receptor function. + Return a factor >1 when astrocyte calcium is high. + """ + # Simple relation: factor = 1 + dserine_max * (A_Ca/(A_Ca+K_D)) + factor = 1.0 + self.p['dserine_max'] * (self.A_Ca / (self.A_Ca + self.p['K_D'])) + return factor + + +# ============================================================================ +# Parameter Sets (with timescales and biological references) +# ============================================================================ + +# All times are in seconds unless noted. + +presyn_params = { + # Calcium dynamics + 'alpha': 0.5, # fast Ca influx per VGCC (µM) + 'beta': 0.02, # slow Ca influx per VGCC (µM) + 'tau_f': 0.020, # fast Ca decay time constant (20 ms) + 'tau_s': 1.0, # slow Ca decay time constant (1 s) + # VGCC + 'N_VGCC0': 10, # baseline number of functional VGCCs + # Release + 'p_max': 0.8, # max release probability + 'Kd': 2.0, # Ca half-activation (µM) + 'n_Hill': 4, # Hill coefficient (cooperativity) + 'gamma': 1.0, # conversion factor from vesicles to NT units + # Vesicle pools + 'R_max': 50, # max RRP size (vesicles) + 'P_max': 200, # max RP size (vesicles) + 'k_fast': 2.0, # fast mobilization rate (s^-1) under high Ca + 'k_med': 0.5, # medium mobilization rate (s^-1) + 'k_slow': 0.1, # slow mobilization rate (s^-1) under low Ca + 'k_rep': 0.03, # RP replenishment rate (s^-1) (time constant ~33 s) + # Calcium trace thresholds (µM) + 'theta_low': 1.0, + 'theta_high': 3.0, +} + +postsyn_params = { + # Glutamate dynamics + 'tau_glu': 0.002, # glutamate decay in cleft (2 ms) + 'K_glu': 1.0, # glutamate half-activation for receptors (a.u.) + # AMPA receptors + 'gA_max': 5.0, # max AMPA conductance (nS) + 'tau_AMPA': 0.002, # AMPA rise/decay time (2 ms, simplified) + 'E_AMPA': 0.0, # reversal potential (mV) + # NMDA receptors + 'gN_max': 2.0, # max NMDA conductance (nS) + 'tau_NMDA': 0.050, # NMDA decay time (50 ms) + 'E_NMDA': 0.0, # reversal potential (mV) + 'Mg': 1.0, # extracellular Mg concentration (mM) + # Membrane parameters + 'Cm': 1.0, # membrane capacitance (µF/cm^2, scaled) + 'gL': 0.1, # leak conductance (nS) + 'E_rest': -70.0, # resting potential (mV) + 'E_Ca': 120.0, # calcium reversal potential (mV) + # Postsynaptic calcium + 'eta_N': 0.02, # NMDA calcium influx factor (µM/ms/nA? simplified) + 'eta_V': 0.01, # VGCC calcium influx factor + 'V_half': -20.0, # half-activation for VGCC (mV) + 'k_vgcc': 5.0, # slope factor for VGCC activation + 'tau_Ca_post': 0.050, # calcium decay time (50 ms) + 'C_post_rest': 0.05, # resting calcium (µM) + # Plasticity thresholds + 'theta_LTD': 0.5, # LTD threshold (µM) + 'theta_LTP': 1.5, # LTP threshold (µM) + 'gamma_LTD': 0.001, # LTD rate (s^-1) + 'gamma_LTP': 0.002, # LTP rate (s^-1) +} + +astro_params = { + # IP3 dynamics + 'beta_IP3': 0.1, # IP3 production rate per glutamate (µM/s) + 'tau_IP3': 2.0, # IP3 decay time constant (2 s) + 'K_mGluR': 0.5, # glutamate half-activation for mGluR (a.u.) + 'n_mGluR': 2, # cooperativity for mGluR + # Astrocyte calcium + 'A_max': 2.0, # max calcium release (µM) + 'K_IP3': 0.5, # IP3 half-activation (µM) + 'n_IP3': 2, # cooperativity for IP3 + 'tau_A_rise': 1.0, # rise time for Ca release (1 s) + 'tau_A_decay': 5.0, # decay time for Ca (5 s) + 'A_Ca_rest': 0.1, # resting calcium (µM) + # Gliotransmitter release + 'theta_S': 0.8, # Ca threshold for release (µM) + 'S_max': 1.0, # max gliotransmitter concentration (a.u.) + 'tau_S': 1.0, # gliotransmitter decay time (1 s) + # Feedback parameters + 'delta_A1': 0.5, # maximum inhibition of VGCCs (fraction) + 'K_A1': 0.3, # half-saturation for adenosine inhibition (a.u.) + 'dserine_max': 0.5, # maximum fractional increase in NMDA conductance + 'K_D': 0.5, # half-saturation for D-serine modulation (µM) +} + +# ============================================================================ +# Simulation Setup +# ============================================================================ + +# Create instances +pre = Presynapse(presyn_params) +post = Postsynapse(postsyn_params) +astro = Astrocyte(astro_params) + +# Simulation parameters +dt = 0.0001 # time step = 0.1 ms +T_total = 30.0 # simulate 30 seconds +time = np.arange(0, T_total, dt) +n_steps = len(time) + +# Storage for plotting +store = { + 't': [], + 'Cf': [], 'Cs': [], 'R': [], 'P': [], + 'G': [], 'gA': [], 'gN': [], 'Vm': [], 'C_post': [], 'w': [], + 'I': [], 'A_Ca': [], 'S': [], + 'AP': [] +} + +# Define spike times (example: two bursts) +spike_times = [] +# 5 spikes at 20 Hz starting at 1.0 s +spike_times.extend(np.linspace(1.0, 1.2, 5, endpoint=False)) +# 10 spikes at 50 Hz starting at 5.0 s +spike_times.extend(np.linspace(5.0, 5.18, 10, endpoint=False)) +# 3 spikes at 10 Hz starting at 15.0 s +spike_times.extend(np.linspace(15.0, 15.2, 3, endpoint=False)) +spike_times.sort() +spike_index = 0 +next_spike = spike_times[spike_index] if spike_times else np.inf + +# Main simulation loop +for i, t in enumerate(time): + # Check for AP + AP = 0 + if t >= next_spike: + AP = 1 + spike_index += 1 + next_spike = spike_times[spike_index] if spike_index < len(spike_times) else np.inf + + # --- Presynapse --- + if AP: + pre.ap_event() + pre.step(dt) + + # --- Astrocyte modulation (feedback) --- + # Astrocyte senses glutamate from cleft (post.G) + astro.step(dt, post.G) + vgcc_factor = astro.get_vgcc_modulation() + pre.set_vgcc_modulation(vgcc_factor) + dserine_factor = astro.get_dserine_modulation() + + # --- Postsynapse --- + # Transfer released NT + if AP: + post.receive_nt(pre.NT_released) + post.step(dt, dserine_factor) + + # Store data every 1 ms to save memory + if i % int(0.001/dt) == 0: + store['t'].append(t) + store['Cf'].append(pre.Cf) + store['Cs'].append(pre.Cs) + store['R'].append(pre.R) + store['P'].append(pre.P) + store['G'].append(post.G) + store['gA'].append(post.gA) + store['gN'].append(post.gN) + store['Vm'].append(post.Vm) + store['C_post'].append(post.C_post) + store['w'].append(post.w) + store['I'].append(astro.I) + store['A_Ca'].append(astro.A_Ca) + store['S'].append(astro.S) + store['AP'].append(AP) + +# ============================================================================ +# Plot Results +# ============================================================================ + +fig, axes = plt.subplots(4, 1, figsize=(10, 12), sharex=True) + +t_plot = np.array(store['t']) + +# Presynapse +axes[0].plot(t_plot, store['Cf'], label='Fast Ca', color='C0') +axes[0].plot(t_plot, store['Cs'], label='Slow Ca', color='C1') +axes[0].set_ylabel('Ca (µM)') +axes[0].legend(loc='upper right') +axes[0].set_title('Presynapse') +ax2 = axes[0].twinx() +ax2.plot(t_plot, store['R'], label='RRP', color='C2', linestyle='--') +ax2.plot(t_plot, store['P'], label='RP', color='C3', linestyle='--') +ax2.set_ylabel('Vesicles') +ax2.legend(loc='lower right') + +# Postsynapse conductances and potential +axes[1].plot(t_plot, store['gA'], label='AMPA', color='C4') +axes[1].plot(t_plot, store['gN'], label='NMDA', color='C5') +axes[1].set_ylabel('Conductance (nS)') +axes[1].legend(loc='upper left') +ax2 = axes[1].twinx() +ax2.plot(t_plot, store['Vm'], label='Vm', color='C6', linestyle='-') +ax2.set_ylabel('Vm (mV)') +ax2.legend(loc='upper right') + +# Postsynaptic calcium and weight +axes[2].plot(t_plot, store['C_post'], label='Post Ca', color='C7') +axes[2].set_ylabel('Ca (µM)') +axes[2].legend(loc='upper left') +ax2 = axes[2].twinx() +ax2.plot(t_plot, store['w'], label='Weight w', color='C8', linestyle='--') +ax2.set_ylabel('Synaptic weight') +ax2.legend(loc='upper right') + +# Astrocyte +axes[3].plot(t_plot, store['I'], label='IP3', color='C9') +axes[3].plot(t_plot, store['A_Ca'], label='Astro Ca', color='C10') +axes[3].plot(t_plot, store['S'], label='Gliotrans.', color='C11', linestyle=':') +axes[3].set_ylabel('Conc. (µM / a.u.)') +axes[3].set_xlabel('Time (s)') +axes[3].legend(loc='upper right') +axes[3].set_title('Astrocyte') + +plt.tight_layout() +plt.tight_layout() +# plt.show() # Comment out or replace +plt.savefig('tripartite_simulation.png', dpi=150, bbox_inches='tight') +print("Simulation finished. Plot saved as 'tripartite_simulation.png'.") +# plt.show() \ No newline at end of file diff --git a/neuron/appunti/traditional_python_models/tripartite_new.py b/neuron/appunti/traditional_python_models/tripartite_new.py new file mode 100644 index 0000000..626cf91 --- /dev/null +++ b/neuron/appunti/traditional_python_models/tripartite_new.py @@ -0,0 +1,648 @@ +# Tripartite Synapse - Multi-Scale Computational Model + +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 — derived once, used wherever a biological rate +# constant is expressed in /s but the loop timestep is in ms. +# Multiplying a /s rate by dt_s gives a dimensionless per-step increment: +# increment = k [/s] * Ca * dt_s [s/step] = dimensionless / step +dt_s = dt / 1000.0 # ms -> s (= 0.001 s per 1 ms step) + +# -- Calcium influx & buffering -- +N_VGCC = 100 # ORIG: number of presynaptic VGCCs + # [CASCADE 6] Sets the absolute ceiling of Ca2+ influx. + # effective_conductance can only reduce from here. +V_pre_voltage = -10.0 # ORIG: mV - membrane potential during AP +k_flux = 0.05 # ORIG: Ca2+ influx per open channel (a.u.) + +B_total = 1.0 # ORIG: total buffer capacity (normalised) + # [CASCADE 4] When B_free -> 0 (buffer saturated), + # all raw_influx enters Ca_micro directly, + # accelerating residual Ca2+ build-up. +tau_buffer_rebind = 200.0 # NEW: ms - time for saturated buffer to recharge + # [CASCADE 4] Slow recharge means buffer stays saturated + # during sustained bursting, removing its + # protective role against Ca2+ accumulation. + +# -- Ca2+ clearance rate constants (per ms) -- +# [CASCADE 3] These three constants define maximum clearance capacity. +# PMCA and SERCA are multiplied by pump_scale (ATP-gated); +# NCX runs independently but cannot fully compensate when they fail. +k_PMCA = 0.03 # NEW: plasma membrane Ca-ATPase - primary high-affinity pump + # [CASCADE 3] First pump to fail under low ATP. Its rate + # constant is the largest ATP-dependent term; + # halving pump_scale loses 0.015/ms of clearance. +k_NCX = 0.10 # NEW: sodium-calcium exchanger - fast, NOT ATP-dependent + # [CASCADE 3 note] NCX keeps clearing during pump failure, + # but k_NCX alone cannot prevent Ca2+ + # accumulation when PMCA + SERCA are both + # impaired - it is a floor, not a rescue. +k_SERCA = 0.01 # NEW: ER pump - slowest, ATP-dependent + # [CASCADE 3] Compounds PMCA failure; also stops loading + # Ca_ER, so the ER store cannot buffer spikes. +ATP_half = 0.3 # NEW: Hill half-saturation for ATP-dependent pumps + # [CASCADE 3] At ATP_level = ATP_half pumps run at 50% speed. + # Below this value pump failure is disproportionately + # severe because the Hill curve is steepest here. + +# -- CDI (calcium-dependent inactivation) -- +k_CDI_rise = 0.8 # ORIG: CDI build rate — units /s per unit Ca_micro + # Biological rate constant expressed in per-second units. + # Applied inside map_calcium_to_inactivation via * dt_s. + # [CASCADE 5] Higher value = VGCCs lock shut sooner at any [Ca2+]. +Ca_micro_saturation = 2.0 # NEW: normalisation ceiling for CDI recovery + # [CASCADE 5] When Ca_micro = Ca_micro_saturation, + # CDI_recovery_rate -> 0 completely. +k_CDI_rec = 0.015 # NEW: CDI de-inactivation rate — units /s + # Biological rate constant expressed in per-second units. + # Applied in the loop via * dt_s. + # [CASCADE 5] Lower value = slower recovery = longer silence + # window after pump failure. + +# -- Vesicle pools -- +Max_RRP = 20 # ORIG: maximum readily-releasable pool size + # [CASCADE 1] Smaller Max_RRP -> exhausted in fewer spikes. +Max_RP = 200 # ORIG: maximum reserve pool size + # [CASCADE 1] Once depleted, only the glutamine shuttle + # (Loop 3, minutes timescale) can refill it. +p_release_base = 0.1 # ORIG: base release probability per vesicle per ms + # [CASCADE 1] Multiplied by Ca_micro inside stochastic_release, + # so higher Ca_micro during a burst draws MORE + # vesicles per spike - a positive feedback that + # accelerates RRP depletion. + +# -- Trace integrator -- +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 = 0.005 # ORIG: fast RP->RRP recruitment rate (/ms) + # [CASCADE 1] Even at maximum, recruitment lags release at + # high firing. When N_RP also falls, the term + # (rate * N_RP * headroom) shrinks on both sides, + # accelerating collapse toward N_RRP = 0. +k_rec_slow = 0.0005 # ORIG: slow RP->RRP recruitment rate (/ms) + +# -- Postsynaptic -- +tau_membrane = 20.0 # ORIG: ms - membrane time constant +tau_desensitization = 500.0 # ORIG: ms - receptor desensitization decay + +# -- eCB retrograde brake -- +tau_eCB_rise = 2000.0 # ORIG: ms - rise time constant + # [CASCADE 6] eCB is a slow-onset brake (~2 s rise). + # Provides early partial conductance + # suppression before CDI lock-out is complete. +tau_eCB_decay = 10_000.0 # ORIG: ms - decay time constant + # [CASCADE 6] Slow decay (10 s) means eCB persists as + # a partial brake even as the burst ends, + # prolonging the silence window. +eCB_threshold = 0.7 # ORIG: V_post level that triggers eCB synthesis + # [CASCADE 6] Requires sustained postsynaptic depolarisation; + # acts as an integrating safety threshold. + +# -- mGluR presynaptic autoreceptor -- +Km_mGluR = 0.5 # NEW: Michaelis-Menten half-saturation for NT_cleft + # [CASCADE 6] Lower Km -> autoreceptor saturates at smaller + # NT_cleft -> faster VGCC suppression onset. +tau_mGluR = 2000.0 # NEW: ms - mGluR activation / desensitization time constant + # [CASCADE 6] Governs how quickly the autoreceptor engages + # and releases as NT_cleft rises and falls. +alpha_mGluR = 0.4 # NEW: max fractional VGCC suppression by mGluR + # [CASCADE 6] Upper bound of this term. Combined with eCB, + # both brakes together can suppress ~70% of + # conductance before CDI becomes significant. + +# -- 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 + # [CASCADE 1] Sets max RP-refill rate. + # Degraded by metabolic fatigue via + # conversion_efficiency in Loop 3. + + +# ----------------------------------------------------------------------- +# HELPER FUNCTIONS +# ----------------------------------------------------------------------- + +def compute_flux(conductance, voltage): + # ORIG: Ca2+ influx into microdomain. + # [CASCADE 4] Raw influx that, after buffer capture, loads Ca_micro. + # Collapses to near zero once CASCADE 6 locks effective_conductance. + return k_flux * conductance * abs(voltage) + + +def stochastic_release(N_RRP, Ca_micro): + # ORIG: Binomial vesicle release from RRP. + # [CASCADE 1] p scales with Ca_micro: higher Ca_micro during the early burst + # draws MORE vesicles per spike than at rest - a positive feedback + # that accelerates RRP depletion beyond what recruitment can match. + p = min(1.0, p_release_base * Ca_micro) + return int(np.random.binomial(int(N_RRP), p)) + + +def map_trace_to_speed(Tr_Ca): + # ORIG: Maps calcium trace to RP->RRP recruitment rate. + # [CASCADE 1] Even at maximum (k_rec_fast = 0.005/ms), recruitment lags + # release at high firing. As N_RP and headroom both fall, + # refill_amount shrinks and N_RRP collapses. + 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 in /s. Multiplying by dt_s converts it to a + # dimensionless per-step increment: + # increment = k_CDI_rise [/s] * Ca_micro * dt_s [s/step] + # Using dt_s here (not dt) makes the unit conversion explicit and + # avoids the implicit /1000 anti-pattern. + # [CASCADE 5] Fires every ms that Ca_micro > 0, including inter-spike + # intervals under pump failure. Result: CDI_factor accumulates + # across spikes instead of resetting, climbing toward 1. + 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] Mechanistic bridge from CASCADE 2 (ATP_level, written by + # Loop 3 once per minute) to CASCADE 4 (Ca_micro staying elevated). + # ATP_level = 1.0 -> pump_scale ~= 0.92 (near full speed) + # ATP_level = 0.3 -> pump_scale = 0.50 (half speed) + # ATP_level = 0.1 -> pump_scale ~= 0.10 (severe failure) + # Hill exponent of 2 gives a sigmoidal response: small drops + # have modest effect; drops below ATP_half are disproportionate. + 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 triggered by sustained high postsynaptic activity. + # [CASCADE 6] Integrates V_post over a 2 s window. Acts as a slow retrograde + # warning signal; its tau_eCB_decay of 10 s means it outlasts + # the burst and keeps conductance partially suppressed during recovery. + 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, step): + # ORIG + NEW ATP output. + # [CASCADE 2] Glucose_level is the root upstream variable of the slow cascade. + # Under sustained high firing, demand outpaces vascular delivery + # and Glucose_level falls. This function converts that into: + # + # ATP_level -> [CASCADE 2->3 bridge] + # Read every ms by compute_pump_atp_factor() in Loop 1. + # This is how minute-scale metabolic state reaches the + # millisecond Ca2+ clearance equations. + # + # conversion_efficiency -> [CASCADE 1, slow refill arm] + # Gates glutamine shuttle throughput; low efficiency starves + # N_RP refill, prolonging vesicle depletion. + # + # Both outputs degrade together as Glucose_level -> 0, + # tightening CASCADE 1 and CASCADE 3-5 simultaneously. + health = np.clip(Glucose_level, 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) +# ----------------------------------------------------------------------- + +# -- Presynaptic Ca2+ -- +Ca_micro = 0.0 # ORIG: free cytosolic [Ca2+] in microdomain + # [CASCADE 4] Central victim variable. + # Normally returns to ~0 between spikes; + # stays elevated when pumps fail. +Ca_ER = 0.5 # NEW: Ca2+ stored in ER (loaded by SERCA) +Ca_buffer_bound = 0.0 # NEW: Ca2+ currently held by buffer proteins +B_free = 1.0 # NEW: free buffer binding sites (= B_total at rest) + # [CASCADE 4] B_free -> 0 during sustained bursting. + # Once saturated, capture_fraction -> 0 and + # raw_influx enters Ca_micro unattenuated. + +# -- CDI -- +CDI_factor = 0.0 # ORIG: fractional VGCC inactivation (0 -> 1) + # [CASCADE 5] Primary silence variable; rises with Ca_micro, + # recovery suppressed by Ca_micro. + # [CASCADE 6] Appears as (1 - CDI_factor) in effective_conductance; + # when CDI_factor -> 1 conductance collapses to 0. + +# -- Vesicle pools -- +N_RRP = 15 # ORIG: readily-releasable pool (docked vesicles) + # [CASCADE 1] Depleted in seconds at high firing rate. +N_RP = 150 # ORIG: reserve pool + # [CASCADE 1] Depleted over tens of seconds; only the glutamine + # shuttle (Loop 3) can refill it on minute timescales. + +# -- Calcium trace -- +Tr_Ca = 0.0 # ORIG: integrative calcium memory; drives RP->RRP recruitment speed. + +# -- NT in cleft -- +NT_cleft = 0.0 # ORIG: normalised NT concentration + # [CASCADE 6] Drives mGluR_activation (fastest conductance brake). + +# -- Postsynaptic -- +V_post = 0.0 # ORIG +receptor_conductance = 0.0 # ORIG +Desensitization_level = 0.0 # ORIG +V_post_history = [] # ORIG: rolling buffer for eCB computation + +# -- Retrograde / autoreceptor -- +eCB_level = 0.0 # ORIG: endocannabinoid retrograde brake (0 -> 1) + # [CASCADE 6] Second multiplicative suppressor of + # effective_conductance. Slow onset (~2 s), + # slow decay (~10 s). +mGluR_activation = 0.0 # NEW: presynaptic autoreceptor occupancy (0 -> 1) + # [CASCADE 6] Third multiplicative suppressor. Fastest + # onset because it reads NT_cleft directly, + # no postsynaptic relay required. + +# -- Astrocyte -- +IP3 = 0.0 # ORIG +Glutamine_pool = 50.0 # ORIG + +# -- Metabolic -- +ATP_level = 1.0 # NEW: normalised ATP (1 = replete) + # [CASCADE 2->3] Bridge variable: written by Loop 3 + # (minutes), read every ms in Loop 1 + # via compute_pump_atp_factor(). +conversion_efficiency = 0.8 # ORIG +Glucose_level = 1.0 # ORIG: external glucose supply (0 -> 1) + # [CASCADE 2] Root input of the slow cascade arm. + # Set < 1.0 to activate metabolic silencing. + # Try 0.2 to observe full pump failure. + +# -- NT cleft decay -- +tau_NT_decay = 5.0 # ms - passive NT diffusion / dilution out of cleft + + +# ----------------------------------------------------------------------- +# 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 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, conversion_efficiency, Glucose_level + + log = {k: [] for k in [ + "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", + ]} + + 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 + + # -- 1A. PRESYNAPTIC SPIKE HANDLING ---------------------------- + if V_pre == 1: + + # [CASCADE 6] OUTCOME of the full cascade. + # All three suppression terms must be near zero + # for significant Ca2+ influx to occur. + # + # Onset timing during a burst: + # mGluR*alpha_mGluR rises first (~seconds, reads NT_cleft) + # eCB_level rises second (~2-10 s, retrograde) + # CDI_factor locks last but is irreversible until + # Ca_micro falls (needs pump recovery) + 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 + ) + + raw_influx = compute_flux(effective_conductance, V_pre_voltage) + + # NEW: Buffer proteins capture a fraction of raw_influx immediately. + # [CASCADE 4] capture_fraction = B_free / B_total. + # Early in a burst B_free ~= B_total so most influx is + # captured and Ca_micro rises slowly (buffer is protective). + # As bursting continues B_free -> 0 (buffer saturates), + # capture_fraction -> 0 and full raw_influx enters Ca_micro, + # accelerating the residual Ca2+ accumulation of CASCADE 4. + 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) + + # ORIG: NT release from RRP. + # [CASCADE 1] Each spike draws from N_RRP via stochastic_release. + # Draw rate (p * N_RRP) exceeds refill rate + # (k_rec_fast * N_RP * headroom) at high firing, + # driving N_RRP -> 0. + if N_RRP > 0: + released_NT = stochastic_release(N_RRP, Ca_micro) + N_RRP = max(0, N_RRP - released_NT) # floor guard + NT_cleft += released_NT + + # -- 1B. Ca2+ CLEARANCE (runs every ms, spike or not) ---------- + + # [CASCADE 3] pump_scale links CASCADE 2 (ATP depletion) to + # CASCADE 4 (residual Ca2+ accumulation). + # Written from ATP_level which is updated once per minute + # in Loop 3; read here every millisecond. + pump_scale = compute_pump_atp_factor(ATP_level) + + # [CASCADE 3->4] PMCA: primary high-affinity Ca2+ extrusion pump. + # ATP-dependent. First pump to be impaired; largest + # single loss of clearance capacity when pump_scale drops. + cleared_PMCA = k_PMCA * Ca_micro * pump_scale + + # [CASCADE 3 note] NCX: fast, NOT ATP-dependent. + # Continues clearing during metabolic failure. + # Acts as a partial floor but cannot alone prevent + # Ca2+ accumulation when PMCA + SERCA are impaired. + # Also enables the auto-reset: once high-frequency + # drive stops, NCX gradually lowers Ca_micro and + # allows CDI recovery even without restored ATP. + cleared_NCX = k_NCX * Ca_micro + + # [CASCADE 3->4] SERCA: ER pump, ATP-dependent. + # Compounds PMCA failure. Also stops loading Ca_ER, + # so the ER store cannot buffer subsequent spikes. + 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 # SERCA fills the ER store when ATP is available + + # NEW: Buffer recharge - bound Ca2+ slowly re-releases back to cytosol. + # [CASCADE 4] During pump failure, released buffer Ca2+ is not promptly + # cleared, so it adds to Ca_micro and sustains its elevation + # between spikes - a secondary positive feedback within 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 + + # -- 1C. CDI - RISE AND RECOVERY -------------------------------- + + # [CASCADE 5] CDI increment: fires every ms proportional to Ca_micro. + # Under normal conditions Ca_micro returns to ~0 between + # spikes and CDI is balanced by recovery below. + # Under pump failure (CASCADE 4) Ca_micro stays > 0 + # between spikes, so this increment fires continuously + # and CDI_factor climbs monotonically toward 1. + CDI_factor += map_calcium_to_inactivation(Ca_micro) + + # [CASCADE 5] CDI recovery: maximised when Ca_micro = 0, suppressed + # as Ca_micro approaches Ca_micro_saturation. + # This is the self-locking mechanism of the cascade: + # pump failure -> Ca_micro stays high + # Ca_micro high -> CDI_recovery_rate -> 0 + # CDI no recovery -> CDI_factor -> 1 + # CDI_factor -> 1 -> effective_conductance -> 0 (CASCADE 6) + # k_CDI_rec is in /s; multiply by dt_s for a per-step decrement. + 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) + + # -- 1D. TRACE INTEGRATOR -------------------------------------- + # ORIG: integrates Ca_micro; drives RP->RRP recruitment speed. + Tr_Ca = Tr_Ca + (Ca_micro - Tr_Ca / tau_Tr_Ca) * dt + + # -- 1E. RP -> RRP RECRUITMENT (with pool guards) -------------- + # [CASCADE 1] Only counter-force to vesicle depletion. + # refill_amount = rate * N_RP * (Max_RRP - N_RRP) + # As N_RP falls (reservoir empties) and N_RRP -> 0 + # (headroom = Max_RRP, maximised), the N_RP term collapses + # and refill becomes negligible - depletion becomes + # self-sustaining once N_RP is exhausted. + current_recruitment_rate = map_trace_to_speed(Tr_Ca) + refill_amount = current_recruitment_rate * N_RP * (Max_RRP - N_RRP) + refill_amount = max(0.0, refill_amount) # never negative + refill_amount = min(refill_amount, N_RP) # cannot exceed RP supply + + N_RRP = min(N_RRP + refill_amount, Max_RRP) # ceiling guard + N_RP = max(0.0, N_RP - refill_amount) # floor guard + + # -- 1F. 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["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) + + # ============================================================== + # LOOP 2 - SLOW / ASTROCYTE (dt_slow = 1 s) + # ============================================================== + + if (step % High_Freq_Multiplier) == 0: + + # ORIG: Astrocyte clears NT from cleft via EAATs. + cleared_NT = NT_cleft * 0.3 + NT_cleft = max(0.0, NT_cleft - cleared_NT) + + # ORIG: IP3 integrates clearance load. + IP3 += cleared_NT - (IP3 / tau_IP3) * dt_slow + IP3 = max(0.0, IP3) + + if IP3 > IP3_threshold: + trigger_slow_astrocyte_calcium_wave() + + # ORIG: eCB retrograde signal. + # [CASCADE 6] Second conductance brake. Onset ~2 s after sustained + # postsynaptic activity crosses eCB_threshold. Its slow + # decay (10 s) means it persists as a partial brake + # during the recovery window, preventing premature + # re-engagement of the synapse after a silence event. + 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) + + # NEW: mGluR presynaptic autoreceptor feedback. + # [CASCADE 6] Fastest conductance brake. Reads NT_cleft directly + # (no postsynaptic relay), so it rises within seconds + # of burst onset when NT_cleft is highest. + # At saturation it suppresses up to alpha_mGluR = 40% + # of conductance, acting as the first partial brake + # and slightly slowing Ca2+ influx before CDI builds. + # Decays rapidly when NT_cleft falls (e.g. after vesicle + # depletion in CASCADE 1), so it does not contribute to + # the irreversible lock-out - that role belongs to CDI. + mGluR_target = NT_cleft / (NT_cleft + Km_mGluR) + mGluR_activation += (mGluR_target - mGluR_activation) * (dt_slow / tau_mGluR) + mGluR_activation = np.clip(mGluR_activation, 0.0, 1.0) + + # ============================================================== + # LOOP 3 - METABOLIC (dt_meta = 1 min) + # ============================================================== + + if (step % Metabolic_Multiplier) == 0: + + # [CASCADE 2] Root of the slow cascade arm. + # Glucose_level is the external input that sets the + # pace of ATP replenishment. Under sustained demand + # it falls, reducing ATP_level. + # + # ATP_level written here is the CASCADE 2->3 bridge: + # it is read every millisecond in Loop 1 by + # compute_pump_atp_factor(), carrying the minute-scale + # metabolic state into the Ca2+ clearance equations. + # + # conversion_efficiency (also written here) is the + # CASCADE 1 slow-refill arm: low efficiency starves + # N_RP, making vesicle depletion more severe and + # prolonged after each burst episode. + conversion_efficiency, ATP_level = compute_astrocyte_metabolic_health( + Glucose_level, step + ) + + # [CASCADE 1 - slow refill] Glutamine shuttle rebuilds N_RP. + # Gated by conversion_efficiency; when low, N_RP cannot + # recover between bursts, so successive episodes deplete + # faster - a cumulative fatigue effect. + 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 before calling run_simulation() to engage + # the full metabolic cascade (CASCADES 2-5). + spike_train = list(range(0, 2000, 50)) + + results = run_simulation(spike_train, total_steps) + + t = np.arange(total_steps) * dt + + fig, axes = plt.subplots(6, 1, figsize=(12, 14), sharex=True) + fig.suptitle("Tripartite Synapse - Metabolic Silencing Cascade", fontsize=13) + + axes[0].plot(t, results["Ca_micro"], color="darkorange", lw=0.8) + axes[0].set_ylabel("[Ca2+] free (a.u.)") + axes[0].set_title("CASCADE 4 - residual Ca2+ stays high under pump failure", + fontsize=9, loc="left") + + axes[1].plot(t, results["CDI_factor"], color="firebrick", lw=0.8, label="CDI factor") + axes[1].plot(t, results["B_free"], color="steelblue", lw=0.8, label="Buffer free") + axes[1].set_ylabel("CDI / Buffer (0-1)") + axes[1].set_title("CASCADE 5 - CDI lock-out | CASCADE 4 secondary - buffer saturation", + fontsize=9, loc="left") + axes[1].legend(fontsize=8) + + axes[2].plot(t, results["N_RRP"], color="teal", lw=0.8, label="RRP") + axes[2].plot(t, results["N_RP"], color="purple", lw=0.8, label="RP") + axes[2].set_ylabel("Vesicles") + axes[2].set_title("CASCADE 1 - vesicle depletion (fast)", fontsize=9, loc="left") + axes[2].legend(fontsize=8) + + axes[3].plot(t, results["NT_cleft"], color="darkgreen", lw=0.8, label="NT cleft") + axes[3].plot(t, results["mGluR_activation"], color="saddlebrown", lw=0.8, label="mGluR (CASCADE 6)") + axes[3].plot(t, results["eCB_level"], color="crimson", lw=0.8, label="eCB (CASCADE 6)") + axes[3].set_ylabel("Cleft / Feedback") + axes[3].set_title("CASCADE 6 - three multiplicative brakes on effective_conductance", + fontsize=9, loc="left") + axes[3].legend(fontsize=8) + + axes[4].plot(t, results["V_post"], color="navy", lw=0.8) + axes[4].set_ylabel("V_post (a.u.)") + axes[4].set_title("CASCADE 6 result - postsynaptic silence", fontsize=9, loc="left") + + axes[5].plot(t, results["ATP_level"], color="goldenrod", lw=0.8) + axes[5].set_ylabel("ATP level (0-1)") + axes[5].set_title("CASCADE 2 - set Glucose_level < 1.0 to activate this arm", + fontsize=9, loc="left") + axes[5].set_xlabel("Time (ms)") + + plt.tight_layout() + plt.savefig("./synapse_simulation.png", dpi=150) + plt.close() + print("Done.") \ No newline at end of file