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

Abstract

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:(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.

[Figure omitted. See PDF.]

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.

[Figure omitted. See PDF.]

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

[Figure omitted. See PDF.]

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.

[Figure omitted. See PDF.]

Molecular docking

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).

[Figure omitted. See PDF.]

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

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.

[Figure omitted. See PDF.]

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.

[Figure omitted. See PDF.]

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.

Molecular mechanics generalized born and surface area (MMGBSA) calculations

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).

[Figure omitted. See PDF.]

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).

[Figure omitted. See PDF.]

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).

[Figure omitted. See PDF.]

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.

Supporting information

S1 File.

https://doi.org/10.1371/journal.pone.0271401.s001

(DOCX)

S1 Graphical abstract.

https://doi.org/10.1371/journal.pone.0271401.s002

(TIF)

Acknowledgments

Authors are also grateful to Department of Science and Technology–Philippine Council for Agriculture, Aquatic, and Natural Resources Research and Development (DOST-PCAARRD); DST-FIST support by Govt. of India to the Department of Microbiology, Ramakrishna Mission Vivekananda Centenary College, West Bengal, India and Department of Botany, Gauhati University, Assam, India. This work was funded by Princess Nourah Bint Abdulrahman University Researchers Supporting Project number (PNURSP2022R62), Princess Abdulrahman University, Riyadh, Saudi Arabia, Researchers Supporting Project number (RSP-2021/26), King Saud University, Riyadh, Saudi Arabia.

Citation: Castrosanto MA, Mukerjee N, Ramos AR, Maitra S, Manuben JJP, Das P, et al. (2022) Abetting host immune response by inhibiting rhipicephalus sanguineus Evasin-1: An in silico approach. PLoS ONE 17(9): e0271401. https://doi.org/10.1371/journal.pone.0271401

About the Authors:

Melvin A. Castrosanto

Roles: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Writing – original draft, Writing – review & editing

Affiliation: Institute of Mathematical Sciences and Physics, College of Arts and Sciences, University of the Philippines Los Baños, Los Baños, Laguna, Philippines

https://orcid.org/0000-0003-2492-2379

Nobendu Mukerjee

Contributed equally to this work with: Nobendu Mukerjee, Swastika Maitra, Padmashree Das

Roles: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Writing – original draft, Writing – review & editing

E-mail: [email protected] (AG); [email protected] (AA); [email protected] (NM); [email protected] (SA)

Affiliations Department of Microbiology, Ramakrishna Mission Vivekananda Centenary College, West Bengal, Kolkata, India, Department of Health Sciences, Novel Global Community Educational Foundation, Hebersham, New South Wales, Australia

https://orcid.org/0000-0002-7129-7003

Ana Rose Ramos

Roles: Writing – original draft, Writing – review & editing

Affiliation: Institute of Chemistry, College of Arts and Sciences, University of the Philippines Los Baños, Los Baños, Laguna, Philippines

https://orcid.org/0000-0001-7094-8911

Swastika Maitra

Contributed equally to this work with: Nobendu Mukerjee, Swastika Maitra, Padmashree Das

Roles: Methodology, Writing – original draft, Writing – review & editing

Affiliation: Adamas University, West Bengal, Kolkata, India

John Julius P. Manuben

Roles: Data curation, Writing – review & editing

Affiliation: National Crop Protection Center, College of Agriculture and Food Science, University of the Philippines Los Baños, Los Baños, Laguna, Philippines

https://orcid.org/0000-0002-9260-7971

Padmashree Das

Contributed equally to this work with: Nobendu Mukerjee, Swastika Maitra, Padmashree Das

Roles: Formal analysis, Writing – original draft, Writing – review & editing

Affiliation: Central Silk Board, Regional Office, Khanapara, Guwahati, Assam, India

Sumira Malik

Roles: Methodology

Affiliation: Department of Biotechnology, Amity University Jharkhand Ranch, Jharkhand, India

https://orcid.org/0000-0001-5077-1493

Mohammad Mehedi Hasan

Roles: Resources, Writing – review & editing

