1、matlab与数值分析糖尿病模型的诊断一、背景知识二、模型的建立三、模型的分析四、总结参考文献附录 糖尿病诊断模型一、背景知识糖尿病 (DiabetesMellitus, DM)是目前病因和发病机理尚未完全认识的常见的内分泌代谢疾病，它是由于胰岛素分泌不足或作用缺陷所引起的以慢性高血糖为主 ,合并脂肪和蛋白质代谢紊乱为特征的综合征。其早期无症状；晚期多尿、多食、多饮，消瘦、无力等临床表现，化验结果血糖、尿糖指标增高。其危害为：易并发感染心脏病、肾病、视网膜病变，引起心肌梗死、中风、肾功能衰竭、白内障失明等后果。糖尿病的产生一般与人体内血糖的调节系统失衡有关。血糖（Blood Glucose）指
2、的是血液中的葡萄糖，它能为体内所有组织和器官提供能源。当人处于安静、空腹状态时，体内正常的血糖含量为 70 100 mg / 100ml，当体内血糖持续过高或过低时称之为为病态。血糖在其自我调节的过程中，受到各种激素及其代谢物质的影响和控制，其中参与血糖调节的激素主要有胰岛素、胰高血糖素、肾上腺素以及其他相关激素等。胰岛素是由胰岛B细胞分泌的一种蛋白质激素，其作用是促进糖、脂肪、蛋白质三大营养物质的合成代谢，它的最主要功能是调节糖代谢，促进全身组织对糖的摄取、储存和利用，从而使血糖浓度降低1。胰岛素是体内唯一能降低血糖浓度的一类激素，它不能直接发挥作用，必须和所要结合的细胞膜上的胰岛素受体紧密
3、结合后，才能产生生理效应。胰岛素受体是一种特殊的蛋白，主要分布在肝脏、肌肉、脂肪等组织的细胞上，它对胰岛素特别敏感，而且识别性极强。如果把胰岛素受体比作是一把锁，那胰岛素就是一把钥匙。胰岛素发挥降血糖的过程就好象是用钥匙打开锁，使细胞的大门打开，血液中的葡萄糖迅速进入细胞内并被利用，从而使血液中的血糖量降低。胰高血糖素是由胰岛A细胞分泌的一种激素，它的作用主要是抑制机体对血糖的吸收，使体内血糖含量升高。肾上腺素对血糖的调节作用和胰高血糖素对机体的调节作用是相同的。体内血糖调节的过程比较复杂，为了便于分析，对其模型进行简化得到人体内血糖调节的简单模型如下图图1.1所示，当体内血糖浓度高于正常水平
4、时，血糖含量的变化会刺激胰岛B细胞分泌胰岛素，胰岛素的增加会促使体内血糖一部分被肝脏吸收转化为肝糖原储存在肝脏内，另一部分血糖则被组织所吸收，从而降低体内血糖的含量；当体内血糖浓度低于正常水平时，血糖含量的变化会刺激体内激素（胰高血糖素和肾上腺素含量的升高），激素含量的增加会促使肝脏中的肝糖原转化为葡萄糖同时组织将不再吸收血糖，此时体内的血糖含量升高2。组织吸收肝脏释放内分泌组织吸收组织拒绝图1.1 人体内血糖调节机制的简单模型二、模型的建立 糖尿病的诊断通常是通过葡萄糖耐量测试（GTT）来检查的，较严重的糖尿病患者很容易被发现，但是对于轻微的糖尿病患者不太容易被诊断出来。二十世纪60年代中期
5、，北爱尔兰马由医院的医生Rosevear和Molnar以及美国明尼苏达大学的Ackeman和Gatewood博士研究了血糖循环系统，建立了一个简单的数学模型，该数学模型是基于图一血糖调节的简单模型而建立的，为轻微糖尿病的诊断提供了较为可靠的依据。葡萄糖耐量测试（GTT）可以作为轻微糖尿病患者诊断的依据，其诊断的主要过程为： 患者前三天正常进食； 前一天晚餐后禁食； 测试当天早上先抽血化验血糖量，再注射75g葡萄糖，在0.5h3h内间隔一定时间对受试者做葡萄糖浓度检测。那么如何对GTT结果进行正确的分析，确定被检测者是否患有糖尿病？为了便于分析和计算，提出了两种假设：H1：假设对血糖浓度有调节作
6、用的各种激素可以用一个综合参数H描述，H升高导致血糖浓度G下降；H下降导致血糖浓度G升高，参数正常值为H0；H2：G和H相互制约相互影响，自动向正常水平G0和H0调节；H3：g=G-G0,h=H-H0相对较小，h代表激素偏离正常水平的变化量，g代表人体内血糖浓度偏离正常水平的变化量。那么，根据以上假设可以建立数学模型，具体如公式1-1所示。式1-1 其中，上式中和为血糖浓度和激素之间的关系函数，为t时葡萄糖的外界摄入量，如果t时刻从外界摄入体内的葡萄糖浓度为零，并且当血糖浓度和激素水平都为正常值时，则有式1-2 其中，G0，H0为正常水平值令g=G-G0,h=H-H0,代入式1-1，并在（G0
7、，H0）处展开成泰勒级数的形式如公式1-3所示：式1-3其中、是g 和h 的高阶无穷小量，并且当、很小时（即检测者至多为轻微病人时），为求解方便，该方程组可近似为：式1-4下面估计式1-4方程组右端各项的符号，（不妨设J(t)=0）：若g 0 且 h = 0，则此人血糖浓度高于正常值，内分泌激素将促使组织吸收葡萄糖，并将其存储进肝脏，此时有,从而应有其激素浓度将增加将抑制血糖浓度的增高，因而又有 若g=0且h0，则说明此人激素浓度高于正常值，血糖浓度及激素浓度均将减少，从而必有:式1-5 则式1-4可以改写为：式1-6 其中，m1、m2、m3、m4均为正数式1-6是关于 g、h的一阶常系数微分
8、方程组，因激素浓度不易测得，对该式再次求导可以化为式1-7：式1-7 式1-8为了简化计算，可令：式1-9 式1-10 式1-11 于是，式1-8可以简写为：式1-12 为了对式1-12方便求解，设在t = 0 时患者开始被测试，被试者需在很短时间内喝下一定数量的外加葡萄糖水，如果忽略这一小段时间，此后方程可写成：式1-13 求得特征方程：式1-14 特征方程的解为：式1-15 则式1-13的解可以写成：式1-16 则血糖浓度最终可表示为：式1-17 对该模型进行求解的方法一般有两种：1、解非线性方程组：线性方程的根可以用解析法求得，但大多数非线性方程不能用解析法求解，而需要用数值方法。非线性
9、方程定义为如果方程中含有变量的非线性函数、或者变量的高次项，则方程为非线性方程3。多项式方程也是非线性方程，表示如下：多项式方程的根就是可以使f(x)等于0的变量x的值。如果f(x)是二次方程，则可求得解析解。三次以及次以上的高次多项式方程，有3个或3个以上的根。除了方程根的数目外，还要重视根的性质。非线性方程的根可为：a.不同的实根；b.重复的实根；c.复数根；d.各种根的组合。 求解非线性方程组的方法常见的有一下几种：逐次代换法：逐次迭代法为单点迭代搜索根：先变换方程f(x)=0,使方程的左边只有一个x ，即x=g(x)，则直线y=x与曲线y=g(x)的交点就是方程的根。选取一个初始近似值
10、x1,可以求得x2：x2=g(x1)；再把x2作为新的初始值，进行下一个迭代。以此类推可以得到xn+1=g(xn)。则其收敛条件为：对于搜索区间上所有x的值，必须满足，在计算 逐次代换法的优点是只需要一个出发点，不需要求函数的导数。 试位法(线性插值法)：试位法寻根的方法是选取根两边的点，在函数中这两点间进行线性插值。优点：无需求函数的导数；计算简单，是求解非线性方程最方便的方法之一。缺点：搜索根时，收敛的准确度和速度受到初始值的限制，初始值成为所有后续代换的关键点。 牛顿一拉弗森法 ：第一步：非线性方程表示为：f(x)=0；第二步：利用泰勒级数在任一点x1展开。总结三种方法：（1）N-R最常
11、用，其搜索过程对于根的初始估计值的选定敏感；（2）逐次代换法确定根的迭代部数最少，计算简单，不需要计算函数的导数（3）线性插值法需要给定两个初始估计值，用于设定搜索范围的起点和终点，所需的迭代次数最多。牛顿法求解非线性方程组：如果模型中包含两个或更多个联立非线性方程，可扩展N-R法同时求解这些方程。该模型可以用非线性方程组的解法进行求解，取t=ti(i=1,2,3,4)各测一次，将测得的值代入式1-14，得到一个方程组，由此可解得相应的参数值。2、最优化方法-最小二乘法： 最小二乘法又称最小平方法，是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得
12、未知的数据，同时满足这些求得的数据与实际数据之间误差的平方和为最小。 为了量化最小误差平方和，借用统计中残差的概念来定义最小二乘法使用过程中的优化判据。残差是指估计或预测数据与实测值之间的差值，它反映的是估计值偏离实测值的程度，是估计误差的一种衡量。最小二乘法的优化判据即是使所有样本的残差平方和最小，残差平方和的定义如公式1-15所示：式1-18 求解该模型也可以借用式1-15最小二乘法的方法进行求解。三、模型的分析以上内容介绍了模型的建立和求解，并利用非线性方程组或者最小二乘法对该模型进行求解，最终可以求得模型的参数。我们知道在进行糖尿病诊断时，通常用葡萄糖耐量测试（GTT），那么就可以利用
13、这些模型参数中的其中一些参数对GTT进行估计。首先，介绍一下参数灵敏度的有关知识。假设模型y=f(x;a)，其中a是参数；如果存在M10，使得|S(y,a)|0，使得|S(y,a)|M2，则称模型y在水平M2下关于参数a是灵敏的。在实际计算中发现，G的微小误差会引起的很大偏差，故任一包含的诊断标准都将是不可靠的，同时也可发现G对并不十分敏感（计算结果与实际值相差较小），故可用的测试结果作为GTT检测值来判断此人是否真的患有轻微的糖尿病。为了判断上的方便，故一般利用自然周期 作为判别标准。根据人们的生活习惯，两餐之间的间隔时间大体为4小时。临床应用显示， 在T 4h? 图1.2 GTT诊断流程图
14、 根据流程图利用matlab进行编程仿真，并根据临床患者血糖的观测数据进行仿真，得到了不同测试者体内血糖的浓度随时间的变化图。下图1.3所示是正常G（t）模型图。从正常的G（t）模型图中，可以看出人体内血糖的浓度随着时间的变化在正常值上下波动，最终基本处于正常水平。这是由于体内激素对血糖的调节作用，当血糖浓度升高时，激素的总体水平使血糖浓度下降；当血糖浓度降低时，激素的总体水平促使血糖升高。总之，由于体内激素对血糖的正常调节作用从而使血糖维持在正常的水平。 图1.3 正常G（t）模型图 在t=0,0.5,1,2,3时刻可观测到患者A体内血糖的值为G0=6.4，G1=8.15，G2=11.43，
15、G3=12.55，G4=9.9，将不同时刻的观测值带入到模型方程中，联立方程组可计算出模型参数的值为： 根据所求的参数值计算出自然周期： 从自然周期T值中可以判断出该病人为轻微糖尿病患者，该患者的G（t）模型图如下图图1.4所示，从该图中可以看出当患者体内血糖含量升高时，激素的总体水平不能使血糖的浓度及时的降到正常水平的范围，由此可导致糖尿病的发生。 图1.4 轻微糖尿病患者的G（t）模型图 在t=0,0.5,1,2,3时刻测到患者B体内血糖的值为G0=6.51，G1=11.41，G2=16.18，G3=12.84，G4=5.33，将不同时刻的观测值带入到模型方程中，联立方程组可计算出模型参数
16、的值为： 根据所求的参数值计算出自然周期： 从自然周期T值中可以判断出该病人不是糖尿病患者，该病人的G（t）模型图如下图图1.5所示，从该图中可以看出病人体内的血糖浓度随着时间在正常值上下波动，最终趋于正常水平。在该病人的体内，激素发挥了正常的作用，当血糖浓度升高时，激素的总体效应使之降低；当血糖浓度降低时，激素的总体效应使之升高，随着时间的变化最终趋于正常水平。 图1.5 非糖尿病病人的G(t)图四、总结在该实验的设计中，首先是对人体内激素对血糖浓度的调节过程做了简化，然后在一些假设前提下，建立激素和血糖之间关系的数学模型，并利用数学中的解方程组的方法对该模型进行求解，得到相应的模型参数，然
17、后对这些模型参数的灵敏度进行研究，最终找到了可以表征轻微糖尿病的参数-自然周期T，当T4h时，可以认为该病人是轻微的糖尿病患者；当T4h时，可以认为该病人是非糖尿病患者，但在该模型中也存在一些缺点。糖尿病是一种新陈代谢疾病，它是由胰岛素缺乏引起的新陈代谢紊乱造成的。糖尿病的诊断是通过葡萄糖容量测试（GTT）来检查的，较严重的糖尿病医生不难发现，较为困难的是轻微糖尿病的诊断。在轻微糖尿病的诊断中，由于内分泌激素浓度不易测量，在上面的建模过程中对各种不同激素未加以一一区别，即对其采用了集中参数法。这样做虽大大简化了模型，但也在一定程度上影响了模型的应用效果。临床应用时发现，在患者饮下葡萄糖水大约1
18、-3小时后，测得的数据有一定的偏差，其原因可能是内分泌激素的作用造成的，因而，要得到更精确的结果，当然要考虑到内分泌激素浓度的变化，建立更精确的模型。罗德岛医院已找到一种测量内分泌浓度的方法，相信在此基础上一定可以设计出诊断轻微糖尿病的更好方法。参考文献1 Pritpal Bansal and Qinghua Wang. Insulin as a physiological modulator of glucagon secretion, Am J Physiol Endocrinol Metab295:E751-E761, 2008; irstpublishedJuly22,2008;doi
19、:10.1152/ajpendo.90295.2008.2 冯海英，丁焕然，史视明. 2型糖尿病患者胰岛素抵抗与胰高血糖素的关系. 中国综合临床2006年8月第22卷第8期:707-708；3 Stanley, M, Dunn, Alkis, Constantinides. 数值方法在生物医学工程中的应用M. 北京:机械工业出版社, 2009. 100-160附录1、参数的输入clear all;clct0=input(shuru t=);%输入参数t1=input(shuru t=);t2=input(shuru t=);t3=input(shuru t=);t4=input(shuru t
20、=);G0=input(shuru G0=);G1=input(shuru G1=);G2=input(shuru G2=);G3=input(shuru G3=);G4=input(shuru G4=);D=t0,t1,t2,t3,t4,G0,G1,G2,G3,G4;save(E:matlab2012bbincanshu.mat,D)2、模型方程组的建立function q=myfun(p) load(canshu.mat);w0=p(1);alpha=p(2);deta=p(3);B=p(4);w2=w02-alpha2;A=0,1; -w02,-2*alpha;q(1)=D(6)+B*e
21、xp(-alpha*D(1)*cos(w2)1/2*D(1)-deta)-D(6);%联立方程组q(2)=D(6)+B*exp(-alpha*D(2)*cos(w2)1/2*D(2)-deta)-D(7);q(3)=D(6)+B*exp(-alpha*D(3)*cos(w2)1/2*D(3)-deta)-D(8);q(4)=D(6)+B*exp(-alpha*D(4)*cos(w2)1/2*D(4)-deta)-D(9);q(5)=D(6)+B*exp(-alpha*D(5)*cos(w2)1/2*D(5)-deta)-D(10);end3、模型的求解clear all;clcload(can
22、shu.mat);x=fsolve(myfun,4,3,2,2,optimset(Display,off)%求解方程组t=0:0.01:3;w2=x(1)2-x(2)2;g=D(6)+x(4)*exp(-x(2)*t).*cos(w2)1/2*t-x(3);plot(t,g)xlabel(时间/h),ylabel(血糖浓度/mmol/L)hold ong0=D(6);g1=find(t=D(2);g2=find(t=D(3);g3=find(t=D(4);g4=find(t=D(5);plot(0,g0,rs)%画出G（t）图text(0,g0,P0( num2str(0) , num2str
23、(g0) )plot(D(2),g(g1),rs)text(D(2),g(g1),P1( num2str(D(2) , num2str(g(g1) )plot(D(3),g(g2),rs)text(D(3),g(g2),P2( num2str(D(3) , num2str(g(g2) )plot(D(4),g(g3),rs)text(D(4),g(g3),P3( num2str(D(4) , num2str(g(g3) )plot(D(5),g(g4),rs)text(D(5),g(g4),P4( num2str(D(5) , num2str(g(g4) )hold onplot(t,g0,r )text(D(2),g(g1),g)T=2*pi/x(1)%求解自然周期Tdisp(end)
