Abetting host immune response by inhibiting rhipicephalus sanguineus Evasin-1: An in silico approach.

The brown dog tick (Rhipicephalus sanguineus) is the most prevalent tick in the world and a well-recognized vector of many pathogens affecting dogs and occasionally humans. Pathogens exploit tick salivary molecules for their survival and multiplication in the vector and transmission to and establishment in the hosts. Tick saliva contains various non-proteinaceous substances and secreted proteins that are differentially produced during feeding and comprise of inhibitors of blood congealing and platelet aggregation, vasodilatory and immunomodulatory substances, and compounds preventing itch and pain. One of these proteins is Evasin-1, which has a high binding affinity to certain types of chemokines. The binding of Evasin-1 to chemokines prevents the detection and immune response of the host to R. sanguineus, which may result in the successful transmission of pathogens. In this study, we screened potential Evasin-1 inhibitor based on the pharmacophore model derived from the binding site residues. Hit ligands were further screened via molecular docking and virtual ADMET prediction, which resulted in ZINC8856727 as the top ligand (binding affinity: -9.1 kcal/mol). Molecular dynamics simulation studies, coupled with MM-GBSA calculations and principal component analysis revealed that ZINC8856727 plays a vital role in the stability of Evasin-1. We recommend continuing the study by developing a formulation that serves as a potential medicine aid immune response during R. sanguineus infestation.

Introduction

The brown dog tick, Rhipicephalus sanguineus, of the Ixodidae family is a cosmopolitan species found widely distributed around the world [1]. The tick parasitizes animals such as dogs, cats, rabbits, bovines, reptiles, and in some cases humans [2]. It is the most studied tick considering its veterinary and public health importance, that is, a vector of disease agents with notable zoonotic concerns such as Coxiella burnetii, Ehrlichia canis, Babesia canis, Rickettsia conorii, and Rickettsia rickettsia [3]. The feeding behavior of R. sanguineus differs at the various stages of its life cycle, with larvae feeding for 2 days while adult females last for several weeks, and adult males can take multiple blood meals [4, 5]. R. sanguineus can have different hosts during its life cycle because of its cryptic behavior of feeding and then dropping off the host for multiple times. Their host seeking activity which appears on active stages enable R. sanguineus to easily find a host due to its increase movement [3].

Once attached to its host, the R. sanguineus uses its mouthpart in penetrating the host’s skin and then inserting the hypostome and mouthpart in the epidermis [6]. While attached, it releases a cement-like fluid that creates a cone on the epidermis up the outermost layer of the epidermis [6]. During feeding, the mouthpart of the R. sanguineus rips capillaries and small blood vessels when probing for blood. This creates hemorrhage, creating a blood pool, where R. sanguineus feeds from the blood and other fluids [7]. There is an alternating phases of blood sucking and salivating with regurgitation occurring normally in the process [8]. The phase of forceful salivation and regurgitation is of great significance for the transmission of pathogens. The R. sanguineus’ saliva components allow the completion of blood feeding through the repression of the host immune and inflammatory response, letting the tick attached for an extended number of time sometimes as long as 2–3 weeks [9, 10]. The saliva contains large number of bioactive molecules including salivary proteins exhibiting anticoagulation, antiplatelet, vasodilatory, anti-inflammatory, and immunomodulatory activities [11, 12]. The discovery that salivary gland extracts from numerous ixodid tick species could neutralize the chemokine CXCL8 and CCL2, CCL3, CCL5, and CCL11 provided the first indications of anti-chemokine activity in ticks.

