AC-DC基因回路的多功能动力学全景解析

来源:Let there be multifunctionality: Uncovering the dynamical zoo of the AC-DC genetic circuit(Cell Systems)| 编译:DNA Lab Space | 原文许可:elsevier-subscription

📎 相关工具:蛋白可溶性优化分析

导读

基因调控网络通常以多稳性或振荡性被单独研究和工程化。本文以AC-DC回路——一个融合了经典双稳态开关与阻遏振荡器的最小三基因网络——为对象,结合热力学形式与贝叶斯推断,系统揭示了单诱导物版本即可呈现振荡与多稳性共存等多功能动力学。研究进一步对涌现行为进行分类、评估鲁棒性并分析临界慢化,发现该回路可产生逾三十种拓扑上不同的分岔图,挑战了"网络拓扑刚性决定动力学结局"的经典观点。

研究背景

基因调控网络(GRN)通过控制基因表达的时空组织来主导生物学特征。借助对外部信号与网络内转录因子的动态解读,GRN调控着多细胞生物的模式形成、环境适应以及代谢与生长的协调等复杂过程。即便仅由少数基因构成的小型GRN,其编排复杂过程的能力也依赖于转录与翻译机器的非线性响应,典型动态行为包括遗传振荡、混沌动力学、多稳基因开关、可兴奋性和完美适应等。

小型GRN虽已被广泛研究,但其可展现的动态行为全貌仍被显著低估。传统上网络拓扑被视为功能的首要决定因素,但这一图景远比想象中微妙。越来越多证据表明,共享同一拓扑的回路仅通过调节生物物理参数即可展现不同的动力学功能——这一现象称为多功能性,已在理论与实验层面得到研究。充分揭示这种动态库对于理解自然回路的能力以及设计合成回路均至关重要。本文即以此为出发点,采用比以往更具机制真实性的生物物理表述,系统探索AC-DC回路的多功能动力学全景。

研究方法

热力学形式体系(主模型)

为模拟基因调控动力学,我们采用热力学形式体系(亦称分数占据模型),以生物物理术语描述DNA、转录因子(TF)与RNA聚合酶的结合事件。在该框架中,基因表达对应于启动子处于结合RNA聚合酶状态的概率。这涉及枚举所有可能的启动子占据状态,并计算包含结合聚合酶的状态所占比例(图S1)。每个状态的相对丰度取决于其结合能,转录因子与信号调节不同构型的能量。最终得到系统中每个启动子(或节点)的转录概率函数。例如,节点i由概率φ_i(X, Y, Z, S)表征,该概率依赖于分子j对启动子i的结合亲和力K_ij,以及协同性参数c_i——描述诱导物信号结合后对聚合酶结合的能量促进。对于单诱导物模型(图S1),概率描述为:

φ_X(S, Z) = (K_XP·P + c_XPS·K_XP·K_XS·P·S) / (K_XP·P + c_XPS·K_XP·K_XS·P·S + (1 + K_XS·S)(1 + K_XZ·Z)^n)

φ_Y(X) = K_YP·P / (K_YP·P + (1 + K_YX·X)^n)

φ_Z(X, Y) = K_ZP·P / (K_ZP·P + (1 + K_ZY·Y)^n(1 + K_ZX·X)^n)

在该表述中,阻遏型TF阻止聚合酶结合,n为结合位点数。本文模拟中固定n = 3。此外,由于模拟中总聚合酶P保持恒定,不失一般性地设P = 1,以聚合酶浓度为单位测量结合亲和力K_iP。类似地,双诱导物情形下唯一变化是额外信号结合到启动子Y:

φ_Y(S, X) = (K_YP·P + c_YPS·K_YP·K_YS·P·S) / (K_YP·P + c_YPS·K_YP·K_YS·P·S + (1 + K_YS·S)(1 + K_YX·X)^n)

所得非线性转录概率φ_i可直接纳入常微分方程(Equation 1),从而模拟系统动力学。此外,在参数搜索中我们将表达水平归一化为X̃ = δ_X·X/α_X,Ỹ = δ_X·Y/α_Y,Z̃ = δ_X·Z/α_Z。进一步重参数化时间t̃ = δ_X·t,得到无量纲降解率δ̃_Y = δ_Y/δ_X和δ̃_Z = δ_Z/δ_X。类似地,结合亲和力相对于相应转录因子的最大表达水平进行重参数化(如K̃_YX = K_YX·α_X/δ_X)。为清晰起见,所有无量纲参数与变量的波浪号均省略。这些重参数化将自由参数数量从17个(单诱导物模型15个)减少至13个(单诱导物模型11个)。

【数据】结合位点数:n=3;聚合酶浓度:P=1(无量纲);自由参数:双诱导物13个(无量纲化后)/单诱导物11个;原始参数:双诱导物17个/单诱导物15个;原文未详述具体亲和力数值范围

