Nontoxic and Naturally Occurring Active Compounds as Potential Inhibitors of Biological Targets in Liriomyza trifolii

In recent years, novel strategies to control insects have been based on protease inhibitors (PIs). In this regard, molecular docking and molecular dynamics simulations have been extensively used to investigate insect gut proteases and the interactions of PIs for the development of resistance against insects. We, herein, report an in silico study of (disodium 5′-inosinate and petunidin 3-glucoside), (calcium 5′-guanylate and chlorogenic acid), chlorogenic acid alone, (kaempferol-3,7-di-O-glucoside with hyperoside and delphinidin 3-glucoside), and (myricetin 3′-glucoside and hyperoside) as potential inhibitors of acetylcholinesterase receptors, actin, α-tubulin, arginine kinase, and histone receptor III subtypes, respectively. The study demonstrated that the inhibitors are capable of forming stable complexes with the corresponding proteins while also showing great potential for inhibitory activity in the proposed protein-inhibitor combinations.


Introduction
A concerning problem that threatens food security around the world is the emergence of insects capable of developing resistance to insecticides [1]. Excessive use of many of these insecticides is associated with various health and environmental issues [2][3][4]. Liriomyza trifolii is a highly polyphagous pest in crop fields and greenhouses that has detrimental economic impacts [5]. Both larvae and adults selectively eat only the layers with the least amount of plant cellulose [6]. Stippling is one example of the damage in crop plants caused by the sap-sucking female fly; internal mining caused by larvae is another such example. These various types of damage allow pathogenic fungi to enter

Binding Site Prediction and Protein-Ligand Interaction
The putative ligand binding sites (both major and minor) for the predicted proteins were identified through Discovery studio software and were visualized ( Figure 1). All target proteins (acetylcholinesterase, α-tubulin, actin, arginine kinase, histone subunit III, Hsp90, elongation factor 1-alpha, and carbomoylphosphate synthase) were docked with the ligands, most of which were phytochemicals derived from the leaves of Phaseolus vulgaris [16,17] and the yeast extract. We evaluated the protein-ligand interaction through SAMSON software [18]. It was found that the tool has discrepancies in results for accurate pose prediction among the various putative docking poses.

Molecular Docking and Binding Free Energy Calculation
The prepared protein structures of (acetylcholinesterase, α-tubulin, actin, arginine kinase, histone subunit III, Hsp90, elongation factor 1-alpha, and carbamoyl phosphate synthase) were docked using SAMSON software with phytochemical compounds and yeast extracted compounds listed in the supplementary data. The results of the docking studies were provided in Table 2, and it was revealed that the phytochemical compounds were superior to the yeast extract compounds based on the docking score. All docking results were monitored by scoring functions that predict how well the ligand binds in a particular docked pose. This scoring function gives the ranking of the ligands. In the present study, the docking score was taken into consideration for the selection of the best ligands. This allowed us to explain the mechanisms of insect death. A mathematical empirical scoring function was used to approximately predict the binding affinity between two molecules after they have been docked by approximating the ligand's binding free energy [20]. It includes various force field interactions such as electrostatic and van der Waals contributions, which influence ligand binding. Subsequently, the docked structures were queried for binding free energy calculation. The results of binding free energy calculation were provided in interaction through SAMSON software [18]. It was found that the tool has discrepanci in results for accurate pose prediction among the various putative docking poses.  [19]. The large red sphere represents th cavities surrounding the active sites and was visualized using the visualization module Discovery Studio 3.0 visualization.

