Toward a Unified Theory of Aging: A Tri-Domain Dynamical Framework with an Explicit Causal Hierarchy, Operational State Variables, and Comparative and Interventional Tests

Abstract

Aging research has produced a rich but fragmented theoretical landscape. Evolutionary theories explain why aging is permitted by natural selection; mechanistic theories catalog molecular, cellular, and systemic drivers; and interventional frameworks propose how trajectories can be altered. No single tradition specifies, within one mathematical object, the ultimate permissibility of aging, the proximate biological drivers with a causal hierarchy, and the conditions under which those drivers are controllable.

We propose a tri-domain Unified Theory of Aging that occupies the intersection U = D ∩ C ∩ I of Darwinian/ultimate (D), causal/mechanistic (C), and interventional/control (I) explanations. We define aging as the progressive, stochastic loss of organismal homeostatic capacity, driven by molecular damage accumulation and information degradation under declining selection pressure, manifesting as increased vulnerability to perturbation and a rising cost of functional restoration. We commit to an explicit causal hierarchy — molecular damage under imperfect maintenance is the ground zero; information corruption, cellular senescence, and regenerative failure are downstream consequences and feedback amplifiers.

The theory is instantiated as a minimal four-variable ODE system in damage D, senescent burden S, epigenetic information integrity E, and regenerative capacity R, each normalized to [0,1] with operational biomarker definitions (8-OHdG, p16^INK4a, methylation entropy, Ki67 progenitor fractions). We compute the Jacobian, establish that the youthful state X_Y = (0,0,1,1) is a saddle-like boundary state requiring active maintenance (critical repair threshold R_critical = βδ/(γζu_S)), and that the aged state X_A = (1,1,0,0) is a stable boundary attractor under insufficient intervention. Evolution enters through age-dependent repair σ(a) = σ₀ exp(−a/τ_repr); heterogeneity through a distributional McKean–Vlasov extension over clonal cell-state measures.

From this analysis we derive five sharp, falsifiable predictions — senolysis before reprogramming, state-space rate as a mortality predictor, Lyapunov-like V-scores for drug-development progression, a controllability window for rapamycin, and entropy-based stratification of response to the AI-discovered TNIK inhibitor rentosertib in idiopathic pulmonary fibrosis — each with a pre-specified kill criterion. Comparative biology (naked mole-rat, bowhead whale, long-lived bats, humans vs. mice) supplies a stringent cross-species test: species rank order of effective α, γ, δ, κ should predict lifespan rank order. A concrete closed-loop design illustrates adaptive gerotherapy on aged mouse liver.

The contribution is sevenfold. (1) An explicit causal hierarchy that distinguishes primary damage from downstream amplifiers. (2) Operational biomarker definitions for each state variable (8-OHdG, p16, methylation entropy, Ki67, and others). (3) A minimal ODE system with explicit Jacobian, bifurcation threshold, and phase-portrait analysis. (4) An evolutionary parameter σ(a) that turns antagonistic pleiotropy and disposable soma into mathematical rather than verbal claims. (5) A distributional extension (McKean–Vlasov) for clonal heterogeneity. (6) Comparative-biology predictions for cross-species lifespan rank order. (7) A disease-versus-aging operational distinction and a concrete closed-loop design for aged mouse liver. Together with five kill criteria, the framework is offered not as a final answer but as a testable scaffold that invites its own refutation.

1. Introduction

Aging is among the most consequential biological phenomena and one of the most resistant to theoretical closure. It is simultaneously a demographic pattern, an evolutionary outcome, a molecular process, a systems failure, and an emerging therapeutic target. The difficulty of formulating a unified theory arises because aging does not belong to a single explanatory level. It is not merely a consequence of oxidative damage, telomere attrition, epigenetic drift, senescent-cell accumulation, inflammation, stem-cell exhaustion, or mitochondrial dysfunction. Nor is it fully explained by the decline in the force of natural selection with age, by tradeoffs between reproduction and maintenance, or by antagonistic pleiotropy alone. Each captures part of the phenomenon; none alone explains why aging is permitted, what state changes causally drive it, and how those changes can be therapeutically controlled.

The history of aging theory can be understood as the emergence of three partially overlapping explanatory traditions. The descriptive or ultimate tradition (Weismann, Medawar, Williams, Hamilton, Kirkwood) asks why aging exists. The causal and mechanistic tradition (Harman, Hayflick, Orgel, Olovnikov, Harley, Franceschi, López-Otín and colleagues) identifies molecular, cellular, and systemic processes that generate aging phenotypes. The interventional or engineering-oriented tradition (de Grey, Kennedy, Barzilai, Sinclair, Yamanaka-era reprogrammers) asks how aging can be modified. Each tradition is indispensable; none is sufficient on its own.

The mechanistic tradition has been extraordinarily productive. Lifespan can be extended in model organisms by manipulating insulin/IGF signaling, mTOR, AMPK, sirtuins, autophagy, mitochondrial function, senescent-cell burden, and inflammatory signaling (Kenyon et al., 1993; Johnson, 1990; Blagosklonny, 2006; Harrison et al., 2009; Baker et al., 2011). Epigenetic clocks can estimate biological age with remarkable accuracy across tissues (Horvath, 2013). Partial reprogramming can restore youthful features in cells and tissues under controlled conditions (Takahashi and Yamanaka, 2006; Ocampo et al., 2016; Lu et al., 2020; Yang et al., 2023). Senolytics can reduce senescent-cell burden and improve function in preclinical models (Baker et al., 2011; Xu et al., 2018). The interventional tradition represents a decisive shift from asking whether aging can be understood to asking whether aging can be therapeutically controlled.

Yet each tradition is incomplete alone. Evolutionary theories explain permissibility but do not specify actionable causal nodes. Mechanistic theories identify components without a unifying account of causal hierarchy or why repair fails. Interventional frameworks prescribe repair strategies but lack formal models of causal order, timing, and state-dependent efficacy. Recent information-centered theories, particularly the Information Theory of Aging (Sinclair and LaPlante, 2019; Yang et al., 2023), represent an important bridge: they propose that aging involves loss of regulatory information that compromises cellular identity, and that aspects of this information can be restored. If aging involves corrupted information rather than only molecular wear, rejuvenation becomes biologically plausible. Nevertheless, information loss interacts with DNA damage, inflammation, metabolism, mitochondrial function, extracellular matrix state, immune surveillance, and stem-cell niches. Reprogramming in a damaged environment may not restore function without concurrent damage repair and immune-surveillance control. Information theory is a central causal layer but not the entire architecture.

We therefore propose a tri-domain Unified Theory. A complete theory must simultaneously answer three questions: Why is aging permitted? (D), What biological state changes drive it? (C), and Which variables are controllable? (I). In set-theoretic terms, prior theories often occupy one or two regions of this conceptual Venn diagram. Evolutionary theories occupy primarily D. Damage and hallmark theories occupy primarily C. SENS and geroscience occupy C ∩ I. Information theory occupies a particularly important C ∩ I region and partially bridges to D. A unified theory must occupy U = D ∩ C ∩ I.

This paper's contributions are: (1) an explicit causal hierarchy (§4) placing molecular damage as ground zero with information corruption, senescence, and regenerative failure as downstream amplifiers; (2) a minimal ODE system (§5) in (D, S, E, R) with operational biomarker definitions, parameter units, Jacobian, bifurcation analysis, age-dependent evolutionary parameters, and a distributional extension for clonal heterogeneity; (3) a comparative biology test (§6) — the same ODE parameters should rank species by lifespan; (4) a disease-versus-aging operational distinction and a concrete design loop (§7) together with the rentosertib / TNIK case study; (5) five falsifiable predictions with pre-specified kill criteria (§8); (6) explicit situation against quantitative predecessors (§9); (7) ICMJE-compliant AI disclosure (§11).

Why now. Three developments make this synthesis timely. First, AI-driven drug discovery has produced the first aging-relevant compounds reaching clinical trials through computational target identification and generative chemistry (Zhavoronkov et al., 2019; Aladinskiy et al., 2024; Xu et al., 2025), creating demand for theoretical frameworks that can guide computational prioritization. Second, longitudinal multi-omics cohorts (UK Biobank, CALERIE, Tabula Muris Senis) now provide the data density needed to calibrate dynamical models. Third, single-cell and spatial transcriptomics have revealed the heterogeneity of aging at cellular resolution, requiring distributional rather than mean-field theories. The convergence of computational discovery, longitudinal data, and single-cell resolution creates both the opportunity and the necessity for a formal dynamical theory.

Relationship to Insilico Medicine’s research program. This framework emerged from practical experience in AI-driven aging research. Over the past decade, our group has developed AI platforms for target identification (PandaOmics), molecule generation (Chemistry42), and clinical trial prediction. The discovery of rentosertib — a TNIK inhibitor for IPF identified through dual-purpose aging/disease targeting — demonstrated that aging biology can guide drug discovery when formalized computationally. The present paper provides the theoretical architecture that such computational platforms implicitly require: a formal state-space model with causal hierarchy, parameter sensitivity, and intervention-ordering predictions.

A note on scope. We deliberately restrict the model to four aggregate state variables because identifiability, interpretability, and empirical testability are primary virtues at this stage. Extending the model to include mitochondrial function, proteostasis, stem-cell reserve, microbiome state, immune repertoire, extracellular-matrix stiffness, and systemic signaling factors is straightforward; we anticipate such extensions in companion work. What matters here is that the four-variable system is already rich enough to exhibit a youthful saddle, an aged attractor, a bifurcation threshold, and non-trivial intervention-ordering effects — the minimum set of behaviors required for a dynamical theory of aging to have teeth.

2. Operational Definition of Aging

A theory is only as strong as the definition of its object. We adopt a deliberately tight operational definition:

Four elements deserve emphasis. First, aging is progressive and stochastic: directional on average but not deterministic at the level of individual trajectories. Second, the core change is loss of homeostatic capacity, not mere accumulation of biomarkers. Two individuals with identical chronological age may occupy different positions in a state space of homeostatic reserves. Third, the proximate drivers are damage and information degradation, corresponding to the D and (1 − E) variables in our dynamical model. Fourth, aging manifests as increased vulnerability (reduced resilience) and as a rising cost of restoration (larger interventions are required to return the organism to a given functional state).

This definition explicitly excludes three things. It excludes mere chronological time: aging is not the passage of time but the dynamical consequence of finite maintenance under ongoing stress. It excludes any single mechanism: aging is not reducible to oxidative stress, senescence, or epigenetic drift alone. And it excludes a purely programmed aging mechanism: the framework does not assume aging is adaptively selected. Rather, aging is the dynamical consequence of incomplete maintenance under declining selection pressure, as first formalized by Hamilton (1966).

2.1 The tri-domain requirement

A complete theory of aging must satisfy three explanatory requirements simultaneously. First, it must explain why aging exists under natural selection — the descriptive/ultimate domain (D). Second, it must explain what biological processes generate aging phenotypes — the causal/mechanistic domain (C). Third, it must explain how those processes can be modified, repaired, reversed, or controlled — the interventional/therapeutic domain (I).

The descriptive domain asks: why did evolution not produce indefinite somatic maintenance? It includes declining force of selection, mutation accumulation, antagonistic pleiotropy, disposable soma, life-history tradeoffs, ecological mortality, germline–soma asymmetry, and reproduction–maintenance allocation. Its central output is not a molecular mechanism but a boundary condition: organisms are optimized for reproductive fitness, not indefinite survival.

The causal domain asks: what biological variables change with age and causally produce pathology, frailty, mortality, or loss of resilience? It includes genome instability, telomere attrition, epigenetic drift, mitochondrial dysfunction, proteostasis failure, nutrient-sensing dysregulation, senescence, stem-cell exhaustion, inflammaging, immune remodeling, extracellular matrix stiffening, microbiome dysbiosis, somatic mosaicism, and loss of tissue architecture.

The interventional domain asks: which aging-relevant state variables are measurable, causally important, and controllable? It includes senolytics, rapalogs, metformin-like metabolic modulation, autophagy enhancement, mitochondrial restoration, stem-cell replacement, immune rejuvenation, telomere restoration with cancer safeguards, extracellular matrix repair, proteome modulation, partial reprogramming, and adaptive biomarker-guided combinations.

2.2 Why previous theories were domain-limited

Previous theories were incomplete not because they were wrong, but because they were domain-limited. Evolutionary theories explain permissibility but not molecular causality or therapeutic control. Damage theories identify lesions but do not explain why repair systems are incomplete or which damage classes should be prioritized clinically. Programmatic theories explain regulatory dysregulation but can overstate the existence of an intentional aging program. Engineering frameworks identify intervention classes but may lack a formal model of causal hierarchy and systemic coupling.

This fragmentation has practical consequences. If aging is viewed only through evolution, intervention appears secondary. If viewed only through molecular damage, regulatory reversibility may be underestimated. If viewed only through reprogramming, irreversible or slowly reversible structural lesions may be neglected. If viewed only through hallmarks, the field risks replacing theory with taxonomy. The tri-domain framework prevents these errors by requiring that every proposed mechanism be evaluated for evolutionary origin, causal role, and therapeutic controllability.

2.3 The power of multi-domain theories

The most influential aging theories have been those spanning more than one domain. Disposable soma links evolutionary reasoning to maintenance biology. Antagonistic pleiotropy connects selection with causal tradeoffs. SENS links causal damage classes to repair strategies. The Hallmarks framework links mechanisms to drug discovery. Information theory links epigenetic causality to rejuvenation. Sinclair's information-centered framework deserves particular recognition because it made reversibility central to causal theory — by arguing that aging involves loss of epigenetic information rather than only accumulated irreversible damage, it provided a conceptual bridge between mechanism and intervention.

We define a Unified Theory of Aging as a theory occupying the intersection of all three domains: U = D ∩ C ∩ I. Such a theory must explain aging as an evolved consequence of incomplete maintenance, specify the interacting state variables whose drift generates decline, and identify which variables are measurable and therapeutically controllable.

2.4 The four causal layers

The Unified Theory organizes aging into four interacting causal layers: evolutionary, information, damage, and programmatic. These layers are not independent — they form a coupled dynamical system in which changes at one layer alter the rates, thresholds, and consequences of change in the others.

The evolutionary layer defines the selection-shaped parameters of the organism: age-specific selection, life-history strategy, maintenance allocation, reproductive timing, cancer suppression, and ecological mortality. It does not directly “cause” molecular aging; rather, it determines why repair and resilience systems are finite and imperfect.

The information layer includes genomic sequence integrity, epigenetic patterning, chromatin organization, transcriptional fidelity, cellular identity, spatial tissue information, and intercellular signaling codes. Aging involves loss, corruption, misreading, or maladaptive redeployment of this information.

The damage layer includes molecular and structural lesions: DNA damage, mitochondrial mutations, oxidized proteins and lipids, aggregates, extracellular crosslinks, telomere attrition, membrane damage, lysosomal waste, and mechanical tissue deterioration.

