21、基于深度学习的新冠潜在药物分子识别模型

基于深度学习的新冠潜在药物分子识别模型

1. 引言

2019 年新型冠状病毒引发的 COVID - 19 是 21 世纪重大医疗挑战之一。多数治疗方法效果不佳，即便接种疫苗，感染风险仍较高。目前缺乏有效治疗手段，因此急需通过药物发现方法开发强效治疗化合物。

药物发现主要有小分子开发和抗体开发两种途径。前者用计算方法创建分子作为配体抑制靶蛋白，后者通过抗体结合病毒蛋白表面阻止其附着宿主细胞受体。近期探索潜在药物多依赖药物再利用，但效果欠佳。而基于计算方法的从头分子设计在新分子生成领域表现出色。

在化学信息学领域，对生成精确逼真化合物的生成模型需求日益增长，常用简化分子线性输入规范（SMILES）字符串表示化合物。然而，巨大的搜索空间给药物设计带来挑战，据估算，有超过(10^{60})种合成可用的类药物化合物。

虚拟筛选可通过相似性度量寻找合适类药物化合物，而从头药物设计则是让创建的独特化合物结合生物靶标。深度学习生成模型在文本、图像生成等多领域取得进展，为提供潜在候选分子打开了大门。

由于 COVID - 19 爆发，新疾病出现可能性增加，且现有治疗手段过度使用导致耐药性，因此需要开发新药物。深度学习和分子对接是受关注的药物开发计算方法，但分子对接存在获取靶蛋白 3D 结构困难、模拟耗时耗钱等缺点。

为此，提出 Deep - DTI 模型用于从头药物开发，采用卷积神经网络（CNN）进行判别，基于长短期记忆网络（LSTM）的生成对抗网络（GAN）完成任务，并使用强化学习更新生成模型权重，使对抗网络学习化合物分布，生成与病毒蛋白生理直接相互作用的新化合物。

2. 文献综述

2.1 分子设计的生成神经网络架构

常见的分子生成架构有循环神经网络（RNN）、变分自编码器（VAE）、生成对抗网络（GAN）和对抗自编码器（AAE）。
- GAN ：由生成模型（G）和判别模型（D）组成，G 从先验数据分布采样生成新数据点，D 区分真实数据和生成数据。二者通过随机梯度下降训练，可看作一个极小极大博弈。
- AAE ：受 VAE 和 GAN 启发，作为传统自编码器（AE）正则化技术，用 KL 正则化使后验分布匹配先验分布。
- VAE ：常通过 KL 正则化对潜在代码施加先验分布。

分子输入通常有 2D 无向分子图和 SMILES 线性表示两种形式，SMILES 字符串是生成模型中常用的分子表示方式。不同研究基于 RNN、VAE、AAE 等构建了多种生成模型，用于不同任务，如生成疟原虫、5 - HT2A 相关分子，抑制 JAK2 激酶等。

虽然已有多种利用深度学习进行药物再利用的方法，但寻找新型类药物分子的治疗方案研究较少。本文提出的 Deep - DTI 是一种基于强化学习开发从头药物样化学物质的对抗模型，可确定药物与靶蛋白的结合亲和力。

3. 材料与方法

3.1 分子表示

采用简化分子输入线输入系统（SMILES）格式以机器可读方式表示分子。其基于原子符号表示原子、方括号表示金属原子、特定符号表示单、双、三键等约定。通过深度学习技术记录序列中的抽象特征和模式，使生成模型理解 SMILES 语言规则并生成合法分子。

3.2 数据集收集

生成模型的分子 SMILES 数据集由 ChEMBL 和 ZINC 两个数据库融合而成，共约两百万个 SMILES。使用 RDKit 工具进行数据预处理，包括功能团归一化、离子酸碱中和等，仅使用长度在 32 到 128 个字符之间的 SMILES 训练生成模型，预处理后得到 130 万个 SMILES 。

