乳腺癌EMT过程的动力学建模与分岔分析

Dynamical Modeling and Bifurcation Analysis of the Epithelial-Mesenchymal Transition Process in Breast Cancer

摘要: 乳腺癌是全球女性最常见的恶性肿瘤之一，其高死亡率主要源于肿瘤转移。上皮–间质转化(EMT)是驱动转移的关键步骤，其核心调控网络涉及miR-200、ZEB1和CDH1之间的复杂相互作用。本研究通过构建该网络的三维动力学模型并运用分岔理论，深入揭示了EMT的动态调控机制。分析表明，ZEB1对miR-200和CDH1的抑制强度可作为关键开关，通过诱导鞍结分岔直接影响细胞状态的稳定性。进一步的双参数分析揭示了更复杂的调控景观，ZEB1与miR-200的相互抑制强度、以及ZEB1对miR-200的抑制与CDH1对ZEB1的抑制这两组参数组合，均能引发导致表型突变的尖点分岔。这些研究结果表明，精准调控关键参数可有效维持上皮状态或抑制EMT进程，从而为逆转乳腺癌转移提供了重要的理论依据。

Abstract: Breast cancer is one of the most common malignant tumors in women worldwide, and its high mortality rate mainly stems from tumor metastasis. The epithelial-mesenchymal transition (EMT) is a key step driving metastasis, and its core regulatory network involves the complex interactions among miR-200, ZEB1, and CDH1. This study, by constructing a three-dimensional dynamic model of this network and applying bifurcation theory, deeply reveals the dynamic regulatory mechanism of EMT. Analysis indicates that the inhibition strengths of ZEB1 on miR-200 and CDH1 can act as critical switches, directly affecting the stability of cellular states by inducing saddle-node bifurcations. Further two-parameter analysis revealed a more complex regulatory landscape: the mutual inhibition strength between ZEB1 and miR-200, as well as the combination of ZEB1’s inhibition on miR-200 and CDH1’s inhibition on ZEB1, can both trigger cusp bifurcations that lead to phenotypic mutations. These research results indicate that precise regulation of the key parameters can effectively maintain the epithelial state or inhibit the EMT process, thereby providing an important theoretical basis for reversing breast cancer metastasis.

1. 引言

乳腺癌是女性最常见的恶性肿瘤之一，其发病率和死亡率居高不下[1]。尽管近年来在诊断和治疗方面取得了显著进展，但乳腺癌的转移和复发仍然是患者生存率低的主要原因。上皮–间质转化(Epithelial-Mesenchymal Transition, EMT)是乳腺癌转移过程中的关键机制之一，其本质是上皮细胞向间充质表型转变的动态过程，见图1。该过程的核心特征是E-cadherin等黏附分子下调导致细胞连接解离，并使细胞获得迁移与侵袭能力[2]，因此与肿瘤侵袭、转移形成及治疗耐药密切相关。该过程通过中间过渡状态逐步完成，表现为上皮标志物与间充质标志物的共表达[3] [4]。在乳腺癌中，EMT引发显著的细胞去分化，形态上由多边形转变为梭形、细胞极性消失并发生骨架重组，功能上则表现为细胞黏附减弱和运动能力增强[5] [6]。

Figure 1. EMT process

图1. EMT过程

值得注意的是，EMT过程在某些情况下是可逆的，这一逆过程被称为间充质–上皮转化(Mesenchymal-Epithelial Transition, MET)。研究表明，MET在转移灶处循环肿瘤细胞的定植中起着关键作用[7]。MET是通过上调E-cadherin等上皮标志物和下调N-cadherin等间充质标志物，使细胞恢复上皮表型[8]。EMT由一个复杂的基因网络控制，这些基因在转录、转录后和翻译后水平上的复杂相互作用精确调控了这一转变过程及其动力学行为。广泛认可的EMT相关基因的例子包括编码ZEB和SNAIL家族等转录因子的基因，以及波形蛋白等效应蛋白[9]。此外，许多microRNA被证明在EMT/MET中发挥关键作用[10]。microRNA (miRNA)是一类单链非编码RNA，主要在转录后水平调控基因表达，已发现的1200多种哺乳动物miRNA可能靶向多达三分之一的蛋白质编码基因，参与发育、分化、代谢、信号转导等多种生物过程[11] [12]。特别是在癌症领域，miRNA表达失调与肿瘤进展、转移及化疗耐药密切相关，其靶标涉及肿瘤抑制基因或癌基因[13] [14]。TGFβ诱导的EMT常伴随CDH1 (E-cadherin)启动子的表观遗传重塑和miRNA介导的E-钙粘蛋白表达下调[15] [16]，其中，miR-200家族是调控EMT的重要因子，通过抑制ZEB1/2维持CDH1表达，阻止细胞转化为间充质状态[17]。该家族包括miR-200a/b/c、miR-429和miR-141 [18]。

