From a molecular state list to a predictive model

A calculation chapter on Markov dynamics, detailed balance, physical weights and fast biological regulation. CCBS · Fall 2026 · Lecture 7 · Lecturer-written exposition · Revised 22 September 2026 The task is to turn a biological description into a calculation you can defend. A promoter switches, an enzyme changes conformation, or a receptor binds a ligand. Which states must be distinguished? What sets the waiting time? When can equilibrium weights replace a kinetic solve? What measurement would reveal that the replacement failed? This chapter builds that competence in four steps. First construct a jump model and calculate its dynamics. Then test whether its stationary distribution is an equilibrium distribution. Next use physical state weights to derive regulatory responses. Finally restore finite pools, changing inputs and unresolved mechanisms. A gated receptor supplies a complete worked test. Optional applications develop the tools further. A historical workshop compares these methods on the original Michaelis–Menten, Hill and MWC examples. The mathematics assumes basic differentiation, integration and probability. Matrix notation is introduced where it is needed. Numerical examples are illustrative unless explicitly identified as reconstructions using published model parameters. Neither kind is presented as new experimental data. The core selects the classroom argument. This chapter supplies the methods needed to repeat and extend it independently.

Part I

Build and solve a molecular jump model

Learn to predict occupancies, count fluctuations and arrival times from a specified set of states and transition hazards.

1Choose a state that retains the needed information

A state is a predictive description, not merely a measured variable. In an autonomous deterministic model, a state

x(t)

obeys

\dot x=f(x)

. Given that state and a well-posed rate function, the future trajectory is fixed. For a spring, position alone is insufficient. The same position with opposite velocities leads to different motion. Position and velocity together form the state.

A stochastic state determines a future probability law. Write $\xi(t)$ for a random molecular configuration. The Markov assumption says that, conditional on the present state and the specified environment, earlier history adds no predictive information. One promoter might have states OFF and ON. A molecule-count model instead has states $N_X(t)=0,1,2,\ldots$ . Random trajectories differ even when they start in the same state.

Constant transition hazards require a fixed environment in addition to the Markov assumption. If an externally controlled ligand changes with time, a Markov model can have time-dependent hazards. This chapter first holds external conditions fixed, giving a time-homogeneous continuous-time Markov chain. Later we let the conditions change and test whether the fast state follows them.

Checkpoint: is “gene OFF” a sufficient state?

Suppose OFF pools two hidden configurations. One activates at $1\,\mathrm{s}^{-1}$ , the other at $10\,\mathrm{s}^{-1}$ , and they do not interconvert during an OFF episode. Equal initial weights give survival $S(t)=\tfrac12e^{-t}+\tfrac12e^{-10t}$ . Its instantaneous activation hazard is $[e^{-t}+10e^{-10t}]/[e^{-t}+e^{-10t}]$ . This falls from $5.5\,\mathrm{s}^{-1}$ toward $1\,\mathrm{s}^{-1}$ . The duration already spent OFF predicts the future. A single fixed-rate OFF state is therefore not exact.

Choose state variables by predictive sufficiency, and test a proposed lumped state for hidden memory.

2Derive the time and identity of the next event

A hazard specifies a short-time conditional probability. At a fixed state, let event channel $r$ have hazard $a_r\geq0$ . The probability that this event occurs in the next $\Delta t$ is $a_r\Delta t+o(\Delta t)$ . Hazards have inverse-time units. Multiple-jump probabilities are of smaller order than $\Delta t$ in this finite-rate model.

Survival converts that local statement into a waiting-time law. Let $T$ be the time to the next event, $S(t)=\Pr(T>t)$ , and $a_0=\sum_r a_r$ . Conditional on no event yet, the state has not changed. For constant hazards,

\begin{aligned} S(t+\Delta t)&=S(t)[1-a_0\Delta t]+o(\Delta t),\\ S'(t)&=-a_0S(t),\qquad S(0)=1,\\ S(t)&=e^{-a_0t},\qquad f_T(t)=a_0e^{-a_0t}. \end{aligned}

The mean waiting time is $\int_0^\infty S(t)\,dt=1/a_0$ . The ratio $S(t+s)/S(t)=e^{-a_0s}$ proves memorylessness: surviving for time $t$ does not change the residual waiting-time law. If $a_0=0$ , the state is absorbing and the next-event time is infinite.

Competing events give a joint density, not a second waiting-time assumption. Write $R$ for the random channel index in this calculation. The probability density that the first event is channel $r$ at time $t$ is the probability of surviving until $t$ , multiplied by that channel's hazard:

f_{T,R}(t,r)=a_re^{-a_0t} =\underbrace{a_0e^{-a_0t}}_{\text{event-time density}} \underbrace{\frac{a_r}{a_0}}_{\text{event probability}}.

This factorization gives Gillespie's two draws. Draw $T=-\ln U/a_0$ with $U$ uniform on $(0,1)$ . Choose event $r$ with probability $a_r/a_0$ . Update the state and recompute all hazards. The same law follows from the minimum of independent exponential clocks, since their survival probabilities multiply. The independent-clock construction applies at the current state, not unchanged throughout the future trajectory.¹

Scroll the figure horizontally to read labels at full size.

Figure 1. Two clocks produce one event-time law. Hazards are 2 and 1 per second. Left: individual survival probabilities and survival of their minimum. Right: the two joint event/time densities, with their sum dashed. Their integrals are 2/3 and 1/3. The combined mean wait is 1/3 s.

Checkpoint: what changes under a time-varying input?

If hazards depend on prescribed time, then $S(t)=\exp[-\int_0^t a_0(u)\,du]$ . The event identity conditional on an event at time $t$ has probability $a_r(t)/a_0(t)$ . The simple constant-rate exponential draw is no longer valid unless the changing input has been resolved appropriately.

Derive the next-event distribution and recognize when Gillespie's constant-hazard step applies.

3Assemble a probability-conserving generator

Each source state supplies one column of the generator. For distinct states $i,j$ , define $q_{ij}$ as the hazard from $j$ to $i$ . The second index is the source. Let $p_i(t)=\Pr[\xi(t)=i]$ and collect these probabilities in a column $p$ . The two-index quantity $q_{ij}$ is a transition hazard. A one-species subscript such as $q_E$ will denote a molecular total, as in earlier lectures.

The diagonal records all ways to leave the source. Set $q_{jj}=-\sum_{i\ne j}q_{ij}$ . Counting arrivals and departures over a short interval gives

\begin{aligned} p_i(t+\Delta t)&=p_i(t)+\Delta t\left[\sum_{j\ne i}q_{ij}p_j(t)+q_{ii}p_i(t)\right]+o(\Delta t),\\ \dot p_i&=\sum_jq_{ij}p_j,\qquad \dot p=Qp,\qquad \mathbf1^\top Q=0. \end{aligned}

Here $\mathbf1$ is the column of ones and $\top$ means transpose. The zero column sums imply $d(\mathbf1^\top p)/dt=0$ . Nonnegative off-diagonal entries prevent probability from flowing out of a state that currently has zero probability. These two checks are the first diagnostics for an assembled generator.

The same matrix acts backward on observables. Give state $j$ a numerical value $g_j$ . Its ensemble mean is $\langle g\rangle=g^\top p$ . Differentiating identifies the observable generator $\mathcal L$ :

\frac{d\langle g\rangle}{dt} =g^\top Qp=\langle\mathcal Lg\rangle,\qquad (\mathcal Lg)_j=(Q^\top g)_j=\sum_{i\ne j}q_{ij}(g_i-g_j).

The forward equation evolves probability columns. The backward operator acts on a function of the starting state. This distinction will determine the transpose in the first-passage equation.

Checkpoint: assemble a three-state chain.

Take transitions $0\to1$ , $1\to0$ , $1\to2$ , and $2\to1$ with hazards 2, 1, 4 and 7 per minute. These are state-transition arrows. The generator is

Q=\begin{pmatrix}-2&1&0\\2&-5&7\\0&4&-7\end{pmatrix}\mathrm{min}^{-1}.

Column 1, the middle column, contains two exits. Their sum is 5 per minute, so its diagonal is −5, not −1 or −4.

Build a column generator from a state graph and obtain the backward operator without changing conventions.

4Separate stationary occupancy from relaxation speed

A two-state promoter can be solved exactly. Label OFF as 0 and ON as 1. Let activation and deactivation hazards be $a>0$ and $b>0$ . Its generator and ON probability equation are

Q=\begin{pmatrix}-a&b\\a&-b\end{pmatrix},\qquad \dot p_{\mathrm{on}}=a(1-p_{\mathrm{on}})-bp_{\mathrm{on}}.

Subtract the stationary value to isolate the decaying part. Setting the derivative to zero gives $p_{\mathrm{on}}^*=a/(a+b)$ . For $e=p_{\mathrm{on}}-p_{\mathrm{on}}^*$ , the equation becomes $\dot e=-(a+b)e$ . Therefore

p_{\mathrm{on}}(t)=\frac{a}{a+b} +\left[p_{\mathrm{on}}(0)-\frac{a}{a+b}\right]e^{-(a+b)t}, \qquad t_{\mathrm{relax}}=\frac1{a+b}.

The ratio fixes occupancy while the sum fixes speed. Multiplying both hazards by ten leaves the stationary fraction unchanged and divides the relaxation time by ten. If both fraction and relaxation time are measured, the two hazards become identifiable within this model: $a=p_{\mathrm{on}}^*/t_{\mathrm{relax}}$ and $b=(1-p_{\mathrm{on}}^*)/t_{\mathrm{relax}}$ .

Scroll horizontally for the full-size figure.

Figure 2. The same stationary fraction can be reached on different clocks. All curves start OFF. Hazard pairs are (0.2, 0.1), (2, 1) and (20, 10) per second. Every curve tends to 2/3. Their relaxation times differ by factors of ten. Time is plotted logarithmically to show all three transients.

Checkpoint: distinguish three times.

For $a=2,b=1\,\mathrm{s}^{-1}$ , the mean OFF holding time is $1/a=0.5\,\mathrm{s}$ . The mean ON holding time is $1/b=1\,\mathrm{s}$ . Ensemble relaxation takes $1/(a+b)=1/3\,\mathrm{s}$ . A residence time and a relaxation time answer different questions.

Infer what occupancy data determine, what they leave undetermined, and which transient measurement supplies the missing clock.

5Recover the count equation and its Poisson law

A chemical master equation is a generator equation on molecule-count states. For species $X_1,\ldots,X_n$ , let $N_X(t)$ be the random count vector, with entries $N_{X_i}(t)$ . Reaction channel $r$ changes the vector by $\gamma_r$ and has propensity $a_r(n)$ in count state $n$ . Insert its hazard into $Q_{n+\gamma_r,n}$ . Sum contributions when distinct channels have the same endpoints. The forward equation becomes

\dot p(n,t)=\sum_r\left[a_r(n-\gamma_r)p(n-\gamma_r,t)-a_r(n)p(n,t)\right].

The first term counts arrivals from a predecessor count state. An impossible predecessor contributes zero. The second term counts departures. The graph's nodes are entire molecular populations, not individual chemical species. This is the connection between a reaction list, the count-state graph and the master equation.

Constant addition and first-order removal give a birth-death chain. The composite chemical accounting is $\varnothing\rightsquigarrow X$ , $X\rightsquigarrow\varnothing$ . Declare the separate propensities $a_+(n)=b_0$ and $a_-(n)=dn$ , with $b_0>0,d>0$ . An event adds or removes one molecule. The probability equation is

\dot p_n=b_0p_{n-1}+d(n+1)p_{n+1}-(b_0+dn)p_n,\qquad p_{-1}=0.

The boundary at zero forces zero stationary edge current. Define $J_n=b_0p_n^*-d(n+1)p_{n+1}^*$ . Stationarity gives $J_n=J_{n-1}$ . At state zero there is no incoming current from a negative count, so $J_0=0$ . Successive ratios and normalization then give

\begin{aligned} p_{n+1}^*&=\frac{\rho}{n+1}p_n^*,\qquad \rho=\frac{b_0}{d},\\ p_n^*&=p_0^*\frac{\rho^n}{n!},\qquad 1=p_0^*\sum_{n=0}^{\infty}\frac{\rho^n}{n!}=p_0^*e^\rho,\\ p_n^*&=e^{-\rho}\frac{\rho^n}{n!}. \end{aligned}

This normalizable stationary law is Poisson. Directly shifting the sum gives $\langle N_X\rangle=\rho$ and $\langle N_X(N_X-1)\rangle=\rho^2$ , hence $\operatorname{Var}N_X=\rho$ . Positive removal is essential. With $d=0,b_0>0$ , counts drift upward and no stationary law exists. A normalizable count law does not establish physical equilibrium of the unresolved synthesis and removal mechanisms.

Scroll horizontally for the full-size figure.

Figure 3. A stationary distribution is more than its mean. Poisson probabilities for mean counts 1, 4 and 10. The coefficient of variation is

1/\sqrt{\rho}

. Increasing both addition and removal speeds leaves each histogram unchanged.

The mean equation follows from the backward operator without closing a hierarchy. For $g(n)=n$ , one addition changes $g$ by +1 and one removal changes it by −1. Thus $\mathcal Lg=b_0-dn$ and

\frac{d\langle N_X\rangle}{dt}=b_0-d\langle N_X\rangle,\qquad \frac{d\langle X\rangle}{dt}=\frac{b_0}{\Omega}-d\langle X\rangle,\qquad X=\frac{N_X}{\Omega},\quad\Omega=N_AV.

Here $V$ is the fixed volume in liters when $X$ is molar, and $N_A$ is Avogadro's constant. The count propensity $b_0$ is an event rate. Because each event adds one molecule, its numerical value also gives molecules added per time. Dividing by $\Omega$ produces a concentration flux.

Keep the deterministic rate constant

For an elementary channel with reactant coefficients $\alpha_{ir}$ , deterministic event flux $v_r(x)=k_r\prod_i x_i^{\alpha_{ir}}$ , and molecularity $m_r=\sum_i\alpha_{ir}$ , the matching count propensity is

a_r(n)=k_r\Omega^{1-m_r}\prod_i(n_i)_{\alpha_{ir}},\qquad (n)_a=n(n-1)\cdots(n-a+1).

This retains the concentration convention for $k_r$ . For elementary $2A\xrightarrow{k}C$ , the event flux is $kA^2$ and $\dot A=-2kA^2$ . The propensity is $kn_A(n_A-1)/\Omega$ . There are $\binom{n_A}{2}$ unordered pairs, so the corresponding per-pair coefficient is $2k/\Omega$ . The factorial belongs in the coefficient if unordered combinations are used. It must not be inserted a second time.

For a general network, the exact mean relation is $d\langle X\rangle/dt=\Gamma\langle a(N_X)\rangle/\Omega$ . Replacing the expected propensity by its value at the mean is an additional approximation for nonlinear reactions. The linear addition-removal example needs no such approximation.

Move consistently between a reaction list, a count-state generator, a stationary law and a concentration-rate equation.

6Calculate a first-arrival time with the backward equation

A first-arrival question needs a target and a starting state. Let $A$ be a target set. Let $\tau_j$ be the expected first time to reach $A$ starting from state $j$ . Set $\tau_j=0$ on the target. Assume the target is reached with finite mean time. A finite irreducible chain and a nonempty target satisfy this condition.

Condition on the first holding interval to derive the equation. The departure hazard from $j$ is $a_j=-q_{jj}$ . Its mean holding time is $1/a_j$ . The next state is $i$ with probability $q_{ij}/a_j$ . Consequently