DTI 模型使用 BindingDB 数据集，包含 52000 个药物 - 靶标配对及其对应的解离常数（Kd）分数。为保持训练稳定，将 Kd 分数转换为对数尺度（PKd），公式如下：
[P_{Kd} = \log_{10}K_d / 9]

3.3 生成方法生成新分子

3.3.1 GANs

生成对抗网络（GANs）可生成与真实数据相似的合成数据。RNN 网络（如 LSTM 和 GRU）可处理序列数据，但存在暴露偏差问题，而 GANs 可减轻该问题。

GANs 由判别器（D）和生成器（G）组成。G 从随机分布（Pr）采样合成数据（S），D 用真实数据分布（Pt）和 Pr 生成的样本训练。G 试图生成让 D 无法区分真假的数据，通过判别器性能优化生成器性能，将无监督问题转化为有监督的极小极大问题，损失函数如下：
[\min_{G}\max_{D} E_{s \sim P_t}[\log D(s)] + E_{s \sim P_r}[\log(1 - D(G(s)))]]

生成器架构由两个堆叠的 LSTM 层组成，判别器模型由两个卷积层和一个密集层组成，使用交叉熵作为优化函数。

3.3.2 强化方法更新生成器权重

由于 GAN 模型在处理序列数据时，生成器无法控制离散输出，因此使用强化学习将离散标记的创建视为随机策略。

在强化学习设置中，当前状态（K）是到时间步（t）为止生成的 SMILES 字符总和。生成器 (G_{\theta}) 根据随机策略生成字符，判别器 (D_{\phi}) 区分生成数据和真实数据。

动作值 (Q[K, C]) 表示在状态 K 应用策略 (G_{\theta}) 生成 SMILES 字符的预期奖励，公式如下：
[Q[K, C] = Q[K, C_{t - 1}] + R_t[C]]

随机策略的预期长期奖励（价值函数）由以下公式提供：
[J(\theta) = E_{C \sim G_{\theta}}[R(K, C)] = \sum_{i = 0}^{N} \frac{1}{N} R(K, C_i) \mid C_i \sim G_{\theta}(C \mid K)]

使用蒙特卡罗搜索获取完整序列，在生成器训练过程中，使用有效性度量 (V(C)) 和判别器分类作为总奖励更新生成器权重，总奖励公式如下：
[R(C) = \frac{1}{2}D_{\phi}(C) + \frac{1}{2}V(C)]

当生成更多有效化合物时，通过 Tanimoto 相似度微调生成器，使其生成与新冠主要蛋白酶（3CLpro）结合的分子。

3.4 DTI 模型

3.4.1 1D 卷积神经网络

1D 卷积用于处理序列数据，通过移动滤波器对编码序列进行卷积操作，滤波器值可学习，卷积层输出经池化层下采样，减少特征图大小和可学习参数数量，加速学习。

3.4.2 长短期记忆（LSTM）网络

LSTM 是一种 RNN，可解决传统 RNN 的梯度消失问题。它通过增加输入门和遗忘门等新门控单元，利用自循环记忆单元管理数据在单元间的流动，存储输入标记为连续值。

3.4.3 DTI 架构

输入为蛋白质和药物分子，固定蛋白质输入长度为 1024 字符，药物分子为 128 字符。输入序列编码后送入嵌入层，维度均设为 32。

蛋白质嵌入向量输入到两个堆叠的卷积层，中间使用最大池化，将输入蛋白质序列维度从 (1024×32) 变为 (256×128) 再变为 (128×128)。药物分子嵌入向量通过 LSTM 层，将维度从 (128×32) 变为 (128×128)。

卷积层和 LSTM 层输出通过拼接层合并，再输入到卷积层和 LSTM 层，形成多视图学习范式。最终拼接表示经扁平层和密集层处理，将解离常数（Kd）预测视为回归问题。

3.4.4 实现细节

使用 Python 的 TensorFlow 库作为后端，在 Tesla P100 GPU 上训练模型。使用 RDKit 完成 SMILES 字符串的验证、清理和理化特性计算等预处理。

