# From child selections to a complete futile-cycle model

Example for manuscript entry 38, **Futile-cycle child-selection determinants and unstable cores**, SHA-256 `370eea21e788907581cc3c677783a0bd1581afe53307babfdc2325e77120f2aa`.

A child selection assigns a distinct consuming reaction to each selected species. Its square stoichiometric matrix records a potential feedback structure. This paper shows how simple enzyme-conversion chemistry constrains every such determinant, yet admits several different kinds of unstable core. The example connects exact structural calculations to a complete saturating kinetic network with conserved enzyme and substrate totals.

## Run

Python 3.11 or newer:

```sh
python -m venv .venv
# Windows: .venv\Scripts\activate
# macOS/Linux: source .venv/bin/activate
python -m pip install -r requirements.txt
python -m unittest -v test_example.py
python example.py --output outputs
```

Inputs at the top of `example.py` control the network size, exhaustive census sizes/budget, kinetic reactivities, perturbation, integration interval, an independent seven-dimensional rate/loss study and asymptotic site counts. The default exhaustive three-site census is the most expensive step. Set `CENSUS_SITES=(2,)` when exploring kinetics repeatedly. Exceeding the census budget raises an error rather than reporting incomplete counts as exhaustive. Explicit networks are limited to 30 sites; scalar spectral calculations support up to 10,000 sites. Full negative-core kinetics require at least three sites. The seven-element flower inputs deliberately describe a separate three-site matrix study.

Outputs include exact certificates and characteristic coefficients in JSON, complete numerical trajectory and sweep CSVs, console output, three PNG/SVG figure pairs and hashes. Tests take a few seconds. No research directory or Lean installation is needed.

## Chemical objects and reusable components

Each `ConversionArm` contains a substrate s, product t, enzyme e and unique intermediate c, with reactions

```text
s + e -> c       binding
c -> s + e       unbinding
c -> t + e       conversion
```

`ConversionSystem` builds literal reactant/product/stoichiometric matrices, validates disjoint substrate/enzyme/intermediate roles, constructs conserved pools, validates a child selection and enumerates all selections. Arms can share substrates and enzymes; a conversion may return the same substrate. A substrate acting as an enzyme is outside this model's role assumptions. Every intermediate belongs to exactly one arm.

`ConversionSystem.futile(n)` builds the n-site phosphorylation/dephosphorylation network. Kinase K converts S_i to S_(i+1) through C_i, and phosphatase F converts S_(i+1) to S_i through D_(i+1). There are 3n+3 species and 6n reactions. Stored species order is substrates, enzymes, then arm intermediates; use names rather than assumed matrix positions.

`ChildSelection` owns a selected matrix and its determinant certificate. Ordering free species before selected intermediates gives a block matrix `[B W; C -I]`, whose determinant is `(-1)^q det(B+WC)`. Changing signs on enzyme rows leaves at most one +1 and one -1 per column. Explicit cofactor reduction proves each determinant belongs to {-1,0,1}. This is a zero-frequency determinant identity, not spectral equivalence after eliminating intermediates.

The complete censuses, including the empty selection, are:

| Sites | Children | det -1 | det 0 | det +1 |
| --- | ---: | ---: | ---: | ---: |
| 2 | 4,455 | 543 | 3,377 | 535 |
| 3 | 204,849 | 15,649 | 173,559 | 15,641 |

## Ordinary instability and instability after scaling

The six-species negative-feedback child uses owners `(S1,S3,K,F,C2,D1)` and reactions `(K1b,F3b,K2b,F1b,K2v,F1v)`. Exact Routh arithmetic gives two roots in the open right half-plane. All 62 nonempty proper principal restrictions are nonunstable at unit rates: all but one reduce to singleton or positive-cycle blocks, and the exceptional five-dimensional B has an explicit positive-definite Lyapunov matrix.

B itself is stable at unit rates, but B*diag(1,1,2,1,2) has two unstable roots. Its 30 proper nonempty principal restrictions remain nonunstable for every positive column scaling. Thus the six-dimensional ordinary unstable core is not minimal for instability under scaling. Exact signs, Lyapunov leading minors, Routh columns and the paper's logarithmic-derivative value are regenerated. A degenerate Routh pivot is reported unresolved rather than guessed.

