key: cord-0069848-eyjy0t28
authors: Alrumaihi, Faris
title: A Comprehensive Computational Screening of Phytochemicals Derived from Saudi Medicinal Plants against Human CC Chemokine Receptor 7 to Identify Potential Anti-Cancer Therapeutics
date: 2021-10-21
journal: Molecules
DOI: 10.3390/molecules26216354
sha: d3df502e84ef21679515812f4e4f7640501f14bf
doc_id: 69848
cord_uid: eyjy0t28

Homeostatic trafficking of immune cells by CC chemokine receptor 7 (CCR7) keeps immune responses and tolerance in a balance. The involvement of this protein in lymph node metastasis in cancer marks CCR7 as a penitential drug target. Using the crystal structure of CCR7, herein, a comprehensive virtual screening study is presented to filter novel strong CCR7 binding phytochemicals from Saudi medicinal plants that have a higher binding affinity for the intracellular allosteric binding pocket. By doing so, three small natural molecules named as Hit-1 (1,8,10-trihydroxy-3-methoxy-6-methylanthracen-9(4H)-one), Hit-2 (4-(3,4-dimethoxybenzyl)-3-(4-hydroxy-3-methoxybenzyl)dihydrofuran-2(3H)-one), and Hit-3 (10-methyl-12,13-dihydro-[1,2]dioxolo[3,4,5-de]furo[3,2-g]isochromeno[4,3-b]chromen-8-ol) are predicted showing strong binding potential for the CC chemokine receptor 7 allosteric pocket. During molecular dynamics simulations, the compounds were observed in the formation of several chemical bonding of short bond distances. Additionally, the molecules remained in strong contact with the active pocket residues and experienced small conformation changes that seemed to be mediated by the CCR7 loops to properly engage the ligands. Two types of binding energy methods (MM/GBPBSA and WaterSwap) were additionally applied to further validate docking and simulation findings. Both analyses complement the good affinity of compounds for CCR7, the electrostatic and van der Waals energies being the most dominant in intermolecular interactions. The active pocket residue’s role in compounds binding was further evaluated via alanine scanning, which highlighted their importance in natural compounds binding. Additionally, the compounds fulfilled all drug-like rules: Lipinski, Ghose, Veber, Egan, and Muegge passed many safety parameters, making them excellent anti-cancer candidates for experimental testing.

The human immune system has the potential to fight against pathogens without harming normal cells and tissues. One of the most important aspects is chemotactic trafficking, which guides immune cells to perform their action at appropriate places and times [1] . Chemotactic trafficking is regulated by 20 G protein-coupled receptors (GPCRs) and more than 40 chemokines. Inflammatory chemokines and receptors are generated in response to inflammatory stimuli, whereas homeostatic chemokines are produced continuously, which direct cells to certain organs [2] . Chemokine ligands such as CCL19 and CCL21 bind the CC chemokine receptor 7 (CCR7) to direct B cells, T cells, and dendritic cells to lymph nodes all across the body. The chemotactic trafficking creates cellular pathways, wherein the inflammatory immune response is induced by inflammation chemokines and target specific cells. CCL19 and CCL21 are the cellular cation receptors that, when attached to

Filtration of the in-house generated drug library was first performed to remove nondrug-like molecules. This step holds great importance, as drug-like molecules enhance the chances to forward only those molecules that could lead to successful drug discovery. Out of 1741, a total of 1080 molecules were found to completely follow the Lipinski rule of five [19] . Only drug-like molecules were selected and used in the virtual screening process.

Target-based virtual screening was performed to identify the best binding from the Saudi medicinal plant library. The filtered compounds from the library done in the previous step were used only in the virtual screening process. The first three Hits that were ranked as best binding molecules considering their lower binding energy value were selected for binding mode and interactions analysis (Table 1) . Hit-3 -7.00

