0917-062-0010

مشاوره رایگان

9 صبح تا 9 شب

شنبه تا پنجشنبه

d8a8db8cd988d981d8a7d8b1d985d8a7d8b3db8c d8afd8b1 d9bed8a7db8cd8aad988d986 d8a8d8aed8b4 dbb1 d8b3db8cd8b3d8aad985 db8cdaa9 d8a8d8ae 6574c0ffd3fd7

در این مطلب یکی از مدل‌های ساده و پرکاربرد برای توزیع داروهای تزریقی در خون را بررسی می‌کنیم که به آن سیستم یک‌بخشی عروقی گفته می‌شود.

فهرست مطالب این نوشته
مدل ریاضی توزیع دارو در خون

مدل‌سازی با داده‌های واقعی در پایتون

جمع‌بندی

faradars mobile

مدل ریاضی توزیع دارو در خون

در شرایطی که دارو به صورت مستقیم وارد خون می‌شود و به صورت درجه یک از خون حذف می‌شود، به صورت زیر خواهیم داشت:

$$ {dC over dt} = -K times C$$

که در این رابطه، $$C$$ نشان‌دهنده‌ غلظت پلاسمایی دارو $$t$$ نشان‌دهنده‌ زمان و عدد $$K$$ نیز برابر با ثابت سرعت حذف دارو است.

آموزش برنامه نویسی پایتون + مثال های عملی در Python
فیلم آموزش برنامه نویسی پایتون + مثال های عملی در Python در تم آف

کلیک کنید

با اندکی تغییر رابطه‌ بالا به فرم زیر می‌رسیم:

$$ {dC over C} = -K cdot dt$$

حال می‌توانیم از طرفین انتگرال بگیریم:

$$ int {1 over C} cdot dC = space – int K cdot dt $$

$$ {ln space C} = {space ln space C_0} {-} {Kt} $$

حال می‌توان مقدار $$C$$ را به شکل زیر به دست آورد:

$$ {ln space C} – {ln space C_0} = {- Kt} $$

$$ {ln} {large(} { {C} over C_0} {large) = -Kt} $$

$$ frac{C}{C_{0}}=e^{-K t}$$

$$ C=C_{0} cdot e^{-K t} $$

به این شکل مدل نهایی به دست آمد. حال می‌خواهیم یک مدل‌سازی روی داده‌های واقعی انجام دهیم.

آموزش فارماکولوژی یا داروشناسی داروهای موثر بر دستگاه گوارش، خون و روماتولوژی
فیلم آموزش فارماکولوژی یا داروشناسی داروهای موثر بر دستگاه گوارش، خون و روماتولوژی در تم آف

کلیک کنید

مدل‌سازی با داده‌های واقعی در پایتون

فرض کنید 50 میلی‌گرم دارو به صورت تک‌دوز به بیماری تجویز شده است و غلظت پلاسمایی دارو در زمان‌های مختلف در جدول زیر آورده شده است. با مدل‌سازی مناسب، ثابت سرعت حذف، سطح زیر نمودار تا بی‌نهایت، نیمه‌عمر و کلیرانس را محاسبه کنید.

Concentration (mg/L) Time (h)
1٫84 ۰
1٫56 ۰٫۵
1٫52 ۱
1٫17 ۲
0٫68 ۴
0٫48 ۶
0٫24 ۸
0٫18 ۱۰

ابتدا وارد محیط برنامه‌نویسی شده و کتابخانه‌های مورد نیاز را فراخوانی می‌کنیم:

import math
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

این کتابخانه‌ها به ترتیب برای موارد زیر استفاده خواهند شد:

  1. محاسبات ریاضی
  2. کار با آرایه و محاسبات برداری
  3. بهینه‌سازی مدل
  4. رسم نمودار
آموزش کتابخانه های NumPy و Matplotlib در پایتون
فیلم آموزش کتابخانه های NumPy و Matplotlib در پایتون در تم آف

کلیک کنید

حال تولید اعداد تصادفی را ثابت کرده و Style نمودارها را تنظیم می‌کنیم:

np.random.seed(0)
plt.style.use('ggplot')

سطر اول باعث می‌شود در صورت استفاده از اعداد تصادفی در برنامه، در دفعات متفاوت از اجرای کد، نتایج یکسانی بگیریم و نتایج قابل بازتولید باشند. در سطر دوم نیز Style معروف ggplot استفاده شده است. حال داده‌ها را به صورت آرایه تعریف می‌کنیم:

T = np.array([0, 0.5, 1, 2, 4, 6, 8, 10])
C = np.array([1.84, 1.56, 1.52, 1.17, 0.68, 0.48, 0.24, 0.18])

