key: cord-0068897-ymsm09on
authors: Allemailem, Khaled S.
title: A Comprehensive Computer Aided Vaccine Design Approach to Propose a Multi-Epitopes Subunit Vaccine against Genus Klebsiella Using Pan-Genomics, Reverse Vaccinology, and Biophysical Techniques
date: 2021-09-27
journal: Vaccines (Basel)
DOI: 10.3390/vaccines9101087
sha: 542cb9e9a3aa242e5f4a4b06db830b709973761b
doc_id: 68897
cord_uid: ymsm09on

Klebsiella is a genus of nosocomial bacterial pathogens and is placed in the most critical list of World Health Organization (WHO) for development of novel therapeutics. The pathogens of the genus are associated with high mortality and morbidity. Owing to their strong resistance profile against different classes of antibiotics and nonavailability of a licensed vaccine, urgent efforts are required to develop a novel vaccine candidate that can tackle all pathogenic species of the Klebsiella genus. The present study aims to design a broad-spectrum vaccine against all species of the Klebsiella genus with objectives to identify the core proteome of pathogen species, prioritize potential core vaccine proteins, analyze immunoinformatics of the vaccine proteins, construct a multi-epitopes vaccine, and provide its biophysical analysis. Herein, we investigated all reference species of the genus to reveal their core proteome. The core proteins were then subjected to multiple reverse vaccinology checks that are mandatory for the prioritization of potential vaccine candidates. Two proteins (TonB-dependent siderophore receptor and siderophore enterobactin receptor FepA) were found to fulfill all vaccine parameters. Both these proteins harbor several potent B-cell-derived T-cell epitopes that are antigenic, nonallergic, nontoxic, virulent, water soluble, IFN-γ producer, and efficient binder of DRB*0101 allele. The selected epitopes were modeled into a multi-epitope peptide comprising linkers and Cholera Toxin B adjuvant. For docking with innate immune and MHC receptors and afterward molecular dynamics simulations and binding free energy analysis, the vaccine structure was modeled for tertiary structure and refined for structural errors. To assess the binding affinity and presentation of the designed vaccine construct, binding mode and interactions analysis were performed using molecular docking and molecular dynamics simulation techniques. These biophysical approaches illustrated the vaccine as a good binder to the immune receptors and revealed robust interactions energies. The vaccine sequence was further translated to nucleotide sequence and cloned into an appropriate vector for expressing it at high rate in Escherichia coli K12 strain. In addition, the vaccine was illustrated to generate a good level of primary, secondary, and tertiary immune responses, proving good immunogenicity of the vaccine. Based on the reported results, the vaccine can be a good candidate to be evaluated for effectiveness in wet laboratory validation studies.

The genus Klebsiella is ubiquitous in nature and comprises Gram-negative, rod shaped, and oxidase negative bacteria, and has a prominent polysaccharide capsule [1] . The genus is named after Edwin Klebs, a German-Swiss microbiologist. Species of this genus are part of animal and humans normal flora found in the mouth, nose, and intestines [2] . They are 0.5 to 5.0 µm long and 0.3 to 1.5 µm wide [3] . From a medical perspective, the polysaccharide and human self-antigens [21] . Genomics based methods predicted 90 previously unknown antigens. Meningococcus B and 29 out of them were found to stimulate antibodies that were able to kill the bacteria. In addition, the reverse vaccinology technique is used for development of protein-based vaccines against Streptococcus pneumoniae and Staphylococcus [21] . This technique can be integrated with several other approaches such as pang-genomics, immunoinformatics, and different biophysical analysis to construct a novel multi-epitope peptide [25, 26] . The core vaccine proteins were identified using pan-proteome analysis and used for broad spectrum epitopes prioritization [23] . Further, it was ensured that the selected epitopes are vital for the pathogen survival that are not human homologs, and are exposed to the host immune system for efficient recognition and processing. Ideally, conserved epitopes that were found to generate specific B-and T-cell responses were selected to enable indication of pathogen specific immune responses.

The methodology of the study can be explained in a stepwise fashion and is illustrated in Figure 1 .

Vaccines 2021, 9, x 3 of 25 pathogen that successfully was addressed by reverse vaccinology. The bacterium surface proteins are extremely variable and homology exists between the pathogen capsular polysaccharide and human self-antigens [21] . Genomics based methods predicted 90 previously unknown antigens. Meningococcus B and 29 out of them were found to stimulate antibodies that were able to kill the bacteria. In addition, the reverse vaccinology technique is used for development of protein-based vaccines against Streptococcus pneumoniae and Staphylococcus [21] . This technique can be integrated with several other approaches such as pang-genomics, immunoinformatics, and different biophysical analysis to construct a novel multi-epitope peptide [25, 26] . The core vaccine proteins were identified using pan-proteome analysis and used for broad spectrum epitopes prioritization [23] . Further, it was ensured that the selected epitopes are vital for the pathogen survival that are not human homologs, and are exposed to the host immune system for efficient recognition and processing. Ideally, conserved epitopes that were found to generate specific B-and Tcell responses were selected to enable indication of pathogen specific immune responses.

The methodology of the study can be explained in a stepwise fashion and is illustrated in Figure 1 . 

