The sequence of monkeypox from Pakistan (GISAID ID: EPI_ISL_17536780) was retrieved from GISAID. The nucleotide sequence is around 200 kb in size and encode 190 proteins. The antigenic potential of all the proteins were determined through VaxiJen web server (Doytchinova and Flower, 2007). The antigenic proteins were further subjected to determine the subcellular localization of proteins using Cell-Ploc-2 web server (Chou and Shen, 2008). The proteins that expressed on the cell membrane and in the extracellular region were further subjected for the identification of signal peptide and their secretory nature using Signal-BLAST web server (Frank and Sippl, 2008).

The shortlisted proteins were subjected for toxicity analysis through Toxinpred web server (Gupta et al., 2015). These proteins were also screened for the allergenicity profile through AllerTop web server (Dimitrov et al., 2013). The proteins were also subjected to BLASTp program for identification of human homologues with percent identity ≤35. The final proteins included in the analysis were antigenic, secretory in nature, non-toxic, non-allergic and non-homologues to humans. The shortlisted proteins were then used for the identification of B- and T-cell epitopes.

For identification of B-cell epitopes, Immune Epitope Database (IEDB) and ABCPred were used. The B-cell prediction tool at IEDB predicts the epitopes on the basis of protein flexibility, antigenicity, and surface accessibility (Fleri et al., 2017). The default threshold of 1.00 was used for the epitope prediction. The ABCPred predict the B-cell epitopes on the basis of recurrent neural network (RNN) (Saha and Raghava, 2006). The 16–18 amino acids long peptides with threshold of 1.0 were selected for further analysis. The epitopes found promising in the two platforms were consider for further investigation.

T cell (cytotoxic and helper T-cell) epitopes play an important role in eliciting immune response. The T-cell epitopes are recognized by histocompatibility complex (MHC) molecules as MHC-I and MHC-II.

The cytotoxic T-cell epitopes recognized by MHC-I were predicted through NetMHCpan4.0. NetMHCI 4.0 predicts the binding affinity through artificial neural network by schooling 81 distinct HLA-A, -B, -C and -E human MHC alleles (Reynisson et al., 2020). The recommended peptides that elicit MHC-I response were of nine residues long hence the 9-mer peptides with binding affinity ≤100 nM were selected. The T-helper cell epitopes were predicted through NetMHCIIpan4.0. It is based on artificial neural network and predict the binding affinity based on three human MHC class II isotypes HLA-DR, HLA-DQ, HLA-DP (Reynisson et al., 2020). The peptide length for binding of epitopes were of 13–17 amino acid in length hence the 15 amino acid long peptides with binding affinity ≤100 nM were selected.

The selected MHC-II peptides were further scrutinized for their potential to elicit IFN-γ response using IFN epitope server ² . The server predicts IFN-γ response on the basis of support vector machine (Dhanda et al., 2013). The peptides with the positive values were considered to be capable of inducing IFN-γ production.

All the shortlisted B- and T-cell epitopes were subjected for further evaluation for their potential antigenicity, toxicity, and allergenicity potential. The peptides that were antigenic (predicted through VaxiJen) were further evaluated for toxicity through Toxinpred and the non-toxic peptides were further evaluated for allergic potential using AllerTop (Dimitrov et al., 2014). The BLAST-p program was used to identify the human homologues. The non-human homologues with identity <30% were selected for final vaccine construction.

In order to determine which population would get benefits from vaccine epitopes, it is important to predict the population coverage of selected epitopes. To predict the population coverage, we used the Population Coverage ³ tool available at IEDB online server was used (Bui et al., 2006). Final selected MHC-I and MHC-II data were used to identify the population coverage across all the regions of the World.