The positive-feedback child has 2n+1 species. A diagonal sign similarity converts it into a matrix `M=-I+adjacency` with a long directed cycle of length 2n and a two-cycle sharing a hub. This sign transformation is a spectral device; it does not turn the chemical network into an autocatalytic mechanism. Proper restrictions contain only isolated vertices and unit positive cycles, proving noninstability under positive column scaling. The full matrix is unstable under every positive scaling. At n=10, its exact Routh calculation gives **three** unstable roots: one positive real root and a nonreal pair. A unique positive real root does not mean a unique unstable root.

## Rates, losses and a scalar spectral solver

`FlowerSpectrum` accepts positive column rates d and nonnegative independent losses eta. For `H=M*diag(d)-diag(eta)`, define a_i=d_i+eta_i. The spectral abscissa beta is the unique root above `-min(a)` of

```text
product_long_cycle d_i/(beta+a_i)
  + d_hub*d_leaf/((beta+a_hub)*(beta+a_leaf)) = 1.
```

The solver evaluates this equation in logarithms, avoiding a large dense eigensolve. It also constructs the positive eigenvector. Exact rational evaluation at beta=0 determines stability: gain less than, equal to or greater than one means negative, zero or positive spectral abscissa. Independent uniform loss delta subtracts delta from the abscissa. With no additional loss, increasing any column rate increases the positive abscissa, bounded between d_min and d_max times the unit-rate value.

At unit rates, alpha solves `(1+alpha)^(-2n)+(1+alpha)^(-2)=1`. Lambert-W lower/upper expressions and the asymptotic approximation W(n-1)/(2(n-1)) show why larger cores have slower growth. The plotted roots and Lambert evaluations are numerical evaluations of exact formulas, not interval certificates. Independent selected-matrix rates and losses are not automatically realizable interventions on the complete conserved chemical system.

## Full kinetics and reuse

`SaturatingKinetics` takes a network, positive reference concentration xbar, positive balanced flux vbar and reactant-supported reactivity matrix R. It requires `S*vbar=0` and `0 < theta_ji=R_ji*xbar_i/vbar_j < 1`. Reaction rates are

```text
v_j(x) = vbar_j * product_reactants z_i/(theta_ji+(1-theta_ji)*z_i),
z_i = x_i/xbar_i.
```

This realizes exactly v(xbar)=vbar and Dv(xbar)=R. It is a saturating kinetic construction, **not mass action**. The default reference has all concentrations one, binding flux four, unbinding/conversion flux two, background reactivity 1/100 and selected negative-child reactivities one. Units are nondimensional; these are mathematical witness inputs, not fitted biochemical parameters.

For three sites, kinase total is 4, phosphatase total is 4 and total free/bound substrate is 10. Exact rank is nine. The full Jacobian has three conservation zeros and exactly three unstable roots on the nine-dimensional stoichiometric subspace. The leading complex pair is about 0.113319 +/- 0.389512i. An unstable eigenvector supplies a small perturbation in the same conserved class. Radau and BDF integrate the nonlinear model independently; a linear prediction supplies a separate local comparison. Integration uses departures from the reference and cancellation-resistant rate differences; it rejects positive-domain exits instead of clipping concentrations.

To reuse, compose new arms, provide a balanced flux, choose a reference state and a supported R within the saturation constraint. Vary these inputs or use `FlowerSpectrum` for inexpensive matrix studies. Structural conclusions require the stated chemical roles; full-network stability must be recomputed for changed kinetics. Conservation rows need not be independent for a disconnected custom network, so calculate their rank rather than assuming three independent pools.

This example establishes no Hopf bifurcation, periodic orbit, global behavior or experimental realism. It does not rerun Lean or exhaustively classify every unstable core of larger networks.

## Provenance and licensing

The reaction construction, matrices and exact constants were checked against the supplied manuscript and accompanying `Futile_Cycle_Child_Selection_arxiv/check_paper.py`. The implementation is self-contained and independently checks its emitted certificates. MIT is proposed for the example code, pending owner approval; this README does not grant a license or change manuscript rights.