\begin{aligned} \tau_j&=\frac1{a_j}+\sum_{i\ne j}\frac{q_{ij}}{a_j}\tau_i,\\ -1&=\sum_iq_{ij}\tau_i=(Q^\top\tau)_j,\qquad j\notin A. \end{aligned}

The reduced matrix retains its original diagonal. Let $B$ be the states outside the target. Remove target rows and columns from $Q$ , then solve $Q_{BB}^{\top}\tau_B=-\mathbf1$ . Do not recompute the diagonal after deleting target transitions. Those exits are exactly how probability is absorbed. The transpose appears because arrival time is a function of the starting state.

A closed/open/ON model gives an explicit biological example. Treat these labels as inaccessible, accessible and transcriptionally active configurations. They are a modelling choice, not a universal three-state description of mammalian promoters. Let closed-to-open, open-to-closed and open-to-ON hazards be $a,b,c$ . The equations are

\begin{aligned} \tau_{\mathrm{closed}}&=\frac1a+\tau_{\mathrm{open}},\\ \tau_{\mathrm{open}}&=\frac1{b+c}+\frac{b}{b+c}\tau_{\mathrm{closed}},\\ \tau_{\mathrm{open}}&=\frac{1+b/a}{c},\qquad \tau_{\mathrm{closed}}=\frac1a+\frac1c+\frac{b}{ac}. \end{aligned}

Substitution of the first equation into the second produces the last line. The extra term $b/(ac)$ is the delay caused by returns to the closed state. With $a=2,b=1,c=4\,\mathrm{min}^{-1}$ , the answers are 0.375 min from open and 0.875 min from closed. Any ON-to-open rate is irrelevant to this first arrival, since the clock stops at ON.

Scroll horizontally for the full-size figure.

Figure 4. A return path increases the first-arrival delay. The arrows are transitions between configurations. Above: the target is ON. Below: mean arrival times as the reclosing hazard varies, with opening 2/min and activation 4/min. Both delays increase linearly with reclosing. Their difference remains the 0.5 min needed for the first opening.

Checkpoint: can removal of the target be treated as a smaller closed chain?

No. With one starting state and one target connected by hazard $k$ , the retained matrix is $[-k]$ . It gives $\tau=1/k$ . Replacing its diagonal by zero would remove the absorption and make the equation unsolvable.

Calculate a first-passage time to one state or a target set using the correct backward generator.

Part II

Test whether stationarity has equilibrium structure

Separate node balance, resolved edge balance and physical reservoir compatibility before using an energy function.

7Solve a stationary ring and its current

Stationarity requires zero net accumulation at each node. A stationary probability column $p^*$ satisfies $Qp^*=0$ and $\mathbf1^\top p^*=1$ . A finite irreducible chain has a unique positive stationary distribution. Numerically, replace one redundant row of $Qp^*=0$ with normalization. Check the residual, nonnegativity and normalization after solving.

Three states reveal what node balance leaves free. Orient the ring $1\to2\to3\to1$ . Call its forward hazards $a,b,c$ and reverse hazards $r,s,u$ , respectively. All are positive and have inverse-time units. Define the one-way flux $\phi_{i\to j}=q_{ji}p_i$ and net current $J_{i\to j}=\phi_{i\to j}-\phi_{j\to i}$ . Node balance gives

ap_1^*-rp_2^*=bp_2^*-sp_3^*=cp_3^*-up_1^*=J_{\mathrm{cycle}}.

Eliminate the probabilities to obtain the current. The first two equalities imply $p_2^*=(ap_1^*-J_{\mathrm{cycle}})/r$ and $p_3^*=[abp_1^*-(b+r)J_{\mathrm{cycle}}]/(rs)$ . Insert the latter into the third equality. This gives $(abc-rsu)p_1^*=(bc+rc+rs)J_{\mathrm{cycle}}$ . Repeating the cyclic relabeling and imposing normalization yields

\begin{aligned} w_1&=rc+bc+rs,&w_2&=as+us+ac,&w_3&=ub+ab+ur,\\ D&=w_1+w_2+w_3,\qquad p_i^*=\frac{w_i}{D},\qquad J_{\mathrm{cycle}}=\frac{abc-rsu}{D}. \end{aligned}

The weights in this kinetic solve have units of inverse time squared. They are not yet Boltzmann factors. Each term is the product of the hazards in a directed spanning tree whose arrows lead to its root. The general matrix-tree formula extends this construction to larger finite graphs.²

Equal occupancies do not imply equal currents. If all six hazards are 1/s, then $p_i^*=1/3$ and the current is zero. If every forward hazard is 2/s and every reverse hazard is 1/s, the same probabilities coexist with $J_{\mathrm{cycle}}=1/3\,\mathrm{s}^{-1}$ . Forward traffic is 2/3 per second on each edge and reverse traffic is 1/3 per second. The probabilities stay constant because the circulating inflow equals circulating outflow.

Scroll horizontally for the full-size figure.

Figure 5. Measure traffic as well as occupancy. A symmetric equilibrium ring and a symmetric driven ring have identical state probabilities. Direction-resolved event frequencies distinguish them. The common edge current counts net cycle completions per time. Adding the three edge currents would count three transitions per completed cycle.

Checkpoint: solve a nonsymmetric ring.

With forward hazards (2, 3, 4)/s and reverse hazards (1, 6, 1)/s, the three weights are (22, 26, 10)/s². Thus $p^*=(11,13,5)^\top/29$ and $J_{\mathrm{cycle}}=9/29\,\mathrm{s}^{-1}$ . Substituting into all three edge currents gives the same answer. This checks both the orientation and the normalization.

Solve a stationary ring and distinguish constant probability from zero microscopic circulation.

8Convert an enzyme mechanism into hazards, flux and work

The ring can represent the states of one supplied enzyme. Resolve the elementary pairs

\begin{aligned} E+S&\xrightleftharpoons[k_{-S}]{k_{+S}}C_{ES},\\ C_{ES}&\xrightleftharpoons[k_{\mathrm{cat}}^-]{k_{\mathrm{cat}}^+}C_{EP},\\ E+P&\xrightleftharpoons[k_{-P}]{k_{+P}}C_{EP}. \end{aligned}

Take free substrate $S$ and product $P$ as maintained reservoirs. The enzyme states are $E,C_{ES},C_{EP}$ . In that order, the three clockwise hazards are $k_{+S}S,k_{\mathrm{cat}}^+,k_{-P}$ . Their reverse hazards are $k_{-S},k_{\mathrm{cat}}^-,k_{+P}P$ . Each association constant has concentration-inverse/time units. Multiplication by its free ligand converts it into an inverse-time hazard.

Scroll horizontally for the full-size figure.

Figure 6. Reservoir conversion drives the enzyme-state cycle. The molecular mechanism has straight arrows labeled by elementary constants. The enzyme configuration graph has state-transition hazards. Supply and removal are composite boundary processes. Sustaining different substrate and product chemical potentials can maintain circulation even though enzyme occupancies are stationary.

One cycle current gives three distinct dimensional outputs. For one enzyme, $J_{\mathrm{cycle}}$ counts net conversions per time. For enzyme concentration $q_E$ , the concentration throughput is $v=q_EJ_{\mathrm{cycle}}$ . For an actual enzyme count $N_E$ , the expected total event current is $N_EJ_{\mathrm{cycle}}$ . None of these quantities is a microscopic association constant.

The rate-product ratio identifies the cycle's dimensionless driving. Define the affinity

\mathcal A=\ln\frac{k_{+S}S\,k_{\mathrm{cat}}^+k_{-P}} {k_{-S}k_{\mathrm{cat}}^-k_{+P}P}.

For a thermodynamically consistent enzyme tightly coupled to one substrate-to-product conversion, with no other work exchange, local detailed balance gives $\mathcal A=(\mu_S-\mu_P)/(k_BT)$ . Here chemical potentials are per molecule. Molar chemical potentials instead require $RT$ in the denominator. Keeping one convention prevents an Avogadro-factor error.³

Current times affinity is a dissipation rate only with those physical assumptions. At steady state, the entropy-production rate per enzyme is $\dot S_{\mathrm{prod}}=k_BJ_{\mathrm{cycle}}\mathcal A$ . The dissipated power is $T\dot S_{\mathrm{prod}}$ . Multiplying all hazards by the same factor leaves affinity and occupancies unchanged, but multiplies current and power. A large driving force does not guarantee a large throughput.

Checkpoint: what maintains the reservoirs?

Supplying $S$ and removing $P$ maintains the chemical free-energy difference. Each boundary event need not hydrolyze ATP directly. The larger metabolic network may supply the free energy. Without maintenance, a closed compatible system changes its substrate/product ratio and eventually removes the driving. Dissipation is not useful work output unless a coupled load has also been specified.

A glycolytic conversion with an explicit energy bill

Triosephosphate isomerase supplies a concrete substrate/product pair. In glycolysis, dihydroxyacetone phosphate (DHAP) interconverts with glyceraldehyde 3-phosphate (GAP). Downstream reactions consume GAP. The overall relation $\mathrm{DHAP}\leftrightsquigarrow\mathrm{GAP}$ does not itself consume ATP. Use the three-state enzyme as a coarse teaching mechanism for this conversion, not as a complete experimentally resolved isomerase mechanism.

A published pathway calculation provides an energy scale. Flamholz and colleagues optimized glycolytic metabolite concentrations between $1\,\mu\mathrm M$ and $10\,\mathrm{mM}$ , with transformed energies at pH 7.5 and ionic strength 0.2 M. They found a bottleneck driving force near $2.9\,\mathrm{kJ/mol}$ for the classical pathway, including the isomerase step. The number is a constrained-model result, not a universal measurement. Our calculation uses that energy scale while choosing the kinetics independently and explicitly.⁹

Convert the energy to a dimensionless affinity before selecting hazards. At the chosen teaching temperature $T=298.15\,\mathrm K$ , the gas constant gives $RT=2.478957\,\mathrm{kJ/mol}$ . With molar reaction free-energy change $\Delta_rG=-2.9\,\mathrm{kJ/mol}$ ,

\mathcal A=-\frac{\Delta_rG}{RT}=1.169847,\qquad \frac{\bar P}{\bar S}=e^{-\mathcal A}=0.310414, \qquad \bar S=\frac{S}{K_{d,S}},\quad\bar P=\frac{P}{K_{d,P}}.

Choose a transparent kinetic family with this affinity. Give both binding pairs the same dissociation frequency $k_{\mathrm{off}}$ and both catalytic directions the same frequency $k_{\mathrm{cat}}$ . The clockwise hazards are $(k_{\mathrm{off}}\bar S,k_{\mathrm{cat}},k_{\mathrm{off}})$ . The reverse hazards are $(k_{\mathrm{off}},k_{\mathrm{cat}},k_{\mathrm{off}}\bar P)$ . Their product ratio is exactly $\bar S/\bar P=e^{\mathcal A}$ . For this chosen family, the concentration equilibrium constant is $K_{\mathrm{eq}}=K_{d,P}/K_{d,S}$ . No unverified absolute affinity or metabolite concentration is required.

The free-enzyme balance first determines its exact probability. Let $p_E,p_{ES},p_{EP}$ denote probabilities, not concentrations. Binding leaves free enzyme at frequency $k_{\mathrm{off}}(\bar S+\bar P)$ and returns it at frequency $k_{\mathrm{off}}(p_{ES}+p_{EP})$ per enzyme. Combining its stationary balance with normalization gives

(\bar S+\bar P)p_E=p_{ES}+p_{EP}=1-p_E,\qquad p_E=\frac1{1+\bar S+\bar P}.

Subtract the two bound-state balances to find the conversion current. Their difference eliminates the common occupancy contribution. The result and its current are

\begin{aligned} 0&=k_{\mathrm{off}}(\bar S-\bar P)p_E -(k_{\mathrm{off}}+2k_{\mathrm{cat}})(p_{ES}-p_{EP}),\\ J&=k_{\mathrm{cat}}(p_{ES}-p_{EP})\\ &=\frac{k_{\mathrm{cat}}k_{\mathrm{off}}}{k_{\mathrm{off}}+2k_{\mathrm{cat}}} \frac{\bar S-\bar P}{1+\bar S+\bar P}. \end{aligned}

This is exact for the declared three-state family. The sum $p_{ES}+p_{EP}=1-p_E$ and the difference above recover both bound probabilities. In the limit $k_{\mathrm{off}}\gg k_{\mathrm{cat}}$ , the probabilities approach weights proportional to $(1,\bar S,\bar P)$ . The current approaches a finite value. Exact equilibrium of the binding-only problem therefore does not assert zero current in the full catalytic problem.

A reproducible energy and current calculation

Choose kinetics, then calculate the output. Set $\bar S=1$ , $k_{\mathrm{off}}=100\,\mathrm{s}^{-1}$ , $k_{\mathrm{cat}}=10\,\mathrm{s}^{-1}$ and $q_E=0.1\,\mu\mathrm M$ . These are illustrative, not fitted isomerase parameters. The exact solution gives

\begin{aligned} p_E&=0.432823,\\ p_{ES}&=0.407950,\\ p_{EP}&=0.159227,\\ J&=2.487236\,\mathrm{s}^{-1}. \end{aligned}

Multiply event current by the cost per event. The entropy-production rate is $k_BJ\mathcal A=2.909685\,k_B/\mathrm s$ per enzyme. Since $k_BT=4.116405\times10^{-21}\,\mathrm J$ , the dissipated power is about $12.0\times10^{-21}\,\mathrm W$ per enzyme. This energy is released by the conversion's maintained chemical-potential difference, not by assigning one ATP hydrolysis to each turn.

Concentration flux gives power per volume. The flux is $q_EJ=0.248724\,\mu\mathrm M/\mathrm s$ . Converting micromoles to moles and kilojoules to joules gives

\mathcal P_V= (0.248724\times10^{-6}\,\mathrm{mol\,L^{-1}s^{-1}}) \times(2.9\times10^3\,\mathrm{J\,mol^{-1}}) \approx721.3\,\mu\mathrm{W/L}.

This power is for the modelled conversion at the chosen enzyme concentration, not for a whole cell or the entire glycolytic pathway.

The edgewise entropy calculation checks the cycle budget. Write directional stationary traffic as $F_i=f_ip_i$ and $B_i=r_ip_{i+1}$ , with cyclic indices. Local detailed balance assigns each resolved edge the contribution $k_B(F_i-B_i)\ln(F_i/B_i)$ . At stationarity every difference equals $J$ . Summing and cancelling the probabilities gives

\frac{\dot{\mathcal S}_{\mathrm{prod}}}{k_B} =\sum_i(F_i-B_i)\ln\frac{F_i}{B_i} =J\ln\frac{\prod_i f_i}{\prod_i r_i} =J\mathcal A.

Do not replace each edge's traffic ratio by the complete cycle's affinity. Some forward binding events simply unbind without completing a conversion. The sum over the physically resolved cycle is what establishes the budget. At equilibrium opposing traffic survives but each difference vanishes. Under reverse driving both $J$ and $\mathcal A$ are negative, leaving nonnegative dissipation.

The larger pathway raises a different accounting question. Its dissipation is the sum $\sum_r v_r(-\Delta_rG_r)$ , with stoichiometric fluxes and consistent molar units. ATP formation must remain explicitly in those reaction energies. It cannot be counted once as conserved chemical output and again as dissipated heat. The core playground keeps the single conversion visible so these distinctions can be explored without claiming a whole-pathway kinetic reconstruction.

Translate a biochemical mechanism into a state generator, solve its current, and turn a stated conversion energy into entropy production and power with consistent units.

9Construct a potential only when all routes agree