The immunogenicity of peptide vaccines can be low, therefore, a multi-epitope peptide vaccine can overcome this problem. Hence, a multiepitope vaccine is created by fusing the overlapping B- and T-cells epitopes with the adjuvant and linkers. Linkers render vaccine constructs stable and permit them to behave independently after being injected into the host while adjuvant augments the immunogenic response (Ismail et al., 2022). Three different vaccine constructs were designed using 50S ribosomal protein L7/L12 adjuvants, TLR4, and β-defensin. The sequence of adjuvants were downloaded from Uniprot. The EAAK linker was used to join the adjuvant at the N-terminus of vaccine construct. The B-cell epiotpes were linked through KK linker while T-cell epitopes were joined through GPGPG linker. The physicochemical properties of the vaccine constructs like molecular weight, GRAVY, PI, stability index and half life were evaluated through ProtParam tool available on ExPASy server (Wilkins et al., 1999). The solubility of vaccine constructs was evaluated through SOLpro web server (Magnan et al., 2009). All the four vaccine constructs were analyzed for their antigenicity, toxicity, allergenicity, and autoimmunity using VacxiJen, Toxinpred, AllerTop and BLAStp Programs.

The secondary structure of the final vaccine construct was predicted via PSIPRED (Jones, 1999). The tertiary structure was predicted through trRosetta and I-TASSER (Du et al., 2021). I-TASEER is an iterative threading program that builds 3D structure by using the hierarchical method. trRosetta is de novo structure prediction tool based on deep neural network. Each method generate five models of each vaccine construct. The models were evaluated using Ramachandran plot, and ERRAT plot (Colovos and Yeates, 1993; Lovell et al., 2003). The protein model with the highest quality factor and accurate geometry according to Ramachandran plot was further undergo refinement using DeepRefiner server (Shuvo et al., 2021).

Molecular dynamics simulations were conducted for all three vaccine constructs using the GROMACS 2024.3 package, with each construct subjected to a 100 ns simulation. The structures were initially minimized using the Amber force field. Energy minimization was performed for 1,000 steps using the Steepest Descent Method, achieving convergence with a maximum force <1,000 (KJ mol⁻¹ nm⁻¹) to eliminate steric clashes. The system underwent equilibration under NVT and NPT ensembles for 100 ps (50,000 steps) and 1,000 ps (1,000,000 steps), respectively, with time steps of 0.2 and 0.1 fs at 300 K, ensuring the systems were fully equilibrated for the production run. The 100 ns production run was carried out at a constant temperature of 300 K and pressure of 1 atm (NPT) using weak coupling velocity-rescaling (modified Berendsen thermostat) and Parrinello-Rahman algorithms. The simulation results were analyzed based on root mean square deviation (RMSD), root mean square fluctuation (RMSF), and radius of gyration (RoG).

In order to generate immune response, the vaccine construct has to interact with the immune cell receptor. Hence, the docking of vaccine construct was studied with the Toll-like receptor 4 (TLR4), Toll-like receptor 3 (TLR3) and Toll-like receptor 7 (TLR7)to evaluate vaccine-immune receptor interactions. Structure of TLR4 (PDB ID: 4G8A), TLR3 (PDB ID: 1ZIW), and TLR7 (PDB ID: 5GMF) were downloaded from protein data bank. All the receptor structures were prepared for docking using Molecular Operating Environment (MOE) software. Briefly, the heteroatoms and co-crystallize ligands were removed followed by energy minimization using Amber12 force field. Molecular docking of all three multi epitope vaccines with TLR4 was performed using ClusPro server (Kozakov et al., 2017). The server performs docking of two interacting proteins. The server uses the rigid body docking approach by exploring billions of initial conformations. After that 1,000 lowest energy structures are shortlisted on the basis of RMSD value in order to find the best cluster that represent the most likely complex. Finally, the refinement of selected complex was performed through energy minimization (Kozakov et al., 2017). The binding flexibility and stability of vaccine-TLR4 complex was predicted through iMODS server. The server utilizes an elastic network model (ENM) to perform Normal Mode Analysis (NMA) in internal (dihedral) coordinates (Lopez-Blanco et al., 2014). For each complex the server provide deformability plots and B-factor, mode variance plot, eigenvalues, elastic network, and covariance map.

To evaluate the immune profile of vaccine construct, computational immune simulation was performed using online C-ImmSim server. This server uses a position-specific scoring matrix (PSSM) and machine learning techniques for the prediction of the immune response (Rapin et al., 2011). The simulation parameters were set at default for a period of 1, 84, and 168 h with default host HLA MHCI and MHCII alleles. The simulation was run for 1,000 steps.

