Dynamical systems driven by ODEs

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} \dA{A}{} = F(A(t)) \end{equation} where $\dA{A}{}$ denotes the vector of derivatives of $A(t)$ with respect to $t$: \begin{equation} \label{ode2_model} \left\{ \begin{array}{lll} \dA{A}{1} & = & F_1(A(t)) \\ \dA{A}{2} & = & F_2(A(t)) \\ \vdots & \vdots & \vdots \\ \dA{A}{L} & = & F_L(A(t)) \end{array} \right. \end{equation}

Notations: \begin{itemize} \item let '"`UNIQ-MathJax7-QINU`"' be the initial condition of the system defined at the initial time '"`UNIQ-MathJax8-QINU`"', \item let '"`UNIQ-MathJax9-QINU`"' be the solution of the system at equilibrium: '"`UNIQ-MathJax10-QINU`"' \end{itemize}

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

\noindent \underline{Example}:  A viral kinetic (VK) model.

\begin{minipage}{8cm}
 In this example, the data file contains the viral load:
\end{minipage}

\hspace*{1cm} \begin{minipage}{8cm} \texttt{ \begin{tabular}{ccc} 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 \\ \end{tabular} } \end{minipage}

Consider a basic VK model with $A=(\nitc,\itc,\vl)$ where $\nitc$ is the number of non infected target cells, $\itc$ the number of infected target cells and $\vl$ 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} \dA{\nitc}{} & = & s - \beta \, \nitc(t) \, \vl(t) - d\,\nitc(t) \\ \dA{\itc}{} & = & \beta \, \nitc(t) \, \vl(t) - \delta \, \itc(t) \\ \dA{\vl}{} & = & p \, \itc(t) - c \, \vl(t) \end{array} \right. \end{equation} The equilibrium state of this system is $\As = (\nitc^\star , \itc^\star , \vl^\star)$. where \begin{equation} \label{eq1} \nitc^\star = \frac{\delta \, c}{ \beta \, p} \quad ; \quad \itc^\star = \frac{s - d\,\nitc^\star}{ \delta} \quad ; \quad \vl^\star = \frac{ p \, \itc^\star }{c}. \end{equation}

Assume that the system has reached the equilibrium state $\As$ 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}
  \dA{\nitc}{} & = & s - \beta(1-\eta) \, \nitc(t) \, \vl(t) - d\,\nitc(t)  \\
  \dA{\itc}{} & = & \beta(1-\eta) \, \nitc(t) \, \vl(t) - \delta \, \itc(t) \\
  \dA{\vl}{} & = & p(1-\varepsilon) \, \itc(t) - c \, \vl(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 (in a block \verb"$EQUATION" with MLXTRAN): \hspace*{4cm} \begin{minipage}[b]{10cm} \begin{verbatim} $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 \end{verbatim} \end{minipage}

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

\noindent{\bf Remark2:} If the initial condition is not given in the model, it is assumed to be 0.

\subsubsection{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 \\ \dA{A}{} &= &F_k(A(t)) \ , \ t_{k-1} \leq t \leq t_{k} \end{eqnarray*}

\noindent \underline{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 \texttt{EVENT} is necessary in the dataset to describe this information (\texttt{EVENT} is an extension of the \texttt{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$:

\hspace*{4cm} \texttt{ \begin{tabular}{cccc} 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 \end{tabular} }

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

\begin{itemize}
 \item  before '"`UNIQ-MathJax32-QINU`"', '"`UNIQ-MathJax33-QINU`"', where '"`UNIQ-MathJax34-QINU`"' is the equilibrium state defined in (\ref{eq1})
\item between '"`UNIQ-MathJax35-QINU`"' and '"`UNIQ-MathJax36-QINU`"',
\begin{equation}
\label{vk3}
\left\{
\begin{array}{lll}
  \dA{\nitc}{} & = & s - \beta \, \nitc(t) \, \vl(t) - d\,\nitc(t)  \\
  \dA{\itc}{} & = & \beta \, \nitc(t) \, \vl(t) - \delta \, \itc(t) \\
  \dA{\vl}{} & = & p(1-\varepsilon) \, \itc(t) - c \, \vl(t)
\end{array}

\right. \end{equation} \item between $T_{Start2}$ and $T_{Stop}$, the system is governed by the ODES described in (\ref{vk2})

$$
\left\{
\begin{array}{lll}
  \dA{\nitc}{} & = & s - \beta(1-\eta) \, \nitc(t) \, \vl(t) - d\,\nitc(t)  \\
  \dA{\itc}{} & = & \beta(1-\eta) \, \nitc(t) \, \vl(t) - \delta \, \itc(t) \\
  \dA{\vl}{} & = & p(1-\varepsilon) \, \itc(t) - c \, \vl(t)
\end{array}
\right.
$$

\item 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} \dA{\nitc}{} & = & s - \beta (1-\eta \, e^{-k_1 (t-T_{Stop})})\, \nitc(t) \, \vl(t) - d\,\nitc(t) \\ \dA{\itc}{} & = & \beta (1-\eta \, e^{-k_1 (t-T_{Stop})}) \, \nitc(t) \, \vl(t) - \delta \, \itc(t) \\ \dA{\vl}{} & = & p(1-\varepsilon\, e^{-k_2 (t-T_{Stop})}) \, \itc(t) - c \, \vl(t) \end{array} \right. \end{equation} \end{itemize} 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 \verb"$EQUATION" and using the statement \verb"SWITCH" with MLXTRAN). We only show the blocks \verb"$VARIABLES" and \verb"$EQUATION" of the code: \hspace*{2cm} \begin{minipage}[b]{10cm} \begin{verbatim} $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 DDT_V = p*(1-epsilon)*I - c*V 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 \end{verbatim} \end{minipage} \noindent{\bf Remark 1:} Here, \verb"EVENT" is a reserved variable name. Then, the information in the column \verb"EVENT" is recognized as a succession of events. Furthermore, the times of the events \verb"Start1", \verb"Start2" and \verb"Stop" are automatically created as \verb"T_Start1", \verb"T_Start2" and \verb"T_Stop". \noindent{\bf 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: \hspace*{2cm} \begin{minipage}[b]{10cm} \begin{verbatim} $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 \end{verbatim} \end{minipage}