Detailed balance cancels flux pair by pair. A positive normalized distribution $p^{\mathrm{eq}}$ satisfies detailed balance if $q_{ji}p_i^{\mathrm{eq}}=q_{ij}p_j^{\mathrm{eq}}$ for every resolved reverse pair. Summing these equalities at each node proves stationarity. The driven ring shows why stationarity does not imply them.

Finite-chain equilibrium criterion

For a finite, irreducible, time-homogeneous Markov chain with a positive reverse hazard for every present edge, the following statements are equivalent: a detailed-balanced probability distribution exists, the forward and reverse hazard products agree on every cycle, and all edge ratios admit a common state potential. Section 10 adds the equivalent stationary time-reversal statement. Physical equilibrium additionally requires that the graph resolves the relevant physical channels.²

Detailed balance implies the cycle condition by cancellation. On an oriented edge, $q_{ji}/q_{ij}=p_j^{\mathrm{eq}}/p_i^{\mathrm{eq}}$ . Multiply around a cycle. Every probability appears once in a numerator and once in a denominator, leaving

\prod_{(i\to j)\in C}\frac{q_{ji}}{q_{ij}}=1,\qquad \mathcal A_C=\sum_{(i\to j)\in C}\ln\frac{q_{ji}}{q_{ij}}=0.

The converse constructs the weights rather than guessing them. Choose a root state and assign it weight 1. For another state, choose any path from the root and multiply the forward/reverse hazard ratios along that path. Two different paths form a closed walk when one is reversed. The cycle condition makes their ratio one, so they assign the same weight. Finite connectivity makes every weight finite and positive. Normalize their sum to obtain $p_i^{\mathrm{eq}}=w_i/\sum_jw_j$ .

A state potential is the logarithm of these compatible weights. Define $G_i=-k_BT\ln w_i$ , up to an additive constant. The construction gives both a probability formula and an edge constraint:

p_i^{\mathrm{eq}}=\frac{e^{-G_i/(k_BT)}}{Z},\qquad Z=\sum_je^{-G_j/(k_BT)},\qquad \frac{q_{ji}}{q_{ij}}=e^{-(G_j-G_i)/(k_BT)}.

For the abstract graph, the temperature factor simply gives the potential energy units. Its identification with physical free energy requires the reservoir assumptions in the next sections. Writing $-k_BT\ln p_i^*$ for an arbitrary stationary distribution is not sufficient. The driven ring also permits that expression, but its edge ratios violate the last equality.

Scroll horizontally for the full-size figure.

Figure 7. Path independence is the compatibility test. The indirect route assigns weight ratio 4. A direct ratio of 4 agrees and gives weights (1, 2, 4). A direct ratio of 1/2 conflicts. Labels are ratios of state-transition hazards, not chemical rate constants. The lower panel states the corresponding stationary time-reversal identity.

A cycle basis makes the test practical. In a connected graph with $V$ vertices and $E$ undirected edge pairs, select a spanning tree. Each of its $E-V+1$ omitted edges closes a fundamental cycle. Checking those cycles suffices because logarithmic affinities add along closed walks. A tree has no such tests and always admits compatible weights for positive reverse pairs.

Checkpoint: what does a tree prove?

At a leaf, stationarity forces its only edge current to vanish. Remove that leaf and repeat. Thus a finite closed tree cannot support stationary circulation in the states and channels retained. This does not prove that unresolved molecular mechanisms consume no fuel. Nor does choosing a spanning tree make all its edges current-free once the omitted edges are restored.

Test independent cycle products and construct equilibrium weights without solving the stationary linear system.

10Reverse a stationary trajectory law

Time reversal compares probabilities of histories, not individual movies. A random equilibrium trajectory usually looks different from its own reversal. Reversibility means that each allowed trajectory and its reversed trajectory have equal probability density when both begin in the stationary ensemble.

Bayes' rule gives the reversed hazard. For a stationary process with probabilities $p^*>0$ , a reversed jump from $j$ to $i$ corresponds to an original jump from $i$ to $j$ . Dividing that stationary joint probability by the probability of the reversed starting state yields

q^{\mathrm{rev}}_{ij} =\lim_{\Delta t\downarrow0} \frac{\Pr[\xi(t-\Delta t)=i\mid\xi(t)=j]}{\Delta t} =q_{ji}\frac{p_i^*}{p_j^*},\qquad i\ne j.

Stationarity ensures that the reversed departure sum equals the original departure sum. Therefore their diagonals also agree. The identity $Q^{\mathrm{rev}}=Q$ holds exactly when detailed balance holds.

The same result can be checked for an entire jump path. Consider states $i_0,\ldots,i_m$ , with a specified duration spent in each visited state during a fixed observation window. The forward path density contains its initial probability, the product of its jump hazards and the exponential survival factors. The reversed path contains the same survival factors in reverse order. Their ratio is

\frac{\Pr[\text{forward path}]}{\Pr[\text{reversed path}]} =\frac{p_{i_0}^*}{p_{i_m}^*} \prod_{\ell=0}^{m-1}\frac{q_{i_{\ell+1},i_\ell}}{q_{i_\ell,i_{\ell+1}}}.

At detailed balance, each hazard ratio is the ratio of the destination and source equilibrium probabilities. The product telescopes and cancels the initial-probability ratio. Conversely, comparing one-jump stationary histories recovers pairwise detailed balance.

Equilibrium still contains molecular traffic. Equal forward and reverse frequencies need not be small. Increasing both hazards can produce many more jumps without changing weights or net current. The theorem here concerns discrete molecular configurations and counts, which are unchanged by physical reversal of time. Mechanical states with velocities require reversing those velocities as well.

Checkpoint: reverse the driven ring.

For uniform stationary probabilities, reversing the (2/s forward, 1/s reverse) ring exchanges its forward and reverse hazards. A clockwise three-jump cycle has hazard-product ratio 8 relative to its reversed cycle. The occupancy histogram cannot reveal this asymmetry, but a direction-resolved trajectory can.

Derive the reversed generator and identify which trajectory measurement distinguishes equilibrium from a circulating steady state.

11Choose a physical free energy and explain relaxation

Equilibrium is constrained by what the environment holds fixed. A system does not generally minimize internal energy and maximize its own entropy simultaneously. The appropriate extremum follows from the allowed exchange. Here $U$ is internal energy, $\mathcal S$ entropy, $T$ temperature and $P_{\mathrm{ext}}$ external pressure.

Held fixed	Equilibrium criterion	Ensemble
Energy, volume, constituent amounts	Maximize total entropy	Microcanonical
Temperature, volume, constituent amounts	Minimize $F=U-T\mathcal S$	Canonical
Temperature, pressure, constituent amounts	Minimize $G=U+P_{\mathrm{ext}}V-T\mathcal S$	Isothermal-isobaric
Temperature, volume, compatible reservoir chemical potentials	Minimize $\Phi=F-\sum_\ell\mu_\ell N_\ell$	Grand canonical

Constituent conservation constrains reacting mixtures. It does not freeze every reacting species count. Reservoirs are compatible when they do not impose a persistent thermal or chemical driving around a cycle. Contact with several maintained reservoirs therefore does not by itself imply equilibrium.

State free energies already contain unresolved internal entropy. For microstates with energies $E_i$ and degeneracies $g_i$ , a canonical weight is $g_i e^{-E_i/(k_BT)}$ . One may instead define a coarse free energy $G_i=E_i-k_BT\ln g_i$ . Use either explicit degeneracy or that entropy contribution, not both. When ligands exchange with a reservoir, the weight additionally includes $e^{\mu N/(k_BT)}$ for the molecules transferred into the subsystem.

The normalized weights minimize an ensemble free-energy functional. For given effective state free energies $G_i$ , define $\mathcal F[p]=\sum_i p_iG_i+k_BT\sum_i p_i\ln p_i$ . Substituting $G_i=-k_BT\ln p_i^{\mathrm{eq}}-k_BT\ln Z$ gives

\mathcal F[p]-\mathcal F[p^{\mathrm{eq}}] =k_BT\sum_i p_i\ln\frac{p_i}{p_i^{\mathrm{eq}}}\geq0.

The final inequality is nonnegativity of relative entropy. It follows, for positive probabilities, from $-\ln z\geq1-z$ , applied with $z=p_i^{\mathrm{eq}}/p_i$ and summed with weights $p_i$ . The limiting convention $0\ln0=0$ covers zero probabilities. This establishes the equilibrium distribution from the specified free energies. It does not turn an arbitrary fitted stationary distribution into a physical equilibrium model.

Use the excess free energy as a Lyapunov function

A minimum alone does not explain why the dynamics approaches it. To establish that link, hold the generator, temperature and compatible reservoirs fixed. Define the dimensionless excess free energy $D(p\Vert p^{\mathrm{eq}})=\sum_i p_i\ln(p_i/p_i^{\mathrm{eq}})$ . We will show that its time derivative is negative away from equilibrium. The probability simplex, not the original molecular configuration space, is the state space for this deterministic relaxation argument.

Detailed balance makes every edge act on an occupancy ratio. Write $u_i=p_i/p_i^{\mathrm{eq}}$ and $c_{ij}=q_{ji}p_i^{\mathrm{eq}}=q_{ij}p_j^{\mathrm{eq}}$ . The symmetric quantity $c_{ij}$ is equilibrium traffic, in inverse-time units. The net current from i to j is $c_{ij}(u_i-u_j)$ . Thus probability flows from a state overrepresented relative to its equilibrium weight toward one less overrepresented.

Differentiate and pair the two ends of each edge. Probability conservation removes the derivative of the normalization term. The current removes probability from i and adds it to j, giving

\begin{aligned} \frac{dD}{dt}&=\sum_i\dot p_i\ln u_i\\ &=-\sum_{i<j}c_{ij}(u_i-u_j)(\ln u_i-\ln u_j)\leq0. \end{aligned}

Each product is nonnegative because logarithm is increasing. For a connected graph with positive reverse hazards, equality requires the ratios to agree across every edge. They must then agree everywhere. Normalization fixes their common value to one. The unique zero of the derivative is the equilibrium distribution.

This is an energy argument for global relaxation. Relative entropy is nonnegative and vanishes only at equilibrium. Its derivative is negative elsewhere in the interior. Finite irreducible continuous-time chains enter that interior for every positive time, even from a boundary initial distribution. This gives the Lyapunov interpretation from Lecture 4 without computing every eigenvalue. It does not determine the speed of relaxation, which still depends on the edge traffic.

Scroll horizontally for the full-size figure.

Figure 7A. Relaxation is not the same as absence of dissipation. The two rings start with probabilities (0.8, 0.1, 0.1). Their stationary probabilities are uniform. Left: relative entropy to that distribution. Right: entropy production from resolved directional traffic. The driven ring continues producing entropy at

k_B\ln2

per second after its distribution has settled. The hazards and their different clocks are stated in the figure.

Why decreasing relative entropy does not prove detailed balance

A general stationary distribution also contracts relative entropy. Let $\pi>0$ satisfy $Q\pi=0$ , without detailed balance. Set $u_i=p_i/\pi_i$ . Using the generator and $\ln z\leq z-1$ gives

\begin{aligned} \frac{d}{dt}D(p\Vert\pi) &=\sum_{i\ne j}q_{ij}\pi_j u_j\ln\frac{u_i}{u_j}\\ &\leq\sum_{i\ne j}q_{ij}\pi_j(u_i-u_j)=0. \end{aligned}

The last cancellation uses stationarity: $\sum_{j\ne i}q_{ij}\pi_j=-q_{ii}\pi_i$ . The first sum in the final expression therefore equals the second. This is an attraction statement, not a claim of physical equilibrium.

Entropy production instead measures irreversible traffic. For physically resolved reverse pairs, define $F_{ij}=q_{ji}p_i$ and $B_{ij}=q_{ij}p_j$ . Under local detailed balance,

\frac{\dot{\mathcal S}_{\mathrm{prod}}}{k_B} =\sum_{i<j}(F_{ij}-B_{ij})\ln\frac{F_{ij}}{B_{ij}}.

At detailed balance this equals $-dD/dt$ . In a driven stationary ring, $D=0$ but the traffic expression stays positive. The stationary histogram cannot supply the missing cycle information. Seifert develops this distinction between stationary probabilities and irreversible currents.¹⁰

Put ligand exchange into the same energy convention

Chemical potentials require dimensionless concentration ratios. For an ideal dilute solute of free concentration $R$ , write $\mu_R=\mu_R^\circ+k_BT\ln(R/c^\circ)$ , where $c^\circ$ is a standard concentration. For binding one repressor to a gene, define the standard binding free-energy change $\Delta G_{\mathrm{bind}}^\circ$ . The relative bound-state weight is

\frac{w_{C_{GR}}}{w_G} =\frac{R}{c^\circ}e^{-\Delta G_{\mathrm{bind}}^\circ/(k_BT)} =\frac{R}{K_d},\qquad K_d=c^\circ e^{\Delta G_{\mathrm{bind}}^\circ/(k_BT)}.

The binding dissociation constant $K_d$ has concentration units. The weight and every logarithm's argument are dimensionless. The binding free energy includes solvent, translational, rotational and conformational contributions at the stated standard state. It is not the energy of one chemical contact.

Checkpoint: change the standard concentration.

Set $c^{\circ\prime}=\eta c^\circ$ . The same physical affinity requires $\Delta G_{\mathrm{bind}}^{\circ\prime}=\Delta G_{\mathrm{bind}}^\circ-k_BT\ln\eta$ . Substitution leaves $K_d$ and $R/K_d$ unchanged. An arbitrary standard-state convention cannot alter a measured binding curve.

Specify the physical free energy, use it to explain passive relaxation, and keep relaxation distinct from the absence of sustained dissipation.

12Check compatibility in a nonlinear reaction network

For nonlinear chemistry, stoichiometric dependencies replace the simple graph-cycle test. Choose one orientation for each positive elementary reversible pair $\sum_i\alpha_{ir}X_i\xrightleftharpoons[k_r^-]{k_r^+}\sum_i\beta_{ir}X_i$ . Its change vector is $\gamma_r=\beta_r-\alpha_r$ . Collect the change vectors in the species-by-reaction matrix $\Gamma$ . The elementary forward and reverse concentration fluxes are $k_r^+\prod_i x_i^{\alpha_{ir}}$ and $k_r^-\prod_i x_i^{\beta_{ir}}$ .

Taking logarithms makes detailed balance a linear consistency problem. Define dimensionless activities $a_i=x_i/c^\circ$ and dimensionless equilibrium constants

K_r=\frac{k_r^+}{k_r^-}(c^\circ)^{|\alpha_r|-|\beta_r|}, \qquad |\alpha_r|=\sum_i\alpha_{ir},\qquad \Gamma^\top\ln a^{\mathrm{eq}}=\ln K.

To obtain the last equation, equate the two fluxes, divide their concentration monomials, and take logarithms. The vector logarithm acts component by component. A positive equilibrium composition corresponds to a real solution for $\ln a^{\mathrm{eq}}$ .

Wegscheider compatibility conditions

A positive detailed-balanced composition exists if and only if

z^\top\ln K=0\quad\text{for every }z\in\ker\Gamma, \qquad\text{equivalently}\qquad\prod_rK_r^{z_r}=1.

This is a statement about positive elementary reversible mass-action pairs. It does not, by itself, prove attraction from every initial condition or select a particular conserved class.⁴

Necessity and sufficiency are ordinary linear algebra. If $\Gamma z=0$ , multiplying the logarithmic balance equation by $z^\top$ gives zero on its left. Conversely, the image of $\Gamma^\top$ is the orthogonal complement of $\ker\Gamma$ . Orthogonality of $\ln K$ to every such $z$ therefore supplies a solution. Exponentiating its components produces a positive composition.

