Theoretical study on the structures and pharmacokinetic evaluation of verticillane-type diterpenes from soft coral Heteroxenia ghardaqensis

Abstract

The density functional theory (DFT) method ωB97XD/6-311++G(2d, p) was applied to calculate and analyze the geometric structures, spectral properties, frontier molecular orbitals, and molecular electrostatic potentials of 14 novel verticillane-type diterpenoids isolated from the soft coral Heteroxenia ghardaqensis. Additionally, reaction index analysis was conducted using conceptual density functional theory, and the drug-likeness of these compounds was evaluated using two different pharmacokinetic prediction platforms. The results showed that the hydroxyl hydrogen, secondary amine hydrogen, carbonyl oxygen, and hydroxyl oxygen in the molecules of these compounds have relatively high reactivity. Compounds 5, 8, and 9 exhibit significant anti-inflammatory activity and have similar electronic delocalization distribution characteristics, showing good stability and excellent biological activity, among which compound 5 demonstrates more significant drug potential. For compounds 2, 8, and 12 with hepatoprotective activity, through the analysis of comprehensive pharmacokinetic parameters and molecular docking data, compound 12 is considered more suitable as a potential hepatoprotective drug.

Graphical abstract

Background

Marine soft corals are considered to be a natural source of bioactive compounds. Their metabolites exhibit diverse biological activities, including anti-inflammatory and antibacterial properties [1, 2], anticancer effects [3], antiviral activity [4], antimalarial and antifouling properties [5, 6], cytotoxicity [7, 8], and enzyme activity inhibition [9]. Studying soft coral metabolites is essential in marine natural product chemistry and provides a rich resource for developing new marine drugs and bioactive substances [2, 10, 11]. For instance, Jiang et al. [12] utilized various chromatographic techniques to separate and purify diterpenoid compounds from the soft coral Sinularia depressa; Pham et al. [13] applied spectroscopic methods such as nuclear magnetic resonance (NMR) and electronic circular dichroism (ECD) to determine the absolute configurations of diterpenoids and steroids in the soft coral S. brassica; Yu et al. [14] conducted structural analysis of four new xanthone-type diterpenoid compounds isolated from the soft coral Sinularia nanolobata using time-dependent density functional theory (TD-DFT); Jin et al. [15] investigated the structure–activity relationship of two new cladiellin-type diterpenoid compounds from the metabolites of the soft coral Cladiella krempfi against epidermal growth factor receptor (EGFR) inhibition through molecular docking studies. These research advancements demonstrate an active trend in studying metabolites from soft corals.

In the extensive investigation of the metabolites of marine soft corals, reports on verticillane-type diterpenoid compounds are relatively rare. Between 2002 and 2010, 26 verticillane-type diterpenoid compounds were isolated from the same soft coral Cespitularia hypotentaculata [16,17,18,19,20,21,22]. In 2012, Chang et al. [23] conducted a study on the metabolites of the soft coral Cespitularia taeniata and reported the isolation of three verticillane-type diterpenoid compounds. In natural product research, verticillane-type diterpenoids have a broad and diverse range of sources. Besides marine soft corals, these compounds are also widely found in various plants, including conifers, dicotyledons, and bryophytes [24]. In early studies, researchers successfully isolated verticillane-type diterpenoids from the frankincense tree Boswellia carteri [25, 26], and verticillol was isolated from the Japanese golden fir Sciadopitys verticillata [27]. In 2005, Fumihiro Nagashima et al. [28] isolated three verticillane-type diterpenoids with enantiomers from Jackiella javanica, a Japanese moss species. In 2020, Yadav et al. [29] discovered ent-verticillane-type diterpenoids in mosses, indicating that these secondary metabolites might play a crucial role in the adaptation process of mosses to their environment. In 2022, Batiha et al. [30] noted in a review work that verticillane-type diterpenoids are among the principal bioactive constituents of Commiphora plants. In 2023, Yuan et al. [31] isolated two new verticillane-type diterpenoids from the gum resin of Boswellia sacra. Recently, Sura and Cheng [32] once again pointed out in a review concerning natural products from medicinal plant resins that verticillane-type diterpenoids primarily originate from plant resins, with Boswellia and Commiphora plants being typical representatives. Thus, it can be seen that verticillane-type diterpenoids from marine soft corals are relatively scarce. The Egyptian Red Sea soft coral Heteroxenia ghardaqensis, belonging to the Xeniidae family [33], is abundant in structurally unique bioactive substances such as ceramides [34], sesquiterpenes and diterpenes [35], steroids [36], and diacylglycerols [37]. In 2023, Han et al. [38] successfully isolated fourteen verticillane-type diterpenoid compounds from Heteroxenia ghardaqensis, including compounds 5, 8, and 9 exhibiting anti-inflammatory activity and compounds 2, 8, and 12 showing hepatoprotective activity.

