Translating ODE Models Across Software

Why ODE Translation Requires More Than Matching Syntax

An ordinary differential equation (ODE) model is specified by its states, the rates of change of those states, its initial conditions, and the outputs calculated from the states. Translating a model between software therefore means preserving all four elements—not merely writing equations that look similar. A change from amount to concentration as the state, for example, changes the equations unless the volume factors are changed consistently.

This tutorial follows one linear two-compartment model with elimination from the central compartment. The structural model is held fixed across NONMEM, Monolix, nlmixr2, and mrgsolve. The purpose is to make their ODE expressions and structural predictions correspond; it is not to compare estimation results or to prescribe a particular data-file setup.

Definition: Let $A_1(t)$ and $A_2(t)$ be the amounts of drug in the central and peripheral compartments at time $t$. Let $CL$ be clearance from the central compartment, $V_1$ and $V_2$ the central and peripheral volumes, and $Q$ the intercompartmental clearance. With linear distribution and elimination, the model is $$ \frac{dA_1}{dt}=-\frac{CL}{V_1}A_1-\frac{Q}{V_1}A_1+\frac{Q}{V_2}A_2,\qquad \frac{dA_2}{dt}=\frac{Q}{V_1}A_1-\frac{Q}{V_2}A_2. $$ The central concentration output is $C_1(t)=A_1(t)/V_1$.

The model assumes that $CL$, $Q$, $V_1$, and $V_2$ are constant over the modeled interval; distribution between compartments and elimination from the central compartment are first-order processes; and amounts and volumes use compatible units. For example, if amounts are in mg, volumes in L, and time in h, clearances must be in L/h, concentrations are in mg/L, and each derivative is in mg/h. The model does not include elimination from the peripheral compartment.

For a bolus dose $D$ placed in the central compartment at time zero, the post-dose initial conditions are $A_1(0)=D$ and $A_2(0)=0$. Without a dose, the corresponding initial conditions might instead be $A_1(0)=A_2(0)=0$. State the conditions explicitly: an ODE expression alone does not determine the predicted trajectory.

Read the Equations as a Model Before Translating Them

The central-compartment equation has three terms. The first is elimination, with rate $CL \times C_1=(CL/V_1)A_1$. The second is distribution out of the central compartment, at rate $Q \times C_1=(Q/V_1)A_1$. The third is distribution into it, at rate $Q \times C_2=(Q/V_2)A_2$, where $C_2=A_2/V_2$. The peripheral equation has the same two distribution rates with the directions reversed.

This rate-based interpretation is a useful translation check. Each clearance multiplied by the concentration in its source compartment gives an amount-per-time rate. Both compartments exchange the same amount per unit time in opposite directions; they do not generally exchange the same fraction of their amounts per unit time, because $V_1$ and $V_2$ may differ.

Worked Example: Calculating the Initial Slopes

Suppose the hypothetical model has $CL=5$ L/h, $Q=8$ L/h, $V_1=20$ L, and $V_2=40$ L. A 100 mg bolus is placed in the central compartment, so immediately after dosing $A_1=100$ mg and $A_2=0$ mg. Calculate the initial derivatives and central concentration.

The rate constants in the amount equations are $CL/V_1=5/20=0.25$ h$^{-1}$, $Q/V_1=8/20=0.40$ h$^{-1}$, and $Q/V_2=8/40=0.20$ h$^{-1}$. Substitution gives $$ \left.\frac{dA_1}{dt}\right|_{0^+}=-(0.25+0.40)(100)+(0.20)(0)=-65\ \text{mg/h}, \qquad \left.\frac{dA_2}{dt}\right|_{0^+}=(0.40)(100)-(0.20)(0)=40\ \text{mg/h}. $$ The post-dose central concentration is $C_1(0^+)=100/20=5$ mg/L. These are useful reference values: a translated implementation using the same parameters and dose should have the same initial states, slopes, and concentration.

The total amount initially changes at $-65+40=-25$ mg/h. This is consistent with central elimination at rate $CL \times C_1=(5\ \text{L/h})(5\ \text{mg/L})=25$ mg/h.