Primarily, the reference proteomic data of all the seven species of Klebsiella were retrieved from the NCBI database [27] . These seven species of Klebsiella include K. aerogenes, K. oxytoca, K. michiganansis, K. pnemoniae, K. quasipneumoniae, K. grimontii, and K. variicola. The proteomes were subjected to pan-proteome analysis to retrieve the core sequence. The Bacterial Pan Genome Analysis (BPGA) tool [28] was used to perform pan-proteome analysis. The core sequence retrieved through BPGA was then clustered through CD-Hit [29] to remove redundant proteins from the core sequence. The threshold was set at 0.5, so all the sequences having 50% similarity were clustered together. The subcellular localization of selected proteins was predicted through PSORTb 3.0 subcellular localization prediction tool [30] . The surface proteins were then filtered using several parameters to prioritize potential vaccine candidates for the pathogen [31] [32] [33] [34] . 

Primarily, the reference proteomic data of all the seven species of Klebsiella were retrieved from the NCBI database [27] . These seven species of Klebsiella include K. aerogenes, K. oxytoca, K. michiganansis, K. pnemoniae, K. quasipneumoniae, K. grimontii, and K. variicola. The proteomes were subjected to pan-proteome analysis to retrieve the core sequence. The Bacterial Pan Genome Analysis (BPGA) tool [28] was used to perform pan-proteome analysis. The core sequence retrieved through BPGA was then clustered through CD-Hit [29] to remove redundant proteins from the core sequence. The threshold was set at 0.5, so all the sequences having 50% similarity were clustered together. The subcellular localization of selected proteins was predicted through PSORTb 3.0 subcellular localization prediction tool [30] . The surface proteins were then filtered using several parameters to prioritize potential vaccine candidates for the pathogen [31] [32] [33] [34] .

Virulent proteins are a primary target for vaccine candidates because they play a significant role in pathogenesis [35, 36] . The nonredundant set of proteins were used in BLASTp search against the core virulent factor database (VFDB) [37] and those having sequence identity ≥30% and bit score >100 were selected. Subsequently, the physiochemical properties of proteins were checked using ProtParam tool [38] . The proteins having molecular weight <110 kDa and thermostability index >40 are considered as ideal vaccine candidates as they can be easily purified [39, 40] . The proteins were then subjected to transmembrane helices prediction by using HMMTOP 2.0 [41] . The proteins predicted to show (1 or 0) outside transmembrane helices were selected [42] . The shortlisted proteins were then BLAST against reference human proteome, which has "taxonomic ID: 9606" in the BLASTp of NCBI and those with E-value < 1, bit score >100 and sequence identity >30% were discarded [43] . All the homologous proteins were discarded to avoid the host's autoimmune reactions. Antigenicity of the proteins was predicted next through VaxiJen 2.0 with the threshold of 0.4 [44] . All proteins depicting an antigenicity score >0.5 were selected as antigenic proteins and the other proteins were discarded. Thus, these proteins that are highly antigenic towards T-cell receptors and antibodies are potential vaccine targets. The next step was to predict the adhesive nature of selected protein candidates as proteins having adhesive property are effective meditators of prompting pathogenesis. The adhesive nature of proteins was checked using SPAAN server keeping all the settings at default [45] . The AllerTOP2.0 [46] was then used to check the allergenicity of proteins that had IC 50 values ≤ 100 nM [43] . All the allergen proteins were removed to avoid allergic reactions, leaving only nonallergen proteins. Another BLASTp search of NCBI against Lactobacillus rhamnosus, Lactobacillus casei, and Lactobacillus johnsonii was performed to check homology of proteins with the mentioned probiotic species [32] . All the proteins that were showing no similarity were selected for epitope mapping analysis.

Successively, the liner B-cells epitopes were then predicted by using Bepipred Linear Epitope Prediction 2.0 [47] , which is accessible at Immune Epitope Database (IEDB) with a threshold of 0.5 [48] . The resulting peptides were then used in T-cell epitopes mapping in IEDB T-cell epitopes prediction tools to predict subsequences that bind to alleles of major histocompatibility complex (MHC) class I and II alleles [49] . The peptides with a percentile score <10 are high affinity binders [40] . The MHCPred 2.0 analysis [50] was then executed to calculate the binding affinity of the predicted epitopes for DRB*0101 allele and only those with IC 50 values ≤ 100 nM were selected [43] . Antigenicity, allergenicity, and virulence of the common peptides finally selected were checked using VaxiJen [44] , AllerTOP 2.0 [51] , and VirulentPred [52] , respectively. The solubility of peptides was then checked using the Innovagen tool, which is used as a peptide solubility calculator. All the peptides showing good water solubility were selected. Subsequently, ToxinPred [53] was used to check the toxicity of antigenic nonallergen virulent epitopes [54] . All the nontoxic peptides were selected and subjected to check whether these epitopes can prompt IFN-gamma or not using the IFN epitope server [55] .

