Methods

A total number of 25 iminoguanidine derivatives reported by PubChem with the accession number AID_1315713 were collected to perform this study. The MIC₅₀ values (minimum inhibitory concentrations) were converted to pMIC50 (-logMIC₅₀) values that were utilized in the QSAR investigation as dependent variables. The compound's pMIC₅₀ values are within the QSAR prerequisite range according to Shirvania and Fassihi [³¹]. The 3D structures of all the compounds including the reference compounds were prepared using MarvinView Europium.6 (http://www.chemaxon.com) and their energies were minimized with the Spartan’14 V1.1.4 PM6. The energy minimized molecules were subjected to optimization via DFT (density functional theory) method with B3LYP function [³²,³³] and 6-31G+ basic set [³⁴] with the GAUSSIAN 09 package [³⁵]. The optimized structures were imported into PaDEL-Descriptors v2.20 [³⁶], which calculated the molecular descriptors for each molecule. The calculated descriptors from the PaDEL software saved in an excel spreadsheet were pre-treated using NEW_V-WSP version 1.2 (DTC) laboratory software to sequester highly correlated descriptors from the distantly related ones, thus reducing redundancy in the calculated values. The pre-treated descriptors saved in the text file were transferred into the QSARINS software for further analysis [³⁷]. The chemical structure of the molecules and their pMIC50 values are available in Table 1. The QSARIN v2.2.4 software was adopted to perform a Genetic algorithm (GA) and model development. GA approach allows the proper selection of some exceptionally good descriptors to improve a model’s quality and predictive accuracy and by doing so, a multiple linear regression (MLR) is curated to obtain the additive relationship between biological endpoints and molecular descriptors [³⁸]. Thus, the technique is referred to as GA-MLR. Before developing the 2D-QSAR model, the data set were divided into 70 % training and 30 % test by applying the random selection implemented in QSARINS software.

Table 1 Chemical structures of iminoguanidine compounds as well as their activity levels. 

Compound No.
	R1	R2	R3	MIC50 (ug/Ml)	pMIC50
1	H	H	Cl	10.7	1.0294
2	H	H	N(CH3)2	32.5	1.5119
3	OH	H	H	58.5	1.7672
4	OCH3	H	H	38.4	1.5843
5	H	H	NO2	61.4	1.7882
6	F	H	H	28.6	1.4564
7	H	Naphc	Naphc	31.7	1.5011
8	H	Cl	H	19.2	1.2833
9	Cl	H	H	16.3	1.2122
10	H	H	CH(CH3)2	20.6	1.3139
11	H	H	OCH3	30.5	1.4843
12	Naphc	Naphc	H	30.4	1.4829
13	H	H	F	87.4	1.9415
14	H	H	Br	19.7	1.2945
15	H	OCH2C6H6	H	25.5	1.4065
16	Br	H	H	62.9	1.7987
17	OH	Br	H	145	2.1614
18	H	H	C6H6	34.6	1.5391
19	H	OH, OHa	H	338	2.5289
20	H	H	OH	443	2.6464
21	H	OH	H	391	2.5922
22	H	OCH3	H	48.6	1.6866
23	H	F	H	70	1.8451
24	H	Br	H	23	1.3617
25	H	Cl	Cl	24.1	1.3820
Reference compounds
Doxycycline (Vibramycin)
Chloroquine
Ruxolitnib
Ampicillin

2D-QSAR Model validation

Internal and external validations of the QSAR model are carried out by using the model to forecast training and test data. To assess the model's validity, we calculate several parameters. To begin, the coefficients of determination (R2) are computed to assess the strength of the relationship between actual and predicted values. The slope value (k) and (r2-r02/r2) are also used to assess the model's performance. The correlation coefficients between actual and expected values with and without intercept are represented by the r02 and r2 parameters, respectively. Ultimately, we compute the modified r2 (rm2) to show the model's predictive performance [³⁹,⁴⁰]. The following equations can be used to calculate the validation parameters:

Where, Y and Y^ denote the actual and predicted value of pMIC50, respectively. The threshold values for all the parameters used in these studies were detailed according to Edache et al [⁴¹].