研究表明，miR-200的表达升高可逆转乳腺癌中的EMT，具有潜在治疗价值[19]。其过表达可阻止EMT并诱导MET，是维持上皮表型的重要因子[20]。miR-200家族与EMT激活因子ZEB1/2之间形成双负反馈回路：ZEB1/2抑制miR-200簇的转录，而miR-200通过结合其3'UTR反向抑制ZEB1/2的表达[21] [22]。这一调控机制在乳腺癌转移中发挥关键作用，使细胞在上皮、部分EMT与完全EMT状态间动态切换，形成表型双稳态[23] [24]，其中部分EMT状态的细胞往往比完全EMT更具侵袭性[25]。已有研究进一步揭示ZEB1可结合miR-200启动子E盒位点抑制其表达[26]，并通过建立了miR200-ZEB1的动力学模型对EMT过程进行定量描述[27]。另有研究构建二维模型模拟乳腺癌细胞在干细胞样、基底和管腔状态间转换[28]。此外，miR-200还能直接抑制PD-L1，而PD-L1可间接下调CDH1表达[29] [30]。然而，已有的miR-200-ZEB1模型[27]往往忽略了CDH1的协同作用和miR-200-PD-L1-CDH1这条间接通路，缺乏对三因子高维耦合动力学的建模，因而难以全面揭示EMT背后的复杂调控机制。

基于此，本文通过构建ZEB1、miR-200、CDH1调控网络的动力学模型，揭示了乳腺癌EMT的多稳态转换机制。首先，基于Sahoo等人[29]提出的由miR-200、ZEB1和CDH1构成的调控网络，构建了一个描述乳腺癌EMT过程的微分方程系统，该系统重点刻画了miR-200、ZEB1和CDH1之间的动态相互作用，以及miR-200通过PD-L1这条路径间接作用于CDH1。其次，对所建立的微分方程系统进行了深入的平衡点稳定性分析，系统推导了系统产生鞍–结点分岔的条件，并提供了完整的理论证明过程。最后，在理论分析的基础上，通过数值模拟方法开展了系统的单参数和双参数分岔分析，绘制了相应的分岔图，并深入探讨了如何通过调控关键参数来控制乳腺癌转移进程，为临床治疗提供了潜在的理论依据。

2. 模型建立

大量实验研究表明，EMT过程受到由SLUG、ZEB1、miR-200、CDH1和PD-L1组成的调控网络的精细调控[29]。见图2，SLUG和ZEB1作为转录因子协同促进上皮细胞向间质细胞转化，并共同抑制miR-200的表达，形成对上皮表型的双重负调控；ZEB1还可通过自我正反馈机制，进一步强化其促EMT功能。miR-200家族作为EMT的关键抑制因子，通过直接抑制ZEB1表达以维持上皮表型，同时也可抑制PD-L1的表达。CDH1作为重要的上皮细胞粘附分子，其表达受ZEB1直接抑制，从而导致细胞极性丧失和迁移侵袭能力增强。PD-L1作为免疫抑制蛋白，除受miR-200负调控外，还可反向抑制CDH1表达，从而使得miR-200能够通过抑制PD-L1来间接促进CDH1的表达。值得注意的是，miR-200、ZEB1与CDH1构成核心调控轴，而miR-200和CDH1对ZEB1的抑制又形成重要的平衡机制。因此，这些相互作用共同构成一个动态平衡的网络，调控上皮与间质状态之间的可逆转换，最终决定EMT的进程。

Figure 2. Cross-talk between EMT regulators. Orange arrows denote promoting interactions, while blue T-shaped lines signify inhibitory interactions

图2. EMT调节因子间的相互作用，橙色箭头表示促进作用；蓝色T形线表示抑制作用

