其他分享
首页 > 其他分享> > 数学建模（四）—— 局部脑血流测定

数学建模（四）—— 局部脑血流测定

作者：互联网

一、题目要求

       用放射性同位素测定大脑局部血流量的方法如下：由受试者吸入含有某种放射性同位素的气体，然后将探测器置于受试者头部某固定处，定时测量该处的放射性记数率（简称记数率），同时测量他呼出气的记数率。
       由于动脉血将肺部的放射性同位素传送至大脑，使脑部同位素增加，而脑血流又将同位素带离，使同位素减少。实验证明由脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比。其比例系数反应该处的脑血流量，被称为脑血流量系数，只要确定该系数即可推算出脑血流量。动脉血从肺输送同位素至大脑引起脑部记数率上升的速率与当时呼出气的记数率成正比。若某受试者的测试数据如表1所示：

表1 某受试者的测试数据
时间 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00 4.25 4.50 4.75 5.00 5.25 5.50 5.75 6.00 6.25 6.50 6.75 7.00 7.25 7.50 7.75 8.00 8.25 8.50 8.75 9.00 9.25 9.50 9.75 10.00
头部计数率 1534 1528 1468 1378 1272 1162 1052 947 848 757 674 599 531 471 417 369 326 288 255 225 199 175 155 137 121 107 94 83 73 65 57 50 44 39 35 31 27
呼出气记数率 2231 1534 1054 724 498 342 235 162 111 76 52 36 25 17 12 8 6 4 3 2 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0

       根据以上题目所给的条件及数据：1.建立确定脑部血流系数的数学模型；2.计算上述受试者的脑血流系数。

二、模型的建立与求解

2.1 模型的建立

       由题意可知脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比，动脉血从肺输送同位素至大脑引起脑部记数率上升的速率与当时呼出气的记数率成正比，因此把脑部计数率 h ( t ) h(t) h(t)作为讨论对象。
       设 t t t 时刻脑部计数率为 h ( t ) h(t) h(t)， △ t △t △t时刻以后脑部计数率为 h ( t + △ t ) h(t+△t) h(t+△t)，由问题假设可知，脑部计数率的增量 △ h △h △h只与以下三个因素有关。
       1.肺动脉血将肺部的放射性同位素送到大脑，使脑部记数率增加 △ h 1 △h_{1} △h1​；
       2.脑血流将同位素带离，脑记数率下降 △ h 2 △ h2 △h2；
       3.放射性元素自身有衰减，设其半衰期为 τ \tau τ，由此引起的记数率下降为 △ h 3 △h3 △h3。
       又由医学实验和假定有 d h 1 d t = k p ( t ) ， d h 2 d t = K h ( t ) \frac{dh_{1}}{dt}=kp(t)，\frac{dh_{2}}{dt}=Kh(t) dtdh1​​=kp(t)，dtdh2​​=Kh(t)

       查阅资料得由放射线元素衰减引起的脑部计数率下降为 Δ h 3 = − h ( t ) ( 1 2 ) Δ t τ × l n 2 × 1 τ × Δ t + o ( Δ t ) \Delta h_{3}=-h(t)(\frac{1}{2})^{\frac{\Delta t}{\tau }}\times ln2\times \frac{1}{\tau }\times \Delta t+o(\Delta t) Δh3​=−h(t)(21​)τΔt​×ln2×τ1​×Δt+o(Δt)

略去高次项，变化率为 d h 3 d t = − l n 2 τ h ( t ) \frac{dh_{3}}{dt}=-{\frac{ln2}{\tau }}h(t) dtdh3​​=−τln2​h(t)

       因此， △ t △t △t时刻内脑部放射性元素计数率有变化，有 Δ h ( t ) = Δ h 1 ( t ) − Δ h 2 ( t ) + Δ h 3 ( t ) \Delta h(t)=\Delta h_{1}(t)-\Delta h_{2}(t)+\Delta h_{3}(t) Δh(t)=Δh1​(t)−Δh2​(t)+Δh3​(t)

