从稳定性角度看细菌群落的噪声特性

裴启明,周祎凡,蒋 龙

(长江大学物理与光电工程学院,湖北 荆州 434023)

0 引言

在现代医学中,抗生素的引入为许多疾病提供了良好的治疗,尤其是由细菌引起的严重感染,从而挽救了无数人的生命.然而,抗生素的广泛使用和滥用导致其耐药性逐渐增强,而发现的新抗生素却越来越少,这给全球健康带来了危机[1].

持久性被认为是抗生素耐药性的根本原因[2],它是一种与异质细菌种群相关的可逆表型转换[3],可以引起抗生素敏感性的非遗传编码和可逆损失[4].目前,消除持久性已成为控制抗生素耐药性的一种新方法.

术语“持久性”由J.W. Bigger[5]于1944年提出.他将金黄色葡萄球菌的遗传同质菌株长时间暴露于杀菌浓度的青霉素中,发现一小部分细菌能逃脱抗生素的杀伤而存活下来,但其后代仍对抗生素敏感,这就是细胞的“持久性”,这一小部分细菌被称为“持留细胞”(persistent cells,简写为PCs).PCs只是暂时对抗生素有抗药性,一旦抗生素压力下降,它们就会恢复为正常细胞(normal cells,简写为NCs)[6-7].PCs不是遗传变异,而是表型转换[8-9].为了响应环境触发因素(如资源紧缺或抗生素)的存在,NCs和PCs可以相互转化[5],从而形成双向表型转换级联.

目前,已在人类病原体(如金黄色葡萄球菌、结核分枝杆菌和铜绿假单胞菌)、真核微生物(如白色念珠菌和酿酒酵母)甚至肿瘤细胞群中发现了PCs[10],由于PCs可以在抗生素消失后重新启动细胞分裂[11],因此它们极大地增加了慢性感染的难治性[7],它们在临床上的相关性已在实验中得到了验证[12].因此,消除PCs的策略会改善感染治疗的结果[13].

一方面,在定量生物学中,在生物分子水平上的随机波动被称为“噪声”或“涨落”[14-16].在细菌群落中,NCs和PCs除了相互转化外,还可以自我增殖、程序性死亡.这些反应均随机发生,从而产生内噪声.微环境的变化将产生外部噪声.由此可见,噪声无处不在.统计物理理论表明,噪声与近独立粒子的个数呈反比.因此,通过NCs和PCs的噪声特性也可以反映细胞个数的变化.

另一方面,稳定性是非线性系统的一个基本问题.对于物理过程的数学模型,只有其稳定解才对应于自然界中可观测到的实际运动[17].稳定性分析是研究非线性系统动态行为的关键步骤,如在病毒模型[18-19]、细胞神经网络[20]以及种群竞争模型[21-22]中,通过稳定性分析可以给出各种动力学行为的参数范围,从而为病毒传播的控制和疾病的临床治疗提供建议和方案[23].

在细菌群落中,如何不破坏细胞平衡态的稳定性能消除PCs呢?为了解决这个问题,本文从稳定性角度对在细菌群落中细胞的噪声特性进行了研究.在随机动力学中,Gillespie算法能够解释在实际生化反应中的噪声及其相关性[24].因此,本文理论公式的模拟数据均与Gillespie算法的模拟结果进行了比较,以验证理论的正确性.

1 在细菌群落中的双向表型转换级联

1.1 模型

在细菌群落中的双向表型转换级联[25]如图1所示.PCs和NCs都可以自我增殖,单位时间的增殖率分别为ai(i=1,2,以下相同,不再重复).由于营养物质的限制,ai会受到环境承载力的影响.另外,为了获得更多的资源,不同细胞间存在剥削竞争,使得一种细胞最大限度地将资源占为己有,从而耗尽另一种细胞的可用资源.因此,实际的自我增殖率为a1=k1(1-(N1+N2)/N0),a2=k2(1-(N1+N2)/N0),其中ki分别为NCs和PCs的最大自我增殖率,Ni为细胞个数,N0为环境可承载的最大细胞数.PCs和NCs可以相互转换,单位时间的转化率分别为r12、r21.其中PCs仅为一小部分,这表明转化率相对较小,r12、r21可取为常数.PCs和NCs都会死亡,单位时间的死亡率分别为b1、b2.通常,这2个死亡率bi也被视为常数.

注:正常细胞(NCs)和持留细胞(PCs)可以自我增殖、死亡和相互转化,相应的概率分别为ai、bi、rij (i, j=1,2;i≠j).

根据质量定律,这2种细胞的个数随时间的演化方程为

dN1/dt=a1N1-b1N1-r12N1+r21N2,

(1)

dN2/dt=a2N2-b2N2-r21N2+r12N1,

