28、基于AMMAS的基因设计优化算法解析

基于AMMAS的基因设计优化算法解析

在合成生物学领域，设计具有特定功能的基因构建体是一项关键任务。为了高效地完成这一任务，研究人员提出了一种基于自适应最大 - 最小蚁群系统（AMMAS）的算法，用于指导用户选择最优的生物部件组合。本文将详细介绍该算法的原理、步骤以及应用实例。

1. 算法核心机制

1.1 自适应信息素浓度更新机制

信息素更新机制旨在将更多的信息素浓度分配给短路径，类似于强化学习模式，更好的解决方案将获得更高的强化。信息素更新公式模拟了蚂蚁在访问边上沉积新信息素和信息素蒸发的过程。为防止算法因最差路径和最佳路径之间信息素密度差异过大而陷入局部最优，设计了动态调整信息素挥发系数的公式：
[
\rho(N + 1) = \frac{1}{\log_2(1 + N_c)}
]
在研究中，信息素挥发系数 $\rho$ 逐渐减小，直到 $\frac{1}{\log_2(1 + N_c)} < \rho_{min}$。

1.2 权重系数的自适应变化

为了增强搜索过程的强度和多样性，使找到的解决方案更具鲁棒性，设计了权重系数 $\alpha$ 和 $\beta$ 的动态变化机制。当当前全局最优路径在 50 轮内未发生变化时，通过降低信息素的相对影响并增加启发式信息的相对影响，来保持搜索新解决方案的高能力，防止算法陷入局部最优。$\beta$ 的计算公式如下：
[
\beta = \frac{sl}{3 \cdot sum_column}
]
其中，$sum_column$ 是构建的晶格的总列数，$sl$ 反映了在一次迭代中通过最优解所有边的蚂蚁总数。

1.3 状态转移规则

为了在利用先验知识和探索新路径之间实现更好的平衡，蚂蚁 $k$ 从节点 $i$ 选择下一个节点 $j$ 时采用伪随机率规则，具体描述如下：
[
j =
\begin{cases}
\arg\max_{u\in allowed_k} {[\tau_{iu}(t)]^{\alpha} \cdot [\eta_{iu}(t)]^{\beta}}, & \text{if } q \leq q_0 \
p_{ij}^k(t), & \text{else}
\end{cases}
]
[
p_{ij}^k(t) =
\begin{cases}
\frac{(\tau_{ij}(t))^{\alpha}(\eta_{ij}(t))^{\beta}}{\sum_{l\in Allowed_k(i)} (\tau_{il}(t))^{\alpha}(\eta_{il}(t))^{\beta}}, & j \in Allowed_k(i) \
0, & \text{else}
\end{cases}
]
其中，$q$ 是在 $[0, 1]$ 上均匀分布的随机数，变量 $q_0$ 定义如下，它决定了利用先验信息和探索新路径之间的相对重要程度。$\alpha$ 和 $\beta$ 分别强调信息素浓度和启发式信息的重要程度，$\eta_j$ 表示指导选择下一个节点的启发式信息。
[
q_0 = \frac{sum_column}{sl}
]
[
\eta_j = P(part_j|part_i) = \frac{P(part_i, part_j)}{P(part_i)}
]

2. AMMAS 算法步骤

