Dynamical systems driven by ODEs
Contents
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: |
|
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} & = & F_1(A(t),u_1(t)) \\ \dot{A_2} & = & F_2(A(t),u_2(t)) \\ \vdots & \vdots & \vdots \\ \dot{1_L} & = & 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$, $ \dot{A_{\ell}} = 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 \begin{equation}\label{inf1} 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}
\label{inf1}
\left\{
\begin{array}{lll}
\dA{A}{1} & = & -k\,A_1(t) - k_{12}A_1(t) + k_{21}A_2(t) + u_1(t)) \\
\dA{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: |
|
Here, \verb"AMT" is the amount and \verb"TINF" the infusion duration. Alternatively, the infusion rate \verb"AMT/TINF" could be provided in a column \verb"RATE".
This is a very basic situation with only one type of administration. Then, there is no need of any additional column \verb"DPT" (depot compartment) in the datafile to distinguish different target depots for different types of administration, and the model reduces to:
\hspace*{2cm} \begin{minipage}[b]{10cm} \begin{verbatim} '"`UNIQ-MathJax70-QINU`"'EQUATION DDT_Ac = -k*Ac - k12*Ac + k21*Ap DDT_Ap = k12*Ac - k21*Ap \end{verbatim} \end{minipage}
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:
\hspace*{2cm} \begin{minipage}[b]{10cm} \begin{verbatim} $EQUATION INPUT(CMT=2) DDT_Ap = k12*Ac - k21*Ap DDT_Ac = -k*Ac - k12*Ac + k21*Ap \end{verbatim} \end{minipage} \subsubsection{Spike inputs} We consider an input-less dynamical system: for any $1\leq \ell \leq L$, '"`UNIQ-MathJax2-QINU`"' 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}$. \noindent \underline{Example}: Consider an IV bolus. \begin{minipage}{6cm} The input is given by the column \verb"AMT" in the dataset: \end{minipage} \hspace*{1cm} \begin{minipage}{8cm} \texttt{ \begin{tabular}{cccc} 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 \end{tabular} } \end{minipage} The model is exactly the same model defined above for an infusion, but without any column <span style="font-family:'Tahorma'">RATE</span> or <span style="font-family:'Tahorma'">TINF</span> in the datafile, spike inputs are assumed. ==='"`UNIQ--h-5--QINU`"'''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 <u>Example:</u> 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}.$ {| align=center; style="width: 700px" cellpadding="20" cellspacing="0" style="border: 1px solid darkgray;" | There is no more information about the input in the dataset: || {| class="wikitable" style="text-align:center; width:80%;" ! 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 <span style="font-family:'Tahorma'">ExpInput</span> 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 \verb"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) . . . . . .
