In silico exploring the mechanisms of action of diterpenoids from Rabdosia serra against lung cancer through inhibition of the anti-apoptotic pathway

Hung Duc Nguyen

Thai Nguyen University of Education, Thai Nguyen, Viet Nam

Corresponding author: Hung Duc Nguyen (hungnd@tnue.edu.vn)

Abstract

Introduction: Apoptosis resistance in non-small cell lung cancer is frequently sustained by pro-survival Bcl-2 family proteins such as Mcl-1, motivating the search for new Mcl-1 inhibitors from natural products, including diterpenoids from Rabdosia serra.

Materials and Methods: An integrated in silico approach was applied to evaluate R. serra diterpenoids as putative Mcl-1 (PDB: 6QFQ) inhibitors, integrating molecular docking, molecular dynamics simulation, MM/GBSA rescoring, pkCSM-based ADMET prediction, and DFT calculations at the B3LYP/6-31G(d,p) level, with Tivantinib as the reference ligand.

Results and Discussion: Docking prioritized CPD1 (-10.31 kcal/mol) over Tivantinib (-9.09 kcal/mol). Molecular dynamics simulation indicated stable complexes, with CPD1-6QFQ showing a more compact, less solvent-exposed ensemble (Rg 1.42-1.45 nm; SASA 83-89 nm²; RMSD 0.09-0.15 nm). MM/GBSA favored CPD1 (ΔTOTAL -26.17 ± 2.79 kcal/mol) versus Tivantinib (-24.25 ± 4.50 kcal/mol). ADMET predicted high intestinal absorption for CPD1 (98.243%), a higher unbound fraction, and fewer liabilities (negative hepatotoxicity; negative hERG II inhibition). DFT supported CPD1 with a smaller ΔE (3.7165 eV) and higher softness (0.5381 eV¹) than Tivantinib.

Conclusion: Convergent computational evidence nominates CPD1 as a leading putative Mcl-1 inhibitor for optimization, while target engagement, efficacy, and safety require experimental validation.

Graphical Abstract

Keywords: anti-apoptosis; diterpenoid; in silico; lung cancer; Mcl-1; Rabdosia serra

Introduction

Cancer imposed a substantial mortality burden in 2022, with approximately 9.7 million deaths worldwide, corresponding to roughly one in six deaths from all causes. Major determinants include modifiable exposures such as tobacco use, alcohol consumption, unhealthy dietary patterns, physical inactivity, and air pollution (Bray et al. 2024). Lung cancer remains a dominant contributor to cancer mortality worldwide, with global estimates for 2022 indicating about 1.8 million deaths (Zhou et al. 2024). Clinicopathological classification typically distinguishes non-small cell lung cancer as the predominant category, comprising about 80% to 85% of cases, and small cell lung cancer as a less frequent but more aggressive entity, comprising about 10% to 15% of cases (Petersen and Warth 2016). Accordingly, mechanistic and therapeutic discussions often focus on non-small cell lung cancer.

Current treatment for lung cancer includes surgery, radiotherapy, chemotherapy, targeted therapy, and immunotherapy. Early-stage cancer is commonly treated with surgery, whereas advanced disease typically requires systemic treatment, with targeted agents used when actionable alterations are present and immunotherapy considered in appropriate clinical contexts (Bouchard and Daaboul 2025; Tahayneh et al. 2025). Despite these advances, resistance remains a significant limitation, and impaired apoptosis is a key survival mechanism in non-small cell lung cancer. This phenotype is frequently linked to dysregulation of the apoptosis regulatory circuitry, including the functional dominance of anti-apoptotic Bcl-2 family members that support the survival of genomically compromised cells and reduce responsiveness to cytotoxic therapy (Chattopadhyay et al. 2023; Kim et al. 2025). Among pro-survival Bcl-2 proteins, Mcl-1 has emerged as a compelling target in apoptosis evasion, and pharmacologic disruption of Mcl-1-mediated survival signaling is therefore of sustained interest as a strategy to restore apoptotic competence and mitigate treatment resistance (Tantawy et al. 2023; Chen et al. 2025).

Natural products remain a central reservoir for drug discovery because of high scaffold diversity and a long record of contribution to approved therapeutics, including oncology indications. Plant-derived metabolites have been widely evaluated as candidate agents or adjuncts in non-small cell lung cancer through effects on proliferation control, apoptosis induction, and modulation of oncogenic signaling and resistance-associated pathways, although the clinical evidence base remains heterogeneous and further trials are commonly recommended (Shah et al. 2013; Asma et al. 2022). Diterpenoids constitute a structurally diverse class of plant-derived secondary metabolites that have attracted sustained attention in oncology-oriented natural product research (Zhou et al. 2012). Within traditional Chinese medicine, Rabdosia serra is used primarily for hepatobiliary and gastrointestinal disorders, including hepatitis, jaundice, cholecystitis, and intestinal ailments. Collectively, the chemical profile of R. serra provides a rationale for continued investigation of its diterpenoid-rich constituents as candidates for mechanistic and translational evaluation in cancer-related research (Liu et al. 2022). A previous study on this species led to the isolation of several diterpenoids, and some of them possess biological activity against the NCI-H661 non-small cell lung cancer cell line (Wang and Xuan 2016). However, the specific molecular targets and pathways mediating the anticancer activities of these diterpenoids remain insufficiently characterized. Accordingly, an integrated in silico strategy using molecular docking, molecular dynamics simulation, MMGBSA, ADMET prediction, and DFT descriptors can provide mechanistic plausibility and candidate prioritization for subsequent experimental validation and structure-guided optimization. At the same time, safety-related interpretation requires confirmation in dedicated in vitro and in vivo systems.

Materials and Methods

Data collection

These selected diterpenoids, including serrin B (CPD1), serrin A (CPD2), isodocarpin (CPD3), and lushanrubescensin J (CPD4), have molecular formulas C20H26O7, C22H30O6, C20H26O5, and C40H52O12, respectively, with molecular weights 378.1679, 390.2042, 346.1780, and 724.3459 m/z. Tivantinib, possessing a molecular formula C23H19N3O13 and a molecular weight 369.1477 m/z, was chosen as the positive control (Fig. 1).

Figure 1. 3D Structures of selected ligands. Note: (A) CPD1, (B) CPD2, (C) CPD3, (D) CPD4, (E) Tivantinib.

Molecular docking analysis

Three-dimensional models of the selected diterpenoids and the reference compound Tivantinib were constructed in PDB format using BIOVIA Discovery Studio Visualizer. Ligand preparation comprised the addition of polar hydrogen atoms, assignment of Gasteiger partial charges, and preservation of conformational flexibility by retaining torsional degrees of freedom for rotatable bonds. The crystallographic structure of the anti-apoptotic protein Mcl-1 (PDB ID: 6QFQ) was obtained from the RCSB Protein Data Bank (Murray et al. 2019). Molecular docking was executed in AutoDock Tools using a cubic grid encompassing the binding region, defined by 60 points along each Cartesian axis (x, y, z) with a spacing of 0.375 Å. The grid center was set at x = 3.934 Å, y = -18.343 Å, and z = 18.680 Å. Pose generation and ranking were performed with the Lamarckian genetic algorithm to identify low-energy binding orientations and recurrent interaction features. The top-ranked protein-ligand complex was subsequently inspected in Discovery Studio Client 2024, and its interaction pattern was evaluated alongside the docking pose of Tivantinib within the same Mcl-1 structure to characterize shared and distinct binding determinants.