4. 结果与讨论

4.1 生成器质量指标

训练后，从生成器采样约 10000 个化合物，用有效性、唯一性和新颖性评估质量：
- 有效性（QMValid） ：衡量生成分子中有效分子的比例，公式为 (QM_{Valid} = \frac{V}{S})，其中 (V) 是有效分子总数，(S) 是采样分子总数。
- 唯一性（QMUnique） ：指有效分子中唯一分子的数量，公式为 (QM_{Unique} = \frac{\vert set(V) \vert}{V})。
- 新颖性（QMNovel） ：表示生成分子中不在训练集中出现的数量，公式为 (QM_{Novel} = 1 - \frac{\vert set(V) \cap N \vert}{\vert set(V) \vert})，其中 (N) 是训练样本总数。

生成器生成有效 SMILES 的比例为 93.12%，其中 94.23% 是唯一的，新颖性为 100%。作为对比，训练的基于 LSTM 的生成模型生成有效 SMILES 的比例为 72.15%，有效序列中 97.38% 是唯一的，新颖性为 90.12%。

不同生成模型的质量指标对比如下表所示：
| Architecture Type | Model Name | Dataset used | Size of Molecule | No. of Trained Molecules | No. of Generated Molecules | Generator Quality Metrics | Task |
| — | — | — | — | — | — | — | — |
| RNN - and AE - Based Architectures | Grammar VAE [25] | ZINC | < 39 heavy atoms | 250,000 | 100,000 | 7.2% (V) | Penalized logP |
| RNN - and AE - Based Architectures | SD VAE [13] | ZINC | < 39 heavy atoms | 250,000 | 100,000 | 43.5% (V) | Penalized logP |
| RNN - and AE - Based Architectures | AAE [8] | ChEMBL | < 121 characters | 1.3 million | No data | 77.4% (V) | Drug analogue generation |
| RNN - and AE - Based Architectures | ECAAE [34] | ZINC | < 58 characters | 1.8 million | 10,000 | No data | Structural analogue |
| RNN - Based Architecture with RL | REINVENT [30] | ChEMBL | 10–50 heavy atoms | 1.5 million | 12,800 | 94%(V), 90% (N) | Drug analogue generation |
| RNN - Based Architecture with RL | ReLeaSE [35] | ChEMBL | No data | 1.5 million | 1 million | 95%(V), 95.3% (N) | Inhibitor of JAK2 |
| RNN - Based Architecture with RL | ChemTS [48] | ZINC | No data | 250,000 | No data | No data | Penalized logP |
| RNN - Based Architecture | Segler et al. [42] | ChEMBL | No data | 1.4 million | 976,327 | 97.7%(V), 89.4% (N) | Plasmodium falciparum,5 - HT2A |
| RNN - Based Architecture | Bjerrum et al. [7] | ZINC | No data | 1,611,889 | 50,000 | 98%(V), 63% (N) | Retrosynthetic route |
| RNN - Based Architecture | Gupta et al. [18] | ChEMBL | 34–74 heavy atoms | 541,555 | 30,107 | 93%(V), 92% (N) | PPARs, Trypsin |
| RNN - Based Architecture | Amilpur et al. [2] | ChEMBL, ZINC | 34–128 characters | 290,000 | 10,000 | 70.5% (V), 98.99% (N), 99.83% (U) | CLPro protease |
| GNN - and RNN - Based Architecture | ATNC [37] | ChemDiv | < 91 characters | 15000 | 157,986 | 72%(V), 77% (N) No. of unique heterocycles | |
| GNN - and RNN - Based Architecture | RANC [36] | ChemDiv | < 91 characters | 15000 | 896,000 | 58%(V), 48% (N) No. of unique heterocycles | |
| Deep - DTI | | ChEMBL, ZINC | 48–132 characters | 2.9 million | 10,000 | 93.12% (V) | 3CLpro main protease of COVID - 19 |

4.2 药物 - 靶标相互作用评估指标