A monomolecular cycle recovers the Markov result. For $A\xrightleftharpoons[k_1^-]{k_1^+}B$ , $B\xrightleftharpoons[k_2^-]{k_2^+}C$ , $C\xrightleftharpoons[k_3^-]{k_3^+}A$ , the three changes sum to zero. Compatibility requires $K_1K_2K_3=1$ . This is exactly equality of forward and reverse rate products.

A stoichiometric dependency need not look like a loop of reaction complexes. Consider the hypothetical elementary pairs $A\xrightleftharpoons[k_1^-]{k_1^+}B$ , $C\xrightleftharpoons[k_2^-]{k_2^+}D$ , and $A+C\xrightleftharpoons[k_3^-]{k_3^+}B+D$ . The third change is the sum of the first two, so $z=(1,1,-1)^\top$ lies in the kernel. Thus $K_3=K_1K_2$ , although the complex graph consists of three separate reversible pairs.

Checkpoint: choose compatible constants.

If $K_1=2$ and $K_2=3$ , choose $K_3=6$ , not 5. One positive detailed-balanced composition is $A=C=c^\circ$ , $B=2c^\circ$ , $D=3c^\circ$ . Its third concentration ratio is $BD/(AC)=6$ . With 5, no positive composition can balance all three pairs.

Test equilibrium constants against the stoichiometric kernel, including dependencies invisible in a drawn reaction-complex graph.

Part III

Calculate regulation from states and weights

Count the allowed configurations, include concentration and internal free energy, normalize, and only then attach the measured activity.

13Derive the concentration factor by counting binding arrangements

Binding competes with the many places a free molecule could occupy. To make this entropy cost explicit, consider one gene site and $N$ identical repressors distributed over $M$ equivalent background locations. Each location accommodates at most one repressor. The free-gene macrostate has $\binom MN$ arrangements. The bound-gene macrostate leaves $N-1$ repressors on those locations and has $\binom M{N-1}$ arrangements.

The multiplicity ratio supplies the missing concentration dependence. If moving one repressor from the background to the gene changes its energy by $\Delta\epsilon$ , then

\begin{aligned} \frac{w_{\mathrm{bound}}}{w_{\mathrm{free}}} &=\frac{\binom M{N-1}}{\binom MN}e^{-\Delta\epsilon/(k_BT)} =\frac{N}{M-N+1}e^{-\Delta\epsilon/(k_BT)}\\ &\simeq\frac NM e^{-\Delta\epsilon/(k_BT)} \qquad (N\ll M). \end{aligned}

The last approximation is dilute occupancy of the background. One extra molecule changes the reservoir little only when the available pool is also sufficiently large. These are separate assumptions. The lattice illustrates the entropy calculation, rather than asserting that cytoplasmic proteins literally occupy a fixed grid. Its reference volume and energy convention are absorbed into the physical standard binding free energy.⁵

A declared reference volume converts the lattice ratio into concentration units. Give each background site volume $v_{\mathrm{site}}$ , so $M=V/v_{\mathrm{site}}$ . In the free-gene state, molar repressor concentration is $R=N/(N_AV)$ . With standard concentration $c^\circ=1/(N_Av_{\mathrm{site}})$ , the ratio becomes $N/M=R/c^\circ$ . Choosing a different reference volume changes the corresponding standard binding free energy, not the physical affinity. This connects the counting model to the dimensionless concentration factor derived from chemical potential in section 11.

The physical weight ratio is free concentration divided by affinity. Using the standard-state relation from section 11 gives $w_G=1$ and $w_{C_{GR}}=R/K_d$ . Thus

Z=1+\frac R{K_d},\qquad p_{\mathrm{bound}}=\frac{R/K_d}{1+R/K_d}=\frac R{K_d+R},\qquad p_{\mathrm{free}}=\frac{K_d}{K_d+R}.

The kinetic derivation agrees when it uses the same assumptions. For elementary $G+R\xrightleftharpoons[k_-]{k_+}C_{GR}$ , detailed balance gives $k_+GR=k_-C_{GR}$ . Hence $C_{GR}/G=R/K_d$ , with $K_d=k_-/k_+$ . Dividing by the gene total $q_G=G+C_{GR}$ reproduces the probability formula in an ensemble description. The weight method needs the compatible state free energies, but not the separate kinetic constants.

Scroll horizontally for the full-size figure.

Figure 8. Entropy and binding free energy determine occupancy. The upper panels compare molecular arrangements. Their ratio becomes the concentration factor in the lower state weights. Bound probability and the activity of a promoter that is silent when bound are different readouts of those same probabilities.

A single gene is never fractionally occupied on one realization. It is free or bound. Its occupancy probability equals the average bound fraction across repeated experiments or an ensemble of equivalent genes at fixed free repressor. Replacing free $R$ by total $q_R$ requires negligible depletion. Small isolated pools need a joint count distribution, developed in section 18.

Checkpoint: compare two binding affinities.

At $R=10\,\mathrm{nM}$ , an affinity $K_d=10\,\mathrm{nM}$ gives bound probability 1/2. A tenfold smaller $K_d$ gives 10/11. At the same standard state, this affinity change corresponds to $\Delta\Delta G^\circ=-k_BT\ln10$ . It does not determine whether association became faster or dissociation became slower.

Derive a binding curve from multiplicity and free energy, while keeping free concentration, total concentration and single-gene probability distinct.

14Attach activities before predicting expression

Occupancy becomes an output only after a readout rule is specified. Let $r_i$ be the transcription-initiation frequency of one gene in configuration $i$ . These frequencies have inverse-time units. With state probabilities $p_i$ , the ensemble mean initiation frequency is

\bar r=\sum_i r_ip_i.

This expectation is exact for the stated conditional readout model. Replacing a rapidly fluctuating promoter by a single averaged initiation hazard over slow times is an additional time-scale approximation. The resulting mean law does not automatically reproduce burst statistics or short-time correlations.

Simple repression is a two-state activity assignment. Give the free state frequency $r_0$ and the bound state frequency $r_b$ . At binding equilibrium,

\bar r(R)=\frac{r_0K_d+r_bR}{K_d+R}. \qquad r_b=0\ \Longrightarrow\ \frac{\bar r(R)}{\bar r(0)}=\frac1{1+R/K_d}.

Complete silencing is one choice, not a consequence of binding. If the bound state retains activity, the high-repressor limit is $r_b$ . If the two frequencies are equal, binding has no effect on this output despite changing occupancy.

Protein abundance requires another dynamical bridge. Let $M$ and $X$ be mRNA and protein concentrations. Hold gene concentration $q_G$ and translation frequency $k_{\mathrm{tr}}$ fixed. Let $\gamma_M,\gamma_X$ be effective first-order removal constants, including dilution if the cell grows. The composite transcription and translation processes have separately specified rate laws:

\begin{aligned} \varnothing&\rightsquigarrow M,\qquad M\rightsquigarrow M+X,\\ \dot M&=q_G\bar r-\gamma_MM,\qquad \dot X=k_{\mathrm{tr}}M-\gamma_XX,\\ M^*&=\frac{q_G\bar r}{\gamma_M},\qquad X^*=\frac{k_{\mathrm{tr}}q_G\bar r}{\gamma_M\gamma_X}. \end{aligned}

Each equation has concentration/time units. A steady protein fold-change equals an initiation fold-change only if the other factors remain unchanged and the reporter has reached its steady state. Growth-dependent dilution, changing translation or slow maturation can break that equality without invalidating the promoter's equilibrium occupancy model. The occupancy-to-expression assumption is explicit in the thermodynamic promoter framework.⁵

Checkpoint: does a bound repressor necessarily repress?

No. Binding determines which state is occupied. Repression means that the occupied state has a lower relevant activity. A ligand can instead activate, sequester an inhibitor or alter a downstream transition. The state ledger and the readout column identify these different models explicitly.

Derive an effective production law and state the additional conditions needed to compare it with a protein-expression measurement.

15Enumerate joint states before combining regulatory inputs

Two single-input curves do not specify a two-input response. Suppose free regulators $A$ and $B$ can bind different sites on one gene. Define the dimensionless free concentrations $\bar A=A/K_A$ and $\bar B=B/K_B$ . If simultaneous binding is allowed, the state ledger is

Gene configuration	Weight	Initiation frequency
Neither bound	1	$r_0$
Only A bound	$\bar A$	$r_A$
Only B bound	$\bar B$	$r_B$
Both bound	$\omega\bar A\bar B$	$r_{AB}$

The positive interaction factor $\omega=e^{-\Delta G_{\mathrm{int}}/(k_BT)}$ accounts for the extra free energy of joint binding. Conditional independence gives $\omega=1$ . Steric exclusion removes the joint state entirely. A very small weight approximates exclusion, but is not the same model.

Normalization and activity averaging give the response directly.

Z=1+\bar A+\bar B+\omega\bar A\bar B,\qquad \bar r=\frac{r_0+r_A\bar A+r_B\bar B+r_{AB}\omega\bar A\bar B}{Z}.

For independent binding with activity only in the jointly bound state, $\bar r/r_{AB}=[\bar A/(1+\bar A)][\bar B/(1+\bar B)]$ . For activity whenever either regulator is bound, with the same activity in all three occupied states, the normalized output is $1-1/[(1+\bar A)(1+\bar B)]$ . These product and OR-like forms follow from different activity columns, not a universal rule for combining Hill curves.

Scroll horizontally for the full-size figure.

Figure 9. The state ledger fixes the mixed response. Along the matched-input path

\bar A=\bar B

, the curves compare independent AND-like activity, independent OR-like activity, cooperative joint-only activity with

\omega=10

, and activity under mutual exclusion. These are declared state and activity choices, not fitted logic functions.

Checkpoint: infer an unmeasured interaction.

Take joint-only activity. If either input is zero, the output is zero for every value of $\omega$ . Single-input measurements therefore cannot identify the joint interaction. At $\bar A=\bar B=1$ , the output is $\omega/(3+\omega)$ . A mixed-input measurement supplies information the isolated input curves cannot contain.

Build a multiregulator response from allowed joint states and activities instead of guessing how isolated response curves combine.

16Sum an allosteric ensemble without specifying every rate

An MWC model couples ligand binding through a shared conformation. Consider a generic protein with active conformation A and inactive conformation I. All $n$ sites switch conformation together. Within one conformation the sites are equivalent and bind independently. Their dissociation constants are $K_A$ and $K_I$ . Let $L>0$ be the unliganded inactive/active weight ratio. Thus $L=e^{-[G_I^0-G_A^0]/(k_BT)}$ . This defines a model class, not a claim that every allosteric protein has this mechanism.⁶

Enumerate the two-site case before taking the binomial shortcut. In either conformation, two sites have patterns 00, 10, 01 and 11. There are eight molecular states. The two singly occupied patterns are distinct but have equal weights. Grouping by occupancy gives six rows:

Conformation	Bound ligands	Multiplicity	Total class weight
A	0	1	1
A	1	2	$2c/K_A$
A	2	1	$(c/K_A)^2$
I	0	1	$L$
I	1	2	$2Lc/K_I$
I	2	1	$L(c/K_I)^2$

Here $c$ is free ligand concentration. A class containing $m$ bound ligands has $\binom nm$ distinct site patterns. Summing their weights gives

\begin{aligned} W_A&=\sum_{m=0}^n\binom nm(c/K_A)^m=(1+c/K_A)^n,\\ W_I&=L\sum_{m=0}^n\binom nm(c/K_I)^m=L(1+c/K_I)^n,\\ Z&=W_A+W_I,\qquad p_A=\frac{W_A}{Z}. \end{aligned}

Activity and bound fraction require different numerators. If every A state is active with the same catalytic frequency and every I state is inactive, $p_A$ is normalized activity. The mean fraction of occupied sites instead weights each state by its bound count divided by $n$ . Differentiating each weight with respect to $\ln c$ brings down that count:

\theta=\frac{\langle N_{\mathrm{bound}}\rangle}{n} =\frac{c}{nZ}\frac{dZ}{dc} =\frac{(c/K_A)(1+c/K_A)^{n-1}+L(c/K_I)(1+c/K_I)^{n-1}}{Z}.

The differentiation shortcut works because each weight is proportional to $c^m$ for a state with $m$ ligands. A fluorescence occupancy assay and an activity assay can therefore report different curves for the same protein.

The limiting behavior can be checked without fitting. At zero ligand, $p_A=1/(1+L)$ and $\theta=0$ . At saturating ligand, $p_A\to1/[1+L(K_A/K_I)^n]$ while $\theta\to1$ . Use the core's two-site example throughout: $L=10,K_A=1\,\mathrm{nM},K_I=100\,\mathrm{nM}$ . Activity ranges from 1/11 to 1000/1001, while site occupancy ranges from zero to one.

One concentration makes the different numerators concrete. At $c=1\,\mathrm{nM}$ , the conformation weights are $W_A=4$ and $W_I=10.201$ . Hence $p_A=4/14.201\approx0.282$ . The occupied-site numerator is $2+0.101=2.101$ , giving $\theta=2.101/14.201\approx0.148$ . The distinction comes from the activity assigned to each state, not from a different partition function.

Scroll horizontally for the full-size figure.

Figure 10. An activity curve, a binding curve and a slope are different observables. Left: active fraction and mean occupied-site fraction for the two-site example in the text. Right: the local slope of active-state log odds, derived in section 17. This slope function is independent of

L

and is not a fitted Hill exponent.

Compose allostery with repression to predict induction

A repressor turns the active fraction into a new biological input. For LacI, “active” means competent to bind the operator and repress. An inducer such as IPTG favors the inactive conformation, so $K_I<K_A$ . Increasing inducer lowers $p_A(c)$ and can raise gene expression. This reverses the ligand preference in the illustrative activating-protein example above.

Count the available active repressors before assigning the DNA weight. Let $N_R$ be the number of repressor dimers and $N_{\mathrm{NS}}$ the number of nonspecific genomic binding sites. Let $\Delta\epsilon_R$ be the active repressor's specific-minus-nonspecific binding free energy. Under dilute background occupancy and negligible depletion by the operator,

w_R^0=\frac{N_R}{N_{\mathrm{NS}}}e^{-\Delta\epsilon_R/(k_BT)}, \qquad w_R(c)=p_A(c)w_R^0.

The count ratio is dimensionless. It is not a molar concentration divided by a site count. The state weights of section 13 supply this background-count formulation. Here we also neglect specific DNA binding by the inactive conformation.

Polymerase competition supplies the final readout. Allow three promoter states: empty, polymerase-bound and active-repressor-bound. Exclude simultaneous polymerase/repressor binding. Their weights are $1,w_P,w_R(c)$ . Divide polymerase occupancy by its value in a matched no-repressor strain:

\mathrm{fold\ change}(c) =\frac{w_P/[1+w_P+w_R(c)]}{w_P/(1+w_P)} =\frac{1+w_P}{1+w_P+w_R(c)} \simeq\frac1{1+p_A(c)w_R^0}.

The final approximation uses $w_P\ll1$ . Identifying this ratio with measured protein expression also requires the unchanged downstream factors discussed in section 14. The allosteric state sum, DNA occupancy and reporter dynamics are separate pieces of the prediction.

Scroll horizontally for the full-size figure.

Figure 10A. From inducer to expression in three explicit steps. Curves reconstruct the central allosteric-induction model of Razo-Mejia and colleagues for O2, using

K_A=139\,\mu\mathrm M

K_I=0.53\,\mu\mathrm M

L=e^{-4.5}

, two sites,

\Delta\epsilon_R=-13.9\,k_BT

and