(2)

(3)

其中ρ=(k1s2-k2s1+((k1s2-k2s1)2+4k1k2r12r21)1/2)/

(2k2r21);s1=k1-b1-r12,s2=k2-b2-r21,它们分别为NCs和PCs的净固有增殖率.

1.2 对数增益

(4)

将式(4)与式(1)、式(2)相结合,可以得到级联中的4个对数增益分别为

H11=(k1〈N1〉/N0+ρr21)/(a1+ρr21),

H21=ρ(k1〈N1〉/N0-r21)/(a1+ρr21),

H12=(k2〈N2〉/N0-r12)/(ρa2+r12),

H22=(r12+ρk2〈N2〉/N0)/(ρa2+r12).

(5)

1.3 平均寿命

τ1=(b1+r12)-1,τ2=(b2+r21)-1.

(6)

1.4 Langevin方程

(7)

求解式(7),可得

xi=exp(-Hit/τi).

由此可见,xi随t呈指数变化.由于τi>0,所以平衡态的稳定性由Hi唯一决定.当Hi>0时,第i种细胞的平衡态是稳定的,否则不稳定.从这个意义上来说,Hi也可以理解为细胞平衡态的稳定因子,可被用来判断细胞受到干扰后能否恢复到平衡态.

根据式(7)中稳定性因子Hi的定义,易得

H1=H11+H21,H2=H12+H22.

(8)

由于NCs和PCs可以相互转化,一种细胞必定受到另一种细胞的影响,因此,一种细胞的稳定性因子不仅与自身有关,还与另一种细胞有关.如H11和H21共同确定NCs平衡态的稳定性,H22和H12共同确定PCs平衡态的稳定性,正如式(8)所示.若Hi>0,则NCs和PCs均处于稳定平衡态,否则,它们处于不稳定平衡态.即一旦NCs或PCs受到干扰,它们就将远离平衡态.

2 噪声公式

在随机动力学中,主方程可以给出每个生化反应的联合概率分布.然而,在大多数情况下,主方程不能精确求解.为了解决该问题,通常将主方程进行Ω展开[27],得到Fokker-Planck方程.而在平衡态下,Fokker-Planck方程的系数满足涨落-耗散关系[26]

MV+(MV)T+D=0,

(9)

其中V是协方差矩阵,其对角元Vii为方差,代表第i种细胞中的噪声;非对角元Vik(k≠i)为协方差,表示第i种细胞的噪声与第k种细胞的噪声之间的关联程度;M是漂移矩阵,矩阵元Mik(i,k=1,2)满足

Mik=-Hki/τi;

(10)

D是扩散矩阵,由式(1)～(2)可得

Dik=-(rik/〈Nk〉+rki/〈Ni〉)(k≠i).

(11)

(12)

(13)

由式(12)可知,任一细胞的噪声可分解为3部分:内噪声、来自另一种细胞的传输噪声以及2种细胞之间的转换噪声.不同于基因表达中的内噪声,NCs和PCs的内噪声不仅有纯内噪声,还有附加内噪声.纯内噪声与细胞的平均数成反比,这与公认的结论一致.附加内噪声由扰动环境产生,反映了NCs和PCs之间相互转化产生的间接影响.第j种细胞的内噪声可以通过级联进行传播,产生第i种细胞的传输噪声.转换噪声取决于2个转换率,反映了NCs和PCs之间相互转化的直接影响.

由式(13)可知,每个细胞的附加内噪声和2个细胞之间的转换噪声是协方差的来源.

3 模拟和结果

由式(3)可知,平衡细胞数取决于在级联中每个生化反应发生的概率.通常,死亡率bi和转化率rij为常数,PCs的最大自我增殖率k2变化也不大.因此,选择参数k1为控制变量.为了确保本文中所有参数的生物学意义,取0.01≤k1≤0.25.其他参数值分别为k2=0.121、b1=b2=0.02、r12=r21=0.001、N0=1 000[28],所有参数均为无量纲.

3.1 稳定性因子的特征

联合式(5)和式(8),可得稳定性因子Hi随NCs的最大自我增殖率k1的变化曲线(见图2)(其他参数值见正文,所有参数均为无量纲,以下各图的参数相同,不再重复).由此可见,Hi>0.因此,NCs和PCs的平衡态是稳定的,这意味着即使NCs或PCs受到ξ1(t)或ξ2(t)的干扰而瞬时偏离平衡态,它们也会通过自我调控而恢复到平衡态.另外,随着k1的增加,H1近似线性增加;H2尽管变化不大,但也略有增加.因此,Hi越大,NCs和PCs的平衡态越难以被破坏.

3.2 稳定性因子对平衡细胞数的影响