Express the Same States in Four Tools

The equations below use the same amount states and central concentration output. Syntax details can depend on software version and model-file conventions, but the mathematical correspondence is explicit. Assume that the dose record places the bolus into the central state and that both states are zero before the dose. The event record, rather than a different ODE, creates the post-dose initial condition $A_1(0^+)=D$, $A_2(0^+)=0$.

Worked Example: Translating the ODE Expressions

Using the parameter values from the previous example, express the two derivatives and central output in each tool. In every case the required numerical model is $$ \frac{dA_1}{dt}=-0.65A_1+0.20A_2,\qquad \frac{dA_2}{dt}=0.40A_1-0.20A_2,\qquad C_1=A_1/20. $$ The derivative coefficients have units h$^{-1}$, and $C_1$ is in mg/L when the state amounts are in mg.

NONMEM. Name the compartments in the model definition, for example with $MODEL COMP=(CENTRAL) COMP=(PERIPHERAL). In $DES, use DADT(1)=-(CL/V1+Q/V1)*A(1)+(Q/V2)*A(2) and DADT(2)=(Q/V1)*A(1)-(Q/V2)*A(2). The structural output in an observation or prediction block is CP=A(1)/V1. Here A(1) and A(2) are the amounts in the two modeled compartments; the dose must be assigned to the central compartment.

Monolix. Define the inputs $CL$, $V1$, $Q$, and $V2$. In the ODE section, use ddt_A1 = -(CL/V1+Q/V1)*A1 + (Q/V2)*A2 and ddt_A2 = (Q/V1)*A1 - (Q/V2)*A2. Define the output as C1 = A1/V1 and expose C1 as the model output. The data and administration specification must place the dose in the state corresponding to $A1$.

nlmixr2. In the model block, the state derivatives can be expressed as d/dt(centr) = -(CL/V1+Q/V1)*centr + (Q/V2)*periph and d/dt(periph) = (Q/V1)*centr - (Q/V2)*periph. Define the structural prediction as cp = centr/V1. The event data must dose the central state, here named centr.

mrgsolve. Declare the amount states as central and peripheral compartments, for example with $CMT CENT PERIPH. In the ODE block use dxdt_CENT = -(CL/V1+Q/V1)*CENT + (Q/V2)*PERIPH; and dxdt_PERIPH = (Q/V1)*CENT - (Q/V2)*PERIPH;. In the table block define CP = CENT/V1; and capture that output as needed. The bolus must be assigned to CENT.

The state names differ, but each equation has the same central loss, central gain, peripheral gain, and peripheral loss. A syntax translation is complete only when the data's dose compartment and the reported output also refer to the intended states.

A Short Justification: Exchange Cancels in the Total Amount

Adding the two state equations provides an independent check on signs and process placement: $$ \frac{d(A_1+A_2)}{dt} = -\frac{CL}{V_1}A_1-\frac{Q}{V_1}A_1+\frac{Q}{V_2}A_2 +\frac{Q}{V_1}A_1-\frac{Q}{V_2}A_2 = -\frac{CL}{V_1}A_1. $$

The two intercompartmental terms in each direction cancel because distribution moves amount between compartments but does not create or remove drug. The total amount therefore decreases only through central elimination. If a purported translation leaves an intercompartmental term in the total-amount derivative, one of the transfer signs or coefficients is inconsistent. If it instead gives $-CL A_2/V_2$, elimination has been attached to the wrong compartment.

Worked Example: Finding an Incorrect Transfer Term

Suppose an implementation uses the correct central equation but mistakenly defines the peripheral equation as $dA_2/dt=(Q/V_1)A_1-(Q/V_1)A_2$. Using the hypothetical parameters above and the same post-dose states, calculate its initial peripheral slope and compare it with the intended model.

