数学建模 · GENEBELLE

一、建模概述与系统描述

1.1 生物学背景

人体腋臭的主要来源是腋下大汗腺分泌的臭味前体分子（如 Cys-Gly-3M3SH、HMHA-Gln 等谷胱甘肽/半胱氨酸-S-缀合物），这些前体分子本身无臭，但经腋下共生菌（主要为棒状杆菌属 Corynebacterium）的 C-S β-裂解酶切割后，释放出挥发性硫醇（如 3M3SH）和短链脂肪酸（如 3-HMHA），产生特征性腋臭。

这些前体分子由肝脏合成后经血液循环输送至腋下大汗腺。因此，如果能在前体分子进入血液循环之前，在肠道内将其截留，即可从源头阻断腋臭的产生。

1.2 工程菌设计策略

本项目设计的内服工程大肠杆菌，其核心除臭模块为 Lpp-OmpA-Linker-ApoD 表面展示系统：

核心逻辑：将人类载脂蛋白D（ApoD）通过 Lpp-OmpA 表面展示系统锚定于工程菌外膜表面。ApoD 作为脂质运载蛋白（lipocalin）家族成员，具有疏水口袋结构，可高亲和力结合臭味前体分子。被捕获的前体分子随工程菌一起通过粪便排出体外。

调控策略：该模块受磷酸化 RpaA 激活的 pKaiBC 启动子调控，实现昼夜节律性表达（模拟白天高表达、夜间低表达），与人体活动周期匹配。

1.3 建模目标

建立从基因表达到全身药代动力学的五层级联动力学模型，量化预测工程菌的除臭效率，并指导关键参数的优化设计。模型结构如下：

$$\underbrace{\text{pKaiBC 启动子}}_{\text{模块一}} \xrightarrow{\text{转录翻译}} \underbrace{\text{表面展示}}_{\text{模块二}} \xrightarrow{\text{结合}} \underbrace{\text{ApoD-前体复合物}}_{\text{模块三}} \xrightarrow{\text{清除}} \underbrace{\text{肠道药代}}_{\text{模块四}} \xrightarrow{\text{效率}} \underbrace{\text{全局模型}}_{\text{模块五}}$$

二、模块一：pKaiBC 启动子昼夜节律表达动力学

2.1 模型假设

pKaiBC 启动子受 KaiABC 振荡器磷酸化的 RpaA（RpaA~P）激活。在工程菌中重建 Kai 振荡器后，RpaA~P 浓度呈现近似正弦的 24 小时振荡周期。我们以 RpaA~P 的振荡驱动 mRNA 转录速率。

2.2 RpaA~P 振荡函数

磷酸化 RpaA 浓度随时间变化近似为：

$$[RpaA{\sim}P](t) = \bar{R} + A_R \cdot \cos\!\left(\frac{2\pi}{T}(t - \phi)\right)$$ (1)

其中：

符号	含义	取值
$\bar{R}$	RpaA~P 平均浓度	0.5 μM
$A_R$	振荡振幅	0.4 μM
$T$	振荡周期	24 h
$\phi$	相位偏移（峰值对应白天正午）	6 h

2.3 mRNA 转录动力学

Lpp-OmpA-ApoD 融合蛋白的 mRNA 转录速率受 RpaA~P 以 Hill 函数形式激活：

$$\frac{d[mRNA]}{dt} = \alpha_m \cdot \frac{[RpaA{\sim}P]^n}{K_d^n + [RpaA{\sim}P]^n} + \alpha_{\text{basal}} - \delta_m \cdot [mRNA]$$ (2)

符号	含义	取值
$\alpha_m$	最大转录速率	50 nM/h
$\alpha_{\text{basal}}$	基础泄漏转录速率	2 nM/h
$K_d$	RpaA~P 半激活浓度	0.3 μM
$n$	Hill 系数	2
$\delta_m$	mRNA 降解速率常数	6.93 h⁻¹（半衰期 ~6 min）

2.4 蛋白翻译动力学

融合蛋白总量（包括胞质内未折叠、折叠中间态和已插入外膜的蛋白）的变化为：

$$\frac{d[P_{\text{total}}]}{dt} = \alpha_p \cdot [mRNA] - \delta_p \cdot [P_{\text{total}}]$$ (3)

符号	含义	取值
$\alpha_p$	翻译速率常数	20 h⁻¹
$\delta_p$	蛋白总降解/稀释速率	0.35 h⁻¹（细胞分裂稀释为主）

三、模块二：Lpp-OmpA-ApoD 表面展示动力学

3.1 表面展示过程建模

翻译后的融合蛋白需经过信号肽切割、跨内膜转运、OmpA 跨膜区折叠插入外膜等步骤，才能成功展示于细胞表面。膜插入速率常数约 0.5 h⁻¹，单细胞最大表面展示容量约 10⁵ molecules/cell，表面蛋白脱落/失活速率约 0.1 h⁻¹。外膜面积有限，当表面展示分子数接近上限时，新蛋白无法继续插入，产生自抑制效应。动力学方程为：

$$\frac{d[S]}{dt} = k_{\text{ins}} \cdot [P_{\text{total}}] \cdot \left(1 - \frac{[S]}{S_{\max}}\right) - \delta_s \cdot [S]$$ (4)

其中 $[S]$ 为单个细菌表面展示的活性 ApoD 分子数（molecules/cell）。

