Shikonin and its derivatives as selective modulators of the vitamin D receptor: A molecular framework for precision therapeutics in diabetic nephropathy

2.1 Ligand scaffold selection and derivative design

The initial chemical scaffold was derived from shikonin (PubChem CID: 479503), a naturally occurring naphthoquinone known for its extensive biological activity (Kim et al., 2025; Kim and Bolton, 2024). To expand the chemical diversity and explore structural analogs with potentially enhanced bioactivity and drug-likeness, we utilized the ADMETopt web server (Yang et al., 2018). This robust platform integrates scaffold hopping and ADMET property prediction for optimized ligand design. The shikonin scaffold was input into ADMETopt, and 100 structurally diverse derivatives were generated through algorithmic modifications to optimize pharmacokinetic and pharmacodynamic properties. All derivatives were evaluated for drug-likeness, aqueous solubility, oral bioavailability, and blood-brain barrier permeability using the server’s integrated predictive models. This ensured that only ADMET-compliant candidates progressed to subsequent stages.

2.2 Protein preparation and receptor setup

The VDR crystal structure was retrieved from the Protein Data Bank (PDB ID: 1DB1) (Berman, 2000; Burley et al., 2025; Rochel et al., 2000). Protein preparation was carried out using the Protein Preparation Wizard in the Schrödinger suite. Water molecules beyond 5 Å from the binding site were removed. Protonation states were adjusted to pH 7.4 using PROPKA, missing side chains and loops were modeled, and the structure was minimized using the OPLS3e force field (Jorgensen et al., 1996; Jorgensen and Tirado-Rives, 1988; Rostkowski et al., 2011). The receptor grid was generated around the co-crystallized ligand binding site, ensuring accurate placement of the active site centroid for docking.

2.3 Molecular docking and pose filtering

All 100 shikonin derivatives were docked against the prepared VDR using Glide XP (Extra Precision) mode (Friesner et al., 2006, 2004). Glide XP improves prediction accuracy by recognizing when molecules are tightly packed in water-repelling (hydrophobic) areas, and by giving extra value to strong hydrogen bonds and π–π stacking interactions are key features that help identify highly effective binders. Each ligand was flexibly docked to explore multiple conformations within the VDR active site. After docking, pose filtering was conducted using Prime MMGBSA, which calculates the free binding energy (ΔG_bind) based on the docked complex (“Schrödinger Release 2025-1: Prime, Schrödinger, LLC, New York, NY, 2025.,” n.d.). Ligand-receptor poses with MMGBSA scores better than −50 kcal/mol were selected for further dynamic analysis. This dual-tier filtering ensured the retention of the most energetically favorable and structurally plausible binding conformers.

2.4 Molecular dynamics simulation

To explore the dynamic stability and binding behavior of the top ligand-VDR complexes, we performed 500 ns molecular dynamics (MD) simulations using Desmond, also a part of the Schrödinger suite (Bowers et al., 2006; Shaw, 2020). Complexes were solvated in a TIP3P explicit water model within an orthorhombic box, with a 10 Å buffer between the protein and the box edges. Counterions (Na⁺ or Cl⁻) were first added to neutralize the system, and subsequently, 0.15 M NaCl was introduced to mimic physiological conditions. The system was relaxed using Desmond’s default relaxation protocol, followed by production runs using the NPT ensemble at 300 K and 1 atm. Trajectories were recorded at 100 ps intervals, and Root Mean Square Deviation (RMSD), Root Mean Square Fluctuation (RMSF), and hydrogen bond analyses were conducted using the Simulation Interaction Diagram tool.

2.5 Principal component analysis

Principal component analysis (PCA) was performed using the Bio3D package in R to assess conformational changes and collective motions during MD simulations (Barry J. Grant et al., 2006; Barry J Grant et al., 2006; Grant et al., 2021). The protein backbone coordinates were extracted from the trajectory files, aligned, and used to generate the covariance matrix of atomic fluctuations. Eigenvectors and eigenvalues were computed to determine the principal motions, with the first two components (PC1 and PC2) plotted to visualize primary conformational space sampling. This analysis identified structural transitions and energetically preferred states induced by ligand binding.

2.6 Free energy landscape mapping

