1.1 Background and Motivation                                                                                                         

Allicin (diallyl thiosulfinate, C₆H₁₀OS₂; CAS 539-86-6), the principal bioactive compound of Allium sativum (garlic), exerts well-documented anticancer, anti-inflammatory, antiviral, and antibacterial activities ^[1,2]. Its principal pharmacological limitation, however, is rapid degradation in aqueous environments (t½ ≈ 3.5–4 days at 37 °C, pH 7.4), which severely restricts therapeutic bioavailability ^[3,4]. In our previous work (Morais et al., IJRPAS 2026 — hereafter referred to as the v2.0 study ^[5]), we presented the first Digital Twin (DT) framework for modeling allicin aqueous stability and encapsulation by Al₁₂N₁₂, B₁₂N₁₂, and C₂₄ nanocages. That framework, parametrized by DFT thermodynamic data from Mozafari et al. (2025) ^[6], demonstrated that the Allicin/Al₁₂N₁₂ complex achieves a Nanocage Protection Factor (NPF) of 5.5 × 10⁶, extending the half-life from 3.6 days to beyond 9,999 days under physiological conditions. Despite its predictive power, the v2.0 framework relied exclusively on analytically generated synthetic data for its Gaussian Process (GP) molecular dynamics surrogate — a deliberate methodological choice that, by design, left open the critical question of whether the surrogate predictions are consistent with real MD trajectory data for allicin–protein complexes.

The present study directly addresses this gap. Li & Cheng (2023) ^[7] performed Discovery Studio CHARMm36 molecular dynamics simulations of alliin–6LU7 and allicin–6LU7 complexes and reported both RMSD and RMSF trajectories for the SARS-CoV-2 main protease (Mpro) active pocket. These represent some of the most directly relevant experimental MD data available for allicin–Mpro interactions. Integrating them into the Digital Twin framework creates a fifth layer of structural calibration previously absent, while simultaneously enabling a new pharmacodynamic analysis: quantifying how much the extended aqueous lifetime conferred by nanocage encapsulation translates into sustained Mpro inhibitory exposure.

1.2 Significance of Molecular Dynamics Integration

 Molecular dynamics simulations predict how every atom in a protein-ligand system evolves over time, revealing conformational changes, fluctuation profiles, and binding stability that static docking scores cannot capture ^[8,9]. The RMSD (Root Mean Square Deviation) metric quantifies structural drift from the initial conformation: lower, stable RMSD indicates a thermally equilibrated, persistent complex. The RMSF (Root Mean Square Fluctuation) characterizes per-residue flexibility: lower RMSF at active-site residues signals tighter, more stable ligand engagement ^[8]. Li & Cheng ^[7] established that alliin–6LU7 equilibrates at approximately 1.10 nm RMSD after 50 conformations, whereas allicin–6LU7 exhibits persistent fluctuations that never fully equilibrate — a qualitative difference in binding stability that our previous GP surrogate could not distinguish because it was not calibrated against real trajectory data.

Generative and physics-informed Digital Twins are increasingly recognized as transformative tools in pharmaceutical sciences ^[10]. Wu & Koelzer (2024) ^[11] demonstrated that generative DTs can recreate spatially resolved representations of biological entities with high fidelity, enabling in silico modeling of therapeutic interventions. Building on this paradigm, the present v3.0 framework demonstrates that integrating real MD data into a multi-layer DT not only improves prediction credibility but also enables emergent pharmacodynamic metrics — NEFI and SBIS — that are impossible to derive from DFT thermodynamics alone.

1.3 Objectives

This study aims to: (i) integrate real MD trajectory data (Li & Cheng 2023) into the Digital Twin v2.0 as Layer 4; (ii) develop Layer 5 — a Mpro binding dynamics model yielding the NEFI and SBIS metrics; (iii) cross-validate the GP surrogate predictions against real CHARMm36 RMSD values; (iv) predict the active-pocket RMSF profile of the Allicin/Al₁₂N₁₂ complex using the calibrated surrogate; and (v) identify open experimental questions that motivate subsequent work.

2 MATERIALS AND METHODS

2.1 Digital Twin v3.0 Architecture