One strategy utilized by R. sanguineus to prevent rejection from its host is the inhibition of immune cell recruitment which is accomplished by a class of chemokine-binding proteins (CKBPs) known as Evasin. During blood feeding, the host’s proper physiological response to injury or tissue damage is through inflammation. Inflamed cells recruit leukocytes to the area which results in inhibition of pathogens and repair the damaged tissues. The recruitment of leukocytes is regulated by chemokines which are secreted in the site of injury which then activates the chemokine receptors [13, 14]. Chemokines are divided into two major families (CC and CXC) and two minor families (CX3CL and XCL) depending on the distance between the two N-terminal Cys residues (Zlotnik and Yoshie, 2000; Zlotnik and Yoshie, 2012). The two classes of Evasin found in R. sanguineus’ saliva binds to the host chemokines which results to inhibition of chemokine receptor activation. The first class of Evasin (Evasin-1 and Evasin-4) consists of proteins with eight conserved Cys residues and completely binds to the CC chemokines. The second class binds to CXC chemokines and consists of six Cys residues (Evasin-3) [15, 16]. Moreover, selectivity profiles of the different Evasins showed that Evasin-1 has highest specificity to several CC chemokines–CCL3, CCL3L1, CCL4, CCL4L1, CCL14, and CCL18 compared to Evasin-4 which binds to around 20 CC chemokines [17]. Their binding profiles also showed that Evasin- 1 and Evasin-3 impede neutrophil migration, which is important in the early stages of the immunological response in mice although expression profile of chemokine receptors on leukocytes in the dog are still unknown. With the parasite successful blood feeding, pathogen transmission via tick saliva is not a simple routine process but it was noted that pathogens were able to exploit the bioactive components of the tick’s saliva for their survival, multiplication, transmission, and establishment to the host [11, 15, 18, 19].

In view of the role of Evasin in the possible transmission of pathogen causing diseases, structure-based virtual screening centered on the investigation of potential inhibitors of Evasin-1 that weaken its binding affinity to chemokines is a good strategy in managing diseases through the salivary-assisted transmission in ticks. Computational approaches such as the molecular-based virtual screening studies are widely used in the screening, identification, and analysis of possible inhibitory compounds [20-23]. Henceforth, this study focused on the use of virtual screenings based on the pharmacophore model derived from the binding site residues, ADMET, and molecular dynamics simulation to identify possible compounds with inhibitory activity against Evasin-1.

Materials and methods

Target protein preparation