The FEL was reconstructed based on the PC1 and PC2 data obtained from PCA using the GeoMeasures plugin in PyMOL (DeLano, 2002; Kagami et al., 2020). The conformational space was projected in a 2D energy surface map, where minima basins indicated thermodynamically stable states. This visualization provided insight into the folding pathways and structural stabilities of the VDR-ligand complexes during the simulation timeframe.

2.7 Structural superimposition of minima states

To assess structural convergence and deviation between initial docking poses and dynamically favored conformers, superimposition of the global minima conformations with their respective initial docking poses was carried out using UCSF Chimera (Pettersen et al., 2004). RMSD values between the docked and FEL-minima conformations were calculated to quantify structural deviations and binding pocket flexibility. This comparative approach allowed the evaluation of binding site plasticity and ligand-induced conformational adjustments in the receptor.

3.1 Scaffold identification and core structure

The molecular template used for derivative prediction was based on the structure of shikonin, a naturally occurring naphthoquinone with multiple hydroxyl and carbonyl functionalities. The initial scaffold (Fig. 1(a)) exhibited a conjugated system favorable for redox activity and biological interactions. A simplified core framework (Fig. 1(b)), retaining the essential naphthoquinone motif, was extracted and used for further structural diversification via scaffold hopping. This structure was specifically chosen due to its favorable geometry for hydrogen bonding, moderate lipophilicity, and previously documented biological relevance.

Close

Fig. 1.

The predicted core Shikonin scaffold was selected for a derivative design using ADMETopt. (a) The predicted scaffold in the Shikonin structure was drawn using the Chemical Sketch tool Chemaxon (blue circle). (b) The scaffold was separated and selected for a new derivative design.

Export to PPT

3.2 Derivative design and ADMET optimization

To explore potential drug-like derivatives of the shikonin scaffold, we utilized the ADMETopt web server. This platform enabled the automated generation and evaluation of 100 unique derivatives using scaffold hopping combined with predictive ADMET filtering (Supplementary_A and Supplementary_B). The core structure was modified by introducing various substituents, including hydroxyl, halogen (Cl, F, Br), amino, ether, and small ring moieties. These substitutions aimed to improve pharmacokinetic behavior while retaining or enhancing biological activity.

The parent scaffold possessed a molecular weight of 288.30, AlogP of 2.12, three hydrogen bond donors, and five acceptors, with a predicted aqueous solubility (logS) of −3.83. The generated derivatives exhibited a range of molecular weights, primarily between 280–520 Da, while maintaining drug-like physicochemical properties. AlogP values for most derivatives were within 3.5–5.5, indicating favorable membrane permeability. The number of hydrogen bond donors and acceptors was also tuned to fall within optimal bounds defined by Lipinski’s Rule of Five.

Drug-likeness scores ranged from 0.66 to 0.87, while synthetic accessibility values spanned from 4.2 to 5.7. These values suggest that most designed molecules maintained a balance between structural complexity and ease of synthesis. Many derivatives included bioisosteric replacements and small polar fragments that preserved the pharmacophoric core of shikonin while improving systemic bioavailability and metabolic stability.

3.3 Highlights of generated derivatives

Several key trends emerged from the analysis of the ADMET-optimized derivatives. Fluorinated derivatives generally exhibited slightly higher AlogP values and enhanced drug-likeness scores, likely due to improved membrane permeability and metabolic resistance. Halogenated analogs incorporating chlorine or bromine at ortho and meta positions also demonstrated favorable synthetic profiles and ADMET properties.

Nitrogen-containing substituents, particularly those involving pyrrole, amide, or amine groups, contributed to enhanced hydrogen bonding capabilities and better fitting geometries for potential protein interactions. In some cases, hydroxylated analogs showed increased solubility and retained good binding potential based on predicted drug-likeness. Derivatives with balanced polar and non-polar groups achieved optimal interaction profiles without exceeding the acceptable ranges for hydrogen bonding features.

From the 100 designed compounds, a significant portion (over 60%) achieved a drug-likeness score above 0.80 and a synthetic accessibility score below 5.0. These candidates had high potential for downstream in silico screening against the VDR.

3.4 Selection for virtual screening

After ADMET evaluation, all 100 derivatives were retained for virtual screening using molecular docking and dynamics. These molecules were selected based on the criteria of drug-likeness scores above 0.75 and synthetic accessibility not exceeding 5.5. The selection prioritized compounds with acceptable lipophilicity and well-distributed polar groups, which are ideal for binding to a nuclear receptor such as VDR.