使用均方根误差（RMSE）和一致性指数（CI）评估 DTI 模型有效性。
- RMSE ：计算预期值和实际值之间的差异，公式为 (RMSE = \sqrt{\frac{1}{N} \sum_{i = 1}^{N} (y_i - \hat{y} i)^2})，其中 (y_i) 是第 (i) 个样本的实际结果，(\hat{y}_i) 是预测结果。
- CI ：比较两个随机选择的药物 - 靶标配对的预期结合亲和力和实际亲和力顺序，公式为 (CI = \frac{1}{Z} \sum {i < j} h(b_i - b_j) \mid \delta_i > \delta_j)，其中 (Z) 是归一化常数，(h[x]) 是阶跃函数：
[h(x) = \begin{cases}
1, & x > 0 \
0.5, & x = 0 \
0, & x < 0
\end{cases}]

不同模型在 BindingDB 数据集上预测亲和力的 RMSE 和 CI 对比如下表所示：
| Model | RMSE | CI |
| — | — | — |
| GanDTI [46] | 0.721 | – |
| DeepAffinity [22] | 0.740 | |
| DeepDTA [31] | 0.782 | 0.812 |
| MONN [27] | 0.764 | – |
| DeepCDA [1] | 0.899 | 0.822 |
| Deep - DTI | 0.698 | 0.852 |

4.3 新型分子生成

训练 15 个周期后，因生成器样本有效性未提高，提前停止训练。生成器生成的分子中，有效 SMILES 比例达 93.12%，其中 94.23% 唯一，新颖性 100%。通过将生成器度量从有效性改为 Tanimoto 相似度，在新数据集上运行模型，生成能与 3CL 蛋白酶结合的化合物。

4.4 生成针对 3CLpro 蛋白酶的活性药物分子的额外训练

训练后生成器能产生合法化合物，通过收集现有抗病毒药物、HIV 抑制剂和有效处方药物创建新训练集，使用 Tanimoto 相似度指数作为奖励函数训练生成器。选择的 10000 个分子中，96.17% 是可接受的 SMILE 序列，其中 86.23% 唯一，99.12% 新颖。使用 DTI 模型确定这些化合物与 SARS - CoV - 2 蛋白酶的平衡解离常数（Kd）。

4.5 评估新创建的结合 3CL 蛋白酶的化合物

4.5.1 与其他基线模型的比较

为与其他模型有效比较，先以 Ki 为目标值训练模型。在 BindingDB 数据集上，不同模型预测亲和力的 RMSE 和 CI 表现不同，Deep - DTI 模型的 RMSE 为 0.698，CI 为 0.852，优于多数基线模型。去除中间拼接层实验中，RMSE 为 0.718，CI 约为 0.812。

4.5.2 前十个潜在药物分子的结合亲和力评估

使用生成模型的药物化合物测试 DTI 模型对 COVID - 19 主要蛋白酶的预测能力，并与市售抗病毒药物比较。结果显示，分子 ID 为 M001 的化合物预测结合亲和力为 15.37 纳摩尔，远高于瑞德西韦的 1557.35 纳摩尔。M001 - M003 分子与现有抗病毒药物相比表现出色，这些潜在药物可进一步进行体外和体内实验。

