Screening of secretory amastigote proteins from experimental studies
Following literature screening, three out of 28 PubMed indexed experimental studies on L. donovani proteome were considered for the study. Two of the studies compared proteomic abundance of promastigotes with that of amastigotes as Leishmania adapts to the changes in conditions resembling the host [57, 58], while the other one compared splenic amastigotes to axenic amastigotes . A total of 118 out of 134 proteins, which had a relative increase of at least 1.5-fold or were reported exclusively in the protein profile of amastigotes and/or splenic amastigotes, were found to have an identity of 90% or above for absolute query coverage with L. donovani proteins of similar functional annotations. After cross-matching of this group of proteins to 151 L. donovani secretory proteins revealed experimentally , 16 proteins were found to have an identity percentage and query coverage of 96.75 ± 1.1% and 99.44 ± 1.13%, respectively. Based on the presence of classical or non-classical secretion signal sequences along with minimum (no more than one) transmembrane helices, one common and an additional 17 secretory proteins were included to the pool. Among 33 amastigote-associated potential secretory proteins, 26 were selected based on their antigenicity probability scores of ≥ 0.5 as predicted by both ANTIGENpro and VaxiJen (Table 1, Additional file 2: Data S1).
Multi-epitope subunit L. donovani vaccine: construction and properties
A total of 79 CTL 9-mer epitopes were initially screened in NetCTL. Among them, only nine epitopes from six proteins were predicted to be non-self, highly immunogenic and high-to-moderate TAP-transporter binder non-B cell epitopes. These epitopes covered theoretically, an average of 66.46 ± 7.88% and a cumulative of 98.57% of the world population. Similarly, HTL 15-mer epitopes were screened to ensure both affinity and coverage. Fourteen selected HTL epitopes from eight proteins were finally predicted to be non-self, highly immunogenic non-B cell epitopes, with a mean theoretical coverage of 96.62 ± 1.35% and a cumulative of 99.52% of the world population. All the CTL and HTL epitopes except for H2-10 and H2-13 were conserved (100% identical) in L. infantum, whereas, less conservancy was found in representative proteins of L. major (13/23) and L. mexicana (10/23). Properties of individual CTL and HTL epitopes are given in Table 2 and Table 3, respectively.
The construct of 397 amino acid residues comprised of the 9 CTL and 14 HTL epitopes, with AAY and GPGPG linkers added in the intra-epitopic positions of CTL and HTL epitopes, respectively. It preceded in N-terminal by TLR4 peptide adjuvant, APPHALS, linked by EAAK linker to the vaccine. The selected rearranged model had the antigenicity score of 0.8 calculated by ANTIGENpro, and 0.74 (bacteria model) and 0.65 (parasite model) by VaxiJen. Furthermore, the construct was found to be non-allergenic for human use. When re-analyzed by the screening tools, all the original CTL and HTL epitopes were found consistent with the pre-screening immunogenicity, cleavage and TAP binding properties in the rearranged model. On the other hand, the arrangement of the construct resulted in generation of only three regions (15-mer overlapping) of IL-10 inducing epitopes and three non-specific CTL epitopes (9-mer) comparable to the potency of target epitopes (Additional file 3: Figure S1).
Prediction on IFN-γ induction capacity revealed a total of 117 epitopes (15-mer) with positive scores. This prediction was consistent with the C-ImmSim simulated immune response in terms of high IFN-γ production after hypothetical immunization (three doses) in a population characterized by a combination of frequent and VL susceptible HLA alleles [59–61]. Since the hypothetical cytokine levels in simulated immune response represent only the outcome of algorithmically set dynamic cellular interactions for a defined time period after antigen priming , it was important to evaluate whether the simulation module can respond dynamically to different constructs . Therefore, we simulated immune responses of two additional peptide vaccine candidates: peptide 1 (L. infantum derived fusion peptide ) and peptide 2 (L. donovani GP63 derived peptide ), which were experimentally found to exhibit varying cytokine response in comparison to soluble Leishmania antigen (SLA). Substantial difference was observed in terms of immunosuppressive IL-10 and TGF-β induction capacity between peptide 1 and peptide 2; however, determination of statistical significance was not possible in the simulation module. Nevertheless, the outcome can be considered consistent with the general trend of in vitro immune response (compared to SLA), with peptide 2 being more prominent IL-10 inducer compared to SLA as reported in . In terms of cytokine induction potential, simulation outcome of our designed construct conformed more closely to that of peptide 1, which did not induce IL-10 level higher than that by SLA in vitro  (Fig. 2).
Tertiary structure of the chimeric protein: physiological properties
Since there was no significant template hit, the choice of 3D model among outputs generated by I-TASSER was based on: (i) cluster size of model replicas; (ii) frequency of model in simulation trajectory; and (iii) C-score. The selected model has the highest C-score of -1.56 which is close to the I-TASSER recommended score (-1.5) for accuracy, and has the highest frequency in the top cluster by size. After further refinement of the protein topology, the PROCHECK  server returned a G-score of -0.04, which indicates that the backbone and side chain of the model correspond to high-probability stereochemical conformations. The model scored 1.73 in X-ray resolution scale by MolProbity , with no poor rotamers and bad bonds, negligible all atom steric overlaps (0.5%) and an increase in Ramachandran-favored residue number from 79.2% (unrefined) to 92.4% (refined) with a subsequent decrease in outliers (Additional file 4: Figure S2). The vaccine construct has a molecular weight of 42.1 kDa, with a basic nature (isoelectric point: 9.16). The score obtained for instability index was 27.26, which implied the stable nature of the vaccine in vitro. The estimated value of aliphatic index was 75.39 which indicated its thermo-stability. The folded structure has a melting temperature of 73.9 °C and folding free energy of -17.7 kcal/mol at neutral pH in humans. Additionally, this model was found to have substantial solubility with a score of 0.38 in folded state in contrast to the unfolded intrinsic score of -3.06, which suggests that hydrophobic residues in this model tend to ideally form the stable core leaving hydrophilic residues much on the solvent accessible surface. The half-life of the construct in mammalian reticulocytes was estimated as 4.4 h in vitro, compared with 20 h and 10 h in yeasts and Escherichia coli in vivo, respectively.
In terms of chimera-specific B cell response, Bepipred predicted six B cell epitopes of 8–12 residues in length above the threshold score, while BCPREDS predicted 11 non-overlapping and linear 20-mer B cell epitopes with specificity scores > 0.99. Residues in those linear epitopes accounted for 41% residues of the 08 non-overlapping conformational epitopes (Fig. 3, Additional file 5: Table S1).
This sequence was used to generate in silico cloning model for E. coli (K12) expression. After optimization of the codon, the codon adaptation index (CAI) value of the chimera was 0.98, while the GC content was 56.09%. For insertion into the E. coli pET28a(+) expression vector, two restriction sites for XhoI and NdeI enzymes were added in the 3'- and 5'-end, respectively, of the vaccine coding strand enclosed by 6-histidine residues at both ends (Additional file 6: Figure S3).
Molecular docking of vaccine in TLR4
Molecular docking of the vaccine construct with TLR4 in ClusPro 2.0 docking server generated 30 models ranked by cluster size of the representative pose. The selected docked complex had the largest cluster size (ClusPro recommended) with second-lowest binding energy score (-1282.3) among the top ten models. The chimeric construct seemed to occupy partially into the lateral concave surface, but not the convex surface, with strong hydrophobic interactions mostly with the beta-sheet adjacent residues at the C-terminal domain of TLR4 ectodomain (ECD) and also with its adapter protein, MD2, with support of several hydrogen bonds, thus establishing ligand mediated cross-link between TLR4 and MD2 (Fig. 4).
Molecular dynamics (MD) simulation of vaccine-TLR4 complex
Molecular dynamics simulation of the docked complex was performed by using OPLS_2005 force field. Using the Simulation Quality Analysis tool of the Desmond software, the mean potential energy for the complex was obtained as -6.4e5 kilocal/mol (Additional file 7: Figure S4). The radius of gyration (Rg) obtained for the docked complex showed that the mean distance in rotating complex from the center of mass is 4.31 nanometers (SD: 0.2 nanometers) about which the model becomes consistent after 4 ns (Fig. 5a). The number of intermolecular hydrogen bond (H-bond) between the side chains of vaccine protein and TLR4 initially fluctuated probably due to solvent effect before matching the trend of Rg in reaching steadiness after 4 ns. This suggests the role of H-bonds in the overall compactness of the complex (Fig. 5b). The trends of Rg and H-bond plots indicate that 6–8 strong H-bonds were persistent over simulation period between vaccine and TLR4, and this might be crucial for stable binding.
The root mean square deviation (RMSD) of the vaccine-TLR4 complex for backbone atoms over the simulation period was 4.0 Å (SD: 0.49 Å), while it was 3.2 Å (SD: 0.35 Å) for ligand-free TLR4 atoms (Fig. 5c), suggesting comparably higher (paired t-test: P < 0.0001) RMSD of the complex backbone. The root mean square fluctuation (RMSF) for side-chain atoms of vaccine-bound TLR4 (1.9 Å, SD: 0.7 Å, range: 0.8–7.0 Å) was higher (Wilcoxon matched-pairs test, P < 0.0001) than unbound TLR4 (1.6 Å, SD: 0.5 Å, range: 0.7–4.2 Å). The RMSF indicates overall less fluctuations for atoms interacting with vaccine residues, while atoms at vaccine unbound regions of the N-terminal and central domain underwent high fluctuations (Fig. 5d). Although the trends toward reaching convergence were very similar, higher RMSD value of the complex than the vaccine-unbound TLR4 indicates structural mobility in the complex due to vaccine interaction and this is likely attributable to the higher RMSD of vaccine protein along the MD simulation time. In congruence, rearrangement of several bonds between the vaccine and TLR4 was observed between pre-simulation and post-simulation models, while the total number of non-covalent bonds increased from 41 in pre-simulation model to 64 in post-simulation model (not shown). Visualization of the interacting residues also indicates that, in comparison to unbound (and also pre-simulation) structure, the post-simulation bonding rearrangement is coupled with increased number of H-bond at the C-terminal domains between TLR4 and TLR4* (second TLR4 ECD) (Additional file 8: Figure S5). This implies likely chance of positive interactions between the TLR4 monomers in physiological condition following vaccine interaction. Overall, the conformation of vaccine-bound receptor supports structural flexibility, which might be in favor of biological response of the receptor.
Immune simulation to predict secondary response
Hypothetical administration of three doses of vaccine construct four weeks apart with 1000 unit/dose was performed to simulate the immune response generated by immunization. In silico immune simulation plots hinted at antigenic recognition and subsequent response in terms of antibody production, and active as well as memory B cell and T cell generation in the population with a VL susceptible HLA profile after hypothetical immunization. The primary response to the proposed chimera can be characterized by a marked increase in chimera-specific IgM and IgG production. After subsequent doses, a corresponding decrease in antigen concentration indicates gradual increase in memory B cell production with persistence. Furthermore, expansion of CD4+ T lymphocytes with memory development following the initial dose was observed. CD8+ T lymphocytes response was also high for the susceptible population reaching at its peak after the second dose. Repeated exposure of 12 doses, on the other hand, did not seem to cause clonal expansion of any epitope-specific T cells- as indicated by the Simpson’s index D, which is inversely related to diversity (Additional file 9: Figure S6).