有些葯物本身不具有葯理活性，進入體內后經代謝產生活性物質，也有一些葯物，其自身與代謝物都具有葯理活性，對於此類葯物，我們需要關注它們的代謝物的濃度。
利用隔室模型進行研究時，步驟是根據隔室模型寫出微分方程組，求出表達式，畫圖。對於常見的一室、二室、三室模型，葯代動力學課本上都給出了血葯濃度公式，但是有了代謝物之后，就不是應用上述模型，換句話說，沒有課本上現成的公式可以套用了。
那么如何畫圖呢？我們一般會想到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個模型的穩態血葯濃度范圍：
