使用 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})$$
有了血药浓度公式，下面就可以画图了。
1 import numpy as np
2 import matplotlib.pyplot as plt
3
4 e=np.e
5 plt.figure(figsize=(6,4))
6 t= np.linspace(0, 20, 1000)
7
8 plt.plot(t,5/(20*(1-0.1))*(e**(-0.1*t)-e**(-1*t)),label='model 1')
9 plt.legend()
10 plt.plot(t,5/(40*(1-0.1))*(e**(-0.1*t)-e**(-1*t)),label='model 2')
11 plt.legend()
12 plt.title('single dose')
13 plt.xlabel('t (h)')
14 plt.ylabel('c (ng/ml)')
15 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}$
分段函数图像分段画，代码如下：
1 import numpy as np
2 import matplotlib.pyplot as plt
3
4 e=np.e
5 t= np.linspace(0,49,1000)
6
7 #intervals
8 i1 = [1 if (i>=0 and i<=5) else 0 for i in t]
9 i2 = [1 if (i>=5 and i<=10) else 0 for i in t]
10 i3 = [1 if (i>=10 and i<=15) else 0 for i in t]
11 i4 = [1 if (i>=15 and i<=20) else 0 for i in t]
12 i5 = [1 if (i>=20 and i<=25) else 0 for i in t]
13 i6 = [1 if (i>=25 and i<=30) else 0 for i in t]
14 i7 = [1 if (i>=30 and i<=35) else 0 for i in t]
15 i8 = [1 if (i>=35 and i<=40) else 0 for i in t]
16 i9 = [1 if (i>=40 and i<=45) else 0 for i in t]
17 i10 = [1 if (i>=45 and i<=49) else 0 for i in t]
18
19 def c(t):
20 return 5/(20*(1-0.1))*(e**(-0.1*t)-e**(-1*t))
21
22 y=c(t)*i1+\
23 (c(t)+c(t-5))*i2+\
24 (c(t)+c(t-5)+c(t-10))*i3+\
25 (c(t)+c(t-5)+c(t-10)+c(t-15))*i4+\
26 (c(t)+c(t-5)+c(t-10)+c(t-15)+c(t-20))*i5+\
27 (c(t)+c(t-5)+c(t-10)+c(t-15)+c(t-20)+c(t-25))*i6+\
28 (c(t)+c(t-5)+c(t-10)+c(t-15)+c(t-20)+c(t-25)+c(t-30))*i7+\
29 (c(t)+c(t-5)+c(t-10)+c(t-15)+c(t-20)+c(t-25)+c(t-30)+c(t-35))*i8+\
30 (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+\
31 (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
32
33 plt.plot(t,y)
34 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}$$
根据此公式画图：
1 import numpy as np
2 import matplotlib.pyplot as plt
3
4 e=np.e
5 t= np.linspace(0,49,1000)
6
7 i1 = [1 if (i>=0 and i<=5) else 0 for i in t]
8 i2 = [1 if (i>=5 and i<=10) else 0 for i in t]
9 i3 = [1 if (i>=10 and i<=15) else 0 for i in t]
10 i4 = [1 if (i>=15 and i<=20) else 0 for i in t]
11 i5 = [1 if (i>=20 and i<=25) else 0 for i in t]
12 i6 = [1 if (i>=25 and i<=30) else 0 for i in t]
13 i7 = [1 if (i>=30 and i<=35) else 0 for i in t]
14 i8 = [1 if (i>=35 and i<=40) else 0 for i in t]
15 i9 = [1 if (i>=40 and i<=45) else 0 for i in t]
16 i10 = [1 if (i>=45 and i<=49) else 0 for i in t]
17
18 def c(t,n):
19 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))
20
21 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
22
23 plt.plot(t,y)
24 plt.show()血药浓度图像如下：
与上一段代码比起来，这段代码在表示y的时候简洁了许多，但仍然需要划分区间，这就导致在改变给药次数$n$后，画图需要改动代码，并且$n$比较大时，代码也比较长。那么如在让我们设置任意的$n$时，图形都可以方便地画出来呢？
之前我们是根据表达式直接画图（虽然Python实现起来是通过生成很多散点），而表达式一旦是分段函数的话，画起来就很麻烦，必须分段（因此可以看到代码中有很长的一段是在分割区间）。现在可以转换一下思路，不通过表达式画图，而是通过表达式生成散点，再通过散点作图。利用分段函数产生点比利用分段函数画图方便多了，尤其是$C(n,t)$的通式已经给出了，只需一个简单的循环语句就能生成每一次给药后的浓度点，再通过Python的列表（list）或数组（array）的append()方法，把不同时间段的浓度点连接起来，就是我们所需要的y（全程的浓度），根据时间散点和浓度散点，就可以把图画出来了。并且给药次数可以任意设置。
给药次数设为10次：
1 import numpy as np
2 import matplotlib.pyplot as plt
3
4 def c(t,n):
5 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))
6
7 e=np.e
8 t1= np.linspace(0,4,100)
9 y1=c(t1,1)
10 t=t1
11 y=y1
12
13 for N in range(2,10):
14 tN=t1+(N-1)*5
15 yN=c(t1,N)
16 t=np.append(t,tN)
17 y=np.append(y,yN)
18
19 plt.plot(t,y)
20 plt.show() 
如果想把给药次数设为50，只需将第13行代码的10改成50即可：
1 import numpy as np
2 import matplotlib.pyplot as plt
3
4 def c(t,n):
5 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))
6
7 e=np.e
8 t1= np.linspace(0,4,100)
9 y1=c(t1,1)
10 t=t1
11 y=y1
12
13 for N in range(2,50):
14 tN=t1+(N-1)*5
15 yN=c(t1,N)
16 t=np.append(t,tN)
17 y=np.append(y,yN)
18
19 plt.plot(t,y)
20 plt.show()
由于开头说比较model1和model2，要把2条曲线画在一起。所以最后再画上model2的药时曲线，添上标题、图例、横纵坐标，画图完成。
1 import numpy as np
2 import matplotlib.pyplot as plt
3
4 plt.figure(figsize=(8,6))
5
6 def c1(t,n):
7 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))
8 def c2(t,n):
9 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))
10
11 e=np.e
12 t0= np.linspace(0,4,100)
13 t=t0
14 y1=c1(t0,1)
15 y2=c2(t0,1)
16
17 for N in range(2,20):
18 tN=t0+(N-1)*5
19 yN1=c1(t0,N)
20 yN2=c2(t0,N)
21 t=np.append(t,tN)
22 y1=np.append(y1,yN1)
23 y2=np.append(y2,yN2)
24
25 plt.plot(t,y1,label='model 1')
26 plt.plot(t,y2,label='model 2')
27 plt.legend()
28 plt.title('multiple dose')
29 plt.xlabel('t (h)')
30 plt.ylabel('c (ng/ml)')
31 plt.show()