N_{\mathrm{NS}}=4.6\times10^6

. They are not measured data or refitted curves. The original uncertainty distributions are not propagated.¹¹

Checkpoint: follow one prediction through its units.

At $c=10\,\mu\mathrm M$ and $N_R=260$ , first compute $p_A=0.2076$ . The unmodulated DNA weight is $w_R^0=\frac{260}{4.6\times10^6}e^{13.9}\approx61.50$ . Thus $w_R=12.77$ and fold-change is $1/(1+12.77)=0.0726$ . Repressor number counts dimers. Substituting tetramer counts without the factor of two would change the prediction.

The physical parameters specify useful perturbations. Changing inducer changes $p_A$ . Changing repressor abundance changes $N_R$ . Changing operator sequence changes $\Delta\epsilon_R$ . A new dose-response curve need not receive a new arbitrary Hill fit. The lineage explains how prior calibration and one fitted induction curve supported tests across other constructs.

Enumerate an allosteric state ensemble and compute both conformation activity and ligand occupancy from the correct weighted sums.

17Interpret a Hill-like slope without inventing a mechanism

A slope needs a named readout and a named normalization. For a probability $p(c)$ strictly between zero and one, define its local log-odds slope $h_p(c)=d\ln[p/(1-p)]/d\ln c$ . A Hill curve $p=c^h/(K^h+c^h)$ has constant slope $h$ . A general molecular-state model need not.

The MWC active-state slope follows from its odds ratio. Using $p_A/(1-p_A)=W_A/W_I$ , take logarithms before differentiating:

\begin{aligned} \ln\frac{p_A}{1-p_A} &=-\ln L+n\ln(1+c/K_A)-n\ln(1+c/K_I),\\ h_A(c)&=n\left[\frac{c}{K_A+c}-\frac{c}{K_I+c}\right]. \end{aligned}

If $K_A<K_I$ , this slope is positive and smaller than $n$ at every finite positive concentration. It tends to zero at both extremes. Its maximum occurs at $c=\sqrt{K_AK_I}$ : differentiating the bracket and setting it to zero gives $K_A/(K_A+c)^2=K_I/(K_I+c)^2$ . For the two-site example, the maximum is $18/11\approx1.64$ , not 2.

Changing the conformation bias shifts the response without changing this slope function. The constant $L$ disappears from the logarithmic derivative, but it determines which concentration corresponds to half activity. A fit restricted to the transition region can therefore report different effective sharpness when the bias changes. Renormalizing activity between its nonzero lower and upper plateaus also changes the odds slope. State the convention before comparing numbers.

A measured Hill coefficient is not a direct count of binding sites. Site number, affinity contrast, conformational bias, hidden states and the selected output all affect the fitted curve. The specific bound above concerns the raw active-state odds in this MWC model. It is not a theorem about all possible biological response functions.⁶

Checkpoint: remove the affinity difference.

If $K_A=K_I=K$ , then $p_A=1/(1+L)$ is independent of ligand, even with many sites. Yet $\theta=c/(K+c)$ still rises from zero to one. Ligand binding alone does not imply allosteric regulation.

Calculate and interpret a response slope while separating site number, observable choice and fitted sharpness.

18Close a binding model when the ligand pool is finite

Equilibrium weights use free ligand, while many experiments control totals. For the binding pair $G+R\xrightleftharpoons[k_-]{k_+}C_{GR}$ , define $q_G=G+C_{GR}$ and $q_R=R+C_{GR}$ . At deterministic binding equilibrium, substituting these exact balances into $K_dC_{GR}=GR$ gives

(q_G-C_{GR})(q_R-C_{GR})=K_dC_{GR}.

The physical interval selects one quadratic root. Complex concentration must lie between zero and $\min(q_G,q_R)$ . The left side minus the right side decreases strictly across this interval, since its derivative is $-(q_G+q_R+K_d-2C_{GR})<0$ . Its endpoint signs bracket zero. The unique physical root, written in a numerically stable form, is

C_{GR}^{\mathrm{eq}}= \frac{2q_Gq_R}{q_G+q_R+K_d+ \sqrt{(q_G+q_R+K_d)^2-4q_Gq_R}}.

The discriminant is $(q_G-q_R)^2+2K_d(q_G+q_R)+K_d^2>0$ . The second root exceeds both totals. With $q_G=10,q_R=1,K_d=1\,\mathrm{nM}$ , the physical complex is approximately 0.901 nM. Only 9.01% of genes are bound. Substituting total repressor into the buffered formula would give 50%, violating the available-pool constraint.

Scroll horizontally for the full-size figure.

Figure 11. Depletion and count fluctuations are separate corrections. Left: bound gene fraction versus total repressor for gene total 10 nM and

K_d=1\,\mathrm{nM}

. The naive buffered curve uses total in place of free repressor. Right: exact count-model mean bound fraction for equal fixed totals

N_G=N_R=m

and

\Omega K_d=m

. The deterministic root is approached as volume and both totals increase together.

A finite-count equilibrium needs probabilities over the bound count. Let fixed constituent counts be $N_G,N_R$ . Let $N_{C_{GR}}=c$ be the random number of bound pairs, with $0\leq c\leq\min(N_G,N_R)$ . The count propensities are

a_+(c)=\frac{k_+}{\Omega}(N_G-c)(N_R-c),\qquad a_-(c)=k_-c,\qquad\Omega=N_AV.

Detailed balance gives an exact recurrence. Write the stationary bound-count distribution as $\pi_c$ . Equating the two fluxes between counts $c$ and $c+1$ yields

\frac{\pi_{c+1}}{\pi_c} =\frac{(N_G-c)(N_R-c)}{\Omega K_d(c+1)},\qquad \pi_c\propto \frac{N_G!N_R!}{(N_G-c)!(N_R-c)!\,c!\,(\Omega K_d)^c}.

Normalize this finite list of weights. The exact mean need not equal the deterministic root. For $N_G=N_R=1$ and $\Omega K_d=1$ , the two weights are equal and $\langle N_{C_{GR}}\rangle=1/2$ . The deterministic equation in count units gives $(3-\sqrt5)/2\approx0.382$ . The discrepancy comes from fluctuations and correlations, not slow equilibration.

Checkpoint: can faster binding remove this discrepancy?

No. Multiplying $k_+$ and $k_-$ together leaves the stationary recurrence unchanged. Increasing system size at fixed concentrations is a different limit. For this example, set $N_G=N_R=m$ and $\Omega K_d=m$ . The exact mean divided by $m$ approaches the deterministic fraction as $m$ grows.

Close an equilibrium model with finite totals and choose a count distribution when molecule-number fluctuations cannot be neglected.

Part IV

Use equilibrium as a tested fast-state reduction

Preserve exact accounting, specify the eliminated clock and check whether the chosen state and channel resolution remains adequate.

19Separate exact total dynamics from approximate reconstruction

Binding can be eliminated from total equations before making any approximation. Return to the reversible enzyme in section 8. Define net binding fluxes $v_{\mathrm{bind},S}=k_{+S}ES-k_{-S}C_{ES}$ and $v_{\mathrm{bind},P}=k_{+P}EP-k_{-P}C_{EP}$ . Define catalytic flux $v_{\mathrm{cat}}=k_{\mathrm{cat}}^+C_{ES}-k_{\mathrm{cat}}^-C_{EP}$ . All have concentration/time units.

Write the full equations to see the cancellation. Let composite supply and removal maintain substrate and drain product, with separately declared concentration fluxes $v_{\mathrm{in}}$ and $v_{\mathrm{out}}$ . At fixed volume,

\begin{aligned} \dot E&=-v_{\mathrm{bind},S}-v_{\mathrm{bind},P},\\ \dot C_{ES}&=v_{\mathrm{bind},S}-v_{\mathrm{cat}},& \dot C_{EP}&=v_{\mathrm{bind},P}+v_{\mathrm{cat}},\\ \dot S&=v_{\mathrm{in}}-v_{\mathrm{bind},S},& \dot P&=-v_{\mathrm{bind},P}-v_{\mathrm{out}}. \end{aligned}

Count the enzyme, unconverted substrate and product constituents. Their totals are $q_E=E+C_{ES}+C_{EP}$ , $q_S=S+C_{ES}$ and $q_P=P+C_{EP}$ . Adding the appropriate species equations gives the exact projection

\dot q_E=0,\qquad \dot q_S=v_{\mathrm{in}}-v_{\mathrm{cat}},\qquad \dot q_P=v_{\mathrm{cat}}-v_{\mathrm{out}}.

Only $q_E$ is conserved by the entire stated network. The substrate and product totals are conserved by binding but changed by catalysis and the boundary. This is why they are natural slow coordinates when binding is fast.

The fast binding problem has a different graph from the complete enzyme cycle. Suppose both binding pairs relax before catalysis or boundary fluxes appreciably change the retained totals. The leading fast graph contains the two binding pairs, but excludes the slow catalytic edge. Its equilibrium ratios are $C_{ES}/E=S/K_{d,S}$ and $C_{EP}/E=P/K_{d,P}$ , where $K_{d,S}=k_{-S}/k_{+S}$ and $K_{d,P}=k_{-P}/k_{+P}$ . Therefore

\begin{aligned} Z_{\mathrm{bind}}&=1+S/K_{d,S}+P/K_{d,P},\\ C_{ES}^{\mathrm{eq}}&=q_E\frac{S/K_{d,S}}{Z_{\mathrm{bind}}},& C_{EP}^{\mathrm{eq}}&=q_E\frac{P/K_{d,P}}{Z_{\mathrm{bind}}},\\ v_{\mathrm{cat}}^{\mathrm{red}} &=q_E\frac{k_{\mathrm{cat}}^+S/K_{d,S}-k_{\mathrm{cat}}^-P/K_{d,P}}{Z_{\mathrm{bind}}}. \end{aligned}

Finite totals close these formulas implicitly. Solve $S=q_S-C_{ES}^{\mathrm{eq}}$ and $P=q_P-C_{EP}^{\mathrm{eq}}$ together with the displayed weights for nonnegative free concentrations. Denote the resulting complex reconstruction by $h_{C_{ES}}(q_S,q_P;q_E,K_{d,S},K_{d,P})$ and likewise for $C_{EP}$ . Replacing the exact complexes by these functions closes the slow equations approximately. Free-ligand dominance can simplify the solve further, but is not part of exact total accounting.

Scroll horizontally for the full-size figure.

Figure 12. Exact projection and equilibrium closure do different jobs. Binding cancels exactly from the total equations. Compatible fast binding supplies an approximate reconstruction of the molecular states. Slow catalysis can maintain throughput while that fast layer remains close to its conditional equilibrium. No equilibrium claim about the whole cell is needed.

Checkpoint: does a reconstructed complex stop changing?

No. Along the slow solution, $d h_{C_{ES}}/dt=(\partial h_{C_{ES}}/\partial q_S)\dot q_S+(\partial h_{C_{ES}}/\partial q_P)\dot q_P$ , with other parameters fixed. At finite binding speed, a small displacement from the algebraic constraint supplies the binding flux needed to follow that motion. Setting the limiting binding balance to zero is not an exact statement that every finite-speed complex derivative vanishes.

Construct exact slow-coordinate equations and identify precisely where equilibrium reconstruction enters as an approximation.

20Bound equilibrium-reconstruction error for a named observable

A relaxation time must be compared with the motion of the distribution being followed. Let a two-state binding model have time-dependent hazards $a(t),b(t)>0$ . Define $\kappa(t)=a(t)+b(t)$ and the instantaneous stationary bound probability $p_{\mathrm{eq}}(t)=a(t)/\kappa(t)$ . The exact occupancy equation is

\dot p=\kappa(t)[p_{\mathrm{eq}}(t)-p].

The error equation separates initial relaxation from ongoing lag. Set $e=p-p_{\mathrm{eq}}$ . Then $\dot e=-\kappa e-\dot p_{\mathrm{eq}}$ . Multiplying by the integrating factor and integrating gives

e(t)=e(0)e^{-\int_0^t\kappa(u)\,du} -\int_0^t e^{-\int_s^t\kappa(u)\,du}\dot p_{\mathrm{eq}}(s)\,ds.

If $\kappa(t)\geq\kappa_{\min}>0$ and $|\dot p_{\mathrm{eq}}(t)|\leq M$ , bound the exponential kernel by $e^{-\kappa_{\min}(t-s)}$ and integrate. This proves

|e(t)|\leq |e(0)|e^{-\kappa_{\min}t} +\frac{M}{\kappa_{\min}}(1-e^{-\kappa_{\min}t}).

The first term is the initial layer. The second is the lag caused by the moving target. Both $M$ and $\kappa_{\min}$ have inverse-time units, so their ratio is an absolute probability-error bound. Small absolute error need not be small relative error for a rare state.

A sinusoid gives a reproducible tracking example. Prescribe $p_{\mathrm{eq}}(t)=1/2+(2/5)\sin(\omega t)$ . Set both hazards through $a=\kappa p_{\mathrm{eq}}$ , $b=\kappa(1-p_{\mathrm{eq}})$ , with constant $\kappa$ , and start at $p(0)=1/2$ . Direct integration gives

p(t)=\frac12+\frac{2}{5(\kappa^2+\omega^2)} \left[\kappa^2\sin(\omega t)-\kappa\omega\cos(\omega t) +\kappa\omega e^{-\kappa t}\right].

This example controls both hazards. It is not fixed-association-constant binding driven solely by a changing ligand concentration. After the transient, the response amplitude is reduced by $\kappa/\sqrt{\kappa^2+\omega^2}$ and its phase lags by $\arctan(\omega/\kappa)$ . The bound uses $M=(2/5)\omega$ .

Scroll horizontally for the full-size figure.

Figure 13. Faster relaxation reduces the lag, not the motion of the target. The instantaneous equilibrium is the dashed sinusoid with angular frequency 1/s. Solid curves show exact solutions for relaxation hazards 0.5/s and 10/s. Both start on the equilibrium curve. Continued motion creates an error even after the initial condition has been forgotten.

Checkpoint: certify a requested tolerance.

Suppose $|\dot p_{\mathrm{eq}}|\leq0.02\,\mathrm{s}^{-1}$ and $\kappa\geq10\,\mathrm{s}^{-1}$ . After the initial transient, the absolute error is at most 0.002. To guarantee this from an arbitrary initial occupancy to within an additional 0.001, wait at least $\ln(1000)/10\approx0.69\,\mathrm{s}$ , using the worst-case initial error of one.

Calculate the error of a rapid-binding enzyme closure

The running enzyme makes the observable dependence explicit. Use the symmetric catalytic family in section 8 at fixed free reservoirs. With $Z=1+\bar S+\bar P$ , its rapid-binding distribution is $p^{\mathrm{bind}}=(1,\bar S,\bar P)^\top/Z$ . This is the equilibrium of the binding-only graph, not of the full catalytic cycle.

The bound-state balance determines the displacement from that graph. The substrate-bound state has zero net accumulation when

0=k_{\mathrm{off}}(\bar S p_E-p_{ES}) -k_{\mathrm{cat}}(p_{ES}-p_{EP}).

Its first term is the binding current. Its second is the catalytic current $J$ . Since the exact free probability is $p_E=1/Z$ , define the probability displacement $\delta=J/k_{\mathrm{off}}$ . Normalization then gives $p_{ES}=\bar S/Z-\delta$ and $p_{EP}=\bar P/Z+\delta$ . Substituting these into the catalytic current solves the displacement:

\delta=\frac{k_{\mathrm{cat}}(\bar S-\bar P)}{(k_{\mathrm{off}}+2k_{\mathrm{cat}})Z}, \qquad J=k_{\mathrm{off}}\delta, \qquad J_{\mathrm{bind}}=k_{\mathrm{cat}}\frac{\bar S-\bar P}{Z}.