基于这些调节机制，本文建立了以下关于miR-200、ZEB1和CDH1三个关键分子的无量纲微分方程：

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\frac{a\cdot {\theta }^{n_1}}{{\theta }^{n_1}+y^{n_1}}-k_{1}x,\\ \frac{\text{d}y}{\text{d}t}=\frac{b\cdot y^{n_2}}{{\theta }^{n_2}+y^{n_2}}+\frac{c{\theta }^{n_3}}{{\theta }^{n_3}+x^{n_3}}+\frac{d{\theta }^{n_4}}{{\theta }^{n_4}+z^{{}^{n_4}}}-k_{2}y,\\ \frac{\text{d}z}{\text{d}t}=\frac{m{\theta }^{n_5}}{{\theta }^{n_5}+y^{n_5}}+\frac{n{x}^{n_6}}{{\theta }^{n_6}+x^{n_6}}-k_{3}z,\end{array}$ (1)

在本模型中，$x$ 、$y$ 和$z$ 分别表示miR-200、ZEB1和CDH1的浓度，参数$a$ 代表SLUG和ZEB1对miR-200的抑制强度，$c$ 为miR-200对ZEB1的抑制强度，$d$ 为CDH1对ZEB1的抑制强度，$n$ 为miR-200对CDH1的激活强度，$k_1$ 为miR-200的自降解率，$k_2$ 为ZEB1的自降解率，$k_3$ 为CDH1的自降解率，$\theta$ 是阈值，$n_1,n_2,n_3,n_4,n_5,n_6$ 是Hill系数。为在有限维度下刻画EMT调控网络的核心动力学特征，本研究对SLUG与ZEB1的部分作用进行了适度简化。具体而言$b$ 为SLUG和ZEB1对ZEB1的激活强度。该处理基于以下动力学与生物学假设：SLUG作为EMT上游诱导因子可上调ZEB1，而ZEB1又能通过自激活维持其高水平，二者对ZEB1启动子的协同激活在模型时间尺度内可视为外源缓变效应，因此归一化为单一等效激活参数b；同时，考虑到CDH1启动子上存在多个EMT转录因子结合位点，SLUG与ZEB1可通过直接结合或招募共抑制复合物实现功能一致的协同抑制作用，因此在建模中将其整合为一个等效“抑制模块”，假设其抑制效应可线性叠加并归一化为等效强度m。通过上述合理简化，模型在保留核心回路非线性动态的同时，显著提升了数学可处理性。

$\frac{a\cdot {\theta }^{n_1}}{{\theta }^{n_1}+y^{n_1}}$ 表示ZEB1对miR-200的抑制过程，$k_{1}x$ 表示miR-200的降解，$\frac{b\cdot y^{n_2}}{{\theta }^{n_2}+y^{n_2}}$ 表示ZEB1的自激活过程，$\frac{c{\theta }^{n_3}}{{\theta }^{n_3}+x^{n_3}}$ 表示miR-200对ZEB1的抑制过程，$\frac{d{\theta }^{n_4}}{{\theta }^{n_4}+z^{{}^{n_4}}}$ 表示CDH1对ZEB1的抑制过程，$k_{\text{2}}x$ 表示ZEB1的降解，$\frac{m{\theta }^{n_5}}{{\theta }^{n_5}+y^{n_5}}$ 表示ZEB1对CDH1的抑制过程，$\frac{n{x}^{n_6}}{{\theta }^{n_6}+x^{n_6}}$ 表示miR-200对CDH1的抑制过程，$k_{\text{3}}x$ 表示CDH1的降解[27]。

上述过程，将SLUG与自身对ZEB1的激活强度设为参数b，并将SLUG与ZEB1对CDH1的抑制强度设为参数m，这一简化不会改变系统的稳态拓扑结构，同时，Hill系数在2~6之间的变化不会破坏双负反馈环的双稳态特性。模型参数的适度扰动亦不会影响EMT状态从E型到M型的过渡趋势。采用一阶降解与Hill激活/抑制函数的组合，不仅符合生物学原理，还表现出优良的鲁棒性，保证了模型的稳定性。

3. 分岔分析

3.1. 平衡点稳定性分析

细胞在上皮型与间充质型表型之间的命运选择，本质上由miR-200、ZEB1与CDH1组成的调控网络所形成的动态平衡状态决定。该系统的稳定性若发生定性改变，即出现动力学分岔，往往预示着细胞表型的切换。为此，本节首先对系统的正平衡点进行局部稳定性分析，并在此基础上，深入探讨诱发鞍–结分岔的条件。

设系统的正平衡点为$E^{*}\left(x^{*},y^{*},z^{*}\right)$ ，系统的在正平衡点$E^{*}$ 处的雅可比矩阵$J\left(E_{*}\right)$ 如下所示：

