基于连续-离散扩展卡尔曼建模的p53-Mdm2调控关系研究

姜 莉 李玉榕*

1(福州大学电气工程与自动化学院，福州 350116)2(福州大学福建省医疗器械和医药技术重点实验室，福州 350116)

基于连续-离散扩展卡尔曼建模的p53-Mdm2调控关系研究

姜 莉1,2李玉榕1,2*

1(福州大学电气工程与自动化学院，福州 350116)2(福州大学福建省医疗器械和医药技术重点实验室，福州 350116)

生物医学研究揭示，肿瘤抑制基因p53的表达水平与肿瘤的形成相关。p53相关的信号转导网络研究为揭示癌症、肿瘤的发病机制和寻找其治疗方法提供新思路，而p53和癌基因Mdm2的负反馈关系是该信号转导网络的核心部分。因此，研究p53-Mdm2的调控关系意义重大。基于电离辐射后人类白血病细胞中测得的时间序列基因表达数据，利用连续-离散的扩展卡尔曼(EKF)算法，建立带有延迟特性的非线性动态连续随机生物网络模型，模拟电离辐射后的p53-Mdm2动态调控过程。同时，对所建立网络的准确性进行验证，可得模型的误差率为0.85%。结果表明：所提出的算法是收敛的，可以预测基因的表达水平，准确地模拟p53-Mdm2调控网络受到电离辐射后的衰减震荡响应过程，该过程和生物实验结果一致。所提出的算法为生物学实验建模提供一种有效的方法，并可以为网络系统动力学特性的研究提供条件。

p53-Mdm2网络；连续-离散扩展卡尔曼(EKF)算法；时间序列数据；延迟

引言

p53是重要的肿瘤抑制因子，在细胞周期抑制、细胞凋亡、生长停滞、抑制血管生成等过程中具有重大的作用[1]，且与肿瘤的形成和治疗有密切的关系。研究表明，许多癌症中都发现了p53功能异常，超过50%的癌症是由p53基因突变引起的[2]。在癌症理疗中，人工恢复p53的活性，可以使癌细胞停止增长或凋亡，来达到治疗的目的[3-4]。因此，越来越多的学者开始关注p53的研究。

p53对其下游多个基因的表达进行调控，其中p53和Mdm2的负反馈调节可以看成p53整个转导通道的核心，p53和Mdm2的表达失衡会引起器官或组织中的病理变化[3]。研究p53和Mdm2蛋白功能及其之间的相互调节机制不仅可以找到新的癌症治疗方法，而且在某些疾病的预防方面也有极其重要的作用。在正常细胞中p53蛋白的含量很少，在细胞受到外界刺激时，如电离辐射(ionization radiation，IR)、紫外线辐射(ultraviolet radiation，UV)等，会引起DNA的损伤，随后激活p53，使得p53表达增多，从而引起一系列复杂的反应。该过程中p53-Mdm2调节过程是一个复杂且精细的过程，而生物实验在这方面存在着很多局限。例如：反应过程中的微观数据很难测量；可供实验的样本数量有限；某些反应过快或过慢，难以实现测量；依靠实验寻找靶点过于盲目等。而数学建模不仅可以模拟得到生物实验中由于技术原因不能测得的数据，而且可以从微观上来分析生物反应过程中的可能响应，为生物实验提供更多的可能靶点。

出于技术、经费、人力的考虑，有学者提出通过建立数学模型来模拟p53-Mdm2负反馈的网络模型[3-7]。对生物网络的建模能够从机制上研究生命系统的一些现象，如细胞凋亡，生长停滞等[8]。这些学者的做法都是先建立能够符合生物特性的模型，然后通过数值分析方法，对模型的动态特性进行分析。这类研究在确定模型参数的时候，一般是比较盲目地去尝试，或者根据经验来取得，而且这些模型只能表示一类动态网络，并不能够针对某一个器官或者组织给出特定的模型参数。而随着基因芯片技术和高通量技术的发展，可以结合生物学实验和数学模型进行研究。通过实验在有限的时间点上检测细胞受到电离辐射后p53和Mdm2的表达量，而后利用这些时间序列的表达数据对模型进行状态和参数辨识，从而确定模型参数。

