ABSTRACT

For the early diagnosis of complex diseases, traditional methods require many patient samples to obtain statistical indicators. In real life, people do not go to hospitals for medical examination frequently, so it is difficult to obtain a large number of patients data, and only the result of individual single sampling has been got, which requires a single sample based algorithm to detect critical points in the development of complex diseases. Based on the theory of dynamic system and dynamic network biomarkers, a new algorithm is developed by using new individual samples combined with local biomolecular networks and the concept of entropy. The algorithm is applied to numerical simulation data and a real disease data. In these applications, the algorithm can detect the critical point in time and accurately. The new algorithm developed in this paper is stable and efficient. It can help medical workers to further understand the dynamic biomolecular mechanism of complex diseases, and can more quickly propose a reasonable treatment plan.

Keywords:Single-Sample, Dynamical System, Dynamic Biomarkers, Local Network Entropy

单样本网络熵与某些复杂疾病恶性 突变的预警

李明亮，刘锐

华南理工大学数学学院，广东 广州

收稿日期：2019年1月2日；录用日期：2019年1月17日；发布日期：2019年1月24日

摘 要

关于复杂疾病的早期诊断，传统方法需要许多病人样本来获取统计指标。而现实生活中人们不会频繁出入医院进行体检，故难以获得大量的患者数据，只有对个体单次取样的结果，这就需要提出基于单样本的探测复杂疾病发展过程中临界点的算法。本文基于动力系统和动态网络生物标志物的理论，利用个体新样本结合局部生物分子网络以及熵的概念开发了一个新算法。将该算法应用于数值模拟数据和一个真实的疾病数据，在这些应用中，该算法能及时准确地探测出临界点。本文开发的新算法具有稳定性和高效性，能帮助医学工作者进一步了解复杂疾病的动态生物分子机制，能更快速地提出合理的治疗方案。

关键词 :单样本，动力系统，动态生物标志物，局部网络熵

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

许多复杂系统的演变过程都存在着临界点，比如种群生物数量系统、金融市场、气候系统、癫痫发作等系统中都存在着在某个很短的时间段内系统从一个稳定状态迅速转变到另一个稳定状态的现象。在临界点时期，复杂系统的鲁棒性急剧降低，状态转移风险大幅升高 [1] 。关于复杂疾病的临界点，相关的研究已经取得了许多有效的成果，比如最近有国内学者提出了利用动态网络生物标志物进行复杂疾病临界点探测的方法 [2] [3] [4] 。一般而言，可以将复杂疾病(比如肝癌、糖尿病、癫痫等)的发展过程划分为三个阶段，分别是健康阶段、临界阶段、疾病阶段，其中临界阶段持续的时间相较于其他阶段而言非常短暂。处于健康阶段的系统比较稳定，不易受到外界环境的扰动而改变，具有较强的恢复力和鲁棒性。而临界阶段的系统恢复力和鲁棒性均很弱，状态很不稳定，经过及时治疗，系统可以从临界阶段恢复到健康阶段，一旦系统进入了疾病阶段，则系统便处于了另一个较稳定状态，即疾病已经发展到重症阶段，再给予多大强度的治疗，系统也难以从疾病阶段恢复到健康阶段。

预测复杂疾病进展过程中的临界点具有非常重要的实际意义，但这项工作却是非常困难的。前人提出的单纯利用基因表达浓度的变化得到相关的不一致性指标的方法是静态的方法，不适用于探测复杂疾病的临界点。而利用传统的临界慢化法则需要病人提供大量的时序列样本数据，现实情况下难以在短时间内获取符合要求的疾病数据，所以利用生物网络对复杂疾病进行早期诊断是大势所趋。最近有学者提出了利用单个样本检测复杂疾病临界点的方法 [5] ，本文则基于利用单个样本与已知样本结合局部网络 [6] [7] [8] [9] [10] 以及熵的概念 [11] [12] [13] [14] 开发出了一个新的探测复杂疾病临界点的算法。本文首先利用希尔方程建立具有11个结点的复杂疾病生物仿真网络，并将我们的算法应用于仿真网络生成的数据，算法成功探测出了仿真数据的临界点。然后对NCBI生物网站上的小鼠肺光气实验(GSE2565)这个真实的高通量时序列数据应用本文的算法进行研究发现，该疾病即将恶化的预警信号出现的时间点与实验观察的结果吻合。

2. 方法