There have been few reports on the theoretical study of verticillane-type natural diterpenes. Cerda-García-Rojas et al. [39] applied the B3LYP method in density functional theory (DFT) at 6-31G (d, p) and DGDZVP levels to perform vibrational circular dichroism (VCD) theoretical calculations, determining the absolute configuration of verticillane-type natural diterpenes isolated from Bursera suntui. The density functional theory method is irreplaceable and crucial in investigating complex marine natural products. This method can precisely predict the molecular geometry, electron density distribution, and spectral properties through theoretical calculations. For instance, DFT can simulate infrared vibrational and nuclear magnetic resonance spectra, offering significant evidence for determining molecular structures; it can also ascertain the chiral features of molecules through calculations, which is paramount for comprehending their biological activities [40, 41]. Based on density functional theory, we can also calculate global and local reaction descriptors, thereby identifying the active sites of molecules. These active sites form the core basis for further exploring the interaction between molecules and biological targets. Take the research of Flores-Holguín et al. as an example. They conducted an in-depth investigation of marine peptides Discodermins A-H and Theopapuamides A-D using the DFT approach. They elaborately explored these molecules' chemical reactivity and biological activity [42, 43]. The combination of ADME/Tox parameter prediction and molecular docking studies enables the efficient screening of compounds with potential pharmacological activities. For instance, Lou et al. [44] screened 52,765 compounds from a vast marine compound database and initially selected 20 compounds through pharmacokinetic tests. Subsequently, molecular docking studies further validated the binding ability of these compounds to target proteins, ultimately identifying compound No. 50843 as a potential CDK4/6 inhibitor. The above-targeted theoretical and computational studies can improve the efficiency and potential success rate of discovering new marine drugs while reducing the randomness and costs of the experimental process.

Further research on specific soft coral metabolites is beneficial for discovering new drug candidates. Therefore, this work computationally analyzed the structures and properties of the fourteen recently reported verticillane-type diterpenoid compounds from Heteroxenia ghardaqensis, comparing differences in geometric and electronic structures and exploring the relationship between structure and their anti-inflammatory and hepatoprotective activities.

Results and discussion

Geometric structure

Figure S1 illustrates the structural formulas of fourteen recently reported verticillane-type diterpenoid compounds (1–14) isolated from the soft coral Heteroxenia ghardaqensis. The 6/12 bicyclic carbon skeleton is designated as the A ring and the C ring, respectively, while the hetero-lactone or hetero-lactam ring is defined as the B ring. These compounds' stable 3D molecular configurations and atomic numbering were obtained at the ωB97XD/6-311++G(2d, p) level, optimized in vacuum, methanol, and water environments. Figure 1 shows only the optimized geometry in a vacuum environment.

Fig. 1

The geometric structures of fourteen verticillane-type diterpenoid compounds were optimized using the ωB97XD/6-311++G(2d, p) method in a vacuum environment

The theoretical values of key bond lengths and bond angles for compounds 1, 2, 8, 9, 12, and 13, optimized in vacuum, are compared with experimental values from crystal structure analysis in Supplementary Material Table S1. Table S2 shows the main bond lengths of 14 compounds after being optimized by the ωB97XD/6-311++G(2d, p) method in vacuum, methanol, and water. The following discussion focuses on gas-phase optimized theoretical bond lengths and dihedral angles. For instance, in compound 1, the carbonyl bond length in the B ring is 1.197 Å, whereas for compounds 2, 12, 13, and 14, it is around 1.210 Å; for compounds 3–11, it ranges between 1.216 and 1.219 Å. These differences may be attributed to the unstable lactone or lactam ring effect.

The carbonyl bond lengths in the C ring of the fourteen compounds show no significant differences, averaging at approximately 1.209 Å. Similarly, the hydroxyl bond lengths in the B ring and the nitrogen atom substituent in the B ring fall within a range of 0.957–0.962 Å across all compounds. In compound 2, the N5-H25 bond length in the imidazole ring is measured at 1.009 Å, slightly longer than that of compounds 7 and 11, which have indole rings with a nitrogen–hydrogen bond length of 1.003 Å by approximately 0.006 Å. The C1–N5 bond lengths for compounds 2–14 range from 1.365 to 1.386 Å, indicating partial double-bond character falling between carbon–nitrogen single and double bonds. C4-C5 bond lengths range from 1.400 to 1.462 Å, closer to the typical carbon–nitrogen single bond length (1.460 Å). Notably, there are distinct differences in dihedral angles (C8–C9–C3–C4) between A/C double-ring carbon skeletons. Compound 1 exhibits − 160.3°, whereas compounds 2–11 display angles ranging from − 130.1° to − 127.4°; and finally, 12–14 exhibit greater torsion angles for A and C rings (− 68.11°, − 55.8°, and − 55.8°), likely because compounds 1–11 possess C2–C3 double bonds while 12–14 have C4–C11 double bonds. From Fig. 1, the key bond lengths of the core carbon skeletons of fourteen compounds in methanol and aqueous environments show minimal variation. The disparities are primarily evident in the main bond lengths in vacuum and polar solvents such as methanol and water, including the C1–O10 and C1–N5, and carbonyl groups in rings B and C for compounds 2–14. Table S1 reveals that the theoretical bond lengths of the C1–O10 carbonyl in the B ring in compounds 1, 2, 4, 8, 9, and 12 differ from the experimental values by only 0.019–0.024 Å. In contrast, for compound 13, this difference is a mere 0.009 Å.

