1.   引　言

癌症是种细胞信号通路的疾病, 严重威胁着人类生命健康^[1]. 抑癌因子p53作为细胞信号转导网络的枢纽, Puszyński等^[2]断言在人类癌症中约有一半会出现p53自身的突变, 而另一半会出现p53调控子的突变. 在静息细胞中, E3泛素连接酶Mdm2 (murine double minute 2) 使得p53多泛素化, 进而被蛋白酶体识别并快速降解; 此时p53的浓度被抑制在一个基底水平^[3]. 而在缺氧、紫外线照射、依托泊苷处理、电离辐射等细胞压力下, p53被稳定并激活 ^[4-7]. 活化的p53与靶基因的绑定亲和力增强, 诱导多种下游蛋白的表达. 这些下游蛋白参与到细胞周期停滞、DNA修复、细胞程序性死亡等生理过程中^[8]. 据实验报道, 细胞命运的信息被p53动力学编码^[9]. 因此, 还有一些问题亟需解决. 例如, 解码p53动力学所关联的细胞命运抉择的细致机制是个极具应用价值的问题; 研究p53蛋白信号网络中本身的动力学控制也是个很有意义的问题, 等等. 本文主要关注的是有关p53蛋白信号网络中的动力学控制问题, 即Mdm2在Ser395位点磷酸化这条信号通路对p53动力学调控的问题.

2000年, Lev Bar-Or等^[10]做了一项关于p53动力学的实验, 发现电离辐射下细胞群体中p53的衰减式脉冲, 也就是阻尼振荡. 后来的实验在单细胞中观测, 发现p53的振荡是无阻尼的, 并且脉冲的个数与辐射剂量正相关^[11,12]. 后续的很多理论模型围绕p53的反复式脉冲动力学开展, 例如, 张丽娟等^[13]用p53-Mdm2负反馈环附加时滞构建的振子模型揭示了群体细胞p53的阻尼振荡可能是单细胞p53振荡持续时间不一、脉冲数量差异的结果; Bottani和Grammaticos^[14]研究了带时滞的p53-Mdm2振子在p53, Mdm2半衰期参数控制下的分岔现象; Ma等^[15]建立了复杂的混合模型 (随机过程对接微分方程) 来模拟p53的数字式脉冲; Xia和Jia^[16]重建了包含p53振子模块和DNA修复模块的简化模型, 用更简单的方式解释了p53脉冲数正比于辐射剂量的可行性机制; Elia${\breve{\rm {s}}}$等^[17]探索了单细胞中Wip1 (wild type p53-induced phosphatase 1) 参与下p53发生振荡的常微分模型和基于反应扩散的偏微分模型; 毕远宏等^[18]研究了PDCD5 (programmed cell death 5) 和Mdm2最大生成速率调控下的p53全局动力学与稳定性等问题, 展示了p53动力系统中存在的丰富动力学现象; Hat等^[19]研究了p53基因网络中常见的正负反馈通路及反馈背后的双稳态、振荡等分岔机制; 王道光等^[20]提取了p53-Mdm2负反馈环和p53-Wip1-ATM (ataxia telangiectasia-mutated)负反馈环共存的模型, 分析了p53信号转导网络中的双节率现象; Ochab和Puszynski^[21]开发了一种分段线性方法, 并以p53振子为例对比了新方法与传统的非线性方法; 刘楠等^[22]重建了p53-Mdm2通路和p53-Wip1-ATM通路耦合的模型, 在Ser395位点磷酸化的Mdm2被考虑在内, 并探索了这种情景下Wip1管理的p53潜在动力学; 等等. 虽然p53-Mdm2之间的正反馈被提及, 但据我们所知, 尚未有细致的研究.

基于上述考虑, 本文重新审视了Mdm2促进p53的作用对p53动态的影响. 首先给出这个基因网络模型对应的动力学方程; 然后利用分岔分析的办法, 找出了关键的分岔点, 系统的动力学特征由分岔点划分的参数区间所决定; 最后根据分岔曲线, 在二维参数平面上找出了振荡的区域. 数值研究结果表明: p53-Mdm2之间正反馈强度的差异可能是p53面对两种都能激活ATM的细胞压力时出现不同动力学响应方式的关键原因; Wip1对这个基因网络的振荡会时而有利时而不利, 取决于Mdm2对p53的正调控强度; 在这个模型中, p53-Mdm2的正负双重反馈, 是振荡出现的根基, 也就是说p53-Mdm2之间的正负反馈对p53振荡而言是同等重要的. 总之, 振荡出现需要正负两种反馈强度接近、相互偶联.

