Translating Reaction Steps Into ODEs

From Reaction Steps to Concentration ODEs

A reaction list describes which transformations are included in a model; a set of ordinary differential equations (ODEs) specifies how the concentrations change over time. The translation is systematic: for each reaction, determine the change in each species per reaction event, multiply that change by the reaction rate, and add the contributions from all reactions. This accounting makes it possible to construct a model without guessing the signs or coefficients in each equation.

The method assumes a well-mixed, constant-volume system in which concentrations are represented by continuous, differentiable state variables. Each listed reaction is treated as a directional reaction channel with a specified rate law. The rate \(v_j\) is the number of reaction events per unit volume per unit time, so it has concentration-per-time units. The rate law and its parameters must be appropriate for the conditions being modeled. In the examples, we use the mass-action rate laws for elementary steps introduced in the earlier material on reaction kinetics.

Definition: For reaction \(j\), the stoichiometric coefficient of species \(i\) is \(\nu_{ij}\), defined as the product coefficient minus the reactant coefficient. Thus \(\nu_{ij}\) is negative for a consumed species, positive for a produced species, and zero if the species does not participate. If the reaction rate is \(v_j\), its contribution to the rate of change of species \(i\)'s concentration is \(\nu_{ij}v_j\).

If the model has \(m\) species and \(r\) reactions, let \(C_i\) denote the concentration of species \(i\), and let \(v_j\) denote the rate of reaction \(j\). The ODE for each species is the sum of its contributions from all reactions:

\frac{dC_i}{dt}=\sum_{j=1}^{r}\nu_{ij}v_j,\qquad i=1,\ldots,m.

In matrix notation, let \(N\) be the stoichiometric matrix whose entry in row \(i\), column \(j\) is \(\nu_{ij}\), and let \(v\) be the vector of reaction rates. Then \(dC/dt=Nv\). The matrix form is compact, but the underlying operation is still reaction-by-reaction addition.

Why the Sum Gives the ODE

Consider a short time interval of length \(\Delta t\). If the rate of reaction \(j\) is approximately \(v_j\) over that interval, then approximately \(v_j\Delta t\) reaction events occur per unit volume. Each event changes the amount per volume of species \(i\) by \(\nu_{ij}\). The reaction’s approximate concentration change is therefore \(\nu_{ij}v_j\Delta t\). With several reactions, their changes add:

\Delta C_i\approx\sum_{j=1}^{r}\nu_{ij}v_j\Delta t.

Dividing by \(\Delta t\) and taking the limit as \(\Delta t\) approaches zero gives the stated ODE. This argument uses the rate of each reaction event and its stoichiometric effect; it does not require the rate law to have the same powers as the stoichiometric coefficients. In particular, a coefficient of 2 means that two units of a species are consumed or produced per event. It does not, by itself, determine the reaction order.

Worked Example: A Simple Conversion

Consider the hypothetical reaction \(A\rightarrow B\), with rate law \(v=k[A]\). Suppose \(k=0.30\ \text{min}^{-1}\), and at the time of interest \([A]=8.0\ \mu\text{M}\). Construct the ODEs and calculate the initial rates of change at these concentrations.

One reaction event consumes one \(A\) and produces one \(B\), so \(\nu_A=-1\) and \(\nu_B=+1\). The reaction rate is

v=(0.30\ \text{min}^{-1})(8.0\ \mu\text{M})=2.4\ \mu\text{M}\,\text{min}^{-1}.

Multiplying the rate by each stoichiometric change gives

\frac{d[A]}{dt}=-v=-k[A],\qquad \frac{d[B]}{dt}=+v=+k[A].

At the stated concentration, \(d[A]/dt=-2.4\ \mu\text{M}\,\text{min}^{-1}\) and \(d[B]/dt=+2.4\ \mu\text{M}\,\text{min}^{-1}\). The negative sign means \(A\) is being depleted; the positive sign means \(B\) is being formed. Since one molecule of \(A\) becomes one molecule of \(B\), the sum \([A]+[B]\) has zero derivative in this closed, one-to-one conversion model.

Association: One Rate, Several Species Equations

For an association reaction, the same reaction rate appears in the ODE for every participating species, but with a sign determined by that species’ stoichiometric coefficient. For example, the elementary step \(A+B\rightarrow C\) has mass-action rate \(v=k_{\mathrm{on}}[A][B]\). One event consumes one \(A\), consumes one \(B\), and produces one \(C\). Thus its contributions to the three ODEs are \(-v\), \(-v\), and \(+v\), respectively.