扩展卡尔曼(extended kalman Filter，EKF)算法因其能辨识非线性模型而被广泛应用。Wang等建立基因调控网络的状态空间模型[9]，通过EKF算法来同时辨识网络的参数和状态，建立动态的基因网络，并对网络的稳定性能进行分析[9]。Xiong等利用EKF线性化的特点，提出一种基于鲁棒性状态估计的算法，来对生物网络的结构进行辨识[10]。Zhang等提出利用一种分数的扩展卡尔曼算法，对于建立的非线性基因调控网络进行建模[11]。以上建模研究中均将生物网络看作离散系统，而生物的响应过程本质上是一个连续的过程，若能够考虑到其连续的性质，则会更加符合生物网络的特性，得到更准确的生化反应动态模型。

本研究首先根据生物调控过程的连续性、延迟性、随机性和非线性特性，建立一个连续的带有延迟的非线性动态随机生物网络模型，来模拟电离辐射后p53-Mdm2的响应过程。而后利用生物学实验得到的电离辐射后测得的时间序列表达数据，基于连续-离散的EKF算法[12]，同时辨识网络的参数和基因表达水平。最后，对网络的准确性进行估算，并对所建立模型的响应过程进行分析。本研究提出的算法具有通用性，可以得到特定网络的更加准确的模型参数，为进一步理解p53调节通道提供了可能，从而为癌症治疗提供更加准确的靶点。

1 材料与方法

基于时间序列的基因表达数据，综合考虑生化反应过程的连续性、时滞性、随机性和非线性，首先提出符合生物特性的微分模型；然后利用连续-离散扩展卡尔曼算法对于模型的状态和参数进行联合辨识；最后利用误差表达式对于已经建立的模型进行评价。

1.1 数据样本

采用的源数据取自文献[13]中：在室温下，对人类白血病细胞(MOTL4)进行5Gy的γ辐射后，分别获得细胞内0、2、4、6、8、10、12 h的基因表达水平；为了减少测量中的误差，独立进行了3次数据测量过程，其均值和标准偏差如图1、2所示。本研究选取3组数据的平均值来进行建模。

图1 p53的平均值和标准差Fig.1 Mean and standard deviation of p53

图2 Mdm2的平均值和标准差Fig.2 Mean and standard deviation of Mdm2

1.2 p53-Mdm2网络模型

p53-Mdm2负反馈网络尚有许多细节需要完善，为了不失一般性，本研究采用简化模型,不考虑mRNA的作用，即把转录和翻译合成一项来建模。生物网络的模型不仅要求能够符合生物反应机理，更要符合生物网络的调节特性[3]。而基因在调控的过程中，即转录和翻译过程中，都会发生延迟的现象[2]。p53-Mdm2动态调控网络中，p53对于Mdm2的调控也存在时滞，若不考虑到延迟的特性，必然会影响基因网络的稳定性，甚至导致模型错误的预测响应过程。因此，本研究在文献[5]中的模型基础上，考虑到生物过程具有延迟的特性[2,7]，把转录和翻译过程中的延迟整合到Mdm2表达式中的一个项中来表示。p53、Mdm2分别表示相应的蛋白水平，其调控网络表示如下：

(1)

式中：第1项r1是p53蛋白本底的合成速率；第2项表示除了Mdm2外其他原因引起p53的降解，d1是本底的降解速率；由于p53转导信号网络具有非常多的正反馈回路[14]，第3项描述了p53系统中的正反馈效应，正反馈可以被认为是一种广义的“自诱导”，模型中用一个Hill函数来描述，vp表示正反馈导致的本底生成速率，k1是Hill函数的常数；第4项表示Mdm2诱导的p53的降解作用；kmax表示Mdm2对p53负调控的降解速率，kp是p53和Mdm2之间的络合物离解常数[3]。

(2)

