From 1fe72a277b170ac4787eedc88426ae61787c8495 Mon Sep 17 00:00:00 2001 From: ocrampal Date: Thu, 11 Jun 2026 18:59:09 +0200 Subject: [PATCH] old --- .../appunti/simple_interpreter/example.snip | 40 - .../appunti/simple_interpreter/interpreter.py | 480 ---------- .../Tripartite-Synapse.py | 690 -------------- .../tripartite-new-1.py | 714 -------------- .../tripartite-soma.py | 883 ------------------ 5 files changed, 2807 deletions(-) delete mode 100644 elements/neuron/appunti/simple_interpreter/example.snip delete mode 100644 elements/neuron/appunti/simple_interpreter/interpreter.py delete mode 100644 elements/neuron/appunti/traditional_python_models/Tripartite-Synapse.py delete mode 100644 elements/neuron/appunti/traditional_python_models/tripartite-new-1.py delete mode 100644 elements/neuron/appunti/traditional_python_models/tripartite-soma.py diff --git a/elements/neuron/appunti/simple_interpreter/example.snip b/elements/neuron/appunti/simple_interpreter/example.snip deleted file mode 100644 index 0099244..0000000 --- a/elements/neuron/appunti/simple_interpreter/example.snip +++ /dev/null @@ -1,40 +0,0 @@ -# example.snip — demonstrates the Snippet language - -snippet init - set x = 5 - set y = 10 - print x - print y - -snippet adjust_x - # If x is small, bump it up; otherwise pull it back down - if x < 8 - increase x - increase x - increase x - if x > 6 - increase x - else - decrease x - else - decrease x - print x - -snippet adjust_y - if y > 8 - decrease y - if y > 5 - decrease y - decrease y - else - increase y - else - increase y - print y - -snippet final - # Both snippets share variables — x and y from above are visible here - if x > 7 - print x - else - print y diff --git a/elements/neuron/appunti/simple_interpreter/interpreter.py b/elements/neuron/appunti/simple_interpreter/interpreter.py deleted file mode 100644 index a6e6bb5..0000000 --- a/elements/neuron/appunti/simple_interpreter/interpreter.py +++ /dev/null @@ -1,480 +0,0 @@ -#!/usr/bin/env python3 -""" -╔══════════════════════════════════════════════════════════════╗ -║ S N I P P E T L A N G U A G E ║ -║ Interpreter v1.0 ║ -╠══════════════════════════════════════════════════════════════╣ -║ Syntax: ║ -║ snippet ║ -║ set = # assign a variable ║ -║ increase # var += 1 ║ -║ decrease # var -= 1 ║ -║ print # print variable value ║ -║ if > # conditional (> or <) ║ -║ ... ║ -║ else ║ -║ ... ║ -║ ║ -║ Rules: ║ -║ • Snippets run top-to-bottom and share all variables ║ -║ • Indentation defines scope (spaces or tabs, be consistent) -║ • No loops ║ -║ • Numbers only (integers and decimals) ║ -║ • Lines starting with # are comments ║ -╚══════════════════════════════════════════════════════════════╝ - -Usage: - python interpreter.py program.snip # run a file - python interpreter.py # REPL (blank line to run) -""" - -import sys -import re -from dataclasses import dataclass -from typing import Any, Dict, List, Optional, Tuple - - -# ══════════════════════════════════════════════════════════════ -# TOKEN TYPES -# ══════════════════════════════════════════════════════════════ - -class TT: - SNIPPET = 'SNIPPET' - IF = 'IF' - ELSE = 'ELSE' - INCREASE = 'INCREASE' - DECREASE = 'DECREASE' - SET = 'SET' - PRINT = 'PRINT' - NAME = 'NAME' - NUMBER = 'NUMBER' - GT = 'GT' - LT = 'LT' - EQ = 'EQ' - -KEYWORDS: Dict[str, str] = { - 'snippet': TT.SNIPPET, - 'if': TT.IF, - 'else': TT.ELSE, - 'increase': TT.INCREASE, - 'decrease': TT.DECREASE, - 'set': TT.SET, - 'print': TT.PRINT, -} - - -@dataclass -class Token: - type: str - value: Any - line: int - - def __repr__(self) -> str: - return f"Token({self.type}, {self.value!r}, line={self.line})" - - -# ══════════════════════════════════════════════════════════════ -# LEXER → list of (indent, [Token …], line_num) -# ══════════════════════════════════════════════════════════════ - -class LexerError(Exception): - pass - - -class Lexer: - def __init__(self, source: str): - self.source = source - - def tokenize(self) -> List[Tuple[int, List[Token], int]]: - result: List[Tuple[int, List[Token], int]] = [] - for line_num, raw in enumerate(self.source.splitlines(), 1): - stripped = raw.rstrip() - if not stripped or stripped.lstrip().startswith('#'): - continue - indent = len(raw) - len(raw.lstrip(' \t')) - tokens = self._lex_words(stripped.lstrip().split(), line_num) - if tokens: - result.append((indent, tokens, line_num)) - return result - - def _lex_words(self, words: List[str], line_num: int) -> List[Token]: - tokens: List[Token] = [] - for word in words: - if word in KEYWORDS: - tokens.append(Token(KEYWORDS[word], word, line_num)) - elif re.fullmatch(r'-?\d+(\.\d+)?', word): - tokens.append(Token(TT.NUMBER, float(word), line_num)) - elif word == '>': - tokens.append(Token(TT.GT, '>', line_num)) - elif word == '<': - tokens.append(Token(TT.LT, '<', line_num)) - elif word == '=': - tokens.append(Token(TT.EQ, '=', line_num)) - elif re.fullmatch(r'[a-zA-Z_]\w*', word): - tokens.append(Token(TT.NAME, word, line_num)) - else: - raise LexerError(f"[Line {line_num}] Unknown token: {word!r}") - return tokens - - -# ══════════════════════════════════════════════════════════════ -# AST NODES -# ══════════════════════════════════════════════════════════════ - -@dataclass -class ProgramNode: - snippets: List - -@dataclass -class SnippetNode: - name: str - body: List - line: int - -@dataclass -class ConditionNode: - left: str # variable name - op: str # '>' or '<' - right: float # literal number - line: int - -@dataclass -class IfNode: - condition: ConditionNode - then_body: List - else_body: List - line: int - -@dataclass -class IncreaseNode: - name: str - line: int - -@dataclass -class DecreaseNode: - name: str - line: int - -@dataclass -class SetNode: - name: str - value: float - line: int - -@dataclass -class PrintNode: - name: str - line: int - - -# ══════════════════════════════════════════════════════════════ -# PARSER -# ══════════════════════════════════════════════════════════════ - -class ParseError(Exception): - pass - - -# Shorthand type for a lexed line -LexLine = Tuple[int, List[Token], int] - - -class Parser: - def __init__(self, lines: List[LexLine]): - self.lines = lines - self.pos = 0 - - # ── low-level helpers ───────────────────────────────────── - - def _peek(self) -> Optional[LexLine]: - return self.lines[self.pos] if self.pos < len(self.lines) else None - - def _consume(self) -> LexLine: - line = self.lines[self.pos] - self.pos += 1 - return line - - # ── top level ───────────────────────────────────────────── - - def parse(self) -> ProgramNode: - snippets = [] - while self._peek(): - indent, tokens, line_num = self._peek() - if tokens[0].type != TT.SNIPPET: - raise ParseError( - f"[Line {line_num}] Expected 'snippet', got {tokens[0].value!r}" - ) - snippets.append(self._parse_snippet()) - return ProgramNode(snippets) - - # ── snippet ─────────────────────────────────────────────── - - def _parse_snippet(self) -> SnippetNode: - indent, tokens, line_num = self._consume() - if len(tokens) < 2 or tokens[1].type != TT.NAME: - raise ParseError(f"[Line {line_num}] Expected a name after 'snippet'") - name = tokens[1].value - body = self._parse_block(parent_indent=indent) - if not body: - raise ParseError(f"[Line {line_num}] Snippet '{name}' has an empty body") - return SnippetNode(name, body, line_num) - - # ── block ───────────────────────────────────────────────── - - def _parse_block(self, parent_indent: int) -> List: - """ - Collect statements whose indentation is strictly greater than - parent_indent and all equal to the first line's indentation. - """ - nxt = self._peek() - if nxt is None or nxt[0] <= parent_indent: - return [] - - block_indent = nxt[0] - stmts = [] - - while True: - nxt = self._peek() - if nxt is None or nxt[0] < block_indent: - break - indent, tokens, line_num = nxt - if indent > block_indent: - raise ParseError( - f"[Line {line_num}] Unexpected indentation " - f"(expected {block_indent}, got {indent})" - ) - if tokens[0].type == TT.ELSE: - break # parent if-parser will handle this - stmts.append(self._parse_statement(block_indent)) - - return stmts - - # ── statements ──────────────────────────────────────────── - - def _parse_statement(self, block_indent: int): - indent, tokens, line_num = self._peek() - tt = tokens[0].type - dispatch = { - TT.IF: self._parse_if, - TT.INCREASE: self._parse_increase, - TT.DECREASE: self._parse_decrease, - TT.SET: self._parse_set, - TT.PRINT: self._parse_print, - } - if tt not in dispatch: - raise ParseError( - f"[Line {line_num}] Unexpected keyword: {tokens[0].value!r}" - ) - return dispatch[tt](block_indent) if tt == TT.IF else dispatch[tt]() - - def _parse_if(self, current_indent: int) -> IfNode: - indent, tokens, line_num = self._consume() - # Expected shape: IF NAME (GT|LT) NUMBER - if len(tokens) < 4: - raise ParseError(f"[Line {line_num}] Incomplete 'if' condition") - if tokens[1].type != TT.NAME: - raise ParseError(f"[Line {line_num}] Expected variable name after 'if'") - if tokens[2].type not in (TT.GT, TT.LT): - raise ParseError(f"[Line {line_num}] Expected '>' or '<' in condition") - if tokens[3].type != TT.NUMBER: - raise ParseError(f"[Line {line_num}] Expected a number in condition") - - cond = ConditionNode( - left = tokens[1].value, - op = tokens[2].value, - right = tokens[3].value, - line = line_num, - ) - then_body = self._parse_block(parent_indent=indent) - - # Optional else at the SAME indent level as the if - else_body: List = [] - nxt = self._peek() - if nxt and nxt[0] == indent and nxt[1][0].type == TT.ELSE: - self._consume() # eat 'else' - else_body = self._parse_block(parent_indent=indent) - - return IfNode(cond, then_body, else_body, line_num) - - def _parse_increase(self) -> IncreaseNode: - indent, tokens, line_num = self._consume() - if len(tokens) < 2 or tokens[1].type != TT.NAME: - raise ParseError(f"[Line {line_num}] Expected variable name after 'increase'") - return IncreaseNode(tokens[1].value, line_num) - - def _parse_decrease(self) -> DecreaseNode: - indent, tokens, line_num = self._consume() - if len(tokens) < 2 or tokens[1].type != TT.NAME: - raise ParseError(f"[Line {line_num}] Expected variable name after 'decrease'") - return DecreaseNode(tokens[1].value, line_num) - - def _parse_set(self) -> SetNode: - indent, tokens, line_num = self._consume() - # set NAME = NUMBER - if len(tokens) < 4: - raise ParseError(f"[Line {line_num}] Syntax: set = ") - if tokens[1].type != TT.NAME: - raise ParseError(f"[Line {line_num}] Expected variable name after 'set'") - if tokens[2].type != TT.EQ: - raise ParseError(f"[Line {line_num}] Expected '=' after variable name") - if tokens[3].type != TT.NUMBER: - raise ParseError(f"[Line {line_num}] Expected a number after '='") - return SetNode(tokens[1].value, tokens[3].value, line_num) - - def _parse_print(self) -> PrintNode: - indent, tokens, line_num = self._consume() - if len(tokens) < 2 or tokens[1].type != TT.NAME: - raise ParseError(f"[Line {line_num}] Expected variable name after 'print'") - return PrintNode(tokens[1].value, line_num) - - -# ══════════════════════════════════════════════════════════════ -# INTERPRETER (tree-walk evaluator) -# ══════════════════════════════════════════════════════════════ - -class SnippetRuntimeError(Exception): - pass - - -class Interpreter: - def __init__(self): - # All snippets share this variable store - self.variables: Dict[str, float] = {} - - def run(self, program: ProgramNode) -> None: - for snippet in program.snippets: - self._exec_snippet(snippet) - - # ── snippet ─────────────────────────────────────────────── - - def _exec_snippet(self, node: SnippetNode) -> None: - print(f"── snippet {node.name} ──") - self._exec_block(node.body) - - # ── block / dispatch ────────────────────────────────────── - - def _exec_block(self, stmts: List) -> None: - for stmt in stmts: - self._exec(stmt) - - def _exec(self, node) -> None: - dispatch = { - IfNode: self._exec_if, - IncreaseNode: self._exec_increase, - DecreaseNode: self._exec_decrease, - SetNode: self._exec_set, - PrintNode: self._exec_print, - } - handler = dispatch.get(type(node)) - if handler is None: - raise SnippetRuntimeError(f"Unknown AST node: {type(node).__name__}") - handler(node) - - # ── statements ──────────────────────────────────────────── - - def _exec_if(self, node: IfNode) -> None: - if self._eval_condition(node.condition): - self._exec_block(node.then_body) - else: - self._exec_block(node.else_body) - - def _eval_condition(self, cond: ConditionNode) -> bool: - val = self._get_var(cond.left, cond.line) - if cond.op == '>': - return val > cond.right - if cond.op == '<': - return val < cond.right - raise SnippetRuntimeError(f"Unknown operator: {cond.op!r}") - - def _exec_increase(self, node: IncreaseNode) -> None: - self.variables[node.name] = self._get_var(node.name, node.line) + 1 - - def _exec_decrease(self, node: DecreaseNode) -> None: - self.variables[node.name] = self._get_var(node.name, node.line) - 1 - - def _exec_set(self, node: SetNode) -> None: - self.variables[node.name] = node.value - - def _exec_print(self, node: PrintNode) -> None: - val = self._get_var(node.name, node.line) - # Show as int when the value is a whole number - display = int(val) if val == int(val) else val - print(f" {node.name} = {display}") - - # ── helpers ─────────────────────────────────────────────── - - def _get_var(self, name: str, line: int) -> float: - if name not in self.variables: - raise SnippetRuntimeError( - f"[Line {line}] Undefined variable: '{name}'" - ) - return self.variables[name] - - -# ══════════════════════════════════════════════════════════════ -# PUBLIC ENTRY POINT -# ══════════════════════════════════════════════════════════════ - -def run_source(source: str) -> None: - """Lex → parse → interpret a source string.""" - try: - lines = Lexer(source).tokenize() - ast = Parser(lines).parse() - Interpreter().run(ast) - except (LexerError, ParseError, SnippetRuntimeError) as err: - print(f"\n ✖ {err}", file=sys.stderr) - sys.exit(1) - - -# ══════════════════════════════════════════════════════════════ -# COMMAND LINE -# ══════════════════════════════════════════════════════════════ - -REPL_BANNER = """\ -╔════════════════════════════════╗ -║ Snippet Language REPL v1.0 ║ -║ blank line → run ║ -║ 'exit' → quit ║ -╚════════════════════════════════╝ -""" - -def _repl() -> None: - print(REPL_BANNER) - while True: - collected: List[str] = [] - prompt = '> ' - try: - while True: - line = input(prompt) - if line.strip().lower() == 'exit': - sys.exit(0) - if line == '' and collected: - break - collected.append(line) - prompt = ' ' - except EOFError: - break - if collected: - run_source('\n'.join(collected)) - print() - - -def main() -> None: - if len(sys.argv) == 1: - _repl() - elif len(sys.argv) == 2: - path = sys.argv[1] - try: - with open(path) as fh: - source = fh.read() - except FileNotFoundError: - print(f"File not found: {path}", file=sys.stderr) - sys.exit(1) - run_source(source) - else: - print(f"Usage: {sys.argv[0]} [file.snip]", file=sys.stderr) - sys.exit(1) - - -if __name__ == '__main__': - main() diff --git a/elements/neuron/appunti/traditional_python_models/Tripartite-Synapse.py b/elements/neuron/appunti/traditional_python_models/Tripartite-Synapse.py deleted file mode 100644 index d16ec92..0000000 --- a/elements/neuron/appunti/traditional_python_models/Tripartite-Synapse.py +++ /dev/null @@ -1,690 +0,0 @@ -# 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 -# 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 - - -# ----------------------------------------------------------------------- -# 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 - -# -- 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, bAP_train=None): - """ - spike_train : list of int - presynaptic AP timestep indices - total_steps : int - bAP_train : list of int - postsynaptic bAP timestep indices (optional) - 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 - - 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", - ]} - - spike_set = set(spike_train) - bAP_set = set(bAP_train) if bAP_train else set() - - for step in range(total_steps): - - # ============================================================== - # LOOP 1 — HIGH-FREQUENCY (dt = 1 ms) - # ============================================================== - - V_pre = 1 if step in spike_set else 0 - bAP = 1 if step in bAP_set else 0 - released_NT = 0.0 - - # -- 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) - # Mg block removal is a sigmoid of V_post. - # bAP (backpropagating AP) boosts V_post further, enabling - # full NMDA opening and larger Ca_post surge. - V_post_effective = V_post + (bAP * 0.5) # bAP adds depolarisation - 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 - - # -- 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) - - # ============================================================== - # 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)) - - # Postsynaptic bAPs coincident with every 5th presynaptic spike - # (simulates partial coincidence for NMDA activation). - bAP_train = list(range(0, 2000, 250)) - - results = run_simulation(spike_train, total_steps, bAP_train) - 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) - - plt.tight_layout() - plt.savefig("./synapse_simulation.png", dpi=150) - plt.close() - print("Done.") diff --git a/elements/neuron/appunti/traditional_python_models/tripartite-new-1.py b/elements/neuron/appunti/traditional_python_models/tripartite-new-1.py deleted file mode 100644 index 254fe55..0000000 --- a/elements/neuron/appunti/traditional_python_models/tripartite-new-1.py +++ /dev/null @@ -1,714 +0,0 @@ -# Tripartite Synapse - Multi-Scale Computational Model -# ===================================================== -# Presynaptic perspective, with all missing behaviors integrated. -# -# Change log: -# ORIG - present from the original document -# NEW - added in the missing-behavior integration pass -# DET - this pass: deterministic Ca2+-driven vesicle release -# NKA - this pass: explicit Na/K-ATPase V_pre decay and ATP cost -# -# Clock structure: -# Loop 1 - dt = 1 ms (Ca2+, vesicle release, short-term traces) -# Loop 2 - dt = 1000 ms (astrocyte clearance, eCB, mGluR feedback) -# Loop 3 - dt = 60000 ms (glutamine shuttle, metabolic health) -# -# ======================================================================= -# METABOLIC SILENCING CASCADE - variable map -# ======================================================================= -# Each step of the cascade is tagged inline with [CASCADE n]. -# -# [CASCADE 1] HIGH FIRING RATE -> VESICLE DEPLETION (fast, ~seconds) -# Driver : spike_train density -> released_NT per ms -# Victim : N_RRP, N_RP -# Bottleneck : RP->RRP recruitment cannot keep up at high rates. -# Outcome : N_RRP -> 0, released_NT -> 0, NT_cleft collapses. -# -# [CASCADE 2] HIGH FIRING RATE -> ATP DEPLETION (slow, ~minutes) -# Driver : Na/K-ATPase recharge cost per AP (NKA_cost_per_AP) -# + PMCA/SERCA pump load + vesicle re-docking -# Victim : ATP_demand accumulator -> ATP_level -# Bottleneck : Glucose_level sets replenishment ceiling. -# Outcome : ATP_level < 1, pump_scale < 1. -# -# [CASCADE 3] LOW ATP -> PUMP FAILURE (PMCA / SERCA slow) -# Driver : pump_scale = Hill(ATP_level) -# Victim : cleared_PMCA, cleared_SERCA -# Outcome : total Ca2+ clearance rate drops. -# -# [CASCADE 4] PUMP FAILURE -> RESIDUAL Ca2+ STAYS HIGH -# Driver : reduced clearance + saturated buffer -# Victim : Ca_micro -# Outcome : Ca_micro persists between spikes. -# -# [CASCADE 5] RESIDUAL Ca2+ -> CDI STAYS ACTIVE (VGCCs lock shut) -# Driver : Ca_micro > 0 between spikes -# Victim : CDI_factor -# Outcome : CDI_factor -> 1, effective_conductance -> 0. -# -# [CASCADE 6] RESULT - SYNAPSE SILENCES (excitotoxicity protection) -# Driver : CDI_factor ~= 1 -# Mechanism : effective_conductance = N_VGCC -# * (1 - eCB_level) -# * (1 - CDI_factor) -# * (1 - mGluR*alpha_mGluR) -# -> raw_influx ~= 0 -> no release. -# Auto-reset : NCX keeps clearing; CDI recovers once drive stops. -# ======================================================================= - -import numpy as np - -# ----------------------------------------------------------------------- -# PARAMETERS -# ----------------------------------------------------------------------- - -dt = 1.0 # ms - high-freq timestep -dt_slow = 1000.0 # ms - astrocyte / slow loop timestep -dt_meta = 60_000.0 # ms - metabolic loop timestep - -High_Freq_Multiplier = int(dt_slow / dt) # 1000 -Metabolic_Multiplier = int(dt_meta / dt) # 60000 - -# Unit-conversion scalar: biological rate constants in /s, timestep in ms. -# increment = k [/s] * signal * dt_s [s/step] -> dimensionless per step -dt_s = dt / 1000.0 # 0.001 s per 1 ms step - -# ----------------------------------------------------------------------- -# -- Voltage / membrane (NKA) -- -# ----------------------------------------------------------------------- -tau_V_pre = 2.0 # NKA: ms - membrane repolarisation time constant after AP - # The AP waveform decays with this time constant. - # Controls how long Ca2+ channels see a depolarised membrane. -V_pre_peak = 1.0 # NKA: normalised peak depolarisation on each AP (dimensionless) -V_rest = 0.0 # NKA: resting membrane potential (normalised) -V_pre_voltage = -10.0 # ORIG: driving force used in compute_flux (mV, kept for continuity) - -# Na/K-ATPase ATP cost -NKA_cost_per_AP = 0.002 # NKA: ATP units consumed per AP for Na/K-ATPase recharge. - # This is the largest single ATP cost per spike. - # [CASCADE 2] Accumulates in ATP_demand each time V_pre = 1; - # at 20 Hz that is 0.002 * 20 = 0.04 ATP units/s, - # dominating pump and docking costs at high rates. - -# ----------------------------------------------------------------------- -# -- Calcium influx & buffering -- -# ----------------------------------------------------------------------- -N_VGCC = 100 # ORIG: number of presynaptic VGCCs - # [CASCADE 6] Absolute ceiling of Ca2+ influx. -k_flux = 0.05 # ORIG: Ca2+ influx per open channel per unit driving force (a.u.) - -B_total = 1.0 # ORIG: total buffer capacity (normalised) - # [CASCADE 4] Saturates during bursting; B_free -> 0. -tau_buffer_rebind = 200.0 # NEW: ms - buffer recharge time constant - -# ----------------------------------------------------------------------- -# -- Ca2+ clearance rate constants (/ms, already per-millisecond) -- -# ----------------------------------------------------------------------- -# [CASCADE 3] All three define maximum clearance; PMCA and SERCA are ATP-gated. -k_PMCA = 0.03 # NEW: plasma membrane Ca-ATPase - primary, ATP-dependent -k_NCX = 0.10 # NEW: sodium-calcium exchanger - fast, NOT ATP-dependent -k_SERCA = 0.01 # NEW: ER pump - slowest, ATP-dependent -ATP_half = 0.3 # NEW: Hill half-saturation for ATP-dependent pumps - -# ATP cost of pumping (per unit Ca2+ cleared per ms) -ATP_cost_PMCA = 0.0005 # NKA: ATP units per unit Ca2+ extruded by PMCA per ms - # [CASCADE 2] Second-largest ongoing ATP drain after NKA. -ATP_cost_SERCA = 0.0002 # NKA: ATP units per unit Ca2+ pumped into ER per ms - -# ----------------------------------------------------------------------- -# -- Deterministic vesicle release (replaces stochastic_release) -- -# ----------------------------------------------------------------------- -# The Ca2+-sensor model (synaptotagmin cooperativity): -# -# release_fraction = k_rel * Ca_micro^n / (Ca_micro^n + KD_rel^n) -# * (1 - NT_suppression_weight * NT_cleft_norm) -# -# Term 1 - Hill equation: -# k_rel : maximum fraction of RRP releasable per spike (saturation ceiling) -# KD_rel : [Ca2+] at half-maximum release (half-saturation constant) -# n_rel : Hill cooperativity exponent (~4 for synaptotagmin-1) -# -# Term 2 - NT suppression modulation: -# NT_suppression_weight : how strongly accumulated cleft NT reduces the -# releasable fraction (autoreceptor-independent brake, -# represents physical occlusion of release sites and -# depletion-sensing). -# NT_cleft_norm : NT_cleft normalised to NT_suppression_sat (0 -> 1). -# -# The product is clamped to [0, 1] before multiplying by N_RRP, -# so released_NT is always a real non-negative number <= N_RRP. - -k_rel = 0.5 # DET: max releasable fraction of RRP per spike (0->1) - # Lower = more reluctant synapse at any [Ca2+]. -KD_rel = 1.0 # DET: half-saturation [Ca2+] (same a.u. as Ca_micro) - # Higher KD = release only at high [Ca2+] peaks. -n_rel = 4 # DET: Hill exponent - cooperativity of Ca2+ sensor - # n=4 matches synaptotagmin-1 (four C2 domain sites). - # Higher n = sharper threshold, more digital release. -NT_suppression_weight = 0.3 # DET: strength of NT_cleft feedback on release fraction - # 0 = no suppression; 1 = full block at saturation. -NT_suppression_sat = 50.0 # DET: NT_cleft level that saturates the suppression term - # (same units as NT_cleft vesicle count) -ATP_cost_docking = 0.001 # NKA: ATP units per vesicle docked (RP -> RRP per step) - # [CASCADE 2] Docking is ATP-dependent (NSF/SNAPs); - # adds to ATP_demand proportionally to refill_amount. - -# ----------------------------------------------------------------------- -# -- CDI (calcium-dependent inactivation) -- -# ----------------------------------------------------------------------- -k_CDI_rise = 0.8 # ORIG: CDI build rate (/s per unit Ca_micro) -Ca_micro_saturation = 2.0 # NEW: normalisation ceiling for CDI recovery -k_CDI_rec = 0.015 # NEW: CDI de-inactivation rate (/s) - # Both expressed in /s; applied with * dt_s in the loop. - -# ----------------------------------------------------------------------- -# -- Vesicle pools -- -# ----------------------------------------------------------------------- -Max_RRP = 20 # ORIG: [CASCADE 1] ceiling of the firing-ready pool -Max_RP = 200 # ORIG: [CASCADE 1] ceiling of the reserve pool - -# ----------------------------------------------------------------------- -# -- Trace integrator (RP->RRP recruitment speed) -- -# ----------------------------------------------------------------------- -tau_Tr_Ca = 1000.0 # ORIG: ms - calcium trace decay -T_high = 0.6 # ORIG: trace threshold -> fast recruitment -T_low = 0.2 # ORIG: trace threshold -> slow recruitment -k_rec_fast = 5.0 # /s — fast RP->RRP recruitment rate - # At high Tr_Ca: refills ~5% of headroom per second. - # Full RRP recovery from empty takes ~4-5 s of sustained activity. -k_rec_slow = 0.5 # /s — slow RP->RRP recruitment rate - # At low Tr_Ca: refills ~0.5% of headroom per second. - # Matches the ~30-60 s recovery seen after deep depletion. - -dt_slow_s = dt_slow / 1000.0 # 1.0 s — the slow loop timestep in seconds - -# ----------------------------------------------------------------------- -# -- Postsynaptic -- -# ----------------------------------------------------------------------- -tau_membrane = 20.0 # ORIG: ms -tau_desensitization = 500.0 # ORIG: ms - -# ----------------------------------------------------------------------- -# -- eCB retrograde brake -- -# ----------------------------------------------------------------------- -# [CASCADE 6] Slow retrograde suppressor of effective_conductance. -tau_eCB_rise = 2000.0 # ORIG: ms -tau_eCB_decay = 10_000.0 # ORIG: ms -eCB_threshold = 0.7 # ORIG: V_post level that triggers eCB synthesis - -# ----------------------------------------------------------------------- -# -- mGluR presynaptic autoreceptor -- -# ----------------------------------------------------------------------- -# [CASCADE 6] Fastest conductance brake; reads NT_cleft directly. -Km_mGluR = 0.5 # NEW: Michaelis-Menten half-saturation for NT_cleft -tau_mGluR = 2000.0 # NEW: ms -alpha_mGluR = 0.4 # NEW: max fractional VGCC suppression - -# ----------------------------------------------------------------------- -# -- Astrocyte / IP3 -- -# ----------------------------------------------------------------------- -tau_IP3 = 3000.0 # ORIG: ms -IP3_threshold = 0.8 # ORIG - -# ----------------------------------------------------------------------- -# -- Glutamine shuttle -- -# ----------------------------------------------------------------------- -conversion_efficiency_base = 0.8 # ORIG: fraction of Gln pool converted per cycle - -# ----------------------------------------------------------------------- -# -- NT cleft -- -# ----------------------------------------------------------------------- -tau_NT_decay = 5.0 # ms - passive NT diffusion / dilution out of cleft - - -# ----------------------------------------------------------------------- -# HELPER FUNCTIONS -# ----------------------------------------------------------------------- - -def compute_flux(conductance, voltage): - # ORIG: Ca2+ influx into microdomain. - # [CASCADE 4] Collapses to near zero once CASCADE 6 locks conductance. - return k_flux * conductance * abs(voltage) - - -def deterministic_release(N_RRP, Ca_micro, NT_cleft): - # DET: Deterministic Ca2+-driven vesicle release. - # - # Replaces stochastic_release (binomial with p = p_base * Ca_micro). - # - # Biology: synaptotagmin-1 is the fast Ca2+ sensor. Its four C2-domain - # Ca2+-binding sites give steep, cooperative Ca2+ sensitivity (Hill n~4). - # Release is not random per vesicle but is driven by the microdomain [Ca2+] - # that the sensor actually sees. At low [Ca2+] essentially no vesicles fuse; - # at high [Ca2+] a saturating fraction fuses within the AP window. - # - # Step 1: Hill equation gives Ca2+-dependent release fraction (0 -> k_rel). - Ca_n = Ca_micro ** n_rel - release_frac = k_rel * Ca_n / (Ca_n + KD_rel ** n_rel) - # - # Step 2: NT suppression modulation (0 -> 1, then inverted). - # Accumulated NT_cleft represents: (a) depletion sensing via presynaptic - # autoreceptors that are faster than the mGluR loop, and (b) physical - # competition for release site access at the active zone. - # This is distinct from the mGluR brake which reduces VGCC conductance; - # this term reduces the fraction of already-docked vesicles that fuse. - NT_norm = min(1.0, NT_cleft / NT_suppression_sat) - suppression = NT_suppression_weight * NT_norm - release_frac = release_frac * (1.0 - suppression) - release_frac = np.clip(release_frac, 0.0, 1.0) - # - # Step 3: Apply to pool size and floor at zero. - released = release_frac * N_RRP - return max(0.0, released) - - -def map_trace_to_speed(Tr_Ca): - # ORIG: Maps calcium trace to RP->RRP recruitment rate. - # [CASCADE 1] k_rec_fast lags release at high firing; collapse is - # self-accelerating as N_RP and headroom both shrink. - if Tr_Ca > T_high: - return k_rec_fast - elif Tr_Ca < T_low: - return k_rec_slow - else: - t = (Tr_Ca - T_low) / (T_high - T_low) - return k_rec_slow + t * (k_rec_fast - k_rec_slow) - - -def map_calcium_to_inactivation(Ca_micro): - # ORIG: Ca2+ drives CDI increment each ms. - # k_CDI_rise is /s; * dt_s gives dimensionless per-step increment. - # [CASCADE 5] Accumulates between spikes under pump failure. - return k_CDI_rise * Ca_micro * dt_s - - -def compute_pump_atp_factor(ATP_level): - # NEW: Hill function - ATP gates PMCA and SERCA speed. - # [CASCADE 3] Bridge from ATP_level (Loop 3) to clearance rate (Loop 1). - # ATP=1.0 -> ~0.92 (near full) | ATP=0.3 -> 0.50 | ATP=0.1 -> ~0.10 - return (ATP_level ** 2) / (ATP_level ** 2 + ATP_half ** 2) - - -def compute_EPSP(receptor_conductance): - # ORIG: postsynaptic potential increment. - return receptor_conductance * 0.1 - - -def compute_postsynaptic_eCB_signal(V_post_history): - # ORIG: eCB synthesis from sustained postsynaptic activity. - # [CASCADE 6] Slow retrograde brake; persists 10 s after burst ends. - recent_mean = (np.mean(V_post_history[-2000:]) - if len(V_post_history) >= 2000 - else np.mean(V_post_history)) - if recent_mean > eCB_threshold: - return recent_mean - eCB_threshold - return 0.0 - - -def compute_astrocyte_metabolic_health(Glucose_level, ATP_demand_accumulated): - # ORIG + NKA: Converts Glucose_level and accumulated ATP demand into: - # ATP_level -> [CASCADE 2->3 bridge] read every ms in Loop 1 - # conversion_efficiency -> [CASCADE 1] gates glutamine shuttle throughput - # - # NKA change: ATP_demand_accumulated (summed in Loop 1 each ms) is now - # subtracted from health before clamping, so high firing rates visibly - # reduce ATP_level within the same metabolic window. - # The demand term is normalised by a scale factor so that a physiological - # 20 Hz burst causes a ~20% ATP drop over one minute. - ATP_demand_scale = 50.0 # normalisation: demand at 20 Hz over 60 s ~ 1.0 - health = np.clip(Glucose_level - ATP_demand_accumulated / ATP_demand_scale, - 0.0, 1.0) - return health, health # (conversion_efficiency, ATP_level) - - -def trigger_slow_astrocyte_calcium_wave(): - # ORIG: placeholder - would release gliotransmitters over ~10 s. - pass - - -# ----------------------------------------------------------------------- -# STATE VARIABLES (initial values) -# ----------------------------------------------------------------------- - -# -- Voltage / membrane -- -V_pre_state = 0.0 # NKA: continuous membrane potential (0=rest, 1=peak AP) - # Decays with tau_V_pre after each spike. - # Controls the effective driving window for Ca2+ influx. - -# -- Presynaptic Ca2+ -- -Ca_micro = 0.0 # ORIG: free cytosolic [Ca2+] in microdomain [CASCADE 4] -Ca_ER = 0.5 # NEW: Ca2+ stored in ER -Ca_buffer_bound = 0.0 # NEW: Ca2+ bound to buffer proteins -B_free = 1.0 # NEW: free buffer sites [CASCADE 4] - -# -- CDI -- -CDI_factor = 0.0 # ORIG [CASCADE 5,6] - -# -- Vesicle pools -- -N_RRP = 15.0 # ORIG: readily-releasable pool [CASCADE 1] (float for deterministic) -N_RP = 150.0 # ORIG: reserve pool [CASCADE 1] - -# -- Calcium trace -- -Tr_Ca = 0.0 # ORIG: integrative Ca2+ memory - -# -- NT in cleft -- -NT_cleft = 0.0 # ORIG [CASCADE 6] - -# -- Postsynaptic -- -V_post = 0.0 -receptor_conductance = 0.0 -Desensitization_level = 0.0 -V_post_history = [] - -# -- Retrograde / autoreceptor -- -eCB_level = 0.0 # ORIG [CASCADE 6] -mGluR_activation = 0.0 # NEW [CASCADE 6] - -# -- Astrocyte -- -IP3 = 0.0 -Glutamine_pool = 50.0 - -# -- Metabolic -- -ATP_level = 1.0 # NEW [CASCADE 2->3] -ATP_demand = 0.0 # NKA: accumulated ATP demand within current metabolic window -conversion_efficiency = 0.8 # ORIG -Glucose_level = 1.0 # ORIG: set < 1.0 to engage metabolic silencing cascade - - -# ----------------------------------------------------------------------- -# MAIN SIMULATION LOOP -# ----------------------------------------------------------------------- - -def run_simulation(spike_train, total_steps): - """ - spike_train : list of int timestep indices at which an AP arrives - total_steps : int number of 1 ms steps to simulate - """ - - global V_pre_state - global Ca_micro, Ca_ER, Ca_buffer_bound, B_free - global CDI_factor - global N_RRP, N_RP - global Tr_Ca, NT_cleft - global V_post, receptor_conductance, Desensitization_level, V_post_history - global eCB_level, mGluR_activation - global IP3, Glutamine_pool - global ATP_level, ATP_demand, conversion_efficiency, Glucose_level - - log = {k: [] for k in [ - "V_pre_state", - "Ca_micro", "Ca_ER", "CDI_factor", "B_free", - "N_RRP", "N_RP", "Tr_Ca", "NT_cleft", - "V_post", "eCB_level", "mGluR_activation", - "released_NT", "ATP_level", "ATP_demand", - ]} - - spike_set = set(spike_train) - - for step in range(total_steps): - - # ============================================================== - # LOOP 1 - HIGH-FREQUENCY (dt = 1 ms) - # ============================================================== - - V_pre = 1 if step in spike_set else 0 - released_NT = 0.0 - - # -- 1A. MEMBRANE VOLTAGE / Na-K-ATPase ------------------------ - # - # NKA: V_pre_state is now a continuous variable. - # On each AP it jumps to V_pre_peak, then decays exponentially - # toward V_rest with time constant tau_V_pre (~2 ms). - # This models the width of the depolarisation window that keeps - # VGCCs open after the peak, giving Ca2+ influx a temporal profile - # rather than a single instantaneous pulse. - if V_pre == 1: - V_pre_state = V_pre_peak # AP fires: membrane jumps to peak - - # Exponential decay toward rest (Na/K-ATPase restores the gradient) - V_pre_state += (V_rest - V_pre_state) * dt / tau_V_pre - - # NKA: ATP cost of Na/K-ATPase recharge on each AP. - # Fired once per spike, not per ms — the cost is per action potential. - # [CASCADE 2] This is the dominant ATP drain at high firing rates. - if V_pre == 1: - ATP_demand += NKA_cost_per_AP - - # -- 1B. PRESYNAPTIC Ca2+ INFLUX & VESICLE RELEASE ------------- - # - # Ca2+ influx now uses V_pre_state (the continuous voltage) instead - # of the binary V_pre flag, so influx has the same temporal profile - # as the depolarisation window and tapers as the membrane repolarises. - # [CASCADE 6] OUTCOME: effective_conductance collapses when any of - # the three suppression terms approaches 1. - effective_conductance = ( - N_VGCC - * (1.0 - eCB_level) # [CASCADE 6] retrograde brake - * (1.0 - CDI_factor) # [CASCADE 5->6] CDI lock-out - * (1.0 - mGluR_activation * alpha_mGluR) # [CASCADE 6] autoreceptor brake - ) - - # Ca2+ influx is gated by V_pre_state: significant only while the - # membrane is depolarised; tapers to zero as V_pre_state -> V_rest. - raw_influx = compute_flux(effective_conductance, V_pre_state) - - # NEW: Buffer proteins capture a fraction of raw_influx immediately. - # [CASCADE 4] B_free -> 0 during sustained bursting (saturation). - capture_fraction = B_free / B_total - captured = raw_influx * capture_fraction - B_free = max(0.0, B_free - captured) - Ca_buffer_bound += captured - Ca_micro += (raw_influx - captured) - - # DET: Deterministic Ca2+-driven vesicle release. - # Released on every ms that Ca_micro > 0 (not only on the spike flag), - # giving a smooth release profile that follows the Ca2+ transient. - # This is more physically accurate than gating release on V_pre == 1: - # in biology, vesicles fuse throughout the Ca2+ microdomain lifetime - # (~1-2 ms), not only at the exact moment of depolarisation. - # [CASCADE 1] released_NT draws from N_RRP; the Ca2+ and NT_cleft - # dependence means release self-limits as both Ca_micro - # falls (clearance) and NT_cleft rises (suppression). - if N_RRP > 0 and Ca_micro > 0: - released_NT = deterministic_release(N_RRP, Ca_micro, NT_cleft) - released_NT = min(released_NT, N_RP) - N_RRP -= released_NT - NT_cleft += released_NT - NT_released_this_window += released_NT # NEW: accumulator for Loop 2 - - # Passive diffusion — this is fine as-is, represents lateral escape - NT_cleft *= (1.0 - dt / tau_NT_decay) - NT_cleft = max(0.0, NT_cleft) - - # -- 1C. Ca2+ CLEARANCE (every ms) ---------------------------- - # - # [CASCADE 3] pump_scale: ATP bridge from Loop 3 to clearance. - pump_scale = compute_pump_atp_factor(ATP_level) - - # [CASCADE 3->4] PMCA: primary, ATP-dependent. - cleared_PMCA = k_PMCA * Ca_micro * pump_scale - # [CASCADE 3 note] NCX: fast, NOT ATP-dependent. Floor, not rescue. - cleared_NCX = k_NCX * Ca_micro - # [CASCADE 3->4] SERCA: slowest, ATP-dependent. - cleared_SERCA = k_SERCA * Ca_micro * pump_scale - - Ca_micro -= (cleared_PMCA + cleared_NCX + cleared_SERCA) - Ca_micro = max(0.0, Ca_micro) - - Ca_ER += cleared_SERCA # ER store loaded while ATP is available - - # NKA: ATP cost of PMCA and SERCA clearing Ca2+ this step. - # [CASCADE 2] Ongoing drain proportional to Ca2+ load; highest - # during the early burst when Ca_micro peaks. - ATP_demand += ATP_cost_PMCA * cleared_PMCA - ATP_demand += ATP_cost_SERCA * cleared_SERCA - - # NEW: Buffer recharge - bound Ca2+ slowly re-releases to cytosol. - # [CASCADE 4] Sustains Ca_micro elevation during pump failure. - rebind_flux = Ca_buffer_bound * dt / tau_buffer_rebind - Ca_micro += rebind_flux - Ca_buffer_bound = max(0.0, Ca_buffer_bound - rebind_flux) - B_free = B_total - Ca_buffer_bound - - # -- 1D. CDI - RISE AND RECOVERY -------------------------------- - # - # [CASCADE 5] Rise: proportional to Ca_micro, fires every ms. - CDI_factor += map_calcium_to_inactivation(Ca_micro) - # [CASCADE 5] Recovery: suppressed when Ca_micro is high. - # Self-locking: pump failure -> Ca_micro high -> - # recovery_rate -> 0 -> CDI_factor -> 1 -> silence. - CDI_recovery_rate = k_CDI_rec * (1.0 - Ca_micro / Ca_micro_saturation) - CDI_factor = np.clip(CDI_factor - CDI_recovery_rate * dt_s, 0.0, 1.0) - - # -- 1E. TRACE INTEGRATOR ------------------------------------- - # ORIG: integrates Ca_micro; drives RP->RRP recruitment speed. - Tr_Ca = Tr_Ca + (Ca_micro - Tr_Ca / tau_Tr_Ca) * dt - - # -- 1G. POSTSYNAPTIC FAST RESPONSE --------------------------- - # ORIG: receptor activation and desensitization. - effective_NT = released_NT * (1.0 - Desensitization_level) - receptor_conductance += effective_NT * 0.05 - receptor_conductance *= (1.0 - dt / tau_membrane) - - V_post += compute_EPSP(receptor_conductance) - (V_post / tau_membrane) * dt - V_post = max(0.0, V_post) - - Desensitization_level += NT_cleft * 0.001 * dt - Desensitization_level -= (Desensitization_level / tau_desensitization) * dt - Desensitization_level = np.clip(Desensitization_level, 0.0, 1.0) - - # ORIG: NT diffuses / dilutes out of cleft each ms. - NT_cleft *= (1.0 - dt / tau_NT_decay) - NT_cleft = max(0.0, NT_cleft) - - V_post_history.append(V_post) - if len(V_post_history) > 5000: - V_post_history.pop(0) - - # -- RECORD --------------------------------------------------- - log["V_pre_state"].append(V_pre_state) - log["Ca_micro"].append(Ca_micro) - log["Ca_ER"].append(Ca_ER) - log["CDI_factor"].append(CDI_factor) - log["B_free"].append(B_free) - log["N_RRP"].append(N_RRP) - log["N_RP"].append(N_RP) - log["Tr_Ca"].append(Tr_Ca) - log["NT_cleft"].append(NT_cleft) - log["V_post"].append(V_post) - log["eCB_level"].append(eCB_level) - log["mGluR_activation"].append(mGluR_activation) - log["released_NT"].append(released_NT) - log["ATP_level"].append(ATP_level) - log["ATP_demand"].append(ATP_demand) - - # ============================================================== - # LOOP 2 - SLOW / ASTROCYTE (dt_slow = 1 s) - # ============================================================== - - if (step % High_Freq_Multiplier) == 0: - - # EAAT clearance: a fraction of what remains in the cleft - # (this is now meaningful because diffusion hasn't zeroed it yet - # at realistic tau_NT_decay values — but we also fix IP3 sourcing) - cleared_NT = NT_cleft * 0.3 - NT_cleft = max(0.0, NT_cleft - cleared_NT) - - # IP3 integrates total release load, not post-diffusion residual. - # This correctly represents the astrocyte sensing cumulative activity. - IP3 += NT_released_this_window - (IP3 / tau_IP3) * dt_slow - IP3 = max(0.0, IP3) - NT_released_this_window = 0.0 # reset accumulator for next window - - if IP3 > IP3_threshold: - trigger_slow_astrocyte_calcium_wave() - - # ORIG: eCB retrograde signal. - # [CASCADE 6] Second conductance brake (~2 s onset, 10 s decay). - eCB_signal = compute_postsynaptic_eCB_signal(V_post_history) - if eCB_signal > 0: - eCB_level += eCB_signal * (dt_slow / tau_eCB_rise) - else: - eCB_level -= eCB_level * (dt_slow / tau_eCB_decay) - eCB_level = np.clip(eCB_level, 0.0, 1.0) - - - - - # ── LOOP 2 — SLOW / ASTROCYTE (dt_slow = 1 s) ────────────────────── - - # -- RP -> RRP RECRUITMENT (with pool guards) -------------------- - # Runs once per second. k_rec_fast and k_rec_slow are in /s, - # so multiplying by dt_slow_s = 1.0 s gives a dimensionless - # per-step fraction — no hidden unit scaling needed. - # - # [CASCADE 1] Recruitment is the only counter-force to depletion. - # Even k_rec_fast fills only ~5% of headroom per second, - # lagging well behind release at high firing rates. - current_recruitment_rate = map_trace_to_speed(Tr_Ca) # /s - refill_amount = current_recruitment_rate * dt_slow_s * N_RP * (Max_RRP - N_RRP) / Max_RRP - refill_amount = max(0.0, refill_amount) - refill_amount = min(refill_amount, N_RP) - - N_RRP = min(N_RRP + refill_amount, Max_RRP) - N_RP = max(0.0, N_RP - refill_amount) - - # ATP cost of docking moves here too — it is per recruitment event, - # not per millisecond. - ATP_demand += ATP_cost_docking * refill_amount - - # ============================================================== - # LOOP 3 - METABOLIC (dt_meta = 1 min) - # ============================================================== - - if (step % Metabolic_Multiplier) == 0: - - # NKA + ORIG: ATP_level is now driven by both Glucose_level - # (supply) and ATP_demand (demand accumulated over this window). - # [CASCADE 2] ATP_demand encodes all three cost sources: - # NKA_cost_per_AP - dominant at high firing rates - # ATP_cost_PMCA/SERCA - proportional to Ca2+ load - # ATP_cost_docking - proportional to recruitment rate - conversion_efficiency, ATP_level = compute_astrocyte_metabolic_health( - Glucose_level, ATP_demand - ) - - # Reset demand accumulator for the next metabolic window. - ATP_demand = 0.0 - - # [CASCADE 1 - slow refill] Glutamine shuttle rebuilds N_RP. - refill_RP = Glutamine_pool * conversion_efficiency - N_RP = min(Max_RP, N_RP + refill_RP) - Glutamine_pool = max(0.0, Glutamine_pool - refill_RP) - - return log - - -# ----------------------------------------------------------------------- -# EXAMPLE USAGE -# ----------------------------------------------------------------------- - -if __name__ == "__main__": - import matplotlib.pyplot as plt - - total_steps = 10_000 # 10 seconds of simulated time - - # 20 Hz burst for 2 s, then silence. - # Set Glucose_level = 0.2 to engage the full metabolic cascade. - spike_train = list(range(0, 2000, 50)) - - results = run_simulation(spike_train, total_steps) - - t = np.arange(total_steps) * dt - - fig, axes = plt.subplots(7, 1, figsize=(12, 16), sharex=True) - fig.suptitle("Tripartite Synapse - Presynaptic Model", fontsize=13) - - axes[0].plot(t, results["V_pre_state"], color="slateblue", lw=0.8) - axes[0].set_ylabel("V_pre (a.u.)") - axes[0].set_title("NKA - membrane repolarises with tau_V_pre after each AP", - fontsize=9, loc="left") - - axes[1].plot(t, results["Ca_micro"], color="darkorange", lw=0.8) - axes[1].set_ylabel("[Ca2+] free (a.u.)") - axes[1].set_title("CASCADE 4 - residual Ca2+ profile follows V_pre decay", - fontsize=9, loc="left") - - axes[2].plot(t, results["CDI_factor"], color="firebrick", lw=0.8, label="CDI factor") - axes[2].plot(t, results["B_free"], color="steelblue", lw=0.8, label="Buffer free") - axes[2].set_ylabel("CDI / Buffer (0-1)") - axes[2].set_title("CASCADE 5 - CDI lock-out | CASCADE 4 - buffer saturation", - fontsize=9, loc="left") - axes[2].legend(fontsize=8) - - axes[3].plot(t, results["N_RRP"], color="teal", lw=0.8, label="RRP") - axes[3].plot(t, results["N_RP"], color="purple", lw=0.8, label="RP") - axes[3].set_ylabel("Vesicles") - axes[3].set_title("CASCADE 1 - deterministic depletion (Hill Ca2+ sensor + NT suppression)", - fontsize=9, loc="left") - axes[3].legend(fontsize=8) - - axes[4].plot(t, results["NT_cleft"], color="darkgreen", lw=0.8, label="NT cleft") - axes[4].plot(t, results["mGluR_activation"], color="saddlebrown", lw=0.8, label="mGluR") - axes[4].plot(t, results["eCB_level"], color="crimson", lw=0.8, label="eCB") - axes[4].set_ylabel("Cleft / Feedback") - axes[4].set_title("CASCADE 6 - three multiplicative brakes on effective_conductance", - fontsize=9, loc="left") - axes[4].legend(fontsize=8) - - axes[5].plot(t, results["V_post"], color="navy", lw=0.8) - axes[5].set_ylabel("V_post (a.u.)") - axes[5].set_title("CASCADE 6 result - postsynaptic silence", fontsize=9, loc="left") - - axes[6].plot(t, results["ATP_level"], color="goldenrod", lw=0.8, label="ATP level") - axes[6].plot(t, results["ATP_demand"], color="tomato", lw=0.8, label="ATP demand (cumul.)") - axes[6].set_ylabel("ATP (a.u.)") - axes[6].set_title("CASCADE 2 - NKA + pump + docking demand drives ATP depletion", - fontsize=9, loc="left") - axes[6].set_xlabel("Time (ms)") - axes[6].legend(fontsize=8) - - plt.tight_layout() - plt.savefig("./synapse_simulation.png", dpi=150) - plt.close() - print("Done.") \ No newline at end of file diff --git a/elements/neuron/appunti/traditional_python_models/tripartite-soma.py b/elements/neuron/appunti/traditional_python_models/tripartite-soma.py deleted file mode 100644 index 450715e..0000000 --- a/elements/neuron/appunti/traditional_python_models/tripartite-soma.py +++ /dev/null @@ -1,883 +0,0 @@ -# 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