Molecular dynamics simulations

Molecular dynamics simulations were performed to characterize the time-dependent stability of the docked Mcl-1 ligand complex, using the Mcl-1 crystal structure (PDB ID 6QFQ) as the structural template. All simulations were carried out in GROMACS 2024.4 with the CHARMM36 force field and a trajectory length of 100 ns (van der Spoel et al. 2005). Protein preparation included completion of missing atoms and residues in Swiss-PdbViewer (Guex and Peitsch 1997). Ligand topology and parameterization compatible with CHARMM were generated with SwissParam (Zoete et al. 2011). The complex was embedded in a triclinic periodic simulation box, solvated using SPC water, and ionized with NaCl to 0.15 M while maintaining overall electroneutrality. Energy minimization was conducted using the steepest descent algorithm for 50000 steps to remove steric conflicts and unfavorable contacts. Equilibration was completed in two stages, consisting of 200 ps in the NVT ensemble at 300 K, followed by 200 ps in the NPT ensemble at 1 bar. The production phase comprised three independent simulations, each spanning 100 ns with a 2 fs integration time step, and coordinate snapshots were saved every 10 ns for downstream analysis. Trajectory metrics were processed in Grace, including root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent-accessible surface area (SASA), and number of hydrogen bonds (Hbonds). Structural similarity among representative conformations was further examined in UCSF Chimera v1.13.3 through coordinate superposition to compare characteristic states across the simulations (Pettersen et al. 2004).

Molecular mechanics generalised born surface area (MM/GBSA) analysis

Binding free energy estimates for the CPD1-6QFQ and Tivantinib-6QFQ complexes were obtained using the gmx_MMPBSA workflow with the CHARMM36 parameter set. The polar solvation contribution was evaluated with a generalized Born implicit solvent model, whereas the nonpolar solvation term was derived from solvent accessible surface area. Energetic components were extracted from the molecular dynamics trajectories by processing 125 uniformly distributed snapshots sampled at 80 ps intervals over an 80 ns analysis window spanning 20 to 100 ns. The resulting snapshot averaged energies enabled direct comparison of ligand protein interaction energetics and supported interpretation of relative binding propensity and complex persistence within the simulated timeframe.

Assay protocol for ADMET prediction

Early assessment of absorption, distribution, metabolism, excretion, and toxicity (ADMET) characteristics provides critical insight into potential pharmacokinetic limitations and safety liabilities, thereby informing compound prioritization and reducing the likelihood of late-stage failure. In the present analysis, the pkCSM web server was used to estimate the ADMET profiles of CPD1 and Tivantinib computationally. pkCSM applies graph-based molecular signature representations to infer predictors associated with absorption, distribution, metabolism, excretion, and toxicity, enabling a structured comparison of developability-relevant attributes across the evaluated compounds (Pires et al. 2015).

Quantum chemistry computation using the density functional theory (DFT) method