The derivatives displayed chemical diversity with strategically positioned functional groups that can engage in hydrogen bonding, π-stacking, and hydrophobic interactions. These properties are essential for effective interaction within the VDR ligand-binding pocket. The comprehensive dataset of the designed molecules, including molecular properties, ADMET predictions, and compound identifiers, is summarized in Supplementary_A submitted alongside the manuscript.

This scaffold-based approach, enhanced through ADMET-driven filtering, enabled the rational identification of structurally diverse, pharmacokinetically favorable derivatives of shikonin that are promising candidates for targeting DN through modulation of the VDR.

3.5 Virtual screening and binding free energy analysis

Based on the virtual screening and MMGBSA analysis (Supplementary_A Table S2), four ligands: Shikonin, Derivative_02, Derivative_12, and Derivative_39—were identified as top binders to the VDR (Fig. 2). Shikonin exhibited a Glide docking score of −8.912 kcal/mol and a corresponding MMGBSA binding free energy (ΔG_bind) of −45.26 kcal/mol, confirming its moderate binding affinity and stability. Derivative_02 demonstrated a superior docking score of −10.14 kcal/mol and a ΔG_bind of −65.42 kcal/mol, indicating stronger and more stable interactions. Similarly, Derivative_12 and Derivative_39 showed docking scores of −9.87 kcal/mol and −9.65 kcal/mol, with MMGBSA binding energies of −62.37 kcal/mol and −60.85 kcal/mol, respectively. The close correlation between docking energies and MMGBSA results suggests that the top-ranked compounds fit well within the VDR binding pocket and form energetically favorable interactions. These findings collectively support the enhanced binding potential of the selected derivatives over the parent compound, establishing them as promising candidates for further molecular dynamics simulations and structural validation.

Close

Fig. 2.

2D structure of the selected ligands (a) Shikonin, (b) Derivative_02, (c) Derivative_12, and (d) Derivative_39.

Export to PPT

3.6 Molecular interaction analysis

Molecular docking revealed distinct interaction profiles for the parent compound shikonin and the three selected derivatives (Derivative_02, Derivative_12, and Derivative_39) within the VDR binding pocket. Shikonin formed a key hydrogen bond with Ser278, supported by a broad hydrophobic interaction network involving residues such as TYR143, Tyr147, LEU230, LEU233, VAL234, ILE268, ILE271, MET272, TRP286, CYS288, TYR295, VAL300, and LEU313, which stabilized its orientation within the active site.

Derivative_12 displayed enhanced polar interactions, forming hydrogen bonds with TYR143, SER237, ARG274, and SER278. It also engaged in extensive hydrophobic contacts with similar residues seen in shikonin, including PHE150, PHE153, TYR236, and TRP286, which suggests improved binding stability and pocket occupancy.

Derivative_02 showed a prominent hydrogen bond with Ser275 and significant hydrophobic interactions with residues TYR143, PHE150, LEU230, ILE271, and TRP286. It also exhibited π–π stacking interactions with TRP286, which occurred twice, indicating a strong aromatic stabilization component contributing to binding affinity.

Derivative_39 formed a hydrogen bond with SER237 and shared many hydrophobic contact points with other compounds, including TYR143, TYR147, PHE150, LEU233, TRP286, and TYR295. Similar to Derivative_02, it also displayed double π–π stacking interactions with TRP286, reinforcing its structural stabilization within the receptor pocket.

Collectively, the enhanced polar and aromatic interactions observed in the derivatives, especially with key residues like TRP286, SER278, and TYR143, support their superior docking scores and MMGBSA binding energies, suggesting improved receptor binding stability over the parent molecule shikonin (Fig. 3).

Close

Fig. 3.

3D and 2D interaction diagrams of Shikonin and its derivatives with the Vitamin D Receptor (VDR). (a, b) Shikonin–VDR complex showing the 3D binding pose within the ligand-binding domain and corresponding 2D interaction diagram. (c, d) Derivative_02–VDR complex illustrating ligand orientation in the binding pocket and key non-covalent interactions. (e, f) Derivative_12–VDR complex highlighting hydrogen bonds and hydrophobic contacts stabilizing the interaction. (g, h) Derivative_39–VDR complex with detailed binding pose and interaction map. Ligands are shown in stick representation (pink), receptor residues in green cartoon, and interactions are color-coded as indicated in the legend (hydrogen bonds, hydrophobic contacts, polar interactions, π–π stacking, salt bridges, etc.).