قبل از شروع، داده‌ها را مصورسازی می‌کنیم. برای این کار می‌توانیم از Plot استفاده کنیم:

plt.plot(T, C, lw=1.2, marker='o', ms=3, label='Plasma Concentration', c='crimson')
plt.title('Data Visualization')
plt.xlabel('Time (h)')
plt.ylabel('Concentration (mg/L)')
plt.axhline(c='k')
plt.axvline(c='k')
plt.legend()
plt.show()

بنابراین، در خروجی شکل زیر را داریم.

شکل خروجی داده‌ها

مشاهده می‌کنیم که نمودار حاصل، به نمودار تابع نمایی نزدیک است، بنابراین سیستم یک‌بخشی عروقی می‌تواند مورد استفاده قرار گیرد.

آموزش کتابخانه SciPy برای محاسبات علمی در پایتون – بخش یکم
فیلم آموزش کتابخانه SciPy برای محاسبات علمی در پایتون – بخش یکم در تم آف

کلیک کنید

اکنون باید مدل را به صورت یک تابع تعریف کنیم. از آن‌جا که محاسبات برای مجموعه‌ای از زمان‌ها انجام می‌شود، باید تابع را متناسب با این نیاز بسازیم. این تابع در ورودی مقدار $$K$$، $$C_0$$ و زمان‌های مورد نظر را خواهد گرفت و در خروجی غلظت در تک تک زمان‌های ورودی را خواهد داد. بنابراین، خواهیم داشت:

def Model (K:float, C0:float, T:np.ndarray):
    C = C0 * np.exp( -K*T )
    return C

توجه داشته باشید که دو ورودی اول تابع اعدادی اعشاری هستند، ولی ورودی سوم یک بردار از نوع numpy.ndarray است. به دلیل اینکه در $$T$$ به جای یک عدد، با مجموعه‌ای از اعداد کار می‌کنیم، باید از تابع numpy.exp به جای math.exp استفاده کنیم. خروجی این تابع نیز یک بردار از نوع numpy.ndarray است.

حالا می‌توانیم با دادن $$K$$ و $$C_0$$ مناسب به تابع Model مقدار غلظت دارو در پلاسما را پیش‌بینی کنیم. ولی مقدار مناسب $$K$$ و $$C_0$$ معلوم نیست و باید آن را تعیین کنیم. برای این کار یک تابع خطا باید تعریف کنیم و آن را بهینه کنیم.

آموزش برنامه نویسی پایتون + مثال های عملی در Python
فیلم آموزش برنامه نویسی پایتون + مثال های عملی در Python در تم آف

کلیک کنید

تابع خطا با گرفتن داده‌های اصلی، $$K$$ و $$C_0$$ خواهد توانست خطای حاصل را برگرداند که یک عدد است. بنابراین داریم:

def Error (Parameters:np.ndarray, T:np.ndarray, C:np.ndarray):
    K = Parameters[0]
    C0 = Parameters[1]
    Prediction = Model(K, C0, T)
    RMSE = np.mean( np.power(C – Prediction, 2) )**0.5
    return RMSE

نکته‌ای که باید به آن توجه کرد این است که $$K$$ و $$C_0$$ پارامترهای مدل هستند و باید به صورت یک آرایه با هم به تابع وارد شوند. بنابراین به صورت قراردادی، عضو اول Parameters را برابر با $$K$$ و عضو دوم آن را برابر با $$C_0$$ در نظر می‌گیریم.

در این تابع به ترتیب عملیات زیر را انجام می‌دهیم:

  1. مقدار $$K$$ را استخراج می‌کنیم.
  2. مقدار $$C_0$$ را استخراج می‌کنیم.
  3. با استفاده از تابع Model پیش‌بینی را انجام می‌دهیم.
  4. میانگین قدرمطلق خطاها را محاسبه می‌کنیم.
    • به این خطا به اختصار RMSE گفته می‌شود که کوتاه‌شده عبارت Root Mean Squared Error است.
    • توجه داشته باشید که چون $$Y$$ و Prediction آرایه‌ای از اعداد هستند، باید از تابع numpy.power به جای pow استفاده کنیم.
  5. مقدار خطا را خروجی می‌گیریم.
آموزش برنامه نویسی پایتون + مثال های عملی در Python
فیلم آموزش برنامه نویسی پایتون + مثال های عملی در Python در تم آف

کلیک کنید

به این صورت تابع خطا نیز تعریف می‌شود.

اکنون باید یک حدس اولیه از پارامترها داشته باشیم. با توجه به اینکه غلظت ثبت شده در زمان $$T=0$$ برابر با $$1.84$$‌است، احتمالاً همین عدد برای غلظت اولیه مناسب خواهد بود. برای ثابت سرعت حذف نیز می‌توان عدد $$0.3$$ را محل شروع در نظر گرفت.