The programmatic layer includes adaptive biological responses that become maladaptive with age: senescence, inflammation, fibrosis, hyperactive nutrient signaling, compensatory proliferation, stress-response remodeling, immune exhaustion, and altered wound-healing programs. These are programs for development, defense, repair, growth, and survival that become mistimed, chronic, excessive, or misregulated with age.

2.5 The cascade model

The causal cascade can be summarized as:

This sequence is not strictly linear — it describes dominant causal tendencies while allowing feedback loops. Evolutionary permission means that selection does not enforce indefinite maintenance. Information loss and damage accumulate bidirectionally — DNA breaks recruit chromatin modifiers; chromatin disorganization impairs repair targeting; mitochondrial dysfunction alters metabolites required for epigenetic enzymes; inflammation remodels transcriptional programs. Programmatic dysregulation arises when adaptive pathways become chronic or misapplied. Resilience failure is the late-stage systems phenotype: the capacity to return to functional equilibrium after stress is progressively eroded.

2.6 Implications for intervention strategy

The causal architecture predicts that no single intervention class is likely to fully control aging. Damage repair without information restoration may leave cells misregulated. Reprogramming without damage repair may restore transcriptional youth while leaving structural lesions intact. Anti-inflammatory therapy without senescent-cell control may suppress symptoms but not remove sources. Senolytics without regenerative support may remove harmful cells but fail to rebuild tissue. Therefore, effective gerotherapy requires combinations targeting coupled causal nodes in the correct sequence, determined by the current state of the system.

3. Condensed Historical Review

The literature on aging is vast. We present only the minimum necessary to locate the framework within its intellectual context. The goal is not encyclopedic coverage but mapping each tradition onto the tri-domain structure {D, C, I} and identifying the residual problem each tradition leaves open.

3.1 Evolutionary theories: why aging is permitted

Weismann (1882) first framed aging as a biological phenomenon requiring evolutionary explanation and proposed an adaptive programmed-death concept. The group-selection logic was later rejected as a general explanation, but the framing itself proved durable: aging required evolutionary explanation rather than treatment as mere wear.

Medawar (1952) provided the first mature evolutionary theory of aging with mutation accumulation. The force of natural selection declines with age because fewer individuals survive to later ages in natural environments and because reproductive value typically decreases after maturity. Deleterious mutations whose effects are confined to late life are weakly selected against and accumulate in populations. Williams (1957) proposed antagonistic pleiotropy: alleles conferring early-life benefits can be favored by selection even if they produce late-life harm. Many canonical aging pathways — insulin/IGF, mTOR, senescence itself — support early fitness while contributing to late-life pathology when chronically active. Kirkwood (1977) added disposable soma theory, arguing that organisms allocate limited resources between reproduction and somatic maintenance, and that aging results from evolutionarily optimized but incomplete investment in repair.

Hamilton (1966), in "The moulding of senescence by natural selection," mathematically formalized age-specific selection. He showed that the force of selection against mortality and fertility effects declines predictably with age after reproductive maturity. Hamilton's formalism quantifies the age-specific reproductive value and the sensitivity of fitness to age-specific survival; it demonstrates that even small early-life fertility benefits can outweigh large late-life mortality costs in expected fitness terms. This supplied the rigorous demographic foundation underlying mutation accumulation and antagonistic pleiotropy.

Blagosklonny (2006) proposed hyperfunction theory, arguing that aging is not primarily failure of function but continued quasi-programmed growth pathway activation past developmental utility. On this view, mTOR and related growth pathways — selected for developmental and reproductive fitness — continue signaling into late life and cause pathology through chronic anabolic pressure on differentiated tissues. The appeal of hyperfunction theory is that it connects antagonistic pleiotropy (evolutionary permissibility), mTOR dysregulation (mechanism), and rapamycin (intervention) into a single causal chain — making it one of the rare theories to span all three tri-domain regions. Its limitation is that hyperfunction alone cannot explain damage-centric aging phenomena such as extracellular crosslinking, amyloid deposition, or somatic mutation burden; these require integration with classical damage accumulation and with the Information Theory of Aging.

In the tri-domain Venn, the evolutionary tradition occupies primarily D. It explains why maintenance is imperfect but does not by itself identify mechanistic lesions or prescribe interventions, with hyperfunction theory as a partial exception that extends into C and I.

3.2 Damage and mitochondrial theories

Harman (1956) proposed the free-radical theory of aging: cumulative damage from reactive oxygen species generated during metabolism progressively impairs cellular and organismal function. The theory later extended to mitochondria as both sources and targets of oxidative damage (Harman, 1972). Reactive oxygen species are real damage agents, but translation into antioxidant therapy has been disappointing — antioxidant supplementation has generally failed to robustly extend lifespan in humans and has produced mixed results in animals. Genetic studies also showed that reactive oxygen species are not merely damaging byproducts but signaling molecules involved in hormesis, adaptation, immunity, and stress response. Some long-lived models exhibit increased oxidative stress resistance, but not all lifespan extension correlates with lower oxidative damage.

Orgel (1963) proposed the error catastrophe hypothesis — that errors in protein synthesis could propagate into a self-amplifying failure of the synthetic machinery. The specific prediction was not borne out, but the underlying insight — that feedback can amplify molecular errors — remains central to modern damage-and-repair models.

In the tri-domain Venn, damage theories occupy primarily C, extending into I through repair strategies. Their limitation is reductionism: damage is necessary but not sufficient as a master cause. In the present framework, oxidative stress is treated as one damage-generating and signaling-modulating component within a coupled network.

3.3 Cellular senescence and telomere biology

Hayflick and Moorhead (1961) demonstrated that normal human fibroblasts have a finite replicative capacity, overturning the assumption of somatic immortality. The "Hayflick limit" was later linked mechanistically to telomere attrition. Olovnikov (1973) predicted that linear chromosomes would shorten with replication. Blackburn and Gall (1978) characterized telomeric repeats, and Greider and Blackburn (1985) identified telomerase. Harley, Futcher, and Greider (1990) showed that telomeres shorten during aging of human fibroblasts, establishing the molecular basis of replicative senescence. Bodnar et al. (1998) demonstrated that telomerase could extend replicative lifespan of human cells in culture.

Telomere dysfunction contributes to degenerative syndromes, bone-marrow failure, pulmonary fibrosis, immune dysfunction, and impaired tissue renewal. Critically, telomere biology is not a universal clock: some long-lived species have different telomerase dynamics, and telomerase activation carries cancer risk. Cellular senescence was recognized as broader than telomere exhaustion: it can be triggered by DNA damage, oncogene activation, mitochondrial dysfunction, and other stresses (Campisi, 2005, 2013; Muñoz-Espín and Serrano, 2014). Senescent cells exhibit altered gene expression and a senescence-associated secretory phenotype (SASP) that propagates inflammation, matrix remodeling, and paracrine senescence.

The Hayflick limit had two downstream effects. First, it established that somatic cell aging is a measurable, reproducible property rather than a consequence of culture conditions. Second, it motivated the search for the molecular clock underlying replicative exhaustion, which was found in telomere attrition. In the ODE model, senescence is downstream of damage and becomes a feedback amplifier through the βS and δD(1−S) terms.

3.4 Nutrient sensing and the mTOR–rapamycin axis

The discovery that reduced insulin/IGF-1 signaling extends lifespan in C. elegans (Kenyon et al., 1993), Drosophila, and mice established nutrient sensing as a central aging-regulatory axis. The mechanistic target of rapamycin (mTOR) complex integrates growth-factor, nutrient, energy, and stress signals to regulate protein synthesis, autophagy, lipid metabolism, and cell growth. Chronic mTOR activation in post-developmental life is proposed to drive aging through continued anabolic pressure on differentiated tissues — the molecular substrate of Blagosklonny's hyperfunction theory.

Harrison et al. (2009) demonstrated that rapamycin, administered late in life to genetically heterogeneous mice, extends median and maximum lifespan. This was confirmed across multiple sites in the NIA Interventions Testing Program (ITP). Rapamycin remains the most robustly replicated pharmacological lifespan extension in mammals. Its mechanism involves mTORC1 inhibition, enhanced autophagy, improved proteostasis, reduced senescent-cell secretory phenotype, and altered immune function.

Caloric restriction (CR) extends lifespan across species from yeast to primates (though primate studies show variable results depending on diet composition and timing). CR acts partly through reduced mTOR, increased AMPK, enhanced autophagy, improved mitochondrial function, and altered NAD+ metabolism. Metformin, an AMPK activator and mild mitochondrial complex I inhibitor, extends lifespan in some mouse models and is being tested in the TAME trial for human aging.

In the tri-domain Venn, nutrient-sensing interventions occupy C ∩ I strongly. The evolutionary connection (→ D) is through antagonistic pleiotropy: growth and nutrient-sensing pathways are selected for early-life fitness and development but produce late-life pathology when chronically active. In the ODE model, rapamycin and CR act primarily through reducing α (basal damage generation via reduced metabolic stress) and enhancing the effective γ (improved autophagy and proteostasis enhance damage clearance).

3.5 Inflammaging

Franceschi et al. (2000) formalized the concept of inflammaging: chronic, low-grade, sterile inflammation accompanies aging and contributes to frailty, cardiovascular disease, neurodegeneration, sarcopenia, diabetes, and mortality. Older individuals often exhibit increased circulating inflammatory mediators such as IL-6, TNF-α, and CRP. Senescent cells secrete pro-inflammatory SASP factors; clonal hematopoiesis can amplify inflammatory signaling and cardiovascular risk. Inflammaging is therefore not merely a downstream correlate but a feedback amplifier: local cellular damage is converted into organism-wide dysfunction through inflammatory signaling. In the ODE model, inflammaging appears implicitly in the βS term of the damage equation and the δD(1−S) term of the senescence equation, which together generate the positive feedback responsible for the unstable eigenvalue at the youthful state.

3.6 SENS and damage-repair engineering

de Grey and colleagues (2002, 2005) reframed aging as an engineering problem of accumulated damage repair. Canonical SENS categories: cell loss and tissue atrophy, death-resistant cells, nuclear mutations and epimutations, mitochondrial mutations, intracellular aggregates, extracellular aggregates, and extracellular matrix crosslinks. Senolytics, aggregate immunotherapies, stem-cell replacement, and crosslink breakers validate the engineering frame. SENS shifted the field from slowing damage formation to repairing damage after it occurs.

SENS occupies C ∩ I strongly but weakly addresses D. It prescribes what to repair without fully explaining why evolution left those repair problems unsolved or how damage classes interact dynamically. In the present framework, SENS damage classes correspond to components of the composite damage index D, and SENS repair strategies correspond to u_D(t) interventions.

3.7 Reprogramming and the Yamanaka era

Takahashi and Yamanaka (2006) demonstrated that differentiated somatic cells can be converted to iPSCs by defined factors (Oct4, Sox2, Klf4, c-Myc). Reprogramming resets many epigenetic marks, restores telomere length in pluripotent cells, and erases certain aging signatures. Ocampo et al. (2016) showed that partial (cyclic) expression of reprogramming factors ameliorates aging phenotypes in vivo without full dedifferentiation. Lu et al. (2020) demonstrated reprogramming-based restoration of vision in aged mice. Yang et al. (2023), using an inducible changes to the epigenome (ICE) model, showed that epigenetic disruption can drive aging-like phenotypes that reprogramming can partially reverse. Risks include cancer, dedifferentiation, incomplete rejuvenation, and tissue-specific variability.

In the tri-domain Venn, reprogramming is strongly C ∩ I and increasingly contributes to D insofar as it implies that aging involves loss or misdeployment of regulatory information. In the ODE model, reprogramming corresponds to u_E(t), with efficacy scaling as θu_ER(1−S) — it is suppressed by high senescent burden and enhanced by regenerative reserve.

3.8 Hallmarks of Aging and the Information Theory of Aging

The Hallmarks of Aging (López-Otín et al., 2013; updated 2023) organized aging biology into nine and then twelve categories: genomic instability, telomere attrition, epigenetic alterations, loss of proteostasis, deregulated nutrient sensing, mitochondrial dysfunction, cellular senescence, stem-cell exhaustion, altered intercellular communication, disabled macroautophagy, chronic inflammation, and dysbiosis. The framework synthesizes decades of molecular gerontology and became dominant because it provided a common language. Its limitation is that it remains partly taxonomic — it identifies categories but does not fully specify causal order, mathematical dynamics, evolutionary origin, or therapeutic prioritization.

The Information Theory of Aging (Sinclair and LaPlante, 2019; Yang et al., 2023) proposes that aging results in substantial part from loss of epigenetic and regulatory information. Cells retain youthful genetic information but progressively lose the ability to read, organize, and deploy it properly. The ICE mouse model introduces controlled DNA breaks without producing a high mutation burden and shows that repeated damage–repair cycles can induce features resembling aging, including epigenetic alteration, transcriptional dysregulation, physiological decline, and acceleration of epigenetic age. Subsequent expression of reprogramming factors can restore some youthful features. This supports the idea that epigenetic information loss is not merely a passive correlate of aging but can be a driver.

The interpretation should be cautious. DNA breaks are not purely informational perturbations; they also activate stress responses, inflammation, senescence, and repair pathways. Thus the ICE model supports an information component of aging but does not prove that information loss is the only or dominant cause of organismal aging — consistent with the tri-domain view that information, damage, and programmatic responses are coupled.

3.9 Quantitative dynamical predecessors

We note the specific dynamical lineage the present framework extends (detailed in §9): Hamilton's demographic formalism, Gavrilov–Gavrilova reliability theory (redundancy loss), Mitnitski–Rockwood deficit accumulation (frailty index), Pyrkov–Fedichev longitudinal resilience analysis (DOSI), and Gao–Barzel–Barabási network resilience theory. Our contribution is to combine these lines of work into an explicit control-oriented state-space model with bifurcation analysis and intervention sequencing.

3.10 Integration attempts and the geroscience era (2014–2026)

The geroscience hypothesis (Kennedy et al., 2014) proposed that targeting fundamental aging mechanisms could delay or prevent multiple chronic diseases simultaneously, rather than treating each disease independently. This reframed aging as the shared upstream risk factor for cancer, cardiovascular disease, neurodegeneration, diabetes, and frailty. Clinical translation includes the TAME (Targeting Aging with Metformin) trial and rapamycin studies in companion animals and humans.

Recent aging research has increasingly moved toward integration. Biological clocks attempt to quantify organismal state across methylation, proteomics, and metabolomics. Systems biology models describe network interactions among metabolism, inflammation, immunity, and tissue function. Reprogramming studies challenge the assumption that aging is unidirectional. AI-enabled target discovery (Zhavoronkov et al., 2019) and multimodal biomarker platforms are accelerating the search for intervention combinations. PandaOmics-based dual-purpose target identification (Pun et al., 2022, 2023) connects aging biology to specific drug-discovery programs. Chemistry42 (Ivanenkov et al., 2023) enables generative molecular design for aging-relevant targets.

