建模案例-药物中毒的施救方案

建模案例-药物中毒的施救方案

书籍：MATLAB数学建模方法与应用。仅作为学习使用，大家理性学习

本内容涉及高等数学计算，建议大家再学一遍，感触颇深

文章目录

• 建模案例-药物中毒的施救方案
• • 一、问题描述
  • 二、问题分析
  • 三、模型假设与符号约定
  • 四、模型建立
  • 五、模型求解
  • 六、结果与分析
  • 七、施救方案

一、问题描述

一名儿童一次性误服 11 片治疗哮喘病的、计量为 100mg/片的氨茶碱片，到达医院就诊时，距服药时刻已过去两个小时，此时孩子已经出现呕吐、头晕等不良症状。按照药品使用说明书，氨茶碱的成人用量是每次 100～200mg，儿童用量是 3～5mg/kg，如果过量服用，可使血药浓度（单位容积血液中的药量）过高，当血药浓度达到 100μg/ml 时，会出现严重中毒，达到 200μg/ml 则可致命。按接诊医生需要根据专业知识迅速判断孩子的中毒情况，选择合理的施救方案。

二、问题分析

人体服用药物后的吸收代谢过程，药物首先进入胃肠道，再由胃肠道的外壁进入血液循化系统，然后随血液循环被人体脏器代谢而排除体外。血液系统系统对药物的吸收率（药物从胃肠道向血液系统转移的转移率）通常与胃肠道中的药量成正比，血液系统对药物的排除率（药物从血液系统向体外转移的转移率）通常与血液中的药物成正比。

接诊医生需要根据药物的吸收代谢规律及时计算出孩子入院时的血药浓度，依次选择接下来的施救方案。如果血药浓度达到危险的水平，可选的临床施救方案包括：口服活性炭和体外血液透析，前者通过口服活性炭来吸附药物，可使血液中药物的排除率增至原来（人体自身）的 2 倍，后者通过体外对血液进行过滤，可使血液中药物的排除率增至原来的 6 倍。

三、模型假设与符号约定

为了判断孩子的中毒状况，需要建立数学模型，求出胃肠道和血液系统中药量随时间变化的规律。这里把$ t$ 时刻胃肠道中的药量记为 x ( t ) x(t) x(t)， 血液系统中的药量记为 y ( t ) y(t) y(t)，以孩子误服药物的时刻为 t = 0 t=0 t=0时刻。在此基础上给出如下模型假设。

1. 由于药物从胃肠道向血液系统转移的转移率与胃肠道中的药量 x ( t ) x(t) x(t)成正比，这里记比例系统为 λ ( λ > 0 ) \lambda (\lambda \gt 0) λ(λ>0) ，并假设总剂量为 1100mg 的药物在 t = 0 t=0 t=0时刻全部进入胃肠道，即 x ( 0 ) = 1100 x(0)=1100 x(0)=1100。
2. 由于血液系统中的药物排除率与血液系统中的药量 y ( t ) y(t) y(t)成正比，这里记比例系数为 μ ( μ > 0 ) \mu(\mu \gt 0) μ(μ>0)，并假设 t = 0 t=0 t=0时刻血液系统中药量为 0，即 y ( 0 ) = 0 y(0)=0 y(0)=0。
3. 假设任意 t t t时刻，血液系统中的药物分布是均匀的。
4. 在药物总量确定的情况下，血药浓度与血液总量有关。人体血液总量约占体重的 7%~8%(L/kg)，即每千克体重约 70～80ml 血液，一个体重为 50～60kg 的成人的血液总量约为 4000ml 左右。假设误服药儿童的血液总量约为成人的一半，约为 2000ml。
5. 比例系数 λ \lambda λ和 μ \mu μ可由氨茶碱吸收和排除的半衰期确定，从药品说明书可知，氨茶碱吸收的半衰期约为 5h，即 x ( 5 ) = x ( 0 ) / 2 x(5)=x(0)/2 x(5)=x(0)/2，排除的半衰期约为 6h，若考虑血液系统对药物的排除，则有 y ( t + 6 ) = y ( t ) / 2 y(t+6)=y(t)/2 y(t+6)=y(t)/2