分岔分析

为获得分岔图,我们求解无量纲系统(Equation 1)在不同S值下的稳态(X, Y, Z)及相应线性稳定性。AC-DC的拓扑允许我们通过数值求解一维方程Z = (1/δ_Z)·φ_Z(X(Z), Y(Z))来找到稳态。其中X(Z)和Y(Z)对应X(Z) = φ_X(Z),Y(Z) = (1/δ_Y)·φ_Y(X(Z*))。极限环通过计算从扰动不稳定点出发的长时间轨迹获得。分岔点定义为分支稳定性改变的点(如Hopf分岔)、稳态数量改变的点(如鞍结分岔),或极限环变为同宿轨道的点(如鞍环分岔、不变圆上的鞍结分岔)。

【数据】结合位点数:n=3;稳态求解:一维方程数值求根;分岔类型:Hopf、鞍结、鞍环、SNIC;原文未详述具体S扫描步长

参数探索

我们使用基于序列蒙特卡洛采样的近似贝叶斯计算(ABC)探索参数空间,先验分布为表S1所示参数范围内的平坦分布。距离函数d可理解为损失函数或得分函数的倒数,用于定义伪似然L,如d = -log L。这引入了基于稳态分析所得动力学的样本拒绝度量。距离函数决定参数集与期望动力学特性的距离,值越大表示未满足期望的动力学标准。我们在S = 1、S = 10²、S = 10³、S = 10⁴处评估相图(见图1E)以获得多功能性得分。每个参数集按距离函数d = d_mono(1) + d_cycle(10²) + d_bistable(10³) + d_mono(10⁴)评分,该函数评估每个参数集的单稳性(d_mono)、振荡性(d_cycle)和双稳性(d_bistable)。振荡行为在S = 10²处评估,要求存在不稳定螺旋。多平衡态以及缺乏特征值复数分量的单稳性均被惩罚。我们还优先选择特征值λ = a ± iω,使排斥更快、振荡频率更高,以筛选生物学上可观测的振荡:

d_cycle = { 3|1 - N_SS|,若N_SS ≠ 1;2.5,若N_SS = 1且λ ∈ ℝ;l(-a) + 2l(-|ω|),若N_SS = 1、λ ∈ ℂ且|ω| > 0 }

其中N_SS为稳态数,l(x) = 1/(1 + e^(-x))为逻辑函数用于归一化距离度量。类似地,S = 1和S = 10⁴处的单稳行为通过惩罚多稳性、偏好非振荡快速节点来选择:

d_mono = { 3|1 - N_SS|,若N_SS ≠ 1;l(max(a)) + 2l(ω) - 1,若N_SS = 1 }

我们还在S = 10³处评估双稳性,选择三个性质未定义稳态的共存:d_bistable = 3|3 - N_SS|。

对于单诱导物与双诱导物分析,我们生成4×10⁴个参数集(每种回路),满足最终距离函数阈值d = 1.35。对于分岔动物园结果(图5),目标不是估计后验分布,而是揭示特定分岔图的存在性。为此,我们采用半自动方法扩展ABC搜索:首先通过修改距离函数放宽初始ABC筛选以允许其他多功能组合,然后在特定参数集周围进行定向手动探索,以锐化观察到的行为或揭示邻近替代方案。例如,改变降解率可诱导与振荡幅度相关的动力学条形码之间的转换;调节信号-聚合酶协同性可显著改变多稳性范围;其他参数变化可微调单个分支并移动多功能性出现的信号范围。导致SS→SS←、SH→SH←、LS→HS←、LH→HH←、LH→JH←、LS→JS←、∅→JH←行为的初始参数集通过以下反转分岔序列的逆距离函数获得:d = d_mono(1) + d_bistable(10) + d_cycle(10²) + d_mono(10⁴)。对于初始HH→HH←、HH→HJ←、HL→HJ←和HS→HS←,我们使用较小输入信号值,距离函数为d = d_mono(10⁻⁶) + d_cycle(10⁻⁴) + d_bistable(10⁻³) + d_mono(10⁻²)。该搜索分类不包括与AC-DC拓扑兼容的平凡情形,如无多稳性的单稳性、振荡。

【数据】参数集数量:4×10⁴/回路;距离阈值:d=1.35;评估信号值:S=1、10²、10³、10⁴;逆距离函数信号值:10、10²、10⁴及10⁻⁶、10⁻⁴、10⁻³、10⁻²;结合位点数:n=3

软件与代码

