Design peptide and multi-epitope protein vaccine candidates against monkeypox virus using reverse vaccinology approach: an in-silico study

Abstract Monkeypox is a zoonotic virus that has recently affected different countries worldwide. On July 23, 2022, the WHO declared the outbreak of monkeypox as a public health emergency of international concern. Surveillance studies conducted in Central Africa in the 1980s and later during outbreaks in the same region showed smallpox vaccines to be clinically somewhat effective against Monkeypox virus. However, there is no specific vaccine against this virus. This research used bioinformatics techniques to establish a novel multi-epitope vaccine candidate against Monkeypox that can induce a strong immune response. Five well-known antigenic proteins (E8L, A30L, A35R, A29L, and B21R) of the virus were picked and assessed as possible immunogenic peptides. Two suitable peptide candidates were selected according to bio-informatics analysis. Based upon in silico evaluation, two multi-epitope vaccine candidates (ALALAR and ALAL) were built with rich-epitope domains consisting of high-ranking T and B-cell epitopes. After predicting and evaluating the 3D structure of the protein candidates, the most efficient 3D models were considered for docking studies with Toll-like receptor 4 (TLR4) and the HLA-A * 11:01, HLA-A*01:01, HLA-A*02:01, HLA-A*03:01, HLA-A*07:02, HLA-A*15:01, HLA-A*30:01 receptors. Subsequently, molecular dynamics (MD) simulation of up to 150 nanoseconds was employed to assess the durability of the interaction of the vaccine candidates with immune receptors. MD studies showed that M5-HLA-A*11:01, ALAL-TLR4, and ALALAR-TLR4 complexes were stable during simulation. Analysis of the in silico outcomes indicates that the M5 peptide and ALAL and ALALAR proteins may be suitable vaccine candidates against the Monkeypox virus. Communicated by Ramaswamy H. Sarma


Introduction
Monkeypox (MPX) is a rare zoonotic viruses that was first discovered in the 1970s in the Democratic Republic of the Congo (DRC) (Orviz et al., 2022).On May 20th, 2022, the cumulative number of reported and suspected cases increased above 100.Since October 28, 2022, a total amount of 109 nations had reported monkeypox cases, with 77,115 laboratory-confirmed cases and 36 fatalities (Silva et al., 2022).Geographically, 50,573 cases were reported in the Americas, 25,297 in Europe, 934 in Africa, 209 in the Western Pacific, 72 in the Eastern Mediterranean, and 30 cases in South-East Asia (Almehmadi et al., 2022;Petersen et al., 2022).Today, monkeypox is therefore thought of as a reemerging contagious virus with the potential to spread out around the world for a long duration, similarly to the coronavirus disease which arised in 2019 , and still is big concern worldwide.
Presently, although monkeypox is not a recently arising infection, a variety of gaps in its understanding doesn't allow for the advancement of efficient disease control.Monkeypox virus (MPXV) transmission to humans can occur from animals or humans by respiratory droplets, body fluids, sexual contact, direct contact, or eating a contaminated animal (Hraib et al., 2022).Furthermore, it has been shown that MPXV can transmit between animals via aerosols (Alakunle & Okeke, 2022;Mart� ınez et al., 2022).Its symptoms are similar to smallpox but more severe with fever, headache, lymphadenopathy and a mortality rate of about 1 to 10% (Doshi et al., 2020;Ogoina & Yinka-Ogunleye, 2022).As per WHO recommendations, nucleic acid amplification test (NAAT) has been performed to detect MPXV virus.Genes commonly used for routine PCR testing include hemagglutinin, the acidophilus gene, and the crmB gene.For most patients, treatment is symptomatic and supportive.There is no specific treatment for MPXV.Antiviral drugs for smallpox, such as brincidofovir, tequovirimate, and cidofovir, may affect MPXV because of their similar genetics (Huang et al., 2022;Li et al., 2017).
Up to now, three vaccines (ACAM 2000, LC16m8 and JYNNEOS (also known by the brand names Imvamune and Imvanex)) are used to protect against the smallpox virus.These vaccines have an 85% efficiency versus human monkeypox infection as a result of resemblances between the two large double-stranded DNA orthopoxviruses (Gong et al., 2022;Sah et al., 2022).
ACAM2000 is a live vaccine, JYNNEOS is based on a live, attenuated non-replicating orthopoxvirus, and LC16m8 is a third-generation vaccine that contains a virus from the Lister strain used in first-generation vaccines.Multiple passages in tissue culture and selection for a weak phenotype resulted in strain LC16m8, which lacks the full-length, functional B5 membrane protein (Forbes et al., 2022;Saijo et al., 2006).The US Food and Drug Administration (FDA) granted provisional approval of the JYNNEOS and ACAM2000 vaccines for the 2022 MPXV outbreak.The results of Meo and colleagues' research showed that both vaccines are useful, but the overall effect of JYNNEOS is better than that of ACAM2000 (Meo et al., 2022).It has already been demonstrated that ACAM 2000 can lower the symptoms of MPX; however, it has adverse effects like atopic dermatitis in immunocompromised patients (Brown & Leggat, 2016).Also, JYNNEOS is capable of preventing smallpox virus and MPXV for use in high-risk adults without showing atopic dermatitis (Petersen et al., 2019).A 3rd authorized vaccinia vaccine is made by KM Biologics and is manufactured and offered in Japan only for USA against bioterrorism (Lozano & Muller, 2023).None of the above mentioned vaccines however, have been approved for the general populace and also they were developed with older technologies (Haider et al., 2022).So, because of the lack of a reliable vaccine for MPXV, there is a need to create a specific vaccine against this infection.In vaccine development, computational methods can decrease the time and expenses of finding the right epitope areas from different parts of the infection pathogen (Jahantigh et al., 2022;Khan & Pichichero, 2014).Using in silico methods, it is possible to find epitopes that can be identified by T and B lymphocytes (Dey, Mahapatra, Lata, et al., 2022).
The MPXV has linear, double-stranded DNA with inverted terminal repeats (ITRs) at its ends (Fantini et al., 2022).MPXV encodes a range of various host proteins that regulate cellspecific antiviral responses to ensure the virus can replicate in some cells, such as A8L, A30L, A35R, A29L, and B21R (Dupuy & Schmaljohn, 2009;Franceschi et al., 2015;McInnes et al., 2006;Meyer et al., 2005;Townsend et al., 2017).For instance, the B21R glycoprotein of MPXV, specifically found only in MPXV, has a remarkable immunogenic feature and is used to differentiate monkeypox from previous vaccination and other orthopoxviruses (Isidro et al., 2022).Also, the MPXV A29L protein is a recombinant protein expressed in mammalian cells containing a C-terminal His-tag, is associated with cell fusion and plays an essential role in virus penetration (Dumont et al., 2014).
The present study intends to utilize a reverse vaccinology technique to find immunogenic epitopes from essential virulence proteins of MPV, particularly E8L, A30L, A35R, A29L, and B21R proteins, to introduce suitable vaccine candidates against MPV.Finally, we have proposed two peptide vaccine candidates and two recombinant protein candidates against MPV.By computational methods, we assess the interaction between the recombinant proteins and peptides with the TLR4 and HLA-A � 11:01, HLA-A � 01:01, HLA-A � 02:01, HLA-A � 03:01, HLA-A � 07:02, HLA-A � 15:01, HLA-A � 30:01 immune system receptors to show the biological activity of the proposed protein and peptide candidates.Finally, using in silico immune simulations, the effect of protein vaccine candidates on the induction of immune responses was analyzed.

