How to construct and test a reduced biochemical model

A standalone worked chapter: exact totals, singular limits, implicit binding constraints, regulatory laws and stochastic averaging.

Lecture 6 · lecturer-written exposition · Companion to the core

This chapter teaches a procedure you can apply to a new reaction network. Start with a specified prediction. Build the full model, identify exact totals, justify the fast-state elimination and test what the smaller model preserves. You need differentiation, elementary linear algebra and the interpretation of a reaction as a state change with a rate. The stochastic sections introduce the generator when it is needed.

Three operations carry different assumptions. Molecular accounting gives exact total balances. A specified singular limit supplies an approximate fast-state closure. Additional dominance assumptions may turn that implicit closure into an explicit Michaelis–Menten or Hill law. Each worked example keeps these operations separate.

The closed enzyme supplies the full calculation first. Supply and drain then expose the role of boundaries. Product competition and self-repression require coupled binding constraints. A binary promoter requires a probability distribution in place of a deterministic fast concentration. The final enzyme–buffer exercise combines the method without assuming attendance at the lecture. All numerical examples in this chapter are model calculations.

1Specify what the model must predict

Suppose the experiment starts with enzyme and substrate mixed in a closed vessel. A detector measures product over several catalytic times. A different detector measures complex formation immediately after mixing. The same reduced law can be excellent for the first measurement and unsuitable for the second. Before doing algebra, name the output, time window, initial conditions, and acceptable error.

For the running calculation we will compare product trajectories and complex reconstruction separately. Concentrations refer to a well-mixed vessel of fixed volume. Enzyme total is constant. The full model explicitly resolves reversible binding and irreversible conversion. It omits product rebinding, enzyme synthesis, compartment transport, and growth. Those omissions define a model, not a statement that the omitted processes never occur.

Write a reduction target that a concrete measurement or simulation could fail.

2Give every channel a flux and a change

Let $E$ be free enzyme, $S$ free substrate, $C_{ES}$ their complex and $P$ free product. Order the species as $x=(E,S,C_{ES},P)^\top$. Here the superscript denotes transpose. The mechanism and channel rates are

Association changes species by $(-1,-1,1,0)$. Dissociation reverses that vector. Conversion changes them by $(1,0,-1,1)$. Multiply each vector by its flux and add. With $v_\mathrm{bind}=v_\mathrm{bind}^{+}-v_\mathrm{bind}^{-}$, this gives

If concentrations are micromolar and time seconds, $k_+$ has units $\mu\mathrm{M}^{-1}\mathrm{s}^{-1}$. The constants $k_-,k_\mathrm{cat}$ have units $\mathrm{s}^{-1}$. All three fluxes have units $\mu\mathrm{M}\,\mathrm{s}^{-1}$. This dimensional check catches comparisons of constants that are not comparable rates.

Construct and dimension-check a full mass-action ODE from state-change vectors.

3Establish the physical region

At a boundary where free substrate is zero, its derivative is $k_-C_{ES}\geq0$. At zero free enzyme, its derivative is $(k_-+k_\mathrm{cat})C_{ES}\geq0$. At zero complex, its derivative is $k_+ES\geq0$. Product cannot decrease in the closed model. Nonnegative initial concentrations therefore remain nonnegative.

Adding the enzyme and complex equations gives $d(E+C_{ES})/dt=0$. Adding substrate, complex, and product gives $d(S+C_{ES}+P)/dt=0$. With fixed initial totals, every species is bounded. These identities give both a physical domain and numerical diagnostics.

Check positivity and conservation in the full model and use them as independent tests of a numerical trajectory.

4Find totals with a counting matrix

Define a row for each constituent: enzyme occurs in free enzyme and complex. Unconverted substrate occurs in free substrate and complex. Product occurs in free product alone. Thus

For the binding column $\gamma_\mathrm{bind}=(-1,-1,1,0)^\top$, direct multiplication gives $L\gamma_\mathrm{bind}=0$. For conversion, $L\gamma_\mathrm{cat}=(0,-1,1)^\top$. Consequently

These are identities at any binding speed. A total is a subnetwork invariant, not necessarily a whole-network conserved quantity. Catalysis changes unconverted substrate into product while preserving enzyme abundance.

For a larger network, the retained rows must describe all independent fast-subnetwork invariants needed to reconstruct its state. Physically named constituent totals are useful rows, but one must check their independence and completeness. Section 27 returns to this coordinate requirement.

Build the exact projection that removes fast reaction fluxes from total dynamics.