2.1. 单样本网络熵的推导

由于生理现象的复杂性，影响生理状态的因素非常多，目前尚无法用一个微分方程来描述。对于复杂疾病的临界点这个特殊时期，疾病的发展依赖于某几个主导因素，于是通过研究低维的动力系统微分方程来代换高维而复杂的疾病生理系统 [2] [3] [4] 。

对于一个复杂疾病临界点时的系统，可通过建模成为离散的动力系统：

$Z\left(t+1\right)=f\left(Z\left(t\right);P\right)$ . (1)

在系统(1)中的 $Z\left(t\right)=\left(z_1\left(t\right),z_2\left(t\right),\cdots ,z_n\left(t\right)\right)$ 代表固定t时刻下n维基因表达浓度的状态变量，参数向量 $P=\left(p_1,p_2,\cdots ,p_s\right)$ 代表生理系统的驱动因子，比如基因拷贝数变异、DNA甲基化等。

系统(1)需要满足一些条件：1) 存在不动点 $\bar{Z}\left(t\right)=\left({\bar{z}}_1\left(t\right),{\bar{z}}_2\left(t\right),\cdots ,{\bar{z}}_n\left(t\right)\right)$ 满足 $\bar{Z}=f\left(\bar{Z};P\right)$ ；2) 存在 $P_c$ 使得雅可比矩阵

${\frac{\partial f\left(Z;P_c\right)}{\partial Z}|}_{Z=\bar{Z}}$ , (2)

有一个或一对特征值的模等于1；3) 当 $P\ne P_c$ 时，系统(1)的特征根的模不总等于1。这三个条件意味着当P达到阈值 $P_c$ 时该系统在 $\bar{Z}$ 处经历了一个相变或一个余维数一的分岔。

在系统(1)的 $\bar{Z}$ 附近，当P到达 $P_c$ 前，按假设系统应该会留在一个稳定的不动点 $\bar{Z}$ 上，因此所有的特征值的模是在(0,1)之间的。在系统处于临界点时的参数值 $P_c$ 称为分岔点值或临界跃迁值。

现在考虑系统(1)的线性化近似方程，具体而言，通过引入新的变量 $Y\left(t\right)=\left(y_1\left(t\right),y_2\left(t\right),\cdots ,y_n\left(t\right)\right)$ 和

一个可逆的转换矩阵S，即

$Y\left(t\right)=S^{-1}\left(Z\left(t\right)-\bar{Z}\right)$ . (3)

结合系统(1)有

$Y\left(t\right)=\Lambda \left(P\right)Y\left(t-1\right)+\xi \left(t-1\right)$ ,(4)

(4)式中 $\xi =\left({\xi }_1,{\xi }_2,\cdots ,{\xi }_n\right)$ 是均值为零的微小高斯噪声， ${\xi }_i\left(i=1,2,\cdots ,n\right)$ 的标准偏差为 ${\sigma }_i\left(i=1,2,\cdots ,n\right)$ ， $\Lambda$ 是关于参数P的对角矩阵。

