The crystal structure of interleukin-8 (IL-8) was obtained from the protein data bank (PDB ID: 5D14) and subjected to energy minimization (Powers et al. 1993). The crystal structure of the protein was examined for missing atoms, alternative locations, the addition of water molecules, the presence of multiple molecules, chain breaks, and other anomalies. Polar hydrogen atoms were added, and Kollman united atom charges were assigned. The chemical structures of the small molecules (ZINC000409414675, ZINC000888088617, ZINC000888090135, ZINC000001893410, ZINC000022339916, and ZINC000012444335) used in the study were retrieved from the ZINC database (https://zinc.docking.org/) (Hikmawati et al. 2022). The preparation process involved the use of ChemDraw Professional 16.0 and Chem3D 16.0 software (Lumbanraja et al. 2023). The three-dimensional molecular structures were designed and refined using the Quantum ESPRESSO v.7.1 program with density functional theory (DFT) methodology. Becke’s potential function theory (B3LYP) with 3 parameters and the 3–21G basis set were employed for calculations (Giannozzi et al. 2020).

The chemical structure identification of the small molecules utilized in the study was conducted using chemical analysis by accessing the prediction of activity spectra for substances (PASS) web server, accessible through the URL http://www.way2drug.com/passonline/ (Lagunin et al. 2000). The PASS platform employs an algorithm capable of identifying interactions between chemical compounds and the target protein under examination. This tool assesses the potential of a chemical compound to act as an inhibitor or agonist for the specified target protein. To perform PASS predictions, the molecular structure data of the compound is input using its SMILES or InChI representation. Subsequently, the platform predicts the compound’s activity spectrum based on the probable activity (Pa) and probable inactivity (Pi) values, aiming to evaluate its potential activity against the target protein (Islam et al. 2022). The Pa value represents the likelihood of the chemical compound demonstrating the desired activity on the target protein. A higher Pa value indicates a greater likelihood of the compound displaying the desired activity. This Pa value serves as an indicator of the potency of the chemical compound as an inhibitor or agonist for the tested target protein. On the other hand, the Pi value signifies the likelihood of the chemical compound exhibiting undesirable activity on the target protein. A higher Pi value suggests an increased possibility of the compound displaying undesired activity. This Pi value serves as an indicator of potential interactions between chemical compounds and unintended target proteins, as well as potential side effects that might arise.

Molecular dynamics simulations and molecular docking were conducted using the Ubuntu 23.04.5 LTS operating system on a Dell workstation equipped with an Intel Core i5–11400 CPU (12 MB cache, up to 4.40 GHz) processor and an NVIDIA GeForce GTX 1080 Ti graphics card. The system had 64 GB of RAM, a 256 GB SATA HDD, and a 256 GB SSD NVME. The docking process involved AutoDock 4.2 from MGLTools (https://ccsb.scripps.edu/mgltools/) for blind docking of all small molecules with IL-8 (Forli et al. 2012). The grid dimensions for X, Y, and Z coordinates were set at 64, 60, and 60, respectively, with a central grid placement at 2.390, 17.146, and –3.842. The grid spacing was maintained at 1.00 Å with an exhaustiveness of 8. Molecular dynamics simulations were carried out using GROMACS 2016.3 (Aragones et al. 2013). To visualize, assess, and analyze molecular dynamics trajectories, various software tools were employed, including PyMOL 3.0 (Delano 2002), QtGrace (Adnan et al. 2022), Visual Molecular Dynamics (VMD) v.1.9.4 (Humphrey et al. 1996), and Discovery Studio Visualizer v.2021 (BIOVIA 2005). The simulations of IL-8-ZINC complex molecular dynamics (200 ns) were performed using the AMBER-FF14SB-ILDN force field within GROMACS 2016.3, employing molecular mechanics at a temperature of 300 K (Gül and Pehlivan 2018). Sampling of conformations occurred at intervals of 10 ps during the 200 ns simulations for the IL-8-ZINC complex. The GROMACS 2016.3 utilities, gmx rms, gmx rmsf, gmx gyrate, gmx sasa, gmx rdf, and gmx hbond, were utilized for trajectory analysis, including root mean squared deviation (RMSD), root mean squared fluctuation (RMSF), radius of gyration (Rg), radial distribution function (RDF), surface area solvent accessible (SASA), secondary structure determination, and hydrogen bond (HBond) analysis. QtGrace was employed to create graphs and figures for visualization and presentation of the obtained data.

The process of calculating interaction-free energy involved employing the Molecular Mechanics/Poisson Boltzmann Surface Area (MM-PBSA) method, which is a well-established approach for analyzing molecular interactions within biological systems (Tuccinardi 2021). The MM-PBSA calculation was facilitated using specific molecular dynamics scripts tailored for this purpose. To compute the binding free energies within the MM-PBSA framework, the g_mmpbsa script provided by the GROMACS 2016.3 software package was utilized. This approach allows for a comprehensive assessment of the energetics of interactions between biomolecules and contributes to a deeper understanding of their behavior and stability.

Molecular docking was executed using IL-8, extracted from the final trajectory of a 200 ns molecular dynamics simulation, in conjunction with cathelicidin (LL-37). The crystal structure of LL-37 was sourced from the Protein Data Bank (PDB ID: 2K6O) and prepared accordingly (Wang 2008). The assessment of protein-peptide interactions involved employing the PatchDock server (http://bioinfo3d.cs.tau.ac.il/PatchDock/) (Schneidman-Duhovny et al. 2005). The grouping procedure adhered to the default root mean square deviation (RMSD) with a maximum threshold of 2.00 Å. The chosen complex type was the protein-peptide configuration. PatchDock utilizes well-suited algorithms to generate surface representations of molecules, presenting them in convex, concave, and planar forms. Subsequent optimization is applied to small-sized molecules within PatchDock, resulting in a collection of stable complexes. During this process, adjustments are made to the relative orientation of molecules, limiting the flexibility of surface side chains involved in interactions while allowing controlled motion for small, rigid entities. The top 10 candidate solutions underwent a re-evaluation based on the PatchDock score and atomic contact energy (ACE) score (Li et al. 2022). The identification of molecular protein-peptide docking outcomes hinged on the scores attained through this approach. The visual representation of interactions and binding within the docked conformation was facilitated using the Discovery Studio Visualizer v.2021 software tool.

Toxicity prediction is a procedure employed to assess the potential toxic effects of the molecular composition of small compounds before embarking on clinical trials involving humans or animals. Among the methodologies employed for toxicity prediction, one notable approach is the preADMET 2.0 (prediction of ADME and TOXicity) server (https://preadmet.webservice.bmdrc.org/) (Eurtivong et al. 2019). PreADMET 2.0 serves to assess the potential toxicity of bioactive compounds by scrutinizing their pharmacokinetic characteristics. These characteristics encompass the absorption, distribution, metabolism, excretion, and interaction of small-molecule compounds with receptor proteins. The preADMET 2.0 method employs statistical techniques and machine learning algorithms to analyze the pharmacokinetic traits of molecular compounds.

The crystal structure of Interleukin-8 (IL-8), obtained from the Protein Data Bank (PDB ID: 5D14), was subjected to energy minimization to optimize its conformational stability. This process was carried out using the GROMACS 2016.3 software, applying the AMBER-FF14SB-ILDN force field parameters. The goal was to mitigate any excess energy present in loops, sheets, and bonds within the protein structure. Subsequently, the optimized structure was superimposed using Discovery Studio Visualizer v.2021 to visualize any deviations in comparison to the initial structure. The root mean squared deviation (RMSD) was calculated to quantify the structural changes, resulting in a calculated RMSD value of 1.260 Å. This value indicates that the structure underwent minimal conformational alterations during the energy minimization process (Fig. 1). This assessment highlights the stability of the structure after energy minimization.

Figure 1. 

The cartoon depiction of IL-8 (PDB ID: 5D14) is presented on the left side, represented in green. On the right side, the superimposed IL-8 structure is displayed, with the original structure shown in green and the energy-minimized structure in red.

Biological activity prediction entails assessing the potential biological effects of bioactive compounds prior to initiating clinical trials in humans or animals. This process aims to identify compounds that might exhibit desired biological activities before embarking on more intricate and costly experimental evaluations. The biological activity of each bioactive compound is forecasted using the Prediction of Activity Spectra for Substances (PASS) web server. Based on the outcomes of the biological activity prediction in Table 1, it is evident that only a select few small-molecule compounds exhibit biological activity against IL-8, including ZINC000022339916 and ZINC000012444335. As such, the ZINC compounds are postulated to possess activity against IL-8. IL-8 is a proinflammatory cytokine with a role in stimulating immune responses and inducing inflammation within the body. The forecasted biological activity serves as an invaluable initial step in guiding the selection and development of potential candidate compounds for medicinal and therapeutic applications.

All the investigated small-molecule compounds in this study also exhibited a likelihood of activity against various interleukins, with their biological activity indicated by the dominant Pa values over the Pi values. Moreover, the Pa values of the compound molecules fall within the range of 0.104 to 0.482. The notion of potential activity is employed to assess the inherent capability of small compound molecules to induce desired biological effects. Conversely, the Pi values of the compound molecules range from 0.003 to 0.065. The concept of inactivity potential pertains to instances where a small compound molecule fails to display activity or does not function effectively.

The molecular docking analysis was employed to assess the binding affinity, binding energy, and potential intermolecular interactions between IL-8 and the remaining amino acids. This analysis aimed to evaluate the distances and bond conformations of these interactions. By combining all the listed small molecules with IL-8, we examined their binding affinity. The comprehensive details of the ZINC compound molecules, including binding energy and inhibition constant, are presented in Table 2. Notably, a high binding affinity was observed for certain small molecules, such as ZINC000001893410 (–7.35 kcal/mol) and ZINC000022339916 (–7.58 kcal/mol), which also displayed significant interactions with conventional amino acids. The investigation revealed that these two ZINC compound molecules acted as robust inhibitors of IL-8, as outlined in Table 2. The molecular structures of ZINC000001893410 and ZINC000022339916 are illustrated in Fig. 2.

Figure 2. 

Structures of small molecules: ZINC000001893410 (on the left) and ZINC000022339916 (on the right).

In Fig. 3, a visual representation illustrates how the chemical structural complexes of several different small compounds (ZINC000409414675, ZINC000888088617, ZINC000888090135, ZINC000001893410, ZINC000022339916, and ZINC000012444335) form. These complexes hinder the binding accessibility of small-molecule compounds through interactions that inhibit binding in that specific area. Subsequently, we conducted a molecular docking analysis involving IL-8 and these small compounds. The primary goal was to identify the most favorable docking position across the entire surface of the protein target (Fig. 3). The results of this analysis uncovered a notable binding pattern between the ZINC group’s small molecules and the central region of the IL-8 molecule. This finding strongly suggests a high potential for interaction between the two entities.

Figure 3. 

Illustration of the IL-8-ZINC compound complex portrayed using a surface depiction. ZINC’s compact molecules are visualized as colored projections. A closer examination of the IL-8 binding region is displayed, highlighting pivotal amino acid constituents responsible for establishing connections with the inhibitory ZINC molecule.

Moreover, an examination was conducted on the bonding arrangement between IL-8 and molecules of the ZINC compound, specifically focusing on ZINC000001893410 and ZINC000022339916. The findings indicate the presence of numerous potent amino acid residues within IL-8, capable of establishing hydrogen linkages with ZINC000001893410 and ZINC000022339916 at optimal bond distances (Fig. 3). Ultimately, it was deduced that ZINC000001893410 predominantly forges the majority of interactions (typical hydrogen bonds) with IL-8’s amino acid residues LYS21, LEU23, VAL25, GLN57, and GLU61. These hydrogen bonds are arranged geometrically around IL-8’s active site. Besides hydrogen bonds, IL-8 also generates hydrophobic associations with ZINC000001893410 through LEU23, VAL25, VAL60, and LEU64.

Nonetheless, the ZINC000022339916 small molecule predominantly establishes the majority of interactions (typical hydrogen bonds) with IL-8’s amino acid residues LEU23, VAL25, and GLU61. These hydrogen bonds trigger the creation of stable molecular linkages encircling IL-8’s active site. Beyond hydrogen bonds, IL-8 also generates hydrophobic partnerships with ZINC000022339916 via LEU23, ARG24, VAL25, and VAL 60. The involvement of hydrogen bonds and hydrophobic interactions contributes to fortifying the IL-8-ZINC complex and offers insights into ZINC000022339916’s potential as an effective inhibitory agent. As such, the synergic merging of these two interaction types could yield significant ramifications in advancing therapeutics or therapeutic strategies centered on the IL-8-ZINC complex, particularly in intervening within biological processes associated with IL-8.

Molecular dynamics simulations spanning 200 ns were employed to comprehend alterations in conformation, stability, and interactions within the free IL8-ZINC complex. Prior to commencing the molecular dynamics analysis, the IL8-ZINC complex’s average potential energy in its free state was determined. Table 3 enlists the values for RMSD, RMSF, SASA, Rg, and RDF. The process of ligand binding to target molecules can introduce modifications to structure and stability and even induce conformational shifts. RMSD values serve as indicators of structural deviations and stability (Ramadhan et al. 2022). For the IL-8-ZINC complex, the observed RMSD values fluctuated between 0.285327 nm and 0.430443 nm (Table 3). The RMSD plot unveiled the attachment of the ZINC small molecule within IL-8’s binding pocket, which led to limited deviations in structure and conformation when compared to the original IL-8 configuration (Fig. 4). However, minor fluctuations were noticeable due to the initial alignment of ZINC with the IL-8 binding pocket. Following these initial shifts, a state of equilibrium emerged, and the system remained stable throughout the simulation. Simulation outcomes conclusively demonstrated that the presence of ZINC compound molecules binding to IL-8 did not induce changes in the stability or structure of IL-8. The RMSD values were consistently low during several segments of the simulation, maintaining equilibrium and indicating the stability of the IL-8-ZINC complex (Fig. 4).

Figure 4. 

Dynamic changes in the structure of IL-8 upon ZINC binding. Plot illustrating the RMSD of the IL-8-ZINC complex over time. Representation of the RMSF of IL-8 post-binding. Evolution of the Rg for the IL-8-ZINC complex over time. Presentation of the SASA of the IL-8-ZINC complex across time. A visual depiction of the RDF for IL-8 after binding. These values have been derived from a 200 ns MD simulation timeframe. Each color signifies the value obtained for the free IL8-ZINC complex.

The variability in the local structure of the IL-8-ZINC complex was assessed through the employment of RMSF values, graphically presented in Fig. 4. A compilation of RMSF values can be found in Table 3, encapsulating the average fluctuation of residues. This comprehensive representation elucidates notable fluctuations within diverse regions of the IL-8 structure. While the simulations exhibited scarce binding of the ZINC small molecule, the binding of ZINC000022339916 to IL-8 induced sporadic, unpredictable residual fluctuations (Fig. 4). In the pursuit of a deeper comprehension of a protein’s tertiary structure volume and overall conformational profile, the calculation of Rg is pursued, serving as a gauge of protein stability within a biological context. A heightened Rg value signifies increased flexibility in protein arrangement (Muchtaridi et al. 2023). The calculated average Rg values for the IL8-ZINC complex were observed to range from 1.25559 nm to 1.28771 nm, respectively (Table 3). Subtle deviations in the Rg plot were noticeable within the IL8 structure upon interaction with the ZINC compound molecule, signifying cohesion within the complex ensemble. Initial Rg plots displayed considerable compactness, potentially attributed to the tightly packed IL-8 configuration, while subsequent simulations demonstrated a balanced Rg equilibrium for IL8 (Fig. 4). Importantly, in the Rg plot, the incorporation of the ZINC small molecule displayed minimal structural variance and no discernible conformational transformation within IL-8. Notably, throughout the simulations, IL-8 maintained structural integrity even after the binding of the ZINC compound on a molecular level.

Subsequently, the assessment of the protein’s interaction surface area with the surrounding solvent was conducted via 200 ns molecular dynamics simulations, focusing on SASA. SASA and a protein’s Rg exhibit a direct correlation. Subsequent observation of the mean SASA values yielded a range between 55.84815 nm² and 57.66862 nm² for the IL-8-ZINC complex, as depicted in Fig. 4. The higher SASA value within the IL-8-ZINC complex might stem from conformational changes that expose a marginal portion of IL-8 residues to the solvent environment. Moreover, the RDF was employed to quantify atom distribution within the simulated protein structure. Illustrated in Fig. 4, the RDF profile displayed minimal differentiation from the Rg chart. The RDF outcome suggests that the small molecule ZINC outperforms azelaic acid throughout the 200 ns MD simulation. RDF fluctuations could potentially signify alterations in protein conformation arising from interactions with ligands, environmental variations, or thermal processes. Furthermore, RDF can discern the nature of interactions occurring within proteins; a decrease in RDF implies hydrophobic interactions, while an increase indicates hydrophilic interactions (Hidayat and Fakih 2021).

Moreover, the kinetic energy of the solvent system, density, and volume, along with enthalpy and potential energy, are discerned and documented in Table 4. Over the course of the 200 ns MD simulation period, a notable dominance in kinetic energy, enthalpy, and potential energy was observed in relation to the compound molecule ZINC000012444335. This occurrence underscores the significant role these compounds play in steering system dynamics. This analytical endeavor was pursued with the intent of unveiling the pivotal contribution offered by the ZINC000012444335 compound to the interactions and transformative processes witnessed throughout the MD simulation experiment.

Table 5 displays the mean values of liberated binding energies along with their corresponding standard deviations. The outcomes of the MM-PBSA analysis suggested that all types of energy facilitate the interaction between IL-8 and the diminutive ZINC molecule. On the whole, the IL-8-ZINC complex displayed a more robust attraction compared to the IL-8 complex coupled with azelaic acid. Notably, the diminutive ZINC000022339916 molecule manifests a liberated binding energy of –165.451 kJ/mol. This discovery portrays that the compound’s molecular arrangement significantly contributes to the heightened bonding and more potent interactions between IL-8 and ZINC, positioning it as a promising contender for a potential inhibitory agent within the pertinent context.

The primary energies governing the interaction between IL-8 and ZINC compound molecules are van der Waals forces, electrostatic interactions, and the accessibility of solvent-exposed surfaces (surface area of solvent accessibility, or SASA). This phenomenon can be attributed to the attractive van der Waals forces operating between neighboring atoms within the IL-8-ZINC complex, as well as the electrostatic interactions that arise between charged groups present in both entities. Furthermore, the extent of surface area accessible to solvent molecules (SASA) assumes a vital role in facilitating contacts and interactions between these molecular constituents (Settanni et al. 2017). The amalgamation of these factors engenders an environment conducive to robust and enduring interactions between IL-8 and ZINC compound molecules within the intricate complex.

While the IL-8-ZINC interaction unfolded, an investigation into hydrogen bond formations within the complex was concurrently conducted. Across all complex systems, hydrogen bonds were identified with an occupancy rate exceeding 10%. Notably, ZINC000888090135 (at 75.48%) and ZINC000001893410 (at 88.17%) exhibited elevated hydrogen bond occupancy rates (Table 6). Remarkably, both demonstrated a distinctive characteristic of engaging in hydrogen bond formations with the same residues as those interacting with azelaic acid, specifically GLN6, LYS9, SER12, LYS13, ARG45, GLU46, and CYS48 located at the IL-8 binding site. This signifies a parallel interaction pattern concerning the development of hydrogen bonds between specific ZINC molecules and IL-8 within their binding pocket.

Protein-peptide docking was conducted to explore the binding of LL-37 to each IL-8-ZINC complex. The atomic contact energy (ACE) score was employed to evaluate the binding strength of LL-37. To enhance the credibility of the docking technique, a validation step was executed using PatchDock. This validation strategy encompassed a re-docking approach, wherein LL-37 was reconstructed to match its original crystal structure. Throughout this re-docking process, close attention was paid to RMSD values and active sites. Additionally, potential binding residues of LL-37 were identified. The outcomes of the method validation phase furnish valuable insights into the precision and efficacy of protein-peptide docking simulations, thereby establishing the dependability of subsequent analyses that assess the interplay between LL-37 and the IL-8-ZINC complex.

This investigation unveiled intriguing outcomes concerning the interplay between complexes of small molecules (ZINC000409414675, ZINC000888090135, ZINC000022339916, and ZINC000012444335) and IL-8 concerning LL-37 (Table 7). In this specific context, the favorable scores generated by the PatchDock analysis are captivating, characterized by binding energies spanning from 94.30 kJ/mol to 366.04 kJ/mol. These findings suggest that the trio of small molecules possesses the capability to effectively impede the binding of LL-37 to the surface of IL-8. Moreover, these results unlock the potential utility of these three compounds as potent inhibitors of IL-8.

To ensure the safety of the ZINC small molecule, a toxicity prediction process is conducted, as it can aid in evaluating the safety and efficacy of a drug even before its administration to humans or experimental animals. The evaluation of absorption, distribution, metabolism, and excretion (ADME) along with toxicity is conducted utilizing an online platform known as preADMET 2.0. For absorption, predictions are made for human intestinal absorption (HIA) and Caco-2 cell permeability. Distribution assessments encompass plasma protein binding (PPB) and blood-brain barrier (BBB) permeability, while toxicity profiles are examined for carcinogenic and mutagenic effects. The HIA value serves as an indicator of the extent to which the active substance is absorbed in the human intestine. Based on the percentage of HIA, a compound is classified as well absorbed (70–100%), moderately absorbed (20–70%), or poorly absorbed (0–20%) (Wadapurkar et al. 2018). As depicted in Table 8, the ZINC compound molecules, on the whole, align with the specified range criteria.

Moreover, the utilization of Caco-2 cell modeling was employed to anticipate the absorption of active compounds through the oral route in a controlled in vitro setting. The notable permeability of Caco-2 cells is reflected by values surpassing 0.9 x 10^–6 cm/s (Mohapatra et al. 2015). Solely the compound molecule ZINC000001893410 exhibits a Caco-2 cell value beneath 0.9 x 10^–6 cm/s, signifying limited permeability in terms of cellular membrane penetration. The plasma protein binding (PPB) value exerts an impact on the pharmacokinetic and pharmacodynamic attributes of drugs. A PPB value surpassing 90% suggests a strong binding affinity between the drug and plasma proteins, whereas a PPB value below 90% indicates a weaker binding propensity with plasma proteins, consequently facilitating proper distribution to the intended target sites. An additional parameter for distribution is the Blood-Brain Barrier (BBB) value, which signifies the concentration of a drug within the brain. This parameter holds significance in gauging a drug’s capability to traverse the blood-brain barrier. BBB values are stratified into two categories: values exceeding 0.3 are interpreted as facile blood-brain barrier penetration, while values below –1 indicate limited distribution into the brain (Vique-Sánchez 2021).