Geometry optimization of CPD1 and Tivantinib was performed with the ORCA quantum chemistry suite version 6.1.0. Starting conformations were constructed in Avogadro, whereas molecular orbital visualization and related electronic examinations were conducted in IboView version 20211019 (http://www.iboview.org). Density functional theory calculations employed the B3LYP exchange correlation functional in combination with the 6-31G(d,p) basis set to locate optimized minima and obtain the corresponding electronic wavefunctions. Frontier orbital energies for the HOMO and LUMO levels were subsequently extracted from the optimized structures, and the associated energy separation ΔE was calculated. Conceptual DFT descriptors, including chemical potential (µ), electronegativity (χ), global hardness (η), softness (σ), and electrophilicity index (ω), were then computed using a Koopmans-type approximation to support the interpretation of electronic structure characteristics and comparative reactivity tendencies of the investigated molecules (Hanwell et al. 2012; Knizia and Klein 2015; Neese 2025).

Results and Discussion

Molecular docking analysis

Molecular docking estimates plausible ligand poses inside a protein cavity and ranks these poses using scoring functions that approximate noncovalent stabilization. In contemporary structure-based workflows, interpretation typically combines the energetic order with recurrence of pocket residues and the distribution of interaction types, because docking scores alone can be sensitive to approximations in solvation, entropy, and protein flexibility (Nivatya et al. 2025). Based on the reported binding energies in Table 1, CPD1 exhibited the most favorable score at -10.31 kcal/mol, followed by CPD3 at -9.32 kcal/mol and CPD2 at -8.64 kcal/mol. In contrast, CPD4 and the reference Tivantinib produced values of -6.24 and -9.09 kcal/mol, respectively, which is consistent with less favorable binding under the scoring convention applied in this docking run. Because docking energies are best treated as relative rankings within the same protocol, subsequent refinement using dynamic sampling and more physics-grounded free-energy estimators remains appropriate for mechanistic claims.

Table 1.

Download as

_XLSX_

_CSV_

The interactions between the docked ligands and the protein 6QFQ

Docked ligands

Binding energy (kcal/mol)

Hydrogen bond interaction

van der Waals interaction

Hydrophobic interaction

CPD1

-10.31

Val249, Phe254, Arg263, Phe270

Phe228, Met231

Met250, Val253, Leu267, Phe270

CPD2

-8.64

Arg263, Thr266

Val249, Phe254

Phe228, Met231, Met250, Val253, Leu267, Phe270

CPD3

-9.32

Leu246, Leu267

Ser247, Val249, Phe254, Thr266, Gly271

Met231, Leu235, Leu246, Met250, Val253, Leu267, Phe270, Val274, Leu290, Ile294

CPD4

-6.24

Arg263, Thr266

Val220, Phe254, Leu267, Phe270

His224, Ala227, Phe228, Met231, Val249, Met250, Val253, Arg263

Tivantinib

-9.09

Phe254, Leu267

Ala227, Leu235, Leu246, Val253, Arg263, Gly271

Phe228, Met231, Val249, Met250, Thr266, Leu267, Phe270, Val274, Leu290, Ile294

CPD1 engaged a compact active site defined by Phe228, Met231, Val249, Met250, Val253, Arg263, Leu267, and Phe270. Hydrogen bonding was assigned to Val249, Phe254, Arg263, and Phe270, indicating polar anchoring distributed across both aromatic and charged environments at the pocket interface. Van der Waals contributions were reported for Phe228 and Met231, consistent with close-range packing against the pocket core. Hydrophobic contacts involved Met250, Val253, Leu267, and Phe270, suggesting that nonpolar burial around residue 250, residue 253, and the Phe270 ring contributes substantially to pose stabilization (Fig. 2A). The combination of directional hydrogen bonding with extensive apolar complementarity matches a common recognition pattern in small molecule binding where dispersion and hydrophobic effects support affinity once a viable polar anchor is established (Ferenczy and Kellermayer 2022; Adediwura et al. 2024).

For concision and emphasis, CPD2-CPD4 can be considered a secondary group with less favorable energetic behavior and more heterogeneous anchoring patterns. CPD2 shared the identical principal active site residues as CPD1, namely Phe228, Met231, Val249, Met250, Val253, Arg263, Leu267, and Phe270. Still, it exhibited a lower docking score and fewer interaction annotations, with hydrogen bonds assigned to Arg263 and Thr266, and van der Waals contacts limited to Val249 and Phe254. CPD3 retained a similar hydrophobic nucleus centered on Met231, Val249, Met250, Val253, Leu267, and Phe270, yet its hydrogen bonding annotation shifted to Leu246 and Leu267, and its hydrophobic interactions expanded toward Val274, Leu290, and Ile294, suggesting sampling of an adjacent apolar subregion and a pose geometry distinct from CPD1 (Sun 2022). CPD4 showed the broadest residue coverage, including Val220, His224, and Ala227 in addition to the shared pocket residues, and it maintained hydrogen bond assignments to Arg263 and Thr266. However, the positive docking energy for CPD4 indicates an unfavorable score under this protocol despite extensive nominal contacts, which may reflect suboptimal shape complementarity or an unfavorable balance between desolvation and intramolecular strain captured by the scoring function. Such divergence between contact counts and docking energy is compatible with known limitations of scoring functions and with the role of interaction geometry and solvation treatment in docking performance.

Figure 2. Molecular docking model and 2D interaction diagram of CPD1 (A) and Tivantinib (B) with 6QFQ protein.

Tivantinib occupied an overlapping region with active site residues Ala227, Phe228, Met231, Val249, Met250, Val253, Leu267, Phe270, and Gly271. Hydrogen bonds were assigned to Phe254 and Leu267, while van der Waals interactions involved Ala227, Leu235, Leu246, Val253, Arg263, and Gly271. Hydrophobic interactions extended across Phe228, Met231, Val249, Met250, Thr266, Leu267, Phe270, Val274, Leu290, and Ile294, indicating extensive apolar complementarity (Fig. 2B). Nevertheless, the positive docking score suggests that, within this docking configuration, the reference compound was not stabilized as effectively as CPD1-CPD3, emphasizing that contact breadth does not guarantee a favorable energy estimate when pose the scoring function considers strain and solvation terms (Zhang et al. 2022; Sahu et al. 2024).

Overall, repeated engagement of Phe228, Met231, Val249, Met250, Val253, Leu267, and Phe270 across ligands indicates a conserved recognition core enriched in aromatic and aliphatic side chains, where hydrophobic effects and dispersion interactions are expected to be prominent contributors to stabilization. Within that framework, CPD1 distinguished itself by the best docking energy and by a balanced profile of hydrogen bonding and hydrophobic packing centered on Arg263, Phe254, and Phe270 together with the recurrent nonpolar core residues. Based on the docking rank order and interaction recurrence within the predicted pocket, CPD1 is prioritized for dynamic evaluation, while Tivantinib is retained as a reference ligand for protocol-matched comparison. To account for receptor plasticity and solvent-mediated rearrangements that are not represented in static docking, all-atom molecular dynamics simulations will be conducted for the CPD1-6QFQ and Tivantinib-6QFQ complexes, using the experimental structure deposited as PDB 6QFQ as the starting receptor model.

Molecular dynamics simulations

Molecular dynamics simulations provide time-resolved atomic trajectories that complement docking by introducing explicit conformational sampling under a defined force field and solvent representation. This approach is routinely used to examine whether an initial docking pose remains stable, to identify persistent interaction motifs, and to quantify structural responses of the receptor to ligand binding through trajectory descriptors including RMSD, RMSF, Rg, SASA, and Hbonds (Uludağ and Tang 2021). The CPD1-6QFQ complex exhibited total and potential energies of -221,657 kJ/mol and -275,625 kJ/mol, respectively, whereas the corresponding values for Tivantinib-6QFQ were -221,462 kJ/mol and -275,381 kJ/mol. Temperature regulation maintained the simulations at 300 K.

RMSD value is a critical analysis tool used to quantify the structural deviation of a molecule from a reference structure over the course of a simulation. In the backbone RMSD, both systems remained within a narrow deviation regime over 0 to 100 ns. CPD1-6QFQ predominantly populated 0.09 to 0.15 nm, with intermittent excursions that reached about 0.17 nm in the second half of the trajectory. Tivantinib-6QFQ occupied 0.07 to 0.16 nm for most frames, with short-lived peaks close to 0.18 nm, most evident after the mid-trajectory region (Fig. 3A). The traces converged near 0.09 to 0.11 nm during the early segment, whereas the interval from roughly 50 to 90 ns showed more frequent sampling of higher RMSD values for Tivantinib than for CPD1, indicating a larger-amplitude global response for the reference complex under the exact fitting definition (Sargsyan et al. 2017).

Figure 3. Results of MD simulation for the bindings of CPD1 (blue) and Tivantinib (red) with 6QFQ protein. Note: (A) RMSD, (B) RMSF, (C) Rg, (D) SASA, (E) Hbonds.

To account for residue-wise flexibility, the analysis was restricted to Asp172-Val321, which encompasses the docking-defined recognition corridor and includes key pocket residues such as His224, Ala227, Phe228, Met231, Val249-Val253, Arg263, Thr266, Leu267, and Phe270. Within this segment, both complexes exhibited low-amplitude fluctuations, with RMSF values largely confined to 0.04-0.07 nm and with substantial overlap between the CPD1-6QFQ and Tivantinib-6QFQ traces (Fig. 3B). No pronounced local maxima were evident in this interval. Fluctuations did not exceed 0.08 nm based on the plotted scale, supporting preservation of a relatively rigid binding-site environment during the 100 ns simulation. RMSF corresponds to the standard deviation of positional fluctuations after optional fitting and is widely used to localize flexible segments while distinguishing them from comparatively restrained regions (Salo-Ahen et al. 2021).

Compactness, as assessed by Rg, indicated a more compact ensemble for CPD1-6QFQ across much of the simulation. CPD1-6QFQ generally occupied 1.42 to 1.45 nm, with intermittent increases that reached about 1.47 nm. Tivantinib-6QFQ more often populated a higher band, commonly 1.44 to 1.46 nm, and displayed repeated excursions toward 1.47 to 1.48 nm, particularly between about 55 and 85 ns. A decrease toward 1.42 to 1.44 nm occurred near the terminal portion of the trajectory for Tivantinib-6QFQ (Fig. 3C). The radius of gyration in GROMACS is mass-weighted and is frequently interpreted as a global measure of compactness, where sustained higher values indicate a less compact ensemble (Sneha and George Priya Doss 2016).

The SASA implementation reports solvent-accessible areas along the trajectory and is commonly used to evaluate changes in exposure associated with conformational expansion or contraction. Surface exposure trends were consistent with the Rg behavior. SASA for CPD1-6QFQ was distributed mainly in the 83 to 89 nm² interval, with occasional departures that reached about 81 nm² at the lower end and near 92 nm² at the upper end. Tivantinib-6QFQ remained shifted upward for extended periods, commonly occupying 86 to 92 nm², with spikes reaching near 96 nm² and late decreases that approached 80 to 85 nm² (Fig. 3D). The average separation between the traces during the early to mid trajectory visually corresponded to several nm², indicating higher solvent exposure for the Tivantinib-bound complex in that period.

Hbonds are defined using a geometric criterion and provide a time-resolved proxy for polar-contact persistence, although hydrogen-bond counts do not necessarily map monotonically onto global stability metrics. CPD1-6QFQ exhibited frequent hydrogen-bond occurrences during the early segment, with many frames at one hydrogen bond and repeated short spikes to two or three between roughly 10 and 25 ns. A single event at four hydrogen bonds was visible at the start of the time axis. After this early phase, the CPD1 trace was dominated by zero to one hydrogen bond with sporadic short increases (Fig. 3E). Tivantinib-6QFQ showed sparse hydrogen bonding through much of the mid-trajectory window, typically zero to one, but displayed a late-stage cluster from about 90 to 100 ns characterized by sustained one to two hydrogen bonds and multiple spikes to three.

Collectively, CPD1-6QFQ showed lower Rg and lower SASA for substantial portions of the simulation, accompanied by a slightly lower RMSD envelope across the mid to late segment, supporting a more compact and less solvent-exposed conformational ensemble relative to Tivantinib-6QFQ. Tivantinib-6QFQ exhibited broader excursions in RMSD, higher Rg and SASA for extended intervals, and a delayed increase in hydrogen-bond formation during the late trajectory, consistent with greater conformational plasticity and a later emergence of stronger polar anchoring. For energetic refinement beyond structural descriptors, endpoint rescoring methods, such as MM/PBSA or MM/GBSA, applied to equilibrated snapshots are widely used, provided that protocol sensitivity and reporting constraints are rigorously addressed.

Free binding energy (MM/GBSA) analysis

MM/GBSA estimates binding using an end-state strategy that evaluates molecular-mechanics interaction energies for the complex in the gas phase and adds an implicit-solvent correction. In this formulation, polar solvation is treated with a generalized Born model, whereas the nonpolar contribution is commonly approximated from solvent-accessible surface area. Because conformational entropy is often excluded due to expense and inconsistent performance across protocols, MM/GBSA results are generally interpreted as relative indicators of binding strength when sampling and model parameters are kept identical across systems (Poli et al. 2020). As summarized in Table 2, both complexes exhibited favorable mean total binding free energies, with CPD1-6QFQ showing ΔTOTAL = -26.17 kcal/mol with SD 2.79 kcal/mol, and Tivantinib-6QFQ showing ΔTOTAL = -24.25 kcal/mol with SD 4.50 kcal/mol. The more negative mean for CPD1 indicates a stronger predicted association under the present MM/GBSA protocol, while the larger dispersion for Tivantinib suggests greater energetic heterogeneity across sampled frames. Such heterogeneity can reflect broader variability in interfacial packing and electrostatic alignment during the trajectory (Liu et al. 2023).

Table 2.

Download as

_XLSX_

_CSV_

Free energy of binding obtained using MMGBSA calculations

Energy Component

Average (kcal/mol)

Standard Deviation

CPD1-6QFQ

Tivantinib-6QFQ

CPD1-6QFQ

Tivantinib-6QFQ

ΔVDWAΔALS

-31.43

-39.74

2.41

6.20

ΔEEL

-3.05

-9.10

4.42

5.56

ΔEGB

12.09

29.68

4.47

4.44

ΔESURF

-3.78

-5.10

0.27

0.63

ΔGGAS

-34.48

-48.84

5.54

6.66

ΔGSOLV

8.31

24.58

4.40

4.32

ΔTOTAL

-26.17

-24.25

2.79

4.50

Energy decomposition indicates distinct balances between gas-phase stabilization and solvation penalty. For CPD1-6QFQ, the gas-phase term was favorable at ΔGGAS = -34.48 kcal/mol with SD 5.54 kcal/mol, dominated by van der Waals attraction ΔVDWAALS = -31.43 kcal/mol with SD 2.41 kcal/mol, while Coulombic electrostatics were weaker at ΔEEL = -3.05 kcal/mol with SD 4.42 kcal/mol. A comparatively small net solvation opposed this interaction gain cost ΔGSOLV = 8.31 kcal/mol with SD 4.40 kcal/mol, driven mainly by the polar component ΔEGB = 12.09 kcal/mol with SD 4.47 kcal/mol and partially compensated by the favorable nonpolar surface term ΔESURF = -3.78 kcal/mol with SD 0.27 kcal/mol. This signature is consistent with a binding mode dominated by close-range packing and dispersion complementarity, with limited electrostatic contribution and a moderate polar desolvation burden.

In contrast, Tivantinib-6QFQ showed substantially stronger gas-phase stabilization at ΔGGAS = -48.84 kcal/mol with SD 6.66 kcal/mol, arising from both enhanced dispersion ΔVDWAALS = -39.74 kcal/mol with SD 6.20 kcal/mol and more favorable electrostatics ΔEEL = -9.10 kcal/mol with SD 5.56 kcal/mol. However, this stronger vacuum interaction was counterbalanced by a markedly larger solvation penalty ΔGSOLV = 24.58 kcal/mol with SD 4.32 kcal/mol, dominated by the polar GB term ΔEGB = 29.68 kcal/mol with SD 4.44 kcal/mol with an additional compensatory nonpolar term ΔESURF = -5.10 kcal/mol with SD 0.63 kcal/mol. The near cancellation between more attractive gas-phase interactions and increased polar desolvation rationalizes why ΔTOTAL for Tivantinib is less favorable than for CPD1, despite a much more negative ΔGGAS. This compensation pattern is a recognized feature of end-state implicit-solvent approaches, in which a higher desolvation cost accompanies greater electrostatic complementarity. Variability analysis further supports CPD1 as the more energetically consistent complex under the sampled ensemble. Tivantinib presented larger SD values for key interaction terms, particularly ΔVDWAALS (SD 6.20 vs 2.41) and ΔEEL (SD 5.56 vs 4.42), consistent with broader fluctuations in short-range packing and electrostatic alignment. In addition, ΔTOTAL dispersion for Tivantinib (SD 4.50) exceeded that of CPD1 (SD 2.79), indicating less stable energetic convergence across frames within the same MM/GBSA settings. Sensitivity of MM/GBSA results to sampling length, dielectric choices, and solvation model details has been emphasized in recent benchmarking and methodological discussions, supporting cautious interpretation of fine rank differences that are comparable to the reported dispersions (Tuccinardi 2021).

Overall, the MM/GBSA decomposition indicates that CPD1-6QFQ achieves a more favorable mean ΔTOTAL through a balance of strong dispersion stabilization and a relatively small solvation penalty, whereas Tivantinib-6QFQ exhibits stronger gas-phase attraction that is substantially offset by a larger polar solvation cost, together with higher frame-to-frame variability. Under the present end-state protocol, these results support CPD1 as the preferable complex for subsequent interpretation and reporting, while maintaining Tivantinib as a mechanistically informative comparator because of its distinct electrostatics-desolvation trade-off.

ADMET prediction analysis

In silico ADMET profiling is commonly applied at the earliest stages of drug discovery to screen developability attributes before committing to resource-intensive in vitro and in vivo testing. Data-driven predictors, including QSAR and machine-learning models, enable rapid triage of chemical series by highlighting likely liabilities in exposure or safety, thereby supporting early deprioritization of candidates with unfavorable pharmacokinetic or toxicological signals and improving the efficiency of subsequent synthesis and experimental validation (Toh et al. 2025). Table 3 summarizes the predicted ADMET properties for CPD1 and the reference Tivantinib.

Table 3.

Download as

_XLSX_

_CSV_

Predicted ADMET properties of CPD1 and Tivantinib

ADMET properties

Unit

CPD1

Tivantinib

Water Solubility

(Log mol/L)

-4.034

-4.405

Caco2 permeability

(Log Papp in 10-6 cm/s)

1.172

1.122

Intestinal absorption (Human)

(% Absorbed)

98.243

94.761

Skin permeability

(Log Kp)

-3.322

-3.221

P-glycoprotein substrate

Yes/No

Yes

Yes

P-glycoprotein I inhibitor

Yes/No

Yes

Yes

P-glycoprotein II inhibitor

Yes/No

No

Yes

VDss

(Log L/kg)

0.201

-0.243

Fraction unbound (human)

(Fu)

0.308

0.046

BBB permeability

(Log BB)

-0.105

0.013

CNS permeability

(Log PS)

-3.156

-1.741

CYP2D6 substrate

Yes/No

No

No

CYP3A4 substrate

Yes/No

Yes

Yes

CYP1A2 inhibitor

Yes/No

No

Yes

CYP2C19 inhibitor

Yes/No

No

Yes

CYP2C9 inhibitor

Yes/No

No

No

CYP2D6 inhibitor

Yes/No

No

No

CYP3A4 inhibitor

Yes/No

No

Yes

Total clearance

(Log mL/min/kg)

0.329

0.458

Renal OCT2 substrate

Yes/No

No

No

AMES toxicity

Yes/No

No

No

Max. tolerated dose (human)

(Log mg/kg/day)

-0.565

0.578

hERG I inhibitor

Yes/No

No

No

hERG II inhibitor

Yes/No

No

Yes

Oral rat acute toxicity (LD50)

(mol/kg)

2.167

2.774

Oral rat chronic toxicity (LOAEL)

(Log mg/kg_bw/day)

2.421

2.22

Hepatotoxicity

Yes/No

No

Yes

Skin sensation

Yes/No

No

No

Tetrahymena pyriformis toxicity

(Log ug/L)

0.286

0.462

Minnow toxicity

(Log mM)

3.398

0.65

CPD1 displayed a higher predicted aqueous solubility than Tivantinib, with water solubility values of -4.034 log mol/L and -4.405 log mol/L, respectively. Epithelial transport surrogates were also slightly higher for CPD1, with Caco-2 permeability of 1.172 versus 1.122 log Papp at 10⁻⁶ cm/s, and higher predicted human intestinal absorption of 98.243% versus 94.761%. Caco-2 permeability is widely used as an intestinal transcellular permeability surrogate, although assay design, transporter expression, and model domain can influence absolute interpretability (Steinbauer et al. 2024). Skin permeability values were close in magnitude, with CPD1 showing -3.322 log Kp and Tivantinib -3.221 log Kp, indicating a slightly lower modeled transdermal permeation tendency for CPD1. Both ligands were classified as P-glycoprotein substrates and P-glycoprotein I inhibitors, while only Tivantinib was predicted as a P-glycoprotein II inhibitor. This combined substrate and inhibitor pattern is relevant because P-glycoprotein can influence oral exposure through efflux in the intestine and can contribute to transporter-mediated drug-drug interaction risk when inhibition co-occurs (Veiga-Matos et al. 2023).

Distribution descriptors favored broader systemic distribution and higher free fraction for CPD1. The predicted steady-state volume of distribution was 0.201 log L/kg for CPD1 compared with -0.243 log L/kg for Tivantinib. The human fraction unbound was markedly higher for CPD1 at Fu 0.308, whereas Tivantinib showed Fu 0.046, consistent with a larger circulating unbound pool for CPD1 under this estimator. Central nervous system indices separated the two ligands in the opposite direction. BBB permeability was -0.105 log BB for CPD1 and 0.013 log BB for Tivantinib, and CNS permeability was -3.156 log PS for CPD1 versus -1.741 log PS for Tivantinib. These values indicate higher modeled CNS exposure surrogates for Tivantinib, while also emphasizing that BBB metrics should be interpreted as comparative indicators, as transporter effects and model assumptions can strongly influence actual brain penetration (Nabi et al. 2025).

Both compounds were predicted as CYP3A4 substrates and non-substrates for CYP2D6, indicating a shared reliance on CYP3A4-associated oxidative metabolism in this prediction set. A divergence emerged in inhibition alerts. CPD1 showed no predicted inhibition for CYP1A2, CYP2C19, or CYP3A4, while Tivantinib was expected to be an inhibitor of CYP1A2, CYP2C19, and CYP3A4. Because CYP3A4 is a major contributor to clinical clearance and interaction liability, the presence of a CYP3A4 inhibition flag for Tivantinib suggests higher perpetrator interaction potential within the limitations of classifier-based predictions (Damoiseaux et al. 2024).

Predicted total clearance was higher for Tivantinib (0.458 log mL/min/kg) than for CPD1 (0.329 log mL/min/kg), consistent with faster modeled systemic elimination for the reference compound. Both ligands were classified as non-substrates for renal OCT2, which reduces but does not eliminate concern about OCT2-mediated renal secretion variability in this output set. OCT2 and related MATE transporters are recognized determinants of renal secretion and clinically relevant transporter-based interactions for susceptible substrates and inhibitors (Ailabouni and Prasad 2025; Asano et al. 2025).

Both ligands were predicted as AMES negative, indicating no modeled bacterial mutagenicity signal in this workflow, while acknowledging that mutagenicity prediction accuracy remains domain-dependent and benefits from experimental confirmation for lead candidates (Uesawa 2024). Cardiac liability flags differed. Both compounds were predicted as hERG I non-inhibitors. In contrast, hERG II inhibition was negative for CPD1 and positive for Tivantinib, consistent with a higher modeled cardiotoxicity concern for the reference compound under the applied classifier. Contemporary in silico hERG screening emphasizes applicability-domain control and orthogonal validation because false negatives and endpoint heterogeneity remain possible. Hepatotoxicity was predicted as No for CPD1 and Yes for Tivantinib, aligning CPD1 with a lower modeled liver-toxicity liability in this panel, while recognizing that hepatotoxicity prediction is a complex endpoint influenced by multiple mechanisms and model assumptions. Additional toxicity descriptors showed quantitative separation. The maximum tolerated dose was -0.565 log mg/kg/day for CPD1 and 0.578 log mg/kg/day for Tivantinib. Acute oral rat LD50 was 2.167 mol/kg for CPD1 and 2.774 mol/kg for Tivantinib, while chronic oral rat LOAEL was 2.421 log mg/kg bw/day for CPD1 and 2.22 log mg/kg bw/day for Tivantinib. Environmental toxicity proxies also differed, with T. pyriformis toxicity 0.286 log µg/L for CPD1 versus 0.462 log µg/L for Tivantinib, and minnow toxicity 3.398 log mM for CPD1 versus 0.65 log mM for Tivantinib, indicating non-uniform organism-dependent trends across endpoints and units (Shin et al. 2023; Falcón-Cano et al. 2025).

In conclusion, CPD1 combined higher predicted solubility and slightly stronger intestinal transport surrogates with a higher unbound fraction and higher predicted volume of distribution. In contrast, Tivantinib showed higher modeled BBB and CNS permeability indices, higher predicted clearance, and multiple metabolic and safety flags, including CYP inhibition, hERG II inhibition, and hepatotoxicity. From a comparative developability perspective, the CPD1 profile suggests more favorable systemic exposure characteristics and fewer predicted liabilities in the metabolism and toxicity panels. In contrast, Tivantinib shows stronger interaction and safety signals that may warrant closer experimental scrutiny upon progression.

Quantum chemistry computation using the DFT method

Electronic-structure characterization was performed using density functional theory to quantify frontier orbital levels and to derive global reactivity indices for CPD1 and the benchmark compound Tivantinib (Table 4). The positions of the highest occupied and lowest unoccupied molecular orbitals provide a compact description of how readily a molecule can participate in charge transfer processes. At the same time, the HOMO-LUMO gap summarizes the energetic cost of promoting electron density into an acceptor manifold. In conceptual DFT, these quantities are routinely interpreted together with chemical potential, electronegativity, hardness, softness, and electrophilicity to frame molecular responsiveness to polarization and external perturbation comparatively (Miranda-Quintana et al. 2022).

Table 4.

Download as

_XLSX_

_CSV_

Quantum descriptors of CPD1 and Tivantinib

Molecule

EHOMO

(eV)

ELUMO

(eV)

ΔE

(eV)

µ

(eV)

χ

(eV)

η

(eV)

σ

(eV-1)

ω

(eV)

CPD1

-7.5695

-3.8530

3.7165

-5.7113

5.7113

1.8583

0.5381

8.7766

Tivantinib

-7.4231

2.9976

10.4207

-2.2128

2.2128

5.2104

0.1919

0.4699

Note: EHOMO (eV) – highest occupied molecular orbitals; ELUMO (eV) –- lowest unoccupied molecular orbitals; ΔE (eV) – energy gap; µ (eV) – chemical potential; χ (eV) – electronegativity; η (eV) – hardness; σ (eV-1) – softness; ω (eV) – electrophilicity index.

Tivantinib has a slightly higher HOMO energy (EHOMO = -7.4231 eV) than CPD1 (EHOMO = -7.5695 eV), consistent with a modestly greater donor tendency for Tivantinib under the same computational protocol. In contrast, CPD1 exhibits a substantially lower LUMO energy (ELUMO = -3.8530 eV) than Tivantinib (2.9976 eV), indicating that CPD1 provides a more accessible acceptor orbital. These orbital trends are reflected in the ordering of the gaps. CPD1 presents a markedly smaller ΔE (3.7165 eV) than Tivantinib (10.4207 eV), suggesting greater electronic pliability for CPD1. At the same time, Tivantinib exhibits a more rigid electronic structure, commonly associated with greater kinetic stability in the conceptual DFT interpretation (Fig. 4).

Figure 4. HOMO and LUMO surface diagrams of CPD1 (A) and Tivantinib (B).

Global descriptors further differentiate the two electronic profiles. CPD1 shows a more negative chemical potential (µ = -5.7113 eV) than Tivantinib (-2.2128 eV), accompanied by a higher electronegativity for CPD1 (χ = 5.7113 eV) relative to Tivantinib (2.2128 eV). This combination is consistent with a stronger overall electron-attracting tendency for CPD1 at the global level. Hardness and softness follow the same direction as the orbital gap. CPD1 exhibits lower hardness (η = 1.8583 eV) and higher softness (σ = 0.5381 eV¹) compared with Tivantinib (η = 5.2104 eV, σ = 0.1919 eV¹), indicating that electron density deformation is energetically more facile for CPD1. In contrast, Tivantinib resists density redistribution more strongly (Geerlings 2022).

The electrophilicity index ω provides an additional global measure of charge acceptance propensity that depends on both µ and η. CPD1 exhibits a high electrophilicity value (ω = 8.7766 eV), while Tivantinib shows a much smaller value (ω = 0.4699 eV). Under the same descriptor definitions, this contrast indicates that CPD1 has a substantially stronger electrophilic signature than Tivantinib in the global conceptual DFT sense, consistent with the combination of a very low LUMO, small ΔE, and high softness observed for CPD1. Conversely, Tivantinib combines a larger gap and higher hardness with a comparatively weak electrophilicity index, which aligns with reduced global susceptibility toward charge acceptance despite the slightly higher HOMO energy.

Conclusion

An integrated in silico approach was applied to evaluate diterpenoids from R. serra as putative Mcl-1 (6QFQ) inhibitors, integrating molecular docking, molecular dynamics, MM/GBSA rescoring, pkCSM-based ADMET prediction, and DFT calculations at the B3LYP 6-31G(d,p) level with Tivantinib as the reference ligand. Docking prioritized CPD1 with a binding energy of -10.31 kcal/mol, whereas Tivantinib yielded -9.09 kcal/mol under the same scoring convention. Molecular dynamics indicated stable complexes overall, with CPD1-6QFQ showing a more compact and less solvent-exposed ensemble, evidenced by Rg of 1.42 to 1.45 nm, SASA of 83 to 89 nm², and RMSD primarily within 0.09 to 0.15 nm. At the same time, Tivantinib sampled higher compactness and exposure metrics and exhibited a late increase in hydrogen bonding. MM/GBSA supported a favorable association for both ligands and favored CPD1 by the mean total binding free energy, ΔTOTAL = -26.17 ± 2.79 kcal/mol versus -24.25 ± 4.50 kcal/mol, consistent with a binding profile dominated by van der Waals stabilization and a minor solvation penalty for CPD1. ADMET estimates suggested high intestinal absorption for CPD1 at 98.243%, a higher predicted unbound fraction than Tivantinib, and fewer predicted liabilities, including negative hepatotoxicity and negative hERG II inhibition. DFT descriptors further indicated greater electronic adaptability for CPD1, with a smaller ΔE of 3.7165 eV and a higher softness of 0.5381 eV¹ than Tivantinib. Collectively, convergence across structure-based, dynamics-based, energetics-based, ADMET, and electronic structure analyses supports CPD1 as the leading candidate for subsequent optimization and experimental validation.

Computational predictions provide a rational basis for candidate selection but cannot substitute for experimental demonstration of target engagement or biological efficacy. Because no biochemical, cellular, or animal studies were performed, the inferred binding strength, dynamic persistence, and pharmacokinetic or safety implications remain provisional and should not be interpreted as validated therapeutic effects. Follow-up work should quantify Mcl-1 binding using direct biophysical assays, evaluate apoptosis-related outcomes in lung cancer models with documented Mcl-1 dependence, and determine exposure and clearance in pharmacokinetic studies to address developability constraints suggested by limited solubility and transporter-associated behavior.

 

Additional Information

Conflict of interest

The authors declare that they have no conflicts of interest.

Funding

No funding was received to conduct this research.

Acknowledgments

The authors have no support to report.

Data availability

All of the data that support the findings of this study are available in the main text.

References

§  Adediwura VA, Koirala K, Do HN, Wang J, Miao Y (2024) Understanding the impact of binding free energy and kinetics calculations in modern drug discovery. Expert Opinion on Drug Discovery 19(6): 671–682. https://doi.org/10.1080/17460441.2024.2349149 [PubMed] [PMC]

§  Ailabouni A, Prasad B (2025) Organic cation transporter 2: Structure, regulation, functions, and clinical implications. Drug Metabolism and Disposition 53(3): 100044. https://doi.org/10.1016/j.dmd.2025.100044 [PubMed]

§  Asano S, Galetin A, Tomita Y, Giacomini KM, Chu X, Yang X, Nakamura T, Kusuhara H, Sugiyama Y (2025) Predicting OCT2/MATEs-mediated drug interactions in healthy volunteers and patients with chronic kidney disease. Clinical Pharmacology & Therapeutics 118(5): 994–1014. https://doi.org/10.1002/cpt.3727 [PubMed] [PMC]

§  Asma ST, Acaroz U, Imre K, Morar A, Shah SR, Hussain SZ, Arslan-Acaroz D, Demirbas H, Hajrulai-Musliu Z, Istanbullugil FR, Soleimanzadeh A, Morozov D, Zhu K, Herman V, Ayad A, Athanassiou C, Ince S (2022) Natural products as a source of anticancer drugs. Cancers 14(24): 6203. https://doi.org/10.3390/cancers14246203 [PubMed] [PMC]

§  Bouchard N, Daaboul N (2025) Lung cancer: Targeted therapy in 2025. Current Oncology 32(3): 146. https://doi.org/10.3390/curroncol32030146 [PubMed] [PMC]

§  Bray F, Laversanne M, Sung H, Ferlay J, Siegel RL, Soerjomataram I, Jemal A (2024) Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. A Cancer Journal for Clinicians 74(3): 229–263. https://doi.org/10.3322/caac.21834 [PubMed]

§  Chattopadhyay S, Sarkar SS, Saproo S, Yadav S, Antil D, Das B, Naidu S (2023) Apoptosis-targeted gene therapy for non-small cell lung cancer. Frontiers in Bioengineering and Biotechnology 11: 1188652. https://doi.org/10.3389/fbioe.2023.1188652 [PubMed] [PMC]

§  Chen H, Yan Z, Xiao C, Tang Z, Wan Y (2025) Advances of small-molecule Mcl-1 inhibitors. European Journal of Medicinal Chemistry 299: 118079. https://doi.org/10.1016/j.ejmech.2025.118079 [PubMed]

§  Damoiseaux D, Beijnen JH, Huitema AD, Dorlo TPC (2024) Early prediction and impact assessment of CYP3A4-related drug-drug interactions for small-molecule anticancer drugs using human-CYP3A4-transgenic mouse models. Drug Metabolism and Disposition 52(11): 1217–1223. https://doi.org/10.1124/dmd.123.001530 [PubMed]

§  Falcón-Cano G, Morales-Helguera A, Lambert H, Cabrera-Pérez MÁ, Molina C (2025) hERG toxicity prediction using machine learning. Scientific Reports 15(1): 15585. https://doi.org/10.1038/s41598-025-99766-3

§  Ferenczy GG, Kellermayer M (2022) Computational and Structural Biotechnology Journal 20: 1946–1956. https://doi.org/10.1016/j.csbj.2022.04.025

§  Geerlings P (2022) From density functional theory to biosystems. Pharmaceuticals 15(9): 1112. https://doi.org/10.3390/ph15091112 [PubMed] [PMC]

§  Guex N, Peitsch MC (1997) SWISS-MODEL and Swiss-Pdb Viewer. Electrophoresis 18(15): 2714–2723. https://doi.org/10.1002/elps.1150181505

§  Hanwell MD, Curtis DE, Lonie DC, Vandermeersch T, Zurek E, Hutchison GR (2012) Avogadro: An advanced semantic chemical editor, visualization, and analysis platform. Journal of Cheminformatics 4(1): 17. https://doi.org/10.1186/1758-2946-4-17 [PubMed] [PMC]

§  Kim H, Jung SH, Jo S, Han JW, Lee JH (2025) Induction of apoptotic cell death in non-small-cell lung cancer cells by MP28 peptide derived from bryopsis plumosa. Marine Drugs 23(12): 481. https://doi.org/10.3390/md23120481 [PubMed] [PMC]

§  Knizia G, Klein JEMN (2015) Electron flow in reaction mechanisms - revealed from first principles. Angewandte Chemie International Edition 54(18): 5518–5522. https://doi.org/10.1002/anie.201410637 [PubMed]

§  Liu F, Zeng Y, Dai P, Huang K, Zhang K, Tao T, Wang M, Zhu C, Lin C (2022) Comparative Pharmacokinetics of Three bioactive diterpenoids of rabdosia serra extract in normal and con a-induced liver injury rats using UPLC-MS/M. Frontiers in Pharmacology 13: 944949. https://doi.org/10.3389/fphar.2022.944949 [PubMed] [PMC]

§  Liu X, Zheng L, Qin C, Cong Y, Zhang JZH, Sun Z (2023) comprehensive evaluation of end-point free energy techniques in carboxylated-pillar[6]arene host-guest binding: III. force-field comparison, three-trajectory realization and further dielectric augmentations. Molecules 28(6): 2767. https://doi.org/10.3390/molecules28062767 [PubMed] [PMC]

§  Miranda-Quintana RA, Heidar-Zadeh F, Fias S, Chapman AEA, Liu S, Morell C, Gómez T, Cárdenas C, Ayers PW (2022) Molecular interactions from the density functional theory for chemical reactivity: Interaction chemical potential, hardness, and reactivity principles. Frontiers in Chemistry 10: 929464. https://doi.org/10.3389/fchem.2022.929464 [PubMed] [PMC]

§  Murray JB, Davidson J, Chen I, Davis B, Dokurno P, Graham CJ, Harris R, Jordan A, Matassova N, Pedder C, Ray S, Roughley SD, Smith J, Walmsley C, Wang Y, Whitehead N, Williamson DS, Casara P, Le Diguarher T, Hickman J, Stark J, Kotschy A, Geneste O, Hubbard RE (2019) Establishing drug discovery and identification of hit series for the anti-apoptotic proteins, Bcl-2 and Mcl-1. ACS Omega 4(5): 8892–8906. https://doi.org/10.1021/acsomega.9b00611 [PubMed] [PMC]

§  Nabi AE, Pouladvand P, Liu L, Hua Y, Ayubcha C (2025) Machine learning in drug development for neurological diseases: A review of blood brain barrier permeability prediction models. Molecular Informatics 44(3): e202400325. https://doi.org/10.1002/minf.202400325 [PubMed] [PMC]

§  Neese F (2025) Software update: the ORCA program system version 6.0. WIREs Computational Molecular Science 15(2): e70019. https://doi.org/10.1002/wcms.70019

§  Nivatya HK, Singh A, Kumar N, Sonam, Sharma L, Singh V, Mishra R, Gaur N, Mishra AK (2025) Assessing molecular docking tools: Understanding drug discovery and design. Future Journal of Pharmaceutical Sciences 11(1): 111. https://doi.org/10.1186/s43094-025-00862-y

§  Petersen I, Warth A (2016) Lung cancer: Developments, concepts, and specific aspects of the new WHO classification. Journal of Cancer Research and Clinical Oncology 142(5): 895–904. https://doi.org/10.1007/s00432-015-2004-4 [PubMed] [PMC]

§  Pettersen EF, Goddard TD, Huang CC, Couch GS, Greenblatt DM, Meng EC, Ferrin TE (2004) UCSF chimera-a visualization system for exploratory research and analysis. Journal of Computational Chemistry 25(13): 1605–1612. https://doi.org/10.1002/jcc.20084 [PubMed]

§  Pires DE, Blundell TL, Ascher DB (2015) pkCSM: Predicting small-molecule pharmacokinetic and toxicity properties using graph-based signatures. Journal of Medicinal Chemistry 58(9): 4066–4072. https://doi.org/10.1021/acs.jmedchem.5b00104 [PubMed] [PMC]

§  Poli G, Granchi C, Rizzolio F, Tuccinardi T (2020) Application of MM-PBSA methods in virtual screening. Molecules 25(8): 1971. https://doi.org/10.3390/molecules25081971 [PubMed] [PMC]

§  Sahu D, Rathor LS, Dwivedi SD, Shah K, Chauhan NS, Singh MR, Singh D (2024) A review on molecular docking as an interpretative tool for molecular targets in disease management. ASSAY and Drug Development Technologies 22(1): 40–50. https://doi.org/10.1089/adt.2023.060 [PubMed]

§  Salo-Ahen OM, Alanko I, Bhadane R, Bonvin AMJJ, Honorato R V, Hossain S, Juffer AH, Kabedev A, Lahtela-Kakkonen M, Larsen AS, Lescrinier E (2021) Molecular dynamics simulations in drug discovery and pharmaceutical development. Processes 9(1): 71. https://doi.org/10.3390/pr9010071

§  Sargsyan K, Grauffel C, Lim C (2017) How molecular size impacts RMSD applications in molecular dynamics simulations. Journal of Chemical Theory and Computation 13(4): 1518–1524. https://doi.org/10.1021/acs.jctc.7b00028 [PubMed]

§  Shah U, Shah R, Acharya S, Acharya N (2013) Novel anticancer agents from plant sources. Chinese Journal of Natural Medicines 11(1): 16–23. https://doi.org/10.1016/S1875-5364(13)60002-3

§  Shin HK, Huang R, Chen M (2023) In silico modeling-based new alternative methods to predict drug and herb-induced liver injury: A review. Food and Chemical Toxicology 179113948. https://doi.org/10.1016/j.fct.2023.113948 [PubMed] [PMC]

§  Sneha P, George Priya Doss C (2016) Molecular dynamics: New frontier in personalized medicine. In: Donev R (Ed) Personalized medicine. Academic Press, Boston, pp. 181-224. https://doi.org/10.1016/bs.apcsb.2015.09.004

§  Steinbauer FF, Lehr T, Reichel A (2024) Exploring the potential of adaptive, local machine learning in comparison to the prediction performance of global models: A case study from Bayer’s Caco-2 permeability database. Journal of Chemical Information and Modeling 64(24): 9163–9172. https://doi.org/10.1021/acs.jcim.4c01083 [PubMed] [PMC]

§  Sun Q (2022) The hydrophobic effects: Our current understanding. Molecules 27(20): 7009. https://doi.org/10.3390/molecules27207009 [PubMed] [PMC]

§  Tahayneh K, Idkedek M, Abu Akar F (2025) NSCLC: Current evidence on its pathogenesis, integrated treatment, and future perspectives. Journal of Clinical Medicine 14(3): 1025. https://doi.org/10.3390/jcm14031025 [PubMed] [PMC]

§  Tantawy SI, Timofeeva N, Sarkar A, Gandhi V (2023) Targeting Mcl-1 protein to treat cancer: opportunities and challenges. Frontiers in Oncology 13: 1226289. https://doi.org/10.3389/fonc.2023.1226289 [PubMed] [PMC]

§  Toh WX, Tang SY, Sim KS, Lim SH, Tan CH (2025) Ethnobotanical and alkaloid composition with their cytotoxicity of Tabernaemontana polyneura: A review with in silico ADMET analysis. Fitoterapia 186: 106851. https://doi.org/10.1016/j.fitote.2025.106851 [PubMed]

§  Tuccinardi T (2021) What is the current value of MM/PBSA and MM/GBSA methods in drug discovery? Expert Opinion on Drug Discovery 16(11): 233–1237. https://doi.org/10.1080/17460441.2021.1942836 [PubMed]

§  Uesawa Y (2024) Progress in predicting Ames test outcomes from chemical structures: An in-depth re-evaluation of models from the 1st and 2nd Ames/QSAR international challenge projects. International Journal of Molecular Sciences 25(23): 1373. https://doi.org/10.3390/ijms25031373 [PubMed] [PMC]

§  Uludağ H, Tang T (2021) How can molecular dynamics simulations assist with gene medicines? Biomaterials and Biosystems 2: 100014. https://doi.org/10.1016/j.bbiosy.2021.100014 [PubMed] [PMC]

§  Van Der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE, Berendsen HJC (2005) GROMACS: Fast, flexible, and free. Journal of Computational Chemistry 26(16): 1701–1718. https://doi.org/10.1002/jcc.20291 [PubMed]

§  Veiga-Matos J, Morales AI, Prieto M, Remião F, Silva R (2023) Study models of drug-drug interactions involving p-glycoprotein: The potential benefit of p-glycoprotein modulation at the kidney and intestinal levels. Molecules 28(22): 7532. https://doi.org/10.3390/molecules28227532 [PubMed] [PMC]

§  Wang W, Xuan L (2016) ent-6,7-Secokaurane diterpenoids from Rabdosia serra and their cytotoxic activities. Phytochemistry 122: 119–125. https://doi.org/10.1016/j.phytochem.2015.11.014 [PubMed]

§  Zhang B, Li H, Yu K, Jin Z (2022) Molecular docking-based computational platform for high-throughput virtual screening. CCF Transactions on High Performance Computing 4(1): 63–74. https://doi.org/10.1007/s42514-021-00086-5 [PubMed] [PMC]

§  Zhou J, Xu Y, Liu J, Feng L, Yu J, Chen D (2024) Global burden of lung cancer in 2022 and projections to 2050: Incidence and mortality estimates from GLOBOCAN. Cancer Epidemiology 93: 102693. https://doi.org/10.1016/j.canep.2024.102693 [PubMed]

§  Zhou L, Wang J, Wang K, Xu J, Zhao J, Shan T, Luo C (2012) Secondary metabolites with antinematodal activity from higher plants. In: Atta-ur-Rahman (Ed) Studies in natural products chemistry. Elsevier, Amsterdam, pp. 67–114. https://doi.org/10.1016/B978-0-444-59514-0.00003-1

§  Zoete V, Cuendet MA, Grosdidier A, Michielin O (2011) SwissParam: A fast force field generation tool for small organic molecules. Journal of Computational Chemistry 32(11): 2359–2368. https://doi.org/10.1002/jcc.21816 [PubMed]

Author Contribution

§  Hung Duc Nguyen, Dr (Pharmaceutical Biology), Associate Professor at the Faculty of Biology, Thai Nguyen University of Education, Viet Nam; e-mail: hungnd@tnue.edu.vn; ORCID ID: https://orcid.org/0000-0002-5764-1242. The sole responsibility for the conception of the study, presented results, and manuscript preparation.