We retrieved the crystal structure of Rhipicephalus sanguineus’ Evasin-1 [24] from the Protein Data Bank (https://www.rcsb.org/structure/3FPR). The Protein Preparation Wizard (PrepWizard) in Maestro [25] was used to prepare Evasin-1 in this study. In PrepWizard, water molecules were removed, bond orders were assigned, and hydrogens were added. Then the protein was refined by optimizing the hydrogen bonds. Lastly, minimization was done via two steps–(1) H-only optimization followed by (2) all-atom minimization with termination based on convergence or reaching a heavy atom root mean square deviation (RMSD) of 0.30 Å using the OPLS3e force field [26]. The Ramachandran plot along with other structure quality information were generated from the SWISS-MODEL structure assessment tool (https://swissmodel.expasy.org/assess) before and after the optimization to check the quality of the protein.

Receptor-based pharmacophore modelling and virtual screening

The probable ligand-binding site of Evasin I (PDB ID: 3FPR) was predicted using the Cavity module of the CavityPlus web server [27]. Then using the cavity with the highest druggability score, receptor-based pharmacophore model was generated through the CavPharmer Module. We noted the important features in the predicted binding site and used the gained information (pharmacophore and coordinates) as input in Pharmit [28] to virtually screen about 120 million conformers of 13 million molecules from the ZINC Purchasable database [29]. To limit the number of hit compounds, exclusive shape was set to receptor with 1 Å tolerance, hit reduction was set to 1 hit per configuration per molecule and a maximum of 50 total hits having the lowest RMSD to the pharmacophore model. We chose the molecular weight range of 180 g/mol to 500 g/mol, and rotational bonds of 5 to 10 as described [30].

Molecular docking

The top 50 ligand hits with were uploaded to the Open Babel [31] tab of the PyRx software [32]. Ligand minimization was done using the mmff94 force field and 1000 steps of conjugate gradient optimization algorithm. Then, file format was converted from sdf to pdbqt simultaneously uploading it as an AutoDock ligand. The optimized and minimized 3FPR was also uploaded as AutoDock macromolecule. Molecular docking was conducted using the Vina Wizard with the grid box of x: 10.7391 Å, y: 11.9125 Å, z: 20.2942 Å and center of x: 33.2765, y: 36.2609, z: 11.6079. Binding affinities (in kcal/mol) of the ligands were recorded and compared. We did virtual ADMET prediction of the top 10 ligand hits with the most favorable binding affinity values using Osiris Property Explorer software [33]. The top ligand that passed the toxicity prediction was further used in molecular dynamics simulation. The protein-ligand interaction diagram of the top ligand was generated in 2D and 3D.

Molecular dynamics (MD) simulation

The MD simulations studies were carried in triplicate on docked complex for Evasin-1 with ZINC8856727 using the Desmond 2020.1 from Schrödinger, LLC. The triplicate samplings were made using same parameters for each MD run in order to obtain reproducibility of the results. The OPLS-2005 force field [34-36] and explicit solvent model with the SPC water molecules were used in this system [37]. Na+ ions were added to neutralize the charge 0.15 M, NaCl solutions were added to the system to simulate the physiological environment. Initially, the system was equilibrated using an NVT ensemble for 100 ns to retrain over the 3FPR-ZINC8856727 complex. Following the previous step, a short run of equilibration and minimization was carried out using an NPT ensemble for 12 ns. The NPT ensemble was set up using the Nose-Hoover chain coupling scheme [38] with the temperature at 37 ºC, the relaxation time of 1.0 ps, and pressure 1 bar maintained in all the simulations. A time step of 2 fs was used. Canonical Ensemble (NVT) collects systems whose thermodynamic state is characterized by a fixed number of atoms, N, a fixed volume, V, and a fixed temperature, T. Whereas, the NPT ensemble describes systems in contact with a thermostat at temperature T and a bariostat at pressure P. Therefore, equilibration at constant volume and pressure along with constant number and temperature will maintain a conformational the sampling from a stable system. The Martyna-Tuckerman–Klein chain coupling scheme [39] barostat method was used for pressure control with a relaxation time of 2 ps. The particle mesh Ewald method [40] was used for calculating long-range electrostatic interactions, and the radius for the coulomb interactions were fixed at 9Å. RESPA integrator was used for a time step of 2 fs for each trajectory to calculate the bonded forces. The root means square deviation (RMSD), radius of gyration (Rg), root mean square fluctuation (RMSF) and number of hydrogen (H-bonds) were calculated to monitor the stability of the MD simulations. The free energy landscape of protein folding on ZINC8856727-bound complex was measured using Geo_measures v 0.8 [41]. Geo_measures include a powerful library of g_sham and form the MD trajectory against RMSD and radius of gyration (Rg) energy profile of folding recorded in a 3D plot using matplotlib python package.

Molecular mechanics generalized born and surface area (MMGBSA) calculations

During MD simulations of Evasin-1 complexed with ZINC8856727, the binding free energy (Gbind) of docked complexes was calculated using the premier molecular mechanics generalized Born surface area (MM-GBSA) module (Schrodinger suite, LLC, New York, NY, 2017–4). The binding free energy was calculated using the OPLS 2005 force field, VSGB solvent model, and rotamer search methods [42]. After the MD run, 10 ns intervals were used to choose the MD trajectories frames. The total free energy binding was calculated using Eq 1: Where, ΔGbind = binding free energy, Gcomplex = free energy of the complex, Gprotein = free energy of the target protein, and Gligand = free energy of the ligand. The MMGBSA outcome trajectories were analyzed further for post dynamics structure modifications.

Principal component analysis (PCA)

During a 100 ns simulation of Evasin-1 complexed with ZINC8856727, PCA analysis was used to recover the global movements of the trajectories. To calculate the PCA, a covariance matrix was created as stated. For conformational analysis of the ZINC8856727 in bound complex, 10 alternative conformational modes of the main component as movements of trajectories were calculated, and a comparison of the mode PC1, PC2, PC3 and the last modes PC9 and PC10 were investigated to understand the convergence of trajectories. The MD trajectory versus PC2 energy profiles of folding was recorded in a 3D plot using the matplotlib python package using Geo measures, which includes a comprehensive library of g_sham.

Results and discussion

Protein preparation

The Evasin-1 that was retrieved from the Protein Data Bank contains water molecules that we removed before we subject it to the structure assessment tool of the SWISS-MODEL web server. The MolProbity score, clash score, and Ramachandran residue outlier were all higher in the unprepared protein compared to the optimized and minimized one (Table 1). Clashes involved in the unprepared protein are A: THR25 and A: ARG45, A: ASN88 and D: LYS20, and A: GLY7 and A: ASP8. These clashes were all removed upon preparing the protein. Similarly, the only Ramachandran outlier, A:ASP8, was also corrected during the protein preparation. MolProbity score improved from 1.33 to 1.01, an indication of the overall improvement of the protein structure’s quality. Superimposition of the unprepared and prepared Evasin-1 is shown in Fig 1 of the S1 File, together with their respective Ramachandran plots.

Table 1

Protein structure quality information as given by the SWISS-MODEL structure assessment tool.

Protein	Ramachandran Outliers	Clash Score	MolProbity Score
Before preparation	0.63% (A:ASP8)	4.21	1.33
After preparation	0.00%	0.42	1.01

Receptor-based pharmacophore modelling

Since no co-crystallized ligand is present in the Evasin-1’s deposited crystal structure, we decided to predict its most probable binding site and determine the binding site residues using the Cavity module of the CavityPlus web server. There was a total of seven cavities predicted (Table 2), but we chose the cavity with the highest drug score and strongest druggability among the results (Fig 1A). The predicted binding site lies in between the two chains of Evasin-1 with chain A having more involved residues. Binding site residues in the chain A includes residues 9–16, 22–26, 28–29, 31–33, 36–46, and 70. On the other hand, binding site residues belonging to chain B are 47–50, 63–67, 75, 77, and 79–83. Using the CavPharmer module of the CavityPlus web server, we gained information on important pharmacophores and their corresponding coordinates (Fig 1B). Selected pharmacophores are classified into positive electrostatic center (blue), hydrogen bond donor center (white), and hydrophobic center (green). The complete set of pharmacophores along with the coordinates were presented at Table 1 of the S1 File.

Table 2

Binding site prediction of the Cavity module.

Cavity No.	Predicted Ave. pKd	Drug Score	Druggability
1	6.81	831.00	Strong
2	6.16	-785.00	Weak
3	5.76	-584.00	Weak
4	5.60	-958.00	Weak
5	5.20	-1046.00	Weak
6	4.94	-1428.00	Weak
7	4.50	-1605.00	Weak

Fig 1

The 3D-surface of Evasin-1’s predicted ligand binding site (left) along with the pharmacophores used in the screening of ligands in Pharmit (right).

Positive electrostatic center is colored blue, hydrophobic centers are green, and hydrogen donor center is white.

Virtual screening was done on ZINC database using the pharmacophores presented in Fig 1B. Molecular weight was restricted from 180 g/mol to 500 g/mol as molecules that are too light tends to be excreted easily, while molecules that are too heavy may be hard to absorbed by the system. The surface of the receptor was set as an excluded space to make sure that the resulting ligand hits will not cause clashes with the side chains of the proteins during docking. Most of the top ligand hits contains nitro groups as the positive electrostatic center, an aromatic ring for the hydrophobic centers, and a nitrogen atom for the hydrogen donor center. Fig 2 shows the superimposition of the top 50 ligand hits with the pharmacophore used in the screening. The top 50 ligand hits according to the RMSD values were retrieved and docked to Evasin-1.

Fig 2

Superimposition of the top 50 ligand hits generated from the ZINC database.

The top 50 ligand hits were docked in the predicted binding site of Evasin-1 to calculate their respective binding energies. The resulting values ranges from -6.5 kcal/mol to -9.3 kcal/mol, having significant binding affinities (Table 3 in S1 File). Exactly 6 have binding energies of -9.0 to -9.3 kcal/mol, 18 have -8.0 to -8.9 kcal/mol, 22 scored -7.0 to -7.9 kcal/mol, and only 4 scored -6.5 to -6.9 kcal/mol. The ligand ZINC9594956 resulted in the most favorable interaction with Evasin-1 based on the binding affinity value. However, virtually evaluated top 10 ligands (Table 2, S1 File), only 3 of them passed the non-toxic, non-irritant and have no effect on reproductive efficiency criterion (ZINC8856727, ZINC13552533, and ZINC585290522). We then proceeded with ZINC8856727 as the ligand of interest for the rest of the study due to its favorable binding affinity and non-toxicity. The 3D structure of Evasin-1-ZINC8856727 complex, and the 2D and 3D interaction diagram (Fig 3).

Fig 3

The 3D rendition, and the 2D and 3D interaction diagram of ZINC8856727 with Evasin-1.

The high affinity of ZINC8856727 to Evasin-1 can be attributed to its various interactions. Two pi-anion interactions are depicted in the 2D diagram, which came from the carboxylic oxygen of D:ASP64. Upon generating the 3D interaction diagram, we have found an internal pi-anion interaction coming from the nitro group of the ligand. A T-shaped pi-pi interaction between D:TYR81 and one aromatic ring of the ligand also contributes to the affinity. Moreover, the alkyl group of D:PRO66 interacts with the pi-electrons of one aromatic ring. Lastly, weak but notable interactions arise from carbon-hydrogen bond of ligand’s nitro group to A:PRO41 and also the presence of hydrophobic alkyl-alkyl interaction.

Molecular dynamics (MD) simulation studies were carried out to determine the stability and convergence of ZINC8856727-bound Evasin-1 complex. Each simulation of 100 ns displayed stable conformation while comparing the root mean square deviation (RMSD) values. The root mean square deviation of apo-Evasin-1 was considerably higher as compared to ZINC8856727-bound Evasin-1, signify the bound state conforms into more stable structure. The Cα-backbone of Evasin-1 bound to ZINC8856727 exhibited a deviation of 0.2 Å (Table 3, Fig 4A; R1, R2, R3) RMSD plots are within the acceptable range signifying the stability of proteins in the ZINC8856727 bound state before and after simulation, and it also suggests that the stability might be due to significant binding of the ligand.

Table 3

The average root mean square deviation, radius of gyration, hydrogen bonding.

	R1 (Å)	R2 (Å)	R3 (Å)	Average value (Å)
RMSD	0.3	0.1	0.2	0.2
Radius of Gyration	0.8	0.7	0.5	0.67
Hydrogen Bonding	5	6	5	5

Fig 4

A. MD simulation trajectory analysis of Root Mean Square Divisions (RMSD) of ZINC8856727 bound with 3FPR, i.e. Evasin-1 100 ns time frame in triplicate displayed: R1 (replicate 1) RMSD plot of ZINC8856727 bound 3FPR (red) with control 3FPR (light green); R2 (replicate 2) RMSD plot of ZINC8856727 bound 3FPR (dark maroon) with control 3FPR (juniper green); R3 (replicate 3) RMSD plot of ZINC8856727 bound 3FPR (yellow) with control 3FPR (cyan) B. MD simulation trajectory analysis of Root Mean Square Fluctuations (RMSF) of ZINC8856727 bound with 3FPR, i.e. Evasin-1 100 ns time frame in triplicate displayed: R1 (replicate 1) RMSF plot of ZINC8856727 bound 3FPR (red) with control 3FPR (black); R2 (replicate 2) RMSF plot of ZINC8856727 bound 3FPR (munsell yellow) with control 3FPR (sap green); R3 (replicate 3) RMSF plot of ZINC8856727 bound 3FPR (munsell yellow) with control 3FPR (cyan blue); C. MD simulation trajectory analysis of Radius of gyration (Rg) of ZINC8856727 bound with 3FPR, i.e. Evasin-1 100 ns time frame in triplicate displayed: R1 (replicate 1) Rg plot of ZINC8856727 bound 3FPR (red) with control 3FPR (light green); R2 (replicate 2) Rg plot of ZINC8856727 bound 3FPR (porcelain) with control 3FPR (dark green); R3 (replicate 3) Rg plot of ZINC8856727 bound 3FPR (black) with control 3FPR (Prussian blue) D. MD simulation trajectory analysis of Hydrogen Bonding (H-Bonds) of ZINC0 = 8856727 bound with 3FPR, i.e. Evasin-1 100 ns time frame in triplicate displayed: R1 (replicate 1) H-Bond plot of ZINC8856727 bound 3FPR (red); R2 (replicate 2) H-Bond plot of ZINC8856727 bound 3FPR (black); R3 (replicate 3) H-Bond plot of ZINC8856727 bound 3FPR (light green) E. Stepwise trajectory analysis for every 25 ns displaying the protein and ligand conformation during 100 ns of simulation of Evasin-1-ZINC8856727 complex. F Free Energy Landscape displaying the achievement of global minima (ΔG, kJ/mol) of 3FPR in presence of ZINC8856727 with respect to their RMSD (nm) and Radius of gyration (Rg, nm). G. MMGBSA trajectory (0 ns, before simulation and 100 ns, after simulation) exhibited conformational changes of ZINC8856727 upon binding with the protein Evasin 1 (PDB I.D. 3FPR). The arrows indicating the overall positional variation (movement and pose) of ZINC08856727 at the binding site cavity.

Radius of gyration is the measure of the compactness of the protein. The ZINC8856727-bound Evasin-1 protein displayed lowering of Radius of Gyration (Fig 4B; R1, R2, R3) as compared to the apo-Evasin-1. Looking at the overall quality of analysis from root mean square deviation (RMSD), and Rg, it can be concluded that ZINC8856727 bound to the protein targets posthumously in the binding cavities and plays a significant role in the stability of the proteins. The plots for root mean square fluctuations (RMSF) displayed a significant high RMSF in the apo-protein at few residues as compared to ligand-bounded Evasin-1 at the specific time function of 100 ns (Fig 4C; R1, R2, R3). Therefore, from RMSF plots, it can be suggested that the protein structures were stable during simulation in ZINC8856727-bound conformation.

The average hydrogen bonds formed between ZINC8856727 and Evasin-1 during the 100 ns simulation were also noted and recorded in Fig 4D. From 0 ns to 100 ns, an average of four hydrogen bonding was found throughout the simulation (Fig 4D; R1, R2, R3). Overall, four hydrogen bonds were formed throughout the simulation and confirmed from 2D ligand binding plot. The number of hydrogen bonds between Evasin-1 and ZINC8856727 has strengthened the binding, helping to make it more stable during the simulation.

The stepwise trajectory analysis for every 25 ns of simulation of Evasin-1 bound to ZINC8856727 displayed the positional alteration with reference to 0 ns structure (Fig 4E). It has been observed that the ligand, ZINC8856727 possessed a structural angular movement at the end frame to achieve its conformational stability and convergence.

The free energy landscape (FEL) of achieving global minima of Cα backbone atoms of protein with respect to the RMSD and radius of gyration (Rg) is displayed in Fig 4F. Evasin-1 bound to ZINC8856727 achieved the global minima (lowest free energy state) at 2.9 Å and Rg 16.17 Å. The FEL envisaged for deterministic behavior of Evasin-1 to lowest energy state owing to its high stability and best conformation at ZINC8856727 bound state. Therefore, FEL is the indicator of the protein folding to attain minimum energy state, and that aptly achieved due to ZINC8856727 bound state.

To assess the binding energy of ligands to protein molecules, the MMGBSA technique is commonly employed. The binding free energy of each Evasin-1-ZINC8856727 complex, as well as the impact of other non-bonded interactions energies, were estimated. With Evasin-1, the ligand ZINC8856727 has a binding energy of -51.69 kcal/mol. Non-bonded interactions like GbindCoulomb, GbindCovalent, GbindHbond, GbindLipo, GbindSolvGB, and GbindvdW govern Gbind. Across all types of interactions, the GbindvdW, GbindLipo, and GbindCoulomb energies contributed the most to the average binding energy. On the other side, the GbindSolvGB and Gbind Covalent energies contributed the least to the final average binding energies. Furthermore, the GbindHbond interaction values of ZINC8856727-Evasin-1 complexes demonstrated stable hydrogen bonds with amino acid residues. In all of the compounds, GbindSolvGB and GbindCovalent exhibited unfavorable energy contributions, and so opposed binding. Fig 4G (left panel) reveals that between pre-simulation (0 ns) and post-simulation (100 ns), ZINC8856727 in the binding pockets of Evasin-1 has undergone a large angular change in the pose (curved to straight). These conformational changes lead to better binding pocket acquisition and interaction with residues, which leads to enhanced stability and binding energy (see Table 4).

Table 4

Binding energy calculation of ZINC8856727 with Evasin-1 and non-bonded interaction energies from MM-GBSA trajectories.

Energies (kcal/mol)	3FPR
ΔG bind	-51.67 ± 7.60
ΔG bind Lipo	-14.75 ± 1.72
ΔG bind vdW	-42.41 ± 4.57
ΔG bind Coulomb	-25.56 ± 7.01
ΔG bind H bond	-2.90 ± 0.71
ΔG bind SolvGB	-29.73 ± 3.53
ΔG bind Covalent	-6.01 ± 1.67

Thus, the MM-GBSA calculations resulted from MD simulation trajectories were well justified with the binding energy obtained from docking results. Moreover, the last frame (100 ns) of MMGBSA displayed the positional change of the ZINC8856727 as compared to 0 ns trajectory signify the better binding pose for best fitting in the binding cavity of the protein (Fig 4G).

Therefore, it can be concluded that the ZINC8856727 molecule has good affinity for the major target Evasin-1.

Principle component analysis (PCA) and energy calculations

Principal component analysis (PCA) determines the relationship between statistically meaningful conformations (major global motions) sampled during the trajectory. PCA of the MD simulation trajectories for Evasin-1 bound to ZINC8856727 molecule were analyzed to interpret the randomized global motion of the atoms of amino acid residues. The internal coordinates mobility into three-dimensional space in the spatial time of 10 ns were recorded in a covariance matrix and rational motion of each trajectory are interpreted in the form of orthogonal sets or Eigen vectors. In the Evasin-1 trajectory, PCA indicates the statistically significant conformations. It is possible to identify the major motions within the trajectory as well as the critical motions required for conformational changes. In Evasin-1 bound to ZINC8856727, two different clusters along the PC1 and PC2 plane have been exhibited, indicating a non-periodic conformational shift (Fig 5A). On the other hand, the global motions are periodic because the groupings along the PC2 and PC3 planes do not totally cluster separately (Fig 5B). Moreover, high periodic global motion was observed along the PC9 and PC10 planes due to the grouping of trajectories in a single cluster at the center of the PCA plot (Fig 5C).

Fig 5

Principal component analysis of the Evasin-1-ZINC8856727 complex showing a stable configuration throughout the 100 ns run.

Centering of the trajectories in a single cluster indicates the periodic motion of MD trajectories due to stable conformational global motion. Therefore, PCA analysis suggested that the Eigen vectors of relative aggregated motion of the trajectories became better at higher mode PC10 into a converted global motion of the trajectories.

The energy profiles of the protein and Evasin-1-ZINC8856727 complex system were determined to display the stability of the entire system. In this regard, the total energy (Total Energy) of the Evasin-1-ZINC8856727 system shown to be very stable with an average total energy -91.00 kcal/mol (pale green). However, Coulombic interactions displayed to be merged over the total energy with an average energy -49.00 kcal/mol (cyan), contemplating as principal contributor to the stability of the Evasin-1-ZINC8856727 complex. In addition, van der Waals energy played an equally important role in the system stability and contributing an average energy -45.00 kcal/mol (red) (Fig 6).

Fig 6

Energy plot of protein Evasin-1 and ZINC8856727 complex system during the entire simulation event of 100 ns.

The total energy (dark green), van der Waal’s energy (cyan) and coulomb energy (red) of the entire system indicating the stability of the individual systems bound to ZINC8856727 molecule.

Conclusion

In this paper, receptor-based pharmacophore screening was done in finding candidate inhibitors of Evasin-1, an important salivary enzyme of the brown dog tick (Rhipicephalus sanguineus). This enzyme is being used by pathogens to their advantage in disease transmission. Evasin-1 binds to a certain chemokine, resulting in the suppression of immune response. The top 50 ligands from virtual screening were filtered using molecular docking and virtual ADMET prediction. It resulted in the discovery of ZINC8856727 as the top ligand. Molecular dynamics simulation followed by energy calculations and structural analyzes proves the stability it contributes to Evasin-1 through important interactions. Till date, this is the first study to use computational methodologies in search for Evasin-1 inhibitors. The results of our study can trigger the discovery of a future medication that can aid or abet the host immune response during Rhipicephalus sanguineus feeding.

(DOCX)

Click here for additional data file.

(TIF)

Click here for additional data file.