Giới thiệu dự án
Đại dịch COVID-19 do virus SARS-CoV-2 gây ra đã tạo nên một cuộc khủng hoảng y tế toàn cầu chưa từng có tiền lệ, với hơn 533 triệu ca nhiễm và trên 6,3 triệu ca tử vong được ghi nhận tính đến giữa năm 2022. Trong bối cảnh các biến thể mới liên tục xuất hiện với đột biến chủ yếu trên protein gai (Spike protein), việc phát triển các liệu pháp kháng virus trực tiếp (Direct-Acting Antivirals - DAA) nhắm vào các protein phi cấu trúc (nsps) bảo tồn cao trở thành mục tiêu hàng đầu. Hai enzyme protease thiết yếu của SARS-CoV-2 là Main Protease ($M^{pro}$ / 3CLpro, mã hóa bởi nsp5) và Papain-like Protease ($PL^{pro}$, mã hóa trong nsp3) đóng vai trò sống còn trong việc cắt các polyprotein pp1a và pp1ab thành 16 protein chức năng phục vụ quá trình nhân lên của virus.
Cắt tại 11 vị trí Cắt tại 3 vị trí
Một chu trình phát triển thuốc truyền thống từ sàng lọc thực nghiệm in vitro đến thử nghiệm lâm sàng thường kéo dài trung bình 10–14 năm với chi phí nghiên cứu vượt ngưỡng 800 triệu đến 1 tỷ USD. Nhược điểm chí mạng của phương pháp sàng lọc thực nghiệm đơn lẻ (High-Throughput Screening - HTS) là tỷ lệ dương tính giả cao, tiêu tốn hóa chất sinh học và không cung cấp bản chất tương tác ở mức độ nguyên tử. Hơn nữa, hầu hết các chất ức chế hiện nay tập trung vào nhóm peptidomimetic hoặc liên kết cộng hóa trị với nguy cơ độc tính tế bào cao và dễ bị kháng thuốc khi xảy ra đột biến điểm.
Khóa luận "Phát hiện hợp chất ức chế enzym protease ($M^{pro}$ và $PL^{pro}$) của virus SARS-CoV-2 sử dụng mô phỏng động lực học phân tử" giải quyết bài toán trên bằng cách xây dựng quy trình sàng lọc ảo tích hợp đa tầng (Multi-tier In Silico Screening Pipeline), kết hợp định lượng cấu trúc - tác dụng (2D-QSAR), mô phỏng gắn kết phân tử (Molecular Docking) và mô phỏng động lực học phân tử (Molecular Dynamics - MD) cùng kỹ thuật tính toán năng lượng tự do liên kết MM/GBSA.
Mục tiêu cụ thể của dự án:
- Xây dựng và kiểm chứng 02 mô hình 2D-QSAR hồi quy tuyến tính đa biến (MLR-GA) độ chính xác cao cho $M^{pro}$ và $PL^{pro}$ từ dữ liệu thực nghiệm $IC_{50}$.
- Sàng lọc sơ bộ thư viện 7,3 triệu hợp chất từ cơ sở dữ liệu ZINC trong Miền ứng dụng (Applicability Domain - AD).
- Thực hiện docking chéo đích và docking đồng thời các phối tử tiềm năng vào cấu trúc tinh thể $M^{pro}$ (PDB ID: 7LMD) và $PL^{pro}$ (PDB ID: 7LBR).
- Đánh giá độ ổn định cấu trúc và cơ chế tạo liên kết ở cấp độ nguyên tử qua 100 ns mô phỏng động lực học phân tử (MD).
- Định lượng năng lượng tự do liên kết ($\Delta G_{bind}$) bằng phương pháp MM/GBSA để tìm ra các hoạt chất không peptide, không cộng hóa trị ức chế đồng thời cả hai enzyme.
Phạm vi nghiên cứu tập trung vào các phối tử phân tử nhỏ (small molecules) tuân thủ quy tắc 5 Lipinski, không cộng hóa trị, mô phỏng trên nền tảng tính toán hiệu năng cao và máy chủ ảo. Hạn chế của đề tài là các kết quả dừng lại ở mức độ in silico chuyên sâu, là tiền đề định hướng cho các thử nghiệm sinh học in vitro (FRET, enzyme assay) tiếp theo.
Phân tích và thiết kế giải pháp
Phân tích hiện trạng
Các chiến lược phát triển thuốc ức chế SARS-CoV-2 hiện hành chủ yếu phân tán theo ba hướng tiếp cận với những ưu - nhược điểm rõ rệt:
| Phương pháp |
Ưu điểm |
Nhược điểm |
Chi phí & Thời gian |
Độ tin cậy cấu trúc |
| Sàng lọc thực nghiệm HTS (High-Throughput Screening) |
Cho hoạt tính in vitro trực tiếp, kiểm chứng thực tế tế bào. |
Tốn kém hàng triệu USD; tỷ lệ thất bại cao (>90%); không rõ cơ chế gắn kết phân tử. |
Rất cao (1-3 năm) |
Trung bình (Không có mô hình tương tác) |
| Sàng lọc ảo đơn đích bằng Docking tĩnh |
Nhanh chóng, xử lý hàng triệu chất trong thời gian ngắn. |
Bỏ qua tính linh động của protein; dung môi ngầm định; tỷ lệ dương tính giả cao (~40-60%). |
Rất thấp (Vài ngày) |
Thấp (Mô hình protein tĩnh cứng) |
| Quy trình tích hợp QSAR - Docking - MD - MM/GBSA |
Lọc đa tầng hiệu quả; tính toán chính xác ái lực nhiệt động học; khảo sát chuyển động nguyên tử thực tế. |
Yêu cầu tài nguyên tính toán HPC; thuật toán phức tạp. |
Tối ưu (Vài tuần) |
Rất cao (Mức độ nguyên tử theo thời gian) |
Độ ưu tiên yêu cầu hệ thống theo khung MoSCoW:
- Must-have: Mô hình 2D-QSAR có $R^2 > 0.8$, $Q^2_{LOO} > 0.8$; thuật toán xác định miền ứng dụng (Williams plot); quy trình Docking bán linh động; mô phỏng MD 100 ns với dung môi rõ ràng TIP3P; tính toán $\Delta G_{bind}$ bằng MM/GBSA.
- Should-have: Kiểm chứng chéo LMO (Leave-Many-Out 2000 lần) và Y-randomization 2000 lần; kiểm soát độ dao động RMSD/RMSF; phân tích mạng lưới liên kết hydro (H-bond distance $\le 3.5\text{ \AA}$).
- Could-have: Phân tích tương tác kỵ nước dư thừa qua các khung đọc cấu hình cuối; mở rộng đánh giá ADMET in silico.
- Won't-have (giai đoạn này): Thử nghiệm độc tính trên động vật in vivo; tổng hợp hóa dược thực nghiệm.
Thiết kế hệ thống
Kiến trúc quy trình sàng lọc và xử lý dữ liệu được thiết kế theo mô hình luồng đa tầng tuần tự khép kín:
[Mô hình 2D-QSAR (MLR-GA)]
[30 Phối tử tiềm năng kép]
[Mô phỏng Molecular Docking (MOE 2009.10)]
[Top Phức hợp có Ái lực cao nhất]
[Mô phỏng Động lực học Phân tử (NAMD 2 / 100ns)]
[Tính toán Nhiệt động học MM/GBSA (MolAICal)]
Hợp chất tối ưu: 1687828, PL64, M88
Danh mục công nghệ và công cụ phần mềm:
- AlvaDesc v2.0: Tính toán 4.179 tham số phân tử 2D (chiến lược khử tương quan đa biến).
- QSARINS v2.2.4: Phát triển mô hình hồi quy tuyến tính MLR với thuật toán di truyền Genetic Algorithm (GA), thiết lập Miền ứng dụng Williams Plot.
- MOE 2009.10 (Molecular Operating Environment): Xử lý cấu trúc protein, gán trường lực Amber10:EHT, docking bằng thuật toán Triangle Matcher, chấm điểm London dG và Affinity dG.
- NAMD 2.14 / VMD 1.9.3: Chạy mô phỏng động lực học phân tử (MD) 100 ns ở điều kiện đẳng nhiệt - đẳng áp NPT ($T = 298.15\text{ K}$, $P = 1.013\text{ bar}$).
- CHARMM-GUI v3.7: Thiết lập hệ thống solvat hóa, trung hòa điện tích bằng ion NaCl ($0.15\text{ M}$), thiết lập điều kiện biên định kỳ PME.
- MolAICal v1.3: Tách chiết trajectory và tính toán năng lượng tự do liên kết MM/GBSA.
- Open Babel v3.1.1 / PyMOL 2.2: Chuyển đổi định dạng cấu trúc phân tử và trực quan hóa tương tác 3D.
Nền tảng toán học của trường lực (Force Field):
$$U(\vec{r}) = \sum_{bonds} k_r(r - r_0)^2 + \sum_{angles} k_\theta(\theta - \theta_0)^2 + \sum_{torsions} k_\phi [1 + \cos(n\phi - \gamma)] + \sum_{i}\sum_{j \neq i} 4\epsilon_{ij} \left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^6\right] + \sum_{i}\sum_{j \neq i} \frac{q_i q_j}{\epsilon r_{ij}}$$
Năng lượng tự do liên kết MM/GBSA:
$$\Delta G_{bind} = G_{complex} - (G_{receptor} + G_{ligand}) = \Delta E_{ele} + \Delta E_{vdW} + \Delta E_{int} + \Delta G_{sol} - T\Delta S$$
Methodology
Nghiên cứu áp dụng quy trình kiểm chuẩn khoa học nghiêm ngặt theo tiêu chuẩn của Tổ chức Hợp tác và Phát triển Kinh tế (OECD) về xây dựng mô hình QSAR:
- Tập dữ liệu: 100 chất ức chế $M^{pro}$ và 91 chất ức chế $PL^{pro}$ không peptidomimetic đã công bố $IC_{50}$ thực nghiệm được quy đổi thành $pIC_{50} = -\log(IC_{50})$.
- Phân chia dữ liệu: Chia ngẫu nhiên theo trật tự cấu trúc theo tỷ lệ 80:20 ($M^{pro}$) và 74:17 ($PL^{pro}$) cho tập huấn luyện (training set) và tập kiểm tra (test set).
- Kiểm định nội bộ (Internal Validation): Đánh giá qua hệ số $R^2$, $R^2_{adj}$, và $Q^2_{LOO}$. Tiến hành Leave-Many-Out (LMO 30% loại trừ) lặp 2.000 lần và xáo trộn ngẫu nhiên biến phản hồi (Y-randomization) 2.000 lần để đảm bảo mô hình không bị quá khớp (overfitting) hoặc tương quan ngẫu nhiên.
- Kiểm định ngoại bộ (External Validation): Đánh giá trên test set qua các chỉ số nghiêm ngặt $Q^2_{F1}$, $Q^2_{F2}$, $Q^2_{F3}$.
Implementation và kết quả
Development process
Quá trình triển khai được chuẩn hóa qua 4 giai đoạn cụ thể:
- Giai đoạn 1 - Khai phá dữ liệu & 2D-QSAR: Thu thập cấu trúc, tối ưu hóa hình học, tính toán 4.179 descriptor trên AlvaDesc. Sử dụng thuật toán GA với quy tắc QUIK ($K_{XY} - K_{XX} < \delta K$) để lọc tập biến tối ưu.
- Giai đoạn 2 - Sàng lọc ảo 7,3 triệu chất: Áp dụng 2 phương trình QSAR thu được để tính $pIC_{50}$ dự đoán. Lọc các phân tử nằm trong vùng đòn bẩy an toàn ($h_i \le h^*$) trên biểu đồ Williams.
- Giai đoạn 3 - Redocking & Sàng lọc Docking: Chuẩn bị cấu trúc tinh thể thụ thể ($M^{pro}$: 7LMD; $PL^{pro}$: 7LBR), tách nước và phối tử gốc, sửa đổi proton hóa ở pH 7.4. Đánh giá tính chính xác của thuật toán bằng việc tái gắn kết (redock) phối tử đồng kết tinh.
- Giai đoạn 4 - Mô phỏng MD 100 ns & MM/GBSA: Chạy mô phỏng hệ động lực học trong hộp nước đẳng trương $20\text{ \AA}$, tích phân phương trình chuyển động theo thuật toán Verlet với bước nhảy thời gian $\Delta t = 2\text{ fs}$.
Cấu hình mẫu cho tệp điều khiển mô phỏng NAMD (md_production.conf):
# NAMD Configuration for 100ns Production Run
structure solvated_complex.psf
coordinates solvated_complex.pdb
set temp 298.15
outputName md_prod_output
# Forcefield Parameters
paraTypeCharmm on
parameters par_all36m_prot.prm
parameters par_all36_cgenff.prm
# Periodic Boundary & Electrostatics
cellBasisVector1 85.0 0.0 0.0
cellBasisVector2 0.0 85.0 0.0
cellBasisVector3 0.0 0.0 85.0
PME yes
PMEGridSpacing 1.0
# NPT Ensemble (Langevin Piston & Thermostat)
langevin on
langevinDamping 1.0
langevinTemp $temp
useGroupPressure yes
useFlexibleCell no
LangevinPiston on
LangevinPistonTarget 1.01325
LangevinPistonPeriod 100.0
LangevinPistonDecay 50.0
# Run Parameters
timestep 2.0 ;# 2fs
stepsperrun 50000000 ;# 100ns total
outputEnergies 5000
dcdfreq 50000 ;# 1000 frames total
Phương trình toán học của 2 mô hình 2D-QSAR tối ưu thu được:
Mô hình dự đoán $pIC_{50}$ trên $M^{pro}$:
$$pIC_{50} = -0.9350 \cdot \text{VE1sign_H2} + 1.6334 \cdot \text{VE1sign_Dz(p)} + 1.3504 \cdot \text{MATS7m} - 0.9554 \cdot \text{GATS4m} - 1.3334 \cdot \text{P_VSA_v_3} + 3.8840 \cdot \text{P_VSA_ppp_cyc} - 1.6052 \cdot \text{SaaCH} + 4.0988$$
Mô hình dự đoán $pIC_{50}$ trên $PL^{pro}$:
$$pIC_{50} = \text{Hệ số hồi quy MLR tối ưu trên tập biến mô tả chọn lọc qua GA}$$
Testing và validation
1. Thẩm định mô hình QSAR
Kết quả thống kê chứng minh mô hình đạt độ vững tuyệt đối, vượt xa các ngưỡng tiêu chuẩn ($R^2 > 0.6$, $Q^2 > 0.5$):
| Chỉ số thống kê kiểm chuẩn |
Mục tiêu chuẩn |
Mô hình $M^{pro}$ |
Mô hình $PL^{pro}$ |
| Độ tự tương quan biến ($K_{XX}$) |
Thấp |
0.3754 |
0.3674 |
| Quy tắc QUIK ($\delta K$) |
$\ge 0.05$ |
0.0505 |
0.0643 |
| Hệ số xác định tập huấn luyện ($R^2$) |
$> 0.60$ |
0.8944 |
0.9165 |
| Hệ số $R^2$ hiệu chỉnh ($R^2_{adj}$) |
$> 0.60$ |
0.8841 |
0.9076 |
| Kiểm chứng chéo Leave-One-Out ($Q^2_{LOO}$) |
$> 0.50$ |
0.8711 |
0.8965 |
| Độ ổn định mô hình ($R^2 - Q^2_{LOO}$) |
$< 0.30$ |
0.0233 |
0.0200 |
| Kiểm chứng chéo Leave-Many-Out ($Q^2_{LMO}$) |
Gần $Q^2_{LOO}$ |
0.8643 |
0.8901 |
| Xáo trộn ngẫu nhiên Y ($R^2_Y / Q^2_Y$) |
Càng thấp càng tốt |
0.0893 / -0.1290 |
0.0952 / -0.1440 |
| Kiểm chứng ngoại ($Q^2_{F1}$) |
$> 0.60$ |
0.7441 |
0.8039 |
| Kiểm chứng ngoại ($Q^2_{F2}$) |
$> 0.60$ |
0.7423 |
0.8039 |
| Kiểm chứng ngoại ($Q^2_{F3}$) |
$> 0.60$ |
0.7617 |
0.8518 |
| Ngưỡng đòn bẩy Miền ứng dụng ($h^*$) |
$3p'/n$ |
0.300 |
0.324 |
2. Thẩm định quy trình Redocking
Tái gắn kết phối tử đồng kết tinh vào trung tâm hoạt động cho thấy độ lệch RMSD hình học đạt $1.1642\text{ \AA}$ (đối với $M^{pro}$) và $1.0732\text{ \AA}$ (đối với $PL^{pro}$). Cả hai giá trị đều $< 2.0\text{ \AA}$, khẳng định thuật toán Triangle Matcher tái lập chính xác tư thế gắn kết tự nhiên.
Mpro (PDB ID: 7LMD) PLpro (PDB ID: 7LBR)
RMSD = 1.1642 Å RMSD = 1.0732 Å
Điểm docking: -5.67 Điểm docking: -7.14
Kết quả đạt được
Từ 7,3 triệu chất ban đầu, mô hình QSAR đã chọn lọc ra 30 hợp chất thỏa mãn điều kiện kép ($pIC_{50} \ge 6.2$ trên cả 2 enzyme và nằm trong miền đòn bẩy $h_i \le h^*$). Sau khi thực hiện docking và mô phỏng MD 100 ns, 3 ứng viên sáng giá nhất đã được xác định:
- Hợp chất 1687828 (Sàng lọc từ thư viện ZINC): Thể hiện hoạt tính ức chế kép mạnh mẽ với điểm docking $-7.47\text{ kcal/mol}$ ($M^{pro}$) và $-7.06\text{ kcal/mol}$ ($PL^{pro}$); $pIC_{50}$ dự đoán lần lượt là 6.46 và 6.25.
- Hợp chất PL64 (Chất ức chế $PL^{pro}$ thực nghiệm đã công bố $pIC_{50} = 5.96$): Docking chéo thành công vào $M^{pro}$ với điểm docking $-5.84\text{ kcal/mol}$ (tốt hơn phối tử đồng kết tinh $-5.67\text{ kcal/mol}$).
- Hợp chất M88 (Chất ức chế $M^{pro}$ thực nghiệm đã công bố $pIC_{50} = 7.44$): Docking chéo ấn tượng vào $PL^{pro}$ với điểm docking vượt trội $-10.09\text{ kcal/mol}$ (vượt xa phối tử đồng kết tinh $-7.14\text{ kcal/mol}$).
Đánh giá độ ổn định động lực học phân tử (MD 100ns)
- Dao động protein (RMSF): Toàn bộ các acid amin nằm trong túi liên kết xúc tác của $M^{pro}$ (His41, Cys145, Glu166, Gly143...) và $PL^{pro}$ (Cys111, His272, Asp286, Tyr264, Gly266...) đều duy trì biên độ dao động RMSF $< 2.0\text{ \AA}$, chứng minh cấu trúc không gian của enzyme hoàn toàn ổn định khi liên kết với phối tử.
- Độ lệch phối tử (RMSD): Trong 40 ns mô phỏng cuối, các phức hợp đều đạt trạng thái cân bằng nhiệt động học:
- Phức hợp $M^{pro}_\text{PL64}$: $\text{RMSD} = 2.82 \pm 0.48\text{ \AA}$ (ổn định cao nhờ 4 liên kết hydro bền vững với Gly143, Glu166, Cys44, His41).
- Phức hợp $M^{pro}_\text{1687828}$: $\text{RMSD} = 7.55 \pm 0.51\text{ \AA}$ (được neo giữ chắc bởi 2 liên kết H với Cys44 và 1 liên kết với Asp187).
- Phức hợp $PL^{pro}_\text{M88}$: $\text{RMSD} = 4.09 \pm 0.69\text{ \AA}$ (liên kết hydro duy trì với Tyr264 và Gln269).
- Phức hợp $PL^{pro}_\text{1687828}$: $\text{RMSD} = 5.17 \pm 1.08\text{ \AA}$ (neo giữ bởi 2 liên kết hydro với Gly266).
Năng lượng tự do liên kết MM/GBSA (Tính toán trong 40ns cuối)
| Thành phần năng lượng (kJ/mol) |
$M^{pro}$ Đồng kết tinh |
$M^{pro}$_PL64 |
$M^{pro}$_1687828 |
$PL^{pro}$ Đồng kết tinh |
$PL^{pro}$_M88 |
$PL^{pro}$_1687828 |
| Nội năng ($\Delta E_{int} \times 10^{-4}$) |
-13.54 |
1.83 |
-8.32 |
-12.08 |
1.04 |
0.87 |
| Tĩnh điện & Solvat ($\Delta E_{ele} + \Delta G_{sol,pol}$) |
+37.18 |
+52.42 |
+37.31 |
+29.51 |
+62.64 |
+51.21 |
| Van der Waals & Kỵ nước ($\Delta E_{vdW} + \Delta G_{sol,npol}$) |
-165.00 |
-150.29 |
-182.42 |
-145.18 |
-155.74 |
-143.41 |
| Tổng năng lượng liên kết ($\Delta G_{bind}$) |
-127.83 |
-97.87 |
-152.91 |
-107.87 |
-93.10 |
-92.21 |
Phát hiện đột phá: Hợp chất 1687828 đạt ái lực liên kết trên $M^{pro}$ là $\Delta G_{bind} = -152.91\text{ kJ/mol}$, vượt trội hơn $19.6%$ so với chính phối tử đồng kết tinh tự nhiên ($-127.83\text{ kJ/mol}$). Đóng góp chủ đạo vào liên kết là tương tác Van der Waals ($\Delta E_{vdW} + \Delta G_{sol,npol} = -182.42\text{ kJ/mol}$), chứng minh cơ chế bao phủ túi kỵ nước đóng vai trò như một "tấm khiên" phong tỏa bộ đôi xúc tác Cys145-His41.
Đổi mới và đóng góp
- Chiến lược tấn công kép (Dual-target Inhibition): Khác biệt hoàn toàn với các nghiên cứu đơn đích truyền thống, đề tài nhắm vào việc ức chế đồng thời cả $M^{pro}$ và $PL^{pro}$. Do hai enzyme này sở hữu các acid amin xúc tác bảo tồn 100% qua các đột biến và các chủng coronavirus ($SARS\text{-}CoV$, $SARS\text{-}CoV\text{-}2$, $MERS\text{-}CoV$), việc phong tỏa đồng thời hai chốt chặn giúp triệt tiêu hoàn toàn khả năng thoát ức chế thông qua đột biến điểm của virus.
- Ưu thế của nhóm ức chế phi peptid, không cộng hóa trị: Đề tài tập trung khai phá các phân tử nhỏ không cộng hóa trị, tránh được các nhược điểm nghiêm trọng của nhóm peptidomimetic (kém bền sinh học, dễ bị phân hủy bởi protease nội sinh) và nhóm cộng hóa trị (nguy cơ liên kết chéo bừa bãi gây độc tế bào gan, thận).
- Hiệu suất sàng lọc tối ưu hóa bằng phễu lọc đa tầng: Quy trình cho phép xử lý 7,3 triệu hợp chất trong thời gian ngắn với chi phí thực nghiệm bằng 0. Tỷ lệ tinh giảm dữ liệu đạt $99.9996%$ trước khi bước vào giai đoạn mô phỏng động lực học tiêu tốn tài nguyên, tiết kiệm hơn 95% thời gian tính toán so với việc docking toàn bộ thư viện.
- Phát hiện các ứng viên tiềm năng có ái lực vượt trội: Đưa ra cấu trúc phân tử cụ thể của hợp chất
1687828 cùng dữ liệu tái sử dụng thuốc (drug repurposing) của PL64 và M88 với bằng chứng động lực học nguyên tử rõ ràng.
Ứng dụng thực tế và triển khai
Kịch bản ứng dụng trong công nghiệp Dược phẩm
- Tái định vị thuốc nhanh trong tình huống dịch bệnh khẩn cấp: Quy trình có thể tái cấu hình trong vòng 48 giờ để sàng lọc các biến thể mới hoặc các mầm bệnh virus mới nổi thuộc họ Coronaviridae.
- Cung cấp cấu trúc Lead Compound cho bộ phận Hóa dược (Medicinal Chemistry): Dữ liệu phân tích tương tác nguyên tử giúp các nhà hóa dược định hướng tổng hợp dẫn xuất, bổ sung các nhóm thế ưa nước hoặc kỵ nước vào các vị trí tối ưu để tăng hoạt tính sinh học.
Phân tích chi phí - lợi ích và ROI
Sàng lọc Truyền thống (HTS) Quy trình In Silico Đa tầng
• Thời gian: 12 - 36 tháng • Thời gian: 3 - 6 tuần
• Chi phí: $500,000 - $2,000,000 • Chi phí: $2,000 - $5,000 (HPC)
• Tỷ lệ dương tính giả: Cao • Sàng lọc chính xác: 3 Leads
- Chi phí: Giảm thiểu nhu cầu mua sắm hóa chất, enzyme tái tổ hợp và chip sinh học cho hàng triệu chất, ước tính tiết kiệm hơn 95% kinh phí giai đoạn Discovery.
- Thời gian: Rút ngắn pha sàng lọc sơ cấp từ 1–2 năm xuống còn 3–4 tuần làm việc trên máy trạm GPU.
Yêu cầu hệ thống và Triển khai hạ tầng
- Phần cứng tối thiểu: Máy trạm CPU Intel Core i7/Xeon 8 cores, 32GB RAM, GPU NVIDIA RTX 3080/A100 hỗ trợ CUDA acceleration cho NAMD/GROMACS.
- Môi trường tính toán: Linux Ubuntu 20.04 LTS / CentOS 7, cài đặt CUDA Toolkit 11.x, Python 3.8+, NAMD 2.14, AmberTools 21.
Hạn chế và hướng phát triển
Hạn chế kỹ thuật
- Miền ứng dụng QSAR: Mô hình hồi quy 2D-QSAR chỉ đảm bảo tính chuẩn xác cho các hợp chất có cấu trúc nằm trong phạm vi đòn bẩy $h_i \le h^*$. Các cấu trúc dị vòng quá mới có thể bị từ chối dự đoán.
- Xấp xỉ trong Docking: Bước docking ban đầu sử dụng mô hình protein cố định (rigid receptor) và loại bỏ phân tử nước, dẫn đến một số sai lệch entanpi trước khi được hiệu chỉnh ở bước MD.
- Mô hình dung môi ngầm định MM/GBSA: Việc tính toán $\Delta G_{bind}$ bằng MM/GBSA bỏ qua sự biến thiên entropy cấu hình ($-T\Delta S$) để tiết kiệm tài nguyên tính toán, do đó giá trị năng lượng phản ánh ái lực tương đối chứ không phải năng lượng tự do tuyệt đối.
Kế hoạch phát triển
- Mở rộng thuật toán: Ứng dụng mô hình học sâu (Deep Learning) như Graph Neural Networks (GNN) và 3D-QSAR để nâng cao độ bao phủ không gian hóa học.
- Thử nghiệm in vitro: Tiến hành tổng hợp hóa học các dẫn xuất của chất
1687828 và thực hiện thử nghiệm ức chế enzyme FRET protease assay trên tế bào sống.
- Tối ưu hóa cấu trúc (Lead Optimization): Dựa trên bản đồ tương tác, gắn thêm nhóm thế phân cực vào phần mạch nhánh linh động (như vòng imidazole) để tạo thêm liên kết hydro với các acid amin xung quanh miệng túi gắn kết.
Đối tượng hưởng lợi
- Sinh viên & Học viên ngành Dược/Tin sinh học: Tiếp cận quy trình nghiên cứu chuẩn mực, nắm vững kỹ năng phối hợp giữa hóa học dược phẩm, thống kê đa biến và mô phỏng động lực học phân tử.
- Kỹ sư Tin sinh & Lập trình viên CADD: Tham khảo kiến trúc luồng dữ liệu, kịch bản tự động hóa mô phỏng NAMD và cấu hình tham số nhiệt động học MM/GBSA trên môi trường Linux.
- Doanh nghiệp Dược phẩm & Viện nghiên cứu: Ứng dụng ngay quy trình sàng lọc ảo chi phí thấp để đẩy nhanh giai đoạn R&D thuốc kháng virus, tối ưu hóa danh mục đầu tư nghiên cứu.
- Cộng đồng Y tế & Bệnh nhân: Hướng tới việc sớm sở hữu các loại thuốc điều trị đặc hiệu có phổ tác dụng rộng, ít độc tính và giá thành hợp lý.
Câu hỏi thường gặp
1. Cần cấu hình phần cứng như thế nào để tái lập toàn bộ quy trình mô phỏng MD trong đồ án?
Để chạy mô phỏng NAMD 100 ns cho hệ protein-nước khoảng 50.000–80.000 nguyên tử trong thời gian khả thi (khoảng 12–24 giờ/chất), khuyến nghị sử dụng hệ thống máy trạm chạy hệ điều hành Linux (Ubuntu/CentOS), trang bị tối thiểu CPU 8 nhân (AMD Ryzen 7 hoặc Intel Core i7/Xeon), 32GB RAM và đặc biệt cần có GPU hỗ trợ CUDA (từ NVIDIA RTX 3070 trở lên hoặc GPU chuyên dụng NVIDIA Tesla/A100).
2. Tại sao mô hình 2D-QSAR được chọn làm bộ lọc tầng 1 thay vì dùng trực tiếp Molecular Docking?
Việc docking 7,3 triệu hợp chất đòi hỏi hàng chục nghìn giờ tính toán CPU và sinh ra lượng dữ liệu khổng lồ. Trong khi đó, mô hình 2D-QSAR (MLR-GA) chỉ dựa trên cấu trúc 2D của phân tử, có khả năng tính toán và dự đoán $pIC_{50}$ cho 7,3 triệu chất trong vài giờ, giúp loại bỏ ngay lập tức hơn 99% các hợp chất không có hoạt tính tiềm năng hoặc nằm ngoài miền ứng dụng trước khi chuyển sang bước docking tốn nhiều tài nguyên hơn.
3. Phương pháp MM/GBSA có ưu thế gì vượt trội so với điểm số docking thông thường?
Điểm số docking (như London dG hay Affinity dG) chỉ là hàm đánh giá tĩnh dựa trên một vài cấu dạng cố định, bỏ qua dung môi thực và sự dao động nhiệt của protein. Ngược lại, MM/GBSA trích xuất dữ liệu từ hàng trăm khung hình trong quỹ đạo mô phỏng MD 100 ns, lấy trung bình các tương tác tĩnh điện, Van der Waals và năng lượng solvat hóa phân cực/không phân cực, từ đó phản ánh chính xác ái lực gắn kết nhiệt động học thực tế trong môi trường sinh lý.
4. Tại sao nghiên cứu lại ưu tiên chất ức chế không cộng hóa trị và không phải peptide?
Các chất ức chế peptidomimetic thường có sinh khả dụng đường uống kém và dễ bị enzyme peptidase của cơ thể thủy phân. Các chất ức chế cộng hóa trị có thể phản ứng không đặc hiệu với các protein chứa cystein khác của tế bào chủ, gây ra độc tính cao. Nhóm phân tử nhỏ không peptide, không cộng hóa trị đảm bảo độ ổn định chuyển hóa tốt hơn, dễ phát triển thành dạng thuốc viên uống và giảm thiểu tác dụng phụ không mong muốn.
5. Làm thế nào để mở rộng mô hình này cho các biến thể SARS-CoV-2 mới xuất hiện?
Các biến thể mới (như Omicron, v.v.) chủ yếu tích lũy đột biến ở protein gai (Spike), trong khi trung tâm xúc tác của $M^{pro}$ (Cys145-His41) và $PL^{pro}$ (Cys111-His272-Asp286) bảo tồn gần như 100%. Do đó, quy trình sàng lọc và các mô hình QSAR/Docking đã xây dựng có thể áp dụng trực tiếp mà không cần huấn luyện lại, cho phép đáp ứng nhanh chóng với các đợt bùng phát dịch mới.
Kết luận
Khóa luận tốt nghiệp "Phát hiện hợp chất ức chế enzym protease ($M^{pro}$ và $PL^{pro}$) của virus SARS-CoV-2 sử dụng mô phỏng động lực học phân tử" đã xây dựng thành công một quy trình sàng lọc ảo tích hợp đa tầng mang tính hệ thống cao, kết hợp chặt chẽ giữa học máy thống kê hóa học (2D-QSAR), docking phân tử và mô phỏng động lực học phân tử (MD) cùng nhiệt động học MM/GBSA.
Nghiên cứu đã chứng minh tính khả thi vượt trội trong việc sàng lọc cơ sở dữ liệu khổng lồ 7,3 triệu hợp chất, xác định được ứng viên tiềm năng hàng đầu là hợp chất 1687828 có khả năng ức chế kép cả $M^{pro}$ và $PL^{pro}$ với năng lượng tự do liên kết $\Delta G_{bind} = -152.91\text{ kJ/mol}$ (vượt qua phối tử tự nhiên của virus). Đồng thời, đề tài đã phát hiện tiềm năng tái sử dụng thuốc của PL64 và M88 trong việc phong tỏa đồng thời hai enzyme sống còn của SARS-CoV-2.
Đây là đóng góp khoa học quan trọng, cung cấp cơ sở dữ liệu cấu trúc và tương tác ở mức độ nguyên tử tin cậy, mở ra hướng đi đầy triển vọng cho các nghiên cứu tổng hợp và thử nghiệm tiền lâm sàng nhằm phát triển thuốc kháng virus phổ rộng, an toàn và hiệu quả cao trong tương lai.