The conventional vaccines using whole organism or large proteins engender nonspecific immune responses as it carries unnecessary antigenic load [56] . The development of a peptide vaccine, as a substitute to conventional vaccines, generates highly specific immune responses as it depends on short peptides that are free from allergenic reactions [34, 57] . However, peptide vaccines are weakly immunogenic [58] . To overcome this problem, all the epitopes selected were fused using GPGPG linkers to induce an immune response sufficiently strong to combat pathogens. This multi-epitope vaccine construct was then joined to an adjuvant cholera toxin B (CTB) [59] using EAAAK linker. The computation of various chemical and physical parameters of the vaccine construct was calculated using the ExPASy Protparam tool [38] . It permits calculation of various parameters such as instability index, amino acid composition, theoretical pI, molecular weight, and grand average of hydropathicity (GRAVY). After this, a file was created in which all the epitopes were added with alleles in order to check worldwide population coverage by additionally using population coverage, which is accessible at Immune Epitope Database (IEDB) [60] . The designed multi-epitopes vaccine construct should show extensive human population coverage. The multi-epitopes vaccine construct was subjected to tertiary structure predictions using 3Dpro of SCRATCH protein predictor [61] . The structure refinement was done through GalaxyRefine of GalaxyWeb [62, 63] .

Blind molecular docking is executed to analyze and observe the binding affinity of multi-epitopes vaccine construct with immune receptors. If the designed vaccine construct is showing interactions with host receptors and immune cells, then good immune responses can be established. The blind docking of multi-epitopes vaccine construct was performed by choosing TLR4, MHC-I, and MHC-II as receptors in PatchDock server [64] . The 3D structure of TLR4, MHC-I, and MHC-II were retrieved from Protein Data Bank (PDB) using code 4G8A, 1I1Y, and 1KG0, respectively. Additionally, the PatchDock complexes were assembled through RMSD, which was set to default at 4.0 Å. The protein-protein docking output generated through PatchDock was further subjected for refinement by using a FireDock (Fast Interaction Refinement in Molecular Docking) [65] . The FireDock helps to provide an effective platform for producing refined PatchDock complexes. After this, all of the complexes that were showing lowest global energy were ranked at the top and subjected to further evaluation. UCSF Chimera 1.13.1 [66] was used to investigate intermolecular interactions among the immune receptors (TLR4, MHC-I, and MHC-II) and the vaccine construct.

The molecular dynamics simulation assay helps to analyze the stability and dynamics of all the docked complexes in three steps: (a) system preparation, (b) preprocessing, and (c) production phase. A 100 ns simulation run for the docked complexes was performed using AMBER20 [67] . The first phase of the simulation was to prepare the complexes parameters using antechamber module. To solvate the complexes into TIP3P solvation box, the leap module of AMBER was used and the input value for size was set at 12 Å. To study the intermolecular and intramolecular interactions in molecular dynamics assay, a force field of ff14SB was applied [68] . Na + ions as counter ions were incorporated into the system for neutralization of charges. In the second round of the preprocessing phase, energies of the systems were minimized: energy minimization of hydrogen atoms (500 cycles), water box (1000 cycles), complete system atoms (1000 cycles), and nonheavy atoms (300 cycles). Subsequently, using NVB ensemble, the heating process was initiated and temperature was maintained at 300 K. To maintain temperature and hydrogen bonds of the systems, the Langevin dynamics [69] and SHAKE algorithm [70] were used, respectively. Then, in the next step, the complexes pressure was maintained for 100 ps to achieve pressure equilibrium using NPT ensemble. The CPPTRAJ module [71] was used to evaluate simulation trajectories.

MMPBSA.py module [72] was employed to calculate the solvation and associated binding free energies produced as an outcome of interactions between the vaccine and all three receptors. The binding free energy was calculated on total of 100 frames taken from the molecular dynamics simulation's trajectories. The equation used to determine the total binding energy for all the docked complexes follows:

Vaccines 2021, 9, 1087 6 of 23

Consequently, the C-Immune Simulation server (or C-ImmSim) [73] was used to check the in silico immune profiling of the multi-epitopes vaccine construct. The algorithm used in this server is mainly based on position-specific matrix (PSSM) to evaluate the immunogenic potential of a given antigen. All the input parameters were set as default except number of steps, which were set to 1000 [74] .

The sequence of vaccine was back-translated to generate DNA sequence and then codon usage was adapted to E. coli K-12 strain to achieve high expression of the vaccine. This was done through Java Codon Adaptation Tool (JCat) server [75] . Expression of the construct was assessed by codon adaptation index (CAI) and GC-content. This was followed by in silico cloning of the vaccine in E. coli pET28a vector by Snapgene tool.

The stability of protein was increased by introduction of disulfide bonds [76] . Design 2.0 [77] was used to perform disulfide engineering.

Advances in sequencing technologies and the field of metagenomics resulted in a paradigm shift in microbial genomics, thus allowing comparison of large scale pangenomics studies. This further facilitates estimation of genomic diversity and determination of highly conserved core genomes. In vaccine designing, this core genome has broad spectrum applications to select the most conserved antigen for a broad spectrum vaccine design [28, 78] . A total number of 973 proteomes were retrieved from the NCBI genome database, the details of which are tabulated in Table S1 , and only reference proteome of each Klebsiella species was used for pan-proteome analysis because the reference proteome is completely annotated. The core-proteome of the referenced species contains~22,393 proteins with an average of 3635 proteins for each strain. The core-pan plot total genes and core genes for each analyzed genome are presented in Figure 2 . The core proteome was then filtered for several parameters that were set to prioritize potential vaccine candidates.