根据假设，对角化矩阵 $\Lambda =diag\left({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n\right)$ 的每个 ${\lambda }_i\left(i=1,2,\cdots ,n\right)$ 的模均在0和1之间。不失一般性，在 $\Lambda$ 的特征值中，假设模最大的一个为 ${\lambda }_1$ ，当参数转换 $P\to P_c$ 发生时， ${\lambda }_1$ 是第一个到达1的。特征值 ${\lambda }_1$ 表征了系统围绕不动点改变的速率， ${\lambda }_1$ 被称为占主导地位的特征值。系统远离临界点时的状态对应于 $\left|{\lambda }_1\right|<1$ ，而系统趋近临界点时的状态对应于 $\left|{\lambda }_1\right|\to 1$ 。设Y中与 ${\lambda }_1$ 相关的变量为 $y_1$ ，又根据方差与均值的性质得到 $y_i\left(i=1,2,\cdots ,n\right)$ 的方差为

$Var\left(y_i\right)=\frac{{\sigma }_i^2}{\left|1-{\lambda }_i^2\right|}$ , (5)

当主导特征值 $\left|{\lambda }_1\right|\to 1$ 时，从(5)式中，可以得到变量 $z_i=s_{i1}{y}_1+s_{i2}{y}_2+\cdots +s_{in}{y}_n+{\bar{z}}_i\left(s_{i1}\ne 0\right)$ 的标准偏差急剧增加的结论。因此，当系统处于远离临界点的状态时， $\left|{\lambda }_1\right|$ 远小于1并且没有大的变化， $Z\left(t\right)$ 和 $Z\left(t-1\right)$ 之间的分布也没有明显的变化。然而，当该系统处于趋近临界点时的状态，即 $\left|{\lambda }_1\right|$ 接近1时， $Z\left(t\right)$ 和 $Z\left(t-1\right)$ 这两者的分布便有差异。

对于以第i个结点为中心的局部网络，假设总共有m个一阶邻居 $i_1,i_2,\cdots ,i_m$ ，现令随机变量 $\mathcal{X}$ 代表局部网络的状态，局部网络中结点的值为0或1，则 $\mathcal{X}$ 总共取值个数为 $2^{m+1}$ ，设 $P_{u,v}\left(t\right)$ 代表t时刻第u个局部网络并且状态为v的概率，则第i个局部网络的熵为

$E_i^{\left(t\right)}\left(\mathcal{X}\right)=-{\displaystyle {\sum }_v{P}_{i,v}\left(t\right)\log P_{i,v}\left(t\right)}$ , (6)

整体网络熵为

$E^{\left(t\right)}\left(\mathcal{X}\right)=-{\displaystyle {\sum }_{u,v}{P}_{u,v}\left(t\right)\log P_{u,v}\left(t\right)}$ . (7)

系统远离临界点时，假设局部网络中某个结点 $i_l$ 为显著改变状态，即值取1的概率为

$P\left(x_{i_l}\left(t\right)=1\right)=P\left(\left|\Delta Z_{i_l}\left(t\right)\right|>d_i\right)=P_1$ , (8)

其中 $\Delta Z_{i_l}\left(t\right)=Z_{i_l}^c\left(t\right)-Z_{i_l}^n\left(t\right)$ ， $Z_{i_l}^c\left(t\right)$ 为t时刻单样本第 $i_l$ 个基因的表达浓度值， $Z_{i_l}^n\left(t\right)$ 为t时刻已知样本第 $i_l$ 个基因的表达浓度值，所以结点 $i_l$ 为不显著改变状态即值取0的概率为

$P\left(x_{i_l}\left(t\right)=0\right)=P\left(\left|\Delta Z_{i_l}\left(t\right)\right|\le d_i\right)=P_0=1-P_1$ , (9)

因为正常状态下基因表达稳定，所以 $P_1\to 0$ ，即

$P\left(x_{i_l}\left(t\right)=1\right)\to 0$ . (10)

令第i个局部网络的状态 $A_1=\left\{x_i=0,x_{i_1}=0,\cdots ,x_{i_m}=0\right\}$ ，则有

$P\left(\mathcal{X}=A_1\right)\to 1$ , (11)

而 $P\left(\left\{\cdots ,x_{i_l}=1,\cdots \right\}\right)\le P\left(x_{i_l}\left(t\right)=1\right)\to 0$ ，所以局部网络i处于 $A_1$ 以外的其他状态的概率是趋于0的，

故在正常状态有

$E^{\left(t\right)}\left(\mathcal{X}\right)=-{\displaystyle {\sum }_{u,v}{P}_{u,v}\left(t\right)\log P_{u,v}\left(t\right)}\to 0$ . (12)

而系统趋近于临界点时，对于第 $i_l$ 个结点，有 $sd\left(Z_{i_l}\right)\to \infty$ ，其中 $Z_{i_l}$ 为第 $i_l$ 个基因的表达浓度值，而 $\Delta Z_{i_l}\left(t\right)=Z_{i_l}^c\left(t\right)-Z_{i_l}^n\left(t\right)$ ，故有

$sd\left(\Delta Z_{i_l}\right)=sd\left(Z_{i_l}\right)$ , (13)

也即说明 $\Delta Z_{i_l}$ 的离散程度非常高，对于第i个局部网络，有

$E_i^{\left(t\right)}\left(\mathcal{X}\right)=-{\displaystyle {\sum }_{v=1}^{2^{m+1}}{P}_{i,v}\left(t\right)\log P_{i,v}\left(t\right)}$ , (14)

和

${\displaystyle {\sum }_{v=1}^{2^{m+1}}{P}_{i,v}\left(t\right)}=1$ . (15)

针对(14)和(15)二式引入拉格朗日乘子 $\lambda$ ，得到

$G_i^{\left(t\right)}\left(\mathcal{X}\right)=-{\displaystyle {\sum }_{v=1}^{2^{m+1}}{P}_{i,v}\left(t\right)\log P_{i,v}\left(t\right)}+\lambda \left({\displaystyle {\sum }_{v=1}^{2^{m+1}}{P}_{i,v}\left(t\right)}-1\right)$ , (16)

令 $G_i^{\left(t\right)}\left(\mathcal{X}\right)$ 对每个 $P_{i,v}$ 及 $\lambda$ 求偏导并令其等于0，可得

, (17)

和

${\displaystyle {\sum }_{v=1}^{2^{m+1}}{P}_{i,v}\left(t\right)}=1$ . (18)

所以 $\lambda -1=\log P_{i,v}\left(t\right)\Rightarrow P_{i,v}\left(t\right)={\text{e}}^{\lambda -1}=\frac{1}{2^{m+1}}$ ，则有

$\lambda =\log \left(\frac{1}{2^{m+1}}\right)+1$ , (19)

故当 $\Delta Z_{i_l}$ 是均匀分布的时候，单样本网络熵能取到最大值，此时有

$E_i^{\left(t\right)}\left(\mathcal{X}\right)=\log \left(2^{m+1}\right)=m+1$ . (20)

2.2. 探测临界点

1) 针对某种复杂疾病，收集该疾病已有的多个样本，并将该疾病对应的全基因组映射到STRING (STRING数据库是一个搜寻已知蛋白质之间和预测蛋白质之间相互作用的系统)的蛋白质–蛋白质相互作用网络。

