cdi traces
This commit is contained in:
@@ -1,493 +0,0 @@
|
||||
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()
|
||||
@@ -1,648 +0,0 @@
|
||||
# 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.")
|
||||
Reference in New Issue
Block a user