Difference between revisions of "Dynamical systems driven by ODEs"

From Popix
Jump to navigation Jump to search
Line 109: Line 109:
  
 
{| 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;"  
  ! EQUATION
+
  | $EQUATION
 
  |-
 
  |-
 
  |T_0 = 0
 
  |T_0 = 0
Line 232: Line 232:
  
 
  $VARIABLES  ID, TIME, VL use=DV, EVENT list=(Start1, Start2, Stop)
 
  $VARIABLES  ID, TIME, VL use=DV, EVENT list=(Start1, Start2, Stop)
 
 
  $EQUATION
 
  $EQUATION
 
 
  SWITCH
 
  SWITCH
 
     CASE T < T_Start1
 
     CASE T < T_Start1
Line 253: Line 251:
 
         DDT_I = beta*(1-eta*exp(-k1*(T-T_Stop)))*N*V - delta*I
 
         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
 
         DDT_V = p*(1-epsilon*exp(-k2*(T-T_Stop)))*I - c*V
 
 
  END
 
  END
 
 
 
 
 
 
 
 
 
  
  
Line 270: Line 258:
 
* ''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:
 
* ''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
! $EQUATION
+
T_0 = T_Start1
|-
+
N_0 = delta*c/(beta*p);
|T_0 = T_Start1
+
I_0 = (s-d*N)/delta
|-
+
V_0 = p*I/c
|N_0 = delta*c/(beta*p);
+
 
|-
+
SWITCH
|I_0 = (s-d*N)/delta
+
  CASE T_Start1 < T < T_Start2
|-
+
    be = beta
|V_0 = p*I/c
+
    pe = p*(1-epsilon)
  
SWITCH
+
  CASE T_Start2 < T < T_Stop
  CASE T_Start1 < T < T_Start2
+
    be = beta*(1-eta)
    be = beta
+
    pe = p*(1-epsilon)
    pe = p*(1-epsilon)
 
  
  CASE T_Start2 < T < T_Stop
+
  CASE T > T_Stop
    be = beta*(1-eta)
+
    be = beta*(1-eta*exp(-k1*(T-T_Stop)))
    pe = p*(1-epsilon)
+
    pe = p*(1-epsilon*exp(-k2*(T-T_Stop)))
 +
END
  
  CASE T > T_Stop
+
DDT_N = s - be*N*V - d*N
    be = beta*(1-eta*exp(-k1*(T-T_Stop)))
+
DDT_I = be*N*V - delta*I
    pe = p*(1-epsilon*exp(-k2*(T-T_Stop)))
+
DDT_V = pe*I - c*V
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'''==
 
=='''Dynamical systems with source terms'''==

Revision as of 16:28, 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, <span style="font-family:'Monospaced';">$T_0 = 0$</span> means that the system is constant and is $A^{\star}$, defined in the script by <span style="font-family:'Monospaced';">$(N_0, I_0, V_0)$</span>, for any $t<0$. * ''Remark2:'' If the initial condition is not given in the model, it is assumed to be 0. ==='"`UNIQ--h-2--QINU`"'''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*} {| 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 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 <span style="font-family:'Monospaced'">EVENT</span> is necessary in the dataset to describe this information EVENT is an extension of the <span style="font-family:'Monospaced'">EVID</span> (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$: {| class="wikitable" style="text-align:center; width:400px;" ! 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 <span style="font-family:'Monospaced'">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: $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, <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:

 $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