Affiliation: Department of Biochemistry and Molecular Biology, Faculty of Life Science, Mawlana Bhashani Science and Technology University, Tangail, Bangladesh

Athanasios Alexiou

Roles: Investigation, Supervision, Writing – review & editing

E-mail: [email protected] (AG); [email protected] (AA); [email protected] (NM); [email protected] (SA)

Affiliation: Department of Science and Engineering, Novel Global Community Educational Foundation, Australia & AFNP Med, Hebersham, New South Wales, Australia

Abhijit Dey

Roles: Supervision

Affiliation: Department of Life Sciences, Presidency University, Kolkata, West Bengal, India

Mohammad Amjad Kamal

Roles: Supervision

Affiliations Institutes for Systems Genetics, Frontiers Science Center for Disease-Related Molecular Network, West China Hospital, Sichuan University, Chengdu, Sichuan, China, King Fahd Medical Research Center, King Abdulaziz University, Jeddah, Saudi Arabia, Department of Pharmacy, Faculty of Allied Health Sciences, Daffodil International University, Dhaka, Bangladesh, Enzymoics, 7 Peterlee Place, Habersham, New South Wales, Australia, Novel Global Community Educational Foundation, Habersham, New South Wales, Australia

Nada H. Aljarba

Roles: Funding acquisition, Validation, Writing – review & editing

Affiliation: Department of Biology, College of Science, Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia

https://orcid.org/0000-0001-9224-8844

Saad Alkahtani

Roles: Funding acquisition, Validation, Writing – review & editing

E-mail: [email protected] (AG); [email protected] (AA); [email protected] (NM); [email protected] (SA)

Affiliation: Department of Zoology, College of Science, King Saud University, Riyadh, Saudi Arabia

Arabinda Ghosh

Roles: Conceptualization, Formal analysis, Investigation, Project administration, Supervision, Validation, Writing – original draft, Writing – review & editing

E-mail: [email protected] (AG); [email protected] (AA); [email protected] (NM); [email protected] (SA)

Affiliation: Microbiology Division, Department of Botany, Gauhati University, Guwahati, Assam, India

https://orcid.org/0000-0002-3891-5949

References

1. Pegram R, Keirans J, Clifford C, Walker J. Clarification of the Rhipicephalus sanguineus group (Acari, Ixodoidea, Ixodidae). II. R. sanguineus (Latreille, 1806) and related species. Systematic Parasitology. 1987; 10:27–44.

2. Walker JB, Keirans JE, Horak IG. The Genus Rhipicephalus (Acari, Ixodidae). The Genus Rhipicephalus (Acari, Ixodidae). Cambridge University Press; 2000.

3. Dantas-Torres F. The brown dog tick, Rhipicephalus sanguineus (Latreille, 1806) (Acari: Ixodidae): From taxonomy to control. Vol. 152, Veterinary Parasitology. 2008. p. 173–85. pmid:18280045

4. Koch HG. Oviposition of the Brown Dog Tick (Acari: Ixodidae) in the Laboratory 1. 1982.

5. Troughton DR, Levin ML. Life Cycles of Seven Ixodid Tick Species (Acari: Ixodidae) Under Standardized Laboratory Conditions. Vol. 44, J. Med. Entomol. 2007. pmid:17915502

6. Szabó MPJ, Szabó S, Bechara GH. Sequential histopathology at the Rhipicephalus sanguineus tick feeding site on dogs and guinea pigs. 1999.

7. Mans BJ, Neitz AWH. Adaptation of ticks to a blood-feeding environment: Evolution from a functional perspective. Vol. 34, Insect Biochemistry and Molecular Biology. Elsevier Ltd; 2004. p. 1–17. pmid:14723893

8. Parola P, Raoult D. Ticks and Tickborne Bacterial Diseases in Humans: An Emerging Infectious Threat [Internet]. 2001. http://cid.oxfordjournals.org/

