[Download the complete reproduction package](/publications/optimal-affinity/versions/v1/reproduction.zip) or [the Python source](/publications/optimal-affinity/versions/v1/optimal_affinity.py). The saved run used Python 3.12.14.

This program makes the paper's main constructions and examples explicit. It starts with reaction counts, constructs rate constants, checks the steady state and its responses, and verifies algebraic identities that prove the stated global maxima. The figures show actual saved calculations, so no installation is needed to inspect the results.

## The main result in numbers

For X ⇌ Y with rate constants (3,2), and 2Y ⇌ 2X+Y with rate constants (7,6), the unique greatest production rate occurs at x=y=1. The forward rates are 3 and 7, the reverse rates are 2 and 6, and both net reaction currents equal 1. Thus Y is stationary and the net production of X is −1+2=1.

The affinity at this maximum is

$$
A^*=\log(3/2)+\log(7/6)=\log(7/4)\approx0.5596157879,
$$

below the conjectured minimum log(2) ≈ 0.6931471806. This comparison concerns a maximum of the production rate, not a minimum of affinity along that one network's entire steady-state curve.

![Production peaks at x=1, where the affinity is below the conjectured bound.](/publications/optimal-affinity/versions/v1/counterexample.png)

The left panel shows production. The right shows affinity at the same steady states. At x=7/4 the system reaches detailed balance, with zero net production and zero affinity. The plotted window is finite; the following identity proves the global conclusion for every positive steady state.

Set E=3x−2y−7y²+6x²y, s=2y, and J=3x−2y. The program checks exactly that

$$
12E=(s-2)^2(4s+3)+(J-1)[4s(J+1)+8s^2+12].
$$

At a steady state E=0 and s>0. If J>1, the right side is positive, a contradiction. If J=1, it forces s=2, hence x=y=1. This is why an algebraic identity, together with its sign argument, establishes a global maximum rather than merely suggesting one from a graph.

## The same molecule count permits very different affinities

The family X ⇌ Y, mY ⇌ 2X+(m−1)Y always has the same overall X count: one consumed and two produced. Choose an integer m≥2 and 1<r<2, and set

$$
k_1^+=r/(r-1),\quad k_1^-=1/(r-1),\quad k_2^-=rm/(2-r),\quad k_2^+=k_2^-+1.
$$

The unique global maximum is still at x=y=1 with J=1, but its affinity satisfies

$$
e^{A^*}=r+(2-r)/m.
$$

Here r is the derivative of log y with respect to log x at the maximum. It measures the relative concentration response, not a reaction-rate constant.

![The optimal affinities differ at fixed overall X count and can approach zero across the family.](/publications/optimal-affinity/versions/v1/recycling_limits.png)

For each fixed m, the lowest limit of exp(A*) is 1+1/m, approached as r decreases toward 1. The endpoint is excluded. Along r=1+1/m the realized value is 1+2/m−1/m², which tends to 1; its logarithm therefore tends to zero. This limit changes m and the rate constants. At the displayed normalization J=1, some rates become large.

![Parameter map of affinity at the global production maximum, with one row per integer m.](/publications/optimal-affinity/versions/v1/affinity_parameter_map.png)

The map evaluates the exact formula on 29 integer values of m and 197 values of r. Color represents A*, not exp(A*). This is a map of a numerical quantity, not a phase diagram of different dynamics, and the model supplies no fractal structure. CSV files contain the plotted numbers.

## How the code constructs an actual local maximum

The generic `reconstruct` function accepts the reactant counts, product counts and a requested forward-response vector f. It normalizes that vector so that log x increases at unit speed, calculates the reverse responses h and the ratios qᵢ=hᵢ/fᵢ, and sets

$$
k_i^- = Jg_i/(q_i-1),\qquad k_i^+=q_i k_i^-.
$$

All monomials equal one at unit concentrations, so these rate constants are also the one-way rates at the constructed state. The code then checks positivity, stationarity, zero first derivatives, the rank of the current Jacobian, regularity of the controlled stationary equations, and the strictly negative second derivative of production.

For the main example the results are f=(1,3), h=(3/2,7/2), q=(3/2,7/6), and d²J/d(log x)²=−33/10. The routing weights give the same curvature as direct second differentiation of the stationary equations. An independent numerical solve at log x=±0.001 gives a finite-difference curvature of approximately −3.30000238; its error relative to −3.3 is about 2.38×10⁻⁶.