符号	含义	取值
$k_{\text{ins}}$	膜插入速率常数	0.5 h⁻¹
$S_{\max}$	单细胞最大表面展示容量	~10⁵ molecules/cell
$\delta_s$	表面蛋白脱落/失活速率	0.1 h⁻¹

饱和项 $(1 - S/S_{\max})$ 的物理意义：外膜面积有限，当表面展示分子数接近上限时，新蛋白无法继续插入，产生自抑制效应。

3.2 群体水平总有效结合位点

考虑肠道中工程菌的总菌量，群体水平的总活性 ApoD 结合位点浓度为：

$$[S_{\text{total}}](t) = [S](t) \times N_{\text{bac}}(t) \times \frac{1}{V_{\text{gut}} \cdot N_A}$$ (5)

其中 $N_{\text{bac}}(t)$ 为肠道内工程菌总数，$V_{\text{gut}}$ 为肠道有效容积，$N_A$ 为阿伏伽德罗常数。建模中假设工程菌已达到稳态定殖水平：

符号	含义	取值
$N_{\text{bac}}$	肠道内工程菌稳态总数	10⁹–10¹⁰ CFU
$V_{\text{gut}}$	肠道有效液体容积	~0.5 L

四、模块三：ApoD 与臭味前体分子的结合动力学

4.1 结合反应方程

表面展示的 ApoD（$S$）与肠道内游离臭味前体分子（$L$，如 Cys-Gly-3M3SH）可逆结合形成表面复合物（$C$）：

$$S + L \underset{k_{\text{off}}}{\overset{k_{\text{on}}}{\rightleftharpoons}} C$$

4.2 结合动力学方程组

基于质量作用定律：

$$\frac{d[C]}{dt} = k_{\text{on}} \cdot [S_{\text{free}}] \cdot [L] - k_{\text{off}} \cdot [C]$$ (6)

$$[S_{\text{free}}] = [S_{\text{total}}] - [C]$$ (7)

其中解离常数为：

$$K_D = \frac{k_{\text{off}}}{k_{\text{on}}}$$ (8)

符号	含义	取值
$k_{\text{on}}$	结合速率常数	10⁶ M⁻¹s⁻¹（典型脂质运载蛋白）
$k_{\text{off}}$	解离速率常数	10⁻³ s⁻¹
$K_D$	平衡解离常数	~1 μM（ApoD 对疏水小分子典型值）

4.3 准稳态近似

由于结合-解离反应速率（秒级）远快于蛋白表达和肠道转运速率（小时级），可对结合反应施加准稳态近似（QSSA）：

$$[C]_{\text{ss}} = \frac{[S_{\text{total}}] \cdot [L]}{K_D + [L]}$$ (9)

这是经典的 Langmuir 等温吸附形式：当 $[L] \gg K_D$ 时结合位点趋于饱和；当 $[L] \ll K_D$ 时捕获率近似线性正比于 $[L]$。

4.4 有效捕获通量

单位时间内被工程菌表面捕获并固定化（随粪便清除）的前体分子通量为：

$$J_{\text{capture}}(t) = k_{\text{off}} \cdot [C]_{\text{ss}} \cdot f_{\text{irr}} + k_{\text{elim}} \cdot [C]_{\text{ss}}$$ (10)

简化后，有效清除速率可表示为：

$$J_{\text{capture}}(t) = k_{\text{eff}} \cdot \frac{[S_{\text{total}}](t) \cdot [L](t)}{K_D + [L](t)}$$ (11)

其中 $k_{\text{eff}}$ 为有效捕获速率常数（包括复合物随菌体排出和不可逆结合的综合效应），取 $k_{\text{eff}} \approx 0.2 \text{ h}^{-1}$。

五、模块四：肠道腔内前体分子药代动力学

5.1 前体分子来源与去向

流入项：① 饮食来源的含硫氨基酸经肝脏代谢后以前体分子形式随胆汁分泌入肠道；② 肠道上皮细胞分泌。

流出项：① 肠壁重吸收进入血液循环（主要经门静脉）；② 工程菌表面 ApoD 捕获清除；③ 随粪便自然排出。

5.2 肠道前体分子动力学方程

$$\frac{d[L]}{dt} = \underbrace{R_{\text{in}}(t)}_{\text{流入}} - \underbrace{k_a \cdot [L]}_{\text{肠壁吸收入血}} - \underbrace{k_{\text{eff}} \cdot \frac{[S_{\text{total}}](t) \cdot [L]}{K_D + [L]}}_{\text{ApoD 捕获}} - \underbrace{k_{\text{fecal}} \cdot [L]}_{\text{粪便排出}}$$ (12)

符号	含义	取值
$R_{\text{in}}(t)$	前体分子流入速率（含餐后波动）	见下文
$k_a$	肠壁吸收速率常数	0.3 h⁻¹
$k_{\text{fecal}}$	粪便清除速率常数	0.04 h⁻¹（肠道转运时间 ~24h）

5.3 前体分子流入函数

假设每日三餐后肝脏代谢产生前体分子，可建模为脉冲式输入叠加基础分泌：

$$R_{\text{in}}(t) = R_0 + \sum_{i=1}^{3} \frac{D_i}{\sigma\sqrt{2\pi}} \exp\!\left(-\frac{((t \bmod 24) - t_i)^2}{2\sigma^2}\right)$$ (13)

