Difference between revisions of "Dynamical systems driven by ODEs"

From Popix
Jump to navigation Jump to search
Line 108: Line 108:
  
  
{| align=left; style="width: 800px"font-style:italic" cellspacing="0" style="border: 1px solid darkgray;"  
+
{| align=left; style="width: 800px"font-style:italic" cellspacing="0" style="border: 1px solid darkgray;"  
  ! EQUATIONS
+
  ! EQUATION
 
  |-
 
  |-
 
  |T_0 = 0
 
  |T_0 = 0
Line 144: Line 144:
 
\end{eqnarray*}
 
\end{eqnarray*}
  
{| align=left; align=center; style="border: 1px solid darkgray;"
+
{| align=left; style="border: 1px solid darkgray;"
  <u>Example:</u>  viral kinetic model. We assume here that a first treatment which blocks the production of virus starts first at  
+
  | <u>Example:</u>  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  
  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}$.
 
  treatments stop at time $T_{Stop}$.
 
|}
 
|}
Line 228: Line 227:
 
We have seen that the information about switching times is given in the data set.
 
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 <span style="font-family:'Monospaced'>$EQUATION<span/> and using the statement <span style="font-family:'Monospaced'>SWITCH<span/> with MLXTRAN).
 
The different dynamical systems are described in the MDL (in a block <span style="font-family:'Monospaced'>$EQUATION<span/> and using the statement <span style="font-family:'Monospaced'>SWITCH<span/> with MLXTRAN).
We only show the blocks <span style="font-family:'Monospaced'>$VARIABLES<span/> and <span style="font-family:'Monospaced'>$EQUATION<span/> of the code:
+
We only show the blocks <span style="font-family:'Monospaced'>VARIABLES<span/> and <span style="font-family:'Monospaced'>EQUATION<span/> of the code:
  
 
{| align=left style="width: 700px" cellspacing="0" style="border: 1px solid darkgray;"
 
{| align=left style="width: 700px" cellspacing="0" style="border: 1px solid darkgray;"
Line 250: Line 249:
 
   |align="right" width="50%"| CASE T_Start1 < T < T_Start2
 
   |align="right" width="50%"| CASE T_Start1 < T < T_Start2
 
     |-
 
     |-
     |DDT_N = s - beta*N*V - d*N
+
     |align="right" width="40%"|DDT_N = s - beta*N*V - d*N
 
     |-
 
     |-
     |DDT_I = beta*N*V - delta*I
+
     |align="right" width="40%"|DDT_I = beta*N*V - delta*I
 
     |-
 
     |-
     |DDT_V = p*(1-epsilon)*I - c*V
+
     |align="right" width="40%"|DDT_V = p*(1-epsilon)*I - c*V
 
     |-
 
     |-
 
     |-
 
     |-
  CASE T_Start2 < T < T_Stop
+
    |align="right" width="40%"|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
+
    |align="right" width="40%"|DDT_N = s - beta*(1-eta)*N*V - d*N
     DDT_V = p*(1-epsilon)*I - c*V
+
     |-
 +
    |align="right" width="40%"|DDT_I = beta*(1-eta)*N*V - delta*I
 +
     |-
 +
    |align="right" width="40%"|DDT_V = p*(1-epsilon)*I - c*V
  
 
   CASE T > T_Stop
 
   CASE T > T_Stop

Revision as of 16:07, 23 January 2013

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<span></span> and using the statement <span style="font-family:'Monospaced'">SWITCH<span></span> with MLXTRAN). We only show the blocks <span style="font-family:'Monospaced'">VARIABLES<span></span> and <span style="font-family:'Monospaced'">EQUATION<span></span> of the code: {| align=left style="width: 700px" cellspacing="0" style="border: 1px solid darkgray;" | $VARIABLES ID, TIME, VL use=DV, EVENT list=(Start1, Start2, Stop) |- |- | $EQUATION |- |- | SWITCH |- |align="right" width="50%"| CASE T < T_Start1 |- | align="right" width="40%"| N = delta*c/(beta*p); |- | align="right" width="40%"| I = (s-d*N)/delta |- | align="right" width="40%"|V = p*I/c |- |- |align="right" width="50%"| CASE T_Start1 < T < T_Start2 |- |align="right" width="40%"|DDT_N = s - beta*N*V - d*N |- |align="right" width="40%"|DDT_I = beta*N*V - delta*I |- |align="right" width="40%"|DDT_V = p*(1-epsilon)*I - c*V |- |- |align="right" width="40%"|CASE T_Start2 < T < T_Stop |- |align="right" width="40%"|DDT_N = s - beta*(1-eta)*N*V - d*N |- |align="right" width="40%"|DDT_I = beta*(1-eta)*N*V - delta*I |- |align="right" width="40%"|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, <span style="font-family:'Monospaced'">EVENT<span></span> is a reserved variable name. Then, the information in the column <span style="font-family:'Monospaced'">EVENT<span></span> is recognized as a succession of events. Furthermore, the times of the events \verb"Start1", \verb"Start2" and <span style="font-family:'Monospaced'">Stop<span></span> are automatically created as <span style="font-family:'Monospaced'">T_Start1<span></span>, <span style="font-family:'Monospaced'">T_Start2<span></span> and <span style="font-family:'Monospaced'">T_Stop<span></span>. * ''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: {| align=left style="width: 700px" cellspacing="0" style="border: 1px solid darkgray;" ! $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