2.   模型与方程

事实上, p53的动力学行为是压力类型依赖的^[23]. 振荡的p53动力学出现在那些能引起DNA双链断链 (DSB) 的细胞压力中, 例如电离辐射和依托泊苷刺激^[12,24]. ATM作为DSB的信号传感器, 接收到DNA损伤信号并迅速传导到p53-Mdm2模块^[25]. 一方面激活p53的转录活性, 另一方面通过催化Mdm2在Ser395位点的磷酸化促使Mdm2自我泛素化并被快速分解来增强p53的稳定性^[26]. 在Ser395位点磷酸化的Mdm2能结合p53的信使RNA, 使其翻译效率被大幅提升^[27]. 此外, Mdm2还有一种在Ser166或Ser168位点的磷酸化方式, 这类Mdm2的作用是促进p53失活和降解^[28]. Mdm2是p53的靶蛋白之一, 实验表明另一种p53的靶蛋白Wip1也对维持形状均匀的p53脉冲至关重要^[29]. 作为一种蛋白激酶, Wip1能促进ATM的去磷酸化, 从而抑制DNA损伤诱导下的ATM激活. 因此, p53-Mdm2之间既有正反馈环又有负反馈环, 而p53-Wip1-ATM仅形成了负反馈环, 如图1所示, 其中红色线形成的回路为负反馈环, 绿色线形成的回路为正反馈环.

图 1  模型示意图. 箭头线表示促进或状态转移; 杠头线表示抑制

Fig. 1.  Schematic diagrams of the model. Arrow-headed lines indicate promotion or state transition; bar-head lines indicate inhibition.

下载: 全尺寸图片 幻灯片

要构建基因网络的动力学模型, 首先要明确网络中的节点. 我们的模型中共有6个节点: 活性的ATM, 活性的p53, Wip1, 在Ser395位点磷酸化的Mdm2 (Mdm2p1), 未修饰的Mdm2 (Mdm2i), 在Ser166或Ser168位点磷酸化的Mdm2 (Mdm2p2). 根据模型所考虑的生化反应 (也就是ATM被DSB激活, 被Wip1灭活; p53的激活需要ATM, 合成受到Mdm2p1影响; Wip1基因的转录依赖于p53; Mdm2基因的转录也依赖p53; 不考虑能终止p53振荡的Akt通路^[30], Mdm2可自由地在Ser166或Ser168位点被磷酸化以及去磷酸化, 导致Mdm2i与Mdm2p2相互转化; Mdm2被ATM催化转变为Mdm2p1; Mdm2p1又能自发地去磷酸化, 变回Mdm2i), 可给出如图1 所示的模型框架. 网络中的每个节点对应着一个方程. 用方括号 “[ ]” 表示无量纲的蛋白质浓度, 在方括号里写了蛋白质的种类, 用一类 “S” 形函数 (广义的希尔函数) 来刻画蛋白质之间的相互调控, 则由这个基因网络“翻译”出的动力学方程组如下:

其中浓度单位为“C”, 时间单位为分钟(min), 希尔系数 $s_0 = s_1 = s_2 = s_3 = s_4 = s_6 = 4$, $s_5 = 2$, 其他的参数值及描述如表1^[22,30,31]所列. 在标准参数下p53的振荡周期在4—7个小时之间 (图2), 与实验数据一致.

表 1  模型中的参数

Table 1.  Parameters of the model