符号	含义	取值
$R_0$	基础分泌速率	0.5 μM/h
$D_i$	第 i 餐产生的前体分子总量	5 μM（早）/ 8 μM（午）/ 10 μM（晚）
$t_i$	第 i 餐后前体到达肠道的时间	9 h / 14 h / 20 h
$\sigma$	高斯展宽（胃肠转运延迟）	1.5 h

本图对比了工程菌处理组与未处理对照组的肠道内臭味前体分子浓度动态，对应肠道药代动力学，灰色阴影为三餐带来的前体脉冲式输入。对照组中，前体浓度随三餐输入出现显著峰值，最高接近 8μM；而工程菌处理组中，ApoD 对前体的高效捕获使肠道内游离前体浓度大幅降低，峰值被显著抑制，全程维持在 2μM 以下，直接减少了可被肠壁吸收进入血液循环的前体总量。

六、模块五：血液-腋下转运与除臭效率全局模型

6.1 血浆前体分子动力学

被肠壁吸收进入血液的前体分子浓度变化为：

$$\frac{d[L_{\text{blood}}]}{dt} = \frac{k_a \cdot [L] \cdot V_{\text{gut}}}{V_{\text{blood}}} - k_{\text{hepatic}} \cdot [L_{\text{blood}}] - k_{\text{axilla}} \cdot [L_{\text{blood}}] - k_{\text{renal}} \cdot [L_{\text{blood}}]$$ (14)

符号	含义	取值
$V_{\text{blood}}$	血浆分布容积	~3 L
$k_{\text{hepatic}}$	肝脏首过代谢清除速率	0.5 h⁻¹
$k_{\text{axilla}}$	腋下大汗腺摄取速率	0.05 h⁻¹
$k_{\text{renal}}$	肾脏清除速率	0.1 h⁻¹

6.2 腋下前体分子蓄积

到达腋下大汗腺区域并被分泌至皮肤表面的前体分子量：

$$\frac{d[L_{\text{axilla}}]}{dt} = k_{\text{axilla}} \cdot [L_{\text{blood}}] \cdot \frac{V_{\text{blood}}}{V_{\text{axilla}}} - k_{\text{secrete}} \cdot [L_{\text{axilla}}]$$ (15)

6.3 臭味强度评估函数

最终腋臭强度与腋下前体分子被细菌裂解产生挥发性硫醇的量成正比：

$$\text{Odor}(t) = \eta \cdot k_{\text{secrete}} \cdot [L_{\text{axilla}}](t)$$ (16)

其中 $\eta$ 为腋下微生物的裂解转化效率。

6.4 除臭效率定义

除臭效率定义为 24 小时内臭味强度的相对降低百分比：

$$\boxed{E_{\text{deodor}} = \left(1 - \frac{\displaystyle\int_0^{24} \text{Odor}_{\text{treated}}(t)\,dt}{\displaystyle\int_0^{24} \text{Odor}_{\text{untreated}}(t)\,dt}\right) \times 100\%}$$ (17)

其中 "treated" 指服用工程菌后的情况，"untreated" 为对照组（无工程菌，$[S_{\text{total}}] = 0$）。

全局模型 ODE 系统总结

核心方程组（共6个耦合常微分方程）：

(1) $\frac{d[mRNA]}{dt}$：式(2)，由 RpaA~P 振荡驱动
(2) $\frac{d[P_{\text{total}}]}{dt}$：式(3)，由 mRNA 驱动翻译
(3) $\frac{d[S]}{dt}$：式(4)，膜插入与饱和展示
(4) $\frac{d[L]}{dt}$：式(12)，肠道前体分子核心方程
(5) $\frac{d[L_{\text{blood}}]}{dt}$：式(14)，血浆动力学
(6) $\frac{d[L_{\text{axilla}}]}{dt}$：式(15)，腋下蓄积

输出指标：式(17)，24h 累积除臭效率 $E_{\text{deodor}}$。

本图是全链路模型的核心输出结果，对应腋下前体蓄积与臭味强度评估。对照组（紫色填充）的腋臭强度随三餐和昼夜节律出现显著峰值，日间活动时段强度最高达 14 a.u.；工程菌处理组（绿色填充）的臭味强度全程被大幅压制，峰值降至 10 a.u. 以下，24 小时累计臭味强度显著降低，直观验证了工程菌的核心除臭效果。

七、参数汇总与取值依据

参数	含义	取值	来源/依据
$\bar{R}$	RpaA~P 平均浓度	0.5 μM	Kai 振荡器文献
$A_R$	振荡振幅	0.4 μM	Kai 振荡器文献
$T$	振荡周期	24 h	设计需求
$\alpha_m$	最大转录速率	50 nM/h	大肠杆菌强启动子文献
$\delta_m$	mRNA 降解速率	6.93 h⁻¹	大肠杆菌 mRNA 半衰期 ~6 min
$\alpha_p$	翻译速率常数	20 h⁻¹	BioNumbers 典型值
$\delta_p$	蛋白降解/稀释速率	0.35 h⁻¹	肠道内 ~2h 倍增时间
$S_{\max}$	单细胞最大展示量	10⁵	Lpp-OmpA 系统文献
$k_{\text{on}}$	ApoD 结合速率常数	10⁶ M⁻¹s⁻¹	脂质运载蛋白典型值
$K_D$	平衡解离常数	1 μM	ApoD 对疏水配体文献估计
$N_{\text{bac}}$	工程菌定殖量	10⁹ CFU	益生菌口服定殖文献
$k_a$	肠壁吸收速率	0.3 h⁻¹	小分子肠道吸收文献
$k_{\text{fecal}}$	粪便清除速率	0.04 h⁻¹	肠道转运时间 ~24h
$k_{\text{hepatic}}$	肝脏清除速率	0.5 h⁻¹	肝脏首过代谢文献

