Immunoinformatics based design and prediction of proteome-wide killer cell epitopes of Leishmania donovani: Potential application in vaccine development

Abstract Despite several extensive and exhaustive efforts, search for potential therapy against leishmaniasis has not made much progress. In the present work, we have employed mining strategy to screen Leishmania donovani proteome for identification of promising vaccine candidate. We have screened 21 potential antigenic proteins from 7960 total protein of L. donovani, based on the presence of signal peptide, GPI anchor, antigenicity prediction and substractive proteomic approach. Secondly, we have also performed comprehensive immunogenic epitope prediction from the screened 21 proteins, using IEDB-AR tools. Out of the 21 antigenic proteins, we obtained 11 immunogenic epitopes from 9 proteins. The final results revealed that four predicted epitopes namely; YPAFAALVF, VAVAATVAY, AAAPTEAAL and MYPLVAVVF, have significantly better binding potential with respective alleles and could elicits immune responses. Docking analysis using PATCHDOCK server and molecular dynamic simulation using GROMACS revealed the potential of the sequences as immunogenic epitopes. In silico studies also suggested that the epitopes occupied almost same binding cleft with the respective alleles, when compared with the reference peptides. It is also suggested from the molecular dynamic simulation data that the peptides were intact in the pocket for longer periods of time. Our study was designed to select MHC class I restricted epitopes for the activation of CD8 T cells using immunoinformatics for the prediction of probable vaccine candidate against L. donovani parasites. Communicated by Ramaswamy H. Sarma


