利用Python繪制血葯濃度-時間曲線——口服吸收一室模型

血葯濃度-時間曲線一般是通過擬合所測定的血葯濃度點而畫出來的，但是在某些時候，如閱讀文獻時，我們需要根據別人報道的PK參數來畫出葯時曲線。Python語法簡單，擁有豐富的開源庫，下面嘗試通過Python來根據已知的血葯濃度公式繪制葯時曲線。

我們選擇口服給葯、符合單室模型的葯物來畫圖，因為這種葯物較為常見，公式也很簡單，並將2個模型的血葯濃度畫在同一張圖上，便於比較。

首先繪制單劑量給葯，對於model1，我們設置給葯劑量$X_0=5mg$，生物利用度$F=100%$，表觀分布容積$V=20L$，一級吸收速率常數$k_a=1$，一級消除速率常數$k=0.1$；對於model2，我們設置給葯劑量$X_0=5mg$，生物利用度$F=100%$，表觀分布容積$V=40L$，一級吸收速率常數$k_a=1$，一級消除速率常數$k=0.1$。

隔室模型示意圖如下：

model1的血葯濃度表達式為：

$$C=\frac{k_a\cdot X_0\cdot F}{V(k_a-k)}(e^{-kt}-e^{-k_at})=\frac{1\times 5 \times 1}{20\times (1-0.1)}(e^{-1}-e^{-0.1})$$

model2的血葯濃度表達式為：

$$C=\frac{k_a\cdot X_0\cdot F}{V(k_a-k)}(e^{-kt}-e^{-k_at})=\frac{1\times 5 \times 1}{40\times (1-0.1)}(e^{-1}-e^{-0.1})$$

有了血葯濃度公式，下面就可以畫圖了。

import numpy as np
import matplotlib.pyplot as plt
 
e=np.e 
plt.figure(figsize=(6,4)) 
t= np.linspace(0, 20, 1000)
 
plt.plot(t,5/(20*(1-0.1))*(e**(-0.1*t)-e**(-1*t)),label='model 1')
plt.legend()
plt.plot(t,5/(40*(1-0.1))*(e**(-0.1*t)-e**(-1*t)),label='model 2')
plt.legend() 
plt.title('single dose')
plt.xlabel('t (h)')
plt.ylabel('c (ng/ml)')
plt.show()

單劑量搞定，下面考慮多劑量給葯的圖形。對於多劑量給葯，我們一般需要使血葯濃度達到穩態。

以model1為例，首先求算其表達式。設置給葯間隔$\tau =5h$，給葯次數$n=10$。在第1個給葯周期內，其血葯濃度為第1次給葯后的濃度；在第2個給葯周期內，其血葯濃度為第1次給葯后5~10 h內的濃度加上第2次給葯后0~5 h內的濃度；以此類推。顯然，多劑量給葯，其體內葯量是一個分段函數：

$\begin{cases}
X_0\cdot e^{-kt} & \text{  } n=1, 0\leq t\leq 5 \\
X_0\cdot e^{-kt}+X_1\cdot e^{-k\cdot(t-5)} & \text{ } n=2, 5\leq t\leq 10 \\
X_0\cdot e^{-kt}+X_1\cdot e^{-k\cdot(t-5)}+X_2\cdot e^{-k\cdot(t-10)} & \text{ } n=3, 5\leq t\leq 10 \\
\cdots \\
X_0\cdot e^{-kt}+X_1\cdot e^{-k\cdot(t-5)}+X_2\cdot e^{-k\cdot(t-10)}+\cdots +X_n\cdot e^{-k\cdot [t-5(n-1)]} & \text{  } n=n, 5(n-1)\leq t\leq 5n
\end{cases}$

分段函數圖像分段畫，代碼如下：

import numpy as np
import matplotlib.pyplot as plt

e=np.e
t= np.linspace(0,49,1000)

