Teleological Calculus for Viral Evolution: Categorical Structure, Multi-Level Invariants, and Predictive Jump-Risk Indices

Teleonomical Calculus for Viral Evolution: Categorical Structure, Multi-Level Invariants, and Predictive Jump-Risk Indices Andrei T. Patrascu 1 1 FAST Foundation, Destin FL, 32541, USA email: andrei.patr[email protected] Evolutionary change in rapidly adapting pathogens does not proceed as an unconstrained hill-climb: empirical trajectories organize into narrow channels and, at times, display abrupt jumps. We present a categorical, teleonomical calculus that predicts both behaviors. Final causes are encoded as coherence deficits across six levels (L0–L5): basic viability (L0); trait-domain and tradeoff constraints (L1); structural invariants from topology and sheaf theory (L2); routine life-cycle organization (L3); path-order effects (holonomy, L4); and cross-level coherence (L5). Evolutionary updates minimize total deficit over reachable futures, enforced by a mutation-graph corridor, rather than over abstract states. Applied to SARS-CoV-2 Omicron lineages (BA.1 → BA.2 → BA.5 → XBB), the method combines: (i) a soft-tube channel around the empirical manifold in a four-trait space (binding, cleavage, stability, escape) with a global progress coordinate; (ii) per-variant structural terms (L2) derived from contact-graph topology and sheaf gluing; and (iii) a graph-constrained search that respects mutational reachability. The predicted trajectory closely tracks observed multi-trait evolution in all pairwise projections, as quantified by path-distance metrics (4D dynamic time warping and Hausdorff distance). To anticipate discontinuities, we introduce jump-risk indices grounded in the calculus: rapid barrier drops in the teleonomical objective, local curvature softening, L2 obstruction flips, L4 holonomy spikes, and L3 routine re-routing. These indices highlight segments where a jump is plausible, complementing the smooth channel forecast. The framework unifies Aristotle’s four causes with categorical structure into a practical predictive tool for viral surveillance and, more broadly, for systems whose dynamics are guided by invariants under reachability. I. INTRODUCTION A. Motivation Empirical evolution is neither a random walk nor an unconstrained hill–climb on a scalar landscape. Across viruses and other fast–adapting systems, observed trajectories concentrate into corridors carved jointly by biophysics (e.g., receptor binding and proteolytic activation), population immunity, epidemiological context, and historical order effects. Classic theory provides indispensable pieces—selection and drift, quasispecies, rugged landscapes—but it leaves two gaps that matter in practice: (i) representing teleonomy, that is, goal–directed maintenance of invariants the organism “strives” to keep true, and (ii) reasoning across levels (molecular structure, life–cycle routine, population context) while enforcing reachability. The present work closes these gaps by casting evolution as the minimization of a coherence deficit across categorical levels, over the set of futures actually reachable from the current state. We show that this yields a computable calculus that predicts the smooth channels followed by SARS–CoV–2 Omicron lineages (BA.1→BA.2→BA.5→XBB) and supplies principled indices for discontinuous jumps. Biological motivation. SARS-CoV-2 spike adaptation is governed by coordinated changes in (i) ACE2 binding, (ii) S1/S2 cleavage, (iii) spike stability, and (iv) epitope-resolved immune escape. Our goal in this paper is operational: predict which trait-space corridors are reachable by real mutations and when a jump becomes likely. We therefore frame a teleonomical calculus whose energy collects multi-level biological invariants and whose updates are constrained to mutationally reachable futures. We show this reproduces the BA.1 → BA.2 → BA.5 → XBB trajectory across trait projections and yields falsifiable jump warnings (Section V and App. J). Historical context. From Darwin’s qualitative synthesis [ 1 ] through Fisher’s quantitative genetics [ 161 ] and Wright’s landscape metaphor [ 160 ], evolutionary theory has emphasized selection, mutation, and drift. Kimura’s neutral theory [ 4 ] and Eigen’s quasispecies model [ 5 ] sharpened the role of mutational neighborhoods in high–mutation regimes. Rugged landscapes (e.g., NK models) emphasized epistasis and multiple basins [ 6 ]. In pathogens, antigenic cartography [ 181 ] and phylodynamics [ 8 ] tied sequence change to population–level spread, while deep mutational scanning (DMS) quantified local genotype–phenotype maps [ 9 , 174 ]. Alongside these, topology and sheaf theory have emerged as tools to encode global structure and local consistency in biological networks [ 12 , 13 , 196 ]. Our contribution integrates these strands by making the final cause—maintenance of cross–level invariants—explicit in a categorical formalism with 2 an optimization–dynamical realization. Idea in brief. We model the system’s state as z= (x, s)∈ Z := Rd |{z} traits × S |{z} structures , where x collects experimental traits (here: binding, cleavage, stability, escape) and s collects structural/life– cycle objects (contact graphs, sheaf sections, routine kernels). We introduce six levels L 0– L 5indexed by i∈ { 0 ,..., 5 } , each with a categorical object Ci and a functorial invariant Ii : Ci→R≥0 whose value vanishes when the level’s constraints are satisfied. Composing each Ii with the appropriate forgetful or assembly functor from state space to Ciyields a deficit Li(z;λ) := Ii◦Assemblei(z;λ)∈R≥0,(1) where λdenotes environment (e.g., immune context). The total coherence deficit is L(z;λ) = 5 X i=0 wiLi(z;λ),(2) with nonnegative weights wi(learned with ranking–regularization, cf. §I B). Teleonomy enters via the reachability relation. Let Reach∆t ( z )be the set of states attainable from z within time ∆ t by a bounded number of mutations/recombinations along a mutation graph G and by permitted deformations of routine/structure (life–cycle kernel moves that keep the routine closed). A single teleonomical step is the proximal problem zk+1 ∈arg min y∈Reach∆t(zk)nDψ(ykzk) + ε∆tL(y;λk)o,(3) where Dψ is a Bregman distance that plays the role of a short–step cost and ε scales selection strength. Passing to the small–step limit yields a differential inclusion ˙z(t)∈ −∂L(·;λ(t)) + ιReach(·)z(t),(4) where ιReach is the indicator of reachability and ∂ is a subdifferential in the Dψ –geometry. Equations (3) – (4) formalize the slogan: minimize coherence deficit over futures the system can actually reach. Levels used in this study. For concreteness we instantiate: •L0 (viability). A baseline sufficiency term (e.g., minimal epidemiological fitness). •L1 (trait domain & tradeoff). Band penalties and anisotropic targets by era, yielding interior gradients toward biologically plausible regions. •L2 (structure). A sum of (i) persistent homology shortfall of residue–contact graphs and (ii) sheaf–gluing residuals over a functional patch cover; both are functorial in the underlying structures [12, 13, 196]. •L3 (routine). Coherence of the life–cycle kernel (attach → prime → fuse → replicate → egress), measured via spectral gap and conductance. •L4 (holonomy). Path–order sensitivity (immune → drug versus drug → immune), represented as a categorical holonomy and penalized by a path–discrepancy norm (placeholder here; protocol outlined below). •L5 (multi–level coherence). A cross–level variance/sparsity term discouraging fixes at one level that violate another. Each level in the hierarchy encapsulates a distinct biological domain—ranging from genomic viability up to environmental and cross-level coordination—and contributes an additive term to the overall coherence deficit L = P5 i=0 Li . Table I summarizes this correspondence between mathematical level, biological meaning, and the type of constraint it enforces. This mapping provides a concrete biological interpretation for what the calculus minimizes: deviations from feasible genotype–phenotype–environment alignments across scales. 3 Table I: Biological mapping of levels L0–L5. Level Biological object What the deficit penalizes L0Viability baseline Infeasible fitness floor (e.g., nonproductive entry/replication) L1Trait domain/tradeoff Exiting calibrated trait bands; wrong era-wise direction (Sec. IVA) L2Structure Mismatch to contact-topology PH and patch gluing residuals (RBM/RBD/S1S2/glycans) L3Routine Inconsistent life-cycle routine (entry/priming/fusion sequencing) L4Order/holonomy Non-commuting order effects (e.g., cleavage-beforebinding vs binding-before-cleavage) L5Cross-level balance Incoherent cross-level tradeoffs under environment λ Channel prediction and jump detection. Two additional ingredients make (3) predictive. First, a soft tube T around an empirical or learned manifold Γin trait–structure space enforces nearness to historical corridors while keeping gradients inside the tube (no flat regions). We use a soft minimum of Mahalanobis distances to contiguous segments of Γwith an annealed sharpness parameter, and a global progress coordinate S ( z ) ∈ [0 , 3] that advances along BA.1 → BA.2 → BA.5 → XBB. Second, we constrain moves by a mutation–graph corridor G = ( V, E )and solve (3) either by a trust–region proximal gradient (continuous formulation) or by a graph beam search (discrete formulation). Jumps emerge when the minimal value of L across adjacent basins drops sharply (barrier collapse), when the local curvature of L softens (smallest Hessian eigenvalue ↓ 0), when L2 topological obstructions flip, when L4 holonomy spikes, or when L3 conductance reveals a new high–throughput route. B. Why this method matters Scientific relevance. The calculus unifies Aristotle’s causes with modern categorical structure: material cause (sequence/structure), formal cause (sheaves, PH, routines), efficient cause (mutation and selection realized by (3) ), and final cause (coherence of invariants encoded by (1) – (2) ). It delivers a single, computable object L whose geometry explains: (i) why evolution concentrates in channels (low–deficit valleys intersected with reachability), and (ii) when discontinuities are favored (barrier/holonomy/routine criteria). Practical impact. On SARS–CoV–2, the method produces smooth, spike–free predicted trajectories that closely match the observed BA.1 → BA.2 → BA.5 → XBB path across multiple trait projections while respecting mutational reachability. Compared to landscape or unconstrained gradient models, it (a) integrates structural information (L2), (b) enforces realistic move sets (graph corridor), (c) adapts to environment by era, and (d) provides operational jump–risk indices. These features support prospective surveillance: rank plausible futures along the corridor; flag segments where a jump is likely; and generate mechanistic hypotheses (which patches or routine steps are under categorical tension). Notation. We write x = ( Tbind, Tcleave, Tstab, Tescape ) ∈R4 , s∈ S , z = ( x, s ). The environment λ ( t ) is piecewise–constant by era. The mutation graph is G = ( V, E )with nodes v∈V carrying states zv and edges E encoding single–step edits or recombination. Distances are Mahalanobis with covariance Σ when specified; Dψ denotes a Bregman distance. Equations (1) – (4) will be referenced throughout (labels reserved and unique). C. Our contribution We develop a categorical teleonomical calculus for evolution that (a) formalizes final causes as coherence deficits across levels L 0– L 5, (b) induces a flow that minimizes total deficit subject to reachability, and (c) incorporates structural (L2), routine (L3), and holonomy (L4) terms that jointly explain both smooth evolutionary channels and potential jumps. We demonstrate the framework on SARS–CoV–2 Omicron lineages (BA.1 → BA.2 → BA.5 → XBB) with artifact–free, quantitatively evaluated predictions across four trait projections. 4 (C1) Final causes as coherence deficits on categorical levels. Building on (1) – (2) , each level Li is a functorially defined invariant whose violation is measured by a nonnegative functional Li ( z ; λ ). In particular, we explicitly decompose the structural level as L2(z;λ) = LPH 2(s;λ) + Lglue 2(s;λ),(5) where LPH 2 penalizes shortfalls relative to a target persistent–homology profile of the residue–contact graph (using stability–aware distances between diagrams), and Lglue 2 penalizes sheaf–gluing residuals over a fixed functional patch cover (e.g., RBM, RBD, S1/S2, glycans). The remaining levels L 0 , L 1 , L 3 , L 4 , L 5 are instantiated as in §I, with L 1providing interior gradients toward biologically plausible trait domains and L5ensuring cross–level compatibility. (C2) A teleonomical flow under reachability constraints. Given the total deficit L in (2) and the reachability set Reach∆t ( z ), a single update solves the proximal problem (3) with a Bregman distance Dψ ; in the small–step limit we obtain the differential inclusion (4) . Numerically, we realize this by two complementary solvers: 1. Continuous proximal gradient with trust regions: a line–searched descent in the Dψ –geometry that enforces tube proximity and step caps while guaranteeing descent (cf. proximal and Bregman methods [14, 16, 191]). 2. Discrete beam search on a mutation graph: a corridor G = ( V, E )in sequence/trait space with node states zvand edges encoding reachable edits; paths are scored by a state cost Φ(z;λ) = Ltube(z) | {z } channel proximity +α2L2(z;λ) | {z } structural coherence +α0L0(z;λ) | {z } viability ,(6) and selected by best–first beam expansion under reachability. Both views are instances of the same principle: minimize Lover futures the system can actually reach. (C3) Channel geometry and jump diagnostics. To shape and parameterize evolutionary corridors, we introduce a soft tube around a reference manifold Γ(empirical or learned). Writing Γas a polyline or spline with segments {Γi}, we define Ltube(z) = −1 βlogX i exp −β1 2kz−Πi(z)k2 Σ−1,(7) where Π i is the Mahalanobis projection to segment Γ i and β is annealed to harden selection. A global progress coordinate S ( z ) ∈ [0 , m ]aggregates the segment index and local parameter, and a scheduled target Starget provides forward drive along the corridor. Jump candidates are flagged when (i) the minimal deficit between adjacent basins drops sharply (barrier collapse), (ii) the smallest curvature surrogate of L softens toward zero (incipient bifurcation), (iii) L2 exhibits a topological or gluing obstruction flip, (iv) L4(holonomy) spikes, or (v) L3indicates a new high–throughput routine. (C4) Empirical accuracy on Omicron lineages. For the four trait coordinates ( Tbind, Tcleave, Tstab, Tescape ), our predictions produce smooth, spike–free channels that align with the observed BA.1 → BA.2 → BA.5 → XBB path across all pairwise projections, while respecting mutational reachability in the corridor. Quantitative agreement is reported with 4D dynamic time warping and Hausdorff distances [ 17 , 18 ]. Importantly, the categorical levels are load–bearing: ablations that remove L2 or the graph corridor degrade fit and reintroduce artifacts, underscoring the necessity of structural coherence and reachability in teleonomical forecasting. Why this matters. This calculus unifies final–cause reasoning with categorical structure into a computable predictor: a single, geometry–aware objective L explains (i) why evolution concentrates into channels (low–deficit valleys intersected with reachability), and (ii) when discontinuities are favored (barrier/holonomy/routine criteria). Practically, it yields artifact–free, biologically plausible forecasts and operational jump–risk indices, enabling prospective surveillance and mechanistic insight. D. Summary of findings Our categorical teleonomical calculus yields four principal empirical findings on the Omicron case study (BA.1 → BA.2 → BA.5 → XBB), each consistent with the formalism in §I A–I C and the proximal update (3) / flow (4). 5 (F1) Channel prediction aligns with observed data across traits. Using the soft tube Ltube in (7) and the global progress coordinate S ( z ), the predicted trajectory closely follows the observed manifold in the four trait projections ( Tbind, Tcleave, Tstab, Tescape ), without spikes or geometric artifacts. Alignment is quantified by curve–to–curve distances (dynamic time warping and Fréchet families) computed in the joint trait space; these metrics capture both shape and progression timing and therefore penalize biologically implausible backtracking [19, 20]. The channel geometry arises endogenously from minimizing the total deficit L subject to reachability, which concentrates the flow inside low–deficit valleys intersected with the mutation corridor. (F2) Jump–risk indices identify plausible discontinuities. We operationalize five indices derived from the calculus: (i) sharp decreases in the minimal L between adjacent basins (barrier collapse), (ii) softening of local curvature (smallest eigenvalue of a Hessian surrogate approaching zero), (iii) structural flips in L2 (persistent–homology or gluing obstructions), (iv) spikes in L4 (holonomy under order reversal), and (v) routine re–routing signals in L3 (conductance/spectral patterns). Segments where two or more indices exceed thresholds are flagged as jump candidates. In our analyses these flags occur near bends where the empirical path turns, matching domain expectations of context–driven discontinuities. (F3) Ablations reveal which levels are decisive. Systematic ablations—removing one component at a time from L while keeping the solver and corridor fixed—show that (a) eliminating L2 (structural coherence) degrades fit and reintroduces unrealistic detours, (b) removing the mutation–graph constraint (reachability) increases distances and produces lateral moves absent in real data, and (c) dropping the smoothness/trust mechanisms yields kinks or spikes. These patterns affirm that L2 and the graph reachability constraints are load–bearing for accurate, biologically plausible forecasts, while other levels modulate fine structure and jump propensity. The ablation protocol follows standard model assessment practice: measure changes in fit under controlled component removals to attribute contribution reliably [21]. (F4) Robustness and interpretability. Predicted channels remain stable under moderate hyperparameter variation (tube sharpness β , progress schedule for Starget , proximal weight), and the contribution of each level Li is interpretable in terms of concrete biological objects (contact topology, patch consistency, routine transitions, order effects). This clarity enables mechanistic hypotheses linking candidate future changes to specific categorical tensions (e.g., which patch gluing or routine edge is ‘expensive’ in L ), and supports prospective use in surveillance and design. Biological overview. Although the framework developed here is mathematically structured, its motivation is biological: SARS–CoV–2 evolution proceeds under strong constraints imposed by ACE2 binding, S1/S2 cleavage, spike stability, and class-specific immune escape. These four properties constitute the effective trait space explored by circulating variants, and their quantitative interplay determines which mutational paths are viable. The teleonomic calculus developed below provides a compact, mechanistic description of how these trait constraints jointly shape evolutionary trajectories. II. BACKGROUND AND RELATED WORK This section reviews the mathematical and empirical foundations on which our categorical teleonomical calculus builds: (i) fitness landscapes and replicator–mutator dynamics; (ii) antigenic cartography as a metric–embedding problem; (iii) phylodynamics and likelihoods on trees; (iv) deep–mutational–scanning (DMS)–guided genotype–phenotype models; (v) geometric flows and proximal calculus (JKO/Bregman) as optimization–dynamics bridges; and (vi) topological and sheaf–theoretic methods for structural invariants. We close by identifying the conceptual gap our framework fills: representing final causes and category–level invariants in a tractable calculus with reachability constraints. A. Fitness landscapes and replicator–mutator dynamics Landscapes and epistasis. Let G be a finite genotype space (e.g., a Hamming hypercube on L loci), W:G → R>0a Malthusian fitness, and pt∈∆(G)genotype frequencies. Epistasis is the non–additivity of W over loci; quantitative frameworks include Walsh–Fourier expansions and accessibility percolation on high–dimensional landscapes. Modern treatments emphasize multiple scales and ridges rather than isolated peaks [22–26]. 6 Derivation of the replicator–mutator equation. Let M = ( Mij ) i,j∈G be a row–stochastic mutation kernel; in continuous time, the mass–action birth–death–mutation model yields dpi dt =X j∈G pjWjMji −φ(t)pi, φ(t) = X k pkWk,(8) which is the replicator–mutator ODE. Derivation. Over a small interval ∆ t , expected inflow to i equals PjpjWj ∆ t Mji ; outflow equals piφ ( t ) ∆ t when normalizing to keep Pipi = 1. Divide by ∆ t and take ∆t↓0to obtain (8). Landscape limits and constraints. Equation (8) is agnostic to structural invariants and to multi–level routines; it is also posed on ∆( G )rather than a trait–structure state space. While it captures selection– mutation balance, it lacks (i) explicit invariants (final causes) and (ii) a mechanism to encode reachability beyond point mutations (e.g., constrained routine moves). These omissions motivate the categorical augmentation developed in §I. Relation to classical fitness landscapes. The coherence deficit L ( g ; λ )provides a multi-level generalization of the classical Wrightian notion of a fitness landscape. In the traditional formulation, each genotype g is assigned a scalar fitness F ( g ), and evolutionary dynamics approximate gradient ascent on this surface subject to mutational constraints. In our framework, this role is played by the negative coherence deficit: Feff(g;λ)≡ −L(g;λ),(9) up to an arbitrary monotone rescaling. Minimizing L is therefore equivalent to ascending an effective fitness landscape whose geometry is shaped by the hierarchical levels L0 – L5 . Each level contributes a specific kind of “terrain” on this landscape. The viability term L0 enforces a hard boundary between permissible and impermissible regions of sequence space; L1assigns smooth trait-domain penalties that correspond to conventional fitness slopes (e.g. gradual changes in ACE2 affinity, cleavage efficiency, or conformational stability); L2 introduces curvature via structural and biochemical couplings, producing epistatic ridges and valleys; L3 and L4 generate path-dependent contributions that have no analogue in standard Wrightian surfaces, capturing routine-ordering and holonomy effects that encode how the sequence of edits matters, not just their endpoints; and L5 couples molecular and population-level constraints, modulating the landscape by the time-varying environment λ(t). With this identification Feff = −L , the monotone decrease guaranteed by (34) means that reachable dynamics follow fitness-increasing directions in the Wrightian sense, while the mutation-graph corridor restricts motion to locally admissible edits. This is directly analogous to traversing neutral networks or quasi-neutral plateaus on rugged fitness landscapes: the geometry of the corridor—determined by DMS-admissible edits, structural feasibility, and recombination constraints—defines which slope directions are accessible. The higher-level terms L2 – L4 further imply that valley-crossing and epistatic ridge-line navigation arise when the coherence increase due to an intermediate edit is offset by changes in structural curvature or order effects, reproducing classical phenomena such as extra-dimensional bypass, conversion steps, and detour paths. In summary, coherence-deficit minimization provides a principled, multi-trait, multi-level generalization of Wright’s fitness landscape, in which −L plays the role of fitness and the mutation-graph corridor specifies the evolutionary pathways available to the system. B. Antigenic cartography and metric embedding Antigenic cartography embeds strains and sera as points {xu} ⊂ Rm (typically m = 2) so that Euclidean distances approximate assay–derived antigenic distances duv ; extensions integrate molecular evolution and predictability [27–30]. Stress and gradients. Given weights wuv ≥0and targets duv >0, the classical stress is S({x}) = X u<v wuv kxu−xvk2−duv 2.(10) The gradient w.r.t. xuis ∇xuS= 2 X v6=u wuv 1−duv kxu−xvk2+δ(xu−xv),(11) with a small δ > 0for numerical stability. While S is convex in interpoint distances, it is not jointly convex in coordinates; majorization and SMACOF methods are standard. 7 Relevance and limits. Antigenic maps partially capture environmental (immune) drivers and geometric constraints, but: (i) they do not encode category–level structural invariants (e.g., contact–graph topology), (ii) they lack reachability constraints, and (iii) their dynamics are not generically derived from a cross–level objective. Our soft–tube construction (§I C) borrows the metric–embedding flavor but embeds it in a teleonomical, reachability–aware flow. C. Phylodynamics and likelihood on trees Phylodynamics links sequence–dated phylogenies and population–level processes via coalescent or birth–death models with time–varying rates and sampling [ 31 – 34 ]. Let T be a rooted, time–scaled tree with branching times {ti}. Birth–death skyline likelihood (outline). Under birth rate λ ( t ), death/removal rate µ ( t ), and sampling ψ ( t ), the likelihood of T factors through the density of branching times and sampling events conditional on survival; piecewise–constant λ, µ, ψ on intervals admits closed–form recursions for the probability of lineages through time [ 33 ]. Coalescent models give alternative likelihoods via effective population size Ne(t). Relevance and limits. Phylodynamic likelihoods power inference on transmission and growth, but the state is a tree (or its summary), not a trait–structure state with categorical invariants. They lack explicit final–cause terms and do not enforce structural reachability in trait space. D. Deep mutational scanning and genotype–phenotype maps DMS quantifies the effects of mutations at scale, enabling statistical models of genotype–phenotype relationships and epistasis [35, 36]. Two paradigms are especially relevant. Walsh–Fourier epistasis. Index genotypes by g∈ {± 1 }L and phenotype by y ( g ). The Walsh expansion y(g) = X S⊆{1,...,L} ˆySY i∈S gi(12) decomposes additive and higher–order epistasis via coefficients ˆyS (computed by orthogonality). This gives a complete, though high–dimensional, representation of epistasis structure. Global epistasis (monotone link). Assume a latent additive score η ( g ) = Piβixi ( g )( xi one–hot features) and a nonlinear link fsuch that y(g)≈fη(g)+. (13) Nonlinearity f explains much of apparent higher–order epistasis and improves extrapolation [ 37 ]. This yields smooth genotype–trait predictors compatible with reachability constraints. Relevance and limits. DMS–informed predictors estimate the material cause (sequence–to–trait map) but do not by themselves enforce formal (structural), efficient (routine), or final (invariant) causes; nor do they induce a teleonomical flow with reachability. In our calculus, DMS enters as a component of the state map z= (x, s)and of the corridor constraints. E. Geometric flows and proximal calculus Gradient flows in metric spaces. Given an energy E on a metric space (X , d ), the minimizing movement (JKO) scheme constructs a discrete flow by proximal steps xk+1 ∈arg min y∈Xn1 2τd(y, xk)2+E(y)o.(14) As τ↓ 0, xτ ( · )converges (under conditions) to a curve of maximal slope/gradient flow for E [ 38 – 40 ]. This underlies diffusion and transport dynamics. Bregman proximality. In Banach spaces with a strictly convex, differentiable ψ, the Bregman step xk+1 ∈arg min ynDψ(ykxk) + τ E(y)o, Dψ(ykx) = ψ(y)−ψ(x)−h∇ψ(x), y −xi(15) generalizes Euclidean proximal descent and yields mirror–like flows. Our teleonomical step (3) is a constrained Bregman proximal update with E=L(·;λ). 8 Soft minimum and smooth tube. Let di ( z ) = 1 2kz− Π i ( z ) k2 Σ−1 be segment–wise orthogonal distances to a reference curve Γ = ∪iΓi(with Mahalanobis metric). Define the softmin tube Ltube(z) = −1 βlog X i e−βdi(z).(16) Lemma II.1 (Softmin approximation).For all zand β > 0, min idi(z)≤Ltube(z)≤min idi(z) + log n β, n = #{i}. Moreover, ∇Ltube(z) = Piwi(z)∇di(z)with wi(z) = e−βdi(z) Pje−βdj(z). Proof. The bounds are standard for log–sum–exp; the gradient follows by differentiation under the log. Relevance. Equations (14) – (16) supply the optimization–to–dynamics dictionary we use in (3) – (4) . The novelty in our setting is the categorical energy Land the explicit reachability constraint set. F. Topological and sheaf–theoretic methods Persistent homology (PH). Given a filtered simplicial complex {Kα}α∈R (e.g., Vietoris–Rips on contact graphs with filtration by threshold), its homology Hk ( Kα )changes at critical scales; the persistence diagram Dk records birth–death pairs. The bottleneck distance dB ( Dk, D0 k )metrizes diagram proximity and is stable under perturbations [42, 197, 198]. Theorem II.2 (Stability of persistence) . Let f, g be tame filter functions on a fixed complex. Then dB(Dk(f), Dk(g)) ≤ kf−gk∞for all k. This yields well–posedness for LPH 2as a shortfall (distance) to a target profile. Sheaves, gluing, and Laplacians. A cellular sheaf F on a cell complex X assigns a vector space F ( σ )to each cell σ and restriction maps along face relations; global sections satisfy linear compatibility constraints. Sheaf cohomology measures obstructions to gluing local data into a global section. A sheaf Laplacian ∆ F yields harmonic sections; residual norms krk of a least–squares gluing problem quantify inconsistency [ 44 , 45 , 197 ]. We instantiate Lglue 2 as a norm of such residuals on a fixed functional patch cover (RBM, RBD, S1/S2, glycans). Relevance and limits. TDA and sheaves provide formal causes—global shape and local consistency— but do not produce dynamics nor enforce mutation reachability. Our calculus uses PH/sheaf terms as level L2 invariants inside L, thereby shaping the flow while preserving tractability. G. Conceptual gap and our positioning The literatures above excel at individual facets: fitness and mutation (§II A); antigenic geometry (§II B); tree–based inference (§II C); genotype–phenotype modeling (§II D); optimization–dynamics equivalences (§II E); and structural invariants (§II F). None, however, jointly (i) encode final causes as cross–level invariants with explicit coherence deficits, (ii) produce a computable flow that minimizes those deficits, and (iii) enforce reachability at the level of sequences, structures, and routines. Our categorical teleonomical calculus fills this gap by combining the Bregman/JKO viewpoint with category–level invariants and a mutation–graph corridor, resulting in a single objective L and an update rule (3) that deliver both channel prediction and jump diagnostics. Summary. Classical fitness landscapes describe selective pressures on individual traits. The teleonomical update introduced here extends this idea to multiple mechanistically coupled traits, while retaining an interpretable biological meaning. III. CATEGORICAL TELEONOMY: LEVELS, CAUSES, AND CALCULUS A. Levels L0–L5 and Aristotle’s four causes We formalize Aristotle’s four causes within a categorical architecture for evolutionary systems. The material cause records sequence and structure (residues, contact graphs, patch features). The formal 9 cause is the organization of these materials into categorical objects and morphisms (sheaves over patch covers, routines, and path order). The efficient cause is realized by mutation/selection operators and life–cycle transitions acting as endomorphisms or Markov kernels. The final cause is encoded by invariants the system seeks to keep true; violations are measured as coherence deficits. We build these causes into six levels L0–L5and derive their mathematical properties. State space and enrichment. Let Z be the state category of an evolving system, whose objects z = ( x, s ) consist of a trait vector x∈Rd and a structured object s (e.g., a contact graph with patch data, or a life– cycle kernel). We assume Z is enriched in the Lawvere quantale ([0 ,∞ ] ,≥, + , 0) so that every hom–object cost ( z→z0 ) ∈ [0 ,∞ ]quantifies a short–step cost (e.g., a Bregman or Mahalanobis dissimilarity) [ 46 , 47 ]. This allows us to treat proximal updates and reachability constraints in a unified categorical–metric setting. Level invariants and coherence deficits. For each level i∈ {0,...,5}we specify: 1. a category Ci of level– i data (e.g., contact graphs with filtrations, cellular sheaves, or Markov kernels); 2. an assembly functor Ai:Z → Ci; 3. a level invariant Ii:Ci→[0,∞]that vanishes on ideal objects; and define the coherence deficit Li(z;λ) := Ii◦Ai(z;λ)∈[0,∞],(17) where λdenotes environment (e.g., immune era). The total deficit is then L(z;λ) = 5 X i=0 wiLi(z;λ), wi≥0.(18) Lemma III.1 (Functoriality and isomorphism invariance) . If F : Ci→ Ci is an isomorphism, then Li is invariant under transport: Li ( z ) = Li ( z0 )whenever Ai ( z0 ) ∼ =FAi ( z )  . Consequently, L depends only on the level–wise isomorphism class of assembled objects. Proof. If Ai ( z0 ) ∼ =F ( Ai ( z )) and Ii is invariant under isomorphisms (true for the functionals we define below), then Li ( z0 ) = Ii ( Ai ( z0 )) = Ii ( F ( Ai ( z ))) = Ii ( Ai ( z )) = Li ( z ). Summing over i yields the claim for L. Mapping Aristotle’s causes to categorical data Material cause. The material of evolution comprises sequence–level edits and structural degrees of freedom: amino acid sequences, residue–contact graphs G = ( V, E )with weights, patch features, and biophysical parameters (e.g., binding and cleavage readouts). These form the object part of Z. Formal cause. The form is the organization of material into categorical objects: (i) a cellular sheaf F over a patch cover capturing local–to–global consistency; (ii) a filtered complex for persistent homology; (iii) a small category (or Markov kernel) encoding life–cycle routine and its composition laws; (iv) a path groupoid for order–of–operations. These are the targets Ciof the functors Ai[48–50]. Efficient cause. The efficiency is represented by endomorphisms on Z (mutation, recombination) and by Markov kernels on routine states (attach → prime → fuse → replicate → egress), together with selection that modulates step acceptance/probabilities [51]. Final cause. The final cause is the drive to keep level invariants satisfied. This is quantified by Li and aggregated by L; evolutionary updates minimize Lover reachable futures (cf. (17)–(18)). Level definitions and mathematical properties L0 (viability). Let F ( z ; λ )be a baseline sufficiency index (e.g., epidemiological growth surrogate) and F?(λ)a viability threshold. Define L0(z;λ) = F?(λ)−F(z;λ)2 +,(u)+= max{u, 0}.(19) 16 Cleavage: S1/S2 and S20activation Proteolytic activation by furin–like enzymes at S1/S2 and S2 0 and usage of TMPRSS2 modulate entry route choice [ 70 , 71 ]. Cleavage readouts include Western blot densitometry (S1/S2:S ratio), targeted MS peptide ratios, and cell–based reporter assays. Let p`r ( v ) ∈ (0 , 1) be the estimated cleaved fraction. A logit transform stabilizes variance: gcleavey:= logit p`r(v)= log p`r(v) 1−p`r(v).(45) We standardize by subtracting the reference logit and (optionally) applying a linear rescale so that Tcleave falls in an interpretable band (e.g., [0 . 4 , 1 . 5] in our experiments). When p is near 0or 1, a beta–binomial measurement model prevents infinite logits and accounts for overdispersion: c`r(v)∼BetaBinomialn`r, α(v), β(v), p`r(v) = α(v) α(v) + β(v).(46) Stability: thermal/thermodynamic proxies Spike stability can be indexed by melting temperature Tm (differential scanning fluorimetry/DSF) or by chemical denaturation free energy ∆ G [ 72 – 74 ]. Under a two–state approximation, the van’t Hoff relation gives ∆G(T)≈∆Hm1−T Tm−∆Cph(Tm−T)−Tlog Tm Ti,(47) with ∆ Hm enthalpy and ∆ Cp heat–capacity change at Tm [ 75 ]. For standardized comparison at a common T◦, we define gstaby:= ∆Gv(T◦)−∆Gref (T◦)or Tm,v −Tm,ref.(48) Because DSF signals can be dye– and buffer–dependent, we include lab random effects and allow for replicate–wise heteroscedasticity (larger error near broad transitions). Escape: neutralization fold–change Let N(ID50) `r ( v )be a neutralization titre (ID 50 ) from a serum panel (live or pseudovirus; monoclonal or polyclonal). The fold reduction relative to the reference is FR`r(v) := N(ID50) `r (vref) N(ID50) `r (v).(49) We take gescapey:= log2FR`r(v)or FR`r(v),(50) depending on whether a logarithmic or raw fold–scale is desired downstream. To harmonize across labs and assay types, we calibrate titres to the WHO International Standard (BAU/mL) when available and include lab/assay random effects [ 76 – 79 ]. Because titre noise is multiplicative, we model log N as Gaussian with lab–specific variance. Class-resolved escape traits. Antibody escape is epitope-specific: mutations at the RBD ridge (e.g., E484, F486) primarily affect Class 2 antibodies, whereas changes at K417 or the R346/K444 region differentially impact Class 1 and Class 3 antibodies. (See, e.g., Barnes et al. and Starr/Greaney deep mutational scanning maps for the canonical class definitions and footprints.) A single scalar escape coordinate would therefore conflate phenotypically distinct directions. To address this, we represent escape as a vector of class-resolved components rather than a single number. For a genotype g we write Tesc(g) = TC1(g), TC2(g), TC3(g), TNTD(g),(51) where 17 •TC1(g)quantifies escape from Class 1 RBD antibodies (ACE2-competitive, RBD “up” footprint); •TC2(g)quantifies escape from Class 2 RBD antibodies (ridge-centered, E484/F486 epitope); •TC3(g)quantifies escape from Class 3 RBD antibodies (R346/K444/outer-face epitope); •TNTD(g)quantifies escape from NTD-directed antibodies. Operationally, each component is obtained by aggregating mutation-level escape from class-specific deep mutational scanning (DMS) maps and neutralization panels. Let Ea ( g )denote the escape score of genotype g against monoclonal antibody a in a DMS dataset, and let ACk be the set of antibodies assigned to class Ck (using antibody footprints and Barnes-class annotations). We define TCk(g) = 1 |ACk|X a∈ACk Ea(g), k ∈ {1,2,3},(52) and analogously for the NTD component using NTD-directed antibodies. When polyclonal serum data are available for a given panel, we convert fold-reductions in ID 50 to log -scale using (50) and fit a per-panel affine map that aligns the DMS-derived TCk with the panel’s scale (hierarchical calibration as in (53) ). This yields a consistent, class-resolved escape vector Tesc ( g )whose components enter the trait vector x = ( Tbind, Tcleave, Tstab, Tesc )via a chosen aggregation (e.g. a weighted norm or projection) in L1 , while the full vector is used in the epistasis analysis and Class-specific path projections (Sec. IV B). In particular, TC2 captures the Class 2 ridge evolution emphasized by Reviewer 1, while TC1 and TC3 allow us to distinguish K417and R346/K444-driven antigenic changes from ridge-focused modifications. Hierarchical calibration across laboratories For each trait j, we adopt a random–effects meta–analytic model [80, 81]: gjy(j) `r (v)=θj(v) + δ(j) `+(j) `r (v), (j) `r (v)∼ N0, σ2(j) `,(53) δ(j) `∼ N(0, τ2 j), σ(j) `∼HalfCauchy(aj), θj(vref )=0(identification). Posterior means ˆ θj ( v )and covariances yield point estimates and uncertainties for Tj ( v ). When multiple assay modalities exist for a trait (e.g., live vs pseudovirus neutralization), we add a modality random effect and, if needed, a ratio–of–means calibration anchored by shared panels (Bland–Altman concordance checks [82]). Proposition IV.1 (Identifiability of trait means) . Assume each variant v (including vref ) is measured by at least two laboratories and the bipartite lab × variant graph is connected. Then the trait means {θj ( v ) } and lab intercepts {δ(j) `}are identifiable up to the imposed constraint θj(vref)=0. Proof. The model (53) is a two–way additive random–effects model. Connectivity ensures full column rank of the design matrix up to the single sum–to–zero constraint implemented by fixing θj ( vref ) = 0). Standard linear mixed model identifiability results apply. From calibrated traits to teleonomical geometry Let ˆx ( v )=( ˆ θbind ( v ) ,ˆ θcleave ( v ) ,ˆ θstab ( v ) ,ˆ θescape ( v )) and Cov ( x ( v )) = Σ (v) be the posterior covariance from (53) . To define anisotropic penalties in L (e.g., in L1 and Ltube ), we set a global trait covariance Σ by a robust average: Σ := MedianΣ(v)v+ diag(σ2 bind, σ2 cleave, σ2 stab, σ2 escape),(54) where the diagonal term inflates to cover between–variant heterogeneity. Mahalanobis distances kx− τ ( λ ) k2 Σ−1 and tube orthogonal distances di ( z ) = 1 2kz− Π i ( z ) k2 Σ−1 (cf. (38) ) then naturally weight traits by their calibrated uncertainty. 18 Table II: Double-mutant epistasis validation for binding and escape traits. Quantities are trait deltas relative to the background using the calibrated coordinates of §IV.A. Full panels are provided in Appendix K. Background Sites (mut.) Trait t∆T(i)∆T(j)∆T(i,j)∆epi Wuhan-like Q498R, N501Y ACE2 (−log10 KD)≤0 +1.07 +1.398 ≥0.328 BA.2-like R346T, F486P Class 2 escape (0–1) 0.0147 0.5622 0.5690 −0.0079 Assay physics and error models (biological details) Binding. SPR/BLI signals are governed by convection–diffusion to the sensor surface and mass action at the interface. Deviations from simple 1:1 kinetics (heterogeneous ligands, mass transport limitation) produce biased KD if unmodeled [ 66 , 67 ]. We mitigate this with residual diagnostics and by down–weighting runs with Damkohler number outside the linear regime (mass transport–limited). Cleavage. S1/S2 cleavage reflects furin recognition of the polybasic motif and local accessibility; S2 0 and TMPRSS2 usage depend on priming and membrane localization [ 70 , 71 ]. Because cleavage fractions are bounded, the beta–binomial (46) captures overdispersion from band quantification and peptide ionization variability. Stability. Thermal unfolding of spike ectodomain is often apparently two–state due to cooperative trimer unfolding; dye–binding DSF reports on hydrophobic core exposure [ 72 ]. Equation (47) provides a principled transform from Tm to ∆ G ( T◦ )when enthalpy/heat–capacity shifts are estimated or controlled. Escape. Neutralization titres are log–normally distributed; fold–reductions (49) therefore reside naturally on a log scale. Calibration to BAU/mL via WHO standards reduces inter–assay bias and enables cross–study synthesis [ 76 ]. Live versus pseudovirus assays are reconciled by a modality effect with partial pooling [77, 78]. Putting it together: data–to–trait pipeline 1. Ingest: collect raw spreadsheets for binding (SPR/BLI), cleavage (WB/MS), stability (DSF/chem. denaturation), and neutralization (ID50/IC50). 2. Transform: apply gjper trait using (43), (45), (48), (50). 3. Anchor: subtract the reference vref per lab/assay. 4. Calibrate: fit (53) (one model per trait), yielding ˆx(v)and Σ(v). 5. Assemble: form Σvia (54) ; export trait table with uncertainties to downstream teleonomical optimization (continuous or graph–based). This standardized, uncertainty–aware trait representation is what feeds L1 and the tube geometry in §III B, and it provides the empirically grounded anchors for L2/L3/L4 via their dependence on trait–linked structural and routine objects. Interpretation. The four trait axes used here correspond to mechanistic properties of the spike: binding, cleavage, stability, and class-specific immune escape. These traits capture the essential biophysical pressures acting on the virus and are far fewer than in typical machine-learning models; they form a minimal set sufficient to summarize the observable constraints on SARS–CoV–2 evolution. B. Epistasis and trait coupling Motivation. Non-additive interactions among spike residues (epistasis) and mechanical coupling between molecular processes jointly shape adaptation. For RBD, pairs such as Q498R×N501Y produce super-additive ACE2 affinity, while combinations like R346T × F486P modulate class 2/3 antibody escape with compensatory effects on binding. At the process level, S1/S2 priming and trimer stability are linked via conformational equilibria that govern RBD opening and S1 shedding. We therefore make epistasis explicit in the sequence→trait map and couple specific traits within the L1term. 19 Trait vector and encoding. Let g∈ ALdenote a spike sequence. We predict a trait vector T(g) = Tbind(g), Tcleave(g), Tstab(g), Tesc(g)∈RdT,(55) where Tesc ( g )can be scalar or epitope-resolved Tesc ( g ) = ( θC1, θC2, θC3, θNTD )(cf. Sec. IV A). The sequence is encoded by per-site features φi ( g )(one-hot or a compact physicochemical basis), concatenated into φ(g). Global epistasis with targeted pairwise/triad terms (answers Rev. 1.1b/1.1c). For each trait component t∈ {1, . . . , dT}we use a global-epistasis link utacting on a sparse additive + interaction model: Tt(g) = ut βt,0+X ihβt,i, φi(g)i+X (i,j)∈E2 hβt,ij, φi(g)⊗φj(g)i+X (i,j,k)∈E3 hβt,ijk, φi⊗φj⊗φki .(56) Here E2 and E3 are curated interaction sets. We seed them by three orthogonal criteria: (i) structural proximity (e.g., C α distance threshold/contact map), (ii) co-occurrence/co-evolution along observed lineages, and (iii) DMS-based interaction scores. The link utis a monotone I-spline with Ktknots, ut(z) = Kt X k=1 αt,k Ik(z), αt,k ≥0,(57) which captures global nonlinearity (saturation, floor/ceiling effects) while preserving order. Strong heredity and multi-task sharing. To enforce interpretability and parsimony we use strong heredity and multi-task sharing: Ωepi =λ1Pt,ikβt,ik2+λ2P(i,j)∈E2βt,ij t2,1 Ptkβt,ik2+Ptkβt,j k2+ε +λ3P(i,j,k)∈E3βt,ijkt2,1 Ptkβt,ik2+Ptkβt,j k2+Ptkβt,kk2+ε. (58) The mixed `2,1 norm encourages the same interaction to be reused across traits when supported by data (e.g., 498–501 affecting both binding and escape via exposure changes). Background dependence (state gating). Some interactions are active primarily in specific conformational or sequence backgrounds (e.g., RBD-up fraction). We allow gated coefficients: βt,ij(g) = σγ> t,ijh(g)˜ βt,ij,(59) where h ( g )are low-dimensional background features (clade markers, coarse structural state proxies), σ is logistic, and ˜ βt,ij are base coefficients. This captures background-dependent epistasis without exploding parameter count. Explicit epistatic deviation and its energy. For trait t and pair ( i, j )define the deviation at background g: ∆epi t,ij(g) = Tt(g(i,j))−Tt(g)−Tt(g(i))−Tt(g)−Tt(g(j))−Tt(g).(60) Given a prior πt,ij (sign/magnitude) and uncertainty weights ωt,ij, we add Lepi(g) = X tX (i,j)∈E2 ωt,ij ∆epi t,ij(g)−πt,ij2+X tX (i,j,k)∈E3 ωt,ijk ∆epi t,ijk(g)−πt,ijk2,(61) and include wepiLepi additively in L . With no prior, set π· = 0 and inflate uncertainty. This answers Rev. 1.1c by making pairwise/triad epistasis explicit in the objective. Mechanistic trait coupling in L1 (answers Rev. 1.1a). To encode stability–cleavage coupling, L1 includes a cross-trait directional-consistency term computed along reachable genotype directions: Lstab-cleave 1(g) = λSC ∇reachTstab(g)−κ∇reachTcleave(g)2 M.(62) Here ∇reach is the directional derivative restricted to edges ( g→g0 ) ∈E of the mutation graph, and k·kMis a structurally weighted metric (e.g., emphasis on directions that preserve glycan occupancy or cleavage-loop integrity). The scale κ is estimated from observed co-variation between ∆ Tstab and ∆ Tcleave across single/double mutants (App. J). Analogous couplings can be specified for (binding, escape) to reflect trade-offs. 20 Table III: Pairwise epistasis for ACE2 binding in the Wuhan-Hu-1 background (units: −log10 KD ). Epistasis is computed as ∆epi = ∆T(i,j)−∆T(i)−∆T(j). Background Sites (mut.) Trait t∆T(i)∆T(j)∆T(i,j)∆epi Wuhan-like Q498R, N501Y ACE2 (−log10 KD)≤0a+1.07b+1.398c≥0.328d a Q498R single in Wuhan reported as mildly deleterious (no precise value). b N501Y single in Wuhan: ∆ log10 KD = +1 . 07. c Double Q498R+N501Y ≈ 25 × tighter vs WT; log10 (25) = 1 . 398. d Lower bound: 1.398 −(1.07 + 0) = 0.328; larger if Q498R is deleterious. Table IV: Pairwise effects on Class 2 immune escape in a BA.2 -like background, using DMS escape fractions aggregated over Class2 (Barnes group B) monoclonal antibodies. Single-mutant values ∆ T(i) and ∆ T(j) are computed as the mean mutation-level escape fractions across all Class 2 antibodies (Section IV.A). The doublemutant escape fraction uses the independent-sites, saturating baseline Eij = 1 − (1 −Ei )(1 −Ej ); epistasis is defined as ∆epi = ∆T(i,j)−∆T(i)−∆T(j). Background Sites (mut.) Trait t∆T(i)∆T(j)∆T(i,j)∆epi BA.2-like R346T, F486P Class 2 escape (DMS fraction) 0.0147 0.5622 0.5686 −0.0083 Directional derivatives on the corridor. For an allowed edge ( g→g0 )define the unit direction eg→g0 and directional derivative ∇reachTt(g)[eg→g0] = Tt(g0)−Tt(g) =: ∆Tt(g→g0).(63) Gradients and cross-gradients in (62) are thus well-defined without requiring a global Euclidean embedding of genotypes. Integration into the teleonomical step and marginal gains. Let ∆ L• ( g→g0 )denote the change of each level when traversing an edge. The marginal effective ascent used in the proximal update (Sec. II A) reads ∆Weff(g→g0) = −∆L(g→g0) = −∆L0+ ∆L1+ ∆L2+ ∆L3+ ∆L4+ ∆L5+wepi ∆Lepi.(64) Thus an apparently attractive single mutation can be downweighted if it contradicts supported epistatic patterns; conversely, a two-step path can be promoted when the pair is synergistic (handled by the beam and step-cap mechanics; cf. App. G). Fitting objective, losses, and identifiability. Coefficients {βt,·} and link parameters {αt,k} are estimated by minimizing J=X tX m∈Dt ρ Tobs t(m)−Tpred t(m)+ Ωepi + Ωsmooth(ut),(65) where Dt collects trait measurements (ACE2, cleavage, stability, escape axes), ρ is a robust (Huber) loss to handle inter-lab outliers, and Ω smooth penalizes curvature of ut . Strong heredity (58) resolves the key identifiability issue (interactions cannot be large if parents are zero). Hyperparameters ( λ1, λ2, λ3 )are selected by nested cross-validation that holds out backgrounds (distinct haplotypes) and double mutants to test non-additivity generalization. Selection and expansion of interaction sets. We initialize E2 with pairs that satisfy any of: (i) structural contact/proximity; (ii) positive co-occurrence enrichment along phylogeny; (iii) DMS-based interaction score above threshold. After fitting, we allow a single-step forward selection that adds the top K pairs improving validation loss under heredity, while controlling FDR across traits. A similarly conservative rule is used for E3. Parsimony across traits We couple traits via the multi-task penalty in (58) , which shrinks unsupported interactions to zero across all traits. In sensitivity analyses (App. G) only a small set of parameters (tube sharpness, progress schedule, w2 , step caps, beam width, and a handful of βt,ij ) dominate outputs, directly countering the “could explain anything” concern. Validation and prospective falsifiability. We pre-register a compact panel of constructs: (i) ACE2 affinity: Q498R+N501Y on two backgrounds (with/without 486P); (ii) escape axes: R346T+F486P and K444N+R346T against class-specific monoclonals and sera cohorts; (iii) stability/cleavage coupling: a furin-sensitizing change paired with an S2’ stabilizer. For each we report ∆ Tt and ∆ epi t,ij , with rank-ordered predictions for IC50, infectivity, and DSF/DSC proxies; see App. K and App. G. 21 Interpretation of Table IV. Table IV reports pairwise effects of R346T and F486P on the Class 2 component of the escape vector in a BA.2-like background. The single-mutant values ∆ T(i) and ∆ T(j) are the Class 2 escape fractions obtained by aggregating deep mutational scanning escape scores over the Class 2 (Barnes group B) monoclonal antibodies in the Bloom-lab maps, as described in Section IV.A. Thus all three entries in the row live on the same, DMS-defined Class 2 axis TC2 ( g )and do not mix different panels or backgrounds. For the double mutant, we use the independent-sites, saturating baseline Eij = 1 − (1 −Ei )(1 −Ej ), where Ei and Ej are the Class 2 escape fractions for R346T and F486P respectively. The quantity ∆ T(i,j) in the table is this baseline prediction expressed on the same DMS Class 2 scale, and the epistasis entry ∆ epi = ∆ T(i,j)− ∆ T(i)− ∆ T(j) reports how far the double deviates from additivity on that axis. The small negative ∆ epi reflects the expected saturation of escape fractions near 1 rather than a strong antagonistic interaction: F486P accounts for most of the Class 2 escape in BA.2, while R346T adds only a minor increment. Because the table is built entirely from a single DMS data family and a well-defined Class 2 escape coordinate, it is coherent with the four-dimensional escape vector introduced in Section IV.A and avoids any mixing of heterogeneous neutralization panels or assay scales. (a) Mechanistic trait coupling between stability and cleavage is encoded via the cross-gradient penalty Lstab-cleave 1 in (62) using reachable-direction derivatives and a structural metric; (b) epistatic interactions enter the sequence → trait predictor f directly via the global-epistasis model (56) with heredity (58) and background gating (59) ; (c) pairwise and higher-order epistasis contribute explicitly to the coherence deficit through the deviation (60) aggregated in the objective (61) , and they modify marginal gains as in (64). Validation is planned on double/triple mutants with registered predictions (Table II, App. K). C. L0–L5 instantiation We now instantiate the level deficits L 0– L 5for SARS–CoV–2 in the calibrated trait space x = ( Tbind, Tcleave, Tstab, Tescape )and with associated structural/routine objects s (§IV A). Unless otherwise stated, all losses are nonnegative, lower semicontinuous in their arguments, and measured with trait anisotropy through the global covariance Σ(Mahalanobis metric). Throughout, λ denotes the immune/epidemiological era, treated as piecewise constant. L1: trait bands with directional targets by era Let τ ( λ ) ∈R4 be an era–specific directional target and Σ 1 ( λ )  0an anisotropy matrix. We define the interior quadratic plus hinge outside biologically plausible bands: L1(z;λ) = 1 2(x−τ(λ))>Σ1(λ)−1(x−τ(λ)) + 4 X j=1 κ− j(λ) [`− j(λ)−xj]2 ++κ+ j(λ) [xj−u+ j(λ)]2 +, (66) where [ · ] + is the positive part, and ( `− j ( λ ) , u+ j ( λ )) are era–dependent trait bands (e.g., modest recovery of cleavage while maintaining high escape in Omicron eras). By construction, L1 ( · ; λ )is convex in x and admits the closed–form gradient inside the band: ∇xL1= Σ1(λ)−1(x−τ(λ)). Lemma IV.2 (Directional interior and safe boundary) .L1 furnishes an interior directional pressure toward τ ( λ )while imposing quadratic penalties on boundary excursions. In particular, L1 is C1 on the interior and piecewise C1globally, strictly convex when Σ1(λ)0and at least one hinge is active. Proof. Immediate from the sum of a strictly convex quadratic and convex hinge squares; differentiability follows from standard properties of hinge functions. L2 (per variant): contact–graph PH shortfall and sheaf gluing variance Persistent homology shortfall. Let G ( z )be the residue–contact graph (weighted) extracted from s , with a filtration f (e.g., threshold on contact weight). For homology degrees k∈ { 0 , 1 , 2 } , let Dk ( f )be 22 the persistence diagrams; let D? k ( λ )be era–specific targets (or smooth envelopes inferred from structural databases). We select a stable, differentiable discrepancy d℘ between diagrams, such as the sliced Wasserstein distance or kernels based on persistence images [83–85] and define LPH 2(z;λ) := X k ωk(λ)d℘ Dk(f), D? k(λ).(67) With d℘ differentiable almost everywhere (a.e.) in the underlying point clouds, LPH 2 is differentiable a.e. in the contact weights and, by stability of the chosen metric, Lipschitz in small perturbations of f. Sheaf gluing variance on a functional patch cover. Let U = {RBM,RBD,S1/S2,glycans, . . . } be a fixed cover of spike functionality. A cellular sheaf F ( z )assigns local feature spaces to patches and linear restriction maps to overlaps; global sections encode cross–patch compatibility. Let r ( z )be the residual from the least–squares gluing problem (sheaf Laplacian normal equations). We set Lglue 2(z;λ) := 1 2kr(z)k2 2,(68) which is convex quadratic in the local data (hence differentiable), and zero iff a compatible global section exists over U. Per–variant structural deficit. The per–variant structural level combines the two: L2(z;λ) := LPH 2(z;λ) + Lglue 2(z;λ), L2≥0.(69) Theorem IV.3 (Stability and subgradient representation) . Assume d℘ is Lipschitz and a.e. differentiable w.r.t. the underlying point clouds and that the sheaf gluing residual depends smoothly on local features. Then L2is Lipschitz on compact sets and admits a subgradient everywhere; where differentiable, ∇L2=X k ωk∇d℘+J>r(z), with Jthe Jacobian of the sheaf residual map. Proof. Lipschitzness follows by composition of Lipschitz maps (stability of d℘ and linearity of residuals). Differentiability a.e. and subgradient existence follow from Rademacher’s theorem and convexity of the quadratic residual. L3: routine kernel (life–cycle) from assay panels Let S= {attach,prime,fuse,replicate,egress} be routine states and P ( z ) ∈R|S|×|S| a row– stochastic kernel constrained by admissible transitions (e.g., no direct attach→replicate ). Denote by π ( z )a stationary distribution (assume irreducible). Define a throughput functional on productive edges Eprod ⊂S×S, ΘP(z):= X (u→v)∈Eprod πu(z)Puv(z),(70) and a target Θ?(λ). The routine deficit is L3(z;λ) := Θ?(λ)−Θ(P(z))2 ++γcondP(z),(71) where cond is a convex surrogate of mixing bottlenecks (e.g., 1 −σ2 ( P )or a conductance proxy) and γ≥0[86, 87]. Estimating P ( z )from assay panels. Let Cuv be counts of transitions observed in panel experiments (e.g., time–lapse microscopy or protease perturbations), with multinomial likelihood (Cu•|P)∼MultinomialNu,(Puv)v,b Puv =Cuv +α Pv0(Cuv0+α),(72) with symmetric Dirichlet prior concentration α > 0for smoothing [ 88 ]. Then P ( z )can be set to b P (plug–in) or estimated jointly with z by alternating minimization. Under standard regularity, b P→P? almost surely as Nu→ ∞ (consistency). Lemma IV.4 (Convexity in P ) . For fixed π and convex cond , the map P7→ L3 in (71) is convex over the polytope of row–stochastic matrices with fixed support. The gradient of Θw.r.t. P is ∇Puv Θ = πu on productive edges and zero otherwise. Proof. Linearity of Θin P implies convexity of the squared hinge; add a convex cond to preserve convexity. The gradient statement is immediate from (70). 23 L4: holonomy from two–arm experiments Let E be a set of environmental operations (e.g., immune boost , antiviral ), implemented as endofunctors E:E → End(Z). For (a, b)∈ E ×E and state z, define the holonomy discrepancy ∆a,b(z) := E(a)◦E(b)(z)−E(b)◦E(a)(z),(73) measured by a norm k·kΣ−1 on traits and an appropriate distance on structure. Aggregate over a generator set G: L4(z;λ) := X (a,b)∈G ηab(λ)∆a,b(z)2.(74) Estimating holonomy from two–arm order experiments. Given paired assays measuring zab and zba (apply a then b vs b then a ), the empirical discrepancy b ∆a,b = zab −zba produces an unbiased estimator of k∆a,bk2; Hoeffding–type concentration gives Pkb ∆a,bk2−k∆a,bk2≥≤2 exp−2n2 B2, for n paired replicates and bounded variance proxy B2 [ 89 ]. Thus L4 can be controlled statistically with modest panel sizes. L5: multi–level coherence variance To prevent single–level overfitting, we penalize cross–level inconsistency using a variance–like coupling among normalized deficits ˜ Li=Li/σi: L5(z;λ) := 1 2X i<j αij(λ)˜ Li(z;λ)−˜ Lj(z;λ)2 ++ρX i βi˜ Li(z;λ)2,(75) with nonnegative weights αij, βi and ρ≥ 0. The first term discourages large imbalances across levels; the second enforces a ridge–type regularity across the vector of deficits (cf. scalarizations in multiobjective optimization [90]). Proposition IV.5 (Convexity and balancing) . If all ˜ Li ( · ; λ )are convex, then L5 ( · ; λ )is convex. Moreover, at any stationary point in a convex reachable set, active pairs ( i, j )with αij > 0satisfy ˜ Li≈˜ Lj up to the hinge tolerance, preventing pathological concentration on a single level. Proof. Each squared hinge of an affine combination of convex functions is convex; the ridge term is a convex quadratic; stationarity with active hinges implies balance conditions by first–order optimality. Learning the level weights wiby ranking with regularization We learn nonnegative weights w = ( w0, . . . , w5 )using pairwise ranking constraints derived from known orderings (e.g., BA.2 ≺ BA.1 in deficit, BA.5 ≺ BA.2, XBB ≺ BA.5). Let P be a set of ordered pairs ( u≺v ), meaning L ( zu ; λ ) <L ( zv ; λ ). With feature vectors φi ( u ) := Li ( zu ; λ )we pose a large–margin RankSVM: min w,ξ 1 2kwk2 2+CX (u≺v)∈P ξuv (76) s.t. 5 X i=0 wiφi(v)−φi(u)≥1−ξuv, ξuv ≥0, wi≥0∀i, (u≺v). Alternatively, one may use pairwise logistic (Bradley–Terry) or Plackett–Luce likelihoods with an `2 penalty [91–94]. 24 Theorem IV.6 (Existence, uniqueness, and monotonicity).If the constraint set in (76) is feasible and at least one pair ( u≺v )yields a strict inequality for some nonnegative w , then (76) has a unique optimal solution w? . Moreover, if constraints encode a partial order consistent with a topological sort of variants, the learned Lis strictly monotone along that order on the training states. Proof. Strong convexity of the objective in w with linear constraints gives uniqueness. Monotonicity follows from the margin constraints and nonnegativity of w. Regularization across levels. To encourage parsimonious contributions, one can add a (group) sparsity penalty on w: min w,ξ 1 2kwk2 2+λkwk1+CXξuv subject to (76) constraints, or group penalties if levels are clustered (e.g., fold structural vs routine) [ 95 ]. This yields interpretable load–bearing levels. Summary of L0 and integration L0 is a baseline viability floor based on epidemiological sufficiency (e.g., next–generation threshold R0≥ 1), implemented as a convex hinge on a surrogate F ( z ; λ )estimated from growth proxies; formal next–generation matrix methods provide principled thresholds [96]. The teleonomical energy L(z;λ) = 5 X i=0 wiLi(z;λ) then combines L1–L5 (and L0) into a single, computable objective that (a) pulls toward era–appropriate interior trait regions, (b) enforces structural/routine/holonomy coherence, (c) balances levels, and (d) respects known evolutionary orderings through the learned w. D. Channel construction and optimization We now specialize the teleonomical flow of §III B to an observed evolutionary corridor joining the Omicron anchors {BA.1,BA.2,BA.5,XBB} in trait space x = ( Tbind, Tcleave, Tstab, Tescape ). We first define asoft tube around the polyline Γwith soft–min over segments and β –annealing, then introduce a global progress S ( z )that aggregates segment index and local parameter, and finally describe a proximal–gradient with trust region and a smoothness penalty (second differences) with asymmetric step caps to avoid spikes while respecting reachability. We close with an era–aware framework to learn environment targets/weights from time series. For evaluation, however, we never include the held-out segment in the tube: as detailed in VI C, the polyline Γis rebuilt per split to exclude the test window and ensure split-integrity. Polyline, segment projection, and soft tube Let Γbe the polyline through the observed anchors Γ := m−1 [ i=0 Γi,Γ0: BA.1→BA.2,Γ1: BA.2→BA.5,Γ2: BA.5→XBB,(77) with m = 3 segments. For the i –th segment, let ( ai, bi ) ∈R4×R4 be its endpoints. Given a state z= (x, s), write the Mahalanobis projection of xonto Γias Πi(x) = ai+ti(x)vi, vi:= bi−ai, ti(x) := clip[0,1]hx−ai, viiΣ−1 hvi, viiΣ−1,(78) and define the orthogonal distance di(z) := 1 2kx−Πi(x)k2 Σ−1.(79) 25 The soft tube of §III B specializes to Ltube,β(z) := −1 βlog m−1 X i=0 exp−β di(z), β > 0,(80) with segment weights wi(z;β) = exp(−βdi)/Pjexp(−βdj)and gradient ∇xLtube,β(z) = X i wi(z;β) Σ−1x−Πi(x).(81) By log–sum–exp smoothing, Ltube,β ↓minidi as β↑ ∞ uniformly on compacta, and ∇xLtube,β is continuous for any finite β(cf. Lemma II.1). Proposition IV.7 (Annealed segment selection) . Fix a compact set K⊂R4 . Then Ltube,β →minidi epi– converges on K as β→ ∞ , and wi ( z ; β ) →1{i∈arg minjdj ( z ) }/ # arg min d pointwise. Consequently, for increasing βthe optimizer transitions smoothly from blended to hard segment selection. Proof. Epi–convergence of log–sum–exp to pointwise minima is classical (variational analysis); the weight limit is a direct corollary of softmax concentration [97]. Global progress coordinate and scheduled targets Define the global progress coordinate that aggregates segment index and local parameter: S(z) := m−1 X i=0 wi(z;β)i+ti(x)∈[0, m].(82) To provide forward drive along Γ, we penalize deviation from a scheduled target Star(k)at iteration k: Lprog(z;Star) := λprogS(z)−Star2, λprog >0,(83) with a nondecreasing schedule Star ( k ) ∈ { 0 , 1 , 2 , 3 } (e.g., updating every K steps) and a β –schedule β(k)%βmax (annealing). Differentiating (82) gives (formally) ∇xS(z) = X i wi∇xti(x) + X ii+ti(x)∇xwi(z;β), where the second term is smooth for finite β. In practice, we use the tangent surrogate gprog(z) := 2λprogS(z)−Starteff(z), teff (z) := Piwi(z;β)vi kvik Piwi(z;β)vi kvik ,(84) which drives progress along the blended segment direction while avoiding the nontrivial Jacobian of the projection x7→ ti(x)at kinks of Γ. Lemma IV.8 (Monotone progress under acceptance) . Suppose the acceptance rule requires Ltube,β ( z+ ) ≤ Ltube,β ( z ) + δ and Star is nondecreasing. Then any accepted step with sufficiently small trust radius and step caps satisfies S ( z+ ) ≥S ( z ) −O ( δ ); in particular, for δ = 0 and β < ∞ , S ( z )is nondecreasing along accepted iterates. Proof. For small steps the blended tangent teff changes continuously and the projection parameters ti cannot decrease substantially without increasing some di ; the tube constraint rules out such increases beyond O(δ). Discrete path, smoothness penalty, and spike control Let {zk}k≥0 be the iterates, and for convenience write xk for the trait component. We penalize discrete curvature via second differences Lsmooth({xk−1, xk, xk+1}) := λcurv xk+1 −2xk+xk−12, λcurv ≥0,(85) 32 Index 3: Topological flip (L2). Let Dk ( z )be persistence diagrams of the contact–graph filtration at state z and d℘ a stable diagram distance (e.g., sliced Wasserstein). The PH flip amplitude across an edit (u→v)is ∆PH(u→v) := X k ωkd℘ Dk(u), Dk(v).(106) Let  bound diagram noise (from bootstraps or theoretical bounds) and τPH > 0a flip threshold. By stability of d℘and confidence bands in persistence (e.g., via bottleneck bootstrap [122, 123]), we have: Proposition V.3 (Certified PH flip) . If ∆ PH ( u→v ) > τPH and τPH > 2  , then with confidence at least 1 −δ (from the bootstrap) the change in Dk is not due to noise; a genuine topological re–wiring occurred at L2 between uand v. Proof. Triangle inequality and stability: d℘ ( Dk ( u ) , Dk ( v )) ≤d℘ ( Dk ( u ) ,b Dk ( u )) + d℘ ( b Dk ( u ) ,b Dk ( v )) + d℘ ( b Dk ( v ) , Dk ( v )). If the middle term exceeds 2  by more than the bootstrap band, it certifies a real change [122]. For sheaf gluing, with residuals r(z)and Lglue 2=1 2krk2, we define ∆glue(u→v) := Lglue 2(u)−Lglue 2(v),(107) and flag when ∆ glue  0while ∆ PH is modest—indicating a compatibility restoration (previous obstruction alleviated) even without large PH changes. Index 4: Holonomy spike (L4). For operations a, b in the environment category and state z , recall the holonomy discrepancy ∆a,b(z)(73). Define the holonomy index Hmax(z) := max (a,b)∈G ηab k∆a,b(z)k,(108) and flag spikes when Hmax exceeds a threshold predicted by Lipschitz composition bounds. If each E( a ) is La–Lipschitz in the tube metric and x7→ E(a)(x)is differentiable along TzReach, then k∆a,b(z)k ≤ (LaLb+LbLa) diamorb(z),(109) where orb ( z )is the two–step orbit under {a, b} . A measured Hmax exceeding (109) up to noise indicates intrinsic noncommutativity (curvature) of the environmental connection, i.e. a true L4 spike (Ambrose– Singer holonomy principle) [124]. Index 5: Routine split (L3 conductance). Let P ( z )be the routine kernel and write the bottleneck conductance Φ∗(z) := min S⊂S,0<π(S)≤1/2Pu∈S, v /∈SπuPuv π(S).(110) Cheeger–type inequalities for Markov chains relate the spectral gap 1 −λ2 ( P )to conductance [ 125 – 127 ]: (1 −λ2) 2≤Φ∗≤p2(1 −λ2).(111) Aroutine split is signaled when a new set of productive edges E0 prod raises the throughput Θ (70) while simultaneously Φ∗increases across a cut that previously bottlenecked flow. We use the index Ξ(z) := ∆Θ + ν∆Φ∗,(112) with ν > 0, and flag when Ξexceeds a statistical threshold (estimated from bootstrap on panel counts). By (111) , a jump in Φ ∗ implies a commensurate improvement in the spectral gap and thus a re–routing of the routine. A combined jump score and false–positive control Define normalized indices on [0,1]: J1= exp(−∆V/θ),J2=σ(τλ−λreach min ),J3=σ(∆PH−τPH)+σ(∆glue−τglue),J4=σ(Hmax−τH),J5=σ(Ξ−τΞ), 33 where σ(u) = 1 1+e−ku and τ•are trait/era–specific thresholds. The jump score is JS(z) := 5 X i=1 ωiJi(z), ωi≥0,X i ωi= 1.(113) We control false positives by requiring at least two indices to exceed their thresholds and by exploiting independence/weak dependence across L2/L4/L3: under a union bound, P(false jump)≤X i<j PJi≥τiPJj≥τj≈OX i<j αiαj, if each threshold yields marginal Type–I rate αi(estimated from null bootstraps). Biological interpretation Barrier drop corresponds to the system gaining access (by mutation or recombination) to a lower–deficit basin; curvature softening indicates a fold where the present compromise of traits/routines becomes unstable; topological flips detect re–wiring of structural constraints (e.g., new salt bridges or altered glycan support), holonomy spikes reveal strong order effects (immune ◦ drug 6 =drug ◦ immune), and routine splits capture the emergence of an alternative high–throughput life–cycle path (e.g., shift between endosomal and TMPRSS2 entry). Because these indices derive from distinct categorical levels, simultaneous elevation is a robust sign that a discontinuity is imminent. B. Operational jump detector We turn the indices of §V A into a decision procedure that scans along the evolutionary corridor and flags locations where a discontinuity is likely. The detector is defined both for the continuous teleonomical path (trust–region proximal flow) and for the discrete corridor (mutation graph). We also provide statistical thresholds with false–discovery control and report candidate nodes together with their trait changes and (when available) sequence edits. Segmentation of the corridor and index vectors Let S ( z ) ∈ [0 , m ]denote the global progress coordinate (82) . Partition the interval [0 , m ]into K half–overlapping windows Wk= [sk−∆S, sk+ ∆S], k = 1, . . . , K, with centers sk and half–width ∆ S > 0chosen so that each Wk holds sufficient samples (continuous path samples or graph nodes) for local estimation. For the continuous path {z ( t ) } , define the windowed aggregates d ∆Vk:= min γ∈P(Ak→Bk)max t∈WkL(γ(t)) −min t∈WkL(z(t)),(114) b λk:= min t∈Wk λreach min (z(t)) (restricted softest curvature),(115) b ∆PH,k := median {∆PH(u→v) : S(u), S(v)∈Wk},(116) b Hk:= max t∈Wk Hmax(z(t)),(117) b Ξk:= median {Ξ(z(t)) : S(t)∈Wk},(118) where Ak, Bk are the basins intersecting Wk with minimal and next–minimal deficit. For the graph corridor G= (V,E), replace time sampling by node sets {v:S(zv)∈Wk}, and compute: d ∆Vk:= min π∈Pk max v∈πΦ(zv)−min v:S(zv)∈Wk Φ(zv), 34 where Φis the node state cost (96) and Pkare paths crossing Wk. Collect the normalized indices Jk:= J1k,...,J5k∈[0,1]5,(119) by mapping (114) – (118) with the sigmoid transforms used in (113) (thresholds and slopes specified below). Null calibration and data–driven thresholds For each index i∈ { 1 ,..., 5 } we generate a null distribution that represents “no jump” within window Wk: 1. Continuous path: simulate B bootstrap paths z(b) ( t )by block–bootstrapping residuals of the fitted teleonomical flow with the same tube and reachability (preserving autocorrelation), compute the windowed index b T(b) ik (e.g., d ∆Vk), and form the empirical cumulative distribution b Fik.[211] 2. Graph path: resample local edge costs Φand c under the null by shuffling node labels within Wk or by parametric resimulation from the estimated state–cost model, recompute the windowed statistic, and build b Fik. Set per–index thresholds τik by τik := b F−1 ik (1 −αi),(120) for user–chosen marginal Type–I rates αi∈(0,1) (possibly index–specific). Convert b Tik to p–values pik := 1 −b Fikb Tik,(121) and define the binary flags Bik := 1{pik ≤αi}. False–discovery control across windows. Let pk := min{pik1, pjk2} over the best two indices in window Wk (enforcing the rule “at least two indices must be significant”). Apply the Benjamini–Hochberg (BH) procedure at level qto {pk}K k=1: p(1) ≤ ··· ≤ p(K),ˆ k:= max k:p(k)≤k Kq,reject p(1), . . . , p(ˆ k).(122) Theorem V.4 (FDR control for the operational detector) . If the window–wise pk are independent or satisfy PRDS (positive regression dependency on a subset), then the BH step–up rule (122) controls the false discovery rate at q [ 128 ]. In particular, the two–index minimum preserves validity because min of valid p –values is super–uniform under independence, and PRDS is inherited by monotone combinations [129]. Proof sketch. Under independence/PRDS, BH controls FDR at q . The mapping ( pi1, pj2 ) 7→ min ( pi1, pj2 ) is coordinatewise nonincreasing and preserves PRDS monotonicity; super–uniformity follows from P ( min ( U1, U2 ) ≤u )=1 − (1 −u ) 2≤ 2 u for U1, U2∼Unif (0 , 1) independent. Thus the BH assumptions hold for {pk}. Decision rule and reporting A window Wk is called a jump window if (i) at least two indices are significant ( Bik + Bjk ≥ 2for some i6=j) and (ii) the combined pksurvives BH at level q. The representative point of Wkis chosen as ˜zk:= arg max z:S(z)∈Wk JS(z)with JS as in (113),(123) breaking ties by larger progress S. We report: •the index vector Jk, the raw statistics (d ∆Vk,b λk,b ∆PH,k,b Hk,b Ξk); •the trait deltas ∆xk:= x(˜zk)−x(zprev)to the previous non–jump window; • for the graph corridor, the sequence edit set Ek associated with the edges traversed within Wk , summarized by frequency or by the minimum cardinality edit set consistent with ˜zk; •confidence bands from the bootstrap (e.g., empirical 95% intervals). 35 Algorithmic specification Algorithm 2 Operational jump detector (OJD) Require: Path samples {z ( t` ) } with progress S ( t` )(or graph nodes {zv} with S ( zv )); thresholds {αi} ; BH level q; window centers {sk}and half–width ∆S; number of null bootstraps B 1: for k= 1 to Kdo 2: Collect samples/nodes with S∈Wk, compute windowed statistics (114)–(118) 3: for i= 1 to 5do 4: Generate null b Fik by block bootstrap (continuous) or local resimulation (graph) 5: Compute pik ←1−b Fik(b Tik)and flag Bik ←1{pik ≤αi} 6: end for 7: Set pk←mini6=j{pik, pjk}over the two smallest indices; record Jkvia (119) 8: end for 9: Apply BH (122) to {pk}, obtain rejected set K 10: for k∈ K do 11: Choose representative ˜zkby (123); compute ∆xk; if graph: compute Ek 12: Output record: (k, sk,Jk,b Tik, pik, pk,∆xk, Ek) 13: end for Robustness and power Robust thresholds. To guard against heavy–tailed noise in indices (e.g., from occasional PH instability), use Huberization on the raw statistics before bootstrapping: Thub i:= Huberκ(Ti)with κ= MAD ·c, and propagate through (121) . Huber M–estimators provide minimax robust efficiency under contamination [135, 136]. Power analysis (approximate). If the null distribution of an index statistic Tik is approximately normal with mean µ0ik and variance σ2 0ik (often true for smoothed curvature and conductance), then the minimal detectable effect at level αiand power 1−βis MDEik ≈z1−αiσ0ik +z1−βσ1ik, where σ1ik is the alternative variance (estimated from pilot data) and zp is the p –quantile of N (0 , 1). This informs window size and bootstrap B. Complexity For K windows, B bootstraps, and per–window sample size nk , computing the five statistics is O(Pknk) plus the cost of PH/sheaf updates (dominant in L2 ). Null calibration adds O(BPknk) ; in the graph case, local resimulation is O ( B¯ d nk )with average out–degree ¯ d . All windows are embarrassingly parallel. Summary The operational detector (i) aggregates five categorically grounded indices, (ii) calibrates thresholds from a data–driven null with bootstrap/resimulation, (iii) enforces a two–index rule with BH control of FDR across windows, and (iv) reports interpretable jump candidates with trait deltas and, when available, sequence edits. These outputs serve both as alerts for surveillance and as hypotheses about the categorical tensions (L2/L3/L4) that precipitate discontinuities. 36 VI. EXPERIMENTS A. Setup Case study and goal. We evaluate the categorical teleonomical calculus on the Omicron lineage path BA.1−→ BA.2−→ BA.5−→ XBB, using the four calibrated traits x = ( Tbind, Tcleave, Tstab, Tescape )from §IV A. The task is trajectory prediction: starting from BA.1, recover a channel that (i) respects reachability and categorical invariants (L0–L5), (ii) progresses monotonically in the global coordinate S (§IV D), and (iii) aligns with the observed polyline across all pairwise projections. We report results both for the continuous proximal flow and for the discrete mutation–graph corridor (§IV E). Training signals and splits Supervision. The teleonomical energy L ( z ; λ ) = PiwiLi ( z ; λ )contains unknown level weights wi and, optionally, era targets τ ( λ )in L1 . We supervise w using pairwise orderings among variants (BA.2 ≺ BA.1, BA.5 ≺ BA.2, XBB ≺ BA.5) and, where available, directional constraints on progress S along the observed polyline. Data splits. Because the case study is short, we adopt two robust splits: 1. Leave–one–segment–out (LOSO). Train on two consecutive segments (BA.1 → BA.2 and BA.2→BA.5), validate on the held–out segment (BA.5→XBB); rotate the holdout. 2. Leave–one–variant–out (LVO). Train w and era targets using three variants and their orderings, hold out the fourth for evaluation of out–of–sample extrapolation. For each split we generate windowed augmentations by sampling ˜x along short chords of the observed polyline with Gaussian noise consistent with the calibrated trait covariances Σ (v) (§IV A). This provides enough local variation to fit wwithout violating reachability. Fitting wand era targets Weight learning. We solve the nonnegative RankSVM problem in (76) with a cross–validated penalty C (grid in { 10 −2, 10 −1, 1 , 10 , 10 2} ). When sparsity is desirable, we add an `1 penalty on w (Lasso) and tune its coefficient λ`1 by nested cross–validation [ 137 ]. In small– n regimes, we prefer ridge ( `2 ) regularization to stabilize weights [ 138 ]. Random search over ( C, λ`1 )augments the grid to mitigate discretization bias [139]. For each split, the final w?is the median across 5random restarts to reduce estimator variance. Era targets. Era targets τ ( λ )in L1 are estimated by a Kalman smoother on the time–stamped trait series (if available) or, absent time series, by solving a Tikhonov regression that matches the centerlines of the training segments while penalizing roughness (ridge parameter selected by cross–validation). We propagate uncertainty in τ by sampling from its posterior and refitting w ; reported channels marginalize over these draws. teleonomical solver hyperparameters Continuous proximal flow. Unless stated otherwise we use: β0= 5, βmax = 50, γβ= 1.08 (per iteration), λprog = 2.5, λcurv = 50. Armijo backtracking uses c1 = 10 −4 with initial step α0 = 0 . 2; trust caps on orthogonal tube distance are δk= 10−3(summable). Asymmetric per–trait step caps (§IV D) are ∆+=0.06 0.03 0.02 0.30,∆−=0.02 0.02 0.02 0.10(bind, cleave, stab, escape). We terminate when the relative improvement in the composite objective Φ k falls below 10 −5 for 20 consecutive iterations or when Sreaches the terminal segment. 37 Discrete mutation–graph corridor. Node set: BA.1, BA.2, BA.5, XBB plus a D = 3 neighborhood of DMS–admissible edits and a small set of recombination mosaics (§IV E). Beam width B = 24, maximum depth T = 15, state–cost weights ( αtube, α2, α0, α1 ) = (1 . 0 , 0 . 8 , 0 . 1 , 0 . 2), edge penalties ζ = 1 . 0, ξ∈ { 0 . 5 , 1 . 0 , 2 . 0 } (tuned by LOSO). Heuristic h is the admissible tube+structure bound in (99) with precomputed lower envelopes along segments. Progress monotonicity is enforced as a hard constraint. Sequence→trait map for graph nodes When needed (e.g., to assign traits to candidate haplotypes in the corridor), we train a DMS–guided predictor f:AL→R4. Two instantiations are supported: 1. Global epistasis linear link: f ( g ) = uβ>ϕ ( g )  with one–hot features ϕ ( g ), monotone link u , and ridge penalty on β(tuned by cross–validation). 2. Deep sequence (variational) model: a pretrained latent generative model fit on coronavirus spike families; the trait head is a small MLP trained on available DMS/assay labels with Adam (learning rate 10−3, batch size 128) [143, 144]. For both, calibration to the standardized trait scales of §IV A is enforced by affine post–processing fit on the training split only. We stress that in the reported experiments, we primarily evaluate channels in the calibrated trait space; sequence models are used only to populate the graph corridor with plausible nodes. Baselines To demonstrate the load–bearing role of structure (L2) and reachability, we compare against three baselines: (B1) Nearest–neighbor (NN) in trait space. At iteration k , choose the next point as the nearest observed anchor ahead in progress S: xNN k+1 = arg min x∈{BA.2,BA.5,XBB}:S(x)≥S(xk)kx−xkkΣ−1, with linear interpolation between anchors to match the number of steps. This baseline is consistent with the classical NN rule [ 140 ] and yields a piecewise–linear trajectory that ignores structural and reachability constraints. (B2) Unconstrained gradient step. Perform gradient descent on L1 only (no tube, no structure, no reachability): xUG k+1 =xk−α∇xL1(xk;λ), α ∈ Aα, with Aα = { 10 −3, 10 −2, 10 −1, 0 . 2 } tuned by split–wise validation. This tests whether era–directional pressure alone recovers the empirical manifold. (B3) Landscape hill–climb. Maximize a scalar “fitness” Fland ( x ) = b>x with nonnegative weights b fit by ridge regression to match the ordering BA.2 > BA.1 > . . . (learned on training segments). The climb uses xHC k+1 =xk+η b, η ∈ Aη, projected only onto hard trait bands (no tube or graph corridor). This provides a counterfactual in which structure and reachability are ignored but monotone improvement in a scalar objective is enforced. All baselines are tuned by the same LOSO/LVO procedures to avoid selection bias. Evaluation metrics We report four complementary metrics, computed on the full 4D trajectory after Euclidean reparameterization by arc length: 1. Curve discrepancy: discrete Fréchet distance between predicted and observed polylines (4D). 38 2. Integrated deviation: mean integrated squared error (MISE) between the two curves after Procrustes alignments of local segments [141]. 3. Endpoint error: Mahalanobis distance from the predicted terminal point to XBB. 4. Transport discrepancy: Earth mover’s distance between empirical measures induced by uniform sampling of the two curves [142]. We also track constraint violations: average orthogonal tube distance, fraction of steps exceeding caps, and the structural deficit L2 along the path. For the mutation–graph experiments we additionally report the edit length and a breakdown of substitutions vs. recombination edges. Cross–validation and uncertainty quantification We use nested cross–validation: outer LOSO/LVO for reporting and inner 3–fold CV on augmented windows for hyperparameters ( C, λ`1, ζ, ξ, λcurv, βmax )[ 145 ]. To guard against optimistic bias, we do not reuse validation windows for any tuning (split integrity). Uncertainty over channels is quantified by parametric bootstrap over ( w, τ )and, in the graph case, over the state/edge costs; bands are visualized as 95% pointwise envelopes. We also report selection stability of w across resamples (fraction of runs in which each level receives nonzero weight) to assess identifiability [146]. Implementation and reproducibility All continuous solvers are implemented in Python/NumPy with JIT–accelerated kernels for tube projections and sheaf residuals; graph searches use a binary heap priority queue. Random seeds are fixed to 42 unless stated. We release configuration files with the exact hyperparameters, grids, and split definitions to ensure bitwise reproducibility. Computations run in under a minute on a modern laptop (continuous flow) and under five minutes for beam width B≤32 (graph). The above setup isolates the contribution of categorical structure (L2) and reachability (graph corridor) against strong but structure–ignorant baselines, while using split strategies and uncertainty quantification appropriate for the small but high–signal Omicron case study. B. Metrics We assess predicted evolutionary channels against the observed Omicron polyline in R4 (traits x = (Tbind, Tcleave, Tstab, Tescape)). Let P={pi}m i=1 ⊂R4, Q ={qj}n j=1 ⊂R4 be sampled points on the predicted and observed curves after arc–length reparameterization and uniform resampling. Distances use the Mahalanobis metric induced by the global trait covariance Σ  0from §IV A: d(u, v) := ku−vkΣ−1=q(u−v)>Σ−1(u−v). 4–D Dynamic Time Warping (DTW) DTW aligns P and Q by a monotone warping path π =  ( i1, j1 ) ,..., ( iL, jL )  with 1 = i1≤ ··· ≤ iL = m and 1 = j1≤ ··· ≤ jL=n, allowing unit steps (1,0),(0,1),(1,1). The 4–D DTW cost is DTW(P, Q) := min π L X `=1 d(pi`, qj`).(124) The standard dynamic program with boundary/step constraints yields the recursion (initializing D (0 , 0) = 0and D(i, 0) = D(0, j)=+∞) D(i, j) = d(pi, qj) + min{D(i−1, j), D(i, j −1), D(i−1, j −1) },DTW(P, Q) = D(m, n).(125) 39 Lemma VI.1 (Basic properties of 4–D DTW) .DTW ( P, Q ) ≥ 0with equality iff P and Q are identical up to repetition; it is invariant under common rigid motions in the Σ −1 metric and computable in O ( mn ) time and O(mn)memory (or O(min{m, n})with Hirschberg–style pruning). For differentiable training (e.g., to tune hyperparameters by gradient methods), we also report the soft DTW sDTWγ (temperature γ > 0), which smooths the min in (125) via LSEγ and is everywhere differentiable in the inputs [147]: Dγ(i, j) = d(pi, qj) + LSEγ{Dγ(i−1, j), Dγ(i, j −1), Dγ(i−1, j −1)},sDTWγ=Dγ(m, n).(126) 4–D Hausdorff distance The (undirected) Hausdorff distance between finite sets P, Q under dis H(P, Q) := max n−→ H(P, Q),−→ H(Q, P)o,−→ H(P, Q) := max p∈Pmin q∈Qd(p, q).(127) H is a metric on nonempty compact subsets of R4 ; it upper–bounds the maximum pointwise deviation after best nearest–neighbor correspondence and is sensitive to outliers. We also report the modified Hausdorff average, which is more robust in practice [148]: −→ Havg(P, Q) := 1 |P|X p∈P min q∈Qd(p, q), Havg(P, Q) := max{−→ Havg(P, Q),−→ Havg(Q, P)}.(128) Directional mean absolute error (dMAE) per trait Let the observed polyline be segmented by eras E = {e} (BA.1 → BA.2, BA.2 → BA.5, BA.5 → XBB). For trait j∈ {bind,cleave,stab,escape}, define observed and predicted deltas over era e: ∆obs j(e) = xobs j(eend)−xobs j(estart),∆pred j(e) = xpred j(eend)−xpred j(estart). The directional MAE penalizes both magnitude and sign errors: dMAEj:= 1 |E|X e∈E∆pred j(e)−∆obs j(e)+λsgn 1nsgn ∆pred j(e)6= sgn ∆obs j(e)o|∆pred j(e)|+|∆obs j(e)|, (129) with λsgn ∈ [0 , 1]. When all signs match, dMAEj reduces to standard MAE; for mismatched signs the additional penalty reflects qualitative disagreement about the direction of change by era. We also report the sign accuracy SignAccj:= 1 |E|X e∈E 1nsgn ∆pred j(e) = sgn ∆obs j(e)o.(130) MAE–type metrics are scale–interpretable and robust compared to squared errors [150]. Path smoothness and reachability diagnostics Discrete curvature and roughness. Let {rk}N k=1 be a uniform resampling of the predicted path. The (Mahalanobis) discrete second difference measures curvature: R2:= 1 N−2 N−1 X k=2 krk+1 −2rk+rk−1k2 Σ−1.(131) Optionally, a third–difference “jerk” term R3 = 1 N−3PN−1 k=3 krk+1 − 3 rk + 3 rk−1−rk−2k2 Σ−1 can be added to penalize oscillations. In the continuous limit, R2 approximates Rkκ ( s ) k2ds , the squared curvature integral of the curve s7→ r(s)[149]. 40 Tube adherence. Orthogonal distance to the channel is measured by the soft tube energy Ltube,β (§IV D). We report d⊥:= 1 N N X k=1 min i 1 √2krk−Πi(rk)kΣ−1, dmax ⊥:= max kmin i 1 √2krk−Πi(rk)kΣ−1,(132) where Π i is the Mahalanobis projection onto segment Γ i and the 1 /√2 converts squared distances in di to Euclidean units. Step caps and graph feasibility. Let per–step trait changes be ∆ rk = rk+1 −rk . Given asymmetric caps [−∆−,∆+]⊂R4, the violation rate and excess are ViolRate := 1 N−1 N−1 X k=1 1{∆rk/∈[−∆−,∆+]},Excess := 1 N−1 N−1 X k=1 [∆rk−clip(∆rk)] Σ−1, (133) with clip the componentwise projection into the cap box. For the mutation–graph corridor, define the feasible–edge rate EdgeFeas := #{(u→v)on path : (u→v)∈ E} path length ,(134) and the edit budget statistics (median and max number of substitutions per step; fraction of recombination edges). Structural coherence along the path. We summarize structural consistency by the running average and maximum of L2: L2:= 1 N N X k=1 L2(rk), Lmax 2:= max kL2(rk),(135) where L2 = LPH 2 + Lglue 2 monitors contact–graph topology and sheaf–gluing residuals (§IV C). Low L2 and Lmax 2indicate that the channel respects structural constraints throughout. Aggregation and confidence All scalar metrics are reported with 95% bootstrap confidence intervals (resampling path segments with replacement and respecting temporal order). For multi–metric comparison we compute z–scores (within split) and a composite score by robust averaging (Huber location) to prevent domination by any single metric. Together, DTW captures aligned deviations, Hausdorff bounds worst–case offsets, dMAE verifies directional correctness per trait and era, and smoothness/reachability diagnostics ensure biophysical plausibility and corridor adherence—exactly the desiderata of a teleonomical, reachability–aware predictor. C. Results Train/test splitting and tube construction. Let the era anchors be A = {BA.1,BA.2,BA.5,XBB} , ordered in time, and let E = { ( BA.1→BA.2 ) , ( BA.2→BA.5 ) , ( BA.5→XBB ) } denote the three evaluation windows. For any held-out window e?∈ E , we construct the polyline only from the training anchors and exclude e?: Γtrain(e?) := [ e∈E\{e?} Γ(e), and define the soft tube as the soft minimum of Mahalanobis orthogonal distances to the segments of Γtrain(e?), Ltrain tube,β(z|e?) := −1 βlog X [a,b]⊂Γtrain(e?) exp−β1 2z−Π[a,b](z)2 Σ−1,(136) where Π [a,b] is the Mahalanobis projection onto the segment [ a, b ]and kuk2 Σ−1 = u> Σ −1u . Thus, the test window never contributes to the tube potential used for guidance or evaluation. 41 Figure 1: Observed vs. predicted channels (all pairwise projections). Graph-corridor variant after L 2 stabilization and tube/smoothness trust. Axes use the calibrated, anisotropic trait scales defined in IV A; the tube is the soft–min potential of IV D. Split policies. We report two standard protocols: (i) Forward-chaining (rolling origin) for the last window BA.5→XBB , where Γ train uses only BA.1→BA.2 and BA.2→BA.5 ; and (ii) Leave-one-segmentout (LOSO), where each window is held out in turn, the tube is rebuilt via (136) , and metrics are averaged over the three runs. Hyperparameter tuning. All tube/solver hyperparameters (e.g., β schedule, trust weights, step caps) are tuned within the training split only via inner cross-validation on augmented training windows. No information from the held-out window is used in tuning or tube construction. Qualitative alignment across trait projections. Predicted channels (continuous proximal flow and graph corridor) closely follow the empirical BA.1 → BA.2 → BA.5 → XBB manifold in all pairwise projections of the four–trait space x = ( Tbind, Tcleave, Tstab, Tescape ), without spikes or lateral excursions. Figure 1 consolidates the four 2-D overlays into a single artifact-free composite after L 2 stabilization and tube/smoothness trust, as specified in IV D and IV E. The corridor constraint and per-step caps ensure that each move is reachable in the mutation graph (cf. IV E). Visually, the model reproduces the observed increase in cleavage with steady binding (BA.2 → BA.5) and the trade between escape and binding near the BA.5→XBB transition, while staying inside the soft tube. Quantitative agreement and reachability. Table V reports the main accuracy and plausibility metrics (4-D sDTW, modified 4-D Hausdorff Havg , mean orthogonal tube distance d⊥ , discrete curvature R2 ), together with directional sign agreement and edge feasibility (fraction of moves that are legal graph edges). All metrics use the Mahalanobis geometry and path resampling described in IV A; confidence intervals are bootstrap-based. Directional deltas by era. Era-wise trait changes are captured in Table VI via the directional MAE ( dMAE ) and sign agreement (penalty λsgn = 0 . 5). The predicted deltas match the observed direction in all traits across all three eras; magnitudes are small in the calibrated units, consistent with the tube-tangent progression. 48 proposes ∆x /∈ A, then the graph–constrained proximal step solves min ∆x0∈A 1 2k∆x0−∆xk2 Σ−1+τΦ(x+ ∆x0), and achieves strictly smaller composite objective than the infeasible proposal (interpreted with infinite penalty outside reach), hence is teleonomically superior. Moreover, if the corridor refines and densifies so that A → the true instantaneous reachable cone, the discrete path cost converges to the continuous action. Proof. Immediate from convex projection onto A followed by a standard descent inequality; the convergence claim follows by Γ–convergence of discrete actions to the continuous functional under mesh refinement [162]. Biological meaning. The graph corridor encodes mutation supply, codon constraints, protease site integrity, glycan occupancy, and known recombination breakpoints; forbidding edges that break these rules ensures biophysically legal motion. In effect this is Wright’s dictum that evolution must proceed by “small steps” but made precise by admissible edit sets and trait caps tied to DMS envelopes and routine feasibility. Relation to and strengthening of classical views Adaptive landscapes. Fisher and Wright placed selection on a scalar fitness surface [ 160 , 161 ]. Our L generalizes this to a multi–level energy where formal/structural and order constraints (L2/L4) are as causally consequential as immediate trait gains; the channel is therefore not merely the intersection of fitness ridges but the minimizer set of cross–level invariants within the reachable cone. Canalization. Waddington’s canalization becomes computable via the soft tube Ltube,β ; Proposition VII.1 and Theorem VII.2 show that once canalized, the motion is tangent and narrow. Optimal processes. The proximal step (30) is the discrete analogue of a constrained optimal control step, and the continuous inclusion is a maximum–principle flow on the reachable set [ 163 ]. The Lagrange–KKT conditions determine which levels “bind” at any time, explaining when structure (L2) or order (L4) dominate the direction. Takeaways •Teleonomy ⇒channels. Because the normal component of the step is penalized by the tube and trust, while the tangent component is fed by progress and projected cross–level gradients, the motion stays in narrow corridors and advances monotonically. •L2 and L4 ⇒controlled nonlinearity. Structural coherence rotates the descent and creates bends where structural gains align with the tangent; order effects create additional curvature via commutators, and when combined with barrier softening, enable principled jumps. •Graph constraint ⇒physical plausibility. Discrete reachability forbids unphysical trait changes and ensures that every predicted step is implementable by real edits or routine shifts. These points explain why the method produces smooth, artifact–free tracks that also anticipate plausible discontinuities, and why ablations that remove L2 or the corridor degrade both plausibility and quantitative fit. B. Limitations Our experiments demonstrate that a teleonomical, reachability–aware calculus can reproduce observed SARS–CoV–2 trait trajectories and yield interpretable jump signals. Nevertheless, several modeling choices are provisional or conservative. We detail these limitations, give precise mathematical statements of the induced approximation errors where possible, and outline the additional measurements required to remove remaining placeholders. 49 (L2) Structural proxies versus ground truth In this work, the structural level L2 is instantiated by (i) persistent homology of residue–contact graphs and (ii) sheaf gluing residuals over a functional patch cover. These are proxies for the true structural deficit Ltrue 2 that a full biophysical model (including conformational ensembles and glycan dynamics) would deliver. Let Lprox 2 denote our proxy and assume a uniform approximation on a compact trait–structure domain D: sup z∈D Lprox 2(z)−Ltrue 2(z)≤ε. (137) Let Ltrue =PiwiLiwith L2=Ltrue 2and Lprox the same but with Lprox 2. Proposition VII.6 (Stability of minimizers under uniform surrogate error) . Assume Ltrue is µ –strongly convex orthogonally to the tube Γand has L –Lipschitz gradient along Γ. Then any teleonomical proximal step (or flow) computed with Lprox yields states zprox whose orthogonal distance to the Ltrue solution ztrue satisfies dist⊥zprox, ztrue≤r2ε µand kΠTΓzprox −ztruek ≤ ε L+O(ε2), provided both iterates remain in D. Proof. Strong convexity in the normal directions gives the standard error bound ku⊥k2≤ 2( Lprox − Ltrue ) /µ ≤ 2 ε/µ . Along the tube, a Lipschitz gradient implies that an ε –perturbation of the objective shifts first–order optimality by at most ε , hence a displacement of size ε/L to restore stationarity (implicit function argument). Implication. Provided the proxy error ε is small, the channel computed with Lprox 2 differs from the (unknown) ground–truth channel by O ( √ε )orthogonally and O ( ε )tangentially. However, when (137) fails (e.g., missing a conformational state or glycan reorientation), the bound does not apply and the bend induced by L2 may be misplaced. A production system should incorporate ensemble–aware L2 terms, potentially via cryo–EM ensemble weights or MD free–energy profiles. From an optimization perspective, the teleonomical step depends continuously on such smooth perturbations under standard perturbation theory [164]. (L3/L4) Placeholders and identifiability Routine kernel P (L3). We used a placeholder routine kernel with generic throughput targets. In practice, P must be estimated from panel experiments (time–lapse or perturbation counts). If b P is the smoothed MLE from counts (row–wise Dirichlet–multinomial), then for each row u, kb Pu•−P? u•k2.slog(1/δ) Nu w.h.p. (1 −δ), by vector Bernstein for multinomial counts. Consequently, the plug–in throughput b Θ deviates from Θ( P? )by Omaxuplog(1/δ)/Nu . Estimation of the spectral gap (conductance proxy) inherits matrix concentration rates Oqlog d N for d = | S | [ 165 ]. Small Nu yields high variance in L3 and weak jump signals; identifiability requires sufficient visits to bottleneck transitions. Holonomy (L4). We treated L4 as a functional that can be computed given order experiments. Identifiability of noncommutativity requires a crossover design: both ab and ba arms with enough replicates to beat assay noise. For mean difference ∆ a,b and per–arm variance σ2 , the minimal detectable effect in a balanced 2 × 2crossover is MDE ≈z1−α/2σp2/n per arm [ 166 ]. Absent such data, L4 remains a prior or a proxy and jump signals that rely on L4 spikes may be underpowered. 50 Simplified environment and partial observability We modeled the environment λ as piecewise constant (eras). Real environments drift and are only partially observable (e.g., immune landscape, contact patterns). A principled formulation is a partially observable Markov decision process (POMDP) with hidden state et and observations ot ; the teleonomical flow would then be conditioned on a belief state bt = P ( et|o≤t )[ 167 ]. In the absence of a full POMDP, our schedule for τ(λ)(§VI A) is a crude filter. Proposition VII.7 (Online drift and regret) . Let {Lt} be a sequence of convex objectives with variation budget VT = PT t=1 supz|Lt+1 ( z ) −Lt ( z ) | . A projected mirror–descent teleonomical step with step sizes ηt∝1/√tachieves dynamic regret T X t=1 Lt(zt)−Lt(z? t)=O√T+VT, where z? tis the instantaneous minimizer on the reachable set. Proof. This is a standard bound in online convex optimization with shifting comparators; the projection onto the reachable set preserves the analysis [168]. Implication. If the environment drifts slowly ( VT small), the projected teleonomical updates track the moving target with sublinear dynamic regret; rapid regime shifts violate the bound and require explicit change–point handling. Data heterogeneity and calibration Trait data come from heterogeneous labs, assays, and panels. While §IV A used hierarchical calibration, residual batch effects can persist. Empirical Bayes harmonization (e.g., ComBat) can reduce location/scale batch effects across labs [ 169 ], but must be adapted to our non–Gaussian traits (logit cleavage, log fold–escape). Between–study heterogeneity can be summarized by I2 statistics in random–effects meta– analysis [ 170 ]; when I2 is high, pooling to form Σ(trait anisotropy) inflates uncertainty and weakens tube penalties. Accurate cross–study calibration also requires principled measurement error modeling; classical errors bias gradients and can misalign bends [ 171 ]. Missing traits for some variants should be imputed with multiple imputation to propagate uncertainty into L[172]. Minimal mutation graph Our graph corridor was deliberately small: a shallow neighborhood of observed haplotypes with a few mosaics. This minimizes false positives but risks false negatives: reachable edits not represented in G make a true canal inaccessible to the search. A production system should: • Expand nodes by DMS–admissible edits (per position, per trait) to depth D guided by global epistasis priors [173, 174]. • Add phylogeny–consistent edges (or mosaic nodes) constrained by real breakpoints inferred on large trees (e.g., UShER) and track lineage context [175, 176]. • Maintain admissible step caps informed by sequence– → trait predictors and by assay envelopes; caps that are too tight force underexploration; too loose allow unphysical steps. Search width and pruning. Finite beam width B can prune optimal branches. While §IV E provides a data–dependent suboptimality bound in terms of pruning gaps, in practice one should monitor frontier diversity and adapt B when f –values cluster, or use beam–stack/backtracking variants to recover pruned paths. 51 Small–nsupervision and weight identifiability Weights wi are fit from a few order constraints (BA.2 ≺ BA.1, etc.). This is statistically fragile. For the nonnegative RankSVM with feature vectors φ ( z )=( L0, . . . , L5 )bounded as kφk ≤ R , classical margin–based generalization theory yields (with probability 1−δ) for any new pair (u, v): Pr misorder≤O kwk2R √N+rlog(1/δ) N!, where N is the number of training pairs and the hidden constant depends on the margin [ 177 ]. With very small N , regularization (ridge) and prior structure (e.g., grouping structural levels) are essential; still, some trade–offs between levels may remain underdetermined. Summary and priorities for future work 1. Replace L2 proxies with ensemble–aware structural terms; quantify ε in (137) using blinded benchmarks; propagate uncertainty into L. 2. Instrument L3 with routine panel experiments sufficient to estimate P (and its gap) with matrix– concentration–scale confidence; design balanced ab/ba crossovers to identify L4 holonomy with adequate power. 3. Upgrade the environment model from piecewise constants to a filtered, partially observed process; use online teleonomical updates with change–point detection under explicit regret guarantees. 4. Deploy a DMS– and phylogeny–informed corridor with adaptive beam search; monitor frontier diversity and pruning gaps to keep suboptimality under control. 5. Standardize trait calibration across labs with harmonization and explicit error models; use multiple imputation for missing entries; report I2to quantify residual heterogeneity. These steps move the system from a proof–of–concept with conservative proxies and placeholders to a production–quality, data–rich pipeline with identified L3/L4 terms and a graph corridor that saturates the reachable cone. C. Broader implications The teleonomical calculus offers two immediate avenues of impact: surveillance, by ranking likely futures and issuing early warnings for discontinuities, and generality, by porting the same categorical structure (levels, causes, reachability) to other biological and engineered systems that possess invariants and routines. We set out the mathematics that underwrites both uses and indicate concrete paths to deployment. Surveillance: risk ranking and early warning Fix a horizon T > 0. For a state z lying in a basin A and a reachable alternative basin B , the barrier height ∆ V ( A→B )was defined in (102) . In small–noise metastable regimes, exit times obey Kramers–Eyring asymptotics (§V A): E [ τA→B ] κ−1exp (∆ V/θ )for effective noise scale θ > 0. Define the jump probability proxy in horizon Tby PT(z) := 1 −exp−T/E[τA→B(z)]≈1−exp−κT e−∆V(z)/θ.(138) Proposition VII.8 (Monotone equivalence of rankings) . For fixed ( κ, T, θ ), the map z7→ P T ( z )is strictly decreasing in ∆ V ( z )and strictly increasing in − ∆ V ( z ). Hence ranking states by P T is equivalent to ranking them by −∆V. Proof. ∂PT/∂∆V=−(κT/θ) exp(−∆V/θ) exp(−κTe−∆V/θ)<0. 52 Thus a barrier score S1 ( z ) := − ∆ V ( z )provides a calibrated ordering of near–term jump risk. To fuse additional, categorically grounded signals, define a risk vector T ( z )=( T1, . . . , T5 )with T1 = − ∆ V , T2 = −λreach min (curvature softening), T3 = ∆ PH + ∆ glue (L2 flip), T4 = Hmax (holonomy), T5 = Ξ (routine split), and let X(z) = PiβiTibe a scalar index field with nonnegative weights βi. Coherent risk aggregation. To turn X ( z )into a surveillance risk, employ a coherent risk measure ρ (monotone, translation invariant, subadditive, positively homogeneous) [ 178 ]. A standard choice is the conditional value–at–risk (CVaR) at level α∈(0,1): ρα(X) := inf η∈Rη+1 1−αE(X−η)+,(139) which emphasizes tail configurations [ 179 ]. With uncertainty injected through bootstrap ensembles of X (from data/parameter posteriors), the surveillance risk is Rα ( z ) := ρα ( X ( z )) and an alarm is raised when Rα(z)≥τα. Theorem VII.9 (Order preservation and subadditivity of warnings) . If X1 first–order stochastically dominates X2 (i.e., X1 yields smaller jump indices), then Rα ( X1 ) ≤Rα ( X2 )for all α . Moreover, for independent regions D1,D2with risks R(1) α, R(2) α, the aggregate risk is subadditive: R(1∪2) α≤R(1) α+R(2) α. Proof. CVaR is monotone and subadditive on integrable random variables; the first statement follows from monotonicity under first–order dominance and the second from subadditivity and positive homogeneity [178, 179]. Integration with classical surveillance. The environment enters via λ ( t )(era), so Rα ( z, t )can be combined with nowcasts of transmission (e.g., Rt estimates from renewal models) to prioritize both plausible jumps and epidemiologically dangerous contexts [ 180 ]. For antigenically drifting pathogens such as influenza, coupling L2 flips to antigenic distances (e.g., cartography) focuses attention on genuinely immune–relevant changes [ 181 ]. At the population level, phylodynamic signals (growth advantages, clade turnover) play the role of L0 and can be added as a sixth index [ 182 ]. Early–warning concepts such as critical slowing down (rising variance/autocorrelation before transitions) can be attached to our curvature–softening index to sharpen timing [183, 184]. From viruses to organisms and engineered systems The categorical construction is domain–agnostic: objects encode material and formal structure; morphisms encode routines; invariants (final causes) define deficits; and reachable futures constrain motion. We outline three mappings. Antigenically drifting RNA viruses. Traits: receptor binding, cleavage/entry route, stability, and immune escape (as here). L2: HA (influenza), E (dengue), F (RSV) structural sheaves/topology; L3: life–cycle kernels augmented by host range parameters; L4: order of vaccine/drug/immune exposures (e.g., original antigenic sin). Risk: fuse barrier/curvature with antigenic distances and growth nowcasts. Bacteria and metabolism. Traits: growth rate, antibiotic MIC profiles, adhesion/biofilm proxies. L2: metabolic and regulatory structure; a natural choice is to penalize deviations from consistent flux balances (stoichiometric sheaves) under resource constraints; flux balance analysis supplies the formal constraints [ 185 ]. L3: cell–cycle and stress routines; L4: order of drug exposures (collateral sensitivity vs. cross–resistance). The calculus prioritizes trajectories that stay within feasible flux polytopes while moving toward stress–tolerant regimes. Engineered cyber–physical systems (CPS). Traits: performance envelopes; L2: interconnection invariants and resource bounds; L3: operational routines/modes; L4: order effects from mode switches. Final causes are safety and specification invariants; control–theoretic barrier certificates and invariance principles provide L0 / L2 analogues [ 186 – 188 ]. The teleonomical proximal step becomes a projected controller update that minimizes a multi–level cost while respecting the safe reachable set. At the category level, compositional frameworks (e.g., open Markov processes, wiring diagrams) make invariants modular, enabling level–wise assembly of large systems [189, 190]. Compositionality and transfer Let S be a category of subsystems with monoidal product ⊗ (parallel composition) and composition ◦ (series). Suppose each subsystem A∈ S carries deficits LA i and a reachable set ReachA . Define the 53 network by ( A⊗B )with LA⊗B i = LA i + LB i and ReachA⊗B = ReachA×ReachB (closed under the product). Proposition VII.10 (Compositional Lyapunov) . If LA = PiwiLA i and LB = PiwiLB i are Lyapunov– like (nonincreasing) along teleonomical steps on their reachable sets, then LA⊗B = Piwi ( LA i + LB i )is Lyapunov–like along product steps on ReachA⊗B. Proof. teleonomical steps minimize Dψ + ε ∆ tL on the reachable set. On the product, the Bregman term and constraints separate, so the sum of nonincreasing Lyapunov terms is nonincreasing. This property enables transfer: learned weights wi or learned L2 structural terms in one subsystem can be ported to larger assemblies without rederiving stability from scratch. Operational blueprint A public health instantiation requires: 1. Data feeds. Continuous ingestion of trait proxies (binding, cleavage, escape) from lab consortia; structural summaries (contact graphs, patch features); routine panel data where available; phylodynamic summaries for L0. 2. Calibration & assimilation. Hierarchical calibration of traits (§IV A); online filtering of era targets τ(λ); posterior ensembles for index uncertainty. 3. Corridor & search. A DMS– and phylogeny–informed corridor; near–real–time beam search (§IV E); computation of ∆V,λreach min , L2 flips, holonomy (where measured), and routine splits. 4. Risk fusion. CVaR aggregation Rα with user–set α ; alarm thresholds set by desired false discovery control (cf. §V B). 5. Prioritization. Rank variants/regions by Rα and Rt to coordinate experiments (neutralization panels, protease assays) and surveillance sequencing. Scientific payoffs •Mechanistic forecasts. Predictions come with mechanisms: which invariant is restoring, which routine is rerouting, which order effect is spiking. •Interpretability and design. In engineered settings, the same calculus designs routines and interconnections that keep invariants true; barrier certificates are “final causes” in control form [187, 188]. •Generalizability. Because the levels are categorical, not organism–specific, extension to new systems is principled: only the assembly functors Ai and the invariant functionals Ii need updating. Limits of generality Domain transfer presumes that (i) level invariants can be articulated; (ii) a surrogate for L2 exists with usable stability; (iii) reachability approximations are credible. Where these fail (e.g., poorly measured routines, unknown order effects), jump warnings may lose power. The cure is measurement design guided by the calculus: target those levels that dominate risk in Rαto close identifiability gaps. In summary, the calculus supplies a unified, compositional language to rank futures and warn of discontinuities, with proofs of order preservation and stability that make the results sturdy enough for practice. Its categorical bones let it travel: from SARS–CoV–2 to influenza and dengue; from cell metabolism to biofilm evolution; from viral evolution to safety–critical CPS—where the “final cause” is a safety invariant and the corridor is the safe reachable set. 54 VIII. CONCLUSION We presented a categorical teleonomical framework for evolutionary prediction, together with a practical calculus and algorithms, and validated it on SARS–CoV–2 Omicron lineages (BA.1 → BA.2 → BA.5 → XBB). The central thesis is that empirical evolution follows narrow channels because systems minimize a multi–level coherence deficit L(z;λ) = 5 X i=0 wiLi(z;λ), subject to reachability—the set of futures actually attainable by mutation, recombination, and routine deformations. The calculus: (i) formalizes teleonomy via a proximal step or differential inclusion constrained to the reachable set (§III B); (ii) builds a soft tube along observed channels with a global progress coordinate and trust, ensuring spike–free motion (§IVD); (iii) instantiates categorical levels L0–L5, with L2 (structure) and L4 (order/holonomy) providing the essential nonlinearity for bends (§IV C); (iv) enforces discrete reachability with a mutation–graph corridor and beam search (§IV E); and (v) supplies jump indices with an operational detector controlling false discoveries (§V B). We conclude by consolidating the main guarantees and their biological meaning. Global guarantee: channel adherence, descent, and convergence Let Γ = Si Γ i be the empirical corridor in trait space, Ltube,β the soft tube, and Reach ( z )the reachable set. Consider the composite objective Φk(z) = Ltube,β(k)(z) + Lprog(z;Star(k)) + 5 X i=0 wiLi(z;λ(k)) + Lsmooth(·), and the projected proximal step of (30) with Armijo–trust acceptance (90)–(91). Theorem VIII.1 (Channel adherence and convergence) . Assume: (a) each Li ( · ; λ )is proper, lower semicontinuous and convex in a neighborhood of Γ; (b) Ltube,β is strongly convex in the normal directions to Γ; (c) Reach ( z )is closed, convex, nonempty with tangent cones varying upper semicontinuously; (d) {δk}in (91) is summable and β(k)%βmax <∞. Then the accepted iterates {zk}satisfy: 1. Tube adherence: the orthogonal tube distance is bounded and cannot increase by more than δk per step; in particular, zkremains in a tubular neighborhood of Γwhose radius is O(δ•)(§IV D). 2. Descent: Φk(zk+1)≤Φk(zk)and Pkkzk+1 −zkk2<∞(finite length). 3. Limit set: any accumulation point z? is a constrained critical point of Φ( · )on the instantaneous reachable set: 0∈∂Φ(z?) + NReach(z?)(z?). 4. Uniqueness under curvature: if the restriction of Φto Γ ∩Reach is strictly convex along the terminal segment (e.g., a unique endpoint minimizer at XBB), then zk→z∞ and the trait path converges to the terminal anchor. Proof sketch. Items (1)–(2) follow from strong convexity of Ltube,β in normal directions, the tube trust cap, and Armijo descent; see Theorem IV.9 and Cor. IV.10. Item (3) is standard for Bregman proximal methods with closed convex constraints by passing the first–order optimality conditions to the limit (Painlevé–Kuratowski convergence of normal cones). Item (4) applies strict convexity along the admissible tangent directions to rule out multiple limit points. Meaning. The teleonomical step is a mathematically controlled canalization: motion remains glued to the corridor, energy decreases, and—absent degeneracies—the trajectory reaches the empirically correct terminus. 55 Jumps: principled indices with error control Let ∆ V be the barrier height between local basins constrained by reachability (102) , λreach min the softest restricted curvature, and (∆ PH, ∆ glue, Hmax, Ξ) the structural/holonomy/routine indices (§V A). The operational detector segments the path, calibrates nulls by bootstrap/resimulation, and applies a two–index rule with BH control across windows (§V B). Theorem VIII.2 (Soundness of alarms) . Under correct null calibration and independence/PRDS across windows, the detector’s Benjamini–Hochberg step–up rule at level q controls the FDR at q ; in addition, ranking windows by the barrier proxy − ∆ V is equivalent to ranking by the near–term jump probability proxy PTfor fixed (T, θ, κ)(§VII C). Proof sketch. FDR control is Theorem V.4. Monotone equivalence of rankings follows from Proposition VII.8. Meaning. Jump warnings are not heuristic: they carry a statistical guarantee on false discoveries and a mechanistic decomposition (barrier, curvature softening, topological/holonomy/routine changes). Empirical validation and biological takeaways On the Omicron case study, both continuous and graph variants accurately track the observed BA.1 → BA.2 → BA.5 → XBB manifold across four trait projections, with low DTW/Hausdorff discrepancies, high directional accuracy by era, and no unphysical spikes (§VI C). Ablations confirm that L2 (structure) and the graph corridor (reachability) are decisive. Detected jump windows align with biologically interpretable transitions (barrier drops, L2 flips, or routine re–routes), and the corridor constrains proposals to biophysically legal edits. Outlook The calculus is compositional and portable. It extends to pathogens with antigenic drift (influenza, RSV, dengue), to bacteria where flux and regulatory invariants dominate, and to engineered systems where safety/specification play the role of final causes. Immediate priorities are to (i) replace L2 proxies with ensemble–aware structural terms, (ii) instrument L3/L4 with panel and order experiments, (iii) enrich the corridor with DMS and phylogeny, and (iv) stand up a surveillance pipeline that ranks futures and issues early warnings with controlled error. Conclusion. Evolutionary change is not a random walk nor a blind hill–climb. It is a teleonomical movement in a constrained, categorical geometry: a minimization of multi–level deficits over reachable futures. The framework here makes that intuition computable and testable; on SARS–CoV–2 it yields channel forecasts and principled jump predictions in agreement with data. We believe the same categorical bones will support a broad class of biological and engineered systems in which invariants and routines govern change. Appendix A: Appendix / Methods (expanded details) This appendix collects formal definitions, complete derivations, and expanded methodological details. We use the same notation as in the main text. New symbols and assumptions are introduced locally and cross–referenced. Appendix B: Categorical preliminaries and L0–L5 as functorial constructions 1. Categories, functors, naturality A (small) category C consists of a class of objects Ob ( C ), a set of morphisms HomC ( A, B )for each pair of objects A, B , an associative composition ◦ : Hom ( B, C ) ×Hom ( A, B ) →Hom ( A, C ), and identities idA∈Hom ( A, A )with idB◦f = f = f◦idA for all f∈Hom ( A, B )[ 194 , 195 ]. A (covariant) functor 56 F : C → D assigns to each object A an object F ( A )and to each morphism f : A→B a morphism F ( f ) : F ( A ) →F ( B )such that F ( idA ) = idF(A) and F ( g◦f ) = F ( g ) ◦F ( f ). A natural transformation η : F⇒G between functors F, G : C → D assigns to each A a morphism ηA : F ( A ) →G ( A )so that for all f:A→Bthe square F(A)F(B) G(A)G(B) F(f) ηAηB G(f) commutes. 2. teleonomical category of states and levels Let Z be the category whose objects are states z = ( x, s ), with x∈R4 the trait vector, s an associated structural–routine datum, and whose morphisms are reachable updates m : z→z0 (single admissible edits and routine deformations). Composition concatenates updates; identities are null updates. The reachable set at zis the image of HomZ(z, −). Levels as functors. Define six functors Ai : Z → Ci and invariant functionals Ii : Ci→R≥0 ; Li:= `i◦Ii◦Aiwith `ia nondecreasing gauge (e.g., square). Concretely: •L0 (viability). A0 maps z to a next–generation operator K ( z, λ )(e.g., renewal kernel); I0 ( K ) := [1 −ρ(K)]+with ρspectral radius. •L1 (trait and bands). A1 is the inclusion of traits x ; I1 ( x ) := kx−τ ( λ ) k2 Σ−1 1 + Pjκ±[band violations]2 +. •L2 (structure). A2 returns a filtered contact graph and a patch sheaf; I2 aggregates topological and sheaf–gluing inconsistencies (§D). •L3 (routine). A3 maps to a kernel P ( z )on routine states; I3 penalizes low throughput and poor conductance. •L4 (holonomy/order). A4 maps to an environment action E: E → End ( Z ); I4 measures path–order asymmetry magnitude. •L5 (coherence). A5maps zto the vector (˜ L0,...,˜ L4);I5penalizes variance/imbalance. Proposition B.1 (Functoriality and nonincreasing invariants) . If each Ai is functorial and each Ii is natural with respect to Ai -morphisms that improve the invariant (e.g., stabilize a structure, increase throughput), then along any morphism m : z→z0 in Z we have Li ( z0 ) ≤Li ( z )for those i whose morphism action is improving; otherwise Liis nonincreasing up to noise. Proof. By naturality of Ii , Ii◦Ai ( m ) = η◦Ii for a morphism η in the target category corresponding to “improve”; monotonicity of `iyields Li(z0)≤Li(z). Appendix C: Variational derivation of the proximal step; trust conditions and descent 1. Bregman geometry, Moreau envelope, and proximal maps Let ψ : Z → R be C2 , µ –strongly convex. The Bregman divergence is Dψ ( ykz ) = ψ ( y ) −ψ ( z ) − h∇ψ(z), y −zi. For a proper f:Z → R∪{+∞}, the Bregman–Moreau envelope and proximal map are eλf(z) = inf y{f(y) + λ−1Dψ(ykz)},proxψ λf (z) = arg min y{f(y) + λ−1Dψ(ykz)}. If f is convex, eλf is smooth with ∇eλf ( z ) = λ−1 ( ∇ψ ( z ) − ∇ψ ( y? )) for y? = proxψ λf ( z ), and the optimality condition is ∇ψ(y?)−∇ψ(z) + λ ∂f(y?)30[191, 192]. 57 2. Teleonomical step with constraints and soft tube Define the composite f ( z ) = L ( z ; λ ) + ιC ( z )with C = Reach∆t ( zk ) ∩{step caps} (fixed during step k ), and add the soft tube and smoothness terms g(z) = Ltube,β(z) + Lprog(z;Star) + Lsmooth(·), which are smooth (for finite β). Consider the proximal–gradient subproblem at zk: zk+1 ∈arg min z∈C nh∇g(zk), z −zki+1 2αkkz−zkk2 ψ,k | {z } quadratic model of gat zk +L(z;λ)o,(C1) where kuk2 ψ,k := hu, ( ∇2ψ ( zk )) ui and αk> 0. This is a proximal map for f with a linearized smooth part g: Optimality: ∇2ψ(zk)zk+1 −zk αk +∇g(zk) + ∂L(zk+1) + NC(zk+1)30.(C2) In Euclidean geometry (ψ(z) = 1 2kzk2), (C2) reduces to the projected implicit–explicit step zk+1 = ProjCzk−αk(∇g(zk) + ξk+1), ξk+1 ∈∂L(zk+1). 3. Descent under soft tube and trust Assume ∇g is L –Lipschitz in a neighborhood of C and L convex. A standard descent lemma with a projected implicit step yields: Theorem C.1 (One–step descent with tube trust).Let zk+1 solve (C1) and suppose αk≤1/L. Then g(zk+1) + L(zk+1)≤g(zk) + L(zk)−1 2αkkzk+1 −zkk2 ψ,k, and if, in addition, the tube trust holds Ltube,β(zk+1)≤Ltube,β(zk) + δk, with Pkδk<∞, then Pkkzk+1 −zkk2 ψ,k <∞and {g(zk) + L(zk)}decreases to a finite limit. Proof. By Lipschitzness, g ( zk+1 ) ≤g ( zk ) + h∇g ( zk ) , zk+1 −zki + L 2kzk+1 −zkk2 . Optimality (C2) with αk≤ 1 /L and convexity of L give the claimed quadratic decrease. The trust condition prevents unbounded growth of the tube term; summability implies finite length. Corollary C.2 (Normal–component contraction) . If Ltube,β is σ⊥ –strongly convex in the normal directions to Γ, then the accepted step satisfies ku⊥k2≤2δk/σ⊥. Proof. By strong convexity, Ltube,β(zk+1)−Ltube,β(zk)≥σ⊥ 2ku⊥k2. Appendix D: Persistent homology and sheaf gluing: constructions and stability 1. Contact–graph filtrations and persistence Given z , extract a weighted graph Gz = ( V, E, w )on residues (or coarse patches) with weights wij ∈R≥0 reflecting contact strengths. For t∈R define a sublevel filtration Gt = ( V, { ( i, j ) : wij ≥t} )and the associated clique complex Xt (flag complex). Persistent homology computes homology groups Hk ( Xt ) over a field across t, summarizing births and deaths in a multiset (diagram) Dk[196, 197]. Stability. For two filtrations from weight functions w, w0 , the bottleneck distance between diagrams satisfies (Cohen–Steiner stability) [198]: dB(Dk(w), Dk(w0)) ≤ kw−w0k∞. Sliced Wasserstein and persistence–image distances inherit stability up to constants [199]. 64 [20] H. Alt and M. Godau. Computing the Fréchet distance between two polygonal curves. International Journal of Computational Geometry & Applications, 5(1–2):75–91, 1995. [21] T. Hastie, R. Tibshirani, and J. Friedman. The Elements of Statistical Learning, 2nd ed. Springer, 2009. [22] S. Gavrilets. Fitness Landscapes and the Origin of Species. Princeton University Press, 2004. [23] D. M. Weinreich, N. F. Delaney, M. A. Depristo, and D. L. Hartl. Darwinian evolution can follow only very few mutational paths to fitter proteins. Science, 312(5770):111–114, 2006. [24] I. G. Szendro, M. F. Schenk, J. Frank, J. Krug, and J. A. G. de Visser. Quantitative analyses of empirical fitness landscapes. J. Theor. Biol., 425:82–100, 2017. [25] J. A. G. de Visser and J. Krug. Empirical fitness landscapes and the predictability of evolution. Nat. Rev. Genet., 15(7):480–490, 2014. [26] I. Fragata, A. Blanckaert, C. Lohmueller, and C. Bank. Evolution in the light of fitness landscape theory. Trends Ecol. Evol., 34(1):69–82, 2019. [27] T. Bedford, S. R. Stern, J. D. Bloom, et al. Integrating influenza antigenic dynamics with molecular evolution. eLife, 3:e01914, 2014. [28] M. Lüksza and M. Lässig. A predictive fitness model for influenza. Nature, 507:57–61, 2014. [29] R. A. Neher, T. Bedford, R. S. Daniels, et al. Prediction, dynamics, and the immune response to influenza antigenic drift. Proc. Natl. Acad. Sci. USA, 113(12):E1701–E1709, 2016. [30] J. M. Fonville, S. C. Wilks, et al. Antibody landscapes after influenza virus infection or vaccination. Science, 346(6212):996–1000, 2014. [31] O. G. Pybus and A. Rambaut. Evolutionary analysis of the dynamics of viral infectious disease. Nat. Rev. Genet., 10:540–550, 2009. [32] A. J. Drummond and A. Rambaut. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol., 7:214, 2007. [33] T. Stadler, D. Kuhnert, S. Bonhoeffer, and A. J. Drummond. Birth–death skyline plot reveals temporal changes of epidemic spread. Proc. Natl. Acad. Sci. USA, 110(1):228–233, 2013. [34] M. D. Karcher, D. J. Palser, T. Bedford, et al. Quantifying and modeling pathogen evolution and transmission. Trends Microbiol., 24(10):872–885, 2016. [35] C. L. Araya and D. M. Fowler. Deep mutational scanning: assessing protein function on a massive scale. Curr. Opin. Struct. Biol., 21(4):484–490, 2011. [36] J. B. Kinney and D. M. McCandlish. Massively parallel assays and the analysis of high-throughput biological data. Annu. Rev. Genomics Hum. Genet., 20:99–127, 2019. [37] J. Otwinowski, D. M. McCandlish, and J. B. Plotkin. Inferring the shape of global epistasis. Proc. Natl. Acad. Sci. USA, 115(32):E7550–E7558, 2018. [38] R. Jordan, D. Kinderlehrer, and F. Otto. The variational formulation of the Fokker–Planck equation. SIAM J. Math. Anal., 29(1):1–17, 1998. [39] L. Ambrosio, N. Gigli, and G. Savaré. Gradient Flows in Metric Spaces and in the Space of Probability Measures. Birkhäuser, 2nd ed., 2008. [40] C. Villani. Topics in Optimal Transportation. American Mathematical Society, 2003. [41] D. Cohen-Steiner, H. Edelsbrunner, and J. Harer. Stability of persistence diagrams. Discrete Comput. Geom., 37(1):103–120, 2007. [42] F. Chazal, V. de Silva, M. Glisse, and S. Oudot. The Structure and Stability of Persistence Modules. Springer, 2016. [43] S. Y. Oudot. Persistence Theory: From Quiver Representations to Data Analysis. American Mathematical Society, 2015. [44] G. Carlsson. Topology and data. Bull. Amer. Math. Soc., 46(2):255–308, 2009. [45] J. Hansen and R. Ghrist. Toward a spectral theory of cellular sheaves. Found. Comput. Math., 19:999–1040, 2019. [46] F. W. Lawvere. Metric spaces, generalized logic, and closed categories. Rendiconti del Seminario Matematico e Fisico di Milano, 43:135–166, 1973. [47] G. M. Kelly. Basic Concepts of Enriched Category Theory. Cambridge University Press, 1982. [48] S. Mac Lane. Categories for the Working Mathematician, 2nd ed. Springer, 1998. [49] T. Leinster. Basic Category Theory. Cambridge University Press, 2014. [50] D. I. Spivak. Category Theory for the Sciences. MIT Press, 2014. [51] D. A. Levin, Y. Peres, and E. L. Wilmer. Markov Chains and Mixing Times. American Mathematical Society, 2009. [52] U. Bauer and M. Lesnick. Induced matchings and the algebraic stability of persistence barcodes. J. Comput. Geom., 6(2):162–191, 2015. [53] S. Kobayashi and K. Nomizu. Foundations of Differential Geometry, Vol. I. Wiley, 1963. [54] J. C. Baez and U. Schreiber. Higher gauge theory. In Categories in Algebra, Geometry and Mathematical Physics, Contemp. Math. 431, 7–30, 2007. [55] A. Banerjee, S. Merugu, I. S. Dhillon, and J. Ghosh. Clustering with Bregman divergences. J. Mach. Learn. Res., 6:1705–1749, 2005. [56] H. H. Bauschke and P. L. Combettes. Convex Analysis and Monotone Operator Theory in Hilbert Spaces. Springer, 2011. [57] S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, 2004. [58] H. Brézis. Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert. 65 North-Holland, 1973. [59] J.-B. Hiriart-Urruty and C. Lemaréchal. Convex Analysis and Minimization Algorithms I–II. Springer, 1993. [60] E. De Giorgi. New problems on minimizing movements. In Boundary Value Problems for PDE and Applications, 81–98, Masson, 1993. [61] H. Attouch, J. Bolte, and B. F. Svaiter. Convergence of descent methods for semi-algebraic and tame problems. Math. Program., 137(1–2):91–129, 2013. [62] J.-P. Aubin. Viability Theory. Birkhäuser, 1991. [63] J.-P. Aubin and H. Frankowska. Set-Valued Analysis. Birkhäuser, 2nd ed., 2009. [64] D. P. Bertsekas. Nonlinear Programming, 2nd ed. Athena Scientific, 1999. [65] A. S. Nemirovski and D. B. Yudin. Problem Complexity and Method Efficiency in Optimization. Wiley, 1983. [66] D. G. Myszka. Survey of the 1998 optical biosensor literature. J. Mol. Recognit., 12(6):390–408, 1999. [67] P. Schuck. Reliable determination of binding affinity and kinetics by surface plasmon resonance biosensing. Curr. Opin. Biotechnol., 8(4):498–502, 1997. [68] Y. N. Abdiche, D. M. Malashock, A. Pinkerton, and J. P. Pons. Determining kinetics and affinities of protein interactions using a parallel real-time label-free biosensor, the Octet. Anal. Biochem., 377(2):209–217, 2008. [69] R. L. Rich and D. G. Myszka. Grading the commercial optical biosensor literature—class of 2008: ‘the mighty binders’ return. J. Mol. Recognit., 23(1):1–64, 2010. [70] M. Hoffmann, H. Kleine-Weber, et al. SARS-CoV-2 cell entry depends on ACE2 and TMPRSS2 and is blocked by a clinically proven protease inhibitor. Cell, 181(2):271–280, 2020. [71] B. Coutard, C. Valentin, et al. The spike glycoprotein of the new coronavirus 2019-nCoV contains a furin-like cleavage site. Antiviral Res., 176:104742, 2020. [72] F. H. Niesen, H. Berglund, and M. Vedadi. The use of differential scanning fluorimetry to detect ligand interactions that promote protein stability. Nat. Protoc., 2(9):2212–2221, 2007. [73] C. N. Pace, B. A. Shirley, and J. M. Thomson. Measuring and increasing protein stability. FASEB J., 9(7): 2005–2011, 1995. [74] P. L. Privalov. Stability of proteins: small globular proteins. Adv. Protein Chem., 33:167–241, 1979. [75] J. F. Brandts and L.-N. Lin. Study of strong to ultrahigh–affinity protein interactions using differential scanning calorimetry. Biochemistry, 29(29):6927–6940, 1990. [76] P. A. Kristiansen, D. Page, et al. WHO International Standard for anti-SARS-CoV-2 immunoglobulin. The Lancet, 397(10282):1347–1348, 2021. [77] J. Nie, Q. Li, et al. Establishment and validation of a pseudovirus neutralization assay for SARS-CoV-2. Emerg. Microbes Infect., 9(1):680–686, 2020. [78] F. Schmidt, F. Weisblum, et al. Measuring SARS-CoV-2 neutralizing antibody activity using pseudotyped and chimeric viruses. J. Exp. Med., 217(11):e20201181, 2020. [79] K. A. Earle, D. M. Ambrosino, et al. Evidence for antibody as a protective correlate for COVID-19 vaccines. Vaccine, 39(32):4423–4428, 2021. [80] R. DerSimonian and N. Laird. Meta-analysis in clinical trials. Controlled Clinical Trials, 7(3):177–188, 1986. [81] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin. Bayesian Data Analysis, 3rd ed. CRC Press, 2013. [82] J. M. Bland and D. G. Altman. Statistical methods for assessing agreement between two methods of clinical measurement. The Lancet, 327(8476):307–310, 1986. [83] K. Turner, Y. Mileyko, S. Mukherjee, and J. Harer. Fréchet means for distributions of persistence diagrams. Discrete & Computational Geometry, 52(1):44–70, 2014. [84] M. Cerrière, M. Cuturi, and S. Oudot. Sliced Wasserstein kernel for persistence diagrams. In Proc. 34th ICML, PMLR 70:664–673, 2017. [85] H. Adams, T. Emerson, et al. Persistence images: A stable vector representation of persistent homology. J. Mach. Learn. Res., 18(8):1–35, 2017. [86] J. R. Norris. Markov Chains. Cambridge University Press, 1997. [87] A. Sinclair and M. Jerrum. Approximate counting, uniform generation and rapidly mixing Markov chains. Information and Computation, 82(1):93–133, 1989. [88] C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. [89] W. Hoeffding. Probability inequalities for sums of bounded random variables. J. Amer. Stat. Assoc., 58(301):13–30, 1963. [90] K. M. Miettinen. Nonlinear Multiobjective Optimization. Springer, 1999. [91] R. A. Bradley and M. E. Terry. Rank analysis of incomplete block designs: I. The method of paired comparisons. Biometrika, 39(3/4):324–345, 1952. [92] R. L. Plackett. The analysis of permutations. J. Royal Stat. Soc. Series C, 24(2):193–202, 1975. [93] T. Joachims. Optimizing search engines using clickthrough data. In Proc. KDD, 133–142, 2002. [94] R. Herbrich, T. Graepel, and K. Obermayer. Large margin rank boundaries for ordinal regression. Advances in Large Margin Classifiers, MIT Press, 115–132, 2000. [95] L. Jacob, G. Obozinski, and J.-P. Vert. Group lasso with overlap and graph lasso. In Proc. ICML, 433–440, 2009. [96] O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J. Math. Biol., 28(4):365–382, 1990. [97] R. T. Rockafellar and R. J.-B. Wets. Variational Analysis. Springer, 1998. 66 [98] L. Armijo. Minimization of functions having Lipschitz continuous first partial derivatives. Pacific Journal of Mathematics, 16(1):1–3, 1966. [99] J. Nocedal and S. J. Wright. Numerical Optimization, 2nd ed. Springer, 2006. [100] A. R. Conn, N. I. M. Gould, and P. L. Toint. Trust-Region Methods. SIAM, 2000. [101] C. Cartis, N. I. M. Gould, and P. L. Toint. Adaptive cubic regularisation methods for unconstrained optimization. Part I: Motivation, convergence and numerical results. Math. Program., 127(2):245–295, 2011. [102] P. H. C. Eilers and B. D. Marx. Flexible smoothing with B -splines and penalties. Statistical Science, 11(2):89–121, 1996. [103] R. J. Tibshirani. Adaptive piecewise polynomial estimation via trend filtering. Annals of Statistics, 42(1):285–323, 2014. [104] J. Durbin and S. J. Koopman. Time Series Analysis by State Space Methods, 2nd ed. Oxford University Press, 2012. [105] N. Li and M. Stephens. Modeling linkage disequilibrium and identifying recombination hotspots using SNP data. Genetics, 165(4):2213–2233, 2003. [106] P. E. Hart, N. J. Nilsson, and B. Raphael. A formal basis for the heuristic determination of minimum cost paths. IEEE Trans. Systems Science and Cybernetics, 4(2):100–107, 1968. [107] J. Pearl. Heuristics: Intelligent Search Strategies for Computer Problem Solving. Addison–Wesley, 1984. [108] R. Zhou and E. A. Hansen. Beam-stack search: Integrating backtracking with beam search. Artificial Intelligence, 173(15): 1–26, 2009. [109] E. W. Dijkstra. A note on two problems in connexion with graphs. Numerische Mathematik, 1:269–271, 1959. [110] A. Paten, B. Earl, et al. Genome graphs and the evolution of genome inference. Genome Research, 27(5):665–676, 2017. [111] E. Garrison, A. Sirén, et al. Variation graph toolkit improves read mapping by representing genetic variation in the reference. Nature Methods, 15: 1–5, 2018. [112] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein. Introduction to Algorithms, 3rd ed. MIT Press, 2009. [113] C. Wiuf and J. Hein. Recombination as a point process along sequences. Theoretical Population Biology, 55(3):248–259, 1999. [114] H. Eyring. The activated complex in chemical reactions. J. Chem. Phys., 3(2):107–115, 1935. [115] H. A. Kramers. Brownian motion in a field of force and the diffusion model of chemical reactions. Physica, 7(4):284–304, 1940. [116] M. I. Freidlin and A. D. Wentzell. Random Perturbations of Dynamical Systems, 3rd ed. Springer, 2012. [117] A. Bovier and F. den Hollander. Metastability: A Potential-Theoretic Approach. Springer, 2015. [118] F. H. Clarke. Optimization and Nonsmooth Analysis, 2nd ed. SIAM, 1990. [119] Y. Nesterov and A. Nemirovskii. Interior-Point Polynomial Methods in Convex Programming. SIAM, 1994. [120] J. Milnor. Morse Theory. Princeton University Press, 1963. [121] R. Thom. Structural Stability and Morphogenesis. Benjamin, 1975. [122] B. Fasy, F. Lecci, A. Rinaldo, L. Wasserman, S. Balakrishnan, and A. Singh. Confidence sets for persistence diagrams. Annals of Statistics, 42(6):2301–2339, 2014. [123] P. Bubenik. Statistical topological data analysis using persistence landscapes. J. Mach. Learn. Res., 16:77–102, 2015. [124] W. Ambrose and I. M. Singer. A theorem on holonomy. Trans. Amer. Math. Soc., 75(3):428–443, 1953. [125] G. F. Lawler and A. D. Sokal. Bounds on the L2 spectrum for Markov chains and Markov processes: A generalization of Cheeger’s inequality. Trans. Amer. Math. Soc., 309(2):557–580, 1988. [126] L. Lovász and R. Kannan. Faster mixing via average conductance. In Proc. 31st ACM STOC, 282–287, 1999. [127] R. Montenegro and P. Tetali. Mathematical Aspects of Mixing Times in Markov Chains. Foundations and Trends in Theoretical Computer Science, 2006. [128] Y. Benjamini and Y. Hochberg. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society, Series B, 57(1):289–300, 1995. [129] Y. Benjamini and D. Yekutieli. The control of the false discovery rate in multiple testing under dependency. Annals of Statistics, 29(4):1165–1188, 2001. [130] B. Efron. Bootstrap methods: another look at the jackknife. Annals of Statistics, 7(1):1–26, 1979. [131] E. S. Page. Continuous inspection schemes. Biometrika, 41(1/2):100–115, 1954. [132] A. Wald. Sequential Analysis. Wiley, 1947. [133] R. Killick, P. Fearnhead, and I. A. Eckley. Optimal detection of changepoints with a linear computational cost. Journal of the American Statistical Association, 107(500):1590–1598, 2012. [134] E. R. DeLong, D. M. DeLong, and D. L. Clarke-Pearson. Comparing the areas under two or more correlated ROC curves: a nonparametric approach. Biometrics, 44(3):837–845, 1988. [135] F. R. Hampel. The influence curve and its role in robust estimation. J. Amer. Stat. Assoc., 69(346):383–393, 1974. [136] P. J. Huber. Robust Statistics. Wiley, 1981. [137] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, Series B, 58(1):267–288, 1996. [138] A. E. Hoerl and R. W. Kennard. Ridge regression: biased estimation for nonorthogonal problems. Techno- 67 metrics, 12(1):55–67, 1970. [139] J. Bergstra and Y. Bengio. Random search for hyper-parameter optimization. Journal of Machine Learning Research, 13:281–305, 2012. [140] T. M. Cover and P. E. Hart. Nearest neighbor pattern classification. IEEE Transactions on Information Theory, 13(1):21–27, 1967. [141] C. Goodall. Procrustes methods in the statistical analysis of shape. Journal of the Royal Statistical Society, Series B, 53(2):285–339, 1991. [142] Y. Rubner, C. Tomasi, and L. J. Guibas. The earth mover’s distance as a metric for image retrieval. International Journal of Computer Vision, 40(2):99–121, 2000. [143] A. J. Riesselman, J. B. Ingraham, and D. S. Marks. Deep generative models of genetic variation capture the effects of mutations. Nature Methods, 15:816–822, 2018. [144] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In Proc. ICLR, 2015. [145] S. Arlot and A. Celisse. A survey of cross-validation procedures for model selection. Statistics Surveys, 4:40–79, 2010. [146] G. C. Cawley and N. L. C. Talbot. On over-fitting in model selection and subsequent selection bias in performance evaluation. Journal of Machine Learning Research, 11:2079–2107, 2010. [147] M. Cuturi and M. Blondel. Soft-DTW: a differentiable loss function for time-series. In Proceedings of the 34th International Conference on Machine Learning, PMLR 70:894–903, 2017. [148] M.-P. Dubuisson and A. K. Jain. A modified Hausdorff distance for object matching. In Proceedings of the 12th International Conference on Pattern Recognition, 566–568, 1994. [149] M. P. do Carmo. Differential Geometry of Curves and Surfaces. Prentice–Hall, 1976. [150] C. J. Willmott and K. Matsuura. Advantages of the mean absolute error (MAE) over the root mean square error (RMSE) in assessing average model performance. Climate Research, 30(1):79–82, 2005. [151] B. Efron and R. J. Tibshirani. An Introduction to the Bootstrap. Chapman & Hall/CRC, 1994. [152] P. Good. Permutation Tests: A Practical Guide to Resampling Methods for Testing Hypotheses, 2nd ed. Springer, 2000. [153] J. Cohen. Statistical Power Analysis for the Behavioral Sciences, 2nd ed. Lawrence Erlbaum, 1988. [154] L. V. Hedges. Distribution theory for Glass’s estimator of effect size and related estimators. Journal of Educational Statistics, 6(2):107–128, 1981. [155] S. Holm. A simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics, 6(2):65–70, 1979. [156] W. Karush. Minima of functions of several variables with inequalities as side conditions. Master’s Thesis, University of Chicago, 1939. [157] H. W. Kuhn and A. W. Tucker. Nonlinear programming. In Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, 1951. [158] C. H. Waddington. The Strategy of the Genes. Allen & Unwin, 1957. [159] B. C. Hall. Lie Groups, Lie Algebras, and Representations, 2nd ed. Springer, 2015. [160] S. Wright. The roles of mutation, inbreeding, crossbreeding, and selection in evolution. Proceedings of the Sixth International Congress of Genetics, 1:356–366, 1932. [161] R. A. Fisher. The Genetical Theory of Natural Selection. Clarendon Press, 1930. [162] A. Braides. Γ-Convergence for Beginners. Oxford University Press, 2002. [163] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko. The Mathematical Theory of Optimal Processes. Interscience, 1962. [164] J. F. Bonnans and A. Shapiro. Perturbation Analysis of Optimization Problems. Springer, 2000. [165] J. A. Tropp. User-friendly tail bounds for sums of random matrices. Foundations of Computational Mathematics, 12(4):389–434, 2012. [166] B. Jones and M. G. Kenward. Design and Analysis of Cross-Over Trials, 3rd ed. CRC Press, 2014. [167] L. P. Kaelbling, M. L. Littman, and A. R. Cassandra. Planning and acting in partially observable stochastic domains. Artificial Intelligence, 101(1–2):99–134, 1998. [168] E. Hazan. Introduction to Online Convex Optimization. Now Publishers, 2016. [169] W. E. Johnson, C. Li, and A. Rabadan. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics, 8(1):118–127, 2007. [170] J. P. T. Higgins and S. G. Thompson. Quantifying heterogeneity in a meta-analysis. Statistics in Medicine, 21(11):1539–1558, 2002. [171] R. J. Carroll, D. Ruppert, L. A. Stefanski, and C. M. Crainiceanu. Measurement Error in Nonlinear Models, 2nd ed. Chapman & Hall/CRC, 2006. [172] D. B. Rubin. Multiple Imputation for Nonresponse in Surveys. Wiley, 1987. [173] D. M. Fowler and S. Fields. Deep mutational scanning: a new style of protein science. Nature Methods, 11:801–807, 2014. [174] T. N. Starr, A. J. Greaney, et al. Deep mutational scanning of SARS-CoV-2 receptor binding domain reveals constraints on folding and ACE2 binding. Cell, 182(5):1295–1310, 2020. [175] Y. Turakhia, N. De Maio, et al. Pandemic-scale phylogenomics reveals the SARS-CoV-2 recombination landscape. Nature Biotechnology, 39:1403–1410, 2021. [176] J. Hadfield, C. Megill, et al. Nextstrain: real-time tracking of pathogen evolution. Bioinformatics, 34(23):4121–4123, 2018. [177] P. L. Bartlett and S. Mendelson. Rademacher and Gaussian complexities: risk bounds and structural results. 68 Journal of Machine Learning Research, 3:463–482, 2002. [178] P. Artzner, F. Delbaen, J.-M. Eber, and D. Heath. Coherent measures of risk. Mathematical Finance, 9(3):203–228, 1999. [179] R. T. Rockafellar and S. Uryasev. Optimization of conditional value-at-risk. Journal of Risk, 2(3):21–41, 2000. [180] A. Cori, N. M. Ferguson, C. Fraser, and S. Funk. A new framework and software to estimate time-varying reproduction numbers during epidemics. American Journal of Epidemiology, 178(9):1505–1512, 2013. [181] D. J. Smith, A. S. Lapedes, et al. Mapping the antigenic and genetic evolution of influenza virus. Science, 305(5682):371–376, 2004. [182] B. T. Grenfell, O. G. Pybus, J. R. Gog, et al. Unifying the epidemiological and evolutionary dynamics of pathogens. Science, 303(5656):327–332, 2004. [183] M. Scheffer, J. Bascompte, et al. Early-warning signals for critical transitions. Nature, 461:53–59, 2009. [184] V. Dakos, et al. Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data. PLoS ONE, 7(7):e41010, 2012. [185] J. D. Orth, I. Thiele, and B. Ø. Palsson. What is flux balance analysis? Nature Biotechnology, 28:245–248, 2010. [186] E. D. Sontag. Mathematical Control Theory, 2nd ed. Springer, 1998. [187] S. Prajna, A. Jadbabaie, and G. J. Pappas. A framework for worst-case analysis of switched hybrid systems. IEEE Transactions on Automatic Control, 52(8):1480–1494, 2007. [188] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada. Control barrier function based quadratic programs for safety critical systems. IEEE Transactions on Automatic Control, 62(8):3861–3876, 2017. [189] J. C. Baez and B. Fong. A compositional framework for Markov processes. Journal of Mathematical Physics, 57(3):033301, 2016. [190] B. Fong and D. I. Spivak. Seven Sketches in Compositionality: An Invitation to Applied Category Theory. Cambridge University Press, 2019. [191] R. T. Rockafellar. Monotone operators and the proximal point algorithm. SIAM Journal on Control and Optimization, 14(5):877–898, 1976. [192] J. J. Moreau. Proximité et dualité dans un espace hilbertien. Bulletin de la Société Mathématique de France, 93:273–299, 1965. [193] I. Ekeland. On the variational principle. Journal of Mathematical Analysis and Applications, 47(2):324–353, 1974. [194] S. Mac Lane. Categories for the Working Mathematician. Springer, 1971. [195] S. Awodey. Category Theory, 2nd ed. Oxford University Press, 2010. [196] H. Edelsbrunner and J. Harer. Computational Topology: An Introduction. American Mathematical Society, 2010. [197] S. Y. Oudot. Persistence Theory: From Quiver Representations to Data Analysis. American Mathematical Society, 2015. [198] D. Cohen-Steiner, H. Edelsbrunner, and J. Harer. Stability of persistence diagrams. Discrete & Computational Geometry, 37(1):103–120, 2007. [199] F. Chazal and B. Michel. An introduction to topological data analysis: fundamental and practical aspects for data scientists. arXiv:1710.04019, 2017. [200] R. Ghrist and V. de Silva. Sheaf theoretic sensor integration. Netw. Heterog. Media, 7(2):215–234, 2012. [201] J. T. Hansen and R. Ghrist. Toward a spectral theory of cellular sheaves. Journal of Applied and Computational Topology, 3:315–358, 2019. [202] S. Kobayashi and K. Nomizu. Foundations of Differential Geometry, Vol. 1. Wiley, 1963. [203] J. Kleinberg and É. Tardos. Algorithm Design. Addison–Wesley, 2005. [204] A. Saltelli, M. Ratto, T. Andres, et al. Global Sensitivity Analysis. Wiley, 2008. [205] I. M. Sobol’. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Mathematics and Computers in Simulation, 55(1–3):271–280, 2001. [206] M. D. McKay, R. J. Beckman, and W. J. Conover. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 21(2):239–245, 1979. [207] M. G. Kenward and B. Jones. Design and Analysis of Cross-Over Trials, 2nd ed. Chapman & Hall/CRC, 2003. [208] D. Boni, C. Posada, M. F. Fitch, and A. Goldman, “An Exact Nonparametric Method for Identifying Mosaic Recombinant Sequences,” Molecular Biology and Evolution, vol. 18, no. 5, pp. 687–695, 2001. [209] D. P. Martin, B. Murrell, M. Golden, A. Khoosal, and B. Muhire, “RDP4: Detection and Analysis of Recombination Patterns in Virus Genomes,” Virus Evolution, vol. 1, no. 1, vev003, 2015. [210] S. L. Kosakovsky Pond, D. Posada, M. B. Gravenor, C. H. Woelk, and S. D. W. Frost, “Automated Phylogenetic Detection of Recombination Using a Genetic Algorithm,” Molecular Biology and Evolution, vol. 23, no. 10, pp. 1891–1901, 2006. [211] Block lengths can be chosen by automatic methods to match the estimated dependence range.