STRC AC1-CREB Phase 3 Bifurcation Analysis

Linear-stability analysis of the 7-variable AC1-CREB ODE across the 20–110 dB SPL range, extending STRC AC1-CREB Parameter Robustness. Three independent findings: (1) global asymptotic stability — max Re(eigenvalue) = −5×10⁻⁶ /s at every SPL, equal to the protein-decay rate; no positive real parts, no Hopf bifurcation, no limit cycles; (2) no true bistability — the apparent 49% hysteresis gap under 96 h dwell collapses to 0.000% under 1000 h dwell from either initial condition (rest SS vs 110-dB SS), so the ramp-observed path-dependence is purely a slow-protein-decay artifact (t½ ≈ 38 h); (3) the apparent sharp switch at 31 dB is a Ca-forcing artifact: the published ca_forcing has an if SPL ≤ 30: Ca = rest branch that causes a discontinuous jump. A smooth-sigmoid Ca-forcing control reduces the Hill coefficient from n=60 (clipped) to n=4.33 with K=19.4 dB (R² = 0.9998), matching the modest steepness expected from 2-subunit Hill kinetics at each step in the cascade.

Method

  1. Imported the published ODE rhs(t, y, SPL_dB, am_freq_hz) from ca_oscillation_ac1_creb_pivot.py by importlib (keeps a single source of truth for the model).
  2. Fine SPL sweep 20→110 dB, 2 dB step, 72 h integration from rest IC per point.
  3. Jacobian by central finite difference (h = 10⁻⁴ relative); eigenvalues via numpy.linalg.eig.
  4. Hill fit: P(SPL) = P_min + (P_max − P_min)·SPL^n / (K^n + SPL^n), bounds n ∈ [0.1, 60].
  5. Hysteresis: pre-converge rest and 110-dB steady states for 400 h each; forward ramp 20→110 dB (from rest SS) and reverse ramp 110→20 dB (from 110-dB SS), 96 h dwell per step.
  6. Closure test: at 30 dB, integrate from both ICs for 1000 h; if protein converges, there is no bistability.
  7. Smooth-Ca control: monkey-patched ca_forcing to remove the if SPL ≤ 30 step, replacing with a single sigmoid across all SPL; re-swept at 2 dB resolution to measure the true Hill slope of the cascade (independent of the model’s Ca boundary condition).

Deterministic. Runtime ≈ 5 min on one CPU.

Results

Linear stability (eigen analysis)

QuantityValue
SPL range tested20–110 dB, 2 dB resolution (46 points)
Max Re(λ) across sweep−5.00×10⁻⁶ /s
Unstable points0 / 46
Dominant (slowest) eigenvalueprotein-decay K_PROT_DECAY_S = 5×10⁻⁶ /s (t½ ≈ 38 h)
Fastest eigenvalueCa²⁺-relaxation, K_EXTRUSION_S = 30 /s
Hopf detectedno (no complex eigenvalues crossing zero)
Saddle-node detectedno (single fixed point per SPL)
Verdictasymptotically stable across the full physiological SPL range

Hysteresis (ramp + closure)

TestProtein at 30 dBGap
Forward ramp (96 h dwell, from rest SS)134.4
Reverse ramp (96 h dwell, from 110-dB SS)200.649.2%
Closure test (1000 h dwell, from rest SS)134.44
Closure test (1000 h dwell, from 110-dB SS)134.440.000%

The 96 h reverse ramp does not fully settle because protein half-life (38 h) is comparable to the dwell time. 1000 h dwell (≈ 26× t½) converges to single basin → no true bistability.

Sharp-switch artifact (smooth-Ca control)

Model variantHill nK_dBNotes
Published (if SPL≤30: Ca=25 nM)60 (clipped at bound)31.0 dB0.995artefactual step at 31 dB
Smooth-sigmoid Ca-forcing4.3319.4 dB0.9998true cascade steepness

Cascade composition: 4 Hill steps with individual n≈2 would predict an aggregate n of ≈ 4–8 under perfect independence (multiplicative Hills). n=4.3 falls inside that range and is the correct biological estimate.

Interpretation for Misha

  • Safety: the construct is provably non-oscillatory and non-bistable under Phase 1 parameter choices. No Hopf bifurcation means no risk of spontaneous oscillations in CREB-P or STRC transcription. No bistability means the system always relaxes to the same steady state at a given sustained SPL — predictable clinical behavior.
  • The self-dosing window is narrow under current parameters (K=19 dB, n=4.3). Any sustained ambient sound ≥ 25–30 dB saturates the response. This matches the earlier NFAT-model finding (STRC Sonogenetics Bifurcation Analysis) and confirms it is a system-level property, not an artifact of the NFAT model specifically.
  • Published “sharp switch at 30 dB” claim is wrong: the apparent step was a model-code artifact. Once the discontinuity is removed, the cascade behaves like a standard 4-Hill bio-amplifier with a threshold near whisper-level sound. Reporting should be updated.
  • To move the switching threshold up to 45–50 dB (where hearing aids actually modulate), the intervention must be upstream of Ca — a Ca-sensor swap, threshold desensitization of AC1/CaM binding, or adding a kinetic filter (CaMKII autoinhibition). No single-parameter perturbation within the published ranges achieves this.

Limitations

  • 1000 h dwell confirms no bistability at 30 dB. Other SPLs between 25 and 35 dB might still harbor very-slow basins; the Jacobian analysis is stronger (proves no local instability anywhere on the sweep).
  • Jacobian was computed via finite differences; condition-number checks were not logged. Eigenvalues near 10⁻⁶ are likely trustworthy because the ODE rhs is analytical and well-conditioned, but a symbolic Jacobian (e.g., sympy) would close that gap.
  • The Ca-forcing step at 30 dB was monkey-patched for the control; the original parameter choice should be revisited upstream — the smooth-sigmoid variant is biologically more defensible than the if/else form.
  • No noise/stochastic terms. A Langevin variant (Ca²⁺ shot noise, transcriptional bursting) could reveal sub-threshold stochastic switching invisible to the deterministic ODE.

Next steps

  1. Update published model (ca_oscillation_ac1_creb_pivot.py) to replace the if SPL ≤ 30 branch with a smooth sigmoid; re-run Phase 1 sweep + Phase 2 sensitivity with the corrected model and note any downstream changes to robustness index.
  2. Parameter-search for therapeutic window: scan K_CA_AC1_NM × AC1_VMAX × K_CAMP_PKA_NM × K_CREB_CRE_HALF to find a combination with n≈4 AND K≈50 dB, if one exists. Expected: upstream Ca-sensor engineering is required; the downstream cascade cannot shift the threshold by more than ~10 dB.
  3. Stochastic variant: add Ca²⁺-shot-noise SDE (Gillespie or Euler-Maruyama) and re-test whether noise can split the stable FP into a bimodal distribution (transcriptional bursting as an alternative “self-dosing” mechanism).

Replication

cd ~/STRC/models
/opt/miniconda3/bin/python3 ca_oscillation_phase3_bifurcation.py
# outputs: ca_oscillation_phase3_bifurcation.json

Files / Models

  • ~/STRC/models/ca_oscillation_phase3_bifurcation.py — full bifurcation sweep + Jacobian + hysteresis-closure + smooth-Ca control
  • ~/STRC/models/ca_oscillation_phase3_bifurcation.json — eigenvalues per SPL, Hill fit, closure-test results, smooth-Ca sweep

Connections