#intervals
i1 = [1 if (i>=0 and i<=5) else 0 for i in t]
i2 = [1 if (i>=5 and i<=10) else 0 for i in t]
i3 = [1 if (i>=10 and i<=15) else 0 for i in t]
i4 = [1 if (i>=15 and i<=20) else 0 for i in t]
i5 = [1 if (i>=20 and i<=25) else 0 for i in t]
i6 = [1 if (i>=25 and i<=30) else 0 for i in t]
i7 = [1 if (i>=30 and i<=35) else 0 for i in t]
i8 = [1 if (i>=35 and i<=40) else 0 for i in t]
i9 = [1 if (i>=40 and i<=45) else 0 for i in t]
i10 = [1 if (i>=45 and i<=49) else 0 for i in t]

def c(t):
    return 5/(20*(1-0.1))*(e**(-0.1*t)-e**(-1*t))

y=c(t)*i1+\
(c(t)+c(t-5))*i2+\
(c(t)+c(t-5)+c(t-10))*i3+\
(c(t)+c(t-5)+c(t-10)+c(t-15))*i4+\
(c(t)+c(t-5)+c(t-10)+c(t-15)+c(t-20))*i5+\
(c(t)+c(t-5)+c(t-10)+c(t-15)+c(t-20)+c(t-25))*i6+\
(c(t)+c(t-5)+c(t-10)+c(t-15)+c(t-20)+c(t-25)+c(t-30))*i7+\
(c(t)+c(t-5)+c(t-10)+c(t-15)+c(t-20)+c(t-25)+c(t-30)+c(t-35))*i8+\
(c(t)+c(t-5)+c(t-10)+c(t-15)+c(t-20)+c(t-25)+c(t-30)+c(t-35)+c(t-40))*i9+\
(c(t)+c(t-5)+c(t-10)+c(t-15)+c(t-20)+c(t-25)+c(t-30)+c(t-35)+c(t-40)+c(t-45))*i10

plt.plot(t,y)
plt.show()

畫出model1給葯次數為10次，給葯間隔為$5h$時的葯時曲線：

由於分段函數需要分段表示，必須要用語句划分區間，並且每一段的表達式都要寫出來，導致上述代碼比較長。在葯代動力學教材中，通過引入多劑量函數（$r=\frac{1-e^{-nk\tau }}{1-e^{-k\tau}}$）這個概念，可利用$n$和$t$把第$n$次給葯的體內葯量$X(n,t)$的通式表示出來：

$$X_n=X_0\cdot \frac{1-e^{-nk\tau }}{1-e^{-k\tau}}\cdot e^{-kt}$$

則第$n$次給葯后，血葯濃度$C_n$為：

$$X_n=\frac{X_0}{V}\cdot \frac{1-e^{-nk\tau }}{1-e^{-k\tau}}\cdot e^{-kt}$$

 根據此公式畫圖：

import numpy as np
import matplotlib.pyplot as plt

e=np.e
t= np.linspace(0,49,1000)

i1 = [1 if (i>=0 and i<=5) else 0 for i in t]
i2 = [1 if (i>=5 and i<=10) else 0 for i in t]
i3 = [1 if (i>=10 and i<=15) else 0 for i in t]
i4 = [1 if (i>=15 and i<=20) else 0 for i in t]
i5 = [1 if (i>=20 and i<=25) else 0 for i in t]
i6 = [1 if (i>=25 and i<=30) else 0 for i in t]
i7 = [1 if (i>=30 and i<=35) else 0 for i in t]
i8 = [1 if (i>=35 and i<=40) else 0 for i in t]
i9 = [1 if (i>=40 and i<=45) else 0 for i in t]
i10 = [1 if (i>=45 and i<=49) else 0 for i in t]

def c(t,n):
    return 5/(20*(1-0.1))*((1-e**(-0.1*5*n))/(1-e**(-5*0.1))*e**(-0.1*t)-(1-e**(-1*5*n))/(1-e**(-1*5))*e**(-1*t))