All three compounds bound deep inside the binding pocket near the TM6 IP linker ( Figure 1A ). The first Hit molecule is 1,8,10-trihydroxy-3-methoxy-6-methylanthracen-9(4H)-one, where major contribution was noticed from the central oxygen atoms and hydroxyl groups of the benzene rings ( Figure 1B ). The compound is strongly engaged by three hydrogen bonds of residues (residue number in crystal CCR7); Lys39 (Lys90), Met41 (Met92), and Asp102 (Asp153), with bond distances of 2.1 A, 1.8 A, and 1.9 A, respectively. Besides hydrogen bonding, several van der Waals bondings were witnessed, for example, Arg37 (Arg88), Leu38 (Leu89), Thr40 (Thr91), Arg103 (Arg154), Ala106 (Ala157), His109 (His160), Arg112 (Arg163), and Arg671 (Arg722). The Hit-2 molecule is 4-(3,4-dimethoxybenzyl)-3-(4-hydroxy-3-methoxybenzyl)dihydrofuran-2(3H)-one. After thorough search in PubChem, it has been found that it is Arctigenin, a potent ingredient of Arctium lappa. This compound is reported to show anti-inflammatory activity by blocking the function of ROS-dependent signal transducers and activators of transcription. It is also reported for antibacterial, antifungal, anti-cancer, and anti-diabetic activities [20, 21] . The PubChem ID for this compound is 230232. The molecule was observed to dock at the same position as Hit-1; however, due to the difference in the binding conformation, the binding interaction pattern of Hit-2 is slightly different ( Figure 1C ). The 1,2-dimethoxy-4-methylbenzene moiety of the compound was seen in maximum interactions with the enzyme. One strong hydrogen bond interaction of the compound with Arg103 (Arg154) was seen. The rest of the compound structure (3-(4-hydroxy-3-methoxybenzyl)dihydrofuran-2(3H)-one) is engaged by many other van der Waals and very weak alky and pi-alkyl bonding. The third molecule is 10-methyl-12,13-dihydro- [1, 2] dioxolo [3,4,5-de] furo [3,2-g] isochromeno [4,3b] chromen-8-ol. The interactions of this molecule with the enzyme active pocket are dominated by very weak van der Waals, carbon hydrogen, pi-cation, pi-anion, pi-sigma, pisulfur, pi-pi stacked, alkyl, and pi-alkyl interactions ( Figure 1D ). The co-crystalized Cmp2105 compound with CCR7 outcompetes the CCL19 (native protein ligand) and has Hit-3 -7.00

All three compounds bound deep inside the binding pocket near the TM6 IP linker ( Figure 1A ). The first Hit molecule is 1,8,10-trihydroxy-3-methoxy-6-methylanthracen-9(4H)-one, where major contribution was noticed from the central oxygen atoms and hydroxyl groups of the benzene rings ( Figure 1B ). The compound is strongly engaged by three hydrogen bonds of residues (residue number in crystal CCR7); Lys39 (Lys90), Met41 (Met92), and Asp102 (Asp153), with bond distances of 2.1 A, 1.8 A, and 1.9 A, respectively. Besides hydrogen bonding, several van der Waals bondings were witnessed, for example, Arg37 (Arg88), Leu38 (Leu89), Thr40 (Thr91), Arg103 (Arg154), Ala106 (Ala157), His109 (His160), Arg112 (Arg163), and Arg671 (Arg722). The Hit-2 molecule is 4-(3,4-dimethoxybenzyl)-3-(4-hydroxy-3-methoxybenzyl)dihydrofuran-2(3H)-one. After thorough search in PubChem, it has been found that it is Arctigenin, a potent ingredient of Arctium lappa. This compound is reported to show anti-inflammatory activity by blocking the function of ROS-dependent signal transducers and activators of transcription. It is also reported for antibacterial, antifungal, anti-cancer, and anti-diabetic activities [20, 21] . The PubChem ID for this compound is 230232. The molecule was observed to dock at the same position as Hit-1; however, due to the difference in the binding conformation, the binding interaction pattern of Hit-2 is slightly different ( Figure 1C ). The 1,2-dimethoxy-4-methylbenzene moiety of the compound was seen in maximum interactions with the enzyme. One strong hydrogen bond interaction of the compound with Arg103 (Arg154) was seen. The rest of the compound structure (3-(4-hydroxy-3-methoxybenzyl)dihydrofuran-2(3H)-one) is engaged by many other van der Waals and very weak alky and pi-alkyl bonding. The third molecule is 10-methyl-12,13-dihydro- [1, 2] dioxolo [3,4,5-de] furo [3,2-g] isochromeno [4,3b] chromen-8-ol. The interactions of this molecule with the enzyme active pocket are dominated by very weak van der Waals, carbon hydrogen, pi-cation, pi-anion, pi-sigma, pisulfur, pi-pi stacked, alkyl, and pi-alkyl interactions ( Figure 1D ). The co-crystalized Cmp2105 compound with CCR7 outcompetes the CCL19 (native protein ligand) and has Hit-3 -7.00