Codon optimization is an important step in reverse vaccinology in-order to avoid the codon bias which in turn enhance the expression of recombinant protein. The Java codon adaptation index (JCAT) ⁴ was used to reverse translate amino acid sequence of vaccine into cDNA for effective expression of vaccine construct in E.coli vector. It provides the Computed Codon Adaptation Index (CAI) values and GC content of the vaccine construct (Grote et al., 2005). These two parameters are important to evaluate the expression of vaccine. The optimal transcriptional and translational efficiency of vaccine achieved with CAI >0.8 and GC content between 30%–70%. The adapted codon sequence was cloned into pET-28a (+) vector using SnapGene ⁵ . This vector enables enhanced protein recovery and purification due to its N-terminally 6 × His-tagged proteins.

On the basis of B-cell epitope prediction server, 16- and 18-mer epitopes were predicted for all the 10 proteins. The epitopes were shortlisted on the basis of their antigenicity and allergenicity. No epitope was predicted for EGF-like domain protein. The epitopes predicted for Envelope protein A28 homolog and Crm-B secreted TNF-alpha-receptor-like protein were non-antigenic and allergen. Hence, not included in the study. 12 epitopes were prioritized from seven proteins based on their high antigenicity, low toxicity, and low allergenic reactions (Table 2).

A total of 122 CTL epitopes were predicted with binding affinity ≤100 nM with MHC-I alleles. Among these epitopes, 33 epitopes were predicted to be antigenic, non-toxic and non-allergen (Table 3). A total of 44 helper T-cell (HTL) epitopes were identified with binding affinity ≤100 nM with MHC-II alleles. Out of 44 HTL epitopes, 12 epitopes satisfied the antigenicity, toxicity, and allergenicity tests (Table 4).

The shortlisted, 12 B-cell, 33 CTL and 12 HTL epitopes were further shortlisted on the basis of overlapping features. 6 B-cell, 8 CTL and 6 HTL overlapping epitopes were used to design vaccine construct (Table 5). The eventual target has been to determine lead epitopes with the potential to prompt humoral and cell mediated immunogenic responses. The selected T-cell epitopes were further subjected to identify the population coverage. The IEDB results revealed that the selected epitopes showed 97.5% worldwide coverage of both the MHC I and MHC II alleles. Population coverage in different regions of world is shown in Figure 1.

To design vaccine construct 6 B-cell, 8 CTL and 6 HTL epitopes were used using EAAK, KK and GPGPG linker. Three different vaccine constructs were designed based on three different adjuvants. To enhance protein purification, a 6xHis tag was added to the C terminal. The adjuvants were added to N-terminal to construct the vaccine. On the basis of adjuvant, three vaccine constructs were generated as V1 (adjuvant: beta defensin), V2 (adjuvant: 50S ribosomal protein L7/L12) V3 (adjuvant: TLR4 agonist).

Immunological profile of all the three vaccine constructs showed that these are antigenic, non-toxic and non-allergic. The antigenic score of all the three-vaccine construct >0.8 indicating the antigenic nature of vaccine constructs. The SOLpred score of all the three vaccine constructs >0.6 indicating the solubility upon overexpression (Table 6). The physiochemical properties of all the three vaccine constructs is shown in Table 7. As all the vaccine constructs have similar protein sequence the only difference is the presence of adjuvant hence the physicochemical properties differ slightly among these three constructs. The molecular weight ranges from ∼20 kDa to ∼25 kDa. Theoretical pI values range from 8.8 to 9.3. The stability index lies in the range of 16–29, that indicate the stability of vaccine constructs. GRAVY scores from −0.425 to −0.583 indicate the hydrophilic nature of the vaccine constructs. Aliphatic index (72.2–86.8) describes the thermostability of vaccine constructs. These results indicate that all the three vaccine constructs may be able to initiate immunological response into the host cell.

In order to understand the interactions of vaccine with the host immune cells, stable three dimensional structure is required. Rosetta was used to generate the structure of all the three vaccine constructs and the resulting models were evaluated for structure stability. Five models were generated for each vaccine construct. Top scoring models were further subjected to structural refinement. The results of structural validation showed that ERRAT score >80. The ramachandran plot showed that 78.1% residues of V1, 85.2% residues of V2, and 85.7% residues of V3 fall in favored region of plot. Hence, these results indicate that the predicted 3D structure of vaccine constructs are of good quality and used for further investigation (Figure 2).