5Change coordinates before taking a limit

Retain $q_S,C_{ES},P$ and treat $q_E$ as a fixed parameter. Reconstruction is $E=q_E-C_{ES}$, $S=q_S-C_{ES}$. The full system becomes

The physical interval for complex is $0\leq C_{ES}\leq\min(q_E,q_S)$. Binding alone can change complex rapidly but cannot change the total coordinate. This is exactly what a clean slow coordinate should do.

Rewrite a mechanism in slow totals and fast occupancy without dropping any full-system term.

6State the general singular perturbation problem

A fast equation can become an algebraic constraint only after a scaling makes the approximation precise. Choose dimensionless state vectors $x^\mathrm{slow}$ and $x^\mathrm{fast}$, with physical reference times $t_\mathrm{slow}$ and $t_\mathrm{fast}$. Define $\varepsilon=t_\mathrm{fast}/t_\mathrm{slow}$. The standard form is

The rate functions $f,g$ are dimensionless in this display. They must remain finite as the selected parameter family sends $\varepsilon$ to zero. State scales are part of the construction. A complex that tends to zero may still have a finite, rapidly changing fractional occupancy.

The fast problem determines which algebraic root the system can follow. Changing the clock to $t/t_\mathrm{fast}$ multiplies the slow rate by $\varepsilon$. The leading fast equations are therefore

At each retained slow state, solve $g(x^\mathrm{slow},x^\mathrm{fast};0)=0$ for a physical root $x^\mathrm{fast}=h(x^\mathrm{slow})$. This graph is the critical manifold. Check that the fast dynamics attracts the intended initial fast state. For a locally attracting smooth branch, a fast Jacobian whose eigenvalues have real parts uniformly below zero supplies the local stability condition.

The slow problem evolves the retained state with the reconstructed fast state:

The reduction applies after the initial layer on a bounded slow-time interval that stays inside the region of attraction. It also needs smooth rate functions, a selected isolated branch and appropriate initial data. These hypotheses connect an algebraic solve to a trajectory approximation.³

The reconstructed fast state keeps moving because the slow state moves. Its derivative is $Dh\,\dot x^\mathrm{slow}$. The term neglected in the fast equation is $\varepsilon t_\mathrm{slow}\dot x^\mathrm{fast}$, rather than a requirement that this derivative vanish on the slow clock. The next four sections calculate every part of this construction for the enzyme.

Specify a singular limit, select its fast branch, and state the conditions needed to use that branch in a slow model.

7Scale the enzyme without renaming its species

For the enzyme, use the fixed concentration scale $K_d=k_-/k_+>0$. Define $\bar q_S=q_S/K_d$, $\bar C_{ES}=C_{ES}/K_d$, $\bar q_E=q_E/K_d$, and $\bar q_P=q_P/K_d$. Physical time remains $t$. The reference times are $t_\mathrm{bind}=1/k_-$ and $t_\mathrm{cat}=1/k_\mathrm{cat}$.

The association term supplies the conversion factor. Since $k_+K_d=k_-$, it becomes $k_-K_d(\bar q_E-\bar C_{ES})(\bar q_S-\bar C_{ES})$. The complex derivative is $K_d\,d\bar C_{ES}/dt$. Dividing the full complex equation by $k_-K_d$ gives

The totals form $x^\mathrm{slow}$, and the complex is $x^\mathrm{fast}$. To shrink the time-scale ratio at fixed affinity, increase $k_+$ and $k_-$ together. Hold catalytic rate, initial totals and retained slow boundary rates fixed. The Michaelis constant $K_M=(k_-+k_\mathrm{cat})/k_+$ then satisfies $K_M/K_d=1+\varepsilon$.

Recover the small parameter by carrying units and derivative transformations through every term.

8Select the physical root and test attraction

The reconstruction $h_{C_{ES}}(q_S;q_E,K)$ returns complex concentration from substrate total, enzyme total and a positive balance constant. Solve $(q_E-C_{ES})(q_S-C_{ES})=KC_{ES}$ and select the smaller root:

There is exactly one root in the physical interval. At $C_{ES}=0$, the quadratic is nonnegative. At $C_{ES}=\min(q_E,q_S)$, it is nonpositive. Its derivative throughout this interval is $-(q_E+q_S+K-2C_{ES})<0$. The larger root exceeds both totals.

The two roots remain distinct for nonnegative totals and positive $K$. The discriminant can be written as $(q_E-q_S)^2+2K(q_E+q_S)+K^2>0$.

