Topological Mismatch and Bifurcation Loss in Autophagy Model Reductions

how a monotonicity-preserving simplification of the AMPK–mTORC1–ULK1 switch quietly deleted the biology's only mechanism for oscillation — and what one real, primary-source rate constant looks like once it's actually checked against its source

WP-31 established a four-stage calibration pipeline specifically to prevent the failure WP-30 documented: an abstract eigenvalue rescaled into a unit-bearing number by a citation that didn't hold up. This paper reports a different failure mode, caught while attempting to complete Stage 2 honestly rather than skip it.

WP-31C built a three-species reduced ODE model of the AMPK–mTORC1–ULK1 switch, verified it computationally, and used it to argue the pathway cannot sustain oscillation for any generic parameter choice. That conclusion is correct — for the model as built. It is not a fact about the actual biology. The published, peer-reviewed model of this exact network (Szymańska et al., PLOS ONE 2015) contains a feedback edge the reduced model omitted, and that edge is precisely what the real system uses to oscillate. The omission was structural, not numerical: no choice of rate constants could have produced oscillation in the reduced topology, because the topology itself forbids it.

1. What the reduced model actually proved

The three-variable active-fraction system

dA/dt = k_a1*S*(1-A) - k_a2*A dM/dt = k_m1*(1-M) - k_m2*M*A - k_m3*M*U dU/dt = k_u1*(1-U)*A - k_u2*U*M

was analyzed in full, and every claim below was checked computationally against the model as actually specified:

Result	Method
CRNT topology: δ = 2	Explicit graph construction (n=9 complexes, l=4 linkage classes, s=3 stoichiometric rank) — an earlier hand count had mis-stated l=3 by missing a shared complex
Cooperativity	K = diag(+1,−1,+1) applied to the actual computed Jacobian at S ∈ {0.1, 1.0, 2.5, 5.0}; all off-diagonal entries confirmed ≥ 0; loop-gain J₂₃J₃₂ confirmed positive (≈0.78–0.84)
No oscillation	Hirsch (1982, 1985) + Markus (1956)/Thieme (1992); verified directly by simulation from five widely separated initial conditions, all converging to the identical fixed point with zero tail variance
No integer clustering	Numerical continuation, S ∈ [0.1, 5.0], 50-point resolution: λ_max(S) ∈ [−0.65, −2.82] min⁻¹, smooth and monotonic

Every one of these results is correct as a statement about the model that was built. The error is upstream of all of them: the model itself is not a reduction of the published biology, it is a different, simpler topology that happens to share three variable names with it.

2. Topological mechanism of the bifurcation loss