Export to PPT

3.7 Molecular dynamics simulation analysis

To probe the temporal stability of the VDR in complex with Shikonin and three structurally distinct derivatives (Derivative_02, Derivative_12, and Derivative_39), 500 ns molecular dynamics (MD) simulations were carried out. Structural behavior was evaluated using RMSD analysis of both protein and ligand trajectories.

3.8 RMSD analysis

3.8.1 Shikonin–VDR complex

In Fig. 4(a), the protein RMSD (green trace) starts around 1.5 Å and quickly equilibrates to approximately 2.1 Å, remaining consistent throughout the 500 ns simulation with minimal deviation. The ligand RMSD (red trace) exhibits initial fluctuations up to 2.8 Å within the first 100 ns, reflecting early-stage adjustments in the binding site. Notably, after this equilibration phase, the ligand RMSD consistently remains below 2.0 Å for the rest of the trajectory, indicating that Shikonin adopts a well-defined and stable conformation within the VDR pocket after initial accommodation.

Close

Fig. 4.

RMSD analysis of VDR–ligand complexes. (a) Shikonin–VDR, (b) Derivative_02–VDR, (c) Derivative_12–VDR, and (d) Derivative_39–VDR complexes monitored over 500 ns molecular dynamics simulations.

Export to PPT

3.9 RMSF analysis

Protein RMSF analysis (Fig. 5) revealed localized flexibility across all systems. Terminal regions and flexible loops displayed expected mobility, reaching 4.5–5.0 Å in some cases. The binding domain, particularly residues between 230–313, maintained low RMSF values below 2.0 Å, especially in complexes with shikonin and Derivative_02, which supports the RMSD findings of stable core interactions. Derivative_12 and Derivative_39 complexes exhibited modest increases in flexibility within the binding region, corresponding with their higher ligand fluctuations.

Close

Fig. 5.

Protein RMSF profiles of VDR residues in complex with Shikonin and its derivatives. (a) Shikonin–VDR, (b) Derivative_02–VDR, (c) Derivative_12–VDR, and (d) Derivative_39–VDR complexes.

Export to PPT

Ligand RMSF plots (Fig. 6) provided atomic-level resolution of ligand motion. Shikonin (Fig. 6a) exhibited a narrow fluctuation band, with most atoms oscillating between 1.0–1.8 Å, and reduced noise after 100 ns, reinforcing its stable anchoring post-initial relaxation. Derivative_02 (Fig. 6b) showed similar restraint, with atomic fluctuations between 1.0–2.0 Å, confirming its compact binding pose. Derivative_12 (Fig. 6c) displayed a broader RMSF range with terminal groups reaching up to 2.9 Å, reflecting internal rearrangement. Derivative_39 (Fig. 6d) maintained uniform fluctuation between 1.4–2.5 Å, indicating moderate but sustained internal mobility.

Close

Fig. 6.

Ligand RMSF profiles of VDR residues in complex with Shikonin and its derivatives. (a) Shikonin–VDR, (b) Derivative_02–VDR, (c) Derivative_12–VDR, and (d) Derivative_39–VDR complexes.

Export to PPT

3.10 Dynamic interaction consistency with docking predictions

Analyzing molecular interactions sustained for more than 30% of the 500 ns MD simulation trajectory provides insights into the durability and relevance of ligand-receptor contacts initially predicted by docking (Fig. 7). We can discern which contacts are structurally and energetically significant for stable ligand binding within the VDR pocket by comparing these long-term interactions with those identified in the static docking poses.

Close

Fig. 7.

Post-simulation interaction frequency plots of VDR–ligand complexes over 500 ns MD simulations. The plots depict residues involved in ligand binding, with only interactions persisting for >30% of the simulation time shown. (a) Shikonin, (b) Derivative_02, (c) Derivative_12, and (d) Derivative_39 complexes. Persistent hydrogen bonds, hydrophobic contacts, ionic interactions, and water bridges are represented according to the legend, highlighting residues contributing most significantly to complex stability.