Attraction requires a separate dynamical check. For rapid binding, take $K=K_d$. Differentiate the net binding flux at fixed totals and evaluate it at $C_{ES}=h_{C_{ES}}$:

The quantity in parentheses is the sum of free enzyme, free substrate and $K_d$. A small complex perturbation therefore decays as $e^{-t/t_\mathrm{relax}}$ when totals are fixed. In affinity units, $t_\mathrm{relax}/t_\mathrm{bind}=K_d/(q_E+q_S+K_d-2h_{C_{ES}})$.

For smooth systems on suitable compact attracting portions of the critical manifold, a uniform negative fast eigenvalue supports a nearby locally invariant slow manifold at small positive $\varepsilon$.³ If a transverse eigenvalue approaches zero in another model, the same argument can fail. Solving an algebraic equation alone does not establish fast attraction.

Check root feasibility and calculate the physical relaxation time needed to justify fast-state reconstruction.

9Solve the initial binding layer

During the leading fast problem, hold substrate total at $q_S(0)$ and neglect catalysis. The physical complex equation is

Use the root just derived. Let $h_{C_{ES}}=h_{C_{ES}}(q_S(0);q_E,K_d)$ denote the physical binding-equilibrium root. The second root is $C_{ES}^{+}=q_E+q_S(0)+K_d-h_{C_{ES}}$. It exceeds the available totals and is not a physical concentration. Factoring the quadratic gives $\dot C_{ES}=k_+(C_{ES}-h_{C_{ES}})(C_{ES}-C_{ES}^{+})$.

For initial complex zero, separation of variables gives an explicit binding transient. Define its local relaxation time $t_\mathrm{relax}=1/[k_+(C_{ES}^{+}-h_{C_{ES}})]$. The partial-fraction identity is

Integrate both sides, multiply by the root difference and impose zero initial complex. Exponentiating gives a ratio of distances from the two roots. Solving that relation for the complex gives

The formula starts at zero and approaches $h_{C_{ES}}$. For $q_E=K_d$ and $q_S(0)=4K_d$, the root is approximately $0.763932K_d$, and $t_\mathrm{relax}=t_\mathrm{bind}/\sqrt{20}$. The reference unbinding time and the actual binding relaxation time are distinct.

The full substrate total changes slightly during this layer because catalysis continues. In affinity units, $d\bar q_S/d(t/t_\mathrm{bind})=-\varepsilon\bar C_{ES}$. Its change over a bounded binding-time interval is of order $\varepsilon$. Only the binding contribution to its derivative cancels exactly.

Compute an initial binding transient and distinguish exact fast-subnetwork conservation from small full-system drift.

10Reconstruct a moving fast state

The reconstruction changes when the substrate total changes. Hold $q_E$ and $K_d$ fixed, and write $h_{C_{ES}}(q_S)=h_{C_{ES}}(q_S;q_E,K_d)$. Implicit differentiation of its binding relation gives

The leading closed-vessel dynamics is $\dot q_S=-k_\mathrm{cat}h_{C_{ES}}$. Its reconstructed complex has derivative $d h_{C_{ES}}(q_S(t))/dt=-k_\mathrm{cat}h_{C_{ES}}\,dh_{C_{ES}}/dq_S<0$ while substrate remains. A fast constraint therefore predicts a changing complex. The nonzero derivative is computed by differentiating the reconstruction along the slow solution.

A more accurate reconstruction must follow an invariant graph of the full finite-speed dynamics. Write this graph as $C_{ES}=h_{C_{ES}}^{\varepsilon}(q_S)$. Its derivative along the total equation must equal its derivative in the full complex equation. Dividing that equality by $k_+$, and using $k_\mathrm{cat}/k_+=\varepsilon K_d$, gives

To find the first displacement from binding equilibrium, expand $h_{C_{ES}}^{\varepsilon}=h_{C_{ES}}+\varepsilon h_{C_{ES}}^{(1)}+O(\varepsilon^2)$. The leading term satisfies binding balance. At first order, differentiating the binding polynomial contributes $-(q_E+q_S+K_d-2h_{C_{ES}})h_{C_{ES}}^{(1)}$. The invariance equation becomes

Solving for the correction gives

The curve where the full complex derivative is zero is a different object, called its nullcline. It uses $K_M=K_d(1+\varepsilon)$ in the quadratic. Expanding that root gives a first correction of $-K_dh_{C_{ES}}/(q_E+q_S+K_d-2h_{C_{ES}})$. The invariant graph also contains the chain-rule contribution. Its small offset from the nullcline supplies the slow derivative.⁴