The RMSD, RMSF, and radius of gyration analyses for all three vaccine constructs indicate that V1 is the most stable (Figure 3). The RMSD of V1 increases during the first 50 ns but stabilizes thereafter. In contrast, the V2 construct remains unstable throughout the 100 ns simulation, with its RMSD continuing to rise. For V3, the RMSD stabilizes after 70 ns, but it remains higher than that of V1.

Regarding the radius of gyration, V1 maintains a value of around 1.8 nm, indicating compactness. V2 has a higher radius of gyration, approximately 1.95 nm, suggesting a more expanded structure. V3 fluctuates between 1.75 and 1.8 nm, showing moderate compactness, but still less stable than V1.

The binding of all the three vaccine constructs with TLR4 was evaluated using ClusPro. For each vaccine complex, Cluspro generated 120 structures. Top ten complexes of all the three vaccine constructs were subjected to protein binding energy prediction using PRODIGY. The server predicts the binding affinity in the complexes. On the basis of binding energy, the V1 vaccine construct showed the lowest binding energy among all the three complexes (Table 8). The interaction analysis revealed that 11 residues of V1 are directly involved in interactions with TLR4 (Figure 4).

All the three complexes were further subjected for deformability analysis. All the three complexes have relatively similar deformability in the residues (Figures 5A, E, I). Furthermore, the force required to deform the complex is measured through Eigenvalue. Lower the Eigenvalue, easier the structure will deform. According to graphs reported in Figures 5B, F, J, V1 complex have high eigenvalue compared to V2 and V3. Hence, the complexes of V2 and V3 with TLR4 are easier to deform. The covariance matrix also showed the correlated movements of atoms in case of V1 compared to V2 and V3 complexes as there is relatively more red region in V1 compared to V2 and V3 (Figures 5C, G, K). The elasticity of vaccine constructs is represented in Figures 5D, H, and L.

The interactions of V1 vaccine construct were further studied with TLR3 and TLR7 (Figure 6). The binding affinity of V1 with TLR3 is lower than with the TLR7 (Table 9). The interaction analysis revealed that there relatively strong interactions of V1 with TLR3 as compared to TLR7. The immune simulation results of V1 with TLR3 and TLR7 is also show a string immune response (Figure 7 ).

Final vaccine construct (V1) was subjected to codon optimization using JCat web server. To obtain high expression in the expression system, the peptide sequence of V1 was reverse translated into DNA. The GC content and predicted CAI value of vaccine construct was 48.4% and 1.0, respectively. Snapgene was used to perform the in silico cloning. The multiepitope vaccine construct exhibits PstI and RsrII restriction sites at the start and end, respectively. The pBAD-DEST49 vector was used for cloning as it exhibits both the PstI and RsrII sites. The multiepiotpe vaccine construct (V1) was cloned into pBAD-DEST49 vector in-order to ensure the expression in the E.coli expression system (Figure 8A). In order to check the insertion of peptide into vector, the vector was digested with PstI and RsrII enzymes. The virtual agarose gel simulation showed the presence of insert along with vector (Figure 8B).

To understand the immune response of vaccine construct, immune simulation was performed. The results showed that there was high level of IgM in the initial phase of vaccine administration followed by IgG1 and IgG2 (Figure 9A). The same trend observed in human body immune responses where the primary immune response is generated by IgM and secondary and tertiary immune response is generated by B-cell, IgG, and IgM. As the antibody response is generated there is decrease in the antigen level as well. The exposure of antigen in the vaccine construct also led to rapid increase in population of memory B-cells (Figure 9B). Along with the memory B-cell population there was an increase in cytotoxic and helper T-cell population (Figures 9C, D). Like real-time immunization, there was constant elevated levels of natural killer cells, dendritic cells, and macrophages (Figures 9E–G). Increased level of cytokines and interleukins were also observed in response to vaccine administration (Figure 9H). These results suggest that the vaccine construct (V1) can generate immune response inside the human host.