药物中毒模型Python实现可视化（数学建模之微分方程模型）

一、模型假设分析：

1.考虑对象：儿童血液总量，胃肠道药物吸收率，血液药物排出率

2.理想条件：

1）儿童血液总量是成人的一半

2）吸收率和转化率均与药量成正比

3.数学描述：列解微分方程，并根据半衰期求出正比的比率，画图分析

二、模型假设

1.胃肠道中药物向血液系统的转移率与药量成正比，且（𝛾 > 0），总剂量𝑘的药物在𝑡 = 0瞬间进入胃肠道。

2.血液系统中药物排除率与药量𝑦(𝑡)成正比，比例系数𝜇（𝜇 > 0），假设𝑡 = 0时血液中无药物。

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

4.孩子的血液总量为2000𝑚𝐿。

三、符号说明

x(t)	t时刻胃肠道中的药量
y(t)	t时刻血液系统中的药量
μ	血液系统中药物排除率与药量的比值
γ	药物向血液系统中转移率与药量的比值
k	药物的总剂量

from sympy import *
import sympy
x = symbols('x',cls=Function)
y = symbols('y',cls=Function)
y1 = symbols('y1',cls=Function)
y2 = symbols('y2',cls=Function)
t = symbols('t')
#λ=0.1386 μ=0.1155

eq_x = Eq(-0.1386*x(t),x(t).diff(t)) #胃肠道中的药量
con_x = {x(0):1100}  #t=0时药量为1100
a=sympy.latex(sympy.Eq(eq_x, x(t), ics=con_x))
x_t = dsolve(eq_x,x(t), ics=con_x).rhs
print('x_t=',x_t)

#x_t = 1100*exp(-0.1386*t)  #胃肠道中的药量函数
eq_y = Eq(0.1386*x_t-0.1155*y(t),y(t).diff(t)) #血液系统中的药量
Ncon_y = {y(0):0} #t=0时肠胃的药量还没有转移到血液系统
y_t = dsolve(eq_y, y(t), ics=con_y).rhs
print('y_t=',y_t)

#y_t = (6600.0 - 6600.0*exp(-0.0231*t))*exp(-0.1155*t)  #血液系统中的药量函数
#到达医院时 t=2h 
y0 = (6600.0 - 6600.0*exp(-0.0231*2))*exp(-0.1155*2)  #送到医院时已转移到血液系统的药量 t=2
eq_y1 = Eq(0.1386*x_t-2*0.1155*y1(t),y1(t).diff(t))  #2倍药物排除率-口服活性炭吸收药物
con_y1 = {y1(2):y0} #t=2时 肠胃的药量转移到血液系统的药量
y1_t = dsolve(eq_y1, y1(t), ics=con_y1).rhs
print(y1_t)

eq_y2 = Eq(0.1386*x_t-6*0.1155*y2(t),y2(t).diff(t)) #6倍药物排除率-体外血液透析
con_y2 = {y2(2):y0} #t=2时 肠胃的药量转移到血液系统的药量
y2_t = dsolve(eq_y2, y2(t), ics=con_y2).rhs
print(y2_t)

import numpy as np
import math
import matplotlib.pyplot as plt
e = math.e
t= np.arange(2,20,0.00001)
x_t = 1100*e**(-0.1386*t)
y_t = (6600.0 - 6600.0*e**(-0.0231*t))*e**(-0.1155*t)
y1_t = (1650.0*e**(0.0924*t) - 1609.49603878015)*e**(-0.231*t)
y2_t = (275.0*e**(0.5544*t) + 112.359036623159)*e**(-0.693*t)
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.title('')
plt.plot(t,x_t,label='胃肠道中的药量')
plt.plot(t,y_t,label='血液系统中的药量')
plt.plot(t,y1_t,label='2倍药物排除率下血液系统的药量')
plt.plot(t,y2_t,label='6倍药物排除率下血液系统的药量')
plt.xlabel('t/h')
plt.ylabel('x,y,y1,y2 /mg')
plt.legend()
plt.show()

from scipy import optimize
#求极值点-不同施救方法下血液系统的药量达到最大时的时间
def y1(t):
    y1 = (1650.0*e**(0.0924*t) - 1609.49603878015)*e**(-0.231*t)
    return y1
def y2(t):
    y2 = (275.0*e**(0.5544*t) + 112.359036623159)*e**(-0.693*t)
    return y2
max_y1 = optimize.fminbound(lambda t: -y1(t), 2,20)
max_y2 = optimize.fminbound(lambda t: -y2(t), 2,20)
print('1.孩子到达医院时刻(t=2)就开始采用口服活性炭吸收药物的办法抢救，\n当血液系统的药量达到最大值的时刻为t={}，此时药量为{}'.format(round(max_y1,2),round(y1(max_y1),2)))
print('2.孩子到达医院时刻(t=2)就开始采用体外血液透析的办法抢救，\n当血液系统的药量达到最大值的时刻为t={}，此时药量为{}'.format(round(max_y2,2),round(y2(max_y2),2)))