八、灵敏度分析

8.1 局部灵敏度分析

对除臭效率 $E_{\text{deodor}}$ 关于各参数 $p_i$ 进行归一化局部灵敏度分析：

$$S_i = \frac{\partial E_{\text{deodor}}}{\partial p_i} \cdot \frac{p_i}{E_{\text{deodor}}}$$ (18)

8.2 关键灵敏度参数预测

基于模型结构分析，以下参数对除臭效率的影响最为敏感：

参数	灵敏度方向	生物学含义
$N_{\text{bac}}$（工程菌定殖量）	正相关	菌量越多，总结合位点越多，捕获效率越高
$K_D$（解离常数）	负相关	亲和力越高（$K_D$ 越小），捕获效率越高
$S_{\max}$（单细胞展示量）	正相关	表面展示密度直接决定捕获能力
$k_a$（肠壁吸收速率）	负相关	吸收越快，ApoD 来不及捕获的前体越多

设计优化建议：模型预测提高除臭效率的最有效策略依次为：① 提高工程菌定殖量（增加口服剂量或添加定殖辅助模块）；② 通过定向进化降低 ApoD 对目标底物的 $K_D$；③ 优化 Lpp-OmpA linker 长度以提升 $S_{\max}$。

8.3 全局灵敏度分析

推荐使用 Sobol' 指数法进行全局灵敏度分析，通过 Monte Carlo 采样在参数空间内计算一阶和总效应灵敏度指数：

$$S_i^{\text{total}} = 1 - \frac{\text{Var}_{X_{\sim i}}(E[E_{\text{deodor}} | X_{\sim i}])}{\text{Var}(E_{\text{deodor}})}$$ (19)

其中 $X_{\sim i}$ 表示除 $p_i$ 外的所有参数。

本图为除臭效率对模型核心参数的归一化局部灵敏度分析结果，对应模型第 8 部分的灵敏度分析公式，量化了各参数对最终除臭效率的影响程度与作用方向，是工程菌设计优化优先级的核心决策依据。

1、第一优先级核心优化靶点：单细胞最大展示容量 $S_{\max}$、工程菌肠道定殖量 $N_{\text{bac}}$、有效捕获速率常数 $k_{\text{eff}}$，三者的归一化灵敏度指数均为 +0.868，是对除臭效率影响最大的核心参数。正灵敏度说明提升这三个参数可显著提高除臭效率，且三者影响权重完全一致，是工程菌优化的第一优先级靶点，与模型预测的核心优化策略完全匹配。

2、第二优先级优化靶点：肠壁吸收速率常数 $k_a$（灵敏度 -0.239）、ApoD 与底物的解离常数 $K_D$（灵敏度 -0.230）为核心负相关参数。负灵敏度说明降低这两个参数可显著提升除臭效率：$K_D$ 越小代表 ApoD 对底物的亲和力越高，捕获能力越强；$k_a$ 越小代表前体分子在肠道内的停留时间越长，ApoD 可捕获的窗口期越久，二者是工程菌优化的第二优先级靶点。

3、非限速低影响参数：膜插入速率 $k_{\text{ins}}$、最大转录速率 $\alpha_m$，灵敏度分别为 +0.068、+0.062，对除臭效率的影响极弱，说明转录与膜插入步骤并非整个除臭系统的限速步骤，无需作为核心优化方向。

4、可忽略参数：粪便清除速率 $k_{\text{fecal}}$，灵敏度仅 -0.032，对除臭效率的影响几乎可以忽略，说明粪便自然清除过程不会成为系统的性能瓶颈，无需针对性优化。

本图为除臭效率的剂量-响应曲面，量化了工程菌肠道定殖量（横轴）与 ApoD 底物亲和力（$K_D$，三条曲线）两个核心参数对除臭效率的协同影响，为工程菌的口服剂量设计、蛋白亲和力定向进化提供了直接的量化参考。

1、定殖量的剂量阈值效应：在所有 $K_D$ 条件下，除臭效率均随工程菌定殖量的提升呈 S 型上升，存在明确的剂量阈值。当定殖量低于 10⁹ CFU 时，无论亲和力高低，除臭效率均低于 10%，几乎无实用除臭效果；当定殖量达到 10¹⁰~10¹¹ CFU 区间时，除臭效率进入快速上升期；定殖量超过 10¹¹ CFU 后，效率逐渐进入平台期，趋近于 100%。

2、亲和力的增效与阈值降低效应：ApoD 的亲和力提升（$K_D$ 降低）可显著降低除臭效率的剂量阈值，大幅提升低定殖量下的工程菌性能。$K_D=0.1$ μM（高亲和力，绿色曲线）时，定殖量 10¹⁰ CFU 即可实现近 70% 的除臭效率，10¹¹ CFU 即可接近 100%；而 $K_D=10$ μM（低亲和力，红色曲线）时，即使定殖量达到 10¹² CFU，除臭效率也仅 95% 左右，10¹⁰ CFU 时效率不足 20%，直观验证了 ApoD 亲和力优化对工程菌性能的关键增效作用。