AMMAS 算法分为两个主要步骤：创建晶格和应用改进的算法求解统计语言模型。具体步骤如下：
1. 初始化 ：
- 初始信息素信息 $\tau_{ij} = \tau_{max}$（常数）。
- 算法迭代计数器 $N_c$。
- 状态转移规则中涉及的可调参数 $q_0$。
- 最大迭代次数 $N_{max}$。
- 每一轮蚂蚁经过的最大节点数 $N_m$。
- 变量 $n$ 表示一只蚂蚁经过的节点数，初始化为 $n = 1$。
- 设定阈值 $U$（$U \leq 10$），用于确定算法早期蚂蚁是否根据公式 (17) 选择下一个节点。
- 将 $m$ 只蚂蚁随机分配到已建立晶格的第一列的所有节点。
2. 节点选择 ：
- 如果 $N_c \leq U$，蚂蚁 $k$ 根据公式 (17) 选择下一个节点，并将该节点放入禁忌列表 $Tabu_k$。
- 如果 $N_c > U$，蚂蚁根据 $q_0$ 和 $q$ 的比较选择合适的节点选择公式。
3. 迭代判断 ：
- 如果蚂蚁 $k$ 在本轮经过的节点数未达到给定的 $N_m$，算法返回步骤 2 决定下一个节点；否则，进入步骤 6。
4. 信息素更新 ：
- 当所有蚂蚁完成一次迭代后，获取本轮所有候选路径中概率最大的路径。
- 根据公式 (14) 改变信息素挥发因子 $\rho$。
- 全局更新所有路径上的信息素信息，并引入奖惩机制，在好的路径上分配更多信息素，在坏的路径上分配更少信息素，以加速收敛速度。
- 信息素更新公式如下：
[
\tau_{ij}(t + 1) = (1 - \rho)\tau_{ij}(t) + \sum_{k=1}^{m} \Delta\tau_{ij}^{good}
]
[
\tau_{ij}(t + 1) = (1 - \rho)\tau_{ij}(t) - \sum_{k=1}^{m} \Delta\tau_{ij}^{bad}
]
其中，$\Delta\tau_{ij}^{good} = w \cdot f(s_{best})$，$\Delta\tau_{ij}^{bad} = w \cdot f(s_{worst})$，$w$（$0 \leq w \leq m$）表示边 $(i, j)$ 在所有候选路径中出现的总数，$f(s_{best})$ 和 $f(s_{worst})$ 分别表示本轮目标函数的最佳和最差值。
5. 全局最优更新 ：
- 将本轮获得的最优路径与全局最优路径进行比较，如果更好，则替换全局最优路径。
- 如果全局最优路径在 50 轮内未发生变化，进一步动态改变 $\beta$ 以跳出局部最优，扩大搜索范围。
- 根据 $T \leftarrow aT$，$t \leftarrow t + 1$ 更新温度 $T$。
- 如果 $T \geq T_{min}$，算法返回步骤 2 开始新的迭代。
- 如果 $T < T_{min}$ 且 $H = H_{max}$，算法输出保存的全局优化路径。
- 如果 $T < T_{min}$ 且 $H < H_{max}$，则 $H \leftarrow H + 1$，$T \leftarrow T_{max}$，并返回步骤 2。

3. 语法规则与应用实例

为了设计基因构建体，研究人员基于基本语法 PRO（启动子） - RBS（核糖体结合位点） - GEN（基因） - TERM（终止子）实现了一组规则，完整的语法模型类似于上下文无关语法（CFG），如下表所示：
| 规则代码 | 规则 | 描述 |
| ---- | ---- | ---- |
| 1 | CAS → 2CAS | 将一个盒式结构转换为两个盒式结构（CAS） |
| 2 | CAS → Pro CIS TERM | 将一个盒式结构转换为一个启动子（Pro）、一个顺反子（CIS）和一个终止子（TERM） |
| 3 | CAS → Pro CIS | 将一个盒式结构转换为一个启动子（Pro）和一个顺反子（CIS） |
| 4 | CIS → 2CIS | 将一个顺反子转换为两个顺反子（CIS） |
| 5 | CIS → RBS GEN | 将一个顺反子转换为一个 RBS 和一个基因（GEN） |
| 6 | TERM → 2TERM | 将一个终止子转换为两个终止子（TERM） |
| 7 | GEN → 2GEN | 将一个基因转换为两个基因（GEN） |

为了说明如何使用该算法将符合 BioBrick 标准的一系列部件组装成一个功能性生物合成系统，研究人员选择了多个实例进行验证，包括冬青气味生物合成系统、香蕉气味生物合成系统等。以冬青气味生物合成系统为例，通过应用上述语法规则，将设计逐步转换为 PRO - RBS - GEN - TERM - TERM - PRO - RBS - GEN - TERM - TERM 的形式。输入部件类别序列后，算法可以自动选择最优的部件组合，模拟结果与验证的部件组合非常相似。

4. 算法流程图

graph TD;
    A[初始化] --> B{是否 Nc <= U};
    B -- 是 --> C[按公式17选节点];
    B -- 否 --> D{q <= q0};
    D -- 是 --> E[按规则选节点];
    D -- 否 --> F[按概率选节点];
    C --> G[放入禁忌列表];
    E --> G;
    F --> G;
    G --> H{节点数达到 Nm};
    H -- 否 --> B;
    H -- 是 --> I[获取最大概率路径];
    I --> J[更新挥发因子 ρ];
    J --> K[全局更新信息素];
    K --> L{本轮最优 > 全局最优};
    L -- 是 --> M[替换全局最优];
    L -- 否 --> N{全局最优 50 轮未变};
    N -- 是 --> O[动态改变 β];
    N -- 否 --> P[更新温度 T];
    M --> P;
    O --> P;
    P --> Q{T >= Tmin};
    Q -- 是 --> B;
    Q -- 否 --> R{H = Hmax};
    R -- 是 --> S[输出全局优化路径];
    R -- 否 --> T[H + 1, T = Tmax];
    T --> B;

综上所述，该算法在基因构建体的设计中表现出了良好的性能，能够快速找到与实际部件组合高度相似的最优解决方案。然而，由于使用的双元统计语言模型存在一定局限性，无法很好地揭示部件类别之间的相互作用关系，因此在未来的研究中，引入更高阶的语法模型可能是进一步改进算法的方向。