$J\left(E_{*}\right)=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right)=\left(\begin{array}{ccc}-k_1& \frac{-a\cdot n_1\cdot {\theta }^{n_1}\cdot {\left(y^{*}\right)}^{n_1-1}}{{\left({\theta }^{n_1}+{\left(y^{*}\right)}^{n_1}\right)}^2}& 0\\ \frac{-c\cdot n_3\cdot {\theta }^{n_3}\cdot {\left(x^{*}\right)}^{n_3-1}}{{\left({\theta }^{n_3}+{\left(x^{*}\right)}^{n_3}\right)}^2}& -k_2+\frac{b\cdot n_2\cdot {\theta }^{n_2}\cdot {\left(y^{*}\right)}^{n_2-1}}{{\left({\theta }^{n_2}+{\left(y^{*}\right)}^{n_2}\right)}^2}& \frac{-d\cdot n_4\cdot {\theta }^{n_4}\cdot {\left(z^{*}\right)}^{n_4-1}}{{\left({\theta }^{n_4}+{\left(z^{*}\right)}^{n_4}\right)}^2}\\ \frac{n\cdot n_6\cdot {\theta }^{n_6}\cdot {\left(x^{*}\right)}^{n_6-1}}{{\left({\theta }^{n_6}+{\left(x^{*}\right)}^{n_6}\right)}^2}& \frac{-m\cdot n_5\cdot {\theta }^{n_5}\cdot {\left(y^{*}\right)}^{n_5-1}}{{\left({\theta }^{n_5}+{\left(y^{*}\right)}^{n_5}\right)}^2}& -k_3\end{array}\right),$

其中，

$\begin{array}{l}{p}_1=\frac{-n_1\cdot {\theta }^{n_1}\cdot {\left(y^{*}\right)}^{n_1-1}}{{\left({\theta }^{n_1}+{\left(y^{*}\right)}^{n_1}\right)}^2},p_2=\frac{-n_3\cdot {\theta }^{n_3}\cdot {\left(x^{*}\right)}^{n_3-1}}{{\left({\theta }^{n_3}+{\left(x^{*}\right)}^{n_3}\right)}^2},\\ p_3=\frac{n_2\cdot {\theta }^{n_2}\cdot {\left(y^{*}\right)}^{n_2-1}}{{\left({\theta }^{n_2}+{\left(y^{*}\right)}^{n_2}\right)}^2},p_4=\frac{-n_4\cdot {\theta }^{n_4}\cdot {\left(z^{*}\right)}^{n_4-1}}{{\left({\theta }^{n_4}+{\left(z^{*}\right)}^{n_4}\right)}^2},\\ p_5=\frac{n_6\cdot {\theta }^{n_6}\cdot {\left(x^{*}\right)}^{n_6-1}}{{\left({\theta }^{n_6}+{\left(x^{*}\right)}^{n_6}\right)}^2},p_6=\frac{-n_5\cdot {\theta }^{n_5}\cdot {\left(y^{*}\right)}^{n_5-1}}{{\left({\theta }^{n_5}+{\left(y^{*}\right)}^{n_5}\right)}^2},\end{array}$

进而得到雅可比矩阵如下：

$J\left(E_{*}\right)=\left(\begin{array}{ccc}-k_1& a\cdot p_1& 0\\ c\cdot p_2& -k_2+b\cdot p_3& d\cdot p_4\\ n\cdot p_5& m\cdot p_6& -k_3\end{array}\right),$

对应的特征方程可表示为：

$\det \left(\lambda E-J\left(E^{*}\right)\right)=0,$

展开得到

${\lambda }^3+q_1{\lambda }^2+q_2\lambda +q_3=0,$ (2)

其中，

$\begin{array}{l}{q}_1=-b{p}_3+k_1+k_2+k_3,\\ q_2=-ac{p}_1{p}_2-dn{p}_4{p}_6-b{k}_1{p}_3-b{k}_3{p}_3+k_1{k}_2+k_1{k}_3+k_2{k}_3,\\ q_3=-adn{p}_1{p}_4{p}_5-ac{p}_1{p}_2{k}_3-dn{p}_4{p}_6{k}_1-b{p}_3{k}_1{k}_3+k_1{k}_2{k}_3.\end{array}$

根据Routh–Hurwitz准则[31]，系统在正平衡点$E^{*}\left(x^{*},y^{*},z^{*}\right)$ 处渐近稳定的条件为：