A small displacement can carry a finite flux. At fixed reservoirs and catalytic frequency, $\delta\to0$ as binding speeds up, but $J\to J_{\mathrm{bind}}$ . Multiplying a small probability imbalance by a large hazard produces the finite current. Setting the leading-order binding balance to zero therefore does not set the finite-speed binding flux to zero.

Different measurements divide the discrepancy by different quantities. Total variation, $\epsilon_{\mathrm{pop}}=\tfrac12\sum_i|p_i-p_i^{\mathrm{bind}}|$ , is an absolute distribution error. For nonzero net driving, define relative current error and relative product-bound-state error using the reduced prediction in their denominators. They are

\begin{aligned} \epsilon_{\mathrm{pop}}&=|\delta|,\\ \epsilon_J&=\frac{|J-J_{\mathrm{bind}}|}{|J_{\mathrm{bind}}|} =\frac{2k_{\mathrm{cat}}}{k_{\mathrm{off}}+2k_{\mathrm{cat}}},\\ \epsilon_{EP}&=\frac{|p_{EP}-p_{EP}^{\mathrm{bind}}|}{p_{EP}^{\mathrm{bind}}} =\frac{|\delta|Z}{\bar P}. \end{aligned}

Scroll horizontally for the full-size figure.

Figure 13A. Name the error before deciding whether binding is fast enough. The maintained inputs are

\bar S=1

and

\bar P=0.01

. Faster binding reduces all three errors at these fixed inputs. The same rate ratio gives very different absolute-population, relative-current and relative-rare-state tolerances. These are exact stationary comparisons, not a finite-time simulation.

Checkpoint: audit a 100-fold separation.

Take $k_{\mathrm{off}}/k_{\mathrm{cat}}=100$ and the plotted reservoirs. Rapid binding gives $p_{EP}^{\mathrm{bind}}=0.004975$ . The displacement is $\delta=0.004829$ , giving exact $p_{EP}=0.009804$ . The population error is below 0.5 percentage points, and the current error is 1.96%. But a readout active only in the product-bound state is 97.1% above its reduced prediction. A small global probability error is not a small relative error for every state-specific output.

The tolerance must include the driving regime. At fixed $\mathcal A=\ln(\bar S/\bar P)$ , every error vanishes in the infinite-separation limit. For a rare product state under strong forward driving, $\epsilon_{EP}=|e^{\mathcal A}-1|/(k_{\mathrm{off}}/k_{\mathrm{cat}}+2)$ . A uniform claim over arbitrary driving would require more than a fixed large rate ratio. At zero driving both currents vanish, so a relative current error is undefined even though the distributions agree exactly.

Turn “fast binding” into a claim with a stated clock, driving regime, observable and error tolerance.

21Distinguish hidden dissipation from hidden memory

Resolving states and resolving transition channels are separate tasks. Two physical mechanisms can connect the same pair of observed states. If the mechanisms exchange different chemical species or different amounts of fuel, their thermodynamic reverse partners must be tracked separately. Summing their hazards may preserve the observed Markov dynamics while concealing physical circulation.³

Two states already suffice to hide a driven cycle when there are two channels. Let channel α have hazards 2/s from 0 to 1 and 1/s from 1 to 0. Let channel β have hazards 1/s and 2/s in the same directions. The aggregate graph has hazards 3/s both ways and stationary probabilities $p_0^*=p_1^*=1/2$ . Its two resolved currents are

J_\alpha=\frac12(2)-\frac12(1)=\frac12\,\mathrm{s}^{-1}, \qquad J_\beta=-\frac12\,\mathrm{s}^{-1}.

The cycle goes forward through α and back through β. Its affinity is $\ln4$ . Under a physical local-detailed-balance assignment, its stationary entropy production is $k_B(\ln2)\,\mathrm{s}^{-1}$ . The aggregate state path has exactly the same law as a genuinely equilibrium single-channel chain with hazards 3/s each way. No statistic of that state path alone can distinguish the two channel assignments.

Scroll horizontally for the full-size figure.

Figure 14. Two different omissions require two different tests. Left: summing physical channels conceals equal and opposite stationary currents while leaving a Markov state path. Right: pooling slow hidden states can produce a non-exponential holding-time law. Channel-resolved chemistry tests the first omission. Dwell-time statistics can test the second.

Pooling hidden states can instead destroy the Markov property. The two-exponential OFF survival from section 1 retains the age of an OFF episode as predictive information. A sufficient condition for exact aggregation at all initial distributions is that every microscopic state in one block has the same total transition hazard into each other block. This is the strong lumpability condition. Fast internal mixing can sometimes justify an approximate aggregate model, but requires a separate time-scale analysis.

Checkpoint: which experiment is needed?

Non-exponential dwell times challenge a fixed-rate state description. Exponential dwell times do not establish physical equilibrium, because hidden parallel channels can preserve them exactly. To test those channels, measure fuel turnover, labeled molecular transfers or mechanism-specific event directions.

Choose a test that matches the missing information: state memory and physical channel dissipation are not the same defect.

22Build a gated-receptor model from the specification

A fresh model tests whether the construction can be transferred. Consider a receptor H that switches between closed $H_{\mathrm c}$ and open $H_{\mathrm o}$ configurations. Only the open configuration binds ligand L. The bound state is $C_{HL}$ . Begin with free ligand concentration $c$ held fixed. The observed outputs are bound fraction, mean activity and time to first binding from a closed receptor.

Translate the specification into a state graph before calculating. The closed-to-open hazard is $\alpha$ , its reverse is $\beta$ , open-to-bound is $k_+c$ , and bound-to-open is $k_-$ . These are state-transition hazards. The corresponding chemical binding pair is elementary $H_{\mathrm o}+L\xrightleftharpoons[k_-]{k_+}C_{HL}$ . In state order closed, open, bound,

Q=\begin{pmatrix} -\alpha&\beta&0\\ \alpha&-(\beta+k_+c)&k_-\\ 0&k_+c&-k_- \end{pmatrix}.

The tree graph supplies compatible weights. Assuming these are the resolved passive channels, choose weights $1,\alpha/\beta,(\alpha/\beta)c/K_d$ , with $K_d=k_-/k_+$ . Normalizing gives

p_{\mathrm{bound}}= \frac{c}{K_{\mathrm{eff}}+c},\qquad K_{\mathrm{eff}}=K_d\left(1+\frac{\beta}{\alpha}\right).

The apparent affinity includes the closed-state penalty. It need not equal the microscopic affinity of the open receptor. If only the bound receptor produces a nonperturbing readout at frequency $r_{\mathrm{bound}}$ , its mean activity is $\bar r=r_{\mathrm{bound}}p_{\mathrm{bound}}$ . This readout rule does not follow from the weights.

The kinetic observable requires the rates that the weight ratio discarded. Apply the first-passage derivation with opening $\alpha$ , reclosing $\beta$ and binding $k_+c$ . The answers are

\tau_{\mathrm{open}}=\frac{1+\beta/\alpha}{k_+c},\qquad \tau_{\mathrm{closed}}=\frac1\alpha+\frac{1+\beta/\alpha}{k_+c}.

Use $\alpha=2,\beta=1,k_-=4\,\mathrm{s}^{-1}$ , $k_+=1\,\mathrm{nM}^{-1}\mathrm{s}^{-1}$ and $c=4\,\mathrm{nM}$ . The weights are (1, 2, 2), so probabilities are (0.2, 0.4, 0.4). Mean first binding takes 0.375 s from open and 0.875 s from closed. With $r_{\mathrm{bound}}=0.1\,\mathrm{s}^{-1}$ , mean activity is 0.04 events/s per receptor.

Scroll horizontally for the full-size figure.

Figure 15. One graph predicts both occupancy and delay. The passive gated-receptor graph is a line. Its effective binding curve is shifted by the closed-state weight. The two first-arrival curves depend on kinetic hazards, not just the apparent affinity. All parameters are those declared in the worked example.

Finite ligand total adds one conservation equation. Let receptor total be $q_H$ and ligand total be $q_L$ . Bound concentration is $C_{HL}=q_Hp_{\mathrm{bound}}$ and free ligand is $c=q_L-C_{HL}$ . Therefore

(q_H-C_{HL})(q_L-C_{HL})=K_{\mathrm{eff}}C_{HL}.

The same physical-root calculation as section 18 now applies with $K_{\mathrm{eff}}$ . For $q_H=q_L=6\,\mathrm{nM}$ and the rates above, $K_{\mathrm{eff}}=6\,\mathrm{nM}$ . The bound fraction is $(3-\sqrt5)/2\approx0.382$ . Treating total ligand as a reservoir would instead predict 1/2.

Checkpoint: speed every transition by ten.

Multiply $\alpha,\beta,k_+,k_-$ by ten, keeping ligand fixed. Both affinity ratios and all occupancies remain unchanged. Both first-arrival times become tenfold shorter. If the measured delay does not change, the experiment has not sped the rate-limiting processes assumed by this model.

Construct a new molecular-state model, separate its static and kinetic predictions, and close it with finite-pool conservation.

Part V · Optional

Use the same tools beyond the equilibrium shortcut

These applications are independent study branches. Each asks a precise question that occupancy weights alone cannot settle.

23Prove a response prohibition, then test a driven counterexample

An equilibrium restriction must specify the allowed state ledger. Consider one transcription-factor binding site and one polymerase binding site. The four states are empty, TF-bound, polymerase-bound and jointly bound. Hold all chemical activities except free TF fixed, and let $\bar R=R/K_{\mathrm{ref}}$ be dimensionless TF concentration. Equilibrium weights have the form $1,a\bar R,b,c\bar R$ , where $a,b,c>0$ collect fixed affinities and interactions.

At most one TF molecule in a state makes both sums affine in its concentration. For any fixed nonnegative state activities, the weighted mean has the form

y_{\mathrm{eq}}(\bar R)=\frac{A_0+A_1\bar R}{B_0+B_1\bar R}, \qquad \frac{dy_{\mathrm{eq}}}{d\bar R} =\frac{A_1B_0-A_0B_1}{(B_0+B_1\bar R)^2}.

The positive denominator cannot change the sign of the constant numerator. An interior peak is impossible under this ledger and these fixed readout assumptions. The argument covers more states if none contains more than one molecule of the varied input. It does not cover changing output activities, hidden varying inputs or states with higher input stoichiometry.

A driven square supplies a constructive alternative on the same four states. Label empty as 0, TF-bound as 1, polymerase-bound as 2 and jointly bound as 3. Use the following state-transition hazards, all in $\mathrm{s}^{-1}$ :

Clockwise edge	Forward hazard	Reverse hazard
0 to 1	$\bar R$	1
1 to 3	1	10
3 to 2	1	$0.01\bar R$
2 to 0	0.1	0.01

For positive input, the forward/reverse product ratio is 100. Thus $\mathcal A=\ln100$ . Choose polymerase occupancy as output, $y=p_2^*+p_3^*$ . Solving the four stationary equations gives the exact illustrative response

y(\bar R)=\frac{600+5501\bar R+50\bar R^2}{6600+11506\bar R+550\bar R^2} =\frac1{11}+\frac{4455\bar R}{6600+11506\bar R+550\bar R^2}.

The peak can be proved rather than inferred from a plot. Differentiating the second term produces a numerator proportional to $6600-550\bar R^2$ . Hence the unique interior maximum is at $\bar R=\sqrt{12}$ , where $y\approx0.3818$ . Both endpoint limits are 1/11. Reducing only the 0-to-1 hazard to $0.01\bar R$ restores detailed balance and gives constant polymerase occupancy 1/11.

Scroll horizontally for the full-size figure.

Figure 16. A bounded equilibrium theorem and an explicit kinetic witness. The four-state graph and eight hazards define the comparison completely. The driven model has an interior maximum. The balanced comparison is flat. This is an illustrative calculation, not a fit to a measured promoter.

Checkpoint: does a measured peak prove fuel consumption?

It rejects the specified one-input equilibrium ledger with fixed activities. It does not reject every equilibrium model. Before concluding that driving is required in the experiment, test additional binding states, concentration-dependent activities, input depletion and hidden changes in other regulators. A complete energetic cost also requires a physical channel and reservoir assignment.

Prove what a specified equilibrium model cannot do and compare it with a reproducible driven model without overgeneralizing the conclusion.

24Calculate error, yield and cost in two-stage proofreading

A second rejection stage can amplify an existing discrimination. Start after an initial selection has produced an incoming wrong/right ratio $f_0<1$ . Each arrival enters a paid intermediate. From a correct intermediate, commitment has hazard $W$ and discard has hazard $f_0\ell$ . From a wrong intermediate, commitment still has hazard $W$ , but discard has hazard $\ell$ . These are exit hazards in a limiting attempt model. They do not resolve a complete thermodynamic mechanism.

Competing clocks give the success probability of each attempt. Integrating the commitment/event density gives $s_{\mathrm{right}}=W/(W+f_0\ell)$ and $s_{\mathrm{wrong}}=W/(W+\ell)$ . Multiply each incoming flux by its success probability. The wrong/right product ratio is

f=\frac{J_{\mathrm{wrong}}}{J_{\mathrm{right}}} =f_0\frac{W+f_0\ell}{W+\ell}.

The squared-error limit requires commitment to be slower than both discard rates. If $W\ll f_0\ell$ , then $f\approx f_0^2$ . The weaker condition $W\ll\ell$ is not enough when $f_0$ is small. At fast commitment, $W\gg\ell$ , nearly every intermediate becomes product and $f\to f_0$ . The product fraction that is wrong is $f/(1+f)$ , not exactly $f$ .

Yield and throughput must be named separately. The probabilities of a correct and wrong arrival are $1/(1+f_0)$ and $f_0/(1+f_0)$ . Yield per intermediate arrival is therefore

Y=\frac{1}{1+f_0} \left[\frac{W}{W+f_0\ell}+f_0\frac{W}{W+\ell}\right]. \qquad J_{\mathrm{product}}=J_{\mathrm{attempt}}Y.

The last equality uses a separately specified attempt flux $J_{\mathrm{attempt}}$ . In an enzyme model, that flux may itself depend on intermediate residence times and enzyme availability. Holding it fixed is a declared comparison, not a theorem about speed.

A fuel-per-product estimate needs an explicit fuel rule. If every intermediate arrival consumes one fuel conversion, including successful attempts, then fuel conversions per product equal $1/Y$ . Counting only discarded attempts would give $(1-Y)/Y$ and omit the successful attempts' cost. Converting either count into dissipated energy requires the fuel chemical potential and a thermodynamically consistent mechanism.

Scroll horizontally for the full-size figure.

Figure 17. Fidelity, yield and fuel conversions trade off in a specified limit. The incoming wrong/right ratio is 0.01 and the wrong-discard hazard is 1/s. The sweep changes commitment

W

. Separate panels show product error ratio, yield per arrival and fuel conversions per product under the one-fuel-per-arrival rule. Neither yield nor fuel count is silently labeled speed or affinity.

Driving is needed to realize the reset mechanism without equilibrium cancellation. In Hopfield's construction, a fuel-coupled step supplies a high-energy intermediate and suppresses the reverse routes that would undo the extra discrimination. The limiting calculation above isolates the second selection, after the first has supplied $f_0$ . It does not establish the same error improvement for arbitrary rates or an equilibrium network. Restoring reverse channels is necessary to calculate finite affinities and dissipation.⁷

Checkpoint: compute one operating point.

