Documenting States, Reactions, and Assumptions

Why a Model Specification Needs More Than a Diagram

A pathway diagram can communicate the intended biology, but it may not show exactly what is calculated. An arrow from one species to another does not, by itself, specify the rate of that process, whether material is conserved, or how much of each species is present at the start of a simulation. A usable model specification makes those choices explicit so another modeler can reconstruct the equations and understand what the equations claim about the biology.

This tutorial develops that practice using a hypothetical three-state pathway. The example is deliberately small: its purpose is to show how to document a model, not to claim that a particular biological pathway follows these kinetics. The same discipline applies when a QSP model contains many more states, compartments, and processes.

Definition: A model specification is a written account of the model’s state variables, processes and their mathematical representations, governing equations, initial conditions, and biological assumptions. It should be detailed enough to identify what is modeled and how the stated assumptions produce the equations.

A Small Pathway and Its States

Consider a hypothetical pathway with three modeled states: \(X\) is an inactive precursor, \(Y\) is an activated intermediate, and \(Z\) is a downstream product. Material moves from \(X\) to \(Y\), then from \(Y\) to \(Z\); \(Z\) is removed. An external input can add material to \(X\). Let \(X(t)\), \(Y(t)\), and \(Z(t)\) denote their concentrations at time \(t\). In this example, concentrations are measured in nM and time in hours.

A state variable is a quantity whose value changes during the simulation and is needed to determine the model’s future behavior. Naming the three states is not enough: the specification should also say what each quantity represents, what is included in it, and what its measurement unit is. Here, each state represents the total concentration of the corresponding modeled species in one well-mixed compartment. No unmodeled forms are included in those totals.

The processes are labeled \(R_1\), \(R_2\), and \(R_3\). The first two transfer material between states; the third removes material from \(Z\). For this baseline model, the corresponding fluxes are first order. The input \(q(t)\) is a rate of addition to \(X\), rather than a state or an internal reaction.

Definition: A reaction or process flux is the amount of material transferred or removed per unit time. In a concentration-based model, its units are concentration per time. A state variable has units of concentration; its time derivative and the fluxes contributing to it have units of concentration per time.

The specification defines \(k_1\), \(k_2\), and \(k_3\) as constant first-order rate constants, each with units of h\(^{-1}\). The flux laws are:

v_1(t)=k_1X(t),\qquad v_2(t)=k_2Y(t),\qquad v_3(t)=k_3Z(t).

Each flux is nonnegative when its state and rate constant are nonnegative. The first-order laws mean, for example, that doubling \(X\) while holding \(k_1\) fixed doubles the rate of the \(X\)-to-\(Y\) process. That is a modeling assumption, not a conclusion established by drawing an arrow.

Deriving the Governing Equations

For each process, record how much of each state changes when one unit of material undergoes that process. In the order \(X,Y,Z\), process \(R_1\) changes the states by \((-1,+1,0)\), \(R_2\) by \((0,-1,+1)\), and \(R_3\) by \((0,0,-1)\). These changes form the stoichiometric matrix \(S\). The state derivatives are the sum of the process contributions, plus the external input.

S= \begin{pmatrix} -1&0&0\\ 1&-1&0\\ 0&1&-1 \end{pmatrix}, \qquad \mathbf{v}= \begin{pmatrix} v_1\\v_2\\v_3 \end{pmatrix}, \qquad \frac{d}{dt} \begin{pmatrix} X\\Y\\Z \end{pmatrix} = S\mathbf{v}+ \begin{pmatrix} q(t)\\0\\0 \end{pmatrix}.

Multiplying the matrix by the flux vector gives one contribution per state: \(X\) loses \(v_1\), \(Y\) gains \(v_1\) and loses \(v_2\), and \(Z\) gains \(v_2\) and loses \(v_3\). Including the input and substituting the flux laws yields the full system:

\begin{aligned} \frac{dX}{dt}&=q(t)-k_1X,\\ \frac{dY}{dt}&=k_1X-k_2Y,\\ \frac{dZ}{dt}&=k_2Y-k_3Z. \end{aligned}

This is a direct derivation: each term follows from a stated process and its effect on a state. It also provides a quick check. Every term in an equation for a concentration must have units of concentration per time. For example, \(k_1X\) has units h\(^{-1}\) times nM, or nM/h, matching \(dX/dt\).

Model assumptions: This baseline specification assumes a single well-mixed compartment; nonnegative concentrations; irreversible transfers in the stated directions; constant first-order rate constants; a specified external input \(q(t)\); and no feedback, delay, spatial gradient, or additional source or sink. These assumptions define the scope of this model. They are not universal properties of biological pathways.