$\begin{array}{l}{\Delta }_1=q_1>0,\\ {\Delta }_2=\left|\begin{array}{cc}{q}_1& q_3\\ 1& q_2\end{array}\right|>0,\\ {\Delta }_3=\left|\begin{array}{ccc}{q}_1& q_3& 0\\ 1& q_2& 0\\ 0& q_1& q_3\end{array}\right|=q_3{\Delta }_2>0.\end{array}$ (3)

3.2. 鞍–结点分岔的存在性

假设系统有一对特征值${\lambda }_1$ 和${\lambda }_{\text{2}}$ ，其中一个是正实部，另一个是负实部。根据特征方程(2)，系统发生鞍结点分岔的条件是在两个特征值合并时，产生一个零特征值。具体地，要求以下条件[32]：

(1) 存在一个参数使得在平衡点处特征方程有一个单零特征值，且其余特征值满足实部不等于0。有

$q_3=0,$

即

$adn{p}_1{p}_4{p}_5+ac{p}_1{p}_2{k}_3+dn{p}_4{p}_6{k}_1+b{p}_3{k}_1{k}_3=k_1{k}_2{k}_3.$

设特征值${\lambda }_1=0$ ，特征方程(2)变为

$\lambda \left({\lambda }^2+q_1\lambda +q_2\right)=0.$

剩余两个特征值满足：

${\lambda }^2+q_1\lambda +q_2=0,$

根据求根公式：

${\lambda }_{2,3}=\frac{-q_1\pm \sqrt{q_1^2-4q_2}}{2}.$

其余特征值满足实部不等于0，即

$q_1\ne 0$

(2) 系统(1)的参数变化导致平衡点的产生或消失，满足横截条件。横截条件可以表示为：

$\frac{\text{d}}{\text{d}\mu }\det \left(J\left(E^{*}\right)\right)\ne 0,$

其中，$\mu$ 表示系统参数中$a,b,c,d,m,n$ 中的一个。

令

$F={\lambda }^3+q_1{\lambda }^2+q_2\lambda +q_3=0,$

F对u求全微分得：

$\frac{\partial F}{\partial \mu }=3{\lambda }^2\frac{\partial \lambda }{\partial \mu }+\frac{\partial q_1}{\partial \mu }{\lambda }^2+2q_1\lambda \frac{\partial \lambda }{\partial \mu }+\frac{\partial q_2}{\partial \mu }\lambda +q_2\frac{\partial \lambda }{\partial \mu }+\frac{\partial q_3}{\partial \mu }=0.$

因为${\lambda }_1=0$ 代入上式得到

$\frac{\partial F}{\partial \mu }=q_2\frac{\partial \lambda }{\partial \mu }+\frac{\partial q_3}{\partial \mu }=0,$

即

$\frac{\partial \lambda }{\partial \mu }=-\frac{\frac{\partial q_3}{\partial \mu }}{q_2},$

要使$\frac{\partial \lambda }{\partial \mu }\ne \text{0}$ ，即

$\frac{\partial q_3}{\partial \mu }\ne 0$

4. 数值分析

本节将根据以上解析结果，从数值角度对系统的动力学行为进行深入探究，通过调节关键参数来改变系统平衡点的稳定性，这种调控将直接影响乳腺癌EMT过程中细胞表型状态的转变。接下来通过开展单参数与双参数分岔分析，系统探究系统(1)的动力学行为，旨在阐明乳腺癌EMT过程的关键参数响应规律。同时参数的选取参考已有研究成果，以保证研究方法的科学性与生物学意义[27] [28]。若没有特殊说明参数值见表1。

Table 1. Parameter values

表1. 参数值

参数	数值	参数	数值	参数	数值	参数	数值	参数	数值	参数	数值
a	0.2	b	0.05	c	0.1	d	0.3	m	0.5	n	0.5
n1	3	n2	2	n3	6	n4	2	n5	2	n6	2
k1	0.1	k2	0.1	k3	0.1	θ	0.5

固定z画出系统(1)的xOy面上的相图见图3，确定性系统(1)存在三个平衡点：两个渐近稳定的焦点，对应乳腺癌EMT过程中的两种表型状态，包括上皮表型(Epithelial Phenotype, E)和间充质表型(Mesenchymal Phenotype, M)，以及一个鞍点，对应乳腺癌EMT过程中的上皮/间充质混合表型(Epithelial/Mesenchymal Hybrid Phenotype, E/M)。

Figure 3. Phase portrait with direction field

图3. 带方向场的相图

4.1. 单参数分岔