Worked Example: Constructing ODEs for Association

Consider the hypothetical association \(A+B\rightarrow C\), with \(v=k_{\mathrm{on}}[A][B]\). Let \(k_{\mathrm{on}}=0.020\ \mu\text{M}^{-1}\text{min}^{-1}\), and let \([A]=5.0\ \mu\text{M}\), \([B]=3.0\ \mu\text{M}\), and \([C]=0.50\ \mu\text{M}\). Find the ODEs and evaluate the derivatives.

The stoichiometric changes are \(\nu_A=-1\), \(\nu_B=-1\), and \(\nu_C=+1\). The reaction rate is

v=(0.020\ \mu\text{M}^{-1}\text{min}^{-1})(5.0\ \mu\text{M})(3.0\ \mu\text{M}) =0.30\ \mu\text{M}\,\text{min}^{-1}.

Therefore,

\frac{d[A]}{dt}=-k_{\mathrm{on}}[A][B],\qquad \frac{d[B]}{dt}=-k_{\mathrm{on}}[A][B],\qquad \frac{d[C]}{dt}=+k_{\mathrm{on}}[A][B].

At the stated concentrations, \(d[A]/dt=-0.30\ \mu\text{M}\,\text{min}^{-1}\), \(d[B]/dt=-0.30\ \mu\text{M}\,\text{min}^{-1}\), and \(d[C]/dt=+0.30\ \mu\text{M}\,\text{min}^{-1}\). The units are consistent: multiplying a rate constant with units \(\mu\text{M}^{-1}\text{min}^{-1}\) by two concentrations gives concentration per time. The value of \([C]\) does not enter this rate law because \(C\) is not a reactant in this directional step.

Branching: Add Every Pathway Separately

A species may be consumed by more than one reaction. Its ODE must then include a negative contribution from each consuming pathway. Products are handled in the same way: contributions from every reaction that forms a product are added. A useful practice is to write one term per reaction before simplifying the equation.

Worked Example: A Species With Two Competing Pathways

Suppose a hypothetical species \(A\) can undergo either \(A\rightarrow B+C\), with \(v_1=k_1[A]\), or \(A\rightarrow D\), with \(v_2=k_2[A]\). Let \(k_1=0.10\ \text{min}^{-1}\), \(k_2=0.050\ \text{min}^{-1}\), and \([A]=4.0\ \mu\text{M}\). Construct the equations and evaluate their derivatives.

For the first pathway, one \(A\) is consumed and one each of \(B\) and \(C\) is produced. For the second, one \(A\) is consumed and one \(D\) is produced. The rates are

v_1=(0.10\ \text{min}^{-1})(4.0\ \mu\text{M})=0.40\ \mu\text{M}\,\text{min}^{-1}, \qquad v_2=(0.050\ \text{min}^{-1})(4.0\ \mu\text{M})=0.20\ \mu\text{M}\,\text{min}^{-1}.

Summing the contributions gives

\frac{d[A]}{dt}=-v_1-v_2=-(k_1+k_2)[A],\qquad \frac{d[B]}{dt}=v_1,\qquad \frac{d[C]}{dt}=v_1,\qquad \frac{d[D]}{dt}=v_2.

At the stated concentration, \(d[A]/dt=-0.40-0.20=-0.60\ \mu\text{M}\,\text{min}^{-1}\); \(d[B]/dt=+0.40\ \mu\text{M}\,\text{min}^{-1}\); \(d[C]/dt=+0.40\ \mu\text{M}\,\text{min}^{-1}\); and \(d[D]/dt=+0.20\ \mu\text{M}\,\text{min}^{-1}\). The two product concentrations \(B\) and \(C\) receive contributions only from the first pathway, whereas \(D\) receives a contribution only from the second. The total rate of loss of \(A\) is the sum of both reaction rates.

Stoichiometric Coefficients and Reliable Checks

If a reaction has a coefficient other than one, use that coefficient in the species balance. For \(2A\rightarrow B\), the reaction rate \(v\) counts reaction events per volume per time; each event consumes two \(A\). The contributions are therefore \(d[A]/dt=-2v\) and \(d[B]/dt=+v\). If the step is elementary with mass-action rate \(v=k[A]^2\), that rate law is specified separately from the stoichiometric accounting. The coefficient 2 appears in the ODE multiplier because of the amount consumed per event, while its appearance in the rate law follows from the stated elementary-step assumption.