These parameters can be briefly discussed as nonredundancy check, subcellular localization, homology check against human host, antigenicity check, allergenicity check, virulent analysis, adhesion analysis, physicochemical characterization, and BLASTp check against human probiotic bacteria. First, redundancy check was performed that predicted 5864 proteins as nonredundant and ensured that only one copy of each protein is present in the core genome sequence file [79] . Furthermore, only 43 protein sequences were found secretory and exoproteome in nature. The surface-localized proteins are good candidates for vaccine design [80] . Virulent proteins have been reported to have an important role in the pathogenesis of the organism and play a vital role as a good vaccine candidate [34] . Virulent peptides can provoke immune response properly, so virulent protein analysis was performed in order to predict virulent proteins. Virulent analysis indicated that 19 proteins are virulent and these proteins were used for further analysis. The 19 virulent proteins are tabulated in Table 1 . Next, transmembrane helices check and physicochemical analysis was carried out that filtered 18 proteins as thermodynamically stable, have few transmembrane helices, and exhibit good physicochemical properties. To avoid cross reactivity and autoimmune reactions inside the host, all core proteins sequences were subjected to comparative homology analysis, among which 15 proteins were predicted as nonhomologous to humans. From 15 human nonhomologous proteins, 34 proteins were homologues to human normal intestinal flora, which were discarded from the study. Antigenic proteins are important for designing a vaccine as they are the prime component to stimulate host immune responses against any antigen. Five proteins were predicted as virulent and forwarded to additional analysis. Similarly, an allergenicity check was performed to discard proteins that are able to generate allergic reactions. As such proteins are discouraged during vaccine design, they were not picked and only nonallergic protein candidates were opted. Only two proteins were picked as nonallergenic and used in downward analysis. Adhesion analysis was important in part to ensure selection of proteins that are key in pathogen attachment to the host cells and significant for initiating infection pathway [81] . In the past, adhesive proteins were unveiled to be important candidates for vaccine development [82] . Herein, both the proteins were found to play adhesive role and were selected. Lastly, to avoid accidental inhibition of the good bacteria, the filtered proteins were also tested for homology against different probiotic species. Proteins with no significant hits were only selected. The BLASTp results concluded both proteins as nonhomologous to the host probiotic bacteria. These two proteins were TonB-dependent siderophore receptor and siderophore enterobactin receptor FepA. These parameters can be briefly discussed as nonredundancy check, subcellular localization, homology check against human host, antigenicity check, allergenicity check, virulent analysis, adhesion analysis, physicochemical characterization, and BLASTp check against human probiotic bacteria. First, redundancy check was performed that predicted 5864 proteins as nonredundant and ensured that only one copy of each protein is present in the core genome sequence file [79] . Furthermore, only 43 protein sequences were found secretory and exoproteome in nature. The surface-localized proteins are good candidates for vaccine design [80] . Virulent proteins have been reported to have an important role in the pathogenesis of the organism and play a vital role as a good vaccine candidate [34] . Virulent peptides can provoke immune response properly, so virulent protein analysis was performed in order to predict virulent proteins. Virulent analysis indicated that 19 proteins are virulent and these proteins were used for further analysis. The 19 virulent proteins are tabulated in Table 1 . Next, transmembrane helices check and physicochemical analysis was carried out that filtered 18 proteins as thermodynamically stable, have few transmembrane helices, and exhibit good physicochemical properties. To avoid cross reactivity and autoimmune reactions inside the host, all core proteins sequences were subjected to comparative homology analysis, among which 15 proteins were predicted as nonhomologous to humans. From 15 human nonhomologous proteins, 34 proteins were homologues to human normal intestinal flora, which were discarded from the study. Antigenic proteins are important for designing a vaccine as they are the 

In B and T cell epitope mapping phase, all the prioritized shortlisted proteins sequence were analyzed for epitope prediction. Consequently, three B cell epitopes were predicted for both proteins and predicted B cell epitopes were further screened for B cell derived T cell epitopes, as shown in Table 2 . Antigenicity of the predicted epitopes was checked. Only those epitopes were selected whose predicted antigenic prediction score was greater than 0.5 and were considered as good candidates for vaccine development. Three epitopes were predicted as good antigenic epitopes for TonB-dependent siderophore receptor while four epitopes were selected for siderophore enterobactin receptor FepA. The epitopes selected fulfill different good antigenic epitope properties, such as antigenic, soluble, nonallergenic, nontoxic, and virulent. The population coverage of the epitopes is 82.86%. Further details on the population coverage can be found in Table S2 . 

A multi-epitopes vaccine construct was engineered that comprises different epitopes from the prioritized vaccine proteins and fused in order to generate substantial immune responses. In total, seven epitopes (noted in Table 2 ) were selected based on their ability to fulfill several different parameters, as presented in Table 2 . The selected epitopes are capable of eliciting both humoral-and cell-mediated immunity as the epitopes contain sequences of both B and T cells. In addition, the epitopes are free from allergic sequence, which ensures avoiding allergic reactions. The epitopes are nontoxic and predicted to not cause any toxicity inside the cells. Similarly, the epitopes are water soluble and have high affinity for the most prevalent DRB*0101 alleles, thus providing them with an efficient ability to binding to the host immune receptors [40] . The epitopes are antigenic and vital in binding with the products of immune system. Further, the epitopes are virulent, which ensures the epitopes to activate infectious pathway and allow the host immune system to respond to the antigen. The epitopes were joined via GPGPG linkers that allow easy presentation of the epitopes to the host immunity and keep the epitopes separated [32] . The final epitopes peptide was then joined to an adjuvant molecule to further boost the immune simulation ability of the designed vaccine [83] . The cholera B subunit was considered as the adjuvant molecule herein as the adjuvant is safe to be used in humans and allows active stimulation of cytotoxic T cells, which is a key player in destroying foreign pathogens inside the human body [59, 84] . The designed vaccine molecule is 222 amino acid long, spanning epitopes from two vaccine proteins. The designed vaccine construct is schematically presented in Figure 3A . 

