作者:欧新宇(Xinyu OU)
当前版本:Release v1.0
开发平台:Python3.11
运行环境:Intel Core i7-7700K CPU 4.2GHz, nVidia GeForce GTX 1080 Ti
本教案所涉及的数据集仅用于教学和交流使用,请勿用作商用。
最后更新:2024年3月10日
微分方程模型
一天夜晚,你作为见习医生正在医院内科急诊室值班,两位家长带着一个孩子急匆匆进来,诉说两个小时前,孩子一口气误吞下了11片治疗哮喘病的、剂量为每片100mg的氨茶碱片,已经出现呕吐、头晕等不良症状。按照药瓶使用说明书,氨茶碱的每次用量成人是 ,儿童是 。如果服用过量,可使血药浓度(单位血液容积中的药量)过高,当血药浓度达到 时,会出现严重中毒,达到 则会致命。
作为一名医生,你清楚地知道,由于孩子服药是在两小时前,现在药物已经从胃进入肠道,无法再用刺激呕吐的办法排除。当前需要作出判断的是,孩子的血药浓度会不会达到 ,甚至是 ,如果会达到,则临床上应采取紧急方案来就救治孩子。
类似问题:
1. 人体血液的基本情况 人体服用一定药物后,血液浓度与人体的血液总量有关。一般来说,血液总量约位体重的 ,即体重 的成年人血液总量约为 4000mL。目测这个孩子的体重约为成年人的一半(提出假设或称量确定),可认为其血液总量约为 。由此,血液系统中的血药浓度与药量之间可以互相转换。
2. 药物转移的过程 药物口服后迅速进入胃肠道,再由胃肠道的外壁进入血液循环系统,被血液吸收。胃肠道中药物的转移率,即血液系统的吸收率,一般与胃肠道中的药量(设为 )成正比。药物在血液吸收的同时,又通过代谢作用由肾脏排出体外,排除率一般与血液中的药量(设为 )成正比。如果认为整个血液系统内药物的分布,即血药浓度是均匀的,可以将血液系统看作是一个房室,建立所谓的一室模型。基于以上思想,我们可以得到药物转移过程的示意图如下:
血液系统对药物的吸收率和排除率可以由半衰期确定,从药品的说明书可知,氨茶碱吸收的半衰期为 ,排除的半衰期为 。(PS:半衰期是指某种物质衰变至原来数量的一半所需的时间。例如,在放射性衰变中,半衰期是指放射性物质衰变至原来数量的一半所需的时间;在医学领域,药物的半衰期指的是药物在血液中浓度降低一半所需的时间,或者药物自体内消除一半所需的时间。这个指标有助于医生确定药物的给药频率和剂量,以及计算多剂量给药中到达稳态浓度所需的时间。)
3. 施救方案 如果血液浓度达到危险水平,临床上施救的一种方法是采用口服活性炭来吸附药物,可使药物的排除率增加到原来(人体自身)的 倍,另一种办法是进行体外血液透析,药物排除率可增加到原来的 倍,但是安全性不能得到充分保证,建议尽量少用。
为了判断孩子的血药浓度会不会达到危险水平,需要寻求胃肠道和血液系统中的药量随时间变化的规律。记胃肠道中的药量为 ,血液系统中的药量为 ,时间 以孩子误服药的时刻为起点 。根据前面的调查与分析可作以下假设:
根据假设对胃肠道中药量 和血液系统中的药量 建立如下模型:
方程(1)(2)中的参数 分别表示药物在胃肠道和血液系统中的转移速率,它们可以由假设3中药物的半衰期确定。
微分方程(1)是可分离变量方程,通过求解可以得到:
其中胃肠道中的药量 随时间 单调减少并趋于0。为了确定超参数 ,我们可以利用药物吸收的半衰期为 5h,即 ,也即 。由此可得 。此时,公式(3)可简化为:
下面给出Python求解微分方程方程(1)得到解(3)的Python代码。此处,我们使用sympy符号计算方法进行化简,然后代入初值 ,以获得公式(3)。
from sympy import symbols, Eq, dsolve, Function
# 1. 定义变量
t = symbols('t')
lamb = symbols('lambda')
# 2. 定义函数
x = Function('x')(t)
# 3. 建立微分方程
diff_eq = Eq(x.diff(t), -lamb*x) # dx/dt = -lambda*x
# 4. 使用dsolve解微分方程,并按照初始条件进行化简
solution = dsolve(diff_eq, x)
solution_with_condition = dsolve(diff_eq, x, ics={x.subs(t, 0): 1100})
# 5. 输出求解结果打印解
print('solution = ', solution)
print('solution_with_condition = ', solution_with_condition)
solution = Eq(x(t), C1*exp(-lambda*t))
solution_with_condition = Eq(x(t), 1100*exp(-lambda*t))
整理一下以上代码的输出,可得:当初值 时,解得公式(3):
将(3)代入方程(2),可以得到描述药物在血液系统中的一阶线性微分方程:
求解可得:
类似的,使用sympy符号计算方法进行化简,然后代入初值 ,求解获得公式(5)。
from sympy import symbols, Eq, dsolve, simplify, Function
from sympy import exp
# 1. 定义变量
t = symbols('t')
lamb = symbols('lambda')
mu = symbols('mu')
# 2. 定义函数
y = Function('y')(t)
# 3. 建立微分方程
diff_eq = Eq(y.diff(t), - mu*y + 1100 * lamb * exp(-lamb * t)) # dy/dt
# 4. 使用dsolve解微分方程,并按照初始条件进行化简
solution = dsolve(diff_eq, y)
solution_with_condition = dsolve(diff_eq, y, ics={y.subs(t, 0): 0})
# solution_with_condition = simplify(solution_with_condition)
# 5. 输出求解结果打印解
print('solution = ', solution)
print('solution_with_condition = ', solution_with_condition)
solution = Eq(y(t), C1*exp(-mu*t) - 1100*lambda*exp(-lambda*t)/(lambda - mu))
solution_with_condition = Eq(y(t), 1100*lambda*exp(-mu*t)/(lambda - mu) - 1100*lambda*exp(-lambda*t)/(lambda - mu))
Note
在上面的例子中,我们可以尝试使用sympy.simplify()对解进行化简。但应注意化简的结果是否符合习惯,如果不符合,可以考虑适当结合手工化简。例如在本例中,若调用simplify函数后,化简结果等于:
这显然不符合习惯,因此我们可以使用直接求解的结果,并作适当的手工化简以得到最终的解:
为了根据药物排除的半衰期为 6h 来确定排除率比例系数 ,可以考虑血液系统只对药物进行排除的情况,这时 满足方程:
设在某时刻 ,有 ,即此时血液中的药物含量为 。则:
利用 ,可得 。
求解方程(7)的Python代码如下所示:
from sympy import symbols, Eq, dsolve, simplify, Function
from sympy import exp
import numpy as np
# 1. 定义变量
t = symbols('t')
lamb = symbols('lambda')
mu = symbols('mu')
tau = symbols('tau')
# 2. 定义函数
y = Function('y')(t)
# 3. 建立微分方程
diff_eq = Eq(y.diff(t), - mu*y) # dy/dt
# 4. 使用dsolve解微分方程,并按照初始条件进行化简
solution = dsolve(diff_eq, y)
solution_with_condition = dsolve(diff_eq, y, ics={y.subs(t, 'tau'): 'a'})
solution_with_condition = simplify(solution_with_condition)
# 5. 输出求解结果打印解
print('solution_with_condition = ', solution_with_condition)
solution_with_condition = Eq(y(t), a*exp(mu*(-t + tau)))
以上根据以上方程可得简化后的血液药物浓度微分方程,适当简化后可得解(8):
将 和 代入方程(3),(6)中,可得:
以上公式可以表明血液系统中的药量 随时间先增后减,并在 小时达到最大值,然后逐渐减小。
首先,我们使用Python根据公式(9)和(10)绘制出胃肠道中的药量 和血液系统中的药量 随时间的变化曲线,如 图1 所示。
import numpy as np
import matplotlib.pylab as plt
lamb = 0.1386
mu = 0.1155
xt = lambda t: 1100 * np.exp(-lamb * t)
yt = lambda t: 6600 * (np.exp(-mu * t) - np.exp(-lamb * t))
t0 = np.linspace(0, 20, 100)
plt.figure(figsize=(10,5))
plt.rc('font', size=16)
plt.plot(t0, xt(t0), '--', label='$x(t)$: Drag in the gut')
plt.plot(t0, yt(t0), label = '$y(t)$: Drug in the blood')
plt.legend()
plt.grid()
plt.xlabel('t/h')
plt.ylabel('x,y/mg')
plt.title('figure 1: No treatment')
plt.show()
根据假设4,孩子的血液总量为 2000mL,出现严重中毒的血液浓度 和致命的血液浓度 分别相当于血液中药量 达到 200mg 和 400mg。由 图1 看出,药量 y(t) 在约2h达到200mg,即孩子到达医院时已经出现严重中毒,如不及时施救,药量 y(t) 将在5小时(到达医院3h)达到致命剂量400mg,即孩子已经死亡。
由公式(10)容易精确地计算出孩子到达医院时血液系统中的药量 y(2) = 236.5mg,而计算药量达到400mg的时间(记作),则需要解非线性方程 ,用Python编程可得到 。
由 图1 还可以看出,血液中药量 达到最大值的时间约在 ,即到达医院后6h,其精确值可也由方程(10)计算获得,记作 :
此时,
from scipy.optimize import fsolve, fmin
y2 = yt(2)
print('2h时血液中的药量:', y2)
eq = lambda t: yt(t) - 400
t1 = fsolve(eq, 2)
print('药量在 {}h 时达到400mg。'.format(t1))
ytm = lambda t: -yt(t) # 定义yt的相反数匿名函数
t2 = fmin(ytm, 2)
print('{} h 时达到最大药量 {}。'.format(t2, yt(t2)))
2h时血液中的药量: 236.52132678290988
药量在 [4.86585945]h 时达到400mg。
Optimization terminated successfully.
Current function value: -442.065329
Iterations: 21
Function evaluations: 42
[7.89267578] h 时达到最大药量 [442.06532921]。
Note
在上面的代码中,我们定义了一个 的相反数匿名函数,用于计算函数 的最大值。这样做的主要原因是在 scipy.iptimize 优化库中没有 fmax
函数,因为优化库主要关注于寻找函数的最小值,而不是最大值。然而,找到函数的最大值是一个常见的问题,通常可以通过将目标函数取负,然后寻找最小值来解决。
利用这个模型,还可以确定对于孩子及成人服用氨茶碱能引起严重中毒和致命的最小剂量问题(拓展练习)。
根据上述分析,我们得出结论,若不即使施救,孩子会有生命危险。根据调查,如采用口服活性炭来吸附药物的办法进行施救,药物的排除率可增加到 的2倍,即 ;而采用体外血液透析的方法,药物排除率可增加到 ,血液中药量下降得更快。虽然血液透析能够更快排除药物,但该方法存在一定的危险性,而口服活性炭则相对简单且安全。在临床中究竟采取哪种方法进行救治需要我们给出响应的数据分析后,征求病人和家属的意见后再确定。
下面,我们给出使用口服活性炭进行施救的方案分析。
设孩子到达医院时(t=2)就立即开始施救。前面已经算出 ,由(2)(3)可得到新模型:
其中, 表示采用口服活性炭后,血液中的药量。不难发现,该模型依然是一个一阶线性微分方程,只不过初始时刻 时,当 ,而 时,该方程(11)的解为:
利用初值 ,可以得到解的简化表达式:
t = symbols('t')
lamb = 0.1386
mu2 = 0.2310
z = Function('z')(t)
diff_eq2 = Eq(z.diff(t), -mu2*z + 1100 * lamb * exp(-lamb * t)) #dz/dt
solution2 = dsolve(diff_eq2, z)
solution2_with_condition = dsolve(diff_eq2, z, ics={z.subs(t, 2): 236.5})
# 5. 输出求解结果打印解
print('solution = ', solution2)
print('solution_with_condition = ', solution2_with_condition)
solution = Eq(z(t), C1*exp(-0.231*t) + 1650.0*exp(-0.1386*t))
solution_with_condition = Eq(z(t), -1609.52988961616*exp(-0.231*t) + 1650.0*exp(-0.1386*t))
使用Python对公式(13)作图得 图2。
zt = lambda t: -1609.5 * np.exp(-0.231 * t) + 1650 * np.exp(-0.1386 * t)
t2 = np.linspace(2, 20, 100)
plt.figure(figsize=(10,5))
plt.rc('font', size=16)
plt.plot(t0, xt(t0), '--', label='$x(t)$: Drag in the gut')
plt.plot(t2, zt(t2), label = '$z(t)$: Drug in the blood')
plt.legend()
plt.grid()
plt.xlabel('t/h')
plt.ylabel('x, z/mg')
plt.title('figure 2: Activated carborn')
plt.show()
由图2可以看出,施救后血液中药量 达到最大值的时间约在 ,即到达医院施救后3h,其精确值可由公式(8)计算获得,记作 。,且 ,远低于 的最大值和致命水平。
图2还表明,虽然采用了口服活性炭吸附药物的办法进行施救,但血液中药量 仍有一段时间在上升,说明用这种方法药物的排除率增加还不够大。不妨计算一下,如果要使 在施救后 () 立即下降,排除率 至少应该多大。
z(t) 在 取得最大值,相当于 时,公式(7)满足:
由(5)算出 ,再利用前面已有的 和 ,立即得到 ,约原来(人体自身 )的 倍。此时,如果采用血液透析法,由于它可使排除率增加到6倍,它确实可以有效抑制药物中毒。但在临床中,血液透析的安全性依然是需要考虑的重要因素之一。
from scipy.optimize import fsolve, fmin
ztm = lambda t: -zt(t) # 定义yt的相反数匿名函数
t3 = fmin(ztm, 2)
print('{} h 时达到最大药量 {} mg。'.format(np.round(t3[0], 4), np.round(zt(t3)[0], 4)))
x2 = xt(2)
mu3 = x2*lamb/y2
n = mu3/mu
print('药量不增加时,排除率达到 mu = {:.4f} 时,才能快速排除危险,\n 此时排除率为人体自身排除率的 {:.4f} 倍。'.format(mu3, mu3/mu))
Optimization terminated successfully.
Current function value: -318.390614
Iterations: 19
Function evaluations: 38
5.2595 h 时达到最大药量 318.3906 mg。
药量不增加时,排除率达到 mu = 0.4885 时,才能快速排除危险,
此时排除率为人体自身排除率的 4.2298 倍。
...请读者自行补充完成