式中：第1项r2表示Mdm2的合成速率；第2项表示Mdm2的本底降解，d2表示本底降解速率；第3项表示p53诱导Mdm2的生成，vm表示p53引起的Mdm2的增长速率，τ表示调控过程存在的滞后，k2是Hill函数的常数；第4项表示Mdm2介导的自降解作用，mdeg表示Mdm2的自降解的速率，km表示Mdm2离解常数。

1.3 连续-离散的EKF算法建模

1.3.1 连续-离散的EKF算法

EKF算法在每个时间步长用一个线性化的过程来把非线性系统近似为时变的线性系统，再用传统的卡尔曼算法对时变的线性系统进行滤波，是卡尔曼算法在非线性系统中的应用。EKF是基于最小方差的次优算法，并不是最优的算法，但其效果良好并成为处理非线性滤波问题中最经典、应用最广泛的方法[9]。详细介绍EKF算法可参考文献[9，15]。EKF算法的收敛性讨论可以参考文献[16-17]。EKF算法已经在很多时候应用在参数识别中[9-11]，一般使用离散的EKF算法。本研究根据模型连续的特性，以及使用的时间序列基因表达数据离散的特点，采用连续-离散的EKF算法来对模型中的参数和p53、Mdm2蛋白的表达水平进行辨识。此外，考虑到基因调控过程会遇到各种噪声，如附近基因表达对于本网络中基因表达信号测量的干扰、由基因芯片放置位置不同而产生的噪声、测量时温度变化引起的噪声、基因芯片使用时产生的测量噪声等，在网络中加入噪声的影响。具体过程如下：

将p53和Mdm2的表达水平分别用x1、x2表示，并在模型中加入噪声，于是式(1)、(2)可以写成

(3)

(4)

式中，w1(t)、w2(t)为互不相关的均值为零的高斯白噪声。

令x(k)=[x1(t)x2(t)]T，则式(3)、(4)可以简化写为

(5)

为了使用EKF算法来辨识参数和基因表达水平，把需要辨识的参数扩展到系统的状态变量中。令

最终用连续-离散EKF算法辨识参数的模型写为

式中，y(t)为基因芯片含噪声的测量值，V(t)、W(t)为互不相关的均值为零的高斯白噪声，其协方差矩阵分别为Q(t)、R(t)，k为采样点。

1.3.2 算法的具体迭代步骤

假设采样点的个数为N。

步骤1：k=0时，设置初始值。

P(0,0)=E[(X(0)-X0)(X(0)-X0)T]=PX0

步骤2：k=1，2，3，4，…，N时，计算如下：

1)时间更新(预测)，其中t∈[tk-1,tk]。

状态时间更新

(9)

(10)

误差协方差时间更新

(11)

P(k|k-1)=P(tk)

(12)

2)测量更新(矫正) 。

计算卡尔曼增益矩阵，有

Kk=P(k|k-1)CkT[CkP(k|k-1)CkT+Rk]-1

(13)

用测量值更新估计值，有

(14)

误差协方差更新，有

P(k|k)=(I-KkCk)P(k|k-1)

(15)

其中

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

1.3.3 算法的具体实现步骤

不断重复步骤1～4，直至采样点全部用完。

1.4 模型评价

为了能够定量化地衡量辨识出的模型质量，引入误差指标，即观测到的值与模型实际预测值之间的误差百分比，模型误差算式[9,18]具体表示为

(25)

式中：l指的是所有基因的数目；s是观测点的数量，即采样的点数；yck是指c基因在k时刻的实际表达水平。

2 结果

本算法是基于Matlab软件。为探索电离辐射后网络的动态过程，实验数据点数越多，则辨识的结果越准确。而本文第1.1节中的原始数据只有7个时间点，使用三次样条差值方法可以得到足够多的点，以避免出现过度拟合的情况[19]。所以，通过三次样条插值方法取25个采样点，每30 min一次。此外，由于实验测得的数据过少，若辨识全部的参数难免会出现误差较大。而文献[5]中指出，r1(p53的合成速率)、vm(负反馈强度)对于系统动力学影响较大，故本研究辨识式(3)、(4)中的参数r1、vm值，其余参数值根据文献[5]设置(τ取为5 min)。

2.1 辨识结果