2) 针对某个未知状态的样本，将该样本分别与每个已知样本建立状态网络。

3) 局部化全基因组状态网络。

4) 计算每个局部化状态网络的单样本网络熵。

5) 将所有局部状态网络的单样本网络熵排序并取前10%的平均值作为差异评价指标。

3. 主要结果

3.1. 仿真网络

本文利用11个结点网络生成疾病的仿真数据，然后利用新开发的算法探测疾病发展过程的临界点。这种类型的基因调控网络常常应用于研究基因转录、翻译等影响基因调控的活动。利用MM (Michaelis-Menten)方程并加入基因自身的降解率，可以得到11个结点的仿真网络，对应的微分方程组为：

(21)

其中P为数值控制参数， ${\zeta }_i\left(t\right)\left(i=1,2,\cdots ,11\right)$ 是均值为0的高斯噪声， $z_i\left(i=1,2,\cdots ,11\right)$ 为第i个mRNA的浓度，11个mRNA的自降解率分别为 $\left(\frac{8+3\ast \left|P\right|}{20},\frac{1}{2},\frac{3}{5},\frac{8+3\ast \left|P\right|}{20},\frac{6}{5},\frac{6}{5},\frac{6}{5},\frac{8}{5},\frac{9}{5},2,\frac{11}{5}\right)$ 。微分方程组(21)完全

体现了仿真网络之中结点间的调控关系，并且有一个稳定点：

$\bar{Z}=\left({\bar{z}}_1,{\bar{z}}_2,{\bar{z}}_3,{\bar{z}}_4,{\bar{z}}_5,{\bar{z}}_6,{\bar{z}}_7,{\bar{z}}_8,{\bar{z}}_9,{\bar{z}}_{10},{\bar{z}}_{11}\right)=\left(0,0,0,0,0,0,0,0,0,0,0\right)$ . (22)

将微分方程组(21)利用一个小 $\Delta t$ 进行差分化，然后用 $J={\frac{\partial f\left(Z\left(k\right);P\right)}{\partial Z}|}_{Z=\bar{Z}}$ 代表差分方程组的雅可比

矩阵，有

$J={\text{e}}^{\Delta t\ast A}$ , (23)

其中A为微分方程组(21)线性主部的系数矩阵，A的11个特征值分别为：

$\left(-\frac{3}{10}\ast \left|P\right|,-\frac{1}{2},-\frac{3}{5},-\frac{4}{5},-1,-\frac{6}{5},-\frac{7}{5},-\frac{8}{5},-\frac{9}{5},-\frac{10}{5},-\frac{11}{5}\right)$ , (24)

所以当 $P\in \left(0,1\right]$ 时平衡点 $\bar{Z}$ 是稳定的，当 $P\to P_c=0$ 时，系统便失去了稳定性，并且系统状态会经