四、模型建立

由假设(1)可知， x ( 0 ) = 1000 x(0)=1000 x(0)=1000，随着药物从胃肠道向血液系统的转移，胃肠道中药量 x ( t ) x(t) x(t)会逐渐下降，其下降速度 d x d t \frac{dx}{dt} dtdx​与 x ( t ) x(t) x(t)成正比(比例系数 λ > 0 \lambda \gt 0 λ>0)，从而建立微分方程模型：
d x d t = − λ x , x ( 0 ) = 1100 \frac{dx}{dt} = - \lambda x, x(0)=1100 dtdx​=−λx,x(0)=1100
有假设(2)可知， y ( 0 ) = 0 y(0)=0 y(0)=0，从 t = 0 t=0 t=0开始，药物从胃肠道向血液系统转移(血液系统吸收药物)，同时血液系统启动对药物的排除，血液系统中的药量 y ( t ) y(t) y(t)先增后减。 y ( t ) y(t) y(t)因吸收作用而增长的速度为 λ x \lambda x λx，因排除作用而减少的速度与 y ( t ) y(t) y(t)成正比(比例系数 μ > 0 \mu \gt 0 μ>0)，所以 y ( t ) y(t) y(t)满足微分方程：
d y d t = λ x − μ y , y ( 0 ) = 0 \frac{dy}{dt}=\lambda x - \mu y, y(0)=0 dtdy​=λx−μy,y(0)=0

五、模型求解

求解式(1)和式(2)的 matlab 代码及结果：

>> syms x(t) y(t) lambda mu %定义符号
>> equ = [diff(x) == - lambda * x, diff(y) == lambda * x - mu * y]; 
>> cond = [x(0) == 1100, y(0) == 0];
cond =
[x(0) == 1100, y(0) == 0]
>> [x(t), y(t)] = dsolve(equ, cond); % 求解函数
x(t) =
1100*exp(-lambda*t)
y(t) =
(1100*lambda*exp(-mu*t))/(lambda - mu) - (1100*lambda*exp(-lambda*t))/(lambda - mu)

由以上代码可知：
x ( t ) = 1100 e − λ t x(t) = 1100 e^{-\lambda t} x(t)=1100e−λt

y ( t ) = 1100 λ λ − μ ( e − μ t − e − λ t ) y(t) = \frac{1100 \lambda}{\lambda - \mu} (e^{-\mu t} - e^{-\lambda t}) y(t)=λ−μ1100λ​(e−μt−e−λt)

接下来由氨茶碱吸收和排除的半衰期确定比例系数 λ \lambda λ和 μ \mu μ。由于氨茶碱的吸收半衰期约为 5h，可知 x ( 5 ) = x ( 0 ) / 2 x(5)=x(0)/2 x(5)=x(0)/2，即 110 0 − 5 λ = 1100 / 2 1100^{-5\lambda}=1100/2 1100−5λ=1100/2，解得 λ = ( ln ⁡ 2 ) / 5 ≈ 0.1386 \lambda = (\ln2)/5 \approx 0.1386 λ=(ln2)/5≈0.1386。

为了根据氨茶碱排查的半衰期确定 μ \mu μ，不考虑血液系统对药物的吸收，只考虑血液系统对药物的排除，此时 d y d t = − μ y \frac{dy}{dt} = -\mu y dtdy​=−μy，假设 τ \tau τ时刻血液系统中的药量 y ( τ ) = a y(\tau)= a y(τ)=a，解得 y ( t ) = a e − μ ( t − τ ) , t ≥ τ y(t) = ae^{-\mu(t-\tau)}, t \geq \tau y(t)=ae−μ(t−τ),t≥τ。由于氨茶碱排除的半衰期为 6h，可知 y ( τ + 6 ) = y ( τ ) / 2 y(\tau + 6) = y(\tau)/2 y(τ+6)=y(τ)/2，即 a e − 6 μ = a / 2 ae^{-6\mu} = a/2 ae−6μ=a/2，解得 μ = ( ln ⁡ 2 ) / 6 ≈ 0.1155 \mu = (\ln 2)/6 \approx 0.1155 μ=(ln2)/6≈0.1155。