Molecular Docking and Binding Free Energy Calculation
The prepared protein structures of (acetylcholinesterase, α-tubulin, actin, arginin kinase, histone subunit III, Hsp90, elongation factor 1-alpha, and carbamoyl phospha synthase) were docked using SAMSON software with phytochemical compounds an yeast extracted compounds listed in the supplementary data. The results of the dockin studies were provided in Table 2, and it was revealed that the phytochemic compounds were superior to the yeast extract compounds based on the docking scor All docking results were monitored by scoring functions that predict how well th ligand binds in a particular docked pose. This scoring function gives the ranking of th ligands. In the present study, the docking score was taken into consideration for th selection of the best ligands. This allowed us to explain the mechanisms of insect deat A mathematical empirical scoring function was used to approximately predict th binding affinity between two molecules after they have been docked by approximatin the ligand's binding free energy [20]. It includes various force field interactions such a electrostatic and van der Waals contributions, which influence ligand bindin Subsequently, the docked structures were queried for binding free energy calculatio The results of binding free energy calculation were provided in Table 2. It was foun that binding energy values supported the docking result well. Hesperidin, Naringin, an Rosmarinic acid have higher binding energies than other compounds. All of the oth values contribute to the ΔG values which reflect the binding energy of the protein-ligan complex.  [19]. The large red sphere represents the cavities surrounding the active sites and was visualized using the visualization module of Discovery Studio 3.0 visualization.

Binding Pose Analysis
The binding mode of the compounds with proteins (acetylcholinesterase, α-tubulin, actin, arginine kinase, histone subunit III, Hsp90, elongation factor 1-alpha, and carbamoyl phosphate synthase) showed the different interactions between the proteins and ligands showed in Table 3. The interactions between the inhibitors and their target proteins, as well as their binding modes and orientations, are shown in Figures 2-9.

Root Mean Square Deviation (RMSD) Analysis
Calculations of the RMSD for the ligand-enzymes complex were used to determine the dynamic stability and conformational perturbation, which occur in each of the simulated systems during the simulation time scale. The RMSD values were calculated for the following protein-inhibitors combinations: acetylcholinesterase with disodium 5′inosinate and petunidin 3-glucoside; actin with calcium 5′-guanylate D and chlorogenic acid; α-tubulin with chlorogenic acid alone; arginine kinase with kaempferol-3,7-di-Oglucoside I, hyperoside D, delphinidin 3-glucoside ID, and histone subunit III complexes with myricetin 3′-glucoside ID and hyperoside D. All the trajectories reached equilibrium state after 20 ns, as shown in Figure 10. The RMSD values for all complexes are observed to be stable during the 50 ns simulation.

Root Mean Square Deviation (RMSD) Analysis
Calculations of the RMSD for the ligand-enzymes complex were used to determine the dynamic stability and conformational perturbation, which occur in each of the simulated systems during the simulation time scale. The RMSD values were calculated for the following protein-inhibitors combinations: acetylcholinesterase with disodium 5 -inosinate and petunidin 3-glucoside; actin with calcium 5 -guanylate D and chlorogenic acid; αtubulin with chlorogenic acid alone; arginine kinase with kaempferol-3,7-di-O-glucoside I, hyperoside D, delphinidin 3-glucoside ID, and histone subunit III complexes with myricetin 3 -glucoside ID and hyperoside D. All the trajectories reached equilibrium state after 20 ns, as shown in Figure 10. The RMSD values for all complexes are observed to be stable during the 50 ns simulation.

Radius of Gyration (Rg) Analysis
The Rg factor is best described for the stability of receptor-ligand complexes during the molecular dynamics simulations. The results demonstrate that the Rg values during different time points for the acetylcholine esterase, actin, α-tubulin, arginine kinase, and histone subunit III complexes to their respective ligands are constant during 50 ns simulation, which indicates the compactness of all of the proteins (Figure 11).

Radius of Gyration (Rg) Analysis
The Rg factor is best described for the stability of receptor-ligand complexes during the molecular dynamics simulations. The results demonstrate that the Rg values during different time points for the acetylcholine esterase, actin, α-tubulin, arginine kinase, and histone subunit III complexes to their respective ligands are constant during 50 ns simulation, which indicates the compactness of all of the proteins (Figure 11).