使用NumPy(Harris等,https://numpy.org)和SciPy(Virtanen等,https://scipy.org)进行模拟与分析。模拟与分析代码见https://zenodo.org/records/19555341。

【数据】软件:NumPy、SciPy;代码仓库:zenodo.org/records/19555341;原文未详述版本号

研究结果

生物物理探索

AC-DC回路的动力学可用常微分方程(ODE)描述,表征三个节点X、Y、Z在外部诱导信号S影响下的表达水平:

(Equation 1) dX/dt = α_X·φ_X(Z; S, θ_X) - δ_X·X,dY/dt = α_Y·φ_Y(X; S, θ_Y) - δ_Y·Y,dZ/dt = α_Z·φ_Z(X, Y; S, θ_Z) - δ_Z·Z

其中参数α和δ分别对应每个转录因子的最大生产率和降解率(见图1B)。调控逻辑编码于调控函数φ(·)的函数形式与参数向量θ中。以往模型使用Hill函数或广义线性模型,而我们采用基于枚举所有可能启动子占据状态及其结合能的热力学描述(见图1B及补充信息),从而计算基因处于活跃转录状态(结合RNA聚合酶)的概率。在此形式体系中,θ_j中每个参数都有生物物理解释,非线性自然产生于预设的分子相互作用,无需强加Hill指数。具体而言,每个基因由分子i对启动子j的结合亲和力集合K_ji表征;例如K_XP表示聚合酶对基因X空启动子的结合亲和力,K_ZY表示蛋白Y对启动子Z的结合亲和力。在该表述中,诱导信号S的作用是通过协同性参数c_i降低RNA聚合酶的结合能。热力学形式体系不仅提供更高的机制真实性,还涵盖更广泛的输入-输出行为类别,Hill函数仅构成其特定极限情形(如Estrada等和Martinez-Corral等)。正文中的调控函数对应于三个阻遏物结合位点且不与信号结合相互作用的特定情形。RNA聚合酶仅在所有阻遏物结合位点未占据时才能结合启动子,而外部信号促进聚合酶招募。函数的显式形式见STAR Methods,其中双诱导物回路包含13个不同参数(无量纲化前17个)。此外,我们分析了纳入不同生物物理假设的替代热力学表述,包括阻遏物结合位点数目变化及信号-阻遏物结合竞争。这允许我们探索AC-DC启动子更广泛的生物物理参数空间。为清晰起见,这些替代架构的大部分结果呈现在补充信息中,并在出现定性差异时于正文引用。

AC-DC是一种多功能基因调控基序。(A)AC-DC电路的拓扑结构。(B)基因表达采用热力学形式体系建模,该体系枚举所有可能的启动子结合状态及其相应的平衡能量。(C)在某一参数区间中的示例轨迹,其中
▲ AC-DC是一种多功能基因调控基序。(A)AC-DC电路的拓扑结构。(B)基因表达采用热力学形式体系建模,该体系枚举所有可能的启动子结合状态及其相应的平衡能量。(C)在某一参数区间中的示例轨迹,其中

模型的初步探索证实其能够重现振荡与稳态之间的双稳性,与先前报道的动力学行为一致(见图1C和1D)。然而,要揭示该表述可及的全部行为范围,需要更系统地探索其高维参数空间。为此,我们采用基于序列蒙特卡洛方案的贝叶斯推断策略,逐步雕琢出与预设动力学行为兼容的参数空间区域。具体而言,我们定义评分函数量化给定参数集与目标多功能行为的距离。例如,在本工作第一部分,多功能AC-DC行为通过不同诱导信号S值下的一组动力学区域定义——即单稳、双稳和振荡区域的组合(见图1E)。相关ABC评分不仅由离散区域分类定义,还纳入连续动力学特征(如特征值谱和振荡质量),以偏好具有生物学意义的吸引子特性。评分函数的显式形式见STAR Methods。该方法的一个关键优势是不试图找到满足动力学的单一参数集点估计,而是识别支持预设多功能性的参数空间区域。这使我们能够评估每种行为的鲁棒性(此处定义为支持给定动力学区域的参数空间相对体积)、促进模型选择,并最终提供具体设计指南。此外,该方法允许通过采用涵盖实验测量范围(跨数量级)的结合亲和力与诱导物浓度的先验概率分布,将探索限制在可行的生物物理参数范围内。

【数据】节点数:3(X、Y、Z);结合位点数:n=3;双诱导物参数:13个(无量纲化后)/17个(原始);单诱导物参数:11个(无量纲化后)/15个(原始);评估信号值:S=1、10²、10³、10⁴;原文未详述具体α、δ、K数值

单诱导物与双诱导物

文献中AC-DC动力学采用双诱导物模型捕获,信号激活作用于回路中两个基因(见图2A),灵感来自指导神经祖细胞分化的原始网络拓扑。问题在于是否存在更简单的系统能展现相同的多功能动力学。为此,我们考虑单诱导物系统(见图2A),信号激活仅作用于基因X。这一简化使我们能够确定双输入调控是否为AC-DC多功能性所必需,还是该行为源于核心AC-DC拓扑,同时为合成实现提供更简单的设计规则。

Comparison of multifunctional parameter spaces for one- and two-inducer circuits (A) One- and two-inducer AC-DC circuits. (B) Posterior probability distributions comparing one-inducer (1i, blue) vs. two-inducer systems (2i, pink), including marginal posterior distribution for the inhibitor binding K
▲ Comparison of multifunctional parameter spaces for one- and two-inducer circuits (A) One- and two-inducer AC-DC circuits. (B) Posterior probability distributions comparing one-inducer (1i, blue) vs. two-inducer systems (2i, pink), including marginal posterior distribution for the inhibitor binding K

利用多功能贝叶斯探索(见STAR Methods和图1E),我们发现单诱导物和双诱导物系统均展现多功能动力学。我们预期双诱导物系统在更广泛参数空间区域展现多功能性,因为单诱导物系统是双诱导物系统的子情形。检查两个模型中共同参数的边缘和成对后验分布证实了这一点(见图2B和S5)。然而,双诱导物模型包含两个额外参数K_YS和c_Y,使其参数空间维度更高。因此,各自超体积共享维度上更广泛的分布本身并不能揭示哪个回路更易出现多功能行为——即在随机采样参数集时哪个回路更可能出现多功能性。为此,我们比较了后验分布相对于先验的相对超体积(类似于贝叶斯因子;详见STAR Methods)。该分析显示,单诱导物模型展现多功能性的可能性是双诱导物模型的1.5倍。这意味着一旦考虑维度,单诱导物模型的有效分布略大于双诱导物模型。成对后验参数相关性进一步支持这一结论(图2C),其中双诱导物模型显示更强的相关性。高相关性(如K_YX与K_XZ之间)表明对参数扰动的脆弱性,需要协调调节多个参数才能使双诱导物回路适用于合成开发。

鉴于与多功能性兼容的广泛参数范围,我们探索参数空间以揭示与该行为相关的结构特征。分析显示启动子强度(聚合酶结合亲和力K_XP、K_YP、K_ZP)的相对排序起关键作用,仅观察到三种主导排序,单诱导物系统中仅两种,对应于K_XP为最弱启动子(见图2D)。这一观察在替代阻遏物结合架构中得以保留(图S8),表明其不依赖于[原文此处截断]。

【数据】单诱导物vs双诱导物多功能性概率比:1.5倍;参数集数量:4×10⁴/回路;距离阈值:d=1.35;启动子强度排序:双诱导物3种主导排序/单诱导物2种;最弱启动子:K_XP;原文未详述具体亲和力数值

讨论与解读

本研究采用比以往更具机制真实性的生物物理表述,探索了AC-DC回路丰富的动力学库。该回路不仅保留了先前报道的多功能性,而且这些行为出现在仅需单一信号输入的简化版本中,无需依赖多诱导物或嵌入更大网络。这对理解AC-DC类基序在自然发育回路中的作用以及合成生物学——从最小组件构建鲁棒多功能回路是核心设计目标——均具有重要意义。ABC的使用为理性回路设计提供了实用途径,返回支持目标动力学行为的结合亲和力与启动子强度的明确范围和排序,使回路更易于用标准合成生物学工具包实现,促进模型驱动推断与可测生物物理参数相结合的设计循环,并解决系统生物学中ODE预测常见的模型不可辨识性问题。

结果还强调调控函数的选择可强烈影响机制预测。例如,先前研究将鞍环分岔确定为连接振荡[原文此处截断]的关键路径,而本工作揭示的多样性表明拓扑约束远比经典观点宽松。

编译者解读

该研究的方法论价值在于用热力学形式体系替代Hill函数,使每个参数具备可测量的生物物理含义,并借ABC将"设计"转化为参数空间体积的定量比较——单诱导物比双诱导物多功能性概率高1.5倍这一结论直接指导合成实现应优先选择更简回路。局限在于全部结论基于确定性ODE与数值分岔分析,未纳入噪声、细胞生长与资源竞争;三十余种分岔图的存在性依赖半自动手动搜索,其生物学可及性仍需实验验证。对合成生物学而言,这提供了一份"动力学条形码"设计目录,但落地仍需与可测启动子强度库对接。

参考来源

期刊:Cell Systems, 2026

DOI: 10.1016/j.cels.2026.101715

DOI: 10.1016/j.cels.2026.101715

延伸阅读:更多代谢通路设计、多基因组装等合成生物学工具,可访问 DNA Lab Space。

📎 相关文章

无抗生素标记蓝细菌底盘的一体化工程 11-05 染色体工程改造酪丁酸梭菌攻克碳分解代谢阻遏 11-05 光遗传学UVA生物传感器平台:高通量防晒阻断剂测试新方案 11-05 敲除限制系统打造转化底盘 11-05