Encoding PBPK Compartments as an ODE State Vector

From Compartments to a State Vector

A PBPK model describes how drug moves among physiological compartments over time. To integrate such a model numerically, the compartment quantities must be represented in a fixed order as a single state vector, and the model must return the rate of change of every state in that same order. A correct set of physiological equations can still produce incorrect predictions if one derivative is assigned to the wrong state or uses inconsistent units.

This tutorial uses amounts as the states and a small, deliberately simplified three-compartment model to make the indexing explicit. Its transfer coefficients are illustrative first-order rate constants; they are not calculations of blood-flow-limited exchange. The purpose is to define a reliable state-vector and derivative-function contract.

For a model with \(m\) compartments, choose an explicit indexing convention, such as \(y_1=A_{\mathrm{central}}, y_2=A_{\mathrm{liver}}, y_3=A_{\mathrm{kidney}}\). Here \(A_i\) denotes drug amount in compartment \(i\). If amount is measured in milligrams and time in hours, each derivative has units of milligrams per hour. A transfer rate constant has units of inverse hours.

Specify the State, Parameters, and Initial Conditions

Before writing the derivative function, make a mapping from each array index to one named quantity. State the units as well as the compartment names. The same mapping must be used when defining initial conditions, calculating derivatives, and extracting model outputs.

An initial condition supplies the state at the start time \(t_0\): \(y(t_0)=y_0\). It is not a derivative or a dose rate. For a dose-free interval that starts just after a bolus placed in the central compartment, one possible initial condition is \(y_0=(100,0,0)^\mathsf{T}\) mg. A different dosing history, or a different starting time, may require different initial amounts.

In the simple model below, transfer from the central compartment to each tissue is first order in the central amount, and return transfer is first order in the tissue amount. Elimination occurs from the liver and kidney states. Assume fixed rate constants, fixed compartment volumes, no input during the interval, and nonnegative initial amounts. These assumptions define this example; more detailed PBPK models may require additional states and processes.

Write the Derivative Function

Let \(k_{12}\) denote transfer from central to liver, \(k_{21}\) return from liver to central, \(k_{13}\) transfer from central to kidney, and \(k_{31}\) return from kidney to central. Let \(k_{20}\) and \(k_{30}\) denote elimination from liver and kidney, respectively. Each rate constant has units \(\mathrm{h}^{-1}\). Under first-order kinetics, for example, the central-to-liver transfer rate is \(k_{12}A_1\), in \(\mathrm{mg/h}\).

Each state’s derivative is the sum of rates entering that compartment minus the rates leaving it. Therefore, for a dose-free interval,

The derivative function returns these three expressions in the exact state order. The central derivative is not just a sum of losses: it also includes return from both tissues. The liver and kidney derivatives include their incoming central transfers, their return transfers, and their elimination losses. All terms in each equation have units \(\mathrm{mg/h}\).

Check the Equations Before Integrating

A useful verification is to add the three amount derivatives. Internal transfers should cancel because each transfer leaving one modeled compartment enters another. In this example, the total drug amount is \(A_{\mathrm{tot}}=A_1+A_2+A_3\). Adding the equations gives

For example, the term \(-k_{12}A_1\) in the central equation cancels with \(+k_{12}A_1\) in the liver equation. The remaining terms are elimination losses. This equality is a direct check on the written equations: with no input, total amount cannot increase in this model, and it decreases only through the specified elimination pathways. At the initial state in the worked example, both eliminating amounts are zero, so the sum of the derivatives is \(-30+20+10=0\ \mathrm{mg/h}\), as expected.

This check depends on the stated model boundaries: all internal transfer destinations are included, and there is no unmodeled input or output. If a process moves drug outside the modeled system, it should appear as an external loss; if drug enters, it should appear as an input. Conservation checks must reflect what the states actually represent.

Amounts, Concentrations, and Model Outputs

An amount state is often convenient for writing transfer balances. Concentrations can be derived from amounts when the compartment volumes are known. With constant volume \(V_i\), concentration in compartment \(i\) is \(C_i=A_i/V_i\). If \(A_i\) is in mg and \(V_i\) is in L, \(C_i\) is in mg/L. The state vector can remain in amounts even when reported outputs are concentrations.

Implementation Choices and Failure Modes

The main alternative is to store concentrations, rather than amounts, as the states. Either representation can be valid, but the derivative equations must match it. For amount states, transfer rates can be expressed directly in amount per time. For concentration states, transfer terms generally require volume factors so that gains and losses are expressed in concentration per time. Mixing an amount-based equation with concentration states, or vice versa, produces dimensional errors even if the numerical values look plausible.

A derivative function should implement the model equations for the current time and state, using a stable parameter-to-rate mapping and returning one derivative per state. It should not silently reorder states, use a different unit for a parameter than the equations assume, or treat an initial dose amount as a continuous input rate. This tutorial’s equations assume a dose-free interval after the initial condition; other input patterns require their own explicit terms and timing.

Common warning signs include a derivative vector of the wrong length, an unexpected sign at a simple test state, a total amount that grows despite having no input, or a concentration that results from dividing by the wrong compartment volume. A mass-balance check may not detect every error: swapping two tissue indices can preserve the total amount while assigning predictions to the wrong organs. Verify both arithmetic and the state-to-compartment mapping.

• Index mismatch: The mistaken reasoning is that all components are interchangeable because they share units. They are not; each position has a defined physiological meaning. Check the index map for the initial state, derivative, volumes, and outputs.
• Inconsistent units: A rate constant in \(\mathrm{h}^{-1}\) multiplied by an amount in mg gives \(\mathrm{mg/h}\). Combining that term with a rate in \(\mathrm{mg/min}\) without conversion is invalid.
• Confusing amount and concentration: An amount of 30 mg is not a concentration of 30 mg/L. Divide by the appropriate volume only when calculating concentration.
• Missing a transfer destination: If a modeled transfer leaves one state but does not enter another modeled state, the total-amount check may reveal an unintended loss. Decide whether that loss is a real external process or an omitted state.

A state-vector representation is an implementation structure, not a guarantee that the physiological model is complete or correct. This reduced example omits many processes and does not establish physiological parameter values. The representation is useful because it makes the relationship between model equations, numerical inputs, and reported compartment quantities explicit.