All three compounds bound deep inside the binding pocket near the TM6 IP linker ( Figure 1A ). The first Hit molecule is 1,8,10-trihydroxy-3-methoxy-6-methylanthracen-9(4H)-one, where major contribution was noticed from the central oxygen atoms and hydroxyl groups of the benzene rings ( Figure 1B ). The compound is strongly engaged by three hydrogen bonds of residues (residue number in crystal CCR7); Lys39 (Lys90), Met41 (Met92), and Asp102 (Asp153), with bond distances of 2.1 A, 1.8 A, and 1.9 A, respectively. Besides hydrogen bonding, several van der Waals bondings were witnessed, for example, Arg37 (Arg88), Leu38 (Leu89), Thr40 (Thr91), Arg103 (Arg154), Ala106 (Ala157), His109 (His160), Arg112 (Arg163), and Arg671 (Arg722). The Hit-2 molecule is 4-(3,4-dimethoxybenzyl)-3-(4-hydroxy-3-methoxybenzyl)dihydrofuran-2(3H)-one. After thorough search in PubChem, it has been found that it is Arctigenin, a potent ingredient of Arctium lappa. This compound is reported to show anti-inflammatory activity by blocking the function of ROS-dependent signal transducers and activators of transcription. It is also reported for antibacterial, antifungal, anti-cancer, and anti-diabetic activities [20, 21] . The PubChem ID for this compound is 230232. The molecule was observed to dock at the same position as Hit-1; however, due to the difference in the binding conformation, the binding interaction pattern of Hit-2 is slightly different ( Figure 1C ). The 1,2-dimethoxy-4-methylbenzene moiety of the compound was seen in maximum interactions with the enzyme. One strong hydrogen bond interaction of the compound with Arg103 (Arg154) was seen. The rest of the compound structure (3-(4-hydroxy-3-methoxybenzyl)dihydrofuran-2(3H)-one) is engaged by many other van der Waals and very weak alky and pi-alkyl bonding. The third molecule is 10-methyl-12,13-dihydro- [1, 2] dioxolo [3,4,5-de] furo [3,2-g] isochromeno [4,3b] chromen-8-ol. The interactions of this molecule with the enzyme active pocket are dominated by very weak van der Waals, carbon hydrogen, pi-cation, pi-anion, pi-sigma, pisulfur, pi-pi stacked, alkyl, and pi-alkyl interactions ( Figure 1D ). The co-crystalized Cmp2105 compound with CCR7 outcompetes the CCL19 (native protein ligand) and has -7.00

All three compounds bound deep inside the binding pocket near the TM6 IP linker ( Figure 1A ). The first Hit molecule is 1,8,10-trihydroxy-3-methoxy-6-methylanthracen-9(4H)-one, where major contribution was noticed from the central oxygen atoms and hydroxyl groups of the benzene rings ( Figure 1B ). The compound is strongly engaged by three hydrogen bonds of residues (residue number in crystal CCR7); Lys39 (Lys90), Met41 (Met92), and Asp102 (Asp153), with bond distances of 2.1 A, 1.8 A, and 1.9 A, respectively. Besides hydrogen bonding, several van der Waals bondings were witnessed, for example, Arg37 (Arg88), Leu38 (Leu89), Thr40 (Thr91), Arg103 (Arg154), Ala106 (Ala157), His109 (His160), Arg112 (Arg163), and Arg671 (Arg722). The Hit-2 molecule is 4-(3,4dimethoxybenzyl)-3-(4-hydroxy-3-methoxybenzyl)dihydrofuran-2(3H)-one. After thorough search in PubChem, it has been found that it is Arctigenin, a potent ingredient of Arctium lappa. This compound is reported to show anti-inflammatory activity by blocking the function of ROS-dependent signal transducers and activators of transcription. It is also reported for antibacterial, antifungal, anti-cancer, and anti-diabetic activities [20, 21] . The PubChem ID for this compound is 230232. The molecule was observed to dock at the same position as Hit-1; however, due to the difference in the binding conformation, the binding interaction pattern of Hit-2 is slightly different ( Figure 1C ). The 1,2-dimethoxy-4-methylbenzene moiety of the compound was seen in maximum interactions with the enzyme. One strong hydrogen bond interaction of the compound with Arg103 (Arg154) was seen. The rest of the compound structure (3-(4-hydroxy-3-methoxybenzyl)dihydrofuran-2(3H)-one) is engaged by many other van der Waals and very weak alky and pi-alkyl bonding. The third molecule is 10-methyl-12,13-dihydro- [1, 2] dioxolo [3,4,5-de] furo [3,2g] isochromeno [4,3-b] chromen-8-ol. The interactions of this molecule with the enzyme active pocket are dominated by very weak van der Waals, carbon hydrogen, pi-cation, pi-anion, pi-sigma, pi-sulfur, pi-pi stacked, alkyl, and pi-alkyl interactions ( Figure 1D ). The co-crystalized Cmp2105 compound with CCR7 outcompetes the CCL19 (native protein ligand) and has IC 50 value of 35 nM in membrane competition experiments. Chemically, Cmp2105 contains thiadiazole-dioxide core motif along with two amine-linked substituents. The CCR7-Cmp2105 crystal structure demonstrated that the compound is bound to the same intracellular part as found for the compounds screened in this study. The binding pocket of the compounds is located between TM1, TM2, TM3, and TM6 end and between the H8 and TM7 loop of the CCR7. The Cmp2105 and all three Hits were reported to interact with the same active pocket residues. The Cmp2105 interacts with Val79, Thr82, Phe86, Thr91, Thr93, Arg154, and Tyr326. All the screened Hits as well as the Cmp2105 interact with highly conserved TM7-H8 patch residues [1] . The mentioned patch is located within the allosteric binding pocket and represents a promising hotspot for targeting chemokine receptors. Thus, it can be speculated that the filtered compounds might serve as broad spectrum blockers of CCR. In short, it can be concluded that the compounds showed favorable interactions with the enzyme active site residues and showed good binding affinity. Hit-1 and Hit-3 seem to be novel, as no such structures were found after a thorough search against PubChem and ChEMBL databases. IC50 value of 35 nM in membrane competition experiments. Chemically, Cmp2105 contains thiadiazole-dioxide core motif along with two amine-linked substituents. The CCR7-Cmp2105 crystal structure demonstrated that the compound is bound to the same intracellular part as found for the compounds screened in this study. The binding pocket of the compounds is located between TM1, TM2, TM3, and TM6 end and between the H8 and TM7 loop of the CCR7. The Cmp2105 and all three Hits were reported to interact with the same active pocket residues. The Cmp2105 interacts with Val79, Thr82, Phe86, Thr91, Thr93, Arg154, and Tyr326. All the screened Hits as well as the Cmp2105 interact with highly conserved TM7-H8 patch residues [1] . The mentioned patch is located within the allosteric binding pocket and represents a promising hotspot for targeting chemokine receptors. Thus, it can be speculated that the filtered compounds might serve as broad spectrum blockers of CCR. In short, it can be concluded that the compounds showed favorable interactions with the enzyme active site residues and showed good binding affinity. Hit-1 and Hit-3 seem to be novel, as no such structures were found after a thorough search against PubChem and ChEMBL databases. Table mentioning residues numbering in the docked CCR7 structure and crystal CCR7 structure is also provided.

