Difference between revisions of "Dynamical systems driven by ODEs"
m (→A basic model) |
m |
||
Line 134: | Line 134: | ||
;Remark2: | ;Remark2: | ||
:If the initial condition is not given in the model, it is assumed to be 0. | :If the initial condition is not given in the model, it is assumed to be 0. | ||
+ | |||
+ | |||
===''Piecewise defined dynamical systems''=== | ===''Piecewise defined dynamical systems''=== | ||
Line 195: | Line 197: | ||
We will consider the same viral kinetics model defined above. This system is now piecewise defined: | 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 [[#eq1|(1.4)]] | + | |
+ | {| cellpadding="5" cellspan="5" | ||
+ | | style="width: 450px;" | | ||
+ | * before $T_{Start1}$, $A(t) = A^{\star}$, where $A^{\star}$ is the equilibrium state defined in [[#eq1|(1.4)]] <br> | ||
* between $T_{Start1}$ and $T_{Start2}$, | * between $T_{Start1}$ and $T_{Start2}$, | ||
− | + | :<div id="vk3" style="text-align: left;font-size: 11pt"><math> \left\{ \begin{array}{lll} | |
− | <div id="vk3" | ||
− | |||
− | |||
− | |||
\dot{N}(t) & = & s - \beta \, N(t) \, V(t) - d\, N(t) \\ | \dot{N}(t) & = & s - \beta \, N(t) \, V(t) - d\, N(t) \\ | ||
\dot{I}(t) & = & \beta \, N(t) \, V(t) - \delta \, I(t) \\ | \dot{I}(t) & = & \beta \, N(t) \, V(t) - \delta \, I(t) \\ | ||
\dot{V}(t) & = & p(1-\varepsilon) \, I(t) - c \, V(t) \\ | \dot{V}(t) & = & p(1-\varepsilon) \, I(t) - c \, V(t) \\ | ||
− | \end{array} \right. </math> | + | \end{array} \right. </math></div> |
− | |||
− | |||
− | |||
− | </div> | ||
* between $T_{Start2}$ and $T_{Stop}$, the system is governed by the ODES described in [[#vk2|(1.5)]] | * between $T_{Start2}$ and $T_{Stop}$, the system is governed by the ODES described in [[#vk2|(1.5)]] | ||
− | <div id="vk4" | + | :<div id="vk4" style="text-align: left;font-size: 11pt"><math> |
− | + | \left\{ \begin{array}{lll} | |
− | |||
− | |||
\dot{N}(t) & = & s - \beta(1-\eta) \, N(t) \, V(t) - d\,N(t) \\ | \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{I}(t) & = & \beta(1-\eta) \, N(t) \, V(t) - \delta \, I(t) \\ | ||
\dot{V}(t) & = & p(1-\varepsilon) \, I(t) - c \, V(t)\\ | \dot{V}(t) & = & p(1-\varepsilon) \, I(t) - c \, V(t)\\ | ||
− | \end{array} \right. </math> | + | \end{array} \right. </math></div> |
− | |||
− | |||
− | |||
− | </div> | ||
* after $T_{Stop}$, the system smoothly returns to its original state governed by the ODES described in [[#vk1|(1.3)]] | * after $T_{Stop}$, the system smoothly returns to its original state governed by the ODES described in [[#vk1|(1.3)]] | ||
− | + | :<div id="Equation8" style="text-align: left;font-size: 11pt"><math>\left\{ \begin{array}{lll} | |
− | |||
− | <div id="Equation8" | ||
− | |||
− | |||
− | |||
\dot{N}(t) & = & s - \beta (1-\eta \, e^{-k_1 (t-T_{Stop})})\, N(t) \, V(t) - d\,N(t) \\ | \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{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) \\ | \dot{V}(t) & = & p(1-\varepsilon\, e^{-k_2 (t-T_{Stop})}) \, I(t) - c \, V(t) \\ | ||
− | \end{array} \right. </math></ | + | \end{array} \right. </math></div> |
− | || | + | |
+ | |style="width: 150px; font-size: 11pt"| (1.6) <br><br><br><br><br> (1.7) <br><br><br><br><br> (1.8) | ||
+ | |style="width: 600px; " | | ||
+ | <pre style=" background-color:#EFEFEF; font-family:'courier new';font-size:10pt;border: 1px solid darkgray; border-radius:1em;"> | ||
+ | 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 | ||
+ | |||
+ | |} | ||
+ | |||
+ | |||
+ | |||
+ | |||
+ | |||
|style="text-align: right; margin-left:20em; font-size: 11pt" | (1.8) | |style="text-align: right; margin-left:20em; font-size: 11pt" | (1.8) | ||
|} | |} | ||
Line 246: | Line 263: | ||
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. | ||
− | + | As you can see, the three different dynamical systems are des | |
+ | cribed in the MDL (in a block <span style="font-family:'courier new';font-size:13pt">EQUATION</span> and using the statement <span style="font-family:'courier new';font-size:13pt">SWITCH</span> with MLXTRAN). We only show the blocks <span style="font-family:'courier new';font-size:13pt">VARIABLES</span> and <span style="font-family:'courier new';font-size:13pt">EQUATION</span> of the code: | ||
{| align=left; style="width: 700px; background-color:#EFEFEF; font-family:'courier new';font-size:12pt;" | {| align=left; style="width: 700px; background-color:#EFEFEF; font-family:'courier new';font-size:12pt;" | ||
− | | | + | | |
− | |||
<code><div style="font-family:'courier new';font-size:12pt;"> $EQUATION </div></code> | <code><div style="font-family:'courier new';font-size:12pt;"> $EQUATION </div></code> |
Revision as of 12:45, 13 February 2013
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)
|
(1.1) |
where $\dot{A}(t)$ denotes the vector of derivatives of $A(t)$ with respect to $t$:
|
(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 there is no 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.
|
|
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:
|
(1.3) |
The equilibrium state of this system is $A^{\star} = (N^{\star} , I^{\star} , V^{\star})$, where
|
(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)
|
(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: viral kinetic model.
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$: |
|
We will consider the same viral kinetics model defined above. This system is now piecewise defined:
|
(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 |} |style="text-align: right; margin-left:20em; font-size: 11pt" | (1.8) |} </div> We have seen that the information about switching times is given in the data set. As you can see, the three different dynamical systems are des cribed in the MDL (in a block <span style="font-family:'courier new';font-size:13pt">EQUATION</span> and using the statement <span style="font-family:'courier new';font-size:13pt">SWITCH</span> with MLXTRAN). We only show the blocks <span style="font-family:'courier new';font-size:13pt">VARIABLES</span> and <span style="font-family:'courier new';font-size:13pt">EQUATION</span> of the code: {| align=left; style="width: 700px; background-color:#EFEFEF; font-family:'courier new';font-size:12pt;" | <code><div style="font-family:'courier new';font-size:12pt;"> $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:
INPUT(CMT=2)
DDT_Ap = k12*Ac - k21*Ap DDT_Ac = -k*Ac - k12*Ac + k21*Ap |
Spike inputs
We consider an input-less dynamical system: for any $1\leq \ell \leq L$, ::
Spike inputs means that there exists a sequence of times $(\tau_{\ell,j})$ and amounts $(a_{\ell,j})$ such that ::
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. 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.
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}.$
There is no more information about the input in the dataset:
|
|
Then, several solutions exist for coding the input. A first solution consists in coding directly the input function in the ODE system:
|
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
|
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:
|
INPUT(DPT='GUT1', CMT=1 , Tlag= TLAG1 , P=F1) |
INPUT(DPT='GUT2', CMT=2 , P=F2) |
|