9. Ribeiro JMC, Makoul GT, Levine J, Robinson DR, Spielman A. Antihemostatic, antiinflammatory, and immunosuppressive properties of the saliva of tick, ixodes dammini a [Internet]. 1985. http://rupress.org/jem/article-pdf/161/2/332/1393852/332.pdf

10. Kazimírová M, Štibrániová I. Tick salivary compounds: Their role in modulation of host defences and pathogen transmission. Vol. 4, Frontiers in Cellular and Infection Microbiology. 2013. pmid:23971008

11. Bowman A, Coons L, Needham G, Sauer J. Tick saliva: recent advances and implications for vector competence. Vol. 11, Medical and Veterinary Entomology. 1997. pmid:9330260

12. Valenzuela JG. Exploring tick saliva: From biochemistry to “sialomes” and functional genomics. Vol. 129, Parasitology. 2004. pmid:15938506

13. Murphy P. Chemokines and chemokine receptors. In: Nihon rinsho Japanese journal of clinical medicine. 2012. p. 212–7.

14. Zlotnik A, Yoshie O. The Chemokine Superfamily Revisited. Vol. 36, Immunity. Cell Press; 2012. p. 705–16. pmid:22633458

15. Déruaz M, Frauenschuh A, Alessandri AL, Dias JM, Coelho FM, Russo RC, et al. Ticks produce highly selective chemokine binding proteins with antiinflammatory activity. Journal of Experimental Medicine. 2008 Sep 1;205(9):2019–31. pmid:18678732

16. Frauenschuh A, Power CA, Déruaz M, Ferreira BR, Silva JS, Teixeira MM, et al. Molecular cloning and characterization of a highly selective chemokine-binding protein from the tick Rhipicephalus sanguineus. Journal of Biological Chemistry. 2007 Sep 14;282(37):27250–8. pmid:17640866

17. Piao L, Chen Z, Li Q, Liu R, Song W, Kong R, et al. Molecular Dynamics Simulations of Wild Type and Mutants of SAPAP in Complexed with Shank3. International Journal of Molecular Sciences. 2019 Jan 8;20(1):224. pmid:30626119

18. Ramamoorthi N, Narasimhan S, Pal U, Bao F, Yang XF, Fish D, et al. The Lyme disease agent exploits a tick protein to infect the mammalian host. Nature. 2005 Jul 28;436(7050):573–7. pmid:16049492

19. Brossard M, Wikel SK. Tick immunobiology. In: Ticks: Biology, Disease and Control. Cambridge University Press; 2008. p. 186–204.

20. Nuttall PA, Labuda M. Saliva-assisted transmission of tick-borne pathogens. In: Ticks: Biology, Disease and Control. Cambridge University Press; 2008. p. 205–19.

21. Kumari M, Subbarao N. Virtual screening to identify novel potential inhibitors for Glutamine synthetase of Mycobacterium tuberculosis. Journal of Biomolecular Structure and Dynamics. 2020 Nov 21;38(17):5062–80. pmid:31755360

22. Marinho EM, Batista de Andrade Neto J, Silva J, Rocha da Silva C, Cavalcanti BC, Marinho ES, et al. Virtual screening based on molecular docking of possible inhibitors of Covid-19 main protease. Microbial Pathogenesis. 2020 Nov 1;148. pmid:32619669

23. Kapusta K, Kar S, Collins JT, Franklin LM, Kolodziejczyk W, Leszczynski J, et al. Protein reliability analysis and virtual screening of natural inhibitors for SARS-CoV-2 main protease (Mpro) through docking, molecular mechanic & dynamic, and ADMET profiling. Journal of Biomolecular Structure and Dynamics. 2021;39(17):6810–27.

24. Opo FADM, Rahman MM, Ahammad F, Ahmed I, Bhuiyan MA, Asiri AM. Structure based pharmacophore modeling, virtual screening, molecular docking and ADMET approaches for identification of natural anti-cancer agents targeting XIAP protein. Scientific Reports. 2021 Dec 1;11(1).