Worked Example: Completing the Baseline Specification

For a hypothetical simulation, let \(k_1=0.40\) h\(^{-1}\), \(k_2=0.25\) h\(^{-1}\), and \(k_3=0.10\) h\(^{-1}\). Let the input initially be \(q(0)=2.0\) nM/h, and set \(X(0)=5.0\) nM, \(Y(0)=1.0\) nM, and \(Z(0)=0.50\) nM. Calculate the initial fluxes and state derivatives.

The flux laws give:

\begin{aligned} v_1(0)&=(0.40\ \mathrm{h}^{-1})(5.0\ \mathrm{nM})=2.0\ \mathrm{nM/h},\\ v_2(0)&=(0.25\ \mathrm{h}^{-1})(1.0\ \mathrm{nM})=0.25\ \mathrm{nM/h},\\ v_3(0)&=(0.10\ \mathrm{h}^{-1})(0.50\ \mathrm{nM})=0.050\ \mathrm{nM/h}. \end{aligned}

Substitute these values into the governing equations:

\begin{aligned} \left.\frac{dX}{dt}\right|_{t=0} &=2.0-2.0=0\ \mathrm{nM/h},\\ \left.\frac{dY}{dt}\right|_{t=0} &=2.0-0.25=1.75\ \mathrm{nM/h},\\ \left.\frac{dZ}{dt}\right|_{t=0} &=0.25-0.050=0.20\ \mathrm{nM/h}. \end{aligned}

At the initial instant, input to \(X\) exactly balances its conversion to \(Y\), so \(X\) has zero initial slope. The other two states are increasing at that instant. This does not mean \(X\) remains constant: as the states change, their fluxes can change as well.

What a Complete Specification Records

A reader should be able to move from biological description to equations without guessing. A concise specification can organize the required information in a table and then state the equations and initial conditions separately. For the pathway above, the process table could be:

Process	State change	Flux law	Interpretation
\(R_1\)	\(X\) to \(Y\)	\(v_1=k_1X\)	First-order conversion of the inactive precursor
\(R_2\)	\(Y\) to \(Z\)	\(v_2=k_2Y\)	First-order conversion of the intermediate
\(R_3\)	Removal from \(Z\)	\(v_3=k_3Z\)	First-order loss from the modeled compartment
Input	Addition to \(X\)	\(q(t)\)	External source, specified as a function of time

The initial condition is the state of the system at the chosen time origin. In this example it is the vector \((5.0,1.0,0.50)\) nM at \(t=0\). If the simulation begins after a prior exposure or equilibration period, those initial values need to represent that modeled starting condition; they should not be silently assumed to be zero.

The assumptions should be written in terms that can be checked against the intended use. “First-order removal” is more informative than “\(Z\) decays,” because it identifies the assumed dependence of the removal rate on \(Z\). “One well-mixed compartment” tells the reader that gradients are not represented. If feedback is omitted, state that omission rather than letting the diagram imply that all biologically plausible feedback has been ruled out.

Worked Example: Checking the Equation at a Later State

Suppose the system is at \(X=8.0\) nM, \(Y=3.0\) nM, and \(Z=2.0\) nM, with the same rate constants as above. For this calculation, the external input is off, so \(q(t)=0\). Find the fluxes and the instantaneous state derivatives.

First calculate each flux:

\begin{aligned} v_1&=(0.40\ \mathrm{h}^{-1})(8.0\ \mathrm{nM})=3.2\ \mathrm{nM/h},\\ v_2&=(0.25\ \mathrm{h}^{-1})(3.0\ \mathrm{nM})=0.75\ \mathrm{nM/h},\\ v_3&=(0.10\ \mathrm{h}^{-1})(2.0\ \mathrm{nM})=0.20\ \mathrm{nM/h}. \end{aligned}

With no input, the equations give:

\begin{aligned} \frac{dX}{dt}&=0-3.2=-3.2\ \mathrm{nM/h},\\ \frac{dY}{dt}&=3.2-0.75=2.45\ \mathrm{nM/h},\\ \frac{dZ}{dt}&=0.75-0.20=0.55\ \mathrm{nM/h}. \end{aligned}

The signs match the process description: \(X\) is being consumed, while the inflow to \(Y\) exceeds its outflow and the inflow to \(Z\) exceeds its removal. This check tests the local consequences of the equations; it does not establish that the chosen rate laws are biologically adequate.

Checking Assumptions and Recognizing Failure Modes