The Digital Twin v3.0 inherits all three layers of the previously published v2.0 framework ^[5] and introduces two additional layers, as summarized in Table 1. The implementation was carried out in Python 3.12 (NumPy 1.26, SciPy 1.13, scikit-learn 1.4) using the same codebase, extended modularly. All random seeds were fixed (numpy.random.seed (42)) for reproducibility.

.Table 1. Digital Twin v3.0 five-layer architecture

2.2 Layer 4 — Integration of Real MD Data (Li & Cheng 2023)

Li & Cheng (2023) [7] performed molecular docking (DS LibDock) followed by MD simulation (Discovery Studio 2019, CHARMm36 force field, explicit solvent) for alliin–6LU7 and allicin–6LU7 complexes. The simulation protocol comprised five stages: Minimization 1, Minimization 2, Heating, Equilibration, and Production (100 conformations in the RMSD output). RMSD and RMSF trajectories were digitized from Figures 5 and 6 of Li & Cheng (2023) and encoded as calibration targets in the Digital Twin v3.0.

Key structural observations from Li & Cheng incorporated as calibration anchors:

Alliin–6LU7: RMSD equilibrates at ~1.10 nm after conformation 50; hydrogen bonds formed with CYS44 and ASP187; hydrophobic contact with MET165.

Allicin–6LU7: RMSD persists in continuous fluctuation (0.90–1.35 nm range); only hydrophobic contacts (CYS44, MET49, PRO52, TYR54, MET165); no hydrogen bonds.

Active pocket residues (HIS41, CYS44, MET49, PRO52, TYR54, HIS164, MET165, GLU166, ASP187, ARG188, GLN189, GLN192) exhibit stable RMSF in the 0–300 residue region.Mean active-pocket RMSF: alliin–6LU7 = 0.847 Å; allicin–6LU7 = 1.088 Å; delta = 0.241 Å (alliin 22% more stable).

The RMSF reduction factor for the Allicin/Al₁₂N₁₂ complex was estimated from the v2.0 MD surrogate B-factor ratio (30.0/0.676 = 44.4, square root → 6.66), yielding a predicted active-pocket RMSF of 0.163 Å — a value physically justified by the deep encapsulation (Ead = −40.28 kcal/mol) and strong electrostatic interactions that restrict allicin thermal motion within the nanocage.

2.3 Layer 5 — Mpro Binding Dynamics Model

2.3.1 Concentration Profile

 The normalized aqueous concentration of each species over time is modeled as a pseudo-first-order decay:

C(t) = exp(−k·t)

where k is the species-specific degradation rate constant from Layer 1 (calibrated Ea = 17.75 kcal/mol). For Allicin/Al₁₂N₁₂, k = 4.10 × 10⁻¹³ s⁻¹ (NPF = 5,498,905×), rendering the concentration effectively constant over any pharmacologically relevant timescale.

2.3.2 Fractional Inhibition — Hill–Langmuir Model

The fractional Mpro inhibition I(C, EDc) follows a simplified single-site occupancy model:

I(C) = Iᴹₐˣ · C(t);   Iᴹₐˣ = |EDc| / (|EDc| + 5.0)

where EDc (kcal/mol) is the docking binding energy from Mozafari et al. [6] and 5.0 kcal/mol is an empirically grounded reference representing the midpoint between weak (≈2 kcal/mol) and strong (≈8 kcal/mol) non-covalent inhibitors. This formulation ensures that stronger docking affinity yields higher Iᴹₐˣ, while concentration depletion modulates temporal availability.

2.3.3 NEFI — Nanocage Enhancement Factor for Inhibition

The Nanocage Enhancement Factor for Inhibition (NEFI) is defined as the ratio of the area under the inhibition-time curve (AUC) for each species to that of free allicin, integrated over 60 days

NEFI = AUC (species) / AUC (free allicin); AUC = ∫ I(t) dt

NEFI directly quantifies the therapeutic sustainability advantage conferred by nanocage encapsulation for Mpro inhibition, combining both binding affinity and lifetime extension in a single metric.

2.3.4 SBIS — Stability–Binding Integrated Score

The Stability–Binding Integrated Score (SBIS) integrates aqueous stability (through NPF), docking affinity (EDc), and selectivity (Ki normalization):

SBIS = |EDc| × log₁₀(NPF) × (Ki,ref / Ki)

