Równanie różniczkowe - wchłanianie lekarstwa

[Marakech]

Witam

Dostałem zadanie: "Należy opracować program komputerowy do rozwiązania układu równań różniczkowych - Model matematyczny wchłaniania lekarstwa przez organizm człowieka". Program ten powinien być wykonany w Derive. Potrzebuje wskazówek jak się do tego zabrać, co należy zrobić ponieważ nie miałem styczności z tym programem.

Treść: Zakładając, że mechanizm wchłaniania lekarstwa przez organizm człowieka odbywa się dwuetapowo. Najpierw po spożyciu, lekarstwo dostaje się do układu żołądkowo-jelitowego a następnie do układu krwionośnego. Równania opisujące ten proces są następujące:

\(\displaystyle{ \frac{dm_{1}}{dt} =−k_{1} m_{1} +y(t),}\)
\(\displaystyle{ \frac{dm _{2} }{dt} =k_{1} m_{1} −k_{2} m_{2} ,}\)
\(\displaystyle{ \frac{dy}{dt}=−ay,}\)

\(\displaystyle{ m_1}\) – masa lekarstwa w układzie żołądkowo - jelitowego,
\(\displaystyle{ m_2}\) – masa lekarstwa we krwi,
\(\displaystyle{ y(t)}\) - szybkość przenikania lekarstwa do układu żołądkowo – jelitowego,
\(\displaystyle{ k_1}\) – współczynnik proporcjonalności charakteryzujący stopień przenikania lekarstwa z układu pokarmowego do krwi,
\(\displaystyle{ k_2}\) – współczynnik proporcjonalności charakteryzujący przemianę materii,
\(\displaystyle{ a}\) – współczynnik proporcjonalności charakteryzujący stopień wchłaniania lekarstwa do układu trawiennego.

Dane do programu to:
- współczynniki \(\displaystyle{ k_1, k_2}\),
- warunki początkowe: \(\displaystyle{ m_1(0), m_2(0), y(0)}\).

[janusz47]

Dzień Dobry

Ostatnia wersja programu Derive - Derive 6 nie uruchamia się w Windows 10. Program nie wytrzymał próby czasu. Pamiętam , że w latach dziewięćdziesiątych ubiegłego stulecia był programem atrakcyjnym.

Jeśli zmieni Pan rozwiązanie tego zadania na program np. Maple, to myślę, że będę mógł w stanie pomóc.

[Marakech]

Dzień dobry,
Derive działa pod Windows 10, wystarczy zainstalować brakujące składniki, fonty. Fakt jest już leciwy dlatego chętnie zobaczę, jak wykonać to zadanie w Maple bądź innym programie.

[janusz47]

Dzień Dobry

Czy warunki początkowe są zerowe ? Czy znane są wartości współczynników \(\displaystyle{ k_{1}, k_{2} ? }\)

[Marakech]

Warunki początkowe są zerowe, współczynniki \(\displaystyle{ k _{1}, k _{2} }\) nie są narzucone.

[janusz47]