变化率为 d h d t = d h 1 ( t ) d t − d h 2 ( t ) d t + d h 3 ( t ) d t \frac{dh}{dt}=\frac{dh_{1}(t)}{dt}-\frac{dh_{2}(t)}{dt}+\frac{dh_{3}(t)}{dt} dtdh​=dtdh1​(t)​−dtdh2​(t)​+dtdh3​(t)​

可以得到 d h d t = k p ( t ) − K h ( t ) − l n 2 τ h 3 \frac{dh}{dt}=kp(t)-Kh(t)-\frac{ln2}{\tau }h_{3} dtdh​=kp(t)−Kh(t)−τln2​h3​

       由于放射性同位素的半衰期一般要远大于两次检测之间的间隔时间，所以由放射线元素衰减引起的脑部计数率下降可忽略不计，因此上式可写为 d h d t = k p ( t ) − K h ( t ) \frac{dh}{dt}=kp(t)-Kh(t) dtdh​=kp(t)−Kh(t)

       要确定脑血流量系数模型必须对 h ( t ) h(t) h(t)和 p ( t ) p(t) p(t)的实验数据进行分析，用MATLAB画出 h ( t ) h(t) h(t)和 p ( t ) p(t) p(t)随时间变化的散点图并观察其变化趋势。 h ( t ) h(t) h(t)和 p ( t ) p(t) p(t)的散点图如图1、图2所示。

在这里插入图片描述

图1  脑部计数率散点图
在这里插入图片描述
图2  呼出计数率散点图

       由图2可以看出，呼出计数率与时间有着近似负指数的关系，而图1中脑部计数率与时间的关系不容易观察得到。因此，设 p ( t ) = a e b t p(t)=ae^{bt} p(t)=aebt，用MATLAB的cftool工具箱拟合可得参数a=10000，b=-1.5，即 p ( t ) = 10000 e − 1.5 t p(t)=10000e^{-1.5t} p(t)=10000e−1.5t。呼出计数率的拟合曲线如图3所示。

在这里插入图片描述
在这里插入图片描述