The secondary and 3D structure of the vaccine is presented in Figure 3A ,B. The 3D structure of the vaccine was modeled ab initio as no appropriate template structure was available for homology-based structure prediction. The plot analysis presented by Ramachandran has shown the residues of the vaccine are in highly preferred zones, with 84.5% in most favored regions, 12.3% in additionally allowed regions, and 2.8% in generously allowed regions ( Figure 3D ). The ERRAT score of the vaccine is 80.24 and the Verify-3D score is 81.98. From the secondary structure point of view, the vaccine has 36.9% of alpha helix, 5.4% 3-10 helix, and 57.7% of beta turns, gamma turns, and helix-helix interactions. Once the vaccine model was achieved, loop modeling was performed to obtain the most optimal structure. In this regard, the residues-range Gln158-Asp160, Ala161-Ser172, and Pro182-Gly199 loop regions were loop modeled using Galaxyloop. Next, the structure was refined for structural errors. The refined structures along with the initial structure are given in Table 3 . The top one structure was selected because of improved galaxy energy of −3648.56 kcal/mol. The structure also has improved residues in the Rama-favored region, i.e., 90.5%. The refined model also has improved rotamer and clash scores. 

The overall weight of the vaccine is 23.89 kDa. This small-size vaccine allows its easy purification during experimental studies. The vaccine construct was found to be positive for antigenicity score of 0.93. In addition, the vaccine construct was further checked for its solubility prediction that revealed good water solubility. The vaccine was further reported to be a nonallergen. Stability index value of the vaccine is 26.23 showing that the vaccine construct is stable. The GRAVY index (−0.546) highlighted that the vaccine hydrophilic. The pI value of the vaccine is 8.41 while its aliphatic index score is 69.10. The physicochemical properties of vaccine are presented in Figure 4 . solubility prediction that revealed good water solubility. The vaccine was further reported to be a nonallergen. Stability index value of the vaccine is 26.23 showing that the vaccine construct is stable. The GRAVY index (−0.546) highlighted that the vaccine hydrophilic. The pI value of the vaccine is 8.41 while its aliphatic index score is 69.10. The physicochemical properties of vaccine are presented in Figure 4 . 

Molecular docking approach was used to dock the vaccine construct with immune receptors in order to check its binding affinity and presentation ability to the host immune system for robust production of immune responses [85] . The PATCHDOCK solutions are given in Tables S3-S5, respectively. The top 10 best solutions of the vaccine with TLR4, MHC-I, and MHC-II are tabulated in Table 4 . 

Molecular docking approach was used to dock the vaccine construct with immune receptors in order to check its binding affinity and presentation ability to the host immune system for robust production of immune responses [85] . The PATCHDOCK solutions are given in Tables S3-S5, respectively. The top 10 best solutions of the vaccine with TLR4, MHC-I, and MHC-II are tabulated in Table 4 .

For each complex, the best solution in term of global energy was selected. For TLR4, solution 4 of PATCHDOCK was selected as the best complex as its global binding energy score is −16.01 kcal/mol, which can be split into −31.00 kcal/mol from attractive van der Waals energy, 8.99 kcal/mol from repulsive van der Waals energy, 10.72 kcal/mol from atomic energy contact, and −4.13 kcal/mol from hydrogen bond energy. The vaccine was disclosed to form multiple hydrogen bonds with the TLR4 receptor via interaction with Gln81, Thr106, Lys130, Lys153, Ser207, His229, Lys351, Ser352, Phe573, Asn575, Val602, Val604, and Glu605 ( Figure 5A) . Similarly, in the case of MHC-I, solution 7 was selected because of global energy of −2.22 kcal/mol, with contribution of −2.06 from attractive van der Waal energy, 0.04 kcal/mol from repulsive van der Waals energy, 0.98 kcal/mol from atomic contact energy, and −0.50 kcal/mol energy from hydrogen bond. The following residues of the MHC-I interact through hydrogen bond with Asp77, Leu81, Gly83, Tyr84, Asn86, Glu89, Arg97, Met98, Tyr99, His114, Gln115, Tyr116, Tyr118, Ile124, Lys127, Thr134, Thr143, Trp147, and Val152 ( Figure 5B ). Solution 7 was selected for MHC-II. The global energy of this complex is −39.32 kcal/mol. The vaccine interacts with Pro16, Asp17, Gln18, Gly20, Asp35, Lys39, Ala64, Leu66, Glu71, Ile72, Lys75, and Ala61. Asn118, Pro127, Thr129, and Glu134 ( Figure 5C ). residues of the MHC-I interact through hydrogen bond with Asp77, Leu81, Gly83, Tyr84 Asn86, Glu89, Arg97, Met98, Tyr99, His114, Gln115, Tyr116, Tyr118, Ile124, Lys127 Thr134, Thr143, Trp147, and Val152 ( Figure 5B ). Solution 7 was selected for MHC-II. The global energy of this complex is −39.32 kcal/mol. The vaccine interacts with Pro16, Asp17 Gln18, Gly20, Asp35, Lys39, Ala64, Leu66, Glu71, Ile72, Lys75, and Ala61. Asn118, Pro127 Thr129, and Glu134 ( Figure 5C ). 

