Deciphering the topologicallandscape of glioma using a network theory framework
用网络理论框架解读胶质瘤的拓扑结构
https://www.nature.com/articles/s41598-024-77856-y
胶质母细胞瘤干细胞已被认为是胶质母细胞瘤复发和治疗抵抗的关键因素，为新型治疗提供了有前景的目标。然而，对胶质母细胞瘤干细胞在胶质母细胞瘤层级中的作用理解有限，引发了争议，并阻碍了研究向治疗的转化。尽管我们对基因调控网络的理解取得了显著进展，但这些网络的动态及其对胶质母细胞瘤的影响仍然难以捉摸。本研究采用系统理论视角，将实验知识整合到胶质母细胞瘤的核心内源性网络模型中，通过网络动力学计算阐明其能量景观。该模型识别出两个稳定状态，分别对应于星形胶质细胞样和少突胶质细胞样肿瘤细胞，并通过一个具有高干细胞特征的过渡状态连接，该过渡状态作为星形胶质细胞样和少突胶质细胞样状态之间的能量屏障之一，表明胶质母细胞瘤干细胞在体内的不稳定性。我们还获得了各种稳定状态，进一步支持了胶质母细胞瘤的多细胞起源，并揭示了一组可能引发肿瘤异质性和治疗抵抗的过渡状态。本研究提出，连接两种胶质母细胞瘤稳定状态的过渡状态是胶质母细胞瘤异质性和治疗抵抗的核心。我们的方法可能通过提供关于胶质母细胞瘤生物学复杂景观的新视角，为胶质母细胞瘤治疗的发展做出贡献。
**关键词**胶质母细胞瘤，胶质母细胞瘤干细胞，基因调控网络，网络动力学，内源性网络理论，能量景观
多年前，具有神经/胶质干细胞特征的胶质母细胞瘤细胞被识别为胶质母细胞瘤中的癌症干细胞或胶质母细胞瘤干细胞（GSCs）[1-3]。GSCs在胶质母细胞瘤的复发和对放疗及化疗的抵抗中起着关键作用[4,5]，因此被认为是胶质母细胞瘤的有前景的治疗靶点。提出了胶质母细胞瘤层级模型，描述了胶质母细胞瘤中的所有细胞形成类似Waddington表观遗传景观的结构，其中GSCs位于顶端[5,6]。然而，对GSCs的有限理解引发了争议[7,8]，并阻碍了研究向治疗的转化[5]。GSCs是否是胶质母细胞瘤的起源细胞仍不清楚[4-6]，体细胞胶质细胞和胶质母细胞瘤细胞之间的中间状态仍然隐藏在迷雾中，甚至GSCs和非干细胞肿瘤细胞（NSTCs）在胶质母细胞瘤层级中的观点也是相互排斥的[5]。
胶质母细胞瘤和正常胶质系统可能遵循相同的分化规则[9]，理解胶质母细胞瘤致癌背后的原理和胶质系统的发展过程，即基因调控网络及其决定的景观[10-13]，可能会启发GSCs在胶质母细胞瘤层级中的作用，从而有助于找到新疗法的见解[5,14,15]。
目前，数据驱动的研究正试图基于高通量数据获得胶质母细胞瘤和胶质系统的基因调控网络和景观[13,16-18]。值得注意的是，Hoang等人[19]从数据中驱动了视网膜再生的基因调控网络，胶质母细胞瘤分子模式在单细胞尺度上的分布已广泛通过时间截面方式获得[20,21]。通过数据重建的景观识别了更多的调控基因[22]。即便如此，这些基于数据的研究目前仍面临数据质量和数量不足的困难[14,23]。
这些数据驱动的研究是过去几十年技术革命带来的数据爆炸的延续。从20世纪末到现在，癌症研究领域的先驱们发现了大量与癌症相关的知识，这些知识被总结为一组层级标志[24]，成为了一个里程碑。然而，这些先驱之一断言，癌症研究目前处于虚假繁荣状态，迫切需要新的理论基础和数学工具[25]。在生命科学的各个领域，建立更好的理论和理论与实验的结合，而不仅仅依赖于无法产生真正知识的数据，也被强烈呼吁[26-28]。
这样的理论早已萌芽。Wright[29]和Waddington[30]分别在进化和发育中引入了适应性/表观遗传景观的隐喻概念。Kauffman[31]引入了布尔动力学来模拟基因调控过程，证明了使用动力学过程来建模基因调控的可行性，并启动了景观的定量化。几十年后，随着必要知识和数据的积累，获得了第一个基于真实数据的定量景观[32]。此后，提出了内源性网络理论（ENT）[33]，该理论认为癌症和其他复杂疾病是进化塑造的细胞-分子网络动力学产生的景观中的内源性状态，作为进化动力学在细胞-分子水平上的构成结构[34,35]，随后解决了癌症和发展中的几个问题[36]。Lord等人[37]提供了连接基因网络动力学理论和实验观察表型的实验证据。基于类似系统视角的研究也在早期脑发育[38,39]和造血系统发育[40,41]中实施。
在此，我们采用上述系统视角，基于实验证据构建胶质母细胞瘤的核心内源性网络。我们希望，通过网络动力学描绘的胶质母细胞瘤景观的主要特征能够为理解胶质母细胞瘤提供新的见解，并为开发新型治疗策略提供理论支持。
结果
胶质母细胞瘤内源性网络的构建
基因调控网络的存在和重要性已无需进一步争论，但如何获得一个动态可行且生物学上合理的网络仍存在争议[14,15,42,43]。作为第一步，我们选择了一种自下而上的基于因果知识的方法，而不是自上而下的数据驱动方法来构建网络。
胶质母细胞瘤的核心功能模块和信号通路已经得到了很好的总结[44,45]。我们选择了以下功能模块和信号通路来构建胶质母细胞瘤的核心内源性网络（图1A）：RTK通路[46-50]、NF-κB通路[51-54]、Ras通路[51-54]、Akt通路[55-58]、HIF通路[59-63]、p53通路[64-68]、细胞周期[69-73]、细胞衰老[74-78]、细胞凋亡[79-86]和胶质分化[87-94]。
该网络由上述10个功能模块和通路组成，包含25个节点和75条边，其中包括44条激活边和31条抑制边（图1A，补充表1）。对于每个节点，我们采用了一个结合激活和抑制效应的方程来构建一个25维方程组（图1B，补充表1）。
网络动力学结果在不同计算框架下具有鲁棒性
布尔动力学由Kauffman在20世纪60年代引入基因网络动力学[31]。作为一种离散方法，布尔动力学可以快速捕捉整体结构。在本研究中，我们使用等效于布尔代数的阈值函数进行计算（图1B，补充表1）。为上述构建的网络获得了一组由25个布尔方程组成的布尔方程组（补充表1）。通过10^7个随机初始向量，我们获得了62个吸引子，包括16个点吸引子和46个线性吸引子，扩展线性吸引子后总共有116个状态（图2A，补充表1）。通过将随机初始向量的数量增加到10^8，我们获得了相同数量和类型的吸引子（补充表1）。此外，通过穷举考虑所有初始值，我们获得了相同数量和类型的吸引子（补充表1）。这些结果表明，我们在布尔动力学下获得了所有吸引子。
虽然布尔动力学作为一种离散方法可以快速产生初步结果，但它无法捕捉过渡状态和模拟小扰动的后果[95]，因此无法确定吸引子之间过渡的具体路径。生物学实验通常涉及时间切片，捕捉死亡或检测时刻的特征，因此不能保证系统已达到平衡。即使系统已达到宏观平衡，微观水平上也可能发生噪声驱动的状态转换。过渡状态对于确定细胞命运如何转换是不可或缺的[96]。体外培养的干细胞也表现出不稳定平衡状态的特征，这些状态在小扰动下容易分化。因此，不包括过渡状态的计算是不完整的。为了解决这个问题，我们需要使用连续计算框架。
朗之万方程或随机微分方程（SDE）自然对应于生物过程[34]。然而，求解一个n维非线性SDE（在本工作中n=25）需要过多的计算能力，更不用说目前还没有关于本工作中随机项的相关知识。此外，应用数学和物理学中常用的两种随机积分，Ito和Stratonovich，在SDE和常微分方程（ODE）之间的极值点上产生不一致的结果[97]。然而，基于先前开发的数学工具，SDE和ODE的极值点是相等的[97]。因此，求解高维非线性SDE可以简化为寻找高维非线性ODE的特殊解，这些特殊解足以描绘系统能量景观的主要特征。
在本工作中，我们选择了Hill方程进行模型构建（图1B；补充表1）。这里，x表示核心节点的浓度/活性，归一化为[0, 1]。降解率在本模型中归一化为1，因此每次迭代的降解量为-x。h和k分别表示Hill系数和表观解离常数的倒数。由于方程已经归一化，为了确保S形曲线的S部分尽可能接近中心（即当x=0.5时，f(x)≈0.5），h和k需要满足关系k≈2^h。基于之前的经验，我们最初选择了h=3和k=10（图1B）。
我们获得了一组25个ODE（补充表1）。我们使用均匀分布的随机初始向量进行数值计算：生成一组25维随机向量，对应于网络中每个节点的浓度/活性，作为每次迭代的起点。同时应用牛顿法和欧拉法以避免偏差。
我们依次将随机初始向量的数量增加十倍，从10^4直到10^8。结果在10^5个随机初始向量下迅速收敛到20个稳定状态（欧拉法），并在10^7个随机初始向量下获得相同的20个稳定状态（牛顿法）（图2B，D；补充表1）。我们还类似地识别了过渡状态，结果在10^7个初始随机向量下收敛到97个过渡状态（图2C；补充表1）。基于这些结果，可以得出结论，我们在ODE下识别了绝大多数（如果不是全部）稳定和过渡状态。
神经胶质景观
在之前的工作中，发现了一种能量函数的通用构建方法[99]。随后，在各种非线性系统中构造性地获得了能量函数[100-102]，并在二维条件下从数学上证明了非线性系统中能量函数的存在[103]。因此，我们相信，对于本研究中可以视为一组运动方程的系统，存在一个能量函数及其相关的能量景观。对于图3A中展示的假设能量景观，运动方程中的固定点对应于其能量景观中的极值点（红色点“峰”为不稳定固定点/过渡状态，蓝色点“谷”为稳定固定点/稳定状态）（图3B）。这些稳定和不稳定固定点之间的拓扑连接足以描述其能量景观的主要特征（图3C），从而绕过了获得这种高维非线性系统的全局能量函数的任务[36,97]。
我们对之前获得的每个过渡状态施加小扰动，以模拟系统从过渡状态到稳定状态的演化，从而获得过渡状态与稳定状态之间的连接，这些连接描述了能量景观的主要特征。我们将结果转化为固定点的拓扑连接图，其中较大的节点代表稳定状态，较小的节点代表过渡状态，节点的形状代表凋亡模块的水平，节点的颜色代表分化类型，节点标签的颜色代表细胞周期的激活状态（图4）。需要强调的是，每个节点的位置并不代表能量。
根据分化类型排列节点，我们获得了七个固定点集群。其中一个固定点集群具有高水平的干细胞特征，仅由过渡状态组成，这些状态将在小扰动下演化为稳定状态，符合体外培养的干细胞的特征，并被归类为胶质干细胞。直接连接到这个过渡状态集群的是两个固定点集群，一个具有高水平的少突胶质细胞和干细胞特征，被归类为少突胶质细胞前体细胞（OPCs），另一个具有高水平的星形胶质细胞和干细胞特征，被归类为星形胶质细胞前体细胞（APCs）。一个固定点集群表现出分化模块中所有基因的低水平，被识别为非（大）胶质细胞。两个固定点集群具有最低水平的干细胞特征和较高水平的星形胶质细胞或少突胶质细胞分化特征，被识别为星形胶质细胞或少突胶质细胞。值得注意的是，有一个过渡状态集群表现出略高的星形胶质细胞相关特征和略低的干细胞和少突胶质细胞相关特征，能够在小扰动下过渡到多种不同的分化类型。我们将这个集群称为星形胶质细胞样不稳定祖细胞，可能对应于NG2+细胞（图4）。
排列后，我们发现获得的胶质细胞的潜在能量景观包含了已知的胶质分化层级[9,104,105]，表明到目前为止在所有提到的建模和计算下的整体结果在生物学上是可接受的。然后，我们基于这些结果进一步探讨了与胶质母细胞瘤相关的含义。
定位胶质母细胞瘤状态
为了定位胶质母细胞瘤状态，我们首先在模块和通路水平上验证了结果。对于单一因素模块或通路，在本工作中，HIF、Ras、细胞衰老和RTK，单一因素本身的值代表模块或通路的水平。对于既有抑制又有激活节点的模块或通路，假设核心网络中的每个节点以及抑制和激活部分同样重要，我们平均激活和抑制水平，并从激活水平中减去抑制水平，然后将结果归一化为[0, 1]，以获得多节点模块的模块水平结果。0表示完全抑制，1表示完全激活，0.5左右的值表示激活/抑制平衡。
稳定状态sta_2、sta_5和过渡状态tra_6在模块和通路水平上表现出RTK、AKT、RAS、细胞周期和HIF的高激活，以及NF-κB的中等激活。p53、细胞衰老和凋亡模块被抑制。这与现有胶质母细胞瘤知识完全一致[44]（图5A，B）。sta_2表现出Sox9、Notch、Pros和NFIA的高水平，表明低分化星形胶质细胞的特征；sta_5表现出Olig、Sox10、Notch和Pros的高水平，表明低分化少突胶质细胞的特征；tra_6表现出Notch和Pros的高水平，以及Olig、Sox10、Sox9和NFIA的中等水平，表明具有胶质转换特征的神经/胶质干细胞，但尚未分化为特定的胶质细胞类型[9,89,93,106-114]。基于这些结果，我们将sta_2、sta_5和tra_6视为胶质母细胞瘤候选状态，其中sta_2对应于星形胶质细胞样肿瘤细胞，sta_5对应于少突胶质细胞样肿瘤细胞，tra_6对应于具有高干细胞特征的胶质母细胞瘤细胞（图5A，B）。
为了进一步确认这三个候选状态与胶质母细胞瘤的对应关系，考虑到低通量方法（如PCR、Western blot、IHC等）仍然是分子生物学中结果验证的黄金标准，我们将低通量实验数据与这三个候选状态在分子水平上进行了比较。我们通过回顾患者样本或细胞系的文献，获得了与模型相关的基因的转录和蛋白质功能水平数据[115-131]。在将我们的计算结果和低通量数据离散化后，我们发现胶质母细胞瘤状态与低通量数据的一致性超过95%，进一步确认了这些候选状态：sta_2、sta_5和tra_6为胶质母细胞瘤状态（图5C，D）。
单个低通量实验很少覆盖模型中的所有相关基因。因此，我们进一步使用高通量数据集验证计算结果。通过GTEx-TCGA数据集验证，一致性率为76.5%（图5F）。为了避免批次效应引起的偏差，我们进一步使用了包含配对肿瘤旁和肿瘤组织的GSE151352数据集，一致性率为77.2%（图5F）。详细信息在补充结果中描述。
由于本工作中的计算结果代表与核心节点相关的基因的功能水平，多组学分析已经揭示了转录组和蛋白质组之间的弱相关性[132]，进一步考虑到组织样本中不可避免的细胞类型异质性[20,133]，因此我们无法期望结果与转录组之间达到完美一致。标准正态分布数据的中位一致性率为64%，一致性率大于77%的概率小于5.9%。高通量数据验证的一致性（76.5%，77.2%）是可以接受的。
使用scRNA-seq数据，可以观察到肿瘤内的细胞状态层级。作为胶质母细胞瘤的“平均场”模型，本研究中获得的稳定和过渡状态应分别对应于平衡系统中分布的中点和集中分布之间的空间。在将sta_2、sta_5和tra_6与Darmanis等人[21]的GBM scRNA-seq数据进行比较时，我们发现所有三个肿瘤状态都落在肿瘤细胞的分布内，其中sta_5更倾向于OPC和少突胶质细胞集群（图S1，A-E，K-M）。
总之，我们可以自信地断言，胶质母细胞瘤候选状态：sta_2、sta_5和tra_6对应于胶质母细胞瘤。
胶质母细胞瘤附近的景观
在前述结果中，我们已经描绘了胶质系统的能量景观并定位了胶质母细胞瘤状态。为了回答开头提到的问题，我们需要仔细观察胶质母细胞瘤附近的景观。通过选择与胶质母细胞瘤状态仅通过一个过渡状态连接的第一邻近稳定状态以及连接这些稳定状态的过渡状态（图6A，B），我们获得了胶质母细胞瘤附近的能量景观，共包含12个稳定状态和22个过渡状态（图6C，D）。
胶质母细胞瘤附近的稳定状态是胶质母细胞瘤的起源细胞
在胶质母细胞瘤附近的能量景观中，所有非胶质母细胞瘤稳定状态只需穿越一个过渡状态即可转化为胶质母细胞瘤状态（图6A，B），这意味着只有一个能量屏障阻止这些非胶质母细胞瘤相关细胞的癌变命运，因此这些状态被认为是全局能量景观中最有可能转化为胶质母细胞瘤状态的细胞，换句话说，与这些稳定状态对应的细胞是胶质母细胞瘤的起源细胞，除了全零的死亡状态sta_17（图6C）。
根据分化模块排序，它们可以分为高干细胞特征的星形胶质细胞或APCs（sta_1）、高干细胞特征的少突胶质细胞OPCs（sta_6）、分化的少突胶质细胞（sta_11、sta_13、sta_18和sta_19），以及无法根据本研究识别的非（大）胶质细胞状态，可能是室管膜细胞、神经元或间皮细胞（sta_14、sta_16、sta_20）（图6A，C，D）。
60%的这些稳定状态表现出激活的凋亡模块特征（sta_1、sta_6、sta_11、sta_13、sta_16和sta_20），表明进入凋亡程序的细胞也可能具有进入癌变进程的潜力（图6C，D）。
先前的研究表明，细胞衰老促进胶质母细胞瘤的进展[134,135]，这在我们的结果中也有所体现。60%的胶质母细胞瘤起源状态（sta_1、sta_6、sta_11、sta_14、sta_19和sta_20）表现出由CDKN2A基因代表的激活的细胞衰老模块，该基因翻译为p16INK4（图6C，D）。
可以诱导形成胶质母细胞瘤的神经胶质细胞，实验表明，神经干细胞（NSCs）、成熟的星形胶质细胞和OPCs都可能是胶质母细胞瘤的起源细胞[136,137]。同样，我们的结果显示，具有高干细胞特征并经历凋亡的星形胶质细胞是一种胶质母细胞瘤起源细胞，而分化的星形胶质细胞不是胶质母细胞瘤的起源细胞。具有高干细胞特征并经历凋亡的少突胶质细胞也是胶质母细胞瘤的起源细胞。分化的少突胶质细胞也可能是胶质母细胞瘤的起源细胞，以及无法根据本研究单独识别分化类型的非（大）胶质细胞。这种一致性进一步验证了我们的理论结果。
理论预测
胶质母细胞瘤干细胞状态是分隔星形胶质细胞样和少突胶质细胞样肿瘤细胞的能量屏障之一
我们获得的所有胶质母细胞瘤状态都具有高干细胞特征（图4；图5C，D），这已经在单细胞测序数据的分析中观察到，大多数低级别胶质瘤和高级别胶质母细胞瘤样本中的高增殖细胞表现出高干细胞特征[20]。
在这些胶质母细胞瘤状态中，存在一个过渡状态，tra_6。它具有激活的细胞周期特征，以及高水平与胶质干细胞相关的基因Notch和Pros，以及中等水平的与星形胶质细胞分化相关的基因NFIA和Sox9，以及与少突胶质细胞分化相关的基因Olig和Sox10，代表已经经历胶质转换但仍保持高干细胞特征的低分化胶质细胞[9,93,108]。tra_6是一个不稳定的固定点，一个过渡状态。在小的扰动下，它将并且只能演化为两个胶质母细胞瘤稳定状态，sta_2和sta_5，这与胶质母细胞瘤干细胞（GSCs）的功能特征知识一致[5]（图7A）。
进一步，景观表明，一个稳定的胶质母细胞瘤状态可以通过tra_6在显著的靶向干预或累积的小扰动下演化为另一个状态。当这种转变大量发生时，过渡状态tra_6可以被清晰观察到。这也在实验中被观察到，部分大鼠胶质母细胞瘤细胞系C6自发表现出干细胞样表型[138]，以及发育细胞类型“外放射状胶质细胞”在胶质母细胞瘤中被重新激活[139]。
一个稳定状态需要足够的时间在累积的小扰动下演化为另一个状态，并且已经跨越能量屏障的状态与未跨越的状态相比呈指数分布，这意味着跨越的细胞远少于未跨越的细胞。Suvà和Tirosh观察到，与组织学分级一致的分化细胞数量远远超过另一种分化细胞类型[20]，这再次与我们的理论结果一致。
基于我们的理论结果与实验观察的比较，整个胶质景观中只有tra_6符合胶质母细胞瘤干细胞的特征，因此我们将tra_6视为胶质母细胞瘤干细胞状态（图7B，C）。
在通过固定点的雅可比矩阵的实特征值评估稳定性后，我们对所有不稳定的固定点引入小扰动，记录它们自发演化为稳定固定点，以映射固定点之间的连接网络。由于它可以自发演化（类似于球在被释放后自然从峰顶滚下），这意味着相应的稳定状态比过渡状态具有更低的能量。能量屏障的定义是驱动反应自发进行所需的最小能量，这可以使用Kramers公式定量计算[140,141]。如图3A-C所示，在能量景观中，连接两个稳定状态的过渡状态可以被视为它们之间的能量屏障。我们发现了多个同时连接两个胶质母细胞瘤稳定状态的过渡状态（图6C，图7I），其中只有一个过渡状态，tra_6，仅连接两个胶质母细胞瘤稳定状态，而其他七个过渡状态连接多个稳定状态，可能充当超级过渡状态。这意味着在我们的结果中，有两种胶质母细胞瘤之间的转变存在八条潜在的“反应路径”。当然，确定超级过渡状态的标准包括在固定点的雅可比矩阵中具有两个或更多不稳定模式，而过渡状态只有一个不稳定模式；我们在本工作中没有强调这两者之间的区别。我们的经验表明，相邻的超级过渡状态通常具有更高的能量，表明tra_6在这一系列过渡状态中更有可能具有最低的能量。虽然真正的挑战在于定量计算能量，这需要实际的动力学参数。但由于从过渡状态到其连接的稳定状态的演化是自发进行的，我们可以定性地断言，这些过渡状态都是它们连接的稳定状态之间的能量屏障。
因此，胶质母细胞瘤干细胞状态是分隔星形胶质细胞样和少突胶质细胞样肿瘤细胞的能量屏障之一，这意味着以下推论和现象：
1. **胶质母细胞瘤干细胞（GSCs）在健康的野生型个体中极不稳定，因此在这些个体中很少被观察到；**
2. **GSCs可以在星形胶质细胞-少突胶质细胞转化过程中被观察到；**
3. **当反应是可逆的时，GSCs总是可以被观察到。**
作为推论1和2的证据，Bhaduri等人发现，一种在成熟个体中很少观察到的放射状胶质样细胞在胶质母细胞瘤中重新出现[139]；
作为推论2和3的证据，Suvà和Tirosh证明，GSCs在所有胶质母细胞瘤样本中与星形胶质细胞样和少突胶质细胞样胶质母细胞瘤细胞连续分布，而在高级别胶质母细胞瘤样本中GSCs的数量远高于低级别胶质母细胞瘤样本[20]；
作为推论3的证据，部分大鼠胶质母细胞瘤细胞系C6自发表现出干细胞样表型，并受到“反应条件”的影响[138]。
目前，我们自信地断言，在我们的理论模型中，GSC状态是分隔星形胶质细胞样和少突胶质细胞样肿瘤细胞的能量屏障之一。这一断言的实验证据仍然是间接和回顾性的，需要进一步的数据分析和进一步的投资来设计和进行实验以验证这一理论预测。
诱导异质性的过渡状态是治疗抵抗和胶质母细胞瘤复发的潜在机制
在胶质母细胞瘤附近的能量景观中，所有过渡状态（除tra_6外）连接胶质母细胞瘤起源状态和胶质母细胞瘤状态，作为潜在的致癌途径，同时也是体细胞胶质细胞和胶质母细胞瘤之间的中间状态（图6C，D）。这些过渡状态可以根据它们与胶质母细胞瘤状态的连接分为两类：
第一类过渡状态仅连接一个胶质母细胞瘤状态，例如仅连接sta_2和体细胞状态的过渡状态（tra_4、tra_13、tra_14、tra_21、tra_23、tra_38、tra_48）以及仅连接sta_5和体细胞状态的过渡状态（tra_7、tra_34、tra_37、tra_39、tra_45、tra_52、tra_54）（图7D，E，F）；
第二类过渡状态同时连接两种胶质母细胞瘤状态，包括tra_6、tra8、tra9、tra15、tra18、tra31、tra33和tra_51（图7G，H，I）。有趣的是，所有诱导异质性的过渡状态在分子水平上表现出高干细胞特征。
当从一种胶质母细胞瘤状态通过第一类过渡状态驱动时，系统只会转移到体细胞稳定状态，这可以被认为是潜在的治愈（图6C，左&右；图7D，E，F）。然而，如果被驱动到第二类过渡状态，系统有可能从一种胶质母细胞瘤稳定状态转移到另一种胶质母细胞瘤稳定状态以及多种体细胞状态，导致肿瘤异质性（图6C，中间，图7G，H，I）。在一个存在噪声的系统中，自发跨越屏障的概率总是非零的，这与屏障的高度负相关，并与噪声的幅度正相关[140,142]。高级别胶质母细胞瘤通常对应于更多的肿瘤细胞，自发跨越屏障事件的频率也随着肿瘤细胞数量的增加而增加。只要通往这些诱导异质性的过渡状态的路径存在，就很难达到真正的治愈，这些治疗行为可能会进一步导致胶质母细胞瘤临床级别的增加。
模型原理与局限性
对于一个随机系统，存在两种等价的描述：从分布角度可以使用偏微分方程（PDE）形式的Fokker-Planck方程，从粒子运动角度可以使用随机微分方程（SDE）或朗之万方程[143,144]。如果采用分布角度，限制将是数据的质量和数量。从这一角度出发，专注于胶质母细胞瘤景观的研究已经取得了一些有价值的结果[22,145]。
我们的方法采用粒子运动角度。我们的一位作者之前描述了一组通用的进化动力学方程[34]。相应的动力学方程类似于牛顿第二定律，是经典的朗之万方程。在之前的工作中，我们已经证明，任何动力学系统都存在能量函数或能量景观[99-101,103]，并且在A型积分下，SDE的固定点与ODE的固定点一致[97,146]。因此，求解SDE可以被求解ODE所取代，以捕捉能量景观的主要特征。使用ODE的前提假设是所有空间维度上的均匀分布，这对于本工作所要回答的问题是足够的。
我们在本工作中选择了自下而上的路径，基于实验知识构建动态方程以获得能量景观。由于没有使用输入数据，获得的能量景观不受数据质量和数量的直接影响，这从根本上将我们的方法与数据驱动方法区分开来。局限性在于动态结果依赖于网络构建，而网络构建依赖于实验知识：基因相互作用网络的细节和参数大多未知。我们的解决方案是基于构成理论[33]构建一个“近似”的网络模型并求解方程，然后将结果与实验知识进行比较、修正和迭代，以获得解释和预测[36]。
至于我们使用的构成方程：Hill方程，当Hill系数趋近于无穷大时，它等价于布尔方程。在我们的工作中，布尔动力学被用作一种快速获得动态结果的方法（即初步判断模型构建是否能获得合理的动态结果），并作为ODE中随机初始向量数量是否足够的并行标准。
如果在上述框架下使用偏微分方程（PDE）模型，它可能有助于解决各向异性问题，例如肿瘤血管生成或沿腹背轴和前后轴的发育过程。最近，本工作的一些作者也尝试使用相关方法基于PDE模型解决一些问题[102,147]。
应用景观回答生物学问题的方法最近引起了越来越多的关注。一些努力采用了自上而下的视角，试图直接从数据中推导出景观[13,14,18,22,145,148]；其他工作，包括我们自己的工作，探索了自下而上的视角，从基本原理出发计算景观[32,33,36,40]。如果现实由同一组规则支配，研究人员观察到相同的现象，即使不同的起点和视角暂时产生不同的场景，我们相信它们最终会在合作和相互增强中收敛于相同的现象。
该模型模拟了野生型生物体中的胶质母细胞瘤景观，表明在没有遗传突变的情况下，胶质母细胞瘤仍然可能发生。然而，在遗传突变下，可能会发生剧烈的景观变化，可能导致肿瘤吸引子的扩大，甚至减少肿瘤和非肿瘤状态之间的过渡状态数量。在未来的研究中，我们将探索胶质母细胞瘤中常见的遗传突变如何改变景观。
进一步设计的实验以验证本研究
本研究中的理论结果验证都是回顾性的。以下前瞻性实验，如调查GSCs在纯化的胶质母细胞瘤细胞系中的存在，或植入来源于不同纯化胶质母细胞瘤细胞的肿瘤，可以用来验证我们的预测，即胶质母细胞瘤干细胞作为连接胶质母细胞瘤稳定状态的过渡状态，以及异质性和治疗抵抗由诱导异质性的过渡状态介导。
理论结果表明，肿瘤间异质性可能由一组过渡状态介导，并且在噪声环境中，系统有可能自发跨越势垒。如果在纯化的次级细胞系中自发发生异质性，并且此过程中的分子模式与理论结果一致，则可以作为理论发现的实验验证。
对于胶质母细胞瘤细胞系或原代胶质母细胞瘤细胞，我们建议引入与分化相关因子连接的敲入荧光蛋白标记，如本工作中所讨论的。我们计划使用流式细胞术分离星形胶质细胞样和少突胶质细胞样肿瘤细胞，在体外培养和传代过程中监测荧光信号。我们预计会观察到随着培养和传代过程中固有背景噪声的增加，荧光信号从均一到异质的转变。此外，可以使用微流控芯片进行长期体外观察，使我们能够如Lord等人[37]所述监测细胞命运切换。在体内实验方面，可以使用纯化细胞进行肿瘤植入。我们期望在植入早期观察到大多数低级别胶质瘤（LGG）表型，低剂量凋亡诱导治疗后异质性增加，以及凋亡诱导+分化传播组中异质性较低。
潜在的临床意义
胶质母细胞瘤的分化治疗已经被呼吁（Laug, Glasgow, and Deneen, 2018; Suvà and Tirosh, 2020）。在这里，我们的理论结果初步表明，连接两种胶质母细胞瘤稳定状态的过渡状态是胶质母细胞瘤异质性和对放化疗抵抗的原因。这些诱导异质性的过渡状态都表现出高干细胞特征和低增殖特征。因此，结合凋亡促进疗法和分化促进疗法可能是胶质母细胞瘤治疗的一个潜在策略。
针对我们模型中核心基因的分子治疗策略已经存在[149]。我们的方法可以作为开发个性化鸡尾酒疗法的体内导航工具，通过微调网络和参数。然而，我们必须承认，我们25节点的模型过于粗粒度，未来的工作应该进一步细化并扩大模型。
本研究中使用的模型是“平均场”模型，适用于“平均人”，并未考虑人群中存在的遗传多态性。网络结构和参数可以调整以进一步适应遗传多态性。该模型是可扩展的，意味着它可以“即插即用”，添加或移除功能模块和信号通路以适应进一步的探索。
结论
在本研究中，通过采用系统理论视角并整合基于实验的知识，我们构建了胶质母细胞瘤的内源性网络，并通过网络动力学计算描绘了胶质母细胞瘤系统的能量景观。通过将理论计算与多层次知识和数据进行比较，我们能够定位胶质母细胞瘤状态。
我们的结果揭示了两个稳定状态，分别对应于星形胶质细胞样和少突胶质细胞样肿瘤细胞，并通过一个对应于GSC的过渡状态相互连接。这表明GSCs在体内是不稳定的，不太可能是胶质母细胞瘤的细胞起源。胶质母细胞瘤状态附近的景观展示了多种分化类型的几个稳定状态，进一步从理论上确认了胶质母细胞瘤的多细胞起源。
在胶质母细胞瘤状态附近的景观中，一组同时连接两种胶质母细胞瘤稳定状态的过渡状态表现出非增殖和高干细胞特征。这些过渡状态，包括GSC状态，都是星形胶质细胞样和少突胶质细胞样肿瘤细胞之间各种细胞命运转换路径上的能量屏障。如果系统通过这些过渡状态驱动，这将导致肿瘤异质性增加。如果在治疗过程中不阻断这些诱导异质性的过渡状态，胶质母细胞瘤的真正治愈将难以实现。
基因调控网络的研究在最近几十年显著推进了我们的理解。然而，将这一系列静态知识整合到一个能够产生动态洞察的系统框架中仍然是一个挑战。在本研究中，我们试图通过整合实验知识来构建胶质母细胞瘤的核心内源性网络来弥合这一差距。通过引入基于生物学原理并在统计物理学中广泛应用的景观概念，我们寻求产生新的见解。我们的理论发现虽然有前景，但代表了一个初步步骤，需要进一步的实验验证。它们提供了一个理论模型与实验研究之间复杂相互作用的初步视角，这种关系以相互增强和持续对话为特征。随着我们继续前进，我们怀着谨慎的乐观态度，期望我们的方法能够为更深入地理解胶质母细胞瘤生物学做出贡献，并最终为开发更有效的治疗干预措施提供支持。
方法
网络动力学
布尔动力学
布尔动力学是一种有用且相对快速的工具，用于获取网络动力学结果的架构信息[95]。典型的布尔函数如图1所示，以Akt为例：每个变量具有离散值0或1，分别表示低/高水平；t表示变量的当前状态，t+1表示下一个状态；NOT、AND和OR是通常的逻辑运算符。在我们的计算中，使用等效的阈值函数在Matlab中进行。
**ODE**
我们使用两种不同的算法来计算稳定状态。
**算法1（欧拉法）**：我们通过生成一个随机起点（初始向量）来启动程序。利用动态系统的规则，我们进入一个迭代计算的循环，逐步使用时间步长乘法因子更新向量。迭代过程设计为继续进行，直到达到预定义的步数，确保向量收敛到稳定状态，这通过涉及预设阈值的特定收敛标准来确定。一旦稳定状态的候选者满足收敛标准，我们将其分类并记录为稳定状态。
**算法2（牛顿法）**：
我们采用数值方法来求解表示系统变量随时间变化率的非线性方程。同样，我们的过程从初始化一个随机向量作为迭代程序的起点开始。使用MATLAB中的fsolve函数，我们应用牛顿法来定位使函数为零的解向量，这表示一个固定点。为了确定该解的稳定性，我们评估了解点处雅可比矩阵的特征值，仅将所有特征值具有负实部的解视为稳定状态，而存在正特征值的解被视为过渡状态。此过程系统地重复进行，以编译一份全面的稳定状态列表。
算法的详细信息已在之前的出版物中描述[36]，本工作中使用的函数见补充表1。
**状态互连分析**
我们的算法旨在生成从每个不稳定状态到其相应稳定状态的轨迹。为此，我们在系统位于已识别的不稳定状态时，通过小的随机向量引入轻微扰动。这使我们能够追踪从不稳定状态到每个连接的稳定状态的路径。详细信息见补充表1。
**ODE结果的模型化**
在我们的研究中，方程中的变量‘x’表示网络中节点的浓度或生物活性。我们将计算结果转换为模块水平，以与已建立的生物学知识对齐。以下方程用于确定模块状态：
功能模块和信号通路中的节点被分类为激活性（xact）或抑制性（xinh），分化模块除外。最终模块状态由阈值确定：当模块状态低于0.33时，表示“关闭”状态，表示抑制；当模块状态超过0.66时，表示“开启”状态，表示激活；当模块状态介于0.33和0.66之间时，分配为“稳定”（STA）状态，反映模块内激活和抑制影响的平衡。在具有内部负反馈的模块中，以NF-κB通路为例：STA状态可能表明类似于NF-κB亚基与其抑制剂pIkB形成复合物的2:1比例的平衡；“开启”状态意味着未与pIkB复合的游离NF-κB亚基占优势，从而能够进入细胞核并发挥其功能；“关闭”状态表示抑制基因占主导地位，例如相对于NF-κB亚基，IkB水平较高。
高通量数据预处理
**TCGA-GTEx**
我们从GTEx皮质样本和TCGA GBMLGG数据集中获取了RNAseq数据（TPM）。从皮质、GBM和LGG集合中随机选择150个样本，我们对数据集进行了两步缩放。首先，我们计算每个样本的z分数以减轻批次效应，然后计算每个基因的z分数。超过一个标准偏差的值被限制在此阈值内，随后归一化到[0, 1]范围。1 SD上限定义了分布两侧约16.5%的数据为异常值，从而防止异常值影响数据分布。更高的阈值可能导致不同组的数据在归一化到[0, 1]区间后集中在0.5附近，这可能会显著影响结果。此外，1 SD阈值下的实验数据与模型结果具有相似的分布形式（图S1，F，G）。基因表达的连续分布由R中的ggplot2包绘制。
**GSE151352**
从GEO数据库中检索GSE151352的数据，转换为指数格式（e），并计算其标准误差。我们通过平均具有前200个最小标准误差的基因的标准误差（≈30%）来确定该数据集的有效信息阈值。原始数据集通过从肿瘤样本值中减去癌旁样本值来推导倍数变化。该数据集的标准误差约为30%，假设为系统误差，因此只有高于此阈值的信息才能与噪声区分开来。高表达阈值设置为ln(1.3)，将高于此阈值的倍数变化限制为ln(1.3)；低表达阈值设置为ln(1/1.3)，将低于此阈值的值调整为ln(1/1.3)。随后，数据归一化到[0, 1]范围。