Hydrogen Bond Analysis
The number of hydrogen bonds for the ligand-enzymes complexes are plotted over a 50-ns MD simulation interval ( Figure 12). Since hydrogen bonds constitute a transient connection that provides stability to the receptor-ligand complex, they constitute an important factor to consider when discussing receptor-ligand stability. These bonds determine the specificity of the binding mode. In this study, we have calculated all of the hydrogen bonds for all of the complexes. The numbers of hydrogen bonds at different time points have been calculated, as shown in Figure 12. The average number of hydrogen bonds calculated for inhibitors (disodium 5′-inosinate and petunidin 3glucoside), (calcium 5′-guanylate D and chlorogenic acid), chlorogenic acid, (kaempferol-3,7-di-O-glucoside I, hyperoside D, delphinidin 3-glucoside ID) and (myricetin 3′-glucoside ID, hyperoside D) are (0-6, 0-5), (0-7, 0-9), 0-9, (0-8, 0-9, 0-10), respectively. All of the predicated ligands have shown continuous hydrogen bonding

Hydrogen Bond Analysis
The number of hydrogen bonds for the ligand-enzymes complexes are plotted over a 50-ns MD simulation interval ( Figure 12). Since hydrogen bonds constitute a transient connection that provides stability to the receptor-ligand complex, they constitute an important factor to consider when discussing receptor-ligand stability. These bonds determine the specificity of the binding mode. In this study, we have calculated all of the hydrogen bonds for all of the complexes. The numbers of hydrogen bonds at different time points have been calculated, as shown in Figure 12. The average number of hydrogen bonds calculated for inhibitors (disodium 5 -inosinate and petunidin 3-glucoside), (calcium 5 -guanylate D and chlorogenic acid), chlorogenic acid, (kaempferol-3,7-di-O-glucoside I, hyperoside D, delphinidin 3-glucoside ID) and (myricetin 3 -glucoside ID, hyperoside D) are (0-6, 0-5), (0-7, 0-9), 0-9, (0-8, 0-9, 0-10), respectively. All of the predicated ligands have shown continuous hydrogen bonding during the 50 ns simulation, which demonstrates the stabil-ity of the complexes. The only exception was chlorogenic acid, which only shows stable hydrogen bonding in the span of 35 ns.
during the 50 ns simulation, which demonstrates the stability of the complexes. The only exception was chlorogenic acid, which only shows stable hydrogen bonding in the span of 35 ns.

Root Mean Square Fluctuation Analysis (RMSF)
The RMSF value refers to the flexibility and mobility of structure-a higher value of RMSF indicates a loosely bonded structure with twists, curves, and coils, while a lower value of RMSF indicates a stable secondary structure, including α-helix and beta-sheets. Our RMSF analysis demonstrates that all of the ligands showed less conformational variations during binding and can act as stable complexes ( Figure 13).

Root Mean Square Fluctuation Analysis (RMSF)
The RMSF value refers to the flexibility and mobility of structure-a higher value of RMSF indicates a loosely bonded structure with twists, curves, and coils, while a lower value of RMSF indicates a stable secondary structure, including α-helix and beta-sheets. Our RMSF analysis demonstrates that all of the ligands showed less conformational variations during binding and can act as stable complexes ( Figure 13).

Molecular Mechanics Poisson-Boltzmann Surface Area Free Energy Calculations
The binding capacity of the ligand towards the receptor is quantitatively estimated by binding free energy analysis. Binding free energy is the summation of all non-bonded interaction energies. The binding free energy of the interactions between acetylcholine esterase, actin, α-tubulin, arginine kinase, and histone subunit III and the docked ligands has been estimated using the molecular mechanics Poisson-Boltzmann surface area tool (G_MMPBSA) [21]. This useful tool allows for efficient and reliable free energy simulation to model protein-ligand interactions. Our binding energy analysis spanning 50 ns MD simulation trajectories show that all ligands have a binding affinity towards enzyme inhibition and form stable complexes. Other different kinds of interaction energies, such as van der Waals energy, electrostatic energy, polar solvation energy, and solvent accessible surface area (SASA) energy, have been also calculated for all the Tools Shapes complexes (Figure 14). Results indicate that van der Waals, electrostatic, and SASA energy negatively contribute to the total interaction energy, while only polar