图3  呼出计数率的拟合曲线

       假设在吸入气体瞬时，脑中放射物记数率为零，即0时刻时，脑部计数率 h ( 0 ) = 0 h(0)=0 h(0)=0，可得到一个带初始条件的微分方程 { d h d t = 10000 k e − 1.5 t − K h ( t ) h ( 0 ) = 0                                                   \left\{\begin{matrix} \frac{dh}{dt}=10000ke^{-1.5t}-Kh(t)\\ h(0)=0\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{}\ _{} \end{matrix}\right. {dtdh​=10000ke−1.5t−Kh(t)h(0)=0 ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​​

求解该微分方程就可以得到脑部血流系数的数学模型为 h ( t ) = 10000 k K − 1.5 ( e − 1.5 t − e − K t ) h(t)=\frac{10000k}{K-1.5}(e^{-1.5t}-e^{-Kt}) h(t)=K−1.510000k​(e−1.5t−e−Kt)

2.2 模型的求解

       为了更好的评价两种求解结果的精确程度，我们定义一个平均绝对误差 e = ∑ i = 1 n ∣ h ( i ) − h 0 ( i ) ∣ n e=\frac{\sum_{i=1}^{n}\left | h(i)-h_{0}(i) \right |}{n} e=n∑i=1n​∣h(i)−h0​(i)∣​

其中 n n n为实验的测试次数，本题目中 n = 37 n=37 n=37， h ( i ) h(i) h(i)和 h 0 ( i ) h_{0}(i) h0​(i)分别为第 i i i次模型的预测值和实际测试值。

2.2.1 离散方程对血流系数求解

       将 d h d t = k p ( t ) − K h ( t ) \frac{dh}{dt}=kp(t)-Kh(t) dtdh​=kp(t)−Kh(t)离散化，记时间间隔为T，由前插公式可得 h n + 1 − h n T = − K h n + k p n \frac{h_{n+1}-h_{n}}{T}=-Kh_{n}+kp_{n} Thn+1​−hn​​=−Khn​+kpn​

或 h n + 1 = ( 1 − K T h n ) + k p n h_{n+1}=(1-KTh_{n})+kp_{n} hn+1​=(1−KThn​)+kpn​

其中 h n = h ( t 0 + n T ) h_{n}=h(t_{0}+nT) hn​=h(t0​+nT)， p n = p ( t 0 + n T ) p_{n}=p(t_{0}+nT) pn​=p(t0​+nT)

       对于离散方程 h n + 1 = ( 1 − K T h n ) + k p n h_{n+1}=(1-KTh_{n})+kp_{n} hn+1​=(1−KThn​)+kpn​，可以通过联立不同时刻的方程组求得一系列 K K K值，但是由于在实际测试中存在随机误差，以及离散化的截断误差，使得这些 K K K值不尽相同。为了充分利用已测数据，可以利用MATLAB的cftool工具箱对参数进行拟合，拟合后的离散方程为 h n + 1 = − 0.8825 h n + 0.078075 p n h_{n+1}=-0.8825h_{n}+0.078075p_{n} hn+1​=−0.8825hn​+0.078075pn​

在这里我们取 t 0 = 1 t_{0}=1 t0​=1， T = 0.25 T=0.25 T=0.25，可求得 K = 0.47 K=0.47 K=0.47， k = 0.3123 k=0.3123 k=0.3123。
       把求得的 K K K和 k k k的值代入 h ( t ) = 10000 k K − 1.5 ( e − 1.5 t − e − K t ) h(t)=\frac{10000k}{K-1.5}(e^{-1.5t}-e^{-Kt}) h(t)=K−1.510000k​(e−1.5t−e−Kt)即可得到 h ( t ) h(t) h(t)的函数，进而求得模型的预测值，再和真实值进行比较，用MATLAB绘出每一个预测值与真实值的差值关系图，如图4所示。在这里插入图片描述

图4  离散方程求解误差图

       离散方程求解得 K = 0.47 K=0.47 K=0.47， k = 0.3123 k=0.3123 k=0.3123，此时模型的平均绝对误差为 e = 77.987362 e=77.987362 e=77.987362。
       用离散方程求取血流系数，该算法的优点是计算简单，且容易理解。缺点是实验中的测量间隔时间，即公式中的步长T，对计算结果会产生比较大的影响，通常步长越小所求的系数越精确，但如果增加实验的测量次数同时也会导致实验的难度和成本同时上升。

2.2.2 微分方程对血流系数求解

       对于参数 K K K和 k k k的求解还可以将脑部血流系数数学模型，即 h ( t ) = 10000 k K − 1.5 ( e − 1.5 t − e − K t ) h(t)=\frac{10000k}{K-1.5}(e^{-1.5t}-e^{-Kt}) h(t)=K−1.510000k​(e−1.5t−e−Kt)，直接利用MATLAB里的cftool工具箱进行拟合，拟合得 K = 0.5 K=0.5 K=0.5， k = 0.4001 k=0.4001 k=0.4001。

在这里插入图片描述
       把求得的和的值代入脑部血流系数数学模型即可得到 h ( t ) h(t) h(t)的函数，可求得模型的预测值，再和真实值进行比较，用MATLAB绘出每一个预测值与真实值的差值关系图，如图5所示。在这里插入图片描述

图5  微分方程求解误差图

       微分方程求解得 K = 0.5 K=0.5 K=0.5， k = 0.4001 k=0.4001 k=0.4001，此时模型的平均绝对误差为 e = 0.221975 e= 0.221975 e=0.221975。

三、结果分析

在这里插入图片描述

图6  离散方程求解模型预测值与真实值对比
在这里插入图片描述
图7  微分方程求解模型预测值与真实值对比

       通过图6与图7的对比我们可以发现，用离散方程对模型进行求解所得到的预测值的收敛的相对较慢，在初期产生的误差较大，从而导致误差的平均值也相对较大。而使用微分方程求解的模型所得到的预测值与实验测试数据基本重合，平均误差也非常小。因此，使用 K = 0.5 K=0.5 K=0.5作为脑血流系数去推算出脑血流量是更为合理的。

标签：1.5,记数,血流,frac,求解,建模,脑部,测定,计数率
来源： https://blog.51cto.com/u_15178976/2790873