>> lambda = solve(x(5) == x(0)/2) % 求解血液对药物的吸收率系数
lambda =
log(2)/5
>> syms y2(t) mu tau a
>> y2(t) = dsolve(diff(y2) == -mu * y2, y2(tau) == a) % 求解微积分方程
y2(t) =
a*exp(-mu*t)*exp(mu*tau)
>> mu = solve(y2(tau + 6) == y2(tau)/2, mu) % 求解血液对药物的排除率
mu =
log(2)/6
>> y(t) = subs(y) % 代入求解函数
y(t) =
6600*exp(-(t*log(2))/6) - 6600*exp(-(t*log(2))/5)
>> x = subs(x) %代入求解函数
x(t) =
1100*exp(-(t*log(2))/5)

将 λ \lambda λ和 μ \mu μ的值代入式(3)和(4)可得：
x ( t ) = 1100 e − ( ln ⁡ 2 ) t / 5 x(t) = 1100e^{-(\ln2)t/5} x(t)=1100e−(ln2)t/5

y ( t ) = 6600 ( e − ( ln ⁡ 2 ) t / 6 − e − ( ln ⁡ 2 ) t / 5 ) y(t) = 6600(e^{-(\ln2)t/6}-e^{-(\ln2)t/5}) y(t)=6600(e−(ln2)t/6−e−(ln2)t/5)

六、结果与分析

下面根据模型求解结果绘制 x ( t ) x(t) x(t)和 y ( t ) y(t) y(t)随时间变化的曲线，并计算孩子就诊时 ( t = 2 ) (t=2) (t=2)的血药含量，从而判断孩子的中毒状况。

根据模型假设(4)和问题描述，孩子的血液总量为 2000ml，当血药浓度达到100μg/ml，即血液中药物含量达到 200mg 时，会出现严重中毒，当血药浓度达到 200μg/ml，即血液中药物含量达到 400mg 时则可致命。将 t = 2 t = 2 t=2 代入式(6)可得 y ( 2 ) = 236.5588 m g y(2) = 236.5588mg y(2)=236.5588mg，由此可知孩子就诊时(A点)已经出现严重中毒。结合图可知，若不及时进行施救，血液中药物含量将在一段时间之后超过 400mg，这将危机孩子生命。

进一步地，还可以根据模型式计算血液中药物含量分别达到 200mg 和 400mg 的具体时刻，以及在不进行施救的情况下血液中药物含量达到峰值的时刻和相应的峰值。记血液中药物含量分别达到 200mg、400mg 和峰值的具体时刻分别为 t B t_B tB​、 t C t_C tC​和 t D t_D tD​，则有
y ( t B ) = 200 ， y ( t C ) = 400 ， y ′ ( t D ) = 0 y(t_B) = 200，y(t_C)=400，y'(t_D)=0 y(tB​)=200，y(tC​)=400，y′(tD​)=0

求解式：

>> y2h = double(subs(y(2)))
y2h =
  236.5588
>> tB = double(solve(y(t) == 200))
tB =
    1.6090
>> tC = double(solve(y(t) == 400))
tC =
    4.8648
>> tD = double(solve(diff(y)==0)) % 不施救情形下，血药量达到峰值的时刻
tD =
    7.8910
>> ymax = double(subs(y,tD))
ymax =
  442.0653