Molecular Mechanics Poisson-Boltzmann Surface Area Free Energy Calculations
The binding capacity of the ligand towards the receptor is quantitatively estimated by binding free energy analysis. Binding free energy is the summation of all non-bonded interaction energies. The binding free energy of the interactions between acetylcholine esterase, actin, α-tubulin, arginine kinase, and histone subunit III and the docked ligands has been estimated using the molecular mechanics Poisson-Boltzmann surface area tool (G_MMPBSA) [21]. This useful tool allows for efficient and reliable free energy simulation to model protein-ligand interactions. Our binding energy analysis spanning 50 ns MD simulation trajectories show that all ligands have a binding affinity towards enzyme inhibition and form stable complexes. Other different kinds of interaction energies, such as van der Waals energy, electrostatic energy, polar solvation energy, and solvent accessible surface area (SASA) energy, have been also calculated for all the Tools Shapes complexes (Figure 14). Results indicate that van der Waals, electrostatic, and SASA energy negatively contribute to the total interaction energy, while only polar solvation energy positively contributes to the total free binding energy. In particular, the contribution of van der Waals interactions is much greater than that of the electrostatic interactions in all cases except the complexes arginine kinase-delphinidin 3-glucoside and histone subunit-myricetin 3glucoside. Furthermore, the contribution of SASA energy is relatively small when compared to the total binding energy. The negative value of van der Waals energy also points to the significant hydrophobic interaction between the ligands and the enzymes [22]. solvation energy positively contributes to the total free binding energy. In particular, the contribution of van der Waals interactions is much greater than that of the electrostatic interactions in all cases except the complexes arginine kinase-delphinidin 3-glucoside and histone subunit-myricetin 3′-glucoside. Furthermore, the contribution of SASA energy is relatively small when compared to the total binding energy. The negative value of van der Waals energy also points to the significant hydrophobic interaction between the ligands and the enzymes [22].

Principal Component Analysis (PCA)
Principal component analysis is a method that utilizes linear combinations of measured variables, which allows for the reduction of the dimensionality of data and helps identify the principal sources of variation. In molecular dynamics simulations, PCA is a popular method to account for the essential dynamics of the system on a lowdimensional free energy landscape [23]. To analyze the collective motion of all complexes, PCA analysis based on C-a atoms has been performed. It was observed that the first few eigenvectors of the principal components (PCs) of the structures play an important role and describe the overall motions of the entire system. These data suggest that kaempferol-3,7-di-O-glucoside ID has formed very stable complexes with arginine kinase and myricetin 3′-glucoside ID with histone subunit III, which can be considered as a lead compound ( Figure 15).

Principal Component Analysis (PCA)
Principal component analysis is a method that utilizes linear combinations of measured variables, which allows for the reduction of the dimensionality of data and helps identify the principal sources of variation. In molecular dynamics simulations, PCA is a popular method to account for the essential dynamics of the system on a low-dimensional free energy landscape [23]. To analyze the collective motion of all complexes, PCA analysis based on C-a atoms has been performed. It was observed that the first few eigenvectors of the principal components (PCs) of the structures play an important role and describe the overall motions of the entire system. These data suggest that kaempferol-3,7-di-O-glucoside ID has formed very stable complexes with arginine kinase and myricetin 3 -glucoside ID with histone subunit III, which can be considered as a lead compound ( Figure 15). solvation energy positively contributes to the total free binding energy. In particular, the contribution of van der Waals interactions is much greater than that of the electrostatic interactions in all cases except the complexes arginine kinase-delphinidin 3-glucoside and histone subunit-myricetin 3′-glucoside. Furthermore, the contribution of SASA energy is relatively small when compared to the total binding energy. The negative value of van der Waals energy also points to the significant hydrophobic interaction between the ligands and the enzymes [22]. Figure 14. Representation of the van der Waals, electrostatic, polar solvation, SASA, and binding energy for docked compounds into: Acetylcholine esterase, Actin, Α-tubulin, Arginine kinase and histone subunit III with inhibitors (disodium 5′-inosinate and petunidin 3-glucoside), (calcium 5′guanylate D and chlorogenic acid), chlorogenic acid, (kaempferol-3,7-di-o-glucoside I, hyperoside D, delphinidin 3-glucoside ID), and (myricetin 3′-glucoside ID, hyperoside D) complexes.