参数	参数描述	值	单位
$v_{\rm {atm}}$	ATM激活速率	1	${\rm{C} }\!\cdot\! {\rm{min} }^{-1}$
$d_{\rm {atm}}$	ATM最大失活速率	0.53	${\rm{min}}^{-1}$
$\rho_ 0$	Wip1依赖的ATM失活反应占据的比例常数	0.9	$\setminus$
$k_ 0$	Wip1催化的ATM失活反应中的酶半饱和浓度	1	${\rm{C}}$
$v_ {\rm{p53}}$	p53最大合成并激活的速率	0.95	${\rm{C} }\!\cdot\! {\rm{min} }^{-1}$
$\rho_ 1$	ATM依赖的p53激活反应占据的比例常数	0.9	$\setminus$
$k_ 1$	ATM催化的p53激活反应中的酶半饱和浓度	0.3	${\rm{C}}$
$\rho_ 2$	Mdm2p1依赖的p53合成反应占据的比例常数	0.98	$\setminus$
$k_ 2$	Mdm2p1催化的p53合成反应中的酶半饱和浓度	0.057	${\rm{C}}$
$d_ {\rm{p53}}$	p53最大降解速率	0.6	${\rm{min}}^{-1}$
$\rho_ 3$	Mdm2p2依赖的p53失活和降解反应占据的比例常数	0.97	$\setminus$
$k_ 3$	Mdm2p2催化的p53失活和降解反应中的酶半饱和浓度	0.09	${\rm{C}}$
$v_ {\rm{mdm2}}$	Mdm2最大生成速率	0.153	${\rm{C} }\!\cdot\! {\rm{min} }^{-1}$
$\rho_ 4$	p53依赖的Mdm2合成反应占据的比例常数	0.98	$\setminus$
$k_ 4$	p53催化的Mdm2合成反应中的酶半饱和浓度	4.43	${\rm{C}}$
$k_ {\rm{p}}$	Mdm2在Ser395位点的最大磷酸化速率	0.65	${\rm{min}}^{-1}$
$\rho_ 5$	ATM依赖的Mdm2磷酸化反应占据的比例常数	0.9	$\setminus$
$k_ 5$	ATM催化的Mdm2磷酸化反应中的酶半饱和浓度	1	${\rm{C}}$
$k_{\rm{ q}}$	Mdm2在Ser395位点的去磷酸化速率	0.24	${\rm{min}}^{-1}$
${k_ {{\rm{p}}*}}$	Mdm2在Ser166, Ser168的磷酸化速率	0.14	${\rm{min}}^{-1}$
${k_ {{\rm{q}}*}}$	Mdm2在Ser166, Ser168的去磷酸化速率	0.01	${\rm{min}}^{-1}$
$d_ {\rm{mdm2}}$	Mdm2本底降解速率	0.034	${\rm{min}}^{-1}$
$v_ {\rm{wip1}}$	Wip1最大生成速率	0.09	${\rm{C} }\!\cdot\! {\rm{min} }^{-1}$
$\rho_ 6$	p53依赖的Wip1合成反应占据的比例常数	0.99	$\setminus$
$k_ 6$	p53催化的Wip1合成反应中的酶半饱和浓度	0.8	${\rm{C}}$
$d_ {\rm{wip1}}$	Wip1降解速率	0.05	${\rm{min}}^{-1}$
f	Mdm2p1的降解速率与Mdm2本底降解速率的比值	3.58	$\setminus$
g	Mdm2p2的降解速率与Mdm2本底降解速率的比值	1.95	$\setminus$

下载: 导出CSV 

| 显示表格

图 2  p53的时间历程图

Fig. 2.  Time occurs of p53.

下载: 全尺寸图片 幻灯片

3.   结果与分析

3.1   关于正反馈强度参数${{k}}_{\bf{p}}$的分岔分析

Mdm2对p53正调控的必要途径是Mdm2在Ser395位点的磷酸化, 与参数$k_{\rm{p}}$对应. 换言之, $k_{\rm{p}}$ 的值反应了p53-Mdm2正反馈的强度. 首先在图3给出了p53平衡态与$k_{\rm{p}}$的函数关系, 即余维1分岔图. 在不动点分支上出现了三个分岔点: 鞍结分岔点 (SN)、鞍结同宿分岔点 (SNIC)、霍普夫分岔点 (HB). 在SN和SNIC之间的蓝线是不稳定的鞍点, 在SN和HB之间的红色线是不稳定的结点和焦点, 在SNIC和HB之外的黑色线代表p53稳定的稳态. 在SNIC和HB之间, 这个动力系统中仅存在稳定的极限环, 也就是p53动态行为仅有稳定振荡. 这里极限环的分支之所以没有被延拓是因为极限环上的折叠分岔点靠近HB, 因此系统的动力学由不动点分支上的分岔点基本确定. 当$k_{\rm{p}}$在(0, SNIC) 参数区间中时, p53最终会稳定在低稳态; 当$k_{\rm{p}}$在(SNIC, HB) 参数区间中时, p53动力学呈现反复的、非衰减的脉冲; 当$k_{\rm{p}}$在(HB, $+\infty$) 参数区间中, p53最终会达到一个较高的稳态.