3、定殖量与亲和力的协同优化效应：定殖量与亲和力存在显著的协同作用，同步提升定殖量与优化亲和力，可在更低的口服剂量下实现超高除臭效率。例如 $K_D=0.1$ μM + 10¹¹ CFU 的组合，即可实现接近 100% 的除臭效率，远优于单一参数优化的效果，为工程菌的最优设计方案提供了量化支撑。

九、模型预测与结论

9.1 定性预测

（1）昼夜节律捕获模式：由于 pKaiBC 启动子在白天（RpaA~P 高水平期）驱动 ApoD 高表达，表面展示的 ApoD 数量在白天达到峰值。这与午餐和晚餐后的前体分子流入高峰时间匹配，实现了「按需捕获」。

（2）剂量-效应关系：在前体分子浓度 $[L] \ll K_D$ 的生理条件下，捕获速率近似为一级动力学 $J_{\text{capture}} \approx k_{\text{eff}} \cdot [S_{\text{total}}] \cdot [L] / K_D$，除臭效率与工程菌数量成近似线性关系。

（3）稳态除臭效率估算：在标准参数下（$N_{\text{bac}} = 10^9$，$K_D = 1$ μM），模型预测工程菌可截留约 30–50% 的肠道前体分子，使腋下臭味前体供应减少约 25–40%，可实现有临床意义的腋臭减轻。

9.2 定量估算：稳态下的捕获分数

在准稳态条件下（$[L] \ll K_D$），定义等效一级捕获速率常数：

$$k_{\text{cap}}^* = k_{\text{eff}} \cdot \frac{[S_{\text{total}}]}{K_D}$$ (20)

则肠道前体分子被工程菌捕获的分数为：

$$\boxed{f_{\text{capture}} = \frac{k_{\text{cap}}^*}{k_a + k_{\text{cap}}^* + k_{\text{fecal}}}}$$ (21)

数值代入估算：

取 $N_{\text{bac}} = 10^9$，$S = 5 \times 10^4$ molecules/cell（展示量均值），则：

$$[S_{\text{total}}] = \frac{10^9 \times 5 \times 10^4}{0.5 \times 6.02 \times 10^{23}} \approx 0.17 \;\mu\text{M}$$

$$k_{\text{cap}}^* = 0.2 \times \frac{0.17}{1} = 0.034 \;\text{h}^{-1}$$

$$f_{\text{capture}} = \frac{0.034}{0.3 + 0.034 + 0.04} = \frac{0.034}{0.374} \approx 9.1\%$$

若将定殖量提升至 $10^{10}$ CFU（高剂量益生菌方案），则 $f_{\text{capture}} \approx 47\%$，达到显著治疗效果。

若同时通过定向进化将 $K_D$ 从 1 μM 降至 0.1 μM，则在 $10^{10}$ 定殖量下：

$$k_{\text{cap}}^* = 0.2 \times \frac{1.7}{0.1} = 3.4 \;\text{h}^{-1}$$

$$f_{\text{capture}} = \frac{3.4}{0.3 + 3.4 + 0.04} \approx 91\%$$

模型核心结论

（1）本数学模型建立了从基因调控到全身药代动力学的完整五层级联 ODE 系统，完整描述了 Lpp-OmpA-ApoD 工程菌在肠道内捕获臭味前体分子的全过程。

（2）模型预测工程菌的除臭效率主要取决于三个关键参数：工程菌定殖量 $N_{\text{bac}}$、ApoD 亲和力 $K_D$ 和表面展示密度 $S_{\max}$。

（3）在标准参数条件下，单次口服剂量可截留约 9% 的前体分子；通过提高定殖量（10¹⁰ CFU）和优化 ApoD 亲和力（$K_D = 0.1$ μM），理论捕获率可达 90% 以上。

（4）pKaiBC 昼夜节律调控实现了与人体代谢高峰的时间匹配，最大化了 ApoD 的捕获效率。

本柱状图量化了不同工程菌设计方案下的稳态前体捕获分数，为工程菌优化提供了明确的量化方向。低剂量组（$N=10^9$ CFU，$K_D=1$ μM）捕获效率仅 8.9%；标准设计组（$N=10^{10}$ CFU，$K_D=1$ μM）效率提升至 49.4%；经定向进化优化 ApoD 亲和力（$K_D=0.1$ μM）后，效率提升至 90.7%；同步提高定殖量与亲和力（High+Evolved，$N=10^{11}$ CFU，$K_D=0.1$ μM），最终捕获效率可达 99.0%。

十、安全模块——「与非门」自杀开关数学建模

1.1 系统描述

安全模块采用基于「与非门」(NAND gate)的自杀开关设计，包含以下核心元件：

元件	功能	响应条件
pfnrS（缺氧启动子）	感知低氧环境（肠道内）	缺氧时激活，表达 SupD
Y38（温敏启动子）	感知体温环境	37°C 时激活，表达 T7ptag
AND Gate（北京 2009）	SupD + T7ptag 共同存在时才输出下游蛋白	双信号同时存在
TetR	抑制 pTet 启动子	抑制 MazF 表达
MazF	毒素蛋白（序列特异性内切酶）	降解 mRNA，导致细菌死亡

1.2 逻辑真值表

缺氧信号 A	37°C 信号 B	A 和 B	TetR 表达	MazF 表达（自杀）	状态
0	0	0	无	有	死亡
0	1	0	无	有	死亡
1	0	0	无	有	死亡
1	1	1	有	无	存活

