数学模型——药物中毒急救模型（基于python）_中毒急救微分方程建模-CSDN博客

本文链接：https://blog.csdn.net/qq_60471758/article/details/123481806

背景

人体服用一定量的药物后，血药浓度与人体的血液总量有关，一般来说，血液总量约为体重的7%~8%，即体重在50~60kg的成年人血液总量约为4000mL。

这里需要考虑，药物在进入胃肠道后，被血液吸收，胃肠道中药物的转移率，即血液系统的吸收率，一般与胃肠道的药量成正比。但药物在被血液吸收时，会通过代谢作用由肾脏排出体外，排除率与血液中药量成正比。

假设整个血液系统中药物的分布，即血药浓度是均匀的，将血液系统看作一个房室，建立一室模型。

首先确定吸收率和排除率，设定本次药物为氨茶碱，其吸收的半衰期为5h，排除的半衰期为6h（不同药物的吸收和排除通过上网查证或者药物说明书获取）

治疗方法：当血药浓度达到危险水平，可以采用口服活性炭来吸附药物，可使药物的排除率增加为原来的2倍；亦或者进行体外血液透析，药物的排除率可达到原来的6倍，但安全性不能得到充分保证。

模型建立

为了判断孩子的血药浓度会不会达到危险的水平，需要寻求胃肠道和血液系统中药量随时间变化的规律。记胃肠道的药量为 $x(t)$ ，血液系统中的药量为 $y(t)$ ，时间t以孩子服药的时刻为起点。

假设

1.胃肠道中药物向血液系统的转移率与药量 $x(t)$ 成正比，比例系数记为 $\lambda$ (>0)，总剂量为1100mg的药物在t=0时刻瞬间进入胃肠道

2.血液系统中药物的排除率与药量 $y(t)$ 成正比，比例系数记为u(>0)，t=0时血液中无药物

3.氨茶碱吸收的半衰期为5h，排除的半衰期为6h

4.孩子血液总量为2000mL

求解模型用到的python库

import matplotlib.pyplot as plt
from scipy.integrate import odeint
import numpy as np
from sympy import *

由假设1可得，x(0)=1100mg，随着药物从胃肠道向血液系统的转移，x(t)下降 $\frac{dy}{dt}=\lambda x-uy,y(0)=0$ 的速度与x(t)本身成正比，故x(t)满足：

$\frac{dx}{dt}=-\lambda x,x(0)=1100$

对应的代码为

def dx_dt(x,t):
    p = 0.1386  # 比例系数
    return np.array(-p*x+0*t)

由假设2可得，y(0)=0,药物从胃肠道向血液系统的转移相当于血液系统对药物的吸收，y(t)由于吸收作用增长的速度为 $\lambda x$ ，由于排除而减少的速度与y(t)本身成正比，故y(t)满足

$\frac{dy}{dt}=\lambda x-uy,y(0)=0$

对应的代码为

def dy_dt(y,t):
    p = 0.1386
    u = 0.1155
    x = np.exp(-p*t)*1100  # x关于t的函数，通过Sympy求出
    return np.array(p*x-u*y+0*t)

求解x关于t的方程代码

t = Symbol('t')
x = Function('x')(t)
end = dsolve(Eq(diff(x,t,1),-0.1386*x))
print(end)

对模型进行求解，得出胃肠道中药量和血液系统中药量的的变化趋势，相应代码为：

t = np.arange(0,25,1)  # 设置时间序列
a = odeint(dx_dt,1100,t)  # 求解方程
b = odeint(dy_dt,0,t)
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.plot(a,c='r')  # 对应x(t)
plt.plot(b,c='b')  # 对应y(t)
plt.ylim(0,1200)
plt.grid(ls='--',color='k',alpha=0.5)
plt.show()

由图可知，根据假设4，出现严重中毒时的血药浓度100ug/mL和致命的血药浓度200ug/mL，药量在约2h达到200mg，即到达医院时已经出现严重中毒，需要立刻进行补救措施

补救措施

使用口服活性炭进行救助。假设孩子在2h就开始施救，此时y(2)=236.5，此时，新的模型为

$\frac{dz}{dt}=\lambda x-uz,t\geqslant 2$

$x=1100e^{-\lambda t}$

对应代码为

def dz_dt(z,t):
    p = 0.1386
    x = np.exp(-p*t)*1100
    u = 0.2310
    return np.array(p*x-u*z)

T = np.arange(0,25,1)
c = odeint(dz_dt,0,T)
plt.plot(c,c='y')
plt.show()