Introduction
Leishmaniasis is a spectrum of disease, prevalent in tropical and subtropical areas and causes several thousand deaths every year . It is transmitted by the bites of sand fly of genus Phlebotomus (Vannier-Santos et al., 2002). Nearly 20 species of Leishmania parasite are responsible for different pathologies which accounts for cutaneous lesions (cutaneous, diffuse cutaneous and mucocutaneous leishmaniasis) as well as fatal liver and bone marrow infection in visceral leishmaniasis (VL) (Chappuis et al., 2007). The countries which are mostly affected by leishmaniasis are Brazil, Afghanistan, Iran, Saudi Arabia, Syria, India and Bangladesh (Kashif et al., , 2019. In these countries, leishmaniasis caused major health problem with significant impact on socio-economic conditions. Several works were undertaken in order to find better chemotherapeutic agent or inhibitor against Leishmania parasites, but with little convincing progress (Kashif et al., 2019). The available chemotherapies have high cost, wide spread side effects, poor tolerability besides showing drug resistance (Hussain et al., 2014). Therefore, it is an urgent need to develop better treatment option apart from drug discovery. Several studies have been made in this direction and few vaccine formulations like live attenuated or killed parasite, DNA or protein have also been made and tested in animal models (Bindu & Rentala, 2004;Kedzierski et al., 2006). However, the success rate is low and hence, development of novel methods for new vaccine candidate is urgently needed (Dumonteil, 2007).
In this context, computation-based antigen prediction and immunoinformatics approach for epitope prediction were performed. Studies were conducted earlier to identify vaccine candidates using pathogen genome (Rappuoli, 2000). Integration of immunology with bioinformatics for the development of vaccine is called as vaccinomics (Poland et al., 2009). These methods have been applied previously to identify vaccine against several diseases like malaria (L opez et al., 2001), dengue (Ali et al., 2017) and influenza (Shahsavandi et al., 2015). According to the literature, effective vaccine would elicit T cell responses; therefore, it is important to design cytolytic T lymphocytes (CTL) epitopes for the vaccine development process. Several algorithms and tools were adopted to predict the T and B cell epitopes for peptide based vaccine development process (He et al., 2010). Thus, the methods of T cell epitope prediction for antigen identification appeared to be an effective strategy, especially for the Leishmania parasite having large genome/proteome (Herrera-Najera et al., 2009). The role played by CD8 þ T cells in controlling primary Leishmania infection in genetically resistant mice remains controversial. Leishmania specific CD8 þ T cells were reported to be observed in higher frequencies in resistant mice compared with susceptible mice (Milon et al., 1986;Titus et al., 1987). Depletion of CD8 þ T cells augmented the severity of footpad lesions in mice inoculated with a highdose of parasites; however, the mice were able to control the infection (Titus et al., 1987). Studies in b2-microglobulin or CD8 deficient mice using a high-dose footpad challenge showed that the mice were able to heal their footpad lesions (Huber et al., 1998;Wang et al., 1993). All these studies provided a sound basis for concluding that CD8 þ T cells have little or no role in the control of primary infection by L. major. Primed CD8 þ T cells derived from acute stage infection of C57BL/6 mice showed both pathology and immunity following transfer (Belkaid et al., 2002). Intradermal challenge with low dose L. major lead to generation of CD8 þ T cells which play a critical role in dermal pathogenesis and immunity to primary infection (Belkaid et al., 2002). Although the role of CD8 cell mediated cytotoxicity in disease course has not been elucidated, elevated numbers of CD8 þ T cells have been observed in blood and lesions of patients infected with L. major and L. Mexicana, mediated via the production of IFN c (Ruiz & Becker, 2007). The numbers of CD8 þ T cells were significantly reduced in lesions of diffuse cutaneous leishmaniasis (DCL) patients infected with L. mexicana, as compared to LCL patients (Salaiza-Suazo, Volkow et al. 1999). CD8 þ T lymphocytes in DCL patients are functionally impaired and stimulation with TLR2 ligands restored the effector functions of CD8 þ T lymphocytes against macrophages infected with L. mexicana (Hernandez- Ruiz et al., 2010). Thus, the prediction of class I restricted immunogenic epitopes from the entire genome of Leishmania parasite as potential vaccine candidate is critically important. CD8 þ T cells are potent effector cells and produces inflammatory cytokines which could directly or indirectly reduce the parasitaemia. Immunotherapeutic role of CD8 þ T cells in leishmaniasis is a subject of intense debate and is likely varies, based on the infection model. Leishmania parasite demonstrates wide variations in tropism which is reflected in clinical manifestations. CD8 þ T cells functions will certainly come within the purview of this clinical multimorphism of the disease. In order to address this issue, we dedicated this work solely to identify the CD8 responsive immunogenic epitopes in human for understanding the possible disease course which could generate immunity against VL.
In the current study, we are exploring the possibility of designing and fabricating epitope-based vaccine protocol as a candidate vaccine protocol against Leishmania donovani (L. donovani) by exploiting large protein network, expressing nearly 7960 proteins (UniProt). We have performed computation-based identification of antigenic proteins responsible for invoking killer immune response in host. CTL mediated immune response is critical for clearing the pathogen elicited virulence for faster recovery and curing of the disease. Our approach of involving the entire protein repertoire and identifying the possible immunodominant epitopes against the disease is immensely important for possible application in designing novel vaccine candidate. This will select a set of peptides which could elicit optimal immune responses in the context of HLA class I molecules following vaccination. We have designed numerous epitopes for vaccine development against L. donovani parasite using vaccinomics method. The present work could be used as a basis for further analysis and designing vaccine against Leishmania parasite for future studies. In essence, we could answer some of the critical unanswered questions like degree and levels of protective immunity, how long it lasts and how to generate long lasting immunity in recurrence of the disease.

Materials and methods
2.1. Total protein sequence retrieval of the parasite A total of 7960 protein sequence from L. donovani was derived from UniProt (www.uniprot.com). All the proteins were in FASTA file format, required for our analysis since; many tools and server only accept the sequences when they are submitted in FASTA format.

Sequence analysis and prediction of candidate antigens
All the protein sequences were screened for the presence of signal peptides and Glycosyl phosphatidyl inositol (GPI) anchors using Signal P 3.0 (http://www.cbs.dtu.dk/services/ SignalP-3.0/) (Dyrløv Bendtsen et al., 2004) and PredGPI tools (http://gpcr.biocomp.unibo.it/predgpi/) (Pierleoni et al., 2008) respectively. SignalP 3.0 server predicts the location of signal peptide cleavage sites in submitted sequences based on HMM (Hidden Markov Model) and ANN (Artificial Neural Network) scheme. The default settings included the truncation of each sequence to maximum 70 residues. PredGPI is based on a support vector machine (SVM) which predicts GPI-anchored proteins by analyzing the anchoring signal. It also uses Hidden Markov Model (HMM) for the prediction of the most probable omega-site in protein sequences. The GPI anchor Predictor (PredGPI) server included the datasets as 145 proteins experimentally annotated as GPI-anchored (positive dataset), 10630 proteins experimentally annotated as non-GPI-anchored (negative dataset) and 26 proteins with an experimentally annotated omega-site (omega-site dataset). After that, the protein sequences were predicted for antigenicity using VaxiJen tool (http://www.ddg-pharmfac. net/vaxijen/VaxiJen/VaxiJen.html). It is a server used for in silico prediction of proteins with immunoprotective ability (Doytchinova & Flower, 2007;. The selected candidate proteins were further analyzed for having transmembrane region, solubility and the secondary structure prediction using TMHMM Server v. 2.0 (http://www.cbs.dtu.dk/services/ TMHMM/) (Krogh et al., 2001), SOLpro server (http://scratch. proteomics.ics.uci.edu/) (Cheng et al., 2005;Magnan et al., 2009) and CFSSP server (http://www.biogem.org/tool/choufasman/) ( Chou & Fasman, 1974;T, 2013) respectively. The TMHMM Server v. 2.0 finds transmembrane regions in protein sequences and helps to discriminate between intracellular and surface proteins. SOLpro is a sequence-based prediction method using two stage SVM (Support Vector Machine) architecture, based on multiple features of about 17000 proteins training data set (Magnan et al., 2009). SOLpro tool was integrated in the SCRATCH suite of predictor (Cheng et al., 2005). The sequences were scrutinized to predict their propensity for solubility. CFSSP (Chou & Fasman Secondary Structure Prediction Server) is an online server which predicts secondary structure regions such as alpha helix, beta sheet, and turns from the protein sequence. The method implemented is Chou-Fasman algorithm, which is based on analyses of the relative frequencies of each residues in alpha helices, beta sheets and turns based on known protein crystal structures ( Chou & Fasman, 1974;T, 2013 ).

Subtractive proteomic approach to exclude the host homologous proteins
BLASTp (Basic Local Alignment Search Tool for proteins) analysis were performed for antigenic proteins against the host (human) (taxid: 9606) proteomes at NCBI database (https:// blast.ncbi.nlm.nih.gov/Blast.cgi). The protein sequences showing hits with E-values less than 10 À4 were removed since these sequences may have certain levels of homology with the human proteome (Hosen et al., 2014;Rana et al., 2015).

Mhc class I restricted T cell epitope prediction
The immune epitope database analysis resource (IEDB-AR) tool (https://www.iedb.org/) was used for the prediction of MHC class I restricted epitopes from the selected set of L. donovani proteins (Kim et al., 2012). The tool is trained to predict the epitopes for human, rodent and other animals. The peptide should bind with the MHC molecules so that it could be recognized by the T cells. The tool calculates the inhibitory concentration (IC 50 ) values for the peptide binding to MHC molecules. On the basis of binding affinity using IC 50 values as a unit of measure, epitopes were classified as; high affinity binders (IC 50 <50 nM), intermediate affinity (IC 50 <500 nM) and low affinity binders (IC 50 <5000 nM) respectively. A lower value of IC 50 indicates the higher binding possibility of the peptides with the host MHC allele. T cell epitopes for human MHC class I alleles were predicted by submitting the L. donovani antigenic protein sequences in FASTA format to IEDB-AR tool. We have used forty-nine human alleles for the peptide binding analysis. ANN method was applied to identify nine-mer MHC class I peptides (Nielsen et al., 2003;Rana et al., 2015). This method was trained on data having continuous binding affinities. It is reported that the performance of ANN method is best among the other theoretical methods which were used for predicting peptide binding affinity (Lundegaard et al., 2008).
Further, the top peptides were subjected for immunogenicity prediction at IEDB (http://tools.iedb.org/immunogenicity), which uses position and properties of the amino acids of the peptide for this analysis (Calis et al., 2013). The IEDB server provides another tool known as MHC-NP (http://tools.iedb. org/mhcnp/) which predicts probability and demonstrates that the submitted peptide sequences were naturally processed and bound to the respective MHC allele (Gigu ere et al., 2013).

HLA allele retrieval and docking
We have downloaded HLA alleles from PDB (Protein Data Bank) (Berman et al., 2000). The probable epitopes were screened using iGEMDOCK tool (Mohammad et al., 2017;Yang & Chen, 2004) and then the top screened peptides were docked to the HLA alleles using PATCHDOCK (http:// bioinfo3d.cs.tau.ac.il/PatchDock/) server which is considered as the reliable server for protein-peptide docking (Bora et al., 2017). The server employs geometry-based docking algorithm to get the best solution. The RMSD value of 4 Å for clustering was used to remove redundant solutions (Bora et al., 2017).

Protein and peptide structure visualization
For visualization and analysis of peptides, proteins and docked complexes, we used PyMol and LIGPLOTv1.4.5 tools Wallace et al., 1995). Ligplot is used for 2 D representation of hydrogen bonding, and hydrophobic contact in docked complexes

Sequence retrieval and analysis using various tools and servers for the prediction of epitopes
We have performed extensive and exhaustive analysis of epitope identification and determination of the antigenic peptide, likely responsible for invoking CD8 þ T cells responses in host ( Figure 1). We have applied similar approaches performed previously for this type of analysis (Ali et al., 2017;Rana et al., 2015). Figure 1 showed overall strategy applied for the identification of MHC Class I epitopes. The screened proteins represent the ideal immunogenic candidates while the predicted epitopes from few of these immunogenic proteins suggested that they could be used for the development of diagnostic markers and vaccines. Herein, we have screened complete protein sequences of L. donovani (7960 proteins) using SignalP 3.0 server which predicted the location of signal peptide cleavage sites in submitted sequences. Out of the 7960 proteins, 1642 proteins were predicted to contain signal peptides/signal anchors by SignalP 3.0 server using default settings. Next, 1642 proteins were subjected to PredGPI server for predicting GPI anchors (data not shown). Further, the identified 51 proteins were analyzed for antigenicity with VaxiJen server by adjusting the threshold at 0.6 for parasites which selected 27 antigenic proteins having VaxiJen score above 0.6 ( Table S1). The sequences which showed similarity with human or mouse were discarded from the set of 27 proteins due to the possibility of unwanted autoimmune responses, associated with these sequences. Finally, we have selected 21 proteins through subtractive proteomic approach for further analysis. These 21 proteins were exclusively present in Leishmania parasite and were further analyzed for predicting transmembrane TM topology with HMM based model using TMHMM Server v. 2.0 (Krogh et al., 2001), solubility by SOLpro and the presence of secondary structure by CFSSP servers (Table 1). We conducted such analysis in order to get information about the proteins. It is documented that the proteins having more TM regions creates difficulty in cloning and purification. Some proteins are soluble and insoluble depending upon solvent availability; hence, these types of information are needed sometimes. It is also shown that many pathogenic proteins have beta-helices which are considered as antigenic part of the proteins and helps in infection (Hasan et al., 2013).
3.2. MHC-I binding predictions by IEDB-AR tool and selection of the epitopes T cell epitopes of human MHC Class I alleles were predicted by submitting the above 21 antigenic protein sequences in FASTA format to IEDB-AR tool. We have employed ANN method and used forty-nine human alleles for the peptide prediction against each submitted protein sequence. These set of alleles (Table S2) were selected by checking the checkbox of "Show only frequently occurring alleles" in MHC class I binding predictions tool of IEDB (Kim et al., 2012). This allowed the selection of only those alleles that occurred in at least 1% of the human population with 1% or higher allele frequencies.
We have predicted the peptide sequences showing biding affinity of IC 50 value <50 nM in order to find high affinity binders (Table S3). We have selected top 5 peptides from each protein and then subjected them to immunogenicity prediction using tool at IEDB which gave immunogenic score for each submitted peptide (Calis et al., 2013). Furthermore, we also used another tool known as MHC-NP which predicted the probability of the submitted peptide sequence whether they are naturally processed and bound to respective MHC allele (Table 2) (Gigu ere et al., 2013). On the basis of the results presented in Table 2, we have selected 30 peptides for docking/screening and enlisted them in Table 3. These peptides were selected based on good enough and optimum score predicted by MHC-NP tool which also predicts the probability of the submitted peptide sequence whether they are bound to respective MHC allele.
3.3. HLA alleles retrieval, peptide modeling and docking for selecting epitopes considered as vaccine candidate  Table 3. The 3D structure of the selected peptides was modeled by PEP-FOLD3 server. The peptide sequences were submitted one by one for modeling. Simulation was set to 200, and coarse-grained force field OPEP (optimized potential for efficient structure prediction) was used for sorting the generated models (Lamiable et al., 2016;Maupetit et al., 2007). First, we have screened 30 modeled peptides with the respective alleles and the interactions were shown with iGEMDOCK tool (Laoye et al., 2016). Only the best conformers, showing energy value < 80 kcal/mol were selected for further analysis (Table S4). Docking of the reference peptides (already present with the PDBs) was also performed to compare with the energy values of the docked 30 peptides. We have shown only 11 modeled structures of peptides in Figure 2 since, these peptides have shown better docked energy values in iGEMDOCK tool and were further selected for analysis. Moreover, it is unnecessary to show the structures of those peptides that are of no use further in the study. As presented in Table S4, Pep3 showed better binding energy with HLA-B Ã 35:01 compared with the reference peptide (KPIVVLHGY, present in PDB; 2H6P). Similarly, Pep7, Pep9, Pep10, Pep18 and Pep20 showed energy values less than 80 kcal/mol which was much better compared with the reference peptide of protein HLA-B Ã 44:03 (EEPTVIKKY, PDB ID; 1SYS). Likewise, other peptides (Pep2, Pep23, Pep25, Pep27 and Pep30) showed better energy values and were considered for further analysis. All the 11 peptides were then docked with the respective alleles using PATCHDOCK, a reliable server for protein-peptide docking (Bora et al., 2017) (Figure 3). We again docked the reference peptides using this server for better understanding of the protein-epitope docked complex (Table S5 & Figure S1). The results of PATCHDOCK based docking of 11 peptides are listed in Table S6. We have plotted the entire eleven protein-peptide complexes using Ligplot tool and observed different interactions as shown in Figure S2.
The results of PATCHDOCK suggested that eight epitopes from six proteins namely; epitope QLLPRLFPA extracted from the protein having uniprot id E9B8P9, epitope YPAFAALVF from protein E9BE51, epitopes VAVAATVAY and AVAATVAYV from protein E9B7Y6, epitope AAAPTEAAL from E9BTY2, epitope MYPLVAVVF from E9BD04, and the epitopes TLFGYVAAV and YVAAVTSAF from protein E9BRS9 were found to have better binding potential with respective alleles and they might elicit immune response. Based on the docking analysis using PATCHDOCK server, it was also suggested that the selected epitopes occupied the same binding pockets as the reference peptides and made several hydrogen bonds and hydrophobic interactions with the surrounding residues of the respective alleles (Figure 4 & S3). For analyzing the docked complexes, we used PyMoland LIGPLOTv1.4.5 tool Wallace et al., 1995). We have also  predicted the secondary structure portions of the selected six proteins from the amino acid sequences using CFSSP server ( Figure S4).
The epitope QLLPRLFPA (Pep2) extracted from uncharacterized protein (E9B8P9) has shown the binding affinity of 3.79 nM and the MHC-NP score of 0.5247 against HLA-A Ã 02:01 allele. It has shown the docking score of 6306 with HLA-A Ã 02:01 protein (PDB ID; 5ENW). This epitope has made several H-Bond and hydrophobic interactions with the docked protein residues like Tyr99, Glu63, Ala69 and Tyr159 which are present in the pocket. The epitope YPAFAALVF (Pep3), from protein E9BE51 (uncharacterized protein) has binding affinity of 3.31 nM and immunogenicity score 0.22019. It showed MHC-NP score of 0.421 against HLA-B Ã 35:01, which is a high value. The protein-peptide complex has made several hydrogen bonds with the pocket residues like Tyr159, Ala150, Arg97, and Asn70, present in the pocket of 2H6P protein. This protein is made of 88 residues with two transmembrane regions and it is soluble upon overexpression. It has 72.7% residues making sheets suggesting that it could be a good antigenic protein and may provide better predicted vaccine candidate.
There were two epitopes VAVAATVAY (Pep9) and AVAATVAYV (Pep10) screened from protein E9B7Y6; Surface antigen-like protein. This protein has only 2 transmembrane regions suggesting that it can be cloned and purified. During the secondary structure prediction, 44.9% of the residues are considered to make sheets. The epitopes have shown MHC-NP score of 0.8622 and 0.7125 against HLA-B Ã 44:03 respectively. The Table S6 has also shown several residues of protein (PDB ID; 1SYS) making interaction with Pep9 and Pe10 suggesting that protein-peptide complex has strong binding. This protein has given two high binding epitopes which could be considered for vaccine development process. The  Table 3. Selected list of peptides based on optimum immunogenicity and MHC-NP score represented in Table 2. These peptides were also enumerated and tabulated in another column, Epitope numbering. epitope AAAPTEAAL (Pep23), coming from protein E9BTY2 (Cyclin-e binding protein 1 like protein) has binding affinity of 2.91 and immunogenicity score was 0.18947. It showed MHC-NP score of 0.7889 against HLA-B Ã 53:01, which is a high value. The protein-peptide complex has made several hydrogen bonds with pocket residues like VAl152, Arg62, Tyr9, Tyr99 and Thr73 which present in the pocket of 1A1M protein. This protein is made of 719 residues with no transmembrane region and it is soluble upon overexpression. It has 51.7% residues making sheets suggesting that it is a good antigen. It has several peptides against submitted set of alleles out of which Pep23 could act as better predicted vaccine candidate. The protein, E9BD04 (Paraflagellar rod protein 2 C) has only one transmembrane TM region and it is also soluble upon overexpression with the score 0.902202. It has 76.7% residues making sheet. The epitope MYPLVAVVF (Pep27) has 0.7399 MHC-NP score against HLA-B Ã 53:01 and made docked complex with pocket residues of protein 1A1M. So, it is possible that the protein (Paraflagellar rod protein 2 C) and its epitope could be a better antigenic candidate for vaccine development. Similarly, two epitopes Pep20 (TLFGYVAAV) and Pep30 (YVAAVTSA) were predicted from protein E9BRS9 (Uncharacterized protein) has shown binding with HLA-B Ã 44:03 and HLA-B Ã 53:01 alleles respectively. This protein has also only one TM region suggesting it is not difficult to clone it. We have predicted several epitopes from this protein out of which two have shown better candidature to be considered as vaccine candidate.

Molecular dynamics simulation of allele-epitopes complex to observer stability
The docked MHC protein-peptide complexes were further studied for its structural stability by performing 100 ns molecular dynamics (   and might be considered as better vaccine candidates ( Figure 5(d)). The backbone RMSF data of each residue from all the complexes and four HLA proteins were plotted. The mobility of the amino acids after binding to each peptide was analyzed which indicated that higher RMSF value means greater flexibility and low binding properties. Pep2 does not  show low RMSF values while Pep3 have lower RMSF values compared with the apo-protein data respectively ( Figure S5A & B). Similarly, Pep9, Pep20, Pep23 and Pep27 have made better interactions during the 100 ns simulation ( Figure S5C & D). We have performed radius of gyration (Rg) analyses of all the systems which provided information on the stability of the proteins during simulation ( Figure S6).

Conservancy Analysis of epitopes in different Leishmania species
We have performed epitope conservancy analysis between the selected epitopes and the target Leishmania proteins sequences in order to find epitope distribution in different strains. For this, IEDB's web based tool (http://tools.iedb.org/ conservancy/) was used for the analysis of conservancy level (Bui et al., 2007). All the predicted four epitopes generated from Immunoinformatics approach and MD simulation had 100% maximum identity against L. donovani. The epitope MYPLVAVVF (Pep27) derived from Paraflagellar rod protein 2 C (E9BD04) has shown 100% maximum identity against L. donovani since, this protein is present only in this species. These epitopes were found to be conserved with conservancy level between 70-100% in L. major. Rests of the species have also shown some degree of conservancy which is enlisted in Table S7.

Discussion
Leishmaniasis has become a bigger threat to large number of populations in different parts of world with manifestation of severe forms of the diseases in patients with immune deficiencies including drug resistant variants. Current drug treatment has no satisfactory outcome and raised the chances of drug resistance. Several vaccine formulations are in various stages of development and trial with the hope and expectations for better treatment options against the disease. Scrutinization of GPI-anchored proteins, which are highly expressed in intracellular and infective stages of protozoans are considered as ideal antigenic targets, eliciting both humoral and cellular immunity (Bhatia et al., 2004).
Promastigotes of nearly all the species of Leishmania parasite highly express leishmanolysin (gp63), (GPI)-anchored glycoprotein which is involved in host-parasite interaction and affecting host defense system. Therefore, it was considered as an attractive vaccine candidate (Connell et al., 1993). We have taken a novel approach of selecting multiple epitopes (four in number) from four different proteins of L. donovani which were screened as antigenic proteins with high immunogenic properties. Low binding affinity is related to higher binding of peptide with the selected alleles, as suggested by IEDB tool. We have enlisted the epitopes with high affinity binders (IC 50 <50 nM) from each of the 21 antigenic proteins (Table S3). We have taken top five epitopes from each protein and submitted for Immunogenicity and MHC-NP analysis and the data is shown in Table 2. We then carefully and logically picked up 30 epitopes and enlisted them by giving names whose values were satisfactory and considered within the acceptable range of all the three tools. These peptides were then screened with the respective alleles with iGEMDOCK tool in order to get best fit peptides. PATCHDOCK tool was explored to demonstrate the binding of each of the selected 11 epitopes with their proteins. We have used peptides GLKEGIPAL, KPIVVLHGY, EEPTVIKKY and TPYDINQML present in the cavity of the selected proteins; 5ENW, 2H6P, ISYS and 1A1M respectively, for validating the docking simulations. The extensive analysis of the docked complexes was performed with the support of Ligplot and the figures were made with PyMol. Further, based on the MD simulation data it was suggested that Pep3, Pep9, Pep23 and Pep27 have better interactions during the simulation period. The binding of the peptides leads to the conformational changes of the proteins (MHC alleles) which proceed for further action. In order to provide the elegant view of protein-peptide interaction (for example; 2H6P-Pep3 complex), we captured the snapshots of the complex at different time interval (0 ns, 10 ns, 50 ns and 100 ns) ( Figure S7). It was observed that the peptide is better fitted and stabilized the protein by intermolecular interaction which persisted during the course of simulation. The figure showed the backbone movement of the allele which changes during interaction with the peptide. This suggested that the allele-peptide interaction can lead to the desired biological action. This also elaborates selectivity of the proteins for peptides. It was observed form RMSD analysis that simulated protein-peptide complexes were reasonably more stable than apo-proteins. These peptides may have the potential to bind with high stability and thus they may be presented by MHC alleles to CTLs which will further trigger an immediate immune response against those non-self-antigens (epitopes). They may be considered for further study in order to develop vaccine against leishmaniasis. For the development of epitope-based vaccine design, the conserved epitopes may provide broader protection against multiple species. In diagnostic kit development, the epitopes derived from a given pathogen may be used to monitor disease responses, caused by a particular virulent strain. For this purpose, crucial information's on the degree of conservancy of epitopes are needed (Bui et al., 2007). The overall strategy applied to extract best epitopes was explained in Figure 1 and various tools/servers and several algorithms are also summarized (Table S8). Although previous studies have shown use of these algorithms in other organisms (Rana et al., 2015;Singh et al., 2015), the proteome-wide analysis performed with L. donovani using these tools/servers for the prediction of vaccine candidate could provide significant understanding against the parasite.
Generation of Leishmania specific CD4 þ TH1 response contributes in elimination of the parasite via macrophage activation by IFN-c. Although CD4 þ T cells role is unquestionably important, activation and expansion of CD8 þ T cells during the healing process cannot be ignored. Evaluation of CD8 þ T cells role in inducing protection against the disease deserve attention. Success in effective vaccine design for a large heterogeneous population like human depends on the perfect tuning with HLA polymorphism. CD8 þ T cells reported to play critical role in controlling parasite growth in animal model of VL infection and reactivation of these responses causes dramatic reduction in parasite burden. Involvement of CD8 þ T cell derived IFN-c and perforin reported to play pivotal role in experimental VL infection model. Therefore, immune interventions that target CD8 þ T cell responses may have great therapeutic potential against VL in human. Obviously, effective peptide based vaccine against VL require multiple epitopes based on various HLA types particularly those common in endemic areas. The concept of polytope vaccine is important in this regard which has been known to induce cytolytic T cell response. Recently, HLA-A2 restricted CTL epitopes has been documented in Leishmania major infected population in Iran which generated CD8 þ T cell responses against peptides, derived from LmsTI1 and LPG-3 (Seyed et al., 2011). These results prompted us to undertake the proteome-wide search for effective CTL epitopes against VL which could provide new therapeutic tools to tackle this deadly parasitic disease.

Conclusion
The development of computer-based design of CD8 þ T cell epitopes and their response to treatment or vaccination can be an important milestone for an effective strategy to fight against leishmaniasis. This could offer the ability to expedite the therapeutic procedure dramatically, and to potentially benefit the time management for vaccine design and development against infections such as leishmaniasis. The CTLs of the immune system expected to discriminate between healthy and infected cells, since only the infected cells are to be eliminated. Quantitative estimation suggests that only 1 out of 200 peptides will bind a given MHC class I allele with adequate strength to score and secure a CTL response. Reliable predictions of immunogenic peptides can minimize the experimental effort needed to identify new epitopes to be used in, for example, vaccine design or for diagnostic purposes. Finding the possible epitopes using computational approach to complex problems such as predicting CTL epitopes in Leishmania potentially involves vast data sets. Our holistic approach of genome wide search for CTL epitopes against Leishmania could provide nearly all the possible predicted epitopes, expected to generate immune responses. Earlier, discrete efforts by others have only produced limited results with narrow window of opportunity. The current approach could eliminate nearly all the gaps in our understanding of designing peptide based vaccine against visceral leishmaniasis and likely to be helpful against cutaneous form as well.