\(\displaystyle{ \begin{cases} m'_{1}(t) = k_{1}\cdot m_{1}(t) + y(t) \ \ (1) \\ m'_{2}(t) = k_{1}\cdot m_{1}(t) - k_{2}\cdot m_{2}(t) \ \ (2) \\ y'(t) = - a\cdot y(t) \ \ (3) \end{cases} }\)

Równanie \(\displaystyle{ (3) }\) :

\(\displaystyle{ y'(t) = -a\cdot y(t) }\)

całkujemy obustronnie:

\(\displaystyle{ \int_{y(0)}^{y(t)} \frac{y'(\tau)}{y(\tau)}dy = -a\int_{t_{0}}^{t} d\tau }\)

\(\displaystyle{ \ln(\tau) \mid_{y(0)}^{y(t)} = -a\cdot \tau \mid_{t_{0}}^{t} }\)

\(\displaystyle{ \ln(y(t)) - \ln(y(0)) = -a\cdot (t - t_{0}) }\)

\(\displaystyle{ \ln\left(\frac{y(t)}{y(0)}\right) = -a\cdot (t- t_{0})) }\)

\(\displaystyle{ y(t) = y(0)\cdot e^{-a\cdot (t - t_{0})} \ \ (4) }\)

Rozwiązanie równania różniczkowego \(\displaystyle{ (3) }\) pokazuje, że szybkość przenikania leku do układu żołądkowo-jelitowego opisana jest malejącą funkcją wykładniczą.

Postać macierzowa układu złożonego z równań: \(\displaystyle{ (1), \ \ (2): }\)

\(\displaystyle{ \left [\begin{matrix} m'_{1}(t) \\ m'_{2}(t) \end{matrix} \right] = \left [\begin{matrix} k_{1} & 0 \\ k_{1} & - k_{2} \end{matrix} \right] \cdot \left [\begin{matrix} m_{1}(t) \\ m_{2}(t) \end{matrix} \right] + \left [\begin{matrix} y(t) \\ 0 \end{matrix} \right] }\)

Niech \(\displaystyle{ A = \left [\begin{matrix} k_{1} & 0 \\ k_{1} & -k_{2} \end{matrix} \right],}\).

rozwiązaniem ogólnym równania jednorodnego

\(\displaystyle{ \left [\begin{matrix} m'_{1}(t) \\ m'_{2}(t) \end{matrix} \right] = A \cdot \left [\begin{matrix} m_{1}(t) \\ m_{2}(t) \end{matrix} \right] }\)

jest trajektoria:

\(\displaystyle{ \left [\begin{matrix} m_{1}(t) \\ m_{2}(t) \end{matrix} \right] = e^{A\cdot t} \left [\begin{matrix} m_{1}(0) \\ m_{2}(0) \end{matrix} \right] }\)

Diagonalizujemy macierz \(\displaystyle{ A: }\)

Wartości własne macierzy :

\(\displaystyle{ \det( A - \lambda I) = \det \left [\begin{matrix} k_{1} -\lambda & 0 \\ k_{1} & - k_{2} -\lambda \end{matrix} \right] = 0 }\)

\(\displaystyle{ w(\lambda) = (k_{1} -\lambda)( -k_{2} -\lambda ) = 0 }\)

\(\displaystyle{ \lambda_{1} = k_{1}, \ \ \lambda_{2} = -k_{2}. }\)

Wektory własne odpowiadające wartością własnym:

\(\displaystyle{ ker(A - k_{1} \cdot I) = ker \left [\begin{matrix} k_{1} - k_{1} & 0 \\ k_{1} & - k_{2} - k_{1} \end{matrix} \right] = lin \left (\left [\begin{matrix} 1 \\ \frac{k_{1}}{k_{1}+k_{2}} \end{matrix} \right] \right ) }\)

\(\displaystyle{ ker(A - k_{2} \cdot I) = ker \left [\begin{matrix} k_{1} - k_{2} & 0 \\ k_{1} & -k_{2} + k_{2} \end{matrix} \right] = lin \left (\left [\begin{matrix} 0\\ 1 \end{matrix} \right] \right ) }\)

Bazą diagonalizującą jest układ wektorów:

\(\displaystyle{ \left\{\left [\begin{matrix} 1 \\ \frac{k_{1}}{k_{1}+k_{2}} \end{matrix} \right] , \left [ \begin{matrix} 0 \\ 1 \end{matrix} \right] \right\} }\)

\(\displaystyle{ \left[\begin{matrix} 1 & 0 \\ \frac{k_{1}}{k_{1}+k_{2}} & 1\end{matrix} \right] ^{-1} =\left[\begin{matrix} 1 & 0 \\ -\frac{k_{1}}{k_{1}+k_{2}} & 1\end{matrix} \right] }\)

Stąd

\(\displaystyle{ e^{A\cdot t} = \left[\begin{matrix} 1 & 0 \\ -\frac{k_{1}}{k_{1}+ k_{2}} & 1 \end{matrix} \right]\cdot \left[\begin{matrix} e^{k_{1}\cdot t} & 0 \\ 0 & e^{-k_{2}\cdot t} \end{matrix} \right] \cdot \left[\begin{matrix} 1 & 0 \\ \frac{k_{1}}{k_{1}+k_{2}} & 1 \end{matrix} \right] =\left[ \begin{matrix} e^{k_{1}\cdot t} & 0 \\ -\frac{k_{1}}{k_{1}+k_{2}} e^{k_{1}\cdot t} + \frac{k_{1}}{k_{1}+k_{2}} e^{-k_{2}\cdot t} & e^{-k_{2}\cdot t} \end{matrix} \right] }\).

Otrzymaliśmy rozwiązanie ogólne układu jednorodnego.

\(\displaystyle{ \left [\begin{matrix} m_{1}(t) \\ m_{2}(t) \end{matrix} \right] = \left[ \begin{matrix} e^{k_{1}\cdot t} & 0 \\ -\frac{k_{1}}{k_{1}+k_{2}} e^{k_{1}\cdot t} + \frac{k_{1}}{k_{1}+k_{2}} e^{-k_{2}\cdot t} & e^{-k_{2}\cdot t} \end{matrix} \right] \cdot \left [\begin{matrix} m_{10} \\ m_{20} \end{matrix} \right] }\)

W celu znalezienia rozwiązania ogólnego układu niejednorodnego - zastosujemy metodę uzmiennienia stałych, przyjmujemy:

\(\displaystyle{ \left [\begin{matrix} m_{1}(t) \\ m_{2}(t) \end{matrix} \right] = \left[ \begin{matrix} e^{k_{1}\cdot t} & 0 \\ -\frac{k_{1}}{k_{1}+k_{2}} e^{k_{1}\cdot t} + \frac{k_{1}}{k_{1}+k_{2}} e^{-k_{2}\cdot t} & e^{-k_{2}\cdot t} \end{matrix} \right] \cdot \left [\begin{matrix} m_{10}(t) \\ m_{20}(t) \end{matrix} \right] }\)

Stąd

\(\displaystyle{ \left [\begin{matrix} m'_{1}(t) \\ m'_{2}(t) \end{matrix} \right] = \left[ \begin{matrix} k_{1}e^{k_{1}\cdot t} & 0 \\ -\frac{k_{1}}{k_{1}+k_{2}} (k_{1}e^{k_{1}\cdot t}+ k_{2}e^{-k_{2}\cdot t} ) & -k_{2} e^{-k_{2}\cdot t} \end{matrix} \right] \cdot \left [\begin{matrix} m'_{10}(t) \\ m'_{20}(t) \end{matrix} \right] }\)

Po wstawieniu do układu niejednorodnego:

\(\displaystyle{ \left [\begin{matrix} m'_{10}(t) \\ m'_{20}(t) \end{matrix} \right ] = \frac{1}{e^{(k_{1}-k_{2})\cdot t}}\cdot \left[ \begin{matrix} \frac{e^{-k_{2}\cdot t}}{k_{1}} & - \frac{(k_{1}\cdot e^{k_{1}\cdot t}+ k_{2}\cdot e^{-k_{2}\cdot t})}{k_{2}\cdot (k_{1}+k_{2})} \\ 0 & -\frac{e^{k_{1}\cdot t}}{k_{2}} \end{matrix} \right]\cdot \left [\begin{matrix} m'_{1}(t) \\ m'_{2}(t) \end{matrix} \right] + \left[\begin{matrix} y(t) \\ 0 \end{matrix}\right]}\)

\(\displaystyle{ \left [\begin{matrix} m_{10}(t) \\ m_{20}(t) \end{matrix} \right] = \int_{1}^{t} \left \{ \frac{1}{e^{(k_{1}-k_{2})\cdot \tau}}\cdot \left[ \begin{matrix} \frac{e^{-k_{2}\cdot \tau}}{k_{1}} & - \frac{(k_{1}\cdot e^{k_{1}\cdot \tau}+ k_{2}\cdot e^{-k_{2}\cdot \tau})}{k_{2}\cdot (k_{1}+k_{2})} \\ 0 & -\frac{e^{k_{1}\cdot \tau}}{k_{2}} \end{matrix} \right]\cdot \left[ \begin{matrix} k_{1}e^{k_{1}\cdot \tau} & 0 \\ -\frac{k_{1}}{k_{1}+k_{2}} (k_{1}e^{k_{1}\cdot \tau}+ k_{2}e^{-k_{2}\cdot \tau} ) & -k_{2} e^{-k_{2}\cdot \tau} \end{matrix} \right] \right\} d\tau + \int_{1}^{t} \left[ \begin{matrix} y(\tau) \\ 0 \end{matrix} \right]d\tau }\)


\(\displaystyle{ \left [\begin{matrix} m_{1}(t) \\ m_{2}(t) \end{matrix} \right] = \left[ \begin{matrix} e^{k_{1}\cdot t} & 0 \\ -\frac{k_{1}}{k_{1}+k_{2}} e^{k_{1}\cdot t} + \frac{k_{1}}{k_{1}+k_{2}} e^{-k_{2}\cdot t} & e^{-k_{2}\cdot t} \end{matrix} \right] \cdot \left\{ \int_{1}^{t} \left \{\frac{1}{e^{(k_{1}-k_{2})\cdot \tau}}\cdot \left[ \begin{matrix} \frac{e^{-k_{2}\cdot \tau}}{k_{1}} & - \frac{(k_{1}\cdot e^{k_{1}\cdot \tau}+ k_{2}\cdot e^{-k_{2}\cdot \tau})}{k_{2}\cdot (k_{1}+k_{2})} \\ 0 & -\frac{e^{k_{1}\cdot \tau}}{k_{2}} \end{matrix} \right]\cdot \left[ \begin{matrix} k_{1}e^{k_{1}\cdot \tau} & 0 \\ -\frac{k_{1}}{k_{1}+k_{2}} (k_{1}e^{k_{1}\cdot \tau}+ k_{2}e^{-k_{2}\cdot \tau} ) & -k_{2} e^{-k_{2}\cdot \tau} \end{matrix} \right] \right\}
d\tau + \int_{1}^{t} \left[ \begin{matrix} y(\tau) \\ 0 \end{matrix} \right]d\tau \right\} }\)

Funkcja \(\displaystyle{ y(t) }\) zostało określona w równaniu \(\displaystyle{ (4). }\)

Pozostało wymnożenie macierzy i obliczenie tych nietrudnych lecz żmudnych całek. Wzory uproszczą się.

Wyniki w programie Mathematica 9

\(\displaystyle{ DSolve[\{m'1[t] ==k1*m1[t] +y0*e^{(a*(t-t0))},m2'[t] == k1*m1[t]-k2*m2[t]\}, \{ m1[t], m2[t]\}, t] }\)

\(\displaystyle{ m1(t) = e^{k1t}\int_{1}^{t} e^{-\xi k1- a (\xi- t0)} y(0 )d\xi + c1 e^{k1 t}. }\)

\(\displaystyle{ m2(t)=\frac{1}{(k1+ k2)} e^{-k2 t}\left [ (k1+k2)\int_{1}^{t} \frac{e^{-\xi k1 - a( \xi -t0)} (-1 + e^{\xi (k1 + k2)})k1 y0 }{(k1 +k2)} d\xi +
k1 (e^{(k1+k2)t} -1) \int_{1}^{t} e^{-\xi k1 - a(\xi - t0)} y0 d \xi \right] +\frac{ c1 k1e^{-k 2 t}( e^{k1t +k2t}-1)}{(k1+k2)} + c2 e^{-k2 t}. }\)

Z tych równań wynika, że proces wchłaniania ponadnaczyniowej jednorazowej dawki leku, w otwartym modelu jednokompartnetowym , opisany zmianą jego masy (stężenia) w czasie \(\displaystyle{ m_{2}(t) }\) we krwi jest bardziej skomplikowany niż w układzie żołądkowo-jelitowym, bo należy uwzględnić zmianę masy leku w czasie zarówno w układzie żołądkowo - jelitowym jak i w układzie krwionośnym.

[Marakech]

Dziękuję za pomoc w zadaniu, pozdrawiam.