图 3  关于参数$k_{\rm{p}}$的一维分岔图

Fig. 3.  One-dimensional bifurcation graph on parameter $k_{\rm{p}}$.

下载: 全尺寸图片 幻灯片

为了验证一维分岔分析的结果, 取了3个$k_{\rm{p}}$值在([p53], [Mdm2])相平面画出轨线(图4), 这里的[Mdm2]为三种形式Mdm2的浓度之和. 当$k_{\rm{p}} = 0.3$时(图3中的黑色虚线), 分岔图显示p53有稳定的低稳态, 在图4中对应着黑色轨道, 从原点开始最终停在了不动点. 当$k_{\rm{p}} = 0.65$时(图3中的红色虚线), 不动点是不稳定的焦点, 在图4中对应着红色轨道, 始于原点的轨迹最终收敛到极限环. 当$k_{\rm{p}} = 1$时(图3中的蓝色虚线), p53最终会稳定在较高的浓度, 在图4中对应着蓝色轨道, 从原点出发沿着渐进的环状轨迹最终到达不动点, 并不再移动. 相平面分析(图4)的结果与分岔分析(图3)的结果保持一致, 证实了分岔分析的正确性. 在正反馈强度很弱的时候, 极限环消失, 说明了Mdm2介导的正反馈对p53振荡发生的必要性. 此外, 正反馈强度过强时候振荡也会消失. 实际上, 正反馈过强时可通过调节一些负反馈的强度来抵消正反馈, 进而使得振荡复现, 详情见3.2和3.3小节.

图 4  相平面分析

Fig. 4.  Phase plane analysis.

下载: 全尺寸图片 幻灯片

上述的分析表明, 要获取二维参数平面上的振荡区域, 只需延拓SNIC和HB两个分岔点. 由于ATM的激活参数$v_{\rm{atm}}$是这个基因网络接收刺激信号的开端, 受到DSB的调控, 所以我们把$v_{\rm{atm}}$作为第二参数张成($k_{\rm{p}}$, $v_{\rm{atm}}$)参数平面(图5). 通过延拓SNIC分岔点(图5中的蓝色线)和HB分岔点(图5中的红色线), 得到二维参数平面上的两条分岔曲线, 由这两条分岔线围成的灰色区域为振荡区域. 图5表明, 当正反馈强度较弱, 即$k_{\rm{p}}$较小时, 例如$k_{\rm{p}} = 0.6$, p53的振荡行为对DNA损伤更鲁棒, 也就是说只要DNA损伤输入的信号超过某个阈值就能触发p53振荡. 这与Lahav等^[12]的实验结果一致, 即电离辐射强度(DNA损伤程度)不能改变p53脉冲式动力学的应答模式. 而正反馈强度较强, 例如$k_{\rm{p}} = 1$时, 过多或过少的DNA损伤都不能使p53振荡, 与Chen等^[24]的实验一致, 即依托泊苷的剂量(DNA损伤程度)需要适中才能出现p53振荡. 因此我们推断, 在引起DSB的两类刺激下, p53出现对DSB相异的鲁棒性有可能是p53-Mdm2之间正反馈强度的差异引起的, 这与Sun和Cui^[32]的观点(ATM诱导的Mdm2磷酸化通路有利于双模式的p53动力学)有相似之处.

图 5  二维参数平面($k_{\rm{p}}$, $v_{\rm{atm}}$)上的振荡区域

Fig. 5.  Oscillation area on the ($k_{\rm{p}}$, $v_{\rm{atm}}$) two-dimensional parameter plane.

下载: 全尺寸图片 幻灯片

