北京师范大学 数学模型与数学建模1 (10).pdf
第六章 医学与遗传的数学模型 6.1 糖尿病的诊断 一.背景和问题 1.糖尿病（Diabetes）：胰岛素相对或绝对不足引起的以糖代谢紊乱，血糖增高为主的慢性疾病。10.症状：早期无症状；晚期多尿、多食、多饮，消瘦、无力等临床表现，化验结果血糖、尿糖指标增高。20.危害：易并发感染心脏病、肾病、视网膜病变，引起心肌梗死、中风、肾功能衰竭、白内障失明等后果。30.血糖调节系统 血糖（Blood Glucose）：血液中的葡萄糖，为体内所有组织和器官的能源。安静、空腹时正常含量 70 100 mg/100ml。持续血糖过高或过低为病态。血糖在其自我调节的过程中，受到各种激素（Hormones）及其代谢物质的影响和控制。其中有 胰岛素（Insulin）：由胰脏的 细胞分泌的一种蛋白质激素。血糖（BG）增高时刺激分泌胰岛素促进组织对葡萄糖的吸收，并将多余的糖变成糖原(Glycogen)存储于肝脏中。高血糖素（Glucagon）：由胰脏的细胞分泌的一种蛋白质激素。血糖降低时刺激其分泌，抑制肌体对糖的吸收，并增加糖原转化为 G 的速度。肾上腺素（Epinephrin）肾上腺髓质分泌的激素。属于应急机制的一部分。当 BG 极低时，加快糖原转化成葡萄糖的速度，抑制分泌胰岛素和细胞吸收葡萄糖，以便很快增加血液中葡萄糖的浓度。其它：糖皮质素、甲状腺素、生长激素等。40.致病机理：人体自身不能提供足够的胰岛素以消耗体内多余的糖分，致使血糖增高而无法控制。50.诊断（1960 年前制定的诊断）方法：葡萄糖耐量检查GTT（Glucose Tolerance Test）：空腹测血糖G0；服用G=75 克葡萄糖，并于 0.5、1、2、3 小时测血糖，描出血糖曲线。比较血糖曲线与标准的偏差 正常人的血糖曲线 糖尿病人的血糖曲线 比较 从复杂的检查结果难以区别正常人与糖尿病患者 2.问题：如何通过血糖曲线的分析给出一个可靠的解释 GTT 检验结果的标准。1965 年左右，美国明尼苏达大学的 Ackerman 和 Gatewood 博士 建立了如下的模型，用以描述 GTT 过程中的血糖调节系统，再从这个模型中选择一个参数来区别正常人与糖尿病患者。二.建模 1.假设：10.只考虑血糖浓度G和以胰岛素为代表的激素H的相互作用。20.没有外界摄入时，人体内有一个平衡状态G0、H0。30.只考虑血糖和激素浓度在平衡状态G0、H0附近的小偏差。2.简化的血糖调节系统 3.简化调节机制 内分泌 血糖降低激素增加 血糖增高 肝脏贮存 组织吸收 内分泌 血糖升高激素减少 血糖降低 肝脏释放 组织拒绝 4.模型：设G(t),H(t)分别表示t时刻血糖与激素的浓度。则有dG/dt=F1(G,H)dH/dt=F2(G,H)其中，F1(G,H)，F2(G,H)分别表示血糖浓度的变化和激素浓度的变化率对G、H的依赖关系。存在G0,H0，使得 F1(G0,H0)=0 和 F2(G0,H0)=0。称 G0,H0 为平衡态，5.线性近似 令g(t)=G(t)-G0,h(t)=H(t)-H0.则 dg/dt=dG/dt=F1(G0+g,H0+h)dh/dt=dH/dt=F2(G0+g,H0+h)当 g(t),h(t)充分小时,将 F1,F2 在(G0,H0)附近展开并略去高阶无穷小,得到在平衡态(G0,H0)附近的局部线性化方程组 dg/dt=-m1g-m2h dh/dt=m3 g-m4 h 其中 m1=-F1/G|(G0,H0)0,m2=-F1/H|(G0,H0)0，m3=F2/G|(G0,H0)0,m4=-F2/H|(G0,H0)0 设在 t=0 时刻，人体内的G，H达到平衡值G0,H0，经外部注射，输入大剂量J的葡萄糖,则 g(0)=J,h(0)=0.消去方程组中的变量 h(t),可以得到 d2g/dt2+2adg/dt+02=0,g(0)=J 其中 2a=m1+m4 0,02=m1 m4+m2 m3 0 三.求近似解.由于血糖曲线应该成阻尼振荡形式,因此设 a2 -02 4)的BG观测值时可以使用非线性最小二乘法的出这四个参数的估计值。eqns:=8.06=4.5+A*exp(-a*0.5)*cos(w*0.5-b),6.7=4.5+A*exp(-a)*cos(w-b),3.91=4.5+A*exp(-a*2)*cos(w*2-b),4.47=4.5+A*exp(-a*3)*cos(w*3-b);eqns:=8.06=4.5+A exp(-.5 a)cos(-.5 w+b),6.7=4.5+A exp(-a)cos(-w+b),3.91=4.5+A exp(-2 a)cos(-2 w+b),4.47=4.5+A exp(-3 a)cos(-3 w+b)vars:=A,a,w,b;vars:=a,b,A,w solve(eqns,vars);w=-2.811151237,a=.6293589754-6.283185307 I,b=-.6078588509,A=-6.983046177,线性逼近模型的拟和曲线 2.参数的灵敏度 模型 G=G(t;a),其中 a=(,)是参数.若存在 M10,使得|DG/Da|0,使得|Da/DG|M2,则称参数 a 在水平M2下关于模型 y 是不灵敏的。模型对参数不灵敏，则参数估计中的误差将不会对模型值产生过大的影响。参数对模型不灵敏，则模型值的测量误差对参数估计值将不会对产生过大的影响。以上模型中主要参数为：BG 向平衡值衰减的速率；BG 震荡的自然频率.不难计算G/G/.、/G 和/G.。由此不难看出 1.在 t=0 附近参数 ，对模型值都很敏感，它表明 BG 的微小观测误差可能导致参数估计值较大的偏差。因此 BG 的观测不宜在 t=0 附近进行。2.无论参数 对模型还是模型对参数 都有反应十分敏感的时间区间。因此包含有参数 的有关诊断指标是不可靠的。3.模型对系统的自然频率 的反应并不十分敏感。因此它可以作为 GTT 的基本指标来进行诊断。3.诊断指标 取系统的自然周期 t=2/作为 GTT 法诊断糖尿病的指标。对于正常人来说将有 T 4 时，如果检验结果的 T 明显地大于 4 时，将诊断为轻微的糖尿病。五.讨论 1.此模型只适用于轻微的糖尿病。通常饭后 2 小时，静脉血浆糖 140mg/dl 诊断为糖尿病,血糖正常浓度应为 70100mg/100ml。2.在摄入的葡萄糖被吸收 35 时，模型可能与数据拟合不好。它表明激素（肾上腺素）应与胰岛素分离开来单独分析讨论，因为肾上腺素的急剧上升会使 BG 迅速下降。这一现象在模型中是没有考虑到的。3.在模型中反应肾上腺素的作用有赖于对人体内肾上腺素浓度的准确的测量。在多年研究的基础上，糖尿病确诊标准得到不断改进。1985 年 WHO 糖尿病诊断暂行标准 糖尿病 静脉血浆糖(mmol/L)空腹血 7.8 GTT 2h 11.1（200mg/dl）G耐量减低 空腹血 k.，当t 充分大时，e-ka t0，先设 C(应t大的测量值(t,C)用最小二乘法确定k和A。然后对=A e-ka t 取对数，则 Ln(Ae-kt-C)=-ka t+ln(A)，取对应t小的测量值(t,C)用最小二乘法确定ka和A。Matlab 程序 T=0.5 1 C=3.71 4.93 5.5 5.7 5.6 5.33 4.8 4.1 3.1 2.2 1.8 1.4;b=polyfit(T(5:end),log(C(5:end),1)b=-0.1362 2.0804 a=polyfit(T(1:4),log(exp(b(2a=-1.4876 2.1485 T1)-exp(e(3).*T1),e,T1)fun=inline(e(1).*(exp(e(2).*ee=nlinfit(T,C,fun,exp(2.1),-0.13,-1.49)ee=.2863 8 -0.1401 -1.3515 C2=ee(1).*(e plot(T,C,*,T1,C1,r,T1,C2)拟合曲线 吸收半衰期T1/2ka=ln(2)/ka和清除半衰期T1/2k=ln(2)/k,时间Tmax 表观分布容积V，清除率 持剂量，给药时间间隔 持剂量，给药时间间隔 最大治疗浓度 Cmax 为稳态血药浓度 cs(t)。n 注 要求 smin 即 D*/Ve-k=D/V(ek-1)/(1-e-k)=VCsmax=VCmax)还可以确定血药浓度的最大峰值Cmax 和达到最大峰值的cl=kV，等参数。(四）给药方案 S=D*,D,负荷剂量，维S=D*,D,负荷剂量，维1 要求首次给药剂量 D*使得血药达到2 当给药次数 n 称血药浓度的极限值要求 cs(t)的极小值 Csmin=最小治疗浓度 Cmics(t)的极大值 Csmin=最大治疗浓度 Cmin 静 knkktnjjtkneeeCCte=11)(010)(0 C用(n-1)+t 取代 t,再令 n得 0(aij 0),则称 M 是正则的,M所,n)T,使得 M=.参见 矩阵论Perro-Frobenious 定理）则链具有性质：正则链的状态向量一定趋向于一个稳定的分布,并且这个分布与状态向量的体模型的转移概率矩阵 Mr Mr0,有唯一的不动点F*=(p2,2pq,q2)T,即有 Mr F*=F*而 中 Ir 为 r 阶单位阵,O为零阵,R 的行和均为正.则称由 M 给出的马氏链S(n+1)=MS(n)状态 S1 称为吸收态,其余的状态 S2 为非吸收态，这里 S=(S1,S2).0 0)吸收前转移到非吸收状态 sj 的平均转移次数.则 F-1-1,就会被吸收状态S1,S2 吸收.RMrRI=2/14/14/1SM2/1,4/14/1,1001=QRI=qqpqppMr2/02/102/=TnnM1),(lim?描述的链为正则链.M具有如下性质。M 有唯一的不动点向量 =(1,?=)0(1)0(lim)(mFFMnFTnnn?li（正初始取值无关.例 随机交配群2 且当 F(n+1)=MrF(n)时,必有limnF(n)=F*20.吸收链.若状态转移矩阵有如下形式，其为吸收链。矩阵 Ir所对应的 吸收链具有性质：T(1)lim S(n)S(n)T=(S S1 1T,(2)称 F=(IQ)-1 为吸收链的基本矩阵.F=(fij),其中fij 为从非吸收状态 si 在被(3)令 B=F R=(bij),则给出了从非吸收状态 fi出发,被吸收状态 fj 吸收的概率.例 自交群体遗传模型的转移矩阵 Ms =(I-Q)=(1-)=2 从非吸收状态S3出发,平均转移 2 次=QOI0lim,=nnnnrnQQOM0010014/1102/1004/10114/1002/1004/11B=FRT=2(1/4，1/4)=(1/2,1/2)从非吸收状态S3出发,被吸收状态S1,S2吸收的概率分别为 1/2.问题:1.乒乓球比赛，每胜一球得一分，先得三分者胜一局。建模描述一局乒乓球的比赛过程。2.在上面的乒乓球比赛中，如果双方打成 2:2 平局，则在后面的比赛中需要一方连续胜两球才能赢得此局。试给出模型描述一局乒乓球比赛的过程。