The general theorem requires the hypotheses listed in the demonstration plan and Theorem 5.15. Checking those hypotheses for one supplied profile does not prove that every profile in an arbitrary source's allowed set is accessible. Nor does the local construction establish global optimality without another argument.

## Additional exact examples and what they mean

| Example | exp(A_w*) | Gross X ratio | What is checked |
|---|---:|---:|---|
| Two-reaction counterexample | 7/4 | 2 | Unique global maximum, detailed balance, responses, rate construction |
| Recycling m=100, r=101/100 | 10199/10000 | 2 | Constructed strict local maximum; the general recycling certificate gives globality |
| Three-reaction circuit | 9/5 | 2 | Zero diagonal in T, detailed balance, exact global certificate |
| Power family, d=2 | 4 | 4/3 | Regular local maximum despite a routing coordinate that cannot be revisited; separate global gap |
| Power family, d=4 | 16/9 | 6/5 | Constructed local maximum and separate global gap |
| Additional interior-minimum teaching example | 3/2+√2 | 3 | Exact minimum of the response objective, matrix balance and local realization |

The additional teaching example is not quoted from the paper. It uses S₊=diag(1,2) and S₋=[[1,2],[1,1]]. If t=f₂/f₁>0, then

$$
\Phi(t)=\frac32+\frac2t+\frac t4,
\qquad
\Phi(t)-\left(\frac32+\sqrt2\right)=\frac{(t-2\sqrt2)^2}{4t}\ge0.
$$

The minimizing t=2√2 satisfies both response ratios greater than one. At this point the matrix-balancing condition Pᵀw=w holds exactly. This illustrates an attained interior minimum, whereas the recycling family has an unattained boundary infimum. The teaching example has N=2 and w=(1,1); A_w already includes those integer weights, so no additional factor N belongs in exp(A_w).

For Section 8's first core example, the supplied positive currents and potential differences are checked exactly. Joint production would require the summed current in each of two branches to lie between c₀ and 2c₀. The given barriers also force one branch sum to be ten times the other, which is impossible under those inequalities.

For the second core example, independently assigned activities (4,3,9/4,11/4) for (A,B,2A,3A) give the required currents. But actual species activities a,b would require a>b>a², hence 0<a<1. Joint production would additionally require a+2a²<3b<a+2a³, hence a>1. The contradiction shows why independently choosing a complex's activity is more permissive than requiring it to be the correct power of a shared species activity. The code checks these witness and contradiction identities; it does not reimplement the general core-minimality formalization.

## Exact calculations, numerical illustrations and proof scope

- SymPy checks matrix and polynomial identities exactly. Negative curvature and nonsingular matrices establish the local-maximum conditions on each checked example.
- For the recycling family, the program checks a coefficient-sum identity by its base and induction step, substitutes that formula into the general global-gap identity, and also checks literal finite sums for m=2,…,8. The accompanying positivity argument gives the global result for all integer m≥2 and 1<r<2.
- For the power family, the finite-sum identity has a checked base and induction step; the literal reaction elimination is checked for d=2,…,8. The paper supplies the general source-family argument.
- The plots and nearby branch solve use floating-point arithmetic. The largest steady-state residual in the main plotted curve is about 1.78×10⁻¹⁴. These numerical samples are not used as a proof of a universal theorem.
- The paper reports Lean verification. No Lean source was supplied or compiled in this run.
- The manuscript contains several inconsistent explanatory formulas. [Website conventions and manuscript notes](/publications/optimal-affinity/versions/v1/website-notes.md) identifies them and gives the conventions used by this script; the original PDF was not edited.

## Run the program

```sh
python -m pip install -r requirements.txt
python optimal_affinity.py --output outputs
```

No external data, network access or randomness is used during the run. The saved `diagnostics.json` records dependency versions, exact fixtures, numerical residuals and a SHA-256 hash of the Python source. `requirements.txt` pins the versions used for these results.

The transcript below is the actual website preparation run. The complete source below is read from the same checksum-verified file offered for download.


The pinned requirements match the versions recorded in the downloadable diagnostics.