3.2   p53-Wip1-ATM负反馈环与p53-Mdm2p正反馈环的联合影响

Batchelor等^[29]的实验发现, 在细胞中加入Wip1信使RNA的抑制剂后, p53振荡消失. 此时的p53出现单个大幅度脉冲, 类似于在紫外线照射实验中观测到的p53脉冲. 因此, Wip1通路作用强度的差异很可能是两种刺激下出现不同的p53动力学应答的根本原因^[23]. 为了探究Wip1和振荡发生的关系, 在($k_{\rm{p}}$, $v_{\rm{wip1}}$)参数平面上延拓出SNIC分岔线和HB分岔线; 同样, 灰色区域为振荡区域(图6). 当$v_{\rm{wip1}}$增大到一定程度时, 两条分岔曲线与x轴接近垂直, 此时Wip1的作用饱和, 也就是Wip1含量虽然充足, 但是反应底物活性ATM的含量较低, 使得ATM失活反应的最大速率被限制. 这点在模型中体现在关于Wip1的希尔函数中, 在Wip1的浓度很高的情形下, 继续增加Wip1浓度不再对这个基因网络的其他节点有影响. 显然, 在标准参数下($k_{\rm{p}} = 0.65$), Wip1是振荡出现的重要条件, 与文献[22]结论一致.

图 6  二维参数平面($k_{\rm{p}}$, $v_{\rm{wip1}}$)上的振荡区域

Fig. 6.  Oscillation area on the ($k_{\rm{p}}$, $v_{\rm{wip1}}$) two-dimensional parameter plane.

下载: 全尺寸图片 幻灯片

有趣的是, 当$k_{\rm{p}}$较小时, $v_{\rm{wip1}} = 0$也会发生振荡. 这意味着Wip1被抑制的情况下, 通过削弱p53-Mdm2之间的正反馈也能使得p53出现振荡, 在生物实验中应注意这一点. 为了进一步说明这些情况, 在两种$k_{\rm{p}}$值下做了[p53] 关于参数$v_{\rm{wip1}}$的余维1分岔图(图7). 当$k_{\rm{p}}$较小时, 如图6黑色小三角所指($k_{\rm{p}} = 0.5$), 此时Wip1可能不利于p53振荡; 在图7(a)中, 随着$v_{\rm{wip1}}$增大, 到超过SNIC分岔点时, 不动点为稳定的结点, 极限环消失. 而当$k_{\rm{p}}$较大时, 如图6红色小三角所 ($k_{\rm{p}} = 0.7$), 此时Wip1有利于p53振荡; 在图7(b)中, 随着$v_{\rm{wip1}}$增大, 到超过HB分岔点时, 不动点为不稳定的焦点, 极限环出现. 因此, 在这个模型中Wip1对系统振荡的影响是由p53-Mdm2之间的正反馈强度决定的. 非时滞模型的振荡需要恰当的正反馈强度, 默认参数下($k_{\rm{p}} = 0.65$)正反馈较强, 附加的p53-Wip1-ATM负反馈回路与p53-Mdm2p1正反馈回路拮抗, 所以能出现振荡.

图 7  关于参数$k_{\rm{wip1}}$的一维分岔图

Fig. 7.  One-dimensional bifurcation graph on parameter $k_{\rm{wip1}}$

下载: 全尺寸图片 幻灯片

3.3   p53-Mdm2正负反馈环的联合影响

最后探究p53-Mdm2之间的正负反馈协作对系统振荡出现的影响. 模型中假设Mdm2要实现对p53的抑制作用, 即破坏p53的稳定性和转录活性, 首先要经过Akt的磷酸化. 事实上, 在文献[30]的模型中, Akt的活性也受到p53的间接调控, 并且这种作用很可能引起p53的双稳态. 在本文考虑的模型中, ${ k_{{\rm{p}}*}}$被视为定常量, 是控制Mdm2抑制p53这条通路的关键参数. 如图8所示, 在($k_{\rm{p}}$, ${ k_{{\rm{p}}*}}$)参数平面中的振荡区域(灰色部分)与x轴和y轴都不相交, 这代表p53-Mdm2p1通路和p53-Mdm2p2通路都是振荡的要素. 负反馈回路是振荡发生的网络结构基础, 在很多p53振子的理论模型中, p53-Mdm2的负反馈都起到了重要作用^[13-16], 本文的模型中起主导作用的负反馈回路同样是p53-Mdm2. 参数$k_{\rm{p}}$的振荡区间(SNIC, HB)随着${ k_{{\rm{p}}*}}$的增大而扩宽, 进一步说明了p53-Mdm2p2负反馈环有利于p53的振荡. 随着负反馈强度的增加, 振荡所需要的正反馈强度也随之增加, 再次说明了振荡的出现需要两者同时参与、相互制约. 事实上, Zhang等^[33]提出, 用正负反馈回路耦合构建的p53振子模型能产生更具鲁棒性的振荡.