在图3(a)中,随着H1的增加,NCs的平衡数曲线呈“S”形,而PCs的平衡数曲线呈反“S”形.在图3(b)中,随H2的增加,NCs的平衡数曲线单调上升,而PCs的平衡数曲线单调下降.比较图3(a)和图3(b)可发现:尽管最初PCs的平衡数较大,但随着Hi的增加,这2种细胞的平衡数会发生交换(数值略有不同).最终,NCs的平衡数较大,而PCs的平衡数相对较小.

图2 稳定性因子Hi(i=1,2)随NCs的最大自我增殖率k1的变化图

注:实线是式(3)的理论预测,离散的空心符号是Gillespie算法的模拟结果,虚线和点划线是辅助线.

3.3 稳定性因子对协方差的影响

联立式(5)～(6)和式(13)可得协方差V12(=V21)随稳定性因子Hi的变化关系图(见图4中的实线).其中,离散的空心符号是Gillespie算法[24]的模拟结果.可见,这2种模拟结果吻合较好.

注:实线是式(13)的理论预测,离散的空心符号是Gillespie算法[24]的模拟结果,点划线是辅助线.

由图4可知,对于任意Hi,总有V12<0.因此,NCs的噪声与PCs的噪声负相关.这是因为在细菌群落中资源是有限的.为了获得足够的营养,这2种细胞之间发生剥削竞争.一种细胞将尽可能耗尽另一种细胞的可用资源,导致另一种细胞的数量减少,而自身的数量大大增加.

在图4(a)中,V12随H1的变化曲线为较对称的、开口向上的抛物线.当H1=4.81(虚线所示)时,V12的绝对值最大,这表明2个噪声之间的关联程度在整个参数范围内最强,细胞间的相互作用最大.在图4(b)中,V12随H2的变化曲线也是开口向上的抛物线,但几乎不对称.当H2=4.81(虚线所示)时,V12的绝对值也最大,细胞间的相互作用也最大.

因此,Hi值越大,细胞的平衡态越稳定,每个细胞也相对越独立.

3.4 稳定性因子对噪声的影响

联立式(5)～(6)和式(12)可得方差Vii(i=1,2)随稳定性因子Hi的变化关系图(见图5中的实线).其中,离散的空心符号是Gillespie算法[24]的模拟结果,这2个模拟结果也吻合得较好.

注:实线是式(12)的理论预测,离散的空心符号是Gillespie算法[24]的模拟结果,虚线和点划线为辅助线.

在图5(a)中,随着H1的增大,V11曲线呈“Z”字形,而V22曲线呈反“Z”字形.在H1=4.81附近(虚线所示)时2个噪声相等.在图5(b)中,随着H2的增加,V11单调下降,最后趋近于较小的常数,而V22单调增加,最后趋近于较大的常数.同样,在H2=4.81附近(虚线所示)时2个噪声相等.

比较图5(a)和图5(b)可发现:随着Hi的增加,2个噪声会发生交换(值略有不同),导致PCs的噪声V22最终远大于NCs的噪声V11.当然,只有当Hi的值达到某个阈值时,这种交换才会发生.一旦交换发生,2种噪声就将迅速超越强关联区(接近Hi=4.81),这2种细胞将达到相对独立的平衡态.

联合图5与图3可知,当Hi较大时,PCs的噪声V22远大于NCs的噪声V11,而PCs的平衡数远小于NCs的平衡数.这些均表明:PCs完全被NCs消除掉,感染即可被治愈.因此,Hi值越大,细胞的平衡态越稳定,感染越容易被治愈.

4 结论

在细菌群落中,为了维持细胞平衡态的稳定性并避免反复感染,本文从稳定性角度对NCs和PCs的噪声特性进行了研究.研究结果表明:当稳定性因子Hi达到一定阈值时,这2种细胞的噪声相互交换并迅速超越强关联区,最终导致PCs的噪声远大于NCs的噪声.此时,PCs几乎完全被NCs消除,从而感染被治愈.而且,稳定性因子Hi越大,NCs和PCs的平衡态越稳定,每种细胞相对越独立,感染越容易被治愈.该结果可为临床上治疗反复感染提供理论解决方案.

尽管笔者的一项工作[29]涉及在细菌群落中的噪声传播,但忽略了细胞平衡态的稳定性.本文从稳定性的角度研究噪声特性,从而得到了治愈反复感染的一种可能方案.然而,这种细菌群落的转换级联仅包含2种细胞表型.若级联中有更多的细胞表型(如乳腺癌细胞系中的管腔细胞、基底细胞和干细胞样细胞[30]),则不同细胞表型之间的相互作用会产生更加复杂的波动环境,稳定性角度下的噪声特性将会更为丰富,这将是我们未来的工作.