The complete genome sequences of the 14 strains of P. melaninogenica were used to retrieve a core genome from the EDGAR tool version 2 (Blom et al., 2016). Clustering techniques were employed to identify duplicate or paralogous proteins via CD-HIT analysis (Li and Godzik, 2006), to remove redundancy. The cutoff value of 60% sequence similarity was kept as a threshold to obtain non-paralogous proteins, which were used for further analysis (Qasim et al., 2023). Alignment coverage was obtained by setting parameters at default.

Human host non-homologous proteins were determined by standalone BLASTp analysis, which was used for scanning against species non-paralogous proteins with E-value 1e-20, bitscore ≥100, percent identity ≥35, and query coverage ≥35 (Shah et al., 2021). The human proteome was downloaded from the UniProt database, comprising 1,076,164 proteins.

Further scanning of pathogen proteins was carried out against the human gut microbiome. The NCBI database was used to acquire a human gut microbiome consisting of 75,176 proteins. BLASTp was performed with criteria as follows: E-value 1e-4, bitscore ≥100, percent identity ≥50, and query coverage ≥35, as in our previous studies (Aslam et al., 2020; Shah et al., 2021; Jaan et al., 2022; Qasim et al., 2023). All the homologs were discarded.

After excluding gut microbiome homologs, essential proteins of the pathogen were determined by subjecting proteins to the database of essential genes (DEG) (http://tubic.tju.edu.cn/deg/) (Luo et al., 2021). DEG contains all the genes required for the survival of the pathogen. BLASTp was performed for scanning against DEG with cutoff parameters of E-value 1e-4 and bitscore ≥100 for the determination of P. melaninogenica essential genes.

The host defense mechanism is destroyed by the virulence factors of bacteria, which help to cause the disease through invasion, colonization, and adhesion. Four categories of virulence factors are included in the virulence factor database (VFDB) (Liu et al., 2019), which is used to identify virulence factors in a pathogen. From 25 pathogenic bacteria, offensive, defensive, non-specific, and virulence-related proteins are assembled in this database. Pathogen proteins were subjected to BLASTp analysis against VFDB with parameters, i.e., E-value 1e-4 and bitscore ≥100 (Qasim et al., 2023).

Databases were employed to identify resistant and host–pathogen-interacting proteins with criteria such as E-value 1e-4 and bitscore ≥100. There is an emerging trend toward developing resistance to available antibiotics among bacterial pathogens. Moreover, antibiotic-resistant genes were identified by the ARG-ANNOT database (Gupta et al., 2014), with selected criteria by subjecting pathogen proteins to BLASTp (Aslam et al., 2020). The host–pathogen interaction database (HPIDB v2.0) (Ammari et al., 2016) was used for the curation of host–pathogen protein–protein interactions, which were essential for the survival of a pathogen. Standalone BLASTp was used for the screening of interacting proteins against the HPIDB repository (Aslam et al., 2020).

Metabolic pathways of living organisms can be retrieved from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Fatoba et al., 2021). Human non-homologous essential proteins were further subjected to KEGG for metabolic pathway determination. Functional annotation of proteins was provided by the KAAS server via these metabolic pathways (Nazir et al., 2018). Metabolic pathways of the human host and pathogen were manually compared to find unique and common pathways in the pathogen. KEGG automatic annotation server (KAAS) (Moriya et al., 2007) employs BLASTp analysis for these metabolic pathways. KEGG orthology (KO) identifiers were provided by the KAAS server for proteins (Fatoba et al., 2021). Unique and common pathway proteins were enlisted and categorized into KEGG-dependent and independent proteins. Unique proteins are those involved only in the pathogen-specific pathways but not in humans.

Screening of essential, non-homologous proteins was further performed by subjecting proteins to BLASTp analysis with a cutoff E-value of 1e-4, against the DrugBank database for identification of druggable proteins (Aslam et al., 2021). All FDA-approved drugs and experimental drugs are contained in the DrugBank database (Law et al., 2014). Druggable targets showed significant matching (80% or more) with FDA-approved drug targets, while proteins not showing any significant matching with already available targets were considered novel targets and used for subsequent analysis.