利用连续-离散扩展卡尔曼算法和图1、2中的实验数据，对p53-Mdm2模型进行辨识，估计出模型参数和状态的值，分别如图3、4所示。图5、6所示分别为连续-离散卡尔曼算法辨识的参数和状态的估计误差方差。

图3 模型参数r1,vm估计值Fig.3 Estimated parameters of r1,vm

图4 p53,Mdm2表达水平的估计值和实验值Fig.4 Actual and estimated expression level of p53,Mdm2

图5 参数r1,vm的估计误差方差Fig.5 Variances of estimated parameters r1,vm

图6 状态p53、Mdm2的估计误差方差Fig.6 Variances of estimated states p53,Mdm2

从图3可以看出，经过迭代滤波，参数r1和vm稳定收敛，确定出r1=0.794 7，vm=0.124 7。从图4中看出，p53和Mdm2的预测值和实验值之间误差较小。同时，图5、6显示参数和状态的估计误差方差都很小，除了开始时期，其余时刻接近于零，这意味着该方法能准确地从时间序列数据中估计出模型的状态和参数(见式(3)、(4))。

2.2 评价结果

表1列出了在不同的纯滞后参数τ下所建立模型对应的误差率。

表1 不同τ的模型误差率Tab.1 Error ratio of the model in different τ

由表1可看出，模型取不同纯滞后时间τ,误差率均较低，利用本算法辨识出的模型是准确的。

3 讨论

在研究中，基于连续-离散卡尔曼算法辨识参数和状态，从而建立调控模型。基于该模型，可以预测12 h之后p53以及Mdm2的表达量。图7表示了p53和Mdm2前12 h的估计值和之后的预测值，左边的坐标轴表示p53的表达水平，右边坐标系表示Mdm2的表达水平。

图7 p53和Mdm2的响应过程Fig.7 The response process of p53 and Mdm2

从p53和Mdm2的响应过程中可以看出，白血病细胞在受到电离辐射后，p53比Mdm2先达到峰值。p53第一个脉冲宽度大约为250 min，达到第一个脉冲峰值为200 min左右；而Mdm2的第一个脉冲宽度大约为600 min，达到第一个脉冲峰值为400 min。生物实验的结果表明[20]，细胞受到电离辐射时，会进行振荡，振荡脉冲的宽度为350±160 min，第一个脉冲峰值对应的时间相差很大，大约在损伤后的360±240 min。可见本模型分析所得到的结果很好地符合了生物学实验结果。

此外，从图7中可以看出，电离辐射后，p53的浓度增加，促进Mdm2的生成，当Mdm2的浓度增加到一定程度后，反过来抑制Mdm2的生成。并且可以看出，p53和Mdm2的响应过程是一个衰减震荡的过程，最后p53和Mdm2的浓度基本不变，达到了一个平衡的状态。已有研究表明[21-23]，电离辐射会引起DNA损伤，在轻度的DNA损伤中，p53表现出脉冲动力学特性，这一模式具有很好的灵活性，既有利于细胞周期停滞期间损伤的有效恢复，又可避免轻度损伤细胞的过度凋亡。在达到平衡状态后，说明此时DNA的损伤已经被修复。此外，震荡的过程即为DNA的修复过程，p53发出第一个脉冲后，等待看损伤是否被正确修复，若没有就会产生第二个脉冲，直到损伤有效修复，p53的浓度不再变化[5,22]。从图7中可以看出，大约800 min细胞完成DNA的修复过程。

本算法去除了离散扩展卡尔曼算法不符合生物过程连续性的缺点，不仅考虑到了实验得到的微阵列数据的离散性，而且考虑到了生物调控的连续性。建立模型后，通过图5、6和表1对辨识结果进行评价，说明建立模型的准确性较高。而图7进一步表明，该辨识的模型可以对之后的基因表达进行预测，且辨识得到的结果和预测结果均具有实际意义，最后得到的调控过程基本符合生物学实验结果，为广大系统动力学研究工作者提供了一种新的建模方式。

4 结论