To examine the stability and dynamics of CCR7-compound complexes, 200 ns of molecular simulation was performed. Different statistical analysis, including root mean square RMSD [22] , root mean square fluctuation (RMSF) [23] , radius of gyration (RoG) [24] , and hydrogen bond analysis, were performed on the simulation trajectories to decipher intermolecular strength of interaction and stable dynamics of complexes. The lower the RMSD, more stable the CCR7 will be in the presence of compounds, and vice versa. The RMSD of all three complexes was observed to be stable until 75 ns and then behaved differently for each complex and stayed constant after 150 ns till the simulation end time Table mentioning residues numbering in the docked CCR7 structure and crystal CCR7 structure is also provided.

To examine the stability and dynamics of CCR7-compound complexes, 200 ns of molecular simulation was performed. Different statistical analysis, including root mean square RMSD [22] , root mean square fluctuation (RMSF) [23] , radius of gyration (RoG) [24] , and hydrogen bond analysis, were performed on the simulation trajectories to decipher intermolecular strength of interaction and stable dynamics of complexes. The lower the RMSD, more stable the CCR7 will be in the presence of compounds, and vice versa. The RMSD of all three complexes was observed to be stable until 75 ns and then behaved differently for each complex and stayed constant after 150 ns till the simulation end time with no substantial variation. Between the period of 75 ns to 150 ns, the complexes experienced some deviations that were due to the protein loops, which allowed better accommodation of the compounds at the docked site. The net RMSD of Hit-1, Hit-2, and Hit-3 is 0.9 A, 1.3 A, and 1.26 A, respectively. The RMSD plot for the systems can be visualized in Figure 2A . Next, RMSF for the complexes was calculated to determine residue specific deviation in the presence of compound at the active pocket of the enzyme. As can be seen in Figure 2B , the majority of the receptor residues are in good stable range, with RMSF < 2 A. Some residues, mainly those that are part of the CCR7 loops, drive some deviations responsible for higher RMSF. The net RMSD of Hit-1, Hit-2, and Hit-3 is 0.8 A, 1.1 A, and 1.3 A, respectively. The CCR7 compactness during simulation time was then calculated through RoG. Higher RoG values indicate a lower compactness of receptor molecule, whereas lower RoG implies a higher compact nature of CCR7. As can be noticed in Figure 2C , all three systems revealed a very stable RoG plot, which is a strong indication of receptor equilibrium in the presence of compounds during simulation time. These results also complement the RMSD analysis, which depicted the stable behavior of the systems. The mean RoG values of the systems are: Hit (22.5 A), Hit-2 (23.4 A), and Hit-3 (24.42 A). The intermolecular interaction strength was determined through hydrogen bond analysis. Hydrogen bonds are key interactions that keep the docked ligands intact at the docked site and play a significant role in long term receptor-ligand stability. All three complexes reported multiple hydrogen bonds during simulation ( Figure 2D ). This shows that hydrogen bonding is significant in making the compounds strongly attached to the CCR7 active pocket residues. with no substantial variation. Between the period of 75 ns to 150 ns, the complexes experienced some deviations that were due to the protein loops, which allowed better accommodation of the compounds at the docked site. The net RMSD of Hit-1, Hit-2, and Hit-3 is 0.9 A, 1.3 A, and 1.26 A, respectively. The RMSD plot for the systems can be visualized in Figure 2A . Next, RMSF for the complexes was calculated to determine residue specific deviation in the presence of compound at the active pocket of the enzyme. As can be seen in Figure 2B , the majority of the receptor residues are in good stable range, with RMSF < 2 A. Some residues, mainly those that are part of the CCR7 loops, drive some deviations responsible for higher RMSF. The net RMSD of Hit-1, Hit-2, and Hit-3 is 0.8 A, 1.1 A, and 1.3 A, respectively. The CCR7 compactness during simulation time was then calculated through RoG. Higher RoG values indicate a lower compactness of receptor molecule, whereas lower RoG implies a higher compact nature of CCR7. As can be noticed in Figure  2C , all three systems revealed a very stable RoG plot, which is a strong indication of receptor equilibrium in the presence of compounds during simulation time. These results also complement the RMSD analysis, which depicted the stable behavior of the systems. Hydrogen bonds are key interactions that keep the docked ligands intact at the docked site and play a significant role in long term receptor-ligand stability. All three complexes reported multiple hydrogen bonds during simulation ( Figure 2D ). This shows that hydrogen bonding is significant in making the compounds strongly attached to the CCR7 active pocket residues. 