由求解结果可知血液中的药量 y = ( t ) y=(t) y=(t)在 t B = 1.6090 t_B=1.6090 tB​=1.6090h达到 200mg，在 t C = 4.8648 t_C=4.8648 tC​=4.8648h（到达医院后 2.8648 h）达到 400mg，在 t D = 7.8910 t_D=7.8910 tD​=7.8910h（到达医院后 5.8910 5.8910 5.8910h）达到峰值 442.0653mg。以上状态分别对应图中 B、C、D 点。

七、施救方案

模型计算结果表明，若不及时进行施救，孩子会有生命危险，因此需要选择合适的施救方案。下面针对口服活性口服碳和体外血液透析两种施救方案，重新计算血液中药量的变化规律。

1. 口服活性炭

假设孩子达到医院时 t = 2 t=2 t=2即开始施救，施救中( t ≥ 2 t\geq2 t≥2)血液里的药量标记为 z ( t ) z(t) z(t)。由于口服活性炭可使血液中的药物排除率增至原来(人体本身)的 2 倍，则血液中药物的排除率与血液中药量的比例系数将变为原来的 2 倍，即 μ = ( ln ⁡ 2 ) / 3 \mu=(\ln2)/3 μ=(ln2)/3，新的模型为
{ d x d t = − λ x , x ( 0 ) = 1100 d z d t = − λ x − μ z , t ≥ 2 , z ( 2 ) = 236.5588 λ = ( ln ⁡ 2 ) 5 , μ = ( ln ⁡ 2 ) 3 \begin{cases} \frac{dx}{dt} = -\lambda x, & x(0)=1100 \\ \frac{dz}{dt} = -\lambda x - \mu z, & t\geq 2,z(2)=236.5588 \\ \lambda = \frac{(\ln 2)}{5}, & \mu= \frac{(\ln2)}{3} \end{cases} ⎩ ⎨ ⎧​dtdx​=−λx,dtdz​=−λx−μz,λ=5(ln2)​,​x(0)=1100t≥2,z(2)=236.5588μ=3(ln2)​​

求解模型式：

>> syms x(t) y(t) z(t) %定义符号函数
>> y(t) = 6600*exp(-(t*log(2))/6) - 6600*exp(-(t*log(2))/5)
>> y = vpa(y,5)
y(t) =
6600.0*exp(-0.11552*t) - 6600.0*exp(-0.13863*t)
>> lambda = log(2) / 5; % 定义 lambda
>> mu = 2 * log(2) / 6; % 定义 mu
>> eq2 = [diff(x) == -lambda * x, diff(z) == lambda * x - mu * z]; % 定义微分方程
>> conf2 = [x(0) == 1100, z(2) == 236.5588]; % 定义初值条件
>> [x(t), z(t)] = dsolve(eq2,conf2); % 求解微分方程
>> z = vpa(z,5)
z(t) =
1650.0*exp(-0.13863*t) - 1609.5*exp(-0.23105*t)

>> figure
>> fplot(x,[0,25],'b--');% 胃肠道中药量变化曲线
>> hold on ;
>> fplot(y,[0,25],'r-'); % 血液中药量变化曲线(不施救)
>> fplot(z,[2,25],'k');%血液中药量变化曲线(施救)
>> grid on;
>> set(gca, 'XMinorGrid', 'on', 'YMinorGrid', 'on')
>> xlabel('时间 t(h)'); ylabel('药量(mg)');
>> legend('x(t)','y(t)','z(t)');

>> t4 = double(solve(diff(z) == 0)); % 施救情况下，血药量达到峰值的时刻
>> t4
t4 =
    5.2582
>> zmax = double(subs(z,t4))
zmax =
  318.3973

由以上代码可知：
z ( t ) = 1650 e − 0.13863 t − 1609.5 e − 0.23105 t , t ≥ 2 z(t) = 1650e^{-0.13863t} - 1609.5e^{-0.23105t}, t \geq 2 z(t)=1650e−0.13863t−1609.5e−0.23105t,t≥2
x ( t ) x(t) x(t)， y ( t ) y(t) y(t)， z ( t ) z(t) z(t)随时间变化的曲线如图，由 z ′ ( t ) = 0 z'(t) = 0 z′(t)=0解得 t = 5.2582 t=5.2582 t=5.2582，相应的 z ( t ) = 318.3973 z(t)=318.3973 z(t)=318.3973。很显然，经口服活性炭施救后，血液中药量 z ( t ) z(t) z(t)在 t = 5.2582 t=5.2582 t=5.2582h 时达到峰值 318.3973mg，远低于 y ( t ) y(t) y(t)的峰值（442.0653mg）和致命水平（400mg）。

2. 体外血液透析

由于血液透析可使血液中药物的排除率增至原来(人体本身)的 6 倍，则血液中药物的排除率与血液中药量的比例系数将变为原来的 6 倍，即 μ = ln ⁡ 2 \mu = \ln 2 μ=ln2，新的模型为：
{ d x d t = − λ x , x ( 0 ) = 1100 d z d t = λ x − μ z , t ≥ 2 , z ( 2 ) = 236.5588 λ = l n 2 5 , μ = ln ⁡ 2 \begin{cases} \frac{dx}{dt} = - \lambda x, & x(0) = 1100 \\ \frac{dz}{dt} = \lambda x - \mu z, & t \geq2,z(2) = 236.5588 \\ \lambda = \frac{ln2}{5}, & \mu = \ln2 \end{cases} ⎩ ⎨ ⎧​dtdx​=−λx,dtdz​=λx−μz,λ=5ln2​,​x(0)=1100t≥2,z(2)=236.5588μ=ln2​
求解模型式：

>> syms x(t) y(t) z(t); % 定义符号
>> lambda = log(2)/5;
>> mu = 6 * log(2)/6;
>> eq2 = [diff(x) == -lambda * x, diff(z) == lambda * x - mu * z]; % 定义微分方程
>> conf2 = [x(0) == 1100, z(2) == 236.5588]; % 定义初值条件
>> [x(t), z(t)] = dsolve(eq2,conf2); % 求解微分方程
>> z = vpa(z,5)
z(t) =
275.0*exp(-0.13863*t) + 112.59*exp(-0.69315*t)
>> y(t) = 6600*exp(-(t*log(2))/6) - 6600*exp(-(t*log(2))/5)
>> y = vpa(y,5)
y(t) =
6600.0*exp(-0.11552*t) - 6600.0*exp(-0.13863*t)

>> figure
>> fplot(x,[0,25],'b--');% 胃肠道中药量变化曲线
>> hold on ;
>> fplot(y,[0,25],'r-'); % 血液中药量变化曲线(不施救)
>> fplot(z,[2,25],'k');%血液中药量变化曲线(施救)
>> grid on;
>> set(gca, 'XMinorGrid', 'on', 'YMinorGrid', 'on')
>> xlabel('时间 t(h)'); ylabel('药量(mg)');
>> legend('x(t)','y(t)','z(t)');

由以上代码可知：
z ( t ) = 275 e − 0.13863 t − 112.59 e − 0.69315 t , t ≥ 2 z(t) = 275e^{-0.13863t} - 112.59e^{-0.69315t}, t \geq 2 z(t)=275e−0.13863t−112.59e−0.69315t,t≥2
x ( t ) x(t) x(t)， y ( t ) y(t) y(t)， z ( t ) z(t) z(t)随时间变化的曲线如图。很显然，经体外血液透析施救后，血液中药量 z ( t ) z(t) z(t)关于时间 t t t单调递减，说明体外血液透析的救治方案是十分有效的。

由以上计算可知两种施救方案均能有效地将血液中的药量降至致命水平以下，虽然体外血液透析可使血液中药量下降更快，但其治疗成本较高，并且存在安全风险，综合考虑，应该选择口服活性炭方案。