25. Dias JM, Losberger C, Déruaz M, Power CA, Proudfoot AEI, Shaw JP. Structural basis of chemokine sequestration by a tick chemokine binding protein: The crystal structure of the complex between Evasin-1 and CCL3. PLoS ONE. 2009;4(12). pmid:20041127

26. Schrodinger Inc. Maestro v12.5. Portland, OR; 2020.

27. Madhavi Sastry G, Adzhigirey M, Day T, Annabhimoju R, Sherman W. Protein and ligand preparation: Parameters, protocols, and influence on virtual screening enrichments. Journal of Computer-Aided Molecular Design. 2013 Mar;27(3):221–34. pmid:23579614

28. Xu Y, Wang S, Hu Q, Gao S, Ma X, Zhang W, et al. CavityPlus: A web server for protein cavity detection with pharmacophore modelling, allosteric site identification and covalent ligand binding ability prediction. Nucleic Acids Research. 2018 Jul 2;46(W1):W374–9. pmid:29750256

29. Sunseri J, Koes DR. Pharmit: interactive exploration of chemical space. Nucleic acids research. 2016 Jul 8;44(W1):W442–8. pmid:27095195

30. Sterling T, Irwin JJ. ZINC 15—Ligand Discovery for Everyone. Journal of Chemical Information and Modeling. 2015 Nov 23;55(11):2324–37. pmid:26479676

31. Kulkarni VM, Bhansali S. Pharmacophore generation, atom-based 3D-QSAR, docking, and virtual screening studies of p38-α mitogen activated protein kinase inhibitors: pyridopyridazin-6-ones (part 2). Research and Reports in Medicinal Chemistry. 2013 Dec;1.

32. O’Boyle NM, Banck M, James CA, Morley C, Vandermeersch T, Hutchison GR. Open Babel: An open chemical toolbox. Journal of Cheminformatics. 2011 Dec 7;3(1):33. pmid:21982300

33. Dallakyan S, Olson AJ. Small-molecule library screening by docking with PyRx. Methods in Molecular Biology. 2015;1263:243–50. pmid:25618350

34. Sander T, Freyss J, von Korff M, Reich JR, Rufener C. OSIRIS, an entirely in-house developed drug discovery informatics system. Journal of Chemical Information and Modeling. 2009 Feb 23;49(2):232–46. pmid:19434825

35. Bowers KJ, Chow DE, Xu H, Dror RO, Eastwood MP, Gregersen BA, et al. Scalable Algorithms for Molecular Dynamics Simulations on Commodity Clusters. In: ACM/IEEE SC 2006 Conference (SC’06). IEEE; 2006. p. 43–43.

36. Chow E, Rendleman CA, Bowers KJ, Dror RO, Hughes DH, Gullingsrud J, et al. Desmond Performance on a Cluster of Multicore Processors Hardware and Operating Environment Benchmark Systems and Simulation Parameters. 2008.

37. Shivakumar D, Williams J, Wu Y, Damm W, Shelley J, Sherman W. Prediction of Absolute Solvation Free Energies using Molecular Dynamics Free Energy Perturbation and the OPLS Force Field. Journal of Chemical Theory and Computation. 2010 May 11;6(5):1509–19. pmid:26615687

38. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of simple potential functions for simulating liquid water. The Journal of Chemical Physics. 1983 Jul 15;79(2):926–35.

39. Martyna GJ, Tobias DJ, Klein ML. Constant pressure molecular dynamics algorithms. The Journal of Chemical Physics. 1994 Sep;101(5):4177–89.

40. Martyna GJ, Klein ML, Tuckerman M. Nosé–Hoover chains: The canonical ensemble via continuous dynamics. The Journal of Chemical Physics. 1992 Aug 15;97(4):2635–43.

41. Toukmaji AY, Board JA. Ewald summation techniques in perspective: a survey. Computer Physics Communications. 1996 Jun;95(2–3):73–92.

42. Kagami LP, das Neves GM, Timmers LFSM, Caceres RA, Eifler-Lima VL. Geo-Measures: A PyMOL plugin for protein structure ensembles analysis. Computational Biology and Chemistry. 2020 Aug;87:107322. pmid:32604028