In silico designing of a novel polyvalent multi-subunit peptide vaccine leveraging cross- immunity against human visceral & cutaneous leishmaniasis: An Immunoinformatics-based approach

2.1 Selection and retrieval of leishmania antigenic protein sequences:

The complete amino acid sequences of twelve leishmania-specific antigenic proteins were retrieved from Uniprot (Universal Protein Resource) database (https://www.uniprot.org/uniprot/) in FASTA format (accessed on 21.03.2022). All these proteins were reported as immunogenic, highly conserved and potential vaccine candidates as found in previous immunological studies involving either animal models or computational in silico prediction methods (Coler et al. 2005; Kedzierski, 2010; Nagill et al. 2011; Brito et al. 2018). The selected protein sequences belonged to any of the five pathogenic species of Leishmania, namely Leishmania major, Leishmania maxicana, Leishmania donovani, Leishmania chagasi and Leishmania amazoensis among which L. donovani is the primary parasite causing visceral leishmaniasis (VL) while the others are mostly responsible for causing cutaneous (CL) or muco-cutaneous leishmaniasis (MCL). Among the proteins under investigation, the LACK (Leishmania-activated C-kinase Antigen) protein is a highly conserved protein across Leishmania species and is considered a viable vaccine candidate against human leishmaniasis (Fernandez et al. 2018). It has been proved to be a key target of immune response in sensitive BAL b/c mice (Kelly et al. 2003). KMP-11 (Kinetoplastid Membrane Protein-11) is present in all kinetoplastid protozoa and is also considered a potential vaccine candidate that is shown to increase IL-10 levels in murine models indicating its immunogenic potential (Mendonka et al. 2015 & Atapour et al. 2021). Leishmanolysin or GP63 is a surface proteinase that has been postulated as a virulence factor involved in direct interaction of the parasite with the host macrophage receptor (Joshi et al. 2001). The membrane proteins LCR1 and GP46 are also found to increase IFN-gamma production and hence, are regarded useful for developing a general vaccine against Leishmania infections (McMahon-Pratt et al. 1993; Joshi et al. 2014). Heat shock proteins (HSPs) are highly conserved molecules which are highly immunogenic too inducing both MHC-I and MHC-II pathways of adaptive immunity and are thought to have a pertinent role in vaccine development against infectious leishmaniasis (Holakuyee et al. 2012). Leishmania Hydrophilic Acylated Surface Protein-B (HASP-B) is expressed only in infective parasites suggesting a role in parasite virulence and therefore, a potential vaccine candidate (MacLean et al. 2016). Kinetoplastid paraflagellar rod (PFR) proteins are also of therapeutic and prophylactic importance due to their restricted evolutionary distribution, high order organization and high immunogenicity (Maharana et al. 2015). Likewise, cysteine proteases and proteophosphoglycans also act as essential virulence factors for the parasites in mammalian hosts and are attractive drug targets for leishmaniasis (Mahmoudzadeh and McKerrow 2004; Mundodi et al. 2005; Rogers, 2012). Moreover, beta-tubulin protein had also previously been characterized as T-cell stimulating antigen from Leishmania although its efficacy as vaccine candidate has not been tried as yet (Bhowmick et al. 2009). Beta-tubulin of L.donovani is as well considered a potential drug target against visceral leishmaniasis (Assis et al. 2014). Lastly, the translation factor protein eIF5A from L. braziliensis had been shown to induce heterologous protection against leishmaniasis in animals when administered as vaccine in combination with another recombinant protein and increases the secretion of parasite-specific cytokines in vaccinated animals. eIF5A protein is also found to be conserved in all eukaryotes (Duarte et al. 2016).

The rationale behind the choice of such diverse types of conserved and immunogenic proteins across different species of Leishmania was to design a general peptide vaccine that would confer broad spectrum cross-immunity against both the forms of leishmaniases (VL and CL) that commonly affect humans. The antigenicity of the selected proteins were computed by employing the ANTIGENpro server (http://scratch.proteomics.ics.uci.edu/) that computes antigenicity of an input protein sequence by exploiting a SVM-based two stage classifier validated by 10-fold cross validation approach (Magnan et al. 2010). The proteins were then subjected to signal peptide analyses using SignalP- 4.1 server (SignalP − 4.1 - Services - DTU Health Tech) in order to discriminate classical secretory and non-secretory (transmembrane) proteins. This method incorporates a prediction of cleavage sites and a signal peptide/non-signal peptide prediction based on a combination of several artificial neural networks (Petersen et al. 2011 & Nielsen 2017). Localization predictions for all the proteins were performed next using DeepLoc (http://www.cbs.dtu.dk/services/DeepLoc/), a template-free algorithm that uses deep neural networks to predict protein subcellular localization exploiting only sequence information, achieving good accuracy (Armenteros et al. 2017). All functional protein sequences as obtained from SignalP 4.1 server were next subjected to predict the presence of antigenic determinants (epitopes) that will be recognized by B and T cell receptors. This was done to ensure that the designed vaccine construct contains only immuno-dominant T-cell and B-cell epitopes and so, would be able to drive an efficient immune response involving both humoral and cell-mediated pathways of immune mechanisms.

2.2 Cytotoxic T lymphocyte (CTL) epitope prediction:

Prediction of Cytotoxic T lymphocyte (CTL) epitope is an essential aspect for designing an ideal peptide vaccine (Ali et al. 2017). A freely accessible web-server namely NetCTL1.2 (https://www.cbs.dtu.dk/services/NetCTL) was utilized to predict the CTL epitopes for all the selected proteins. The method integrates prediction of peptide MHC class I binding, proteasomal C terminal cleavage and TAP transport efficiency (Larsen MV et al. 2007) all of which are essential features that a sequence of CTL binding epitope should possess. The server allows for predictions of CTL epitopes restricted to 12 MHC class I super type among which only A1 super type was used for the present study. MHC class I binding and proteasomal cleavage is performed using artificial neural networks. TAP transport efficiency is predicted using weight matrix (Larsen MV et al. 2007). The threshold value for epitope prediction was set at 0.75 (default) that corresponds to higher sensitivity (0.70) and specificity (0.970) for more accurate predictions. The server predicts motifs of nine amino acid residues from each input protein sequences and represents them with a particular score that reflects their MHC-I binding potency. Sequences with higher binding scores are considered potent CTL epitopes.

In addition to NetCTL 1.2, another publicly available web-server based approach namely SYFPEITHI (http://www.syfpeithi.de/bin/MHCServer.dll/EpitopePrediction.htm) was also availed for prediction of CTL epitopes from all the selected leishmania proteins. SYFPEITHI is a database comprising of more than 7000 peptide sequences known to bind class I and class II MHC molecules. The entries are compiled from published reports only. The prediction is based on published motifs (pool sequencing, natural ligands) and takes into consideration the amino acids in the anchor and auxiliary anchor positions, as well as other frequent amino acids (Rammensee et al. 1999 & Schuler et al. 2017). From the list of MHC types provided in this server, HLA-A*01, HLA-A*02:01 and HLA-A*03 were chosen for the prediction. The epitope search was specified for nonamer sequences (9 amino acid sequence motifs) for uniform comparison with the results given by the NetCTL 1.2 server. Only the common epitopes from each protein (if any) which were predicted by both of the NetCTL 1.2 and SYFPEITHI servers with compatible higher scores were selected for the final vaccine construct. This strategy was devised to improve the accuracy and sensitivity of the prediction.

2.3 Helper T lymphocyte (HTL) epitope prediction:

Helper T lymphocytes are probably the most important immune cells as they are required in all types of adaptive immune responses. They function essentially to stimulate antibody production by B cells, to activate macrophages for phagocytosis of the invading pathogens and also for the activation of the cytotoxic T cells to kill the infected target cells. Prediction of HTL epitopes (mostly of the Th-1 type), therefore, becomes an essential element for designing an effective peptide vaccine. In the current study, MHC-II binding analysis tool in the IEDB ( Immuno Epitope Data Base and Analysis) server (http://tools.iedb.org/mhcii/) was employed for prediction of 15-mer HTL epitopes for a set of three human alleles namely HLA-DRB1*01:01, HLA- DRB1*01:02 and HLA-DRB1*01:03. Human HLA alleles (as opposed to mouse alleles) were chosen for more realistic estimation of the MHC-II binding affinity of the epitopes as the vaccine was being designed against human leishmaniases. The predicted output is given in units of IC₅₀nM for combinatorial library and SMM align. Therefore, a lower number indicates higher affinity. As a rough guideline, peptides with IC₅₀values < 50 nM are considered high affinity, < 500 nM intermediate affinity and < 5000 nM low affinity (Ali et al. 2017). Most known epitopes have high or intermediate affinity. IEDB recommended prediction method was used in this approach in which low adjusted rank is indicative of good MHC-II binders and hence, can be defined as potent HTL epitopes.

2.4 Identification and Selection of cytokine inducing HTL epitopes:

Cytokines are of immense importance for proper functioning of the immune system and a battery of cytokines like IL-2, IL-4, IL-6 and interferron-γ have the ability to induce both CTL-mediated and humoral immune response (Pandey et al. 2018). Hence, the HTL epitopes that can trigger cytokine production are preferable for designing an immune-prophylactic vaccine. With this view, IFNepitope Server ((http://crdd.osdd.net/raghava/ifnepitope/)) was used for identification of HTL epitopes having cytokine inducing abilities. IFNepitope is an online prediction server that aims to predict the peptides from protein sequences having the capacity to induce IFN-gamma released from CD4⁺ T cells. The web server has been developed on the basis of a dataset which comprises of IFN-gamma inducing and non-inducing MHC class II binders (Dhanda et al. 2013). The HTL epitopes predicted by IEDB server with higher MHC-II binding affinities (lower adjusted ranks) were provided as input to check for their cytokine producing ability by the IFNepitope server. The epitopes thus identified were then subjected to next level of analyses by IL4pred server (IFNepitope (iiitd.edu.in) and also by IL10Pred server (https://webs.iiitd.edu.in/raghava/il10pred/predict3.php) in order to ensure that the selected HTL epitopes would as well evoke the secretion of IL-4 and IL-10, respectively. IL4pred and IL10Pred servers exploit motif-based as well as SVM based methods trained on experimentally validated positive and negative datasets. They also consider positional conservation of amino acid residues and amino acid composition of the peptides for more accurate prediction (Dhanda et al. 2013; Nagpal et al. 2016). The MHC-II binding HTL epitopes that were predicted to be positive by all the three servers were selected for incorporation in the final vaccine construct. This strategy was adopted so as to improve the efficacy of the vaccine since it would be able to instigate an intense response involving a wide variety of cytokines (IFN-gamma, IL-4 and IL-10) leading to an accelerated and long-lasting immunity.

2.5 Epitope conservation analysis:

Epitope conservation analysis becomes important to check whether the chosen epitope(s) is conserved across different species of Leishmania proteins. Designing a subunit vaccine consisting of conserved epitopes offers a promising scope for broad spectrum protective immunity against leishmaniases. Conservation analysis for each of the selected CTL and HTL epitopes was done by using IEDB “Epitope Conservation Analysis” tool available within the IEDB portal (http://tools.iedb.org/conservancy/). This tool calculates the degree of conservation of the target epitope (e) within a set of homologous protein sequences (P) as the fraction of {p} that matched the aligned e above the chosen identity level while considering that the target epitope is sequential in nature (Bui et al. 2007). Firstly, the corresponding source proteins from which the CTL and HTL epitopes were finally selected were subjected to a BLAST search for which we used the Uniprot BLAST tool (https://www.uniprot.org/blast/). From the BLAST result, the aligned protein sequences (belonging to the genus Leishmania) that showed significantly high similarity with the query sequence (> 40%) were selected. Such sequences were compiled to constitute an epitope-specific dataset and the same was provided as the input set of proteins (P) to be used for conservation analysis by the tool. Next, the sequences for selected CTL and HTL epitopes were screened individually against their respective datasets that revealed the relative conservation of the epitopes across a set of related leishmania proteins. The threshold identity level was set at > 100% (default).

2.6 Designing of multi epitope vaccine sequence:

On the basis of the immunoinformatic analyses, a primary sequence of the subunit vaccine was constructed by incorporating the selected CTL and HTL epitope sequences. These epitopes were linked together by AAY and GPGPG linkers, respectively (Ali et al. 2017; Khatoon et al. 2017; Pandey et al. 2018). Linkers are inserted for proper separation of the individual epitopes that is required for efficient functioning of the vaccine construct (Nezafat et al. 2014; Ali et al. 2017). Secondly, addition of an appropriate adjuvant is crucial in a subunit vaccine in order to boost up the immune response (Pandey et al. 2018). The choice of a suitable adjuvant is cardinal in designing vaccines for human trials. IL-12 produced by various immune cells is important in developing cell-mediated immunity against leishmaniasis and the potential of IL-12 as an adjuvant in vaccines against leishmaniasis has been reported in murine models (Afonso et al. 1994; Handman, 2001; Mutiso et al. 2010). In reference to this, IL-12 was selected as adjuvant for designing the vaccine and the sequence of the same was included in the N terminal of the vaccine construct by using EAAK linker (Ali et al. 2017; Khatoon et al. 2017; Pandey et al. 2018). The sequence of human IL-12 alpha chain was retrieved from Uniprot database (www.uniprot.org; Accession No- P29459). Finally, the peptide vaccine construct was obtained having adjuvant, linker, CTL epitopes and HTL epitopes (with intra-epitopic AAY and GPGPG linkers) in a sequence moving from N terminal to C terminal.

2.7 B cell epitope prediction:

Humoral immunity involves synthesis and secretion of specific antibodies by activated B cells which requires the presence of antigenic determinants that are recognized by B cell receptors present on the surface of B lymphocytes. An ideal immune-protective peptide vaccine should have such components in order to interact with and thereby, stimulate the B cells. Therefore, two servers, namely ABCpred and BepiPred 2.0 were applied for authentic prediction of B cell epitopes in the primary sequence of the designed vaccine. ABCpred server (https://webs.iiitd.edu.in/raghava/abcpred/ABC_submission.html ) is an artificial neural network (ANN) based approach which predicts B cell epitopes in a given antigen sequence using fixed length parameters (Saha et al. 2006). The training dataset used in this method contains B cell epitopes from B cell epitope database (BCIPEP). The server is able to predict epitopes with 65.93% accuracy using recurrent neural network (Saha et al. 2006). The primary sequence of the final vaccine construct was provided as the input in a plain format and the ABCpred server default parameters including a threshold value of 0.51 with a window length of 16 amino acid residues were used for the prediction.

BepiPred 2.0 server (BepiPred − 2.0 - Services - DTU Health Tech), on the other hand predicts sequential B cell epitopes from an input sequence by random forest algorithm trained on epitopes and non-epitope amino acids determined from crystal structures followed by a sequential prediction smoothing (Jespersen MC et al. 2017). The primary sequence of the vaccine was provided in a FASTA format as the input. The application of two different servers that predict B cell epitopes on the basis of two different algorithms makes the prediction more realistic.

2.8 Profiling of antigenicity, allergenicity and toxicity of the vaccine construct:

Antigenicity is an important criterion for an immune-prophylactic vaccine and to ensure that the designed vaccine construct would be able to elicit long-lasting immune response by interacting with the B and T cell receptors, antigenicity of the same was evaluated by using ANTIGENpro and VaxiJen 2.0 servers. ANTIGENpro (http://scratch.proteomics.ics.uci.edu/) server correctly classifies 82% of the known protective antigens when trained using only the protein microarray datasets. The accuracy on the combined dataset is estimated at 76% by 10-fold cross-validation experiments that allows a significantly better recognition of antigenic peptides. It runs on a two stage architecture including a SVM-based second stage classifier. For a new input protein sequence, the probability computed by the second stage SVM predictor is the final ANTIGENpro prediction score (Magnan et al. 2010). VaxiJen2.0 (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html) is another publicly accessible web-server that evaluates antigenicity of a protein sequence by an alignment-free approach which is based on auto cross covariance (ACC) transformation of protein sequences into uniform vectors of principal amino acid properties and allows antigen classification solely based on various physicochemical properties of the protein without referring to sequence alignment (Doytchinova et al. 2007).

Secondly, in order to rule out any possibility that the vaccine construct could trigger an allergic response in an individual, AlgPred (https://webs.iiitd.edu.in/raghava/algpred/submission.html) server was used to determine its allergenicty. The server tool provides a variety of approaches including mapping of IgE epitopes, MEME/MAST motif, SVM method based on dipeptide composition and BLAST ARPs all of which are exploited for a meticulous prediction. The threshold was set at -0.4 (default) for higher precision of the prediction. In addition, AllerTop v.2.0 (https://www.ddg-pharmfac.net/AllerTOP/index.html), a freely accessible web-server based tool was also applied to confirm the non-allergic nature of the vaccine construct. This method is based on auto cross covariance (ACC) transformation of protein sequences into uniform equal-length vectors. It has been applied to quantitative structure-activity relationships (QSAR) studies of peptides with different length and the proteins are classified by k-nearest neighbor algorithm (Dimitrov et al. 2014).

Toxicity, if any, was predicted by utilizing the ToxinPred (https://webs.iiitd.edu.in/raghava/toxinpred/protein.php) server. This prediction tool works on the basis of both machine learning technique and quantitative method that is trained on a dataset of various toxic and non-toxic peptides obtained from SwissProt and TrEMBL (Gupta et al. 2013). Toxicity profiling is crucial to assess the safety of the vaccine and other such immunotherapeutic peptides that are intended for human trial. The analysis was performed by choosing SVM (TrEMBL) + Motif based approach at an E-value threshold of 10.0 for motif based method and at a SVM threshold of 0.1.

2.9 Evaluation of physicochemical properties:

Describing the physical and chemical features of a chimeric peptide becomes important since when administered as a vaccine, it should be able to induce a proper immune response. The ProtParam server available with the ExPASY (Expert Protein Analysis System) portal (https://web.expasy.org/protparam/) was used to compute various physicochemical parameters of the primary vaccine construct that included features like molecular weight, theoretical pI, instability index, aliphatic index, in-vivo and in-vitro half-life, Grand Average of Hydrophobicity (GRAVY) and others.

2.10 Secondary structure prediction of the vaccine peptide:

Secondary structural elements of the designed vaccine construct were predicted by applying PSIPRED server (http://bioinf.cs.ucl.ac.uk/psipred/) and RaptorX Property server (http://raptorx.uchicago.edu/StructurePropertyPred/predict/) both of which are freely accessible on line protein structure analysis servers. The primary sequence of the vaccine peptide was provided as the input. PSIPRED predicts the secondary structure of the input protein sequence by exploiting two feed-forward neural networks that function on the basis of position-specific scoring matrices (PSSM) generated by PSI-BLAST (Jones et al. 1999). It incorporates identification of sequences that are homologous to our vaccine construct by PSI-BLAST followed by generation of PSSM. PSIPRED 3.2 and higher versions were found to achieve an average Q₃ score of 81.6% as evaluated by a stringent three-fold cross validation method which makes this approach significantly precise and accurate. RaptorX Property web server engages a very recent machine learning technique called DeepCNF (Deep Convolutional Neural fields) for predicting various features like secondary structure, solvent accessibility and disordered regions, simultaneously. DeepCNF is an integrated approach combining both Conditional Random Fields (CRFs) and shallow neural networks that models complex sequence-structure relationship by a deep hierarchical architecture and can obtain approximately 84% Q3 accuracy for 3-class secondary structure (Wang et al. 2016).

2.11 Tertiary structure prediction:

The tertiary structure of the final vaccine peptide was predicted by using the on line web server I-TASSER (https://zhanglab.ccmb.med.umich.edu/I-TASSER/). I-TASSER is a hierarchical protein structure modeling approach based on Profile-Profile Threading Alignment (PPA) and iterative structure assembly simulations followed by atomic level structure refinement (Zhang, 2008; Xu & Zhang, 2011; Yang & Zhang, 2015).It works on three stages: structural template identification, iterative structure assembly and structure based function annotation (Yang & Zhang, 2015). The top ranked structure models with global and local accuracy estimations are returned with their representative C-scores and TM-scores which are correlated with the quality of the predicted model. The I-TASSER program represents one of the most successful methods demonstrated in CASP for automated predictions of protein structure and function (Yang & Zhang, 2015).

2.12 Refinement of vaccine tertiary structure:

Conventional computational modeling of a protein structure alone does not guarantee the authenticity and accuracy of the predicted model since such modeling strategies largely depend on the degree of likeliness of the input (target) with the available template structures. Enhancement of the quality of the template based model beyond the accuracy was, therefore, thought to be necessary and this could be achieved by refining the whole protein structure. With this perspective, the tridimensional model output from I-TASSER server with the best C-score was subjected to further refinement by GalaxyRefine server (http://galaxy.seoklab.org/cgi-bin/submit.cgi?type=REFINE). It is a CASP10 tested refinement approach in which protein side chains are rebuilt first and then side chain repacking and overall structure relaxation are performed by molecular dynamics simulation. This leads to precise improvement of both local and global tertiary structures of target proteins. (Heo et al. 2013).

2.13 Validation of vaccine tertiary structure:

The refined tertiary structure of the vaccine needed to be validated to spot potential errors that might occur in the predicted 3D model. ProSA-web server, (https://prosa.services.came.sbg.ac.at/prosa.php) an interactive web-based platform based on the standard ProSa programme was used for this purpose. It uses knowledge based potentials of mean force to evaluate model quality and subsequently generates z-scores that indicate model quality as well as overall deviation of total energy of the predicted model from energy obtained from random conformations (Wiederstein et al. 2007). The refined model obtained from GalaxyRefine server was provided as the input structure in a .pdb format. Besides, Ramachandran Plot is another very useful approach for visualizing energetically allowed and disallowed regions of phi (ϕ) and psi (ψ) dihedral angles of the amino acid side chains that are considered crucial components for validation of a protein structure. Hence, Molprobity server (Main page - MolProbity (duke.edu))was employed that generated Ramachandran plot for the given model of the vaccine based on python based molprobity style phenix validation (Williams et al. 2018). We further used ERRAT server (SAVESv6.0 - Structure Validation Server (ucla.edu)) that analyzes the statistics of non-bonded interactions between different atoms in the model and plots the value of error function versus a 9-residue sliding (Colovos et al. 1993). This step was performed to get more sophisticated and precise validation of the modeled structure.

2.14 Prediction of continuous and discontinuous B cell epitopes in the vaccine peptide:

Most of the antigenic determinants recognized by the B cell receptor and antibodies are discontinuous meaning they contain amino acid residues that are located distantly in the primary structure of the immunogen but are brought adjacent to each other during the folding of the protein (Barlow et al. 1986; Van-Regenmortel, 1996). Prediction of such discontinuous B cell epitopes in the validated tertiary structure of the vaccine was performed by utilizing ElliPro server (http://tools.iedb.org/ellipro/) incorporated within the IEDB portal. ElliPro is a freely accessible web tool that implements modified Thronton’s method and takes into account each residue’s center of mass rather than the Cα atoms. By doing so, it relates each of the epitopes to an assigned score defined as PI (Protrusion Index) averaged over epitope residues. Discontinuous epitopes are defined based on their respective PI values and clustered based on the distance between residue’s centers of mass(R). Larger values of R correlate to larger discontinuous epitopes being predicted (Ponomarenko et al. 2008). At the same time, linear B cell epitopes in the input protein model can also be predicted and visualized by using ElliPro. Detection of both linear and discontinuous B cell epitopes in the refined 3D model of the vaccine ensures a successful vaccine designing strategy.

2.15 Disulfide engineering of the vaccine peptide:

Disulfide bonds contribute to protein stability by lowering the conformational entropy and also by increasing the free energy of the denatured state. Introduction of novel disulfide bonds has been considered as a key biotechnological tool to improve the thermostability of native, folded proteins. Disulfide by Design 2 (DbD 2) v 2.12 server was used to enhance the overall structural stability of the designed vaccine by introducing disulfide linkages. This method utilizes native geometry and calculates an energy value for each potential disulfide, thereby providing a means to rank the candidate disulfide bonds. The disulfide bonds that confer maximum stability were the candidates with high B factors (Craig et al. 2013). The flexible regions of the peptide were selected based on high B factor values and corresponding stabilizing mutations were created by forming disulfide linkages between the selected residue pairs (Khatoon et al. 2017).

2.16 Molecular docking of vaccine construct with immune receptors (Human TLRs-2, 4, 5, 8 and mouse TLR-9):

Toll like receptors (TLR) are major attributes of cellular immune responses which recognize pathogen-associated molecular patterns (PAMPs) and evoke innate immune responses against infections (Faria et al. 2011). The direct activation of TLR-2 with leishmania components was subsequently reported (Faria et al, 2011). In other related experimental studies, lack of TLR-4 was resulted in increased parasitic growth and delayed healing of cutaneous lesions caused by L. major infections indicating the plausible role of TLR-4 in inducing immunity against leishmania (Kropf et al. 2014). On the other hand, TLR-2, 4 and 9 were found to be involved in immuno-pathologic spectrum of CL caused by L. braziliensis and L. amazonensis (Campos et al, 2018). TLR-9 was also evaluated as a key player in activation of dendritic cells in pathogenesis of VL in humans (Tuon et al. 2008). Some of the leishmania-derived components have been shown to activate TLR-2, 4 and 9 in majority of studies conducted in this field (Faria et al. 2011). These observations emphasize the fact that various TLRs are involved in triggering protective immunity against Leishmania. Therefore, an interaction between the TLRs and the designed vaccine construct was considered to be necessary for eliciting effective immunity against leishmania infections. This was checked by performing molecular docking of the multi epitopic vaccine peptide with human TLRs- 2, 4, 5, 8 and mouse TLR- 9 using PatchDock server (PatchDock Server (tau.ac.il)). The PatchDock algorithm has three main stages namely, molecular shape representation; surface patch matching and filtering and scoring. It is a geometry-based molecular docking algorithm that principally operates on molecular shape complementarity between the receptor and the ligand. The complementary patches are matched and candidate transformations are generated. Each of these transformations is further ranked by assigning a scoring function that considers both atomic desolvation energy and geometric fit (Schneidman et al, 2003 & Schneidman et al 2005). In our study, human TLR-2 (PDB id: 4G8A), TLR-4(PDB id: 6NIG), TLR-5 (PDB id: 3J0A) and TLR-8 (PDB id: 4QC0) and mouse TLR-9 (PDB id: 3WPF) were selected as receptors and their respective PDB files were retrieved from Protein DataBank (www.rcsb.org). The refined 3D model of the vaccine peptide was used as the ligand for all the docking simulations. The PatchDock server was set at default parameters (clustering RMSD: 4.0).

We further analyzed best docking model using Cluspro server which is based on possible bond rotation algorithm & possible interactions between amino acids were analyzed using Pymol software (Kozakov et al. 2017; Vajda et al. 2017; Desta et al.2020).

2.17 Codon optimization, mRNA structure prediction and in silico cloning for expression of vaccine protein:

Optimization of codon is required to achieve maximum expression of foreign genes in the host organism when the usage of codon by the host differs significantly from that of the native host from which the original sequences for the final vaccine construct had been culled. JCat (Java Codon Adaptation Tool) available publicly at (http://www.jcat.de/) was used for reverse translation and codon optimization by providing primary sequence of the final vaccine construct as the input. The codon usage was optimized to the most sequenced prokaryotic organism, E.coli K12 (Grote et al. 2005 & Khatoon et al. 2017). Three additional options provided by the tool were selected to avoid rho-independent transcription terminators, prokaryotic ribosome binding sites and unwanted restriction cleavage sites. The codon adaptation index (CAI) and GC content as calculated by the JCat tool are indicative of how good the optimization was. Better optimization ensures higher expression of the vaccine protein. Further, restriction sites for the enzymes XhoI and NdeI were incorporated into the optimized cDNA sequence provided by the JCat tool to facilitate cloning. The modified sequence (with restriction sites) was then inserted into the E.coli pET 28-a(+) plasmid vector by using the restriction cloning module of SnapGene tool (Khatoon et al. 2017; Aathmanathan et al. 2018; Pandey et al. 2018; Khan et al. 2019). SnapGene software is available on the web at (https://www.snapgene.com/).

Besides, we had also used the mfold web server (http://unafold.rna.albany.edu/?q=mfold/RNA-Folding-Form) for prediction of secondary structure of the mRNA encoded by the genes of the chimeric peptide. The reverse translated optimized sequence obtained from JCat tool was provided as the input. The mfold program is based on a core algorithm that predicts a minimum free energy, ∆ G as well as minimum free energies for folding that must contain any particular base pair. Various folding constraints are also applied for accurate prediction of folding energies. The folding energy for RNA folding is fixed at 37⁰C while the ionic conditions are fixed at [Na⁺] = 1M and [Mg⁺⁺] = 0 (Zuker, 2003).

2.18 C-immsim based immune simulation:

C-immsim ( https://kraken.iac.rm.cnr.it/C-IMMSIM/index.php?page=1 ) is server based immune response prediction tool which utilized agent based class & it’s prediction rely on neural networks. It can predict both humoral & cell-mediated immune responses. C-immsim utilizes Celda-Seiden model, bit-string polyclonal lattice model & Simpson index etc. to characterize prediction of active, resting & memory cells and their longevity. It helps to predict cytokine profile of given reaction to evaluate inflammatory responses (Rapin et al. 2010 & Castiglione et al. 2021).

2.19 Molecular Dynamic simulation:

Molecular dynamic simulation was performed to evaluate stability of our best complex i.e TLR with proposed vaccine using iMOD (iMod Server home page (csic.es) server. This web based server uses improved normal mode (NMA) analysis for characterization of various dynamic properties of biomolecules by measuring internal coordinates, protein flexibility & stability using intra-atomic force field. In this current study, we analyzed the main chain deformities, B-factor, egienvalues, covariance factor & elastic network model as standard NMA parameters for evaluation of protein-protein complexes obtained by in silico approach (Alexandrov et al. 2005; Lopéz-Blanco et al. 2011; López-Blanco et al. 2014 ; Bauer & Bauerová-Hlinková 2020).

3.1 Selection and retrieval of leishmania antigenic protein sequences:

The amino acid sequences of twelve leishmania proteins were retrieved from Uniprot database in FASTA format and subjected to further analyses for the purpose of designing a multi subunit peptide vaccine against leishmaniases. These proteins were selected based on previous studies that reported them as potential vaccine candidates as evaluated in experimental animal models or by in silico prediction tools. Antigenicity analysis by the ANTIGENpro server showed all the chosen proteins as probable antigens with varying degrees of antigenicty that is in keeping with the information obtained from literature mining. However, ten among the twelve proteins were presented with relatively higher prediction scores (> 0.51) reflecting their higher potency for vaccine design (Supplementary table S1). The proteins were discriminated as secretory and non-secretory based on the prediction of signal peptide cleavage sites by the SignalP 1.4 server. Two of the selected leishmania proteins viz. L. mexicana cysteine proteinase A (Uniprot id: P25775) and L. amazoensis surface glycoprotein GP-46/M2 (Uniprot id: P21978) were found to contain signal peptides and hence, predicted as secretory proteins. Rest of the ten leishmania antigens were not found to have any signal peptide indicating their non-secretory nature. Localization analyses predicted that the protein cysteine proteinase A of L. mexicana (Uniprot id: P25775) to be localized in lysosome/vacuole, proteins leishmanolysin C1 and proteophosphoglycan of L.mexicana (Uniprot ids: P43150 and Q9TW13; respectively) to be membrane proteins, LACK protective antigen of L. donovani (Uniprot id: Q9BIJ5) to be in the nucleus and Heat shock 70 related protein 1 of L. major (Uniprot id: P12076) to be localized in mitochondria. All the other proteins were predicted as cytoplasmic. In addition to this, the amino acid sequence of alpha chain of human IL-12 was also obtained from Uniprot database (Uniprot id: P29459) to be used as adjuvant in the final vaccine construct.

3.2 Prediction of CTL and HTL epitopes:

From all the proteins being investigated, a total number of 59 common nonameric CTL epitopes were identified jointly by NetCTL 1.2 and SYFPEITHI servers. Each of the predicted CTL epitopes was assigned a particular score computed by the respective servers. Only 26 of these common 9-mer CTL motifs were selected due to their higher scores (indicating higher MHC I binding affinity). These twenty six CTL epitopes identified commonly by both the servers with compatible higher scores were included in the final vaccine construct (supplementary table S2A). Likewise, IEDB MHC-II binding analysis tool predicted 15-mer HTL epitopes from all the selected proteins against a set of three human HLA alleles (HLA- DRB1*01:01, HLA- DRB1*01:02 and HLA- DRB1*01:03). The predicted HTL epitopes with lowest adjusted rank for each one of these alleles were documented as in supplementary table S2B. From this pool, a total number of nine sequences were chosen on the basis of their low adjusted ranks (Table 1) that make these HTL epitopes the most potent MHC-II binders for any of the three HLA alleles. The selected HTL epitopes were next subjected to further computational analyses in order to predict their cytokine producing abilities.

3.3 Screening for HTL epitopes with cytokine inducing abilities:

All of the nine chosen HTL epitopes (Table 1) were evaluated for their ability to produce interferon-γ by employing IFNepitope server. The server predicted (IFN-γ Vs. Non IFNγ, using SVM & Motif based hybrid analysis) only three such peptides as FRRLYKTLGQLVYKK from paraflagellar rod protein 2C of L. major, DNGAYLGMEPSSVAA from proteophosphoglycan of L. mexicana and RLHFFMMGFAPLTSR from tubulin beta chain of L. mexicana as “positive” with scores lesser than 1.0 (Table 2).

These three peptides were further checked for their IL-4 and IL-10 producing abilities by using IL4pred server and IL10pred server, respectively. Both the servers utilize a SVM classifier and consider various SVM input parameters like amino acid composition, dipeptide composition, amino acid propensity and physico-chemical properties for the prediction (Dhanda et al. 2013; Nagpal et al. 2017). The respective servers predicted that these three 15-mer IFN- γ producing epitopes were inducers for IL-4 and IL-10 as well (Supplementary table S3).

Hence, these three MHC-II binding 15-mer HTL epitopes were finally selected for construction of the vaccine considering their ability to elicit the production of cytokines IFN- γ, IL-4 and IL-10 and thereby, mediating a heightened immune response.

3.4 Epitope conservation analysis:

IEDB module “conservation across antigen” available within the IEDB server was availed to find the degree of conservation of the selected CTL and HTL epitopes that would be incorporated in the final vaccine construct. Each of the twenty six CTL epitopes and three HTL epitopes were screened individually against a specific set of homologous protein sequences obtained by Uniprot BLAST. The BLAST results showed high degree of conservation of the selected proteins in various species of Leishmania and also in other related parasitic protozoa, mostly Trypanosoma. All the HTL epitopes were found to be highly conserved across the set of related leishmania proteins used in the analysis while the range of conservation for selected CTL epitopes varied from 0%-100% (Supplementary table S4). The incorporation of such conserved CTL and HTL epitopes in the final vaccine construct would ensure broad spectrum cross-immunity against multiple species of pathogenic Leishmania that commonly infect humans.

3.5 Construction of multi-epitope vaccine sequence:

The primary sequence of the multi-subunit peptide vaccine was constructed by merging the sequences of twenty six CTL epitopes and three HTL epitopes selected on the basis of their antigenicity, MHC binding affinity and cytokine inducing ability. The individual CTL and HTL epitopes were amalgamated by using AAY and GPGPG linkers, respectively and the overall length of the resulting vaccine peptide was 369 amino acids. The adjuvant, 219 amino acid long alpha chain of human IL-12 (Uniprot id: P29459) was added to the N terminal of the designed peptide by EAAK linker. Addition of IL-12 as a natural adjuvant was supposed to enhance the affectivity of the vaccine by accelerating the immune response as already validated in animal studies (Afonso et al. 1994; Mutiso et al. 2010; Raman et al. 2012). In addition, a 6XHis tag was also added to the C terminal of the peptide to facilitate the identification and purification of the protein following its production by recombinant DNA technology. After the addition of linkers, adjuvant and histidine tag, the length of the final vaccine peptide was found to be 598 amino acid residues (Fig.-1).

3.6 Prediction of B cell epitopes:

The presence of antigenic determinants in the primary sequence of the designed vaccine was identified by employing two different servers namely, ABCpred and BepiPred 2.0 for more convincing prediction. ABCpred server predicts B cell epitopes by using trained recurrent neural network and represents the candidate sequences as 16-mer peptides ranked according to their binding scores (Supplementary table S5). Higher score of the peptide means higher probability to be an epitope. The designed vaccine peptide was found to have seven non-overlapping 16-mer peptides with predicted scores 0.90 and above at a default threshold value of 0.51 reflecting their very high potency to bind and stimulate B cells. The other server, BepiPred 2.0 predicts the probability of each amino acid residue in the input sequence to be a part of the B cell epitope and depicts the output as an illustration with an orange color gradation indicating the degree of probability. The predicted residues that belong to B cell epitopes are marked as “E”. The server identified a large number of residues to be linear B cell epitopes with high probability within a residue stretch of 380–530 position of the peptide (Fig.-2) when the epitope threshold was set at 0.5 (default). Thus, the identification of antigenic determinants for B cell receptors by two independent servers established the fact that the peptide, when administered as vaccine, would be able to stimulate the production of antibodies by B cells to protect individuals from leishmania infections.

3.7 Antigenicity, Allergenicity and Toxicity profiling of the vaccine construct:

The antigenicity of the final vaccine sequence (with adjuvant) was evaluated by ANTIGENpro and VaxiJen 2.0 servers. ANTIGENpro server predicted the same to be a probable antigen with antigenicity score of 0.714582. VaxiJen 2.0 server evaluated the vaccine peptide as a probable antigen with a score of 0.4723 at a threshold of 0.4 when target organism was bacteria and 0.5325 at a threshold of 0.5 when target organism was selected to be parasites.

However, when the original vaccine sequence without the adjuvant was analysed by VaxiJen 2.0, it too was predicted as probable antigen with a score of 0.5768 and 0.4648 in parasite and bacterial models, respectively. The ANTIGENpro prediction score for the original sequence was 0.745396. Hence, both the constructed vaccine sequences (with and without adjuvant) were found to be antigenic in nature supporting the fact that the peptide components of the vaccine were antigenic by themselves even in absence of the added adjuvant.

As far as allergenicity is concerned, the final vaccine sequence was found to be non-allergen by AlgPred with a predicted score of 0.4018 at a threshold of -0.4. In addition, the protein sequence was not found to contain any experimentally proven IgE epitope and no hits were found for BLAST results for allergen-representative peptides (ARPs) indicating the non-allergen nature of the vaccine. AllerTOP v.2.0 server also defined the protein as probable non allergen. Toxicity analysis by ToxinPred server predicted the peptide to be non-toxin when SVM (TrEMBL) + Motif based approach was adopted at a SVM threshold of 0.1 and E-value cut off for motif based was set at 10.0. Hence, the generated vaccine sequence was predicted as a probable antigen, non-allergen and non-toxin all of which are important criteria for an ideal multi-subunit peptide vaccine.

3.8 Evaluation of physico-chemical parameters:

The physico-chemical parameters of the final vaccine construct were computed by Protparam server on the basis of the amino acid residues present in the given sequence. The molecular weight (MW) of the designed peptide was calculated to be 65.68 kDa that makes it suitable to induce an immunogenic response. The theoretical pI (Isoelectric Point) was computed as 5.98 that designates the peptide as slightly acidic in nature. The total numbers of positively charged and negatively charged residues in the peptide were 48 and 55, respectively. The half-life was calculated as 30 hours in mammalian reticulocytes in vitro; more than 20 hours in yeast and more than 10 hours in E.coli, in vivo. Persistence of the protein for 30 hours in mammalian cells (including humans) was sufficient for eliciting desirable immune response since the mechanism requires processing and presentation of the protein by the immune cells. The aliphatic index (relative volume occupied by the aliphatic side chains) was predicted to be 78.34 indicating the protein as thermostable in nature as aliphatic index may be considered as a positive factor contributing to the increase in thermo-stability of globular proteins (Wilkins et al. 2005). The predicted Grand Average of Hydropathy (GRAVY) value of the vaccine construct was predicted to be -0.140; the negative value defines it as hydrophilic in nature and will interact with water molecules (Kyte et al. 1982; Ali et al. 2017; Shey et al. 2019). To conclude, the immunoinformatic analyses predicted the designed peptide construct as antigenic, thermostable, persistent in mammalian cells and hydrophilic. All of these factors impart to the effectiveness of the construct as a probable vaccine candidate.

3.9 Prediction of secondary structure:

Prediction of secondary structure for the 592 residue long chimeric peptide by PSIPRED server revealed the presence of 45.77% alpha helix, 4.89% beta sheet and 49.32% coils as depicted in the following Fig. 3. In addition, RaptorX structure prediction server analysed the relevant solvent accessibility for the peptide in which 36% residues were predicted to be exposed, 28% were predicted as medium buried and 35% residues were predicted as buried indicating that the protein will have an overall fair contact with water (solvent). Only 32 residues (5% of the peptide) were predicted to be in disordered domain of the protein by RaptorX server.

3.10 Tertiary structure prediction:

I-TASSER web server was utilized for the prediction of 3-D structure of the designed vaccine. The server predicted five models for the said peptide on the basis of multiple threading alignments using ten different templates among which the template with PDB id 7e2cl showed the best alignment with a TM-score of 0.931 and a RMSD value of 2.75. The top-ranked model exhibiting the highest confidence score (C-score) of – 1.11 and an estimated RMSD of 10.3 ± 4.6Å (Fig. 4a) was selected for further refinement. A higher C-score indicates higher level of confidence with which the server predicts the tertiary structure of the target protein based on the predictions obtained from modeling simulations (Zhang, 2008; Shey et al, 2019). The C-score has a strong correlation with the overall quality of the tertiary structure and it has been used widely for quantitative estimation of the RMSD and TM-scores of predicted models in comparison with the native models. Further, the estimated TM-score of the chosen model was found to be 0.58 ± 0.14 which is another strong indicator of good modeling strategy since a TM score of > 0.5 indicates a model with correct fold. TM score is a sequence-length independent parameter for measuring structural similarity of the predicted model with other identical proteins in the same SCOP/CATH fold family (Yang et al. 2015).

3.11 Refinement and validation of vaccine 3-D model:

The “crude” tertiary model for the peptide vaccine was subsequently refined by processing it through GalaxyRefine server. GalaxyRefine yielded five models out of which model 1 was found to be the best after the evaluation and comparison of the associated parameters such as GDT-HA (0.8826); RMSD (0.608); MolProbity (2.398); clash score (21.1); poor rotamers (0.6) and Rama-favoured region (88.8) with the rest of the models (Fig. 4b).

The Ramachandran plot generated by MolProbity server showed that the refined 3-D model has 88.8% residues in the favoured region, 98.6% residues in the allowed region and only 0.013% outliers. In contrast, the initial crude model was found to have 77.1%, 94.2% and 0.057% residues in the favored, allowed and outlier regions, respectively. These results clearly conclude that the quality of the refined tertiary model was significantly improved. ERRAT server were also exploited for additional validation of the refined model. ERRAT server analyzed the model with an overall quality factor of 43.542 that actually represents the percentage of the input protein model which falls below the 95% rejection limit. Thus, we had successfully obtained an appreciably refined and validated tertiary structure for the vaccine construct. The refined tertiary structure of the vaccine along with the plots obtained from MolProbity server and ERRAT server are represented in Fig. 5. ProSA-web server had shown a Z score of -1.08 (supplementary image M1) in which the negative value indicates that the model falls outside the scores’ range of similar proteins with comparable sizes whose structures had been determined by X-ray diffraction or by NMR (Wiederstein et al. 2007; Khatoon et al. 2017).

3.12 Prediction of conformational and linear B cell epitopes:

Ellipro suit (available within IEDB portal) predicted a total of 323 amino acid residues located in eight discontinuous B cell epitopes in the refined tertiary structure of the peptide vaccine with scores ranging from 0.509–0.802. Among them, the largest conformational epitope contained 107 amino acid residues with a predicted score of 0.772. The default parameters of the server were used for the analysis. Ellipro also predicted twelve linear B cell epitopes with a total number of 324 residues within a score range of 0.535–0.885 (Supplementary table S6 & supplementary image M2).

3.13 Disulfide engineering of the vaccine peptide:

Disulfide by Design 2 server retuned 44 pairs amino acid residues that are probable sites for disulfide bond formation.

3.14 Analyses of docking interactions of the vaccine peptide with the TLRs:

The interactions of the chimeric peptide with the immune receptors (human TLRs-2, 4, 5, 8 and mouse TLR-9) were studied by performing docking analyses using PatchDock server in which the binding affinities of the vaccine (ligand) with various TLRs (receptors) were evaluated. Top ten models generated by PatchDock for each docking were considered among which the best ranking models with highest docking scores were selected for each case. The scores represented the shape and electrostatic complementarities between the receptors and ligand. The best complementarity was found in case of human TLR-2 (PDB id: 6NIG) with a score of 18714; surface area 3963.10 and atomic contact energy (ACE) 217.99 (Fig. 6a). The vaccine also exhibited efficient binding with TLR-5 (PDB id: 3J0A) and TLR-9 (PDB id: 3WPF) with a docking score of 18258; surface area 2879.80 and ACE of -666.94 and docking score 18164; surface area 4090.70 and ACE of 495.74 respectively. Docking with other TLRs showed comparatively weaker binding affinities as reflected by lower docking scores generated by PatchDock (16708 for TLR-4(PDB id: 4G8A), 17800 for TLR-8(PDB id: 4QC0). Thus, our vaccine construct was able to bind with different immune receptors with varying binding affinities. This would be highly desirable for an effective vaccine in terms of inducing a proper and heightened immunity to prevent leishmania infections. Further, FireDock web tool was employed for refinement and rescoring of the best complex (TLR-2-vaccine complex) which refined the same with a minimum global energy of -19.01; Van Der Waals associations (-25.74); atomic contact energy (-0.03) and binding free energy (-24.26). These results corroborate that fact that the vaccine peptide interacts most efficiently with human TLR-2 and thereby, form a fairly stable complex that might eventuate a marked immune response against leishmania parasites. In addition to this we further investigate docking between TLR-2 & vaccine using Cluspro web server, top model (Fig. 6b) is selected based on 5D rotational FFT and CAPRI (Critical assessment of predicted interactions) with balanced weighted score − 1192.5 having largest cluster size of 42 members.

3.15 Codon optimization and in silico cloning for expression of the vaccine protein:

With an attempt to overcome the problem of codon bias, JCat tool was utilized for codon optimization that returned the input amino acid sequence of the vaccine peptide in the form of a reverse translated cDNA sequence with an optimized CAI value of 1. This value lies within the accepted CAI value range of 0.8-1.0 (Ali et al. 2017) ensuring high rate of expression of the gene in the host organism, E.coli K12. The higher is the value of CAI, the better will be the expression rate of the foreign gene. In addition, the GC content of the improved sequence was found to be 51.85% indicating a good codon optimization. For better gene expression, 35%-70% GC contents is required (Khan et al. 2019). The restriction sites for the enzymes XhoI and NdeI were created at the 5^/ and 3^/ ends of the sequence followed by insertion of the modified sequence into pET28 a(+) vector. The cloned insert, represented in red color lies in between the specified restriction sites in the vector. A stretch of six histidine residues is also located at both the ends of the cloned sequence to facilitate the purification process after the production of the recombinant vaccine protein (Khatoon et al. 2017; Pandey et al. 2017). The size of the final plasmid vector construct with the inserted gene fragment was about 7 kilo bases. The mfold web server predicted a total number of fifty folding conformations of the mRNA among which the top structure was chosen due to its lowest energy score (ΔG = -492.50). This particular folded conformation was predicted to contain a single external loop and four stacks with minimum free energies indicating the stability of the conformation. The loop was formed with 21 single-stranded bases and 2 closing helices. The energy dot plot calculated the optimal energy (δG) to be -494.4 kcal/mol. The different colored dots in the plot refer to the different energy-grades of a particular base pair present in the mRNA considering the δG calculated by the server (Zuker, 2003). According to the plot, the free energies of most of the base pairs were found to be within − 487.5 < δG <= -485.8 kcal/mol range. Conclusively, a thermodynamically-favored and energetically stable secondary structure of the mRNA was predicted. All the folding constraints were set at their default values and the RNA was considered to be linear in nature.

3.16 Immune simulation:

In-silico immune simulation using C-immsim web server predicted antigenic count reaching up-to 6.9×10⁵ / ml within two days of initial dose of injection & antigen level dropped to 0 on fourth day onwards. This increment of antigen count is correlated with remarkably increase in the concentration of both IgM + IgG level and IgM level alone during early infection period & reaching peak level approximately within 12day of infection cycle on an arbitrary scale. Other subtypes of IgG1 + IgG2 level & IgG1 increase markedly as a result of secondary immune response (Fig. 8a). Cytokine profiles shows inflammatory mediators like INF-γ reaching level upto 4×10⁶ ng/ ml on 12th day while other cytokine levels remains quite low (Fig. 8b). Cell mediated immunological responses suggested by presence of active & resting Cytotoxic-T lymphocytes increases after 2–3 days after initial injection & presence of proliferative memory helper-T cell was also prominent (Fig. 8c & d). Total B cell population also increases within 2–3 days along with generation of memory B-cell after one week of administration of primary dose (Fig. 8e).

3.17 Interpretation of molecular dynamic simulation:

Best docking score achieved by TLR2-vaccine complex was further analyzed by NMA assisted molecular dynamic simulation. Deformability of main chain represented by high hinges (Fig. 9a) which means high deformability. B-factor (Fig. 9b) was calculated by NMA which signify protein flexibility, based on atomic displacement parameters (Sun et al. 2019). The eigenvalue (Fig. 9c) delineate about required energy to deform the structure. TLR2-vaccine complex shows eigenvalue of about 7.778×10^− 7. Variance (Fig. 9d) is inversely related to eigenvalue. Covariance matrix (Fig. 9e) illustrate whether pairs of residues are correlated (red), uncorrelated (white) or anti-correlated (blue). The elastic network model (Fig. 9f) computes springs between corresponding pairs of atom, dots in graph represented by color gradient indicating rigidity or flexibility of spring, the more darker color reflects more stiffer spring. Thus results of molecular dynamic simulation suggested that our proposed peptide vaccine is stable.