Proteins are located in different positions, including the cytoplasm, inner membrane, periplasmic membrane, and outer membrane, by subcellular localization prediction. Proteins are categorized as drug or vaccine targets according to their location. Cytoplasmic and outer membrane proteins are considered suitable drug and vaccine targets, respectively (Jaan et al., 2022). PSORTb v3.0 (http://www.psort.org/psortb/) was employed for the prediction of the subcellular location of the resulting target proteins. Further validation of proteins was needed to ensure the location, which was carried out by a server named subcellular localization predictor (CELLO v2.5) (http://cello.life.nctu.edu.tw/) (Yu et al., 2014). The STRING database v10.5 (http://string-db.org) (Szklarczyk et al., 2016) was employed for the analysis of protein–protein interactions by setting parameters at default. The number of interactions was presented by hub proteins showing node degree (K ≥ 5). The transmembrane topology of proteins was predicted by the TMHMM v0.2 server (Krogh et al., 2001). Proteins spanning the lipid membrane are termed transmembrane proteins and are prioritized as potent drug targets (Shah et al., 2021).

Proteins prioritized in the previous step were employed in SWISS-MODEL for the generation of 3D structures. The server uses homology modeling techniques for the prediction of protein structures when experimental structures are not accessible. Reliable protein templates with expected accuracy were generated (Jaan et al., 2022). A protein verification tool called “ERRAT” was used to verify the quality of the 3D structure of proteins. ERRAT signifies the quality of protein by a parameter called an overall quality factor, where a high value (≥50) specifies good quality for non-bonded atomic interactions (Aslam et al., 2020). The molecular weight and theoretical PI of proteins were determined by the ExPasy server.

The estimation of protein pocket ability to bind small molecules with high affinity is a valuable step in the identification of drug targets. PockDrug, an online server (http://pockdrug.rpbs.univ-paris-diderot.fr.), was used to calculate druggable scores for the prioritized targets. Protein structure information and ligand proximity assessment were used to predict the pocket druggability of proteins by using different estimation methods via the PockDrug server. By using pocket estimation methods, consistent druggability results were obtained, which differentiated more druggable proteins from less druggable proteins (Shah et al., 2021).

For the evaluation of the antigenic nature of proteins, VaxiJen server v2.0 was employed for outer membrane and extracellular proteins with criteria, i.e., accuracy rate 70–89% and probability score >0.5 (Doytchinova and VaxiJen, 2007). The allergenicity of the shortlisted proteins was evaluated by the AllergenFP server to distinguish allergens from non-allergens. Allergic proteins that are harmful to the host were discarded in this evaluation (Dimitrov et al., 2014). Different parameters of proteins were assessed by the ExPasy ProtParam tool, calculating molecular weight, theoretical PI, instability index, aliphatic index, half-life, number of amino acids, and GRAVY (Grand Average of Hydropathicity) (Qasim et al., 2023).

For a specific pathogen, vaccine targets are identified with the help of T-cell epitopes involved in the elicitation of immune responses. An antigen is processed by an MHC molecule to bind with either cytotoxic or helper T cells (Shah et al., 2021). CD4 or CD8 T cells are stimulated by an epitope defined as the shortest immunogenic peptide having the ability to stimulate immune responses. Two types of major histocompatibility complex (MHC) classes I and II recognized by CD4 and CD8, respectively, are represented by T-cell epitopes (Aslam et al., 2020). MHC prediction was carried out by the Immune Epitope Database (IEDB) server. For the prediction of MHC class-I epitopes, vaccine candidates were subjected to the NetMHCpan 4.1 web server with a percentile rank <0.2, recommended by the IEDB server (Fleri et al., 2017). MHC class-II epitopes were predicted by the IEDB recommendation of 2.22 with an adjusted rank of <0.5. Further evaluation of preferred epitopes was performed.

The solvent-exposed regions in the antigen can be identified by the prediction of B-cell epitopes, by which antibodies can be recognized (Nazir et al., 2018). For the prediction of B-cell epitopes, different servers were used. BCPreds, FBCPreds, AAP, and BepiPred online servers were used to predict linear B-cell epitopes. BCPred uses the SVM method, which has 5-fold cross-validation and utilizes five different kernel methods. Linear-length B-cell epitopes were predicted by subsequent kernel followed by FBCPred. Biochemical properties such as amino acid composition, hydrophobicity, hydrophilicity, secondary structure, and surface accessibility of peptides were used to predict linear epitopes by BepiPred (IEDB) (Solanki and Tiwari, 2018).

A designated vaccine candidate must be able to act as an antigen. For this purpose, the VaxiJen v2.0 server was employed with a threshold value of >0.4. The sequence alignment method is used in the VaxiJen server, and the physicochemical properties of peptides are used to confirm their antigenicity (Doytchinova and VaxiJen, 2007). The allergenic properties of vaccine candidates were measured by the AllergenFP server. Another server, ToxinPred, was employed to check the toxic and non-toxic peptides, where non-toxic peptides were selected for further evaluation (Jaan et al., 2022). The conservancy level of epitopes was evaluated by the IEDB server by setting parameters at default, where conserved or variable epitopes could be assessed (Jaan et al., 2022). These properties of epitopes were evaluated for the construction of chimeric vaccines.

Multi-epitope vaccine constructs were generated by joining overlapped epitopes of MHC-I, MHC-II, and B cells with the help of amino acid linkers (i.e., EAAAK, GGGS, and KK). To increase the immunogenicity of vaccine constructs, adjuvants were added to the sequence. Four adjuvants, i.e., HBHA protein, L7/L12 ribosomal proteins, beta-defensin, and HBHA conserved sequences, were added to improve the efficacy of vaccine models. PADRE peptide sequences, containing 13 amino acids (i.e., AKFVAAWTLKAAA), were also added, which helped induce CD4+ T cells (Aslam et al., 2020). PADRE sequences help to overcome the problems of polymorphism and elicit better immune responses. Different combinations of vaccine constructs were made and checked for immunogenicity, toxicity, and allergenicity. The SOLPro server was employed to check the solubility of vaccine constructs (Magnan et al., 2009).

Intrinsic physical and chemical characteristics of a substance are defined as physicochemical properties. For the estimation of physicochemical properties, the ProtParam tool was used, which calculated the number of amino acids, half-life, theoretical PI, instability index, aliphatic index, molecular weight, and GRAVY of vaccine models (Jaan et al., 2022).

For the prediction of secondary structure, two online servers, i.e., PSIPRED and SOPMA, were used. Properties such as transmembrane helices, bend regions, random coils, and beta sheets can be assessed through these self-optimized prediction methods. A graphic presentation of proteins is obtained as output (Jaan et al., 2022). 3D modeling of chimeric constructs was generated by Phyre2, increasing the accuracy of alignment by using the alignment of the Hidden Markov model (Raza et al., 2021). GalaxyRefine and 3DRefine servers were employed for the refinement of predicted 3D models, in which side chains are rebuilt and repacked. The probable errors of 3D models were checked by the ProSA-web server (Jaan et al., 2022). The Ramachandran plot was obtained by an online server called PROCHECK, which explains the stereochemical properties of the protein (Shah et al., 2021). Phi/Psi angles were assessed by the Ramachandran plot for a comprehensive understanding of the protein backbone.

Vaccine constructs were disulfide-engineered by DbD2, to check the stability and conformational entropy of the protein. Increased stability and decreased conformational entropy of protein can be confirmed by the accessibility of protein, and it can add novel disulfide bonds (Craig and Dombkowski, 2013).

An energy-minimizing step used to obtain a stable structural configuration of vaccine targets with the ligands is determined by molecular docking analysis. Immune receptor–peptide interactions were inferred by docking vaccine constructs with six various human HLA alleles using the PatchDock server. More refinement and re-scoring of docked complexes were performed by the Fast Interaction Refinement in Molecular Docking (FireDock) server, which gave the 10 best models as a result. Based on the lowest binding energy, the vaccine construct (V5) was prioritized as the best-docked complex (Jaan et al., 2022). Furthermore, TLR4 (toll-like receptor) acquired from PDB (2Z65) was docked with vaccine construct V5 in the same way by the PatchDock server. The top 10 generated models from the FireDock server were redirected to PatchDock for refinement. Assumed refined solutions were recognized by the binding score and global binding energies of docked refined models.

The stability of the docked complex and molecules' behavior was assessed by the molecular dynamic approach. The interaction between the designed vaccine and receptor was estimated by the iMOD server. Four main factors are measured by this tool by estimating the direction and range of basic movements of the docked complex. These factors include Eigenvalue, B-factors, covariance, and deformability. A high Eigenvalue indicates the value of much harder deformation (Jaan et al., 2022).

The possible expression level of vaccine constructs was tested by using the Java Codon Adaptation tool (JCAT), which used back-translated amino acid sequences. E. coli (K12 strain) was selected for the expression of the protein. The high expression level in E. coli was verified by the consequential GC content and codon adaptation index (CAI) (Jaan et al., 2022). The rho-independent transcription terminators, cleavage sites of some restriction enzymes, and prokaryotic ribosome binding sites were dodged. SnapGene software was used to introduce codon-adapted sequences into plasmid vector pET30a (+), to construct recombinant plasmid sequences (Solanki and Tiwari, 2018). C-ImmSim server was employed for immune simulation to predict the production of interferons, antibodies, and cytokines in response to a foreign particle. Parameters were set to default (Puzone et al., 2002).

The core genome method uses conserved sequences of available complete genome strains in accordance with the reference genome. A total of 1,415 core proteins were retrieved from 14 complete genome strains of P. melaninogenica. Some proteins are produced in response to duplication events taking place in organisms within one species, known as paralogous proteins. To reduce redundancy in sequences, proteins were subjected to CD-HIT analysis (Li and Godzik, 2006). A total of 1,414 non-paralogous proteins were obtained by subtracting one duplicate protein from the core protein sequences. The threshold was set to 0.6 (60%) sequence similarity.

For the identification of novel targets, bacterial proteins were subjected to homology filters as an essential step for the screening of human host homologs. After removing the proteins homologous to the host proteins, 1,295 candidates were obtained by screening against the human proteome database. Different metabolic reactions such as homeostasis, development, defense system, and physiological functions are performed by the gut flora, which colonizes the human body throughout life (Kho and Lal, 2018). The scanning of core proteins against the human gut microbiota was performed by standalone BLASTp, to prevent the likelihood of serious complications in the host. The resulting 615 proteins were used for downstream analysis.

A minimal set of genes required to perform basic cellular functions and important for the survival of an organism are defined as essential genes. A database of essential genes (DEG) based on BLASTp scanning was employed to dig out essential proteins of the pathogen. As a result, 123 proteins were recognized as essential for the maintenance of cellular functions in P. melaninogenica. The virulence factor database (VFDB) was used to identify the virulence proteins of the pathogen. The pathogenic capacity of bacteria to infect the host is increased by virulence factors (Aslam et al., 2020). By subjecting proteins to BLASTp scanning, VFDB calculated 19 proteins as virulence factors.

Antibiotic resistance proteins were determined by Antibiotic-Resistant Gene ANNOTation (ARG-ANNOT) (Gupta et al., 2014). This analysis determined three proteins as resistant to antibiotics, which might act as important drug targets in future studies. Host–pathogen protein–protein interactions (HP-PPI) mediate infectious diseases in response to molecular cross-conferences between host and pathogens. For treating infectious diseases, such proteins need to be identified for the discovery of potential drug targets. Proteins were subjected to BLASTp against HPIDB v2.0, which identified only one interacting protein (Ammari et al., 2016). The proteins from this analysis have the potential to be selected as alternative drug targets and deserve further experimental investigation.

Metabolic pathway analysis is used to find out the preferable drug target as involved in pathogen-specific unique pathways (Dorella et al., 2013). The KEGG database was used to predict the metabolic pathways of Homo sapiens (humans) and P. melaninogenica, which were 345 and 94, respectively. A manual comparison of human and bacterial metabolic pathways showed that 25 unique pathways were present but absent in the human host, while the remaining 69 pathways were common to both humans and bacteria. The KAAS server was employed for the functional pathway analysis of 114 essential, non-homologous pathogen proteins via BLASTp, which assigned them a KO (KEGG orthology) identifier (Fatoba et al., 2021). Out of 114 proteins, only 5 were involved in pathogen-specific unique pathways called KEGG-dependent (Table 1), while the rest of the 109 proteins were categorized as KEGG-independent.

The potential therapeutic targets were represented by the metabolic pathway analysis of non-homologous, essential human proteins. Suitable and least-targeted were identified by applying filters to reduce time, labor, and resources. Druggability screening of the prioritized proteins was performed to find novel druggable targets, which demonstrated that in KEGG-dependent, one protein is a novel target, while in KEGG-independent proteins, 90 proteins are novel and have the potential to be explored as potential drug targets.

Suitable and effective targets can be identified by predicting subcellular location for a comprehensive understanding of the function and mechanisms of proteins. Purification and assay of membrane-localized proteins are critical, so cytoplasmic proteins are favored as potent drug targets (Mondal et al., 2015). For suitable vaccine targets, outer membrane and extracellular proteins were prioritized for elicitation of better immune responses as membrane proteins were exposed to host cells (Nogueira et al., 2021). Subcellular location prediction by CELLO and PSORTb tools showed three proteins to be cytoplasmic and only one protein to be outer membrane in the case of KEGG-dependent proteins. In the case of KEGG-independent proteins, 49 proteins were found to be cytoplasmic and 23 proteins to be outer membrane in location. These proteins were prioritized as drug or vaccine targets for further analysis in accordance with their location.

Further characterization of selected targets was performed by STRING v10.5 databases for the estimation of protein–protein interactions. A significant number of interactions were represented by node degree (K) ≥5 as hub proteins (Szklarczyk et al., 2016). This characterization prioritized 20 proteins as suitable targets (Figure 2). The TMHMM server was employed for the prediction of transmembrane helices set to be less than or equal to 1. Due to multiple transmembrane helices, there may be complications in the cloning, expression, and purification of proteins (Solanki and Tiwari, 2018). For this purpose, proteins with more than one helix were excluded, and 18 proteins were selected (Table 2). The molecular weight of selected targets was calculated by the ExPasy tool, which prioritized proteins with < 110 KDa molecular weight. Proteins submitted to SWISS-MODEL were checked for model quality estimation by different parameters, i.e., percentage identity, Q-mean, and GMQE scores, which prioritized three drug targets. Protein pocket ability prediction is one of the crucial steps in the development of therapeutic drugs to ensure the affinity of proteins with small drug-like molecules. The PockDrug server was employed for this purpose, which identified significant druggable targets with a drug score of >0.5 (Borrel et al., 2015). Furthermore, the ERRAT server was used to evaluate the quality of models, which was considered significant (85%) for all models.

Surface-exposed or outer membrane proteins which were exposed to host cells were prioritized for the reverse vaccinology, in the case of gram-negative bacteria. Certain characteristics (i.e., instability index, aliphatic index, PI, GRAVY, and molecular weight) of 24 outer membrane proteins were analyzed by the ProtParam tool. The AlgPred server was employed to check the allergenicity of proteins by using the hybrid method for the best vaccine candidate selection. Antigenicity was estimated by the VaxiJen server based on the threshold of >0.5. Two proteins with IDs ADK95685.1 and ADK97014.1 were selected as the best vaccine targets used for subsequent analysis. By this analysis, it was ensured that induced immune responses were aimed to target the pathogen, not the host cells.

Targeted immune responses are induced by the short peptide fragments of a peptide vaccine; accordingly, allergenic sequences are avoided. A variety of peptides termed epitopes can bind to MHC molecules with high affinity as mutations are possible in MHC-binding epitopes caused by pathogens (Shah et al., 2021). One of the major keystones in the development of vaccines is the B- and T-cell epitope prediction (Li et al., 2014). MHC-I and MHC-II molecules provide both types of T cells, i.e., (i) cytotoxic T cells and (ii) helper T cells, with epitopes. High-binding affinity epitopes were predicted by multiple prediction sources to find potent vaccine-candidate epitopes. For the prediction of MHC-binding epitopes in two prioritized vaccine candidates, an IEDB server was employed. For MHC-I binding, the NetMHCpan 4.1 server was used with criteria of <0.2 percentile rank. Hydrogen binding can attach peptides with only 8 to 10 amino acids due to the closed binding site cleft of MHC-I molecules (Nazir et al., 2018). Criteria for binding with human alleles were set to 9-mer length, and predicted 19 and 16 epitopes for both vaccine candidates were used for further analysis. These epitopes were further filtered by toxicity, antigenicity, conservancy, and immunogenicity analyses, which resulted in 5 and 6 final epitopes for chimera vaccine construction.

Due to open and shallow pockets, prediction by MHC-II molecules is less accurate. The IEDB-recommended 2.22 method was employed for the prediction of MHC-II epitopes with a percentile rank <0.5, and the 15-mer length of amino acids was selected due to the higher binding affinity of these molecules as compared with MHC-I (Shah et al., 2021). After screening, the final MHC-II and MHC-I epitopes were selected (Supplementary Tables S1, S2). For the prediction of linear B-cell epitopes, different online servers, i.e., BCPred, BepiPred, FBCPred, and AAP, were employed. The predicted epitopes of T and B cells were compared for the determination of overlapping epitopes (Table 3). These overlapping epitopes may be useful in the development of a vaccine against P. melaninogenica in future. Conformational B-cell epitopes have an important role in antigen and antibody association, but unfortunately, these cannot be predicted due to the inaccessibility of the 3D structures of all the proteins.

For the construction of chimera vaccine constructs, different combinations of non-toxic, antigenic, and conserved lead epitopes were chosen. On the N-terminal of prioritized epitopes, adjuvants were added to enhance the immunogenic nature of the multi-epitope vaccine. Linkers, i.e., EAAAK, GGGS, and KK, were used to join the adjuvants. Allergenicity and antigenicity of the constructs affected by adjuvants were analyzed by the different combinations of adjuvants (Solanki and Tiwari, 2018). The conformation of designed constructs was not altered by adding linkers. Around the global population, the problems triggered by polymorphisms in HLA-DR molecules can be overcome by adding PADRE sequences (synthetic peptides containing 13 amino acids) to vaccine constructs. Better immune protection and cytotoxic T-lymphocyte (CTL) responses were provided by PADRE sequences (Wu et al., 2010). Different combinations of lead epitopes, adjuvants, and PADRE sequences provided eight vaccine constructs, i.e., V1, V2, V3, V4, V5, V6, V7, and V8 (Table 4).

The antigenic behavior of vaccine constructs was examined by the AntigenPro server with a cutoff value of >0.90. The VaxiJen server was employed to estimate the antigenicity, with a threshold of 0.75 showing good antigenic behavior. The allergic nature of constructs was checked by the AllergenFp, and allergic constructs were excluded as being harmful to the human host. The SOLPro server with a probability score of >0.5 was employed to analyze the solubility of the designed chimeric constructs. GC content (30–70%) and CAI value (0.90–1.0) were considered significant as inferred from the results. The ProtParam tool was used to evaluate properties such as the number of amino acids, aliphatic index, instability index, PI, molecular weight, and GRAVY. Eventually, vaccine constructs V5 and V8 fulfilled the standard criteria and were selected for further investigations (Supplementary Tables S3, S4).

A 2D diagram of the vaccine represents the residues involved in the formation of coils, b-sheets, and a-helices. These particular secondary structures of vaccine are involved in its overall structural stability. The 2D diagram indicates a big picture of the complete structure of the vaccine construct. Therefore, prioritized vaccine constructs (V5 and V8) were subjected to PSIPRED and SOPMA servers for the prediction of secondary structure to evaluate the stability of peptides (Raza et al., 2021). The results revealed by servers calculated the percentage of α-helix (36.90), β-strands (10.69), extended strands (21.80), and random coils (30.61) for V5 (Supplementary Figure S1A). Similarly, percentages for V8 were 29.43, 24.44, 11.22, and 34.91, respectively (Supplementary Figure S1B). The function and stability of proteins are critically analyzed by tertiary structure prediction. The Phyre2 server was employed for the prediction of the tertiary structure of prioritized vaccine constructs (V5 and V8). The vaccine construct V8 was deselected due to limited data for three-dimensional structure formation. Hence, V5 was used for tertiary structure, and distortions in protein structure were overcome by refinement of structure. Closeness to the native structure of proteins was acquired by the refinement of models by GalaxyRefine and 3DRefine. Efficient protein structure refinement is obtained by applying knowledge-based force fields and composite physics. Energy minimization steps and repetitive optimization of hydrogen bonding at the atomic level are employed by these online servers for refinement (Raza et al., 2021).

Validation of the refined 3D model was executed by PROCHECK and PROSAweb, which were used to generate the Ramachandran plot. This plot assessment elaborates on the combination and orientation of dihedral angles falling in the disallowed region on the basis of steric hindrance (Maxwell and Popelier, 2017). The Z-score (assessed by PROSA Web) improved from −2.98 to −5.37, which ensured the better quality of models. More than 90% of residues in the most favored regions of the Ramachandran plot indicate the better folding of the protein structure and validate its quality. The results of the Ramachandran plot analysis of our vaccine showed 95.2% of residues in the most favored regions, which indicated the good quality of our vaccine. Visualization of the tertiary structures was made with Pymol v2.5.5 software for better representation (Figure 3).

The designed vaccine construct was stabilized by performing disulfide engineering, where 13 pairs of amino acids were found to be able to make disulfide bonds predicted by the DbD2 server. However, considering the parameters such as energy and chi3, five pairs were considered favorable. The value of energy was calculated to be <4 kcal/mole for 5 pairs, with a chi3 value between −105 and +109.

Favorable interactions between vaccine constructs and HLA alleles were predicted by molecular docking analysis, an energy minimization step. The vaccine construct (V5) was docked with six different HLA alleles, i.e., 1A6A, 1AQD, 2Q6W, 3C5J, 4MD4, and 5NI9, acquired from PDB by using the PatchDock server. For refinement of models, the results were subjected to the FireDock server. Docked complex V5-1AQD fulfilled the criteria by showing the highest docking score (14814) as compared with the docking score of other docked complexes. Vigorous immune-stimulatory and CTL response effects can be elicited by interactions between TLRs (Rana and Akhter, 2016). The vaccine construct (V5) was further docked with TLR4 (toll-like receptor) to enhance immune responses (Table 5). Profound interactions of the docked complex V5-TLR4 were confirmed through binding energy.

In the cellular environment, the stability of docked complexes was inferred by molecular dynamic simulation studies. The iMOD server was employed to explore the movement of molecules and atoms. The relation of the docked complex between the NMA and the PDB sector was clearly visualized by the B-factor graph (Figure 4A). The coupling between residues was indicated by a covariance matrix, where the motions, i.e., correlated, uncorrelated, and anti-correlated, of different pairs were indicated by red, blue, and white colors, respectively (Figure 4B). The distortion of each residue was counted as the deformability of the complex, as shown in the graphs obtained by the server with higher peaks indicating the greater deformability (Figure 4C). A low Eigenvalue indicates that the complex can be deformed easily. The motion stiffness of the complex was represented by an Eigenvalue (9.549413e-05) (Figure 4D). The Eigenvalue and variance related to each normal mode were found to be in an inverse relation (Figure 4F). The elastic network model indicated the stiffness of the residues of the vaccine in the form of springs, with darker grays indicating stiffer springs (Figure 4E). The interpretation obtained from the results confirmed the stability of the docked complex (V5-TLR4) in the cellular environment.

For the sake of maximal protein expression, the vaccine construct's codon was optimized by JCat (Java Codon Adaptation Tool) in E. coli (strain K12). The JCat resource provided the percentage of average GC content up to 70% and CAI value (0.95–1.0) of adapted sequences, which ensured a high expression rate of vaccine constructs in E. coli. The adapted codon sequences were introduced in the plasmid vector pET 28a (+) for the construction of recombinant plasmid sequences, which ensured expression and heterologous cloning in E. coli. Cloning of the vaccine DNA sequence provided a recombinant vector of 6,804 bp (Supplementary Figure S2).

Immune simulations play a role in the development of an immune profile of the vaccine construct. The memory of immune cells was observed under C-ImmSim analysis to increase their half-life and successive immune responses by the cell (Bibi et al., 2021). Steadiness and real immune reactions were confirmed by the outcomes of this server. High IgM values indicated a primary immune response. Furthermore, the immunoglobulin expression level (IgG1+IgG2, IgM, and IgG+IgM) was increased, showing a decrease in antigen concentration (Figure 5A), which confirmed an increase in the active B-cell population (Figure 5D). Immunization led to an increase in the concentration of the active cytotoxic and helper T lymphocytes (Figures 5B, C). The cytokinin levels and IFN-gamma-inducing properties were also increased after every dose of the vaccine. Furthermore, the inset plot indicated that the chances of the danger caused by the vaccine are too low (Figure 5E). All the immune responses indicated that our vaccine is capable of causing a potent immune response.