Protein generation and preparation

Pseudomonas aeruginosa and the 2019 novel coronavirus protein sequence were downloaded from the NCBI database (https://www.ncbi.nlm.nih.gov/) with accession number KJJ19201 and PDB: 7E23, respectively. A BLAST exploration was carried out to find appropriate templates. The UPGMA algorithm [⁴²] was used to infer the evolutionary history. The evolutionary distances were calculated using the Poisson correction method [⁴³] and are measured in amino acid substitutions per site. Each descending clade's proportion of sites with at least one unambiguous base in at least one sequence is shown next to each internal node in the tree. This study looked at five and six different amino acid sequences of P. aeruginosa and SARS-CoV-2 spike proteins, respectively. For each sequence pair, all unclear locations were deleted (pairwise deletion option). In total, there were 215 and 239 positions in the final dataset, respectively. MEGA X was used to undertake evolutionary analyses [⁴⁴]. Pair-wise alignment of the template and query sequences were constructed using MODELLER v10.1 to check the results of the phylogenetic tree analysis [⁴⁵-⁴⁷]. For query and template sequences, BioEdit [⁴⁸] was used to build a sequence alignment file in PIR format. Python Scripts were used to construct the Multiple Template Modeling (MTM) protocol from module 5 structures. With the help of Swiss-Pdb Viewer [⁴⁹], for energy minimization, the model with the most reduced DOPE score (Discrete Optimized Protein Energy) was chosen. For P. aeruginosa and SARS-CoV-2 spike, the template and the predicted modelled structures were superimposed with RMSD values of 0.176 and 0.261, respectively (Fig. 1).

Fig. 1 The modelled protein structures of (A) P. aeruginosa and (B) SARS-CoV-2 spike protein from Modeller v10.1 using multiple template modeling. 

Equipment and programming

Molegro Virtual Docker 6.0 2013 [⁵⁰] and Autodock Vina with PyRx v0.8 [⁵¹] were used for docking simulations with core i3 windows processors. The molecular docking simulations of all compounds were performed with the predicted modelled crystal structure of P. aeruginosa and SARS-CoV-2. Molegro Virtual Docker (MVD) which was considered as one of the best and accurate docking packages were used to study the interactions energies between the stationary modelled proteins and the ligands. The design of the demonstrated proteins was created by pressing icons, and the amino acid flaws that presented as a warning were repaired. The protein surface and binding cavity have both been identified. After resetting the icon, the design of the demonstrated proteins is ready for docking simulation. The docking wizard was used to perform interactions of the protein-ligand and the parameters for the docking interactions were set as ten (10) runs, 50 population size, 1500 iterations, and lastly, the maximum number of poses is 4. Next was the PyRx [⁵¹] docking simulations setup. The modelled proteins were prepared using PyRx v0.8 software and the docking experiments were subsequently accomplished using the AutoDock Vina package [⁵²]. The network place for docking of the ligands against the displayed protein of P. aeruginosa was set as exhaustiveness = 8, center_x = 16.8669, center_y = 53.3921, center_z = -4.3156, size_x = 55.3409311867, size_y = 155.016245728, size_z = 141.969700775, and against the modelled SARS-CoV-2 spike protein was set as center (x = -1.5643, y = 12.818, and, z = 12.3015), size (x = 49.2203297997, y = 44.6869613934, and, z = 51.2521103764), respectively. The Discovery studio 2020 client was used to visualize the outcomes of the docking simulation.

Molecular Dynamics Simulations

Utilizing the CHARMM-GUI web server [⁵³], a cubic box of TIP3 water particles with an edge distance of 10 Å was utilized to submerge the docked complex. A 0.15M KCl counterion was utilized to neutralize the system. To deliver the inner strain energies of the whole system, the whole framework or system was exposed to 5,000 run steps of energy minimization utilizing the steepest drop integrator. The NVT equilibration was then performed at 303.15 K for 0.01 ns. 500,000 runs i.e., 1 ns of unrestrained isothermal isobar group (NPT) were performed for the production and were completed at 303.15 K, utilizing 1 atm (1.01325 bar) with 2 fs time steps. All computations were performed using the CHARMM36m force field [⁵⁴], NAMD [⁵⁵], and VMD [⁵⁶] for representation.

Phylogenetic tree and Homology modeling

The sequence of P. aeruginosa and SARS-CoV-2 spike glycoproteins were obtained from the NCBI (National Center for Biotechnology Information) database (http://www.ncbi.nlm.nih.gov/) with accession number KJJ19201.1 and PDB 7E23, respectively, were used to BLASTp search against the protein data bank to select the template structures. The query sequence and the BLAST template structures were first analyzed through the MEGA-X [⁶³] software to investigate the phylogenetic relationship based on the test UPGMA (Unweighted pair-group method with arithmetic) method as shown in Fig. 2. The phylogenetic tree which is obtained by MUSCLE [⁶⁴] shows the evolutionary relationship between the proteins sequences as it was a rooted tree. Phylogenetic analysis demonstrated clear sequestration of P. aeruginosa and SARS-CoV-2 in two different biological groups and these biological groups were also divided into subgroups. All the biological groups are grouped according to their protein sequence similarity. The P. aeruginosa query sequence (accession number KJJ19201.1) shows high similarity to PDB 4PS6 and 1UUZ, respectively as shown in Fig. 2(A). The phylogenetic sequence analysis of the spike protein of SARS-CoV-2 (PDB: 7E23) and the other five templates shows a close similarity (Fig. 2(B)). Concerning overall sequence similarity PDB 7MZF spike protein S1 has the highest similarity to PDB 7E23. The chain A inhibitors of vertebrate lysozyme (PDB 4PS6 and PDB 1UUZ) and spike protein S1 (PDB 7MZF) have very close relationships to the query sequence, respectively.

Fig. 2 Phylogenetic analysis of (A) P. aeruginosa and (B) SARS-CoV-2 spike protein homologs were constructed by the UPGMA method using the MEGA X program.  

Further analysis using MODELLER 10.1 was performed on the query sequence and the template sequences structures to confirm the phylogenetic tree analysis and a better insight into their protein structure similarity. According to the result malign3D command implemented in MODELLER 10.1, the tree diagram calculates a weight pair group average clustering on a distance matrix (Fig. 3). Out of the template structure, the best homologous structure belonged to an inhibitor of vertebrate lysozyme (PDB: 4PS6) (Fig. 3(A)) and spike protein S1 (PDB: 7mzf) (Fig. 3(B)) because of its similarity and high resolution with the query structure was selected for homology modeling. The results obtained thereof confound to the phylogenetic analysis. The selected PDB was used as a template structure for the construction of the P. aeruginosa and SARS-CoV-2 spike protein, respectively. Pseudomonas aeruginosa, with the accession number of KJJ19201.1, has 199 amino acids in its protein sequence, while SARS-CoV-2 spike protein S1with PDB code 7E23 has 196. The Pseudomonas aeruginosa (PDB 4PS6) and SARS-CoV-2 spike protein S1 (PDB code 7MZF) were selected as the template based on the phylogenetic tree analysis and dendrogram command using the align2d script implemented in MODELLER. An alignment of Pseudomonas aeruginosa (accession number of KJJ19201.1) and SARS-CoV-2 spike protein S1with PDB code 7E23 with the template sequences are presented in Fig. 4. The MODELLER constructed five (5) different structure modelled proteins. The receptor generated by MODELLER was ranked and scored using discrete optimized protein energy (DOPE) score (Data presented in Table 3 and 4). The top one out of 5 models with the lowest DOPE scores were selected and assessed using ERRAT, VERIFY3D, and Ramachandran plots.

Fig. 3 The weighted pair-group matrix of (A) P. aeruginosa and (B) SARS-CoV-2 spike protein homologs were constructed by the modeller program. 

Fig. 4 Pairwise sequence alignment between the template and the query (A) and putative inhibitor of vertebrate lysozyme [Pseudomonas aeruginosa], (B) chain A, Spike protein S1 [severe acute respiratory syndrome coronavirus 2]. 

We perform energy minimization to the assessed receptors to remove side chains and later use it in the molecular docking and MD simulations.

Model number two (PA-model02.pdb) and model number one (CoV-model01.pdb) having lower DOPE scores as presented in Table 4 and Table 5, respectively were considered as the best model of P. aeruginosa and SARS-CoV-2.

The crystal structures of Pseudomonas aeruginosa and SARS-CoV-2 were predicted, polished, and validated. The final predicted structure of P. aeruginosa protein has a modeller objective (molpdf) score of 886.59521, genetic algorithm 342 (GA341) score of 1.0000, and the discrete optimized protein energy (DOPE) score of -15917.49902 was picked out. The modelled crystal structures SARS-CoV-2 spike protein S1 has a molpdf score of 1122.96790, GA341 score of 1.0000, and a DOPE score of -22118.33008 was also picked out. Following the sequence alignment, the DOPE score profile was generated using MODELLER's "assess dope" function to give an insight into the quality of input alignment (Fig. 5). Fig. 5 shows that the DOPE score profile of the modelled (query) and template proteins meet over a wide range of residues, and the results suggest no problematic locations.

Fig. 5 The DOPE score visibility of (A) P. aeruginosa and (B) SARS-CoV-2 spike protein incurred from Modeller software.  

Validity of the structures analyzed for the P. aeruginosa receptor by the Ramachandran plot shows 94 % residues residing in most favored regions, 4.8 % in the additional allowed region, and 0.6 % in the generously allowed region. Only one residue was settled in the disallowed region, which constituted 0.6 % of the total protein, which constituted a good quality of the predicted model (Fig. 6). Other tests, such as Verify 3D and ERRAT also known as “overall quality factor”, were used to ensure the quality of the revised models. For non-bonded atomic interactions, the overall quality factor with larger scores indicates higher quality [⁶⁵]. In general, a range of greater than 50 (> 50) indicates a high-quality model [⁶⁶]. In this research, the ERRAT score for predicted modelled P. aeruginosa protein has an overall quality factor of 83.582 % (Fig. 7). Therefore, this confirmed that modelled putative inhibitor of vertebrate lysozyme pseudomonas aeruginosa protein has a high resolution and quality protein models. Verify-3D score profile access the model's quality. The Verify 3D profile of the modelled protein is shown in Fig. 8, and residues with an averaged 3D-1D score larger than zero are regarded as reliable. Some of the residues in the modelled protein have a computability score that is greater than zero. It was discovered that the structurally and functionally significant residues in the modelled potential inhibitor of vertebrate lysozyme pseudomonas aeruginosa protein had a score ranging from 0.30-0.60, indicating that the model's quality is comparable to that of high-resolution crystal structures.

Fig. 6 Ramachandran plot calculation for 3D model of P. aeruginosa. 

Fig. 7 ERRAT gave a quality factor of 83.582 % to the modelled P. aeruginosa protein structure. 

Fig. 8 Verify 3D score diagram validating the modelled P. aeruginosa protein. 

Ramachandran plot analysis was carried out on the modelled SARS-CoV-2 spike protein showed 91.1 % residues residing in highly favored regions, 8.9 % in additional allowed regions, 0.0 % residues in generously allowed and disallowed regions, respectively (Fig. 9). This also validates that the modelled crystal structure is a good quality model. The quantity of the structural error at each amino acid residue in the 3D structural model was given by the ERRAT plot. The overall quality factor of the computed model was 79.032 % (Fig. 10). Also, the Verify3D plot of the SARS-CoV-2 modelled spike protein (Fig. 11) was obtained for the structure validation and it showed as PASS. Verify-3D ratings near 1.0 correspond to scores that would be expected for a legitimate protein of similar size. It was also discovered that more than 90 % of the residues had a score of more than 0.2, implying that more than 90 % of the residues were complemented by the 1D-3D model, implying that the predicted model's quality is adequate.

Fig. 9 Ramachandran plot computation for a 3D model of SARS-CoV-2 spike protein S1. 

Fig. 10 ERRAT gave a quality factor of 79.032 % to the modelled SARS-CoV-2 spike protein structure. 

Fig. 11 Verify 3D score diagram validating the modelled SARS-CoV-2 spike protein. 

Rating of docking simulation methodological analysis

At first, the Molegro virtual docker (MVD) and AutoDock Vina programs estimated with PyRx were selected to carry out a docking-based virtual screening approach. Both Scoring Functions were supplied to select the best ligand against P. aeruginosa and SARS-CoV-2.

Molecular docking simulation is of two types; namely the oriented docking simulation and the blind docking simulation [⁶⁷,⁶⁸]. Blind molecular docking simulation is a docking procedure in which the position of the enzyme's active site is unknown, and the grid center of attention used during grid box determination is not unique to a specific area but is applied throughout the entire enzyme area. The grid box was specifically resolved after orienting docking by establishing the position of the enzyme active site. Because the active site region has been established, oriented docking was used in this study (using Melogro Virtual Docker). In this case, an oriental docking simulation was carried out. The molecular docking simulation of 25 iminoguanidine derivatives, against modelled proteins Pseudomonas aeruginosa and SARS-CoV-2 Spike protein, were performed using MVD software. The modelled protein was taken as the substrate for molecular docking. The substrate had only three cavities. From the three (3) predicted cavities the one with the highest volume was picked for molecular docking simulations. The cavity with volumes 48.64 ų (Fig. 12) and 11.78 ų (Fig. 13) are used for docking against P. aeruginosa and SARS-CoV-2, respectively. The MolDock score [⁶⁹] was picked out for rating the inhibitor poses, and for all the modelled proteins docking executed here, the poses picked out as the best were those with the most prominent MolDock score. The MolDock score of all the ligands including the reference drugs ranges from -125.553 kcal/mol to 7320.88 kcal/mol as presented in Table 6. The reference drug (Ruxolitinib) has the highest negative MolDock score of -125.553 kcal/mol and -117.237 kcal/mol, followed by compound 15 with -117.489 and -102.85 kcal/mol against P. aeruginosa and SARS-CoV-2 Spike protein, respectively.

Fig. 12 Modelled protein of potential inhibitor of vertebrate lysozyme [Pseudomonas aeruginosa] with electrostatic surface and active site cavity (green grid). 

Fig. 13 The electrostatic surface and active site cavity (green grid) of modelled protein of SARS-CoV-2 Spike protein S1. 

The reference drug (Ruxolitinib) against the protein of putative inhibitor of vertebrate lysozyme (Pseudomonas aeruginosa) formed two hydrogen bond interactions with Cys132 and Asp193 (Fig. 14(A)). Next is compound 15 which also formed two hydrogen bond interactions with the modelled protein first with Asn125, secondly with Gln189 (Fig. 14(B)). Against the modelled protein of SARS-CoV-2 Spike protein S1, the reference drug Ruxolitinib with the highest negative MolDock score of -117.237 kcal/mol formed hydrogen bond interactions with Lys126, Asn149, and Cys148, respectively (Fig. 15(A)). In Fig. 15(B), compound 15 formed five (5) hydrogen bond interactions with the same modelled proteins. The five hydrogen bond interactions of compound 15 are with Arg125, Lys126, and Ser127, respectively.

Fig. 14 Hydrogen bond interactions between (A) Ruxolitinib and (B) compound 15 with modelled protein of putative inhibitor of vertebrate lysozyme [Pseudomonas aeruginosa]. 

Fig. 15 Hydrogen bond interactions between (A) Ruxolitinib and (B) compound 15 with modelled protein of SARS-CoV-2 Spike protein S1. 

Again, blind molecular docking simulations were carried out for all 25 compounds including the reference compounds with PyRx v0.8 software. In the molecular docking simulations with PyRx, we make use of the modelled proteins as the target. The binding affinity values obtained as a result of the molecular docking analyses are shown in Table 7. The negative value of binding affinity as presented in Table 6 suggests the stability of the protein-ligand complex. From the table, the reference compound “Vibramycin” has the highest binding affinity value of -7.4 kcal/mol for P. aeruginosa and SARS-CoV-2, respectively. The reference compound (Vibramycin) interacted with P. aeruginosa modelled protein structure (Fig. 16) showed the corresponding amino acids such as aspartic acid, phenylalanine, glycine, glutamine, tyrosine, and so on. The docking of Vibramycin forms some hydrogen bonds and electrostatic interactions with the modelled protein. Vibramycin is observed to show three (3) hydrogen bond interactions with Ser192 (3.87 Å), Gln189 (5.33 Å), and Gln189 (4.27 Å) against P. aeruginosa modelled protein (Fig. 16(A)). Vibramycin exhibits three electrostatic interactions, one unfavorable acceptor-acceptor interaction, one carbon-hydrogen bond with Pro131 (6.23 Å), and two conventional hydrogen bond interactions with Ser182 (3.87 Å) and Glu184 (4.79 Å), respectively against SARS-CoV-2 protein (Fig. 16(B)).

Fig. 16 The 2D illustrations of interactions between Vibramycin and the modelled (A) P. aeruginosa and (B) SARS-CoV-2 spike protein receptors. 

Compounds 7 and 15 were selected based on their remarkable binding affinity of -6.9 and -7.2 kcal/mol for P. aeruginosa and SARS-CoV-2 spike protein receptors respectively compared to the reference compound’s -7.4 kcal/ mol. The molecular docking contact modes for compounds 7 and 15 are shown in Figures 17 and 18, respectively. The molecular docking interaction of compound 7 with the active site of modelled P. aeruginosa protein reveals two conventional hydrogen bonds between the hydrogen atoms of compound 7 and the amino acid Ser109 with distances of 3.45 and 4.04, respectively, and one carbon-hydrogen bond between the carbon atom and the amino acid Pro107 with a distance of 5.47. The hydrophobic interaction was performed between one of the naphthalene ring and amino acids Arg136 (5.84 Å) and Cys127 (5.04 Å). Moreover, the electrostatic interaction is observed between one of the naphthalene rings and the amino acids Asp193 (5.94 Å) and Arg136 (5.48 Å).

Fig. 17 The 3D and 2D representations of interactions between compound 7 and modelled P. aeruginosa receptor. 

Fig. 18 The 3D and 2D representations of interactions between compound 15 and modelled SARS-CoV-2 spike protein receptor. 

Nevertheless, compound 15 forms two conventional hydrogen bonds with amino acid Asn11 (3.60 and 3.95 Å), one pi-donor hydrogen bond with Trp104 (3.95 Å) with the active site of modelled SARS-CoV-2 spike protein. The hydrophobic interaction was performed between one of the benzyloxy and toluene rings with amino acids Leu36 (5.13 Å) and Val35 (4.78 Å), respectively. Furthermore, the electrostatic interaction is observed between hydrogen and the amino acid Trp104 (4.13 Å). The hydrogen bonds, the hydrophobic interactions, and the electrostatic interactions in the complexes showed that compounds 7 and 15 of iminoguanidine derivatives are most active against modelled P. aeruginosa receptor and modelled SARS-CoV-2 spike protein, respectively.

Lipinski's rule of five aids [⁷⁰] in the differentiation of drug-like and non-drug-like compounds. It forecasts a high likelihood of success or failure due to drug-likeness for compounds that meet two or more of the rules stated in Table 7. The chemical structure of compounds 7 and 15 was submitted to the MolSoft server in SMILE format, to guess the in silico pharmacokinetic parameters. Based on these properties, we investigated how our lead ligands met the Lipinski rule of five (5). The table (Table 8) shows the properties of compounds 7 and 15 in which both of the compounds fulfilled the Lipinski rule of 5. The MollogS value is less than or equal to -4.0. These compounds have a low MolLogS value, indicating that they are difficult to dissolve in water [⁷¹]. The aromatic and non-polar fragments around the chemical structures are most likely to blame. Despite this, their intestinal absorption is high, exceeding 90 %.

The leads in Fig. 19(A) and 19(B) showed that these molecules display ADME profile within the desired range as presented in Table 7, and therefore the values of these parameters shed light on their drug-likeness features, which could set up an essential starting point to encourage experimental procedures.

Fig. 19 The drug-likeness profile of (A) compound 7 (B) compound 15. 

Molecular dynamics (MD) simulations

The study of molecular dynamics simulation in the modelled SARS-CoV-2 spike protein was being used to figure out how flexible is docked compound 15, and this evaluation was lasted up to 1 ns. To compute these RMSD, RMSF, and SASA values we used the VMD plugin with a probe radius set to 1.4 Å. Other thermodynamic parameters investigated are the total energy, potential energy, and kinetic energy of the docked protein-ligand complex [⁷¹]. To understand the structural aberrations across the simulation track, the RMSD of the c-alpha atoms in the protein-ligand complex was analyzed. The RMSD value of the complex increased at the start of the simulation, as seen in Fig. 20(A). The protein-ligand complex then stabilized at 400 and 500 ps, before rapidly increasing from 600 to 800 ps. The protein-ligand complex remained stable until the simulation's final phase, 1 ns. To estimate the amino acid residue flexibility, the RMSF of the protein-ligand complex was calculated. The majority of residues had RMSF values < 2.0 Å, as seen in Fig. 20(B). The presence of peaks indicates the areas fluctuate the most during the simulation. The fluctuations observed are around 1.9 Å (600 ps) and less around 0.4 Å (1 ns) which signifies the stability of the system. This RMSF profile indicates a complex that is less flexible and inflexible. The surface area of a group that is accessible to a solvent probe is known as the solvent-accessible surface area (SASA) [or temporarily accessible surface area (ASA)]. Protein folding and stability have always been determined by the solvent-accessible surface area (SASA). The SASA trend in the predicted protein-ligand complex was higher, indicating an increase in surface volume. Nonetheless, this complex had a high number of SASA deviations, indicating that the protein's surface area was changing significantly. The SASA value for the un-simulated protein-ligand complex has a lower value in comparison with the simulated protein-ligand (Table 9) which indicates that the simulated protein-ligand complex is more exposed to water and also increases in cavities. The SASA profile in Fig. 20(C) confirmed the claim.

Fig. 20 The molecular dynamics simulation of the complexes. (A) Root mean square deviation (RMSD), (B) Root mean square fluctuations of the amino acid residues (RMSF), and (C) Solvent accessible surface area (SASA). 

The MD simulations of the 10 frames and the hydrogen bond information of these 10 confirmations are listed in Fig. S6. From the figures, conformation 3, 5, 8, and 9 has 4, 2, 2, and 2 hydrogen bonds, respectively. The amino acids that appear more as conventional hydrogen bond and carbon-hydrogen bond interactions are Ser41 and Asp 32, respectively. From conformation 1 to 10, amino acid residue Ser41 appears in all the interactions either as hydrogen bond or van der Waals interactions to the ligand. After 1 ns MD simulation (the last confirmation “10”), the residue interactions and protein-ligand 2D contact contour map show that van der Waals and hydrophobic interactions were crucial contributions for the stability of the hit drug.

Fig. S7 depicts a change in total energy versus time in 1 ns of molecular dynamics simulation, which displays a slight fluctuation throughout the 1 ns. It also suggested that when the protein was simulated up to or beyond 1 ns, total energy (at or above -137713. 02 KJ/mol), and kinetic energy (at or beyond 31498.99 kJ/mol), compound 15 produced a stable combination with expected protein structure. Moreover, the examined system's potential energy in 1 ns of simulation at 303.15 K was -169212.01 KJ/mol (Fig. S7(C)). As a result, under stable temperature settings, the derived equilibrium structure of the expected spike protein S1 structure was obtained at 303.15 K. These findings reveal that in the MD simulation, energy conservation was met. At the end of the MD simulation, the coordinates of the ligand in the inserted binding site were changed (Fig. S6), demonstrating the utility of MD simulations after ligand docking in the binding site. Except for Asp32, Leu36, Phe42, Phe10, Phe6, Asn11, and Val35, the active site residues determined by docking were changed, and new residues such as Cys29, Val192, Tyr33, Val30, Ala31, Ile26, Leu181, and Val63 are positioned near ligands and could involve in the interaction, according to exploratory molecular dynamics simulations on the complex between the ligand.