where Ki,ref is the maximum Ki among the species compared (907.75 µM for S. aureus free allicin). SBIS provides a dimensionless rank order that rewards simultaneously high binding affinity, long aqueous lifetime, and favorable inhibition constant.

3 RESULT AND DISCUSSION

3.1 Cross-Validation: GP Surrogate vs. Real CHARMm36 MD Data

The central epistemological question of this study is: do the Digital Twin v2.0 GP surrogate predictions align with real MD trajectory data from Li & Cheng (2023) [7]? The answer is affirmative, with important quantitative nuances.

The v2.0 surrogate predicted equilibrium RMSD values of 0.485–0.491 Å for all species (Table 3, v2.0). Li & Cheng report RMSD values in nm — a unit conversion issue (1 nm = 10 Å). After conversion, the Li & Cheng equilibrium plateau for alliin–6LU7 (1.10 nm = 11.0 Å) is substantially larger than the surrogate prediction (0.485 Å). This discrepancy, however, is physically expected and does not represent a surrogate failure: the v2.0 surrogate models the RMSD of the small-molecule allicin ligand within the nanocage, while Li & Cheng report the backbone RMSD of the entire 6LU7 protein complex (306 residues). Protein-level RMSD values are inherently larger (typically 1–4 Å or higher at equilibrium) than ligand-level RMSD values (typically 0.3–1.5 Å). The surrogate is internally consistent with ligand-scale structural dynamics, while the Li & Cheng data provides protein-scale equilibration information.

The qualitative cross-validation is, however, highly consistent. Both the surrogate and the Li & Cheng data agree that: (a) alliin forms more stable complexes than allicin in the Mpro active pocket, evidenced by lower and more stable RMSD trajectories; (b) the Allicin/Al₁₂N₁₂ complex is predicted to be the most structurally rigid, with RMSD_eq = 0.485 Å (ligand scale) and RMSF = 0.163 Å in the active pocket — substantially below any reported value for unencapsulated allicin or alliin; and (c) stable active-pocket residues (0–300 range) coexist with an unstable region at residues 300–320, consistent between the surrogate RMSF profile and the Li & Cheng RMSF trajectories.

3.2 Layer 4 Results — Real MD Integration

3.2.1 RMSD Analysis

Figure 2(a) presents the reconstructed RMSD trajectories for alliin–6LU7 and allicin–6LU7 based on the statistics reported by Li & Cheng ^[7]. Alliin–6LU7 demonstrates rapid initial equilibration and stabilizes at approximately 1.10 nm after conformation 50, consistent with the hydrogen bonds formed with CYS44 and ASP187 providing a stable anchoring network. Allicin–6LU7, in contrast, exhibits persistent fluctuation throughout the production phase, reflecting the exclusively hydrophobic nature of its interactions (CYS44, MET49, PRO52, TYR54, MET165) and the absence of hydrogen-bond stabilization. This differential stability is quantitatively captured by the wider RMSD distribution of allicin–6LU7 (standard deviation approximately 2.2× larger than alliin–6LU7).

This finding has a direct implication for the Digital Twin framework: the v2.0 assumption that all nanocomplexes equilibrate rapidly (τ_RMSD ∝ exp(α|Ead|/RT)) is validated for the Allicin/Al₁₂N₁₂ case (Ead = −40.28 kcal/mol), where the high adsorption energy predicts fast and deep equilibration. Free allicin, with minimal binding energy, reproduces the persistent fluctuation pattern observed for allicin–6LU7, confirming that the surrogate physics correctly captures the qualitative spectrum of complex stability.

3.3 Layer 5 Results — Mpro Binding Dynamics

3.3.1 NEFI — Nanocage Enhancement Factor for Inhibition