Materials and methods
In this study, we designed peptide and multi-epitope protein vaccine candidates against MPVX, by exploiting the programs of reverse vaccinology (Figure 1).

Retrieval of the virus protein sequences and antigen prediction
The experiment was performed on the MPVX proteome, focusing on some of the main virulence proteins of MPVX, including E8L, A30L, A35R, A29L, and B21R, for which the immune action versus these viral proteins is highly likely to show neutralizing effects (Kugelman et al., 2014).Entire amino acid (aa) sequences of E8L, A30L, A35R, A29L, and B21R were obtained from NCBI (https://www.ncbi.nlm.nih.gov).The VaxiJen v2.0 server based on the automatic crosscovariance (ACC) was used for predicate antigenicity (Doytchinova & Flower, 2007, 2008).To avoid false positive results and increase accuracy, a threshold of 0.4 was considered (Doytchinova & Flower, 2007).

Prediction of MHC class I epitopes of CD8 þ T cells
To predict CD8þ T cell-restricted MHC class I epitopes, the protein sequences of E8L, A30L, A35R, A29L, and B21R in FASTA format were submitted to Immune Epitope Database Analysis Resource (IEDB-AR).This database uses an artificial neural network (ANN) method to identify epitopes for B and T cells (Narang et al., 2021;Nielsen et al., 2003).Based on previous reports, 7 MHC class I alleles which cover 85-90% of the world's population, were analyzed to predict epitopes (Friend et al., 2009).

Prediction of CD4 þ T cell MHC class II epitopes
The ProPred and IEDB-AR servers were used to predict CD4þ T cell MHC class II epitopes (Manici et al., 1999;Singh & Raghava, 2001).At first, epitope identification of the 15-mer sequence was performed by IEDB-AR using the consensus method (Bui et al., 2005;Nielsen et al., 2007).ProPred was also used for MHC class II 9-mer epitope prediction.ProPred is a visual server that predicts MHC class II binding regions using a matrixbased prediction formula (Dey, Mahapatra, Raj, et al., 2022;Nielsen et al., 2007).Finally, we selected the overlapping 9-mer peptides predicted by both different servers.

Prediction of B cell epitopes
Prediction of B epitopes was made using IEDB-AR and ABCpred servers.The ABCpred server uses the Jordan network and features a surprise layer with optional sizes in a single hidden layer for B epitope prediction (Saha & Raghava, 2006).
The IEDB-AR server using the BepiPred method was used to predict direct B-cell epitopes (Dey, Mahapatra, Patnaik, et al., 2022;Ponomarenko & Bourne, 2007).In the following step, overlapping peptides were kept for further analysis.

3D structure modelling and characterization of overlapping epitopes discovered on both B and T cells
A peptide solubility calculator (https://pepcalc.com/)was used to evaluate the solubility of peptides.The antigenicity and allergenicity of overlapping epitopes were predicted using VaxiJen and AllerTOP servers, respectively (Dimitrov et al., 2014;Doytchinova & Flower, 2008).Epitope toxicity prediction was performed using the ToxinPred server (https://webs.iiitd.edu.in/raghava/toxinpred/multi/submit.php)(Gupta et al., 2013).The sequence of the epitopes was entered and other settings remained at their default values.

Homology modelling and peptide-allele docking
Overlapping epitopes of MPXV proteins recognized by T and B lymphocytes were selected for evaluating the interaction with immune system receptors.FASTA formats of selected epitopes were submitted to the PEP-FOLD server for three-dimensional (3D) structure prediction (Beaufays et al., 2012).Only structures with the ideal points, i.e. the most points, were selected.Energy efficient sOPEP (i.e.highest value) was used as the theme (Berman et al., 2003).The molecules' crystal structure of HLA-A � 11:01 (PDB ID: 4uq2), HLA-A � 01:01 (PDB ID: 4nqx), HLA-A � 02:01 (PDB ID: 4u6y), HLA-A � 03:01 (PDB ID: 7l1c), HLA-A � 07:02 (PDB ID: 6at5), HLA-A � 15:01 (PDB ID: 6vb3), HLA-A � 30:01 (PDB ID: 6j1w) receptors were retrieved from the Protein Database Banking (PDB).The overlapping epitopes and alleles were prepared for docking using Cluspro 2 (Desta et al., 2020;Mahapatra et al., 2022).The binding affinity and dissociation constant (Kd) of peptide-protein and protein-protein complexes were predicted using the protein binding energy prediction web server (PRODIGY).Finally, the best complexes will be selected for molecular dynamics simulation based on binding energy, number of hydrogen and hydrophobic bonds, and binding affinity.The PyMOL program (www.pymol.org/)was utilized to illustrate the predicted epitope structure, the setting of these epitopes and accessory regions of HLA-I for the predicted peptides (DeLano, 2002;DeLano & Bromberg, 2004).

Design of multi-epitope protein vaccine candidate
After analysing the results obtained from B and T cell epitopes, epitope-rich regions were identified and linked to each other by appropriate linkers.(GGGS) 3 , (GGGS) 2 EAAAK, and (EAAAK) 2 linkers were evaluated for the desired structures.We examined the epitope-rich domains with different linkers and different states in terms of structural conformation and finally selected the best and most stable structure as the final constructs.

Reviewing solubility and other physicochemical specifications
The protein-sol software was used to analyze the candidate protein's solubility (Nielsen et al., 2007).On the data available for the solubility of the Escherichia coli protein, the protein-sol server uses 35 sequence-based features to evaluate the solubility of a protein (Doytchinova & Flower, 2007).

Secondary and tertiary structure prediction
The second structure prediction of the candidate proteins was performed using the Garnier-Osguthorpe-Robson (GOR) server based on the aa sequence (Garnier et al., 1978).The 3D structure of the candidate proteins was predicted using the Galaxy Web, I-TASSER and SWISS-MODEL servers (Bharti et al., 2019;Shin et al., 2014).Finally the best structure was chosen.
Recognizing errors in experimental and theoretical models of candidate proteins structures is one of the essential problems.The designed structures should be close to the structure of naturally expressed and functional proteins in cells (Wiederstein & Sippl, 2007).The three-demential structures were evaluated using ProSA-web to check possible errors in the designed 3D models.This server evaluates the structure based on Z-score and energy plots.An out-of-range Z-score for native proteins indicates an error in the 3D structure (Wiederstein & Sippl, 2007).One of the most critical programs for evaluating the tertiary structure is calculating the Ramachandran diagram.In this study, the Ramachandran plot obtained from PROCHECK was used to calculate the torsion angles of the residues in the candidate proteins and indicate whether the residues are in outlier, allowed, and favoured regions (Gopalakrishnan et al., 2007).Also, MolProbity analysis was used to validate the quality of three-dimensional structures of macromolecules such as proteins, nucleic acids, and complexes (Chen et al., 2010).MolProbity analysis results include rotamer analysis, MolProbity score, clash score, and Ramachandran plot.

Conformational B cell epitope prediction
Structural epitopes of the candidate proteins were predicted using the ElliPro server (Ponomarenko et al., 2008).This web server uses a modified version of Thornton's method, MODELLER program, and Jmol visitor to predict direct antibody epitopes.The ElliPro server predicts conformational Bcell epitopes using designed proteins structure.

Molecular docking analysis with immune receptors of the developed protein
ClusPro 2.0 has been used to predict the interactions between the TRL4 (PDB ID: 3FXI) and the designed proteins (Desta et al., 2020).The TLR4 receptor was picked and downloaded from the PDB server.At first, the receptors were cleaned, and crystallographic water and other chemicals were erased.All these procedures were carried out using the PyMOL v2.3.4 software program.The ClusPro 2.0 server is based upon PIPER, an FFT correlation strategy for protein docking with pairwise communication capacity.Protein binding energy prediction (PRODIGY) webserver was used to predict binding affinity and dissociation constant (Kd) of receptor-ligand complexes.Finally, interactions were evaluated to analyze hydrogen bonds and hydrophobicity using LigPlot þ v. 4.5.3(Wallace et al., 1995).

Simulation of molecular dynamics
To stabilise and relax the structures in an aqueous physiological environment, the TLR4 and HLA-A � 1101 receptors, free ALAL and ALALAR proteins were simulated for 150 nanoseconds (ns).The peptide-receptor and protein-receptor docked complexes were picked for MD simulations done with the GROMACS 2018 bundle and OPLS-AA force field (Abraham et al., 2014).The M5-HLA-A � 1101 and M6-HLA-A � 1101 peptide complexes and ALAL-TLR4 and ALALAR-TLR4 protein complexes were inserted in a 10A˚ solvent box full of simple point charge (SPC) water molecules.The system charge became neutralized by introducing enough Na þ and Cl.The energy minimization method for all MD simulations consisted of two parts: systems were equilibrated for 100 picoseconds (ps) using a constant Number of particles, Volume, and Temperature (NVT) and a constant Number of particles, Pressure, and Temperature (NPT).The system was equilibrated using a Parrinello-Rahman barostat to achieve a constant pressure and temperature of 1.0 bar and 300K.The Particle Mesh Ewald (PME) method was used to treat longrange electrostatics with a 0.16 nm grid spacing and a 10 Å cut-off distance.With a cut-off of 1 nm, the Van der Waals (VDW) interactions were calculated.Finally, Molecular Dynamics simulations were used for each simulation system (for a duration of 150 ns) with two femtoseconds (fs) time steps).To examine the stability of the peptide-protein and protein-protein complexes, RMSD, hydrogen bonds (H-bonds), the distance of gyration (RG), RMSF, and solventaccessible surface (SASA) were determined (Enayatkhani et al., 2021;Mobini et al., 2020;Sahoo et al., 2022).

Molecular technicians Poisson-Boltzmann solventaccessible area (MMPBSA) binding-free power
The relative free binding energy in the M5-HLA-A � 1101 and M6-HLA-A � 1101 peptide complexes were calculated using the Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) method developed by Rashmi kumari (Kumari et al., 2014).In this method, the free binding energy was calculated by MMPBSA according to the following formula:

