局部脑血流的测定

局部脑血流的测定摘要随着医学的快速发展，脑部的研究越来越得到人们的重视。而对脑血流系数的精确测定可以帮助人们快速得到相关指标，所以对医学领域来说很有意义。问题1，首先我们根据题设中的信息及放射性元素衰减性，列出头部记数率和呼出记数率关于脑部血流系数的微分表达式；再分析试验得到的呼出气记数率的数据，用Matlab的拟合工具箱拟合得到关于呼出气记数率的函数；最后代入上述微分表达式求解微分方程即可得到关于脑部血流系数的表达式，表达式为)(5.11000)(5.1KtteeKktN。问题2，先用Matlab拟合工具箱求出头部记数率的函数，用对比系数法可得到K和k的近似值，即3977.0,5015.0kK，但这种方法是不精确的，只是用于后面方法得到参数的验证；再利用最小二乘法求解拟合后曲线参数的函数以及由问题1得到的关于脑部血流系数的表达式，即可得到脑血流系数：0.5000K，4001.0k，最后对得到的值进行误差分析，可知脑血流的预测值和实际值很吻合，比较符合题意。关键词：对比系数法放射性元素衰减曲线拟合2一问题重述用放射性同位素测定大脑局部脑血流量的方法如下：由受试者吸入含有某种放射性同位素的气体，然后将探测器置于受试者头部某固定处，定时测量该处的放射性记数率(简称记数率)，同时测量他呼出气的记数率。由于动脉血将肺部的放射性同位素传送至大脑，使脑部同位素增加，而脑血流又将同位素带离，使同位素减少。实验证明由脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比，其比例系数反映了该处的脑血流量，被称为脑血流量系数，只要确定该系数即可推算出脑血流量。动脉血从肺输送同位素至大脑引起脑部记数率上升的速度与当时呼出气的记数率成正比。若某受试者的测试数据如附表1所示：根据以上题目所给的条件及数据，回答以下问题：1、建立确定脑部血流系数的数学模型；2、计算上述受试者的脑血流系数。二模型假设1.脑部记数率的上升只与从肺部输送的放射性同位素有关；2.脑部记数率的下降只与当时该处的脑血流量有关；3.脑血流量在测定期间恒定，心脏博动、被测试者大脑活动、情感波动等带来的变化忽略不予考虑；4.每次仪器测量为相互独立事件，各测量值无记忆相关；5.在吸入气体瞬时，脑中放射物记数率为零；6.脑血流量与脑血流量系数成单值函数关系，求得后者即可确定前者。三符号说明符号表示意义3t时间k脑部记数率上升的速率与呼出气记数率的比例系数K脑部记数率下降的速率与当时该处脑部记数率的比例系数，即脑血流量系数)(tN脑部记数率)(tP呼出气记数率1N脑部记数率的增量2N脑部记数率的减量3N放射性衰减引起的脑部记数率的减量四问题分析问题1，首先根据题设可知：一方面，由脑血流引起局部地区记数率下降的速率与当时该处的记数率成正比,并且其比例系数反映了该处的脑血流量；另一方面，动脉血从肺输送同位素至大脑引起脑部记数率上升的速率与当时呼出气的记数率成正比，由上述两方面可得到头部记数率关于脑部血流系数和呼出气记数率的表达式；再分析试验得到的呼出气记数率的数据，用MATLAB的拟合工具箱拟合得到关于呼出气记数率的函数；最后代入上述表达式求解微分方程即可得到关于脑部血流系数的表达式。问题2，利用Matlab里专门求解拟合后曲线参数的函数，再由问题1得到的关于脑部血流系数的表达式即可得到脑血流系数。最后再对得到的值进行误差分析。五模型的建立与求解5.1求解脑血流系数的数学模型设某时刻0t时，脑部记数率为)(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，说明拟合的非常精确。拟合后的图像如下图所示：1234567891005001000150020002500Pvs.xfit1图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.0kK，但这种方法是不精确的，只是用于后面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]佚名[2]王家文，王皓，刘海.MATLAB7.0编程基础[M],北京：机械工业出版社，2005.7[3]刘志平，石林英.最小二乘法原理及其MATLAB实现[J],中国科技西部，2008，17(7):33-349附录附录1附表1某受试者的测试数据时间头部记数率呼出气记数率时间头部记数率呼出气记数率1.00153422315.7522521.25152815346.0019911.50146810546.2517511.7513787246.5015512.0012724986.7513712.2511623427.0012102.5010522357.2510702.759471627.509403.008481117.758303.25757768.007303.50674528.256503.75599368.505704.00531258.755004.25471179.004404.50417129.253904.7536989.503505.0032669.753105.25288410.002705.502553附录2x=[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=[00];a=lsqcurvefit(@(a,x)(10000*a(1)*(exp(-1.5*x)-exp(-a(2)*x)))/(a(2)-1.5),a0,x,y)附录3x=[1.00:0.25:10.0];fori=1:37y=10000*0.4001*(exp(-1.5*x(i))-exp(-a(2)*x(i)))/(0.5000-1.5)end附录410clccleara=3996;b=-0.4998;c=-3998;d=-1.502;N=[15341528146813781272116210529478487576745995314714173693262882552251991751551371211079483736557504439353127];P=[2231153410547244983422351621117652362517128643211110000000000000];forx=1:0.25:10;i=4*x-3;Nff(i)=(a*exp(b*x)+c*exp(d*x));plot(x,Nff(i)-N(i),'*');holdon;axis([010-2020])endNff-N