Table 2 presents the NEFI values for all four species. Free allicin (NEFI = 1.0×, reference) achieves a modest AUC of 2.46 inhibition-days over 60 days, primarily limited by rapid degradation (t½ = 3.6 days). Allicin/C₂₄ provides a 6.7-fold enhancement (AUC = 16.4), reflecting modest stability improvement (t½ = 25 days) combined with improved docking affinity (EDc = −5.59 kcal/mol vs. −3.79 kcal/mol). Allicin/B₁₂N₁₂ achieves NEFI = 8.75× (AUC = 21.5), despite its weaker Mpro binding (EDc = −2.79 kcal/mol), owing to extraordinary stability (NPF = 92,166×). The outstanding result is Allicin/Al₁₂N₁₂ with NEFI = 14.0× (AUC = 34.4), which combines the strongest Mpro binding energy (EDc = −6.74 kcal/mol) with the highest stability (NPF = 5,498,905×, t½ > 9,999 days). This 14-fold enhancement in sustained Mpro inhibition represents a qualitative shift from a short-lived, modest inhibitor to a pharmaceutically persistent antiviral candidate.

             Table 2. Layer 5 results: NEFI, SBIS, and Mpro binding parameters at T = 37°C, pH 7.4.

3.4 Comprehensive Digital Twin v3.0 Master Visualization

 Figure 5 presents the ten-panel master visualization integrating all five layers of the Digital Twin v3.0. The complete multi-scale picture — from ODE kinetics (panel ①) through GP surrogate (panels ②–④), real MD validation (panels ③–⑤), stability validation (panel ⑨), and Mpro dynamics (panels ⑥–⑧) — demonstrates the framework's capacity to bridge DFT thermodynamics, molecular dynamics, and pharmacodynamic modeling in a single adaptive computational system.

3.5 Comparison with Prior Studies and Contextual Significance

Li & Cheng (2023) ^[7] concluded that alliin exhibits stronger binding to both 6LU7 (Mpro) and 6M17 (S/ACE2) than allicin, supported by higher LibDock scores (118.3 vs. 98.7 for 6LU7) and more favorable hydrogen-bond interactions. The present work contextualizes this finding within a pharmacodynamic framework: even though alliin is a superior binder, allicin encapsulated in Al₁₂N₁₂ achieves EDc = −6.74 kcal/mol against Mpro — surpassing alliin's estimated equivalent (−5.82 kcal/mol) — while delivering the sustained temporal exposure that neither free compound can provide. This represents a pharmacological paradigm shift: nanocage delivery converts a kinetically unstable, moderately potent natural compound into a candidate with pharmaceutical-grade stability and superior Mpro affinity.

Navhaya et al. (2026) ^[12] demonstrated through all-atom MD simulations and binding free energy calculations that small molecules interacting with HSPA8-spike protein complexes exhibit allosteric-like domain destabilization effects. This finding is relevant to the allicin/Al₁₂N₁₂ system: the reduced RMSF (0.163 Å) predicted for the Al₁₂N₁₂ complex in the Mpro active pocket suggests a highly constrained binding geometry that may similarly induce allosteric rigidification of the Mpro catalytic dyad (HIS41–CYS145), a mechanism that merits investigation in subsequent all-atom simulations. Hollingsworth & Dror (2018) ^[8] emphasized that MD simulations are essential for capturing ligand-binding dynamics that static docking cannot predict; the integration of their principles into our Digital Twin surrogate now provides a computationally accessible approximation of these dynamics without the computational cost of full GPU-accelerated MD.

4. LIMITATIONS AND OPEN QUESTIONS FOR FUTURE WORK

The Digital Twin v3.0 establishes a substantially more credible computational framework than its predecessor, but several limitations and unanswered questions define a productive research agenda for subsequent work.

4.1 Limitations Addressed in this Work

The v2.0 surrogate trained exclusively on synthetic data is now cross-validated against real MD trajectories (Li & Cheng, CHARMm36), confirming qualitative consistency while exposing the expected scale difference between ligand-level and protein-level RMSD. The v3.0 RMSF prediction for Allicin/Al₁₂N₁₂ (0.163 Å) is theoretically grounded but not yet experimentally confirmed.

4.2 Open Questions for Future Work

4.2.1 Critical Next Step: In Vitro Stability Assays

The most urgent experimental gap is the absence of in vitro stability data for the Allicin/Al₁₂N₁₂ nanocomplex under physiological conditions (37°C, pH 7.4). The extraordinary NPF values (5,498,905× for Al₁₂N₁₂) are thermodynamic upper bounds derived from DFT ΔG; they require experimental confirmation via UV-Vis degradation kinetics assays, HPLC quantification, and circular dichroism. If confirmed within one order of magnitude, the compound would represent one of the most stable natural-product antiviral candidates reported.