A written specification makes several common errors easier to detect. First, check the direction and sign of every process. If \(R_1\) converts \(X\) into \(Y\), its contribution must be negative in the \(X\) equation and positive in the \(Y\) equation. Second, check dimensions: a term added to a state derivative must have units of concentration per time. Third, check whether every process that changes a state appears in that state’s equation. A diagram with an arrow into \(Y\) but no corresponding positive term in \(dY/dt\) is inconsistent.

The reaction representation also permits a conservation check. If \(q(t)=0\), add the three equations:

\frac{d}{dt}(X+Y+Z)=-k_3Z.

The internal transfers cancel, as expected. The total modeled concentration decreases only through the specified removal from \(Z\). With input present, the corresponding result is \(d(X+Y+Z)/dt=q(t)-k_3Z\). If the model instead showed a loss of total material during \(R_1\) or \(R_2\), that would conflict with the stated transfer assumptions unless an additional loss process had been defined.

A failed assumption may show up as a mismatch between intended biology and model behavior: for example, data or mechanistic knowledge may indicate saturation, feedback, or multiple compartments where the specification assumes simple first-order conversion in one well-mixed compartment. That mismatch is a reason to reconsider the formulation, not a reason to hide an extra term in the code without documenting it.

Worked Example: Documenting a Saturable Alternative

Suppose the \(X\)-to-\(Y\) process is believed to become saturated at high \(X\). Replace only its first-order flux law with \(v_1=V_{\max}X/(K_m+X)\), where \(V_{\max}=4.0\) nM/h and \(K_m=5.0\) nM. At \(X=20\) nM, calculate the saturable flux and compare it with the baseline first-order flux using \(k_1=0.40\) h\(^{-1}\).

The saturable flux is:

v_1=\frac{(4.0\ \mathrm{nM/h})(20\ \mathrm{nM})} {5.0\ \mathrm{nM}+20\ \mathrm{nM}} =\frac{80}{25}\ \mathrm{nM/h} =3.2\ \mathrm{nM/h}.

The baseline flux at the same state is:

v_1=k_1X=(0.40\ \mathrm{h}^{-1})(20\ \mathrm{nM})=8.0\ \mathrm{nM/h}.

The alternative therefore predicts 3.2 nM/h rather than 8.0 nM/h at this concentration. Its documentation must identify the saturable law and the meanings and units of \(V_{\max}\) and \(K_m\); retaining the first-order equation in the written specification would make the model ambiguous. This calculation illustrates the consequence of a different assumption, not evidence that saturation is the correct biological description.

Choosing a Useful Level of Detail

An informal diagram is useful for communicating pathway structure, but does not fully specify quantitative behavior. A list of differential equations is more precise mathematically, but can still leave the biological meaning of states, the direction of fluxes, or initial conditions unclear. A reaction table paired with equations and assumptions links those perspectives: the table explains what each term represents, and the equations state how the model evolves.

A specification is not a guarantee that the model is biologically correct. It cannot resolve uncertainty about whether the pathway is complete, whether the assumed rate law is appropriate, or whether the chosen initial state represents the situation of interest. It does, however, make those choices visible and enables a reviewer to identify what would need to change if an assumption is not suitable.

A practical final review can ask whether every state has a definition and initial value, every process has a direction and mathematical law, every equation follows from the listed processes, and every term has compatible units. Also check that stated omissions—such as no feedback or no spatial structure—are clear enough for a reader to judge the model’s intended scope.

Check Your Understanding

Answer each question using the baseline pathway specification. Short solutions follow the questions.

If \(X=6.0\) nM and \(k_1=0.40\) h\(^{-1}\), what is \(v_1\), and what are its units?
At \(X=4.0\) nM, \(Y=2.0\) nM, and \(Z=1.0\) nM, with \(q=0\), calculate \(dY/dt\).
Why must \(R_1\) appear with opposite signs in the \(X\) and \(Y\) equations under the stated transfer assumption?
What does the assumption of a well-mixed compartment exclude from this model?
With no input, what is the equation for the rate of change of \(X+Y+Z\), and which process explains that change?

Solutions. (1) \(v_1=(0.40\ \mathrm{h}^{-1})(6.0\ \mathrm{nM})=2.4\) nM/h. (2) \(dY/dt=k_1X-k_2Y=(0.40)(4.0)-(0.25)(2.0)=1.6-0.50=1.1\) nM/h. (3) \(R_1\) removes material from \(X\) and adds the same transferred amount to \(Y\), so its contributions are \(-v_1\) and \(+v_1\). (4) It excludes modeled spatial gradients or separate locations with different concentrations. (5) Adding the equations gives \(d(X+Y+Z)/dt=-k_3Z\); removal through \(R_3\) is the only internal process that reduces the modeled total.