利用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()