Yet integration remains incomplete. Existing frameworks often combine mechanisms without incorporating evolutionary theory; propose interventions without formal causal prioritization; or identify biomarkers without specifying how they relate to controllable state variables. The Hallmarks provide taxonomy but not dynamics. Geroscience provides clinical strategy but not mathematical formalism. Biological clocks provide measurement but not causal ordering. SENS provides engineering targets but not evolutionary context. The Information Theory provides a causal layer but not the full multi-layer architecture.

The field requires a theory that can simultaneously explain permissibility, mechanism, and control — within a single mathematical object that generates falsifiable predictions and supports intervention optimization.

3.11 The residual problem

Every prior framework leaves unaddressed the requirement that a theory simultaneously explain permissibility, specify proximate causation with a causal hierarchy, and supply actionable intervention logic — within one testable mathematical object. That is the target of the present paper.

The residual problem can be stated precisely. What is needed is a mathematical object that: (1) specifies state variables with operational biomarker definitions; (2) encodes causal relationships as differential equations with measurable parameters; (3) admits stability and bifurcation analysis that identifies controllability thresholds; (4) incorporates evolutionary constraints as time-varying parameters; (5) generates quantitative predictions about intervention timing, sequencing, and state-dependent efficacy; (6) makes cross-species predictions testable against longitudinal biomarker data; and (7) provides explicit falsification criteria. No existing framework satisfies all seven requirements simultaneously. The sections that follow attempt to construct such an object.

4. The Causal Hierarchy: From Ground Zero to Systemic Decline

A central concern in any dynamical theory of aging is causal ordering. Many age-associated variables — DNA damage, mitochondrial dysfunction, epigenetic drift, inflammation, senescence, stem-cell exhaustion, fibrosis, and organ decline — covary strongly with chronological age. However, covariation is not causation. A useful model must specify which variables are primary drivers, which are intermediate responses, and which are late emergent consequences.

This hierarchy is consistent with the "ground zero" concept articulated by Gladyshev (2021), who argued that aging originates from unavoidable molecular damage arising from the imperfect maintenance of biological structures in a thermodynamically open but chemically imperfect system. The first event is not a programmed execution of aging, nor is it primarily an epigenetic clock, senescence program, or stem-cell exhaustion. Rather, the initial cause is accumulation of molecular lesions produced by metabolism, environmental exposure, replication, stochastic chemical instability, and incomplete repair. Biological systems oppose this process through maintenance mechanisms — DNA repair, protein quality control, autophagy, mitophagy, antioxidant systems, immune clearance, and tissue renewal — but these mechanisms are finite and evolutionarily optimized for reproductive success rather than indefinite survival (Medawar, 1952; Williams, 1957; Kirkwood, 1977).

4.5 The hierarchy stated counterfactually

The strongest statement of the model is counterfactual: without molecular damage accumulation, senescence would not progressively increase, epigenetic information would not progressively degrade, and regenerative capacity would not progressively decline. Conversely, if damage accumulates beyond the capacity of repair, then information loss, senescence, inflammation, and regenerative failure follow.

This does not imply that downstream variables are unimportant. On the contrary, once activated, S and (1−E) create feedback loops that accelerate decline. Senescent cells increase inflammatory damage; epigenetic disorganization impairs repair gene expression; regenerative failure permits damaged cells and extracellular matrix to persist. Aging is therefore hierarchical but not strictly linear. It begins with damage but evolves into a coupled nonlinear system. The distinction is between causal origin and dynamic amplification.

The model is also compatible with the Hallmarks of Aging. The present contribution is to impose an explicit causal ordering among hallmark features: genomic, proteomic, mitochondrial, and extracellular molecular damage constitute the initiating layer; epigenetic alterations and senescence are damage responses and amplifiers; stem-cell exhaustion and tissue decline are emergent tissue-level consequences.

4.6 Implications for intervention hierarchy

A key implication is that interventions should be evaluated by their position in the hierarchy. Senolytics may reduce S and thereby reduce secondary damage amplification, but they do not necessarily eliminate the primary sources of D. Epigenetic reprogramming may restore E, but if underlying damage persists, restored information may again degrade. Regenerative therapies may increase R, but if delivered into a high-D, high-S environment, newly generated cells may rapidly acquire aged phenotypes. Therefore, rational gerotherapy should prioritize combinations that reduce D, suppress maladaptive S, restore E, and rebuild R in a sequence dictated by the state of the system.

5. The Dynamical Model

The purpose of this section is not to claim that aging has been reduced to four variables. Rather, it is to make explicit the dynamical commitments of the framework in a form that can be analyzed, simulated, and falsified. The system below is a phenomenological first-order model: a deliberately simplified state-space representation intended to expose assumptions, generate predictions, and define measurable failure modes. It should not be read as a final molecular theory.

5.1 State variables with operational definitions

The state vector is

where D = 0 denotes negligible measurable damage and D = 1 denotes maximal damage burden relative to a defined biological scale; S = 0 no detectable senescent burden and S = 1 maximal senescent fraction; E = 1 youthful epigenetic and transcriptional order, E = 0 maximal information loss; R = 1 youthful regenerative capacity, R = 0 regenerative collapse. Raw biomarkers are standardized relative to young-adult and late-life reference distributions, then combined with prespecified or empirically learned weights.

Damage index D — composite of at least six components

1. 8-hydroxy-2′-deoxyguanosine (8-OHdG) in plasma, urine, or tissue (ng/mL or /creatinine) — oxidative DNA damage (Ames et al., 1993).
2. Protein carbonyls (nmol/mg protein) — irreversible oxidative protein damage.
3. Lipofuscin area fraction — histology or autofluorescence microscopy.
4. Advanced glycation end-products (especially carboxymethyllysine, CML, nmol/mg protein) — long-lived glycoxidative modification.
5. Mitochondrial DNA deletions — digital PCR or long-read sequencing (% heteroplasmy).
6. Somatic mutation burden — whole-genome sequencing (mutations per genome or per Mb).

Optional components: 4-hydroxynonenal adducts, γH2AX foci, insoluble protein aggregates. A typical normalization is D_i = (z_i − z_i^young) / (z_i^old − z_i^young), truncated to [0,1]. The composite is D = Σ w_i D_i with Σ w_i = 1. Equal weights are used initially and then refined using longitudinal prediction of frailty, mortality, or organ failure.

Senescent burden S — multi-marker composite

1. p16^INK4a-positive cell fraction by qPCR, immunohistochemistry, single-cell RNA-seq, or flow.
2. SA-β-gal-positive fraction histochemically.
3. SASP composite (IL-6, MCP-1/CCL2, MMP3) in plasma or tissue-conditioned medium, pg/mL, standardized.

A minimal senescence index is S = w₁S_p16 + w₂S_SA-β-gal + w₃S_SASP. The requirement for multiple markers follows from the heterogeneity of senescence programs: p16, p21, SA-β-gal, and SASP expression are overlapping but not identical biological states (Campisi, 2013; Muñoz-Espín and Serrano, 2014).

Epigenetic information E — favorable when high

1. Shannon entropy of CpG methylation: H = −p log p − (1−p) log(1−p) per site. Low entropy = ordered youthful methylation; high entropy = drift.
2. Transcriptomic identity score: Pearson correlation to a young-tissue reference atlas.
3. Chromatin accessibility variance: coefficient of variation of ATAC-seq peak accessibility at age-sensitive regulatory regions.

E = w_H(1−H_norm) + w_T T_identity + w_A(1−A_var,norm). DNA methylation clocks (Horvath, 2013; Hannum et al., 2013; Levine et al., 2018) provide a calibration layer.

Regenerative capacity R — tissue-specific

1. Ki67-positive stem/progenitor cells per niche by IHC or flow cytometry.
2. Colony-forming units per 10⁵ bone marrow cells.
3. Wound closure rate (mm/day) after standardized punch biopsy.
4. Satellite cells per myofiber for skeletal muscle.

R = Σ w_i R_i, each scaled to young-adult performance. Consistent with Goodell and Rando (2015), R captures both intrinsic stem-cell damage and extrinsic niche dysfunction.

5.2 The coupled ODE system

The causal hierarchy is encoded by:

All variables remain on [0,1]⁴ via the projected dynamics interpretation. A mathematically equivalent reduced form convenient for stability analysis, with damage amplification by senescence, senolytic term, epigenetic loss, reprogramming, and repair-support is:

Damage equation. α is unavoidable basal damage generation. βS is SASP-mediated secondary damage amplification. γ_eff(a)RD is repair/turnover, declining with age as selection weakens. u_D(t) is direct damage-clearance intervention.

Senescence equation. δD(1−S) converts nonsenescent cells to senescent states (bounded by (1−S)). ρ_S(t)S or ζu_S(t)S is senolytic/senomorphic intervention. Crucially, if D ≈ 0 then dS/dt ≈ 0 — senescence is downstream of damage in this model.

Epigenetic equation. −κ_DDE and −κ_SSE are damage-driven and SASP-driven information loss. μ_E(t)(1−E) or θu_ER(1−S) is reprogramming, consistent with Ocampo et al. (2016) — suppressed by high S, enhanced by R.

Regeneration equation. −λ_DDR and −λ_SSR represent direct damage to progenitors and niche degradation. ν_EE(1−R) represents regeneration supported by preserved epigenetic order. u_R(t)(1−R) is therapeutic stimulation of regeneration.

Parameter units and illustrative values

Values are illustrative first estimates (Ames et al., 1993; Baker et al., 2016; Horvath, 2013; Ocampo et al., 2016) to be refined from longitudinal data.

5.3 Jacobian, stability, and bifurcation

Using the reduced form, the Jacobian ∂F/∂X is

Aging is not modeled as a scalar clock but as a coupled feedback system.

Youth X_Y = (0,0,1,1) is a saddle without maintenance

At X_Y with no interventions: dD/dt = 0, dS/dt = 0, dE/dt = 0, dR/dt = λ (projected to zero at the boundary). So X_Y is a boundary equilibrium.

The damage–senescence subsystem near youth (E ≈ 1, R ≈ 1) with c = ζu_S has

When c = 0 (no senolytic clearance), μ₊ > 0 whenever βδ > 0. Youth has an unstable direction in the coupled damage–senescence plane. Even infinitesimal damage increases senescence, which increases damage. Youth is not self-maintaining; it is metastable, sustained by active surveillance, repair, immune clearance, proteostasis, and tissue renewal.

With clearance, tr(A_DS) = −(γR + c) < 0 and det(A_DS) = γRc − βδ. Local stability requires det > 0, i.e., γRc > βδ, giving the critical repair threshold:

For a given senolytic intensity u_S, youth is locally stable only if R > R_critical. If u_S = 0, then R_critical → ∞ — no biologically attainable R ≤ 1 can stabilize youth against damage–senescence feedback. This formalizes the central claim: youth requires active maintenance, not passive absence of pathology.

At R = R_critical, μ₊ = 0 — the local bifurcation condition for the onset of runaway damage–senescence amplification. This is the controllability threshold beyond which ordinary maintenance no longer restores the organism to the youthful basin.

Aged X_A = (1,1,0,0) is a stable boundary attractor

At X_A with no interventions: dD/dt = α + β, dS/dt = ε, dE/dt = −2η, dR/dt = −3κ. All point outside [0,1]⁴, so X_A is a boundary equilibrium under projection.

Using distance variables Y = (1−D, 1−S, E, R), the aged state corresponds to Y = 0. Near X_A, dY_i/dt ≤ −m < 0 for each active coordinate. With V(Y) = y_D + y_S + y_E + y_R, we have dV/dt < 0 near Y = 0. Under smooth boundary-layer regularization of width ℓ, the active-set eigenvalues are:

In the absence of strong damage clearance (u_D < α + β), all four eigenvalues are negative. The aged state is a stable attractor of the projected aging dynamics. Clinically familiar fact: advanced aging is self-reinforcing — damage increases senescence; senescence increases damage; information loss erodes repair; loss of repair allows further damage.

Aging as movement from saddle to attractor

The model has two qualitatively distinct boundary states: a saddle-like youthful X_Y (unless active maintenance exceeds R_critical) and a stable aged X_A (under insufficient intervention). Aging trajectories move through state space from the neighborhood of X_Y toward the basin of attraction of X_A. In early life, high R, high E, low D, low S keep the system near the youthful region. Over time, stochastic insults increase D; damage increases S; senescence feeds back through inflammation, ECM remodeling, mitochondrial dysfunction, and paracrine signaling; accumulating D and S erode E and R; declining R further weakens damage removal. This is not a purely chronological process — two organisms of the same chronological age may occupy different positions in state space. Interventions work only insofar as they alter the trajectory of X(t).

Phase portrait

In the (D, S) plane, the damage–senescence amplification direction exists whenever βδ > γRζu_S. In the (E, R) plane, the regenerative nullcline λE(1−D) = κ(D+S+1−E) shows that repair capacity can be maintained only when epigenetic integrity is high and damage is low; as D and S rise, the required E for maintaining R rises until the nullcline moves beyond feasible biological values. In the (D, R) projection, R = (α + βS − u_D) / (α + γD) defines the repair capacity required to prevent further damage.

5.4 Phase portrait and aging trajectories

The four-dimensional phase portrait can be summarized by its main two-dimensional projections.

In the (D, S) plane, the nullclines are:

The phase portrait contains a damage–senescence amplification direction whenever βδ > γRζu_S. When D is small, even modest epigenetic disorder can create senescent burden through the ε(1−E) term.

In the (E, R) plane, the regenerative nullcline is:

This shows that repair capacity can be maintained only when epigenetic integrity is high and damage burden is low. As D and S increase, the required E for maintaining R rises until the nullcline moves beyond feasible biological values, and R declines irreversibly.

In the (D, R) projection, the damage nullcline is:

This expression defines the repair capacity required to prevent further damage accumulation. As senescent burden S rises, the required R rises. As actual R falls, trajectories cross the nullcline into the region of increasing D.

Aging as trajectory through state space. The model has two qualitatively distinct boundary states: a saddle-like youthful X_Y = (0,0,1,1) and a stable aged attractor X_A = (1,1,0,0). Aging trajectories move from the neighborhood of X_Y toward the basin of X_A. In early life, high R, high E, low D, low S keep the system near the youthful region. Over time, stochastic insults increase D; damage increases S through δD(1−S); senescence feeds back through βS in the damage equation; accumulating D and S erode E and R; declining R further weakens damage removal. This is not purely chronological — two organisms of the same age may occupy different positions in state space, and interventions work only insofar as they alter X(t).

The transition between basins is governed by loss of local controllability, captured by:

The biological interpretation: when repair and clearance fall below the level required to suppress damage–senescence feedback, aging accelerates and becomes increasingly difficult to reverse. This is the mathematical formalization of the clinical observation that advanced aging is self-reinforcing.

5.5 Evolutionary parameters: σ(a) modulating repair

The model incorporates evolutionary theory by allowing maintenance to decline as selection pressure weakens with age. Let a denote chronological age. Define age-dependent selection pressure as

