Screening of secretory amastigote proteins from experimental studies
Following a literature screening, 03 out of 28 pubmed indexed experimental studies on L. donovani proteome were considered for the study. Two of the studies compared proteomic abundance of promastiogtes with that of amastigotes as Leishmania adapts to the changes in conditions resembling host [49, 50], while the other one compared splenic amastigotes to axenic amastigotes [17]. A total of 118 out of 134 proteins that had fold increase of at least 1.5 times or found specifically in amastigotes/splenic amastigotes had 90% or above identity with L. donovani proteins for absolute query cover and with similar functional annotations. After cross-matching of this group of proteins to 151 L. donovani secretory proteins found experimentally [19], 16 proteins were found to have a percent identity and query cover of 96.75±1.1% and 99.44±1.13%, respectively. Prediction of classical and non-classical secretion signal sequences with minimum (no more than one) transmembrane helices revealed 01 common and an additional 17 secretory proteins. 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).
Multi-epitope subunit L. donovani vaccine: construct and properties
A total of 79 CTL 9-mer epitopes were initially screened in NetCTL. Among them only 9 epitopes from 6 proteins were predicted to be non-self, highly immunogenic and high-to-moderate TAP-transporter binder non B-cell epitopes, and covered theoretically, an average of 66.46±7.88% and a cumulative of 98.57% world population. Similarly, HTL 15-mer epitopes were screened to ensure both affinity and coverage. 14 selected HTL epitopes from 08 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 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 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, and 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. All the original CTL and HTL epitopes were consistent in immunogenicity, cleavage and TAP binding properties in the rearranged model, when re-analyzed by the screening tools. However, arrangement of the construct resulted in generation of three non-specific CTL epitopes (9-mer) comparable to the potency of target epitopes according to our original prediction method, and three regions (15-mer overlapping) of IL–10 inducing epitopes (Additional file 3: Figure S1).
IFN-γ epitopes
Prediction on IFN-γ induction capacity revealed a total of 117 epitopes (15-mer) with positive scores. This prediction was consistent with the C-ImmSim simulation predicted immune response in terms of high IFN-γ production following hypothetical immunization (three doses) in population having frequent and VL susceptible HLA alleles [51–53]. Since the hypothetical cytokine levels in simulated immune response represent only the outcome of dynamic cellular interactions following antigen priming and for a defined time period [37], it was thus important to evaluate whether the simulation module can respond dynamically to different constructs. Therefore, we simulated immune response of two additional peptide vaccine candidates, which were experimentally found to exhibit varying cytokine response in comparison to soluble leishmania antigen (SLA) (Fig. 2). 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 since it is not the purpose of simulation. Nevertheless, the outcome can be considered consistent with the general trend of in vitro immune response (compared to SLA) generated by the peptides: peptide 2 being more prominent IL–10 inducer compared to SLA as reported in [39]. In terms of cytokine induction potential, simulation outcome of our designed construct, however, conformed more closely to that of peptide 1, which did not produce more IL–10 than SLA did in vitro [38].
Tertiary structure of chimera: physiological properties
Since there was no significant template hit, the choice of 3D model among outputs generated by I-TASSER was based on- 1) cluster size of model replicas, 2) frequency of model in simulation trajectory, and 3) 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 (Figure 4), the PROCHECK [54] 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 [55], 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, showing the stable nature of vaccine in vitro. The estimated value of aliphatic index was 75.39, indicating 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 human. Besides, 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, suggesting 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 hours in vitro; while 20 and 10 hours in yeast and E. coli in vivo.
In terms of chimera specific B-cell response, Bepipred predicted 06 B-cell epitopes 8–12 residues in length above the threshold, 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 adaptive index (CAI) value of the chimera was 0.98, while the GC content was 56.09%. Two restriction sites, XhoI and NdeI were added in the 3’ and 5’, respectively, of vaccine coding strand enclosed by 6-histidine residues on 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-most 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 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, suggesting 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-bond were persistent over simulation period between vaccine and TLR4, and this might be crucial for binding stability.
The mean 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-pair 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 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 is likely attributable to the higher RMSD of vaccine protein along the MD simulation time. In congruence, rearrangement of several bonds between vaccine and TLR4 was observed between pre-simulation and post-simulation models, while 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 suggest post-simulation bonding rearrangement is coupled with increased number of H-bond at the C-terminal domains between TLR4 and TLR4* (second TLR4 ECD) in comparison to unbound (and also pre-simulation) structure (Additional file 8: Figure S5). This suggests likely chance of positive interactions between the TLR4 monomers in physiological condition following vaccine interaction. Overall, the conformation of vaccine bound receptor suggests structural flexibility, which might be in favor of biological response of receptor.
Immune simulation to predict secondary response
In silico immune simulation plots (Additional file 9: Figure S6) hinted on antigenic recognition and subsequent response in terms of antibody production and active as well as memory B cell and T cell generation in population with VL susceptible HLA profile following hypothetical immunization. Primary response to the proposed chimera can be characterized by 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 initial dose was observed. CD8 T lymphocytes response was also high for the susceptible population reaching at its peak after second dose. Repeated exposure of 12 doses on the other hand did not seem to cause any clonal expansion of epitope specific T-cells.