Dynamical systems driven by ODEs

Autonomous dynamical systems

Consider a time-varying system $A(t)=(A_1(t),A_2(t),\ldots A_J(t))$ defined by a system of Ordinary Differential Equations (ODE)

\begin{equation} \label{ode1_model} \dot{A} = F(A(t)) \end{equation}

where $\dot{A}(t)$ denotes the vector of derivatives of $A(t)$ with respect to $t$: \begin{equation} \label{ode2_model} \left\{ \begin{array}{lll} \dot{A}_1(t) & = & F_1(A(t)) \\ \dot{A}_2(t) & = & F_2(A(t)) \\ \vdots & \vdots & \vdots \\ \dot{A}_L(t) & = & F_L(A(t)) \end{array} \right. \end{equation}

Notations:

• let $A_0 = A(t_0)$ be the initial condition of the system defined at the initial time $t_0$,
• let $A^{\star}$ be the solution of the system at equilibrium: $F(A^{\star}) =0$

A basic model

We assume here that there is no input: \begin{eqnarray*} A(t_0) &= &A_0 \\ \dot{A}(t) &= &F(A(t)) \, \ t \geq t_0 \end{eqnarray*}

Example: A viral kinetic (VK) model. In this example, the data file contains the viral load:

ID TIME VL 1 -5 6.5 1 -2 7.1 1 1 6.3 1 5 4.2 1 12 2.1 1 20 0.9 $\vdots$ $\vdots$ $\vdots$

Consider a basic VK model with $A=(N,I,V)$ where $N$ is the number of non infected target cells, $C$ the number of infected target cells and $V$ the number of virus.

After infection and before treatment, the dynamics of the system is described by this ODE system: \begin{equation} \label{vk1} \left\{ \begin{array}{lll} \dot{N}(t) & = & s - \beta \, N(t) V(t) - d N(t) \\ \dot{I}(t) & = & \beta \, N(t)\, V(t) - \delta \, I(t) \\ \dot{V}(t) & = & p I(t) - c \, V(t) \\ \end{array} \right. \end{equation}

The equilibrium state of this system is $A^{\star} = (N^{\star} , I^{\star} , V^{\star})$, where

\begin{equation} \label{eq1} N^{\star} = \frac{\delta \, c}{ \beta \, p} \quad ; \quad I^{\star} = \frac{s - d\, N^{\star}}{ \delta} \quad ; \quad V^{\star} = \frac{ p \, I^{\star} }{c}. \end{equation}

Assume that the system has reached the equilibrium state $A^{\star}$ when the treatment starts at time $t_0=0$. The treatment inhibits the infection of the target cells and blocks the production of virus. The dynamics of the new system is described with this new ODE system:

\begin{equation} \label{vk2} \left\{ \begin{array}{lll} \dot{N}(t) & = & s - \beta(1-\eta) \, N(t) \, V(t) - d\, N(t) \\ \dot{I}(t) & = & \beta(1-\eta) \, N(t) \, V(t) - \delta \, I(t) \\ \dot{V}(t) & = & p(1-\varepsilon) \, I(t) - c \, V(t) \\ \end{array} \right. \end{equation}

where $0<\varepsilon <1$ and $0 < \eta < 1$.

The initial condition and the dynamical system are described in the MDL EQUATION with MLXTRAN):

EQUATION

T_0 = 0

N_0 = delta*c/(beta*p);

I_0 = (s-d*N)/delta

V_0 = p*I/c

DDT_N = s - beta*(1-eta)*N*V - d*N

DDT_I = beta*(1-eta)*N*V - delta*I

DDT_V = p*(1-epsilon)*I - c*V

Remark1: Here, $T_0 = 0$ means that the system is constant and is $A^{\star}$, defined in the script by $(N_0, I_0, V_0)$, for any $t<0$.

Remark2: If the initial condition is not given in the model, it is assumed to be 0.

Piecewise defined dynamical systems

More generally, we can consider input-less systems which are piecewise defined: there exists a sequence of times $t_0< t_1< ...<t_K$ and functions $F^{(1)}, F^{(2)},\ldots,F^{(K)}$ such that \begin{eqnarray*} A(t_0) &= &A_0 \\ \dot{A}(t) &= &F_k(A(t)) \ , \ t_{k-1} \leq t \leq t_{k} \\ \end{eqnarray*}

Example: viral kinetic model. We assume here that a first treatment which blocks the production of virus starts first at time $T_{Start1}$, then a second treatment which inhibits the infection of the target cells starts at time $T_{Start2}$. Both treatments stop at time $T_{Stop}$. The values of the switching times $(T_{Start1}, T_{Start2}, T_{Stop})$ are part of the data and then should be contained in the datafile itself. Using the NONMEM format for example, a column EVENT is necessary in the dataset to describe this information EVENT is an extension of the EVID (Event Identification) column used by NONMEM and which is limited to some very specific events). In the following example, $T_{Start1}=0$ is used as the reference time, $T_{Start2}=20$ and $T_{Stop}=200$:

ID TIME VL EVENT 1 -5 6.5 . 1 -2 7.1 . 1 0 . Start1 1 5 5.2 . $\vdots$ $\vdots$ $\vdots$ $\vdots$ 1 18 4.6 . 1 20 . Start2 1 25 2.3 . $\vdots$ $\vdots$ $\vdots$ $\vdots$ 1 175 1.4 . 1 200 . Stop 1 250 2.8 . $\vdots$ $\vdots$ $\vdots$ $\vdots$

We will consider the same viral kinetics model defined above. This system is now piecewise defined:

• before $T_{Start1}$, $A(t) = A^{\star}$, where $A^{\star}$ is the equilibrium state defined in \ref{eq1}
• between $T_{Start1}$ and $T_{Start2}$,

\begin{equation} \label{vk3} \left\{ \begin{array}{lll} \dot{N}(t) & = & s - \beta \, N(t) \, V(t) - d\, N(t) \\ \dot{I}(t) & = & \beta \, N(t) \, V(t) - \delta \, I(t) \\ \dot{V}(t) & = & p(1-\varepsilon) \, I(t) - c \, V(t) \\ \end{array} \right. \end{equation}

• between $T_{Start2}$ and $T_{Stop}$, the system is governed by the ODES described in \ref{vk2}

\begin{equation} \label{vk3bis} \left\{ \begin{array}{lll} \dot{N}(t) & = & s - \beta(1-\eta) \, N(t) \, V(t) - d\,N(t) \\ \dot{I}(t) & = & \beta(1-\eta) \, N(t) \, V(t) - \delta \, I(t) \\ \dot{V}(t) & = & p(1-\varepsilon) \, I(t) - c \, V(t)\\ \end{array} \right. \end{equation}

• after $T_{Stop}$, the system smoothly returns to its original state governed by the ODES described in \ref{vk1}

\begin{equation} \label{vk4} \left\{ \begin{array}{lll} \dot{N}(t) & = & s - \beta (1-\eta \, e^{-k_1 (t-T_{Stop})})\, N(t) \, V(t) - d\,N(t) \\ \dot{I}(t) & = & \beta (1-\eta \, e^{-k_1 (t-T_{Stop})}) \, N(t) \, V(t) - \delta \, I(t) \\ \dot{V}(t) & = & p(1-\varepsilon\, e^{-k_2 (t-T_{Stop})}) \, I(t) - c \, V(t) \\ \end{array} \right. \end{equation}

We have seen that the information about switching times is given in the data set. The different dynamical systems are described in the MDL (in a block EQUATION and using the statement SWITCH with MLXTRAN). We only show the blocks VARIABLES and EQUATION of the code:

VARIABLES  ID, TIME, VL use=DV, EVENT list=(Start1, Start2, Stop) 
---
EQUATION
---
SWITCH
   CASE T < T_Start1
       N = delta*c/(beta*p);
       I = (s-d*N)/delta
       V = p*I/c
   CASE T_Start1 < T < T_Start2
       DDT_N = s - beta*N*V - d*N
       DDT_I = beta*N*V - delta*I
   CASE T_Start2 < T < T_Stop
       DDT_N = s - beta*(1-eta)*N*V - d*N
       DDT_I = beta*(1-eta)*N*V - delta*I
       DDT_V = p*(1-epsilon)*I - c*V
   CASE T > T_Stop
       DDT_N = s - beta*(1-eta*exp(-k1*(T-T_Stop)))*N*V - d*N
       DDT_I = beta*(1-eta*exp(-k1*(T-T_Stop)))*N*V - delta*I
       DDT_V = p*(1-epsilon*exp(-k2*(T-T_Stop)))*I - c*V
END

Remark 1: Here, EVENT is a reserved variable name. Then, the information in the column EVENT is recognized as a succession of events. Furthermore, the times of the events Start1,Start2 and Stop are automatically created as T_Start1, T_Start2 and T_Stop.

Remark 2: In this particular example, the dynamical system is described by parameters $\beta$ and $p$ whose definition switches. Then, the same model could be encoded as follows:

EQUATION
T_0 = T_Start1
N_0 = delta*c/(beta*p);
I_0 = (s-d*N)/delta
V_0 = p*I/c
---
SWITCH
  CASE T_Start1 < T < T_Start2
    be = beta
    pe = p*(1-epsilon)
  CASE T_Start2 < T < T_Stop
    be = beta*(1-eta)
    pe = p*(1-epsilon)
  CASE T > T_Stop
    be = beta*(1-eta*exp(-k1*(T-T_Stop)))
    pe = p*(1-epsilon*exp(-k2*(T-T_Stop)))
END
--- 
DDT_N = s - be*N*V - d*N
DDT_I = be*N*V - delta*I
DDT_V = pe*I - c*V

Dynamical systems with source terms

Consider now the system (\ref{ode1_model}) with inputs: \begin{equation} \label{ode3_model} \dot{A}(t) = F(A(t),u(t)) \end{equation} where: \begin{equation} \label{ode4_model} \left\{ \begin{array}{lll} \dot{A_1}(t) & = & F_1(A(t),u_1(t)) \\ \dot{A_2}(t) & = & F_2(A(t),u_2(t)) \\ \vdots & \vdots & \vdots \\ \dot{A_L}(t) & = & F_J(A(t),u_L(t)) \end{array} \right. \end{equation} The input $u(t)=(u_1(t),u_2(t),\ldots,u_J(t))$ of the system can be defined either in the datafile (doses in a PK model for instance) or in the model.

Piecewise constant inputs

We assume an additive model for the inputs: for any $1\leq \ell \leq L$ \begin{equation} \dot{A_{\ell}}(t) = F_\ell(A(t)) + u_\ell(t) \end{equation} Then, there exists a sequence of times $(\tau_{\ell,j})$, durations $(d_{\ell,j})$ and amounts $(a_{\ell,j})$ such that \begin{equation}\label{inf1bis} u_\ell(t) = \left\{\begin{array}{ll} \frac{a_{\ell,j}}{d_{\ell,j}} & \textrm{if } \tau_{\ell,j} \leq t \leq \tau_{\ell,j}+d_{\ell,j} \\ 0 & \textrm{otherwise} \end{array} \right. \end{equation}

Example: Consider an IV infusion with two compartments \begin{equation} \left\{ \begin{array}{lll} \dot{A_1} & = & -k\,A_1(t) - k_{12}A_1(t) + k_{21}A_2(t) + u_1(t)) \\ \dot{A_2} & = & k_{12}A_1(t) - k_{21}A_2(t) \end{array} \right. \end{equation} input $u_1$ is the infusion rate defined in the dataset:

ID TIME AMT TINF DV 1 0 10 3 . 1 2 . . 21 1 4 . . 63 1 6 . . 48 1 12 15 2 . 1 15 . . 72 1 18 . . 39 $\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$

Here, AMT is the amount and TINF the infusion duration. Alternatively, the infusion rate AMT/TINF could be provided in a column RATE. This is a very basic situation with only one type of administration. Then, there is no need of any additional column DPT (depot compartment) in the datafile to distinguish different target depots for different types of administration, and the model reduces to:

VARIABLES  ID, TIME, AMT, TINF, DV
---
EQUATION
DDT_Ac = -k*Ac - k12*Ac + k21*Ap
DDT_Ap = k12*Ac - k21*Ap

If no additional information about the input is given in the model, the default is to assume that the input goes to the first component of the ODE system.

Assume now that the target compartment is not the first component of the system. Then, it is mandatory to associate in the model a component with the target compartment. As an example, consider the same datafile as before, but assume that the ODE system has been permuted, then the model should define the second component of the system as the target component:

EQUATION
-----
INPUT(CMT=2)
DDT_Ap = k12*Ac - k21*Ap
DDT_Ac = -k*Ac - k12*Ac + k21*Ap

Spike inputs

We consider an input-less dynamical system: for any $1\leq \ell \leq L$, \begin{equation} \dot{A_\ell}(t) = F_\ell(A(t)) \end{equation}

Spike inputs means that there exists a sequence of times $(\tau_{\ell,j})$ and amounts $(a_{\ell,j})$ such that \begin{equation} \label{bol1} A_\ell(\tau_{\ell,j}) = A_\ell(\tau_{\ell,j}^{-}) + a_{\ell,j} \end{equation} In other words, the amount $a_{\ell,j}$ is added to the component $A_\ell$ at time $\tau_{\ell,j}$.

Example: Consider an IV bolus. The input is given by the column AMT in the dataset:

ID TIME AMT DV 1 0 10 . 1 2 . 81 1 4 . 63 1 6 . 48 1 12 15 . 1 15 . 72 1 18 . 39 $\vdots$ $\vdots$ $\vdots$ $\vdots$

The model is exactly the same model defined above for an infusion, but without any column RATE or TINF in the datafile, spike inputs are assumed.

Inputs defined in the model

Only some very basic inputs can be directly derived from the information in the datafile. More complex inputs should be defined in the model, or using some external forcing function

Example: Consider now a 2 compartments model \begin{equation} \label{inf1} \left\{ \begin{array}{lll} \dot{A_1} & = & -k\,A_1(t) - k_{12}A_1(t) + k_{21}A_2(t) + u_1(t)) \\ \dot{A_2} & = & k_{12}A_1(t) - k_{21}A_2(t) \end{array} \right. \end{equation} where the input $u_1$ in the central compartment is defined as $u_1(t) = a \, e^{-b \, t}.$ There is no more information about the input in the dataset:

ID TIME DV 1 2 51 1 4 63 1 6 48 1 15 23 1 18 16 $\vdots$ $\vdots$ $\vdots$

Then, several solutions exist for coding the input. A first solution consists in coding directly the input function in the ODE system:

EQUATION

DDT_Ac = -k*Ac - k12*Ac + k21*Ap + a*exp(-b*T)

DDT_Ap = k12*Ac - k21*Ap

This is clearly the simplest solution, but the input function is not defined as such. Then, it would be impossible to encode the input in the MML.

The input function can be defined by an external forcing function (called ExpInput in the following example)

EQUATION

INPUT(CMT=1,RATE=ExpInput(a,b))

DDT_Ac = -k*Ac - k12*Ac + k21*Ap

DDT_Ap = k12*Ac - k21*Ap

or directly in the model

EQUATION

INPUT(CMT=1,RATE=inline('a*exp(-b*T)'))

DDT_Ac = -k*Ac - k12*Ac + k21*Ap

DDT_Ap = k12*Ac - k21*Ap

Multiple inputs

Different inputs in different components of the ODE system can easily be combined.

Example: We consider a multiple administration with two different oral administrations using two different formulations with different release profiles, one subcutaneous injection and one skin patch:

ID TIME AMT TINF DPT DV 1 0 10 3 SC . 1 2 . . . 21 1 4 . . . 63 1 6 . . . 48 1 8 4 . GUT1 . 1 12 15 24 SP . 1 15 . . . 72 1 16 10 . GUT2 . 1 18 . . . 39 $\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$ $\vdots$

Here, AMT is the amount, TINF the infusion duration for the subcutaneous injection and for the skin patch and DPT is the ``depot compartment: GUT1 and GUT2 are the two oral administrations, SC holds for subcutaneous and SP holds for skin patch.

Important Remark: The datafile only contains information about the administration of the drugs, not about the PK model. The link between the depot compartments and the components of the ODES is defined in the model:

EQUATION
---
INPUT(DPT='GUT1', CMT=1)
INPUT(DPT='GUT2', CMT=2)
INPUT(DPT='SP', CMT=4)
INPUT(DPT='SC', CMT=3)
DDT_Ad1 = -ka1*Ad1
DDT_Ad2 = -ka2*Ad2
DDT_Ac =  ka1*Ad1 + ka2*Ad2 - k*Ac
DDT_As = -ks*As
  .         .
  .         .
  .         .

Remark: Bioavailability or lag-times can easily be taken into account in the model:

EQUATION
---
INPUT(DPT='GUT1', CMT=1 , Tlag= TLAG1 , P=F1)
INPUT(DPT='GUT2', CMT=2 , P=F2)
  .         .
  .         .
  .         .