Principal Component Analysis (PCA)
Principal component analysis is a method that utilizes linear combinations of measured variables, which allows for the reduction of the dimensionality of data and helps identify the principal sources of variation. In molecular dynamics simulations, PCA is a popular method to account for the essential dynamics of the system on a lowdimensional free energy landscape [23]. To analyze the collective motion of all complexes, PCA analysis based on C-a atoms has been performed. It was observed that the first few eigenvectors of the principal components (PCs) of the structures play an important role and describe the overall motions of the entire system. These data suggest that kaempferol-3,7-di-O-glucoside ID has formed very stable complexes with arginine kinase and myricetin 3′-glucoside ID with histone subunit III, which can be considered as a lead compound ( Figure 15). Figure 15. The PCA analysis, the plot of eigenvalues vs. eigenvectors have been considered. Figure 15. The PCA analysis, the plot of eigenvalues vs. eigenvectors have been considered.
Since it was previously found that the first five eigenvectors constitute the majority portion of the total dynamics of the whole system, we plotted only the first two eigenvectors against each other, where each dot represents correlated motions ( Figure 16). The well-stable clustered dots signify the more stable structure, and low-clustered dots represent the weaker stable structure.
Since it was previously found that the first five eigenvectors constitute the majority portion of the total dynamics of the whole system, we plotted only the first two eigenvectors against each other, where each dot represents correlated motions ( Figure  16). The well-stable clustered dots signify the more stable structure, and low-clustered dots represent the weaker stable structure.