y=c(t,1)*i1+c(t-5,2)*i2+c(t-10,3)*i3+c(t-15,4)*i4+c(t-20,5)*i5+c(t-25,6)*i6+c(t-30,7)*i7+c(t-35,8)*i8+c(t-40,9)*i9+c(t-45,10)*i10

plt.plot(t,y)
plt.show()

 血葯濃度圖像如下：

與上一段代碼比起來，這段代碼在表示y的時候簡潔了許多，但仍然需要划分區間，這就導致在改變給葯次數$n$后，畫圖需要改動代碼，並且$n$比較大時，代碼也比較長。那么如在讓我們設置任意的$n$時，圖形都可以方便地畫出來呢？

之前我們是根據表達式直接畫圖（雖然Python實現起來是通過生成很多散點），而表達式一旦是分段函數的話，畫起來就很麻煩，必須分段（因此可以看到代碼中有很長的一段是在分割區間）。現在可以轉換一下思路，不通過表達式畫圖，而是通過表達式生成散點，再通過散點作圖。利用分段函數產生點比利用分段函數畫圖方便多了，尤其是$C(n,t)$的通式已經給出了，只需一個簡單的循環語句就能生成每一次給葯后的濃度點，再通過Python的列表（list）或數組（array）的append()方法，把不同時間段的濃度點連接起來，就是我們所需要的y（全程的濃度），根據時間散點和濃度散點，就可以把圖畫出來了。並且給葯次數可以任意設置。

給葯次數設為10次：

import numpy as np
import matplotlib.pyplot as plt

def c(t,n):
    return 5/(20*(1-0.1))*((1-e**(-0.1*5*n))/(1-e**(-5*0.1))*e**(-0.1*t)-(1-e**(-1*5*n))/(1-e**(-1*5))*e**(-1*t))

e=np.e
t1= np.linspace(0,4,100)
y1=c(t1,1)
t=t1
y=y1

for N in range(2,10):
    tN=t1+(N-1)*5
    yN=c(t1,N)
    t=np.append(t,tN)
    y=np.append(y,yN)

plt.plot(t,y)
plt.show()

 

如果想把給葯次數設為50，只需將第13行代碼的10改成50即可：

import numpy as np
import matplotlib.pyplot as plt
 
def c(t,n):
    return 5/(20*(1-0.1))*((1-e**(-0.1*5*n))/(1-e**(-5*0.1))*e**(-0.1*t)-(1-e**(-1*5*n))/(1-e**(-1*5))*e**(-1*t))

e=np.e
t1= np.linspace(0,4,100)
y1=c(t1,1)
t=t1
y=y1

for N in range(2,50):
    tN=t1+(N-1)*5
    yN=c(t1,N)
    t=np.append(t,tN)
    y=np.append(y,yN)
 
plt.plot(t,y)
plt.show()

 

由於開頭說比較model1和model2，要把2條曲線畫在一起。所以最后再畫上model2的葯時曲線，添上標題、圖例、橫縱坐標，畫圖完成。

import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(8,6)) 

def c1(t,n):
    return 5/(20*(1-0.1))*((1-e**(-0.1*5*n))/(1-e**(-5*0.1))*e**(-0.1*t)-(1-e**(-1*5*n))/(1-e**(-1*5))*e**(-1*t))
def c2(t,n):
    return 5/(20*(1-0.2))*((1-e**(-0.2*5*n))/(1-e**(-5*0.2))*e**(-0.2*t)-(1-e**(-1*5*n))/(1-e**(-1*5))*e**(-1*t))

e=np.e
t0= np.linspace(0,4,100)
t=t0
y1=c1(t0,1)
y2=c2(t0,1)

for N in range(2,20):
    tN=t0+(N-1)*5
    yN1=c1(t0,N)
    yN2=c2(t0,N)
    t=np.append(t,tN)
    y1=np.append(y1,yN1)
    y2=np.append(y2,yN2)

plt.plot(t,y1,label='model 1')
plt.plot(t,y2,label='model 2')
plt.legend()
plt.title('multiple dose')
plt.xlabel('t (h)')
plt.ylabel('c (ng/ml)')
plt.show()