图 8  二维参数平面($k_{\rm{p}}$, ${ k_{{\rm{p}}*}}$)上的振荡区域

Fig. 8.  Oscillation area on the ($k_{\rm{p}}$, ${ k_{{\rm{p}}*}}$) two-dimensional parameter plane.

下载: 全尺寸图片 幻灯片

4.   结论与讨论

总而言之, 本文全面探讨了Mdm2对p53的正调控(即Mdm2p1促进p53的翻译)和p53振荡发生之间的关联. 首先利用生物事实给出一个p53基因网络的数学模型; 然后利用数值分岔分析的方法, 研究了刻画p53-Mdm2正反馈强度的参数$k_{\rm{p}}$: 1) k_p与DSB的协作关系, 提出了实验发现的p53对两类能引起DSB的刺激(电离辐射、依托泊苷处理)出现不同动力学响应模式的原因可能是$k_{\rm{p}}$的不同; 2) k_p与Wip1的协作关系, 指出了Wip1一方面在高强度正反馈下能促进振荡, 另一方面在低强度正反馈下能抑制振荡; 3) k_p与Mdm2p2的协作关系, 再次证实了p53-Mdm2的负反馈回路是振荡发生的结构基础. 本文强调的是正负反馈对p53基因网络发生振荡的同等重要性. 值得注意的是, 希尔函数和微分方程刻画的是宏观层次的模型, 今后的研究还需建立更为细致的多态模型^[34]. p53的振荡与细胞周期停滞或细胞凋亡等过程密不可分^[35], 希望本研究能为今后的实验设计和医学应用提供思路.

此外, 本文的结果是在非线性程度(协同性)很高的调节关系上得出的, 这里的非线性程度很高主要体现在数值较大的希尔系数上. 为了使结论更加可靠, 检测希尔系数较小时振荡发生的条件很有必要. 把$s_0$调到1; $s_1$, $s_2$和$s_3$下调到2; $s_4$, $s_6$的值不变 (因为p53以4聚体的形式绑定靶基因^[36], 其希尔系数理论上应取为4). 为了使p53出现振荡, $k_1$的值改为0.09, $k_3$的值改为0.05, 其他参数保持默认值, 在图9中作出关于参数$k_{\rm{p}}$的分岔图组. 从图9可以看出, 不同于强非线性的分岔现象是分岔点SN以及SNIC消失, 振荡区间变为两个HB之间的参数区间(HB1, HB2); 其他的振荡控制现象未发生变化. 也就是($k_{\rm{p}}$, $v_{\rm{atm}}$)参数平面上的振荡区域仍由两条 “$\lfloor$” 形分岔曲线所围成, 当$k_{\rm{p}}$较小时p53能在ATM过度激活的情况下保持振荡; ($k_{\rm{p}}$, $v_{\rm{wip1}}$)参数平面上的振荡区域由两条 “$\rfloor$” 形分岔曲线所围成, 存在一部分$k_{\rm{p}}$的值使得Wip1抑制振荡; ($k_{\rm{p}}$, $k_{{\rm{p}}*}$)参数平面上的振荡区域由两条“$/$”形直线所围成, p53-Mdm2的正负反馈偶联依然是这个基因网络振荡的“引擎”. 因此, 本文的定性结论可能仅取决于这个基因网络的结构, 对希尔函数的非线性程度并不敏感. 结论的可靠性得到了进一步的验证.

图 9  分岔图组

Fig. 9.  Bifurcation graph group.

下载: 全尺寸图片 幻灯片