#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Neuronal Membrane Circadian Oscillator Model This code implements the neuronal membrane-associated circadian oscillator described in: "Autonomous circadian oscillators intrinsic to cell membranes" Forlino et al., 2026 The model couples two timescales: 1. Fast timescale (milliseconds): Hodgkin-Huxley type action potential generation 2. Slow timescale (hours): Membrane-based circadian oscillator The membrane oscillator operates through: - Clock Channels (CC) mediate Ca2+ influx - Ca2+ activates cAMP production (positive feedback) - Ca2+ produces inhibitor Y (delayed negative feedback) - Inhibitor Y inactivates CCs, closing the feedback loop - This generates circadian modulation of neuronal excitability and firing rate The simulation runs for 4 days and outputs time series data for analysis. """ import numpy as np from scipy.integrate import solve_ivp from scipy.signal import find_peaks import os np.seterr(all='ignore') # Final packed dataset DATA_FILE = 'neuron_data.npz' # Prefix for the incremental (crash-safe) plain-text dumps written during the # run. They are packed into DATA_FILE at the end and then deleted; if the run # is interrupted, the partial data survives in these files. RAW_PREFIX = '_raw_' # Spike times are stored as integer multiples of this quantum (in ms), as a # first time plus uint32 inter-spike intervals. Quantising the *absolute* times # before differencing means the reconstruction has no cumulative drift. SPIKE_QUANTUM_MS = 1e-3 # ============================================================================ # SLOW SYSTEM PARAMETERS (Membrane Circadian Oscillator) # ============================================================================ # Clock Channel (CC) dynamics - units: h^-1 alpha = 0.7 # CC recruitment/activation rate beta = 0.9 # CC inactivation rate (modulated by inhibitor Y) # cAMP dynamics (ligand that activates CCs) - units: h^-1 alpha_cAMP = 2.87 # cAMP production rate (Ca2+-dependent) gamma_cAMP = 2.87 # cAMP degradation rate # Inhibitor Y dynamics (provides delayed negative feedback) - units: h^-1 alpha_Y = 0.2 # Inhibitor production rate (Ca2+-dependent) gamma_Y = 0.2 # Inhibitor degradation rate # Clock Channel properties gs = 0.9 # Single channel conductance (mS/cm²) Ncc_Tot = 2.22 # Total number of clock channels (a.u.) # ============================================================================ # FAST SYSTEM PARAMETERS (Hodgkin-Huxley Neuronal Dynamics) # ============================================================================ # --- Reversal Potentials (mV) --- vNa = 45 # Sodium reversal potential vK = -105 # Potassium reversal potential vCa = 120 # Calcium reversal potential (HVA channels) vcc = 10 # Clock channel reversal potential vL = -62.5 # Leak reversal potential # --- Maximum Conductances (mS/cm²) --- # Spiking currents gNa = 29 # Transient sodium (action potential upstroke) gK = 13 # Delayed rectifier potassium (action potential repolarization) gNaP = 20 # Persistent sodium (spontaneous firing) gHVA = 0.2 # High-voltage-activated calcium gL = 3 # Leak current # --- Gating Variable Parameters --- # Transient Sodium (Na) - activation only, inactivation via h = 1-nK theta_mNa = -25 # Half-activation voltage (mV) beta_mNa = -6.5 # Slope factor (mV) # Delayed Rectifier Potassium (K) theta_nK = -26 # Half-activation voltage (mV) beta_nK = -9 # Slope factor (mV) tau_nK = 10 # Base time constant (ms) # Persistent Sodium (NaP) theta_mNaP = -40 # Activation half-voltage (mV) beta_mNaP = -4 # Activation slope (mV) theta_hNaP = -54 # Inactivation half-voltage (mV) beta_hNaP = 5 # Inactivation slope (mV) tau_hNaP = 500 # Base inactivation time constant (ms) # High-Voltage-Activated Calcium (HVA) theta_mHVA = -10.0 # Half-activation voltage (mV) beta_mHVA = -6.5 # Slope factor (mV) # --- Other Electrophysiological Parameters --- C = 2 # Membrane capacitance (μF/cm²) f = 0.0039 # Current-to-concentration conversion factor (cm²/nC) tau = 40 # Ca2+ clearance time constant (ms) Ca0 = 0.37 # Baseline Ca2+ concentration (a.u.) # ============================================================================ # GATING FUNCTIONS (Steady-state activation/inactivation) # ============================================================================ def nK_inf(V): """Potassium channel activation (steady-state).""" return 1 / (1 + np.exp((V - theta_nK) / beta_nK)) def mNa_inf(V): """Sodium channel activation (steady-state).""" return 1 / (1 + np.exp((V - theta_mNa) / beta_mNa)) def mNaP_inf(V): """Persistent sodium activation (steady-state).""" return 1 / (1 + np.exp((V - theta_mNaP) / beta_mNaP)) def hNaP_inf(V): """Persistent sodium inactivation (steady-state).""" return 1 / (1 + np.exp((V - theta_hNaP) / beta_hNaP)) def mHVA_inf(V): """High-voltage-activated calcium activation (steady-state).""" return 1 / (1 + np.exp((V - theta_mHVA) / beta_mHVA)) def Hill(cAMP): """ Hill function for cAMP-dependent CC activation. Parameters ---------- cAMP : float Intracellular cAMP concentration (a.u.) Returns ------- float CC open probability (0 to 1) Notes ----- Uses Hill coefficient n=3 and K_1/2 = 1.47 This introduces the nonlinearity necessary for oscillations. """ return 1 / (1 + (1.47 / cAMP)**3) # ============================================================================ # TIME CONSTANT FUNCTIONS # ============================================================================ def nK_tau(V): """Voltage-dependent time constant for potassium activation.""" return tau_nK / np.cosh((V - theta_nK) / (2.0 * beta_nK)) def hNaP_tau(V): """Voltage-dependent time constant for persistent sodium inactivation.""" return tau_hNaP / np.cosh((V - theta_hNaP) / (2.0 * beta_hNaP)) # ============================================================================ # CONVERT RATE CONSTANTS FROM HOURS TO MILLISECONDS # ============================================================================ # The slow oscillator operates on hour timescale, but ODEs are solved in ms alpha_ms = alpha / 3600000 beta_ms = beta / 3600000 alpha_cAMP_ms = alpha_cAMP / 3600000 gamma_cAMP_ms = gamma_cAMP / 3600000 alpha_Y_ms = alpha_Y / 3600000 gamma_Y_ms = gamma_Y / 3600000 # ============================================================================ # COUPLED SYSTEM OF DIFFERENTIAL EQUATIONS # ============================================================================ def compute_derivatives(t0, y): """ Compute time derivatives for the coupled neuron-oscillator system. This function integrates the fast electrophysiological dynamics (millisecond timescale) with the slow membrane oscillator (hour timescale). Parameters ---------- y : array_like State vector containing: y[0] : Ncc - Number of active clock channels (a.u.) y[1] : Y - Inhibitor concentration (a.u.) y[2] : cAMP - Cyclic AMP concentration (a.u.) y[3] : V - Membrane potential (mV) y[4] : Ca - Intracellular calcium concentration (a.u.) y[5] : nK - Potassium channel activation (0-1) y[6] : hNaP - Persistent sodium inactivation (0-1) t0 : float Current time (ms) - not explicitly used but required by integrator Returns ------- dy : ndarray Time derivatives of all state variables Notes ----- The membrane oscillator (Ncc, Y, cAMP) operates on hour timescale. The electrical dynamics (V, Ca, nK, hNaP) operate on millisecond timescale. Ca2+ couples both timescales by modulating both spiking and oscillator dynamics. """ dy = np.zeros((7,)) # Extract state variables Ncc = y[0] # Active clock channels Y = y[1] # Inhibitor cAMP = y[2] # Cyclic AMP V = y[3] # Membrane potential Ca = y[4] # Calcium concentration nK = y[5] # Potassium activation hNaP = y[6] # Persistent sodium inactivation # ======================================================================== # COMPUTE IONIC CURRENTS # ======================================================================== # Clock Channel current (modulated by cAMP) n = Hill(cAMP) Icc = Ncc * gs * n * (V - vcc) # Leak current IL = gL * (V - vL) # Transient sodium current (fast Na+ for action potentials) # Uses approximation: hNa = 1 - nK (inactivation tied to K activation) INa = gNa * (1 - nK) * (mNa_inf(V)**3) * (V - vNa) # Delayed rectifier potassium current IK = gK * (nK**4) * (V - vK) # Persistent sodium current (supports spontaneous firing) INaP = gNaP * mNaP_inf(V) * hNaP * (V - vNa) # High-voltage-activated calcium current IHVA = gHVA * mHVA_inf(V) * (V - vCa) # ======================================================================== # DIFFERENTIAL EQUATIONS - SLOW SYSTEM (Membrane Oscillator) # ======================================================================== # Clock Channel dynamics # Channels are recruited at rate alpha, inactivated by inhibitor Y dy[0] = alpha_ms * (Ncc_Tot - Ncc) - beta_ms * Ncc * Y # Inhibitor dynamics (delayed negative feedback) # Produced by Ca2+, provides delayed suppression of CC activity dy[1] = Ca * alpha_Y_ms - Y * gamma_Y_ms # cAMP dynamics (positive feedback modulator) # Produced by Ca2+, activates CCs via Hill function dy[2] = Ca * alpha_cAMP_ms - cAMP * gamma_cAMP_ms # ======================================================================== # DIFFERENTIAL EQUATIONS - FAST SYSTEM (Electrophysiology) # ======================================================================== # Membrane potential (current balance equation) dy[3] = -(Icc + IL + INa + IK + INaP + IHVA) / C # Calcium dynamics # Influx through CC (fraction 0.3) and HVA channels, with clearance dy[4] = -f * (0.3 * Icc + IHVA) + (Ca0 - Ca) / tau # Potassium activation (first-order kinetics) dy[5] = (nK_inf(V) - nK) / nK_tau(V) # Persistent sodium inactivation (first-order kinetics) dy[6] = (hNaP_inf(V) - hNaP) / hNaP_tau(V) return dy # ============================================================================ # INITIAL CONDITIONS # ============================================================================ # These initial conditions represent a point on the limit cycle init_cond = np.array([0.815165095152994, 1.40770110623839, 1.90609996371345, -5.09219970e+01, 1.84736123e+00, 5.85939821e-02, 3.53005281e-02]) # ============================================================================ # MAIN SIMULATION LOOP - 4 DAYS, 1-HOUR CHUNKS # # Each solve_ivp call covers 1 hour instead of 24, reducing peak memory # from ~1.4 GB/day to ~60 MB/chunk. Spike detection runs on each chunk's # full-resolution V trace before the solution is discarded. Slow variables # are interpolated onto a regular 60-second grid per chunk. # ============================================================================ CHUNK_HOURS = 1 CHUNK_MS = CHUNK_HOURS * 3_600_000 # milliseconds per chunk SLOW_GRID_MS = 60_000 # slow-variable output resolution: 60 s N_CHUNKS = 4 * 24 // CHUNK_HOURS # total chunks (4 days) LAST_DAY_START_CHUNK = N_CHUNKS - 24 // CHUNK_HOURS # index of first chunk of last day print("Starting simulation of neuronal membrane oscillator...") print(f"Total duration: {N_CHUNKS * CHUNK_HOURS} h | " f"Chunk size: {CHUNK_HOURS} h | " f"Slow-variable resolution: {SLOW_GRID_MS//1000} s") print("=" * 60) # ----------------------------------------------------------------------- # Output files — one per variable, written incrementally after each chunk # ----------------------------------------------------------------------- OUTPUT_FILES = ['time.txt', 'Ncc.txt', 'Y.txt', 'cAMP.txt', 'spike_train.txt'] for fname in OUTPUT_FILES: open(RAW_PREFIX + fname, 'w').close() # clear / create def append_array(filename, arr): """Append a 1-D array to a scratch text file, one value per line.""" with open(RAW_PREFIX + filename, 'ab') as f: np.savetxt(f, arr) # last_day_chunk_inits[k]: initial condition at the start of chunk # (LAST_DAY_START_CHUNK + k) — needed to re-run a specific chunk for # the 1-second window extraction after the main loop. last_day_chunk_inits = {} for chunk_idx in range(N_CHUNKS): t0 = chunk_idx * CHUNK_MS t1 = t0 + CHUNK_MS if chunk_idx >= LAST_DAY_START_CHUNK: last_day_chunk_inits[chunk_idx - LAST_DAY_START_CHUNK] = init_cond.copy() # ------------------------------------------------------------------ # Solve 1-hour chunk # ------------------------------------------------------------------ chunk = solve_ivp(compute_derivatives, [t0, t1], init_cond) # Spike detection on the full-resolution voltage trace peaks = find_peaks(chunk.y[3], height=-20)[0] chunk_spikes = chunk.t[peaks] # Slow variables on a regular 60-second grid — append to files t_grid = np.arange(t0, t1, SLOW_GRID_MS) append_array('time.txt', t_grid) append_array('Ncc.txt', np.interp(t_grid, chunk.t, chunk.y[0])) append_array('Y.txt', np.interp(t_grid, chunk.t, chunk.y[1])) append_array('cAMP.txt', np.interp(t_grid, chunk.t, chunk.y[2])) append_array('spike_train.txt', chunk_spikes) # Advance and free init_cond = chunk.y[:, -1] del chunk print(f" chunk {chunk_idx+1:3d}/{N_CHUNKS} " f"(t = {t0/3_600_000:.0f}–{t1/3_600_000:.0f} h)", end='\r') print("\n" + "=" * 60) print(f"Simulation complete! ({N_CHUNKS * CHUNK_HOURS} h total)") # ============================================================================ # EXTRACT 1-SECOND WINDOWS FOR DETAILED ANALYSIS (Figures 2E and 2F) # # Load the last day's spike times from spike_train.txt, compute firing # rate, then re-run only the relevant 1-hour chunk with t_eval. # ============================================================================ print("\nExtracting 1-second time windows at min/max firing rates...") # Load all spike times and isolate the last day all_spikes = np.loadtxt(RAW_PREFIX + 'spike_train.txt') last_day_t_start = LAST_DAY_START_CHUNK * CHUNK_MS last_day_spikes = all_spikes[all_spikes >= last_day_t_start] ISI_last = np.diff(last_day_spikes) last_day_frequency = 1000 / ISI_last # Hz argmin = np.argmin(last_day_frequency) argmax = np.argmax(last_day_frequency) Tmin = last_day_spikes[argmin + 1] Tmax = last_day_spikes[argmax + 1] def extract_window(T_window_start): """ Re-run the 1-hour chunk containing T_window_start and return a 1-second window of V and Ca starting at T_window_start. No t_eval: the adaptive solver keeps its internal sub-ms steps, which is essential for resolving action potential peaks accurately. Only 1 hour is re-run, so memory is not a concern. """ k = min(int((T_window_start - last_day_t_start) / CHUNK_MS), 24 // CHUNK_HOURS - 1) t_chunk_start = last_day_t_start + k * CHUNK_MS t_window_end = T_window_start + 1000 # 1 second = 1000 ms sol = solve_ivp( compute_derivatives, [t_chunk_start, t_window_end], last_day_chunk_inits[k], dense_output=False, ) # Trim to the 1-second window of interest mask = sol.t >= T_window_start t_rel = sol.t[mask] - T_window_start return t_rel, sol.y[3][mask], sol.y[4][mask] # relative time, V, Ca T_series_min, V_series_min, Ca_series_min = extract_window(Tmin) T_series_max, V_series_max, Ca_series_max = extract_window(Tmax) print(f" Min-rate window: ZT {Tmin/3_600_000:.2f} h") print(f" Max-rate window: ZT {Tmax/3_600_000:.2f} h") # ============================================================================ # PACK EVERYTHING INTO A SINGLE COMPRESSED ARCHIVE # ============================================================================ print(f"\nPacking results into {DATA_FILE} ...") # Spike times: quantise the absolute times first, then store the differences. # Differencing quantised times (rather than quantising the intervals) keeps the # reconstruction exact to one quantum, with no drift accumulating over the run. spike_q = np.round(all_spikes / SPIKE_QUANTUM_MS).astype(np.int64) spike_t0 = spike_q[0] spike_isi = np.diff(spike_q) assert spike_isi.min() > 0, "spike train is not strictly increasing" assert spike_isi.max() < np.iinfo(np.uint32).max, "ISI overflows uint32" spike_isi = spike_isi.astype(np.uint32) np.savez_compressed( DATA_FILE, # Slow variables, on the regular SLOW_GRID_MS grid (ms) time=np.loadtxt(RAW_PREFIX + 'time.txt'), Ncc=np.loadtxt(RAW_PREFIX + 'Ncc.txt'), Y=np.loadtxt(RAW_PREFIX + 'Y.txt'), cAMP=np.loadtxt(RAW_PREFIX + 'cAMP.txt'), # Spike train, quantised (see reconstruct_spike_train in NEURON_plot.py) spike_t0=spike_t0, spike_isi=spike_isi, spike_quantum_ms=SPIKE_QUANTUM_MS, # 1-second windows at the daily firing-rate extremes T_series_min=T_series_min, V_series_min=V_series_min, Ca_series_min=Ca_series_min, T_series_max=T_series_max, V_series_max=V_series_max, Ca_series_max=Ca_series_max, ) # Verify the round-trip before discarding the scratch files _check = (spike_t0 + np.concatenate([[0], np.cumsum(spike_isi, dtype=np.int64)])) \ * SPIKE_QUANTUM_MS _err = np.max(np.abs(_check - all_spikes)) print(f" Spike-train round-trip error: {_err:.2e} ms " f"({len(all_spikes)} spikes, {os.path.getsize(DATA_FILE)/1e6:.2f} MB)") # Remove the scratch files now that the packed archive is written for fname in OUTPUT_FILES: os.remove(RAW_PREFIX + fname) print("\nAll data saved successfully!") print("=" * 60) # ============================================================================ # END OF SIMULATION # ============================================================================