考虑增殖效应的细胞分化能量景观分解
发布时间：
2026-09-28 21:02
考虑增殖效应的细胞分化能量景观分解
史继凡 :h1,2，相良和幸 ，李铁军 ，陈洛南 :h4,5,6,7,*
摘要
基因间的复杂相互作用决定了细胞的发育与分化。我们建立了考虑增殖效应的细胞分化景观理论，将发育过程建模为带有生死项的随机动力系统。研究发现，两种不同的能量景观(记为 和 )共同促成了非平衡稳态分化的形成。势能 是导致稳态分布的能量景观，其亚稳态代表细胞类型；而 则指示从多能细胞到分化细胞的分化方向。这种细胞分化的解释不同于以往未考虑增殖效应的景观理论。我们提出了构建景观 和 的可行数值方法及平均场近似。将其成功应用于典型生物模型，验证了能量景观分解的有效性，并揭示了所研究过程的生物学见解。
关键词:能量景观，细胞分化，随机系统
引言
细胞分化使生物世界丰富多彩。对基因调控网络(GRN)的动力学进行建模和理解，是探究细胞分化潜在机制的关键[1]。自Waddington提出具有开创性的表观遗传景观隐喻[2,3]以来，分化就被直观地描述为一个沿曲面滚动的小球，即能量景观。
构建生物系统能量景观的方法主要有两种。一种是基于数据的方法，该方法试图从实验数据，尤其是单细胞RNA测序(scRNA-seq)数据中识别细胞簇/细胞类型、分化轨迹、拟时间和细胞多能性[4-8]。目前已从图 、熵 和动力学系统理论中提出了多种算法[13-16]。在这些基于数据的方法中，种群平衡分析(PBA)[13-15]和分化动力学景观(LDD)[16]尤其考虑了细胞增殖率和死亡率，它们将细胞分化的假设从平衡态推广到非平衡稳态过程。另一种方法是基于模型的方法，该方法试图构建
© 2022作者。牛津大学出版社代表中国科学出版传媒股份有限公司出版。本文遵循CC BY 4.0许可协议(https://creativecommons.org/licenses/by/4.0/)，在适当引用的前提下，允许在任何媒介中自由使用、分发和复制。
知识共享署名许可协议(https://creativecommons.org/licenses/by/4.0/)允许在适当引用原作的情况下，在任何媒介中自由重复使用、分发和复制。 基因调控网络的动力学方程，分析系统的动力学行为，然后通过数值模拟构建能量景观。Wang等人提出了一种实用框架，用于构建不考虑细胞生死效应的生物系统能量景观[17-19]。他们将景观定义为 ，其中 是基因表达向量， 是系统的稳态概率密度函数(PDF)， 与内在小噪声的振幅有关。通过改变参数值来控制生物过程，景观模式会发生变化，从而直观地描述细胞类型的位置和细胞簇之间的转移概率。该框架已广泛应用于芽殖酵母细胞周期建模[17]、人类干细胞命运研究[19]以及秀丽隐杆线虫衰老研究[20]。然而，在Wang的景观中，由于未考虑细胞增殖效应，分化方向通常并非由内在因素决定，而是通过手动设置参数来控制。除上述主题外，还有许多其他与能量景观理论相关的有趣研究[21-27]。本研究中，我们采用基于模型的方-
方法，而是考虑分化过程中细胞的出生率和死亡率(BDR)，这一思路受PBA和LDD的启发。我们证明，存在两种重要的能量景观可用于描述分化动力学，分别记为 和 ，且可通过该模型计算。 中的亚稳态代表细胞类型， 的值与多能性相关。 的负梯度显示分化方向。考虑BDR对于构建多能性景观 至关重要。在此，我们阐述了考虑增殖效应的细胞分化能量景观分解(ELD)理论。还提出了用于构建能量景观的数值算法，以及高维情况下的平均场近似(MFA)。我们用三个例子展示ELD的应用，其中 和 直观地解释了细胞分化过程。
大学，上海200433，中国； 东京大学高等研究院神经智能国际研究中心，东京大学，东京113-0033，日本；3数学及其应用重点实验室(LMAM)与数学科学学院，北京大学，北京100871，中国； 系统生物学重点实验室，中国科学院分子细胞科学卓越创新中心，上海生物化学与细胞生物学研究所，中国科学院，上海200031，中国；5浙江省系统健康科学重点实验室，中国科学院大学杭州高等研究院，中国科学院，杭州310024，中国； 上海科技大学生命科学与技术学院，上海201210，中国；以及 广东省智能科学与技术研究院，珠海
，珠海
519031，中国
*通讯作者。电子邮箱:kaihara@g.ecc.utokyo.ac.jp；tieli@pku.edu.cn；Inchen@sibs.ac.cn
2022年1月6日收到；2022年5月22日修订；2022年5月25日接受
结果
在群体水平上建模细胞分化
受Weinreb等人[13]和Briggs等人[15]的启发，我们通过带有生死项的随机动力系统对细胞发育进行建模。我们用 表示基因表达水平。在细胞分化过程中，细胞群体用概率密度函数(PDF) 量化，其演化遵循著名的福克-普朗克方程(FPE)的广义形式，即
在 处具有初始概率密度函数(PDF) ，其中 表示基因间的生物相互作用， 是代表噪声幅度的参数， 表示 处细胞的净生死率(BDR)。在实际应用中，根据所考虑的基因调控网络(GRN)， 通常用激活或抑制的希尔函数建模。此处 表示细胞增殖， 表示细胞死亡。当 简化为不包含任何细胞增殖效应的标准福克-普朗克方程(FPE)时。为确保(1)存在非爆炸且非退化的稳态概率密度函数(PDF) ，我们要求以下条件作为 的约束:
需要存在非平凡的 来表征稳态下细胞生死的群体平衡。
单细胞水平的细胞分化建模
(1)中的种群水平动态也可视为单细胞水平的概率解释。我们考虑一个基因表达为 的细胞 ，其从 出发，于 时刻开始，提出伊藤意义下的加权随机动力学模型如下:
其中 的分布遵循 ，是具有独立分量的标准布朗运动，且 代表细胞 的时变权重。(1) 与 (3) 之间的关联可表示为
其中 为狄拉克δ函数，期望是对所有可能的轨迹 求取(详见在线补充材料)。当 >0时，权重 在 时刻增大，对应细胞增殖；当 时，权重 在 时刻减小，对应细胞死亡。若 ，则权重 对每个细胞均为常数，系统简化为无任何细胞增殖的情况，该情况由Wang等人[17,18]讨论。
能量景观分解
基于上述细胞发育的动力学模型，我们采用基于模型的方法构建景观，即假设 和 为已知先验。
我们将 表示为具有已知 BDR 的稳态 PDF，将 表示为具有 的稳态 PDF，对应式 (1)。注意，符号 与初始分布 不同。随后，构建两个能量景观如下:
驱动系统或细胞趋向稳态分布，且
图 1. ELD 框架示意图。(a) 调控细胞分化动力学的 GRN。动力学流 可通过 GRN 建模，其向量场绘制于示意空间 (X, Y) 中。(b) 细胞可按自然出生率和死亡率 R(x) 增殖或死亡。增殖过程(对应 )以红色显示，细胞死亡过程(对应 )以蓝色显示(采用 jet 配色图)。(c) 本研究中的 ELD 理论。动力学 ( , ) 可通过两个势项(能量景观) 和 以及一个旋度部分项 表征。我们使用两个稳态 PDF 定义 和 ，记为 和 。(d) 动力学流 (箭头)与出生率和死亡率 (jet 配色背景)的组合。(e–g) 分别为细胞类型景观 、多能性景观 和旋度部分项 的示意图。(e) 中 的亚稳态代表不同细胞类型。(f) 中 的值表示多能性，其负梯度描绘分化方向。(f) 中的色条通过 沿 轴的平均值计算得出，显示直观的全局分化方向。(g) 中，景观 作为背景绘制，旋度部分项 由向量表示。
量化由 引起的势变化，即细胞增殖和死亡的影响。景观 中的亚稳态盆地代表细胞类型，每个盆地的深度表征其稳定性。 的值描绘多能性，其负梯度场描述分化方向。从干细胞状态到分化细胞状态， 逐渐降低。势 为细胞类型景观， 为多能性景观。除这两个势函数外，剩余项 定义为
即描述所考虑动力学非梯度性的旋度部分。 项满足
对应于 中定义的旋度通量 的无散条件。
综上，我们在ELD框架中提出了针对以 对为特征的分化动力学的能量分解，如下所示:
势函数 、 及旋度部分 有助于深入理解沿分化通路的细胞分化与亚型。图1展示了ELD流程示意图。
ELD框架与现有景观理论的联系如下:在细胞无生死的情况下，即 ，式(9)简化为Wang的景观分解[17-19]，其中包含 和 ；若系统为梯度型，当 时，简化为PBA和LDD 中讨论的景观分解，其中包含 。
在仅已知 的情况下，PBA和LDD通过求解方程为梯度系统定义了能量景观
(10)并满足适当的边界条件。然而，对于非梯度系统，我们证明式(10)并非总是成立，除非旋度部分 几乎处处满足 (详见在线补充材料)。就此而言，我们利用式(6)给出了更通用的 定义。关于ELD的更详细讨论可在在线补充材料中找到。
能量景观的数值构建
根据式 (5) 中的势 和式 (6) 中的势 的定义，一旦获得稳态分布 和 ，我们就可以计算能量景观。接下来，我们分别讨论低维情况和高维情况的数值方法。
求解低维情形下的FPE
对于低维情形(通常维度小于3)，我们采用有限差分法或有限元法对(1)进行数值求解，直至达到稳态。对于小参数 ，传统方法中网格的空间尺度 需设为 ，才能解析这类对流主导方程的动力学行为，这并不实用。我们采用流线扩散法[28]，即使在 的情况下也能避免数值振荡。流线扩散法的详细过程在二维情形下的在线补充材料中有介绍，该方法也应用于我们“结果”部分的数值算例中。
平均场近似
高维情形
对于高维情形，直接求解(1)不可行。因此，我们从单细胞动力学(3)出发，提出一种MFA方法，将计算复杂度从指数级降低到多项式级。总体而言，在MFA中，势函数 和 由高斯混合模型近似，即
其中 为确定性动力学中通过稳定状态计数得到的组分数量， 是第 个高斯组分，均值为 ，协方差为 ， 为混合权重。这些近似由(3)关于小参数 的渐近分析得出。
短时渐近分析。我们将 记为(3)中 动力学的第 个吸引子。在 时间尺度(或 中的跃迁时间尺度)下，不存在状态跳变，我们推导出近似式 ，其中当初始 时 。均值和协方差满足以下方程
其中 和 是 在 附近泰勒展开至 项所得函数，详见在线补充材料。参数 为式(12)稳态。当 近似下 (即 )时，可得经典MFA [17,29]。
长时间渐近性。混合权重 的确定旨在考虑更长时间尺度下的 basin 跃迁，如 所示。在该机制下， 的扩散动力学被升尺度为连续时间马尔可夫链，其中 Arrhenius 型跃迁速率取决于相应吸引 basin 之间的能垒 [30,31]。我们假设升尺度后的跃迁速率矩阵为 ，并将 定义为平均 BDR。权重 在有无生死项 情况下的演化展现了该渐近性
(13)
其中 为 的混合权重， 为 的混合权重。由于难以获得精确的 ，我们利用该方程对稳态 进行蒙特卡洛近似
其中轨迹在 处的有限域上以均匀初始分布进行模拟，直至合适的有限终止时间 。这种选择实际上给出了与 中提出的无 项的 MFA 相同的结果。在遍历性假设下， 状态 的稳态混合权重 推导如下
其中 为状态 的退出速率。在约束方程 (2) 及 的附加假设下，还可得到 的近似，进而得到 的近似(详见在线补充材料)。
利用平均场近似 和 ，由 (5) 和 (6) 分别估计 landscape 与 。
需要说明几点详细说明。(i) 对于(1)也存在一种MFA，但对于多势系统，当分别近似 和 时，针对(1)的MFA可能具有不同数量的分量和不同的参数。这可能会导致出现奇怪的 。因此，针对(3)的MFA的优点之一在于，使用(12)计算的分量数 和 与BDR 无关，仅在(14)和(15)中估计的混合权重有所不同。(ii) 针对(3)的MFA不适用于单稳态系统，因为在唯一吸引域内轨迹的混合权重应相同。对于单分量系统，我们使用针对(3)的MFA，在线补充材料中也对此进行了详细描述。对于实际问题，大多数系统是具有多个不同分量的多稳态系统。(iii) 在(12)中，MFA在函数 和 中被展开为 。在Wang的框架 (其中 )中，仅被展开为 以估计能量景观 。在线补充材料中，我们指出，为了获得 的恰当近似，有必要采用针对 的MFA。还使用了两个分析实例来验证这种必要性。
数值示例
在本节中，我们通过三个示例展示如何应用ELD框架来描述细胞分化。
示例1:二维漂移-扩散过程
在第一个示例中，我们采用二维漂移-扩散过程作为模拟细胞分化的基础示例 。在本示例中，我们将 设为具有以下形式的梯度场
图2 二维漂移-扩散过程的势能景观。(a) 势能景观 。三个亚稳态对应三种细胞类型。(b) 势能景观 。细胞从 值较高的状态分化至值较低的状态。
并设定速率函数 ，其中 。噪声参数为 。势能 存在三个稳定平衡点，即 和 ，它们分别对应三种细胞类型。根据 ，大多数细胞在 附近增殖，并在 和 附近死亡。因此，位于 附近的细胞类型 代表具有高多能性的干细胞，而位于 附近的细胞类型 和位于 附近的细胞类型 则代表分化状态。采用流线扩散法求解含有和不含有 的方程(1)，分别得到分布 和 。此外，构建了如图2所示的能量景观 和 。在 中可以清晰识别出三种细胞类型。出生率和死亡率 的作用决定了该过程的分化方向，而 的值则表征从高到低的干性。在线补充材料中的图S1还显示了原始势能 ，以及 和 的二维投影。
示例2:双基因命运决定系统
第二个示例模拟由两个基因间相互作用控制的二元细胞命运决定过程[18,32,33]。图3a显示了两个基因之间的调控关系。该模型采用
图3 (a) 命运决定系统中两个基因 和 的相互作用。观察到自激活和相互抑制现象。(b) 当增殖-死亡率的幅值设为 时的势能景观 。(c) 当 从30降至3时，势能景观 的变化。(d) 当 从30降至3时，势能景观 的变化。
非梯度Hill动力学，其表达式为
其中参数设置为 、 和 。通过自激活和相互抑制，基因 和 (如 GATA1 和 PU.1)在多能干细胞中共表达，当系统定向为两种不同谱系时，其中一个基因会逐渐占据主导地位。在 Wang 的框架中， 和 发生变化以控制分化 [18]。然而，根据当前细胞分化的观点，ELD 理论认为，即使这些参数固定，细胞命运决定也可根据 BDR 的变化进行调整。为构建图谱 和 ，我们设置噪声幅度 和生死速率 。速率幅度 从 30 变为 0。图 和 展示了当 减小时图谱的变化，图 显示了 的图谱，其中 (即图谱 )。我们通过计算得到以下结果:(i) 当 BDR 较高( 为 30)时，仅存在一种以 为特征的细胞类型，其中 和 共表达；(ii) 随着 BDR 幅度降低，该单一状态分裂为两种，即发生分化；(iii) 当 时，两种分离的细胞类型中，一种基因占据主导地位；(iv) 的值表征细胞的多能性， 指示分化方向。总体而言，除了改变基因间的相互作用强度外，BDR 项 也可能对解释细胞分化起作用，它控制着细胞命运决定。
示例3:T细胞分化
第三个例子用于描述高维空间中的T细胞分化[34]。四个基因 (TCF-1)、 (PU.1)、 (GATA3)和 (BLC11B)通过激活或抑制相互作用(见图4a)。动力学项 采用希尔函数建模，参数列于在线补充材料中。我们将BDR设为 。利用MFA方程(11)，构建的景观如图4b和c所示。势函数 对应带有生死项 的稳态景观。四个亚稳态代表T细胞发育的四个阶段(ETP/DN1、DN2a、DN2b和DN3)。图 中的势函数 通过其值显示细胞的多能性，通过其负梯度场显示分化方向。应用于MFA时，BDR也针对每个高斯分量取平均值，因此每个亚稳态/细胞类型内 的变化远小于 所示的变化。因此，如图4c所示，相邻细胞类型之间的 存在急剧变化。
图4 (a)T细胞分化过程中的基因调控网络。四个基因 和 分别为TCF-1、PU.1、GATA3和BLC11B。(b)景观 和(c)景观 。亚稳态代表T细胞阶段，即ETP/DN1、DN2a、DN2b和DN3。势函数V表示多能性和分化方向。
除了针对T细胞分化的四维示例，我们还对在线补充材料第IV(D)节和图S2所示的含10个变量的高维示例进行了ELD分析。这些结果支持以下结论:ELD和MFA可通过能量景观 研究亚稳态，通过能量景观 研究多能性。
结论与讨论
本文提出了ELD框架，用于描述具有增殖效应的细胞分化过程。两个能量景观 和 可解释分化过程中系统的动力学行为。势函数 代表对应细胞类型的吸引子， 则表征多能性与分化方向。BDR的考虑对构建 至关重要。引入额外的 后，ELD理论是对传统能量景观理论的推广，也是带有生死项的Waddington表观遗传景观的自然实现。ELD的数值构建，尤其是其在 高维情况下的平均场近似，已被介绍。模拟实例证明了ELD的适用性。然而，仍有若干问题需要进一步讨论和研究。
首先，ELD框架与其他用于构建细胞分化景观的模型不同。PBA和LDD是数据驱动的，基于scRNA-seq数据，难以处理非梯度系统。Wang的景观未考虑分化过程中细胞的生死。ELD是基于模型的，分析由细胞增殖和死亡产生的景观 ，以显示分化方向。
其次，我们模型中的BDR 是先验已知的。对于实际情况， 可通过实验估计(如PBA，[13])，或针对每种细胞类型进行近似(如LDD，[16])。若仅关注细胞类型对应的 ，即 对每个亚稳态为常数，则对应的 在单个状态内变化很小，而在两个相邻状态间变化剧烈。
第三，根据吸引域的大小(落入某一亚稳态的轨迹百分比)，(14)中估计的混合权重 仅为粗略近似。更合理的方法是获取不同状态间的转移率矩阵 [31,35-37]。然而，在有限的模拟时间内难以捕捉罕见的转移。快速且准确地估计混合权重仍是一个未解决的问题。
第四，若给定合适的BDR，也可为具有极限环的系统构建景观。根据[17]中的近似， 可用作(11)中的第 个分量，其中 为极限环的周期， 为通过模拟(12)得到的 处概率密度函数(PDF)的高斯平均场近似。为具有极限环行为的真实系统构建景观将是我们未来的工作。
最后，本文采用带有BDR项的FPE研究正常细胞分化过程的动力学行为。对于由正常细胞发展而来的癌细胞，多能性和增殖速率可能通过外部因素恢复。我们将癌症或肿瘤细胞的能量景观研究留待未来开展。
总体而言，能量景观是表征系统动力学行为的通用概念，本研究提出的ELD可帮助理解除纯反应外还涉及增殖与死亡的系统。该工作也可应用于各类生物系统的基于模型和基于数据的动力学分析[38-42]。
材料与方法
详细方法可参见在线补充材料。
补充数据
补充数据可在 NSR 在线获取。
基金
本研究获日本科学技术厅登月型研发项目(JPMJMS2021)、日本医疗研究开发机构(JP21dm0307009)、日本学术振兴会科学研究费补助金(JP20H05921)及东京大学AI与超越研究所资助；并获中国国家重点研发计划(2017YFA0505500)、中科院战略性先导科技专项(XDB38040400)、国家自然科学基金(11825102、12131020、31930022、12026608)、北京人工智能研究院(BAAI)、广东省科技创新战略专项资金(2021B0909050004、2021B0909060002)及鹏城实验室重大重点项目(PCL2021A12)支持。
作者贡献
J.S.、T.L.、L.C. 和 K.A. 设计研究；J.S. 和 T.L. 执行研究；J.S. 分析数据；J.S.、T.L.、L.C. 和 K.A. 撰写论文。
利益冲突声明:无。
参考文献
1. Davidson EH. The Regulatory Genome: Gene Regulatory Networks in Development and Evolution. New York: Academic Press, 2006.
2. Waddington CH. Principles of Development and Differentiation. New York: Macmillan, 1966.
3. Waddington CH. The Strategy of the Genes. London: Routledge, 2014.
4. Tang F, Barbacioru C and Wang Y et al. mRNA-Seq whole-transcriptome analysis of a single cell. Nat Methods 2009; 6: 377-82.
5. Saliba AE, Westermann AJ and Gorski SA et al. Single-cell RNA-seq: advances and future challenges. Nucleic Acids Res 2014; 42: 8845-60.
6. Shi J, Teschendorff AE and Chen W et al. Quantifying Waddington’s epigenetic landscape: a comparison of single-cell potency measures. Brief Bioinformatics 2020; 21: 248-61.
7. Cannoodt R, Saelens W and Saeys Y. Computational methods for trajectory inference from single-cell transcriptomics. Eur J Immuno/2016; 46: 2496-506.
8. Saelens W, Cannoodt R and Todorov H et al. A comparison of single-cell trajectory inference methods. Nat Biotechnol 2019; 37: 547-54.
9. Trapnell C, Cacchiarelli D and Grimsby J et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol 2014; 32: 381-6.
10. Haghverdi L, Buettner F and Theis FJ. Diffusion maps for high-dimensional single-cell analysis of differentiation data. Bioin-formatics 2015; 31: 2989-98.
11. Teschendorff AE and Enver T. Single-cell entropy for accurate estimation of differentiation potency from a cell’s transcriptome. Nat Commun 2017; 8: 15599.
12. Jin S, MacLean AL and Peng T et al. scEpath: energy landscape-based inference of transition probabilities and cellular trajectories from single-cell transcriptomic data. Bioinformatics 2018; 34: 2077-86.
13. Weinreb C, Wolock S and Tusi BK et al. Fundamental limits on dynamic inference from single-cell snapshots. Proc Natl Acad Sci USA 2018; 115: E2467-76.
14. Tusi BK, Wolock SL and Weinreb C et al. Population snapshots predict early haematopoietic and erythroid hierarchies. Nature 2018; 555: 54-60.
15. Briggs JA, Weinreb C and Wagner DE et al. The dynamics of gene expression in vertebrate embryogenesis at single-cell resolution. Science 2018; 360: eaar5780.
16. Shi J, Li T and Chen L et al. Quantifying pluripotency landscape of cell differentiation from scRNA-seq data by continuous birth-death process. PLoS Comput Bio/ 2019; 15: e1007488.
17. Wang J, Li C and Wang E. Potential and flux landscapes quantify the stability and robustness of budding yeast cell cycle network. Proc Natl Acad Sci USA 2010; 107: 8195-200.
18. Wang J, Zhang K and Xu L et al. Quantifying the Waddington landscape and biological paths for development and differentiation. Proc Natl Acad Sci USA 2011; 108: 8257-62.
19. Li C and Wang J. Quantifying cell fate decisions for differentiation and reprogramming of a human stem cell network: landscape and biological paths. PLoS Comput Biol 2013; 9: e1003165.
20. Zhao L and Wang J. Uncovering the mechanisms of Caenorhab-ditis elegans ageing from global quantification of the underlying landscape. J R Soc Interface 2016; 13: 20160421.
21. Zhou JX, Aliyu M and Aurell E et al. Quasi-potential landscape in complex multi-stable systems. J R Soc Interface 2012; 9: 3539- 53.
ELD能量景观,细胞分化,细胞增殖,死亡效应,无散度卷曲项,非梯度动力学,平均场
表观信息
2026-10-09
2026-10-08
2026-10-03
纤维蛋白降解产物通过肌动蛋白细胞骨架调控创伤性脑损伤中的内皮屏障功能
2026-10-01
E3 SUMO连接酶PIASγ是调控糖尿病相关肝细胞核因子-1α活性的新型相互作用伴侣
2026-10-01
艾维缔官网
艾德官网
B站IVDSHOW
抖音军哥聊表观
视频号艾维缔
小红书艾维缔
快手表观盒子
表观遗传学
联系我们