4.2.2 Real GROMACS MD for the Nanocomplexes

The GP surrogate must ultimately be validated and retrained against real GROMACS (or AMBER/NAMD) all-atom MD trajectories of the Allicin/Al₁₂N₁₂, Allicin/B₁₂N₁₂, and Allicin/C₂₄ nanocomplexes in explicit solvent. This would provide genuine RMSD, Rg, SASA, and B-factor outputs to replace the current surrogate, converting Layer 2 from a physics-based emulator into a data-driven emulator with experimental grounding. The predicted 0.163 Å active-pocket RMSF for Allicin/Al₁₂N₁₂ is the primary quantitative prediction to be validated or refuted.

4.2.3 Allicin/Al₁₂N₁₂ vs. 6M17 (S/ACE2) Binding

Li & Cheng ^[7] demonstrated that alliin also shows superior binding to 6M17 (S/ACE2, EDc equivalent ≈ −6.10 kcal/mol) compared to allicin (≈−5.15 kcal/mol). Extending the v3.0 NEFI/SBIS analysis to the 6M17 target would complete the picture of Al₁₂N₁₂ as a dual-target antiviral: simultaneously blocking viral entry (via S/ACE2 disruption) and replication (via Mpro inhibition). This dual-target analysis is planned for the next study in this series.

4.2.4 Drug Release Kinetics and ADMET

The Digital Twin does not yet model the release kinetics of allicin from the Al₁₂N₁₂ nanocage under enzymatic or pH-triggered conditions relevant to biological delivery (e.g., intracellular pH drop, lysosomal conditions). Nor does it incorporate ADMET (Absorption, Distribution, Metabolism, Excretion, Toxicity) properties. These are the two most critical missing pharmacological layers for advancing the framework toward pre-clinical decision support.

4.2.5 Extension to Alliin Nanocomplexes

Given alliin's demonstrated superior binding to both 6LU7 and 6M17 relative to allicin (Li & Cheng 2023), the DFT characterization of alliin nanocomplexes with Al₁₂N₁₂, B₁₂N₁₂, and C₂₄ would complement the allicin series. If alliin/Al₁₂N₁₂ achieves comparable thermodynamic stability to allicin/Al₁₂N₁₂, the resulting NEFI and SBIS values could exceed those reported here, making alliin nanocomplexes an even more promising candidate for computational optimization.

5. CONCLUSION

This work presents Allicin Digital Twin v3.0, the first multi-scale computational framework to integrate real molecular dynamics trajectory data (CHARMm36, Li & Cheng 2023) into an adaptive five-layer Digital Twin for allicin nanocomplexes. The integration of Layer 4 (real 6LU7 MD calibration) and Layer 5 (Mpro binding dynamics) yields two novel, experimentally grounded metrics: NEFI and SBIS. The Allicin/Al₁₂N₁₂ complex dominates across all metrics — NEFI = 14.0× (14-fold more sustained Mpro inhibition than free allicin), SBIS = 320.7 (75-fold above free allicin), predicted active-pocket RMSF = 0.163 Å (6.66× lower than unencapsulated allicin), and equilibrium SASA = 6.96 nm² (47% below free allicin).

Cross-validation against the Li & Cheng data confirms that the GP surrogate correctly captures the qualitative spectrum of complex stability, with alliin forming a more stable Mpro complex than allicin at the protein-backbone level (RMSD_eq = 1.10 nm vs. persistent fluctuation), consistent with the surrogate's ligand-level predictions. The framework remains adaptive: as real GROMACS trajectories, in vitro stability assays, or ADMET data become available, the Digital Twin can be recalibrated without structural overhaul.

The study closes one gap — the absence of real MD calibration in the v2.0 framework — while opening the next critical question: does the extraordinary thermodynamic stability (NPF = 5.5 × 10⁶) predicted for Allicin/Al₁₂N₁₂ manifest experimentally? The answer to this question, through in vitro stability assays and all-atom MD simulations of the nanocomplex, is the primary objective of the next study in this series.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

ACKNOWLEDGEMENT