Because $A_2(0^+)=0$, this particular error does not change the initial peripheral slope: both implementations give $(8/20)(100)-(8/20)(0)=40$ mg/h. An initial-slope check alone would miss it. At a later state such as $A_1=50$ mg and $A_2=20$ mg, the correct peripheral slope is $$ (8/20)(50)-(8/40)(20)=20-4=16\ \text{mg/h}. $$ The incorrect slope is $(8/20)(50)-(8/20)(20)=20-8=12$ mg/h. It predicts too little return from the peripheral compartment because it uses $V_1$, rather than $V_2$, for the peripheral concentration.

The total-amount check also identifies the error. Adding the incorrect equations leaves an extra term $(Q/V_2-Q/V_1)A_2$, which is nonzero when $V_1\ne V_2$ and $A_2\ne0$. The correction is to use the same peripheral transfer rate $Q A_2/V_2$ as the loss from the peripheral state.

Diagnostics, Alternatives, and Limits

A practical structural comparison should use identical parameter values, initial conditions, dose placement, time units, and output definitions. Compare predicted central concentrations at selected times, including immediately after a bolus and at later times when both compartments contain drug. Do not expect two tools to agree if one reports amount and another reports concentration, or if a dose is placed in a different state.

Check units. Each coefficient multiplying an amount in the derivatives must have units of inverse time. The output $A_1/V_1$ must have concentration units.
Check the bolus limit. For a central bolus, $A_1(0^+)=D$, $A_2(0^+)=0$, and $C_1(0^+)=D/V_1$. A wrong initial concentration often indicates a volume or state mismatch.
Check directional rates. Evaluate both derivatives at a state where $A_1$ and $A_2$ are nonzero. This can reveal a wrong volume that an initial state with $A_2=0$ conceals.
Check total amount balance. The derivative of $A_1+A_2$ must equal $-CL A_1/V_1$ under this model.
Check numerical integration separately. If the equations and events match but small differences remain, examine solver settings and tolerances before changing the model. Large or systematic differences more often signal a structural, unit, event, or output mismatch.

For this linear, constant-parameter model, an analytic solution is an alternative to numerical ODE integration and can provide an additional benchmark. Software ODE solvers are useful for implementing this model in a broader modeling workflow, but the particular equations here do not require a numerical solver in principle. Neither approach can compensate for incorrect states, equations, or initial conditions. Conversely, small discrepancies due to numerical tolerances do not by themselves establish that the structural models differ.

The translation shown here is limited to two amount states, linear exchange, and elimination from the central compartment. More complex structures, different state choices, or time-dependent inputs require their own carefully matched equations and outputs. Preserve the model’s meaning first; software-specific convenience should not silently change it.

Common pitfall: Seeing $Q$ in both compartment equations and copying the same coefficient for both losses. The clearance $Q$ is shared, but the rate from each compartment is $Q$ times that compartment’s concentration. Thus the central loss coefficient is $Q/V_1$, while the peripheral loss coefficient is $Q/V_2$.

Check Your Understanding

Answer each question using the model and conventions in this tutorial. Brief solutions follow the questions.

With $CL=6$ L/h, $V_1=30$ L, $Q=9$ L/h, and $V_2=45$ L, what is the coefficient of $A_2$ in the central derivative?
A 120 mg bolus is given centrally when $V_1=24$ L. What are the post-dose initial conditions and central concentration?
Why must the intercompartmental terms cancel when the two amount derivatives are added?
An implementation reports $A_1$ as its output, while the requested observation is central concentration. What transformation is needed, and what are its units?

Solutions. (1) The coefficient is $Q/V_2=9/45=0.20$ h$^{-1}$, multiplying $A_2$ as the return term in the central equation. (2) $A_1(0^+)=120$ mg and $A_2(0^+)=0$ mg; $C_1(0^+)=120/24=5$ mg/L. (3) Distribution transfers drug between compartments without changing total amount: the central loss $-QA_1/V_1$ cancels the peripheral gain $+QA_1/V_1$, and the peripheral loss $-QA_2/V_2$ cancels the central gain $+QA_2/V_2$. (4) Divide the central amount by its volume: $C_1=A_1/V_1$, giving concentration units such as mg/L when amount is in mg and volume in L.