本节中将讨论$a,m$ 的单参数分岔，在乳腺癌EMT过程中，a促进ZEB1的积累，从而驱动细胞向间质表型转化。然而，m会削弱这一平衡，导致CDH1表达下降，进一步加速EMT并增强肿瘤细胞的迁移和侵袭能力。在本节中，用MATCONT画出了系统(1)的分岔图，以说明miR-200、ZEB1和CDH1调控乳腺癌EMT过程。

4.1.1. a作为分岔参数

系统(1)关于参数a的单参数分岔图见图4，其中图4(a)、图4(b)、图4(c)分别是变量x，y，z关于参数a的分岔图。从图上可以看出，当参数a从0变化到5时，系统发生鞍结分岔，LP1和LP2对应的参数值分别为a = 0.07和a = 4.43，详细的分岔图如下。

(a) (b)

(c)

Figure 4. One-parameter bifurcation diagram with respect to a

图4. 关于a的单参数分岔图

图4展示了系统关于参数a的分岔图，揭示了EMT过程中的多稳态行为及其生物学意义。当$0<a<0.07$ 时，系统处于上皮状态，miR-200表达高，ZEB1低，CDH1表达高，细胞具有良好的连接性和低迁移性，维持典型的上皮表型。随着a增大，在$0.07<a<4.43$ 时，系统处于上皮/间充质混合表型，表征miR-200、ZEB1和CDH1处于不同水平，显示出细胞表型的高度可塑性，是EMT转化的关键窗口。尤其在靠近$a\approx 0.07$ 的区域，系统可稳定于部分EMT状态，此时miR-200和ZEB1处于中等水平，CDH1表达下降但未完全丧失，细胞具备较强的迁移性与侵袭性。当$a>4.43$ 时，系统转入间充质状态，miR-200表达低，ZEB1高，CDH1被抑制，细胞丧失极性，迁移性增强，表现出典型的EMT完成状态。因此，参数a的调控可决定细胞是否启动EMT及其程度，降低a可抑制ZEB1表达、恢复CDH1，从而阻止细胞向间充质状态转化。基于此机制，临床可通过抑制SLUG或上调miR-200等手段干预EMT，有望实现肿瘤转移的抑制与精准治疗。

4.1.2. m作为分岔参数

系统(1)关于参数m的单参数分岔图见图5，其中图5(a)、图5(b)、图5(c)分别是变量x，y，z关于参数m的分岔图。从图上可以看出，当参数m从0.8变化到2时，系统发生鞍结分岔，LP1和LP2对应的参数值分别为m = 1.78和m = 1.83，详细的分岔图如下。

图5展示了系统关于参数m的单参数分岔图，随着m从0.8增加到2，系统经历两个鞍结点分岔点对应不同的细胞状态转变。当$0<m<1.78$ 时，系统稳定于上皮状态，CDH1表达高，miR-200水平较高，ZEB1受限，细胞结构完整，迁移性弱；在$1.78<m<1.83$ 区间，系统处于上皮/间充质混合表型，即同时存在上皮、部分EMT和间充质状态，细胞表型高度可塑，是EMT动态调控的关键窗口期，尤其在接近$m\approx 1.78$ 时，细胞可处于部分EMT状态，表现出较强的侵袭与转移能力。当$m>1.83$ 时，系统进入间充质状态，CDH1表达极低，ZEB1显著上调，miR-200受抑，表征EMT完全激活，细胞高度侵袭。生物学上，LP1和LP2分别标志着系统从混合状态向间充质状态或上皮状态的跃迁。临床上，降低参数m的值意味着减弱SLUG/ZEB1对CDH1的抑制，有助于恢复CDH1表达、稳定上皮状态，从而抑制EMT或诱导其逆转为MET。

(a) (b)

(c)

Figure 5. One-parameter bifurcation diagram with respect to m

图5. 关于m的单参数分岔图

4.2. 双参数分岔

本节将展示参数对a-c，a-d，b-m及c-m的双参数分岔分析结果，在乳腺癌EMT调控网络中，参数a与c形成相互拮抗的作用关系，a促进间质转化而c维持上皮特性。通过MATCONT绘制的双参数分岔图显示，当c/a比值超过临界阈值时，系统将稳定在上皮态；反之则转向间质态。类似地，a-d组合中d的增强可抵消a的促EMT效应，b-m组合中b与m的相互作用揭示了表型转换的能量势垒，而c-m组合则体现了miR-200通过双重通路协同维持上皮稳态的机制。这些双参数分析不仅揭示了EMT调控网络的非线性特征，更为开发组合靶向治疗策略提供了重要理论依据。