历一个临界转变。

现在将雅可比矩阵J进行对角化，先引入可逆矩阵S：

$S=\left[\begin{array}{ccccccccccc}1& 0& 0& -1& 0& 0& 0& 0& 0& 0& 0\\ 1& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ -1& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 1& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 1& 0& -1& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& -1& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 1& 0& -1& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 1\end{array}\right]$ , (25)

使得 $S^{-1}JS=\Lambda$ ，对于 $\Lambda$ 有：

$\Lambda =diag\left({0.7408}^{\left|P\right|},0.6065,0.5488,0.4493,0.3679,0.3012,0.2466,0.2019,0.1653,0.14,0.11\right)$ ,(26)

即 $\Lambda$ 是个对角矩阵。现假设 $\Lambda$ 的特征值 ${\lambda }_i\left(i=1,2,\cdots ,11\right)$ 对应的特征向量分别为：

$\left(y_1,0,0,0,0,0,0,0,0,0,0\right),\left(0,y_2,0,0,0,0,0,0,0,0,0\right),\cdots ,\left(0,0,0,0,0,0,0,0,0,0,y_{11}\right)$ .(27)

注意到变量 $y_1$ 与最大的特征值 ${0.7408}^{\left|P\right|}$ 相关并且最接近1，因此，当P趋于0时，主特征值 ${0.7408}^{\left|P\right|}$ 趋于1，系统发生临界转变。又根据：

$Y\left(k\right)=S^{-1}\left(Z\left(k\right)-\bar{Z}\right)$ , (28)

可知 $\left(z_1,z_2,\cdots ,z_{11}\right)$ 中的四个变量 $z_1,z_2,z_3$ 和 $z_4$ 是直接与主特征值对应的 $y_1$ 相关的，根据2.1节的理论推导可知， $\left\{z_1,z_2,z_3,z_4\right\}$ 构成了11个结点仿真网络系统的主要部分，即相当于真实数据中的动态网络生

物标志物。

基于本文新开发的算法，收集11个仿真结点的时序列数据，并将算法应用于数值模拟数据得到的结果如下图1：

图1. 算法应用于仿真网络数据结果曲线

3.2. 真实数据

实验数据是从NCBI的GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)下载得到，数据编号为GSE2565。氯化羰(光气)是一种有毒的工业化合物，人或小鼠暴露于光气之中会产生肺水肿和不可逆急性肺损伤，令CD-1雄性小鼠全身暴露于空气(控制组)和光气(实验组)之中，并分别在暴露0.5、1、4、8、12、24、48、72小时后提取小鼠肺部RNA以获得基因表达浓度，实验中观测到在GSH还原转化进入氧化型GSSG的过程起非常重要作用的GPX-2在暴露4到12小时的时候基因表达显著增加，表现出显著性 [15] 。根据小鼠肺光气实验数据，利用本文基于单样本开发的新算法，得到单样本网络熵在第4小时显著增加并达到峰值，即探测到的临界点为第4个时间点，得到的结果与实验观测吻合，真实数据的结果如下图2：

图2. 算法应用于数据GSE2565结果曲线

4. 讨论

本文思考传统方法的局限性，创新性地将单样本地方法与熵进行结合，提出了一种新的有效稳定的预警复杂疾病临界点的方法。基于小鼠肺光气实验数据，本文开发的算法成功探测出了第四个时间点为临界点。复杂疾病比如前列腺癌、肝癌、肿瘤等等对人类健康造成巨大的威胁，而能够成功探测出复杂疾病即将恶化的信号，便能帮助人类及时预防甚至治疗复杂疾病，本文开发出的新算法是传统方法的一个突破。

本文仍然有不足之处，比如为了使模型简单，只讨论了线性化系统系数矩阵的特征值全不相同的这一种情况，并没有考虑有重根和共轭复根时算法的效果，这是以后需要改进的一个主要方向。

基金项目

广州市科技计划珠江科技新星专项项目资助(No. 201610010029)。

文章引用

李明亮,刘 锐. 单样本网络熵与某些复杂疾病恶性突变的预警
Single-Sample Network Entropy and Early Warning of Sudden Deterioration for Some Complex Diseases[J]. 应用数学进展, 2019, 08(01): 152-159. https://doi.org/10.12677/AAM.2019.81017

参考文献