The radial distribution function analysis was further performed to highlight intermolecular bonds that keep the compounds in contact with enzyme active pocket residues. The intermolecular interactions were determined through an in-house script that filtered all the enzyme residues present within 3 A at the compounds bounded site. Only interactions that have high density during simulation time were used in radial distribution plotting. As can be seen in Figure 3A , for Hit-1, three interactions were reported to show high radii density. All three interactions are hydrogen bonds and are also predicted by docking and simulation studies to play role in CCR7-Hit-1 stable docking. A higher interaction density was seen between Asp102 (Asp153) and Hit-1 hydrogen atom with g(r) value > 3 at a bond distance length of~2 A. The Met41 (Met92) and Hit-O atom interaction has a g(r) value of 1.5 and was at 2.1 A. Similarity, for Lys39 and Hit-1, hydrogen atom maximum density plotted is 0.8 at a distance of 2.3 A. In the case of the Hit-2 molecule, Asp102 (Asp153) interaction with compound hydrogen is reported to reveal high density distribution (>1.12) at a distance of 1.9 A ( Figure 3B ). For Hit-3, the intermolecular interaction between Asp102 (Asp153) and Hit-3 hydrogen shows a maximum g(r) value at 2 A ( Figure 3C ). 

The radial distribution function analysis was further performed to highlight intermolecular bonds that keep the compounds in contact with enzyme active pocket residues. The intermolecular interactions were determined through an in-house script that filtered all the enzyme residues present within 3 A at the compounds bounded site. Only interactions that have high density during simulation time were used in radial distribution plotting. As can be seen in Figure 3A , for Hit-1, three interactions were reported to show high radii density. All three interactions are hydrogen bonds and are also predicted by docking and simulation studies to play role in CCR7-Hit-1 stable docking. A higher interaction density was seen between Asp102 (Asp153) and Hit-1 hydrogen atom with g(r) value > 3 at a bond distance length of ~2 A. The Met41 (Met92) and Hit-O atom interaction has a g(r) value of 1.5 and was at 2.1 A. Similarity, for Lys39 and Hit-1, hydrogen atom maximum density plotted is 0.8 at a distance of 2.3 A. In the case of the Hit-2 molecule, Asp102 (Asp153) interaction with compound hydrogen is reported to reveal high density distribution (>1.12) at a distance of 1.9 A ( Figure 3B ). For Hit-3, the intermolecular interaction between Asp102 (Asp153) and Hit-3 hydrogen shows a maximum g(r) value at 2 A ( Figure 3C ). 

Both MM/GBSA and MM/PBSA methods are frequently applied in drug design works to examine drug binding affinity towards a given biological macromolecule [25] . In particular, the benefit of the methods is the use of modest computational power compared to alchemical binding free energy methods. Additionally, the methods are end point methods to test the different binding energies produced by the interaction of the compounds with the receptor. All three virtually screened compounds reported robust electrostatic and van der Waals energy; in particular, the former contributes more in the case of Hit-1, while the latter contributes significantly for Hit-2 and Hit-3. The different energy values of the compounds can be seen in Table 2 . The electrostatic and van der Waal energy combine to make the gas phase energy component of the MM/GBPBSA. The net gas phase energy of Hit-1, Hit-2, and Hit-3 is -56.25 kcal/mol, −56.49 kcal/mol, and 22.39 kcal/mol, respectively. The solvation part of MM/GBPBSA was revealed to contribute negatively to overall binding energy of the complexes. In the case of MM/GBSA, the net solvation energy of Hit-1 is 14.89 kcal/mol, Hit-2 is 11.58 kcal/mol, and Hit-3 is 13.45 kcal/mol, while in MM/PBSA, Hit-1 net solvation energy is 15 kcal/mol, Hit-2 is 12.83 kcal/mol, and Hit-3 is 9.97 kcal/mol. The analysis indicated that all three compounds achieved a stable mode inside the pocket of the enzyme and are engaged by a strong chemical bonding network. However, further structure modification of these compounds can show an increased interaction pattern, which might increase their affinity towards the CCC7 receptor. 