Database Search, Structural Modeling, and Model Validation
All protein sequences were obtained from NCBI (https://www.ncbi.nlm.nih. accessed on 10 August 2022) in FASTA format and are mentioned by their Gen accession number in Table 2. The Liriomyza trifolii NCBI taxonomy (tax ID: 32 proteins were selected by searching all of the sequential homolog and orthologs u NCBI Blast server [24] with the default values, and against the nonredundant pr sequences. The sequences were retrieved in the FASTA format as an amino sequence. The initial atomic structures of the proteins, based on homology mode were built using the Swissmodel server (https://Swissmodel.expasy.org/, accessed o August 2022). In this study, a sequence of Blast-P similarities for recognition of cl related structural homologs in Liriomyza trifolii was queried against a PDB database The first hit on the annotation Blast-p was obtained to identify the templates base PDB template ID. The Protein Data Bank collected the PDB file of the templates and aligned using BLAST. The Swissmodel server used the target sequence file alignment file, the PDB file for the prototype, and all the template proteins to build homology model. Homology models with a score of <−4 were chosen. The optim models (acetylcholinesterase, α-tubulin, actin, arginine kinase, histone subunit III, shock protein 90 (Hsp90), and elongation factor 1-alpha) were found to be suitable b on several qualitative background checks, including the PROCHECK (PDBSum) Swissmodel server (https://saves.mbi.ucla.edu/, accessed on 10 August 2022). Ramachandran plot evaluated that the predicted models were closer to the tem with (99.1%, 92.6%, 86.7%, 84.4%, 88.6%, 88.4%) residues lying in the favored reg The ERRAT score values of 99.1304, 89.7527, 96.4539, 82.0707, 96.5217, 90.9774, QMEAN score indicated that the predicted models were reliable and satisfactor they are higher than the ideal values of the QMEAN score <−4, and ERRAT around for a model with a satisfactory resolution [24].

Database Search, Structural Modeling, and Model Validation
All protein sequences were obtained from NCBI (https://www.ncbi.nlm.nih.gov/, accessed on 10 August 2022) in FASTA format and are mentioned by their Gen Bank accession number in Table 2. The Liriomyza trifolii NCBI taxonomy (tax ID: 32264) proteins were selected by searching all of the sequential homolog and orthologs using NCBI Blast server [24] with the default values, and against the nonredundant protein sequences. The sequences were retrieved in the FASTA format as an amino-acid sequence. The initial atomic structures of the proteins, based on homology modeling, were built using the Swissmodel server (https://Swissmodel.expasy.org/, accessed on 10 August 2022). In this study, a sequence of Blast-P similarities for recognition of closely related structural homologs in Liriomyza trifolii was queried against a PDB database [18]. The first hit on the annotation Blast-p was obtained to identify the templates based on PDB template ID. The Protein Data Bank collected the PDB file of the templates and was aligned using BLAST. The Swissmodel server used the target sequence file, the alignment file, the PDB file for the prototype, and all the template proteins to build the homology model. Homology models with a score of <−4 were chosen. The optimized models (acetylcholinesterase, α-tubulin, actin, arginine kinase, histone subunit III, heat shock protein 90 (Hsp90), and elongation factor 1-alpha) were found to be suitable based on several qualitative background checks, including the PROCHECK (PDBSum) and Swissmodel server (https://saves.mbi.ucla.edu/, accessed on 10 August 2022). The Ramachandran plot evaluated that the predicted models were closer to the template with (99.1%, 92.6%, 86.7%, 84.4%, 88.6%, 88.4%) residues lying in the favored regions. The ERRAT score values of 99.1304, 89.7527, 96.4539, 82.0707, 96.5217, 90.9774, and QMEAN score indicated that the predicted models were reliable and satisfactory, as they are higher than the ideal values of the QMEAN score <−4, and ERRAT around 95% for a model with a satisfactory resolution [24].

Preparation of Proteins and Ligands
The sequences of the Liriomyza trifolii proteins (acetylcholinesterase, actin, α-tubulin, arginine kinase, elongation factor 1-alpha, Hsp90, and histone subunit III) with GenBank accession no. number (CAI30732.1, ARQ84036.1, ARQ84030.1, ARQ84038.1, ARQ84034.1, AGI19327.1, ARQ84032.1, ABL07756.1, respectively) were obtained from NCBI. The protein sequences were retrieved in the FASTA format. The 3-D structures of proteins were built using the Swissmodel server (https://Swissmodel.expasy.org/, accessed on 10 August 2022). Here, proteins were selected as target receptor proteins and were imported to the 3-D refine server to perform energy minimization for the six proteins (http://sysbio. rnet.missouri.edu/3Drefine/, accessed on 10 August 2022). During docking studies, all water molecules and ligands were removed, and hydrogen atoms were added to the target proteins. The docking system was built using SAMSON software 2020 (https: //www.samson-connect.net/, accessed on 10 August 2022). The structures were prepared using the protein preparation wizard of the Autodock Vina extension of SAMSON 2020 software. The X, Y, and Z dimensions of the receptor grid, used for the blind docking of ligands to proteins, are reported in Table 3. The ligands were retrieved from the PubChem database in SDF format. Subsequently, each ligand was converted into MOL2 format using OpenBabel software (http://openbabel.org/wiki/Main_Page, accessed on 10 August 2022), followed by an energy minimization at pH 7.0 ± 2.0 in SAMSON software.

Binding Site Prediction and Protein-Ligand Docking
Discovery studio software and SAMSON software were used for binding site prediction. SAMSON software uses Autodock Vina as an extension to maximize the accuracy of these predictions while minimizing computer run-time [25]. It uses the interaction energy between the protein and a simple van der Waals probe to locate energetically favorable binding sites. The program is based on quantum mechanics, and it predicts the potential affinity, molecular structure, geometry optimization of the structure, vibration frequencies of coordinates of atoms, bond length, and bond angle. Following an exhaustive search, 100 poses were analyzed, and the best scoring poses were used to calculate the binding affinity of the ligands. The ligands that tightly bind to a target protein with high scores were selected in Table 3. The proteins were docked against a variety of bioactive compounds that are phytochemical components from the HPLC of leaves of Phaseolus vulgaris (ref) and yeast extract using SAMSON software [21]. The 2-D interaction was carried out to find favorable binding geometries of the ligand with the proteins using Discovery Studio software. Thus, the 2-D interaction images of the docked protein-ligand complexes with high scores to the predicted active sites were obtained.

Protein Ligand Interaction Using SAMSON and Discovery Studio Software
The ligands were docked with the target proteins (acetylcholinesterase, actin, αtubulin, arginine kinase, elongation factor 1-alpha, Hsp90, and histone subunit III), and the best docking poses were identified. Figures 1 and 2 show the 2-D and 3-D structures of the binding poses of the compounds, respectively.

Protein-Protein Interaction Network
The Liriomyza trifolii proteins were submitted to the server for functional interaction associated network between partners for the STRING (Research Online of Interacting Genes/Proteins Data Basis version 10.0)13 [24]. The interactions were examined at medium and high confidences.

Molecular Dynamics Simulation
The molecular dynamic approach is widely used to assess atoms' behavior and structural stability, and to study the conformational changes at an atomic level. Understanding the stability of protein upon ligand binding is significantly improved by molecular dynamics simulation studies. Gromacs 4.6.2 [26] with GROMOS96 54a7 force field [27] was used for MD simulation studies of two systems, at 50 ns each. The ProdrG2 Server was used to generate the topology of the analysis of enzyme-ligand complexes. Each system was placed in the center of the cubic box, with a distance of 1.0 nm between the enzyme and the edge of the simulation box. Each system was solvated with explicit water molecules.
Before proceeding towards energy minimization, all systems were neutralized by adding Na + and Cl − ions, accordingly. The steepest descent method was used for the energy minimization of each system. MD simulations with NVT (isochoric-isothermal) and NPT (isobaric-isothermal) ensembles (N 1 4 constant particle number, V 1 4 Volume, P 1 4 Pressure, T 1 4 Temperature) were performed for 1 ns, each, to equilibrate the enzyme-ligand system for constant volume, pressure (1 atm), and temperature (300 K). To calculate electrostatic interaction, the Particle Mesh Ewald (PME) algorithm [25] was used with a grid spacing of 1.6 Å and a cutoff of 10 Å, and the LINCS method was used to restrain the bond length. Finally, the trajectories were saved at every 2-fs time step and the production MD simulation of the enzyme-ligand complex was performed for 50 ns [28].

Conclusions
This study presented an array of naturally occurring, nontoxic, easily extractable, low-cost ligands that show great potential as inhibitors of a variety of proteins found in the gut of the polyphagous pest L. trifolii that is known to attack a myriad of crops. The target proteins are acetylcholinesterase, actin, α-tubulin, arginine kinase, and histone receptor III subtypes. The proposed inhibitors or inhibitor combinations are (disodium 5 -inosinate and praliciguat), (calcium 5 -guanylate and chlorogenic acid), chlorogenic acid alone, (kaempferol-3,7-di-O-glucoside with hyperoside and delphinidin 3-glucoside), and (myricetin 3 -glucoside and hyperoside), respectively. In lieu of an experimentally available structure of the target proteins, the initial models of the protein of L. trifolii origin were constructed using homology modeling. The analyses used in this investigation included structural modeling, binding site interaction prediction, molecular docking free energy calculations, binding pose analysis, dynamic stability and conformational perturbation analysis, radius of gyration analysis, hydrogen bond analysis, and molecular mechanics PBSA free energy calculations. The results demonstrated that the proposed inhibitors formed stable complexes with their target proteins while also having great potentials for inhibitory activity. All five ligand-protein complexes have favorable parameter values in RMSD, RMSF, RoG, intermolecular hydrogen bonding, and binding free energy for 50 ns. Trajectories analysis showed that the studied complexes displayed structural stability during the MD runs.