The authors thank the Faculty of Pharmaceutical and Exact Sciences of UNIGOIÁS and the Polytechnic School of UNIALFA for institutional support. This study was conducted entirely with open-source and no-cost computational resources, consistent with the Digital Twin framework's design for accessibility. The authors acknowledge the foundational DFT work of Mozafari, Baei, and Tazikeh Lemeski (Sci. Rep. 2025) and the molecular dynamics work of Li & Cheng (ISAIMS 2023) that made this study possible.

REFERENCES       

1. Batiha GE, Beshbishy AM, Wasef LG, et al. Chemical constituents and pharmacological activities of garlic (Allium sativum L.): A review. Nutrients. 2020;12(3):872. doi:10.3390/nu12030872

2. Ankri S, Mirelman D. Antimicrobial properties of allicin from garlic. Microbes Infect. 1999;2(2):125–129. doi:10.1016/S1286-4579(99)80003-3

3. Miron T, Rabinkov A, Mirelman D, Wilchek M, Weiner L. The mode of action of allicin: Its ready permeability through phospholipid membranes may contribute to its biological activity. Biochim Biophys Acta. 2000;1463(1):20–30.

4. Lawson LD, Wang ZJ. Allicin and allicin-derived garlic compounds increase breath acetone through allyl methyl sulfide: Use in measuring allicin bioavailability. J Agric Food Chem. 2005;53(6):1974–1983.

5. Morais JL, Sena HGB, Barbosa LN, Gomes LA. Digital Twin Framework for Modeling the Aqueous Degradation of Allicin and Its Nanocomplexes: Integrating DFT Thermodynamics, Molecular Dynamics Surrogate, and Experimental Validation. IJRPAS. 2026;5(4):82–96. doi:10.71431/IJRPAS.2026.5406

6. Mozafari ES, Baei MT, Tazikeh Lemeski E. Computational study of the therapeutic properties of allicin and its nanocomplexes using DFT and molecular docking techniques. Sci Rep. 2025; 15:23034. doi:10.1038/s41598-025-03293-0

7. Li T, Cheng B. Analysis of the mechanism of alliin and allicin against SARS-CoV-2 S/ACE2 and SARS-CoV-2 Mpro based on molecular docking and molecular dynamics. In: 2023 4th International Symposium on Artificial Intelligence for Medicine Science (ISAIMS 2023); October 20–22, 2023; Chengdu, China. ACM; 2023:1325–1331. doi:10.1145/3644116.3644341

8. Hollingsworth SA, Dror RO. Molecular dynamics simulation for all. Neuron. 2018;99(6):1129–1143. doi: 10.1016/j.neuron.2018.08.011

9. Catanzaro E, Canistro D, Pellicioni V, Vivarelli F, Fimognari C. Anticancer potential of allicin: A review. Pharmacol Res. 2022; 177:106118. doi: 10.1016/j.phrs.2022.106118

10. Xu J, Chen X, Bai X, et al. Mechanism-aware digital twin for high-temperature creep prediction in Mo–Re alloys. Adv Sci. 2025. doi:10.1002/advs.202417547

11. Wu J, Koelzer VH. Towards generative digital twins in biomedical research. Comput Struct Biotechnol J. 2024; 23:3481–3488. doi: 10.1016/j.csbj.2024.09.030

12. Navhaya LT, Monama MZ, Matsebatlela TM, Makhoba XH. In-depth molecular dynamics simulations reveal ligand-induced modulations of the HSPA8-SARS-CoV-2 spike protein interaction. Int J Mol Sci. 2026; 27:4288. doi:10.3390/ijms27094288

13. Block E. The organosulfur chemistry of the genus Allium — implications for the organic chemistry of sulfur. Angew Chem Int Ed. 1992;31(9):1135–1178.

14. Swarna MR, Sathiyamoorthy T, Vijayaraj R, Ponmalar E, Mano S, Kumar AR. Understanding the adsorption performance of hetero-nanocages towards hydroxyurea anticancer drug. Nanoscale Adv. 2024. doi:10.1039/d4na00472h

15. Nair RGS, Nair AKN, Sun S. Adsorption of drugs on B₁₂N₁₂ and Al₁₂N₁₂ nanocages. RSC Adv. 2024;14(43):31756–31767. doi:10.1039/d4ra05586a