The reference proteomes of K. pneumoniae (strain ATCC 700721/MGH 78578) and P. aeruginosa (strain ATCC 15692/DSM 22644/CIP 104116/JCM 14847/LMG 12228/1C/PRS 101/PAO1) were downloaded from the UniProt webserver (https://www.uniprot.org/) under the proteome ID of UP000000265 and UP000002438, respectively. As mentioned in the introduction section, the current study aims to design an epitope-based vaccine through the filtration of protein candidates belonging to the outer membrane and iron uptake proteins. Therefore, we selected nine K. pneumoniae protein candidates namely FepA, FepB, FepC, FhuA, FhuF, FuR (iron uptake proteins), OmpA, OmpC, and OmpF (outer membrane proteins), and filtered them through their antigenicity score estimated by VaxiJen v2.0 (Doytchinova and Flower, 2007) with the cutoff score of 0.4 (the threshold value of bacterial antigenic proteins). The assessment of the antigenicity score revealed that there were 8 antigenic proteins, out of the selected 9 ones therefore we selected the top 2 proteins (one protein from each category) based on their antigenicity score where the final 2 protein candidates of K. pneumoniae were FepA and OmpF with antigenicity scores of 0.76 and 0.81 respectively. Moving to P. aeruginosa, we followed the same approach where six protein candidates namely FoxA, FpvA, HasR, HitA (iron uptake proteins), OprF, and OprH (outer membrane proteins) were filtered and 2 proteins (also one from each category) namely HasR and OprF with the antigenicity scores of 0.59 and 0.8 respectively were selected as our final candidates for P. aeruginosa.

The filtered 4 proteins from the previous step were uploaded to SignalP- 5.0 server (Almagro Armenteros et al., 2019) to predict the location of signal peptides. Following that, the mature polypeptides were analyzed for their T and B cell epitopes through the Immune Epitope Database (IEDB) (Dhanda et al., 2019). Firstly, we mapped CTLs for the protein candidates using the HLA allele reference set, which provided more than 97% in terms of population coverage (Weiskopf et al., 2013), and the NetMHCpan EL 4.1 prediction tool (that was recommended by the IEDB database). Secondly, we mapped for HTLs against the HLA reference set to cover more than 99% in terms of population coverage (Greenbaum et al., 2011) and used IEDB recommended 2.22 as a prediction method. Furthermore, HTL peptides were assessed for their ability to induce several cytokines such as IFN-gamma (Dhanda et al., 2013b), IL-4 (Dhanda et al., 2013a), IL-10 (Nagpal et al., 2017), IL-6, and IL-13 (Jain et al., 2022). The last analysis for HTLs and CTLs was the conservancy prediction where multiple sequence alignment against the corresponding proteins in other reference sequences was employed to validate the conservancy of the selected epitopes. The last set of epitopes; namely BCLs were finally estimated through IEBD using the BepiPred-2.0 prediction method (Jespersen et al., 2017). Following prediction, the estimated epitopes were filtered based on the consideration of several characteristics such as the number of reacting alleles (to achieve a high population coverage percentage), conservancy percentage, and antigenicity score.

Following the detection of B and T cell epitopes, we selected top-ranking single epitopes to initiate a multitope vaccine construct representing the epitopes of 4 proteins of our 2 bacteria of interest. Therefore, the best 2 candidates of CTL, HTL, and BCL epitopes per each protein were linked through GGGS, GPGPG, and KK amino acid linkers, respectively, in order to apply in vivo separation of the joined epitopes (Hajighahramani et al., 2017). In addition to the single epitopes with their linkers, PADRE sequence and β-defensin adjuvant were incorporated to complete the potential multitope vaccine construct. Finally, we assessed the multitope vaccine construct for its antigenicity, allergenicity, toxicity, and percentage of population coverage through VaxiJen v2.0, AllerTop v2.0 (Dimitrov et al., 2014), ToxinPred (Gupta et al., 2013), and IEDB (Bui et al., 2006) webservers.

In this stage, we utilized ProtParam, a tool available on the Expasy server (Gasteiger et al., 2005) SOLpro server (Magnan et al., 2009), and PSIPRED 4.0 (Buchan and Jones, 2019) web server in order to anticipate the physicochemical properties, the propensity upon overexpression in E. coli, and the protein secondary structure of the generated potential vaccine construct from the previous step.

After the design with physicochemical and immunological properties assessment of the multitope vaccine construct, we aimed to predict its 3D structure to be used for a docking study with the human immune receptor. For this purpose, the Robetta server (Kim et al., 2004) was employed. Robetta server utilizes a unique approach for protein structure prediction where if a confident match to a protein of known structure is found using BLAST, PSI-BLAST, FFAS03, or 3D-Jury, this protein is employed for the modeling process. Alternatively, if no match is found, the modeling process is performed through the de novo Rosetta fragment insertion method. Following that, we utilized the GalaxyRefine server (Heo et al., 2013) to refine the 3D protein structure estimated by Robetta and evaluated this refinement through the generated scores of Ramachandran plot analysis (Laskowski et al., 1993) and ProSA (Wiederstein and Sippl, 2007).

While continuous B cell epitopes, which were predicted in the above sections, are estimated through the primary amino acid sequence of a protein, another type of epitopes, which is conformational (or discontinuous) B cell epitopes, are predicted based on the 3D structure of the antigenic protein. For this purpose, the ElliPro Server (http://tools.iedb.org/ellipro) was utilized.

Prior to molecular docking simulations, proteins such as the human TLR2 (PDB; 2Z7X) and TLR4 (PDB: 3FXI) and constructed hybrid multi-epitope vaccine were independently prepared through partial charge assignment, 3D-protein protonation and/or removal of any bound ligand(s), crystalized solvent, and ionic metals/salts (Helal et al., 2020; Elhady et al., 2021; Elmaaty et al., 2022). For identifying the active and surrounding residues within the structures of the investigated proteins, the vaccine and TLR targets were submitted at Consensus Prediction Of interface Residues in Transient complexes (CPORT; https://alcazar.science.uu.nl/services/CPORT/) (de Vries and Bonvin, 2011). Docking the constructed vaccine on the TLRs was performed using an on-line docking server; ClusPro v2.0 (Boston and Brook Universities; https://cluspro.org/) (Kozakov et al., 2017; Porter et al., 2017). Relying on the Fast Fourier Transform correlation protocol, ClusPro predicted the vaccine/TLR4 complex through a multi-stage process; PIPER-based rigid docking, interaction energy-based conformation filtration, ranking based on clustering properties, and finally refinement through minimization (Kozakov et al., 2006). Interaction energy adopted by ClusPro included energy terms for van der Waals (EVDW = Eatt + Erep), electrostatic (Eelec), and pairwise structure-dependent potentials (EDARS) resulting from Decoys as reference state, however, lacked entropic energy terms (Chuang et al., 2008). It was suggested to utilize the cluster ranking, in terms of cluster populations, rather than the obtained ClusPro interaction energy for ranking and identifying the best-clustered structure complex (Zhang et al., 2022).

Evaluation of the best docking pose proceeded using the online PDBePISA v1.5.2 server tool (European Bioinformatics Institute/EMBL-EBI; https://www.ebi.ac.uk/msd-srv/prot_int/cgi-bin/piserver) for macromolecular interface analysis (Evgeny, 2010; Schlee et al., 2019). This tool provided descriptions for the sole and bound protein interface including interface residues, total solvent-accessible surface area (Å2), numbers/types of binding interactions, as well as the gained solvation-free energy (ΔiG; Kcal/mol), and its p-value (ΔiG p-value). The last two descriptors are indices for higher interface hydrophobicity/protein affinity (high negative ΔiG values) and to how many degrees the protein-protein interface can be interaction-specific (p < 0.5) (Krissinel and Henrick, 2007). Estimating the vaccine binding affinity for the predicted docked complex was further investigated using the Molecular Mechanics energy-guided Generalized Born and Surface Area (MM/GBSA; Kcal/mol) calculations implemented at the HawkDock server (Zhejiang University; http://cadd.zju.edu.cn/hawkdock/) (Weng et al., 2019). HawkDock MM/GBSA calculation permitted an estimation of energy term components including van der Waal, electrostatic, Generalized Born-predicted polar solvation free energy, and empirical model-predicted non-polar solvation contribution, besides dissection them down to the protein’s per-residue energy contributions (Hou et al., 2011; Zhong et al., 2020).

The best-docked vaccine-TLR complex was subjected to 100 ns all-atom molecular dynamics simulations under CHARMM36m forcefield and GROMACS program (Páll et al., 2015). Protein complexes were solvated at TIP3P cubic box under periodic boundary conditions maintaining a minimum distance of 10 Å between the protein atoms and box boundaries. The net charge of the system was neutralized via sufficient 0.15 M sodium and chloride ions. Systems were subjected to the steepest-descent minimization at 0.05 ns, followed by two-staged equilibration at standard thermo- and barostats (Berendsen-temp for NVT ensembles; 1 ns at 310 K followed by Parrinello-Rahmann barostat for NPT ensemble; 1 ns at 1 atm. and 310 K. Molecular dynamics were run for 100 ns under NPT ensemble and Particle-Mesh-Ewald algorithms for computing long-range electrostatic interaction (Andricioaei and Karplus, 2001). Trajectory analysis was performed using root-mean-square deviations (RMSDs; Å), radii of gyration (ROG; Å), solvent-accessible surface area (SASA; Ų), and RMS-fluctuations (RMSFs; Å) relying on the protein’s backbone alpha-carbon atoms. Conformational analysis and visualization of the simulated complexes at specified timeframes were done using PyMOL software (Schrödinger; v2.0.6, United States).

The final stage of our study was the prediction of the generated immune response upon the administration of the potential vaccine construct where The C-ImmSim server (Rapin et al., 2010) was employed for this purpose. In the estimation, we utilized the prime-booster-booster technique, an approach that was carried out via the injection of the potential vaccine three times at 4 weeks intervals. This approach has been followed to get a long-lasting immune response.

The last stage of the current study was the codon optimization for the designed potential vaccine where we employed the JCAT server for this purpose (Grote et al., 2005). Here, we selected E. coli k-12 strain as the expression organism as it is frequently used in gene cloning experiments (the first stage for wet lab validation of the current potential vaccine). The codon adaptation index (CAI), a value that is calculated by the server, gives an estimation for the constructed potential vaccine to be expressed in E. coli k-12.

Regarding K. pneumoniae, 9 protein candidates namely FepA, FepB, FepC, FhuA, FhuF, FuR (iron uptake proteins), OmpA, OmpC, and OmpF (outer membrane proteins) were filtered through their antigenicity score estimated by VaxiJen v2.0 with the cutoff score of 0.4 (the threshold value of bacterial antigenic proteins) and the final 2 protein candidates of K. pneumoniae were FepA and OmpF with antigenicity scores of 0.76 and 0.81 respectively (highest scores). Moving to P. aeruginosa, the same approach was followed where six protein candidates namely FoxA, FpvA, HasR, HitA (iron uptake proteins), OprF, and OprH (outer membrane proteins) were filtered and 2 proteins (also one from each category) namely HasR and OprF with the antigenicity scores of 0.59 and 0.8 respectively were selected as our final candidates for P. aeruginosa. Collectively, 4 proteins were selected as final targets for epitope prediction (2 proteins per microorganism).

At this stage, the 4 filtered protein candidates were mapped for their T and B cell epitopes. For T cell epitopes, the mapping of MHC-I epitopes resulted in a large number of epitopes for each protein with a percentile rank ranging from 0.01 to 100 where a small percentile rank represents a better binding in comparison to the large percentile rank score, and for this reason, we filtered only the epitopes with percentile score less than 1 (to guarantee for good binding) and the best candidates with a large antigenicity score were listed in Table 1 (that demonstrates the top five epitopes identified for each protein). In addition to that, the same 4 protein candidates were mapped for MHC-II epitopes and the output results were arranged based on the binding affinity. For this type of epitope, we selected only the top 1% of the results and selected the best candidates based on antigenicity score and the ability of the epitope to induce INF-γ, IL-4, IL-10, IL-6, and IL-13 cytokines. Table 2 combines the top five epitopes identified for each of the 4 selected proteins. Lastly, we mapped for B cell epitopes using the BepiPred-2.0 prediction method where the peptides with a length between 9:20 amino acids were submitted to VaxiJen for antigenicity score estimation (Table 3) and the ones that were highly antigenic were selected for the multitope vaccine construction.

As the next stage after epitope mapping, we planned to select the top 2 single epitopes (for FepA, OmpF, HasR, and OprF) from our final list of MHC-I, MHC-II, and B cell epitopes and after selecting the most promising B and T cell epitopes, we get a final count of 8 CTL, 8 HTL, and 8 BCL epitopes. These epitopes were linked together using GGGS, GPGPG, and KK linkers respectively. Following that, 2 important components, namely PADRE peptide and β-defensin, were incorporated to finalize the sequence of the multitope construct which lastly constituted 476 amino acids and sequenced as the following:

“EAAAKGIINTLQKYYCRVRGGRCAVLSCLPKEEQIGKCSTRGRKCCRRKKEAAAKAKFVAAWTLKAAAGGGSKQEPKKYNYGGGSQTNPNYILYGGGSRTNGVATYRGGGSDTSSDDTTYGGGSSVDGMRQNYGGGSRVKHSPVAYGGGSKVHGNLTSLGGGSKSKVKENSYGPGPGLVRWEFAPMQSLEFEGPGPGIIPEYTLNSTLSWQVGPGPGFFGLVDGLSFALQYQGPGPGAVVQYQFDFGLRPSIGPGPGVSQDDLVQMSPSVISGPGPGDLSTLRANYGLEFFYGPGPGGAGLKYYFTENFFAKGPGPGGVGLRPYVSAGLAHQKKGNSTSGQRGNNKKAPEHKDEGSTKRTKKDASNVEGSQTTKKKSNRSVDQKADGNGDKAKKDSASQQQTALRRVRLDIKKSGLGGDANGTKKNITNINSDSQGRQQMTKKSVGTDAYNQKLSERKKAKFVAAWTLKAAAGGGS”.

As mentioned in the methodology section, T cell epitopes mapping was performed against the HLA allele reference set to achieve a high percentage of population coverage, consequently, assessment of the percentage of population coverage revealed the scores of 92.56%, 99.98%, and 100% for world MHCI, MHCII, and combined population coverage respectively (Supplementary Figure S1). Prior to the tertiary structure prediction and analysis of the docking properties of the potential vaccine, it was essential to test and assess its physicochemical properties. The potential vaccine sequence was submitted to AllerTop, VaxiJen, and ToxinPred servers to assess for its allergenicity, antigenicity, and toxicity and the server’s output revealed that our potential vaccine is non-allergen, non-toxic, and antigenic with an antigenicity score of 1.38. Moreover, the potential vaccine sequence was analyzed for its solubility upon overexpression and obtained a score of 0.85 which was a satisfactory one (a score that exceeds 0.5 indicating solubility upon overexpression). Another essential analysis for the potential vaccine sequence was the assessment through Blastp where the results demonstrated that the submitted multitope sequence did not significantly resemble human protein sequences therefore this potential vaccine would not elicit autoimmune reactions during human usage. Furthermore, ProtParam online tools were employed to assess the multitope construct for its other physicochemical properties where the satisfactory output scores were summarized in Supplementary Table S1. Lastly, the secondary structure prediction for the multitope construct was performed using PSIPRED 4.0 webserver and the output predicted 28% helix, 13.4% strand, and 58.6% coil of the potential vaccine construct secondary structure (Supplementary Figure S2).

In order to assess the chemical properties of the potential multitope vaccine it was essential to predict its tertiary structure and for this purpose, we employed the Robetta webserver. The generated 3D model was then uploaded to the GlaxyRefine webserver where Ramachan plot analysis and ProSa Z-score were utilized to assess the quality of the original model generated by Robetta and the refined one. For the current study primary model (output of Robetta prediction), 85.3%, 13.4%, and 1.4% of residues were located in favored, allowed, and outlier regions, respectively. After refinement on the GlaxyRefine server, the refined model demonstrated 89.9%, 10%, and 1.1% of residues located in favored, allowed, and outlier regions, respectively with a Z-score of −4.59. The predicted tertiary structure of the multitope construct and the output of its structural validation are shown in Figure 2.

The generated 3D structure of the multitope construct was uploaded to EIIiPro webserver for the prediction of conformational B cell epitopes. The server estimated 9 discontinuous epitopes (Figure 3) with a score ranging between 0.578 and 0.977. Table 4 mentions the residues of each predicted epitope and the corresponding score.

The employed docking server illustrated significant binding of the constructed multiepitope vaccine at both TLR isotypes, TLR4 and TLR2, where these innate immunity receptors are highly reported for their important role within the host’s defense throughout infections by K. pneumonia (Wieland et al., 2011; Paczosa and Mecsas, 2016; Jeon et al., 2017; Dar et al., 2019; Allemailem, 2021). Selection of the relevant binding mode was guided by furnishing high docking scores as well as achieving contacts with the CPORT-suggested hot-spot residues for relevant protein-protein binding (Supplementary Figure S3). Binding modes were quite comparable at both TLRs, where the vaccine was anchored at the inner concave surface of the TLRs by the virtue of the latter extended conformation shaping as a paddle (Figure 4). Relevant vaccine/TLR binding was mediated by the vaccine’s anti-parallel β-sheets near its N-terminus. Notably, the vaccine N-terminus was further extended towards the other side, yet it was kept at relatively small distances. Representing the handle of a paddle, the α-helices at the vaccine’s carboxy end depicted extended conformation towards the solvent side being quite far from the TLR binding interface. Analyzing the complex interfaces via PDBePISA illustrated a relatively larger interface solvent-accessible area as well as interaction contacts (combined hydrogen bonds and salt bridges) for the obtained vaccine towards TLR4 over those at TLR2 (Table 5). Anchoring preferentiality at TLR4 was also translated into higher negative Δ1G scores as well as lower Δ1G p-value implying profound interaction-specificity for the vaccine and TLR4 interface surface being of higher hydrophobicity than would be averaged for given structures. Based on the above findings, the obtained vaccine-TLR complexes were considered significant for further interface evaluation and to be used for subsequent molecular dynamic simulation studies.

Key binding residues at the interface between the vaccine and each TLR target were also highlighted via PDBePISA interface analysis (Table 6; vaccine residue will be represented in italics starting from here forward). At the TLR4 bound complex, the vaccine Arg17 was depicted as significant for mediating highly network polar contacts, both hydrogen bonds, and salt bridges, with the corresponding TLR4 interface anionic residues (aspartates; Asp60 and Asp84) (Figure 4A). Additionally, the cationic Lys13 and Lys65 at the vaccine’s loop region mediated both hydrogen bond and salt bridge with the respective carboxylate sidechain of Glu154 and Asp60 at the TLR4 side. Several other polar vaccine residues including; Tyr231, Gln187, Asn421, as well as the mainchain of Phe244 predicted relevant hydrophilic contacts with TLR4 amino acids at distinct proximal distances (1.19 to 3.48 Å). Besides polar interacting residues, hydrophobic interface residues of the multiepitope vaccine showed significant closeness and relevant non-polar contacts with neighboring TLR4 amino acids. Vaccine hydrophobic residues including; Trp181, Phe183, Tyr202, Leu208, Phe227, Tyr231, Tyr242, Phe244, Tyr286, Leu288, Tyr292, Tyr303, Phe305 depicted ≤5.0 Å distances from Met41, His229, His256, Phe377, Tyr403, Tyr451, His426, Ile450, Phe408, Phe500, and/or Val602 of the TLR4 interface, respectively.

Moving towards the vaccine-TLR2 docking complex, several polar contacts were depicted majorly through hydrogen bond contacts without relevant salt bridges (Figure 4B). The vaccine’s acidic residue, Glu182, was the most frequent amino acid contributing to the vaccine’s parallel β-sheet stability towards the TLR2 concave interface. Other polar vaccine residues mediated significant hydrogen bonding with corresponding TLR2 amino acids at close distances and relevant angles (1.65 to 2.35 Å and 141 to 179 Å). Interestingly, the peptide backbones of several vaccine’s hydrophobic residues including; Trp62, Pro185, Leu189, Gly195, Ile198, Pro200, Leu288, and Tyr324, depicted relevant polar contacts with close TLR2 residues at the binding interface. In addition to hydrophilic residue pairing, hydrophobic interface residues of the multiepitope vaccine showed significant closeness and relevant non-polar contacts with neighboring TLR2 amino acids. Vaccine’s hydrophobic residues such as; Phe58, Val59, Trp181, Phe183, Met186, Phe191, Ile199, Leu208, Trp210, Leu225, Phe227, Leu229, Ala238, Val240, Tyr242, Phe246, Tyr286, Phe290, Tyr292, Leu301, Tyr302, Phe305, Val325, Ala327, and Leu329 showed distances less than 5.0 Å from Tyr109, His202, Ile204, Leu250, Ile251, Leu280, Tyr364, Leu392, Ile393, Tyr440, Ile461, Tyr483, Met500, Leu502, Val503, Trp529, His 531, and/or Pro568 of the TLR4 side.

Putting all the above preferential per-residue interaction within energy terms, the MM/GBSA binding energy calculations were estimated for both docked vaccine/TLR complexes. As expected, higher residue-wise energy contributions were assigned for the key interacting residues of both vaccine and TLR proteins (Table 6). The higher negative total binding energy has been correlated to the preferential vaccine affinity towards the TLR4 interface (−198.21 Kcal/mol) as compared to that of TLR2 (−127.72 Kcal/mol). Electrostatic energy contributions (−8,323.49 Kcal/mol and −5,374.32 Kcal/mol) were higher than those of hydrophobic van der Waal potentials (−270.25 Kcal/mol and −251.28 Kcal/mol) for TLR4 and TLR2 complexes, respectively. A penalty for polar solvation-free energy was estimated at 8,431.56 Kcal/mol, while a non-polar contribution at −35.93 Kcal/mol corresponds to a relatively large TLR4 hydrophobic surface available for vaccine binding. A similar pattern of polar and non-polar solvation-free energies was depicted at the TLR2 interface, yet at lower values (5,528.13 Kcal/mol and −30.25 Kcal/mol) than those at TLR4. The electrostatic potential preferentiality was further illustrated since several key binding polar residues showed the top-five binding energy contribution values (Arg17–11.06 Kcal/mol, Tyr231–9.01, Asp60–7.18 Kcal/mol, Asp84–7.01 Kcal/mol, and Glu154–5.71 Kcal/mol) conferring their significant role for vaccine binding and TRL4-associated complex stability. Similarly, the vaccine’s polar residues; Glu182–8.49 Kcal/mol, Gln187–7.01 Kcal/mol, and Ser326–6.82 Kcal/mol, as well as those of TLR2 interface; Cys226–5.51 Kcal/mol, Asp366–4.28 Kcal/mol, and Glu310–3.79 Kcal/mol, were assigned with the high binding energy contributions.

The RMSD trajectories of the multitope vaccine, TLR proteins, and their respective complex were monitored across the 100 ns all-atom simulation runs in reference to the alpha-carbon atoms (αC) of their respective initial structures. All monitored RMSD tones showed a gradual increase for the first 25 ns, while after that RMSDs started to attain steadier trajectories till the end of the simulation runs (Figure 5A). The latter corresponds to optimum thermodynamic behavior since all the constraints were removed with the simulation start and the proteins began to relax, till they converged at their equilibration plateau having their RMSDs being fluctuated across their respective averages. General trends of higher RMSD tones and fluctuations were assigned for the vaccine/TLR2 complex in relation to that with TLR4 (6.92 ± 1.16 Å and 11.34 ± 1.83 Å for TLR4 and TLR2, respectively). Regarding the comparative TLR protomeric units, significantly lower RMSD values were depicted for TLR4 (3.65 ± 0.62 Å) than for TLR2 (7.39 ± 1.33 Å). Similarly, the vaccine in complex to TLR4 depicted steadier and less fluctuation RMSD tones (9.77 ± 1.73 Å) as compared to its RMSD tones in complex with TLR2 (15.23 ± 2.56 Å). It is worth noting that the monitored RMSDs of the vaccine/TLR complex were more influenced by the TLR thermodynamic behaviors over those of the vaccine showing values close to the receptor alpha-carbon atoms.

Dissecting the protein fluctuation pattern down to their respective residues, the αC-RMSF trajectories were monitored across the simulation run for each simulated protein, including the vaccine and TLR protomeric unit (Figure 5B). Notably, vaccine residues depicted higher fluctuation patterns as compared to their corresponding bounded TLR, the thing that was further recapitulated by the above-depicted RMSD tones. Additionally, general trends of higher RMSF values were assigned to the vaccine-TLR2 complex over those for vaccine-TLR4 ones (3.89 ± 2.80 Å/1.72 ± 0.55 Å for vaccine/TLR4 versus 6.16 ± 3.91 Å/2.73 ± 0.76 Å for vaccine/TLR4). The highest vaccine fluctuations up to 10.05–12.72 Å were assigned to its C-terminal residues and vicinal amino acid region (extending across Ala327-to-Pro349 at both TLR4 and TLR2 systems; Glu369 to Ala303 only at TLR2 system). On the other hand, the vaccine’s N-terminal residues were of relevant stability and lower RMSF values. However, single residue-range Thr100-Gly162 depicted the highest fluctuations at the TLR4-bounded system (RMSF up to 13.03 Å) and even much higher flexibility at the TLR2-bounded one (RMSF up to 19.87 Å). It is worth noting that both vaccine’s side high-fluctuating residue regions were parts of the vaccine that are greatly extended towards the solvent side, quite distal from the TLR binding interface. Regarding the inherited fluctuations of the simulated TLR proteins, RMSF tones were at regular fluctuation patterns resembling repetitive camel humps. Interestingly, typical protein dynamic behaviors were depicted for simulated TLRs since their core residue ranges were assigned with lower RMSFs than those at the terminal sides. The observed differential fluctuation patterns between the vaccine and/or TLRs would question the influence of the protein’s inherited ternary structure, dimerization interface, and vaccine binding on protein simulation patterns.

Two other trajectory-based stability parameters, ROG and SASA, were also monitored across the simulation runs. In relation to their respective central masses, both simulated multitope vaccines showed differential ROG values across different simulation times (Figure 5C). Bound to the TLR4 protein target, the vaccine’s ROG tones gradually decreased across an average of 32.50 ± 0.64 Å from the simulation start till halfway of the simulation run (∼55 ns). Afterward, the vaccine’s ROGs were maintained across a plateau of 30.13 ± 0.29 Å which was kept steady with minimal fluctuations till the end of the simulation. Similar dynamic behavior was depicted for the TLR2-bound vaccine, yet higher ROG values (33.15 ± 0.74 Å) were obtained while being accompanied by higher fluctuations, particularly at the first half of the simulation run. This was reflected in the total ROG of the vaccine-TLR2 complex showing significant fluctuations before the 50 ns frames. Notably, both described vaccine’s ROG findings conferred differential conformation changes being attained by the vaccines beyond the first 50 ns time frames. Regarding TLR data, much steadier ROG tones were depicted for the target protein as compared to those of the bound vaccines owing to the TLR inherited stability and higher central masses. Limited ROG fluctuations were depicted for TLR4 in relation to those TLR2 (32.65 ± 0.25 Å versus 31.37 ± 0.33 Å, respectively). Complex ROGs were maintained steady over an average of 43.64 ± 0.36 Å for TLR4 and 35.45 ± 0.41 Å for TLR2, with relevant fluctuations below 50 ns time frames. It is worth noting that ROG data should be carefully interpreted since higher ROG of TLR4 and related complex are more reasoned for the larger size (more constituting residues) of TLR4 over TLR2, rather than just reflecting structure compactness and stability (Tanner, 2016).

Regarding the SASA trajectory-based analysis, buried surface area (BSA; Ų) for each simulated complex was estimated by adopting the sole SASAs of each bounded protein (SASA_vaccine and SASA_TLR) in addition to the entire bounded complexes; BSA = 0.5*(SASA_vaccine + SASA_TLR – SASA_complex) (Zhang et al., 2022). The simulated vaccine-TLR4 complex showed higher BSA values around 2,666.09 ± 162.77 Å² as compared to the investigated vaccine-TLR2 model (2,276.82 ± 237.73 Å²) (Figure 5D). Notably, the vaccine-TLR4 system illustrated increased BSA beyond the 50 ns time window while greater fluctuations were demonstrated at the initial frames. The latter dynamic behavior indicated conformational stability, as well as more surface protein areas, which were covered at the second half of the simulation run and till its end. On the contrary, the vaccine-TLR2 complex depicted lower BSAs across the second half of the simulation run conferring fewer areas being covered. Nevertheless, the vaccine-TLR2 BSAs depicted steadier tones following the 50 ns time frame as compared to its initial times. Further tracking of the vaccine’s thermodynamic behavior was done by exploring the time-evolution conformational changes of the vaccine/TLR complex following the molecular dynamics run.

Conformational analysis through aligning the starting and final complex structures illustrated differential orientations for the bound ligand at the TLR binding interface (Figure 5E). Observed visually and correlated to the above obtained RMSD and RMSF trajectories, limited conformational and orientation changes were observed for the TLR proteins. A slight conformational shift was depicted for the TLR protomers showing minimal movement in its flexible loops (aligned RMSD 1.56 Å and 2.22 Å for TLR4 and TLR2, respectively). On the contrary, more dramatic conformational and orientational shifts were assigned to the vaccine protein with higher RMSD between its initial and final frame (5.33 Å and 11.57 in bound to TLR4 and TLR2, respectively). The latter RMSD values showed greater conformational/orientational shifts for the vaccine at the TLR2 interface. Significant drift (∼40 Å) was depicted for the vaccine’s carboxy-terminal and vicinal residue range (Ala327-to-Pro349) as well as the solvent-exposed amino acids (Thr100-Gly162) near the N-terminus. In the TLR4-bound complex, the latter vaccine’s high-flexible regions adopted more compacted conformations as they became more directed toward the lateral interface of the TLR4 protein at the end of the simulation time. On the contrary, the vaccine’s flexible regions at the TLR2 complex could not manage to adopt similar compactness the thing that was translated into higher distortions for the vaccine’s C-terminal side. Such dynamic behavior could explain the higher RMSF-related fluctuation patterns associated with these solvent-exposed residue regions as well as the higher RMSD/Rg tones for TLR2-associated vaccine over that at the TLR4 interface. The rest of the vaccine ternary structures at both TLR complexes showed limited conformational changes, particularly those being anchored deep at the TLR binding interface (inner concave surface).

Moving towards the vaccine-TLR binding affinity across the simulated thermodynamic trajectories, the MM_PBSA calculations revealed great concordance with the preliminary HawkDock-MM/GBSA docking binding scores. The trajectory-based MM_PBSA revealed higher negative average free binding interaction energies for the vaccine towards the TLR4 receptor interface as compared to those at the TLR2 side (ΔG_Total = −10804.19 ± 138.19 kJ/mol versus −2,454.162 ± 410.95 kJ/mol) (Figure 6A). At both vaccine-TLR complexes, pronounced dominance of the electrostatic potentials (ΔG_{Electrostatic}) was depicted over the van der Waal interactions (ΔG_vdw). The electrostatic potential energy contributions were more than 3-folds and 12-fold that of the van der Waal interactions at TLR4 and TLR2 interfaces, respectively. Notably, the polar solvation penalties (ΔG_{Polar solvation}) at both TLR2 and TLR4 systems were almost comparable and were translated as repulsive forces against vaccine binding as well as compromised complex stabilities since binding is a solvent-displacement process (Barillari et al., 2007). Nevertheless, the extremely higher attractive electrostatic potentials as well as the total non-polar interactions (summation of ΔG_SASA and ΔG_vdw) for the TLR4-associated complex over the TLR2 would have managed to overcompensate these unfavored repulsive polar solvation forces. All these furnished binding energy terms were translated into better vaccine affinity towards TLR4 over the TLR2 side. Finally, the relatively high total non-polar potentials could be correlated to the large hydrophobic receptor interface available for the multitope vaccine binding.

Exploring the residue-wise free binding energy contributions for the simulated complex were highlighted in Figures 6B, C. Higher negative-value energy contributions were assigned for the vaccine N-terminal residues as well as those constituting the extensive anti-parallel sheets. Repulsive positive-value energy contributions were limited to the residue area being exposed towards the solvent side, particularly C-terminal residues and vicinal amino acid region. Findings were consistent with the above-depicted RMSF values showing residues with high positive-value energy contributions to be correlated with high respective RMSF trajectories across the simulation time. Concerning the protein receptor, higher negative-value residue-wise energy contributions and lower repulsive positive contributions were assigned for TLR4 over the TLR2 unit. On the other hand, preferentially higher negative-value energy contributions were assigned for the N-terminal half of each TLR structure.

The computational assessment through the C-ImmSim server predicted an overall satisfactory response with the successive doses of the potential multitope vaccine injection (Figure 7). It is clear that a high level of IgM + IgG antibodies was induced upon the multitope vaccine injection. Moreover, as we considered the induction of cytokines during T cell epitope filtration, it is not a surprise that the chimeric construct was predicted to stimulate several cytokines with INF-γ coming at the top of the stimulated cytokines list. Finally, the server output revealed that there was an obvious increase in both B and Th cell levels upon the potential vaccine injection where the highest level was observed after the second booster dose.

The output of the JCAT servers showed that the GC content of the improved sequence was 51.6% (the accepted range is between 30% and 70%. In addition to that, the calculated CAI value was 0.95, a value that indicates the suitability of the improved sequence to be easily expressed in E. coli k-12 (the CAI value ranges between 0 and 1 and the accepted range is between 0.8 and 1).