生成和现有药物的结合亲和力得分对比如下表所示：
| ID | SMILE | Kd (nM) |
| — | — | — |
| M001 | CNC(=O)c1ccc(OCCN2CCC(O)CC2)cc1 | 15.37 |
| M002 | COc1c c(C(= O)NC( C)c2c cc(OC )c(OC )c2)c cc1c1 ccco1 | 33.13 |
| M003 | COc1c c(C=C c2ccc (C(=O )NC3C C3)cc 2)cc( OC)c1 OC | 40.92 |
| M004 | CC(C)C(=O)Nc1nc2ccc(C#CC(=O)O)cc2cn1 | 46.32 |
| M005 | CC(C)NC(=O)c1cccc(C(CC(=O)O)c2ccccc2)c1 | 50.78 |
| M006 | Cc1cc c(C(= O)Nc2 ccc(F )c(C3 (C)C= C(N)C 3)n2) cc1N( C)C | 75.73 |
| M007 | COC(= O)c1c cc2nc (N3CC N(C(= O)C(C )C)CC 3)[nH ]c2c1 | 78.01 |
| M008 | O=C(O )C1=C (C2CC N(CCc 3ncc4 ncncc 4n3)C C2)C1 =O | 102.08 |
| M009 | CC(C)C(=O)NC(Cc1ccccc1)C(O)CNCCc1ccccc1 | 168.75 |
| M010 | CC1(C)OC(=O)C(CCCC2CCCCC2)Cc2ccccc21 | 225.46 |
| Dolutegravir | CC1CC OC2N1 C(=O) C3=C( C(=O) C(=CN 3C2)C (=O)N CC4=C (C=C( C=C4) F)F)O | 2458.33 |
| Ritonavir | CC(C) C1=NC (=CS1 )CN(C )C(=O )NC(C (C)C) C(=O) NC(CC 2=CC= CC=C) CC(C(CC3=CC=CC=C3)NC(=O)OCC4=CN=CS4)O | 2247.78 |
| Entecavir | C=C1C(CC(C1CO)O)N2C=NC3=C2N=C(NC3=O)N | 2215.04 |
| Remdesivir | CCC(C C)COC (=O)C (C)NP (=O)( OCC1C (C(C( O1)(C #N)2 =CC=C 3N2N= CN=C3 N)O)O )OC4= CC4 | 1557.35 |
| Efavirenz | C1CC1 C#CC 2(C3= C(C=C C(=C3 )Cl)N C(=O) O2)C( F)(F) F | 1507.03 |
| Atazanavir | CC(C) (C)C( C(=O) NC(CC 1=CC= CC=C1 )C(CN (CC2= CC=C( C=C2) C3=CC =CC=N 3)NC( =O)C( C(C)( C)C)N C(=O) OC)O) NC(=O )OC | 48.68 |
| Abacavir* | C1CC1 NC2=C 3C(=N C(=N2 )N)N( C=N3) C4CC( C=C4) CO | 43.51 |

注：* 表示 HIV 药物。

综上所述，提出的基于生成对抗网络的强化学习技术，结合分子度量作为目标，能产生合法化合物，调整后可生成选择性结合新冠主要蛋白酶的物质。基于深度学习的药物 - 靶标相互作用模型在预测平衡常数 [Ki] 值方面优于其他现有模型，为开发新治疗化合物提供了潜在途径。

5. 总结与展望

5.1 研究成果总结

本研究提出了一种基于生成对抗网络的强化学习技术，以分子度量为目标，成功地生成了大量合法的化合物。该模型深入掌握了分子语法规则，能够生成符合化学规则的分子结构。进一步对模型进行调整后，使其能够产出选择性结合新冠主要蛋白酶的物质，为抗击新冠疫情提供了新的药物研发方向。

同时，我们还提出了一种创新的基于深度学习的药物 - 靶标相互作用（DTI）模型。该模型在 BindingDB 数据集上进行了测试，结果表明它在预测平衡常数 [Ki] 值方面显著优于其他现有模型。经过训练，该模型能够准确预测解离常数（Kd）分数，为评估药物与靶标之间的结合亲和力提供了可靠的方法。

5.2 生成化合物的优势

生成的药物分子展现出了与新冠主要蛋白酶良好的结合能力。通过对生成化合物的结合亲和力进行评估，我们发现部分分子的结合亲和力远高于现有的抗病毒药物，如瑞德西韦。例如，分子 ID 为 M001 的化合物预测结合亲和力为 15.37 纳摩尔，而瑞德西韦的预测结合亲和力仅为 1557.35 纳摩尔。这表明我们生成的化合物具有成为新型抗病毒药物的潜力，有望为新冠治疗带来新的突破。

5.3 研究的意义与价值

本研究的成果对于药物研发领域具有重要的意义。在新冠疫情持续蔓延的背景下，寻找有效的治疗药物迫在眉睫。我们的研究提供了一种新的药物研发思路，通过生成对抗网络和强化学习技术，能够快速、高效地生成大量潜在的药物分子，并通过 DTI 模型准确评估其与靶标的结合亲和力，大大缩短了药物研发的周期和成本。

此外，我们的研究也为解决其他疾病的药物研发问题提供了借鉴。随着科技的不断发展，越来越多的疾病需要新的治疗方法和药物。我们的方法可以应用于不同疾病的药物研发中，为攻克各种疑难病症提供有力的支持。

5.4 未来研究方向

尽管本研究取得了一定的成果，但仍有许多方面需要进一步探索和完善。以下是一些未来的研究方向：
1. 优化生成模型 ：进一步改进生成对抗网络和强化学习算法，提高生成化合物的质量和多样性。可以尝试引入更多的分子特征和约束条件，使生成的化合物更符合药物研发的要求。
2. 拓展 DTI 模型 ：扩大 DTI 模型的应用范围，不仅可以预测药物与新冠主要蛋白酶的结合亲和力，还可以预测与其他病毒蛋白或疾病相关靶标的结合亲和力。同时，提高 DTI 模型的准确性和泛化能力，使其能够更好地适应不同的数据集和任务。
3. 实验验证 ：对生成的潜在药物分子进行体外和体内实验验证，评估其抗病毒活性和安全性。通过实验验证，筛选出具有实际应用价值的药物分子，并进一步优化其结构和性能。
4. 结合其他技术 ：将本研究的方法与其他先进的技术相结合，如人工智能、机器学习、量子计算等，进一步提高药物研发的效率和准确性。例如，可以利用人工智能技术对大量的分子数据进行分析和挖掘，发现新的药物靶点和作用机制。

5.5 研究流程回顾

为了更清晰地展示本研究的流程，以下是一个 mermaid 格式的流程图：

graph LR
    A[问题提出] --> B[数据收集与预处理]
    B --> C[生成模型训练]
    C --> D[生成化合物评估]
    D --> E[DTI 模型训练与验证]
    E --> F[结合亲和力预测]
    F --> G[潜在药物筛选]
    G --> H[实验验证与优化]

5.6 总结表格

为了方便对比不同模型的性能，以下是对生成模型和 DTI 模型的总结表格：
| 模型类型 | 模型名称 | 数据集 | 生成有效 SMILES 比例 | 唯一性 | 新颖性 | RMSE | CI |
| — | — | — | — | — | — | — | — |
| 生成模型 | Deep - DTI | ChEMBL, ZINC | 93.12% | 94.23% | 100% | - | - |
| 生成模型 | LSTM - based | ChEMBL, ZINC | 72.15% | 97.38% | 90.12% | - | - |
| DTI 模型 | Deep - DTI | BindingDB | - | - | - | 0.698 | 0.852 |
| DTI 模型 | GanDTI | BindingDB | - | - | - | 0.721 | - |
| DTI 模型 | DeepAffinity | BindingDB | - | - | - | 0.740 | - |
| DTI 模型 | DeepDTA | BindingDB | - | - | - | 0.782 | 0.812 |
| DTI 模型 | MONN | BindingDB | - | - | - | 0.764 | - |
| DTI 模型 | DeepCDA | BindingDB | - | - | - | 0.899 | 0.822 |

通过以上表格可以直观地看到不同模型在生成化合物质量和预测结合亲和力方面的差异，为进一步的研究和应用提供了参考。

总之，本研究为新冠药物研发提供了一种新的方法和思路，通过生成对抗网络和强化学习技术生成潜在的药物分子，并通过 DTI 模型评估其与靶标的结合亲和力。未来的研究将继续深化和拓展这些方法，为解决新冠疫情和其他疾病的药物研发问题做出更大的贡献。