Taking into account the effectivity of molecular dynamics simulation in validating complexes dynamics and stability, 200 ns long all atom MD simulations was performed for the best docked vaccine and immune receptors complex. This analysis was essential for understanding the valuable information regarding the system's dynamics and to shed light on key vaccine binding interactions with the receptors. The trajectories of MD simulations were evaluated for different structural parameters such as root mean square deviations (RMSD) [86] , root mean square fluctuations (RMSF) [87] , and radius of gyrations (RoG) [88] . All these analyses were carried out considering carbon alpha atoms of the systems. RMSD was calculated first for all systems to decipher the extent of deviations the complexes can experience, considering the initial complex conformation as a reference. Higher RMSD is the output of significant instability and can be correlated to conformational changes of the investigated system. As depicted in Figure 6A , the systems are undergoing consistent conformational changes as the time proceeds; however, the changes are not sudden and minor regarding the docked conformation. These changes are within the limits and are not affecting overall intermolecular binding and interactions. The inspection of the steady RMSD increase reveal large and complexity of the systems and the presence of large number of loops that is flexible in nature and more dynamic. The TLR4 vaccine complex was found to show high RMSD compared to the MHC-I and MHC-II vaccine complexes. The maximum RMSD reached by the systems were as follows: TLR4 vaccine complex (8 Å), MHC-I vaccine complex (5 Å), and MHC-II vaccine complex (4 Å). Further insights were gained by knowing residue stability of the docked site. For this purpose, the systems RMSDs were split as per residue. The majority of the interacting receptor residues with the vaccine molecule unraveled a highly stable nature, as can be interpreted by the lower RMSF value shown in Figure 6B . The terminal residues in contrast to the core residues were found to have more fluctuations, which are expected due to the flexible nature of the biomolecule terminals. Overall, the average RMSF of the systems is <3 Å, which indicates formation of highly stable complexes and good affinity of the vaccine molecule for the receptors. Lastly, conformational equilibrium of the complexes was evaluated by RoG analysis. Lower RoG value implies a tight packing of the complexes while a higher value implies and corresponds to a loose packing of system's atoms. The RoG analysis results were in line with that of RMSD ( Figure 6C ). The TLR4 vaccine complex was reported to produce a higher RoG value, which is obvious because of the heavy nature of the complex. The highest RoG for the TLR4 vaccine complex reaches 90 Å, while for the MHC-I and MHC-II vaccine complexes the maximum RoG value is 60 and 45 Å, respectively.

The MMGB/PBSA method is a popular technique and commonly employed to accurately predict the binding free energies of complexes [89] . This method is less computationally expensive and more productive than most molecular docking approaches. MMGB/PBSA infers net binding free energy, where negative binding free energy is an indication of high receptor-ligand affinity, while positive net binding energy demonstrates low-docked stability. It was estimated that the TLR4 and MHC-I vaccine complexes report that both gas phase and solvation energy highly dominate the intermolecular interactions. In particular, the van der Waals energy, followed by electrostatic energy, was favorable in complex formation. Likewise, the polar solvation energy showed significant domination while nonpolar solvation energy is also supporting complex formation. The net binding free energy of complexes greatly supports intermolecular affinity by securing energy values of −714. 

Codon optimization defines genetic engineering methodologies that optimize codon to assure a maximum level of targeted protein production and expression. Before codon optimization, the peptide sequence of the final vaccine construct was reverse transcribed to DNA sequence for codon optimization by using the java codon optimization tool (JCat tool), and the expression system was investigated by CAI (codon adaptation index) and GC content of construct. The CAI value of the vaccine is 0.96 and GC content is 52.4%. The CAI value is considered good and the vaccine sequence can be inferred to have good expression rate. Lastly, the optimized vaccine construct was cloned into pET-28a(+) and the expression vector, as shown in Figure 7 .

to DNA sequence for codon optimization by using the java codon optimization tool (JCat tool), and the expression system was investigated by CAI (codon adaptation index) and GC content of construct. The CAI value of the vaccine is 0.96 and GC content is 52.4%. The CAI value is considered good and the vaccine sequence can be inferred to have good expression rate. Lastly, the optimized vaccine construct was cloned into pET-28a(+) and the expression vector, as shown in Figure 7 . 

C-Immune simulation was conducted in order to predict putative immune responses of the model vaccine construct. Both primary and secondary responses with respect to days and major immune players are shown in Figure 8A . Robust level of IgM and IgG can be seen, and combined IgM and IgG ratios reached 25,000 counts/mL. Among the interleukins and cytokines, the IFN-g was found in high concentration in response to the antigen ( Figure 8B ). 