The binding energy of MMGBSA was further decomposed to dissect the energy contribution of enzyme active pocket in recognition and stabilizing compounds [26] . The binding energy of each active pocket residue in MMGBSA is tabulated in Table 3 . The findings are in line with the docking predictions that highlighted several residues in compounds binding. In particular, Lys39, Met41 (Met92), Asp102 (Asp153), Arg103 (Arg154), and Arg112 (Arg163) are reported to show highly stable energy in case of Hit-1. In the case of Hit-2, Lys39 (Lys90), Thr40 (Thr91), Thr42 (Thr93), Asp102 (Asp153), Arg103 (Arg154), Arg671 (Arg722), and Lys735 (Lys786) are high-energy contributing residues. The Hit-3 molecule is strongly engaged by Lys39 (Lys90), Thr40 (Thr91), Met41 (Met92), His109 (His160), Arg103 (Arg154), Glu661 (Glu712), and Arg671 (Arg722).

The MM/GBPBSA method, which uses an implicit water system, picked snapshots at regular intervals of molecular dynamics simulations and did not take into account the enzyme-water and drug-water interaction details. On the other hand, WaterSwap employed an explicit water model and circumvented the limitations of MM/GBPBSA [27] . Overall, the net binding free energy value for all the three systems was recorded to be <−33 kcal/mol (Figure 4) . This signifies the formation of strong intermolecular affinity of the interacting enzyme and drug molecules. Additionally, it can be identified that the difference between the energy values for each complex is <1 kcal/mol, which describes good convergence of the systems. 

The MM/GBPBSA method, which uses an implicit water system, picked snapshots at regular intervals of molecular dynamics simulations and did not take into account the enzyme-water and drug-water interaction details. On the other hand, WaterSwap employed an explicit water model and circumvented the limitations of MM/GBPBSA [27] . Overall, the net binding free energy value for all the three systems was recorded to be < −33 kcal/mol (Figure 4) . This signifies the formation of strong intermolecular affinity of the interacting enzyme and drug molecules. Additionally, it can be identified that the difference between the energy values for each complex is < 1 kcal/mol, which describes good convergence of the systems. 

Entropy energy contribution to the systems net energy was determined through AM-BER normal mode method, which predicted the net binding energy as: Hit-1 (−9.678 kcal/mol), Hit-2 (−6.894 kcal/mol), and Hit-3 (−2.3697 kcal/mol). These values suggested that, upon compounds binding, the targeted macromolecule shift the system to stable microstate level and efficient binding of the compounds to enzyme during simulation. 

Entropy energy contribution to the systems net energy was determined through AMBER normal mode method, which predicted the net binding energy as: Hit-1 (−9.678 kcal/mol), Hit-2 (−6.894 kcal/mol), and Hit-3 (−2.3697 kcal/mol). These values suggested that, upon compounds binding, the targeted macromolecule shift the system to stable microstate level and efficient binding of the compounds to enzyme during simulation.

Computational predictions of compound pharmacokinetics holds significant importance in drug designing process, as the selection of compounds that result in poor clinical outcomes could lead to high economic lost [28, 29] . Detailed pharmacokinetics of the compounds are tabulated in Table 4 . All three compounds are predicted to be good druglike molecules and follow eminent drug rules, such as Lipinski, Ghose, Veber, Egan, and Muegge [30] . All these rules determine the success of compounds in pre-clinical investigation and clinical trials and ensure the drug molecule reaches the market successfully. The compounds also showed good water solubility, which allows maximum drug concentration availability at the target site. Additionally, the compounds have high gastrointestinal absorption. Further, the compounds have zero alerts for pan-assay interference compounds (PAINS) that limited their cross interaction with non-specific biological targets and minimize the chances of false positive results. The topological surface area (TPSA) value of the compounds is within acceptable range, which affirms the compounds' good cell membrane permeability. The compounds showed good synthetic accessibility score, which predicts that the compounds can be easily synthesized for experimental studies. From toxicity perspectives, the compounds in different toxicity tests are found to be non-toxic. They are non-hepatotoxic, non-carcinogenic in mouse, and have low T. pyriformis and minnow toxicity. From an excretion point of view, the compounds are good candidates to be readily excreted from the body. All the predictions indicated the compounds to be good molecules and have high chances to serve as lead molecules. 

