Abstract
Cervical cancer, induced by human papillomavirus (HPV), ranks as the fourth most prevalent malignancy among women globally. Unfortunately, existing prophylactic vaccines lack therapeutic efficacy. This study aimed to design a multi-epitope vaccine targeting the L1 and E7 proteins of HPV 16, 18, 33, and 45, with both preventive and therapeutic potential. Epitopes predicted using Immune Epitope Database (IEDB) and ABCpred were screened via immunoinformatics tools for antigenicity, immunogenicity, safety, conservancy, population coverage, and homology, and appropriate epitopes were assembled into a vaccine with suitable linkers and a 50-S L7/L12 adjuvant. The modeled and optimized vaccine was immunogenic, antigenic, safe, and displayed favorable physicochemical and solubility properties. Docking studies using ClusPro 2.0 and HDOCK indicated robust interactions between the vaccine and toll-like receptors TLR2/TLR4, and molecular dynamics simulations with Desmond validated the structural stability. Furthermore, molecular mechanics with generalized born and surface area solvation (MM/GBSA) analysis employing HawkDock showed favorable binding free energies of −82.86 and −76.72 kcal/mol, respectively. The vaccine’s potential efficacy was demonstrated by C-IMMSIM immune simulations, which revealed robust and long-lasting cellular and humoral responses, and also strong cytokine production. Finally, codon optimization for Escherichia coli K12 using JCat yielded a guanine-cytosine content of 50.69% and a Codon Adaptation Index of 0.97, and in silico cloning into pET28a(+) using SnapGene confirmed high expression potential. Our results indicate that the designed vaccine is a viable candidate for both preventive and therapeutic measures against high-risk HPV, requiring additional laboratory and animal studies.
Keywords
Introduction
Cervical cancer (CC) has been recognized as the fourth most prevalent malignancy among women globally, with approximately 0.6 million new cases and 0.3 million deaths annually. 1 In 2024, it was anticipated that there would be 13820 new cases and 4360 associated deaths in the United States; meanwhile, the European Union was expected to report 58 169 cases (with 56% originating from Central and Eastern Europe) and 22 989 related deaths. 2 Several factors contribute to the risk of CC, including early sexual activity, multiple sexual partners, limited access to screening, and a history of human papillomavirus (HPV) infection, cervical dysplasia, or abnormal Papanicolaou test results. 3 The HPV, mainly spread through sexual contact, stands as the leading cause of CC and is associated with more than 95% of cancers related to HPV, accounting for around 5% of all cancers globally.4-6 While the majority of infections resolve without clinical manifestations, chronic HPV infection may lead to cervical, other anogenital, and oropharyngeal malignancies, as well as anogenital warts. 7 Human papillomavirus infection induces papillomatous, hyperplastic, and wart-like lesions involving squamous cells on skin and several mucosal surfaces of humans, utilizing the host’s cellular machinery for replication. 8 Human papillomavirus genotypes are classified based on their potential to cause cancer: high-risk types (16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, 59, 68) are linked with epithelial and genital cancers, whereas low-risk types (6, 8, 11, 40, 42, 43, 44, 53, 54, 61, 72, 73) are linked to benign lesions and genital warts.9,10 According to a recent systematic review, of the more than 200 HPV genotypes, HPV 16, 18, 45, and 33 are most strongly associated with invasive CC, showing worldwide population attributable fractions of 61.7%, 15.3%, 4.8%, and 3.8%, respectively. 11 Region-specific studies confirm this global pattern, with HPV 16 consistently predominant (52% in Asia to 58% in Europe), followed by HPV 18 (13% in South/Central America to 22% in North America). HPV 33 and HPV 45, though less common, are detected across all regions at 3% to 7%. These data collectively justify targeting HPV 16, 18, 33, and 45 in CC prevention strategies worldwide. 12 Furthermore, HPV 16 and 18 are considered the primary causative factors for various malignancies affecting epithelial cells. 13
The HPV genome is composed of circular double-stranded DNA (dsDNA), roughly 8 kb in size, and lacks an envelope. 14 The HPV genome is divided into three regions: an early region encoding proteins E1 to E7, a late region encoding structural proteins L1 and L2, and a long control region (LCR) containing cis-elements crucial for viral replication and transcription. 15 In high-risk HPVs, E7 oncoprotein is one of the key drivers of HPV-induced carcinogenesis, enhancing epithelial cell transformation, proliferation, and immortalization.16,17 The E7 oncoproteins contribute to malignant transformation by targeting tumor suppressor proteins, particularly retinoblastoma protein (pRB). This interaction releases E2 F transcription factors, thereby promoting uncontrolled cell proliferation.18,19 Consequently, this oncoprotein is contemplated as a promising target for therapeutic HPV vaccine development.20,21 The major capsid protein, L1 (55 kDa), is highly abundant, constituting 80% of the capsid proteins (with an L1/L2 ratio of 10:1). 22 It plays a key role in antigenicity, receptor interaction, and the induction of neutralizing antibodies.18,23
As of now, the HPV vaccine landscape offers six authorized alternatives, all of which are classified as virus-like particle (VLP) vaccines. These include Cervavac®, GARDASIL9®, Cervarix®, GARDASIL®, WalrinvaxV, and Cecolin®. 24 These vaccines are prophylactic in nature, specifically designed to prevent HPV infection; however, they do not have therapeutic effects. Although existing vaccines have decreased the incidence of HPV-positive cases, it is crucial to enhance their efficacy and address adverse effects, including the potential for stronger immune responses compared to natural infection. 23 Recent advancements have focused on therapeutic HPV vaccines that stimulate cellular immunity, particularly CD4+ and CD8+ T cells. Unlike prophylactic vaccines, which primarily induce antibody responses, therapeutic vaccines aim to eliminate infected cells. 25 Considering the high incidence and mortality of HPV-related CC, there is a critical need for a potent therapeutic vaccine or a dual-action vaccine with both therapeutic and prophylactic effects.
To address this gap, we developed a multi-epitope chimeric vaccine targeting HPV 16, 18, 33, and 45. Recent studies have demonstrated the immunogenic potential of the multi-epitope vaccination strategy. For instance, a computationally designed dengue vaccine, which includes conserved CD4+ T-cell epitopes from all four serotypes, elicited robust antibody responses in vivo in rabbits. In another study, a multi-epitope Leishmania donovani vaccine enhanced cytokine and nitric oxide production in vitro, proving the potential for immunoinformatics-driven vaccine development.26,27 We selected the E7 oncoprotein for its role in carcinogenesis and the L1 capsid protein for its ability to form highly immunogenic VLPs. Traditional vaccine development is time-consuming and resource-intensive. 28 On the other hand, immunoinformatics and reverse vaccinology offer a more efficient and cost-effective approach. Reverse vaccinology is a genome-based method that systematically screens pathogen proteins to identify promising candidates for vaccines. This strategy is particularly advantageous for pathogens that are difficult to culture or antigens that are challenging to express in vitro.28-30 The polymorphism observed in host genetics may influence the immunological reaction to an infectious agent within the target population. Considering the diverse characteristics of major histocompatibility complex (MHC) alleles, their ability to attach to specific repertoires of epitopes (short peptides) derived from processed pathogens is crucial. The delivery of these complexes to T cells plays a significant role in developing a robust cytotoxic T-cell (Tc) and helper T-cell (Th) response targeting the infectious agent. 31
The objective of this research was to construct a multi-epitope chimeric vaccine with the potential to elicit both prophylactic and therapeutic immune responses against the four high-risk human papillomaviruses (HPV 16, 18, 33, and 45) utilizing a reverse vaccinology approach. Initially, we identified the most effective Tc-cell, Th-cell, and linear B lymphocyte epitopes of the E7 oncoprotein and L1 major capsid protein that could elicit robust immune responses, including CD4+ Th-cell activation and CD8+ cytotoxic T (Tc) cell activation, interferon-alpha (IFN-alpha) production, and humoral responses (including activation of B-lymphocytes and production of antibody) by utilizing a range of algorithms. Prioritized epitopes assessed for antigenicity, allergenicity, toxicity, cytokine-inducing capability, and human homology were connected with suitable linkers and an immune-enhancer adjuvant for formulating a chimeric multi-epitope vaccine (MEV) targeting HPV infection and CC. Furthermore, an array of analyses was performed, including the examination of physicochemical properties, structural prediction, validation, molecular docking and dynamics simulations, free binding energy calculations, disulfide engineering, in silico cloning, immune simulations, and post-translational modifications, to assess the vaccine’s stability, functionality, and suitability. In summary, this study presents a scientifically rigorous in silico design of a multi-epitope chimeric vaccine with both prophylactic and therapeutic potential against high-risk HPV genotypes associated with CC. A graphical overview of our vaccine design and methodology is presented in Figure 1.
Schematic workflow of the multi-epitope vaccine (MEV) design strategy against high-risk human papillomavirus (HPV) types 16, 18, 33, and 45. The process involved sequential steps beginning with L1 and E7 protein sequence retrieval, followed by epitope prediction, filtering, and MHC clustering. Prioritized CTL epitopes were docked with MHC alleles and used for multi-epitope vaccine construction. Physicochemical and immunological properties, secondary and tertiary structure modeling, conformational B-cell epitope prediction, and receptor (TLR2/4) docking were performed. Structural stability was assessed using molecular dynamics (MD) simulation and molecular mechanics with generalized born and surface area solvation (MM/GBSA) free energy analysis, along with disulfide engineering for stability enhancement. Immune response simulation, codon optimization, in silico cloning, and post-translational modification analysis led to the final MEV candidate.
Methods
Protein sequence retrieval
Sequence retrieval is an important first step in computational vaccine design, acting as the basis for identifying, analyzing, and selecting potential vaccine candidates. These sequences are required for identifying potential antigens that may trigger an immunological response. In this study, we retrieved the amino acid sequences of the major capsid protein L1 and the E7 oncoprotein of HPV 16, 18, 33, and 45 from the UniProt database in FASTA format (https://www.uniprot.org/, accessed on February 6, 2025). 32 The specific UniProt IDs and corresponding protein lengths for each genotype are as follows: for HPV 16, L1 (P03101, 505 AA) and E7 (P03129, 98 AA); for HPV 18, L1 (P06794, 568 AA) and E7 (P06788, 105 AA); for HPV 33, L1 (P06416, 499 AA) and E7 (P06429, 97 AA); and for HPV 45, L1 (P36741, 539 AA) and E7 (P21736, 106 AA).
Prediction of cytotoxic T cell epitopes
Peptides displayed on MHC I molecules at the host cell surface are recognized by the CD8+ T cell population to identify oddities, particularly pathogens. The effectiveness of the peptide-MHC (p-MHC) complex interaction with T-cell receptors (TCRs) is influenced by both the MHC I molecule and the bound peptide. Tc cells rapidly eliminate infected cells and viruses from the host organism. 33 Immune Epitope Database’s (IEDB) MHC-I Binding Predictions v2.24 server (http://tools.iedb.org/mhci/, accessed on February 7, 2025) was employed to forecast the Tc-cell epitopes of the E7 oncoprotein and L1 major capsid protein of HPV 16, 18, 33, and 45. 34 The method for prediction was configured according to the IEDB recommended 2023.09 (NetMHCpan 4.1 EL), with the MHC originating species designated as human, the MHC allele(s) aligned with the human leukocyte antigen (HLA) allele reference set, and the length of the epitope set to 9. 35 Predicted Tc-cell epitopes were retrieved based on their percentile rank. Those with values below 0.5 were prioritized for further evaluation, as lower percentile ranks indicate stronger binding affinity. The VaxiJen v2.0 (https://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html, accessed on February 11, 2025) server was utilized to evaluate the antigenic properties of the Tc-cell epitopes, establishing a threshold score of 0.4 and choosing the virus as the target. 36 VaxiJen v2.0 has a misclassification rate of 11% to 30%. Furthermore, the immunogenicity of these epitopes was evaluated through the MHC Class I Immunogenicity server (http://tools.iedb.org/immunogenicity/, accessed on February 14, 2025), retaining those with positive scores. The VaxiJen v3.0 server (https://www.ddg-pharmfac.net/vaxijen3/home/, accessed on February 16, 2025) was also utilized, selecting epitopes predicted to be “probable immunogen.” 37 The Tc-cell epitope allergenicity was forecasted utilizing two distinct servers: AllerTOP v2.1 (https://www.ddg-pharmfac.net/allertop_test/, accessed on February 18, 2025) 38 and AllergenOnline (http://www.allergenonline.org/databasefasta.shtml, accessed on August 30, 2025). 39 For AllergenOnline, an E-value cutoff of 1 and a maximum of 20 alignments were applied as default parameters, and sequence matches showing greater than 35% identity over an 80-amino acid sliding window were considered potential allergens. The toxicity was predicted through ToxinPred (http://crdd.osdd.net/raghava/toxinpred/, accessed on February 20, 2025) employing the default parameters and the support vector machine (SVM) method. 40 The predicted error rates were approximately 14.7% for AllerTOP v2.1 and 5.5% for ToxinPred, respectively. Efficient CTL epitopes are produced through proteasomal degradation in infected cells and subsequently translocated into the endoplasmic reticulum via a transporter associated with antigen processing (TAP) for presentation on MHC I molecules. 41 Hence, the selected Tc-cell epitopes were further assessed for proteasomal cleavage and TAP transport efficiency via the IEDB MHC Class I Processing tool (http://tools.iedb.org/processing/, accessed on August 30, 2025). 42 The analysis was performed using the IEDB recommended prediction method (version 2013-02-22), specifying human MHC and selecting the immuno-type proteasomal cleavage prediction. Transporter associated with antigen processing transport predictions were made with a maximum precursor extension of 1 and an alpha factor of 0.2 to determine the potential for natural presentation.
Prediction of helper T-cell epitopes
The Th cells are vital to adaptive immunity, sensing exogenous antigens and activating Tc cells and B cells to eliminate them, which is crucial for developing an effective vaccine. 43 IEDB’s MHC-II Binding Predictions server (http://tools.iedb.org/mhcii/, accessed on February 25, 2025) was employed to forecast the Th-cell epitopes (15-mer) from selected L1 and E7 proteins of HPV. The IEDB prediction method was configured to utilize the IEDB-recommended version 2023.09 (NetMHCIIPan 4.1 EL), with humans chosen as the originating species for MHC and the seven-allele HLA reference set designated as the MHC allele(s). A consensus percentile rank of <0.5 was implemented to filter the Th-cell epitopes to identify high-affinity epitopes associated with MHC II alleles. The evaluation of Th-cell epitopes’ antigenicity was performed using the VaxiJen v2.0 server. Toxicity assessments were carried out with ToxinPred, while allergenicity evaluation was conducted using the AllerTOP v2.1 and AllergenOnline servers. The Th cells secrete cytokines, including interleukin 4 (IL-4) and interferon-gamma (IFN-γ), which can trigger the immune cells within the body. Moreover, cytokines can persist past inflammatory responses and prevent the damage of tissue. Hence, the IFNepitope server (webs.iiitd.edu.in/raghava/ifnepitope/index.php, accessed on March 5, 2025) was employed to predict the interferon-gamma (IFN-γ) inducibility of Th-cell epitopes, configuring the model as IFN-γ versus non-IFN-γ and utilizing the Motif and SVM hybrid method. 44 The hybrid model demonstrated performance metrics of 86.97% for sensitivity, 68.38% for specificity, and 80.93% for accuracy. Only a positive IFN-γ score displaying Th-cell epitopes was retained. Finally, the peptides were assessed to determine their capacity to induce IL-4 through the utilization of the IL4pred server (http://crdd.osdd.net/raghava/il4pred/, accessed on March 5, 2025), with an estimated error rate of ~24.24%. 45
Linear B-cell epitope prediction
B-cell epitopes are crucial for eliciting antibody-mediated (humoral) immunity. They are composed of clusters of amino acid residues that interact with secreted antibodies, activating the body’s defense system to combat intruders. 46 In this study, the ABCpred server (http://crdd.osdd.net/raghava/abcpred/, accessed on March 7, 2025) was employed for predicting, evaluating, and screening immune-enhancing linear B-cell epitopes, owing to its notable accuracy. 47 The screening threshold was maintained at 0.51, and the epitope length was restricted to 16 amino acids. The prediction performance of ABCpred was 65.93%, with an anticipated misclassification of ~34.07%. Subsequently, VaxiJen v2.0, VaxiJen v3.0, AllerTOP v2.1, AllergenOnline, and ToxinPred servers were used to test the antigenicity, immunogenicity, allergenicity, and toxicity of these candidate B-cell epitopes.
Human homology analysis of selected epitopes
The epitopes were analyzed for homology to eliminate any potential for autoimmune or lack of response. The BLASTP service, an online sequence alignment tool provided by NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE=Proteins, accessed on March 13, 2025), was utilized to ascertain the human homology of the chosen epitopes. The epitopes were supplied separately in FASTA format for “Homo sapiens (taxid:9606)” as the organism, while all other parameters remained at their default settings. Epitopes were deemed non-homologous if their E-value exceeded 0.05.48,49
Estimation of conservancy and population coverage
Selecting epitopes from conserved regions is crucial for developing effective epitope-based vaccines, ensuring optimal immunity against diverse pathogen strains. The IEDB Conservancy Analysis tool (http://tools.iedb.org/conservancy/, accessed on March 14, 2025) was utilized to evaluate the conservation levels of the screened T-cell and B-cell epitopes obtained from the L1 and E7 protein sequences across four high-risk HPV genotypes (HPV 16, 18, 33, and 45). Moreover, vaccination initiatives must focus on achieving coverage across different communities and regions while considering the significant differences in HLA profiles to safeguard a diverse population. In computational vaccine design, this involves assessing the distribution of HLA alleles linked to predicted epitopes. Consequently, the IEDB Population Coverage Analysis tool (http://tools.iedb.org/population/, accessed on March 16, 2025) was employed to evaluate the overall population coverage of potential MHC I and MHC II epitopes across various countries and ethnic groups. This comprehensive methodology offers an integrated approach to evaluating and addressing the factors that impact vaccine efficacy.
Tc-cell epitope modeling and molecular docking studies
The interaction of the selected Tc-cell epitopes with the HLA binding allele was analyzed by molecular docking to assess their binding affinity and examine structural compatibility, as well as the immunologic response generated, which is crucial for vaccine development endeavors. The three-dimensional structures of the MHC I epitopes were predicted utilizing the sOPEP sorting scheme, which involved running 200 simulations on the PEP-FOLD v3.5 online platform (https://bioserv.rpbs.univ-paris-diderot.fr/services/PEP-FOLD3/, accessed on March 17, 2025). 50 This server utilizes a novel method to predict the conformations of small peptides, ranging from 5 to 50 amino acid residues. It employs the Taboo/Backtrack sampling algorithm to forecast five probable structures for each peptide sequence. Evaluating PEP-FOLD3 with 56 peptides in aqueous solution resulted in experimental-like conformations for 80% of the targets. The human alleles HLA-B*44:03 (PDB code: 4JQX), HLA-A*24:02 (PDB code: 7JYV), HLA-A*68:01 (PDB code: 6PBH), HLA-B*35:01 (PDB code: 8EMJ), and HLA-A*02:03 (PDB code: 3OX8) were chosen for selected MHC I epitopes, with the co-crystallized peptides serving as positive controls to validate our docking methodology. The Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB PDB) server (https://www.rcsb.org/, accessed on March 18, 2025) was utilized for retrieving the three-dimensional crystal structures of HLA alleles, including co-crystallized peptides. 51 Initially, Discovery Studio 2024 was utilized in the process of protein preparation by detaching the ligands from the structures and incorporating polar hydrogens, followed by energy minimization utilizing SWISS PDB Viewer (SPDV) v4.1. (http://www.expasy.org/spdbv/, accessed on March 18, 2025).52,53 Molecular docking simulations were subsequently performed utilizing Cluspro 2.0 (https://cluspro.org/home.php, accessed on March 19, 2025) and HDOCK (http://hdock.phys.hust.edu.cn/, accessed on March 19, 2025) protein-protein docking servers.54,55 ClusPro 2.0 employs a rigid-body docking algorithm that produces a wide array of potential docking poses (decoys) and subsequently organizes them through clustering based on energy minimization, thereby pinpointing the optimal binding orientation patterns. On the other hand, HDOCK operates by sampling and calculating the atomic shape representation; coordinating surface adjustments while distinguishing the surface of the potential binding modes between the two proteins through a Fast Fourier Transform-based global search method and evaluating the sampled binding modes with an improved iterative template-based scoring function for protein–protein interaction. The molecular docking outcomes were visualized utilizing Discovery Studio 2024 software, with figures generated through UCSF ChimeraX. 56 Finally, we employed the PRODIGY server (https://rascar.science.uu.nl/prodigy/, accessed on March 19, 2025) to determine the binding strength of these epitope-MHC I allele complexes. 57
Analysis of MHC clusters
The MHC genomic region exhibits significant polymorphism across a wide range of species. The human MHC genomic region (HLA), characterized by its high diversity, contains thousands of alleles. In most instances, the potentially unique specificity of MHC alleles remains uncharacterized. 58 Cluster analysis of MHC alleles can reveal both MHC I and MHC II molecules that share comparable binding specificities. Therefore, we used the MHCcluster-2.0 online server (http://www.cbs.dtu.dk/services/MHCcluster/, accessed on March 20, 2025) to create phylogenetic trees and highly intuitive heatmaps of the functional clusters among MHC variants. For both MHC class I and II analyses, 50000 peptides were included, with 100 bootstrap calculations and a 10% fraction of peptides used for correlation analysis. For MHC class I, the NetMHCpan-2.8 prediction method was selected, and the allele set HLA Prevalent and Characterized was used. For MHC class II, the allele set HLA-DR representatives was applied.58,59
Construction of a multi-epitope vaccine
In constructing the MEV, we strategically incorporated high-priority, immunodominant Tc-cell, Th-cell, and B-cell epitopes with immunological adjuvant and specified linker sequences. To enhance immunogenicity, AAY, GPGPG, KK, and EAAAK linkers were strategically utilized. Specifically, AAY linkers were deployed for joining the Tc-cell epitopes, GPGPG linkers were utilized to link the Th-cell epitopes, and KK linkers combined B-cell epitopes. The use of these linkers was crucial for achieving effective separation of individual epitopes while also preventing the formation of junctional or neoepitopes, thereby enhancing the processing and presentation of the epitopes. The junctional epitopes significantly impede the effectiveness of the chosen epitopes. 60 Furthermore, the 50S ribosomal protein L7/L12 (UniProt accession: P0A7 K2), a recognized TLR4 agonist, was selected as an adjuvant based on its established use in prior HPV vaccine designs and seamlessly fused to the amino-terminal end of the MEV sequence using an EAAAK linker to improve its immunogenic potential.1,61 This structured arrangement ensures optimal epitope presentation, minimizes junctional immunogenicity and enhances vaccine stability. Furthermore, intrinsic disorder and unstructured regions within the vaccine protein were predicted using DisEMBL 1.5 (http://dis.embl.de/, accessed on September 14, 2025) and IUPred3 (https://iupred3.elte.hu/, accessed on September 14, 2025), both widely used computational tools for disorder analysis in proteins.62,63 This analysis was performed to evaluate the flexibility of linker regions and their potential impact on epitope accessibility and structural integrity of the multi-epitope construct.
Evaluation of immunological, physicochemical, and solubility characteristics
The MEV was thoroughly evaluated regarding its allergenicity, antigenicity, immunogenicity, toxicity, and physicochemical properties. The prediction of antigenic properties, a crucial element in vaccine development, was conducted utilizing the VaxiJen v2.0 server with a threshold value established at 0.4. Furthermore, the assessment of immunogenicity, allergenicity, and toxicity was conducted utilizing the VaxiJen v3.0, AllerTop v.2.1, and ToxinPred servers, respectively. In addition, potential cross-reactivity with known allergens was evaluated utilizing the AllergenOnline database. The Expasy ProtParam online server (https://web.expasy.org/protparam/, accessed on March 21, 2025) was employed for predicting several physicochemical features of the MEV, encompassing the composition of amino acids, theoretical pI, molecular weight, half-life, instability index, the grand average of hydropathy (GRAVY), and aliphatic index. 64 The isoelectric point refers to the pH at which an amino acid or protein becomes neutral, resulting in the loss of its charge and its inability to migrate in a direct current electric field. Determining the theoretical pI of a protein can be significantly helpful in the selection and optimization of protein purification methods, such as isoelectric focusing electrophoresis and ion-exchange chromatography. 65 Half-life is the predicted time it takes for half of the protein molecules to degrade after synthesis in a given expression system. 66 The instability index of proteins reflects their stability in vitro, with proteins exhibiting an instability index below 40 classified as stable. 67 The aliphatic index of a protein denotes the proportion of the protein’s volume attributed to aliphatic side chains, serving as an indicator of the protein’s thermal stability. 68 The GRAVY is calculated by dividing the combined hydropathic values of all amino acids by the total number of residues in the sequence. The GRAVY illustrates the proteins’ amphipathic characteristics, with negative and positive values indicating the hydrophilic and hydrophobic properties of the designed multi-epitope vaccine, respectively. 69 To forecast the solubility of the MEV construct, the SOLpro server (http://scratch.proteomics.ics.uci.edu/, accessed on March 21, 2025) was utilized, where a threshold probability score of ⩾ 0.5 indicates solubility and scores <0.5 indicate insoluble characteristics of the protein. 70 The server demonstrated an overall prediction accuracy exceeding 74% based on multiple 10-fold cross-validation runs.
Secondary structure prediction of the multi-epitope vaccine construct
The secondary structural elements of the designed MEV construct, such as alpha-helical regions, extended strands, and random coils, were predicted by the SOPMA (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=/NPSA/npsa_sopma.html, accessed on March 22, 2025) 71 and PSIPRED 4.0 (http://bioinf.cs.ucl.ac.uk/psipred/, accessed on March 22, 2025) online platforms. 72 Each of the parameters was set to its default value. PSIPRED is an effective and straightforward method for predicting 2D structures. It utilizes two feed-forward neural networks to analyze results generated from PSI-BLAST (Position-Specific Iterated-BLAST). These two servers execute the prediction with an accuracy exceeding 80%. The structural features (two-dimensional) were obtained and assessed to comprehend the structural quality of the designed multi-epitope vaccine.
Tertiary structural modeling, refinement, and validation
We sought to predict and model the three-dimensional structure of the MEV construct to facilitate docking studies with human immune receptors, following an evaluation of its physicochemical and immunological properties. In pursuit of this aim, the vaccine construct was modeled utilizing the Robetta (robetta.bakerlab.org, accessed on March 23, 2025) server. 73 Furthermore, the GalaxyRefine module available on the GalaxyWEB server (https://galaxy.seoklab.org/, accessed on March 24, 2025) was utilized to enhance the quality of the crude 3D model of the MEV construct. 74 Five refined models were produced in the results, integrating a range of parameters such as root mean square deviation (RMSD), poor rotamers, clash score, MolProbity, GDT-HA, and Rama-favored. Based on the criterion that a superior model is indicated by a lower MolProbity value, we selected the optimal best-refined 3D model and submitted it to the PROCHECK tool within the SAVES v6.1 structure validation service (https://saves.mbi.ucla.edu/, accessed on March 25, 2025) 75 and the ProSA-Web server (https://prosa.services.came.sbg.ac.at/prosa.php, accessed on March 25, 2025) 76 for quality assessments. The refined models were evaluated, and the most suitable one was selected for further analysis.
Prediction of conformational B-lymphocyte epitopes
The linear and conformational B-lymphocyte epitopes in the refined MEV construct were anticipated utilizing the Ellipro tool from the IEDB online server (http://tools.iedb.org/ellipro/, accessed on March 26, 2025). 77 This method forecasts conformational B-lymphocyte epitopes by analyzing the protein structure’s geometrical characteristics, flexibility, and solvent accessibility, resulting in epitopes that are longer than others. The default value for the minimum score was established at 0.5, whereas the maximum distance was predetermined to be 6 Å. ElliPro serves as a trustworthy tool for recognizing antibody epitopes and forecasting B-cell conformational epitopes, achieving an AUC score of 0.732. 77
Molecular docking
Evaluating and forecasting the effectiveness of the vaccine construct is a crucial step in reverse vaccinology, as it is vital for generating an appropriate immune response. The molecular docking process elucidates the preferred binding orientations and patterns of interactions of vaccines with immunological receptors when in stable complex form. In this study, we utilized the ClusPro 2.0 (https://cluspro.org/help.php, accessed on March 27, 2025) and HDOCK (http://hdock.phys.hust.edu.cn/, accessed on March 27, 2025) protein-protein docking servers to anticipate the binding affinity of the MEV construct with two immunological receptors, TLR2 (PDB entry code: 2Z7X) and TLR4 (PDB entry code: 3FXI).54,55 Notably, both ClusPro 2.0 and HDOCK rely on rigid-body docking, which may limit accuracy for proteins undergoing substantial conformational changes. In addition, ClusPro 2.0’s results depended on the chosen scoring schemes, while HDOCK’s predictions were influenced by the quality of homology-modeled structures, particularly for low-identity or multi-chain proteins. Initially, the aforementioned TLR receptors were acquired from the RCSB PDB database in PDB format, and water molecules, bound ligands, and additional chains were eliminated using the Discovery Studio 2024 software package. 52 The best-docked complexes exhibiting the most robust interactions among residues and the minimal interaction energy were selected, as a lower energy score suggests a higher binding affinity, and these complexes were utilized as the starting point for molecular dynamics simulation. 78 Finally, all the molecular interactions within the interacting plane between the MEV construct and TLRs were visualized and scrutinized utilizing UCSF Chimera 1.19 software 79 and the PDBsum online server (https://www.ebi.ac.uk/thornton-srv/databases/pdbsum/, accessed on March 28, 2025). 80
Molecular dynamics simulation
Inspecting the intermolecular binding interaction patterns between the designed MEV construct and the prioritized toll-like receptors (TLR) requires a detailed analysis of their dynamic behavior. Therefore, MD simulations lasting 100 ns were performed using Desmond 2023.2 software by Schrödinger LLC on a Linux platform to scrutinize the structural integrity of the MEV-apo (unbound protein state), MEV-TLR2, and MEV-TLR4 complexes in a simulated physiological environment. 81 The docked complexes underwent preprocessing using the Protein Preparation Wizard of Maestro, which involved incorporating hydrogens, creating disulfide bonds, filling in missing side chains, and removing water. The protonation states were optimized utilizing Epik at pH 7.0 ± 2.0, followed by PROPKA for final pH-based adjustments at pH 7.0. 82 The System Builder tool was used to construct the system, with complexes positioned within a TIP3P predefined water solvent model, confined by orthorhombic periodic boundary conditions. The process of electrical neutralization involved the incorporation of sodium and chloride ions, with a 0.15 M NaCl solution being utilized to replicate physiological conditions. The system was minimized and equilibrated under the OPLS3e force field, including energy minimization and equilibration with positional restraints. 83 The simulations were conducted in the Isothermal-Isobaric (NPT) ensemble, using a Nose-Hoover thermostat to control temperature at 300 K, and Martyna-Tobias-Klein barostats to maintain a pressure of 1.01 bar. 84 The simulation interaction diagram (SID) from the Schrödinger package was utilized to interpret the results. Trajectory data was recorded at 100 picosecond intervals, resulting in a total of 1000 frames throughout the 100 ns simulation. Conformational fidelity of the MEV, MEV-TLR2, and MEV-TLR4 complexes was examined through RMSD and root mean square fluctuation (RMSF), offering a thorough insight into the dynamics of the vaccine-receptor complexes.
Molecular mechanics with generalized born and surface area solvation binding free energy calculation
The molecular mechanics with generalized born and surface area solvation (MM/GBSA) program on the HawkDock web server (http://cadd.zju.edu.cn/hawkdock/, accessed on April 3, 2025) was employed to estimate the free binding energy of interactions between the MEV and immune receptors TLR2 and TLR4. 85 This approach integrates molecular mechanics with the Generalized Born model to analyze key intermolecular interactions, including van der Waals forces (VDW), electrostatic interactions (ELE), and solvation effects—both polar (GB) and nonpolar (SA). Recognized for its reliability, the HawkDock server offers an accuracy of 80% to 95% in MM/GBSA calculations for both crystallized and predicted structures. 86
Disulfide engineering
Disulfide engineering is a cutting-edge technique involving mutating cysteine residues to establish disulfide bonds in a specific, highly flexible region of a protein. Disulfide linkages are strategically incorporated to enhance stability in the folded conformation of the protein. These linkages increase the free energy of the denatured state and decrease conformational entropy. 87 Consequently, we used the Disulfide by Design 2.0 v2.13 online platform (http://cptweb.cpt.wayne.edu/DbD2/, accessed on April 4, 2025) to locate potential residue pairs within the vaccine construct that could be altered to cysteine, thereby promoting the establishment of new disulfide bonds. 88 The χ3 value and the Cα-Cβ-Sγ angle were maintained at their default parameters. Based on earlier studies, for disulfide bridging, the χ3 angle should be maintained within the −87 to +97° range, while the energy score should not surpass 2.2 kcal/mol. The configuration of the disulfide bond at a higher angle deviates from the optimal condition, impacting its stability. Furthermore, a disulfide bond with significantly high bond energy is prone to hydrolysis or reduction, influencing its stability. Furthermore, epitope mapping and solvent-accessible surface area (SASA) calculations were performed using PyMOL version 3.1 to assess the immunogenic exposure of each B-cell epitope before and after disulfide bond introduction. 89
Immune simulation
Immunological simulations were performed employing the C-IMMSIM server (https://kraken.iac.rm.cnr.it/C-IMMSIM/index.php, accessed on April 5, 2025) to assess the vaccine’s immunogenicity and analyze its real-life immunogenicity profile. 90 This server functions as an agent-based simulator for immunological responses. It serves as a versatile online tool which employs position-specific scoring matrices (PSSM) techniques for predicting immune epitopes and their interactions with immune receptors. Simultaneously, it replicates three unique anatomical regions present in mammals, namely the tertiary lymphoid organs, thymus, and bone marrow. 91 The vaccine was delivered via 3 injections aimed at mimicking immunity, with each dose comprising 1000 molecules and a 4-week gap between each administration. The three injections were given at intervals of 1, 84, and 168, each roughly corresponding to an 8-hour duration in real life. 92 The random seed and simulation volume parameters were established at 12345 and 10 μl, respectively, with a total of 1095 simulation steps performed. The recommended dosage for the vaccine involves administering three injections at 4-week intervals, which aligns with the standard timing between doses for all commercial vaccines to inspect the impact of vaccine exposure.
Codon optimization and in silico cloning of the final vaccine construct
Producing recombinant vaccine proteins is a crucial step in developing MEV. To prepare a vaccine sequence for cloning and expression in a suitable vector, reverse translation and codon optimization were performed utilizing the Java Codon Adaptation Tool (JCat) (http://www.jcat.de/, accessed on April 7, 2025). 93 Consequently, the vaccine amino acid sequence was uploaded to JCat, with Escherichia coli (K12 strain) chosen as the host for the vaccine protein expression. The whole operation was performed while circumventing each of the 3 parameters: (1) restriction enzyme recognition sites, (2) binding sites for the prokaryotic ribosome, and (3) rho-independent termination of transcription. The adjusted sequence was assessed using the Codon Adaptation Index (CAI) value and guanine-cytosine (GC) content, both of which are crucial in determining the extent of protein expression. After optimization, the sequence was inserted into the pET28a (+) vector utilizing EcoRI and BamHI restriction sites at the 5′ and 3′ ends, employing SnapGene software for the procedure (https://www.snapgene.com/, accessed on April 11, 2025).
Forecasting post-translational modifications
The MusiteDeep platform facilitated the analysis of the post-translational modifications (PTM) of the MEV (https://www.musite.net/, accessed on April 8, 2025). 94 The emphasis was placed on phosphorylation (p), acetylation (ac), ubiquitination (ub), glycosylation (gl), hydroxylation (hy), and methylation (me). The selection of these processes was deliberate to facilitate optimal protein folding, thereby enhancing the vaccine’s stability, solubility, and biological effectiveness. 95 To improve robustness, further PTM evaluations were carried out utilizing specialized servers accessible at (http://www.cbs.dtu.dk/services/, accessed on August 31, 2025), specifically NetAcet-1.0 for acetylation, NetPhos-3.1 for phosphorylation, and NetNGlyc-1.0 for N-linked glycosylation.96-98
Result
Prediction of cytotoxic T-cell epitopes
Major histocompatibility complex I-binding cytotoxic T (Tc) cell epitope prediction was conducted for the E7 and L1 proteins of HPV genotypes 16, 18, 33, and 45. For the L1 protein, 98, 119, 111, and 112 epitopes were identified, respectively, while the E7 protein yielded 13, 11, 14, and 10 epitopes with a percentile rank below 0.5 (Supplementary Table S1). After a comprehensive evaluation of the shortlisted epitopes based on antigenicity, immunogenicity, toxicity, allergenicity, MHC I immunogenicity, and human homology, the most suitable Tc-cell epitopes were selected for MEV formulation. The immunological profile analysis of the L1 protein identified a total of 5, 8, 7, and 6 Tc-cell epitopes for HPV types 16, 18, 33, and 45, respectively. On the other hand, the analysis of the E7 protein showed 3 epitopes for HPV 16 and 1 epitope for each of HPV types 18, 33, and 45, respectively. These epitopes demonstrated strong immunogenic and antigenic properties while lacking allergenicity, toxicity, or similarity to the human proteome. Based on these criteria, four of the most promising Tc-cell epitopes from both proteins were considered optimal candidates for MEV incorporation (Supplementary Table S2). Processing analysis of the shortlisted Tc-cell epitopes was performed using proteasomal cleavage, TAP transport, and MHC binding predictions (Supplementary Table S3). Total scores, which integrate these three components, reflect the overall likelihood of a peptide being naturally processed, transported, and presented on MHC I molecules, with higher scores indicating greater likelihood. LTVGNPYFR (HLA-A*68:01), EEYDLQFIF (HLA-B*44:03), and QAQPATADY (HLA-B*35:01) exhibited the highest total scores, suggesting they are the most likely to be naturally presented. Some epitopes, such as HVEEYDLQF, showed allele-dependent variability. These immune-enhancing cytotoxic T cell epitopes of E7 and L1 proteins were prioritized for further conservancy analysis.
Prediction of helper T-cell epitopes
For the L1 protein of HPV genotypes 16, 18, 33, and 45, a total of 13, 14, 14, and 21 MHC II-binding Th-cell epitopes were identified, respectively (Supplementary Table S1). However, no predictions were found for the E7 protein. After assessing a range of parameters, including antigenic properties, potential allergenic reactions, toxicity assessments, human homology, and the ability to stimulate cytokine production such as interferon-γ (IFN-γ) and IL-4 in Th cells, a selection process was conducted to pinpoint the most appropriate Th-cell epitopes. Finally, two Th-cell epitopes were chosen: GRKFLLQAGLKAKPK (from HPV 16 and 33) and RDNVSVDYKQTQLCI (from HPV 18 and 45) for MEV incorporation (Supplementary Table S4). These high-priority immune-enhancing Th-cell epitopes were prioritized for subsequent analysis regarding their conservancy.
Evaluation of linear B lymphocyte epitopes
A total of 53, 58, 47, and 54 B-cell epitopes were identified in the L1 protein of HPV genotypes 16, 18, 33, and 45, respectively. The E7 protein yielded 8, 8, 10, and 9 B-cell epitopes for the same genotypes (Supplementary Table S1). After a comprehensive evaluation of the shortlisted epitopes based on immunogenicity, antigenicity, toxicity, allergenicity, and human homology, the most suitable B-cell epitopes were selected for MEV formulation. For L1 protein, a total of 8, 9, 9, and 7 epitopes were identified for HPV genotypes 16, 18, 33, and 45, respectively. These epitopes exhibited immunogenicity and antigenicity while being devoid of allergenicity, toxicity, and homology to the human proteome. Through this assessment, five ideal B-cell epitopes of the L1 protein were recognized for MEV design. In the case of the E7 protein, a total of 3, 2, 1, and 2 epitopes were identified for HPV types 16, 18, 33, and 45, respectively. One epitope (LQPETTDLYCYEQLND) was predicted as allergenic by AllerTOP v2.0 but classified as non-allergenic by the AllergenOnline server, while another epitope (TIQQLLMGTVNIVCPT) was predicted as non-immunogenic. Consequently, the four most promising B-cell epitopes of E7 were recognized as the optimal candidates for vaccine development (Supplementary Table S5). These immune-enhancing B-cell epitopes of E7 and L1 proteins were prioritized for further conservancy analysis.
Epitope conservancy and population coverage analysis
A total of 8 Tc-cell epitopes, 2 Th-cell epitopes, and 9 B-cell epitopes were identified with strong immunological characteristics after screening high-priority immunodominant epitopes. These epitopes were then evaluated for population coverage and conservancy across HPV genotypes 16, 18, 33, and 45 and incorporated into the MEV. Subsequently, these epitopes were sequentially labeled from Ep1 to Ep19 to represent conservancy values graphically (Figure 2). The epitopes predicted from the L1 protein showed remarkable conservation across all four genotypes, underscoring their potential to induce cross-protective immunity. However, the oncogenic E7 protein epitopes displayed considerable variability in their amino acid sequences (Figure 2). The epitopes Ep1 and Ep2 showed complete conservation across the three genotypes (HPV 16, 18, and 45) and approximately 90% in HPV 33. The epitopes Ep3, Ep4, Ep9, Ep10, and Ep11 exhibited 100% conservation in at least two of the four HPV genotypes, while their conservation in the remaining genotypes ranged from 60% to 82%. Overall, their conservation across all four genotypes ranged from 60% to 100%. While Ep12 (from HPV 18), Ep13 (from HPV 16), Ep14 (from HPV 33), and Ep15 (from HPV 45) are specific to individual genotypes, they showed conservation levels ranging from 74% to 94% across the other three genotypes. The E7-derived epitopes demonstrated moderate conservation across the four HPV genotypes. Consequently, we meticulously selected 4 Tc-cell (Ep5–Ep8) and B-cell (Ep16–Ep19) epitopes of the E7 oncoprotein, ensuring that these epitopes are associated with at least one of the four HPV genotypes, thereby provoking an immune response that specifically targets the E7 oncogenic protein. Human papillomavirus, as a DNA virus, exhibits relatively slow mutation rates compared to RNA viruses. 99 The E7 protein contains short, conserved motifs, with mutational analyses indicating that many residues are intolerant to substitution without loss of oncogenic function. 100 The chosen E7 and L1 epitopes exhibit a significant degree of conservation among high-risk genotypes. However, infrequent mutations or shifts in genotype may arise over time, highlighting the importance of ongoing surveillance.
Percentage conservation of predicted epitopes (Ep1–Ep19), which include Tc-, Th-, and B-cell epitopes sourced from L1 and E7 proteins, across four high-risk HPV genotypes: HPV 16, 18, 33, and 45. The conservation of L1-derived epitopes across all genotypes highlights their significant potential for cross-protection. E7-derived epitopes exhibited increased variability but were selected to ensure representation in at least one genotype.
The distribution of HLA alleles varies across geographical regions and ethnic populations worldwide. 101 Therefore, the chosen Tc and Th epitopes and their corresponding HLA binding alleles were analyzed using data from the IEDB database to determine the extent of population coverage, focusing on the most effective binders for the predicted epitopes. Population coverage ranged from 65.08% to 90.31% across various countries and 60.77% to 82.52% across different ethnic groups, based on highly interacting MHC alleles (Supplementary Tables S6 and S7). The combined Tc-cell and Th-cell epitopes demonstrated significant population coverage across various countries, such as Japan (85.25%), Italy (90.3%), the United States (81.62%), France (82.92%), Russia (81.66%), England (83.68%), Morocco (85.82%), Argentina (81.99%), Taiwan (80.78%), and Germany (81.29%) (Figure 3A). The analysis of population coverage among various ethnic groups revealed significant representation, with East Asia at 81.54%, Southeast Asia at 71.82%, Europe at 79.13%, West Africa at 79.13%, the West Indies at 82.52%, and North America at 81.64% (Fig. 3B). In our study, certain populations in underrepresented regions such as Northeast Asia, Southwest Asia, East Africa, Central Africa, South America, and Southeast Asia may harbor high-frequency alleles that were not fully captured by our selected allele panel. These regions possess distinctive alleles that could influence local coverage estimates. Nonetheless, this targeted selection strategy ensures that population coverage calculations reflect biologically relevant immune recognition. Overall, the selected vaccine peptide demonstrated extensive population coverage, exceeding 75% globally, with strong representation across diverse ethnic groups and geographical regions, underscoring its potential for broad-spectrum vaccine efficacy.
Global population coverage analysis of chosen Tc- and Th-cell epitopes based on HLA alleles utilizing the IEDB web server. (A) Coverage across countries. (B) Coverage across ethnic groups.
Modeling Tc-cell epitopes and docking with MHC I alleles
The molecular docking analysis using ClusPro 2.0 demonstrated strong binding interactions between the Tc-cell epitopes and the matching HLA alleles. The HLA-B*44:03 (PDB code: 4JQX) allele was selected for the EEYDLQFIF and NELDPVDLL epitopes (Figure 4A and B). The HLA-B*35:01 allele (PDB code: (was found to bind with the HVEEYDLQF and QAQPATADY epitopes (Figure 4C and D). Furthermore, the alleles HLA-A*68:01 (PDB code: 6PBH), HLA-A*24:02 (PDB code: 7JYV), and HLA-A*02:03 (PDB code: 3OX8) were selected for the epitopes LTVGNPYFR, KFGFPDTSF, and LLMGTLGIV, respectively (Figure 4F, G, and I). Subsequently, molecular docking-based binding analysis was validated using co-crystallized experimental peptides as positive controls (Figure 4E and H). The Tc-cell epitopes exhibited lower binding affinities compared to the positive controls for the respective alleles. However, the Tc-cell epitopes exhibited a higher number of hydrogen bonds (H-bonds), which contribute to the stability of the epitope-HLA interactions, except for HVEEYDLQF, which formed three conventional H-bonds. The QAQPATADY epitope most significantly formed eight hydrogen bonds in its interaction with the HLA-B*35:01 allele, surpassing the six hydrogen bonds established by its positive control. Although QAQPATADY exhibited a lower binding affinity in comparison to the positive control, the higher number of H-bonds formed by QAQPATADY suggests that the epitope could still engage effectively with the HLA-B*35:01 allele, potentially leading to a stable immune response despite the reduced binding affinity. Furthermore, the Tc-cell epitopes, along with the positive controls, demonstrated notable hydrophobic interactions with the HLA receptors (Table 1 and Figure 4). To enhance the validation of the docking results, the HDOCK server was employed, yielding binding affinity scores for the Tc-cell epitope-HLA complexes that ranged from −274.26 to −159.16 (Table 1). Finally, the Gibbs free energy values obtained from the PRODIGY server corroborated the robust interactions observed between the selected Tc-cell epitopes and MHC class I alleles. Among the analyzed complexes, HLA-B*35:01–HVEEYDLQF exhibited the highest binding affinity with a ΔG of −14.1 kcal/mol, while HLA-B*44:03–EEYDLQFIF showed the lowest at −7.5 kcal/mol. Negative ΔG values, measured in kcal/mol, signify thermodynamically favorable interactions, whereas lower ΔG values indicate stronger binding affinities. Taken together, the findings demonstrate the reliability and appropriateness of the chosen epitopes in generating immunological responses via HLA-mediated antigen presentation.
Analysis of molecular docking simulation involving the predicted epitope and its corresponding alleles: (A) EEYDLQFIF and (B) NELDPVDLL to the HLA-B*44:03 allele groove, (C) HVEEYDLQF, (D) QAQPATADY, and (E) Control to the HLA-B*35:01 allele groove, (F) LTVGNPYFR to the HLA-A*68:01 allele groove, (G) KFGFPDTSF and (H) Control to the HLA-A*24:02 allele groove, and (I) LLMGTLGIV to the HLA-A*02:03 allele groove. The epitopes are illustrated using ball-and-stick representations, while interacting residues are shown as sticks. Conventional H-bonds are represented by green-colored lines, pi–pi/pi–alkyl stacking bonds are indicated with pink-colored lines, pi–sigma bonds are depicted as violet-colored lines, carbon-H-bonds are shown as white-colored lines, attractive charges are represented by gold-colored lines, and unfavorable bumps are marked with red-colored lines.
Interaction of Tc-cell epitopes with MHC class I alleles.
Analysis of MHC clusters
Cluster analysis of interacting MHC alleles revealed functional relationships among 25 MHC I and 3 MHC II alleles associated with the selected epitopes. Figure 5A and C illustrates the cluster analysis of MHC I and MHC II alleles as heatmaps. Furthermore, Figure 5B and D represents a sophisticated tree map emphasizing the cluster analysis of MHC I and MHC II.
Outcomes of the MHC allele cluster analysis. (A) Heatmap illustrating MHC I cluster analysis. (B) Advanced tree map depicting MHC I cluster analysis. (C) Heatmap representing MHC II cluster analysis. (D) Advanced tree map showcasing MHC II cluster analysis. In the heatmap, red areas indicate stronger interactions among the clustered HLA alleles, while yellow areas represent weaker interactions.
Multi-epitope vaccine formulation
We chose a total of 19 epitopes that had the best immunological properties for designing an MEV. These epitopes included cytotoxic T lymphocyte (CTL), helper T lymphocyte (HTL), and B-cell epitopes. The eight CTL epitopes (9-mer) were linked using AAY linkers, while the two HTL epitopes (15-mer) were joined with GPGPG linkers, and the nine B-cell epitopes (16-mer) were joined via KK linkers. Furthermore, the TLR4 agonist, the 50S ribosomal L7/L12 protein (130 amino acids), was integrated at the amino terminus of the vaccine construct using an EAAAK linker (Figure 6). Intrinsic disorder analysis revealed several unstructured regions within the vaccine protein. According to IUPred3, residues 250 to 260 and 380 to 395 were predicted as disordered (Supplementary Figure S1), while DisEMBL identified residues 45 to 60, 111 to 119, and 377 to 389 as flexible regions. Notably, the disordered regions overlap with the linker sequences (GPGPG, KK), indicating that the linkers provide sufficient local flexibility to maintain proper spatial separation of epitopes while preserving overall protein stability.
Graphical illustration of the final multi-epitope vaccine assembly, depicting CTL, HTL, and B-cell epitopes interconnected by suitable linkers and an adjuvant.
Physicochemical and immunological properties of vaccine
The antigenicity assessment confirmed that the MEV construct is antigenic in both adjuvant-containing and adjuvant-free forms. The construct was also predicted to be non-toxic, non-allergenic, and a probable immunogen. Cross-reactivity analysis using the AllergenOnline database revealed a maximum of 36.2% identity over an 80-aa window, which is below the 50% to 70% threshold, indicating low allergenic cross-reactivity. Physicochemical analysis indicated that the MEV is acidic, stable, and thermostable, with a favorable in vivo half-life and aliphatic index. It also displayed hydrophilic and soluble properties in E. coli (Table 2). Collectively, these features support its potential as a vaccine candidate.
Evaluation of physicochemical properties, immunogenicity, antigenicity, toxicity, allergenicity, and solubility of the designed multi-epitope vaccine construct.
Prediction of secondary structure
The secondary structural characteristics of the MEV were analyzed using PSIPRED, indicating an architecture with a greater proportion of coil (50.47%) than helix (36.28%) and strand (13.25%) (Figure 7A). Furthermore, we conducted the secondary structure analysis using the SOPMA tool, providing a comparable outcome where random coils (42.79%) were more abundant than alpha-helixes (41.16%) and extended strands (16.05%) (Figure 7B).
Diagrammatic representation of the two-dimensional structural elements of the designed multi-epitope vaccine construct utilizing (A) PSIPRED and (B) SOPMA online server.
Tertiary structural modeling, refinement, and validation
The Robetta server generated five three-dimensional models for the designed MEV construct. Figure 8 highlights the final MEV construct, with epitopes and the adjuvant represented in ribbon style and linkers depicted using the ball-and-stick model, each shown in distinct colors.
Epitope mapping and structural modeling of the MEV construct. The adjuvant is presented in olive color, and the adjuvant connected by the EAAAK linker is displayed in purple color, while the CTL epitopes are highlighted in blue color linked by the AAY linkers in green color. Moreover, the cyan color depicts HTL epitopes connected by GPGPG linkers in red, whereas the B-cell epitopes are rendered in black, joined by KK linkers in orange.
In the pre-refinement phase, the PROCHECK tool facilitated the structural evaluation of five crude 3D models, with model 4 emerging as the best, demonstrating the highest percentage (83.8%) of residues situated in the most favored region of the Ramachandran plot. Furthermore, the additionally allowed, generously allowed, and disallowed regions contained 12.9%, 1.6%, and 1.6% of the residues, respectively (Figure 9A). Consequently, we refined the selected crude model using the GalaxyRefine server and retrieved 5 models with various parameters (Supplementary Table S8). The post-refinement quality analysis identified model 2 as the best, exhibiting the highest percentages of residues in the most favored region of the Ramachandran plot (91.5%) (Figure 9B). The additional allowed, generously allowed, and disallowed regions accounted for 4.7%, 2.5%, and 1.4% of the residues, respectively. Based on the established criteria for Ramachandran plot analysis, a model qualifies to have good quality if 90% of its residues fall within the most favored regions. Moreover, we validated the overall structural quality of the MEV construct using the ProSA-web server, where the crude model’s Z-score was found to be −9.59, which improved to −8.91 following the refinement process (Figure 9C and D). Therefore, the findings suggested a high level of accuracy in the molecular conformation of the refined model compared to the initial model. Figure 10 depicts the pre- and post-refinement vaccine models’ superimposition to compare their quality.
Structural evaluation of the three-dimensional MEV construct utilizing the Ramachandran plot and ProSA-web server before and following the refinement procedure. (A) Ramachandran plot of the initial model, showing 83.8% of residues in the most favored region. (B) Ramachandran plot of the refined model, showing 91.5% of residues in the most favored region. (C) ProSA-web Z-score of the initial model: −9.59. (D) ProSA-web Z-score of the refined model: −8.91.
The crude and refined three-dimensional vaccine models are represented in dark cyan and dark blue, respectively, to facilitate their structural comparison. The models were superimposed using Discovery Studio 2024.
Conformational B-lymphocyte epitope prediction
The continuous (linear) and discontinuous (conformational) epitopes of B cells were forecasted using the default parameters of the Ellipro server. Twelve conformational epitopes were identified, with sequence locations and residues spanning from 3 to 75 in the vaccine construct and scores varying from 0.50 to 0.78 (Table 3). Interestingly, the graphical representation of these epitopes indicated that their presence was observed at the terminal ends, which demonstrates less constraint and greater exposure in comparison to the tightly folded core (Figure 11A to L). This heightened exposure enables immune cells, including B-lymphocytes, to access these areas more readily, promoting antibody binding and activating the immune response. In addition, 14 linear B-cell epitopes comprising residues from 4 to 33 were identified in the vaccine, exhibiting scores between 0.522 and 0.803 (Supplementary Table S9 and Figure S2). A linear epitope within the vaccine construct was found to overlap with the discontinuous epitopes found in the developed vaccine.
The information set includes residues, size, and scores for anticipated conformational B-cell epitopes on the multi-epitope vaccine construct.
The yellow-colored region (A-L) of the vaccine construct illustrates the mapping of conformational B-cell epitopes, highlighting each epitope that encompasses residues 3 to 75, with score values between 0.50 and 0.78.
Protein-protein docking of vaccine-immunological receptors
Protein-protein docking studies using the ClusPro 2.0 and HDOCK platforms assessed the interaction energy and dynamics between the MEV construct and the human immune receptors TLR2 and TLR4. In each docking simulation, the Cluspro 2.0 server produced 30 docked complexes involving many cluster members with the lowest energy configurations (Supplementary Table S10). Among these, we prioritized the models that successfully engaged the receptor, showing the lowest energy scores (a sign of high binding affinity) and optimal functional interactions. The docking results reported that the MEV construct exhibited the lowest energy scores in its interactions with TLR2 and TLR4, with recorded values of −868.2 (69 cluster members) and −894.4 (cluster member 31), respectively. Despite the incorporation of 50S ribosomal protein L7/L12 as an adjuvant, an acknowledged TLR4 agonist, docking analysis indicated that the MEV construct also effectively interacted with TLR2, showing comparable binding affinity and suggesting a robust interaction with both immunological receptors. 102 Furthermore, the HDOCK server produced 100 structures for each docking process, with the top 10 presented in Supplementary Table S11. The vaccine construct had the lowest binding scores when it interacted with TLR2 and TLR4, with values of −247.46 and −245.97, respectively. Graphical representations of the top docked complexes were generated using UCSF Chimera (Figure 12A and B), and PDBsum showed the molecular interactions and the residues involved between the MEV construct and TLRs (Figure 12C and D). For each complex, information regarding the quantity of interacting residues, interface area, H-bonds, salt bridges, and non-bonding contacts formed between MEV and TLRS is shown in Figure 12E and F. In contrast to the TLR2-MEV complex, the TLR4-MEV complex exhibited the most significant interactions, with a higher number of H-bonds, salt bridges, and nonbonded contacts. Salt bridges between the protein chains of the dock complexes are represented in red, H-bonds are depicted in blue, disulfide bonds are shown in yellow, and unbound interactions are illustrated with dotted orange lines.
Molecular docking-based binding strength and interaction dynamics between the vaccine construct and toll-like receptors. (A) Interaction of TLR2 (orange) with MEV (cyan). (B) Interaction of TLR4 (dark yellow) with MEV (cyan). (C) Amino acid residues participating in the interaction between MEV and TLR2. (D) Amino acid residues participating in the interaction between MEV and TLR4. (E) Quantification of interaction parameters between MEV and TLR2, including the quantity of interacting residues, interface area, H-bonds, salt bridges, and nonbonding contacts. (F) Corresponding interaction parameters (number of salt bridges, H-bonds, and nonbonding contacts) between MEV and TLR4. These analyses offer comprehensive insights into the structural dynamics and intermolecular interactions facilitating MEV binding to TLR2 and TLR4.
Evaluation of biophysical characteristics through molecular dynamics simulation
To assess the conformational stability, flexibility, and compactness of the MEV construct, MEV-TLR2, and MEV-TLR4 complexes, 100-ns simulations were conducted, and the RMSD of alpha-carbon and RMSF were analyzed using the simulation trajectories.
RMSD analysis
The RMSD value, reflecting atomic displacement relative to a reference frame, revealed distinct stability patterns across the apo MEV, MEV–TLR2, and MEV–TLR4 systems. The MEV construct fluctuated between 2.20 and 6.69 Å during the first 15 ns, reaching a peak of 6.69 Å before rapidly declining to 5.1 Å at 16 ns. Thereafter, it stabilized after ~16 ns and remained steady with only very minor fluctuations, converging to a mean RMSD of 5.55 ± 0.44 Å. The MEV–TLR4 complex exhibited a similar pattern, with an average RMSD of 5.19 ± 0.81 Å, and a range spanning from 1.58 to 8.39 Å. It achieved stabilization after ~16 ns and remained stable until ~84 ns, followed by a transient increase that peaked at 8.39 Å near 97 ns, before returning to stability by the final frame. Conversely, the MEV–TLR2 complex exhibited a higher mean RMSD of 7.70 ± 1.09 Å, with values ranging from 2.40 to 10.46 Å. Its RMSD increased progressively during the first 31 ns, reaching 9.43 Å with some fluctuations, and then entered a transient stabilization phase from 31 to 84 ns, maintaining a similar range but showing minor downward fluctuations. Shortly after 84 ns, the RMSD reached its highest peak of 10.46 Å, followed by a slight downward adjustment, and then continued an overall upward trend until the end of the simulation. Unlike apo MEV and MEV–TLR4, MEV–TLR2 did not achieve long-term stabilization, indicating higher structural flexibility throughout the simulation. Overall, these observations suggest that TLR4 binding stabilizes the MEV construct, whereas TLR2 interaction allows more dynamic behavior.
RMSF analysis
The RMSF profiles of the analyzed vaccine–immune receptor complexes revealed distinct residue-specific fluctuation patterns (Figure 13B). For the MEV construct, the average RMSF was 2.34 ± 1.04 Å. Residues from Lys_74 to Asp_98 showed moderate flexibility, representing a relatively stable yet dynamic portion of the structure, whereas residues surrounding Lys_248 showed the highest fluctuations, reaching an RMSF peak of 6.75 Å. This highly flexible region corresponds to a loop/coil segment overlapping with the Th-cell epitope region. Loop regions are inherently more dynamic, and such flexibility is consistent with previous reports that MHC class II ligands often occupy solvent-exposed regions with higher coil content. 103 This structural adaptability may facilitate efficient antigen processing, presentation, and recognition by Th-cell receptors. In addition, several residues within B-cell epitope regions, including Gly_322, Ser_343, and Ser_381, together with their surrounding residues, displayed moderate flexibility (RMSF values around 4 Å). Such flexibility may enhance solvent exposure and epitope accessibility, 104 which are essential for effective B-cell receptor recognition and antibody binding, while the overall MEV construct remained largely stable during simulation. Most residues remain below 4 Å, indicating restricted local flexibility and stable core regions. On the other hand, the RMSF profile of the MEV–TLR2 complex displayed moderate flexibility throughout, with values spanning from 0.814 to 5.603 Å and an average RMSF of 2.32 ± 0.97 Å. Most residues up to position 535, except for Lys_218, displayed fluctuations below 4 Å, indicating structural stability with minimal mobility. Residues beyond position 535 showed some flexibility but stabilized by the end of the simulation, reflecting a region with localized flexibility and overall stability in the receptor. In contrast, the MEV–TLR4 complex had a slightly higher average RMSF of 2.68 ± 1.69 Å, indicating comparatively greater flexibility. While residues up to 529 remained stable (<4 Å), pronounced fluctuations were observed between residues 535 and 724, likely corresponding to regions of TLR4 undergoing conformational adjustments during vaccine binding. Beyond residue 724, the complex stabilized, suggesting that a stable receptor–vaccine configuration followed initial dynamic fluctuations. Overall, these results indicate that both complexes maintained structural stability, with the MEV–TLR4 complex showing greater localized flexibility, suggesting receptor adaptability during ligand engagement.
Evaluation of the conformational stability, flexibility, and binding interactions of TLR2 and TLR4 with the multi-epitope vaccine through dynamics simulations. (A) The RMSD of the TLR2 receptor (blue) and the TLR4 receptor (green) when bound to the vaccine (ligand), as well as the vaccine alone (gray), indicates enhanced stability upon binding. (B) The RMSF of TLR2 (blue) and TLR4 (green) bound to the vaccine, compared to the vaccine alone (gray) illustrates structural flexibility upon binding.
MM/GBSA binding free energy calculation
The free binding energy (MM/GBSA) calculation of the vaccine-immune receptor complexes estimated a binding energy of −82.86 kcal/mol for the MEV-TLR2 complex. The computed VDW, ELE, GB, and SA values for this complex were −104.75, −101.19, 137.68, and −14.59 kcal/mol, respectively (Figure 14). Conversely, the “MEV-TLR4” complex demonstrated a total free binding energy of −76.72 kcal/mol. The values obtained for the VDW, ELE, GB, and SA were −202.82, 823.64, −670.65, and −26.89 kcal/mol, respectively.
The MM/GBSA free binding energy analysis for the MEV-TLR2 and MEV-TLR4 complexes.
Disulfide engineering
The DbD2 server listed 42 eligible pairs of amino acids that could potentially form disulfide bonds (Supplementary Table S12). In the residue screening process, we identified only six pairs that are appropriate for substitution with cysteine to form a disulfide bridge (Figure 15 and Table 4). Moreover, SASA values showed minimal changes in the immunogenic exposure of B-cell epitopes after disulfide engineering, and structural mapping confirmed that all B-cell epitopes remained accessible on the protein surface in both constructs (Supplementary Table S13 and Figure S3).
Mutant model for the multi-epitope vaccine construct designed by introducing cysteine residues.
Pairs of residues mutated as cysteine in the mutant vaccine model.
Immune simulation
The C-IMMSIM tool effectively forecasted immunological responses in humans following three vaccination doses, given at 4-week intervals. The second and third vaccinations elicited an enormous spike in immune response when compared with those observed following the initial vaccination. The immune response demonstrated an enormous spike in IgM, and the antigen count reached a peak of 670000 per cell following the initial injection; however, it subsequently showed a notable decline. The concentration of IgM and IgG reached its peak at 50 days, reaching approximately 650000 counts per cell, before declining by day 350 (Figure 16A). Furthermore, following the administration of the second and third vaccine doses, there was a significant increase in the overall count of B cells. The injection of the vaccine led to an increase in the population of active B cells (Figure 16B and C). Vaccination also induced an upsurge in the number of helper T cells (TH) and active TH cells throughout the immune stimulation (Figure 16D and E). Furthermore, during the simulation of immunologic response, cytotoxic T cells (TC) were identified, reaching a peak count of over 1150 cells per mm³ (Figure 16F). Following the administration of the chimeric antigen, there was an immediate rise in the count of active cytotoxic T cells (TC), while a sharp reduction in the count of resting TC cells was observed, as illustrated in Figure 16G. Vaccination stimulated strong innate immune responses in dendritic cells (DC), and the observation that these cells maintained normal function for a year straight suggests that the immunity remained unaffected (Figure 16H). The vaccination prompted cytokine responses, notably including interferon-g, which peaked at around 420000 ng/ml, subsequently declining to approximately 400000 ng/ml after 50 days. Interleukin-12 (IL-12) and interleukin-10 (IL-10) exhibited the lowest levels in this simulation, approaching 50000 ng/ml, while tumor growth factor-b measured approximately 130000 ng/ml (Figure 16I). Initially, cytokines were detected but disappeared 60 days after the vaccination. Moreover, the developed MEV reliably triggered elevated levels of IL-2 following each administration, reaching a peak of around 650000 ng/mL. Enhanced immune responses were noted after the second and third doses. The findings demonstrate that the vaccine generates a strong immune response to our targeted HPV genotypes.
Simulation of immunological response for the multi-epitope-based subunit vaccine. (A) Variations in IgM and IgG antibody titers following the administration of the vaccine. (B) B cell enumeration. (C) B cell enumeration by state. (D) TH cell enumeration. (E) TH cell enumeration by state. (F) TC cell enumeration. (G) TC cell enumeration by state. (H) Immune response of dendritic cells (DC) following vaccination. (I) Levels of interleukins and cytokines.
In silico cloning of the final vaccine construct
Back translation and optimization of the codon usage of the vaccine amino acid sequence in E. coli (strain K12) were performed by JCat for the maximal level of vaccine protein expression. The multi-epitope vaccine protein (430 amino acid residues) produced 1302 nucleotide sequences. The CAI and GC content of the MEV construct were 0.97% and 50.69%, respectively, both falling within the commonly reported ranges for efficient protein expression in E. coli (CAI 0.8–1.0, GC 30%–70%), 105 indicating a high likelihood of successful expression. SnapGene software was utilized to introduce adapted codon sequences into the pET28a (+) vector by assisting EcoRI and BamHI restriction enzymes as the start and end cut points, respectively, and obtained a cloned vaccine of 6667 base pairs (bp) (Figure 17).
The in silico cloning of the formulated multi-epitope vaccine construct into the pET-28a (+) vector. The red region illustrates the gene sequence of our developed vaccine construct, while the black section indicates the backbone of the vector. All colored arrows show the location and direction of gene expression: green represents the kanamycin resistance gene, and yellow symbolizes the origin of replication.
Forecasting post-translational modifications
PTM analysis of the vaccine construct was performed to identify potential sites that could affect stability, function, or immunogenicity. MusiteDeep analysis highlighted regions requiring modifications, including phosphorylation (10), ubiquitination (10), glycosylation (2), acetylation (3), methylation (4), and hydroxylation (1). Comprehensive predictions of PTMs were carried out utilizing specialized servers, with the NetNGlyc-1.0 server identifying only one N-glycosylation site (N: 256). The analysis of phosphorylation modifications through the NetPhos 3.1 web tool revealed that 25 specific locations within the construct (Ser: 8, Thr: 12, Tyr: 5) were identified as phosphorylated. No N-acetylation sites were identified through the NetAcet-1.0 server analysis.
Discussion
Despite widespread prophylactic vaccination, HPV-associated CC remains a global health challenge with high mortality rates. 106 Current prevention strategies rely on six authorized VLP-based vaccines; however, no therapeutic option exists.4,61 Therefore, a vaccine with both prophylactic and therapeutic properties would avert HPV infections and provide a treatment for established infections, greatly enhancing disease management. 8 In contrast to traditional vaccines, peptide-based vaccines offer enhanced safety, greater stability, and simplified production processes, as they do not depend on the cultivation of pathogens. These are notably adaptable, focusing on conserved viral regions to maintain efficacy against emerging strains, thereby minimizing the likelihood of immune evasion.107,108 This study provides valuable insights into the potential development of peptide-based vaccines targeting four high-risk HPV genotypes (HPV 16, 18, 33, and 45), employing an immunoinformatics-based approach. The L1 major capsid protein and E7 oncoprotein were selected as key targets due to their roles in viral entry and cancer progression. L1, forming the viral capsid with 72 pentameric capsomeres, is crucial for host cell interaction and immune recognition, induces strong neutralizing antibody responses and forms the basis of current prophylactic vaccines (Cervarix, Gardasil), but lacks therapeutic efficacy. 109 On the other hand, E7 is consistently expressed in nearly all CC cells and has become a prime target for immunotherapy in HPV-related cancers 8 Unlike existing L1-only prophylactic vaccines or the experimental E7-targeted therapeutic vaccines, 110 our strategy that employs a combined L1 and E7 peptide-based approach integrates both prophylactic and therapeutic immunological benefits. L1 facilitates VLP-mediated antibody responses to prevent infection, while E7 activates cytotoxic T-cell responses to destroy infected or malignant cells, successfully addressing the limitations of vaccines that concentrate only on L1 or E7 separately.
Predicting proper T-cell and B-cell epitopes is an essential and critical step in the development of MEVs, as T-cells recognize pathogen-derived antigens and carry out direct effector functions like cytotoxicity, while also aiding B-cells in regulating the development and maturation of antibody responses. 111 The ability of the vaccine to generate a strong and lasting immune response is primarily due to the meticulous choice of epitopes applying various immune filters. Strong immunogenicity was achieved through the selection of highly antigenic epitopes, while safety issues were addressed by omitting allergenic, toxic, and human-homologous epitopes. The ability to induce cytokines and the conservation of epitopes further enhanced immune activation and provided broad coverage across four HPV genotypes. Unlike previous studies that applied immunoinformatics for HPV vaccine design without addressing population coverage, this study demonstrates epitope-binding alleles with over 75% global population coverage.1,4,8
The Tc-cell epitopes from the L1 and E7 proteins showed high binding affinity to MHC class I molecules, as indicated by low docking scores. This strong interaction was primarily driven by hydrogen bonding at the HLA binding sites. However, combined strong–weak H-bond pairings can sometimes reduce the ligand binding affinity because of interactions with bulk water, which helps explain why the random strengthening of ligand-receptor H-bonds does not always correlate well with experimental binding affinity.112,113 Despite this, the Tc-cell epitopes exhibited moderate to high numbers of strong H-bonds and more negative docking scores, indicating that these H-bonds enhanced the binding strength.
Immunodominant multi-epitope-based vaccines are considered a promising therapeutic approach for addressing pathogenic and viral infections. The strategic incorporation of specific linkers and an adjuvant in the final MEV construct played a critical role in enhancing its immunogenic potential and structural stability. Linkers were carefully selected to ensure independent epitope presentation, proper protein folding, and minimal immunological interference between adjacent epitopes. The use of the EAAAK linker to fuse the 50S ribosomal protein L7/L12 adjuvant at the N-terminus contributed to a more potent immune response by promoting structural rigidity through its α-helical conformation, thereby efficiently separating functional domains.114,115 The inclusion of the 50S ribosomal protein L7/L12 adjuvant as a TLR4 agonist is supported by previous immunoinformatics-based HPV vaccine studies.1,61 In this study, Tc epitopes were joined using the AAY (Ala-Ala-Tyr) linker, as described by Bhatnager et al. 116 This linker was selected for its capacity to promote proteasomal cleavage in mammalian cells, thereby ensuring precise epitope processing. This approach reduces junctional immunogenicity while boosting MEV immunogenicity through the facilitation of effective T-cell recognition. Similarly, the GPGPG linker promotes Th-cell responses and mitigates junctional immunogenicity, as demonstrated in mouse models by Livingston et al. 117 In addition, the KK linker aids B-cell epitope presentation by serving as a target for Cathepsin B, a lysosomal protease crucial for MHC-II-restricted antigen processing. This mechanism inhibits the formation of antibodies against linearly connected peptide sequences, further improving immunogenicity. 118 Altogether, the careful integration of these linkers contributed to a rational MEV design, promoting targeted immune activation while minimizing immunogenic interference. Furthermore, intrinsic disorder analysis (IUPred3: 250–260, and 380–395; DisEMBL: 45–60, 111–119, and 377–389) highlights flexible linker regions (GPGPG, KK) that maintain epitope separation and accessibility, while the overall construct remains stable.
While evaluating epitopes is important, concentrating exclusively on them neglects other significant elements that affect the effectiveness and safety of vaccines. The interplay between vaccine components influences the immune response as a whole, underscoring the importance of thorough physicochemical and immunological assessments. The vaccine effectively stimulates the immune system without any adverse reactions, as demonstrated by its antigenic, immunogenic, non-toxic, and non-allergenic characteristics. Our constructed MEV demonstrated strong antigenicity (0.8865), exceeding previously reported HPV vaccine antigenicity scores, which range from 0.522 to 0.745.1,8,25,61 The high antigenicity observed in the vaccine is due to the optimal design of its epitopes, which are well-suited to interact with the immune system, triggering a strong immune response. While earlier studies on HPV MEV focused solely on antigenicity predictions through VaxiJen v2.0, our analysis expanded this scope by integrating both antigenicity (VaxiJen v2.0) and immunogenicity (VaxiJen v3.0) evaluation, thus providing a more comprehensive assessment of the construct’s immunogenic potential. The molecular weight of 46.16 kDa, which falls within an ideal range, supports efficient expression and purification processes, facilitating the production of the vaccine. 119 Its acidic nature and predicted half-life across systems indicate prolonged in vivo stability. 64 The vaccine construct is predicted to be stable, as reflected by an instability index of 25.27, which falls well below the threshold value of 40. 64 Thermostability and solubility are key to vaccine efficacy, particularly during storage and administration. A high aliphatic index suggests structural stability at varying temperatures, while a negative GRAVY score reflects hydrophilicity, supporting favorable solubility in physiological conditions. 120
The determination of the two-dimensional and three-dimensional structure of an MEV construct is essential in vaccine formulation, and a key challenge in structural biology is to pinpoint differences between theoretical predictions and experimental findings of protein structures. 121 The outcomes of secondary structure analysis indicate a greater proportion of coil than helix and strand. The crude vaccine model has undergone substantial optimization following meticulous fine-tuning and the presentation of relevant statistics for the Ramachandran plot. The overall quality evaluated by ProSA Z-score analysis indicated that the vaccine construct aligns with the Z-scores of previously established constructs derived from experimental X-ray crystallography and nuclear magnetic resonance data. 76
Toll-like receptors play a pivotal role in bridging innate and adaptive immunity by activating antigen-presenting cells (APCs) and promoting both cellular and humoral responses. In HPV infection, immune evasion is often mediated by the downregulation of TLR2, whereas TLR4 expression is significantly upregulated. TLR4 activates multiple intracellular pathways, including inflammasome-mediated signaling, which may restrict viral DNA integration while also contributing to cellular immunity and antibody isotype switching against HPV antigens.122,123 The strong binding affinity observed in molecular docking, corroborated by MD simulations, highlights the vaccine construct’s robust and stable interaction with TLRs, suggesting effective immune receptor engagement. The stability analysis demonstrated that the vaccine maintained its conformational integrity under physiological conditions, particularly in complex with TLR4, due to reduced structural deviations and fluctuations compared to TLR2. This suggests greater stability and potentially a more effective interaction with TLR4. The localized flexibility noted near Lys_248 in the TLR2 complex may have an important impact on receptor engagement. The observed dynamic behavior and positive interaction patterns suggest that the vaccine construct demonstrates structural competence for eliciting robust immune responses.4,44,115
The binding free energy calculations revealed that the MEV-TLR2 complex exhibited slightly more favorable binding energy due to the optimal interplay of van der Waals and electrostatic interactions. In contrast, the MEV-TLR4 complex demonstrated noteworthy stability due to van der Waals interactions, whereas pronounced electrostatic destabilization played a counteractive role in the interaction dynamics. The MM/GBSA analysis showed that the association process between the vaccine and receptors is spontaneous and thermodynamically favorable because the ΔTOTAL values were negative. Furthermore, the vaccine construct was subjected to disulfide engineering for creating a mutant variant. Six amino acid pairings in the MEV construct met the selection requirements of bond energies below 2.2 kcal/mol, making them appropriate for the generation of mutant proteins. Disulfide engineering has been shown to enhance vaccine design in several studies.111,124,125 The SASA analysis and epitope mapping indicated that the B-cell epitopes were accessible after disulfide engineering, confirming that the construct preserved its immunogenic potential while enhancing stability.
The vaccine effectively triggered both innate and sustained adaptive immunity, resulting in the production of antibodies, B cells, T cells, dendritic cells, and a variety of cytokines, such as interleukin-2. Memory T cells and B cells were induced, with B-cell-mediated immunity persisting for 1 year. Distinct characteristics were observed in the activation of TH, leading to the subsequent production of IFN-γ and IL-2 following the initial injection, with their levels remaining elevated through repeated exposure to antigen. This observation indicates an activation of the humoral immune response, as there is an increased synthesis of TH cells and immunoglobulins (IgM and IgG). The production of recombinant proteins requires the expression of the vaccine construct in an appropriate vector [116]. In the present study, the sequence of the vaccine was reverse transcribed and optimized for the E. coli (K12 strain) before being cloned into the pET-28a(+) vector to enhance translation efficiency. The experiment showed a satisfactory expression level of the vaccine, suggested by an average GC content of 50.69% and a CAI value of 0.97. Additional in vitro and in vivo research on expression vectors is necessary for the cloning of our designed vaccine, as they may be necessary for large-scale production. Analysis of post-translational modifications revealed that phosphorylation and ubiquitination were the predominant modifications in the MEV construct. Phosphorylation may augment epitope presentation by MHC Class I molecules, facilitating CTL recognition and amplifying immunological responses. Ubiquitination promotes the proteasomal breakdown of the vaccine antigen, producing peptides that are effectively presented to CD8+ T cells, thus enhancing CTL-mediated immunity and overall vaccination immunogenicity. 95 Together, these alterations may improve epitope accessibility and the effectiveness of the vaccine.
This study pioneers an immunoinformatics-driven vaccine targeting the four most prevalent high-risk HPV genotypes (16, 18, 33, 45) with both prophylactic and therapeutic potential. Our findings demonstrate the capability of the vaccine construct to generate both cell-mediated and antibody-mediated responses, integrate all critical immunogenic, physicochemical, and structural features, and successfully trigger an immunological response to tackle HPV infection and CC progression. This vaccine will be subjected to clinical trials in both animal and human subjects to assess its effectiveness. Finally, this study enhances the importance of computational methods in modern biomedical studies by deepening our insights into HPV and CC vaccine development and highlighting their broader potential in combating other challenging infections.
Limitations
While these findings show potential, it is crucial to recognize the particular limitations inherent in our study. Initially, the computational predictions and analyses encompass various aspects, including antigenicity, allergenicity, immunogenicity, molecular docking, molecular dynamics simulations, MM/GBSA binding energy calculations, disulfide engineering, immune simulation, and post-translational modification analyses. These aspects rely on in silico tools that contain significant information but require experimental validation. Further in vitro and in vivo studies are crucial to confirm the immunogenicity and protective efficacy of the proposed multi-epitope vaccine construct. Second, while we employed a range of computational techniques to ensure the stability of the vaccine and its interaction with immune receptors, it is important to acknowledge that variations in real biological systems may affect the expected outcomes. Finally, while the multi-epitope vaccine construct has been optimized for bacterial expression, it may be required to explore alternative expression systems, like yeast or mammalian cells, to enhance expression efficiency and immunogenic performance.126,127
Conclusion
The rising incidence of CC associated with HPV shows the urgent need for effective treatment approaches since existing vaccines do not target prevailing infections. The present study highlighted the advantages of peptide-based vaccines, including their safety, specificity, and ability to adapt to emerging variants. With the help of immunoinformatics, we created a multi-epitope vaccine that targets the L1 major capsid protein and the E7 oncoprotein of HPV genotypes 16, 18, 33, and 45. The vaccine contained Tc-cell, Th-cell, and B-lymphocyte epitopes, along with an immunostimulatory adjuvant that made immune responses prolonged and stronger. Computational evaluations validated its antigenicity, immunogenicity, non-allergenicity, and advantageous physicochemical properties. Docking analysis demonstrated significant interactions between the vaccine and TLRs, which are crucial for immune activation. In addition, the conformational stability, flexibility, and binding strength of vaccine-TLR complexes were demonstrated by the MD simulations and MM/GBSA binding free energy calculations. The analysis of binding interactions revealed a positive correlation between their strength and the anticipated outcomes of immune responses. Following vaccination, there were noteworthy sustained responses from both B-cells and T-cells, accompanied by the activation of the innate immune system. Codon optimization facilitated successful expression in E. coli vectors. The outcomes establish a strong basis for the continued development of vaccines and the initiation of clinical trials, offering significant potential for the prevention and treatment of HPV-associated CC on a global scale. Further trials are essential to convert these computational predictions into viable immunotherapeutic solutions.
Supplemental Material
sj-docx-1-bbi-10.1177_11779322251391076 – Supplemental material for Designing a Multi-Epitope Vaccine Against HPV 16, 18, 33, and 45 Targeting L1 and E7 Proteins: An Immunoinformatics Approach for Cervical Cancer Prevention and Therapy
Supplemental material, sj-docx-1-bbi-10.1177_11779322251391076 for Designing a Multi-Epitope Vaccine Against HPV 16, 18, 33, and 45 Targeting L1 and E7 Proteins: An Immunoinformatics Approach for Cervical Cancer Prevention and Therapy by Md. Touki Tahamid Tusar, Niamul Haq, Hafizur Rahman Gazi, Raduyan Farazi, Mamun Bhuya, Md. Enamul Haque, Md. Golzar Hossain and Abdullah-Al-Jubayer in Bioinformatics and Biology Insights
Supplemental Material
sj-xlsx-2-bbi-10.1177_11779322251391076 – Supplemental material for Designing a Multi-Epitope Vaccine Against HPV 16, 18, 33, and 45 Targeting L1 and E7 Proteins: An Immunoinformatics Approach for Cervical Cancer Prevention and Therapy
Supplemental material, sj-xlsx-2-bbi-10.1177_11779322251391076 for Designing a Multi-Epitope Vaccine Against HPV 16, 18, 33, and 45 Targeting L1 and E7 Proteins: An Immunoinformatics Approach for Cervical Cancer Prevention and Therapy by Md. Touki Tahamid Tusar, Niamul Haq, Hafizur Rahman Gazi, Raduyan Farazi, Mamun Bhuya, Md. Enamul Haque, Md. Golzar Hossain and Abdullah-Al-Jubayer in Bioinformatics and Biology Insights
Footnotes
Acknowledgements
The authors express their gratitude to the Department of Biotechnology and Genetic Engineering, Faculty of Life Sciences, Gopalganj Science and Technology University, Gopalganj-8100, Bangladesh, for permitting the execution of computational studies and for supplying all the resources required for the study.
ORCID iDs
Ethical considerations
This article does not contain any studies with human or animal participants.
Consent to participate
Not applicable.
Consent for publication
Not applicable.
Author contributions
Funding
The authors received no financial support for the research, authorship, and/or publication of this article.
Declaration of conflicting interests
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Data availability statement
All data generated or analyzed during this study are included in this published article.
Supplemental material
Supplemental material for this article is available online.
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
