organism/neuron/appunti/traditional_python_models/tripartite.py

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()