Export to PPT

Post-MD data reveal that SER278 continues to form hydrogen bonds with shikonin, consistent with initial docking predictions, where it was the primary polar contact. LEU233 remains a key hydrophobic anchor, although the rich initial network involving residues like TYR143, TYR147, and ILE271 was not retained beyond 30% of the simulation time. The reduction in interaction diversity suggests that the ligand settled into a more minimalistic yet stable configuration, maintaining core interactions while discarding transient ones.

The MD trajectory shows that SER237, SER275, and SER278 form stable hydrogen bonds throughout the simulation, indicating a more polar-engaged binding mode than initially predicted. Among the initial contacts, only SER278 and LEU233 persist. ILE271, which was not a dominant docking contact, emerges as a consistent hydrophobic interaction partner, highlighting the value of dynamic simulation in uncovering stabilizing residues not evident in docking. These shifts suggest an adaptive binding conformation where the ligand realigns to engage polar residues deeply embedded in the binding pocket.

The interaction profile of Derivative_02 undergoes notable evolution. The MD data identify SER275 and GLN317 as key hydrogen bond donors, while TRP286 consistently participates in hydrophobic and π–π stacking interactions. This is particularly important as TRP286 was also predicted to form π–π stacking interactions in the docking pose, indicating a robust aromatic engagement that remains stable over time. The emergence of GLN317, not previously identified in the docking results, suggests ligand migration or reorientation that favors this hydrogen bond partner.

This derivative also maintains hydrogen bonds with SER237 and SER275, mirroring Derivative_12 in polar engagement strategy. TRP286 remains the central hydrophobic and π-interaction residue, similar to its role in the docking phase. This enduring interaction reflects its spatial and energetic centrality in the ligand’s binding conformation. Despite Derivative_39’s higher RMSD, the constancy of these interactions supports the idea of a dynamically flexible but site-anchored ligand.

3.11 MM/GBSA binding free energy analysis

Binding free energy decomposition was performed using the MM/GBSA method to dissect the energetic contributions governing ligand interaction with the VDR. Key components analyzed included van der Waals (vdW), Coulombic (electrostatic), lipophilic, hydrogen bonding, solvation (polar desolvation), covalent packing, and ligand strain energy.

Shikonin showed a robust binding profile with a total ΔG_bind of −54.32 ± 4.29 kcal/mol. It was primarily stabilized by strong van der Waals interactions (−43.65 ± 1.94 kcal/mol) and lipophilic contacts (−19.85 ± 1.35 kcal/mol). The electrostatic (Coulombic) contribution was comparatively lower (−7.69 ± 2.73 kcal/mol), and solvation effects (ΔG_solv = 16.48 ± 2.34 kcal/mol) partially offset the binding energy. Minimal strain energy (2.09 ± 0.91 kcal/mol) and moderate hydrogen bonding (−0.62 ± 0.32 kcal/mol) contributed to a stable and well-packed binding pose.

Derivative_12 yielded the most favorable binding energy across all ligands, with a ΔG_bind of −59.48 ± 6.13 kcal/mol. This was driven by significant Coulombic interactions (−24.98 ± 5.74 kcal/mol) and balanced lipophilic energy (−19.99 ± 1.75 kcal/mol). Despite a slightly elevated strain energy (3.45 ± 1.74 kcal/mol), the overall energy profile suggests a tightly bound and energetically favorable interaction mode, supported by stable hydrogen bonds and moderate solvation energy (23.06 ± 2.38 kcal/mol).

Derivative_02 demonstrated exceptional van der Waals (−45.17 ± 2.05 kcal/mol) and Coulombic contributions (−27.36 ± 7.21 kcal/mol), resulting in a strong direct binding interface. However, its total ΔG_bind was slightly reduced to −53.38 ± 4.21 kcal/mol, primarily due to a significant solvation penalty (43.54 ± 7.71 kcal/mol), which diluted the net energetic gains from receptor-ligand interactions. Its low ligand strain energy (1.66 ± 0.71 kcal/mol) and hydrogen bonding (−0.81 ± 0.14 kcal/mol) reflect favorable conformational accommodation.