With $f_0=0.01$ , $\ell=1\,\mathrm{s}^{-1}$ and $W=0.001\,\mathrm{s}^{-1}$ , the error ratio is approximately $1.099\times10^{-4}$ . Yield is about 0.0900 and fuel conversions per product are about 11.11. These numbers do not determine product throughput until the attempt flux has been specified.

Calculate proofreading discrimination and its costs while keeping the limiting regime, observable and fuel-accounting assumptions explicit.

25Derive sensing precision from temporal correlations

Time averaging gains information only as the sensor forgets earlier samples. Consider a two-state receptor at stationary occupancy $p=a/(a+b)$ . Let $B(t)$ equal one when bound and zero when free. Write $\kappa=a+b$ . The ON/OFF solution gives $\Pr[B(t)=1\mid B(0)=1]=p+(1-p)e^{-\kappa t}$ . Multiplying by the initial bound probability gives the covariance

\operatorname{Cov}[B(0),B(t)] =p(1-p)e^{-\kappa|t|}.

Integrate this covariance over the observation window. Define the time-averaged occupancy $\bar B_T=T^{-1}\int_0^T B(t)\,dt$ . Stationarity makes its variance a double integral depending only on the time difference. Folding the square integration region about its diagonal gives

\begin{aligned} \operatorname{Var}\bar B_T &=\frac{2p(1-p)}{T^2}\int_0^T(T-u)e^{-\kappa u}\,du\\ &=2p(1-p)\frac{\kappa T-1+e^{-\kappa T}}{(\kappa T)^2}. \end{aligned}

At short times this tends to the Bernoulli variance $p(1-p)$ . At long times it is approximately $2p(1-p)/(\kappa T)$ . Comparing with the variance of a mean of independent Bernoulli samples gives an effective sample count $N_{\mathrm{eff}}\approx\kappa T/2$ . The factor of two belongs to this continuous-time occupancy estimator.

Convert occupancy error into concentration error through the response slope. For $a=k_+c$ , $b=k_-$ and known $K_d=k_-/k_+$ , concentration is $c=K_dp/(1-p)$ . For sufficiently small estimation error, linear propagation gives

\frac{\operatorname{Var}\widehat c}{c^2} \simeq\frac{\operatorname{Var}\bar B_T}{p^2(1-p)^2} \simeq\frac{2}{\kappa T p(1-p)} \qquad(\kappa T\gg1).

This is the performance of occupancy averaging in the stated well-mixed Markov model. It is not a bound on every estimator of a fully resolved binding trajectory. Near complete saturation, a small occupancy error corresponds to a large concentration error.

Scroll horizontally for the full-size figure.

Figure 18. A long record contains correlated samples. Left: stationary receptor autocovariance for hazards 2/s and 1/s. Right: the exact variance of average occupancy, divided by its instantaneous variance, compared with the long-time approximation. At short times the approximation fails. At long times the normalized variance falls as

2/(\kappa T)

Diffusion introduces another clock and another measurement model. For a spherical sampling region of radius $a_{\mathrm{sens}}$ , renewal takes roughly $a_{\mathrm{sens}}^2/D$ . The region contains of order $\rho_La_{\mathrm{sens}}^3$ independent dilute ligands. Combining population fluctuations with repeated renewal gives fractional error of order $(Da_{\mathrm{sens}}\rho_LT)^{-1/2}$ , where $\rho_L$ is molecular number density, not molar concentration.

The precise benchmark depends on the detector

Berg and Purcell's transparent-volume instrument observes the instantaneous number of molecules inside a sphere and averages that count. In the long-time limit its rms concentration error is

\frac{\sigma_c}{c}= \sqrt{\frac{3}{5\pi Da_{\mathrm{sens}}\rho_LT}}.

This is their equation (41), not a universal numerical prefactor for all sensors. An ideal perfectly absorbing sphere with Poisson capture events instead has mean arrival count $4\pi Da_{\mathrm{sens}}\rho_LT$ and fractional counting variance $1/(4\pi Da_{\mathrm{sens}}\rho_LT)$ . Absorbing a ligand and repeatedly observing a freely diffusing ligand are different experiments.⁸

Changing the correlation time can improve precision without proving a need for fuel. Speeding both passive binding hazards at fixed ratio already increases the effective sample count in the simple Markov model. Active resetting may do something similar, but may also change occupancy, rebinding, receptor availability and the estimator. A claim that driving improves sensing must specify that causal comparison and its physical cost.

Checkpoint: check the concentration units.

For a molar concentration $c_{\mathrm M}$ , number density in molecules per cubic meter is $\rho_L=1000N_Ac_{\mathrm M}$ . In micrometer units, 1 nM is approximately 0.6 molecules/μm³. Inserting molarity directly into a formula with $D$ in μm²/s and radius in μm would be dimensionally wrong.

Derive the precision of a stated time-average estimator and distinguish kinetic sampling limits from detector-specific diffusion benchmarks.

26Work a new model independently and test its assumptions

A finished model should expose what an experiment could disprove. State the subsystem, state variables, resolved channels, controlled inputs, output and time window. Assign either kinetic hazards or compatible physical weights. Close finite pools. Then choose an occupancy, transient or channel measurement that interrogates the approximation actually being used.

Problem A: a competitive sensor.

A receptor binds either A or B, never both. The free concentrations are maintained. Derive its A-bound probability and explain what a dose-response fit would infer if B were ignored.

Solution. The three weights are $1,A/K_A,B/K_B$ . Hence $p_A=(A/K_A)/(1+A/K_A+B/K_B)$ . At fixed B, half of the asymptotic A-bound fraction occurs at $A=K_A(1+B/K_B)$ . A single-input fit sees an apparent affinity shifted by the competing ligand. It does not measure the isolated microscopic $K_A$ .

Problem B: an equilibrium constant and a first-passage time.

A passive two-state receptor has bound/free weight ratio 4. Does this determine its mean binding time from free?

Solution. No. Hazards $a=4k_0,b=k_0$ for any positive rate scale $k_0$ give bound fraction 4/5. The mean first binding time is $1/(4k_0)$ . Measure a holding time or relaxation transient to identify the speed.

Problem C: repair a cycle.

A three-state ring has forward hazards (2, 3, 4)/s and reverse hazards (1, 6, 1)/s. Change only the third reverse hazard to make the graph detailed balanced.

Solution. The forward product is 24/s³. The reverse product is $6u$ /s², so choose $u=4\,\mathrm{s}^{-1}$ . Starting from weight 1 at state 1 gives state-2 weight 2 and state-3 weight 1. Normalization yields $p^{\mathrm{eq}}=(1,2,1)^\top/4$ . Check all three pairs, not just the product.

Problem D: a moving rare state.

A fast-state approximation has an absolute probability-error bound of 0.002. The activity is produced only by a state whose equilibrium probability is 0.001. Is the reduced activity necessarily accurate to 1%?

Solution. No. The stated bound exceeds the predicted rare-state probability. Activity-weighted error or a relative bound for that particular state is needed. A good approximation for the dominant populations can be poor for a rare functional state.

Problem E: infer equilibrium from a histogram.

All three enzyme states have stationary probability 1/3. Does that show the enzyme performs no net conversion?

Solution. No. The two rings in section 7 have this same histogram and different currents. Count directional transitions or independently measure substrate-to-product throughput. If only a two-state projection is observed, also examine whether physical channels have been hidden.

The method ends with separate claims, not one global equilibrium label. State sufficiency justifies the dynamical variables. Fixed conditions justify time-homogeneous hazards. Cycle compatibility justifies equilibrium weights on resolved channels. Fast attraction and slowly moving inputs justify a quasi-equilibrium reconstruction. Conservation connects free species to totals. Activities connect states to function. Each claim has its own calculation and possible failure.

Construct, solve and challenge a molecular-state model on a biological system not worked in the lecture.

Optional workshop: reanalyze the original saturation experiments

The goal is to distinguish a measured response from the mechanism used to explain it. We will revisit three original examples with the methods of Lectures 3–7. Michaelis–Menten concerns catalytic speed. Hill and the haemoglobin example of MWC concern oxygen occupancy. Each calculation must identify its states, its input and its readout before comparing with the data.

1. Invertase: one hyperbola, two reductions

Begin with a resolved binding step and a lumped conversion. Let $S$ be substrate, $E$ free enzyme, $C_{ES}$ the complex and $P$ product. Use elementary reversible binding but keep catalysis composite:

E+S\xrightleftharpoons[k_-]{k_+}C_{ES}, \qquad C_{ES}\rightsquigarrow E+P, \qquad v_{\mathrm{cat}}=k_{\mathrm{cat}}C_{ES}.

The first-order conversion law is an effective assumption, not evidence that catalysis has only one microscopic step.

The constants and flux carry different units. The association constant $k_+$ has concentration⁻¹ time⁻¹ units. The constants $k_-$ and $k_{\mathrm{cat}}$ have inverse-time units. Here $k_{\mathrm{cat}}$ is a conversion frequency per bound enzyme, while $v_{\mathrm{cat}}$ is a concentration flux.

The rate equation keeps both binding and conversion. For $x=(S,E,C_{ES},P)^\top$ , the Lecture 3 form is

\dot x=\Gamma v(x),\qquad \Gamma=\begin{pmatrix}-1&1&0\\-1&1&1\\1&-1&-1\\0&0&1\end{pmatrix}, \qquad v=\begin{pmatrix}k_+ES\\k_-C_{ES}\\k_{\mathrm{cat}}C_{ES}\end{pmatrix}.

Enzyme conservation gives $q_E=E+C_{ES}$ . The substrate-bearing total is $q_S=S+C_{ES}$ , so $\dot q_S=-v_{\mathrm{cat}}$ and $\dot P=v_{\mathrm{cat}}$ exactly in this model. These are Lecture 6's slow variables. Reconstructing $C_{ES}$ from them is the separate approximation.

Rapid binding equilibrium balances only the binding pair. If binding relaxes before appreciable conversion and free substrate changes slowly, set $k_+ES\simeq k_-C_{ES}$ . Define $K_d=k_-/k_+$ and use the enzyme total:

C_{ES}^{\mathrm{eq}}=q_E\frac{S}{K_d+S}, \qquad v_{\mathrm{cat}}^{\mathrm{red}}=k_{\mathrm{cat}}q_E\frac{S}{K_d+S}.

The same result follows from state weights $w_E=1$ , $w_{C_{ES}}=S/K_d$ and activities $0,k_{\mathrm{cat}}$ . The dimensionless binding free-energy difference is $-\ln(S/K_d)$ , including the substrate reservoir. This route predicts an allocation. It does not supply a relaxation time.

Quasi-steady complex balance includes the conversion exit. The exact complex equation is $\dot C_{ES}=k_+S(q_E-C_{ES})-(k_-+k_{\mathrm{cat}})C_{ES}$ . Neglecting its derivative at leading order on a justified slow timescale gives

C_{ES}^{\mathrm{qss}}=q_E\frac{S}{K_M+S}, \qquad K_M=\frac{k_-+k_{\mathrm{cat}}}{k_+}, \qquad v_{\mathrm{cat}}^{\mathrm{red}}=V_{\max}\frac{S}{K_M+S}, \quad V_{\max}=k_{\mathrm{cat}}q_E.

This is the Briggs–Haldane derivation. It need not assume $k_{\mathrm{cat}}\ll k_-$ . That extra condition makes $K_M\simeq K_d$ . As Lecture 6 emphasized, a complex following the slow manifold still changes as substrate is depleted. Its physical derivative is not identically zero.^BH

2. The same model as a stochastic enzyme, with a visible clock

A fixed substrate bath turns one enzyme into a two-state Markov model. Its states are free and substrate-bound. Association has hazard $k_+S$ . Unbinding and conversion are distinct exit channels from the complex. With probability column $p=(p_E,p_{C_{ES}})^\top$ ,

\dot p=Qp,\qquad Q=\begin{pmatrix}-k_+S&k_-+k_{\mathrm{cat}}\\ k_+S&-(k_-+k_{\mathrm{cat}})\end{pmatrix}. \qquad p_{C_{ES}}^*=\frac{S}{K_M+S}.

The stationary mean conversion frequency is a state activity average. It equals $k_{\mathrm{cat}}p_{C_{ES}}^*$ per enzyme. The instantaneous bound probability solves

p_{C_{ES}}(t)=p_{C_{ES}}^* +[p_{C_{ES}}(0)-p_{C_{ES}}^*]e^{-t/t_{\mathrm{relax}}}, \qquad t_{\mathrm{relax}}=\frac1{k_+S+k_-+k_{\mathrm{cat}}}.

This is also the linear stability calculation: the nonzero eigenvalue is $-1/t_{\mathrm{relax}}$ . Multiplying all three rates by the same factor leaves the occupied fraction unchanged while changing both the clock and the conversion frequency. A static binding fraction cannot determine either.

Even the same catalytic rate curve can hide different binding kinetics.

As an illustrative identifiability test, hold $k_{\mathrm{cat}}=1\,\mathrm{s}^{-1}$ and $q_E$ fixed. Compare $(k_+,k_-)=(10,99)$ with $(0.2,1)$ , using $(\mu\mathrm M)^{-1}\mathrm{s}^{-1}$ and $\mathrm{s}^{-1}$ , respectively. Both have $K_M=10\,\mu\mathrm M$ and the same $V_{\max}$ . Their binding constants are 9.9 and 5 μM. At $S=10\,\mu\mathrm M$ , their relaxation times are 0.005 and 0.25 s. A time-resolved experiment separates what the rate curves cannot. These numbers are not estimates for invertase.

Two-state balance is not a thermodynamic certificate

The aggregate free/bound chain balances its only edge at stationarity. Yet the catalytic exit releases product and the binding entry consumes substrate. Erasing channel identity and reservoir changes erases the chemical turnover. Equilibrium weights apply physically to the binding-only subsystem only when that reduction is justified. For a thermodynamic current calculation, restore product states and reverse chemistry, as in the core's enzyme cycle.

3. What the 1913 data test beyond an initial-rate curve

Michaelis and Menten measured sucrose hydrolysis through optical rotation. They controlled hydrogen-ion concentration and temperature, stopped samples, and waited for product mutarotation before reading the signal. Initial slopes reduced the immediate complications of product inhibition. Their table gives seven initial slopes at different sucrose concentrations.^MM

Our initial-rate fit is a reanalysis of the printed table. Unweighted least squares gives $K^{\mathrm{fit}}=17.89\,\mathrm{mM}$ . Their graphical analysis reported about 16.7 mM. The two procedures are not the same estimate. Both data and line below use our fitted amplitude for normalization.

Scroll horizontally to read both panels at full size.

Historical Figure H1. Fit a response, then challenge the mechanism over time. Left: the seven initial-rate values in the translation's p. 5 table, normalized by our fitted amplitude, 3.911 in the reported ordinate units. No conversion of those graphical units to molar flux is inferred. Right: series IV of the p. 33 progress table. Both curves use the original integrated-equation constants, not a fit to these seven points. The table's final 1497-minute point is outside the 0–180-minute plot.

Product inhibition becomes another equilibrium state-allocation problem. For actual invertase, one sucrose yields one fructose $F$ and one glucose $G$ . Add nonproductive enzyme states bound to each product. The four weights are $1,S/K_S,F/K_F,G/K_G$ . Only sucrose-bound enzyme contributes to hydrolysis. Assuming rapid competitive binding and negligible enzyme-bound solute compared with the free pools gives

v=\frac{V_{\max}(S/K_S)} {1+S/K_S+F/K_F+G/K_G}.