4.2.1. a-c作为分岔参数

系统(1)关于参数a和c的双参数分岔图见图6，其中图6中上部分图是虚线框的放大部分。从图上可以看出，当参数a从0变化到0.12，c从0变化到0.1时，系统发生分岔从左到右对应的a = 0.05，0.10，c = 0.02，0.03，详细的分岔图如下。

Figure 6. Two-parameter bifurcation diagram in the (a, c)-plane

图6. 关于a和c的双参数分岔图

图6揭示了由SLUG-ZEB1-CDH1调控网络驱动的EMT过程具有明显的非线性多稳态特征，其中参数a与参数c共同决定细胞表型的稳定状态。在绿色封闭区域内，系统处于混合状态，表现出显著的表型可塑性，可作为EMT与MET共存及动态转换的“可逆窗口”；而在白色区域内，系统呈现单稳态：上皮状态通常与较低的a值和较高的c值相关联，间充质状态则与较高的a值和较低的c值相对应，此时细胞表型趋于稳定且不易发生逆转。生物学上，CP1和CP2标志着EMT从可逆向不可逆转化的边界，对识别治疗干预的关键时机具有重要意义。临床上可通过降低参数a、提高参数c，如抑制SLUG、上调miR-200、增强CDH1表达，阻断肿瘤细胞向间充质状态转化，降低转移风险；对混合态细胞的早期干预可促其逆转为上皮态，提升治疗效果。

4.2.2. a-d作为分岔参数

系统(1)关于参数a和d的双参数分岔图见图7，其中图7中上部分图是虚线框的放大部分。从图上可以看出，当参数a从0.5变化到2，d从0变化到0.5时，系统发生分岔，系统发生分岔从左到右对应的a = 0.68，0.9，d = 0.11，0.18，详细的分岔图如下。

图7所示的以参数a和d为坐标的双参数分岔图揭示了EMT过程中的多稳态动力学特征，在绿色封闭区域内，系统处于混合状态，表现出显著的表型可塑性；而在白色区域内，系统呈现单稳态：上皮状态通常与较低的a值和较高的d值相关联，间充质状态则与较高的a值和较低的d值相对应，此时细胞表型趋于稳定且不易发生逆转。尖点分岔CP1和CP2划定了这一可逆性与不可逆性之间的边界，具有重要的生物学与临床意义。临床上，若希望抑制EMT、维持上皮特征，可通过降低a和d (如抑制SLUG/ZEB1、增强CDH1功能或上调miR-200)实现干预；而对处于混合态的细胞群体，若能在绿色区域内及时识别并调控，有望逆转其向上皮态转化，恢复药物敏感性，提升治疗效果。因此，该图不仅体现了EMT调控的复杂性，也为精准治疗提出了一个可检验的动力学假说和参数干预策略。

Figure 7. Two-parameter bifurcation diagram in the (a, d)-plane

图7. 关于a和d的双参数分岔图

4.2.3. b-m作为分岔参数

系统(1)关于参数b和m的双参数分岔图见图8，从图上可以看出，当参数b从0变化到0.075，m从0变化到2.2时，系统发生分岔从左到右对应的b = 0.07，m = 2.04，详细的分岔图如下。

Figure 8. Two-parameter bifurcation diagram in the (b, m)-plane

图8. 关于b和m的双参数分岔图

图8所示的以参数c和m为坐标的双参数分岔图揭示了EMT过程中的多稳态调控机制。在绿色封闭区域内，系统处于混合状态，表现出显著的表型可塑性；而在白色区域内，系统呈现单稳态：当系统处于高c值与低m值的参数组合时稳定于上皮状态，处于低c值与高m值的参数组合时则稳定于间充质状态，此时细胞表型趋于稳定且不易发生逆转。尖点分岔CP1与CP2划定了细胞状态由可逆向不可逆转化的临界边界，具有重要的生物学和临床指导意义。临床上，通过增强CDH1对ZEB1的抑制作用和减弱ZEB1对CDH1的抑制作用，可有效阻断或逆转EMT过程，从而抑制肿瘤细胞的转移潜能，为精准治疗提供理论支持和干预策略。

4.2.4. c-m作为分岔参数

系统(1)关于参数c和m的双参数分岔图见图9，从图上可以看出，当参数c从0变化到0.11，m从0变化到2.1时，系统发生分岔从左到右对应的c = 0.02，0.11，d = 0.45，2.03，详细的分岔图如下。

