Local energetic coupling enhances the expressivity of chemical computation

Marco Tuccio; Jason W. Rocks; Joshua E. Goldford

Source-linked AI summary

Local energetic coupling enhances the expressivity of chemical computation

Marco Tuccio, Jason W. Rocks, Joshua E. Goldford

arXiv:2609.17339v1Sep 15, 2026q-bio.MNcond-mat.stat-mechcs.ETphysics.chem-ph

TL;DR

The paper asks which physicochemical features determine computational expressivity in chemical reaction networks. It inverse-designs thermodynamically consistent networks through steady-state implicit differentiation and finds that nonequilibrium thermodynamic drive is the most effective single resource for steady-state computation.

Problem
The study addresses which physicochemical features control the computational expressivity of chemical systems.
Method
The authors inverse-design fixed-topology polymerization networks by implicitly differentiating through steady-state solutions while training standard chemical potentials, transition-state energies, and thermodynamic drives.
Results
Thermodynamic drives alone reach 60% training success versus 8% for standard chemical potentials, while the fully trainable control reaches 73%; expressivity correlates most strongly with reaction count (Spearman ρ = 1.00).
Takeaways & Limitations
Nonequilibrium drive is the most effective single resource for steady-state computational expressivity in chemical reaction networks.
Takeaways & Limitations
The parameter-freezing result may primarily reflect steady-state concentration targets; transient objectives could depend more strongly on kinetic parameters.

Abstract

from arXiv · show

Living systems compute with chemistry by mapping environmental signals onto specific internal chemical states. Despite recent advances in molecular programming, it remains unclear which physicochemical features control the computational expressivity of chemical systems. Here we inverse-design thermodynamically consistent chemical reaction networks whose steady-state response to an environmental input computes a target nonlinear function. Using implicit differentiation we train the free-energy landscape directly: standard chemical potentials, transition-state energies and thermodynamic drives. Increasingly large networks generated by elementary ligation and cleavage steps fit increasingly complex nonmonotonic polynomial functions, with expressivity scaling logarithmically with network size, predicted primarily by the number of reactions. Training individual energetic parameter classes reveals that internal thermodynamic drives, capable of breaking detailed balance, dominate trainability, with comparable performances achieved only by pairs of parameter classes. These results identify nonequilibrium drive as the most effective single resource for steady-state computational expressivity in chemical reaction networks.

I. SETUP

The paper frames chemical computation as matching an environmental concentration to a target steady-state concentration, using physically parameterized reversible CRNs in a CSTR. It trains energetic parameters through steady-state root finding and implicit differentiation, demonstrating fits to varied nonlinear responses.

CRN computing is formulated as realizing a target function y = f(x), where x is an environmental concentration and y is a chosen species’ steady-state concentration.
Polymerization Networks: Polymerization networks model reversible condensation–hydrolysis chemistry through ligation and cleavage reactions over monomer alphabets and bounded polymer lengths.
The networks use mass-action kinetics in a continuously stirred tank reactor, with monomer inflow and uniform dilution coupling environmental concentrations to internal chemistry.
Reaction rates are decomposed into transition-state energies, standard chemical potentials, and reaction-specific thermodynamic drives, giving each learned parameter an energetic interpretation.Transition-state energies control intrinsic speed, while thermodynamic drives bias reactions directionally and can represent energetic coupling such as ATP hydrolysis.
Training sweeps the input across a grid, minimizes mean squared error with Adam, and obtains exact steady-state gradients by implicit differentiation rather than backpropagating through stiff ODE integration.A stability regularizer biases solutions toward dynamically stable fixed points, and successful fits include diverse forms such as Gaussian bumps, step functions, and nonmonotonic targets.
Figure 2 evaluates success against target complexity using nine networks of up to ∼500 reactions, with success requiring R2 ≥0.98 and quasistatic dynamical verification.Input is the reservoir concentration of monomer a and output is the steady-state concentration of ab; H50 marks the complexity where success reaches 50%.

II. EXPRESSIVITY VERSUS NETWORK SIZE

Larger polymerization networks fit more complex nonmonotonic target functions, but topology-specific compatibility still matters. Expressivity correlates most strongly with reaction count, although the coupled topology set prevents isolating a unique structural determinant.

Some networks fit targets that others cannot: ab4 nearly fits a quadratic tent that ab3 fails to reproduce, showing dependence on topology–target compatibility.
The study compares nine polymerization topologies while fixing monomer a as input and polymer ab as output, using shared random squared trigonometric-polynomial targets across complexity classes.Each complexity class contains 50 targets, producing 4500 training runs; success requires R2 ≥0.98 plus quasistatic verification.
Expressivity tracks network size most strongly by reaction count, with exact ranking by R (Spearman ρ = 1.00), followed by trainable parameters (ρ = 0.97) and species count (ρ = 0.88).

III. FREEZING PARAMETER CLASSES

