Translating ODE Systems Into Right-Hand-Side Functions

From a Mechanistic Balance to a Function

An ordinary differential equation model specifies how its states change at an instant. To simulate the model, its equations must be expressed as a function that takes the current time, current state, and parameter values, then returns the state derivatives. The function does not itself move the system forward in time: it supplies the rate information that a numerical integration method uses to do so.

This translation is a useful test of whether the model has been written unambiguously. Every production, conversion, or loss term in the biological description must appear with the correct sign and units in the derivative. A function can be syntactically valid and still encode the wrong balance.

A One-Compartment Turnover Model

Consider a hypothetical drug amount \(A\) in a single, well-mixed compartment and a response state \(R\) whose production is inhibited by the drug concentration. The compartment has fixed volume \(V\). Drug is eliminated by first-order clearance \(CL\); during the interval considered, there is no drug input. Drug effect is immediate and depends on the current concentration \(C=A/V\). The response is produced at baseline rate \(k_{\mathrm{in}}\), with a fraction of production inhibited by the drug, and is lost at first order with rate constant \(k_{\mathrm{out}}\).

Assume \(A\) and \(C\) are nonnegative, \(V\), \(CL\), \(IC_{50}\), \(k_{\mathrm{in}}\), and \(k_{\mathrm{out}}\) are positive, and \(0\leq I_{\max}\leq1\). The response \(R\) is also taken to be nonnegative. The fractional inhibition is an Emax-type relationship:

Here \(IC_{50}\) has concentration units, and \(I(C)\) is dimensionless. With these assumptions, \(I(C)\) lies between zero and \(I_{\max}\); therefore the remaining production fraction, \(1-I(C)\), is nonnegative.

For a unit check, if \(A\) is measured in mg, \(V\) in L, time in h, and \(CL\) in L/h, then \(C\) is in mg/L and \(CL\,C\) is in mg/h. If \(R\) is measured in response units, \(k_{\mathrm{in}}\) has units of response units/h and \(k_{\mathrm{out}}\) has units of \(1/\mathrm{h}\), so both terms in \(dR/dt\) are response units/h.

Deriving the RHS and Its Implementation

The first equation follows directly from the amount balance: elimination removes drug at rate \(CL\,C\). Since \(C=A/V\) and \(V\) is fixed, this rate can also be written \(CL\,A/V\). The minus sign indicates loss from the compartment. The response equation adds its two balance terms: the drug-modified production rate is \(k_{\mathrm{in}}(1-I(C))\), and first-order loss removes response at rate \(k_{\mathrm{out}}R\). Substituting the inhibition function gives the form used by the function:

This substitution is also a short check on how the model should be coded: compute concentration from amount, use that concentration in the inhibition term, and subtract the response loss. The derivative function should use the current values of \(A\) and \(R\) each time it is called. It should return both derivatives together, keeping their meaning and units explicit.

In Python-like notation, a named-state implementation can be written as follows. The function receives time \(t\), a state record with fields “A” and “R,” and a parameter record with fields “V,” “CL,” “kin,” “kout,” “Imax,” and “IC50.” It calculates \(C=\text{state["A"]}/\text{p["V"]}\), then calculates \(\text{inhibition}=\text{p["Imax"]}*C/(\text{p["IC50"]}+C)\). It sets \(\text{dA}=-\text{p["CL"]}*C\) and \(\text{dR}=\text{p["kin"]}*(1-\text{inhibition})-\text{p["kout"]}*\text{state["R"]}\), and returns the named record \(\{"A":\text{dA},"R":\text{dR}\}\).

Time \(t\) is part of the general function interface because models may depend explicitly on time. In this particular model, the right-hand side has no explicit time dependence: if the state and parameters are unchanged, its returned derivatives are unchanged. The function still accepts \(t\), but does not need to use it. The returned fields represent rates, not updated values; for example, returning \(A+\text{dA}\) would not be an RHS implementation.

Checks That Connect the Code to the Biology

Before using an RHS function in a larger model, check that each expression corresponds to the stated balance. A practical review follows the terms from model to implementation rather than relying on whether the code runs:

The model also supplies simple limiting checks. When \(A=0\), concentration and inhibition are zero, so the drug derivative is zero and response production is at baseline. When \(C\) is large relative to \(IC_{50}\), inhibition approaches \(I_{\max}\), rather than exceeding it. At any fixed concentration, the response derivative is zero when production balances loss:

In the absence of inhibition, this balance reduces to \(R=k_{\mathrm{in}}/k_{\mathrm{out}}\). These checks do not prove that every implementation detail is correct, but they can reveal a misplaced sign, a missing denominator, or use of the wrong state in a term.

What the RHS Does—and Does Not Do

An RHS function provides an instantaneous rate. For example, \(dA/dt=-12\) mg/h at a particular state does not mean that the amount after one hour is exactly \(108\) mg. The rate can change as the state changes. A numerical integration method repeatedly evaluates the RHS and combines those evaluations to approximate the state trajectory. That integration procedure is distinct from the derivative function itself.

This distinction also clarifies the main alternative to an RHS implementation. One could try to update the state directly with a finite time step, such as \(A_{\mathrm{new}}=A+\Delta t(dA/dt)\). That is an integration step, not a replacement definition of the model’s derivative; its accuracy depends on the step size and the changing rate. Keeping the mechanistic RHS separate makes it possible to inspect and test the balance independently of a particular time-stepping method.

The function shown here is deliberately limited to a fixed-volume compartment, first-order elimination, immediate concentration-dependent inhibition, and first-order response loss. It does not represent changing volume, delayed drug effect, saturable elimination, or drug input. Such mechanisms require different balance terms and, when appropriate, additional states or explicit time dependence. Do not add them implicitly through a code shortcut: state the biological assumption and derive the corresponding term first.

Common Mistakes and Failure Modes

• Returning an updated state instead of a derivative. The RHS returns \(dA/dt\) and \(dR/dt\), with units per time. A state plus a rate is not the requested output.
• Using amount where concentration is required. The inhibition law uses \(C=A/V\). Substituting \(A\) directly makes the expression dimensionally inconsistent and changes the modeled effect.
• Omitting the loss sign. First-order response loss subtracts \(k_{\mathrm{out}}R\). A positive sign reverses the biological effect.
• Mixing units or scales. \(IC_{50}\) and \(C\) must use the same concentration units, and the clearance and volume units must make \(CL/V\) an inverse time.
• Assuming successful execution verifies the model. A function may run while using a wrong rate constant or state. Hand calculations, units, signs, and boundary checks test meaning as well as syntax.

A RHS function is most useful when it is a direct, readable expression of the model’s balances. For this turnover example, every returned derivative can be traced to a named biological term, and each can be checked at a specific state. That traceability is more valuable than compressing the code into a single opaque expression.