C-Immune simulation was conducted in order to predict putative immune responses of the model vaccine construct. Both primary and secondary responses with respect to days and major immune players are shown in Figure 8A . Robust level of IgM and IgG can be seen, and combined IgM and IgG ratios reached 25,000 counts/mL. Among the interleukins and cytokines, the IFN-g was found in high concentration in response to the antigen ( Figure 8B ). C-Immune simulation was conducted in order to predict putative immune responses of the model vaccine construct. Both primary and secondary responses with respect to days and major immune players are shown in Figure 8A . Robust level of IgM and IgG can be seen, and combined IgM and IgG ratios reached 25,000 counts/mL. Among the interleukins and cytokines, the IFN-g was found in high concentration in response to the antigen ( Figure 8B ). 

Disulfide engineering is considered an important tool and used extensively in improving protein stability and modifying functional characteristics. The designed vaccine construct was subjected to disulfide engineering to highlight less stable regions of the vaccine and improve it for future applications. Concerning the nine residue pairs (Lys44-Arg56, Gln70-Ser76, Leu98-Lys102, Val103-Ile120, Asn144-Thr149, Gly155-Asp163, Gln158-Asp160, Ala161-Ser172, and Pro182-Gly199), the binding energy is an important contributor to vaccine instability and was high for all the residues (>2 kcal/mol). The wild and mutant structure of the vaccine is presented in Figure 9 . 

Disulfide engineering is considered an important tool and used extensively in improving protein stability and modifying functional characteristics. The designed vaccine construct was subjected to disulfide engineering to highlight less stable regions of the vaccine and improve it for future applications. Concerning the nine residue pairs (Lys44-Arg56, Gln70-Ser76, Leu98-Lys102, Val103-Ile120, Asn144-Thr149, Gly155-Asp163, Gln158-Asp160, Ala161-Ser172, and Pro182-Gly199), the binding energy is an important contributor to vaccine instability and was high for all the residues (>2 kcal/mol). The wild and mutant structure of the vaccine is presented in Figure 9 . 

Klebsiella species are responsible for variety of infections, including life threatening respiratory, bloodstream, and liver infections. Since no FDA licensed vaccine is available for preventing Klebsiella infections, there is an urgent need to propose novel vaccine candidates to be experimentally tested for immune protection efficacy. The traditional vaccine design strategy is expensive and very slow. Computational approaches in this regard can accelerate vaccine development process by providing ideal vaccine candidates against Klebsiella. In this present study, we applied a computational vaccine design approach on all reference genome of Klebsiella genus species and as such identified two highly conserved and potential vaccine candidates that fulfilled several vaccine candidacy parameters. These proteins are TonB-dependent siderophore receptor and siderophore enterobactin receptor FepA. Both proteins were subjected to epitope mapping phase where only seven epitopes were highly ranked and used in a multi-epitope peptide vaccine construction. Further, the designed vaccine ensemble revealed robust interaction energy and has a stable binding conformation with the tested immune receptor. The application of computer aided vaccine design can significantly reduce the vaccine development cost and in less time. The process can also aid in identifying the correct vaccine candidates for experimental evaluations. The results of the study are promising and must be evaluated experimentally to validate biological effectiveness of the designed vaccine construct. 

The data presented in this study are available within the article.

Klebsiella pneumoniae and Pseudomonas aeruginosa

Synthesis of Silver Nanoparticles from Mirabilis Jalapa and Evaluation of Antioxidant Activity

as Nosocomial Pathogens: Epidemiology, Taxonomy, Typing Methods, and Pathogenicity Factors

Healthcare-Associated Infections

The role of Klebsiella in Crohn's disease with a potential for the use of antimicrobial measures

The emerging problems of Klebsiella pneumoniae infections: Carbapenem resistance and biofilm formation

Multidrug resistant gram-negative bacteria in community-acquired pneumonia

Antibiofilm Peptides Increase the Susceptibility of Carbapenemase-Producing Klebsiella pneumoniae Clinical Isolates to β-Lactam Antibiotics

Extended spectrum β-lactamases (ESBL)-An emerging threat to clinical therapeutics

Treatment of ESBL-producing Klebsiella pneumoniae bacteraemia with carbapenems or flomoxef: A retrospective study and laboratory analysis of the isolates

Mechanisms of Resistance and Clinical Relevance of Resistance to β-Lactams, Glycopeptides, and Fluoroquinolones

Genetics and Worldwide Spread of New Delhi Metallo-β-lactamase (NDM): A threat to public health

Detection of New Delhi Metallo-Beta-Lactamase-1 (NDM-1) in carbapenem-resistant Klebsiella pneumoniae isolated from a university hospital in Iran

The Epidemiology of Carbapenem-Resistant Enterobacteriaceae: The Impact and Evolution of a Global Menace

Carbapenem-non-susceptible Enterobacteriaceae in Europe: Conclusions from a meeting of national experts

Erratum for Diago-Navarro et al. Novel, Broadly Reactive Anticapsular Antibodies against Carbapenem-Resistant Klebsiella pneumoniae Protect from Infection

Choice of therapeutic interventions and outcomes for the treatment of infections caused by multidrug-resistant gram-negative pathogens: A systematic review

Colistin resistance in Klebsiella pneumoniae

Progress towards the development of Klebsiella vaccines

Vaccinology in the genome era

Pulling the Brakes on Fast and Furious Multiple Drug-Resistant (MDR) Bacteria

Pangenome-reverse vaccinology approach for identifications of potential vaccine candidates in microbial pangenome

