Abstract
The current research investigated the development of a multi-epitope mRNA vaccine against the rabies virus on the basis of viral proteomes via the use of bioinformatic tools and reverse vaccinology. The aim of this study was to address the limitations of the currently available rabies vaccine by eliciting strong and long-lasting humoral and cellular immune responses. The cytotoxic T lymphocytes (CTLs), helper T lymphocytes (HTLs), and linear B-cell epitopes (LBLs) were mapped and prioritized from four top-ranking vaccine targets (nucleoprotein, phosphoprotein, matrix, and glycoprotein) that were highly antigenic, nonallergenic, nontoxic, and nonhuman homologs. The selected epitopes exhibited strong binding affinity to high-frequency HLA alleles, as evidenced by highly negative ΔG values and low dissociation constants, predicting efficient T-cell recognition and broad population coverage (96.01% globally). A single mRNA construct encompassing 21 shortlisted epitopes (four CTL, four HTL, and thirteen LBL epitopes) was designed with appropriate linkers and the immunostimulatory 50 S ribosomal protein L7/L12 adjuvant. Physicochemical analysis revealed stable, soluble, and hydrophobic properties, with an overall Ramachandran score of 93.2%, an ERRAT quality factor of 94.724%, and a Z score of -5.39. Additionally, molecular docking and normal mode analysis demonstrated the strong binding affinity of the vaccine construct-TLR-4 complex, with a minimum energy of -1655.0 kcal/mol, which was maintained by 23 hydrogen bonds and 2 salt bridge interactions, indicating significant structural stability and stiffness. The structural integrity and stable interaction of the complex were validated through 200 ns molecular dynamics simulations, as evidenced by stable RMSD and radius of gyration values, minimal fluctuations in RMSF, consistent solvent-accessible surface area (SASA), and well-defined conformational transitions observed in principal component analysis (PCA). In silico immune simulation revealed the capacity of the vaccine to stimulate the release of high levels of immunoglobulin, TH, and TC and the release of cytokines. It also has the ability to produce long-lasting memory cells, induce macrophage activity, and promote natural killer cell and neutrophil production. Moreover, further validation, including codon optimization and mRNA secondary structure prediction, confirmed the stable structure and high level of expression in the host. Overall, this study proposed a promising multi-epitope-based mRNA vaccine as an innovative therapeutic candidate against rabies. However, experimental validations are needed with systemic animal studies.
Similar content being viewed by others
Introduction
Rabies is a major public health concern that causes severe suffering, death and economic disruption. Annually, rabies causes approximately 59,000 human deaths across over 150 countries. The highest burden occurs in Africa, where dog rabies remains endemic or epizootic and accounts for most human exposures1,2.
This disproportionate burden reflects systemic challenges in veterinary infrastructure and public health systems, as exemplified by Tunisia’s recent outbreak series documenting 25 human fatalities between 2021 and 20243. The continued prevalence of rabies in the modern era underscores the urgent need for innovative prevention strategies.
This zoonosis is caused by negative stranded RNA rabies virus (RV), a rhabdovirus of the genus Lyssavirus4. The viral genome encodes five proteins: nucleoprotein (N), matrix protein (M), glycoprotein (G), viral RNA-dependent RNA polymerase (L) and its cofactor phosphoprotein (P)5,6. Glycoprotein G is the major antigen responsible for the induction of protective immunity and pathogenicity7,8, and nucleoprotein N triggers rabies virus-specific T cells, facilitating the production of neutralizing antibodies and other immune mechanisms9,10. The structural M protein plays an important role in RV pathogenesis by regulating virus replication and facilitating cell-to-cell spread11. The IFN antagonist P protein is indispensable not only for viral replication but also for evasion of host innate immunity12.
The rabies virus infects the central nervous system of mammals, leading to encephalitis and ultimately death13. The transmission of RV occurs primarily through the saliva of infected animals, which can enter the body through bites, scratches, or contact with mucosa or fresh skin wounds. Once symptoms of rabies appear, the disease is always fatal, emphasizing the critical importance of prompt medical attention and vaccination to prevent infection14.
Collaborative initiatives involving organizations such as the World Health Organization (WHO) and the Global Alliance for Rabies Control are being implemented worldwide to eradicate rabies-related deaths by 2030. In fact, rabies prevention and eradication necessitate coordinated efforts by animal and human health authorities, with an emphasis on controlling rabies in dogs and wild animals and providing prompt prophylaxis to exposed humans. One challenge in the control of rabies is that the disease affects humans and livestock and can re-emerge via transmission out of wild animal reservoirs15.
Because rabies is of great interest and priority for the Tunisian government, especially after the last outbreak occurred, a cooperative framework has been strengthened to create key conditions to support successful improvements in public health. In fact, low vaccination rates among livestock and stray dogs, coupled with limited vaccine availability in rural areas, hinder effective control. Attempts to reduce rabies transmission through fox population control have proven counterproductive, disrupting ecological balance and increasing disease spread.
The prevention and treatment of rabies heavily rely on the development of effective vaccines. Commonly, the WHO recommends two main immunization regimens, namely, pre-exposure prophylaxis (PrEP) and post-exposure prophylaxis (PEP), as crucial measures for preventing human rabies. At present, safe Rabies cell culture vaccines (CCV) are widely used for the prevention of rabies in humans and animals as a post-exposure prophylaxis measure, superseding the nerve tissue vaccine (NTV) developed by Pasteur due its severe neurological side effects, the Live-attenuated vaccines that carry reversion risks, and adjuvanted vaccines (such as those based on aluminum) that may trigger autoimmune reactions16,17. These vaccines, although effective, encounter limitations such as incomplete immune responses, side effects, multiple dose vaccination requirement, expense, insufficient vaccine supply17, and cold chain requirements18,19. Nowadays, no single vaccination strategy is likely to be effective or feasible. Therefore, alternative approaches must be deployed to reduce the number of doses and increase the mass coverage of vaccination. The development of next-generation vaccines is underway, aiming to overcome these limitations20.
Recent advances in computational immunology and reverse vaccinology, focused on in silico and bioinformatics approaches for designing multi-epitope subunit vaccines, have revolutionized antigen selection and vaccine development strategies. These approaches have enabled the precise identification of broadly protective epitopes for SARS-CoV-221, rabies22, influenza viruses23, Zika24, tuberculosis25, malaria26, and various cancers (melanoma, lung, breast)27,28, while offering early evaluation of potential candidates targeting conserved immunogenic sites capable of eliciting both humoral and cellular immunity in a way that is more convenient, cost-effective, and time-efficient. Both in vitro29,30 and in vivo31,32,33 studies have validated the efficacy and prophylactic potential of multiepitope vaccines, with several related candidates already advancing to clinical trials34,35,36,37.
Encoding viral antigens into mRNA vaccines has shown effectiveness in preclinical models for infectious disease threats like Hendra virus (HeV)38 and HIV39. Such vaccines utilize optimized epitopes that trigger strong cellular and humoral immune responses. Alongside these developments, comparable mRNA-based methods for rabies virus glycoprotein G are incorporating computational methods for identifying conserved immunogenic epitopes. Nonetheless, significant obstacles still need solving: (1) the need for several booster shots to generate long-lasting immunity, (2) temporary side effects, (3) encapsulation and delivery systems, and (4) need for ultra-low cold storage.
The current research aimed to apply innovative immunoinformatic approaches for predicting consensus epitopes and developing broad-spectrum rabies vaccine candidates for prophylactic vaccination.
Methodology
Retrieval of sequences and consensus sequence generation
Rabies virus protein (G, P, NP, M and L) sequences were downloaded from the National Center for Biotechnology Information Database (NCBI) and then subjected to multiple sequence alignment with reference strain sequences via Clustal omega V1.2.4 (https://www.ebi.ac.uk/jdispatcher/msa/clustalo). Consensus sequences were then generated for each RV protein via the EMBOSS Cons CLI V6.6.0 tool (https://www.ebi.ac.uk/Tools/). These target sequences were used as a reference for prediction and epitope mapping analyses, and their antigenicity was evaluated via the VaxiJen v2.0 server (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html, accessed on 17 February 2024). This server system predicts the antigenicity of a protein via alignment-independent prediction, with a cutoff score of > 0.45. VaxiJen server shows the antigenic nature of a protein-based on the physicochemical properties of amino acids. This server’s approach for antigen prediction is based on auto cross covariance (ACC) transformation of protein sequences into uniform vectors with amino acid features40.
Prediction of linear B-cell epitopes and analysis
Linear B-cell epitopes of the subunit vaccine candidates were separately predicted via several combined servers to increase the accuracy of epitope prediction.
The BCPREDS server (http://ailab.cs.iastate.edu/bcpreds/) uses three developed methods, AAP, BCPred and FBCPred, to predict B cell epitopes. BCpred employs a subsequent kernel-based support vector machine (SVM) classifier with an accuracy AUC value of 0.758 and a cut-off score of 0.8 to predict linear B-cell epitopes41. However, ABCpred (http://crdd.osdd.net/raghava/abcpred/) applies the recurrent artificial neural network ANN method to classify epitopes and nonepitopes42. The threshold value was defined as 0.8, with a corresponding window length of 16mer.
BepiPred2.0 (http://services.healthtech.dtu.dk/service.php?BepiPred-2.0) uses a random forest algorithm trained on the epitopes annotated from antibody-antigen protein structures. BepiPred-2.0 utilizes sequence-based epitope prediction both on epitope data derived from 3D structures, and on a large collection of linear epitopes downloaded from the IEDB database43. BcePred (http://crdd.osdd.net/raghava/bcepred/), which accepts the primary sequence of the protein as an input, predicts epitopes on the basis of physico-chemical properties of amino acids such as (hydrophilicity > 1.9, flexibility > 2.0, accessibility polarity > 1.8), and exposed surface to forecast B cell epitopes44.
Epitope 1D (http://biosig.lab.uq.edu.au/epitope1d/) uses a random forest algorithm while employing a graph-based signature representation of protein sequences and integrating organism ontology information45.
Default parameters were set in all tools, and the partial or complete overlapping sequences of the top predicted epitopes from at least three software programs were chosen as B-cell dominant epitopes. The antigenicity, allergenicity, toxicity, and immunogenicity of B-cell epitopes were evaluated using computational tools. Antigenicity was predicted using VaxiJen v2.0, while AllerTop v2.0 (http://www.ddg-pharmfac.net/AllerTOP) assessed allergenicity by converting protein sequences into uniform amino acid composition (ACC) vectors of fixed length, applying a threshold of −0.4 for specificity46.
ToxinPred (http://webs.iiitd.edu.in/raghava/toxinpred) determined toxicity by analyzing physicochemical properties against a dataset of 1,805 toxic peptides47. Finally, IEDB’s immunogenicity tool (http://tools.iedb.org/immunogenicity/) predicted immunogenic potential based on physicochemical features, including side chain composition and amino acid positioning48.
Cytotoxic CD8 + T lymphocyte (CTL) epitope prediction
The NetCTL1.2 server (http://www.cbs.dtu.dk/services/NetCTL/)49 was utilized to predict 9mer-length CTL epitopes for the G, P, N, and M RV proteins. The threshold value of 0.75 was set against class 1, 12 MHC supertypes (A2, A3, A24, A26, B7, B8, B27, B39, B44, B58, and B62). The server employed default values to predict these CTL epitopes, which involved the combination of three features: TAP transport efficiency = 0.15, C-terminal cleavage = 0.05, and MHC-I binding peptides.
Artificial neural networks (ANNs) were used to estimate major histocompatibility complex (MHC) class I binding and proteasomal C-terminal cleavage, whereas a weight matrix was used to calculate the efficacy of the TAP transporter. The output derived from the neural network predicting MHC/peptide binding is a log-transformed value related to the IC50 values in nM units. Epitopes that confirmed interactions with ≥ 4 human leukocyte antigen (HLA) I supertypes were subjected to downstream investigation.
The recognized epitopes were further detected via immune epitope consensus (IEDB) (http://tools.iedb.org/mhci/)50 method. Peptides with percentile ranks < 0.2 and IC50 values < 50 nM were considered strong binders against a complete HLA reference set.
The predicted CTL epitopes were then screened for antigenicity, allergenicity, and immunogenicity. Only epitopes that showed a positive value for immunogenicity were retained for the next stage of evaluation.
Helper T lymphocyte epitope prediction
NetMHCIIpan–4.0 (http://services.healthtech.dtu.dk/service.php? NetMHCIIpan-4.0) was used to predict 15-mer-long peptide binding to any MHC-II molecule of the known sequence via ANNs51. The most frequent HLA reference set molecules in all populations, DRB1∗01:01, DRB1∗03:01, DRB1∗04:01, DRB1∗07:01, DRB1∗08:01, DRB1∗11:01, DRB1∗13:01 and DRB1∗15:01, were chosen for epitope hunting. The predicted epitopes were ranked via the IEDB SMM approach, applying a cutoff of 500 nM for the IC50. The shortlisted epitopes that interacted with ≥ 3 supertype alleles were assessed on the basis of their antigenicity, nonallergenicity, nontoxicity, and immunogenicity.
Peptide-protein Docking between the CTL and HTL epitopes and MHC alleles
The PepDock (https://galaxy.seoklab.org/cgi-bin/submit.cgi?type=PEPDOCK) available on GalaxyWEB (https://galaxy.seoklab.org/)52 was used for molecular docking of CTL and HTL epitopes against their relevant MHC alleles to determine their likelihood of being presented on the cell surface for T-cell recognition. the PDB structures of HLA alleles such as HLA-A*01:01 (PDB: 6MPP), HLA-A*02:01 (PDB: 2GIT), HLA-A*02:06 (PDB: 3OXR), HLA-A*08:01 (PDB: 7NVI), HLA-B*15:01 (PDB: 8ELG), HLA-B*35:01 (PDB: 4LNR), HLA-DRB1*01:01 (PDB: 5V4N), HLA-DRB1*04:01 (PDB: 5LAX), HLA-DRB1*11:01 (PDB: 6CPL), and HLA-DRB1*15:01 (PDB: 5ELH) were retrieved from RCSB Protein Data Bank server (https://www.rcsb.org/). Any additional peptides, ions, water, and ligands were eliminated using ChimeraX v1.8. Finally, the docking score between epitope and alleles was computed using the PROtein BinDIng enerGY prediction PRODIGY tool of the HADDOCK server (https://wenmr.science.uu.nl/prodigy/)53. This server incorporates the following computational parameters: the binding affinity (ΔG) and the dissociation constant (Kd) of the complex at room temperature.
Population coverage analysis
T-cell epitopes with their respective HLA alleles were subjected to the IEDB population coverage analysis tool (http://tools.iedb.org/population/) for population coverage calculations across the global population. This tool is utilized to estimate the proportion of individuals likely to respond to a defined epitope set, based on HLA genotypic frequencies and data related to MHC binding or T cell restrictions54.
Nonhuman homologous epitope selection
To avoid the occurrence of autoimmunity and screen nonhuman homologous epitopes, BLASTp (https://blast.ncbi.nlm.nih.gov/Blast.cgiPAGE=Proteins) was used. Epitopes without hits below the E-value inclusion threshold of 0.005 and identities less than 35% were chosen as nonhost/nonhomologous proteins.
Multi-epitope mRNA vaccine construct design
The final selected LBL, HTL, CTL and epitopes were joined with KK, GPGPG and AAY linkers to design vaccine constructs. To create an efficient vaccine and enhance immune responses, the 50 S ribosomal protein L7/L12 adjuvant sequence, a potent TLR4 agonist, was sequentially added to the shortlisted epitopes. Additionally, a tissue Plasminogen Activator (tPA) secretory signal sequence (tPA) was added at the N-terminus to help transport the vaccine out of the cell, and an MHC I-targeting domain (MITD) at the C-terminus was used to direct CTL epitopes to the MHC-I compartment of the endoplasmic reticulum.
A 5′ cap, a Kozak sequence, a poly(A) tail of 120–150 bases, and the 3′ and 5′ untranslated regions (UTRs), which are essential for translation and stability, were also added to the mRNA vaccine construct.
Prediction of the antigenicity, immunogenicity, allergenicity and toxicity of the vaccine construct
The antigenicity, allergenicity, toxicity, and immunogenicity of the designed mRNA vaccine candidates were estimated via online servers as previously described.
Physicochemical characteristics of the vaccine construct
The physical and chemical parameters, including the theoretical isoelectric point, molecular weight, total number of positive and negative residues, extinction coefficient, instability index, aliphatic index, and grand average hydropathy of the vaccine construct, were computed via the ExPASy ProtParam server (http://us.expasy.org/tools/protparam.htm)55.
Prediction of protein solubility
The SOLpro server (http://scratch.proteomics.ics.uci.edu) was used to predict vaccine construct solubility. SOLpro predicts the propensity of a protein to be soluble upon overexpression in Escherichia coli via a two-stage SVM architecture that is based on multiple representations of the primary sequence56.
Prediction of the secondary structure
The secondary structure of the overall vaccine construct was predicted via PSIPRED (version 4.0) (http://bioinf.cs.ucl.ac.uk/psipred/). PSIPRED is a simple and accurate secondary structure prediction method that incorporates two feed-forward neural networks that perform an analysis on outputs obtained from PSI-BLAST (position-specific iterative-BLAST)57. The solvent accessibility (ACC) and disordered regions (DISO) are predicted by the RaptorX Property server58.
3D structure modeling and refinement
The tertiary structure of the vaccine subunit was predicted via the transform restrained trRosetta server (https://yanglab.nankai.edu.cn/trRosetta/)59. This web-based platform combines deep learning and Rosetta to achieve fast and accurate protein structure prediction. By employing constrained Rosetta and direct energy minimizations, along with the implementation of a deep neural network, the protein structure was successfully built while simultaneously determining the locations and orientations of the residues.
To improve the quality of the 3D model obtained from the Rosetta server, the 3Drefine web server (http://sysbio.rnet.missouri.edu/3Drefine/) was used60. The refining protocol is based on a two-step process: optimization of the hydrogen-bonding network and atomic-level energy minimization. Hence, the quality of the improved model was assessed via the GDT-HA, GTD-TS, MolProbity, clash and RMSD scores.
Vaccine 3D structure validation
Refined 3D structure models were validated for possible structure errors through Ramachan- dran Plots and Z scores obtained from the PROCHECK (http://saves.mbi.ucla.edu/)61 and ProSA (http://prosa.services.came.sbg.ac.at/prosa.php) servers62. ERRAT was used to examine nonbonded atomic interactions63.
Discontinuous (Conformational) B-cell epitope screening
The refined 3D structure of the peptide vaccine model was submitted to the IEDB Ellipro tool (http://tools.iedb.org/ellipro/) to identify the conformational B-cell epitopes. ElliPro provides a score for each output epitope on the basis of neighboring residues, estimation of shape (ellipsoid), and protrusion index parameters64.
Molecular Docking between the multiepitope subunit vaccine and toll-like receptor (TLR4)
The docking was executed with the ClusPro 2.0 protein‒protein docking online server (http://cluspro.bu.edu/)65 to determine the binding affinity of the vaccine construct with target immune cell receptors. The TLR4 (PDB 3FXI) 3D structure was retrieved from the RCSB Protein Data Bank (RCSB PDB) (https://www.rcsb.org/). This server, by running rigid body docking with the PIPER docking program, can generate different models on the basis of different scoring schemes. The structures of the peptides were visualized via PyMOL66. Docked complexes with the largest cluster size and high binding affinities (lowest energy structure) were further subjected to molecular dynamics simulations.
Finally, PDBsum67 was used to perform a detailed analysis of protein–protein interactions in the complexes. The binding strength of the docked complexes was explored with PRODIGY.
Vaccine-TLR4 complex molecular dynamic (MD) simulation
To evaluate the dynamic behavior and stability of the docked TLR4–vaccine complex, an explicit all-atom Molecular Dynamics Simulation (MDS) was performed using GROMACS 2023.168, compiled with CUDA support to leverage GPU acceleration. The simulation was executed on the high-performance computing cluster at the Institut Pasteur de Tunis (termed omics), equipped with NVIDIA GPUs and 200 CPU cores, enabling efficient parallel computation69.
The initial 3D structure of the complex was first processed using the pdb2gmx module with the CHARMM36 force field, selecting appropriate protonation states for physiological conditions and employing the SPC/E water model for solvation. The processed structure was centered in a cubic simulation box with a 1.0 nm margin from the edges and subsequently solvated using the spc216 predefined water configuration. To neutralize the system and mimic physiological ionic conditions, sodium (Na+) and chloride (Cl-) ions were added by replacing solvent molecules using the genion tool. Energy minimization was performed using the steepest descent algorithm to remove steric clashes and relax the system. The minimized structure then underwent a two-step equilibration process: first under constant volume and temperature (NVT ensemble), and then under constant pressure and temperature (NPT ensemble), each for 100 ps while restraining heavy atoms. The production molecular dynamics simulation was executed for a total of 200 ns in the NPT ensemble, with a time step of 2 fs, using 50 MPI threads and 4 OpenMP threads per MPI process, and leveraging GPU acceleration via CUDA. Electrostatics were treated using the Particle Mesh Ewald (PME) method, and all bond lengths were constrained using the LINCS algorithm to ensure numerical stability. Trajectory files were saved every 10 ps. Post-simulation analyses were performed to evaluate system stability and conformational changes: the root mean square deviation (RMSD) of the backbone atoms (Cα) was computed to assess structural drift over time; the average temperature was extracted from the energy profile to ensure thermal stability; hydrogen bonds between main-chain atoms were quantified across the trajectory to evaluate intra- and intermolecular cohesion; and representative structural snapshots were extracted every 10 ns to visualize the time-resolved evolution of the complex70. This comprehensive protocol ensured accurate modeling of the protein–vaccine interactions in explicit solvent conditions over an extended timescale.
Normal (NMA) of the TLR4-vaccine complex
Normal (NMA) of the TLR4-vaccine structure was performed the iMODS server (http://imods.chaconlab.org/)71 to study the structural dynamics and flexibility of the vaccine-TLR4 complex through structural deformation analysis, the eigenvalue of the interacting residues (energy required to deform the structure), the B-factor value scores (mobility profiles), the covariance map among individual residues and the elastic network model. The basic interface was used with the CA option for thecoarse grain model atomic representations.
In silico immune simulation
The C-IMMSIM server (https://kraken.iac.rm.cnr.it/C-IMMSIM/)72 (accessed on 17 August 2024) was used to perform an in silico immune simulation that elucidates its potential to trigger the host immune system toward the designed peptide vaccine. The server simultaneously simulates three compartments of immune cells in mammals—the bone marrow, the thymus, and a lymph node—and makes use of position-specific scoring matrices (PSSMs) derived from machine learning techniques for predicting immune interactions.
At intervals of four weeks, three doses of injection, each containing 1000 antigen proteins and no LPS, were administered, and all simulation parameters were set at default values, with time steps set at 1, 84, and 168 simulation steps. A single time step is equivalent to 8 h of daily life, and the first dose was given at time = 0. The input parameters for the immune simulations included a standard volume of 10 µl, 12,345 random seeds and 1050 simulation steps.
Codon optimization and secondary structure prediction of the mRNA vaccine
The optimization of the vaccine construct codons was performed via GenScript (http://www.genscript.com/tools/gensmart-codon-optimization) (accessed on 17 August 2024) to optimize the codon usage of the vaccine in the human host. The GC content of the optimized DNA sequence, codon adaptive index (CAI) and codon frequency distribution (CFD) were determined via the GenScript Rare Codon Analysis Tool (http://www.genscript.com/tools/rare-codon-analysis). The ideal percentage range of the GC content for full expression in the host is between 30% and 70%. The CAI is a measure of synonymous codon usage bias for a target sequence. It ranges from 0 to 1, with values higher than the cutoff (> 0.8), indicating that genes that are more likely to be highly expressed. The occurrence of unusual tandem codons is denoted by the codon frequency distribution (CFD). The optimized DNA sequence was then converted to an RNA sequence via a translation tool, and the prediction of the secondary structure was carried out via the RNAFold web server (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi)73 for thermodynamic analysis and minimal free energy score.
Results
A flow chart showing a graphical outline of the approaches used for designing mRNA-based vaccines against RV is depicted in Fig. 1.
Proteome screening and protein prioritization
The available sequences of the G, P, NP, M and L proteins from RV were retrieved from NCBI and subjected to multiple sequence alignment. EMBOSS software generated a consensus sequence for each protein, and the conservancy was analyzed with the respective reference sequences from the Pasteur strain.
The results revealed that four (G, P, N and M) of the five proteins were antigenic (Table 1).
These antigens were selected as vaccine targets for the in silico analysis of mRNA vaccine development against RV.
Linear Bcell epitope prediction
ABCpred, BCEpred, BCPREDS, Bepipred and epitope 1D servers were used for linear B-cell prediction (Fig. 2, Table S1). The predicted epitopes that overlapped at least three high-score servers and satisfied the criteria of antigenicity, toxicity, allergenicity and immunogenicity were retained for the final vaccine construction. A total of 13 epitopes were predicted for all the tested protein components (Table 2).
Workflow and tools used in this study for the in silico design of rabies mRNA-based multiepitope vaccines.
Cytotoxic T lymphocyte (CTL) epitope prediction
Among the four chosen RV proteins, a total of 13 CTL epitopes with a length of 9 aa exhibited strong binding to at least four MHC class I supertypes (Table S2). These epitopes were then assessed for antigenicity, immunogenicity, nonallergenicity, and nontoxicity. Ultimately, four epitopes (one from N, two from G, one from M) were predicted to be antigenic, immunogenic, nonallergenic, nontoxic, and CTL epitopes (Table 3). The IEDB server was subsequently confirmed through a combined approach in which epitopes were selected as potent CTL epitopes.
Helper T lymphocyte (HTL) epitope prediction
According to the predictions made via NetMHCIIpan version 4.0, a total of 14 MHC-II epitopes (15-mers) in the antigenic proteins were identified to have binding affinity with at least three unique HLA DRB alleles (Table S3). After thorough assessment for strong binding affinity, antigenicity, nonallergenicity, and nontoxicity, two epitopes for the G protein, as well as a single epitope for the NP and M proteins, were selected for the final vaccine design (Table 4). The majority of these HLA epitopes demonstrated the ability to induce cytokines such as IFN-γ and IL-10.
Phosphoprotein
B-cell epitope mapping results from different server tools. The black squares represent epitopes predicted by BCPREDS, the red squares represent epitopes predicted by ABCpred, the green frames denote the epitopes predicted by BCEpred, the dotted squares represent epitopes predicted by Bepipred, and the purple squares denote epitopes predicted by 1Depitope. The black lines with numbers on both ends represent the length of the consensus-generated sequences used in this study.
Molecular Docking of CTL and HTL epitopes with HLA alleles
The GalaxyPepDock server was utilized to perform molecular docking of CTL and HTL epitopes against their relevant MHC alleles, which revealed a strong binding interaction and thermodynamic stability of these epitopes with their respective HLA alleles, as indicated by the negative Gibbs free energy values (ΔG <−7.6 Kcal mol−1, average =−17.5) and low dissociation constant values. Figure 3 showed that the epitopes are docked within the binding groove of their corresponding HLA alleles.
Molecular docking of N protein epitopes (NPE), M protein epitopes (NPE) and G protei epitopes (GPE) with respective HLA alleles. For each protein, the most significant representative epitopes were identified, and their binding with the alleles exhibiting the highest affinity was presented. The docking of protein-peptide was carried out using the GalaxyPepDock server. The free energy (ΔG) values for each binding demonstrate the affinity between the epitopes and alleles, determined through the PRODIGY server. Ribbon structures illustrate HLA alleles, while stick structures depict the corresponding epitopes.
Population coverage analysis
The worldwide combined population coverage of the shortlisted MHC-I and MHC-II epitopes reached 96.01% (Fig. 3), with a pc90 average of 1.5 and an average hit per HLA antigen of 4.08. Europe and North America had the highest rates at 98.42% and 97.62%, respectively, followed by North Africa at 92.54% (Fig. 3). Among the African regions identified as most affected by rabies, Guinea Bissau, Cap Verde, Morocco, Tunisia, and Zimbabwe had over 90% population coverage for the chosen epitopes. Additionally, population coverage rates of over 80% were observed in Sudan, Cameroon, Senegal, South Africa, Uganda, and Zambia (Fig. 4 and Fig. S1).
IEDB global population coverage analysis of the selected T lymphocyte epitopes used in the vaccine construct.
.
Construction of a multiepitope mRNA vaccine candidate sequence
The final construct of our mRNA vaccine was designed from the N- to C-terminus, as shown in Fig. 5.
Schematic representation of the designed mRNA-based multi-epitope construct.CTL: cytotoxic T cell, HTL: helper T cell, HTL: helper T cell, tPa: tissue plasminogen activator, MITD: MCH I- trafficking domain.
The coding region of the vaccine consists of three sets of antigenic, nontoxic, and nonallergenic CTL, HTL and B-cell shortlisted epitopes linked sequentially with AAY, GPGPG, and KK linkers, respectively. To increase the immunological response, the adjuvant 50 ribosomal protein L7/L12 was added via the GPGPG linker. Crucial mRNA elements, such as the 5′ m7G cap sequence, the Kozak sequence including the AUG codon, the globin 5′ and 3′ untranslated regions (UTRs) flanking the ORF of mRNA, and the poly (A) tail of 120–150 bases, have been integrated into the vaccine to improve translation efficiency. Moreover, to facilitate the release of translated epitopes, a tissue plasminogen activator (tPA) secretory signal sequence (tPA) was added at the N-terminal extremity through the GPGPG linker. Additionally, an MHC I-targeting domain (MITD) was attached via an EAAK linker at the C-terminus to guide CTL epitopes toward the MHC-I compartment of the endoplasmic reticulum.
Evaluation of antigenicity, allergenicity, toxicity, and immunogenicity
Upon further analysis, the chimeric construct was observed to be antigenic, with high probability scores of 0.6006 and 0.869, as predicted by the VaxiJen v2.0 and ANTIGENpro servers, respectively. The multi-epitope vaccine was also detected as a nonallergen, nontoxic, soluble (0.828496) and immunogen (0.29221).
Secondary structure prediction
The secondary structure of the final vaccine was estimated to contain 32% alpha helices, 14% extended strands and 53% coils, as predicted by the method of Garnier-Robson. The solvent accessibility (ACC) predicted by the RaptorX Property server revealed that 56% of the residues are predicted to be exposed, 21% are medium exposed, and 21% are buried exposed. Only 11% of the residues are predicted to be located in disordered regions (DISO). The secondary structure of the RV is presented in Fig. 6.
Physicochemical property analysis
The results of the physicochemical property analysis are presented in Table 5. The designed construct is approximately 52.10 kDa (less than 110 kDa) and consists of 72 positively and 70 negatively charged residues. The isoelectric point value is 8.13 (basic). Furthermore, the SOLpro server, which predicts solubility when overexpressed, showed good predictive value (0.88). The instability index is computed to be 32.76, which classifies the protein as stable. The aliphatic index, GRAVY (grand average of hydropathicity) and extinction coefficient of the construct are 81.98, − 0.430 and 58,455 M − 1 cm − 1 at 280 nm, respectively.
Tertiary structure modeling and refinement
Five tertiary structures of the designed multi-epitope vaccine were successfully modeled via the trRosetta server. The estimated TM score provides insight into the structural similarity between two structures, with values ranging from 0 to 1. A TM score above − 0.5 signifies correct topology, and a score below − 0.17 implies random similarity. The selected model exhibited an estimated TM score of −0.331.
The refinement result is reported as a 3Drefine score. The retained model showed great overall quality, attaining a global distance test-high accuracy (GDT-HA) of 0.9686 (a score ranging between [0, 1] indicates a similar score with respect to the initial model), an RMSD of 0.366 Å, which is the root-mean-square deviation score of Cα atoms (Cα-RMSD) of the refined model with respect to the initial model, a low molProbity of 1.277, a clash score of 5.2, and a poor rotamer score of 0.0. Figure 7 A and B illustrate the retained refined model and its sequence.
Secondary structure prediction of the vaccine construct via the position-specific iterated prediction (PSIPRED) server.
Vaccine 3D structure validation
The refined tertiary structure was validated via ProSA-web, RAMPAGE, and ERRAT severs. The refined model revealed that 93.2%, 5.4%, 0.7% and 0.7% of the residues were located in the most favored, favored, allowed, generously allowed and outlier regions, respectively (Fig. 7 C). In addition, ProSA has a Z score of −5.39 (Fig. 7D), which lies inside the range of acceptable scores. The local model quality analyzed by ERRAT was 94.72 (Fig. 7E). All these results indicated that the predicted structure of the construct was reasonable and reliable.
Multiple-epitope vaccine sequence structure prediction and validation. (A) Sequential arrangements of the multiple epitope vaccine (adjuvant, EAAK, AAY, GPGPG, KK, and epitopes, labeled with blue, orange, red, brown, green, and black, respectively). (B) Tertiary structure of the refined construct indicating the α-helix, β-strand, and random coil. (C) The ProsA Z score (− 5.39) indicates that the model quality is closer to the protein structure determined by X-ray crystallography; hence, the positively refined (D) ERRAT plot has a score of 94.72. (E) Ramachandran plot of the refined model showing 93.2%, 5.4%, 0.7% and 0.7% of residues in the most favored (red) and favored (yellow), generously allowed (pale yellow), and disallowed regions (white), respectively.
Screening for conformational B-cell epitopes
ElliPro was used to predict the validated 3D structure of discontinuous (conformational) B-cell epitopes. In total, eight BCL conformational epitopes involved 242 residues, with scores ranging from 0.529 to 0.871 (Fig. 8 and Table S4).
Conformational B-cell epitopes predicted in vaccine candidate by the ElliPro tool of IEDB server. The yellow spheres represent the conformational B-cell epitopes. (A) 43 residues (AA 1–43) with a score of 0.871. (B) 69 residues (AA 382–450) with a score of 0.845. (C) 16 residues (AA 251–266) with a score of 0.669. (D) 14 residues (178–191) with a score of 0.663. (E) 52 residues (AA 269–320) with a score of 0.657. (F) 6 residues (AA 376–381) with a score of 0.638. (G) 36 residues (AA 72–129) with a score of 0.583. (H) 5 residues (AA192–196) with a score of 0.529.
Molecular Docking of the vaccine protein and TLR4
The molecular docking between the vaccine construct and TLR4 was accomplished via the ClusPro 2.0 server, which generated 10 docked complexes ranked by energy and cluster members (Table S5). Cluster 7, which had the lowest energy score of −1600.5 kcal/mol with 21 cluster members, was chosen as the best complex (Fig. 9A).
PDBsum interaction analysis revealed 23 hydrogen bonds between chain B of TLR4 and the vaccine (chain A) and 2 salt bridges, while 302 were nonbonded contacts (Fig. 9B).
The amino acids involved in hydrogen bonds, as well as their lengths, are listed in Table S6.
Advanced analysis conducted with the PRODIGY web server to evaluate Gibbs energy (ΔG) and the dissociation constant of the docked proteins indicated that the TLR4 vaccine complex was energetically viable, exhibiting negative ΔG values of −18.4 Kcal mol-1 and a lower dissociation constant of 3.4e-14.
Molecular docking of the vaccine construct with TLR4. (A) Docked complex of the vaccine (chain A) and TLR4 (chain B). The vaccine (ligand) is shown in cyan, and chain B from TLR4 (receptor) is shown in green. (B) Amino acid interactions between the docked vaccine (chain A) and the TLR4 receptor (chain B) obtained from PDBsum analysis. Residues involved in hydrophobic interactions are shown as red dotted lines, and those involved in hydrogen bonding are shown as green dotted lines. (C) The interacting residues of the vaccine construct (chain A) and the receptor (chain B) visualized by Pymol.
Molecular dynamic simulation
To evaluate the overall structural stability of the TLR4–vaccine complex during the molecular dynamics simulation, the Root Mean Square Deviation (RMSD) of the backbone atoms was calculated over a period of 200 nanoseconds (ns). As shown in Fig. 10 A and Supplementary Fig.S2 the RMSD profile exhibited the following characteristics. Initial increase (0–50 ns): The RMSD rose sharply from ~ 0.3 nm to ~ 2.5 nm, indicating a rapid relaxation of initial structural constraints imposed during docking. This phase typically corresponds to equilibration and early conformational adjustment. Intermediate fluctuations (50–150 ns): The RMSD values fluctuated between 1.5 nm and 3.5 nm, suggesting conformational flexibility of the complex. This dynamic region may involve loop motions or domain rearrangements, common in immune receptors such as TLR4. Stabilization phase (after 150 ns): From around 150 ns onwards, the RMSD values plateaued around 3.2–3.4 nm, indicating that the system reached a relatively stable conformational state or local energy minimum. This plateau reflects a convergence of the simulation and supports the reliability of the sampled conformations for further structural and energetic analyses. The high RMSD amplitude (> 2.0 nm) suggests that the complex exhibits significant structural plasticity, which may be biologically relevant for receptor recognition or signal transduction. Despite these fluctuations, the system demonstrates progressive stabilization, a hallmark of successful long-timescale molecular dynamics simulations.
This RMSD trend reinforces the quality of the docking pose and confirms that the TLR4–vaccine interaction remains structurally stable over 200 ns, providing a reliable foundation for downstream analyses including binding free energy calculations, hydrogen bond stability, and essential dynamics.
To further validate the structural stability of the TLR4–vaccine complex observed in the global RMSD analysis, we performed a focused evaluation of the interaction interface between the receptor and the vaccine construct. As shown in Fig. 10A and B, the Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) were calculated specifically for the interfacial residues involved in binding. The interface RMSD remained consistently low (~ 0.12–0.15 nm) over the 200 ns simulation, indicating stable interactions with minimal deviation. RMSF profile reveals moderate flexibility across the protein, with values ranging from 0.04 to 0.14 nm. Higher fluctuations are observed around residues 190 (~ 0.14 nm), 300 (~ 0.135 nm), and 340 (~ 0.145 nm), corresponding to loop or terminal regions. Most residues remain within 0.04–0.11 nm, suggesting limited local flexibility and stable core regions.
Molecular dynamics simulation analyses of the TLR4–vaccine complex with 200 ns. This figure illustrates the dynamic and structural behavior of the TLR4–vaccine complex: (A) RMSD shows an initial rise to ~ 0.15 nm, followed by moderate fluctuations and stabilization between 0.10–0.125 nm, indicating structural convergence. (B) RMSF highlights flexible regions near residues 190, 300, and 340, while most residues remain stable (0.04–0.11 nm), reflecting limited local mobility. (C) Radius of Gyration (Rg) decreases from ~ 6.0 nm to below 4.25 nm, suggesting progressive compaction of the complex. (D) SASA stabilizes between 614–634 nm² after an initial peak, indicating consistent solvent exposure and structural folding. (E) Intermolecular Hydrogen Bonds (F) PCA-Based Conformational Landscape of the TLR4–Vaccine complex derived from a 200 ns molecular dynamics simulation of a TLR4–vaccine complex, mapping the system’s conformational states onto the first two principal components (PC1 and PC2), which represent the dominant motions and energetically favorable conformational states.
In addition, Fig. 10C and D illustrate the Radius of Gyration (Rg) and the Solvent Accessible Surface Area (SASA) of the entire complex, respectively. Radius of Gyration (Rg) begins at ~ 6.0 nm and decreases steadily to ~ 5.2 nm at 25,000 ps, with transient fluctuations observed throughout the simulation. A major drop to ~ 4.25 nm at 100,000 ps indicates system compaction. The Rg values fluctuate slightly between 4.25 and 4.5 nm from 100,000 to 150,000 ps, and drop below 4.25 nm by the end of the simulation, suggesting increased compactness.
The SASA exhibited a steady decline from ~ 637 to ~ 620 nm², indicating reduced exposure to solvent, likely due to stabilization of the tertiary structure and maintenance of the binding interface. Together, these structural descriptors provide strong evidence for the conformational stability and compactness of the TLR4–vaccine complex, both globally and at the molecular interface, supporting its suitability for downstream analyses such as binding free energy and essential dynamics. The PCA heatmap complements the free energy landscape, highlighting how the TLR4–vaccine complex progressively stabilizes into a narrow ensemble of low-energy conformations. The dominance of PC1 (explaining maximal variance) suggests concerted structural rearrangements, possibly linked to functional motions (e.g., binding-site adjustments). Principal component analysis (PCA) was conducted to capture dominant collective motions within the TLR4–vaccine complex throughout the MD simulation. Figure 10F shows the projection of the simulation trajectory onto PC1 and PC2, the first two principal components representing the dominant conformational changes. The plot displays the trajectory progression over time, with a color gradient from dark red (start of simulation) to gray (end of simulation), illustrating the continuous yet confined exploration of conformational space. Multiple regions are revisited along the trajectory, suggesting the presence of stable and semi-stable conformations with persistent intermolecular interactions at the protein–vaccine interface. This observation is consistent with the pattern revealed by the Intermolecular Hydrogen Bonds analysis, which showed a consistent number of hydrogen bonds throughout the simulation (Fig. 10E). The plot represents simulation snapshots indicating the presence of stable and frequently sampled conformational states of the complex. Conversely, higher-energy regions correspond to less stable, transient conformations, illustrating that the TLR4–vaccine complex explores a diverse yet confined conformational space, consistent with a stable protein-ligand system and corroborating the findings from other analyses such as RMSD, Rg, and hydrogen bond analysis, which collectively affirm the structural robustness and biologically favorable interaction of the vaccine candidate with the TLR4 receptor. The presence of dense red clusters reflects energetically favorable and frequently visited conformational basins, supporting the results from the FEL analysis. These stable states suggest a compact, robust structural behavior of the complex, corroborating findings from RMSD, Rg, SASA, and hydrogen bond analyses. Overall, this dynamic yet confined motion underscores the biologically relevant and persistent interaction between the vaccine construct and the TLR4 receptor reinforcing the structural integrity of the vaccine construct. In conclusion, these results demonstrate that the TLR4–vaccine complex maintains strong structural stability, consistent intermolecular interactions, and a compact folded conformation during the 200 ns MD simulation. These findings validate the vaccine candidate’s potential for productive receptor engagement under physiological conditions.
NMA of the vaccineTLR4 complex
To further analyze the vaccine-TLR4 complex, NMA was performed via the iMODS server. iMODS uses a normal mode analysis approach to explore the collective motions of macromolecules using normal mode analysis (NMA). The results are shown in Fig. 11.
Normal mode analysis (NMA) of the docked complex (Fig. 11A). The deformability of the complex is based on the specific inaccuracy of each residue, represented by chain hinges. The peaks in the graphs are correlated with the deformed regions of the protein (Fig. 11B). The B-factor value scores corroborated the NMA results (Fig. 11C). The eigenvalue associated with each normal mode represents the motion stiffness. Its value is directly related to the energy required to deform the structure. The lower the eigenvalue is, the easier the deformation. As shown in Fig. 11D, the eigenvalue of the complex was 1.81729 e–08, which signifies that the complex is stable. Fig. 11E illustrates the variance associated with each normal mode. The colored bars show the individual an cumulative variances. The variance is inversely related to the eigenvalue. The covariance analysis is depicted in Fig. 11F to describe the correlated, uncorrelated (white), or anti-correlated (blue) atomic movements in the dynamical regions of the complex molecule. The elastic network model defines the pairs of atoms connected by springs. Each dot in the graph represents (Fig. 11G) one spring between the corresponding pair of atoms. The dots are colored according to their stiffness; the darker grays indicate stiffer springs. Taken together, these results predict the stability of the RAB-VAC TLR4 complex.
Normal mode analysis of the docked complex. (A) Normal mode analysis mobility (NMA), (B) deformability of complex, (C) B-factor revealing the association between the NMA and PDB fields in the complex, (D) Eigenvalue representing the energy required to deform the structure, (E) variance inversely related to the eigenvalue cumulative variance shown in green and variance in red, (F) covariance map with correlated (red), uncorrelated (white) or anticorrelated (blue) motions, and (G) elastic network where darker gray areas indicate greater stiffness.
In Silico immune response
The effective immune response profile of the designed vaccine was studied using the C-ImmSim server. The immune simulation analysis predicts a strong humoral immune response. The primary immune response showed high titers of immunoglobulin IgM. Following subsequent exposures, secondary and tertiary immune responses were more active, as evidenced by increased levels of IgG1 + IgG2 and IgM compared to lower levels of IgG + IgM antibodies (Fig. 12A). High titers of (IgG, IgM, IgG1, and IgG1 + IgG2) antibodies were recognized, resulting in a decrease in antigen concentration. Furthermore, several simulated B cell isotypes and memory B cell production were detected, which suggests a productive long- lasting immune reaction induced by the vaccine construct and also a greater antigen clearance (Fig. 12B and C). The TH (T helper) and TC (T cytotoxic) cell populations increased significantly (reaching 11000 and 900 cells/mm3 around 50 days after exposure, respectively) and remarkably responded with the corresponding memory development (Fig. 12D and E).
The Th0 type immune reaction had a lower number (less than 20000 cells/mm3) than the Th1 type immune reaction peaking to 120,000 cells/mm3 around day 50 (Fig. 12F). In fact, sustained significant levels of Th1 cells were achieved after a series of exposures (Fig. 12F).
Meanwhile, the activity of natural killer cells (NK cells) and dendritic cells (DCs) increased concomitantly after exposure together with an improved proliferation of macrophages (Fig. 12G, H and I), suggesting an efficient antigen presentation to the T cells. Remarkable induction of different cytokines such as of IFN-γ, TGF-β, IL-23, IL-10, and IL-12 was also evident over the injections (Fig. 12J). According to the immune simulation results, the designed chimeric vaccine is predicted to trigger a robust immune response and to ensure immunogenicity and stability.
In silico immune simulation to predict the immunological potential of the designed vaccine RDV-V chimeric peptide as predicted by the C-ImmSim server. (A) Increased level of immunoglobin antibodies with a decrease in antigen levels upon vaccine injection. (B) Increased levels of B-cell populations with a decrease in antigen levels upon vaccine injection. (C) Rising B-cell populations after repeated exposure to antigen. (D, E) Increase in the population of T-cytotoxic and T-helper cells upon repeated antigen exposure. (F) Switching between different classes of T helper cells. (G–I) The populations of dendritic cells, macrophages, and natural killer cells increased during the immunization period. (J) Increased concentrations of cytokines and interleukins after repeated antigen exposure. The inset plot shows the danger signal together with the leukocyte growth factor IL-2.
Codon optimization and secondary structure prediction of the mRNA vaccine
The codon usage of the vaccine construct was optimized to ensure its expression efficiency in human hosts via an in silico program. The optimized codon sequence has a length of 1788 nucleotides and a CAI of 0.92 (Fig. 13A). The average GC content and the CFD value of the improved construct were 58.33% and 0%, respectively. Overall, these results clearly demonstrate the great potential for proficient expression in the target expression host.
The prediction of the secondary structure of mRNAs is essential, as it is a key factor influencing translation initiation, elongation, and mRNA biogenesis. The free energy related to the whole mRNA structure was obtained via the RNAfold online service. The minimum free energy of the secondary RNA structure was ΔG = −634.70 kcal/mol (Fig. 13B). The minimum free energy of the centroid secondary structure was − 533.34 kcal/mol (Fig. 13C). Mountain plots of the MFE structure, the thermodynamic ensemble of RNA structures and the centroid secondary structure reveals a significant overlap (Fig. 13D). The positional entropy for each position is also represented in Fig. 13E. These scores reflect the stability and efficiency of protein translation within hosts. Enhanced mRNA stability is directly related to increases in the expression rate.
Codon optimization and secondary structure of the designed mRNA vaccine construct via the RNAFold webserver. (A) CAI value. (B) Minimum free energy (MFE) structure of the mRNA vaccine construct. (C) Centroid free energy (CFE) of the mRNA vaccine construct. (D) Mountain plot representation of the MFE structure, the thermodynamic ensemble of RNA structures, and the centroid structure. The red, green, and blue lines represent the MFE structure, the thermodynamic ensemble of RNA structures, and the centroid structure, respectively. (E) Positional entropy plot for each position.
Discussion
Rabies is a life-threatening disease with no effective antiviral treatment and complex clinical symptoms. Hence, vaccines are the most effective measure for preventing rabies disease. Despite their effectiveness, current rabies vaccines do not confer lifelong protection, and booster doses are needed. New approaches to vaccination employing computational methods and bioinformatics tools for antigen mapping are advantageous and promising. Interestingly, a multi-epitope-based vaccine, by incorporating different immunogenic components in a single vaccine construct, vigorously stimulates both adaptive and innate markers of the immune system74 and thus can be broadly utilized both as prophylactic and therapeutic measures against a variety of pathogens and customized for specific populations or individuals.
To improve vaccine regimens and enhance the potency of the subunit vaccine, we designed a bioinformatics multi-epitope mRNA candidate vaccine that relies on in silico predictions and simulations. To our knowledge, this is the first study to design an mRNA multi-epitope RV vaccine by screening the most antigenic epitopes in the whole rabies virus proteome.
The key to developing a vaccine that can target a broad spectrum of viral strains lies in aligning their different protein sequences and identifying conserved sequences to select the most effective candidate epitopes. Consequently, consensus sequences were generated from multiple alignments of four proteins, G, P, NP and M, from different RV strains. The resulting sequences were used to predict linear B cell, CTL and HTL epitopes, which were then prioritized on the basis of their antigenicity, immunogenicity, nonallergenicity, nontoxicity and nonhomology, to increase their efficacy, safety and accuracy. Finally, B cell, HTL and CTL epitopes that were antigenic, nonhomologous, immunogenic, nontoxic, and nonallergenic and able to induce IFNγ and/or IL-4 or IL-10 production (for HTL epitopes) were shortlisted.
Molecular docking assessments determined the immunogenic potential of selected epitopes by measuring their binding affinity with their respective HLA alleles. The results demonstrated strong interactions, as reflected by consistently negative Gibbs free energy ΔG change and low dissociation constant values, which correlate with energetically favorable binding75. The strong binding affinity remains essential for proper epitope display on cell surfaces, thereby facilitating recognition by cytotoxic T lymphocytes (CTLs) and helper T lymphocytes (HTLs).
Population coverage analysis revealed that the combined CTL and HTL epitopes included in the vaccine construct encompassed 96.01% of the global population, with particularly high representation in regions historically affected by outbreaks.
The final peptide vaccine construct was designed by joining the B cell, HTL and CTL epitopes together via KK, GPGPG and AAY linkers, respectively. These linkers play a key role in the stability and folding of protein fusion as well as distinct functional domains to avoid overlapping and ensure effective antigen presentation76.
The presence of key components such as the 5‘cap, 5’ and 3’ UTR elements flanking the coding sequence, along with the translation initiation site “Kozak sequence” and poly-A tail with 120 bases, is crucial for the effective stability and translation of mRNA into protein77. Additionally, the tPA signal sequence, transcribed but not translated, is included prior to the coding sequence being identified by the translational machinery. Finally, the MHC class I trafficking signal domain (MITD) is also integrated to increase antigen presentation efficiency78. In addition, the 50 S ribosomal protein L7/L12 (RPLL) adjuvant was joined to the N-terminus of the construct via GPGPG linkers to increase the duration of the immune response. The RPLL from Mycobacterium tuberculosis Rv0652 is recognized by Toll-like receptor 4 (TLR4) to induce DC maturation and proinflammatory cytokine production and plays a crucial role in the innate immune response79. Coupled with the CTL, HTL and linear B-cell epitopes, the final rabies construct was 474 amino acids in length.
Analyses of the secondary structure revealed that the predicted peptide was composed predominantly of 53% coils, 32% α-helices, and 11% disordered residues. These regions have been identified as important ‘structural antigen’ types80.
High solubility of the vaccine in the aqueous environment of the host is essential for successful development. The GRAVY and solubility properties of these vaccines also indicate their potential for future use. The trRosetta-modeled structures were deemed satisfactory because of their rapid and accurate results for predicting the 3D structure of the vaccine constructs with the TM value of the best model; a higher TM value is indicative of the best modeled structure. This model was submitted to 3Drefine for further refinement. The top five final models were generated and ranked on the basis of the 3Drefine score, GDT-TS, GDTHA, RMSD, MolProbity, and RWPlus. The refined model with lower 3D refine score, RWplus, and MolProbity values, along with higher GDT-TS, GDT-HA, and RMSD values, indicates superior quality. The quality of the refined model was assigned by the Ramachandran plot, the ERRAT quality factor, and the ProSA z score. Ramachandran plot analysis revealed that 93.2%, 5.4%, and 0.7% of the residues were favored, allowed, and allowed regions, respectively, whereas only 0.7% were disallowed, suggesting minimal steric clashes between the side chain atoms and main chain atoms81. Similarly, the Z score and ERRAT quality factor were analyzed to identify possible errors in the structure. A negative Z score of −5.65 suggests that the 3D structure of the vaccine candidate is of high quality and reliability62. ERRAT, by analyzing the statistics of nonbonded interactions between different atom types, generates an overall quality factor of > 50, i.e., 86%, which signifies a high-quality model63.
Toll-like receptors (TLRs), important mediators of inflammatory pathways, are well known for their ability to recognize various pathogen-associated molecular patterns (PAMPs) and play a major role in mediating both innate and adaptive immune responses82. Interestingly, by incorporating a TLR4 agonist (50 S ribosomal protein L7/L12) into its sequence, our multi-epitope vaccine may effectively enhance potential receptor interactions. Thus, molecular docking and dynamic simulations of TLR4 and the vaccine complex were performed to assess the molecular connection, binding affinity, dynamic stability efficiency and structural integrity of the complex. Cluspro analysis revealed a strong affinity (lowest energy) score of −1655.0 kcal/mol with 28 members, which was confirmed by the binding free energy value (− 19.1 kcal/mol) and the dissociation constant (9.8 × 10 − 15 M). The remarkable stability of the complex was maintained by a network of 23 hydrogen bonds and 2 salt bridge interactions.
To evaluate the overall structural stability of the TLR4–vaccine complex during the molecular dynamics simulation, the RMSD of the backbone atoms was analyzed over a 200-nanosecond (ns) trajectory. The RMSD profile (Figure S2) showed an initial sharp increase from ~ 0.3 nm to ~ 2.5 nm within the first 50 ns, indicating rapid relaxation of structural constraints imposed during docking, consistent with equilibration and early conformational adjustments. Between 50 and 150 ns, the RMSD fluctuated between 1.5 nm and 3.5 nm, reflecting the complex’s conformational flexibility, likely due to loop motions or domain rearrangements typical of immune receptors like TLR4. After 150 ns, the system stabilized, with RMSD values plateauing around 3.2–3.4 nm, suggesting convergence to a relatively stable conformational state or local energy minimum. This stabilization supports the reliability of the sampled conformations for further analysis. The high RMSD amplitude (> 2.0 nm) underscores the complex’s structural plasticity, which may be biologically relevant for receptor recognition or signal transduction. Despite these fluctuations, the progressive stabilization observed is indicative of a well-equilibrated, long-timescale simulation.
The surface interaction profile of the TLR4–vaccine complex, derived from the 200-ns simulation, provided insights into the dynamic binding interface, capturing key conformational changes and intermolecular interactions. The resulting interaction map highlights stable binding modes and critical residues involved in molecular recognition and complex stability under physiological conditions. Further RMSD analysis (Fig. 10A) revealed an initial equilibration phase followed by stable fluctuations (~ 0.10–0.125 nm), confirming structural convergence and suggesting that the complex maintains a stable conformation beyond the early relaxation period.
The RMSF profile indicated minimal atomic fluctuations throughout the protein, with only slight peaks observed around residues 190, 300, and 34, likely associated with flexible loops or terminal regions83,84. The generally low fluctuation values reinforce the rigidity and dynamic stability of the receptor–vaccine interface. The structural compaction of the complex was further evaluated using the Rg, which showed a gradual decline from approximately 6.0 to below 4.25 nm during the simulation. This pattern signifies enhanced compactness and structural tightening of the TLR4–vaccine construct85. Likewise, the SASA stabilized between 614 and 634 nm² following an initial peak, indicating consistent folding and reduced solvent exposure. These findings suggest effective maintenance of tertiary structure and conformational folding. Hydrogen bonding analysis revealed a consistently high number of intermolecular hydrogen bonds between the TLR4 receptor and the vaccine construct, ranging from 420 to 500, which reinforces the strong and stable nature of the interaction interface. Additionally, Principal Component Analysis (PCA) was performed to investigate the dominant conformational motions of the complex throughout the 200 ns simulation. The time-resolved PCA projection revealed a progressive and directional exploration of conformational space, with snapshots transitioning smoothly along the first two principal components (PC1 and PC2). The presence of revisited regions suggests intermediate metastable conformations and persistent interactions within a dynamically confined space. This dynamic yet ordered behavior supports a flexible but stable binding profile, which is essential for receptor engagement and immune stimulation. Altogether, the RMSD, RMSF, Rg, SASA, hydrogen bond, and PCA results consistently demonstrate that the TLR4–vaccine complex remains structurally stable, dynamically compact, and functionally competent throughout the entire 200 ns MD simulation.
These findings validate the strength and resilience of the molecular interaction, highlighting the vaccine construct’s potential to persistently engage TLR4 under physiological condition, an essential property for effective immune activation.
According to the NMA simulations carried out by iMODS, the minimal fluctuation observed suggests that the RabVac-TLR4 complex is stable and exhibits proper stiffness motion. Overall, these data suggest that our candidate vaccine might trigger an adaptive immune response via the engagement of TLR4 signaling pathways.
Following in silico immunization, high levels of neutralizing antibodies, CD4 + T helper cells, memory cells, and cytokines were detected after repetitive vaccine doses. During the second and tertiary immune responses, the levels of Igs and their isotypes are increasingly greater than those of the IgM secreted after the primary reaction, with a consistent decrease in the antigen titer.
Moreover, vaccines may trigger the proliferation of TC, TH and long-lasting memory B cells, induce macrophage activity and stimulate the production of natural killer cells and neutrophils. High levels of IFN- and IL-2 release were also observed despite the irrelevant Simpson index (D).
Ensuring efficient expression in the host relies heavily on the optimization of codons in the final construct. Consequently, the codons were optimized to attain an ideal GC content of 58% and a codon adaptation index of 0.92. The prediction of the secondary structure of mRNAs is also essential, as it is a key factor influencing translation initiation, elongation, and mRNA biogenesis86,87. The predicted structure showed high mRNA stability, which was directly related to an improved expression rate.
Conclusion
Taken together, the proposed vaccine demonstrates significant antigenicity, a strong safety record, no allergenic or toxic effects, no human homology, extensive predicted population coverage, stable, and efficacy in inducing a protective immune response, imply that it could be a promising candidate in the fight against rabies infections.
However, subsequent wet laboratory assessments, both in vitro and in vivo, are required to ascertain the vaccine efficacy.
Limitations
While in silico-designed multi-epitope rabies vaccines show great potential, they have significant limitations. One major concern is the accuracy of epitope prediction algorithms, which often fail to fully reflect the multi-level complexity of antigen processing and presentation in vivo. Computational models alone may not perfectly predict in vivo immune dynamics, such as the interaction between innate and adaptive responses, the intricate cytokine signaling pathways involved, and the inducing of lasting immune memory – all critical for vaccine efficacy.
The proposed chimeric vaccine model demonstrates strong stability and immunogenicity in simulations, yet its true protective efficacy against rabies infection remains uncertain without further testing. Over all, future investigations should prioritize in vitro and in vivo studies, including animal challenge experiments to assess immunogenicity and protection efficacy.
Data availability
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.
References
Hampson, K. et al. Global alliance for rabies control partners for rabies prevention. Estimating the global burden of endemic canine rabies. PLoS Negl. Trop. Dis. 9, e0003709 (2015).
Khairullah, A. R. et al. Tracking lethal threat: in-depth review of rabies. Open. Vet. J. 13, 1385–1399 (2023).
Kalthoum, S. et al. Factors associated with the Spatiotemporal distribution of dog rabies in Tunisia. PLoS Negl. Trop Dis. 18, e0012296 (2024).
Afonso, C. L. et al. Taxonomy of the order mononegavirales: update 2016. Arch. Virol. 161, 2351–2360 (2016).
Dietzgen, R. G., Kondo, H., Goodin, M. M., Kurath, G. & Vasilakis, N. The family rhabdoviridae: mono- and bipartite negative-sense RNA viruses with diverse genome organization and common evolutionary origins. Virus Res. 227, 158–170 (2017).
Nishizono, A. & Yamada, K. Rhabdoviruses Uirusu 62, 183–196 (2012).
Coulon, P., Rollin, P. E. & Flamand, A. Molecular basis of rabies virus virulence. II. Identification of a site on the CVS glycoprotein associated with virulence. J. Gen. Virol. 64, 693–696 (1983).
Pulmanausahakul, R., Li, J., Schnell, M. J. & Dietzschold, B. The glycoprotein and the matrix protein of rabies virus affect pathogenicity by regulating viral replication and facilitating cell-to-cell spread. J. Virol. 82, 2330–2338 (2008).
Morales-Martínez, M. E., Rico-Rosillo, G. & Gómez-Olivares, J. L. Aguilar-Setién, A. Immunologic importance of the N protein in the rabies virus infection. Vet. Mex. 37, 351–367 (2006).
Masatani, T. et al. Rabies virus nucleoprotein functions to evade activation of the RIG-I-mediated antiviral response. J. Virol. 84, 4002–40128 (2010).
Finke, S., Mueller-Waldeck, R. & Conzelmann, K. K. Rabies virus matrix protein regulates the balance of virus transcription and replication. J Gen Virol. 84, 1613–1621 (2003).
Okada, K. Roles of the rabies virus phosphoprotein isoforms in pathogenesis. J. Virol. 90, 8226–8237 (2016).
Faber, M. et al. Identification of viral genomic elements responsible for rabies virus neuroinvasiveness. Proc. Natl. Acad. Sci. U S A. 101, 16328–1633246 (2004).
Brunker, K. & Mollentze, N. Rabies virus. Trends Microbiol. 26, 886–887 (2018).
Shepherd, J. G., Davis, C., Streicker, D. G. & Thomson, E. C. Emerging rhabdoviruses and human infection. Biology (Basel). 12, 878 (2023).
Chen, S. J., Rai, C. I., Wang, S. C. & Chen, Y. C. Infection and prevention of rabies viruses. Microorganisms 13, 380 (2025).
Kaye, A. D. et al. Rabies vaccine for prophylaxis and treatment of rabies: A narrative review. Cureus 16, e62429 (2024).
Tarantola, A. et al. Rabies vaccine and rabies Immunoglobulin in cambodia: use and Obstacles to use. J. Travel Med. 22, 348–352 (2015).
Stitz, L. et al. A thermostable messenger RNA based vaccine against rabies. PLoS Negl. Trop. Dis. 11, e0006108 (2017).
Sautto, G. A., Kirchenbaum, G. A., Diotti, R. A., Criscuolo, E. & Ferrara, F. Next generation vaccines for infectious diseases. J. Immunol. Res. 2019, 5890962 (2019).
Deepthi, V., Sasikumar, A., Mohanakumar, K. P. & Rajamma, U. Computationally designed multi-epitope vaccine construct targeting the SARS-CoV-2 Spike protein elicits robust immune responses in Silico. Sci Rep. 15, 9562 (2025).
Abdulhameed Odhar, H., Hashim, A. F., Humadi, S. S. & Ahjel, S. W. Design and construction of multi epitope- peptide vaccine candidate for rabies virus. Bioinformation 19, 167–177 (2023).
Momajadi, L., Khanahmad, H. & Mahnam, K. Designing a multi-epitope influenza vaccine: an immunoinformatics approach. Sci Rep. 14, 25382 (2024).
Antonelli, A. C. B. et al. In Silico construction of a multiepitope Zika virus vaccine using immunoinformatics tools. Sci. Rep. 12, 53 (2022).
Nayak, S. S., Sethi, G. & Ramadas, K. Design of multi-epitope based vaccine against Mycobacterium tuberculosis: a subtractive proteomics and reverse vaccinology based immunoinformatics approach. J Biomol Struct Dyn. 41, 14116–14134 (2023).
Maharaj, L. et al. Immunoinformatics approach for multi-epitope vaccine design against P. falciparum malaria. Infect. Genet. Evol. 92, 104875 (2021).
Yazdani, Z., Rafiei, A., Irannejad, H., Yazdani, M. & Valadan, R. Designing a novel multiepitope peptide vaccine against melanoma using immunoinformatics approach. J Biomol Struct Dyn. 40, 3312–3324 (2022).
Dhanushkumar, T. et al. Structural immunoinformatics approach for rational design of a multi-epitope vaccine against triple negative breast cancer. Int J Biol Macromol. 243, 125209 (2023).
Jiang, F. et al. A comprehensive approach to developing a multi-epitope vaccine against Mycobacterium tuberculosis: from in Silico design to in vitro immunization evaluation. Front Immunol. 14, 1280299 (2023).
Rojas-Caraballo, J. et al. In vitro and in vivo studies for assessing the immune response and protection-inducing ability conferred by fasciola hepatica-derived synthetic peptides containing B- and T-cell epitopes. PLoS One. 9, e105323 (2014).
Jiang, P. et al. Evaluation of tandem chlamydia trachomatis MOMP multi-epitopes vaccine in balb/c mice model. Vaccine 35, 3096–3103 (2017).
Rezaei, M. et al. In Silico design and in vivo evaluation of two multi-epitope vaccines containing build-in adjuvant with Chitosan nanoparticles against uropathogenic Escherichia coli. Int Immunopharmacol. 117, 109999 (2023).
Qiu, J. et al. Integrated in-silico design and in vivo validation of multi-epitope vaccines for Norovirus. Virol. J. 22, 166 (2025).
Doorn, E. et al. Evaluation of the immunogenicity and safety of different doses and formulations of a broad spectrum influenza vaccine (FLU-v) developed by SEEK: study protocol for a single-center, randomized, double-blind and placebo-controlled clinical phase IIb trial. BMC Infec Dis. 17, 1–24 (2017).
Safavi, A. et al. Purification, and in vivo evaluation of a novel multiepitope peptide vaccine consisted of immunodominant epitopes of SYCP1 and ACRBP antigens as a prophylactic melanoma vaccine. Int Immunopharmacol. 76, 105872 (2019).
Wang, C. Y. et al. A multitope SARS-CoV-2 vaccine provides long-lasting B cell and T cell immunity against delta and Omicron variants. J Clin Invest. 132, e157707 (2022).
Kovalenko, A., Ryabchevskaya, E., Evtushenko, E., Nikitin, N. & Karpova, O. Recombinant protein vaccines against human betacoronaviruses: strategies, approaches and progress. Int J Mol Sci. 24, 1701 (2023).
Mahdeen, A. A., Hossain, I., Masum, M. H. U., Islam, S. & Rabbi, T. M. F. Designing novel multiepitope mRNA vaccine targeting Hendra virus (HeV): an integrative approach utilizing immunoinformatics, reverse vaccinology, and molecular dynamics simulation. PloS One. 19, e0312239 (2024).
Ahmed, M. Z. et al. Immunoinformatic-driven design and evaluation of multi-epitope mRNA vaccine targeting HIV-1 gp120. Front. Immunol. 16, 1480025 (2025).
Doytchinova, I. A. & Flower, D. R. VaxiJen: a server for prediction of protective antigens, tumor antigens and subunit vaccines. BMC Bioinform. 8, 4 (2007).
El-Manzalawy, Y., Dobbs, D. & Honavar, V. Predicting linear B-cell epitopes using string kernels. J. Mol. Recognit. 21, 243–255 (2008).
Saha, S. & Raghava, G. P. S. Prediction of continuous B-cell epitopes in an antigen using recurrent neural network. Proteins 65, 40–48 (2006).
Jespersen, M. C., Peters, B., Nielsen, M. & Marcatili, P. BepiPred-2.0: improving sequence-based B-cell epitope prediction using conformational epitopes. Nucleic Acids Res. 45, W24–W29 (2017).
Saha, S., Raghava, G. P. S. & BcePred, Prediction of continuous B-cell epitopes in antigenic sequences using physico-chemical properties. In Artificial Immune Systems (eds Nicosia, G. et al.) 197–204 (Springer Berlin Heidelberg, 2004).
da Silva, B. M., Ascher, D. B. & Pires, D. E. V. epitope1D: accurate taxonomy-aware B-cell linear epitope prediction. Brief. Bioinform. 24, bbad114 (2023).
Dimitrov, I., Bangov, I., Flower, D. R. & Doytchinova, I. AllerTOP v. 2—a server for in Silico prediction of allergens. J Mol Model. 20, 2278 (2014).
Gupta, S. et al. Peptide toxicity prediction. Methods Mol. Biol. 1268, 143–157 (2015).
Calis, J. J. et al. Properties of MHC class I presented peptides that enhance immunogenicity. PLoS Comput Biol. 9, e1003266 (2013).
Larsen, M. V. et al. Large-scale validation of methods for cytotoxic T-lymphocyte epitope prediction. BMC Bioinform. 8, 424 (2007).
Karosiene, E., Lundegaard, C., Lund, O. & Nielsen, M. NetMHCcons: a consensus method for the major histocompatibility complex class I predictions. Immunogenetics 64, 177–186 (2012).
Reynisson, B., Alvarez, B., Paul, S., Peters, B. & Nielsen, M. NetMHCpan-4.1 and NetMHCIIpan-4.0: improved predictions of MHC antigen presentation by concurrent motif Deconvolution and integration of MS MHC eluted ligand data. Nucleic Acids Res. 48, W449–W454 (2020).
Lee, H., Heo, L., Lee, M. S. & Seok, C. GalaxyPepDock: a protein-peptide Docking tool based on interaction similarity and energy optimization. Nucleic Acids Res. 43, W431–W435 (2015).
Xue, L. C., Rodrigues, J. P., Kastritis, P. L., Bonvin, A. M. & Vangone A. PRODIGY: a web server for predicting the binding affinity of protein–protein complexes. Bioinformatics 32, 3676–3678 (2016).
Bui, H. H. et al. Predicting population coverage of T-cell epitope-based diagnostics and vaccines. BMC Bioinform. 7, 153 (2006).
Gasteiger, E. et al. Protein identification and analysis tools on the expasy server. In The Proteomics Protocols Handbook (eds Walker, J. M. et al.) 571–607 (Humana, 2005).
Magnan, C. N., Randall, A. & Baldi, P. SOLpro: accurate sequence-based prediction of protein solubility. Bioinformatics 25, 2200–2207 (2009).
Buchan, D. W. A., Moffat, L., Lau, A., Kandathil, S. M. & Jones, D. T. Deep learning for the PSIPRED protein analysis workbench. Nucleic Acids Res. 52, W287–W293 (2024).
Wang, S., Li, W., Liu, S. & Xu, J. RaptorX-Property: a web server for protein structure property prediction. Nucleic Acids Res. 44, W430–435 (2016).
Du, Z. et al. The trrosetta server for fast and accurate protein structure prediction. Nat. Protoc. 16, 5634–5651 (2021).
Bhattacharya, D., Nowotny, J., Cao, R. & Cheng, J. 3Drefine: an interactive web server for efficient protein structure refinement. Nucleic Acids Res. 44, W406–409 (2016).
Laskowski, R. A., Mac Arthur, M. W., Moss, D. S. & Thornton, J. M. PROCHECK: a program to check the stereochemical quality of protein structures. J Appl Crystallogr. 26, 283–291 (1993).
Wiederstein, M. & Sippl, M. J. ProSA-web: interactive web service for the recognition of errors in three-dimensional structures of proteins. Nucleic Acids Res. 35, W407–410 (2007).
Colovos, C. & Yeates, T. O. Verification of protein structures: patterns of nonbonded atomic interactions. Protein Sci. 2, 1511–1519 (1993).
Ponomarenko, J. et al. ElliPro: a new structure-based tool for the prediction of antibody epitopes. BMC Bioinform. 9, 514 (2008).
Kozakov, D. et al. The cluspro web server for protein–protein Docking. Nat. Protoc. 12, 255–278 (2017).
DeLano, W. L. & Pymol An open-source molecular graphics tool. CCP4 Newsl. Protein Crystallogr. 40, 82–92 (2002).
Laskowski, R. A. et al. Structural summaries of PDB entries. Protein Sci. 27, 129–134 (2018).
Abraham, M. J. et al. High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 1, 19–25 (2015).
Ghedira, K. et al. Design and implementation of a scalable high-performance computing (HPC) cluster for omics data analysis: achievements, challenges and recommendations in LMICs, GigaScience 13, giae060 (2024).
Khamessi, O., Ben Mabrouk, H., ElFessi-Magouri, R. & Kharrat, R. RK1, the first very short peptide from Buthus occitanus Tunetanus inhibits tumor cell migration, proliferation and angiogenesis. Biochem. Biophys. Res. Commun. 499, 1–7 (2018).
López-Blanco, J., Aliaga, R., Quintana-Ortí, J. I., Chacón, P. & E.S., & iMODS: internal coordinates normal mode analysis server. Nucleic Acids Res. 42, W271–276 (2014).
Rapin, N., Lund, O., Bernaschi, M. & Castiglione, F. Computational immunology Meets bioinformatics: the use of prediction tools for molecular binding in the simulation of the immune system. PLoS One. 5, e9862 (2010).
Kazanskii, M. A. et al. RNAfold: RNA tertiary structure prediction using variational autoencoder. BioRxiv 2024.06 (2024).
Parvizpour, S., Pourseif, M. M., Razmara, J., Rafi, M. A. & Omidi, Y. Epitope-based vaccine design: a comprehensive overview of bioinformatics approaches. Drug Discov Today. 25, 1034–1042 (2020).
Hata, H., Phuoc Tran, D., Marzouk Sobeh, M. & Kitao, A. Binding free energy of protein/ligand complexes calculated using dissociation parallel cascade selection molecular dynamics and Markov state model. Biophys. Physicobiol. 18, 305–316 (2021).
Chen, X., Zaro, J. L. & Shen, W. C. Fusion protein linkers: property, design and functionality. Adv. Drug Deliv Rev. 65, 1357–1369 (2013).
Schlake, T., Thess, A., Fotin-Mleczek, M. & Kallen, K. J. Developing mRNA-vaccine technologies. RNA Biol. 9, 1319–1330 (2012).
Kreiter, S. et al. Increased antigen presentation efficiency by coupling antigens to MHC class I trafficking signals. J. Immunol. 180, 309–318 (2008).
Lee, S. J. et al. A potential protein adjuvant derived from Mycobacterium tuberculosis Rv0652 enhances dendritic cells-based tumor immunotherapy. PLoS One. 9, e104351 (2014).
Corradin, G., Villard, V. & Kajava, A. V. Protein structure based strategies for antigen discovery and vaccine development against malaria and other pathogens. Endocr. Metab. Immune Disord Drug Targets. 7, 259–265 (2007).
Porter, L. L. & Rose, G. D. Redrawing the Ramachandran plot after inclusion of hydrogen-bonding constraints. Proc. Natl. Acad. Sci. U S A. 108, 109–1131 (2011).
Sameer, A. S. & S Nissar. Toll-Like Receptors (TLRs): Structure, functions, signaling, and role of their polymorphisms in colorectal cancer susceptibility. Biomed. Res. Int. 2021, 1157023 (2021).
Shah, M. et al. Deciphering the immunogenicity of Monkeypox proteins for designing the potential mRNA vaccine. ACS Omega. 8, 43341–4335530 (2023).
Akhtar, N. et al. Immunoinformatics-Aided design of a peptide based multiepitope vaccine targeting glycoproteins and membrane proteins against Monkeypox virus. Viruses 14, 2374 (2022).
Musa, M. S. et al. Structure-based virtual screening of trachyspermum Ammi metabolites targeting acetylcholinesterase for alzheimer’s disease treatment. PloS One. 19, e0311401 (2024).
Faure, G., Ogurtsov, A. Y., Shabalina, S. A. & Koonin, E. V. Role of mRNA structure in the control of protein folding. Nucleic Acids Res. 44, 10898–10911 (2016).
Georgakopoulos-Soares, I., Parada, G. E. & Hemberg, M. Secondary structures in RNA synthesis, splicing and translation. Comput Struct. Biotechnol J. 20, 2871–2884 (2022).
Author information
Authors and Affiliations
Contributions
W.T. contributed to conceptualization, study design, methodology, investigation, data collection, formal analysis, validation, visualization, and writing—original draft; O.K. conducted molecular dynamics simulations and performed data analysis; H.O. and O.KA. provided specialized scientific expertise, validated results, and contributed to critical manuscript revisions; R.M. assisted in molecular dynamics experiments and analysis; K.G. coordinated and supervised molecular dynamics experiment; A.T. supervised the study, and reviewed the manuscript. All authors reviewed, edited, and approved the final version of the manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Supplementary Information: the online version contains supplementary material available at.
Additional information
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Below is the link to the electronic supplementary material.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
About this article
Cite this article
Tombari, W., Khamessi, O., Othman, H. et al. A novel mRNA-based multi-epitope vaccine for rabies virus computationally designed via reverse vaccinology and immunoinformatics. Sci Rep 15, 30355 (2025). https://doi.org/10.1038/s41598-025-16143-w
Received:
Accepted:
Published:
Version of record:
DOI: https://doi.org/10.1038/s41598-025-16143-w