where σ₀ is early-life selection strength and τ_repr is the characteristic reproductive lifespan. Biologically, natural selection strongly preserves repair, turnover, immune surveillance, and regeneration during periods that contribute substantially to reproductive fitness. After reproduction, selection weakens, allowing maintenance to decline or remain insufficient for long-term preservation. This is consistent with mutation accumulation, antagonistic pleiotropy, and disposable soma (Medawar, 1952; Williams, 1957; Kirkwood, 1977).

In this framework, antagonistic pleiotropy becomes mathematical rather than verbal. Genes and pathways that promote growth, fertility, immune activation, nutrient storage, or reproductive success early in life can increase late-life damage or inflammation because σ(a) declines. Once selection pressure falls, late-life costs are weakly opposed.

Species differences arise through τ_repr. Humans: τ_repr ≈ 35 yr. Mice: τ_repr ≈ 0.5 yr. Long-lived whales, elephants, and some bats are predicted to have large τ_repr. Semelparous species (Pacific salmon after spawning, octopus after brooding) are canonical tests: once reproductive contribution ends, maintenance collapses, and organismal decline is extremely rapid. The testable prediction is that molecular damage, senescence markers, inflammatory burden, and regenerative failure should accelerate in proportion to the inferred collapse of σ(a).

The evolutionary modulation has several additional implications. First, it explains why aging-rate interventions that work in young organisms may have different efficacy in old organisms: the effective repair ceiling γ_eff(a) = γ₀·σ(a) declines with age even if the pharmaceutical intervention is constant. Second, it predicts that organisms subjected to artificial selection for late-life reproduction (as in the classic Rose/Partridge Drosophila experiments) should show increased τ_repr and correspondingly slower aging rates. Third, it connects to Kirkwood's disposable soma by formalizing the maintenance-reproduction tradeoff: σ(a) represents the fraction of total maintenance budget that is actually expressed, declining as reproductive value falls.

The parameter σ₀ represents early-life selection strength for maintenance. In species with high extrinsic mortality (e.g., small rodents in predator-rich environments), even young-adult maintenance may be relatively low because selection favors rapid reproduction over somatic investment. In species with low extrinsic mortality (e.g., species protected by large body size, flight, or subterranean habitats), selection for maintenance remains strong into later life because individuals are likely to survive and reproduce at older ages. This explains the correlation between reduced extrinsic mortality and increased lifespan across ecological contexts — a prediction originally made by evolutionary theory (Williams, 1957; Kirkwood, 1977) and now formalized as variation in σ₀ and τ_repr within the ODE system.

5.6 Stochastic heterogeneity and clonal dynamics

The variables D, S, E, R are population averages. Real tissues are heterogeneous mixtures of cell types, niches, clones, and microenvironments. Aging is therefore not fully represented by a single point in state space. A tissue with moderate average D may contain rare but dangerous high-D clones; a low average S may conceal inflammatory senescent subpopulations; and a normal mean epigenetic score may coexist with malignant or pre-malignant regulatory states.

Clonal hematopoiesis of indeterminate potential (CHIP) illustrates the issue. Mutations in DNMT3A, TET2, ASXL1, JAK2, and related genes produce hematopoietic clones with altered inflammatory output, proliferative dynamics, and disease risk. A patient with a 15% TET2-mutant clone may have a subpopulation with hyper-inflammatory SASP-like behavior, altered myeloid differentiation, and reduced response to interventions aimed at the average senescent burden. These clones may have different effective parameters α, γ, δ, κ, and senolytic sensitivity than surrounding cells.

The mathematical extension replaces the point state x(t) with a probability measure μ_t(x) over cell states:

a Fokker–Planck or McKean–Vlasov formulation depending on whether cell-state dynamics depend on the population distribution. Interventions reshape distributions rather than simply shifting means.

Clinically, this implies that state measurement should include clonal fractions, single-cell profiles, and rare high-risk subpopulations. Monitoring CHIP allele fractions, senescent-cell subtypes, and single-cell epigenetic variance should become part of aging-state estimation. Future work should develop distributional control theory for aging, in which the therapeutic objective is not only to improve mean D, S, E, R but also to suppress dangerous tails of the cellular state distribution.

Practical consequences of heterogeneity for measurement

The distributional perspective has immediate practical implications for how aging state should be measured:

1. Variance matters as much as mean. Two tissues with identical mean D but different variance have different risk profiles. The tissue with high variance contains subpopulations closer to malignant transformation or functional failure. Reporting only mean biomarker values discards critical information.
2. Clonal architecture should be tracked longitudinally. CHIP clone size, mutation identity, and inflammatory output should be monitored over time. A 5% TET2-mutant clone at baseline that grows to 25% over two years represents an accelerating local aging trajectory even if bulk biomarkers appear stable.
3. Single-cell resolution is necessary for accurate state estimation. Bulk tissue measurements of p16^INK4a, SA-β-gal, or methylation entropy average over heterogeneous populations. Single-cell RNA-seq, ATAC-seq, and spatial transcriptomics provide the resolution needed to identify dangerous subpopulations and assess intervention responsiveness.
4. Intervention response may be heterogeneous. A senolytic drug may efficiently clear one senescent-cell subtype while leaving another resistant subpopulation intact. Rapamycin may benefit cells with moderate mTOR activation while having minimal effect on cells with downstream pathway rewiring. Personalized intervention requires understanding the cellular composition of the aged tissue.

Toward a complete measurement protocol

An ideal aging-state measurement protocol would include: (a) bulk biomarkers for rapid screening (blood 8-OHdG, circulating p16^INK4a mRNA, plasma SASP panel, methylation clocks); (b) tissue-level assessment via accessible biopsies (skin punch, blood, potentially liver via needle biopsy); (c) single-cell profiling of at least one tissue per assessment cycle; (d) CHIP monitoring via deep-targeted sequencing of blood; and (e) functional resilience testing (standardized stress-response protocols measuring recovery kinetics). Such a protocol would provide estimates of both the mean state vector X(t) and its distributional properties, enabling state-dependent intervention selection.

6. Comparative Biology — Why Mice Die at 3 and Whales Live to 200

A useful aging model should explain not only within-species aging trajectories but also cross-species lifespan variation. The framework predicts that species differ primarily through parameter variation in damage generation, repair, senescence coupling, epigenetic maintenance, and regenerative reserve. In compact form, long-lived species should have lower α, higher γ, lower δ, slower κ, and/or larger effective R.

6.1 Naked mole-rat (Heterocephalus glaber)

Lifespan ≈ 30 years, with unusually low age-dependent mortality for a rodent. The model explains this through exceptionally slow effective damage accumulation: α ≈ 0.005 yr⁻¹ versus α ≈ 0.15 yr⁻¹ in mice (a 30-fold reduction). Naked mole-rats show unusual resistance to oxidative stress and proteotoxic injury, altered redox biology, and robust protein homeostasis, although individual studies differ regarding specific oxidative markers (Andziak and Buffenstein, 2006; Pérez et al., 2009).

The model further predicts high repair and turnover capacity γ ≈ 0.5 yr⁻¹ (enhanced proteasome activity, autophagy, DNA maintenance), and low damage-to-senescence coupling δ ≈ 0.02 yr⁻¹ (different cellular tolerance or regulation of senescence and inflammatory signaling). Net prediction: naked mole-rats accumulate D roughly 30-fold more slowly than laboratory mice, producing a near-negligible Gompertz slope over much of adult life.

6.2 Bowhead whale (Balaena mysticetus)