Freezing studies show that thermodynamic drives dominate trainability in a fixed polymerization topology. Adding other parameter classes helps little when drives are trainable, whereas the fully trainable model performs best.

The three parameter classes differ structurally: chemical potentials are shared per species, transition-state energies set bidirectional reaction speed, and drives are per-reaction and antisymmetric.Only thermodynamic drives can break detailed balance on an individual reaction.
Thermodynamic drives alone achieve a 60% success rate, compared with 8% for standard chemical potentials and no successes for transition-state energies.Success requires both R2 ≥0.98 and passing ODE quasistatic verification.
Adding a second parameter class to trainable drives provides no further benefit, while the fully trainable control reaches the highest success rate of 73%.

IV. DISCUSSION

The study shows that physically parameterized polymerization networks can compute diverse steady-state responses, with expressivity scaling logarithmically with network size and relying most strongly on nonequilibrium drive.

IV. DISCUSSION: Physically parameterized, fixed-topology polymerization networks reproduce diverse smooth steady-state responses without catalysis or autocatalysis by adjusting reaction energies.The networks are trained through implicit differentiation of steady-state solutions, avoiding backpropagation through stiff ODE integration.
IV. DISCUSSION: The capacity to fit increasingly complex responses scales logarithmically with network size, although specific structural causes remain unresolved because topology attributes covary with size.Different network-generating primitives could help decouple the topological features responsible for expressivity.
IV. DISCUSSION: Nonequilibrium thermodynamic drive is the most effective single resource for training networks to compute steady-state responses.Freezing studies indicate that driving reactions out of equilibrium matters more than fine-tuning individual transition-state energies.
IV. DISCUSSION: The conclusions primarily concern steady-state concentration targets; transient time-dependent objectives may depend more strongly on kinetic parameters.This scope boundary qualifies how broadly the dominance of thermodynamic drive should be generalized.
IV. DISCUSSION: The framework may support experimental CRN design, including biosensing and computation using enzyme-mediated amino-acid ligation coupled to ATP hydrolysis.The authors identify enzyme-family and substrate-dependent kinetic diversity as a possible route to tuning dissipation and barrier heights.
IV. DISCUSSION: The reaction model assigns standard chemical potentials to species, transition-state energies to reactions, and thermodynamic drives with opposite effects on forward and reverse barriers.Transition-state energies alter reaction timescales, while chemical potentials and thermodynamic drives determine directional bias; nonzero drives can break detailed balance.

Appendix B: Steady-State Training and Dynamical Verification

Training evaluates steady-state responses across input points, updates physicochemical parameters iteratively, and regularizes solutions for numerical accuracy, stability, and conditioning.

Appendix B: Steady-State Training and Dynamical Verification: Training repeatedly evaluates the steady-state response on the input grid, differentiates the fitting objective with respect to physicochemical parameters, and updates them.After each update, the network reaches a new steady state before the next training loop.
Appendix B: Steady-State Training and Dynamical Verification: Regularization favors steady states that are accurately resolved, dynamically stable, and numerically well conditioned.

1. Training Objective

The training objective matches the readout species’ steady-state concentration to a target function across a grid of environmental inputs and updates parameters with Adam.

1. Training Objective: For each environmental input x_k, the steady-state readout concentration y(x_k) is compared with the target function f.
1. Training Objective: The fitting loss and regularization terms update the parameter vector θ using the Adam optimizer.

a. Regularization terms

Training uses residual, spectral, and conditioning regularization to obtain accurate, stable, nonsingular steady states, while quasistatic verification separately tests dynamical reachability.

a. Regularization terms: Three regularizers target low-residual, dynamically stable, and well-conditioned steady-state solutions.Residual penalties improve root quality, spectral penalties discourage unstable Jacobian eigenvalues, and conditioning penalties address unreliable root finding and implicit differentiation near singular Jacobians.
a. Regularization terms: The spectral penalty does not guarantee that the response lies on one dynamically accessible steady-state branch, so quasistatic verification is required after training.Verification follows the ODE dynamics along the ordered input grid rather than checking independently solved roots alone.
a. Regularization terms: Implicit differentiation computes steady-state parameter gradients without storing and backpropagating through the full root-finder trajectory.Adjoint solves are performed independently at each input point and combined across the input grid.
a. Regularization terms: Training discards updates when more than 50% of input points have untrusted roots and terminates at numerical, convergence, trust, or epoch-budget limits.A run's exit condition alone does not determine success; fit quality and post-training dynamical verification are also required.
a. Regularization terms: Quasistatic verification follows converged ODE states across the input grid and rejects any run with a failed grid point or mismatched readout.This excludes unstable roots, disconnected steady-state branches, and trajectories that do not settle to the predicted fixed points.
a. Regularization terms: The target ensemble contains smooth, strictly positive profiles whose complexity knob H controls oscillations, with outputs rescaled to [0.2, 0.6].Squared harmonic sums ensure nonnegative concentration-compatible targets, while increasing H permits more peaks, valleys, and curvature changes.

