有些药物本身不具有药理活性，进入体内后经代谢产生活性物质，也有一些药物，其自身与代谢物都具有药理活性，对于此类药物，我们需要关注它们的代谢物的浓度。
利用隔室模型进行研究时，步骤是根据隔室模型写出微分方程组，求出表达式，画图。对于常见的一室、二室、三室模型，药代动力学课本上都给出了血药浓度公式，但是有了代谢物之后，就不是应用上述模型，换句话说，没有课本上现成的公式可以套用了。
那么如何画图呢？我们一般会想到2种方法：
- 根据微分方程组解出方程，通过方程画图；
- 直接根据微分方程画图。
对于以上2种方法，可以综合使用。
下面以抗精神分裂药利培酮（Risperidone）为例，画出其母药、代谢物以及二者之和的浓度曲线。利培酮在体内代谢为9-羟利培酮，它们都具有药理活性，近年来国内外有不少文献使用群体建模方法对该药物进行了研究。本文画图所使用的的PK参数来自如下文献：
- model1：Feng2008
- model2：Sherwin2012
- model3：Yoo2012
- model4：JSM2016
先解决绘制model1的问题。
文献中给出了隔室模型的示意图：
根据上图，写出微分方程组：
$$\frac{\mathrm{d} X_a}{\mathrm{d} t}=-K_a\cdot X_a \tag{1}$$
$$\frac{\mathrm{d} X_2}{\mathrm{d} t}=K_a\cdot X_a-K_{20}\cdot X_2-K_{23}\cdot X_2 \tag{2}$$
$$\frac{\mathrm{d} X_3}{\mathrm{d} t}=K_{23}\cdot X_2-K_{30}\cdot X_3 \tag{3}$$
其中（1）代表吸收部位药量，（2）代表利培酮药量，（3）代表9-羟利培酮药量。
从文献中找到PK参数的值：$K_a=1.7 h^{-1}$，$KF=0.595$，$CL=65.4L/h$，$V=444 L$，$VM=444 L$，$CLM=8.8 L/h$，我们设置给药剂量为2mg。
将参数带入微分方程组，可以算出方程（1）表达式：
$$X_a=2\cdot e^{-1.7t} \tag{4}$$
将参数和式（4）带入方程（2），方程可写为：
$$\frac{\mathrm{d} X_2}{\mathrm{d} t}=1.7\cdot 2\cdot e^{-1.7t}-0.148\cdot X_2 \tag{5}$$
可求出表达式：
$$X_2=2.19\cdot (e^{-0.148t}-e^{-1.7t}) \tag{6}$$
将参数和式（6）带入方程（3），方程可写为：
$$\frac{\mathrm{d} X_3}{\mathrm{d} t}=0.088\cdot 2.19\cdot (e^{-0.148t}-e^{-1.7t})-0.02\cdot X_3 \tag{7}$$
方程（7）暂时先不解。
有了以上的式子，可以着手用Python画图了。scipy这个科学计算库可以根据微分方程表达式直接画出图形而不必解方程。上面之所以解出了（1）和（2）的表达式，是因为（2）和（3）的方程里分别含有$X_a$和$X_2$。
1 # Feng2008 单剂量 2 3 import numpy as np 4 from scipy.integrate import odeint 5 import matplotlib.pyplot as plt 6 7 e=np.e 8 9 def diff_x2(x,t): 10 return 1.7*2*e**(-1.7*t)-0.148*x 11 def diff_x3(x,t): 12 return 0.088*2.19*(e**(-0.148*t)-e**(-1.7*t))-0.02*x 13 14 t=np.arange(0,100,0.01) 15 x2=odeint(diff_x2,0,t) 16 x3=odeint(diff_x3,0,t) 17 c2=x2/444*1000 18 c3=x3/444*1000 19 c_risperidone=c2+c3 20 21 plt.plot(t,c2,label='RISP') 22 plt.plot(t,c3,label='9-OH') 23 plt.plot(t,c_risperidone,label='RISP+9-OH') 24 plt.legend() 25 plt.title('Feng2008 single dose') 26 plt.xlabel('t (h)') 27 plt.ylabel('c (ng/ml)') 28 plt.show()
这段代码画出了利培酮（RISP）、9-羟利培酮（9-OH）以及二者浓度之和（RISP+9-OH）：
图画出来了，还有一些我们值得思考的地方：
第一，如上所述，解微分方程时，我们是笔算的，能否通过Python来求解微分方程？
第二，我们写出的其实一个微分方程组，但实际画图时是把它们当做独立的微分方程来画图，能否直接通过线性（微分）方程组画图？
对于问题一，答案是Python是可以解微分方程的。科学计算库sympy是一种方法，解式（1）：
1 import sympy as sy 2 3 def differential_equation(x,f): 4 return sy.diff(f(x),x,1)+1.7*f(x) 5 6 x=sy.symbols('t') 7 f=sy.Function('Ca') 8 9 print(sy.dsolve(differential_equation(x,f),f(x)))
解得表达式如下：
Eq(Ca(t), C1*exp(-1.7*t))
这个结果与笔算的结果一致，是正确的。
解式（2） ：
1 import numpy as np 2 import sympy as sy 3 4 e=np.e 5 6 def differential_equation(x,f): 7 return sy.diff(f(x),x,1)-3.4*e**(-1.7*x)+0.148*f(x) 8 9 x=sy.symbols('x') 10 f=sy.Function('f') 11 12 print(sy.dsolve(differential_equation(x,f),f(x)))
结果如下：
Eq(f(x), (C1 - 2.19072164948454*exp(-1.552*x))*exp(-0.148*x))
这个结果是错误的，改变一下写法：
1 # -- coding:utf-8 -- 2 3 from sympy import * 4 5 #用sympy符号运算解方程 6 x=symbols('x',real = True) # real 保证全是实数，自变量 7 y=symbols('y',function = True) # 全部为函数变量 8 eq=y(x).diff(x,1)-3.4*exp(-1.7*x)+0.148*y(x) 9 10 print(dsolve(Eq(eq,0),y(x)))
结果依然如此。错误的原因我还没有找到，此问题悬而未决。
对于问题二，目前我也没找到相关的方法，也有待解决。
画出了model1的单剂量给药的血药浓度图，多剂量该如何画呢？参考画一室模型多剂量给药的方法，先生成散点，由于没有借助通式，代码非常长。
1 # Feng2008 多剂量 2 3 import numpy as np 4 from scipy.integrate import odeint 5 import matplotlib.pyplot as plt 6 7 plt.figure(figsize=(8,4)) 8 e=np.e 9 10 def diff_x2(x,t): 11 return 1.7*2*e**(-1.7*t)-0.148*x 12 def diff_x3(x,t): 13 return 0.088*2.19*(e**(-0.148*t)-e**(-1.7*t))-0.02*x 14 15 t1=np.arange(0,120,0.4) 16 x2_1=odeint(diff_x2,0,t1) 17 x2_11=x2_1[range(0,30)] 18 x2_12=x2_1[range(30,60)] 19 x2_13=x2_1[range(60,90)] 20 x2_14=x2_1[range(90,120)] 21 x2_15=x2_1[range(120,150)] 22 x2_16=x2_1[range(150,180)] 23 x2_17=x2_1[range(180,210)] 24 x2_18=x2_1[range(210,240)] 25 x2_19=x2_1[range(240,270)] 26 x2_110=x2_1[range(270,300)] 27 28 x3_1=odeint(diff_x3,0,t1) 29 x3_11=x3_1[range(0,30)] 30 x3_12=x3_1[range(30,60)] 31 x3_13=x3_1[range(60,90)] 32 x3_14=x3_1[range(90,120)] 33 x3_15=x3_1[range(120,150)] 34 x3_16=x3_1[range(150,180)] 35 x3_17=x3_1[range(180,210)] 36 x3_18=x3_1[range(210,240)] 37 x3_19=x3_1[range(240,270)] 38 x3_110=x3_1[range(270,300)] 39 40 t2=np.arange(12,120,0.4)-12 41 x2_2=odeint(diff_x2,0,t2) 42 x2_22=x2_2[range(0,30)] 43 x2_23=x2_2[range(30,60)] 44 x2_24=x2_2[range(60,90)] 45 x2_25=x2_2[range(90,120)] 46 x2_26=x2_2[range(120,150)] 47 x2_27=x2_2[range(150,180)] 48 x2_28=x2_2[range(180,210)] 49 x2_29=x2_2[range(210,240)] 50 x2_210=x2_2[range(240,270)] 51 52 x3_2=odeint(diff_x3,0,t2) 53 x3_22=x3_2[range(0,30)] 54 x3_23=x3_2[range(30,60)] 55 x3_24=x3_2[range(60,90)] 56 x3_25=x3_2[range(90,120)] 57 x3_26=x3_2[range(120,150)] 58 x3_27=x3_2[range(150,180)] 59 x3_28=x3_2[range(180,210)] 60 x3_29=x3_2[range(210,240)] 61 x3_210=x3_2[range(240,270)] 62 63 t3=np.arange(24,120,0.4)-24 64 x2_3=odeint(diff_x2,0,t3) 65 x2_33=x2_3[range(0,30)] 66 x2_34=x2_3[range(30,60)] 67 x2_35=x2_3[range(60,90)] 68 x2_36=x2_3[range(90,120)] 69 x2_37=x2_3[range(120,150)] 70 x2_38=x2_3[range(150,180)] 71 x2_39=x2_3[range(180,210)] 72 x2_310=x2_3[range(210,240)] 73 74 x3_3=odeint(diff_x3,0,t3) 75 x3_33=x3_3[range(0,30)] 76 x3_34=x3_3[range(30,60)] 77 x3_35=x3_3[range(60,90)] 78 x3_36=x3_3[range(90,120)] 79 x3_37=x3_3[range(120,150)] 80 x3_38=x3_3[range(150,180)] 81 x3_39=x3_3[range(180,210)] 82 x3_310=x3_3[range(210,240)] 83 84 t4=np.arange(36,120,0.4)-36 85 x2_4=odeint(diff_x2,0,t4) 86 x2_44=x2_4[(range(0,30))] 87 x2_45=x2_4[range(30,60)] 88 x2_46=x2_4[range(60,90)] 89 x2_47=x2_4[range(90,120)] 90 x2_48=x2_4[range(120,150)] 91 x2_49=x2_4[range(150,180)] 92 x2_410=x2_4[range(180,210)] 93 94 x3_4=odeint(diff_x3,0,t4) 95 x3_44=x3_4[(range(0,30))] 96 x3_45=x3_4[range(30,60)] 97 x3_46=x3_4[range(60,90)] 98 x3_47=x3_4[range(90,120)] 99 x3_48=x3_4[range(120,150)] 100 x3_49=x3_4[range(150,180)] 101 x3_410=x3_4[range(180,210)] 102 103 t5=np.arange(48,120,0.4)-48 104 x2_5=odeint(diff_x2,0,t5) 105 x2_55=x2_5[range(0,30)] 106 x2_56=x2_5[range(30,60)] 107 x2_57=x2_5[range(60,90)] 108 x2_58=x2_5[range(90,120)] 109 x2_59=x2_5[range(120,150)] 110 x2_510=x2_5[range(150,180)] 111 112 x3_5=odeint(diff_x3,0,t5) 113 x3_55=x3_5[range(0,30)] 114 x3_56=x3_5[range(30,60)] 115 x3_57=x3_5[range(60,90)] 116 x3_58=x3_5[range(90,120)] 117 x3_59=x3_5[range(120,150)] 118 x3_510=x3_5[range(150,180)] 119 120 t6=np.arange(60,120,0.4)-60 121 x2_6=odeint(diff_x2,0,t6) 122 x2_66=x2_6[range(0,30)] 123 x2_67=x2_6[range(30,60)] 124 x2_68=x2_6[range(60,90)] 125 x2_69=x2_6[range(90,120)] 126 x2_610=x2_6[range(120,150)] 127 128 x3_6=odeint(diff_x3,0,t6) 129 x3_66=x3_6[range(0,30)] 130 x3_67=x3_6[range(30,60)] 131 x3_68=x3_6[range(60,90)] 132 x3_69=x3_6[range(90,120)] 133 x3_610=x3_6[range(120,150)] 134 135 t7=np.arange(72,120,0.4)-72 136 x2_7=odeint(diff_x2,0,t7) 137 x2_77=x2_7[range(0,30)] 138 x2_78=x2_7[range(30,60)] 139 x2_79=x2_7[range(60,90)] 140 x2_710=x2_7[range(90,120)] 141 142 x3_7=odeint(diff_x3,0,t7) 143 x3_77=x3_7[range(0,30)] 144 x3_78=x3_7[range(30,60)] 145 x3_79=x3_7[range(60,90)] 146 x3_710=x3_7[range(90,120)] 147 148 t8=np.arange(84,120,0.4)-84 149 x2_8=odeint(diff_x2,0,t8) 150 x2_88=x2_8[range(0,30)] 151 x2_89=x2_8[range(30,60)] 152 x2_810=x2_8[range(60,90)] 153 154 x3_8=odeint(diff_x3,0,t8) 155 x3_88=x3_8[range(0,30)] 156 x3_89=x3_8[range(30,60)] 157 x3_810=x3_8[range(60,90)] 158 159 t9=np.arange(96,120,0.4)-96 160 x2_9=odeint(diff_x2,0,t9) 161 x2_99=x2_9[range(0,30)] 162 x2_910=x2_9[range(30,60)] 163 164 x3_9=odeint(diff_x3,0,t9) 165 x3_99=x3_9[range(0,30)] 166 x3_910=x3_9[range(30,60)] 167 168 t10=np.arange(108,120,0.4)-108 169 x2_10=odeint(diff_x2,0,t7) 170 x2_1010=x2_10[range(0,30)] 171 172 x3_10=odeint(diff_x3,0,t10) 173 x3_1010=x3_10[range(0,30)] 174 175 x21=x2_11 176 x22=x2_12+x2_22 177 x23=x2_13+x2_23+x2_33 178 x24=x2_14+x2_24+x2_34+x2_44 179 x25=x2_15+x2_25+x2_35+x2_45+x2_55 180 x26=x2_16+x2_26+x2_36+x2_46+x2_56+x2_66 181 x27=x2_17+x2_27+x2_37+x2_47+x2_57+x2_67+x2_77 182 x28=x2_18+x2_28+x2_38+x2_48+x2_58+x2_68+x2_78+x2_88 183 x29=x2_19+x2_29+x2_39+x2_49+x2_59+x2_69+x2_79+x2_89+x2_99 184 x210=x2_110+x2_210+x2_310+x2_410+x2_510+x2_610+x2_710+x2_810+x2_910+x2_1010 185 x2=np.concatenate((x21,x22,x23,x24,x25,x26,x27,x28,x29,x210)) 186 187 x31=x3_11 188 x32=x3_12+x3_22 189 x33=x3_13+x3_23+x3_33 190 x34=x3_14+x3_24+x3_34+x3_44 191 x35=x3_15+x3_25+x3_35+x3_45+x3_55 192 x36=x3_16+x3_26+x3_36+x3_46+x3_56+x3_66 193 x37=x3_17+x3_27+x3_37+x3_47+x3_57+x3_67+x3_77 194 x38=x3_18+x3_28+x3_38+x3_48+x3_58+x3_68+x3_78+x3_88 195 x39=x3_19+x3_29+x3_39+x3_49+x3_59+x3_69+x3_79+x3_89+x3_99 196 x310=x3_110+x3_210+x3_310+x3_410+x3_510+x3_610+x3_710+x3_810+x3_910+x3_1010 197 x3=np.concatenate((x31,x32,x33,x34,x35,x36,x37,x38,x39,x310)) 198 199 c2=x2/444*1000 200 c3=x3/444*1000 201 c_Feng2008=c2+c3 202 203 plt.plot(t1,c2,label='RISP') 204 plt.plot(t1,c3,label='9-OH') 205 plt.plot(t1,c_Feng2008,label='RISP+9-OH') 206 plt.legend() 207 plt.title('Feng2008') 208 plt.xlabel('t (h)') 209 plt.ylabel('c (ng/ml)') 210 plt.show()
此段代码过于冗长，有待优化。
图形如下：
model1完成，下面画model2、3、4，步骤几乎完全一样，只是模型稍有不同。
model2的模型和model1完全一致；
model3中，利培酮多了一个周边室，是一个三室模型，所以方程更复杂一点；
model4中吸收部位有一个向9-羟利培酮转化的过程。
分别画出model2、3、4的多剂量给药图形如下：
将4篇文献中的总浓度画在同一张图上，可以比较4个模型的稳态血药浓度范围：