Immune simulation
Analysis of the immune response against the designed protein candidates were performed by C-ImmSim online simulation server (http://150.146.2.1/C-IMMSIM/index.php)(Rapin et al., 2010).The C-ImmSim server uses Miyazawa and Jernigan's protein-protein dimensions to analyze the molecular association of vaccine candidates with immune system receptors in animals' humoral and cellular immune complications.The simulation was set for three injections with time steps of 1, 84, and 168 and system defaults, each time step equalling 8 h (Castiglione et al., 2012;Singh et al., 2020).The simulation steps were 1000, the simulation quantity was 10 mL, and the random seed was 12,345.

In silico cloning of the candidate vaccine
The SnapGene device was used to insert the designed codon sequences of the developed proteins in the pET28a (þ) vector for limitation cloning and informatics conformation of the vaccine candidates by double digestion.NcoI and XhoI restriction sites were incorporated into the designed vaccine sequence at the N-and C-terminal sites, respectively.

CD8þ T cell MHC I epitope identification
The IEDB-AR server predicted 114 peptides for the four studied MPX proteins.Finally, 36 peptides with the most efficient capacity based on high affinity score and PIR results were selected for further study, which showed that they did not induce autoimmune response and were soluble in water (Supplementary Table S1).

CD4þ T cell MHC II epitope identification
ProPred and IEDB-AR servers were chosen to predict MHC class II epitopes.In the initial evaluation, 180 epitopes were selected.Finally, 111 epitopes were selected for further evaluations based on high affinity scores and PIR results, which indicated that they might not activate an autoimmune activity (Supplementary Table S2).

B cell epitope prediction
It is also necessary to stimulate the production of the humoral immune response against the virus.B cell epitopes contribute to immunity and antibody production against infections.Supplementary Table S3 shows the 62 overlapping B-cell epitopes identified using IEDB-AR and ABCpred.

3D structure modelling and characterization of B and T cells overlapping epitopes
The 7 epitopes with overlapping MHC class I and class II T and B cells were forecasted (Table 2).The evaluation of the outcomes reveals that M5 and M6 epitopes have the greatest antigenicity, are not allergenic and are soluble.

Molecular docking studies for HLA and peptide interaction analysis
The crystal structure of HLA-A � 11:01, HLA-A � 01:01, HLA-A � 02:01, HLA-A � 03:01, HLA-A � 07:02, HLA-A � 15:01, HLA-A � 30:01 receptors were taken from the PDB data source.The most ideal binding energy and the variety of residues associated with hydrogen and hydrophobic bonds of the complexes were considered as the major standards for picking the most ideal complexes (Figure 2).The results of weighted score of molecular docking, Kd and affinity binding obtained from PRODIGY are shown in Table 3. Analysis of results  shows that M5 and M6 peptides have strong interaction with the HLA-A � 11:01, HLA-A � 01:01, HLA-A � 02:01, HLA-A � 03:01, HLA-A � 07:02, HLA-A � 15:01, HLA-A � 30:01 receptors.

Design of a multi-domain protein as vaccine candidates
The epitope-rich domains were chosen to build the last constructs as protein vaccine candidates.These domains are strongly immunogenic and cover both B and T cell epitopes (Table 4).The most ideal three-dimensional structures were predicted by (GGGGS) 3 linker (Figure 3A).A schematic diagram of the multi-epitope proteins including the immunodominant domains of A30L, A35R, and A29L proteins connected by (GGGGS) 3 linker is displayed in Figure 3(B).The structure of ALALAR protein was predicted by the I TASSER server and ALAL protein by the SWISS-MODEL server and they had the best score in terms of evaluating the tertiary structure.

Antigenicity and allergenicity evaluation
The antigenicity of candidate proteins was calculated by the VaxiJen v2.0 server.The ALALAR and ALAL proteins had an antigenicity of 0.5926 and .5921,respectively, for the viral model, indicating that our designed candidates are antigenic and can induce immune responses.We also analyzed the allergenicity of the desired proteins by the AllerTOP server.
The result of our evaluation showed that the designed proteins are not allergenic and do not cause IgE-dependent inflammatory reactions.

Evaluating the solubility and other physicochemical parameters
The ProtParam tool was used to evaluate the physicochemical parameters of the vaccine candidates.The results of this evaluation showed that the ALALAR multi-epitope protein with 283 aa has a molecular weight of 29957.38Da.The aliphatic index was predicted as 72.76 and GRAVY was À 0.270.
The theoretical IP was determined as 7.1 for ALALAR protein.
Also, ALAL protein had 152 aa and a molecular weight of 16731.80Da.The aliphatic index was predicted as 69.93 and GRAVY was À 0.520.The theoretical IP was determined as 8.98 for this protein.We evaluated the solubility of the multi-epitope candidates using the Protein-sol software tool.
The solubility of the ALALAR and ALAL protein candidates were determined to be 0.434 and 0.54 that indicate ALAL protein have a higher solubility than the average soluble E.coli protein (0.45) from the experimental solubility dataset and ALALAR is predicted to be less soluble (Niwa et al., 2009).

Secondary and tertiary structure prediction and validation
Prediction of the secondary structure showed that the ALALAR vaccine candidate was made up of an alpha-helix (21.2%), an extended strand (31.1%) and a random coil (47.7%) (Figure 4A), while the ALAL protein was made up of an alpha-helix (27.63%), an extended strand (25.00%), and a random coil (47.37%) (Figure 4B).The predicted structures were verified by ProSA-web, MolProbity, and Ramachandran plot to pick the most ideal and most consistent 3D structure closest to the native protein structure (Table 5).MolProbity evaluation showed the clash score of ALALAR and ALAL proteins were 10.14 (71% of the best among structures) and 11.99 (63% of the best among structures) respectively.The ProSA-web reviewed a 3D model of the vaccine candidate utilizing an energy plot and Z-score.The Z-score of the ALALAR (Figure 4C) and ALAL (Figure 4D) proteins was À 4.74 and À 3.59, respectively, which remains within the variation of the native protein structure.The Ramachandran plot evaluation for ALALAR protein showed 98.7% residues in the favoured and allowed regions and 1.3% in the disallowed regions (Figure 4E).100% of residues of ALAL protein located in favoured and allowed regions (Figure 4F).

Defining discontinuous B cell epitopes
The 3D structure of the ALALAR and ALAL proteins as vaccine candidates were evaluated by the ElliPro server.The results of this analysis revealed regions of the protein as structural epitopes for ALALAR (Figure 5) and ALAL proteins (Figure 6).B cell conformational epitopes were predicted with scores varying from 0.722 to 0.871 for ALALAR protein (Table 6) and from .755 to .904 for ALAL protein (Table 7).

Protein-protein molecular docking
Protein-protein molecular docking was conducted between the ALALAR and ALAL proteins with the TLR4 receptor.We used cluster 2.0 to pick the best interaction based on the weighted score and members.Additionally, hydrogenic and hydrophobic bonds between vaccine candidates with TLR4 were pictured utilizing the LIGPLOT program.We considered energy, the highest clusters, the most negative weighted result, the lowest  PRODIGY (Kd) score, and the greatest number of residues involved in hydrogen and hydrophobic intermolecular binding as the crucial standards for choosing the most potent complexes.The results suggested a strong interaction between the ALALAR and ALAL proteins with the TLR4 (Table 8).The interactions between the TLR4 binding and the docked ALALAR candidate (Figure 7) and ALAL (Figure 8) were visualized using the PyMOL and LIGPLOT programs.

Simulation of molecular dynamics
The most ideal interactions of the peptide-protein HLA-A � 1101 receptor that is the most common human allele (Habel et al., 2022) (M5-HLA-A � 1101 and M6-HLA-A � 1101) and protein-protein (ALALAR protein-TLR4 and ALAL protein-TLR4) docked complexes were selected for MD simulation.To stabilize and relax the structures in an aqueous physiological environment, each of the unbonded receptors and designed proteins were simulated for 150 ns.After 150 ns of simulation, the RMSD values of the unbound receptors showed that the backbone structures are stable (Figure 9).HLA-A � 1101 receptor had relatively strong fluctuations up to 70 ns, but it became stable after about 70 ns.The ALALAR and ALAL protein structures also reached stability after a slight conformational change of about 8 ns after the start of the simulation.The fluctuations of ALALAR and ALAL proteins after 8 ns of the simulation were less than .2nm, which indicates the stability of these two proteins (Figure 9).Unbonded HLA-A � 1101, TLR4 and free protein structures were visually reviewed and considered sufficiently stable for the next phase of the project.The RMSD values of complexes of HLA-A � 1101-M5 and HLA-A � 1101-M6 had an average of .99 and 1.09 nm, respectively (Figure 10A).The RMSD graph shows more stability of HLA-A � 1101-M5 complex compared to HLA-A � 1101-M6 complex.
The ALALAR-TLR4 and ALAL-TLR4 complexes had average RMSD values of .83 and .81nm respectively.The ALAL-TLR4 complex reached stability after about 60 ns.ALALAR-TLR4 complex also reached stability after about 25 ns (Figure 10A).Hbonds are essential for the interactions between the proteinpeptide and protein-protein complexes.The folding and compressibility of the complexes were reviewed making use of Rg.
The values of Rg in all four complexes also revealed that the structures of the complexes are stable throughout the simulation and the receptor and ligand are not divided (Figure 10B).H-bond formation in the protein plays an essential role in its stability.The H-bond diagram of the complexes showed that throughout the MD simulation the bond between the HLA-A � 1101 and TLR4 receptors and the ligands is relatively stable and indicates the binding stability for the complexes (Figure 10C).Additionally, we decided to carry out an RMSF analysis to examine the variations of backbone atoms of unbonded receptors and complexes.In this analysis, the average variation value of each residue throughout the simulation was plotted.As shown in Figure 11, the RMSF values show slight fluctuations for the majority of residues in complex styles compared to unbonded receptors.Compared to the TLR4 and HLA-A � 1101 unbonded receptors, the peptide and protein complexes are more stable and have a lower average RMSF.

MMPBSA binding-free energy study
Binding energies were analyzed using the MMPBSA tool to evaluate the molecular interactions of HLA-A � 1101-M5 and HLA-A � 1101-M6 complexes.Electrostatic energies, polar solubility, van der Waals, SASA, streptavidin (SAV) and binding energies of both complexes were calculated by MMPBSA technique.As can be seen in Table 9, the van der Waals and electrostatic energies of the complexes are solid enough to maintain the M5 and M6 peptides in contact with the HLA-A � 1101 receptor.However, considering that the binding energy of M5-HLA-A � 1101 complex is negative, it indicates a stronger binding of this peptide to the HLA-A � 1101 receptor.

Immune simulation
The C-ImmSim server was used to simulate the immune system response to the designed vaccine candidates.In silico  In addition, increased levels of IgG1, IgG2, IgG, and T and B cell populations indicate a secondary response to the ALALAR and ALAL antigens.The focus of cytokines and interleukins were also considerably enhanced after immunization.The resulting evaluation indicated that the ALALAR protein (Figure 12) stimulated stronger immune response compared to the ALAL protein (Figure 13).

In silico cloning of the candidate vaccine
We plan to clone the ALALAR and ALAL constructs into the pET-28a (þ) vector to use for protein expression.For this reason, we cloned the designed sequences in pET-28a (þ) vector using SnapGene software (Figure 10).For virtual confirmation of the clone, SnapGene agarose gel cloning method was performed.The presence of vaccine candidate genes (849 bp for ALALAR protein and 468 bp for ALAL protein) along with pET-28a (þ) vector (5231 bp) after double digestion with NcoI and XhoI enzymes confirms the protein candidate gene clone in the pET-28a (þ) vector (Figure 14).

Discussion
The search for creating unique vaccines against emerging diseases requires significant theoretical and technical expertise (Poland et al., 2018).The first generation of vaccines mainly relied upon live attenuated or inactive pathogens that require high-level biosafety regulations, stringent cold chains and packaging lines (Poland et al., 2018).Newgeneration vaccines, particularly synthetic and mRNA vaccines, are becoming unique approaches for controlling contagious diseases (Soleymani et al., 2021).Developing a vaccine with in silico can help reduce time and cost (Likitlersuang et al., 2018).Currently, there are no approved vaccines that protect against MPXV in the general population (Keckler et al., 2020).This study intended to identify new immunogenic epitopes by T and B cells to find the most optimal peptides and multiepitope proteins from the MPXV proteome for developing vaccine candidates (Van Regenmortel, 2016).To confirm peptides binding with high affinity, docking analysis was carried out to validate that the target peptides might bind to a particular area.In the following step, by using different methods, 62 linear B cell epitopes were identified (Table S3).At the  end of the process, only the peptides that could be detected by both T and B lymphocytes were kept because they can generate both humoral and cellular immunity.
To avoid autoimmune reactions, the anticipated B and T cell epitopic areas were checked for sequence similarity to the human proteome (Alais et al., 2015;Liljeroos et al., 2015).allergenic response and with high solubility.Together with the peptide vaccine candidates, two multi-epitope protein candidates against MPV were also made by incorporating the immunodominant regions of the A35R, A29L, and A30L proteins.
Based upon the evaluation of the predicted epitopes, five epitope-rich domains, consisting of high-scoring epitopes and those shared between T-and B-cell epitopes were chosen and linked with various linkers.The most ideal 3D structures were predicted by (GGGGS) 3 linker.The ALALAR protein consisted of three domains from A29L, A30L, and A35R, while the ALAL protein was made from A29L and A30L.The final ALALAR protein had 283 aa and a molecular weight of 29957.38Da.The ALAL protein had 152 aa with a molecular weight of 16731.80Da.Protein tertiary structure analysis was performed using MolProbity, ProSA-web and Ramachandran servers.The outcomes of these analyses showed that our designed candidates had an appropriate structure and were close to the structure of natural proteins.A more detailed analyses indicated that the developed proteins were antigens and not allergens.To assess the performance of the vaccine candidate in generating an immune response, ALALAR and ALAL protein interactions with TLR4 and M5 and M6 selected peptides with HLA-A � 1101, HLA-A � 01:01, HLA-A � 02:01, HLA-A � 03:01, HLA-A � 07:02, HLA-A � 15:01 and HLA-A � 30:01 immune receptors were carried out based on molecular docking.Examination of the binding results revealed that the designed proteins efficiently bind to TLR4, and that the selected peptides bind to the HLA-A � 1101, HLA-A � 01:01, HLA-A � 02:01, HLA-A � 03:01, HLA-A � 07:02, HLA-A � 15:01 and HLA-A � 30:01 receptors.The interaction energy of peptides M5 and M6 with these receptors   The results showed that the fluctuations for the HLA-A � 1101-M5 complex after about 40 ns were less than .3nm, which indicates the stability of the complex after this simulation time.The HLA-A � 1101-M6 complex had strong fluctuations after about 40 ns and reached relative stability after about 60 ns, but because the results of MMPBSA of this complex was positive, it indicated that the HLA-A � 1101-M6 complex is unstable in general.The ALAL-TLR4 complex reached stability after about 60 ns.Considering that the fluctuations of this complex after this simulation period were less than .3nm, it can be concluded that the ALAL-TLR4 complex is stable (Figure 6).The ALALAR-TLR4 complex reached stability after about 25 ns, with fluctuations less than .2nm.Therefore, we concluded that it is the most stable complex among the four peptide and protein complexes.MD simulation, the bond between the HLA-A � 1101 and TLR4 receptors and the ligands was relatively stable, indicating the complexes' binding stability.In total, the highest number of hydrogen bonds, and as a result the stronger bond, is related to the ALALAR-TLR4 complex.Additionally, the Rg specification was computed to examine the folding and compressibility of the complexes during the molecular dynamics simulation.The values of Rg in all three complexes also showed that the structure of the complexes was maintained throughout the simulation and that the receptor and ligand were not separated.The outcomes of the free energy prediction utilizing MMPBSA revealed that the M5 peptide had a stronger bond with the HLA-A � 1101 receptor than the M6 peptide.Furthermore, the C-ImmSim server showed that MPXV could simulate a strong immune reaction.
To validate our work, the IEDB-AR web server was used to confirm the results.The results showed that Liusong Yin et al. ( 2013) recognized YSIYENYGNI as a part of A28L that can be recognized by HLA-DR1 to detect vaccine viruses (Wallace et al., 1995).In addition, Kennedy et al. (2010) revealed that EQTSVFSATVYGDKIQ, which contains the SATVYGDKIQ sequence, and IRISMVISLLSMITMS, which contains ISLLSMITMS in the A33R of MPV might elicit IFN-c responses (Abraham et al., 2014).KQRLTNLEKK is one more peptide that can induce CTL actions to vaccinia virus in infected-HLA transgenic mice and creates an immune response against the vaccinia and variola virus (Enayatkhani et al., 2021).Furthermore, it has been shown that the ETLKQRLTNLEKKIT peptide, which includes the KQRLTNLEKK sequence as a part of the A27 protein, can bind to HLA-II and induces a protective immune reaction (Mobini et al., 2020;Sahoo et al., 2022).In this research, the predicted M5 peptide, ALAL, and ALALAR protein candidates revealed an appropriate lead in generating an immune response by immunoinformatic analysis.

Conclusion
In conclusion, the five MPXV proteins, E8L, A30L, A35R, A29L, and B21R were analyzed to find the best immunogenic peptide to induce a robust immune response.Docking and molecular dynamics results indicated that the M5 peptide (ISLLSMITMS) and ALAL and ALALAR proteins show suitable immune response against the MPX virus and that this response is stable during molecular dynamics simulation.In addition,

Figure 1 .
Figure 1.Strategies employed in the overall study.

Figure 3 .
Figure 3. (A) The tertiary structure of the ALALAR multi-epitope protein as a vaccine candidate against the MPX virus.(B) The tertiary structure of the ALAL multiepitope protein as a vaccine candidate against the MPX virus (C) The Schematic diagram of epitope-rich domains along with linkers for the ALALAR protein vaccine construct.(D) The Schematic diagram of epitope-rich domains along with linkers for the ALAL protein vaccine construct.

Figure 4 .
Figure 4. Secondary and tertiary structures of protein vaccine candidates were evaluated.Results confirmed the structures to be reliable and accurate.(A) Structure validation of ALALAR protein using Ramachandran plot.(B) Structure validation of ALAL protein using Ramachandran plot (C) Secondary structure of ALALAR protein constructed using GOR tool.(D) Secondary structure of ALAL protein constructed using GOR tool.(E) Structure validation of ALALAR protein using ProSA-web.(F) Structure validation of ALAL protein using ProSA-web.

Figure 7 .
Figure 7. (A) Graphical illustration of the ALALAR protein and the TLR4 complex.The residues of the protein are presented in red, purple, and yellow, with the TLR4 receptor in green.(B) LIGPLOT representation of the interaction between protein vaccine candidate and TLR4.Hydrogen bonds between TLR4 (blue) and the protein vaccine candidate (green) and hydrophobic interactions with TLR4 (black) and the protein vaccine candidate (blue) are indicated by dark green lines.

Figure 8 .
Figure 8. (A) Graphical illustration of the ALAL protein and the TLR4 complex.The residues of the protein are presented in red, purple, and yellow, with the TLR4 receptor in green.(B) LIGPLOT representation of the interaction between protein vaccine candidate and TLR4.Hydrogen bonds between TLR4 (blue) and the protein vaccine candidate (green) and hydrophobic interactions with TLR4 (black) and the protein vaccine candidate (blue) are indicated by dark green lines.

Figure 9 .
Figure 9.The RMSD values of the HLA-A � 1101 and TLR4 unbonded receptors, free ALALAR and ALAL proteins during the 150 ns of MD simulation.

Figure 11 .
Figure 11.The RMSF values of each receptor in the simulated complexes compared to the simulated unbonded forms of the receptors throughout the 150 ns simulation.(A) Comparison of the RMSF values of the HLA-A � 1101-M5 and HLA-A � 1101-M6 complexes with the unbonded HLA-A � 1101form.(B) Comparison of the RMSF values of the ALALAR-TLR4 and ALAL-TLR4 complexes with the unbonded TLR4 form.
Furthermore, RMSF analysis to examine the backbone atoms of unbonded receptors and complexes showed that fluctuations for M5-HLA-A � 1101 and M6-HLA-A � 1101 peptide complexes and designed protein-TLR4 complexes were more stable compared to the unbonded receptors.H-bond formation in proteins plays an essential role in their stability.The H-bond diagram of the complexes showed that during the

Figure 12 .
Figure 12.C-ImmSim is an in silico immune simulation against ALALAR protein as antigen.Simulations after the next three injections at steps 1, 84 and 168 are presented.Each time step is equal to 8 h.(a) Antigen and immunoglobins.(b) Cytokine production.(c) TH cell population.(d) B cell population.

Figure 13 .
Figure 13.C-ImmSim is an in silico immune simulation against ALAL protein antigen.Simulations after the next three injections at steps 1, 84 and 168 are presented.Each time step is equal to 8 h.(a) Antigen and immunoglobins.(b) Cytokine production.(c) TH cell population.(d) B cell population.

Figure 14 .
Figure 14.(A) Clone of the ALALAR protein gene in the pET-28a vector (highlighted in orange).(a) Informatics validation of the ALALAR protein clone by double digestion.(B) Clone of the ALAL protein gene in the pET-28a vector (highlighted in blue).(a) Informatics validation of the ALAL protein clone by double digestion.

Table 2 .
Overlapping epitopes of E8L, A30L, A35R and B21R proteins recognized by T and B lymphocytes with their immunogenic, allergenicity, solubility and toxicity properties to design multi-epitope vaccine construct.

Table 4 .
The epitope-rich regions of A29L, A30L, and A35R proteins selected for the final construct of vaccine candidate against MPX virus.

Table 3 .
Results of molecular docking evaluations between M5 and M6 peptides predicted as vaccine candidates and immune receptors.

Table 5 .
The features of ALALAR and ALAL multi-epitope proteins 3D structure designed as vaccine candidates against MPX virus.

Table 8 .
Evaluation of molecular docking results between protein vaccine candidates and TLR4.

Table 9 .
The van der Waals, electrostatic, polar solvation, SASA, SAV, and binding energies of the peptide-receptor complexes (kJ/mol), calculated by the MMPBSA method.