基于AMMAS的基因设计优化算法解析

5. 多实例验证结果

除了冬青气味生物合成系统，研究还选取了多个其他实例来验证算法的有效性，具体如下：
- 香蕉气味生物合成系统 ：该系统由iGEM团队在2006年设计和实现，包含两个表达盒。通过应用语法规则，将设计逐步转换为特定形式。输入部件类别序列后，算法自动挑选出的最优部件组合，其模拟结果与验证的部件组合十分接近。
- RBS.GFP + PBad CFP设计 ：同样是2006年iGEM团队的设计，在输入特定的设计序列后，算法推荐的部件系列与实际组合非常相似。
- I0500.Q04400.E0430实例 ：由2004年参与iGEM的团队设计，算法推荐的部件组合与实际结果相比，只有两个基本部件不同。
- http://parts.igem.org/Part:BBa_S01664设计 ：2004年团队的设计，算法模拟结果与验证结果相似。
- “AHL - 依赖逆变器”设计 ：2006年团队的设计，算法模拟结果与实际部件组合仅有一个不同的基本部件。
- QPI Test Construct Intermediate实例 ：2004年团队的设计，算法输出的模拟结果与有效部件系列非常相似。

这些实例的模拟结果表明，该算法能够快速有效地为输入的部件类别序列找到高度相似的最优部件组合。与参考文献中用于解决3 - 克或4 - 克模型的动态规划算法相比，该双元模型算法使用更少的变量和计算量，也能得出与设计实际部件组合高度相似的装配方案，证明了该方案的实际意义。

6. 算法优势与局限性分析

• 优势 ：
  • 高效性 ：能够在短时间内为基因构建体找到最优的部件组合，大大减少了生物实验中的冗余操作、时间和成本。
  • 适应性 ：用户可以根据已知信息，在输入部件系列时预先确定某些部件或给予更高的评估，还可以排除某些部件并迭代算法，以满足设计的其他需求。
  • 科学性与可靠性 ：利用提取的统计参数和自适应最大 - 最小蚁群系统（AMMAS）求解2 - 克数学模型，得到的最优部件组合科学可靠。
• 局限性 ：
  • 数据库稀疏问题 ：在一些实例中，如R0010总是与B0012同时出现，部分原因是使用的数据库稀疏。
  • 模型局限性 ：双元统计语言模型只能简单地表示一个部件是否能在设计中启用与前一个部件有关，无法很好地揭示DNA组装过程的真实机制。例如，基因的高效表达不仅与启动子有关，还与RBS、质粒骨架和其他调控序列有关。

7. 改进方向与应用展望

为了克服双元统计语言模型的局限性，需要引入更高阶的语法模型来描述部件类别之间的相互作用关系。更高阶的模型意味着设计中涉及的一个部件与多个先前的部件相关，但这些模型涉及的变量众多，计算复杂度较高。

该算法不仅可以在GenoCAD等机器人平台的设计最后一步选择优化的部件组合，还可以用于自动化DNA组装过程。通过输入设计的部件类别序列，算法可以根据iGEM网站上先前的成功组装案例，自动生成一组合适的部件来形成基因构建体。这对于合成生物学家设计新项目具有重要的指导意义，能够帮助他们更高效地完成基因构建体的设计任务。

8. 总结

本文介绍的自适应最大 - 最小蚁群系统（AMMAS）算法，通过自适应信息素浓度更新机制、权重系数的自适应变化和状态转移规则等核心机制，结合特定的算法步骤，能够有效地解决生物部件选择问题，为基因构建体的设计提供了一种高效的解决方案。多个实例的验证结果表明，该算法在实际应用中表现良好，但也存在一定的局限性。未来，引入更高阶的语法模型有望进一步提升算法的性能，使其更好地模拟真实的DNA组装过程，为合成生物学的发展提供更有力的支持。

实例名称	设计团队年份	模拟结果与实际结果差异
冬青气味生物合成系统	2006年MIT iGEM	非常相似
香蕉气味生物合成系统	2006年iGEM团队	非常接近
RBS.GFP + PBad CFP设计	2006年iGEM团队	十分相似
I0500.Q04400.E0430实例	2004年iGEM团队	两个基本部件不同
http://parts.igem.org/Part:BBa_S01664设计	2004年iGEM团队	相似
“AHL - 依赖逆变器”设计	2006年iGEM团队	一个基本部件不同
QPI Test Construct Intermediate实例	2004年iGEM团队	非常相似

graph LR;
    A[算法优势] --> B[高效性];
    A --> C[适应性];
    A --> D[科学性与可靠性];
    E[算法局限性] --> F[数据库稀疏];
    E --> G[模型局限性];
    H[改进方向] --> I[引入高阶语法模型];
    J[应用展望] --> K[自动化DNA组装];
    J --> L[指导新项目设计];