Derive the invariance equation and its first correction instead of replacing a moving reconstruction by a zero-derivative curve.

11Match the lost initial condition

For a separate two-variable example, let $\tau=t/t_\mathrm{slow}$ and let primes denote differentiation with respect to this dimensionless time. Consider $x'=-x,\varepsilon y'=x-y$, with $x(0)=1,y(0)=0$. First solve $x=e^{-\tau}$. Multiplying the second equation by its integrating factor $e^{\tau/\varepsilon}$ gives

Integrate from zero and apply the initial condition:

The exact slow graph is $y=x/(1-\varepsilon)$. Its expansion begins $x+\varepsilon x$. The initial-layer term restores the original value zero. To suppress an order-one fast mismatch to order epsilon, solve $e^{-\tau/\varepsilon}\sim\varepsilon$, giving $\tau\sim\varepsilon|\log\varepsilon|$. An e-folding time is shorter than this matching time.

Construct a matched solution that retains the fast initial condition while approaching a moving slow graph.

12Integrate the slow enzyme conversion

Rapid equilibrium gives $C_{ES}=q_ES/(K_d+S)$. Since $q_S=S+C_{ES}$, differentiation produces

Let $S_a$ be the free substrate after the leading binding projection, and set that reduced initial time to zero. Separate variables. The integrand reduces to $1+(K_d+q_E)/S-q_E/(K_d+S)$, divided by $k_\mathrm{cat}q_E$. Thus

This implicit solution predicts the full reduced progress curve. For $q_E=K_d,q_S(0)=4K_d$, the projection gives $S_a/K_d=4-0.763932=3.236068$. Using four as the projected free value would ignore initial sequestration.

Near depletion the total-coordinate root is $C_{ES}\approx q_Eq_S/(K_d+q_E)$, so total substrate has an exponential tail with rate $k_\mathrm{cat}q_E/(K_d+q_E)$. Depletion does not itself create a fold or invalidate this rapid-binding limit.

Integrate a closed reduced trajectory and initialize it with the correct binding projection.

13Derive the small-enzyme criterion

If free substrate changes little while complex forms, freeze it temporarily at $S_\mathrm{ref}$. The complex equation becomes linear with relaxation rate $k_+(S_\mathrm{ref}+K_M)$ and plateau $C_{ES}^{\mathrm{plateau}}=q_ES_\mathrm{ref}/(S_\mathrm{ref}+K_M)$, where $K_M=(k_-+k_\mathrm{cat})/k_+$.

The amount initially captured relative to substrate is $C_{ES}^{\mathrm{plateau}}/S_\mathrm{ref}=q_E/(S_\mathrm{ref}+K_M)$. Requiring this to be small gives the familiar initial reactant-stationarity criterion.

The catalytic depletion time estimates how long the initial substrate would last at the plateau conversion rate. Compare this with the complex relaxation time:¹

The ratio $t_\mathrm{complex}/t_\mathrm{substrate}=k_\mathrm{cat}q_E/[k_+(S_\mathrm{ref}+K_M)^2]$ can be small without guaranteeing negligible binding sequestration. Conversely, large affinity-scale denominator can make the initial criterion small even without a tiny enzyme/substrate ratio. These estimates apply under the stated hypotheses and over the initial adjustment interval. They do not supply a uniform error bound for the full trajectory.

Derive the small-enzyme criterion from initial sequestration and distinguish it from a bare timescale ratio.

14Restore supply and drain in both limits

Substrate supply adds $v_\mathrm{in}$ to the free-substrate equation. Product drain subtracts $k_\mathrm{out}P$ from the product equation. Binding is described by $v_\mathrm{bind}=k_+ES-k_-C_{ES}$. The full equations are now $\dot E=-v_\mathrm{bind}+k_\mathrm{cat}C_{ES}$, $\dot S=v_\mathrm{in}-v_\mathrm{bind}$, $\dot C_{ES}=v_\mathrm{bind}-k_\mathrm{cat}C_{ES}$ and $\dot P=k_\mathrm{cat}C_{ES}-k_\mathrm{out}P$. Adding free and bound contributions gives the total balances below.

Accelerate binding at fixed affinity