1.3 建模目标

量化描述工程菌在不同环境条件下的存活/死亡动力学，验证安全模块的可靠性和逃逸风险。

2. 模型变量定义

环境输入变量：O₂ 环境氧浓度（μM），肠道内约 0–10 μM，大气中约 200 μM；T 环境温度（℃），体内 37℃，体外通常 20–25℃。

分子浓度变量：[SupD] SupD 抑制性 tRNA 浓度；[T7ptag] T7ptag 浓度；[T7] 功能性 T7 RNA 聚合酶浓度（AND 门输出）；[TetR] TetR 阻遏蛋白浓度；[MazF] MazF 毒素蛋白浓度；N 工程菌细胞数（CFU）。

2.3 关键参数

参数	含义	典型值
$\alpha_{\text{SupD}}$	pfnrS 最大转录速率	50 nM/h
$K_{O_2}$	pfnrS 氧半抑制常数	5 μM
$\alpha_{\text{T7p}}$	Y38 最大转录速率	40 nM/h
$T_0$	Y38 激活阈值温度	37°C
$\alpha_{\text{TetR}}$	T7 驱动 TetR 最大转录速率	100 nM/h
$\alpha_{\text{MazF}}$	pTet 驱动 MazF 最大转录速率	80 nM/h
$K_{\text{TetR}}$	TetR 对 pTet 半抑制常数	15 nM
$\delta_p$	蛋白降解/稀释速率	0.3 h⁻¹
$\mu_{\max}$	最大生长速率	0.7 h⁻¹
$K_{\text{MazF}}$	MazF 致死半饱和常数	10 nM
$d_{\max}$	MazF 引起的最大致死速率	2.0 h⁻¹

十一、常微分方程 (ODE) 模型

3.1 环境感知层

pfnrS 缺氧启动子：氧浓度越低，表达越强（抑制型 Hill 函数）：

$$\frac{d[\text{SupD}]}{dt} = \alpha_{\text{SupD}} \cdot \frac{K_{O_2}^n}{K_{O_2}^n + [O_2]^n} - \delta_p \cdot [\text{SupD}]$$

Y38 温敏启动子：以 Sigmoid 函数模拟温度开关特性：

$$\frac{d[\text{T7ptag}]}{dt} = \alpha_{\text{T7p}} \cdot \frac{1}{1 + e^{(T-T_0)/\sigma_T}} - \delta_p \cdot [\text{T7ptag}]$$

3.2 AND 门逻辑层

AND 门核心机制：T7ptag 含有琥珀终止密码子，只有在 SupD（琥珀抑制 tRNA）存在时才能通读翻译为功能性 T7 RNA 聚合酶。用乘法项建模双信号协同：

$$\frac{d[T7]}{dt} = k_{\text{AND}} \cdot \frac{[\text{SupD}]}{K_S+[\text{SupD}]} \cdot \frac{[\text{T7ptag}]}{K_P+[\text{T7ptag}]} - \delta_p \cdot [T7]$$

典型参数：$k_{\text{AND}} = 60$ nM/h，$K_S = 10$ nM，$K_P = 10$ nM。

3.3 TetR 抑制层与 3.4 MazF 毒素层

功能性 T7 驱动 T7 promoter 下游的 TetR 表达。pTet 为组成型表达启动子，被 TetR 抑制；当 TetR 不足时，MazF 大量表达。MazF 表达式含 TetR 抑制项及基础泄漏 $\beta_{\text{MazF}}$（约 0.5–2 nM/h）。

3.5 菌群动力学

工程菌的净增长由正常生长和 MazF 致死效应共同决定：

$$\frac{dN}{dt} = \left[\mu_{\max} - d_{\max} \cdot \frac{[\text{MazF}]^{n_3}}{K_{\text{MazF}}^{n_3} + [\text{MazF}]^{n_3}}\right] \times N$$

当 $[\text{MazF}] \gg K_{\text{MazF}}$ 时，死亡速率趋近 $d_{\max}$，净增长率为负；当 $[\text{MazF}] \approx 0$ 时，菌群正常生长。

本图对应模型 3.2 节 AND 门逻辑层的核心设计，量化展示了自杀开关核心双输入信号（SupD、T7ptag）对下游功能性 T7 RNA 聚合酶输出的协同调控特性，验证了 AND 门的逻辑可靠性。

坐标轴定义：X 轴为琥珀抑制 tRNA SupD 的浓度，Y 轴为含琥珀终止密码子的 T7ptag 浓度，Z 轴与色阶为功能性 T7 聚合酶的稳态浓度，单位均为 nM。

核心逻辑验证：曲面呈现典型的「双输入协同激活」特性，仅当 SupD 和 T7ptag 两个输入同时处于高浓度区间时（对应肠道内缺氧 + 37℃ 的双信号激活场景），T7 功能酶的输出才会显著升高，进入平台期并达到 180 nM 左右的高稳态水平。

安全开关特性验证：当任意一个输入浓度极低时（对应体外单信号缺失/无信号场景），无论另一个输入浓度多高，T7 输出均维持在接近 0 的极低水平，完美契合 AND 门「全 1 出 1，有 0 出 0」的逻辑，从根源上避免了单信号泄漏导致的安全失效。

鲁棒性验证：曲面的梯度变化体现了 Michaelis-Menten 动力学的饱和特性，输入浓度超过半饱和常数后，输出进入稳定平台期，保证了肠道内双信号同时存在时，下游通路输出的稳定性。