The block structure of the surrogate Jacobian makes the loss mechanical, not just qualitative. With A ordered first, the matrix is lower block triangular: the zero entries sit above the diagonal (A's equation has no dependence on M or U), while the coupling sits below it (A still drives M and U). For any block-triangular matrix, the spectrum decomposes as the union of the diagonal blocks' spectra:

sigma(J_surrogate) = {-1/tau_A} U sigma(J_(M,U))

and J_(M,U), the 2×2 cooperative block, is forced to have real eigenvalues unconditionally — not merely generically. For any 2×2 matrix (a b; c d) with b,c ≥ 0 (the cooperative sign pattern), the eigenvalues are [(a+d) ± √((a−d)²+4bc)] / 2, and the discriminant (a−d)²+4bc is a sum of two nonnegative terms, hence never negative. No complex-conjugate pair can appear — a more elementary and more airtight statement than the general Metzler/Perron–Frobenius result (which only guarantees a real leading eigenvalue for n≥3, not reality of the whole spectrum). Since a Hopf bifurcation requires exactly a complex-conjugate pair crossing the imaginary axis, the decoupled block cannot produce one at any parameter value, by this discriminant alone.

Surrogate Graph (WP-31C): Primary Graph (Szymanska 2015): [ AMPK (A) ] [ AMPK (A) ] | ^ | v u2| | +-----+-----+ | v v v +-+-----+---+ [ mTORC1 ] --> [ ULK1 ] [ mTORC1 ] <--> [ ULK1 ] (feedforward, acyclic at A) (closed loop through A) J_01 = J_02 = 0 -> real spectrum only J_02 != 0 -> complex pairs possible

Restoring the feedback edge (J₀₂ ∝ u₂, the ULK1→AMPK dephosphorylation/phosphorylation pair u₂, p₆) breaks the block-triangular structure entirely. The result is a genuine 3×3 coupled system with no forced block decomposition, and the discriminant argument no longer applies — complex-conjugate eigenvalue pairs are permitted, which is exactly what a Hopf bifurcation requires and exactly what Szymańska et al. report. The mechanism is not subtle: one nonzero entry, J₀₂, is the entire difference between a spectrum that is provably always real and one that is free to cross into the complex plane.

3. What the published model actually contains

Szymańska, Martin, MacKeigan, Hlavacek & Lipniacki (PLOS ONE 10(3): e0116550, 2015) built a rule-based (BioNetGen/BNGL) mechanistic model of the same three-node circuit, fit and analyzed in the original paper. Two structural features are stated directly in the text and are absent from the reduced model above.

A closed feedback loop, not a cascade. The published network includes "negative feedback from ULK1 to AMPK," mediated by ULK1-dependent phosphorylation of sites in the AMPK α subunit, governed by rate constants the authors name u₂ (dephosphorylation of those sites) and p₆ (the phosphorylation itself). With the feedback loop active, the system is a genuine relaxation oscillator, transitioning through a saddle-node-on-invariant-circle (SNIC) bifurcation and both subcritical and supercritical Hopf bifurcations as AMPK* or rapamycin* level varies — real limit cycles, sustained across a broad region of input space (their Fig. 5). Removing the feedback edge, as the reduced model does, removes the only mechanism the real system has for oscillating. It does not reveal that the real system doesn't oscillate; it constructs a different system that structurally can't.

A nonlinearity the reduced model also lacks. The authors state that simple mutual inhibition is not sufficient to produce switching at all: "multi-site phosphorylation of RPTOR by ULK1 provides the nonlinearity needed for mutual inhibition to give rise to hysteretic switching." The reduced model's M–U coupling is plain bilinear mass action with no multisite cooperativity.

Diagnostic verdict Forcing Szymańska's fitted rate constants into the reduced 3-node Jacobian would not have been a real test of anything. The surrogate cannot oscillate regardless of which numbers are supplied, because the missing feedback edge — not the kinetic rates — is what the real oscillation depends on. Any Monte Carlo uncertainty propagation run on that Jacobian, however carefully executed, would have produced a confident, precise, and meaningless answer.

4. The one number that did clear

Buried in the published discussion of parameter sensitivity, stated as plain text rather than locked in a table image, is this:

Szymańska et al. 2015, in text "the value of u₀ sets the system's characteristic response time (1/u₀ = 100 s)"

u₀ is the paper's lumped dephosphorylation rate constant — the single parameter, by the authors' own sensitivity analysis, that the network's overall response timescale depends on most strongly. That gives u₀ ≈ 0.01 s⁻¹, directly stated by the primary source, for the dominant relaxation rate of the actual AMPK–mTORC1–ULK1 switch.

This is logged here as the first rate constant in this corpus's autophagy material that traces cleanly from claim to source with no arithmetic gap and no ambiguity about which paper it came from — the failure mode that sank WP-30's τ_mTOR figures (two files, two numbers, neither traceable) does not reproduce here.

It is not being proposed as a replacement anchor for μ = −2 or −3. 1/u₀ = 100 s is a real system response time on a real switch; it is three to four orders of magnitude away from a single-digit second, consistent with the hours-scale timescales the broader ULK1/mTORC1 literature survey already indicated. It does not connect to the dimensionless −2 or −3 by any relation established in this corpus, and no such relation is claimed. The other 21 parameters in the model's Table 1 remain unrecovered — PLOS renders that table as a raster image, not machine-readable text.

5. Standing status

Item	Status
Reduced 3-ODE model: δ=2, no limit cycles, continuous non-integer spectrum	Verified, for the model as built
Reduced 3-ODE model as a stand-in for the published network	False. Missing the ULK1→AMPK feedback edge (u₂, p₆) that produces the real SNIC/Hopf oscillations, and the multisite RPTOR nonlinearity required for switching at all
u₀ ≈ 0.01 s⁻¹ (τ=100 s), Szymańska et al. 2015	Established as a primary-source, traceable, unit-bearing rate constant. Not connected to μ=−2/−3
Full 22-parameter set (Table 1)	Not yet recovered — locked in a raster table image; S1 File (BNGL source) is the correct target
μ_max physiological calibration	Still not established, and this paper does not attempt to re-establish it

Path forward: pull S1 File (the BNGL model specification) to recover the full parameter table in machine-readable form and confirm u₀, u₂, p₆, and the multisite RPTOR rules exactly as specified; then rebuild the reduced ODE model to actually contain the feedback edge and the nonlinearity before fitting anything to it — a corrected topology first, a parameter fit second, in that order. Fitting real numbers to the wrong graph is not a smaller error than fitting wrong numbers to the right graph; it produces a result that looks quantitative while answering a question nobody asked.

All numerical and structural claims in §1–2 were executed directly (Python 3, NumPy/SciPy/NetworkX) against the model as specified, not derived by hand. All claims in §3–4 were checked against the primary source text of Szymańska et al. (2015), retrieved directly, rather than taken from a secondary summary. Method per WP-31B.

← WP-31C · Executing the Pipeline WP-32 · The Forced Urgency Gap →