Ordinary differential equations (ODEs)

Introduction

Ordinary differential equations (ODEs) can be implemented in the EQUATION: block of the [LONGITUDINAL] section. The ODEs describe a dynamical system and are defined by a set of equations for the derivative of each variable, the initial conditions, the starting time and the parameters.

 

Derivatives

The keyword ddt_ in front of a variable name, as ddt_x and ddt_y, defines the derivative of the ODE system with respect to time. The variable names denote the components of the solution. These variables are defined at the whole [LONGITUDINAL] section level through their derivatives. The derivatives themselves aren’t variables but keywords, and cannot be referenced by other equations nor be defined under a conditional statement.

 

Initial conditions

The keyword t_0 defines the initial time of the differential equation, while the variable names with suffix _0 define the  initial conditions of the system. More precisely, in a differential equation for the variable x, x is defined at t = t0  by the initial condition x_0, i.e. x(t0)=x_0 and by the differential equation after t0.





 

Linking administration with ODEs

Mlxtran provides the PK macro depot() that can be used to link the administration from a data set with the model. Details and examples how to use depot() for the administration are given here.

 

Example

We consider in this example a first order differential equation on x defined by:

\left\{\begin{array}{l} \dot{x} = -kx \\ x(t)=V,~~\textrm{for}~~t\leq10\end{array}\right.

where (V,k)=(10,.05), this ODE is typically a linear elimination process of a volume. This model is implemented in the file mlxt_ode.txt:

DESCRIPTION: Model description with a simple ODE representing a linear elimination precess of a volume
[LONGITUDINAL]
input = {V, k}

EQUATION:
t_0 = 10
x_0 = V
ddt_x = -k*x

 

Stiff ODEs

A numerical algorithm providing an efficient and stable way to solve the ODEs and DDEs is used to compute the solution of the dynamical models. Different numerical algorithms can be used to solve the ODE depending on the properties of the ODE system (Adams methods for non stiff ODEs, and Backward Differentiation Formulas methods for stiff ODEs). The right numerical algorithm must be selected in order to get an accurate and stable solution within a reasonable computational time. The stiffness of an ODE systems is one of the main properties that influences the choice of the numerical algorithms. Stiff ODEs can lead to very long computational times if an inappropriate algorithm is used. Stiff ODEs are typical for models that include processes taking place at different time scales. For example models with reaction and diffusion processes are known to be stiff. Reactions can occur at the seconds time scale while diffusion processes can take hours. For such cases a dedicated numerical algorithm for stiff ODEs is available.

A solver for stiff ODEs can be selected in Mlxtran using the keyword odeType. There are two possible choices for odeType:

; ode is considered as stiff
odeType = stiff 
; ode is considered as non-stiff (default)
odeType = nonStiff

This command has to be defined in the [LONGITUDINAL] section. When odeType is not specified then the nonStiff ODE solver is used by default.

Note that:

  • For the DDEs, the keyword odeType is ineffective.
  • For many ODEs the nonStiff solver will provide the best solution. However, it is good practice to always test both solvers to see which one is providing the shorter computational time.

ODE solver tolerances

The stiff and non-stiff solvers compute the solution of the ODE system up to a certain precision, which is defined through the tolerance. The values by default are the following:

  • non-stiff solver: absolute tolerance 1e-6, relative tolerance 1e-3
  • stiff solver: absolute tolerance 1e-9, relative tolerance 1e-6

The tolerances can be modified by specifying the desired value in the structural model file using the keywords odeAbsTol and odeRelTol, for instance:

EQUATION:
odeAbsTol = 1e-12
odeRelTol = 1e-9

The smaller the tolerance, the closer the computed prediction will be from the true solution. The absolute tolerance indicates how much the absolute difference between the computed prediction \(y_{solver}\) and the true solution \(y_{true}\) can be, i.e \(|y_{solver}-y_{true}|<\textrm{odeAbsTol}\). The relative tolerance is related to the number of digits which are correct, i.e  \(\frac{|y_{solver}-y_{true}|}{|y_{true}|}<\textrm{odeRelTol}\).

Note that the right hand side needs to be a number (such as odeAbsTol=1e-12 or odeAbsTol=0.000001) and cannot be an expression (odeAbsTol=10^(-9) is not accepted).

Rules and Best Practices





  • We encourage the users to define the initial conditions explicitly, although the default initial value is x_0=0.
  • When no starting time t0 is defined in the Mlxtran model for Monolix then by default t0 is selected to be equal to the first dose or the first observation, whatever comes first.
  • If t0 is defined, a differential equation needs to be defined.
  • Initial conditions can depend on
    • time
    • an input parameter. Therefore, initial condition values can also be estimated in Monolix.
    • a regressor value. In that case, again, a differential equation needs to be defined.
  • Only first order differential equation can be used with Mlxtran. If you want to use a differential equation of the form \frac{d^{n}x}{dt^n}=f(x,\phi), you need to transform the higher order differential equation into a first order differential equations by defining n successive differential equations as proposed in the following “Best Practices Ex. 1: Second order differential equation”.
  • The derivatives of a ODE system cannot be conditional, i.e. ddt_ cannot be used within an if statement. Instead, conditional intermediate variables for the values of the derivatives should be used as shown in “Best Practices Ex. 2: IF statement”. A conditional statement can be built by combining the keywords if, elseif,  else and end. Several elseif keywords can be chained, and the conditions are exclusive in sequence. A default value can be provided using the keyword else, but also as a simple definition preceding the conditional structure. An unspecified conditional value is null. The enclosed variables are not local variables, which means that a variable with remaining conditional definitions is still incomplete and cannot be referenced into another definition, even under the same condition.

Best Practices Ex. 1: Second order differential equation
If you want to use a second order differential equation of the type \ddot{x} +\dot{x}+x = 0, with (x(0),\dot{x}(0))=(a,b), the structural model in Mlxtran writes

[LONGITUDINAL]
input = {a, b}

EQUATION:
t0 = 0
x_0 = a
dx_0 = b
ddt_x = dx
ddt_dx = -x-dx

Best Practices Ex. 2: IF statement
If a parameter value depends on a condition, it can be written with in Mlxtran with the usual if-else-end syntax. For example, if a parameter c should be 1 if t<=10, 2 if t<20 and 3 otherwise, write

if t<=10
   c = 1
elseif t<20
   c = 2
else
   c = 3 
end

Derivatives cannot be written within an if-else-end structure. Instead, intermediate variables should be used. For example, if a derivative of x over time should be 1 if t<=10 and 2 otherwise, write

if t<=10
   dx = 1
else
   dx = 2
end
ddt_x = dx