A single extent variable closes this dilute-enzyme progress model. Starting without products, let $P$ denote the concentration of sucrose hydrolyzed, not the sum of both products. Then $F\simeq G\simeq P$ and $S\simeq S_0-P$ . Bound-pool corrections have been neglected in these free-concentration relations. Define $\beta=1/K_F+1/K_G$ . Separating variables in $\dot P=v$ yields

V_{\max}t=(1-K_S\beta)P +K_S(1+\beta S_0)\ln\frac{S_0}{S_0-P}.

The argument of the logarithm is dimensionless. The right-hand side has concentration units. The original reciprocal estimates $K_S=1/60\,\mathrm M$ , $K_F=1/17\,\mathrm M$ and $K_G=1/11\,\mathrm M$ give the rounded coefficients in the paper's integrated equation:

\underbrace{\frac{V_{\max}}{K_S}}_{0.0454\,\mathrm{min}^{-1}}t =(1+28\,\mathrm M^{-1}S_0)\ln\frac{S_0}{S_0-P} +32\,\mathrm M^{-1}P.

The progress comparison reuses parameters rather than fitting each time point. At $S_0=41.6\,\mathrm{mM}$ and 90 minutes, this equation predicts $P/S_0=0.7585$ . Setting $\beta=0$ while retaining the same $V_{\max},K_S$ predicts 0.8572. The historical observation at 90.75 minutes is 0.747. Product inhibition matters on this timescale even when the initial-rate curve alone looks simple.

The enzyme amount was not known in molar units. The historical data can constrain an amplitude proportional to enzyme amount, but they cannot separately determine $k_{\mathrm{cat}}$ and $q_E$ . Johnson and Goody's translation and reanalysis explain both this limitation and the original full-progress analysis. The constants were not measurements of every microscopic rate.^MM

4. Hill's oxygen data: a slope is not a molecularity

The measured quantity is now a bound-site fraction. Write $\theta$ for fractional oxygen saturation and $p_{\mathrm{O_2}}$ for oxygen pressure. For one independent binding site, a pressure-form dissociation constant $K_p$ gives $\theta=p_{\mathrm{O_2}}/(K_p+p_{\mathrm{O_2}})$ . At fixed temperature, converting to dissolved oxygen concentration uses a solubility coefficient. We retain the historical pressure units rather than silently treating pressure as concentration.

Hill's 1910 paper proposed two different descriptions of one data series. For haemoglobin in Na₂HPO₄ solution, it gives $n=1.670$ and $K=0.0453\,(\mathrm{mmHg})^{-1.670}$ in $\theta=Kp_{\mathrm{O_2}}^n/(1+Kp_{\mathrm{O_2}}^n)$ . The product in that formula is dimensionless. Its midpoint is $p_{50}=K^{-1/n}=6.379\,\mathrm{mmHg}$ , so a cleaner modern form is

\theta_H(p_{\mathrm{O_2}})= \frac{(p_{\mathrm{O_2}}/p_{50})^n}{1+(p_{\mathrm{O_2}}/p_{50})^n}.

Source excerpt H2a. Hill's observed and calculated rows are different evidence. Cropped from Hill (1910), printed p. vi, series I. Hill credits the observations to Barcroft and Camis. Our points use the observed row. The calculated row is not experimental data. The scan is unretouched.

His aggregation description used a different state hypothesis. Let $u=p_{\mathrm{O_2}}/(1\,\mathrm{mmHg})$ denote the numerical pressure in those units. The published mixture is

\theta_{\mathrm{mix}} =0.848\frac{0.0236u^2}{1+0.0236u^2} +0.152\frac{0.125u}{1+0.125u}.

The historical hypothesis allocated binding capacity to monomer and dimer aggregates. This is not the modern structure of haemoglobin. Yet its curve nearly overlaps the simpler Hill curve on the reported data. The observation constrains a response more directly than it selects a state hypothesis.^Hill

Scroll horizontally to read the curve comparison.

Historical Figure H2. Two published explanations, one measured table. Solid and dashed curves reconstruct Hill's parameters without refitting. The dotted one-site curve has the same midpoint as the Hill curve. The source reports no replicate uncertainty here, so no confidence band or formal model-selection score is asserted.

The local Hill slope can be defined without assuming a Hill model. For any differentiable binding curve with $0<\theta<1$ , define

n_H(p)=\frac{d\ln[\theta/(1-\theta)]}{d\ln p} =\frac{p\,\theta'(p)}{\theta(1-\theta)}.

For the exact Hill law this equals $n$ everywhere. For a general equilibrium state sum it varies with input. A fitted constant exponent, a midpoint derivative and the number of binding sites are therefore distinct quantities. None is the molecularity of an elementary reaction merely because it appears as an exponent in an effective formula.

Sequential states show how a Hill-like law can emerge. For a two-site molecule, let the aggregated weights for zero, one and two bound ligands be $1,a p,b p^2$ . Multiplicity factors are already included in $a,b$ . The occupied-site fraction is

\theta=\frac{a p+2b p^2}{2(1+a p+b p^2)}.

If the singly occupied weight is negligible over the transition range, then $\theta\simeq b p^2/(1+b p^2)$ . A Hill exponent of two emerges without requiring simultaneous two-ligand binding. Hill's 1913 discussion explicitly considers this intermediate and its suppression. The approximation concerns equilibrium populations, not the absence of a sequential route.^{Hill 1913}

5. MWC: retain partial occupancy and add a shared conformation

The original haemoglobin example uses four binding sites and two conformations. Keep the historical names $R,T$ . They do not automatically mean catalytically active and inactive. Sites are independent and equivalent conditional on a conformation. The concerted model excludes mixed-conformation oligomers but permits every occupancy from zero to four. Real haemoglobin has nonidentical alpha and beta chains, so equivalence is an idealization of this model, not a literal description of four identical chains.^MWC

Enumerate occupancy before normalizing. Set $\alpha=p_{\mathrm{O_2}}/K_R$ , $c=K_R/K_T$ and $L=w_{T,0}/w_{R,0}$ . Here $c$ is the original dimensionless affinity ratio, not the ligand concentration used elsewhere in this chapter. For $j=0,\ldots,4$ ,

w_{R,j}=\binom4j\alpha^j,\qquad w_{T,j}=L\binom4j(c\alpha)^j,\qquad Z=(1+\alpha)^4+L(1+c\alpha)^4.

The readout determines which weighted sum is needed. Averaging $j/4$ gives occupied-site fraction. Averaging the indicator of conformation $R$ gives conformation fraction:

\theta= \frac{\alpha(1+\alpha)^3+Lc\alpha(1+c\alpha)^3}{Z}, \qquad p_R=\frac{(1+\alpha)^4}{Z}.

MWC's Figure 3 uses $L=9054$ and $c=0.014$ for Lyster's horse-haemoglobin oxygen data. It reports 4.6% haemoglobin, 0.6 M phosphate, pH 7 and 19°C. These are not the conditions of Hill's earlier table. The two data sets must not be treated as measurements of one unchanged curve.

Source excerpt H3a. Data and fitted line in the original MWC example. Cropped from Monod, Wyman and Changeux (1965), Figure 3, p. 94. The paper credits the points to R. W. J. Lyster, then unpublished. The solid line is the authors' model fit. The excerpt preserves that distinction and has not been retouched or digitized.

The printed figure does not give the absolute affinity scale. Its two printed constants determine the dimensionless shape, not $K_R$ . Solve $\theta(\alpha_{50})=1/2$ to obtain $\alpha_{50}=9.8968$ . We can then compare shapes on $p_{\mathrm{O_2}}/p_{50}=\alpha/\alpha_{50}$ without inventing the missing pressure calibration.

Scroll horizontally to read the distinct readouts.

Historical Figure H3. Reconstruct the model, not unreported measurements. All lines are calculated. The MWC occupancy and conformation curves share the same occupancy midpoint. The equal-affinity hyperbola is separately normalized by its own midpoint. At half occupancy, the MWC conformation fraction is 0.4810 and the occupancy Hill slope is 2.853.

The one-site-shaped hyperbola is a limit of the four-site model. If $c=1$ , then $Z=(1+L)(1+\alpha)^4$ and $\theta=\alpha/(1+\alpha)$ . Also $p_R=1/(1+L)$ : ligand no longer selects a conformation. Four independent sites have the same fractional occupancy curve as one independent site. Taking $L\to0$ likewise leaves only the $R$ branch and its hyperbolic occupancy.

A nearly all-or-none occupancy distribution gives a different limit. If the empty $T$ state and fully occupied $R$ state dominate over a transition range, then $Z\simeq L+\alpha^4$ and $\theta\simeq\alpha^4/(L+\alpha^4)$ . This is a Hill-like limiting regime, not the general MWC formula. The original Figure 3 parameter pair gives a midpoint slope of 2.853 because intermediate states still matter.

An effector tests the conformation hypothesis across curves. In a simple extension with $m$ independent effector sites, an effector concentration $a$ multiplies the conformation weights by $(1+a/K_{R,a})^m$ and $(1+a/K_{T,a})^m$ . Its effect on oxygen binding is therefore captured by

L_{\mathrm{eff}}(a)=L \left(\frac{1+a/K_{T,a}}{1+a/K_{R,a}}\right)^m.

An effector favoring $R$ lowers $L_{\mathrm{eff}}$ . The prediction is a related family of curves, not an arbitrary new Hill exponent for every condition. MWC's original paper develops such heterotropic regulation and cautions that a successful saturation fit alone does not establish its assumptions.

6. Choose the method by the question, not by the curve's name

Course method	What it calculates here	Information still required
Lecture 3: $\dot x=\Gamma v(x)$	Binding, substrate depletion and product formation together	Rate laws, initial state and kinetic constants
Lecture 4: linear stability	Relaxation of a perturbed fixed-bath enzyme occupancy	Absolute hazards, not only their ratios
Lecture 5: Markov dynamics	Single-enzyme occupancy and random event timing	Channel identity for conversion counts and thermodynamics
Lecture 6: time-scale separation	QSSA or rapid-binding reconstruction with evolving totals	The relevant scale hierarchy and conservation constraints
Lecture 7: states and weights	Binding fractions and state-weighted activities	Compatible passive fast states, free input and a readout rule

Try a discriminating experiment.

For invertase, measure binding relaxation or challenge the product-inhibition prediction. For haemoglobin, measure partial occupancy, conformation or responses to an independently varied effector. State which alternative each observation would reject. A better fit to the same curve is not automatically the stronger experiment.

Use the historical data to test a model's declared prediction. Do not infer a binding constant, a site count, a kinetic clock or an energy budget from the name of a fitted saturation law.

Binding playground: change affinity, speed and activity separately

This passive two-state example starts with a free gene. Free repressor is a fixed reservoir. Changing the dissociation hazard at fixed affinity also changes association, since $k_+=k_-/K_d$ . The plot spans five relaxation times, so its physical time axis changes when the clock changes. Activity is normalized to the free-state initiation rate.

One state distribution, two observables

Move any control. The occupancy and activity curves update immediately. Numbers are illustrative.

Free repressor

R

(nM): Dissociation constant

K_d

(nM): Dissociation hazard

k_-

(per second): Bound-state relative activity

r_b/r_0

Scroll horizontally for the full-size plot.

The horizontal range is five relaxation times, with ticks in seconds. Matching the ligand/affinity ratio preserves the stationary fraction. Changing the clock still changes the physical time scale.

References

These sources support the construction and its assumptions. Figures identify illustrative parameter choices and reconstructions of published models separately. Historical development and priority are separate questions from reproducing a calculation.

D. T. Gillespie (1977). Exact stochastic simulation of coupled chemical reactions. Journal of Physical Chemistry 81, 2340–2361. Article. The event/time density and reaction-specific propensities underlie sections 2 and 5.
K. M. Nam, R. Martinez-Corral and J. Gunawardena (2022). The linear framework: using graph theory to reveal the algebra and thermodynamics of biomolecular systems. Interface Focus 12, 20220013. Article. Sections 2 and 4 develop column graph dynamics, tree formulas, reservoirs and equilibrium cycle compatibility.
C. Maes (2021). Local detailed balance. SciPost Physics Lecture Notes 32. Article. The physical interpretation of rate ratios requires a reservoir model. Section V.A demonstrates the danger of unresolved channels.
M. Feinberg (2019). Foundations of Chemical Reaction Network Theory. Springer. Book. Chapter 14 distinguishes graph-cycle conditions from the additional stoichiometric conditions needed for detailed balance.
L. Bintu et al. (2005). Transcriptional regulation by the numbers: models. Current Opinion in Genetics & Development 15, 116–124. Article. Molecular multiplicities produce occupancy weights. Equilibrium and occupancy-to-expression are separate assumptions.
S. Marzen, H. G. Garcia and R. Phillips (2013). Statistical mechanics of Monod–Wyman–Changeux models. Journal of Molecular Biology 425, 1433–1460. Article. The state-counting treatment distinguishes activity, binding and effective sharpness.
J. J. Hopfield (1974). Kinetic proofreading: a new mechanism for reducing errors in biosynthetic processes requiring high specificity. PNAS 71, 4135–4139. Article. The driven intermediate and rate inequalities are essential to the additional discrimination.
H. C. Berg and E. M. Purcell (1977). Physics of chemoreception. Biophysical Journal 20, 193–219. Article. Equations (31)–(41) concern transparent-volume count averaging. Detector geometry and correlation time are part of the prediction.
A. Flamholz, E. Noor, A. Bar-Even, W. Liebermeister and R. Milo (2013). Glycolytic strategy as a tradeoff between energy yield and protein cost. PNAS 110, 10039–10044. Article. Pages 10040–10041 and Figure 3 report the constrained energy benchmark. The isomerase teaching calculation uses that energy scale, not the paper's fitted enzyme kinetics.
U. Seifert (2012). Stochastic thermodynamics, fluctuation theorems and molecular machines. Reports on Progress in Physics 75, 126001. Article. Sections 6.1–6.3 distinguish stationary currents, ensemble entropy production and physically meaningful energy assignments. The Lyapunov calculation in section 11 is derived explicitly here.
M. Razo-Mejia et al. (2018). Tuning transcriptional regulation through signaling: a predictive theory of allosteric induction. Cell Systems 6, 456–469.e10. Article. In-press full text pp. 3–6, equations (1)–(5), Figure 2 and the parameter-inference section. The reconstructed curves use central estimates, not experimental measurements.
L. Michaelis and M. L. Menten. “Die Kinetik der Invertinwirkung.” Biochemische Zeitschrift 49 (1913), 333–369. Inspected in K. A. Johnson and R. S. Goody's full English translation, supplementary to Biochemistry 50 (2011), 8264–8269. Translation pp. 1–6, 10–16, 27, 31–33. The accompanying article explains the historical units and progress analysis. Source
G. E. Briggs and J. B. S. Haldane. “A note on the kinetics of enzyme action.” Biochemical Journal 19 (1925), 338–339. Both pages: complex balance, the kinetic constant and the inability to identify unbinding versus conversion from the rate curve. Source
A. V. Hill. “The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves.” Journal of Physiology 40, Proceedings (1910), iv–vii. Original Na₂HPO₄ observations and two published model descriptions, pp. v–vi. Observations credited there to Barcroft and Camis. Source
A. V. Hill. “The combinations of haemoglobin with oxygen and with carbon monoxide. I.” Biochemical Journal 7 (1913), 471–480. Pages 471–474 explicitly consider a singly occupied intermediate and its small-population limit. Source
J. Monod, J. Wyman and J.-P. Changeux. “On the nature of allosteric transitions: a plausible model.” Journal of Molecular Biology 12 (1965), 88–118. State sums pp. 90–93, haemoglobin Figure 3 and qualifications p. 94, effectors pp. 95–103 and haemoglobin structural discussion pp. 110–112. Source