本研究讨论了p53-MDM2网络在受到电离辐射后的响应过程。首先，根据生物过程延迟和非线性的特性，并依据生物学实验，建立生物网络的非线性随机动态模型。一般模型中的参数依据经验试凑确定，而且无法针对某个个体或某个组织、器官的特异性以及不同疾病模型来确定。本研究根据生物响应连续的特性，提出连续-离散的EKF算法，利用实验获得的电离辐射后测的时间序列基因表达数据，对模型中的参数进行辨识，进而得到该个体或者该组织、器官的模型参数。然后，根据得到的模型，不仅可以得到p53-Mdm2网络在测量时间内的表达水平及动态响应，而且可以定量地预测之后任意时刻的基因表达量及动态响应。

本算法可以较准确地估计出模型的参数，这样在研究某一疾病或者肿瘤治疗的时候，就可以准确得到该疾病患者癌细胞或肿瘤细胞中p53的含量以及响应过程，为进一步决定电离辐射治疗的量提供参考，为医疗工作者提供指导。存在的问题是，由于技术的限制，研究中使用的原始数据的点数太少。这样，为了考虑参数的准确性，不得不只选择比较重要的参数进行辨识，其他参数用经验取值。随着基因芯片技术的不断发展，肯定可以得到更多时间点的数据，这样，把模型中的所有参数辨识出来成为可能，并且提高其准确性，为研究网络的动力学特性提供基础，进而可更精确地选择特定癌症化疗中辐射使用量以及更准确的治疗靶点。

[1] Taberlay PC, Jones PA. DNA methylation and cancer [J]. Oncogene, 2012, 21(35): 89-95.

[2] Yang Y, Lee KS, Xiang C, et al. Biological mechanisms revealed by a mathematical model for p53-Mdm2 core regulation [J]. Iet Systems Biology, 2009, 3(4): 229-238.

[3] Voropaeva OF, Shokin YI, Nepomnyaschchikh LM, et al. Mathematical modeling of the tumor markers network [C] // International Conference on Biomedical Engineering and Computational Technologies. Novosibirsk: IEEE, 2015: 225-229.

[4] Wang Shaomeng, Sun Wei, Zhao Yujun, et al. SAR405838: an optimized inhibitor of MDM2-p53 interaction that induces complete and durable tumor regression [J]. Cancer Research, 2014, 74(20): 5855-5865.

[5] 孙延哲. P53信号转导网络动力学模型研究 [D]. 南京: 南京大学, 2014.

[6] 毕远宏, 杨卓琴, 何小燕. Mdm2生成速率调控的p53-Mdm2振子的全局动力学和稳定性 [J]. 物理学报, 2016, 65(2): 360-367.

[7] Mihalas GI, Simon Z, Balea G, et al. Possible oscillatory behavior in p53-Mdm2 interaction computer simulation [J]. Journal of Biological Systems, 2011, 08(1): 21-29.

[8] Lavoie H, Hogues H, Mallick J, et al. Evolutionary tinkering with conserved components of a transcriptional regulatory network [J]. PLoS Biology, 2010, 8(3): e1000329.

[9] Wang Zidong, Liu Xiaohui, Liu Yurong, et al. An extended Kalman filtering approach to modeling nonlinear dynamic gene regulatory networks via short gene expression time series [J]. IEEE/ACM Transactions on Computational Biology & Bioinformatics, 2009, 6(3): 410-419.

[10] Xiong Jie, Zhou Tong. Structure identification for gene regulatory networks via linearization and robust state estimation [J]. Automatica, 2013, 50(11): 2765-2776.

[11] Zhang Yongqing, Pu Yifei, Zhang Haisen, et al. An extended fractional Kalman filter for inferring gene regulatory networks using time-series data [J]. Chemometrics & Intelligent Laboratory Systems, 2014, 138: 57-63.

[12] Kulikova MV, Kulikov GY. Adaptive ODE solvers in extended Kalman filtering algorithms [J]. Journal of Computational & Applied Mathematics, 2014, 262(10): 205-216.

[13] Barenco M, Tomescu D, Brewer D, et al. Ranked prediction of p53 targets using hidden variable dynamic modeling [J]. Genome Biology, 2006, 7(3): R25.