Alanine scanning mutagenesis was used to get a better understanding of the contribution of active site amino acid residues to the free energy of interaction between the enzyme and compounds. The analysis suggested that a little alteration in the enzyme active pocket can affect overall binding affinity of enzyme towards ligand binding. With alanine scanning analysis, a decrease in the MMGBSA binding free energy was found for selected mutants (Lys39 (Lys90), Met41 (Met92), and Asp102 (Asp153) in the case of Hit-1; Ly39 (Lys90), Thr40 (Thr91), Thr42 (Thr93), and Arg103 (Arg154) in case of Hit-2; and Arg103 (Arg154) and Glu661 (Glu712) in case of Hit-3. The findings are congruent with the findings of the per residue energy decomposition. The substitution of alanine for active site amino acid residues resulted in a decrease in binding affinity for each active residue.

In total, a library of 1741 compounds was created by collecting phytochemicals from 63 Saudi based medicinal plants. From chemical diversity point of view, the compounds belong to alkaloids, anthraquinones, coumarins, terpenes, flavonoids, phenolic, steroids, saponins, cardenolides, tannins, phenanthrenes, glucosinolates, phenylpropanoids, polyketides, limonoids, etc. The phytochemicals database is in a process of development and will be made available soon on the web. The library was transported to LigandScout4.4 [31] , where the library was filtered based on the Lipinski rule of five [19] to filter out non-druglike molecules. The remaining molecules were then transferred to PyRx 0.8 software [32] , where they were energy minimized and converted to pdbqt format.

Structure-based virtual screening of filtered drug-like molecules was performed against the allosteric binding pocket of CCR7. The CCR7 crystal structure was retrieved from protein data bank in UCSF Chimera version 1.15 using PDB id of 6qzh [1] . The structure was first subjected to a short preparatory phase, where co-associated ligands were removed except for some water molecules at the active pocket, which are relevant from a function perspective. The structure was then energy minimized via steepest descent and conjugate gradient methods. The energy-minimized structure was then used in PyRx software [32] . The grid box was set at the allosteric pocket with dimension set at 15 A along XYZ dimensions. The screening process was then run, and each molecule was ranked based on binding energy in kcal/mol. For each molecule, 20 docked solutions were generated, and the one represented frequently was selected as the best binding mode. The top three best Hits were selected for further molecular dynamics simulation analysis.

Molecular dynamics simulations were accomplished for selected complexes to study the receptor/ligand stability and estimating binding free energy calculation based on the simulation trajectories [33] . Molecular dynamics simulations were carried out using AMBER20 [34] . The force field used for compounds is GAFF [35] , while FF14SB [36] was employed for the CCR7 receptor. The systems were solvated into a TIP3P solvation box, large enough in size to cover the complexes. Energy minimization was done for 1000 steps through the steepest decent and conjugate gradient algorithms. In the heating step, temperature was gradually increased to 310 K through canonical ensemble. The systems were equilibrated for 500 ps, where periodic boundary conditions were applied at constant pressure, and Langevin thermostat is utilized [37] . The production run was accomplished for 200 ns for each system through isothermal and isobaric ensemble. The long-range electrostatic interactions are treated via the particle mesh Ewald method [38] . The SHAKE algorithm [39] was utilized to constrain hydrogen bonds. The simulation trajectories were then analyzed via the CPPTRAJ module [40] to investigate system structure deviations throughout the length of simulation time. The number of hydrogen bonds between the CCR7 and compounds was measured through the hydrogen bond plugin in visual molecular dynamics (VMD) version 1.9.3 software [41] . Further, the intermolecular key residue interactions were plotted in radial distribution function (g(r)) analysis [42] .

The different binding free energies between the CCR7 and compounds were estimated using an MM/GBPBSA method. The analysis was performed through MMPBSA.py module [43] . The number of frames used from simulation trajectories was 100. Entropy of the systems was calculated using AMBER normal mode analysis over 10 frames. Additionally, the key residues were mutated to alanine using alanine scanning to examine the importance of such residues in ligand binding and stability. The detail methodology used for performing MM/GBPBSA and alanine scanning was followed from Asama et al., 2017 [26] .

WaterSwap is a binding free energy predicting method that used Monte Carlo simulation to estimate binding free energies for default 1000 iterations [27, 44] . The water molecule role is usually skipped in MM/GBPBSA, which is important in cases where water molecules act as birding molecules between protein-ligand interactions. Three algorithms were run during the analysis: free energy perturbation (FEP), thermodynamic integration (TI), and Bennett's acceptance ratio (BAR). The difference of 1 kcal/mol among the mentioned methods reflects on a values convergence [45] .

To assess the feasibility of the compounds for experimentations and biological testing, in silico pharmacokinetics, medicinal chemistry, and toxicity were evaluated using pkCSM [46] and SWISSADME [47] servers.

Computer-aided drug design is gaining popularity to virtually screen small drug molecule libraries against a given biological time, thus saving time and economic losses associated with experimental drug design. In this present study, medicinal plant-derived phytochemicals were used in structure-based virtual screenings against an allosteric binding pocket of CCR7 receptor. This filtered three drug-like molecules that showed excellent binding affinity for the CCR7. Dynamically, the compounds' binding was observed to be stable, and the interacting network involved short-distance hydrogen and van der Waals bindings. The binding energies by MM/GBPBSA indicated the gas phase energies (electrostatic and van der Waals) to dominate complex formation, which is additionally validated by the explicit WaterSwap method. The compounds were also predicted to follow druglike rules and have good pharmacokinetics, medicinal chemistry, and low toxicity. These findings suggest the compounds are worthy to be utilized in experimental studies to examine their potency against the CCR7 receptor.

Author Contributions: Conceptualization, methodology, formal analysis, data curation, writingoriginal draft preparation, visualization, project administration, funding acquisition, F.A. All authors have read and agreed to the published version of the manuscript.

Funding: This research received no external funding.

Informed Consent Statement: Not applicable.

The data presented in this study are available within the article.

Structural basis for allosteric ligand recognition in the human CC chemokine receptor 7

Chemotaxis, chemokine receptors and human disease

Rot, A. CCR7 and its ligands: Balancing immunity and tolerance

Cancer and the chemokine network

Dendritic cells are early cellular targets of Listeria monocytogenes after intestinal delivery and are involved in bacterial spread in the host

The chemokine system and cancer

Prediction of lymph node metastasis in colorectal carcinoma by expressionof chemokine receptor CCR7

Involvement of chemokine receptors in breast cancer metastasis

A deep learning model based on sparse auto-encoder for prioritizing cancer-related genes and drug target combinations

Homeostatic chemokine receptors and organ-specific metastasis

Expression of the CC chemokine receptor 7 mediates metastasis of breast cancer to the lymph nodes in mice

Cancer incidence and mortality worldwide: Sources, methods and major patterns in GLOBOCAN 2012

CXCR4 and CCR7: Two eligible targets in targeted cancer therapy

Protein Data Bank (PDB): Database of threedimensional structural information of biological macromolecules

Computer-aided drug design methods

Structural insight into the binding pattern and interaction mechanism of chemotherapeutic agents with Sorcin by docking and molecular dynamic simulation

Molecular Docking and Simulation Studies of Antidiabetic Agents Devised from Hypoglycemic Polypeptide-P of Momordica charantia

Structure-based virtual screening and molecular dynamics of phytochemicals derived from Saudi medicinal plants to identify potential COVID-19 therapeutics

Lead-and drug-like compounds: The rule-of-five revolution

Molecular mechanisms of the action of Arctigenin in cancer

Determination of ensemble-average pairwise root mean-square deviation from experimental B-factors

Molecular dynamics simulation studies of novel β-lactamase inhibitor

Radius of gyration as an indicator of protein structure compactness

The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities

Binding free energy based analysis of arsenic (+3 oxidation state) methyltransferase with S-adenosylmethionine

Rapid decomposition and visualisation of protein-ligand binding free energies by residue and by water

Anticancer screening of medicinal plant phytochemicals against Cyclin-Dependent Kinase-2 (CDK2): An in-silico approach

Comparative Studies of the Dynamics Effects of BAY60-2770 and BAY58-2667 Binding with Human and Bacterial H-NOX Domains

Screening of medicinal plant phytochemicals as natural antagonists of p53-MDM2 interaction to reactivate p53 functioning

LigandScout: 3-D pharmacophores derived from protein-bound ligands and their use as virtual screening filters

Small-molecule library screening by docking with PyRx

Mutational Landscape of Pirin and Elucidation of the Impact of Most Detrimental Missense Variants That Accelerate the Breast Cancer Pathways: A

A General Amber Force Field for the accurate molecular dynamics simulation of phospholipid

The FF14SB force field

Langevin stabilization of molecular dynamics

Accuracy and efficiency of the particle mesh Ewald method

A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations

Software for processing and analysis of molecular dynamics trajectory data

VMD: Visual molecular dynamics

Radial Distribution Functions of Some Structures of the Polypeptide Chain

py: An efficient program for end-state free energy calculations

A water-swap reaction coordinate for the calculation of absolute protein-ligand binding free energies

Toward novel inhibitors against KdsB: A highly specific and selective broad-spectrum bacterial enzyme

Predicting small-molecule pharmacokinetic and toxicity properties using graph-based signatures

SwissADME: A free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules

The researcher would like to thank the Deanship of Scientific Research, Qassim University for funding the publication of this project.

The author declares no conflict of interest.Sample Availability: Not applicable.