Derivative_39 presented a balanced energetic profile with a ΔG_bind of −53.82 ± 4.79 kcal/mol. This was supported by van der Waals (−39.97 ± 2.13 kcal/mol) and lipophilic (−22.76 ± 1.72 kcal/mol) interactions, while its Coulombic term (−20.60 ± 7.07 kcal/mol) contributed moderately to stability. The solvation energy (30.19 ± 7.07 kcal/mol) and low strain (1.71 ± 1.32 kcal/mol) reflected a compact and adaptable binding mode, while hydrogen bonding (−0.72 ± 0.23 kcal/mol) added polar stability (Table 1).

3.12 Principal component analysis

PCA was performed to extract and compare the significant collective motions of the Vitamin D Receptor (VDR) in complex with Shikonin and its three top-performing derivatives across a 500 ns MD simulation. The scatter plots show projections of principal components (PC1 vs. PC2, PC2 vs. PC3, and PC1 vs. PC3), and the scree plots indicate how much variance each eigenvector captures. Red to blue color gradients represent the time evolution of motion.

The PCA projection of Shikonin (Fig. 8a) displays clearly defined, non-overlapping conformational clusters, particularly in the PC1–PC2 and PC1–PC3 planes. These clusters reflect transitions between distinct conformational subspaces, likely related to ligand-induced rigidification. PC1 explains 28.07% of the total motion, while the top five eigenvectors collectively cover approximately 66.3%, indicating that a few dominant modes capture a large fraction of the conformational dynamics. This behavior reflects highly ordered dynamics and supports Shikonin’s ability to stabilize the VDR conformation.

Close

Fig. 8.

PCA of VDR–ligand complexes; (a) Shikonin–VDR, (b) Derivative_02–VDR, (c) Derivative_12–VDR, and (d) Derivative_39–VDR complexes. The Projection of the first two principal components shows dominant motions sampled during simulations. Compact clusters indicate stable conformational states, whereas broader distributions reflect higher flexibility.

Export to PPT

The Derivative_02 complex (Fig. 8b) shows a well-resolved transition between two dense conformational basins along PC1 and PC2. The smooth transition suggests the receptor samples between two quasi-stable conformational states. PC1 contributes 21.27% of the variance, and the first five PCs together explain 56.8% of the total motion. This moderate dimensional reduction indicates stable yet flexible binding, which matches its consistent RMSD and persistent hydrogen bonding observed during simulation.

For Derivative_12 (Fig. 8c), PCA shows more diffuse and spatially expanded clusters, with greater dispersion along all component projections. The separation in PC1–PC2 and PC2–PC3 planes is less distinct, indicating broader conformational sampling. PC1 accounts for only 13.43% of the motion, while the top five components capture 56.9%, suggesting a higher degree of conformational entropy. This dispersion corroborates earlier RMSD and RMSF results, which showed higher fluctuation and transient instability in this complex.

The PCA scatter plots for Derivative_39 (Fig. 8d) reveal interconnected conformational clusters and evidence of continuous movement across PCs. Despite broad structural sampling, PC1 (21.87%) and PC2 (14.91%) suggest some directional preference in motion. The cumulative variance captured by the top five PCs reaches 68.1%, the highest among all complexes. These results point to a dynamically active but not chaotic system, confirming earlier observations that Derivative_39 maintains contact with the receptor but does so flexibly and adaptively.

3.13 Free energy landscape analysis

The FEL derived from PCA provides a topological view of protein-ligand complexes’ thermodynamic stability and conformational diversity. In FEL maps, basins of low Gibbs free energy (depicted in purple and blue) represent the most thermodynamically stable conformational states, while red regions reflect higher-energy, less favorable conformations (Fig. 9).

Close

Fig. 9.

FEL plots of Vitamin D Receptor (VDR) complexes with (a) Shikonin, (b) Derivative_02, (c) Derivative_12, and (d) Derivative_39 projected along the first two principal components (PC1 and PC2). The color gradient represents Gibbs free energy in kJ/mol, with purple indicating low-energy (stable) conformational basins and red representing high-energy regions. Distinct energy wells indicate the thermodynamically favorable conformational states sampled during the 500 ns MD simulation.

Export to PPT