Parameters0 = np.array([0.3, 1.84])

حال می‌توانیم بهینه‌سازی را انجام دهیم:

Result = opt.minimize(Error, Parameters0, args=(T, C))

این خط کد، تابع Error را با تغییر پارامترها بهینه می‌کند. با توجه به اینکه تابع خطا به جز پارامترها، دو آرایه‌ دیگر مربوط به زمان و غلظت پلاسمایی را نیز در ورودی می‌گیرد، این دو مورد را در بخش args به صورت یک tuple تعریف می‌کنیم. خروجی این بخش کد یک دیکشنری با اسم Result خواهد بود که کمترین مقدار تابع، بهترین جواب، نتیجه‌ بهینه‌سازی را در خود جای داده است.

آموزش برنامه نویسی پایتون + مثال های عملی در Python
فیلم آموزش برنامه نویسی پایتون + مثال های عملی در Python در تم آف

کلیک کنید

از بین این موارد باید بهترین جواب را استخراج کنیم که با کلید ‘x’ در دسترس است. بنابراین می‌توان نوشت:

Parameters = Result['x']

حال از روی Parameters می‌توانیم مقدار $$K$$ و $$C_0$$ را استخراج کنیم:

K = Parameters[0]
C0 = Parameters[1]

برای مشاهده‌ مقادیر پارامترها نیز، اعداد را تا 4 رقم اعشار در خروجی نمایش می‌دهیم:

print(f'K: {round(K, 4)} 1/h')
print(f'C0: {round(C0, 4)} mg/L')

پس از اجرا، نتایج زیر در خروجی ظاهر می‌شود.

K: 0.2352 1/h
C0: 1.8359 mg/L

همان‌طور که انتظار داشتیم، مقدار $$C_ 0$$ عددی نزدیک به $$1.84$$ ظاهر شد.

آموزش برنامه نویسی پایتون + مثال های عملی در Python
فیلم آموزش برنامه نویسی پایتون + مثال های عملی در Python در تم آف

کلیک کنید

حال می‌توانیم کمترین مقدار خطا را نیز از Result استخراج کرده و نمایش دهیم:

RMSE = Result['fun']
print(f'RMSE: {round(RMSE, 4)} mg/L')

پس از اجرا خواهیم داشت:

RMSE: 0.0426 mg/L

به این صورت مشاهده می‌کنیم که به خطای مناسبی دست یافته‌ایم. توجه داشته باشید که واحد خطای RMSE برابر با داده‌های اصلی یعنی mg/L است.

آموزش برنامه نویسی پایتون + مثال های عملی در Python
فیلم آموزش برنامه نویسی پایتون + مثال های عملی در Python در تم آف

کلیک کنید

حال می‌توانیم سایر مقادیر خواسته شده را محاسبه کنیم:

سطح زیر نمودار تا بینهایت

برای به دست آوردن این مقدار، باید از رابطه‌ اصلی به دست آمده انتگرال بگیریم:

$$begin{aligned}
int_{0}^{infty} C cdot d t=& int_{0}^{infty} C_{0} cdot e^{-K t} cdot d t=C_{0} int_{0}^{infty} e^{-K t} cdot d t=C_{0}left(frac{e^{-infty}}{-K}-frac{e^{0}}{-K}right) \
&=C_{0}left(frac{0}{-K}-frac{1}{-K}right)=frac{C_{0}}{K}
end{aligned}$$

حال، با وارد کردن رابطه در برنامه داریم:

AUC = C0/K
print(f'AUC: {round(AUC, 4)} mg.h/L')

که پس از اجرا به نتیجه‌ زیر می‌رسیم:

AUC: 7.8035 mg.h/L

نیمه عمر

بر اساس تعریف، نیمه عمر برابر است با مدت زمانی که طول می‌کشد تا غلظت دارو نصف شود. بنابراین داریم:

$$ begin{align*}
frac{1}{2} C_{0} & =C_{0} cdot e^{-K t_{h l}} \ Rightarrow 2^{-1}&=e^{-K t_{h l}} \ Rightarrow 2 &=e^{K t_{h l}} \ Rightarrow ln 2 &=K t_{h l} \ Rightarrow t_{h l}
&=frac{ln 2}{K}
end{align*} $$

حال خواهیم داشت:

HL = math.log(2)/K
print(f'HL: {round(HL, 4)} h')

که در خروجی به این حالت داریم:

HL: 2.9482 h

این عدد با داده‌های واقعی نیز سازگار است و می‌توانیم روی نمودار به صورت تقریبی این نقطه را مشاهده کنیم.

کلیرانس