Reverse Vaccinology, and Bacterial Pathogenesis. Cold Spring Harb

Pan-vaccinomics approach towards a universal vaccine candidate against WHO priority pathogens to address growing global antibiotic resistance

Development of a Novel Multi-Epitope Vaccine against Crimean-Congo Hemorrhagic Fever Virus: An Integrated Reverse Vaccinology

Database Resources of the National Center for Biotechnology Information

BPGA-an ultra-fast pan-proteome analysis pipeline

Accelerated for clustering the next-generation sequencing data

PSORTb 3.0: Improved protein subcellular localization prediction with refined localization subcategories and predictive capabilities for all prokaryotes

Integrated Core Proteomics, Subtractive Proteomics, and Immunoinformatics Investigation to Unveil a Potential Multi-Epitope Vaccine against Schistosomiasis

Combating tigecycline resistant Acinetobacter baumannii: A leap forward towards multi-epitope based vaccine discovery

Towards a peptide-based vaccine against Shigella sonnei: A subtractive reverse vaccinology based approach

Designing multi-epitope vaccine against Staphylococcus aureus by employing subtractive proteomics, reverse vaccinology and immuno-informatics approaches

An In Silico Identification of Common Putative Vaccine Candidates against Treponema pallidum: A Reverse Vaccinology and Subtractive Genomics Based Approach

Immuno-Informatics Analysis of Pakistan-Based HCV Subtype-3a for Chimeric Polypeptide Vaccine Design

VFDB: A reference database for bacterial virulence factors

ExPASy-ProtParam Tool

Design of a Novel Multi Epitope-Based Vaccine for Pandemic Coronavirus Disease (COVID-19) by Vaccinomics and Probable Prevention Strategy against Avenging Zoonotics

Pangenome and immuno-proteomics analysis of Acinetobacter baumannii strains revealed the core peptide vaccine targets

The HMMTOP transmembrane topology prediction server

VacSol: A high throughput in silico pipeline to predict potential therapeutic targets in prokaryotic pathogens using subtractive reverse vaccinology

Identification of putative vaccine candidates against Helicobacter pylori exploiting exoproteome and secretome: A reverse vaccinology based approach

A server for prediction of protective antigens, tumour antigens and subunit vaccines

A software program for prediction of adhesins and adhesin-like proteins using neural networks

AllerTOP-A server for in silico prediction of allergens

BepiPred-2.0: Improving sequence-based B-cell epitope prediction using conformational epitopes

The Immune Epitope Database (IEDB): 2018 update

Major Histocompatibility Complex (MHC) Class I and MHC Class II Proteins: Conformational Plasticity in Antigen Presentation

MHCPred: A server for quantitative prediction of peptide-MHC binding

AllerTOP v. 2-A server for in silico prediction of allergens

VirulentPred: A SVM based prediction method for virulent proteins in bacterial pathogens

In silico approach for predicting toxicity of peptides and proteins

Peptide Toxicity Prediction

Designing of interferon-gamma inducing MHC class-II binders

Peptide-Based Vaccines: Current Progress and Future Challenges

Design and development of synthetic peptide vaccines: Past, present and future

Peptide-based immunotherapeutics and vaccines

Cholera Toxin Subunit B as Adjuvant--An Accelerator in Protective Immunity and a Break in Autoimmunity

Predicting population coverage of T-cell epitope-based diagnostics and vaccines

SCRATCH: A protein structure and structural feature prediction server

Protein structure refinement driven by side-chain repacking

Galaxy: A platform for interactive large-scale genome analysis

PatchDock and SymmDock: Servers for rigid and symmetric docking

A web server for fast interaction refinement in molecular docking

UCSF Chimera-A visualization system for exploratory research and analysis

Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB

Langevin stabilization of molecular dynamics

A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations

Software for processing and analysis of molecular dynamics trajectory data

py: An Efficient Program for End-State Free Energy Calculations

An information theoretic approach to interacting multiple model estimation

JCat: A novel tool to adapt codon usage of a target gene to its potential expression host

Protein disulfide engineering

Disulfide by Design 2.0: A web-based tool for disulfide engineering in proteins

Reverse vaccinology assisted designing of multiepitope-based subunit vaccine against SARS-CoV-2

Identification of plausible drug targets by investigating the druggable genome of MDR Staphylococcus epidermidis

Exoproteome and Secretome Derived Broad Spectrum Novel Drug and Vaccine Candidates in Vibrio cholerae Targeted by Piper betel Derived Compounds

Development of an anti-adhesive vaccine for Pseudomonas aeruginosa targeting the C-terminal region of the pilin structural protein

Adhesins as Targets for Vaccine Development

Application of built-in adjuvants for epitope-based vaccines

Multiepitope-Based Subunit Vaccine Design and Evaluation against Respiratory Syncytial Virus Using Reverse Vaccinology Approach

Peptide vaccine against chikungunya virus: Immuno-informatics combined with molecular docking approach

Binding mode analysis, dynamic simulation and binding free energy calculations of the MurF ligase from Acinetobacter baumannii

Comparative subtractive proteomics based ranking for antibiotic targets against the dirtiest superbug: Acinetobacter baumannii

Radius of gyration as an indicator of protein structure compactness

The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities

The researcher would like to thank the Deanship of Scientific Research, Qassim University for funding the publication of this project.

The author declares no conflict of interest.