局部脑血流的测定
摘要
随着医学的快速发展，脑部的研究越来越得到人们的重视。而对脑血流系数的精确测定可以帮助人们快速得到相关指标，所以对医学领域来说很有意义。
问题1，首先我们根据题设中的信息及放射性元素衰减性，列出头部记数率和呼出记数率关于脑部血流系数的微分表达式；再分析试验得到的呼出气记数率的数据，用Matlab的拟合工具箱拟合得到关于呼出气记数率的函数；最后代入上述微分表达式求解微分方程即可得到关于脑部血流系数的表达式，表达式为)(5.11000)(5.1KtteeKktN。
问题2，先用Matlab拟合工具箱求出头部记数率的函数，用对比系数法可得到K和k的近似值，即3977.0,5015.0kK，但这种方法是不精确的，只是用于后面方法得到参数的验证；再利用最小二乘法求解拟合后曲线参数的函数以及由问题1得到的关于脑部血流系数的表达式，即可得到脑血流系数：0.5000K，4001.0k，最后对得到的值进行误差分析，可知脑血流的预测值和实际值很吻合，比较符合题意。
关键词：对比系数法 放射性元素衰减 曲线拟合
2
一 问题重述
用放射性同位素测定大脑局部脑血流量的方法如下：由受试者吸入含有某种放射性同位素的气体，然后将探测器置于受试者头部某固定处，定时测量该处的放射性记数率(简称记数率)，同时测量他呼出气的记数率。
由于动脉血将肺部的放射性同位素传送至大脑，使脑部同位素增加，而脑血流又将同位素带离，使同位素减少。实验证明由脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比，其比例系数反映了该处的脑血流量，被称为脑血流量系数，只要确定该系数即可推算出脑血流量。动脉血从肺输送同位素至大脑引起脑部记数率上升的速度与当时呼出气的记数率成正比。若某受试者的测试数据如附表1所示：
根据以上题目所给的条件及数据，回答以下问题：
1、建立确定脑部血流系数的数学模型；
2、计算上述受试者的脑血流系数。
二 模型假设
1.脑部记数率的上升只与从肺部输送的放射性同位素有关；
2.脑部记数率的下降只与当时该处的脑血流量有关；
3.脑血流量在测定期间恒定，心脏博动、被测试者大脑活动、情感波动等带来的变化忽略不予考虑；
4.每次仪器测量为相互独立事件，各测量值无记忆相关；
5.在吸入气体瞬时，脑中放射物记数率为零；
6.脑血流量与脑血流量系数成单值函数关系，求得后者即可确定前者。
三 符号说明
符号 表示意义 3 t 时间
k 脑部记数率上升的速率与呼出气记数率的比例系数
K 脑部记数率下降的速率与当时该处脑部记数率的比例
系数，即脑血流量系数
)(tN 脑部记数率
)(tP 呼出气记数率
1N 脑部记数率的增量
2N 脑部记数率的减量
3N 放射性衰减引起的脑部记数率的减量
四 问题分析
问题1，首先根据题设可知：一方面，由脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比, 并且其比例系数反映了该处的脑血流量；另一方面，动脉血从肺输送同位素至大脑引起脑部记数率上升的速率与当时呼出气的记数率成正比，由上述两方面可得到头部记数率关于脑部血流系数和呼出气记数率的表达式；再分析试验得到的呼出气记数率的数据，用MATLAB的拟合工具箱拟合得到关于呼出气记数率的函数；最后代入上述表达式求解微分方程即可得到关于脑部血流系数的表达式。
问题2，利用Matlab里专门求解拟合后曲线参数的函数，再由问题1得到的关于脑部血流系数的表达式即可得到脑血流系数。最后再对得到的值进行误差分析。
五 模型的建立与求解
5.1求解脑血流系数的数学模型
设某时刻0t时，脑部记数率为)(tN，在t时刻后记数率为)(ttN， 4 由题设及基本假设1和2可知，脑部记数率的增量)()(tNttNN只与下面两个因素有关：
(1)动脉血从肺部输送放射性同位素至大脑引起脑部记数率的增量为1N；
(2)脑血流将放射性同位素带离使得脑部记数率的减量为2N。
(3)根据文献[1]，放射性元素自身有衰减，设其半衰期为，由此引起的记数率下降为3N
又由医学实验和假定有：
)()(21tKNdtdNtkPdtdN，
又
311(t)ln2(t)2tNNt
3ln2(t)dNNdt
所以考虑△: 时刻内头部放射性元素记数率变化, 有
123(t)N(t)(t)(t)NNN
其变化率为：
312(t)(t)dNdNdNdNdtdtdtdt
于是得到
3ln2()()dNkPtKNtNdt (1)
由于在测试时放射性同位素的半衰期一般很大, 这样会给测量和试验带来严重影响,
因此假定: 于是( l) 变为
()()dNkPtKNtdt （2）
分析式(2)，要确定脑血流系数的模型，必须分析)(tP和)(tN的实验数据，观察其变化趋势。首先用Matlab绘出)(tP和)(tN的散点图并观察其变化 5 趋势，)(tP和)(tN的散点图如下：
1234567891005001000150020002500tP(t)
图5-1 呼出气记数率的散点图
1234567891002004006008001000120014001600tN(t)
图5-2 脑部记数率的散点图
由上图可知，t和)(tP有近似于btae的关系，而t和)(tN的关系暂时不能 6 直接观察出，设btaetP)(，用MATLAB的拟合工具箱可得到参数a=1000，b-1.5，而其相关系数为1，说明拟合的非常精确。拟合后的图像如下图所示：
1234567891005001000150020002500
P vs. xfit 1
图5-2 拟合函数)(tP的图像
由基本假设5，即0)0(N和式(1)联立可得一带有初值的微分方程：
0)0()(10005.1NtKNkedtdNt (3)
解此微分方程得到
)(5.11000)(5.1KtteeKktN (4)
5.2 求解脑血流系数的算法模型
以下为求解)(tN的算法模型：
算法模型I：再分析)(tN的试验数据，同样用Matlab的拟合工具箱里现有的指数函数可精确拟合出)(tN的函数，表达式为：)(tN=ttee479.15015.039834032，可近似认为5.1479.1,39834032，则对比系数可得到3977.0,5015.0kK，但这种方法是不精确的，只是用于后面 7 方法得到参数的验证。
对于参数K和k的求解，可以利用Matlab里专门求解拟合后曲线参数的函数lsqcurvefit求得，0.4001k0.5000, K。
最后对得到的数据进行误差分析，把求得的K和k的值代入（4）式即可得到)(tN的函数，求出预测值，再和实际值进行作差，用Matlab绘出作差后的值和时间t的关系图，如下所示：
012345678910-20-15-10-505101520
图5-3 误差分析比较图
由图可知实际值与预测值的差值很小，在0上下波动，比较符合题意。
六 模型的评价与推广
1在建模时忽略了动脉血从肺部到脑部所需要的时间，如在模型考虑这些因素后，只须在测试中测得这些因素的数值，用上述方法仍是容易实现的。
2本题模型比较简单，利用解微分方程和最小二乘法求得的脑部血流系 8 数比较精确，拟合出的)(tN和)(tP的函数精度很高，第I种算法模型比较稳定。
七 参考文献
[1] 佚名/HotWord/1269774.htm 2013.08.22
[2]王家文，王皓，刘海.MATLAB7.0编程基础[M],北京：机械工业出版社，2005.7
[3]刘志平，石林英.最小二乘法原理及其MATLAB实现[J],中国科技西部，2008，17(7):33-34
9
附录
附录1
附表1 某受试者的测试数据
时间 头部记数率 呼出气记数率 时间 头部记数率 呼出气记数率
1.00 1534 2231 5.75 225 2
1.25 1528 1534 6.00 199 1
1.50 1468 1054 6.25 175 1
1.75 1378 724 6.50 155 1
2.00 1272 498 6.75 137 1
2.25 1162 342 7.00 121 0
2.50 1052 235 7.25 107 0
2.75 947 162 7.50 94 0
3.00 848 111 7.75 83 0
3.25 757 76 8.00 73 0
3.50 674 52 8.25 65 0
3.75 599 36 8.50 57 0
4.00 531 25 8.75 50 0
4.25 471 17 9.00 44 0
4.50 417 12 9.25 39 0
4.75 369 8 9.50 35 0
5.00 326 6 9.75 31 0
5.25 288 4 10.00 27 0
5.50 255 3
附录2
x=[1.00:0.25:10.0];
y=[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];
a0=[0 0];
a=lsqcurvefit(@(a,x)(10000*a(1)*(exp(-1.5*x)-exp(-a(2)*x)))/(a(2)-1.5),a0,x,y)
附录3
x=[1.00:0.25:10.0];
for i=1:37
y=10000*0.4001*(exp(-1.5*x(i))-exp(-a(2)*x(i)))/(0.5000-1.5)
end
附录4