十二、稳态分析

4.1 肠道内环境（缺氧+37°C）——设计存活态

令 $[O_2] \approx 0$，$T = 37°C$：$[\text{SupD}]_{\text{ss}} \approx 167$ nM（高表达），$[\text{T7ptag}]_{\text{ss}} \approx 67$ nM（高表达）。AND 门输出：$[T7]_{\text{ss}}$ 较高 → $[\text{TetR}]_{\text{ss}}$ 高 → MazF 被完全抑制 → 菌群正常生长。

4.2 体外有氧环境（有氧+25°C）——设计死亡态

令 $[O_2] = 200$ μM，$T = 25°C$：$[\text{SupD}]_{\text{ss}} \approx 0.1$ nM（几乎不表达），$[\text{T7ptag}]_{\text{ss}} \approx 0.045$ nM。AND 门输出：$[T7]_{\text{ss}} \approx 0$ → $[\text{TetR}]_{\text{ss}} \approx 0$ → MazF 高表达：$[\text{MazF}]_{\text{ss}} \approx 267$ nM。净增长率 $r_{\text{net}} \approx -1.3$ h⁻¹，菌群半衰期 $t_{1/2} \approx 0.53$ h ≈ 32 分钟。

4.3 单信号缺失场景

场景 1：有氧 + 37℃（如血液中）$[\text{SupD}] \approx 0$，$[\text{T7ptag}]$ 高 → AND 门输出 ≈ 0 → MazF 高表达 → 死亡√。

场景 2：缺氧 + 25℃（如厌氧培养箱低温）$[\text{SupD}]$ 高，$[\text{T7ptag}] \approx 0$ → AND 门输出 ≈ 0 → MazF 高表达 → 死亡√。

本图以环境氧浓度、环境温度两个核心输入变量为维度，绘制了系统的稳态相图，完整验证了自杀开关在不同环境下的功能有效性与鲁棒性。

菌群命运图：基于菌群净增长率绘制的生死相图，绿色区域净增长率 > 0，菌群可存活增殖；红色区域净增长率 < 0，菌群快速死亡；黑色曲线为生死边界（净增长率 = 0）。图中标注了两个核心设计场景：肠道内环境（低氧 + 37℃）处于绿色存活区，符合「肠道定殖存活」的设计需求；体外环境（有氧 200 μM + 25℃）处于红色死亡区，符合「逃逸自杀」的安全设计需求。

单信号缺失场景验证：相图同时验证了单信号缺失时的安全性——即使温度达到 37℃（如血液环境），只要氧浓度升高，或即使处于低氧环境，只要温度低于 35℃，菌群都会进入死亡区，避免了非肠道环境的异常定殖风险，验证了双信号校验设计的高鲁棒性。

十三、逃逸风险分析

5.1 突变逃逸概率模型

工程菌可通过突变逃逸自杀开关。主要逃逸途径及概率：MazF 基因失活突变 p₁ ≈ 10⁻⁶/代；pTet 启动子失活突变 p₂ ≈ 10⁻⁶/代；TetR 组成型突变 p₃ ≈ 10⁻⁷/代。单代逃逸概率 P_escape/代 ≈ p₁ + p₂ + p₃ ≈ 2.1×10⁻⁶。经过 G 代后至少出现一个逃逸突变体的概率：P_escape(G, N₀) = 1 − (1 − P_escape/代)^(G×N₀)。

5.2 逃逸菌的竞争适应度

即使突变逃逸，逃逸菌需在肠道内与野生型菌群竞争。适应度成本 s 表示突变造成的适应度损失（通常 0.01–0.1）。逃逸菌能稳定定殖的条件：μ_max(1−s) > w（w 为肠道洗出率约 0.1–0.3 h⁻¹）。

5.3 多重安全设计建议

为降低逃逸风险，建议增加冗余：

$$P_{\text{escape,redundant}} \approx (P_{\text{single}})^m$$

其中 $m$ 为独立自杀模块数。若 $m = 2$：$P_{\text{escape}}/\text{代} \approx (2.1\times10^{-6})^2 \approx 4.4\times10^{-12}$，高度安全。

十四、动态响应时间分析

6.1 从存活态到自杀态的切换

当工程菌从肠道（存活）排出体外（死亡）时，需分析自杀启动的延迟时间。TetR 衰减：

$$[\text{TetR}](t) = [\text{TetR}]_{\text{ss}} \cdot e^{-\delta_p \cdot t}$$

TetR 降到抑制失效水平（$[\text{TetR}] < K_{\text{TetR}}$）的时间：

$$t_{\text{switch}} = \frac{1}{\delta_p} \cdot \ln\!\left(\frac{[\text{TetR}]_{\text{ss}}}{K_{\text{TetR}}}\right)$$

若 $[\text{TetR}]_{\text{ss}} = 300$ nM，$K_{\text{TetR}} = 15$ nM：$t_{\text{switch}} \approx 10$ h。

本图完整模拟了工程菌从「肠道内存活态」到「体外环境死亡态」的全时域动态响应过程，直观呈现了自杀开关完整的时序调控逻辑。

0–24h（肠道内环境，缺氧 + 37℃）：双环境信号激活 AND 门，TetR 抑制蛋白快速表达并维持在 300 nM 左右的高稳态水平，强效抑制 pTet 启动子的转录活性，因此 MazF 毒素始终维持在极低水平；菌群净增长率为正（绿色生长区），活菌数正常增殖，完全符合设计的「肠道内正常存活定殖」的核心需求。