Figure 9. Two-parameter bifurcation diagram in the (c,m)-plane

图9. 关于c和m的双参数分岔图

图9所示的以参数c和m为坐标的双参数分岔图揭示了EMT过程中的多稳态调控机制。在绿色封闭区域内，系统处于混合状态，表现出显著的表型可塑性；而在白色区域内，系统呈现单稳态：在高c值与低m值的参数组合下系统稳定于上皮状态，在低c值与高m值的参数组合下则稳定于间充质状态，此时细胞表型趋于稳定且不易发生逆转。尖点分岔CP1与CP2划定了细胞状态由可逆向不可逆转化的临界边界，具有重要的生物学和临床指导意义。临床上，通过增强CDH1对ZEB1的抑制作用和减弱ZEB1对CDH1的抑制作用，可有效阻断或逆转EMT过程，从而抑制肿瘤细胞的转移潜能，为精准治疗提供理论支持和干预策略。

5. 结语

本文通过构建一个包含miR-200、ZEB1和CDH1关键调控回路的三维微分方程模型，定量地模拟并刻画了乳腺癌细胞在EMT过程中的分子动态演变。研究结合单参数与双参数分岔分析，系统性地揭示了该网络在不同调控强度下的多稳态与转换阈值，描绘出系统的全局稳定性图谱。

研究表明，降低参数a与m可显著抑制系统的分岔行为，从而使细胞更稳定地维持于高上皮特性的状态。这一发现从动力学机制上阐释了保持上皮表型有助于抑制肿瘤的侵袭与转移潜力。与之相对，关键参数组合a-c、a-d、b-m、c-m的协同变化会引发分岔点的分割与移动，导致系统状态发生急剧跃迁，模拟了细胞在EMT过程中从上皮态向间质态的不可逆转换。

在分子机制层面，本研究清晰地揭示了EMT的核心调控逻辑：miR-200通过抑制转录因子ZEB1的表达，间接上调黏附蛋白CDH1，从而维持细胞的上皮特性和细胞间连接；反之，ZEB1的高表达则会形成一个双负反馈回路，强力抑制miR-200与CDH1，驱动细胞获得间质细胞的迁移与侵袭能力。基于此，理论模型预测，通过靶向干预(如施用miR-200模拟物或ZEB1小分子抑制剂)来调控模型中的关键参数，有望规避由“尖点分岔”所触发的系统状态突变，从而为阻断乳腺癌的转移进程提供新颖的治疗策略。

此外，本研究提出了一个可检验的动力学假说：miR-200/ZEB1/CDH1回路构成了一个决定细胞命运的双稳态开关，其调控机制具有明确的可操作性与可验证性。这不仅为开发针对miR-200或ZEB1的靶向药物提供了坚实的理论依据，更指向了一条精准医疗的应用路径。在未来临床实践中，通过检测患者肿瘤样本中这些标志分子的表达水平，可以精准判别其肿瘤细胞所处的功能状态。在此基础上，联合调控多个关键参数，制定个性化的、旨在将细胞稳定锁定于上皮状态的综合治疗方案，将有望更有效地遏制乳腺癌的侵袭与转移。

6. 讨论

本研究构建了miR-200-ZEB1-CDH1核心调控回路的动力学模型，系统揭示了该网络在乳腺癌EMT过程中表现出的双稳态与分岔特性，为理解细胞命运转换的动力学机制提供了新的理论框架。这一简化模型成功捕捉了EMT调控的核心特征。

在此基础上，未来的工作可从以下几个方向进一步深化：第一，在现有确定性模型基础上引入随机微分方程，探索噪声环境下EMT转换的概率特性；第二，整合转录后修饰、表观遗传调控等多层次调控因素，考察核心回路在复杂生理环境中的稳健性；第三，将本模型作为核心调控模块嵌入更复杂的多尺度模型中，实现从分子机制到组织水平的跨尺度模拟。

这些拓展研究将有助于构建更接近生物实况的EMT动态模型，为发展靶向EMT过程的抗乳腺癌转移策略提供更精准的理论指导，同时也为复杂生物网络的建模研究提供新的思路和方法。

基金项目

云南省教育科学规划(高等学校教师教育联盟)教师教育专项课题(GJZ2308)；云南省2024年本科高校教育教学改革研究项目(JG2024024)；云南师范大学2024年课程思政建设项目《程序设计与算法语言》。

NOTES

^*通讯作者。