At this theoretical level, the calculated values of the carbonyl in the C ring align more closely with experimental values, differing by only 0.004–0.009 Å. When comparing the optimized bond lengths and angles to X-ray crystallographic results, the primary bond lengths deviate from -0.309 to 0.293 Å, with an average absolute deviation of 0.042 Å and an error of 3.33%. The primary bond angles exhibit a deviation range of − 6.1° to 4.8°, with an average absolute deviation of 1.4° and an error of 1.19%. Hence, it can be inferred that the geometric structural parameters obtained from theoretical calculations agree with the crystal structural parameters. Minor discrepancies may arise from differences in compound environments. The experimental results pertain to the solid phase, while the theoretical calculations were conducted in the gas phase [45, 46].

Infrared absorption spectroscopy (IR)

Infrared absorption spectroscopy can be used to analyze molecular structures, and density functional theory offers an economically efficient method for calculating the vibrational spectra of organic molecules [47, 48]. Figure 2 presents the vacuum IR absorption spectra (a–c) as well as the linear regression fitting of theoretical vibrational frequencies and experimental vibrational frequencies (d) for fourteen compounds calculated using the ωB97XD/6-311++G(2d, p) method. Table S3 provides the assignment of the IR peaks in both theoretical and experimental data and lists the absorption intensity of each IR absorption peak. By referencing a comparative computational chemistry database and benchmark database (http://cccbdb.nist.gov/vibscale.asp), we selected a frequency correction factor of 0.957 for comparable calculation methods and similar basis sets. Figure 2 shows that the range 3700–3800 cm⁻¹ absorption peak corresponds to O–H stretching vibration, and compounds 2, 7, and 11 contain amino groups in the 3400–3500 cm⁻¹ range.

Fig. 2

The IR absorption spectra (a–c) of fourteen verticillane-type diterpenoid compounds calculated in the gas phase using the ωB97XD/6-311++G(2d, p) method, along with the linear regression fitting of theoretical vibrational frequencies and experimental vibrational frequencies (d)

In the 2800–3100 cm⁻¹ range, multiple absorption peaks with similar intensities are primarily attributed to the stretching vibration of C-H bonds. Compounds 1–11, 13, and 14 exhibit C=O stretching vibrations in the C ring at approximately 1726 cm⁻¹, while compound 12 shows this vibration at 1733 cm⁻¹. The C=O stretching vibration in the B ring is maximal in compound 1 (1780 cm⁻¹) and minimal in compound 7 (1705 cm⁻¹), ranging from 1710 cm⁻¹ for compounds 3–6. Compounds 2 and 12–14 display this position's absorption peaks around 1735 cm⁻¹. All compounds contain a C=C bond in either the B or C ring, with absorption peaks around 1670 cm⁻¹ for compounds 1–12 and 1659 cm⁻¹ for compounds 13 and 14. Absorption peaks below 1000 cm⁻¹ are attributed to bending vibrations of C-H bonds.

Upon comparing compounds 1 and 2, it was observed that the C=O stretching in the B ring exhibited differences, with compound 2 displaying a blueshift of approximately 40 cm⁻¹. However, the C=O stretching vibration and the O–H stretching vibration of the B ring were broadly consistent between both compounds. This analysis suggests compound 2 contains p–π conjugation in its lactam ring, reducing the double bond properties of C=O and causing an absorption peak shift towards lower wavenumbers. Comparison of compounds 12, 13, and 14 reveals that the O–H stretching vibrations (3735 cm⁻¹, 3733 cm⁻¹, and 3734 cm⁻¹) and C=O stretching vibrations (1740 cm⁻¹, 1739 cm⁻¹, and 1739 cm⁻¹) in the B ring are nearly identical. However, the C=O stretching vibrations in the twelve-membered ring of compounds 13 and 14 are blue-shifted by 9 cm⁻¹ compared to compound 12.

Based on Table S3 and Fig. 2d, it can be observed that R² = 0.988, which is greater than 0.950, indicating a high degree of agreement between the theoretical and experimental results using this method. The characteristic absorption peaks of the theoretical data and experimental data show poor alignment in the high-frequency range (2500–3700 cm⁻¹), particularly within the range of 3400–3800 cm⁻¹; however, they exhibit good alignment in the low-frequency range (below 1700 cm⁻¹). This may be attributed to the theoretical calculation simulating individual molecular behavior. At the same time, the experimental results represent the collective behaviors of molecules in solid powder form involving various intermolecular interactions.

Ultraviolet–visible absorption spectrum (UV–Vis)

The ωB97XD method, which incorporates dispersion and long-range corrections, provides a more accurate calculation of the ultraviolet absorption spectra of organic molecules and offers valuable spectral information [49]. The UV–Vis absorption spectra curves for fourteen compounds calculated using the ωB97XD/6-311++G(2d, p) method in methanol and water are presented in Fig. 3. In the range of 150–300 nm, compounds 9, 10, and 11 exhibit a single absorption peak, while compounds 1, 3, and 4 display three peaks; the remaining compounds show two peaks. Notably, compounds 7 and 11 lack an absorption peak in the far ultraviolet region (150–200 nm); compound 1 exhibits two peaks within this range, while others have one peak. In the middle ultraviolet region (200–300 nm), compounds 9 and 10 do not exhibit any absorption peaks; however, compounds 3, 4, 7, and 11 show two peaks each, with others displaying only one. Solvent effects have minimal impact on the intensity and wavelengths of absorption peaks across all studied compounds. The nearly overlapping nature of the absorption spectra curves in methanol and water may be attributed to limitations arising from a simplified computational solvent model that fails to provide clear differentiation. Referring to Table 1, the primary maximum absorbance peak for compound 7 is attributed to electron transition from HOMO orbital to LUMO + 2, with a contribution rate reaching 93.3%. In contrast, compound 3 is due to an electron transition from the HOMO-2 orbital to the LUMO with a contribution rate of 88.7%.

Fig. 3

Ultraviolet–visible absorption spectra of fourteen compounds in methanol and water environments calculated using the ωB97XD/6311++G(2d, p) method

Table 1 Ultraviolet–visible absorption spectra of fourteen compounds in methanol calculated with ωB97XD/6-311++G(2d, p) method

Upon comparison, it was observed that compound 1 (222 nm) exhibited a red shift in UV wavelength compared to compound 2 (202 nm), which aligned with the experimental observations for both compounds. This can be attributed to the presence of ester groups containing polar bonds in compound 1, which act as electron-withdrawing groups. Compounds 3 and 4 displayed UV absorption wavelengths of 203 nm, closely matching their respective experimental data (both at 213 nm), indicating good agreement. The substitution at the terminal end of the nitrogen atom in the lactam ring with either a hydroxyl group or an isopropyl group has no impact on the absorption wavelength. Compounds 3 (203 nm) and 8 (204 nm) showed similar UV absorption wavelengths, suggesting that substituting a hydroxyl group adjacent to the nitrogen atom in the lactam ring does not affect the absorption wavelength. Similar phenomena were observed between compounds 5 and 9 and 7 and 11. However, it was found that compound 10 exhibited a red shift of 23 nm compared to compound 6; moreover, experimental data reported only a 9 nm red shift for compound 10. Compounds 12, 13, and 14 exhibit absorption peaks at wavelengths of 220 nm, 209 nm, and 218 nm in the middle-ultraviolet region. Compound 14 displays a characteristic absorption band resulting from the π–π* transition of the benzene ring, which is influenced by the phenolic hydroxyl group, leading to p-π conjugation and the disappearance of fine structure. Consequently, compared to compound 13, there is an increase in absorption intensity and a redshift in the absorption peak. Compounds 5 and 6, as well as 9 and 10, exhibit similar trends between theoretically calculated absorption spectra and experimental spectra.

Molecular surface electrostatic potential maps (MSESPMs)

The molecular surface electrostatic potential (MSESP) is a method used to characterize the reactive sites of molecules by identifying potential sites for nucleophilic and electrophilic attacks [50, 51]. Yu et al. [52] conducted molecular docking and electrostatic potential analysis on eight non-prodrug third-generation cephalosporins, and the results indicated a very high degree of alignment between them. Figure 4 illustrates the distribution of extreme molecular surface electrostatic potential points and the area distribution for six verticillane-type diterpenoid compounds, and the electrostatic potential diagram of the remaining eight compounds is shown in Supplementary Materials Fig. S2. The blue region represents areas with a negative charge, making them susceptible to electrophilic reagents; the red region represents areas with a positive charge, making them susceptible to nucleophilic reagents; and the white indicates regions with zero electric potential. Among these regions, yellow spheres correspond to maximum electrostatic potential points, while green spheres correspond to minimum values.

Fig. 4

Electrostatic potential maps of six compounds (compounds 1, 3, 5, 6, 10, and 14) calculated using the ωB97XD/6-311++G(2d, p) method in the gas phase

Figure 4 and Fig. S2 show that most of the maximum points of electrostatic potential for fourteen compounds are located near the hydroxyl and secondary amine hydrogen atoms. In contrast, the minimum points are mainly near the carbonyl and hydroxyl oxygen atoms. All compounds are susceptible to attack by electrophilic reagents at positions adjacent to the carbonyl oxygen in ring B, with compound 1 exhibiting a minimum value of -39.71 kcal/mol at this position, which is smaller than the minimum values of other compounds (− 45.12 to − 40.02 kcal/mol). The maximum value (48.46 kcal/mol) of compound 1 near the epoxide atom in the lactone ring is greater than that of compound 2 (44.61 kcal/mol) near the nitrogen atom in the succinimide ring. Compound 3 exhibits a higher maximum value (46.50 kcal/mol) for hydroxyl hydrogen in the ring B compared to substituted hydroxyl hydrogen in the ring B–H with a maximum value of 44.74 kcal/mol, making it more susceptible to nucleophilic reagent attacks. Compounds 3 and 5–11 have minimum values near carbonyl oxygen in ring B, ranging from − 40.96 to − 45.00 kcal/mol; compared with compounds where there is a hydroxyl substitution adjacent to the nitrogen atom in ring B, those without such substitution are more vulnerable to attack by electrophilic reagents. The maximum values of the hydrogen atoms near the phenolic hydroxyl groups in compounds 6, 10, and 14 are the highest, at 52.67 kcal/mol, 52.41 kcal/mol, and 51.37 kcal/mol, respectively.

Frontier molecular orbital analysis (FMO)

By calculating and analyzing the frontier molecular orbitals, it is possible to predict molecules' electron-donating and accepting abilities. Using density functional theory ωB97XD/6-311++G(2d, p), a frontier molecular orbital analysis was performed on fourteen verticillane-type diterpenoid compounds. However, Fig. 5 only shows the six compounds; the remaining eight are shown in Fig. S3 in the Supplementary Materials. Compounds 7 and 11 exhibit relatively low energy gap differences compared to other molecules, with values of 8.21 eV and 8.18 eV, respectively, indicating high occupied orbitals and solid electron-donating capabilities. Compounds 1–5, 8, and 9 display larger energy gap differences ranging from 9.23 eV to 9.62 eV, suggesting poor electron acceptance ability and good molecular stability [53]. Compounds 6, 10, and 12–14 have molecular energy gap differences ranging from 8.66 eV to 8.90 eV.

Fig. 5

Frontier molecular orbital diagram of six verticillane-type diterpenoid compounds (4–9)

Compound 1's lowest unoccupied molecular orbital (LUMO) is primarily distributed on the lactone and hexagonal rings (B-ring). At the same time, the highest occupied molecular orbital (HOMO) is mainly localized on the carbonyl group and two methyl groups of the dodecagonal ring (C-ring). The frontier molecular orbital distribution of compound 2 resembles that of compound 3, with both LUMO and HOMO predominantly located on the lactam ring and B-ring. For compounds 4–7 and 12–13, the LUMO consists mainly of the carbonyl group in the C-ring. In contrast, for compounds 8–11, the LUMO is primarily delocalized on a methyl group near the lactam ring and on a carbonyl group in the C-ring. The HOMO orbitals of compounds 7 and 11 are mainly delocalized onto a nitrogen atom substituent at one end of an indole moiety attached to a C-ring. Compounds 4, 5, 8, and 9 exhibit HOMOs that are delocalized over both lactam rings, hexagonal rings, and nitrogen atom substituents. The HOMO orbitals of compounds 6, 12, and 14 are mainly delocalized over the lactam ring and the nitrogen-substituted ring, particularly the phenol substituent. For compound 10, its HOMO orbitals are primarily delocalized over the phenol substituent at the end of the nitrogen-substituted ring, with a smaller portion delocalized over the lactam ring. As for compound 13, its HOMO orbitals are mainly delocalized over the lactam ring and the methyl group near this ring.

From the perspective of frontier molecular orbitals, compounds 5, 8, and 9 exhibit good molecular stability with anti-inflammatory activity. Further analysis of the electron delocalization distribution in the molecular orbitals of the fourteen compounds reveals that compounds 5, 8, and 9 share similar distributions of LUMO and HOMO.

Global reactivity descriptors

The global reactivity indices for fourteen verticillane-type diterpenoid compounds are shown in Table 2. The ionization potentials of these compounds range from 7.30 eV to 8.90 eV, and their electron affinities in vacuum are all negative, indicating a weak overall binding ability to electrons and a limited tendency to gain electrons. The parameters μ and η can reflect the stability of molecules; compounds 1 and 2 exhibit smaller μ values (− 4.36 eV and − 4.13 eV) and larger η values (4.51 eV and 4.60 eV), suggesting that these two compounds are more stable than the other twelve compounds. While the S values of the fourteen compounds show minimal differences, compound 1 has higher χ and ω values, indicating a relatively stronger electron affinity with high electrophilic reagent binding capacity, resulting in elevated reactivity levels. Compound 11 displays the lowest η value (4.12 eV) and ω value (1.36 eV), signifying poor stability and reactivity. Compounds 5 and 9 differ only in their C4 substituents, and the analysis shows that these two compounds' η and S values are consistent, but the μ value of 5 is smaller. The ω value is more significant, indicating that 5 is more stable and has higher activity than 9. Upon comparing compounds 8, 9, 10, and 11, it was observed that compound 8, featuring a nitrogen atom connected to a hydroxyethyl group, exhibited superior stability and activity compared to the other three compounds. Similarly, upon comparing compounds 2–7, compound 2 demonstrated better stability and activity than the rest. Compound 7, with a nitrogen atom substituted with an indole group, displayed the minimum stability and activity among all compounds. Further analysis of compounds 12, 13, and 14 further revealed that their S and ω values were identical. However, the μ value of compound 13 was smaller while its η value was more considerable, suggesting more excellent stability compared to the other two compounds. According to the global reactivity descriptors, compound 1 is more stable and has higher reactivity; compounds 5, 8, and 9, which have anti-inflammatory activity, are located in the middle of the fourteen compounds in terms of activity and stability; compounds 2, 8, and 12 with hepatoprotective activity, 2 has better stability and activity, and 12 is poor. According to Table 2, introducing a hydroxyl group into the heterocyclic imide C4 structure of structurally similar compounds may increase their activity; introducing an indole group may lead to a decrease in the stability and activity of the compounds.

Table 2 Global reactivity indices of fourteen verticillane-type diterpenoid compounds at the level of ωB97XD/6-311++G(2d, p)

Local reactivity descriptors

The maximum value of the local reactivity descriptor represents the site on the molecule that is more susceptible to nucleophilic or electrophilic attack than other atoms [54, 55]. Further analysis was conducted on the condensed Fukui functions (f⁺, f ⁻, f ⁰), softness (s⁺, s⁻), and electrophilic index (ω⁺, ω⁻) of the fourteen compounds, and their extrema are listed in Table S4. A comparison between compounds 1 and 2 reveals that the f⁻, s⁻, and ω⁻ indices of heterocyclolactone O5 are more significant than those of heterocyclonitrile N5. In contrast, the f⁺, s⁺, and ω⁺ indices of heterocyclonitrile H31 are greater than those of heterocyclolactone H30. According to Table S4, compounds 2–5, 8, 9, and 11 exhibit the highest f⁻, s⁻, and ω⁻ indices at positions N5/N27. Additionally, the f⁻, s⁻, and ω⁻ indices of the carbonyl oxygen (O10/O11/O22) in rings B and C are relatively large, indicating that these two sites are most susceptible to electrophilic attack. Compounds 6 and 10 show relatively large f⁻, s⁻, and ω⁻ indices for the hydroxyl oxygen (O33/O32) in the benzene ring. Meanwhile, compounds 12 and 13 display relatively large f⁻, s⁻, and ω⁻ indices at position O10 in ring B. Compounds 2 and 3 demonstrate the most significant f⁺, s⁺, and ω⁺ indices for hydrogen atoms attached to hydroxyl groups in the ring B (H31/H33), suggesting susceptibility to nucleophilic attack. Compounds 4–6, 8–10, and 12–14 have maximum f⁺, s⁺, and ω⁺ indices for hydrogen atoms adjacent to carbonyl carbon in ring C, while compound 11 exhibits maximum f⁺, s⁺, and ω⁺ indices at position H69. Upon comparing the hydrogen atoms (H37/H39/H42 on rings 3, 4, and 5) adjacent to the carbonyl group in the twelve-membered ring, it was found that H42 on compound 5 was more susceptible to nucleophilic attack. Upon comparing the hydrogen atoms (H38/H40/H43 in rings 12, 13, and 14) adjacent to the carbonyl group in the twelve-membered ring, it was found that H40 on compound 13 was more susceptible to nucleophilic attack. Compounds 5 and 13 feature a phenyl group connected at the terminal end of a nitrogen-substituted alkyl group. Comparing compounds 3 and 8, 5 and 9, 6 and 10, it was found that the nitrogen atom (N5) at the imidazole ring with hydroxyl substituent in 3, 5, and 6 had smaller f⁻, s⁻, and ω⁻ indices than those in 8, 9, and 10 without hydroxyl substituent. Compounds 7 and 11 differ only in the substituent at C4, and it was found that the f⁺, s⁺, and ω⁺ indices at H69 in 11 were more significant than those at the corresponding atom in 7, while the f⁻, s⁻, and ω⁻ indices at N27 in 11 were smaller than those in 7.

Pharmacokinetic evaluation

The ACD/Percepta software is based on quantitative structure–activity relationship (QSAR) models for predicting the absorption (A), distribution (D), metabolism (M), excretion (E), toxicity (T), and physicochemical properties of compounds. According to Table S5, all verticillane-type diterpenoid compounds comply with Lipinski's rule, indicating good drug-likeness. However, except for compounds 1, 2, and 3, the solubility of the remaining compounds is poor. Pharmacokinetic evaluation was performed for these fourteen compounds, and the predicted results of relevant parameters are presented in Table 3.

Table 3 ACD/Percepta software predicts pharmacokinetic parameters for fourteen verticillane-type diterpenoid compounds

The fourteen compounds are not substrates of P-glycoprotein (P-gp). They exhibit minimal inhibition of the major liver metabolic enzyme CYP1A2 subtype, with compounds 10 and 11 showing slightly more potent inhibition than the others. Compound 11 demonstrates the most potent inhibition of the CYP2C9 subtype, followed by compounds 10 and 7. Compounds 3 and 12 show the weakest inhibition of the CYP2C19 subtype. Compounds 7, 11, and 14 significantly inhibit the CYP3A4 subtype. Results from the Ames test and human Ether-à-go-go-Related Gene (hERG) assay suggest that these compounds may have mild genotoxicity and cardiotoxicity, with compound 11 exhibiting the highest cardiotoxicity. The carcinoma of colon adenocarcinoma-2 (Caco-2) data indicate good absorptive properties for all fourteen compounds, with compound 11 showing lower absorption (33.3%) than the others, while the compound is also poorly absorbed. Compounds 1, 8, and 12 demonstrate excellent absorptive properties. Plasma protein binding rate (PPB) significantly influences drug therapeutic effects; among these compounds, compound 1 exhibits weaker plasma protein binding ability than the others do, whereas compound 12 comes in second place. All fourteen compounds show poor permeability across the blood–brain barrier, resulting in low concentrations within the central nervous system (CNS). These fourteen compounds demonstrate similar metabolic stability profiles and human intestinal absorption (HIA) performance.

Among the eighteen available network platforms for free prediction of physicochemical and pharmacokinetic parameters, the ADMETLab platform offers the broadest prediction range and the most precise results [56]. ADMETLab 3.0 (https://admetmesh.scbdd.com/) represents the latest enhanced version of ADMET Lab 2.0 [57]. Therefore, we utilized the ADMETLab 3.0 platform to predict the ADMET properties of fourteen verticillane-type diterpenoid compounds, with specific results detailed in Table S6. The results indicate no significant difference in the absorption of fourteen compounds in the human intestine. These compounds do not inhibit CYP1A2 and CYP2D6 subtypes or act as hERG blockers. Compounds 2 and 7 demonstrate poor Caco-2 permeability. Compound 6 has the highest probability of being a P-gp substrate, followed by compound 2. Drugs with high protein binding rates may have a low therapeutic index, with compounds 4–7, 9–11, 13, and 14 exhibiting binding rates exceeding 90%. Compounds 3, 4, and 7–11 show substantial liver toxicity, while compounds 2 and 12 exhibit the weakest liver toxicity. All compounds have high nephrotoxicity, with 3 potentially causing the most tremendous damage to the kidneys, while 4, 7, and 11 have strong ototoxicity. The genetic toxicity of 6, 7, and 14 is substantial, while 1 and 2 are relatively safe. 1 has the most prolonged half-life, followed by 2 and 14. Compounds 4, 7, 9, and 11 significantly inhibit CYP2C19, CYP2C9, CYP3A4, CYP2B6, and CYP2C8 subtypes. Compound 8 has the best stability in human liver microsomes, while the rest are unstable, with a 100% probability of instability for 4–7 and 9–11.

Combining the pharmacokinetic parameter prediction results from two different platforms, compound 11 exhibits solid inhibitory effects on specific subtypes of cytochrome P450. Long-term abuse of liver-protecting drugs carries the risk of drug-induced liver injury, which may increase the liver's metabolic burden [58]. Compared to 8, 2 and 12 have lower liver toxicity. Liver metabolic stability is crucial for drug discovery [59], and the human liver microsomal stability of 8 is the best. At the same time, 2 has low genetic toxicity but poor absorption, and 12 has good Caco-2 permeability.

Molecular docking

The potential anti-inflammatory mechanism of active diterpenoid compounds may involve their binding to inducible nitric oxide synthase (iNOS)/cyclooxygenase-2 (COX-2) [60, 61]. In the CuSO₄-induced transgenic fluorescent zebrafish inflammation model experiment, compounds 1–14 were evaluated for their anti-inflammatory activities with 20 μM indomethacin as the positive control. The experimental results indicated that compounds 5, 8, and 9 exhibited moderate anti-inflammatory activities [38]. To further elucidate the possible binding modes of compounds 5, 8, and 9 in exerting anti-inflammatory activity, molecular docking was conducted by retrieving iNOS (PDB: 3E6T) and COX-2 (PDB: 1PXX) proteins from the PDB protein database (https://www.rcsb.org). The results presented in Table 4 demonstrate that compound 5 can effectively bind within the binding pockets of iNOS and COX-2, with respective binding free energies of − 8.58 kcal/mol and − 3.62 kcal/mol, confirming its strong affinity. In the binding pocket of iNOS, compound 5 primarily forms stable hydrogen bonds with Gln257 and Glu371 while also engaging in hydrophobic interactions with four other amino acid residues; its interaction with COX-2 is predominantly driven by hydrophobic forces involving a total of sixteen amino acid residues.

Table 4 The binding free energy (kcal/mol) and corresponding targeting residue distance (Å) between five compounds and protein active sites

Zong et al. [62] discovered that tumor necrosis factor TNF-α (PDB: 6X82) was a crucial target protein for Hendersin B methyl ester in combating various liver injuries caused by viral hepatitis. In the zebrafish liver injury model induced by isoniazid and alcohol, compounds 1–14 were evaluated for their liver injury protective effects with 50 μM S-adenosyl-l-methionine (SAM) as the positive control. The experimental results indicated that compounds 2, 8, and 12 exhibited moderate hepatoprotective activities, similar to the positive control SAM [38]. Therefore, molecular docking was conducted with tumor necrosis factor-α (TNF-α) as the target protein to further elucidate the potential binding modes of compounds 2, 8, and 12 in exerting hepatoprotective activity. The results in Table 4 demonstrate that compound 12 exhibits the lowest binding free energy at − 11.18 kcal/mol and forms interactions through hydrophobic forces with amino acid residues Leu57, Tyr59, Tyr119, Tyr151, Gln61, Ile155, and Leu157. Detailed three-dimensional binding modes of the ligands with the macromolecular protein are shown in Fig. 6, and part of the binding model is shown in supplementary data Fig. S4.

Fig. 6

Schematic representation of the binding modes of the verticillane-type diterpenoid compounds with sizeable molecular protein (a iNOS; b COX-2; c TNF-α)

Conclusions

In this work, density functional theory was employed, with the ωB97XD functional and 6-311++G(2d, p) basis set, to optimize the molecular structures of 14 novel verticillane-type diterpenoids isolated from soft corals Heteroxenia ghardaqensis, thereby determining their stable conformations. These compounds' spectroscopic properties, frontier molecular orbitals, and electrostatic potential distributions were thoroughly calculated and analyzed. Simultaneously, in conjunction with conceptual density functional theory, the reactive sites of the compounds were precisely predicted. Through the ACD/Percepta software and ADMETLab 3.0 platform, the pharmacokinetic (ADME) and toxicological (Tox) properties and the drug-likeness of these compounds were systematically predicted. Comprehensive analysis indicated that compounds 5 and 12 exhibited excellent molecular stability and biological activity, presenting remarkable drug potential.

Methods

Theoretical calculations were performed using the ωB97XD method of density functional theory at the 6-311++G(2d, p) basis set level to optimize the geometric and electronic structures of fourteen verticillane-type diterpenoid compounds depicted in Fig. S1 within a vacuum environment. Subsequently, infrared (IR) vibrational spectroscopy calculations were carried out. There have been reports in the literature that show the ωB97XD method obtains more stable geometric configurations [63,64,65], exhibits minimal error in computing ultraviolet–visible (UV–Vis) spectra [66], and achieves high accuracy in calculating nuclear magnetic resonance (NMR) spectra [67]. Verticillane-type diterpenoids possess complex polycyclic structures and a wide variety of functional groups. There may be some intramolecular hydrogen bonds and hydrophobic interactions within the molecules. The 6-311++G(2d, p) basis set, which includes polarization and diffuse functions, can accurately describe the electronic structure and various weak interactions within the molecules. Cramer elaborated on the significant advantages of the 6-311++G(2d, p) basis set in enhancing computational accuracy in his book "Essentials of computational chemistry: Theories and models" [68]. This basis set has been widely used in existing computational research literature. For instance, Castillo et al. [69] employed the 6-311++G(2d, p) basis set in their study of the conformational isomerism of trans-3-methoxycinnamic acid, successfully and accurately describing the molecular conformation and precisely calculating the spectral parameters, providing crucial data support for a deeper understanding of the molecule's properties. Aiswarya et al. [70], during their exploration of boron-nitrogen nanostructures as carriers for melphalan drugs, employed this basis set to precisely optimize the drug molecule structure and accurately calculate the energies of frontier molecular orbitals and thermodynamic parameters, significantly promoting the research on the performance of drug carriers. Furthermore, the 6-311++G(2d, p) basis set has exhibited relatively high accuracy in calculating nuclear magnetic resonance spectra, among others [71,72,73,74], furnishing more reliable outcomes for related studies.

Based on the optimized structure, the conductor-like polarizable continuum model (CPCM) within the self-consistent reaction field approach was employed to simulate solvent environments, including water (Wat) and methanol (Met). This was followed by UV–Vis absorption spectroscopy calculations using time-dependent density functional theory at the same basis set level. Combining the ωB97XD method with CPCM solvent modeling effectively analyzed molecular interactions between polyacrylamide hydrogels and antimalarial drugs [75]. Xu et al. [76] applied ωB97XD/6-311++G(2d, p)/CPCM for systematic computational studies on the effects of brominated imidazole-based ionic liquids on cytosine structure and properties.

The wave function of the gas-phase stable structure was adopted for electronic structure analysis and conceptual density functional theory (CDFT) series index analysis. This includes predictions of molecular electrophilicity and nucleophilicity, as well as calculations of overall reaction properties descriptors, including ionization potential (IP), electron affinity (EA), electronegativity (χ), overall hardness (η), chemical potential (μ), overall electrophilicity index (ω) and overall softness (S) based on defined formulas in reference literature [77, 78]. Pharmacokinetic calculations were conducted using ACD/Labs Percepta software and the ADMETLab 3.0 platform. All quantitative calculations and analyses were performed using the Gaussian 16 program [79], Multiwfn 3.8 software [80], and VMD software [81]. Molecular docking was carried out using AutoDock software, and the results were displayed using the protein–ligand interaction analysis platform PLIP 2021 [82]. The results were further visualized using Pymol. In molecular docking studies, PrankWeb (https://prankweb.cz/) was applied to rapidly and accurately predict protein ligand-binding sites [83].