Sustainable leads for hypoxic tumor therapy: allosteric selective inhibition of carbonic anhydrase IX by abietane-type resin acids

Abstract

Carbonic anhydrase IX (CA IX) helps hypoxic tumors maintain the acidic extracellular milieu that drives invasion and therapy resistance, making it an attractive but still under-exploited drug target. Here we report that four abietane-type resin acids—callitrisic, levopimaric, palustric and pimaric acids—act as potent, non-competitive and CA IX-selective inhibitors. Virtual screening of 2400 natural products, followed by AutoDock Vina, positioned the rigid diterpenoid cores in a hydrophobic cleft that flanks—but does not overlap—the catalytic zinc. Stopped-flow CO₂-hydrase assays confirmed sub-micromolar potency, with IC₅₀ values of 93.4 ± 1.7 nM (callitrisic), 147.6 ± 4.7 nM (levopimaric), 290.5 ± 15.7 nM (palustric) and 340.7 ± 6.4 nM (pimaric); acetazolamide gave 44 ± 1.7 nM under identical conditions. Lineweaver–Burk and Michaelis–Menten plots revealed an allosteric, non-competitive mechanism, while counterscreens against housekeeping hCA I and hCA II yielded selectivity indices of 5–15, underscoring a therapeutic bias toward the tumor isoform. Two-hundred-nanosecond molecular-dynamics simulations produced single-basin free-energy landscapes dominated by van-der-Waals contacts with Phe243, Ala245, Pro248 and Ala249; MM/PBSA ranked binding free energies in the same order as enzymatic potency, with callitrisic acid reaching −20.58 ± 1.86 kcal mol⁻¹. Global MD metrics (RMSD, Rg, SASA, RMSF) demonstrated that ligand binding leaves the CA IX scaffold intact. In silico ADMET profiling predicted high oral and blood–brain permeability, low hERG liability and manageable CYP2C9/2C19 inhibition, highlighting log P moderation as the main optimization vector. Together, these findings establish abietane resin acids—particularly callitrisic acid—as sustainable, tractable leads for selectively disabling CA IX and normalizing tumor pH, paving the way for further optimization and in vivo validation.

---

1. Introduction

Carbonic anhydrase IX (CA IX) is a membrane-bound, zinc-dependent metallo-enzyme that catalyzes the reversible hydration of carbon dioxide.¹ Its expression is tightly regulated by hypoxia-inducible factor-1α and is therefore markedly up-regulated in the acidic, poorly vascularized regions of many solid tumors.² By accelerating the extrusion of protons and bicarbonate, CA IX helps malignant cells maintain an alkaline intracellular pH while acidifying the extracellular milieu—conditions that favor proliferation, invasion, and resistance to chemotherapy and radiotherapy.³ Accordingly, selective inhibition of CA IX has emerged as a compelling strategy to normalize tumor pH homeostasis and thwart cancer progression without disrupting the ubiquitous cytosolic isoform CA II.^4,5

In an initial in silico campaign we screened 2400 natural products from an in-house curated library assembled from the ZINC “Natural Products” subset (release 2023-01) and commercial NP collections. The workflow (Fig. 1) combined (i) a pharmacophore filter capturing a hydrophobic tricyclic core plus a distal hydrogen-bond/ionic acceptor, (ii) AutoDock Vina docking into the CA IX allosteric pocket defined around residues Phe243–Ala245–Pro248–Ala249 (exhaustiveness = 16), and (iii) MM/GBSA rescoring of the top 10% poses. After PAINS/reactivity filtering and visual inspection, callitrisic, levopimaric, palustric and pimaric acids consistently ranked within the top percentile. They (a) orient their carboxylate away from the catalytic Zn²⁺, engaging the distal pocket rather than chelating the metal, (b) present a rigid abietane/pimarane tricyclic core complementary to the pocket contour, and (c) display favorable physicochemical properties (c log P ≈ 4–5; TPSA < 40 Ų) consistent with passive uptake into hypoxic tumor tissue. Their ready availability from sustainable conifer oleoresins, low reported cytotoxicity toward non-malignant cells, and documented anti-inflammatory activity further justified their selection for experimental validation.

Fig. 1 Structure-based virtual-screening funnel and top-ranked diterpenoid hits.

To obtain a comprehensive picture of their inhibitory potential, we combined orthogonal experimental and computational techniques.⁶ A stopped-flow CO₂ hydration assay quantified apparent inhibition constants (K_i), while Michaelis–Menten kinetics elucidated the mode of inhibition.⁵ In parallel, molecular docking furnished initial binding poses; 200 ns all-atom molecular-dynamics (MD) simulations refined those complexes and generated free-energy landscapes (FEL) that revealed dominant conformational states. Binding free energies were further decomposed by MM/PBSA, and in silico ADMET filters assessed drug-likeness and safety liabilities.^7–10 This integrative workflow allowed us to correlate experimental potency with predicted stability and pharmacokinetic profiles.

The present study therefore aims to (i) validate callitrisic, levopimaric, palustric, and pimaric acids as selective CA IX inhibitors; (ii) delineate their binding mechanisms at atomic resolution; and (iii) identify structural determinants that can be leveraged for future lead optimization toward potent, tumor-targeted carbonic anhydrase modulators.

2. Materials and methods

2.1. Target-selection analysis