24h（红色虚线，排出体外，环境切换为有氧 + 室温）：缺氧、37℃ 双输入信号同时消失，AND 门输出关闭，TetR 蛋白停止合成并按一级动力学衰减；当 TetR 浓度跌破 15 nM 的半抑制阈值后，对 pTet 的抑制作用完全解除，MazF 毒素开始快速积累，迅速突破 10 nM 的致死浓度阈值。

24–72h（体外死亡过程）：MazF 毒素高表达后，菌群净增长率由正转负（红色死亡区），并快速降至 -1.3 h⁻¹ 左右的最大死亡速率，活菌数呈指数级下降，最终在 60h 左右降至检测限以下，实现了工程菌排出体外后的彻底清除。

设计优化提示：图中清晰呈现了系统的延迟特性，工程菌排出体外后存在约 10h 的开关切换窗口期，与模型计算的未优化响应时间完全一致，为后续提升 TetR 降解速率、缩短响应时间的优化方案提供了直观的时序依据。

6.2 MazF 达到致死浓度的时间

MazF 积累后达到 $K_{\text{kill}}$ 的额外时间：

$$t_{\text{kill}} \approx \frac{K_{\text{kill}}}{\alpha_{\text{MazF}}} = \frac{10}{80} \approx 0.125\;\text{h} \approx 7.5\;\text{min}$$

6.3 总响应时间与优化

总响应时间：$t_{\text{total}} = t_{\text{switch}} + t_{\text{kill}} \approx 10.1$ h。注意：此响应时间较长，工程菌排出体外后可能在约 10 小时内仍具有活性。优化建议：① 降低 TetR 蛋白稳定性（增加 ssrA 降解标签），将 $\delta_p$ 提升至 1–2 h⁻¹；② 增加 MazF 的基础泄漏表达。添加 ssrA 降解标签后 $\delta_{\text{TetR}} = 1.5$ h⁻¹：$t_{\text{switch,opt}} \approx 2.0$ h，$t_{\text{total,opt}} \approx 2.1$ h。

本图对应模型量化了 TetR 蛋白降解速率 $\delta_{\text{TetR}}$ 对自杀开关启动速度、毒素积累效率和菌群清除效果的影响，验证了添加 ssrA 降解标签的优化收益。

• 左图（TetR 衰减）：展示工程菌排出体外后，不同降解速率下 TetR 蛋白的衰减动力学。初始稳态 TetR 浓度一致，降解速率越高（$\delta$ 从 0.3 h⁻¹ 提升至 2.0 h⁻¹），TetR 浓度下降越快，能更快跌破 15 nM 的半抑制常数阈值（虚线），解除对 MazF 毒素的转录抑制，是决定自杀开关启动速度的核心环节。

• 中图（MazF 积累）：与 TetR 衰减过程完全对应，TetR 降解速率越高，MazF 毒素的启动积累时间越早，浓度上升速度越快，可快速突破 10 nM 的致死半饱和常数。未优化的 $\delta=0.3$ h⁻¹ 组，MazF 在 20 h 后才达到高浓度；而优化后 $\delta=2.0$ h⁻¹ 组，在 5 h 左右就进入快速积累期，大幅缩短了杀伤启动时间。

• 右图（菌群清除速度）：直观呈现了优化后的安全收益，TetR 降解速率越高，工程菌活菌数下降速度越快，半衰期越短。$\delta=2.0$ h⁻¹ 组在 10 h 内活菌数就降至检测限以下，而未优化组在 25 h 仍有 10⁸ CFU 以上的活菌，验证了提升 TetR 降解速率可大幅缩短自杀开关的响应时间，显著降低工程菌的环境逃逸风险。

十五、安全模块灵敏度与模型总结

7. 灵敏度分析

定义灵敏度系数 $S_i = \partial\ln([\text{MazF}]_{\text{ss}}) / \partial\ln(p_i)$。关键参数：$\alpha_{\text{MazF}}$ 增大→杀伤更快；$K_{\text{TetR}}$ 增大→泄漏表达增加；$n_2$ 越大开关性越好；$\delta_p$ 增大→切换更快；$\beta_{\text{MazF}}$ 适度泄漏可提高安全性。

7.2 鲁棒性评估

系统在以下参数波动范围内仍维持安全性（体外必死）：α_MazF 波动±50%，体外 MazF 稳态仍 >100 nM，安全√；K_TetR 波动±50%，影响切换时间但不改变最终结果，安全√；温度 20–30℃ 范围，Y38 启动子输出接近 0，安全√。

8.1 核心结论

指标	数值	评价
体外自杀启动时间（未优化）	~10 h	偏长，需优化
体外自杀启动时间（优化后）	~2 h	可接受
体外菌群半衰期	~32 min	高效杀伤
单代逃逸概率	~2×10⁻⁶	需冗余设计
双冗余逃逸概率	~4×10⁻¹²	高度安全

8.2 模型局限性与改进方向

随机性：菌量少时确定性 ODE 不准确，应引入 Gillespie 算法进行随机模拟。
空间异质性：肠道内氧浓度梯度，应考虑 PDE 模型。
代谢负担：安全模块表达对宿主菌的代谢负担未纳入模型。
实验校准：参数多基于文献估计，需通过实验数据（荧光报告基因等）进行参数拟合。