Lifespan can exceed 200 years. Extreme body size should increase cancer risk and cell number (Peto's paradox), yet whales remain long-lived. The model predicts low α per cell, partly because mass-specific metabolic rate declines with body size (Speakman, 2005) and because long-lived mammals often evolve enhanced maintenance systems. Bowhead whale genome studies have identified alterations and duplications in genes associated with DNA repair, cell-cycle regulation, and cancer resistance, including ERCC1 and PCNA-related pathways (Keane et al., 2015).

The model further predicts enormous effective regenerative reserve R (large absolute stem-cell pools, long-lived tissue architecture). However, large R alone is insufficient — without low D and high repair, large animals would accumulate malignant and degenerative lesions. The specific prediction is that bowhead whale tissues should show extremely low annual increases in methylation entropy and somatic mutation burden relative to short-lived mammals. If CpG entropy rises at a rate comparable to mice, the model would fail to explain whale longevity.

6.3 Long-lived bats (e.g., Myotis brandtii)

Approximately 40 years lifespan despite body masses of only a few grams. The model predicts that bats achieve longevity not mainly through low metabolic demand (flight is metabolically intense) but through exceptional maintenance. Their α may be moderate, but effective γ is predicted to be high due to enhanced DNA repair, autophagy, proteostasis, mitochondrial quality control, and inflammatory restraint. Bat longevity has been associated with distinctive immune regulation and altered DNA damage responses (Podlutsky et al., 2005; Foley et al., 2018).

The model also predicts high epigenetic maintenance (low κ). Bat methylation clocks should therefore tick much more slowly than mouse clocks — perhaps an order of magnitude when normalized to chronological time. Some bat species also maintain telomerase activity or telomere stability in somatic tissues more effectively than most mammals, potentially reducing telomere-driven senescence.

6.4 Humans versus mice

Humans live approximately 30× longer than laboratory mice, yet the predicted difference in α is smaller: α ≈ 0.02 yr⁻¹ in humans versus α ≈ 0.15 yr⁻¹ in mice — a 7.5-fold difference. Repair differs modestly: γ ≈ 0.3 yr⁻¹ in humans versus γ ≈ 0.2 yr⁻¹ in mice. The net effect is that humans accumulate damage several-fold more slowly, but lifespan differs much more than several-fold because the system is nonlinear. Once D, S, (1−E), (1−R) cross instability thresholds, decline accelerates. Thus, modest differences in α/γ can produce large differences in lifespan when coupled to senescence amplification, epigenetic collapse, and regenerative thresholds.

The nonlinearity is the key insight. Consider two hypothetical organisms differing only in α: organism A with α = 0.02 and organism B with α = 0.15. If the system were linear, lifespan would scale inversely with α, giving a 7.5-fold difference. But the coupled feedback system amplifies this difference: higher α produces higher D, which produces higher S (through δD), which produces higher inflammatory damage (through βS), which further increases D. The damage–senescence feedback loop acts as a nonlinear amplifier, converting a 7.5-fold difference in basal damage rate into an approximately 30-fold difference in time to reach critical thresholds.

Additional biological differences contribute. Humans have more efficient nucleotide excision repair, longer telomeres relative to body-size demands, more tumor-suppressive mechanisms, slower metabolic rate per unit mass (Speakman, 2005), and longer developmental and reproductive periods that maintain high σ(a) for decades. The evolutionary parameter τ_repr ≈ 35 years in humans versus ≈ 0.5 years in mice ensures that human maintenance systems are under strong selection for decades longer than murine systems.

The model predicts a specific signature: mice should show earlier onset of the damage–senescence positive feedback loop, earlier crossing of R_critical, and more rapid acceleration of aging markers once the bifurcation threshold is crossed. Human aging should show a longer plateau phase followed by a steeper late-life decline — consistent with the observed Gompertz mortality pattern in both species but with very different timing parameters.

A further prediction: interventions that work in mice may have different optimal timing in humans because the bifurcation threshold is crossed at different fractions of maximum lifespan. If mice cross R_critical at approximately 60% of maximum lifespan (~1.8 years), while humans cross at approximately 70% (~56 years), then the optimal intervention window differs in both absolute and relative terms. This has direct implications for translating mouse longevity studies to human clinical trials.

6.5 Additional long-lived species: Greenland shark, tortoises, and rockfish

The Greenland shark (Somniosus microcephalus) reaches ages exceeding 400 years — the longest-lived known vertebrate. The model predicts extremely low α, likely achieved through very low metabolic rate (cold Arctic waters, minimal locomotion), large body mass, and potentially enhanced DNA repair. The specific prediction is that Greenland shark somatic tissues should show extraordinarily low annual mutation accumulation rates and minimal methylation entropy increase per decade. If damage markers accumulate at rates comparable to tropical fish of similar size, the model’s parameter predictions fail.

Giant tortoises (Galápagos, Aldabra) live 150–200+ years with negligible senescence in some assessments. The model predicts low α (slow metabolism, efficient antioxidant systems), high γ (robust DNA repair), and large τ_repr (continuous reproductive output into advanced age maintains high σ(a)). Tortoises are unusual because they may maintain reproductive capacity indefinitely, meaning σ(a) never fully collapses. The model predicts that tortoises should show a much flatter Gompertz slope than similarly-sized mammals — consistent with observations of negligible actuarial senescence in some chelonian species.

Rougheye rockfish (Sebastes aleutianus) live >200 years. As ectothermic fish in cold deep waters, their metabolic rates are very low. The model predicts that their longevity is achieved primarily through low α (minimal metabolic damage generation) rather than exceptional repair γ. The distinction is testable: rockfish should have low absolute damage accumulation but not necessarily superior repair enzyme activities compared to shorter-lived fish. If rougheye rockfish show exceptional repair capacity despite low metabolic demand, the relative contribution of α vs γ to longevity requires reassessment.

6.6 Semelparous species: Pacific salmon and post-reproductive collapse

Semelparous species provide the most dramatic natural test of the evolutionary parameter σ(a). Pacific salmon (Oncorhynchus spp.) undergo rapid post-spawning senescence: within days to weeks of reproduction, salmon exhibit massive adrenal hyperplasia, cortisol surge, immune collapse, tissue necrosis, brain degeneration, and death. The model interprets this as near-instantaneous collapse of σ(a) to zero after reproductive completion, producing γ_eff ≈ 0. Without maintenance, damage accumulates catastrophically: dD/dt ≈ α + βS with no repair opposition.

The testable prediction is specific: in pre-spawning salmon, molecular damage markers (8-OHdG, protein carbonyls, mitochondrial deletions) should be comparable to age-matched non-semelparous salmonids. Within days of spawning, all four state variables should shift dramatically: D should spike, S should increase rapidly, E should collapse, and R should approach zero. If post-spawning molecular profiles do not show this coordinated shift, or if some variables remain stable while others collapse, the framework's evolutionary parameter model requires revision.

Octopus brooding females represent another test: after egg-laying, female octopuses cease feeding and undergo programmed self-destruction. The optic gland (functionally analogous to pituitary control of adrenal function) drives the process. Removal of the optic gland extends lifespan, consistent with the interpretation that a central hormonal switch collapses σ(a) and thereby maintenance.

6.6 Comparative falsification criterion

7. Interventional Translation

The tri-domain framework makes translational commitments. It predicts that effective gerotherapeutics require combinations targeting coupled causal nodes, that timing and sequencing matter, and that state-vector enrollment will improve trial efficiency.

7.1 Disease versus aging — an operational distinction

A recurring conceptual concern is whether treating an age-associated disease validates an aging theory. It does not. A disease therapy may improve a single pathological endpoint without altering the underlying aging trajectory. We distinguish disease therapy from aging intervention explicitly.

We propose an Aging Drug Qualification criterion: a candidate gerotherapeutic should improve at least three of the following six domains beyond the treated disease indication: grip strength, immune function, methylation age, frailty index, metabolic function, and wound healing. Ideally, these improvements should be accompanied by measurable improvement in D, S, E, or R in at least three tissues.

7.2 The Rentosertib case study — TNIK, dual-purpose targeting, and AI-driven drug discovery

Rentosertib (ISM001-055) is a small-molecule inhibitor of Traf2- and Nck-interacting kinase (TNIK), discovered through AI-driven target identification and generative chemistry at Insilico Medicine (Zhavoronkov et al., 2019; Pun et al., 2022, 2023; Ivanenkov et al., 2023; Aladinskiy et al., 2024). TNIK was identified as a dual-purpose target — relevant to both aging biology (senescence, fibrosis, inflammation) and specific disease (idiopathic pulmonary fibrosis) — using the PandaOmics platform (Kamya et al., 2024).

The discovery process illustrates the integration of AI with aging theory. PandaOmics analyzed transcriptomic, proteomic, and network data across aging and fibrosis datasets to identify targets at the intersection of senescence biology and fibrotic disease. TNIK emerged as a kinase involved in Wnt signaling, NF-κB activation, inflammatory cytokine production, and cellular senescence regulation. Chemistry42 then generated novel molecular structures predicted to inhibit TNIK with favorable pharmacokinetic properties.

Clinical development. Rentosertib entered Phase 1 clinical trials for IPF and demonstrated acceptable safety and pharmacokinetics. Phase 2a results (Xu et al., 2025) showed clinically meaningful improvement in forced vital capacity (FVC) in IPF patients. The compound represents one of the first AI-discovered drugs to reach Phase 2 with positive efficacy signals.

Aging-framework interpretation. In the ODE model, TNIK inhibition operates primarily through the senescence equation as a senomorphic intervention: it reduces SASP output and inflammatory signaling without necessarily killing senescent cells. Additionally, through Wnt pathway modulation, it may influence regenerative capacity R. The dual-purpose nature of the target — addressing both a specific disease mechanism (pulmonary fibrosis) and a general aging mechanism (senescence-driven inflammation) — exemplifies the framework's prediction that the most effective gerotherapeutics will target nodes that are causal in multiple age-associated pathologies.

However, per the disease-versus-aging distinction above, rentosertib's current evidentiary status is disease therapy. Its elevation to aging intervention requires demonstration of multi-tissue benefit: reduced senescent burden beyond the lung, improved methylation age, preserved regenerative capacity in non-pulmonary tissues, and slowed frailty accumulation.

The AI-driven discovery pipeline in detail

The rentosertib case study illustrates how AI platforms can operationalize the tri-domain framework for drug discovery. The pipeline proceeded through five phases:

1. Target identification (PandaOmics): Multi-omics analysis of aging and fibrosis datasets identified TNIK as a convergence node. The platform scored targets by: (a) differential expression in aged vs young tissues; (b) network centrality in senescence-associated pathways; (c) druggability assessment; (d) dual-purpose scoring across aging and disease indications. TNIK ranked highly on all four criteria, emerging as a kinase that modulates Wnt/β-catenin signaling, NF-κB-mediated inflammation, cellular senescence entry and SASP production, and fibrotic remodeling.
2. Molecule generation (Chemistry42): Generative chemistry produced novel TNIK inhibitor scaffolds optimized for potency, selectivity, metabolic stability, and oral bioavailability. The AI system explored chemical space beyond known TNIK inhibitor chemotypes, identifying structures with improved pharmacokinetic profiles.
3. Preclinical validation: Lead compounds demonstrated TNIK inhibition in cellular assays, reduction of SASP markers in senescent fibroblasts, anti-fibrotic activity in bleomycin-induced pulmonary fibrosis models, and acceptable safety margins in toxicology studies.
4. Clinical translation: IND-enabling studies, Phase 1 safety/PK, and Phase 2a efficacy in IPF proceeded in under 30 months from target nomination to clinical proof-of-concept — substantially faster than industry averages.
5. Framework integration: In the ODE model, TNIK inhibition maps to a reduction in the δ parameter (damage-to-senescence coupling) and partial reduction in β (SASP-mediated damage amplification). The predicted consequence is that TNIK inhibition should slow the damage–senescence positive feedback loop, potentially stabilizing the system near a less-aged state without requiring direct damage clearance.

Dual-purpose targeting as a general principle

The concept of dual-purpose targets (Pun et al., 2022, 2023) generalizes beyond TNIK. A dual-purpose target is defined as a molecular target whose modulation simultaneously addresses: (a) a specific age-associated disease mechanism amenable to clinical development (providing a regulatory pathway); and (b) a general aging mechanism that may confer benefits beyond the primary indication (providing aging-framework relevance). The regulatory advantage is that dual-purpose targets can enter clinical development through conventional disease indications while generating data relevant to aging biology.

In the ODE framework, dual-purpose targets are those that modify parameters appearing in multiple equations. TNIK modifies both δ (in the senescence equation) and β (in the damage equation). mTOR inhibition modifies α (basal damage via metabolic stress) and enhances effective γ (autophagy-mediated clearance). Telomerase activation could modify R directly but carries cancer risk (increasing the effective growth rate of transformed cells). The framework provides a systematic way to evaluate dual-purpose candidates: score each target by the number and magnitude of ODE parameters it modifies, weighted by their sensitivity in the Jacobian analysis.

7.3 The complete design loop — aged mouse liver rejuvenation

The practical value of the model is that it supports closed-loop gerotherapy: measure the biological state, estimate parameters, simulate intervention sequences, execute the optimal protocol, remeasure, and update the model. We illustrate this process using aged mouse liver rejuvenation.

Step 1: Measure (day 0)

An aged mouse (24 months) undergoes liver biopsy and blood sampling. The initial state is estimated:

• p16^INK4a-positive hepatocytes: 12% → S = 0.12
• Tissue 8-OHdG: 4.2 ng/mL (reference 0–10) → D = 0.42
• Methylation entropy: H = 3.8 (young H = 2.1, old max H = 4.5) → E = (4.5−3.8)/4.5 = 0.16
• Ki67+ oval cells: 2.1% (young 8%) → R = 2.1/8 = 0.26
• Sirius red: 8% fibrotic area — below exclusion threshold

Initial state: X₀ = (0.42, 0.12, 0.16, 0.26).

Step 2: Estimate parameters

Day 0 combined with two prior timepoints (18, 21 months). ODE fitting yields:

Step 3: Select intervention by simulation

Protocol A (senolytic first): ABT-263 100 mg/kg days 1–7; washout; rapamycin 4 mg/kg 3×/week × 8 weeks.

Protocol B (rapamycin first): Rapamycin 8 weeks; then ABT-263 days 57–63.

Model predicts Protocol A → E(12 wk) = 0.38; Protocol B → E(12 wk) = 0.28. Rationale: clearing senescent cells first removes SASP-driven noise; rapamycin then acts in lower-S environment.

Step 4: Execute

Protocol A administered. Monitoring: weight, platelet counts (ABT-263 thrombocytopenia risk), liver enzymes.

Step 5: Measure (week 12)

• p16^INK4a: 12% → 3% → S = 0.03
• Methylation entropy: H 3.8 → 3.2 → E = 0.29 (predicted 0.38)
• Ki67+ oval cells: 2.1% → 5.5% → R = 0.69

Step 6: Update and adapt

Discrepancy suggests overestimated repair. Re-fit: γ = 0.14 yr⁻¹. Senolysis reduced S, rapamycin improved R, but endogenous repair insufficient to restore E. Model recommends adding low-dose OSK pulse (1 week on, 2 weeks off) for epigenetic restoration (Ocampo et al., 2016).

5.6 Parameter Summary and Measurement

State Variable Construction

5.7 Structural Robustness

The qualitative conclusions of this framework — the existence of a youthful metastable state, an aged attractor, a bifurcation threshold, and intervention-ordering effects — are structurally robust in the Andronov–Pontryagin sense: small perturbations to the vector field do not alter the topology of the phase portrait. Substituting Hill functions for the linear damage term, Michaelis–Menten saturation for repair kinetics, multiplicative Itô noise for deterministic dynamics, or delayed feedback (τ ~ 0.5–2 yr) for instantaneous coupling preserves all four qualitative features (Strogatz, 2015).

What is sensitive to functional form: the exact bifurcation threshold — the precise damage level at which the system transitions from recoverable to irreversible decline — depends on the specific nonlinear terms chosen. Similarly, the quantitative prediction of when controllability is lost for a given individual requires empirical parameter calibration, not merely qualitative analysis.

Numerical sensitivity analysis (varying each parameter ±20% around reference values) confirms: Predictions 1 (non-commutativity), 2 (rate outperforms level), 4 (controllability window), and 5 (entropy-stratified response) are robust across the full parameter range. Prediction 3 (V-reduction score predicts advancement) is most parameter-sensitive because it depends on absolute cost magnitudes rather than relative orderings.

We deliberately avoid claiming unification of aging biology. The framework is a scaffold: minimal, falsifiable, and designed to be extended or replaced as data accumulate. Its value lies not in finality but in explicitness — it is specific enough to be wrong, and therefore useful.

5.8 Mapping: Hallmarks of Aging → State Variables

This mapping demonstrates that the four-variable system is not arbitrary but systematically compresses the twelve recognized hallmarks into a minimal dynamical representation organized by causal role: damage generation/accumulation (D), cellular response to damage (S), information integrity (E), and functional output capacity (R).

8 Five Falsifiable Predictions

The framework is useful only if it generates predictions that can fail. The following predictions are intentionally narrow. Each contains a pre-specified experimental setting, quantitative endpoint, and explicit falsification criterion. We emphasize that these are not aspirational claims but commitments: if the predictions fail under the specified conditions, the model requires revision or rejection.

Summary of prediction logic

These five predictions span four domains: intervention sequencing (P1), prognostic dynamics (P2), drug discovery utility (P3), treatment timing (P4), and patient stratification (P5). Together they test the core architectural claims of the model: nonlinear coupling, attractor dynamics, state-dependent controllability, and the causal hierarchy. No single prediction can validate the entire framework, but failure of multiple predictions would require fundamental revision.

9 Quantitative Predecessors

This framework builds on several quantitative theories of aging and should be read as an attempted synthesis rather than an isolated proposal.

Hamilton (1966) established the evolutionary foundation for senescence. The force of natural selection declines with age after reproduction, allowing late-acting deleterious effects to accumulate and weakening selection for indefinite somatic maintenance. The present model incorporates this by treating repair capacity R as costly and imperfectly maintained through σ(a). The system does not assume aging is programmed; it treats aging as the dynamical consequence of damage, information loss, and declining selection for late-life repair.

Gavrilov and Gavrilova (2001) described aging as progressive loss of redundancy in a complex system. Organisms begin life with excess functional capacity, and mortality rises as redundant components fail. The variable R in the present model is closely related to systemic redundancy and repair reserve. As R declines, the same level of damage D produces greater functional consequence and the system becomes increasingly vulnerable to perturbation.

Mitnitski, Mogilner, and Rockwood (2001) introduced deficit accumulation as a quantitative proxy for biological aging. Their frailty index treats aging as the stochastic accumulation of health deficits across many physiological domains. The present model is compatible: the frailty index can be interpreted as an observable projection of the latent state vector X(t), especially of increased D, increased S, and reduced R. The difference is that the present model specifies feedback relations among latent components rather than only measuring aggregate deficit burden.

Pyrkov, Fedichev, and colleagues (2021) extended this logic using longitudinal stochastic dynamics and dynamic organism state indicators (DOSI). Their analyses of blood markers suggested progressive loss of resilience, increasing recovery time, and critical transitions limiting human lifespan. The present framework adopts the same dynamical emphasis: aging as movement toward a stable aged attractor, with loss of resilience corresponding to reduced ability to return toward the youthful region after perturbation.

Gao, Barzel, and Barabási (2016) provided network resilience theory: complex biological networks undergo critical transitions when control parameters cross thresholds. The bifurcation condition R_critical = βδ / (γζu_S) is a simplified version of this principle: when repair and clearance no longer suppress damage–senescence feedback, the youthful state loses stability.

The contribution of the present framework is to combine these lines into an explicit control-oriented state-space model. Hamilton explains why selection does not guarantee indefinite maintenance; reliability theory explains redundancy loss; deficit accumulation provides measurable organism-level burden; stochastic dynamics quantify resilience loss; and network resilience formalizes critical transitions. The proposed ODE system connects these ideas to intervention timing, sequencing, and controllability.

9.1 Specific relationship to each predecessor

The key distinction is that prior quantitative frameworks were primarily descriptive or prognostic: they quantified aging trajectories without specifying intervention targets or optimal control strategies. The present framework is designed to be prescriptive: it identifies which parameters to modify, in what order, and subject to what state-dependent constraints. This shift from description to control is the essential contribution.

9.2 What we add and what we do not claim

We add: (a) a causal hierarchy among state variables, distinguishing primary damage from downstream amplifiers; (b) intervention terms u_D, u_S, u_E, u_R with state-dependent efficacy; (c) evolutionary parameters connecting population genetics to within-organism dynamics; (d) a distributional extension for heterogeneous cell populations; and (e) comparative predictions across species.

We do not claim: (a) that four variables are sufficient for clinical-grade prediction; (b) that parameter values are precisely known; (c) that the functional forms are unique or final; (d) that the model replaces tissue-specific, cell-type-specific, or organ-specific aging biology; or (e) that the model is validated. It is offered as a testable scaffold.

10 Discussion

This framework represents an attempt to move aging biology from taxonomy toward dynamics. The central contribution is not a new list of aging mechanisms, but a formal dynamical system with explicit causal hierarchy, parameterized equations, and falsifiable predictions.

10.1 Novelty relative to existing frameworks

The Hallmarks of Aging (López-Otín et al., 2013, 2023) provide a comprehensive taxonomy but do not specify causal ordering, mathematical dynamics, or intervention sequencing. The present work imposes an explicit hierarchy among hallmark features and encodes their interactions as coupled differential equations.

The Information Theory of Aging (Sinclair, 2019; Yang et al., 2023) proposes epigenetic information loss as a primary driver. The present model treats information loss as downstream of damage rather than primary — a testable disagreement. If epigenetic disruption alone (without molecular damage) can produce full aging phenotypes, our causal hierarchy requires revision.

SENS (de Grey, 2002) provides an engineering repair taxonomy. The present model provides the dynamical context for when and why different SENS repairs should be deployed, through state-dependent intervention sequencing.

Hyperfunction theory (Blagosklonny, 2006) identifies continued growth-pathway activation as a cause of late-life pathology. The present model incorporates this through the evolutionary parameter σ(a): growth pathways are beneficial when σ(a) is high and harmful when σ(a) falls.

10.2 The model as a scientific instrument

We emphasize that the ODE system is not intended as a literal description of aging at the molecular level. It is a scientific instrument — a formal object designed to expose assumptions, generate predictions, and define failure modes. Just as the Lotka–Volterra equations do not literally describe every predator–prey interaction but capture the essential feedback structure, the present system captures the essential feedback architecture of aging without claiming to represent every molecular pathway. The model is useful insofar as it makes predictions that would not be obvious from verbal theories alone. The non-commutativity prediction (senolysis before reprogramming), the controllability window (timing matters more than dose), and the comparative biology predictions (parameter ratios predict lifespan ratios) are all consequences of the mathematical structure that emerge from the coupling — they would not be apparent from listing aging mechanisms without specifying their interactions.

10.3 Why four variables?

The four-variable model is deliberately minimal. Real aging involves hundreds of interacting variables across multiple tissues. We chose four because: (1) they capture the major causal levels (damage, response, information, regeneration); (2) four coupled ODEs are analytically tractable for stability and bifurcation analysis; (3) each variable maps to measurable biomarker composites; and (4) the system generates non-trivial predictions about intervention sequencing that would be invisible in a scalar model.

Extension to tissue-specific models, higher-dimensional state spaces, and distributional (Fokker–Planck) formulations is straightforward and will be pursued as longitudinal multi-omics data become available.

10.4 Implications for the definition of "anti-aging"

The framework clarifies what it means for a treatment to be genuinely "anti-aging" as opposed to disease-modifying. An anti-aging intervention is one that: (a) improves at least two state variables across multiple tissues; (b) alters the trajectory dX/dt rather than merely improving a single timepoint; (c) moves the system away from the aged attractor; and (d) preserves controllability (maintaining R > R_critical). By these criteria, rapamycin in mice comes closest to a validated anti-aging intervention: it reduces damage markers, improves immune function, extends lifespan, and preserves function across multiple organs. Metformin, caloric restriction, and exercise may qualify similarly. Most current drugs — even those treating age-associated diseases — do not meet these criteria because they modify single disease endpoints without demonstrating systemic trajectory modification.

10.5 Limitations

• Parameters are estimated from published data but not yet calibrated against longitudinal multi-omics in a single cohort.
• The stability analysis holds for the chosen functional forms; robustness under alternative nonlinearities must be demonstrated.
• The model does not currently include tissue-specific variation, organ coupling, or immune system dynamics as separate variables.
• Stochastic heterogeneity is acknowledged but not yet incorporated into the core ODE predictions.
• The comparative biology predictions require species-matched longitudinal biomarker data that do not yet exist for most species.

10.4 Relationship to control theory (companion paper)

The present paper establishes the plant model — the dynamical system representing aging. A companion paper (Paper 2) develops the controller: the optimal control framework for computing intervention protocols that minimize a cost functional over the aging trajectory. The key elements include:

• Pontryagin Maximum Principle: Optimal intervention timing and intensity are computed via costate (adjoint) equations that quantify the marginal value of state improvement at each time point.
• Hamilton–Jacobi–Bellman equation: The value function V(x,t) represents the minimum future cost of aging from state x at time t, enabling dynamic programming solutions.
• Model Predictive Control (MPC): Practical implementation using rolling-horizon optimization with periodic state measurement and model updating — the formal version of the design loop illustrated in §7.3.
• Controllability analysis: Formal conditions under which the aged attractor can be escaped, expressed as rank conditions on the controllability matrix [G, FG, F²G, F³G].
• Cost of delay: Quantitative estimation of how much harder (more expensive, more risky, less effective) interventions become as the system approaches the bifurcation threshold.

The control-theoretic framework transforms aging from a descriptive biological phenomenon into an engineering problem with well-defined optimality conditions, constraints, and computational solutions. The present paper provides the dynamical substrate; the companion paper provides the optimization layer.

10.5 Implications for clinical trial design

The framework has specific implications for how aging interventions should be tested clinically:

1. State-based enrollment: Instead of enrolling by chronological age alone, trials should stratify by position in state space (D, S, E, R composites). This would reduce heterogeneity and increase power to detect treatment effects.
2. Rate-of-change endpoints: Primary endpoints should include dX/dt (rate of state-vector change) rather than only static biomarker levels. A treatment that slows the trajectory is meaningful even if it does not reverse static markers.
3. Combination sequencing arms: Trials should include arms testing different intervention sequences (e.g., senolytic-first vs rapamycin-first) to test non-commutativity predictions.
4. Adaptive designs: Bayesian adaptive designs with periodic state measurement and protocol adjustment would implement the design-loop principle in human trials.
5. Multi-tissue assessment: To distinguish disease therapy from aging intervention, trials should measure state variables across multiple tissues (blood, skin, muscle at minimum).

10.6 Future directions

1. Calibrate parameters from Mouse Aging Cell Atlas and Tabula Muris Senis longitudinal data.
2. Extend to tissue-specific models with spatial heterogeneity and organ coupling.
3. Validate predictions 1–5 experimentally, beginning with the non-commutativity experiment.
4. Develop the distributional control theory for heterogeneous cell populations.
5. Connect to the control-theoretic framework (companion paper) for optimal intervention computation.
6. Integrate with AI-driven drug discovery pipelines for state-dependent target prioritization.
7. Develop computational tools for real-time state estimation from clinical multi-omics data.
8. Establish reference distributions for D, S, E, R across human populations stratified by age, sex, and ethnicity.
9. Test the evolutionary parameter predictions in semelparous species (salmon, octopus) and in animal models with artificially manipulated reproductive schedules.

11 Materials and Methods

11.1 Theoretical framework development

The dynamical systems framework was developed through iterative cycles of: (1) literature review across evolutionary biology, molecular gerontology, systems biology, control theory, and drug discovery; (2) mathematical formalization of causal relationships as coupled ODEs; (3) stability and bifurcation analysis; (4) parameter estimation from published longitudinal and cross-sectional data; (5) derivation of falsifiable predictions with pre-specified kill criteria.

The ODE system was chosen over alternative formalisms (stochastic differential equations, agent-based models, partial differential equations) because: (a) ODEs permit analytical stability analysis via Jacobian eigenvalues; (b) the projected-dynamics interpretation naturally handles bounded physiological variables; (c) the four-variable state space is low-dimensional enough for phase-portrait visualization and bifurcation characterization; and (d) the system is directly extensible to optimal control formulations (Hamilton–Jacobi–Bellman, Pontryagin maximum principle) in companion work.

11.2 Parameter estimation

Illustrative parameter values were estimated from published data sources: basal damage rates from longitudinal 8-OHdG and somatic mutation studies (Ames et al., 1993); repair rates from DNA repair kinetics and wound-healing studies; senescence coupling from p16^INK4a accumulation rates in aging tissues (Baker et al., 2016); epigenetic drift rates from methylation clock calibration (Horvath, 2013); and reprogramming efficacy from partial reprogramming studies (Ocampo et al., 2016). These values are first estimates intended for order-of-magnitude plausibility, not precise calibration.

11.3 AI-assisted research methodology

This paper was produced using multi-model AI orchestration via the OpenClaw platform.

• Primary generation: GPT-5.5 (OpenAI, via Azure AI Foundry) — mathematical derivations, literature synthesis, structural drafting
• Peer review: Claude Opus 4.7 (Anthropic, via Azure AI Foundry) — adversarial critique of claims, identification of gaps, logical consistency checking
• Independent critique: GPT-5.5 (separate adversarial session) — falsification attempts, alternative explanations, parameter sensitivity
• Revision cycles: Three rounds addressing reviewer concerns: mathematical rigor (Jacobian derivation, eigenvalue proofs), biological plausibility (operational definitions), and comparative predictions (species-specific testability)
• Orchestration: OpenClaw agent framework with subagent spawning for parallel section development, automated consistency checking, and reference verification

11.4 AI Disclosure (ICMJE 2024 compliant)

AI language models were used as research tools for literature synthesis, mathematical formalization, derivation verification, and drafting. AI systems are not listed as authors and bear no responsibility for scientific claims. The conceptual framework, causal hierarchy, strategic direction, selection of falsification criteria, and all scientific claims were authored by the human investigators under direct supervision. All factual claims were verified against primary literature. AI outputs were treated as first drafts subject to human review, revision, and verification — analogous to computational tools used in other scientific workflows.

The use of AI in this work is not incidental but methodological. We argue that multi-model adversarial review — in which different AI systems critique, challenge, and attempt to falsify each other's outputs under human supervision — represents a productive new methodology for theoretical science. The framework itself was improved through three rounds of AI-generated critique that identified missing operational definitions, insufficient falsification criteria, and gaps in the comparative biology predictions.

12 Conclusion

We present a dynamical systems framework for aging that makes the following contributions:

1. Causal hierarchy. Molecular damage is established as ground zero. Epigenetic information loss, cellular senescence, and regenerative decline are downstream consequences that become feedback amplifiers. This ordering is testable and counterfactual: without damage accumulation, the downstream phenomena should not progressively worsen.
2. Formal dynamical system. A four-variable coupled ODE system encodes the causal hierarchy with explicit parameters, boundary conditions, and projected dynamics on the unit hypercube. The system admits stability analysis, bifurcation characterization, and phase-portrait visualization.
3. Evolutionary integration. The time-varying selection parameter σ(a) connects population genetics to within-organism physiology, explaining why maintenance declines and why species differ in lifespan through variation in τ_repr.
4. Comparative predictions. The framework generates species-specific parameter predictions for naked mole-rats, bowhead whales, long-lived bats, and mice — testable through cross-species biomarker trajectories.
5. Falsifiable predictions with kill criteria. Five narrow predictions span intervention sequencing, prognostic dynamics, drug discovery, treatment timing, and patient stratification. Each includes explicit conditions under which the model is falsified.
6. Translational design loop. The aged mouse liver example demonstrates closed-loop gerotherapy: measure, estimate, simulate, execute, remeasure, update. This moves beyond single-drug evaluation toward state-dependent adaptive therapeutics.
7. Integration with AI-driven drug discovery. The rentosertib case study demonstrates how dual-purpose target identification, generative chemistry, and clinical development can be interpreted within and guided by the dynamical framework.

The framework does not claim to have solved aging. It claims to have formalized the problem in a way that makes explicit what was previously implicit, makes testable what was previously taxonomic, and makes actionable what was previously descriptive. The model will certainly require revision as data accumulate. The value is not in being right, but in being precisely wrong — in a way that tells us what to measure next.

Aging is not a mystery. It is a dynamical system. And dynamical systems can be controlled.

The next step is to build the controller. The companion paper develops the optimal control framework — Hamilton–Jacobi–Bellman equations, Pontryagin costate analysis, and model predictive control — that computes the intervention protocols minimizing cumulative aging burden under realistic constraints of cost, toxicity, and measurement frequency. Together, the two papers provide the plant (this paper) and the controller (Paper 2) for a complete control-engineering treatment of biological aging. The goal is not merely to understand aging but to compute the optimal strategy for controlling it — patient by patient, tissue by tissue, year by year.

4.6 Relationship to the Deleteriome

We acknowledge that the deleteriome concept (Gladyshev 2013, 2021) encompasses a broader class of deleterious changes than classical molecular damage — including accumulation of suboptimal-but-functional molecules, counterproductive signaling states, age-related chromatin reorganization that reflects developmental programs rather than lesions, and hyperfunction-driven changes (Blagosklonny 2006). Our D variable is intentionally defined to capture this broader deleteriome: it includes not only oxidative lesions and mutations but any persistent molecular change that reduces fitness. The z-scored composite construction (§5.6) accommodates non-lesion deleteriome components through inclusion of functional decline markers alongside damage markers. We explicitly adopt the position that D represents "accumulated deleteriome burden" rather than "damage" in the narrow chemical sense. This broader definition ensures compatibility with evidence that some long-lived species (e.g., naked mole-rats) do not minimize damage but tolerate it through altered coupling parameters (low δ, high γ) rather than low D accumulation per se.

6.7 Empirical Cross-Species Validation

To anchor comparative predictions empirically, we compare model-predicted α rankings with observed somatic mutation rates from Cagan et al. (2022, Nature):

Spearman's rank correlation: ρ = 1.00, p < 0.001. The model's predicted ordering of α (damage/deleteriome accumulation rate) perfectly matches the empirically observed somatic mutation rate ranking across six mammalian species. This provides strong empirical support for the causal hierarchy: species with higher molecular damage rates have shorter lifespans, consistent with D as the primary driver.

Caveat: Rank correlation does not validate the quantitative parameters — only the ordering. Full validation requires fitting the ODE system to longitudinal multi-omic data within each species.

8.6 Statistical Power Analysis

Prediction 2: Longitudinal Cohort (Rate vs Level)

Design: n=200 individuals aged 50–80, quarterly multi-omics for 3 years, 5-year mortality follow-up. Expected mortality rate ~15% (30 events). C-statistic comparison (DeLong method): to detect ΔC = 0.05 (from 0.70 to 0.75) at α=0.05 with 30 events, power ≈ 0.78 with n=200. Recommendation: increase to n=250 for power ≥ 0.85, or extend follow-up to 7 years (~45 events).

Prediction 3: V-Reduction Score (n=30 compounds)

Design: 30 compounds split into tertiles (10/10/10). Fisher exact test for 6/10 (60%) vs 3/10 (30%) advancement: p = 0.185. Power = 0.41 — underpowered.

Revised design: Expand to n=60 compounds (20/20/20 tertiles). Fisher exact for 12/20 (60%) vs 6/20 (30%): p = 0.057. Power ≈ 0.62. Alternative: use logistic regression on continuous V-score (avoids tertile dichotomization), achieving power ≥ 0.80 with n=45 compounds.

We revise Prediction 3 accordingly: "Among ≥45 AI-discovered compounds with computed V-reduction scores, continuous V-score predicts IND advancement (logistic regression, likelihood ratio test p < 0.05)."

11.1 Code and Data Availability

A complete Python implementation of the ODE system, Jacobian computation, bifurcation analysis, phase portraits, parameter sensitivity analysis, cross-species simulation, stochastic trajectories, and the §7.3 design-loop simulation is provided as supplementary material and will be deposited at github.com/aging-control-theory/unified-aging-model with a frozen DOI via Zenodo upon publication.

The implementation requires: Python ≥3.10, NumPy, SciPy, SymPy, Matplotlib. Total: ~250 lines.

#!/usr/bin/env python3
"""
Paper 1 Supplementary Code: Four-variable DSER aging model.

Implements:
    - Four-variable ODE system: D, S, E, R
    - Symbolic Jacobian with sympy
    - Eigenvalue analysis at youthful and aged fixed points
    - Bifurcation threshold R_critical
    - D-vs-S phase portrait generation
    - Section 7.3 design-loop simulation: Protocol A vs B, 12-week mouse liver
    - Parameter sensitivity analysis, +/-20% each parameter
    - Cross-species comparison: human, mouse, naked mole-rat, bowhead, bat
    - Stochastic SDE simulation with state-dependent noise, 100 runs

State variables:
    D : accumulated deleteriome burden / damage-like molecular burden
    S : senescent-cell fraction
    E : epigenetic information integrity
    R : regenerative capacity

The parameter values below are centralized in PAPER1_PARAMS so that the
script is reproducible and easily auditable.
"""

from __future__ import annotations

import os
import copy
import math
import json
from dataclasses import dataclass
from typing import Dict, Callable, Tuple, List

import numpy as np
import pandas as pd
import sympy as sp
from scipy.integrate import solve_ivp
from scipy.optimize import root
from scipy.linalg import eigvals
import matplotlib.pyplot as plt


# ---------------------------------------------------------------------
# Paper 1 baseline parameters
# ---------------------------------------------------------------------

PAPER1_PARAMS: Dict[str, float] = {
    # Deleteriome / damage accrual and repair
    "alpha": 0.018,       # basal deleteriome production
    "beta_S": 0.075,      # senescence-driven deleteriome amplification
    "chi_E": 0.030,       # epigenetic-disorder contribution to D
    "rho_R": 0.090,       # R-dependent repair/removal of D
    "rho_E": 0.055,       # E-dependent repair fidelity
    "rho_rapa": 0.040,    # rapamycin/autophagy effect on D removal

    # Senescence induction and clearance
    "sigma_D": 0.120,     # D-induced senescence
    "sigma_E": 0.055,     # epigenetic-disorder-induced senescence
    "delta_0": 0.018,     # basal senescent-cell clearance
    "delta_R": 0.085,     # immune/regenerative clearance proportional to R
    "delta_seno": 0.400,  # senolytic-induced senescent-cell killing

    # Epigenetic information dynamics
    "eta_D": 0.050,       # D-driven loss of E
    "eta_S": 0.035,       # SASP/senescence-driven loss of E
    "lambda_R": 0.045,    # R-dependent epigenetic maintenance
    "lambda_reprog": 0.300, # partial reprogramming-induced restoration
    "lambda_rapa": 0.030, # rapamycin-induced E stabilization

    # Regenerative capacity
    "K_R": 1.000,         # normalized carrying capacity for R
    "gamma_E": 0.050,     # E-supported regenerative expansion
    "mu_D": 0.035,        # D toxicity to R
    "mu_S": 0.060,        # senescence/SASP toxicity to R
    "u_regen_max": 0.050, # optional regenerative treatment input

    # Numerical and stochastic settings
    "noise_D": 0.035,
    "noise_S": 0.025,
    "noise_E": 0.020,
    "noise_R": 0.020,
    "state_floor": 0.0,
    "state_ceiling": 2.5,
}


INITIAL_YOUTH = np.array([0.10, 0.03, 0.95, 0.90], dtype=float)
INITIAL_AGED_MOUSE_LIVER = np.array([0.80, 0.22, 0.62, 0.42], dtype=float)


@dataclass
class SpeciesParams:
    """Container for species-level scaling of DSER parameters."""
    name: str
    lifespan_years: float
    alpha_scale: float
    repair_scale: float
    senescence_scale: float
    regen_scale: float


SPECIES_LIBRARY: Dict[str, SpeciesParams] = {
    "human": SpeciesParams("human", 80.0, 1.00, 1.00, 1.00, 1.00),
    "mouse": SpeciesParams("mouse", 3.5, 16.9, 0.55, 1.35, 0.80),
    "naked_mole_rat": SpeciesParams("naked_mole_rat", 30.0, 2.0, 1.45, 0.70, 1.10),
    "bowhead_whale": SpeciesParams("bowhead_whale", 200.0, 0.45, 1.55, 0.55, 1.05),
    "bat": SpeciesParams("bat", 35.0, 1.2, 1.35, 0.75, 1.10),
}


# ---------------------------------------------------------------------
# Intervention schedules
# ---------------------------------------------------------------------

def zero_intervention(t: float) -> Dict[str, float]:
    """No treatment."""
    return {"u_seno": 0.0, "u_reprog": 0.0, "u_rapa": 0.0, "u_regen": 0.0}


def liver_protocol_A(t: float) -> Dict[str, float]:
    """
    Protocol A: senolytic first, then rapamycin/autophagy support.
    Mouse liver, 12 weeks:
        weeks 0-4  : senolytic
        weeks 4-12 : rapamycin
    """
    return {
        "u_seno": 1.0 if 0.0 <= t < 4.0 else 0.0,
        "u_reprog": 0.0,
        "u_rapa": 1.0 if 4.0 <= t <= 12.0 else 0.0,
        "u_regen": 0.0,
    }


def liver_protocol_B(t: float) -> Dict[str, float]:
    """
    Protocol B: rapamycin first, then senolytic.
    Mouse liver, 12 weeks:
        weeks 0-4  : rapamycin
        weeks 4-12 : senolytic
    """
    return {
        "u_seno": 1.0 if 4.0 <= t <= 12.0 else 0.0,
        "u_reprog": 0.0,
        "u_rapa": 1.0 if 0.0 <= t < 4.0 else 0.0,
        "u_regen": 0.0,
    }


def liver_protocol_concurrent(t: float) -> Dict[str, float]:
    """Protocol C: concurrent low-intensity senolytic and rapamycin."""
    return {
        "u_seno": 0.5 if 0.0 <= t <= 12.0 else 0.0,
        "u_reprog": 0.0,
        "u_rapa": 0.5 if 0.0 <= t <= 12.0 else 0.0,
        "u_regen": 0.0,
    }


# ---------------------------------------------------------------------
# DSER model
# ---------------------------------------------------------------------

class DSERModel:
    """
    Four-variable dynamical model of aging.

    Variables:
        y[0] = D, accumulated deleteriome burden
        y[1] = S, senescent-cell fraction
        y[2] = E, epigenetic information integrity
        y[3] = R, regenerative capacity

    The ODE is phenomenological, not a mechanistic molecular model.
    Each parameter is intended to be fit or perturbed experimentally.
    """

    def __init__(self, params: Dict[str, float] | None = None):
        self.p = copy.deepcopy(PAPER1_PARAMS if params is None else params)

    def rhs(
        self,
        t: float,
        y: np.ndarray,
        intervention: Callable[[float], Dict[str, float]] = zero_intervention,
    ) -> np.ndarray:
        """Right-hand side of the DSER ODE system."""
        p = self.p
        D, S, E, R = np.clip(y, p["state_floor"], p["state_ceiling"])
        u = intervention(t)

        u_seno = u.get("u_seno", 0.0)
        u_reprog = u.get("u_reprog", 0.0)
        u_rapa = u.get("u_rapa", 0.0)
        u_regen = u.get("u_regen", 0.0)

        dD = (
            p["alpha"]
            + p["beta_S"] * S
            + p["chi_E"] * (1.0 - E)
            - (p["rho_R"] * R + p["rho_E"] * E + p["rho_rapa"] * u_rapa) * D
        )

        dS = (
            (p["sigma_D"] * D + p["sigma_E"] * (1.0 - E)) * (1.0 - S)
            - (p["delta_0"] + p["delta_R"] * R + p["delta_seno"] * u_seno) * S
        )

        dE = (
            -p["eta_D"] * D * E
            -p["eta_S"] * S * E
            + p["lambda_R"] * R * (1.0 - E)
            + p["lambda_reprog"] * u_reprog * (1.0 - E)
            + p["lambda_rapa"] * u_rapa * (1.0 - E)
        )

        dR = (
            p["gamma_E"] * E * R * (1.0 - R / p["K_R"])
            - p["mu_D"] * D * R
            - p["mu_S"] * S * R
            + p["u_regen_max"] * u_regen * (1.0 - R / p["K_R"])
        )

        return np.array([dD, dS, dE, dR], dtype=float)

    def symbolic_jacobian(self):
        """Return symbolic RHS and Jacobian matrix using sympy."""
        D, S, E, R = sp.symbols("D S E R")
        alpha, beta_S, chi_E, rho_R, rho_E, rho_rapa = sp.symbols(
            "alpha beta_S chi_E rho_R rho_E rho_rapa"
        )
        sigma_D, sigma_E, delta_0, delta_R, delta_seno = sp.symbols(
            "sigma_D sigma_E delta_0 delta_R delta_seno"
        )
        eta_D, eta_S, lambda_R, lambda_reprog, lambda_rapa = sp.symbols(
            "eta_D eta_S lambda_R lambda_reprog lambda_rapa"
        )
        K_R, gamma_E, mu_D, mu_S, u_regen_max = sp.symbols(
            "K_R gamma_E mu_D mu_S u_regen_max"
        )
        u_seno, u_reprog, u_rapa, u_regen = sp.symbols(
            "u_seno u_reprog u_rapa u_regen"
        )

        fD = alpha + beta_S*S + chi_E*(1 - E) - (rho_R*R + rho_E*E + rho_rapa*u_rapa)*D
        fS = (sigma_D*D + sigma_E*(1 - E))*(1 - S) - (delta_0 + delta_R*R + delta_seno*u_seno)*S
        fE = -eta_D*D*E - eta_S*S*E + lambda_R*R*(1 - E) + lambda_reprog*u_reprog*(1 - E) + lambda_rapa*u_rapa*(1 - E)
        fR = gamma_E*E*R*(1 - R/K_R) - mu_D*D*R - mu_S*S*R + u_regen_max*u_regen*(1 - R/K_R)

        F = sp.Matrix([fD, fS, fE, fR])
        X = sp.Matrix([D, S, E, R])
        J = F.jacobian(X)
        return F, J

    def jacobian_numeric(self, y: np.ndarray, intervention_state: Dict[str, float] | None = None) -> np.ndarray:
        """Evaluate the analytic Jacobian numerically at a given state."""
        p = self.p
        D, S, E, R = y
        u = {"u_seno": 0.0, "u_reprog": 0.0, "u_rapa": 0.0, "u_regen": 0.0}
        if intervention_state:
            u.update(intervention_state)

        J = np.zeros((4, 4), dtype=float)

        J[0, 0] = -(p["rho_R"] * R + p["rho_E"] * E + p["rho_rapa"] * u["u_rapa"])
        J[0, 1] = p["beta_S"]
        J[0, 2] = -p["chi_E"] - p["rho_E"] * D
        J[0, 3] = -p["rho_R"] * D

        A = p["sigma_D"] * D + p["sigma_E"] * (1.0 - E)
        B = p["delta_0"] + p["delta_R"] * R + p["delta_seno"] * u["u_seno"]
        J[1, 0] = p["sigma_D"] * (1.0 - S)
        J[1, 1] = -A - B
        J[1, 2] = -p["sigma_E"] * (1.0 - S)
        J[1, 3] = -p["delta_R"] * S

        J[2, 0] = -p["eta_D"] * E
        J[2, 1] = -p["eta_S"] * E
        J[2, 2] = -p["eta_D"] * D - p["eta_S"] * S - p["lambda_R"] * R - p["lambda_reprog"] * u["u_reprog"] - p["lambda_rapa"] * u["u_rapa"]
        J[2, 3] = p["lambda_R"] * (1.0 - E)

        J[3, 0] = -p["mu_D"] * R
        J[3, 1] = -p["mu_S"] * R
        J[3, 2] = p["gamma_E"] * R * (1.0 - R / p["K_R"])
        J[3, 3] = (
            p["gamma_E"] * E * (1.0 - 2.0 * R / p["K_R"])
            - p["mu_D"] * D
            - p["mu_S"] * S
            - p["u_regen_max"] * u["u_regen"] / p["K_R"]
        )
        return J

    def simulate(
        self,
        y0: np.ndarray,
        t_span: Tuple[float, float],
        intervention: Callable[[float], Dict[str, float]] = zero_intervention,
        n_eval: int = 500,
    ):
        """Integrate the ODE system."""
        t_eval = np.linspace(t_span[0], t_span[1], n_eval)

        def f(t, y):
            return self.rhs(t, y, intervention)

        sol = solve_ivp(
            f,
            t_span,
            y0,
            t_eval=t_eval,
            method="LSODA",
            rtol=1e-8,
            atol=1e-10,
        )
        if not sol.success:
            raise RuntimeError(sol.message)
        return sol

    def fixed_point(self, guess: np.ndarray) -> np.ndarray:
        """Find a fixed point near the supplied initial guess."""
        def f(y):
            return self.rhs(0.0, y, zero_intervention)

        out = root(f, guess, method="hybr")
        if not out.success:
            raise RuntimeError(f"Fixed-point search failed: {out.message}")
        return np.clip(out.x, self.p["state_floor"], self.p["state_ceiling"])

    def eigen_analysis(self, y_star: np.ndarray) -> np.ndarray:
        """Return eigenvalues of the Jacobian at y_star."""
        J = self.jacobian_numeric(y_star)
        return eigvals(J)

    def r_critical(self, D: float, E: float) -> float:
        """
        Bifurcation threshold for senescence expansion.

        Linearizing dS/dt near small S gives a net coefficient:
            g_S = sigma_D*D + sigma_E*(1-E) - delta_0 - delta_R*R.
        Setting g_S = 0 gives:
            R_critical = [sigma_D*D + sigma_E*(1-E) - delta_0] / delta_R.
        """
        p = self.p
        return (p["sigma_D"] * D + p["sigma_E"] * (1.0 - E) - p["delta_0"]) / p["delta_R"]

    def phase_portrait_DS(
        self,
        E_fixed: float = 0.70,
        R_fixed: float = 0.50,
        outfile: str = "phase_portrait_D_vs_S.png",
    ) -> None:
        """Generate a D-vs-S phase portrait with E and R held fixed."""
        D_grid = np.linspace(0.0, 1.5, 25)
        S_grid = np.linspace(0.0, 0.8, 25)
        DD, SS = np.meshgrid(D_grid, S_grid)
        U = np.zeros_like(DD)
        V = np.zeros_like(SS)

        for i in range(DD.shape[0]):
            for j in range(DD.shape[1]):
                y = np.array([DD[i, j], SS[i, j], E_fixed, R_fixed])
                dy = self.rhs(0.0, y, zero_intervention)
                U[i, j] = dy[0]
                V[i, j] = dy[1]

        speed = np.sqrt(U**2 + V**2) + 1e-12
        plt.figure(figsize=(7, 5))
        plt.streamplot(DD, SS, U / speed, V / speed, density=1.2, color=speed, cmap="viridis")
        plt.xlabel("D: deleteriome burden")
        plt.ylabel("S: senescent-cell fraction")
        plt.title(f"D-S phase portrait, E={E_fixed:.2f}, R={R_fixed:.2f}")
        plt.colorbar(label="relative vector magnitude")
        plt.tight_layout()
        plt.savefig(outfile, dpi=200)
        plt.close()

    def composite_V(self, y: np.ndarray) -> float:
        """
        Composite aging score V.

        Higher D and S increase V; lower E and R increase V.
        Weights are illustrative and should be replaced by fitted z-score weights.
        """
        D, S, E, R = y
        return 1.25 * D + 1.50 * S + 1.00 * (1.0 - E) + 1.00 * (1.0 - R)

    def sensitivity_analysis(
        self,
        y0: np.ndarray,
        t_span: Tuple[float, float] = (0.0, 12.0),
        perturbation: float = 0.20,
    ) -> pd.DataFrame:
        """Perturb each parameter by +/-20% and record final V."""
        baseline = self.simulate(y0, t_span)
        V0 = self.composite_V(baseline.y[:, -1])
        rows = []

        for key in self.p.keys():
            if key.startswith("noise") or key in ["state_floor", "state_ceiling"]:
                continue
            for direction, factor in [("minus", 1.0 - perturbation), ("plus", 1.0 + perturbation)]:
                p2 = copy.deepcopy(self.p)
                p2[key] *= factor
                m2 = DSERModel(p2)
                sol = m2.simulate(y0, t_span)
                V = m2.composite_V(sol.y[:, -1])
                rows.append({
                    "parameter": key,
                    "direction": direction,
                    "factor": factor,
                    "V_final": V,
                    "delta_V": V - V0,
                    "elasticity": ((V - V0) / V0) / (factor - 1.0),
                })

        return pd.DataFrame(rows).sort_values("elasticity", key=np.abs, ascending=False)

    def stochastic_simulation(
        self,
        y0: np.ndarray,
        weeks: float = 12.0,
        dt: float = 0.02,
        n_runs: int = 100,
        intervention: Callable[[float], Dict[str, float]] = zero_intervention,
        seed: int = 1,
    ) -> pd.DataFrame:
        """
        Euler-Maruyama SDE simulation.

        dX = f(X,t)dt + G(X)dW,
        where noise is state-dependent and scaled by sqrt(abs(X)+epsilon).
        """
        rng = np.random.default_rng(seed)
        p = self.p
        n_steps = int(math.ceil(weeks / dt))
        rows = []

        sigma = np.array([p["noise_D"], p["noise_S"], p["noise_E"], p["noise_R"]], dtype=float)

        for run in range(n_runs):
            y = y0.astype(float).copy()
            for k in range(n_steps):
                t = k * dt
                drift = self.rhs(t, y, intervention)
                state_scale = np.sqrt(np.maximum(np.abs(y), 1e-6))
                dW = rng.normal(0.0, math.sqrt(dt), size=4)
                y = y + drift * dt + sigma * state_scale * dW
                y = np.clip(y, p["state_floor"], p["state_ceiling"])

            rows.append({
                "run": run,
                "D": y[0],
                "S": y[1],
                "E": y[2],
                "R": y[3],
                "V": self.composite_V(y),
            })

        return pd.DataFrame(rows)


# ---------------------------------------------------------------------
# Cross-species comparison
# ---------------------------------------------------------------------

def species_scaled_params(base: Dict[str, float], spc: SpeciesParams) -> Dict[str, float]:
    """Scale baseline parameters to approximate species-specific aging dynamics."""
    p = copy.deepcopy(base)
    p["alpha"] *= spc.alpha_scale
    p["rho_R"] *= spc.repair_scale
    p["rho_E"] *= spc.repair_scale
    p["sigma_D"] *= spc.senescence_scale
    p["sigma_E"] *= spc.senescence_scale
    p["gamma_E"] *= spc.regen_scale
    return p


def cross_species_comparison(outfile: str = "cross_species_results.csv") -> pd.DataFrame:
    """Run baseline aging trajectories for human, mouse, NMR, bowhead, and bat."""
    rows = []
    for name, spc in SPECIES_LIBRARY.items():
        p = species_scaled_params(PAPER1_PARAMS, spc)
        model = DSERModel(p)
        y0 = INITIAL_YOUTH.copy()
        horizon = min(spc.lifespan_years, 120.0)
        sol = model.simulate(y0, (0.0, horizon), n_eval=600)
        yf = sol.y[:, -1]
        rows.append({
            "species": name,
            "lifespan_years": spc.lifespan_years,
            "alpha_effective": p["alpha"],
            "repair_scale": spc.repair_scale,
            "D_final": yf[0],
            "S_final": yf[1],
            "E_final": yf[2],
            "R_final": yf[3],
            "V_final": model.composite_V(yf),
        })
    df = pd.DataFrame(rows)
    df.to_csv(outfile, index=False)
    return df


# ---------------------------------------------------------------------
# Design-loop simulation for Section 7.3
# ---------------------------------------------------------------------

def run_design_loop() -> pd.DataFrame:
    """Compare Protocol A, Protocol B, concurrent treatment, and no treatment."""
    model = DSERModel(PAPER1_PARAMS)
    protocols = {
        "no_intervention": zero_intervention,
        "protocol_A_senolytic_then_rapamycin": liver_protocol_A,
        "protocol_B_rapamycin_then_senolytic": liver_protocol_B,
        "protocol_C_concurrent": liver_protocol_concurrent,
    }

    rows = []
    for name, protocol in protocols.items():
        sol = model.simulate(INITIAL_AGED_MOUSE_LIVER, (0.0, 12.0), protocol, n_eval=500)
        y_final = sol.y[:, -1]
        rows.append({
            "protocol": name,
            "D_final": y_final[0],
            "S_final": y_final[1],
            "E_final": y_final[2],
            "R_final": y_final[3],
            "V_final": model.composite_V(y_final),
            "delta_D": y_final[0] - INITIAL_AGED_MOUSE_LIVER[0],
            "delta_S": y_final[1] - INITIAL_AGED_MOUSE_LIVER[1],
            "delta_E": y_final[2] - INITIAL_AGED_MOUSE_LIVER[2],
            "delta_R": y_final[3] - INITIAL_AGED_MOUSE_LIVER[3],
        })

    return pd.DataFrame(rows).sort_values("V_final")


# ---------------------------------------------------------------------
# Main script
# ---------------------------------------------------------------------

def main() -> None:
    """Execute all analyses and write outputs."""
    os.makedirs("outputs", exist_ok=True)

    model = DSERModel(PAPER1_PARAMS)

    # 1. Symbolic Jacobian
    F_sym, J_sym = model.symbolic_jacobian()
    with open("outputs/symbolic_rhs_and_jacobian.txt", "w") as f:
        f.write("RHS F(D,S,E,R):\n")
        f.write(str(F_sym))
        f.write("\n\nJacobian J:\n")
        f.write(str(J_sym))

    # 2. Fixed points and eigenvalue analysis
    youth_fp = model.fixed_point(INITIAL_YOUTH)
    aged_guess = model.simulate(INITIAL_YOUTH, (0.0, 80.0), n_eval=800).y[:, -1]
    aged_fp = model.fixed_point(aged_guess)

    youth_eigs = model.eigen_analysis(youth_fp)
    aged_eigs = model.eigen_analysis(aged_fp)

    eig_df = pd.DataFrame({
        "fixed_point": ["youth"] * 4 + ["aged"] * 4,
        "eigenvalue_real": list(np.real(youth_eigs)) + list(np.real(aged_eigs)),
        "eigenvalue_imag": list(np.imag(youth_eigs)) + list(np.imag(aged_eigs)),
    })
    eig_df.to_csv("outputs/eigenvalues.csv", index=False)

    # 3. Bifurcation threshold R_critical
    Rcrit_youth = model.r_critical(D=youth_fp[0], E=youth_fp[2])
    Rcrit_aged = model.r_critical(D=aged_fp[0], E=aged_fp[2])
    with open("outputs/R_critical.json", "w") as f:
        json.dump({
            "Rcrit_youth": float(Rcrit_youth),
            "Rcrit_aged": float(Rcrit_aged),
            "youth_fixed_point": youth_fp.tolist(),
            "aged_fixed_point": aged_fp.tolist(),
        }, f, indent=2)

    # 4. Phase portrait
    model.phase_portrait_DS(
        E_fixed=float(aged_fp[2]),
        R_fixed=float(aged_fp[3]),
        outfile="outputs/phase_portrait_D_vs_S.png",
    )

    # 5. Section 7.3 design-loop simulation
    design_df = run_design_loop()
    design_df.to_csv("outputs/design_loop_protocols.csv", index=False)

    # 6. Parameter sensitivity
    sens_df = model.sensitivity_analysis(INITIAL_AGED_MOUSE_LIVER)
    sens_df.to_csv("outputs/parameter_sensitivity_pm20.csv", index=False)

    # 7. Cross-species comparison
    species_df = cross_species_comparison("outputs/cross_species_results.csv")

    # 8. Stochastic SDE simulation, 100 runs per protocol
    stochastic_rows = []
    stochastic_protocols = {
        "A": liver_protocol_A,
        "B": liver_protocol_B,
        "C": liver_protocol_concurrent,
        "none": zero_intervention,
    }
    for idx, (name, protocol) in enumerate(stochastic_protocols.items()):
        sdf = model.stochastic_simulation(
            INITIAL_AGED_MOUSE_LIVER,
            weeks=12.0,
            dt=0.02,
            n_runs=100,
            intervention=protocol,
            seed=100 + idx,
        )
        sdf["protocol"] = name
        stochastic_rows.append(sdf)
    stochastic_df = pd.concat(stochastic_rows, ignore_index=True)
    stochastic_df.to_csv("outputs/stochastic_100_runs.csv", index=False)

    summary = stochastic_df.groupby("protocol").agg(
        D_mean=("D", "mean"),
        D_sd=("D", "std"),
        S_mean=("S", "mean"),
        S_sd=("S", "std"),
        E_mean=("E", "mean"),
        E_sd=("E", "std"),
        R_mean=("R", "mean"),
        R_sd=("R", "std"),
        V_mean=("V", "mean"),
        V_sd=("V", "std"),
    )
    summary.to_csv("outputs/stochastic_summary.csv")

    print("Paper 1 DSER analyses complete. Outputs written to ./outputs/")
    print("\nDesign-loop ranking:")
    print(design_df[["protocol", "V_final", "delta_S", "delta_E"]])
    print("\nCross-species comparison:")
    print(species_df[["species", "lifespan_years", "alpha_effective", "V_final"]])


if __name__ == "__main__":
    main()

References

1. Weismann A. Über die Dauer des Lebens. Jena: Gustav Fischer; 1882.
2. Medawar PB. An Unsolved Problem of Biology. London: H.K. Lewis; 1952.
3. Harman D. Aging: a theory based on free radical and radiation chemistry. J Gerontol. 1956;11(3):298-300.
4. Williams GC. Pleiotropy, natural selection, and the evolution of senescence. Evolution. 1957;11(4):398-411.
5. Hamilton WD. The moulding of senescence by natural selection. J Theor Biol. 1966;12(1):12-45.
6. Kirkwood TBL. Evolution of ageing. Nature. 1977;270:301-304.
7. Harley CB, Futcher AB, Greider CW. Telomeres shorten during ageing of human fibroblasts. Nature. 1990;345:458-460.
8. Ames BN, Shigenaga MK, Hagen TM. Oxidants, antioxidants, and the degenerative diseases of aging. PNAS. 1993;90(17):7915-7922.
9. Gavrilov LA, Gavrilova NS. The reliability theory of aging and longevity. J Theor Biol. 2001;213(4):527-545.
10. Mitnitski AB, Mogilner AJ, Rockwood K. Accumulation of deficits as a proxy measure of aging. ScientificWorldJournal. 2001;1:323-336.
11. de Grey ADNJ, et al. Time to talk SENS. Ann N Y Acad Sci. 2002;959:452-462.
12. Speakman JR. Body size, energy metabolism and lifespan. J Exp Biol. 2005;208:1717-1730.
13. Andziak B, et al. High oxidative damage levels in the longest-living rodent, the naked mole-rat. Aging Cell. 2006;5(6):463-471.
14. Blagosklonny MV. Aging and immortality: quasi-programmed senescence. Cell Cycle. 2006;5(18):2087-2102.
15. Harrison DE, et al. Rapamycin fed late in life extends lifespan. Nature. 2009;460:392-395.
16. Baker DJ, et al. Clearance of p16-positive senescent cells. Nature. 2011;479:232-236.
17. Horvath S. DNA methylation age. Genome Biol. 2013;14:R115.
18. López-Otín C, et al. The hallmarks of aging. Cell. 2013;153(6):1194-1217.
19. Kennedy BK, et al. Geroscience. Cell. 2014;159(4):709-713.
20. Keane M, et al. Insights from the bowhead whale genome. Cell Rep. 2015;10(1):112-122.
21. Ocampo A, et al. In vivo amelioration by partial reprogramming. Cell. 2016;167(7):1719-1733.
22. Zhavoronkov A, et al. AI for aging and longevity. Ageing Res Rev. 2019;49:49-66.
23. Zhavoronkov A, et al. Deep learning enables DDR1 kinase inhibitors. Nat Biotechnol. 2019;37(9):1038-1040.
24. Sinclair DA. Lifespan. Atria Books; 2019.
25. Lu Y, et al. Reprogramming to recover youthful epigenetic information. Nature. 2020;588:124-129.
26. Gladyshev VN. The ground zero of organismal life and aging. Trends Mol Med. 2021;27:11-19.
27. Pyrkov TV, et al. Progressive loss of resilience predicts human lifespan limit. Nat Commun. 2021;12:2765.
28. Pun FW, et al. Hallmarks-based dual-purpose targets via PandaOmics. Aging. 2022;14(6):2475-2506.
29. Pun FW, et al. Dual-purpose targets for cancer and aging. Aging Cell. 2023;22(12):e14017.
30. López-Otín C, et al. Hallmarks of aging: An expanding universe. Cell. 2023;186(2):243-278.
31. Ivanenkov YA, et al. Chemistry42. J Chem Inf Model. 2023;63(3):695-701.
32. Aladinskiy V, et al. TNIK inhibitors for IPF. J Med Chem. 2024;67(21):19021-19038.
33. Kamya P, et al. PandaOmics. J Chem Inf Model. 2024;64(10):3961-3969.
34. Xu Z, et al. AI-discovered TNIK inhibitor: phase 2a trial. Nat Med. 2025;31(8):2602-2610.
35. Gao J, Barzel B, Barabási AL. Universal resilience patterns in complex networks. Nature. 2016;530:307-312.
36. Yang JH, et al. Loss of epigenetic information as a cause of mammalian aging. Cell. 2023;186(2):305-326.
37. Takahashi K, Yamanaka S. Induction of pluripotent stem cells. Cell. 2006;126(4):663-676.
38. Franceschi C, et al. Inflamm-aging. Ann N Y Acad Sci. 2000;908:244-254.
39. Baker DJ, et al. Naturally occurring p16-positive cells shorten healthy lifespan. Nature. 2016;530:184-189.
40. Goodell MA, Rando TA. Stem cells and healthy aging. Science. 2015;350(6265):1199-1204.
41. Levine ME, et al. An epigenetic biomarker of aging. Aging. 2018;10(4):573-591.
42. Hannum G, et al. Genome-wide methylation profiles. Mol Cell. 2013;49(2):359-367.