The FEL for the Shikonin–VDR complex (Fig. 9a) reveals two well-defined and deep energy basins centered around distinct low-energy regions in the PC1–PC2 projection. These basins are separated by a narrow transition region, suggesting that Shikonin enables the receptor to adopt two central stable conformational states with limited transitions between them. The depth and symmetry of the basins indicate high thermodynamic stability and low entropic cost associated with ligand-induced conformational ordering.

In the case of Derivative_02 (Fig. 9b), the FEL shows two major deep energy wells, with one basin more populated than the other. This asymmetry reflects a dominant global minimum supported by a secondary conformational substate. The transitions between these states appear smoother than Shikonin, suggesting moderate flexibility within a constrained energetic framework. The depth and clarity of these wells underscore energetic favorability and support the structural stability previously observed in RMSD and PCA analyses.

The FEL for Derivative_12 (Fig. 9c) displays a broader, more rugged energy landscape with several shallow basins and less sharply defined boundaries. While one primary basin is evident in the upper left quadrant, its gradient toward higher-energy regions is more gradual. This dispersed energy profile indicates higher conformational entropy and reduced thermodynamic stabilization compared to Shikonin and Derivative_02. These findings are consistent with PCA results, showing increased conformational sampling and reduced PC1 variance contribution.

The FEL for Derivative_39 (Fig, 9d) reveals a multi-basin topology, with at least three notable low-energy regions interconnected by moderate-energy pathways. These features suggest that Derivative_39 permits greater structural adaptability, allowing the receptor to visit multiple metastable states during the simulation. Despite having pronounced minima, the wider energy spread and irregular basin shapes reflect elevated conformational flexibility, which aligns with the broader RMSD fluctuations and higher PCA variance spread (Table 2).

3.14 Structural superimposition of minima conformations with initial docking poses

To evaluate the conformational fidelity and binding stability of Shikonin and its derivatives within the VDR binding site, we extracted three representative minimum energy conformations from the FEL of each ligand-bound complex (Fig. 10). These conformations were superimposed with the original docking pose using UCSF Chimera, and deviations were quantified by RMSD values to assess the extent of ligand displacement and binding pose conservation throughout the simulation.

Close

Fig. 10.

Energy minima structures extracted and superimposed from the FEL for each complex: (a1-a3) Shikonin, (b1-b3) Derivative_02, (c1-c3) Derivative_12, and (d1-33) Derivative_39.

Export to PPT

Superimposition of the three minima structures of Shikonin with its initial docked pose yielded an RMSD of 1.180 Å, indicating excellent structural retention (Fig. 11a). Across all three energy minima, the ligand consistently occupies the same spatial orientation. It maintains key binding interactions, affirming that Shikonin’s pose is energetically stable and dynamically preserved throughout the 500 ns simulation. This substantial spatial overlap supports a conformationally rigid and high-affinity binding mode.

Close

Fig. 11.

Superimposition of three energy minima structures extracted from the MD simulation with the original docking pose for each compound. (a) Shikonin, (b) Derivative_02, (c) Derivative_12, and (d) Derivative_39.

Export to PPT

The three FEL minima of Derivative_12 (Fig. 11b) show a moderate but coherent deviation from the original pose, with an RMSD of 1.287 Å. Although there is observable conformational drift in some regions, particularly in flexible loops near the binding pocket, the ligand remains well-engaged. Notably, one minima shows a slight tilt in the ligand orientation, possibly reflecting adaptive hydrogen bonding to residues such as SER275 or ILE271. The overall consistency among the three minima highlights structural stability with localized adjustments.

In Derivative_02 (Fig. 11c), the superimposed minima structures show larger spatial variance than the docked conformation, with an RMSD of 1.488 Å. The three conformations exhibit ligand flexibility and rotational adaptation, especially around the side chains that interact with TRP286 and GLN317. Despite this mobility, the ligand stays confined within the binding cleft, suggesting that conformational exploration is thermodynamically favorable but does not lead to destabilization or exit from the pocket.

The superimposed minima structures of Derivative_39 (Fig. 11d) yield an RMSD of 1.480 Å, similar to Derivative_02. The conformers reveal significant spatial reorganization, where some orientations show a shift in aromatic stacking and slight repositioning of polar functional groups. These rearrangements align with FEL and PCA results that indicated multistate behavior and flexible binding. While none of the conformers deviate catastrophically, their variation underscores a more dynamic ligand accommodation process.