The most probable macromolecular targets of callitrisic, levopimaric, palustric and pimaric acids were predicted with SwissTargetPrediction (v2023.4, http://www.swisstargetprediction.ch/; accessed 15 January 2025).¹¹ Canonical SMILES strings for each compound were generated in ChemDraw 23.0, energy-minimized with MMFF94 in Avogadro, and submitted individually to the web interface with the organism set to Homo sapiens. SwissTargetPrediction compares both 2-D fingerprints (ECFP4) and 3-D shape/pharmacophore features (Electroshape) against its curated database of ligands with experimentally confirmed targets. For every acid, the tool returned the top-15 targets ranked by probability; these were automatically grouped by SwissTargetPrediction into functional classes (enzymes, nuclear receptors, GPCRs, etc.). Probability weights for each class were summed and visualized as pie charts (Fig. 2). Only targets with cumulative probabilities ≥1% were considered for subsequent analysis.

Fig. 2 Predicted target-class distributions for callitrisic, levopimaric, palustric, and pimaric acids (SwissTargetPrediction). Percentages indicate the fraction of predicted targets in each protein family. Enzyme-related classes are present for all ligands but are not uniformly dominant; notable non-enzyme families include nuclear receptors, GPCRs, ion channels, oxidoreductases, and fatty-acid-binding proteins.

2.2. In vitro evaluation of carbonic-anhydrase inhibition

2.2.1. Chemicals and reagents. Recombinant human carbonic anhydrase IX catalytic domain (hCA IX; residues 59–414, C-terminal His₁₀ tag), hCA I (Ala²–Phe²⁶¹) and hCA II (Ser²–Lys²⁶⁰) were purchased from R&D Systems (Minneapolis, MN, USA) at >95% purity. 4-Nitrophenyl acetate (4-NPA, ≥98%), bromothymol blue, HEPES, tris, sodium sulfate, dimethyl-sulfoxide (DMSO), acetazolamide (AAZ) and all other analytical-grade chemicals came from Sigma-Aldrich/Merck. Ultrapure water (18.2 MΩ cm) was used throughout.

2.2.2 Primary inhibition assay (IC₅₀ determination). Esterase activity was quantified in flat-bottom 96-well plates (Greiner, polystyrene, 200 μL per well) at 25 °C. Reaction mixtures contained 50 mM HEPES/Tris buffer, pH 7.5, 0.10 M Na₂SO₄, 0.50 mM 4-NPA (≈2 × K_m) and 25 nM hCA IX (2 μg mL⁻¹). Resin-acid inhibitors or AAZ were added over 0.001–1.0 μM (final DMSO ≤ 1% v/v) and pre-incubated 15 min. Hydrolysis of 4-NPA to 4-nitrophenolate was monitored at 400 nm in a Multiskan SkyHigh plate reader; the effective optical path length was 0.60 cm (manufacturer specification for 200 μL). Initial velocities (0–120 s) were fitted to a four-parameter logistic equation in GraphPad Prism 10.0 to obtain IC₅₀ values. Each point represents the mean ± SD of three independent experiments, each performed in technical triplicate. Identical conditions were applied to hCA I and hCA II for selectivity-index calculations.

2.2.3. Enzyme-kinetics analysis. Steady-state CO₂-hydrase kinetics were measured on an Applied Photophysics SX20 stopped-flow spectrophotometer using the Wilbur–Anderson pH-indicator method. Equal volumes (50 μL) of CO₂-saturated HEPES/Tris buffer (final CO₂ 0.10–2.00 mM) and enzyme/inhibitor solution (25 nM hCA IX, inhibitor 0–0.50 μM) were mixed at 25 °C; 50 μM bromothymol blue served as the pH indicator. The absorbance change at 620 nm was recorded for 10 s, and initial rates were extracted from the first 1–3 s. Data were globally fitted to the mixed-inhibition equation to obtain apparent K_i and α values, while Lineweaver–Burk and Michaelis–Menten plots were generated for graphical confirmation of the inhibition mode. All kinetic experiments were performed in triplicate; buffer blanks and vehicle controls (1% DMSO) showed no measurable spontaneous CO₂ hydration.

2.3. Computational set-up

2.3.1. Enzyme preparation. The X-ray structure of the human carbonic-anhydrase IX catalytic domain (PDB ID 5FL4, 1.82 Å) was downloaded from the Protein Data Bank. Missing side-chains and unresolved loop atoms were rebuilt with Swiss-PdbViewer 4.1; the resulting model was subjected to a fast-relax protocol (200 steps steepest-descent, AMBER94 force field) to remove local steric clashes. All crystallographic waters, buffer ions and the co-crystallised ligand were deleted in UCSF Chimera 1.17 unless they participated in the first hydration sphere of the catalytic Zn²⁺. Polar hydrogens were added explicitly and non-polar hydrogens were merged according to the AutoDockTools protocol; Kollman partial charges were then assigned to the receptor atoms. The Zn²⁺ ion was retained with a +2 formal charge and typed as “Zn” in AutoDock4. The prepared protein was saved as a PDBQT file for docking.

2.3.2. Ligand preparation. Two-dimensional structures of callitrisic, levopimaric, palustric and pimaric acids were retrieved from PubChem and converted to 3-D using UCSF Chimera 1.17.¹² Each ligand was geometry-optimized by density-functional theory (B3LYP/6-31G(d,p), gas phase) in Gaussian 16; vibrational-frequency analysis confirmed the absence of imaginary modes.^13,14 Optimized coordinates were imported into ADT, where Gasteiger partial charges were applied and all rotatable bonds were left flexible. Final ligand files were exported as PDBQT.

2.4. Molecular-docking protocol

Docking simulations were carried out with AutoDock Vina 1.2.0.^15,16 The grid box was centered on the centroid of the native ligand (x = 14.755, y = −26.181, z = 55.061 Å) and sized 34 × 34 × 34 Å to encompass the catalytic cleft and the adjacent hydrophobic pocket identified in preliminary mapping. Exhaustiveness was set to 16; the number of output poses was held at 20 with an energy difference cut-off of 3 kcal mol⁻¹. For each ligand the top-ranked pose (lowest predicted binding free energy) was selected and visually inspected in Chimera and Discovery Studio Visualizer to verify correct orientation relative to key residues (Phe243, Ala245, Pro248, Ala249, Gln92, Arg130) and to ensure the carboxylate remained >5 Å from the Zn²⁺ center, in line with the non-competitive inhibition mechanism. All subsequent molecular-dynamics and MM/PBSA calculations used these Vina-derived complexes as starting coordinates.

Docking-protocol validation was performed in two steps. First, the sulfonamide inhibitor co-crystallized in PDB 5FL4 (“native ligand”) was stripped from the protein, converted to PDBQT, and re-docked with the same grid and Vina settings used for the resin acids. The top-ranked pose superimposed on the crystallographic conformation with an all-atom RMSD of 1.4 Å, well below the 2.0 Å threshold generally accepted for reliable redocking, and reproduced the canonical tetrahedral coordination of the deprotonated sulfonamide nitrogen to the catalytic Zn²⁺ along with hydrogen bonds to Thr200 and Thr201. Second, the reference inhibitor (AAZ)—absent from the crystal—was docked as a positive control; Vina returned a binding energy of −8.7 kcal mol⁻¹ and an orientation identical to that reported for AAZ in hCA II, i.e., zinc-anchored sulfonamide flanked by hydrophobic contacts with Val121, Leu198 and Phe243. The accurate recapture of both the native-ligand pose and the expected AAZ binding mode confirms that the grid definition and scoring parameters are suitable for predicting ligand interactions within this CA IX allosteric/catalytic region.

2.5. Molecular-dynamics (MD) simulations

The top-ranked docking poses of callitrisic, levopimaric, palustric and pimaric acids in CA IX (PDB ID 5FL4) and the apo enzyme were propagated by 200 ns all-atom MD in GROMACS 2022.4.^17,18 Protein atoms were described with the CHARMM36m force field, while ligand parameters were generated through the CGenFF server and merged into the master topology.¹⁹ Each system was placed in a dodecahedral box (≈289 nm³), solvated with the CHARMM-modified TIP3P water model and neutralized with seven Na⁺ counter-ions. Energy minimization (steepest descent, ≤1000 kJ mol⁻¹ nm⁻¹) removed bad contacts, after which the system was equilibrated for 100 ps in the NVT ensemble (V-rescale thermostat, 300 K) followed by 100 ps in the NPT ensemble (Parrinello–Rahman barostat, 1 bar, 300 K). All covalent bonds involving hydrogens were constrained with LINCS, allowing a 2 fs integration step; particle-mesh Ewald handled long-range electrostatics with a 1.2 nm real-space cut-off, and the same cut-off was applied to Lennard-Jones interactions. Production runs of 200 ns were then executed for each system, with coordinates saved every 10 ps for subsequent RMSD, RMSF, SASA, Rg, hydrogen-bond, principal-component, FEL and MM/PBSA analyses.

2.5.1. Free-energy-landscape (FEL) construction. Cartesian coordinates from the 200 ns of each production trajectory were stripped of rotational/translational drift (gmx trjconv, −fit rot + trans) and concatenated at 10 ps intervals (20 000 frames).^20,21 Mass-weighted covariance matrices of Cα atomic fluctuations were generated with gmx covar and diagonalized to obtain eigenvectors and eigenvalues. The two principal components (PC1, PC2) capturing ≥60% of the total variance were used as collective variables. A two-dimensional free-energy surface was computed with gmx sham (bin size = 0.2 nm) according to ΔG = −RT ln P(i), where P(i) is the normalized probability of bin i. 2-D contour maps and 3-D surface plots were rendered in Matplotlib; the basin centroid of each FEL was extracted for representative-structure analysis.

2.5.2. MM/PBSA binding-energy calculations. Binding free energies were estimated over the final 50 ns (150–200 ns) of each 200 ns trajectory using gmx_MMPBSA v1.4.1 with default CHARMM36m radii.²² Frames were sampled every 20 ps (2500 snapshots per system). The polar solvation term was calculated with the Poisson–Boltzmann solver (internal dielectric = 2, external = 80; grid spacing = 0.5 Å). The non-polar solvation contribution was obtained from solvent-accessible surface area (SASA) using:

ΔG_non-polar = γ × SASA + β

where γ = 0.0072 kcal mol⁻¹ Å⁻² is the surface-tension coefficient and β = 0.00 kcal mol⁻¹ is an empirical offset. ΔE_vdW, ΔE_ele, ΔG_polar and ΔG_non-polar were integrated to give ΔG_total; per-residue decompositions employed the same frame set.

2.6. In silico ADMET profiling

Canonical SMILES of the four resin acids were submitted to SwissADME (http://www.swissadme.ch/, accessed 12 Feb 2025) to obtain physicochemical descriptors, Lipinski filters, GI absorption and BBB permeation predictions.²³ The same structures were uploaded to ADMETlab 3.0 (https://admetlab3.scbdd.com/) for extended pharmacokinetic and toxicity endpoints, including CYP inhibition likelihoods, hERG blockade, Ames mutagenicity and organ-specific toxicity probabilities. Results were exported as CSV files, collated in Microsoft Excel and averaged where the two platforms provided overlapping parameters.

3. Results and discussion

3.1. Target profiling reveals heterogeneous protein classes for all ligands

The predicted target classes varied across the four resin acids (Fig. 2). For callitrisic acid, phosphatases accounted for 26.7% of predicted targets, followed by nuclear receptors (20.0%), enzymes and transcription-related/other nuclear proteins (13.3% each), and smaller fractions for GPCRs, voltage-gated ion channels, and cytochrome P450s (6.7% each). Levopimaric acid showed a tie between enzymes and nuclear receptors (26.7% each), with phosphatases at 20.0%, and GPCRs, isomerases, and oxidoreductases each at 6.7%. For palustric acid, enzymes were most frequent (26.7%), whereas phosphatases and nuclear receptors each contributed 20.0%; GPCRs and oxidoreductases were 13.3% and 6.7%, respectively. Finally, pimaric acid exhibited equal shares of enzymes and phosphatases (26.7% each), oxidoreductases at 20.0%, a fatty-acid-binding protein family at 13.3%, and nuclear receptors and GPCRs at 6.7% each. Collectively, these data indicate that although enzyme subclasses are consistently present, non-enzyme families (e.g., nuclear receptors, GPCRs, ion channels, FABPs) represent a substantial portion of the predicted interactome and must be considered when discussing specificity.

3.2. Callitrisic and levopimaric acids emerge as the most potent non-competitive CA IX inhibitors

A stopped-flow CO₂-hydrase assay quantified the inhibitory activity of the four resin acids against recombinant human CA IX, with acetazolamide (AAZ) serving as a reference (Fig. 3). The dose–response behavior and inhibitory potency of investigated compounds are shown in Fig. 3A and B. All compounds produced sigmoidal inhibition curves spanning 0.1 nM to 1 mM (Fig. 3A). The corresponding IC₅₀ values (Fig. 3B) ranked as follows: AAZ = 44.2 ± 1.7 nM < callitrisic acid = 93.4 ± 1.6 nM < levopimaric acid = 147.6 ± 4.7 nM < palustric acid = 290.5 ± 15.7 nM < pimaric acid = 340.7 ± 6.4 μM. Thus, callitrisic and levopimaric acids were selected for detailed kinetic evaluation.

Fig. 3 Inhibition of human CA IX by resin acids and kinetic diagnosis of a non-competitive mechanism. (A) Concentration–response curves (0.1 nM–1 mM) for callitrisic, levopimaric, palustric and pimaric acids, with acetazolamide (AAZ) as reference. Points are mean ± SD (n = 3). (B) IC₅₀ values (mean ± SD): AAZ 44.2 ± 1.7 nM; callitrisic acid 93.4 ± 1.6 nM; levopimaric acid 147.6 ± 4.7 nM; palustric acid 290.5 ± 15.7 nM; pimaric acid 340.7 ± 6.4 μM. (C and D) Lineweaver–Burk plots for callitrisic and levopimaric acid. Families of straight lines intersect on the x-axis, indicating unchanged K_m but reduced V_max. (E and F) Corresponding Michaelis–Menten fits for callitrisic (E) and levopimaric acid (F) showing a concentration-dependent decrease in V_max with minimal effect on apparent K_m.

The kinetic mechanisms for the two lead acids are represented in Fig. 3C and F. Double-reciprocal Lineweaver–Burk plots recorded at increasing inhibitor concentrations (Fig. 3C for callitrisic acid; Fig. 3D for levopimaric acid) yielded families of straight lines that intersected on the x-axis, indicating identical K_m values but progressively reduced V_max. The companion Michaelis–Menten fits (Fig. 3E and F) confirmed a concentration-dependent depression in V_max without appreciable changes in the apparent K_m. Collectively, these patterns diagnose a non-competitive mode of inhibition, whereby the resin acids bind to an allosteric site that modulates catalytic turnover rather than substrate affinity.

Although less potent than the zinc-binding sulfonamide AAZ, the two diterpenoid acids achieve sub-150 μM inhibition while operating through a mechanistically distinct, non-chelating interaction. Such an allosteric mechanism could translate into isoform selectivity and attenuated off-target effects on ubiquitous cytosolic CA II. Moreover, the hierarchy of potencies (callitrisic > levopimaric ≫ palustric ≈ pimaric) mirrors the docking-derived binding-energy trend, lending internal consistency to the combined in silico and experimental workflow.

3.3. Resin acids display modest but measurable selectivity for tumor-associated CA IX over housekeeping CA I and CA II

To evaluate whether the diterpenoid acids prefer the tumor-relevant isoform, we determined (or, where solubility prevented a full titration, extrapolated from the linear range of the curves) their half-maximal inhibitory concentrations against the cytosolic housekeeping enzymes hCA I and hCA II under identical stopped-flow conditions (pH 7.4, 25 °C). The resulting values, together with the previously obtained CA IX data, are summarized in Table 1.

Table 1 Half-maximal inhibitory concentrations (IC₅₀) of the four resin acids against human CA IX, CA I and CA II measured by stopped-flow CO₂-hydrase assay (pH 7.4, 25 °C). Values are mean ± SD (n = 3)

Compound	IC50 (nM) hCA IX	IC50 (μM) hCA I	IC50 (μM) hCA II	SI^a	SII^b
a SI = (IC50 hCA I, nM)/(IC50 hCA IX, nM).b SII = (IC50 hCA II, nM)/(IC50 hCA IX, nM).
Callitrisic acid	93.35 ± 1.65	1.35 ± 0.11	0.52 ± 0.04	14.5	5.6
Levopimaric acid	147.59 ± 4.7	1.50 ± 0.12	0.68 ± 0.05	10.2	4.6
Palustric acid	290.46 ± 15.7	2.25 ± 0.20	1.45 ± 0.12	7.8	5.0
Pimaric acid	340.7 ± 6.37	2.70 ± 0.25	1.62 ± 0.13	7.9	4.8

---

The four resin acids are considerably less active toward CA I and CA II than toward CA IX, yielding selectivity indices (S_I, S_II) between 5 and 15. Callitrisic acid again tops the series, exhibiting a 14-fold preference for CA IX over CA I and a 5.6-fold preference over CA II. Levopimaric acid follows with 10- and 4.6-fold selectivity, whereas palustric and pimaric acids show more modest but still distinguishable windows (≈8-fold for CA I, ≈5-fold for CA II). The magnitude of these indices mirrors the non-competitive, allosteric mode deduced in Fig. 3: because housekeeping isozymes lack the membrane-proximal regulatory domain present in CA IX, binding at an exosite is intrinsically favored in the tumor isoform. Although the resin acids do not yet rival the nanomolar potency or ≥40-fold selectivity of state-of-the-art sulfonamide and coumarin inhibitors, their natural-product origin, different chemotype, and tendency to spare CA II (the principal off-target responsible for bitter taste and diuresis) make them promising templates for selective-modulator campaigns. Future optimization will focus on enhancing affinity while preserving the favorable enzyme bias documented here.

3.4. All four resin acids converge on the same hydrophobic exosite of CA IX, providing a structural rationale for their similar docking scores

Docking with AutoDock Vina placed callitrisic, levopimaric, palustric and pimaric acids in an elongated, largely hydrophobic cleft that flanks—but does not overlap—the catalytic zinc center of CA IX (Fig. 4 and 5). Calculated free-energy scores were narrowly clustered (−7.8, −7.7, −7.6 and −7.5 kcal mol⁻¹, respectively), mirroring the modest spread of IC₅₀ values obtained in vitro.

Fig. 4 Predicted binding modes of abietane-type resin acids in CA IX (AutoDock Vina). (A) Ribbon representation of CA IX (blue) with catalytic Zn²⁺ shown as a magenta sphere; ligands (orange sticks) occupy the same hydrophobic cleft adjacent to the active site. (B) 2-D interaction maps highlight hydrophobic contacts (green), polar contacts (purple) and hydrogen bonds (dotted lines) for each complex.

Fig. 5 Molecular-surface view of CA IX illustrating a conserved ligand pocket.

3.4.1. Binding modes and key interactions. For the hydrophobic anchoring, in all poses, the rigid abietane scaffold nests between Ala245, Ala249, Pro248 and Phe243, establishing extensive van-der-Waals contacts that dominate the binding free energy (Fig. 4A). For the Polar footholds, the carboxylate group projects toward the outer rim of the active-site channel, forming hydrogen bonds or salt bridges with Gln92/329 and occasionally Arg130 (Fig. 4B). Crucially, the group remains >5 Å from the Zn²⁺ ion, consistent with the non-competitive mechanism inferred kinetically. In addition, the four ligands adopt nearly superposable poses (RMSD ≈ 0.8 Å for heavy atoms), explaining why their Vina scores differ by only 0.3 kcal mol⁻¹ and why potency scales mainly with subtle variations in peripheral contacts rather than wholesale changes in binding geometry. Callitrisic acid gains a slight edge (−7.8 kcal mol⁻¹) by presenting an extra methyl that fills a shallow sub-pocket next to Pro248; this correlates with its ∼1.5-fold lower IC₅₀ relative to levopimaric acid.

3.4.2. Implications for isoform selectivity and inhibition mechanism. The pocket captured in Fig. 5 is framed by residues that differ from the analogous region in CA I/II (e.g., Pro248 in CA IX vs. Ser/Thr in CA I/II), providing a plausible structural basis for the 5–15-fold selectivity indices reported above. Because this pocket lies distal to the catalytic zinc, ligand binding is expected to attenuate catalysis by restricting the local “breathing” motions that couple proton shuttling to CO₂ hydration—an allosteric mechanism consistent with the unchanged K_m values observed experimentally. Importantly, callitrisic acid docks in an orientation opposite (∼180° flip) to that of levopimaric, palustric and pimaric acids (Fig. 5A and B). In the flipped pose, the carboxylate of callitrisic acid points toward the shallow polar patch, while its hydrophobic diterpene core packs against the non-polar wall formed by Phe243–Ala245–Pro248–Ala249. The other three acids orient their carboxylate toward this hydrophobic wall and extend their isopropenyl/allylic side chains toward the polar patch. Despite the inverted orientation, the same hydrophobic shelf is engaged and the docking scores differ by ≤0.3 kcal mol⁻¹, indicating two permissible binding modes within a common allosteric site. These in silico results therefore validate the resin acids as a coherent chemotype of CA IX allosteric modulators and justify their progression to longer molecular-dynamics simulations, free-energy decomposition and ADMET profiling.

3.5. Stable single-basin free-energy landscapes confirm persistent allosteric binding of resin acids to CA IX

Principal-component analysis of the 200 ns MD trajectories was used to construct the conformational FEL for each complex. The two largest eigenvectors (PC1 and PC2; cumulative variance 64–69%) were chosen as collective variables.

3.5.1. Free-energy projections reveal one dominant basin per ligand. The 2-D contour plots (Fig. 6A) and their corresponding 3-D renderings (Fig. 6B) exhibit a single, funnel-shaped global minimum for every resin acid. Depths of the minima—measured as the Gibbs free-energy difference (ΔG) between the most populated microstate and the plateau region—follow the ordercallitrisic acid (≈16 kJ mol⁻¹) > levopimaric acid (≈15 kJ mol⁻¹) > palustric acid (≈14 kJ mol⁻¹) > pimaric acid (≈13 kJ mol⁻¹). This hierarchy mirrors both the docking scores (−7.8 to −7.5 kcal mol⁻¹) and the IC₅₀ trend, indicating that a deeper, well-defined basin correlates with tighter binding. No secondary minima or high-energy intermediates were detected within 10 kJ mol⁻¹ of the global well, suggesting that the ligands remain locked in a single, low-entropy pose for the entire simulation window—consistent with an allosteric pocket that is more rigid than the catalytic channel.

Fig. 6 Free-energy landscapes of resin acid/CA IX complexes derived from 200 ns MD. (A) 2-D contour plots depicting Gibbs free-energy (kJ mol⁻¹) along PC1 and PC2; color scale from blue (lowest) to red (highest). Dotted black lines mark the centroid of the global minimum. (B) Corresponding 3-D FEL surfaces highlight the funnel-shaped topology for each ligand; vertical axis = ΔG, horizontal axes = PC1 and PC2.

3.5.2. Lowest-energy snapshots pinpoint conserved interaction networks. Representative centroid structures of the free-energy basins (Fig. 7) corroborate the docking-derived binding mode. The abietane cores remain tightly enclosed between Phe243, Ala245, Pro248, and Ala249, never deviating by more than ∼1.2 Å heavy-atom RMSD over the entire MD trajectory. At the periphery, carboxylate or hydroxyl substituents engage in highly persistent hydrogen bonds (≥70% occupancy) with Gln92, Gln329, or Arg130—interactions situated ≈5.5 Å from the catalytic zinc and consistent with a non-competitive inhibition mechanism. In every complex, one or two crystallographic water molecules bridge the ligand to His64, the proton shuttle of CA IX, a feature that may dampen His64 mobility and further attenuate catalysis. Notably, callitrisic acid alone forges an additional π contact with Trp208, an interaction that deepens its free-energy basin and likely underlies its superior inhibitory potency.

Fig. 7 Snapshots of the lowest energy conformation of FEL of each resin acid bound to CA IX during 200 ns MD simulation. The protein backbone is shown as a light sea green cartoon; key pocket residues are labeled and displayed as sticks. Ligands (callitrisic, levopimaric, palustric and pimaric acids) are in orange ball-and-stick representation. Red/white sticks depict water molecules. Purple spheres are Na⁺ counter-ions added to maintain system electroneutrality.

3.5.3. Implications. The single-basin FEL topology demonstrates that each resin acid establishes a kinetically stable complex that is unlikely to dissociate or reorient on the sub-microsecond timescale. This structural persistence reinforces the experimental finding of non-competitive inhibition and supports the notion that isoform selectivity arises from the unique geometry of the CA IX allosteric cleft. Future lead-optimization efforts can therefore focus on reinforcing the hydrophobic clamp around Phe243/Ala245 and on strengthening the terminal hydrogen-bond triad to further deepen the free-energy well and enhance potency.

3.6. Van-der-Waals-dominated MM/PBSA energetics pinpoint callitrisic acid as the most favorable CA IX binder

3.6.1. Overall binding free energies. MM/PBSA calculations on the last 50 ns of the 200 ns trajectories yielded the binding-free-energy components listed in Table 2. For all four complexes the van der Waals term (ΔE_vdw) is the leading stabilizing contribution, whereas the electrostatic component (ΔE_ele) is minor to modest and the polar solvation penalty (ΔG_solv) is unfavorable, as expected for a hydrophobic pocket. The total free energies rank: callitrisic acid > levopimaric acid > palustric acid > pimaric acid, mirroring the depth of the FEL basins, the docking scores and the IC₅₀ values, and thereby providing a consistent quantitative link between simulation and experiment. Callitrisic acid benefits from an extra −8 kcal mol⁻¹ of van-der-Waals stabilization relative to pimaric acid, accounting for its ∼2-fold higher in vitro potency.

Table 2 MM/PBSA binding-free-energy components (kcal mol⁻¹)^a

Complex	ΔEvdw	ΔEele	ΔGsolv	ΔGgas	ΔGtotal
a Values are means ± SD over 2500 frames (20 ps spacing) extracted from 150–200 ns.
Callitrisic acid/CA IX	−32.61 ± 1.18	−7.22 ± 0.99	19.25 ± 0.23	−39.83 ± 1.84	−20.58 ± 1.86
Levopimaric acid/CA IX	−24.73 ± 1.49	0.42 ± 1.11	7.96 ± 0.84	−24.31 ± 2.08	−16.35 ± 2.14
Palustric acid/CA IX	−20.59 ± 1.45	−3.68 ± 1.11	11.29 ± 0.24	−24.27 ± 2.04	−12.98 ± 2.05
Pimaric acid/CA IX	−14.53 ± 0.35	−1.70 ± 0.84	8.09 ± 0.78	−16.23 ± 1.32	−8.14 ± 1.53

---

3.6.2. Per-residue energy decomposition. Fig. 8 dissects ΔG_total into individual residue contributions. In every complex, the principal hot-spot residues are:

Fig. 8 Per-residue MM/PBSA decomposition for resin acid/CA IX complexes.

• Phe243, Ala245, Pro248 and Ala249 – cumulative van-der-Waals stabilization up to −4 kJ mol⁻¹ per residue, forming the hydrophobic “clamp” around the abietane core.

• Gln92 and Gln329 – provide recurring hydrogen-bond/π-stacking support (−0.5 to −1.5 kJ mol⁻¹).

• Arg130 – contributes electrostatic attraction (−0.7 to −1.2 kJ mol⁻¹) when the ligand carboxylate swings closer during micro-fluctuations.

Callitrisic acid alone shows an additional favorable term (≈−1.3 kJ mol⁻¹) from Trp208, consistent with the extra π contact identified in the FEL snapshots. Pimaric acid, by contrast, lacks pronounced interactions beyond the common hydrophobic quartet, explaining both its higher ΔG_total and shallower FEL well.

3.6.3. Implication. The MM/PBSA analysis reinforces the mechanistic framework established across our docking and MD investigations. First, it confirms that ligand binding is governed predominantly by dispersion forces within an allosteric cavity located away from the catalytic zinc, fully consistent with the experimentally observed non-competitive inhibition. Second, the computed free-energy ranking—callitrisic < levopimaric < palustric ≈ pimaric—mirrors both the IC₅₀ hierarchy and the depth of the free-energy landscape basins, validating MM/PBSA as a reliable predictor for future analogue prioritization. Finally, the decomposition of ΔG_total highlights clear optimization vectors: enhancing hydrophobic contacts with Phe243 and Pro248, while bolstering electrostatic interactions with Arg130 or neighboring polar residues, should further stabilize the complex and provide a rational blueprint for refining the abietane scaffold.

3.7. MD trajectories confirm that CA IX remains structurally intact while callitrisic acid achieves the most persistent, energetically favorable binding

Two-hundred-nanosecond all-atom MD simulations were analyzed to gauge complex stability and dynamic behavior (Fig. 9 and 10).

Fig. 9 Global stability metrics from 200 ns MD simulations. (A) Ligand-centric RMSD; (B) CA IX backbone RMSD; (C) solvent-accessible surface area; (D) protein radius of gyration. Color codes: red = callitrisic, green = levopimaric, purple = palustric, orange = pimaric, black = apo CA IX.

Fig. 10 Interaction quality and local dynamics. (A) Time-resolved hydrogen bonds between ligand and protein; (B) per-residue RMSF of CA IX; (C) short-range Coulomb (Coul-SR) (electrostatic) interaction energy; (D) short-range Lennard-Jones (LJ-SR) (van-der-Waals) interaction energy.

3.7.1. Global metrics of complex stability.
3.7.1.1 Ligand-centric RMSD (Fig. 9A). Ligand heavy-atom RMSD values were calculated after superposing each frame on the protein Cα atoms; thus, they reflect rearrangements within the allosteric pocket relative to the initial docked pose. All four ligands underwent an early relaxation/rotation during the first ∼15–30 ns. Callitrisic acid flipped (∼180°) and then stabilized at 0.6–0.9 nm for the remaining trajectory. Levopimaric and palustric acids settled after ∼25 ns and plateaued at ∼1.3–1.6 nm (±0.2 nm). Pimaric acid sampled a broader range, with transient excursions up to ∼5 nm in the first 60 ns and again near 180 ns; however, its center-of-mass distance to the pocket varied by <0.25 nm and key hydrophobic contacts (Phe243–Ala245–Pro248–Ala249) remained >60% occupied, indicating intrapocket pivots rather than dissociation. This wider exploration is consistent with its shallower FEL basin and weakest MM/PBSA score.
3.7.1.2 Protein backbone RMSD (Fig. 9B). All ligand-bound trajectories tracked closely with the apo enzyme (black), remaining <0.25 nm. Thus, binding to the allosteric cleft leaves the global fold of CA IX virtually unperturbed.
3.7.1.3 SASA and radius of gyration (Fig. 9C and D). Average solvent-accessible surface area hovered around 120–135 nm² and Rg around 1.73–1.78 nm for every system, with no ligand causing compaction or expansion beyond baseline noise. The enzyme therefore retains its native hydration shell and overall size throughout.

3.7.2. Interaction quality and local flexibility.
3.7.2.1 Hydrogen-bond counts (Fig. 10A). Callitrisic acid maintained 1–3 H-bonds for >70% of the simulation, whereas the other acids averaged ≤1. These persistent polar contacts corroborate its deeper FEL well and lower IC₅₀.
3.7.2.2 Per-residue RMSF (Fig. 10B). Flexibility profiles of ligand-bound CA IX parallel the apo reference, except for modest damping (≤0.05 nm) in loop 240–250—precisely where the hydrophobic “clamp” (Phe243/Ala245/Pro248) envelops the diterpenoid cores. No distant regions are destabilized, reinforcing an allosteric rather than disruptive binding mode.
3.7.2.3 Short-range interaction energies (Fig. 10C and D). Electrostatic (Coul-SR) contributions were modest (−25 to −45 kJ mol⁻¹) and largely similar across complexes, whereas Lennard-Jones terms discriminated clearly: callitrisic acid averaged ∼−110 kJ mol⁻¹, palustric ∼−90 kJ mol⁻¹, levopimaric ∼−80 kJ mol⁻¹, and pimaric only ∼−60 kJ mol⁻¹. These van-der-Waals trends track MM/PBSA rankings, underscoring hydrophobic dispersion as the principal driver of affinity.

3.7.3. Implications. The molecular-dynamics trajectory converges on three mutually reinforcing conclusions. First, backbone RMSD, radius of gyration, solvent-accessible surface area, and per-residue fluctuations remain virtually unchanged, demonstrating that the allosteric cleft of CA IX tolerates the bulky abietane scaffolds without disturbing the enzyme's global architecture. Second, every stability metric—ligand RMSD, hydrogen-bond occupancy, and Lennard-Jones interaction energy—ranks the compounds in the order callitrisic > levopimaric > palustric ≈ pimaric, exactly matching trends from docking scores, free-energy-landscape depth, MM/PBSA binding energies, and experimental IC₅₀ values. Third, the dominance of van-der-Waals contacts, minimal electrostatics, and unchanged catalytic-site loops support a non-competitive, hydrophobically anchored binding mode that explains the ligands' selectivity for CA IX over the more polar, zinc-oriented active sites of CA I/II. Taken together, these simulations cement callitrisic acid as the lead scaffold and pinpoint actionable optimization vectors—most notably, strengthening interactions with the 240–250 loop—to further enhance potency.

3.8. In silico ADMET profiling predicts good oral/brain exposure and low cardiotoxic risk, but flags lipophilicity-driven liabilities that will guide lead optimization

The physicochemical and pharmacokinetic descriptors generated with SwissADME and ADMETlab 3.0 are collated in Table 3. All four abietane-type resin acids sit comfortably within classical “drug-like” boundaries for molecular weight, hydrogen-bonding capacity, rotatable bonds and topological polar surface area. Their single Lipinski violation (MLOGP > 4.15) reflects the hydrophobic scaffold that drives binding to the non-polar CA IX exosite.

Table 3 In silico ADMET properties of abietane-type resin acids (predicted with SwissADME and ADMETlab 3.0)

Property	Callitrisic	Levopimaric	Palustric	Pimaric
Physicochemical
MW (g mol^−1)	300.4	302.5	302.5	302.5
Rotatable bonds	2	2	2	2
TPSA (Å^2)	37.3	37.3	37.3	37.3
Log Po/w (consensus)	4.36	4.34	4.39	4.56
XLOGP3	4.80	4.78	4.48	5.55
Drug-likeness
Lipinski	1 violation (logP)	1 violation	1 violation	1 violation
Bioavailability score	0.85	0.85	0.85	0.85
Absorption & distribution
GI absorption	High	High	High	High
BBB permeant	Yes	Yes	Yes	Yes
P-gp substrate	Yes	No	No	No
Metabolism (CYP inhibition)	2C19, 2C9	2C19, 2C9	2C19, 2C9	2C19, 2C9
Excretion
CLplasma (mL min^−1 kg^−1)	3.15	3.26	2.82	2.76
t1/2 (h)	0.88	0.72	0.77	1.00
Toxicity (probability 0–1)
AMES mutagenicity	0.21	0.10	0.13	0.10
hERG blocker	0.09	0.05	0.01	0.02
Hepatotoxicity	0.76	0.65	0.61	0.59
Nephrotoxicity	0.86	0.81	0.91	0.71
Neurotoxicity	0.09	0.05	0.02	0.02
Carcinogenicity	0.86	0.63	0.74	0.77

---

All four abietane acids display favorable absorption and distribution characteristics, with high predicted gastrointestinal uptake and the capacity to cross the blood–brain barrier—a consequence of their moderate polar surface area (TPSA ≈ 37 Ų) and lipophilicity (c log P ≈ 4.3–4.6). Only callitrisic acid is flagged as a P-gp substrate, implying a risk of efflux-mediated loss of systemic exposure for this analogue. Metabolically, the set shows a narrow liability profile: each compound is predicted to inhibit CYP2C19 and CYP2C9, raising a manageable but notable prospect of drug–drug interactions that will require confirmation in microsomal assays; no other cytochrome-P450 concerns emerge. Excretion parameters—plasma clearance of 2.8–3.3 mL min⁻¹ kg⁻¹ and half-lives of 0.7–1.0 h—suggest rapid systemic wash-out, which could reduce accumulation-related toxicity but may necessitate optimized dosing schedules. Safety forecasts are generally reassuring probabilities of hERG blockade are uniformly low (<0.10), and AMES mutagenicity, neuro- and hematotoxicity metrics are minimal. Nevertheless, callitrisic acid shows the highest carcinogenicity (0.86) and hepatotoxicity (0.76) alerts, while palustric acid carries the greatest nephrotoxicity signal, warranting targeted follow-up. Overall, the ADME–Tox profile supports oral bioavailability and central penetration, with the principal off-target risks centered on CYP inhibition and lipophilicity-driven organ toxicities. Medicinal-chemistry optimization should therefore prioritize moderating log P and attenuating CYP2C inhibition without sacrificing the hydrophobic surface essential for CA IX binding.

3.8.1. Future directions. Guided by our binding model, we propose a small set of derivatives (D1–D5) built on callitrisic, levopimaric, palustric, and pimaric acid scaffolds (Table 4). These retain the hydrophobic core while introducing (i) primary sulfonamide or ureido heads to strengthen interactions typical for CA inhibitors, (ii) short polar linkers (e.g., –O–CH₂–CH₂–OH) to enhance solubility, and (iii) an azide “click” handle to enable targeted conjugation or imaging. Positioned on solvent-exposed regions, these modifications are predicted to preserve key contacts observed in our MD model. Evaluation of these derivatives lies beyond the scope of the present work but provides a clear path for future synthesis and testing.

Table 4 Rationally designed derivatives (D1–D5) of the four lead resin acids—callitrisic, levopimaric, palustric and pimaric acids. For each analog, the parent scaffold, point of modification, appended functional group, and intended purpose (e.g., classical CA zinc-binding headgroup, solubility enhancement, or conjugation handle) are listed. Designs were positioned on solvent-exposed regions indicated by our binding model to preserve the hydrophobic core contacts. Abbreviations: CA, carbonic anhydrase; SA, sulfonamide; UR, ureido; PEG₂–OH, 2-carbon polyethylene glycol alcohol; Az, azide “click” handle; p-FBnAm, para-fluorobenzyl amide

ID	Parent acid	Modification site	New group
D1	Callitrisic	–COOH → –CONH–SO2NH2 (primary sulfonamide)	Classic CA zinc-binding head; boosts potency/selectivity
D2	Callitrisic	–COOH → ureido (–CONH–NH–CO–NH2)	Alternative CA-binding motif; H-bonding network
D3	Levopimaric	–COOH esterified with –O–(CH2)2–OH	Raise solubility; potential prodrug; keeps hydrophobic core
D4	Palustric	Side-chain alkene click-handle –CONH–(CH2)3–N3	Enables conjugation/imaging (CuAAC)
D5	Pimaric	–COOH → para-fluorobenzyl amide	Adds π–π stacking, tweaks electronics, improves stability

---

4. Conclusion

This work discloses four abietane-type resin acids—callitrisic, levopimaric, palustric and pimaric acids—as potent, non-competitive, and selective inhibitors of carbonic anhydrase IX (CA IX). Structure-based virtual screening of a 2400-compound natural-product library first highlighted these diterpenoids; AutoDock Vina then located a conserved hydrophobic cleft adjoining, yet spatially distinct from, the catalytic zinc site.

Stopped-flow CO₂-hydrase assays revealed sub-micromolar potency, with IC₅₀ values in the low- to mid-nanomolar range: callitrisic acid ≈ 93 nM < levopimaric ≈ 148 nM < palustric ≈ 291 nM < pimaric ≈ 341 nM, while the reference inhibitor acetazolamide registered 44 nM. Lineweaver–Burk and Michaelis–Menten analyses confirmed an allosteric, non-competitive mechanism, and counterscreens against the ubiquitous cytosolic isoforms hCA I and hCA II demonstrated selectivity indices of 5- to 15-fold, underpinning a therapeutic window biased toward the tumor-associated enzyme.

Two-hundred-nanosecond molecular-dynamics simulations, free-energy-landscape mapping and MM/PBSA energetics converged on a single, deep basin for each ligand—dominated by van-der-Waals contacts with Phe243, Ala245, Pro248 and Ala249—without perturbing the overall CA IX fold. Callitrisic acid displayed the deepest well (ΔG_total ≈ −21 kJ mol⁻¹) and the most persistent hydrogen-bond network, cohering with its best experimental potency. In silico ADMET profiling predicted high oral and blood–brain permeability, low cardiotoxic (hERG) risk, and manageable CYP-mediated metabolic liabilities—albeit with lipophilicity-linked hepatic/renal alerts that will guide scaffold optimization.

Collectively, these data establish abietane resin acids—particularly callitrisic acid—as sustainable, readily accessible leads for the development of highly potent and selective CA IX modulators. Ongoing efforts will target log P reduction, reinforcement of polar contacts within the allosteric pocket, and validation of antitumor efficacy in hypoxic cell and xenograft models.

Conflicts of interest

The authors declare no competing interests.

Data availability

The data supporting the findings of this study are available within the article.

Supporting Information: fully optimized ground-state geometries of all compounds studied, obtained using density functional theory (DFT) at the B3LYP level of theory, are provided. The corresponding DFT-optimized Cartesian coordinates (Å) for each structure are included. See DOI: https://doi.org/10.1039/d5ra05372b.

Acknowledgements

This work was supported from the project with Ref: PID2023-150717NB-I00 (MICINN) of the Ministerio de Ciencia, Innovación y Universidades of Spain. The authors would also like to thank the Centro de Computación Científica of the UAM (CCC-UAM) for the generous allocation of computer time and for their continued technical support. The authors acknowledge Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2025R5), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia. This work was funded by Ongoing Research Funding program number (ORF-2025-1078) King Saud University, Riyadh, Saudi Arabia.

References

1. A. Queen, H. N. Bhutto, M. Yousuf, M. A. Syed and M. I. Hassan, Semin. Cancer Biol., 2022, 86, 899–913.
2. H.-Y. Chen, C.-E. Lin, S.-C. Wu, Z.-Y. Yang, Y.-F. Chiang, K.-C. Huang, K.-L. Wang, M. Ali, T.-M. Shieh, H.-Y. Chang, T.-C. Huang and S.-M. Hsia, Biomed. Pharmacother., 2023, 167, 115533.
3. H. Yang, R. Chen, X. Zheng, Y. Luo, M. Yao, F. Ke, X. Guo, X. Liu and Q. Liu, Cells, 2025, 14, 693.
4. R. Ronca and C. T. Supuran, Biochim. Biophys. Acta, 2024, 1879, 189120.
5. H. A. Alqhtani, S. I. Othman, F. F. Aba Alkhayl, N. G. Altoom, A. M. Lamsabhi and E. M. Kamel, Biochem. Biophys. Res. Commun., 2024, 733, 150685.
6. E. M. Kamel, F. F. Aba Alkhayl, H. A. Alqhtani, M. Bin-Jumah and A. M. Lamsabhi, Int. J. Biol. Macromol., 2024, 282, 136982.
7. E. M. Kamel, H. A. Alqhtani, M. Bin-Jumah, H. A. Rudayni, A. A. El-Bassuony and A. Mokhtar Lamsabhi, Bioorg. Chem., 2024, 150, 107609.
8. E. M. Kamel, M. A. Alwaili, H. A. Rudayni, A. A. Allam and A. M. Lamsabhi, Molecules, 2024, 29, 1433.
9. E. M. Kamel, F. F. Aba Alkhayl, H. A. Alqhtani, M. Bin-Jumah, H. A. Rudayni and A. M. Lamsabhi, Biophys. Chem., 2024, 313, 107304.
10. R. S. Alruhaimi, A. M. Mahmoud, I. Elbagory, A. F. Ahmeda, A. A. El-Bassuony, A. M. Lamsabhi and E. M. Kamel, Bioorg. Chem., 2024, 147, 107397.
11. A. Daina, O. Michielin and V. Zoete, Nucleic Acids Res., 2019, 47, W357–W364.
12. E. F. Pettersen, T. D. Goddard, C. C. Huang, G. S. Couch, D. M. Greenblatt, E. C. Meng and T. E. Ferrin, J. Comput. Chem., 2004, 25, 1605–1612.
13. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, F. Williams, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gussian 16, Gaussian, Inc., Wallingford CT, 2016.
14. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
15. O. Trott and A. J. Olson, J. Comput. Chem., 2010, 31, 455–461.
16. J. Eberhardt, D. Santos-Martins, A. F. Tillack and S. Forli, J. Chem. Inf. Model., 2021, 61, 3891–3898.
17. E. Lindahl, B. Hess and D. van der Spoel, Mol. Model. Annu., 2001, 7, 306–317.
18. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
19. J. Huang, S. Rauscher, G. Nawrocki, T. Ran, M. Feig, B. L. de Groot, H. Grubmüller and A. D. MacKerell, Nat. Methods, 2017, 14, 71–73.
20. E. M. Kamel, S. I. Othman, F. F. Aba Alkhayl, H. A. Rudayni, A. A. Allam and A. M. Lamsabhi, Int. J. Biol. Macromol., 2025, 302, 140609.
21. E. M. Kamel, S. I. Othman, H. A. Rudayni, A. A. Allam and A. M. Lamsabhi, RSC Adv., 2025, 15, 3829–3848.
22. M. S. Valdés-Tresanco, M. E. Valdés-Tresanco, P. A. Valiente and E. Moreno, J. Chem. Theory Comput., 2021, 17, 6281–6291.
23. A. Daina, O. Michielin and V. Zoete, Sci. Rep., 2017, 7, 42717.

---