با توانایی بدن در پاک‌سازی حجم توزیع از دارو متناسب است و می‌توان از فرمول زیر مقدار آن را به دست آورد:

$$C l=frac{D}{A U C}$$

توجه داشته باشید که در این فرمول، $$D$$ مقدار داروی تزریقی است، پس خواهیم داشت:

D = 50
Cl = D/AUC
print(f'Cl: {round(Cl, 4)} L/h')

و در خروجی نیز نتیجه به این صورت خواهد بود:

Cl: 6.4074 L/h

حال می‌توانیم نمودار حاصل از مدل‌سازی را با داده‌های اصلی مقایسه کنیم. برای این کار ابتدا با استفاده از پارامترهای به دست آمده، پیش‌بینی انجام می‌دهیم:

Prediction = Model(K, C0, T)

اکنون می‌توانیم دو نمودار را در کنار هم رسم کنیم:

plt.plot(T, C, lw=1.3, marker='o', ms=3, label='Real', c='crimson')
plt.plot(T, Prediction, lw=1.1,marker='x', ms=4, label='Predicted', c='teal')
plt.title('Model Result')
plt.xlabel('Time (h)')
plt.ylabel('Concentration (mg/L)')
plt.axhline(c='k')
plt.axvline(c='k')
plt.legend()
plt.show()

و در خروجی شکل زیر را داشته باشیم.

نمودار کلیرانس

مشاهده می‌کنیم که مدل دقت خوبی از خود نشان داده است.

آموزش برنامه نویسی پایتون + مثال های عملی در Python
فیلم آموزش برنامه نویسی پایتون + مثال های عملی در Python در تم آف

کلیک کنید

می‌توان خطای مدل و رابطه‌ نهایی را نیز وارد نمودار کرد و به این صورت نوشت:

plt.plot(T, C, lw=1.3, marker='o', ms=3, label='Real', c='crimson')
plt.plot(T, Prediction, lw=1.1,marker='x', ms=4, label='Predicted', c='teal')
plt.text(2, 1.5, f'C = {round(C0, 4)} * exp(-{round(K, 4)} * t)')
plt.text(4, 1.25, f'RMSE = {round(RMSE, 4)}')
plt.title('Model Result')
plt.xlabel('Time (h)')
plt.ylabel('Concentration (mg/L)')
plt.axhline(c='k')
plt.axvline(c='k')
plt.legend()
plt.show()

تا در خروجی شکل زیر را داشته باشیم.

نمودار خطا

به این صورت می‌توانیم نتیجه‌ کار را به صورت جامع و کامل در یک نمودار نشان دهیم.

آموزش برنامه نویسی پایتون + مثال های عملی در Python
فیلم آموزش برنامه نویسی پایتون + مثال های عملی در Python در تم آف

کلیک کنید

بهتر است در این‌گونه مسائل، کمترین و بیشترین درصد خطا نیز آورده شود. برای این منظور، می‌توانیم درصد قدرمطلق خطا (APE یا Absolute Percentage Error) را با استفاده از فرمول زیر محاسبه کنیم:

$$
A P E=100 timesleft|frac{y-hat{y}}{y}right|
$$

درون کد می‌توانیم با استفاده از numpy این بخش را پیاده‌سازی کنیم:

APE = 100*np.abs( np.divide(C-Prediction, C) )

و در نهایت کمترین، بیشترین و میانگین درصد قدرمطلق خطاها را بیابیم و در خروجی نمایش دهیم:

MinAPE = np.min(APE)
MeanAPE = np.mean(APE)
MaxAPE = np.max(APE)

print(f'Minimum APE: {round(MinAPE, 4)} %')
print(f'Mean    APE: {round(MeanAPE, 4)} %')
print(f'Maximum APE: {round(MaxAPE, 4)} %')

که خواهیم داشت:

Minimum APE: 0.2894 %
Mean APE: 5.3739 %
Maximum APE: 16.5446 %

جمع‌بندی

در این مطلب یاد گرفتیم که مدل یک‌بخشی عروقی را روی داده‌های واقعی برازش (Fit) کنیم، پارامترهای مورد نیاز را محاسبه و مدل را ارزیابی کنیم.

مجموعه آموزش برنامه نویسی پایتون (Python)
فیلم مجموعه آموزش برنامه نویسی پایتون (Python) در تم آف

کلیک کنید

نتایج را به صورت نمودار مشاهده کنیم. مدل ارائه‌شده حالت ساده‌ای از سیستم و ارتباط موجود بین دارو و بدن است. در مطالب بعدی به مدل‌های برای مدل‌سازی شرایط پیچیده‌تر خواهیم پرداخت.

ارسال پاسخ

This site is protected by reCAPTCHA and the Google Privacy Policy and Terms of Service apply.