For rapid binding, the normalized substrate equation is $t_\mathrm{cat}\dot{\bar q}_S=t_\mathrm{cat}v_\mathrm{in}/K_d-\bar C_{ES}$. The complex equation is unchanged from Section 7. Keep $t_\mathrm{cat}v_\mathrm{in}/K_d$ finite as $\varepsilon_\mathrm{bind}=t_\mathrm{bind}/t_\mathrm{cat}\to0$. The total then changes negligibly during binding, and the fast constraint uses $K_d$.

Slow substrate change by reducing enzyme

For small-enzyme QSSA, the fast state is better scaled by enzyme availability. Keep the ratio $C_{ES}/q_E$ explicit rather than giving it a new letter. Fix a positive substrate reference $S_\mathrm{ref}$. Define $t_\mathrm{complex}=1/[k_+(K_M+S_\mathrm{ref})]$ and $t_\mathrm{substrate}=(K_M+S_\mathrm{ref})/(k_\mathrm{cat}q_E)$. The exact equations become

These equations expose a second singular perturbation family. Let $q_E\to0$ at fixed microscopic constants and substrate scale, while keeping $v_\mathrm{in}/(k_\mathrm{cat}q_E)$ finite. Then $t_\mathrm{complex}/t_\mathrm{substrate}\to0$. The fast equation is attracting, and its stationary occupancy is $C_{ES}/q_E=S/(K_M+S)$. Also $S=q_S+O(q_E)$ on this family.

The fast fraction can remain finite while the absolute complex tends to zero. This is why concentration scale matters when identifying slow and fast variables. In this family the ratio $\varepsilon_\mathrm{QSSA}=q_E/(K_M+S_\mathrm{ref})$ controls sequestration, and

Retaining $S=q_S-C_{ES}$ in the algebraic balance gives $C_{ES}\approx h_{C_{ES}}(q_S;q_E,K_M)$. This total-QSSA reconstruction agrees with standard QSSA at leading order in the stated small-enzyme family. It does not automatically acquire a validity theorem at arbitrary enzyme concentration.

Check the boundary processes on the fast interval

Boundary forcing must respect the selected limit. Check the supplied amount during fast adjustment against the substrate scale: $v_\mathrm{in}t_\mathrm{bind}/K_d\ll1$ in the first family, or $v_\mathrm{in}t_\mathrm{complex}/S_\mathrm{ref}\ll1$ in the second. Fixed finite input during a limit that sends enzyme capacity to zero does not satisfy the latter family.⁵

Product removal can be retained in its exact downstream equation $\dot q_P=k_\mathrm{cat}C_{ES}-k_\mathrm{out}q_P$. If product is counted among the slow variables, its drain must be slow on the fast-adjustment interval. If drain is faster, product has its own faster response. Because product does not feed back in this model, that response does not alter the substrate–complex reduction.

Put both open-enzyme limits into singular perturbation form and test their state scales, source scaling and fast attraction.

15Separate total closure from Michaelis–Menten dominance

In total coordinates the full complex nullcline satisfies $(q_E-C_{ES})(q_S-C_{ES})=K_MC_{ES}$. Its physical root is $h_{C_{ES}}(q_S;q_E,K_M)$. Rationalizing the numerator avoids subtracting nearly equal numbers:

Integrate $\dot q_S=-k_\mathrm{cat}h_{C_{ES}}(q_S;q_E,K_M)$ and reconstruct free species. Replacing $K_M$ by $K_d$ produces the rapid-equilibrium closure. Its different balance requires different assumptions.

The explicit Michaelis–Menten law makes an extra substitution

The algebraic relation gives $C_{ES}=q_ES/(K+S)$ in free substrate, with the constant selected by the justified limit. A law in total substrate replaces $S=q_S-C_{ES}$ by $q_S$. The bound fraction must be small for this replacement:

Thus either $q_E\ll K$ or $q_E\ll q_S$ suffices. The resulting explicit flux is $k_\mathrm{cat}q_Eq_S/(K+q_S)$. The inequalities must hold over the concentration range being predicted. A large initial substrate pool need not remain large during depletion.

Fast binding alone does not enforce either dominance condition. With $q_E=q_S=K$, the physical root is $0.381966K$, whereas the free-to-total substitution gives $0.5K$. The difference survives when both binding directions accelerate together. In the tight-binding limit, complex approaches $\min(q_E,q_S)$, so a substantial part of the substrate total can remain bound.

Measure the two approximation errors separately