Appendix D: Illustrative Fits

The framework fits diverse smooth one-dimensional steady-state responses, including localized peaks, sharp transitions, oscillatory targets, and a two-dimensional parametric curve.

Appendix D: Illustrative Fits: A network with N = 258 species and R = 462 reactions fits a highly oscillatory H = 5 trigonometric target during training.The response is shown at intermediate epochs 800, 8000, and 12000 for the abcdef3 topology, using monomer a as input and dimer ab as output.
Appendix D: Illustrative Fits: Illustrative fits include Gaussian bumps, step functions, random Legendre polynomials, and random trigonometric polynomials with six or eight harmonics.These examples span localized peaks, sharp transitions, and oscillatory response shapes across ab4, ab5, and ab6 topologies.
Appendix D: Illustrative Fits: The procedure also reproduces a von Neumann elephant parametric curve using one input and two readout species.The ab6 network uses readouts aba and bab, with the dimer ab as input and 64 sampled input values.
Appendix D: Illustrative Fits: Comparing ab3 and ab4 on the same quadratic tent target isolates the effect of increasing maximum polymer length.The input species a and output species ab are held fixed in this comparison.

1. Multiple-Readout Training

Multiple-readout training extends the framework from scalar responses to jointly fitted coordinates, while the topology set varies alphabet size and polymer length to change network size.

1. Multiple-Readout Training: One input can train two readout species jointly against two target coordinates, enabling a parametric response rather than a single scalar map.The two readout concentrations at each input-grid point are optimized together with a joint loss.
1. Topology Set and Network-Size Measures: The topology family is generated by sweeping alphabet size A and maximum polymer length L, producing networks with specified species and reversible-reaction counts.The study includes fixed-A length sweeps, fixed-L alphabet sweeps, and an intermediate A = 3, L = 4 configuration.
1. Topology Set and Network-Size Measures: Each topology has N standard chemical-potential parameters, 2R energetic reaction parameters, and one flow-rate parameter, giving P = N + 2R + 1 trainable parameters.The two reaction parameters are a transition-state energy and a thermodynamic drive for each reaction.
1. Topology Set and Network-Size Measures: Because changing A or L increases N, R, and P together, the experiment cannot identify one size measure as the unique structural determinant of expressivity.The available polymerization topologies therefore strongly covary across all three measures.

2. Experimental Design

The scaling study evaluates common-input, common-readout polymerization networks across random target complexities using verified fit success, while separately tracking optimization and dynamical failures.

2. Experimental Design: Quasistatic verification is necessary because independently root-found responses may use unstable points or splice coexisting branches.A failed grid point discards the run, separating optimization failures from responses that nonlinear ODE dynamics cannot follow.
2. Experimental Design: Readout choice affects fit quality: varying all 80 non-monomer readouts in abcd3 shows that topology alone does not determine performance.The scaling study fixes ab as the readout to remove this variation, but does not vary the input species with the readout fixed.

Appendix F: Parameter-Freezing Study

The parameter-freezing study evaluates how individual energetic parameter classes contribute to CRN trainability under fixed topology, input, readout, and flow conditions. It compares fit quality across parameter subsets using shared H = 2 targets and quasistatic verification.

Experimental design: The study varies energetic parameter classes while holding topology, target, input, readout, and flow rate fixed.Excluded classes are initialized to zero and frozen; the flow rate is fixed at γ = 1.
Readout dependence: Figure 5 measures readout-dependent fit quality in abcd3 using the base-10 logarithm of the best loss for each run and readout choice.The figure identifies a non-trivial readout-species dependence that the study does not further investigate.
Fit-quality comparison: Figure 6 illustrates how trained response curves differ across parameter-freezing bins for a common target, complementing the aggregate success rates.The examples are visualization only and are not additional trials beyond the 128 targets used for success rates.

Appendix G: Implementation and Reproducibility

The implementation uses 64-bit JAX and Diffrax computations with GPU-parallelized independent runs. Reproducibility is supported by separate random streams, deterministic seed derivation, and shared target draws across comparisons.

Implementation: CRN dynamics, steady-state solving, implicit gradients, optimization, and ODE verification are implemented in JAX and Diffrax using 64-bit precision.Scaling and parameter-freezing studies run on T4 NVIDIA GPUs with independent runs parallelized across targets and parameter bins.
Reproducibility: Target generation and parameter initialization use separate random streams so either process can be reproduced independently.Randomness is used for both target functions and trainable-parameter initialization.
Reproducibility: Run seeds are derived deterministically from experiment-grid positions, with coprime offsets and strides preventing collisions.Both target and initialization seeds are recorded in released configuration files.
Experimental comparability: Target draws are shared across topologies and parameter bins, while initialization seeds differ to provide independent starting points.This makes architecture and parameter-subset comparisons use identical target ensembles.