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)

\( \dot{A} = F(A(t)) \)

(1.1)

where $\dot{A}(t)$ denotes the vector of derivatives of $A(t)$ with respect to $t$:

\( \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. \)

(1.2)

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 the system doesn't have any input:

\(\begin{align} A(t_0) & = &A_0 \\ \dot{A}(t) & = &F(A(t)) \, \ t \geq t_0 \end{align}\)

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:

\( \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. \)

(1.3)

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

\( N^{\star} = \displaystyle{\frac{\delta \, c}{ \beta \, p}} \quad ; \quad I^{\star} = \displaystyle{\frac{s - d\, N^{\star}}{ \delta}} \quad ; \quad V^{\star} = \displaystyle{\frac{ p \, I^{\star} }{c}}. \)

(1.4)

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 the new ODE system (1.5). This system, and the initial conditions, can be exactly reproduced using the MLXTRAN macros (block EQUATION)

\( \left\{ \begin{array} \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. \) where $0<\varepsilon <1$ and $0 < \eta < 1$.

(1.5)

$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{align} A(t_0) &= A_0 \\ \dot{A}(t) &= F_k(A(t)) \ , \ t_{k-1} \leq t \leq t_{k} \end{align}\)

Example: A 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 (the switching times is still given in the data set). This system is now piecewise defined:

before $T_{Start1}$, $A(t) = A^{\star}$, where $A^{\star}$ is the equilibrium state defined in (1.4) between $T_{Start1}$ and $T_{Start2}$, \( \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. \) between $T_{Start2}$ and $T_{Stop}$, the system is governed by the ODES described in (1.5) \( \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. \) after $T_{Stop}$, the system smoothly returns to its original state governed by the ODES described in (1.3) \(\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. \)

(1.6) (1.7) (1.8)

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

As you can see, the three 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.

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:

  
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 (1.1) with inputs:

\( \dot{A}(t) = F(A(t),u(t)) \)

(2.1)

where:

\( \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. \)

(2.2)

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$

\( \dot{A_{\ell}}(t) = F_\ell(A(t)) + u_\ell(t) \)

Then, there exists a sequence of times $(\tau_{\ell,j})$, durations $(d_{\ell,j})$ and amounts $(a_{\ell,j})$ such that

\( u_\ell(t) = \left\{\begin{array}{ll} \displaystyle{\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. \)

(2.3)

Example:

Consider an IV infusion with two compartments, i.e. a very basic situation with only one type of administration. Such infusion can be described by the following system of equations:

\( \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. \) The input $u_1$ is the rate of infusion and some of its possible values are indicated 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$ AMT is the infusion amount and TINF the infusion duration. Alternatively, the infusion rate AMT/TINF could be provided in a column RATE.

If there is no need of any additional column DPT (depot compartment) in the datafile to distinguish different target depots for different types of administration, 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 <br> DDT_Ac = -k*Ac - k12*Ac + k21*Ap

Spike inputs

We consider an input-less dynamical system: for any $1\leq \ell \leq L$,

\( \dot{A_\ell}(t) = F_\ell(A(t)) \)

Spike inputs means that there exists a sequence of times $(\tau_{\ell,j})$ and amounts $(a_{\ell,j})$ such that

\( A_\ell(\tau_{\ell,j}) = A_\ell(\tau_{\ell,j}^{-}) + a_{\ell,j} \)

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. In this case, the input is given by the column AMT in the dataset. 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.

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$

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 \( \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. \) where the input $u_1$ in the central compartment is defined as $u_1(t) = a \, e^{-b \, t}$ and 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: 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. 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$

$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 ... ... $EQUATION INPUT(DPT='GUT1', CMT=1 , Tlag= TLAG1 , P=F1) INPUT(DPT='GUT2', CMT=2 , P=F2) ...

Remark 1: 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:

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