[14] Vilborg A, Glahder JA, Wilhelm MT, et al. The p53 target Wig-1 regulates p53 mRNA stability through an AU-rich element [J]. Proceedings of the National Academy of Sciences of the United States of America, 2009, 106(37): 15756-15761.

[15] Hatami E, Salarieh H, Vosoughi N. Design of a fault tolerated intelligent control system for a nuclear reactor power control: using extended Kalman filter [J]. Journal of Process Control, 2014, 24(7): 1076-1084.

[16] Huang S, Dissanayake G. Convergence and consistency analysis for extended Kalman filter based SLAM [J]. IEEE Transactions on Robotics, 2007, 23(5): 1036-1049.

[17] Krener AJ. The convergence of the extended Kalman filter [J]. Mathematics, 2003, 286: 173-182.

[18] Zeng Nianyin, Wang Zidong, Li Yurong, et al. Inference of nonlinear state-space models for Sandwich-Type lateral flow immunoassay using extended Kalman filtering [J]. IEEE Transactions on Biomedical Engineering, 2011, 58(7): 1959-1966.

[19] Tu CT, Chen BS. On the increase in network robustness and decrease in network response ability during the aging process: a systems biology approach via microarray data [J]. IEEE/ACM Transactions on Computational Biology & Bioinformatics, 2013, 10(2): 468-480.

[20] Lahav G, Rosenfeld N, Sigal A, et al. Dynamics of the p53-Mdm2 feedback loop in individual cells [J]. Nature Genetics, 2004, 36(2):147-150.

[21] Bar-Or RL, Maya R, Segel LA, et al. Generation of oscillations by the p53-Mdm2 feedback loop: a theoretical and experimental study [J]. Proceedings of the National Academy of Sciences of the United States of America, 2000, 97(21): 11250-11255.

[22] 张小鹏, 刘锋, 王炜. p53信号网络的非线性动力学研究 [J]. 中国科学:物理学力学天文学, 2014(12): 1311-1318.

[23] Lahav G. The strength of indecisiveness: oscillatory behavior for better cell fate determination [J]. Journal of Health Services Research & Policy, 2014, 19(1): 1-2.

Research of p53-Mdm2 Regulatory Networks Based on Modeling by Continuous-Discrete Extended Kalman Filter

Jiang Li1,2Li Yurong1,2*

1(CollegeofElectricalEngineeringandAutomation,FuzhouUniversity,Fuzhou350116,China)2(FujianKeyLaboratoryofMedicalInstrumentation&PharmaceuticalTechnology,FuzhouUniversity,Fuzhou350116,China)

The studies of biomedicine show that the expression level of p53 (tumor suppressor gene) is related with the tumor forming. Researches of the p53 related signal transduction networks would provide new ideas for revealing the pathogenesis of cancer or tumor and looking for the treatment method. As the p53-Mdm2 negative feedback network plays core roles, the study of the p53-Mdm2 regulation net work is of great significance. In this paper, based on the gene expression time series data of human leukemia cells after ionizing radiation, the nonlinear dynamic continuous random biological networks with time delay were established, then continuous-discrete extended Kalman filter (EKF) algorithm was used to simulate the dynamic regulation of p53 and Mdm2. Meanwhile, the accuracy of the model established in this paper was validated and we found that the error rate of the model was only 0.85%. Results showed that the algorithm thatis convergent could predict gene expression level at any time and simulate accurately the damped oscillation process of p53-Mdm2 regulatory networks after ionizing radiation. The response process was consistent with biological experiment results. The algorithm proposed in this paper provided an effective method for biological experiment modeling, and supplied the foundation for the research of system dynamics.

p53-Mdm2 network; continuous-discrete extended Kalman filter (EKF) algorithm; time-series data; delay

10.3969/j.issn.0258-8021. 2017. 02.008

2016-08-05， 录用日期:2016-10-13

国家自然科学基金(61403319);福建省科技厅国际合作项目(2015I0003);福建省教育厅科技项目(JK2014001)

R318

A

0258-8021(2017) 02-0180-07

*通信作者(Corresponding author)， E-mail: liyurong@fzu.edu.cn