Integrate three models from the same substrate total: the full species dynamics, the total DAE and the explicit Michaelis–Menten law. Let their reconstructed fractional occupancies be $C^\mathrm{full}_{ES}/q_E$, $C^\mathrm{DAE}_{ES}/q_E$ and $C^\mathrm{MM}_{ES}/q_E$. Define the error between any two models on the stated interval as

The full-to-DAE comparison tests finite-speed elimination. The DAE-to-MM comparison tests the extra dominance substitution and its accumulated effect on the total trajectory. The full-to-MM comparison includes both. These maximum errors are not additive because their peaks and signs can differ.

Both families start at $q_S=1\,\mu\mathrm M$. Supply changes from $0.3k_\mathrm{cat}q_E$ to $0.8k_\mathrm{cat}q_E$ at three reference times. Full-model complex starts at the chosen total-closure root. In the rapid-binding family, $q_E=K_d=1\,\mu\mathrm M$, $k_\mathrm{cat}=0.5\,\mathrm{s}^{-1}$, $k_-=k_\mathrm{cat}/\varepsilon_\mathrm{bind}$ and $t_\mathrm{ref}=t_\mathrm{cat}$. The DAE uses $K_d$.

In the small-enzyme family, $k_+=10\,\mu\mathrm M^{-1}\mathrm{s}^{-1}$, $k_-=1\,\mathrm{s}^{-1}$, $k_\mathrm{cat}=9\,\mathrm{s}^{-1}$ and $S_\mathrm{ref}=1\,\mu\mathrm M$. Thus $K_M=1\,\mu\mathrm M$ and $q_E=2\varepsilon_\mathrm{QSSA}\,\mu\mathrm M$. The reference is $t_\mathrm{substrate}$, and the DAE uses $K_M$.

The first row gives a useful counterexample: accurate fast-state elimination can coexist with a poor explicit law in total substrate. The second row shows the joint improvement used in core Figure 6. Neither row is a uniform theorem for arbitrary forcing or initial conditions.

Choose a physically justified closure and test free-to-total substitution as a separate approximation.

16Calculate capacity and recovery

At a full steady state, total balance requires $C_{ES}^*=v_\mathrm{in}/k_\mathrm{cat}$. Here $C_{ES}^*$ is the steady complex concentration. The complex equation then gives $S^*=K_MC_{ES}^*/(q_E-C_{ES}^*)$. Product balance gives $P^*=v_\mathrm{in}/k_\mathrm{out}$. A finite positive state needs $0<v_\mathrm{in}<k_\mathrm{cat}q_E$ and positive drain.

For a reduced increasing flux $v(q_S)=k_\mathrm{cat}h_{C_{ES}}(q_S;q_E,K)$, define perturbations $\xi=q_S-q_S^*$, $\eta=q_P-q_P^*$. Linearization is explicitly

The triangular matrix has negative eigenvalues, so both modes decay locally. Close to enzyme saturation, $v'$ is small: the slow recovery clock lengthens even while binding remains fast. Above capacity, the exact inequality $\dot q_S\geq v_\mathrm{in}-k_\mathrm{cat}q_E>0$ proves accumulation without any closure.

Find the processing steady state, derive its linearized system, and distinguish slow stability from an exact overload bound.

17Choose a numerical test that can reject the reduction

The reproducible calculation integrates the full total-coordinate enzyme and two reduced ODEs. RK4 steps are chosen small relative to the fastest binding relaxation and hit output times exactly. No negative-state clipping is used. A second run at half the step tests discretization. Conserved totals test the equations independently.

Separate graph reconstruction from trajectory error. The former compares full complex with $h_{C_{ES}}(q_S^{\mathrm{full}};q_E,K_d)$ at the actual total. The latter compares independent full and reduced total trajectories. A closure can have a small instantaneous graph error and still accumulate a noticeable trajectory discrepancy over a long time.

For stiff large networks, use a tested implicit solver with tolerances checked against invariants. The explicit solver used for this bounded example resolves the fast relaxation. The coupled product-competition calculation below uses an implicit solver and an independent species-coordinate implementation.

Separate discretization error, closure error, and accumulated trajectory error in a reproducible comparison.

18Work the product-competition DAE from the full model

A product-bound enzyme adds a coupled fast constraint. Let free enzyme $E$ bind substrate $S$ or product $P$, forming $C_{ES}$ or $C_{EP}$. The complexes interconvert through slower reversible catalysis. Substrate enters at flux $v_\mathrm{in}$. Only free product leaves, at flux $k_\mathrm{out}P$.

Define the three net fluxes, then apply each stoichiometric change: