《基因调控网络模型毕业论文设计.doc》由会员分享，可在线阅读，更多相关《基因调控网络模型毕业论文设计.doc（41页珍藏版）》请在三一办公上搜索。
1、绪 论11课题背景随着基因组学的发展，在短时间内可获得生物体基因表达的海量数据，这为研究和揭示基因及其产物之间的相互关系，特别是基因表达的时空调控机制奠定了基础。基因表达的调控不是单一的，孤立的，而是彼此联系，相互制约的，构成了复杂的基因表达调控网络。几乎所有的细胞活动和功能都受基因网络调控。孤立的研究单个基因及其表达几乎完全不能确切地反映生命现象本身和内在规律。因此，科学家必须从系统的观点研究多基因的调节网络，才能阐明生命的本质和疾病发生的机理。因此，基因调控网络是后基因组时代研究的重要课题。基因网络研究的目的是通过建立基因转录调控网络模型对某一个物种或组织中的全部基因的表达关系进行整体的模
2、拟分析和研究，在系统的框架下认识生命现象，特别是信息流动的规律。调控可在分子水平上分为三个层次：DNA水平、RNA水平和蛋白质水平。DNA水平主要是研究基因在空间上的关系影响基因的表达；RNA水平上，也就是转录水平上的调控，主要研究代谢或者是信号传导过程决定转录因子浓度的调控过程；蛋白质水平主要研究蛋白质翻译后修饰加工，从而影响基因表达的活性和功能。12 基因调控网络模型研究现状基因调控网络研究，离不开生物信息学和系统生物学。运用生物信息学的方法和技术通过数据采集、分析、建模、模拟和推断等手段研究复杂的网络关系，揭示有关的作用机理，是当前生命科学的热点之一。从一个初始提供的基因表达数据和调控机
3、制先验知识建立起的模型开始，系统行为可以在各种各样的实验条件下进行模拟。将模拟表达数据和观察得到的真实表达数据做以比较，就可以给出模型适应性评价。如果模拟数据和真实数据不匹配，并且真实数据可靠，则有必要继续对模型进行修正。不断地重复重构和修正网络模型，模拟生物系统行为，测试的预测结果的工作，直到获得一个理想的 模型，其可以很好的与真实基因表达数据相吻合。一般来说基因调控网络有如下特点：(1) 复杂的高维结构 网络中节点和边的数目庞大。染色质构象，结构基因上游各种各样的调节序列，反式作用因子，RNA 聚合酶活性等因素都决定了基因调控网络是一个非常高维的复杂结构。(2) 网络结构动态性 两个基因间
4、的那条边是否存在、作用的方向在不同时期是可能不一样的。 (3) 网络中节点间作用类型复杂多变 目前的研究表明，基因间的相互作用可能是一种非线形的作用关系。在多因子调控模式中还要考虑不同的调控因子对同一个目标调控基因产生作用时的某种逻辑关系，这种逻辑关系是由调控模式中各调控因子的相互关系决定。 (4) 节点类型多种多样 基因调控网络中节点的元素是多种多样。可以是DNA、mRNA、蛋白、分子、大分子、外界环境等等。 (5) 各节点状态是不断变化的 当细胞受到不同的外界刺激或处于不同的发育阶段时，参与表达的基因是不同的，基因的表达量的变化会影响到相互作用的变化，会引起网络结构的变化。(6) 一般存在
5、反馈循环回路 在生物体中各种生理上的周期现象，容易理解生物体中的相互作用存在周期性，至少在网络的局部上是循环的，存在反馈循环回路。要利用数学模型建立基因调控网络时，我们必须考虑以上网络的特性。但同时也需对这些特性进行抽象化、简单化，如对随机性的描述有时只能用确定性模型。在很大程度上，基因的表达信息反映了生物体的阶段功能和状态，这就驱使人们对大规模的基因表达数据进行研究，建立基因调控和反馈网络。目前，许多实验室和研究者应用芯片技术探索基因调控网络，并取得了一些有效的结果。目前除了本文正文使用的微分方程模型以外，常用的基因调控网络模型有布尔网络模型、线性组合模型、加权矩阵模型、贝叶斯网络模型、微分
6、方程模型等等。121 布尔网络模型布尔网络模型最早由Kauffman于1969年引入的，奠定了使用布尔网络研究基因调控网路的基础。在布尔网络中,每个基因所处的状态或者是“开”，或者是“关”。状态“开”表示一个基因转录表达,形成基因产物，而状态“关”则代表一个基因未转录。基因之间的相互作用关系由布尔表达式来表示，例如: A and not BC 表示“如果A基因表达，且B基因不表达，则C基因表达”。以有向图G= (V , H)表示布尔网络，其中 V是图的节点集合(图中的A、B、C)，每个节点代表1条基因，或者代表1个环境刺激，H表示转录表达路径。在应用方面，Yuh 等人综合以往的研究结果，详细分
7、析了海胆 Stronglocentrotus Purpuratus基因Endol16，研究了如何对这一基因转录水平的基因调控网络进行了精确的逻辑描述，他们基于布尔原理描述了基因的顺式调控系统，并使其能够模拟Endol16给定的转录条件下的表达情况。Arnone和Davidson回顾了一些基因调控系统在结构和功能上的研究结果，他们认为对网络各水平下顺式调控系统的直接分析是非常重要的。122 线性组合模型线性组合模型是一种连续网络模型，在这种模型中，一个基因的表达值是若干个其它基因表达值的加权和。基本表示形式为： 其中，是基因在时刻的表达水平, 是基因在时刻的表达水平,而代表基因的表达水平对基因的
8、影响。在这种基因相互关系表示形式中，还可以增加其它数据项，以逼近基因调控的实际情况。例如，可以增加一个常数项，反映一个基因在没有其它调控输入下的活化水平：。 这样，在给定一系列基因表达水平的实验数据之后，即给定每个基因的时间序列，就可以利用最小二乘法或者多重分析法求解整个系统的差分方程组，从而确定方程中的所有参数，即确定。最终，利用差分方程分析各个基因的表达行为。实验结果表明，该模型能够较好地拟合基因表达实验数据。Dhaeseleer 等人用这种方法分析了大鼠脊髓和海马回的基因表达数据，建立了一个包含有65个基因的调控网络模型。并通过在不同时点上对微分方程的迭代运算，精确地复制出网络调控轨迹，
9、包括脊髓发育、海马回发育到海马回损伤。123 加权矩阵模型加权矩阵模型与线性组合模型相似，在该模型中，一个基因的表达值是其他基因表达值的函数。含有n条基因的基因表达状态用n维空间中的向量表示，的每一个元素代表1条基因在时刻的表达水平。以加权矩阵表示基因之间的相互调控作用，的每一行代表1条基因的所有调控输入，代表基因的表达水平对基因的影响。 在时刻，基因对基因的净调控输入为的表达水平即乘以对的调控影响程度。基因的总调控输入为：。这一形式与线性组合模型相似，若为正值，则基因激发基因的表达,负值表示基因抑制基因的表达，0表示基因对基因没有作用。与线性组合模型不同的是，基因最终表达响应还需要经过一次非
10、线性映射： 该函数是神经网络中常用的Sigmoid函数，其中和是2个常数，规定非线性映射函数曲线的位置和曲度。通过上式，计算出时刻基因的表达水平。在最初阶段，加权矩阵的值是未知的。但是可以利用机器学习方法，根据基因表达数据估计加权矩阵中各个元素的值。对于这样的模型，可以利用线性代数方法和神经网络方法进行分析。实验 表明,该模型具有稳定基因表达水平,与实际生物 系统相一致。在这种模型中还可以加入新的变量, 模拟环境条件变化对基因表达水平的影响。Reinitz 和Sharp利用加权矩阵模型构造了果蝇基因调控网络,以此用来描述果蝇基因在果蝇条纹形成过程中的机制,并找到了在果蝇分节中发挥 重要作用的基
11、因。这个基因网络是一种权重矩阵系统:通过模拟退火优化算法得到相互作用连接参数, 同时空间参数也可由细胞间影响因素确定。124 贝叶斯网络模型贝叶斯网络模型本质上是一种概率图模型。Friedman等人于2000年提出了用贝叶斯网络模型分析基因表达数据的方法。它的基本思想是使用简单的局部概率乘积来近似复杂的高维概率分布。贝叶斯网络引入有向无圈图模型和隐马尔可夫链来描述变量间的联系与相互作用，构建调控网络模型，通常贝叶斯网络可用数对表示。其中，为一有向无圈图，图中结点对应随机变量，在微阵列数据中表示基因的表达向量。中另一部分，表示一组条件概率分布。根据马尔可夫假设：每个变量在给定中的父结点前提下各变
12、量之间相互独立，于是得到随机变量的联合概率分布：其中表示的父结点集合。为了确定的联合概率分布，需要确定上式中出现的各个条件概率，所有这些条件概率构成了。贝叶斯网络的核心就是通过将这种条件独立关系解释为因果关系，并用来表示基因间的因果调控关系。 给定一组基因表达谱数据，一般是通过一个打分函数，利用打分函数在条件独立性的条件下可分解性，采用局部搜索的方法寻找使得得分增加的路径，最后得到得分最大的网络结构的基因调控网络结构和参数：Smith 等人还提出了动态贝叶斯网络模型(DBNs)，这种模型和普通贝叶斯网络模型不同之处在于它将两个节点来表示同一基因前后时间点的表达向量，这种模型的优势在于可以将调控
13、的负反馈和延时因素考虑进去，克服了普通贝叶斯网络是一个无环图带来的不足。Husmeier和Ong等人将这种模型用于微阵列数据，进行基因调控网络的构建工作。125 微分方程模型在24中，假设单基因自调控基因网络为 模型等效图：其中，和分别表示在mRNA和蛋白质在时间时的浓度；和分别是mRNA和蛋白质降解速率；是翻译速率；函数代表转录时蛋白质的反馈调节。考虑到转录的时间延迟，Monk16提出了以下模型： 模型等效图：并显示所观察到三种蛋白质的振荡表达和活性很有可能受到转录的延迟驱动。一个基因调控网络由多个基因的相互作用和调节其他基因的表达蛋白（基因衍生物）组成。刺激和抑制蛋白质在转录，翻译，翻译后
14、的过程控制中基因表达的变化。在本文中，我们考虑以下差分方程5,19所描述的GRNS： 其中分别是第个节点的mRNA和蛋白质浓度。参数和分别是mRNA和蛋白质的衰变率；是翻译速率，函数代表转录蛋白质的反馈调节，通常是一个非线性函数，但每个变量都具有单调性。人们可以发现，网络中的任何单个节点或基因中，有到其他节点输出或基因和来自其他节点的多输入。这些基因网络的结构和调节机理可以作为参考。作为一个单调递增或递减的调控函数，是通常的Michaelis-Menten或Hill形式。在本文中，我们所采取的就是这种所谓的SUM逻辑。也就是说，每个转录因子是基因的一个激励，于是如果转录因子是基因的一种阻遏，于
15、是其中是Hill系数，是一个正的常数，是有界常数(无量纲的转录因子到的转录速率)。因此，我们可以改写系统如下： 其中和是集所有的，这是一个基因阻遏，定义如下：系统改写成紧凑的矩阵形式，我们得到 126 模型比较 同其他基因网络比较，布尔网络是一个相对粗糙的模拟，不过，布尔网络可提供一个框架，它描述了基因之间复杂的相互作用，并且从生物学意义上展示基因网络特征，如整体复杂性、自组织性、冗余性等等，而这些自然特征在权重矩阵模型中是不予考虑的。线性模型将基因间的相互作用都近似地看成一种线性关系，这与实际中基因间的相互作用关系往往是非线性的不相符，因此模型尚需进行改 进。布尔网络模型是一种粗放的定性方法
16、，而加权矩阵模型又是通过精细的数学分析来量化描述生物过程，贝叶斯网络模型则可以看成是这两种模型的一种折衷，贝叶斯网络方法也有其局限性：有向无环结构的假设与生物体的生命周期现象并不符合。其次，在贝叶斯网络的学习过程中，网络的结构会非常的复杂，计算量非常之大。微分方程模型是以其它基因表达水平和外部环境的因素组成的函数来描述基因表达的变化，可以充分模拟基因调控网络的动态行为。相比较其他模型，微分方程模型非常强大灵活，利于研究基因网络中的复杂关系。基因的网络分析是生命信息挖掘的重要手段之一,但目前在许多方面尚处于尝试和探索段。大量模型不断涌现，各种数学工具不断引进，这为网 络调控模型的构建创造了良好的
17、数学理论基础。随着基因数据的不断扩展以及数据质量的进一步提高，基因调控网络建模的准确性将得到进一步提高。随着后基因组学的不断发展，基因网络必然会在生命科学的研究中发挥巨大的作用。13 无源性研究现状无源性概念作为耗散性的特例广泛存在于物理学，应用数学以及力学等领域。它们在控制领域里的应用起源于Kalmna、Popov、Yakubocich和Willems等先驱们在超稳定性，正实性等方面开创性的工作。后经众多控制工作者的共同努力形成了系统的耗散性和无源性理论。刻画耗散系统、无源系统最重要的工具就是KYP引理，它是判断系统是否耗散、无源的充分必要条件。KYP引理最早是由Kalman，Popov以及
18、Yakubovieh提出来的，后由Anderson及其合作者将其推广到多变量系统。线性系统的KYP引理也称正实引理，它建立了频域框架下系统的正实性条件和时域里状态空间中相应正实性条件的等价关系。131 无源性、稳定性与最优性(1)内联系统的稳定性运用不同的方法，Popov和Zmaes分别提出了无源性理论中具有深远意义的定理:两个无源系统的并行内联或负反馈内联仍然是无源的。在一定的可探测条件下，无源系统是稳定的。(2)无源性与稳定性无源性与稳定性间的紧密联系可以从两方面得到体现:1、无源系统的KYP引理；2、无源系统的正定存储函数可以作为李雅普诺夫函数。无源系统并不是稳定系统的子集，无源系统加上
19、可探测性条件才会是稳定的。(3)无源性与最优性Kalman最早研究了线性系统最优性和无源性的关系，指出在可探测条件下，最优性和稳定性等价。Myolan将其结论推广到仿射非线性系统。最优控制可以给系统带来鲁棒性。通常最优控制问题归结为求解HJB方程，而这是件非常困难的事，目前仍没有解决该问题系统可行的方法。为了避免求解HJB方程，Freeman和Kokotovic提出了逆最优控制方法。在逆最优控制过程中，首先要给出镇定控制律和相应的Lypaunov函数，然后再设计性能指标使得其最优，即根据控制律反过来设计性能指标。无源性理论主要用来构造控制李雅普诺夫函数。保证系统稳定的控制李雅普诺夫函数正是HJ
20、B方程的解，这样保证系统的稳定性和性能指标的最优性。132 无源性理论应用由于利用了物理系统的结构特点，无源性方法在和其它控制技巧结合使用时，可以简化相应的控制方法。用无源化方法设计的非线性观测器可以减少观测器的参数，而且它也可以简化自适应控制，鲁棒控制，滑模控制，神经网络和模糊控制等。近年来，无源性理论在机械系统，电力系统，化工过程和机器人控制等方面获得了广泛的应用。Ortgea针对机器人系统讨论了基于无源性的控制问题；SiraRamirez针对电力系统及化工过程应用无源性理论研究系统的镇定问题。国内清华大学的研究人员也成功地将无源性理论应用于电力系统的控制研究中。相信这些实际系统研究的深入
21、和成熟，必将为无源性理论的研究注入新的活力。问题描述21系统模型建立变时滞基因调控网络： (1)其中，；，分别是时刻节点的mRNA和蛋白质的浓度；，表示mRNA和蛋白质的降解或稀释率；，是耦合矩阵。被定义如下:非线性函数表示转录蛋白的反馈调节，用HIll形式的单调函数，即，是Hill系数随时间变化的延迟和假定满足，并且，；，其中并是受节点阻碍的节点的集合。设是(1)平衡点，是下列方程的解： (2)为合宜起见，我们替换系统(1)预期的平衡点为原点。定义， ，可得 (3)其中，,。由于是单调饱和的递增函数，从的定义我们知道满足下面的条件： (4)考虑到基因的外部控制输入，可得GRNs为 (5) 其
22、中，表示时间时外部基因控制输入。，并认为为基因调控网络的输入，为基因调控网络的输出。22考虑不确定性的时滞系统模型考虑基因调控网络不确定性 (6)这里, , 和是参数不确定性，是分别满足和的随时间变化延迟。这里和是正的常数。，和。随时间变化的不确定矩阵，和定义如下： (7) 其中和是已知的适当维数的常数实矩阵。和是未知的时变矩阵满足 (8) 23无源性定义定义1 基因调控网络(5)被认为是无源的，如果存在一个标量，所有的即 (9) 在零初始条件下。24引理引理1.对于任何常数矩阵,标量，向量函数和有关的积分的明确界定，得.引理2.对于任何矩阵和向量下面等式成立引理3.设，和是适当尺寸的实矩阵且
23、满足。同时，任何标量24本章小结考虑到无源性，将原来的基因网络(1)加入外部输入得到基因网络(5)，之后再本章的第二节，考虑了基因网络的不确定性，使得讨论的基因网络鲁棒性更具体，更完善。无源性条件31变时滞基因网络的无源性条件这节中，我们分析了不确定的基因网络的无源性，随时间变化的延误(6)通过使用Lyapunov稳定性定理和无源性定义(9)。首先，我们考虑了系统的无源性(5)，也就是说，案件, , , 和(6)。在下面，我们用*表示对称的一部分。311 定理一定理1. 模型(5)在任何和时渐近鲁棒稳定且无源，如果分别存在正定矩阵，对角矩阵，满足： 312 定理一证明使用的GRNs模型：给出以
24、下的Lyapunov-Krasovskii泛函：为了显示无源性考虑到零初始条件，可推导通过我们使用的GRNs模型下面讨论泛函下面将问题分解证明各部分下面，推导要用的不等式根据引理1 考虑和关系，可以分类讨论(a) (b) 设S和T是正定矩阵，则因此，由引理2 注意到边界条件，任何对角正定矩阵U，可以得到 综合上述推导得到下面的不等式其中这是很容易看到，当且仅当和，从Lyapunov-Krasovskii泛稳定性定理和无源性定义(9)，变时滞GRNs(5)渐近稳定和无源。32不确定变时滞基因网络的无源性条件321 定理二无源网络(6)在任何和时渐近鲁棒稳定，如果存在正定矩阵，对角矩阵，标量，这样
25、得到下面的线性矩阵不等式： 322 定理二证明以相同的Lyapunov泛函，证明定理1，分别替换，和为和。写为 然后使用上述不等式 这里这是很容易看到了当且仅当和。因此，不确定变时滞基因网络GRNs(6)是渐近稳定和无源。33本章小结本章使用引理1和引理2的到了基因调控网络的渐进鲁棒稳定性和无源性的条件，如果给定初始条件和常数，通过Matlab的LMI工具箱进行检验可以得到具体的基因调控网络是否渐进鲁棒稳定性和无源性。三节点变时滞基因网络研究仿真repressilator动态理论上已经预测，并在大肠杆菌中的实验调查7。在该repressilater是一个循环的负反馈回路，包括三个阻遏基因（la
26、cl, tetR和cl）和它们的催化剂。系统的动力学描述如下： 这里。和是第三阻遏mRNA和蛋白的浓度，表示蛋白质的衰变率，mRNA降解率的比例。41 三节点变时滞基因网络模型考虑到转录的时间延迟，我们从一些参数调整成矢量改写上述方程： (27)其中 ,在上述方程中,，其中是Hill系数，且。这里很容易检验的导数的最大值小于和，。42 用LMI工具箱检验鲁棒稳定性程序421 程序内容zzjNetwork1Stab.mfunction zzjNetwork1Stab%用LMI工具箱检验变时滞基因调控网络1的鲁棒稳定性clear allclose all%常数A=3 0 0;0 3 0;0 0 3
27、;C=2.5 0 0;0 2.5 0;0 0 2.5;D=0.8 0 0;0 0.8 0;0 0 0.8;B=2.5,2.5,2.5;W=0 0 -2.5;-2.5 0 0;0 -2.5 0;delt=1/8;tau0=3/8;sigma0=3/16;eta=1/16;K=0.65 0 0;0 0.65 0;0 0 0.65;%初始化setlmis()%定义变量Q1=lmivar(1,3,1);Q2=lmivar(1,3,1);Q3=lmivar(1,3,1);Q4=lmivar(1,3,1);R1=lmivar(1,3,1);R2=lmivar(1,3,1);R3=lmivar(1,3,1)
28、;R4=lmivar(1,3,1);X=lmivar(1,3,1);Y=lmivar(1,3,1);U=lmivar(1,3,1);%添加项%不等式1lmiterm(1 1 1 Q1,-2,A); %-2*Q1*Almiterm(1 1 1 Q2,1,1); %Q2lmiterm(1 1 1 Q3,-1/tau0,1); %-(1/tau0)*Q3lmiterm(1 1 1 Q1,1,1); %Q1 lmiterm(1 1 2 X,-A,1,s); %-AX lmiterm(1 1 3 Q3,1/tau0,1,s); %(1/tau0)*Q3 lmiterm(1 2 2 Q3,tau0,1);
29、 %tau0*Q3lmiterm(1 2 2 Q4,2*delt,1); %2*delt*Q4lmiterm(1 2 2 X,-1,1); %-X lmiterm(1 3 3 Q2,-1,1); %-Q2lmiterm(1 3 3 Q4,-1/delt,1); %-(1/delt)Q4lmiterm(1 3 3 Q3,-1/tau0,1); %-(1/tau0)*Q3 lmiterm(1 3 4 Q4,1/delt,1,s); %(1/delt)*Q4 lmiterm(1 4 4 Q4,-1/delt,1); %-(1/delt)*Q4lmiterm(1 4 4 R1,D,D); %D*R1*
30、Dlmiterm(1 4 4 Y,D,D); %D*Y*D%不等式2lmiterm(2 1 1 R1,-2,C); %-2*R1*Clmiterm(2 1 1 R2,1,1); %R2lmiterm(2 1 1 R3,-1/sigma0,1); %-(1/sigma0)*R3lmiterm(2 1 1 R1,1,1); %R1 lmiterm(2 1 2 Y,-C,1,s); %-C*Y lmiterm(2 1 3 R3,1/sigma0,1,s); %(1/sigma0)*R3 lmiterm(2 2 2 R3,sigma0,1); %sigma0*R3lmiterm(2 2 2 R4,2*
31、eta,1); %2*eta*R4lmiterm(2 2 2 Y,-1,1); %-Y lmiterm(2 3 3 R2,-1,1); %-R2lmiterm(2 3 3 R3,-1/sigma0,1); %-(1/sigma0)*R3lmiterm(2 3 3 R4,-1/eta,1); %-(1/eta)*R4 lmiterm(2 3 4 R4,1/eta,1,s); %(1/eta)*R4 lmiterm(2 4 4 R4,-1/eta,1); %-(1/eta)*R4 lmiterm(2 4 5 U,1,1,s); %U lmiterm(2 5 5 Q1,W,W); %W*Q1*Wlm
32、iterm(2 5 5 X,W,W); %W*X*Wlmiterm(2 5 5 U,-2,inv(K); %-2*U*inv(K)%获取LMI系统描述lmisys=getlmis;tmin,xfeas=feasp(lmisys)%tmin =-0.00390，表示可行，用pmat=dec2mat(lmisys,xfeas,P)来得到一个可行的P矩阵Qd1=dec2mat(lmisys,xfeas,Q1)Qd2=dec2mat(lmisys,xfeas,Q2)Qd3=dec2mat(lmisys,xfeas,Q3)Qd4=dec2mat(lmisys,xfeas,Q4)Rd1=dec2mat(l
33、misys,xfeas,R1)Rd2=dec2mat(lmisys,xfeas,R2)Rd3=dec2mat(lmisys,xfeas,R3)Rd4=dec2mat(lmisys,xfeas,R4)Xd=dec2mat(lmisys,xfeas,X)Yd=dec2mat(lmisys,xfeas,Y)Ud=dec2mat(lmisys,xfeas,U)422 结果根据定理1，通过Matlab LMI工具箱求解线性矩阵不等式(13)得到和的值%Network1稳定检验可行性和可行的矩阵tmin =-0.00390 %可行Qd1 =4.2872,0,0；0,4.2872,0；0,0,4.2872;
34、Qd2 =14.6935,0,0；0,14.6935,0；0,0,14.6935;Qd3 =-0.8856,0,0；0,-0.8856,0；0,0,-0.8856;Qd4 =1.4802,0,0；0,1.4802,0；0,0,1.4802;Rd1 =7.5094,0,0；0,7.5094,0；0,0,7.5094;Rd2 =20.9242,0,0；0,20.9242,0；0,0,20.9242;Rd3 =0.3054,0,0；0,0.3054,0；0,0,0.3054;Rd4 =1.9380,0,0；0,1.9380,0；0,0,1.9380;Xd =0.1229,0,0；0,0.1229,0
35、；0,0,0.1229;Yd =0.6773,0,0；0,0.6773,0；0,0,0.6773;Ud =17.2935,0,0；0,17.2935,0；0,0,17.2935;43 用ddesd算法仿真基因调控网络1过程ddesd算法是仿真变时滞偏分方程的一种算法431 展开矩阵方程得到微分方程组程序zzjgetNetwork1Eq.mfunction zzjgetNetwork1Eq%获得基因网络1微分方程组clear allclose all%常数A=3 0 0;0 3 0;0 0 3;C=2.5 0 0;0 2.5 0;0 0 2.5;D=0.8 0 0;0 0.8 0;0 0 0.8
36、;B=2.5,2.5,2.5;W=0 0 -2.5;-2.5 0 0;0 -2.5 0;%定义变量syms t y1 y2 y3 y4 y5 y6;syms yy1 yy2 yy3 yy4 yy5 yy6;M=y1,y2,y3;P=y4,y5,y6;Mdelt=yy1,yy2,yy3;fP=yy4*y4/(1+yy4*y4),yy5*yy5/(1+yy5*yy5),yy6*yy6/(1+yy6*yy6);-A*M+W*fP+B-C*P+D*M结果：其中conj（）是共轭复数函数%Network1微分方程dy(1)=5/2 - (5*conj(yy6)2)/(2*(conj(yy6)2 + 1)
37、 - 3*conj(y1)dy(2)=5/2 - (5*conj(y4)*conj(yy4)/(2*(conj(y4)*conj(yy4) + 1) - 3*conj(y2)dy(3)=5/2 - (5*conj(yy5)2)/(2*(conj(yy5)2 + 1) - 3*conj(y3)dy(4)=(4*conj(y1)/5 - (5*conj(y4)/2dy(5)=(4*conj(y2)/5 - (5*conj(y5)/2dy(6)=(4*conj(y3)/5 - (5*conj(y6)/2432 仿真基因调控网络1仿真程序zzjNetwork1.mfunction zzjNetwork
38、1%用ddesd算法仿真变时滞影响下的基因调控网络1过程 t0 = 0; %设置时间初始值tfinal = 20; %设置时间终值tspan = t0, tfinal; options=odeset(RelTol,0.00001); %精度设置%-sol = ddesd(zzj1Eq,zzj1delay,zzj1hist,tspan,options); %作图1figure(1);plot(sol.x,sol.y(1,:),g,sol.x,sol.y(2,:),b,sol.x,sol.y(3,:),r);xlabel(time t);ylabel(solution y);title(mRNA浓
39、度);hold all;figure(2);plot(sol.x,sol.y(4,:),g,sol.x,sol.y(5,:),b,sol.x,sol.y(6,:),r);xlabel(time t);ylabel(solution y);title(蛋白质浓度);hold all; %-sol = ddesd(zzj1Eq,zzj1delay,zzj2hist,tspan,options); %作图2figure(1);plot(sol.x,sol.y(1,:),g,sol.x,sol.y(2,:),b,sol.x,sol.y(3,:),r);figure(2);plot(sol.x,sol.
40、y(4,:),g,sol.x,sol.y(5,:),b,sol.x,sol.y(6,:),r);hold off %-函数集-function v = zzj1hist(t)v=0 0.2 0.4 0 0.1 0.2; %给出参数y初始值 function v = zzj2hist(t)v=1.5 1.1 0.8 0.5 0.4 0.3; %给出参数y初始值 function d = zzj1delay(t,y)%延时函数矩阵d = t-(abs(sin(t)+1)/4;t-(abs(sin(t)+1)/4;t-(abs(sin(t)+1)/4;t-(abs(cos(t)+1)/8;t-(ab
41、s(cos(t)+1)/8;t-(abs(cos(t)+1)/8; function dy=zzj1Eq(t,y,z)dy=zeros(6,1);%延迟设置yy1=z(:,1);yy2=z(:,2);yy3=z(:,3);yy4=z(:,4);yy5=z(:,5);yy6=z(:,6);%微分方程组dy(1)=-3*y(1)-2.5*yy6(6)*yy6(6)/(1+yy6(6)*yy6(6)+2.5;dy(4)=-2.5*y(4)+0.8*yy1(1);dy(2)=-3*y(2)-2.5*yy4(4)*yy4(4)/(1+yy4(4)*yy4(4)+2.5;dy(5)=-2.5*y(5)+0.8*yy2(2);dy(3)=-3*y(3)-2.5*yy5(5)*yy5(5)/(1+yy5(5)*yy5(5)+2.5;dy(6)=-2.5*y(6)+0.8*yy3(3);44 用ddesd算法仿真基因调控网络1过程结果下图图显示了变量和分别在2个不同的初始函数下的轨迹。因此，伴随有时变延迟的GRN（27）是渐近稳定的。431 mRNA浓度时间t432 蛋白质浓度时间t45 本章小结本章使用定理1，说明了定理的实用性。五节点变时滞基因网络研究仿真51 三节点变时滞基因网络模型考虑下面范数有界变时滞基因调控网络的不确定性：
