diff --git a/neuron/BEH-SOMA.md b/neuron/BEH-SOMA.md index e2ee370..b331654 100644 --- a/neuron/BEH-SOMA.md +++ b/neuron/BEH-SOMA.md @@ -33,13 +33,10 @@ The soma is therefore not a simple threshold device. It is a dynamic integrator In this model we decide to simplify: - We do not model the axon hillock as a separate compartment — threshold crossing is computed directly from V_soma -- We do not model channel kinetics — the AP is treated as an instantaneous threshold event with no rise time or repolarisation dynamics -- We do not model the refractory period — the soma can fire on every ms if input is sufficient - We do not model neuromodulatory inputs — threshold and gain are fixed parameters - We do not model subthreshold oscillations — V_soma is a simple leaky integrator - We do not model the f-I curve explicitly — firing rate emerges from the threshold crossings of V_soma across the simulation - We do not model somatic ATP separately — the soma shares the postsynaptic ATP pool (`ATP_level_post`) drawn from the same astrocyte glucose supply -- Soma firing is driven by an external `soma_spike_train` rather than emerging from V_soma threshold crossings — V_soma is computed for reference but does not itself trigger firing in this version The simplifications imply that: diff --git a/neuron/README.md b/neuron/README.md index 1f05ddd..5599e4d 100644 --- a/neuron/README.md +++ b/neuron/README.md @@ -25,3 +25,62 @@ bAP ## Flusso da SOMA a PRE AP + + +Based on the computational model provided, here is the complete breakdown of all simulated behaviors, categorized by functional compartment. + +## 1. Presynaptic Behaviors + +* **Action Potential Arrival (`V_pre`):** When a spike occurs, the membrane potential (`V_pre_state`) jumps to a peak and decays based on `tau_V_pre`. This profile determines the duration of ion channel opening. +* **Calcium Influx (`VGCC`):** Voltage-Gated Calcium Channels open based on `V_pre_state`. The flow is regulated by three "brakes": **eCB** (retrograde), **CDI** (inactivation), and **mGluR** (autoreceptor). +* **Intracellular Buffering:** Free calcium (`Ca_micro`) is immediately captured by buffers (`B_free`). As activity increases and buffers saturate, the effective calcium concentration rises faster (**Metabolic Cascade 4**). +* **Vesicle Release (NT):** Neurotransmitter release is **deterministic** and follows a Hill equation (simulating Synaptotagmin-1 cooperativity). It is limited by the number of vesicles in the Prontly Releasable Pool (`N_RRP`) and suppressed by high existing levels of NT in the cleft. +* **Vesicle Recycling:** Vesicles move from the Reserve Pool (`N_RP`) to the `N_RRP` at a rate determined by the calcium trace (`Tr_Ca`). Fast recruitment occurs during high activity; slow recruitment occurs at rest. +* **Calcium-Dependent Inactivation (CDI):** Local calcium entering through channels causes them to close (`CDI_factor`). If calcium clearance fails due to low ATP, the CDI "locks" the synapse into a silent state to prevent damage. + +## 2. Postsynaptic Behaviors + +* **AMPA Activation:** Released NT opens AMPA receptors, allowing **Na+** influx. This generates the local excitatory post-synaptic potential (EPSP). +* **Receptor Desensitization:** Continuous exposure to NT reduces the sensitivity of the receptors (`Desensitization_level`), mimicking the presynaptic CDI behavior to prevent over-stimulation. +* **NMDA Coincidence Detection:** NMDA channels open only if **NT is present** AND the **membrane is depolarized** (removing the Mg2+ block). Depolarization is achieved via local AMPA drive plus the back-propagating action potential (**bAP**) from the soma. +* **eCB Synthesis:** When postsynaptic calcium (`Ca_post`) crosses a specific threshold, **Endocannabinoids** are synthesized and sent back to the presynapse to suppress further NT release. + +## 3. Dendritic Behaviors + +* **EPSP Summation:** The dendritic branch (`DB`) acts as a passive integrator. It collects `receptor_conductance` from all active spines and sums them into `V_dend`. +* **Passive Decay:** `V_dend` decays over time according to `tau_dend`, determining the temporal window in which multiple inputs can summate to trigger a somatic spike. +* **bAP Distribution:** When the soma fires, a **back-propagating Action Potential** (`V_bAP`) is broadcasted instantly through the dendrite to all spines to enable NMDA coincidence detection. + +## 4. Somatic Behaviors + +* **Leaky Integration:** The soma integrates the signal from the dendrite (`V_dend`) scaled by `soma_weight`. It acts as a leaky integrator with a time constant of `tau_soma`. +* **Action Potential (AP) Generation:** If `V_soma` crosses the threshold, a multi-phase AP is triggered: + 1. **Rising Phase:** Simulated Na+ channel opening (reaches `V_AP_peak`). + 2. **Falling Phase:** Simulated K+ channel opening (drops to `V_AHP`). + 3. **AHP Phase:** Recovery from hyperpolarization back to rest. +* **Refractory Periods:** After firing, the soma enters an **Absolute Refractory Period** (cannot fire) followed by a **Relative Refractory Period** (threshold is temporarily much higher). + +## 5. Astrocytic Behaviors + +* **Neurotransmitter Clearance:** The astrocyte actively removes NT from the synaptic cleft, governed by the `tau_NT_decay` and metabolic cycles. +* **Glutamine Shuttle:** Cleared NT is converted and recycled back to the presynaptic Reserve Pool (`RP`) with a specific `conversion_efficiency`. +* **IP3 Signaling & Calcium Wave:** Accumulated NT triggers **IP3** production in the astrocyte. If it crosses a threshold, an **astrocytic calcium wave** is triggered. +* **Metabolic Support:** The calcium wave provides a "boost" to the `conversion_efficiency`, helping the synapse recover vesicles faster during high demand. + +## 6. Metabolic & Shared Behaviors (ATP Loop) + +* **ATP Consumption:** Every Action Potential (Pre and Post) and every calcium pumping action (`PMCA`, `SERCA`) drains a shared **Glucose/ATP** budget. +* **Pump Scaling:** The speed of ion pumps is determined by a Hill function of available `ATP_level`. Low energy leads to **Pump Failure**. +* **Metabolic Silencing:** A 6-stage cascade where high firing leads to ATP depletion, which causes pump failure, leading to residual calcium, which triggers CDI, finally **silencing the synapse** to protect against excitotoxicity. + +--- + +### Logic Summary Table + +| Input | Process | Output | +| :--- | :--- | :--- | +| **NT in Cleft** | AMPA / NMDA Opening | **V_post** (Postsynaptic Potential) | +| **V_post** | Dendritic Summation | **V_dend** (Dendritic Potential) | +| **V_dend** | Somatic Integration | **V_soma** (Somatic Potential) | +| **V_soma > Threshold** | Spike Kinetics | **Forward AP** & **Retrograde bAP** | +| **Low ATP** | Pump Failure | **Synaptic Silencing** (Protection) | diff --git a/neuron/appunti/traditional_python_models/synapse_simulation.png b/neuron/appunti/traditional_python_models/synapse_simulation.png deleted file mode 100644 index 974c1a8..0000000 Binary files a/neuron/appunti/traditional_python_models/synapse_simulation.png and /dev/null differ diff --git a/neuron/appunti/traditional_python_models/tripartite-soma.py b/neuron/appunti/traditional_python_models/tripartite-soma.py new file mode 100644 index 0000000..450715e --- /dev/null +++ b/neuron/appunti/traditional_python_models/tripartite-soma.py @@ -0,0 +1,883 @@ +# Tripartite Synapse - Multi-Scale Computational Model +# ===================================================== +# Presynaptic + Postsynaptic perspectives, fully integrated. +# +# Change log: +# ORIG - present from the original document +# NEW - added in the missing-behavior integration pass +# DET - deterministic Ca2+-driven vesicle release +# NKA - explicit Na/K-ATPase V_pre decay and ATP cost +# POST-ATP - postsynaptic Ca2+ dynamics and ATP loop +# DEND - dendritic branch: EPSP summation, V_dend, V_bAP +# SOMA - somatic integration: V_soma, AP threshold, refractory, +# channel kinetics, emergent bAP replacing external bAP_train +# FIX - corrections applied in this pass: +# * NT_released_this_window accumulator (was missing entirely) +# * k_rec_fast / k_rec_slow converted to /s, recruitment moved to Loop 2 +# * dt_slow_s added +# * mGluR now reads NT_released_this_window (not NT_cleft snapshot) +# * IP3 now reads NT_released_this_window (not cleared_NT residual) +# * wave_active flag + conversion_efficiency boost on astrocyte wave +# * CDI rise gated to spike window only +# +# Clock structure: +# Loop 1 - dt = 1 ms (Ca2+, vesicle release, traces, postsynaptic fast) +# Loop 2 - dt = 1000 ms (astrocyte clearance, eCB, mGluR, recruitment) +# Loop 3 - dt = 60000 ms (glutamine shuttle, metabolic health) +# +# ======================================================================= +# THREE CLOSED LOOPS +# ======================================================================= +# +# PRESYNAPTIC: +# NT loop : release (ms) -> cleft -> astrocyte clearance (s) -> +# glutamine shuttle (min) -> RP refill -> RRP -> release +# Ca2+ loop : VGCC influx (ms) -> Tr_Ca -> recruitment speed (s) -> +# eCB retrograde from post (s) -> VGCC suppression +# ATP loop : NKA + pump costs (ms) -> ATP_demand (min) -> ATP_level -> +# pump_scale -> Ca2+ clearance rate -> CDI recovery +# +# POSTSYNAPTIC: +# NT detection loop : NT_cleft -> AMPA -> V_post -> desensitization -> +# reduces next response +# Ca2+ coincidence : NMDA (NT + V_post) -> Ca_post -> eCB -> pre brake +# ATP loop : NKA + PMCA costs (ms) -> ATP_demand_post (min) -> +# ATP_level_post -> pump_scale_post -> Ca_post clearance +# +# SHARED: +# eCB_level : post synthesises -> pre reads (retrograde brake) +# NT_cleft : pre releases -> post detects -> astrocyte clears +# Glucose_level : astrocyte supplies both sides from same budget +# +# ======================================================================= +# METABOLIC SILENCING CASCADE (presynaptic) +# ======================================================================= +# [CASCADE 1] HIGH FIRING -> VESICLE DEPLETION (~seconds) +# release rate >> recruitment rate -> N_RRP -> 0 +# [CASCADE 2] HIGH FIRING -> ATP DEPLETION (~minutes) +# NKA + PMCA + docking demand > glucose-driven supply +# [CASCADE 3] LOW ATP -> PUMP FAILURE +# pump_scale = Hill(ATP_level) -> cleared_PMCA/SERCA fall +# [CASCADE 4] PUMP FAILURE -> RESIDUAL Ca2+ STAYS HIGH +# Ca_micro persists between spikes +# [CASCADE 5] RESIDUAL Ca2+ -> CDI LOCKS VGCCs SHUT +# CDI rise (spike only) + recovery blocked by Ca2+ -> CDI -> 1 +# [CASCADE 6] SYNAPSE SILENCES (excitotoxicity protection) +# effective_conductance = N_VGCC*(1-eCB)*(1-CDI)*(1-mGluR*alpha) +# -> 0; NCX auto-reset when drive stops +# +# POSTSYNAPTIC ATP CASCADE (no CDI equivalent -> dangerous): +# [POST-ATP 1] HIGH V_post + NMDA -> ATP_demand_post rises +# [POST-ATP 2] ATP_level_post falls -> pump_scale_post falls +# [POST-ATP 3] Ca_post clearance slows -> Ca_post stays elevated +# [POST-ATP 4] Ca_post > eCB_threshold without real coincidence +# -> false retrograde signal suppresses presynapse +# [POST-ATP 5] Critically low ATP_post -> runaway Ca_post -> excitotoxicity +# ======================================================================= + +import numpy as np + +# ----------------------------------------------------------------------- +# CLOCK +# ----------------------------------------------------------------------- +dt = 1.0 # ms +dt_slow = 1000.0 # ms +dt_meta = 60_000.0 # ms + +High_Freq_Multiplier = int(dt_slow / dt) # 1000 +Metabolic_Multiplier = int(dt_meta / dt) # 60000 + +dt_s = dt / 1000.0 # 0.001 s/step - for /s rate constants in Loop 1 +dt_slow_s = dt_slow / 1000.0 # 1.0 s/step - for /s rate constants in Loop 2 + +# ----------------------------------------------------------------------- +# PRESYNAPTIC PARAMETERS +# ----------------------------------------------------------------------- + +# -- Voltage / membrane -- +tau_V_pre = 2.0 # ms - AP waveform decay (Na/K-ATPase recharge) +V_pre_peak = 1.0 # a.u. - normalised AP peak +V_rest = 0.0 # a.u. - resting potential +V_pre_voltage = -10.0 # mV - driving force for compute_flux + +NKA_cost_per_AP = 0.002 # ATP units per AP (dominant drain at high rates) + +# -- Ca2+ influx & buffering -- +N_VGCC = 100 # number of VGCCs (ceiling of effective_conductance) +k_flux = 0.05 # Ca2+ influx per open channel per unit driving force +B_total = 1.0 # total buffer capacity (normalised) +tau_buffer_rebind = 200.0 # ms - buffer recharge time constant + +# -- Ca2+ clearance (/ms constants) -- +k_PMCA = 0.03 # ATP-dependent primary pump +k_NCX = 0.10 # ATP-independent floor +k_SERCA = 0.01 # ATP-dependent ER pump +ATP_half = 0.3 # Hill half-saturation for presynaptic pumps + +ATP_cost_PMCA = 0.0005 # ATP per unit Ca2+ extruded by PMCA +ATP_cost_SERCA = 0.0002 # ATP per unit Ca2+ pumped into ER +ATP_cost_docking = 0.001 # ATP per vesicle docked (RP->RRP) + +# -- Deterministic release (Hill + NT suppression) -- +k_rel = 0.5 # max releasable fraction of RRP per spike +KD_rel = 1.0 # half-saturation [Ca2+] +n_rel = 4 # Hill cooperativity (synaptotagmin-1) +NT_suppression_weight = 0.3 # max NT_cleft brake on release fraction +NT_suppression_sat = 50.0 # NT_cleft level that saturates suppression + +# -- CDI -- +k_CDI_rise = 0.8 # /s - CDI build rate (applied * dt_s, spike only) +Ca_micro_saturation = 2.0 # normalisation ceiling for CDI recovery +k_CDI_rec = 0.015 # /s - CDI de-inactivation rate (applied * dt_s) + +# -- Vesicle pools -- +Max_RRP = 20 +Max_RP = 200 + +# -- Calcium trace -- +tau_Tr_Ca = 1000.0 # ms +T_high = 0.6 # Tr_Ca threshold -> fast recruitment +T_low = 0.2 # Tr_Ca threshold -> slow recruitment + +# -- RP->RRP recruitment (/s, runs in Loop 2) -- +k_rec_fast = 5.0 # /s - fast recruitment (at Tr_Ca > T_high) +k_rec_slow = 0.5 # /s - slow recruitment (at Tr_Ca < T_low) + +# -- NT accumulator for Loop 2 signals -- +NT_window_sat = 40.0 # vesicles/s that saturates mGluR and IP3 + # at 20 Hz releasing ~2/spike = 40/s + +# -- eCB retrograde brake -- +tau_eCB_rise = 2000.0 +tau_eCB_decay = 10_000.0 +eCB_threshold = 0.7 # Ca_post level that triggers eCB synthesis + +# -- mGluR presynaptic autoreceptor -- +Km_mGluR = 0.5 +tau_mGluR = 2000.0 # ms +alpha_mGluR = 0.4 # max fractional VGCC suppression + +# -- Astrocyte / IP3 -- +tau_IP3 = 3000.0 # ms +IP3_threshold = 0.8 +wave_boost = 0.2 # conversion_efficiency boost when wave fires +tau_wave_decay = 2 # metabolic cycles before boost decays back + +# -- Glutamine shuttle -- +conversion_efficiency_base = 0.8 + +# -- NT cleft -- +tau_NT_decay = 5.0 # ms + +# ----------------------------------------------------------------------- +# POSTSYNAPTIC PARAMETERS +# ----------------------------------------------------------------------- + +# -- NMDA coincidence detection -- +k_NMDA = 0.08 # Ca_post influx per unit NT * (1 - Mg_block) per ms +V_NMDA_half = 0.3 # V_post at which Mg block is 50% lifted + +# -- Ca_post clearance -- +k_Ca_post_clear = 0.05 # /ms - ATP-dependent PMCA in spine +k_Ca_post_NCX = 0.02 # /ms - ATP-independent NCX floor +ATP_half_post = 0.3 # Hill half-saturation for postsynaptic pumps + +# -- Postsynaptic ATP costs -- +NKA_cost_per_bAP_post = 0.002 # ATP per unit V_post per s (continuous) +ATP_cost_Ca_post_pump = 0.0005 # ATP per unit Ca_post cleared +ATP_demand_scale_post = 50.0 # normalisation (same as presynaptic) + +# -- Receptor desensitization -- +tau_membrane = 20.0 # ms +tau_desensitization = 500.0 # ms + + +# ----------------------------------------------------------------------- +# DENDRITE PARAMETERS +# ----------------------------------------------------------------------- + +# DEND: Single passive dendritic branch connecting postsynaptic spines to soma. +# No active conductances, no spine-neck attenuation, no bAP distance decay. +# The branch sums EPSPs from all active spines (one spine in current model) +# and passes V_dend to the soma each ms. + +tau_dend = 20.0 # DEND: ms - dendritic membrane time constant + # controls how long EPSPs persist before decaying + # longer tau -> broader temporal summation window +AMPA_weight = 0.1 # DEND: scales receptor_conductance -> EPSP contribution + # to V_dend; shared across all spines on the branch + +# bAP: back-propagating AP from soma to all spines (no distance attenuation). +# Generated internally when V_soma crosses threshold (replaces external bAP_train). +V_bAP_peak = 1.0 # DEND: normalised bAP amplitude at all spines +tau_bAP = 3.0 # DEND: ms - bAP decay time constant + # controls width of coincidence window: + # longer tau_bAP -> NT arriving slightly after + # bAP can still achieve NMDA coincidence + +# ----------------------------------------------------------------------- +# SOMA PARAMETERS +# ----------------------------------------------------------------------- + +# SOMA: Leaky integrator with threshold crossing, channel kinetics, and +# refractory period. Firing emerges from V_soma dynamics — not driven by +# an external spike train. Each AP generates a bAP (sent to dendrite) +# and a forward AP (available as output for the next neuron's presynapse). + +tau_soma = 20.0 # SOMA: ms - somatic membrane time constant +soma_weight = 0.5 # SOMA: scales V_dend contribution to V_soma + # reflects electrical coupling efficiency + +V_soma_threshold = 0.5 # SOMA: normalised firing threshold (0->1) + # when V_soma crosses this, AP fires +V_soma_reset = 0.0 # SOMA: V_soma after AP (instantaneous reset + # after repolarisation completes) + +# Channel kinetics — AP waveform profile +# SOMA: The AP is not instantaneous. After threshold crossing: +# (1) Na+ channels open -> V_soma rises to V_AP_peak (depolarisation) +# (2) K+ channels open -> V_soma falls past rest to V_AHP (repolarisation) +# (3) K+ channels close -> V_soma recovers to rest (V_soma_reset) +# tau_AP_rise and tau_AP_fall control the width and shape of the AP waveform. +V_AP_peak = 1.0 # SOMA: normalised AP peak amplitude +V_AHP = -0.1 # SOMA: after-hyperpolarisation trough (below rest) + # negative value: V_soma briefly goes below 0 +tau_AP_rise = 0.5 # SOMA: ms - Na+ channel opening (rising phase) +tau_AP_fall = 1.5 # SOMA: ms - K+ channel repolarisation (falling phase) +tau_AHP = 5.0 # SOMA: ms - recovery from AHP back to rest + +# Refractory period +# SOMA: After an AP fires, the soma cannot fire again until the membrane +# has recovered from inactivation and AHP. +# Absolute refractory: no firing possible regardless of input +# Relative refractory: firing possible but requires stronger input +t_refractory_abs = 2.0 # SOMA: ms - absolute refractory period +t_refractory_rel = 8.0 # SOMA: ms - relative refractory period (total from AP) + # during relative period threshold is elevated + +# ----------------------------------------------------------------------- +# HELPER FUNCTIONS +# ----------------------------------------------------------------------- + +def compute_flux(conductance, voltage): + return k_flux * conductance * abs(voltage) + + +def deterministic_release(N_RRP, Ca_micro, NT_cleft): + # Hill equation: Ca2+ sensor cooperativity (synaptotagmin-1, n=4) + Ca_n = Ca_micro ** n_rel + release_frac = k_rel * Ca_n / (Ca_n + KD_rel ** n_rel) + # NT suppression: physical crowding + fast local autoreceptors + NT_norm = min(1.0, NT_cleft / NT_suppression_sat) + release_frac = release_frac * (1.0 - NT_suppression_weight * NT_norm) + release_frac = np.clip(release_frac, 0.0, 1.0) + return max(0.0, release_frac * N_RRP) + + +def map_trace_to_speed(Tr_Ca): + # Returns /s recruitment rate based on Tr_Ca level + if Tr_Ca > T_high: + return k_rec_fast + elif Tr_Ca < T_low: + return k_rec_slow + else: + t = (Tr_Ca - T_low) / (T_high - T_low) + return k_rec_slow + t * (k_rec_fast - k_rec_slow) + + +def compute_pump_atp_factor(atp, atp_half): + # Hill function: ATP gates pump speed (shared by pre and post) + return (atp ** 2) / (atp ** 2 + atp_half ** 2) + + +def compute_EPSP(receptor_conductance): + return receptor_conductance * 0.1 + + +def compute_astrocyte_metabolic_health(Glucose_level, ATP_demand_accumulated, + demand_scale=50.0): + # Converts glucose supply and accumulated demand into ATP_level (0->1) + # and conversion_efficiency (0->1). Both sides use this function with + # their own demand accumulators but the same Glucose_level — shared + # metabolic vulnerability. + health = np.clip(Glucose_level - ATP_demand_accumulated / demand_scale, + 0.0, 1.0) + return health, health # (conversion_efficiency, ATP_level) + + +def trigger_slow_astrocyte_calcium_wave(): + # Placeholder - gliotransmitter release over ~10 s + pass + + +# ----------------------------------------------------------------------- +# STATE VARIABLES +# ----------------------------------------------------------------------- + +# -- Presynaptic membrane -- +V_pre_state = 0.0 + +# -- Presynaptic Ca2+ -- +Ca_micro = 0.0 +Ca_ER = 0.5 +Ca_buffer_bound = 0.0 +B_free = B_total + +# -- CDI -- +CDI_factor = 0.0 + +# -- Vesicle pools -- +N_RRP = 15.0 +N_RP = 150.0 + +# -- Calcium trace -- +Tr_Ca = 0.0 + +# -- NT cleft -- +NT_cleft = 0.0 + +# -- NT accumulator for slow signals -- +# FIX: this was missing. Accumulates every ms in Loop 1, +# consumed by mGluR and IP3 in Loop 2, reset each second. +NT_released_this_window = 0.0 + +# -- Postsynaptic membrane + receptors -- +V_post = 0.0 +receptor_conductance = 0.0 +Desensitization_level = 0.0 +V_post_history = [] + +# -- Postsynaptic Ca2+ (spine compartment) -- +Ca_post = 0.0 +# Driven by NMDA coincidence (NT + V_post). Cleared by PMCA (ATP-gated) +# and NCX (always). Drives eCB synthesis. No CDI equivalent -> +# elevated Ca_post under ATP failure has no self-limiting mechanism. + +# -- Retrograde / autoreceptor -- +eCB_level = 0.0 +mGluR_activation = 0.0 + +# -- Astrocyte -- +IP3 = 0.0 +wave_active = 0 # countdown: cycles remaining of wave boost +Glutamine_pool = 50.0 + +# -- Presynaptic ATP -- +ATP_level = 1.0 +ATP_demand = 0.0 +conversion_efficiency = conversion_efficiency_base +Glucose_level = 1.0 # set < 1.0 to engage metabolic silencing + +# -- Dendrite -- +V_dend = 0.0 # DEND: dendritic membrane potential (normalised, 0->1) + # sum of attenuated spine EPSPs, decaying each ms + # passed to soma each ms as the integration input +V_bAP = 0.0 # DEND: back-propagating AP amplitude at all spines (0->1) + # set to V_bAP_peak when soma fires + # decays with tau_bAP each ms + # replaces external bAP_train input + +# -- Soma -- +V_soma = 0.0 # SOMA: somatic membrane potential (normalised, 0->1) + # integrates V_dend, decays with tau_soma + # triggers AP when crosses V_soma_threshold +AP_phase = 'rest' # SOMA: current AP waveform phase + # 'rest' | 'rising' | 'falling' | 'ahp' +AP_phase_t = 0.0 # SOMA: ms elapsed in current AP phase +refractory_t = 0.0 # SOMA: ms remaining in refractory period (0 = not refractory) + # absolute refractory if refractory_t > t_refractory_rel - t_refractory_abs + # relative refractory if 0 < refractory_t <= t_refractory_rel - t_refractory_abs +soma_fired = False # SOMA: flag — soma fired this ms + # read by dendrite to trigger V_bAP + # read by simulation output as forward AP signal + +# -- Postsynaptic ATP -- +ATP_level_post = 1.0 # separate pool; same glucose budget as presynaptic +ATP_demand_post = 0.0 # accumulates from NKA (V_post) and PMCA (Ca_post) + + +# ----------------------------------------------------------------------- +# MAIN SIMULATION LOOP +# ----------------------------------------------------------------------- + +def run_simulation(spike_train, total_steps): + """ + spike_train : list of int - presynaptic AP timestep indices + total_steps : int + if None, no bAPs are delivered + """ + + global V_pre_state + global Ca_micro, Ca_ER, Ca_buffer_bound, B_free + global CDI_factor + global N_RRP, N_RP, Tr_Ca, NT_cleft, NT_released_this_window + global V_post, receptor_conductance, Desensitization_level, V_post_history + global Ca_post + global eCB_level, mGluR_activation + global IP3, wave_active, Glutamine_pool + global ATP_level, ATP_demand, conversion_efficiency, Glucose_level + global ATP_level_post, ATP_demand_post + global V_dend, V_bAP + global V_soma, AP_phase, AP_phase_t, refractory_t, soma_fired + + log = {k: [] for k in [ + "V_pre_state", "Ca_micro", "Ca_ER", "CDI_factor", "B_free", + "N_RRP", "N_RP", "Tr_Ca", "NT_cleft", + "V_post", "Ca_post", "eCB_level", "mGluR_activation", + "released_NT", "ATP_level", "ATP_demand", + "ATP_level_post", "ATP_demand_post", + "V_dend", "V_bAP", "V_soma", "soma_fired", + ]} + + 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 + soma_fired = False + + # -- 1A. PRESYNAPTIC MEMBRANE / Na-K-ATPase ------------------- + # AP fires: membrane jumps to peak, then decays with tau_V_pre. + # Ca2+ influx uses V_pre_state (continuous) not binary V_pre, + # giving a temporal influx profile that tapers as membrane repolarises. + if V_pre == 1: + V_pre_state = V_pre_peak + ATP_demand += NKA_cost_per_AP # dominant presynaptic ATP cost + + V_pre_state += (V_rest - V_pre_state) * dt / tau_V_pre + + # -- 1B. PRESYNAPTIC Ca2+ INFLUX ------------------------------ + # Three multiplicative brakes on effective_conductance: + # eCB_level : retrograde brake from postsynapse (Loop 2) + # CDI_factor : Ca2+-dependent inactivation (below) + # mGluR_activation : autoreceptor brake (Loop 2) + effective_conductance = ( + N_VGCC + * (1.0 - eCB_level) + * (1.0 - CDI_factor) + * (1.0 - mGluR_activation * alpha_mGluR) + ) + raw_influx = compute_flux(effective_conductance, V_pre_state) + + # Buffer proteins capture a fraction immediately (fast sponge). + # B_free -> 0 during sustained bursting -> capture_fraction -> 0 + # -> full raw_influx enters Ca_micro (CASCADE 4 acceleration). + capture_fraction = B_free / B_total + captured = raw_influx * capture_fraction + B_free = max(0.0, B_free - captured) + Ca_buffer_bound += captured + Ca_micro += (raw_influx - captured) + + # -- 1C. VESICLE RELEASE -------------------------------------- + # Deterministic: Hill Ca2+ sensor * NT suppression * N_RRP. + # Runs every ms that Ca_micro > 0 (release profile follows Ca2+ + # transient, not locked to spike flag). + if N_RRP > 0 and Ca_micro > 0: + released_NT = deterministic_release(N_RRP, Ca_micro, NT_cleft) + released_NT = min(released_NT, N_RRP) + N_RRP -= released_NT + NT_cleft += released_NT + # FIX: accumulate for Loop 2 mGluR and IP3 signals. + # This is the only correct way to feed slow signals from fast + # events — snapshot of NT_cleft at Loop 2 time would be ~0 + # because passive diffusion has already cleared it. + NT_released_this_window += released_NT + + # Passive NT diffusion out of cleft each ms. + NT_cleft *= (1.0 - dt / tau_NT_decay) + NT_cleft = max(0.0, NT_cleft) + + # -- 1D. PRESYNAPTIC Ca2+ CLEARANCE --------------------------- + # pump_scale: Hill(ATP_level) — bridges Loop 3 ATP to Loop 1 clearance. + # NCX is ATP-independent (floor); PMCA and SERCA are ATP-gated. + pump_scale = compute_pump_atp_factor(ATP_level, ATP_half) + cleared_PMCA = k_PMCA * Ca_micro * pump_scale + cleared_NCX = k_NCX * Ca_micro + cleared_SERCA = k_SERCA * Ca_micro * pump_scale + + Ca_micro -= (cleared_PMCA + cleared_NCX + cleared_SERCA) + Ca_micro = max(0.0, Ca_micro) + Ca_ER += cleared_SERCA + + ATP_demand += ATP_cost_PMCA * cleared_PMCA + ATP_demand += ATP_cost_SERCA * cleared_SERCA + + # Buffer recharge: bound Ca2+ slowly re-releases back to cytosol. + # During pump failure this sustains Ca_micro elevation (CASCADE 4). + rebind_flux = Ca_buffer_bound * dt / tau_buffer_rebind + Ca_micro += rebind_flux + Ca_buffer_bound = max(0.0, Ca_buffer_bound - rebind_flux) + B_free = B_total - Ca_buffer_bound + + # -- 1E. CDI — RISE (spike only) AND RECOVERY (every ms) ------ + # RISE: Ca2+ entering through open channels inactivates them locally. + # Gated to spike window — requires channels to be open. + # (Running every ms was wrong: CDI needs Ca2+ flowing through + # the channel, not ambient cytosolic Ca2+.) + if V_pre == 1: + CDI_factor += k_CDI_rise * Ca_micro * dt_s + + # RECOVERY: continuous, suppressed when Ca_micro is high. + # Self-locking: pump failure -> Ca_micro high -> recovery ~0 + # -> CDI_factor -> 1 -> effective_conductance -> 0 (CASCADE 5-6). + CDI_recovery_rate = k_CDI_rec * (1.0 - Ca_micro / Ca_micro_saturation) + CDI_factor = np.clip(CDI_factor - CDI_recovery_rate * dt_s, 0.0, 1.0) + + # -- 1F. CALCIUM TRACE ---------------------------------------- + # Leaky integrator — integrates full Ca2+ waveform every ms + # including inter-spike clearance. Drives Loop 2 recruitment speed. + Tr_Ca = Tr_Ca + (Ca_micro - Tr_Ca / tau_Tr_Ca) * dt + + # -- 1G. POSTSYNAPTIC: NT DETECTION & AMPA -------------------- + # Desensitization reduces effective NT — sustained NT exposure + # progressively silences receptors (postsynaptic equivalent of CDI). + effective_NT = released_NT * (1.0 - Desensitization_level) + receptor_conductance += effective_NT * 0.05 + receptor_conductance *= (1.0 - dt / tau_membrane) + + V_post += compute_EPSP(receptor_conductance) - (V_post / tau_membrane) * dt + V_post = max(0.0, V_post) + + Desensitization_level += NT_cleft * 0.001 * dt + Desensitization_level -= (Desensitization_level / tau_desensitization) * dt + Desensitization_level = np.clip(Desensitization_level, 0.0, 1.0) + + V_post_history.append(V_post) + if len(V_post_history) > 5000: + V_post_history.pop(0) + + # -- 1H. POSTSYNAPTIC: NMDA COINCIDENCE DETECTION ------------- + # Ca_post enters only when BOTH conditions hold simultaneously: + # (1) NT_cleft > 0 — ligand gate (glutamate present) + # (2) V_post elevated — voltage gate (Mg2+ block lifted) + # V_bAP (from dendrite, generated by soma firing) adds to V_post, + # enabling full Mg block removal only on true pre+post coincidence. + # DEND: V_bAP replaces the old external bAP * 0.5 placeholder. + V_post_effective = V_post + V_bAP # AMPA drive + bAP boost + Mg_block_removal = V_post_effective / (V_post_effective + V_NMDA_half) + NMDA_Ca_influx = k_NMDA * NT_cleft * Mg_block_removal + Ca_post += NMDA_Ca_influx + + # Postsynaptic NKA: membrane recharge cost proportional to V_post. + # [POST-ATP 1] Dominant postsynaptic ATP drain at high activity. + ATP_demand_post += NKA_cost_per_bAP_post * V_post * dt_s + + # -- 1I. POSTSYNAPTIC: Ca_post CLEARANCE ---------------------- + # pump_scale_post: Hill(ATP_level_post) — same structure as presynaptic. + # NCX is ATP-independent floor (enables auto-reset after ATP recovery). + # [POST-ATP 3] When pump_scale_post falls, Ca_post stays elevated -> + # eCB threshold crossed without genuine coincidence -> false retrograde. + pump_scale_post = compute_pump_atp_factor(ATP_level_post, ATP_half_post) + cleared_Ca_post_pump = k_Ca_post_clear * Ca_post * pump_scale_post + cleared_Ca_post_NCX = k_Ca_post_NCX * Ca_post + Ca_post -= (cleared_Ca_post_pump + cleared_Ca_post_NCX) + Ca_post = max(0.0, Ca_post) + + # [POST-ATP 2] ATP cost of postsynaptic PMCA. + ATP_demand_post += ATP_cost_Ca_post_pump * cleared_Ca_post_pump + + # -- 1J. DENDRITE: EPSP SUMMATION & bAP DISTRIBUTION ---------- + # DEND: The dendritic branch collects the EPSP from this spine + # (receptor_conductance * AMPA_weight) and adds it to V_dend. + # V_dend then decays passively with tau_dend. + # No spine-neck attenuation in this simplified model — + # all spines contribute equally regardless of position. + V_dend += receptor_conductance * AMPA_weight + V_dend *= (1.0 - dt / tau_dend) + V_dend = max(0.0, V_dend) + + # DEND: bAP distribution — set by soma firing (section 1K below). + # Decays each ms with tau_bAP. All spines receive the same amplitude + # (no distance attenuation in this simplified model). + V_bAP += (0.0 - V_bAP) * dt / tau_bAP + V_bAP = max(0.0, V_bAP) + + # -- 1K. SOMA: INTEGRATION, AP KINETICS, REFRACTORY -------------- + # SOMA: V_soma integrates V_dend as a leaky integrator. + # When V_soma crosses V_soma_threshold (and not refractory), + # an AP fires. The AP has a three-phase waveform: + # rising : Na+ channels open -> V_soma climbs to V_AP_peak + # falling : K+ channels open -> V_soma falls to V_AHP + # ahp : K+ channels close -> V_soma recovers toward rest + # After the waveform completes, the soma enters the refractory period. + # Absolute refractory: no firing possible (Na+ channels inactivated). + # Relative refractory: threshold is effectively elevated. + + # Step 1: integrate dendritic input (only when not in AP waveform) + if AP_phase == 'rest': + V_soma += V_dend * soma_weight + V_soma *= (1.0 - dt / tau_soma) + V_soma = max(V_AHP, V_soma) + + # Threshold check — blocked during refractory period. + # During relative refractory (0 < refractory_t <= t_refractory_rel): + # effective threshold is raised proportionally to remaining time. + abs_ref_remaining = refractory_t - (t_refractory_rel - t_refractory_abs) + in_absolute = abs_ref_remaining > 0 + effective_threshold = V_soma_threshold + if refractory_t > 0 and not in_absolute: + # Linear threshold elevation during relative refractory + rel_fraction = refractory_t / t_refractory_rel + effective_threshold = V_soma_threshold * (1.0 + rel_fraction) + + if V_soma >= effective_threshold and not in_absolute: + # AP fires: enter rising phase + AP_phase = 'rising' + AP_phase_t = 0.0 + soma_fired = True + refractory_t = t_refractory_rel # start refractory countdown + + # DEND: bAP generated — broadcast to all spines immediately + V_bAP = V_bAP_peak + + # Step 2: AP waveform phases + elif AP_phase == 'rising': + AP_phase_t += dt + # V_soma rises exponentially toward V_AP_peak + V_soma += (V_AP_peak - V_soma) * dt / tau_AP_rise + if AP_phase_t >= tau_AP_rise * 3: # ~3 time constants = near peak + AP_phase = 'falling' + AP_phase_t = 0.0 + + elif AP_phase == 'falling': + AP_phase_t += dt + # V_soma falls exponentially toward V_AHP (after-hyperpolarisation) + V_soma += (V_AHP - V_soma) * dt / tau_AP_fall + if AP_phase_t >= tau_AP_fall * 3: + AP_phase = 'ahp' + AP_phase_t = 0.0 + + elif AP_phase == 'ahp': + AP_phase_t += dt + # V_soma recovers from AHP toward rest (V_soma_reset) + V_soma += (V_soma_reset - V_soma) * dt / tau_AHP + if AP_phase_t >= tau_AHP * 3: + AP_phase = 'rest' + AP_phase_t = 0.0 + V_soma = V_soma_reset + + # Step 3: refractory countdown (runs every ms regardless of phase) + if refractory_t > 0: + refractory_t = max(0.0, refractory_t - dt) + + # -- RECORD --------------------------------------------------- + log["V_pre_state"].append(V_pre_state) + log["Ca_micro"].append(Ca_micro) + log["Ca_ER"].append(Ca_ER) + log["CDI_factor"].append(CDI_factor) + log["B_free"].append(B_free) + log["N_RRP"].append(N_RRP) + log["N_RP"].append(N_RP) + log["Tr_Ca"].append(Tr_Ca) + log["NT_cleft"].append(NT_cleft) + log["V_post"].append(V_post) + log["Ca_post"].append(Ca_post) + log["eCB_level"].append(eCB_level) + log["mGluR_activation"].append(mGluR_activation) + log["released_NT"].append(released_NT) + log["ATP_level"].append(ATP_level) + log["ATP_demand"].append(ATP_demand) + log["ATP_level_post"].append(ATP_level_post) + log["ATP_demand_post"].append(ATP_demand_post) + log["V_dend"].append(V_dend) + log["V_bAP"].append(V_bAP) + log["V_soma"].append(V_soma) + log["soma_fired"].append(float(soma_fired)) + + # ============================================================== + # LOOP 2 — SLOW / ASTROCYTE (dt_slow = 1 s) + # ============================================================== + + if (step % High_Freq_Multiplier) == 0: + + # Astrocyte EAAT clearance — active NT removal from cleft. + cleared_NT = NT_cleft * 0.3 + NT_cleft = max(0.0, NT_cleft - cleared_NT) + + # FIX: IP3 integrates NT_released_this_window (total release + # since last Loop 2), not the post-diffusion NT_cleft residual + # which is ~0 by the time Loop 2 runs. + IP3 += NT_released_this_window - (IP3 / tau_IP3) * dt_slow + IP3 = max(0.0, IP3) + + if IP3 > IP3_threshold: + trigger_slow_astrocyte_calcium_wave() + # FIX: wave boosts conversion_efficiency in the next mins cycle. + # The astrocyte responds to heavy load by upregulating its + # recycling machinery — shipping more glutamine back to the + # presynapse. Boost decays over tau_wave_decay metabolic cycles. + wave_active = tau_wave_decay + + # FIX: mGluR reads NT_released_this_window (accumulated release + # load), not NT_cleft snapshot. NT_cleft is ~0 at Loop 2 time + # due to diffusion; the accumulator correctly represents the + # burst load the autoreceptor has sensed during this window. + NT_window_norm = min(1.0, NT_released_this_window / NT_window_sat) + mGluR_target = NT_window_norm + mGluR_activation += (mGluR_target - mGluR_activation) * (dt_slow / tau_mGluR) + mGluR_activation = np.clip(mGluR_activation, 0.0, 1.0) + + # FIX: reset accumulator for next window. + NT_released_this_window = 0.0 + + # eCB retrograde synthesis: now driven by Ca_post (spine Ca2+), + # not V_post_history. The actual eCB synthesis in the spine is + # triggered by Ca2+-dependent enzymes (DAGL, PLC), not voltage. + # Under normal conditions Ca_post only rises with coincidence. + # Under POST-ATP failure Ca_post stays elevated without genuine + # coincidence -> false retrograde signal (POST-ATP 4). + recent_Ca_post = (np.mean(log["Ca_post"][-2000:]) + if len(log["Ca_post"]) >= 2000 + else (np.mean(log["Ca_post"]) if log["Ca_post"] else 0.0)) + eCB_signal = max(0.0, recent_Ca_post - eCB_threshold) + if eCB_signal > 0: + eCB_level += eCB_signal * (dt_slow / tau_eCB_rise) + else: + eCB_level -= eCB_level * (dt_slow / tau_eCB_decay) + eCB_level = np.clip(eCB_level, 0.0, 1.0) + + # FIX: RP->RRP recruitment moved here from Loop 1. + # Biological timescale: vesicle docking and priming take seconds, + # not milliseconds. k_rec_fast/slow are /s; * dt_slow_s = 1.0 s + # gives dimensionless per-step fraction — no hidden unit scaling. + current_recruitment_rate = map_trace_to_speed(Tr_Ca) # /s + refill_amount = (current_recruitment_rate * dt_slow_s + * N_RP * (Max_RRP - N_RRP) / Max_RRP) + refill_amount = max(0.0, refill_amount) + refill_amount = min(refill_amount, N_RP) + + N_RRP = min(N_RRP + refill_amount, Max_RRP) + N_RP = max(0.0, N_RP - refill_amount) + ATP_demand += ATP_cost_docking * refill_amount + + # ============================================================== + # LOOP 3 — METABOLIC (dt_meta = 1 min) + # ============================================================== + + if (step % Metabolic_Multiplier) == 0: + + # Presynaptic ATP: glucose supply minus accumulated demand. + conversion_efficiency, ATP_level = compute_astrocyte_metabolic_health( + Glucose_level, ATP_demand + ) + ATP_demand = 0.0 + + # FIX: wave boost applied to conversion_efficiency. + # Astrocyte calcium wave (triggered by high IP3) upregulates + # glutamine synthetase -> faster NT recycling -> more RP refill. + # Boost decays over tau_wave_decay cycles. + if wave_active > 0: + conversion_efficiency = min(1.0, conversion_efficiency + wave_boost) + wave_active -= 1 + + # Glutamine shuttle: astrocyte converts cleared NT to glutamine, + # presynapse repackages it into vesicles -> N_RP replenished. + refill_RP = Glutamine_pool * conversion_efficiency + N_RP = min(Max_RP, N_RP + refill_RP) + Glutamine_pool = max(0.0, Glutamine_pool - refill_RP) + + # Postsynaptic ATP: same glucose budget, own demand accumulator. + # Both sides draw from Glucose_level -> shared metabolic vulnerability. + # Presynaptic silence reduces NT -> less NMDA -> less Ca_post -> + # less ATP_demand_post: presynaptic protection indirectly + # protects the postsynapse. + _, ATP_level_post = compute_astrocyte_metabolic_health( + Glucose_level, ATP_demand_post, ATP_demand_scale_post + ) + ATP_demand_post = 0.0 + + return log + + +# ----------------------------------------------------------------------- +# EXAMPLE USAGE +# ----------------------------------------------------------------------- + +if __name__ == "__main__": + import matplotlib.pyplot as plt + + total_steps = 10_000 # 10 seconds + + # Presynaptic 20 Hz burst for 2 s. + spike_train = list(range(0, 2000, 50)) + + # Soma firing emerges from V_soma threshold crossings — no external bAP_train. + results = run_simulation(spike_train, total_steps) + t = np.arange(total_steps) * dt + + fig, axes = plt.subplots(8, 1, figsize=(12, 18), sharex=True) + fig.suptitle("Tripartite Synapse — Presynaptic + Postsynaptic", fontsize=13) + + axes[0].plot(t, results["V_pre_state"], color="slateblue", lw=0.8) + axes[0].set_ylabel("V_pre") + axes[0].set_title("Presynaptic membrane (AP waveform)", fontsize=9, loc="left") + + axes[1].plot(t, results["Ca_micro"], color="darkorange", lw=0.8) + axes[1].set_ylabel("[Ca2+] pre") + axes[1].set_title("CASCADE 4 — presynaptic Ca2+", fontsize=9, loc="left") + + axes[2].plot(t, results["CDI_factor"], color="firebrick", lw=0.8, label="CDI") + axes[2].plot(t, results["B_free"], color="steelblue", lw=0.8, label="Buffer free") + axes[2].set_ylabel("CDI / Buffer") + axes[2].set_title("CASCADE 5 — CDI lock-out", fontsize=9, loc="left") + axes[2].legend(fontsize=8) + + axes[3].plot(t, results["N_RRP"], color="teal", lw=0.8, label="RRP") + axes[3].plot(t, results["N_RP"], color="purple", lw=0.8, label="RP") + axes[3].set_ylabel("Vesicles") + axes[3].set_title("CASCADE 1 — vesicle depletion", fontsize=9, loc="left") + axes[3].legend(fontsize=8) + + axes[4].plot(t, results["NT_cleft"], color="darkgreen", lw=0.8, label="NT cleft") + axes[4].plot(t, results["mGluR_activation"], color="saddlebrown", lw=0.8, label="mGluR") + axes[4].plot(t, results["eCB_level"], color="crimson", lw=0.8, label="eCB") + axes[4].set_ylabel("Cleft / Feedback") + axes[4].set_title("CASCADE 6 — three brakes on conductance", fontsize=9, loc="left") + axes[4].legend(fontsize=8) + + axes[5].plot(t, results["V_post"], color="navy", lw=0.8, label="V_post") + axes[5].plot(t, results["Ca_post"], color="coral", lw=0.8, label="Ca_post (spine)") + axes[5].set_ylabel("Postsynaptic") + axes[5].set_title("Postsynaptic potential + NMDA spine Ca2+", fontsize=9, loc="left") + axes[5].legend(fontsize=8) + + axes[6].plot(t, results["ATP_level"], color="goldenrod", lw=0.8, label="ATP pre") + axes[6].plot(t, results["ATP_level_post"], color="darkorange", lw=0.8, label="ATP post") + axes[6].set_ylabel("ATP level") + axes[6].set_title("CASCADE 2 / POST-ATP — presynaptic and postsynaptic ATP", fontsize=9, loc="left") + axes[6].legend(fontsize=8) + + axes[7].plot(t, results["ATP_demand"], color="tomato", lw=0.8, label="demand pre") + axes[7].plot(t, results["ATP_demand_post"], color="orangered", lw=0.8, label="demand post") + axes[7].set_ylabel("ATP demand") + axes[7].set_title("Accumulated ATP demand (resets each min cycle)", fontsize=9, loc="left") + axes[7].set_xlabel("Time (ms)") + axes[7].legend(fontsize=8) + + fig2, ax2 = plt.subplots(3, 1, figsize=(12, 8), sharex=True) + fig2.suptitle("Dendrite + Soma", fontsize=13) + ax2[0].plot(t, results["V_dend"], color="mediumblue", lw=0.8) + ax2[0].set_ylabel("V_dend") + ax2[0].set_title("DEND — summed EPSPs (leaky integrator)", fontsize=9, loc="left") + ax2[1].plot(t, results["V_soma"], color="darkgreen", lw=0.8) + ax2[1].axhline(V_soma_threshold, color="red", lw=0.5, ls="--", label="threshold") + ax2[1].set_ylabel("V_soma") + ax2[1].set_title("SOMA — membrane potential + threshold (dashed)", fontsize=9, loc="left") + ax2[1].legend(fontsize=8) + ax2[2].plot(t, results["V_bAP"], color="darkorchid", lw=0.8) + ax2[2].plot(t, results["soma_fired"], color="crimson", lw=0.5, alpha=0.5, label="fired") + ax2[2].set_ylabel("V_bAP / fired") + ax2[2].set_title("DEND — bAP distributed to spines on soma firing", fontsize=9, loc="left") + ax2[2].set_xlabel("Time (ms)") + ax2[2].legend(fontsize=8) + fig2.tight_layout() + fig2.savefig("./dendrite_soma.png", dpi=150) + plt.tight_layout() + plt.savefig("./synapse_simulation.png", dpi=150) + plt.close() + print("Done.") \ No newline at end of file