Several checks can expose errors before the ODEs are used in a simulation:

Check signs and participation. Reactants have negative stoichiometric coefficients and products have positive coefficients. A species absent from a reaction gets no contribution from it.
Check multiplicities. A coefficient of 2 in a reaction must lead to a factor of 2 in that species’ contribution, even if the reaction rate itself is defined per reaction event.
Check dimensions. Every term in a concentration derivative must have concentration-per-time units. All reaction contributions to a given ODE must share those units.
Check pathway completeness. For a species involved in several reactions, verify that each relevant pathway contributes exactly once.
Check conservation where appropriate. A closed reaction set may conserve a weighted total of species. For the one-to-one conversion \(A\rightarrow B\), \(d([A]+[B])/dt=0\). The appropriate conserved total depends on the reaction stoichiometry and on what quantities the concentrations represent.

These checks do not establish that the mechanism or rate laws are biologically correct. They test whether the equations faithfully implement the chosen reaction list. Numerical solutions can also flag problems: negative concentrations, for example, may indicate unsuitable numerical settings, an inconsistent model formulation, or a state that has been driven beyond the physical domain. They are a warning to investigate, not by themselves proof of which cause is responsible.

Common Pitfall and Alternative Representations

Wrong reasoning: “The reaction is \(2A\rightarrow B\), so the rate law must be \(k[A]^2\), and the \(A\) ODE is just minus the reaction rate.” This combines two different steps and omits a stoichiometric factor. Correction: first use the given or justified rate law to calculate the event rate \(v\); then apply the stoichiometric changes, giving \(d[A]/dt=-2v\) and \(d[B]/dt=+v\). For an elementary reaction with the stated mass-action law \(v=k[A]^2\), this becomes \(d[A]/dt=-2k[A]^2\), not \(-k[A]^2\). The earlier material on reaction order versus stoichiometric coefficients addresses why the rate-law powers and the ODE multipliers must not be conflated.

For a small reaction list, writing each species equation directly is transparent. A stoichiometric matrix is a more compact alternative for a larger network: place each reaction in a column, each species in a row, and enter the signed coefficients. Multiplying this matrix by the vector of reaction rates automatically performs the same sum. Neither representation replaces the need to specify valid rate laws or check the reaction list. The direct method can be easier to inspect by hand; the matrix method can make omissions and repeated patterns easier to manage computationally.

Key takeaway: For every reaction, multiply its rate by each species’ signed stoichiometric coefficient, then sum those products across reactions to obtain that species’ concentration derivative.

Check Your Understanding

For each problem, identify the reaction contributions and check the signs before simplifying.

For the hypothetical reaction \(X\rightarrow Y\) with \(v=0.40[X]\), write both ODEs. If \([X]=2.0\ \mu\text{M}\), what are the derivatives?
For \(2M\rightarrow N\) with reaction rate \(v=0.10\ \mu\text{M}\,\text{min}^{-1}\), write the contributions to \(d[M]/dt\) and \(d[N]/dt\).
For \(P+Q\rightarrow R\) with \(v=k[P][Q]\), which concentration derivatives are negative, and which is positive?
Suppose \(Z\rightarrow U\) has rate \(v_1\) and \(Z\rightarrow V\) has rate \(v_2\). Write the ODE for \(Z\) and explain why both terms appear.

Solutions. (1) One \(X\) is consumed and one \(Y\) is formed, so \(d[X]/dt=-0.40[X]\) and \(d[Y]/dt=+0.40[X]\). At \([X]=2.0\ \mu\text{M}\), the rate is \(0.80\ \mu\text{M}\,\text{min}^{-1}\), giving derivatives \(-0.80\) and \(+0.80\ \mu\text{M}\,\text{min}^{-1}\). (2) Each event consumes two \(M\) and forms one \(N\), so \(d[M]/dt=-2v=-0.20\ \mu\text{M}\,\text{min}^{-1}\) and \(d[N]/dt=+v=+0.10\ \mu\text{M}\,\text{min}^{-1}\). (3) \(P\) and \(Q\) are reactants, so \(d[P]/dt=-v\) and \(d[Q]/dt=-v\); \(R\) is the product, so \(d[R]/dt=+v\). (4) \(Z\) is consumed by both reaction channels, so \(d[Z]/dt=-v_1-v_2\). Both terms are required because each pathway removes \(Z\).