Designing an Epitope-Based Peptide Vaccine Derived from RNA-Dependent RNA Polymerase (RdRp) against Dengue Virus Serotype 2

Dengue fever (DF) continues to be one of the tropical and subtropical health concerns. Its prevalence tends to increase in some places in these regions. This disease is caused by the dengue virus (DENV), which is transmitted through the mosquitoes Aedes aegypti and A. albopictus. The treatment of DF to date is only supportive and there is no definitive vaccine to prevent this disease. The non-structural DENV protein, RNA-dependent RNA Polymerase (RdRp), is involved in viral replication. The RdRp-derived peptides can be used in the construction of a universal dengue vaccine. These peptides can be utilized as epitopes to induce immunity. This study was an in silico evaluation of the affinity of the potential epitope for the universal dengue vaccine to dendritic cells and the bonds between the epitope and the dendritic cell receptor. The peptide sequence MGKREKKLGEFGKAKG generated from dengue virus subtype 2 (DENV-2) RdRp was antigenic, did not produce allergies, was non-toxic, and had no homology with the human genome. The potential epitope-based vaccine MGKREKKLGEFGKAKG binds stably to dendritic cell receptors with a binding free energy of −474,4 kcal/mol. This epitope is anticipated to induce an immunological response and has the potential to serve as a universal dengue virus vaccine candidate.


Introduction
Dengue fever (DF) is a vector-borne infectious disease that is transmitted by the dengue virus (DENV) through mosquito bites [1].The prevalence of DF has increased 30-fold in the last 50 years, making it a global health concern.An estimated 390 million DENV infections occur each year, with 96 million exhibiting symptoms.Asia accounts for 70% of these cases [2].The increase in dengue cases demands the development of a vaccine that can induce an immunological response to prevent infections, specifically the development of vaccinations that can induce an immunological response [3].Until now, DF treatment has only been supportive with the aim of providing comfort to the patient and reducing the patient's symptoms [4].Efforts to reduce the high mortality rate caused by DENV require preventive measures.One of them is through vaccination.The development of a vaccine based on a peptide is one of the possible approaches to the problem.Peptide vaccines are made from in vitro synthesized peptides containing 8 to 30 amino acids, which are highly immunogenic and able to elicit an immune response [5][6][7].
DENV is a member of the genus Flavivirus of the family Flaviviridae [8].This virus causes a wide range of diseases in humans, from the acute febrile illness dengue fever (DF) to life-threatening dengue hemorrhagic fever/dengue shock syndrome (DHF/DSS) [9,10].Currently, only a small number of FDA-approved vaccines and antiviral treatments are available for the prevention or treatment of DENV infection.Dengvaxia, which is considered as a safe and effective dengue vaccine candidate, is currently undergoing clinical trials and is the only vaccine currently approved by the FDA [11].
The DENV comprises distinct serotypes [12], with serotype 2 (DENV-2) being the most prevalent and malignant of the four [13,14], especially in the tropics.In comparison to DENV-1 and DENV-4 cases, DENV-2 cases in Brazil had a higher percentage of severe dengue cases [14].However, this does not rule out the possibility that the development of severe manifestations may involve host, epidemiological, and viral factors [15].The virus encodes seven nonstructural proteins including NS1, NS2a, NS2b, NS3, NS4a, NS4b, and NS5 [16].The non-structural protein NS5 of Flavivirus is about 900 amino acids long and comprises a methyltransferase domain at its N-terminus and an RNA-dependent RNA polymerase (RdRp) domain at its C-terminal end, which play a role in the viral replication process [17,18].This protein can serve as promising antiviral development targets [19,20], in addition to its potential as antigen for a dengue vaccine formulation [21].The generated peptides also may have the potential to be employed as vaccine candidates based on epitopes.
Utilizing immunoinformatics permits the design of potential peptide-based epitopes for a DENV vaccine [22,23].The principle of designing a peptide-based epitope vaccine is to predict the interaction between RdRp DENV and dendritic cells [24].This vaccine functions as an antigen that can be targeted to one of the antigen-presenting cells (APCs) to stimulate an immune response.[25].Dendritic cells are immune system cells that control adaptive immune responses and can identify infections that activate T cells and B cells [26].Antigens bound by dendritic cells will be ingested by the cytosol and cut into peptides, which will subsequently be presented in MHC for T cells and B cells, inducing a specific immune response that can be recognized by all organs of the body when utilized as vaccine candidates [26][27][28].It is essential that epitope vaccines interact with dendritic cell receptors in order to elicit an immune response [29].According to studies in mice, targeting antigens to dendritic cells (DCs) can elicit robust CD4+T cell responses [30].Additionally, DCs can process antigens, such as peptides.Several studies involving DCs for vaccine development have been reported [31].Furthermore, directing antigens to DCs, the primary antigen presentation cells and orchestrators of the adaptive immune response, can boost immunogenicity by activating T and B cells [26].
In order to observe the interactions between dendritic cells and epitope-based vaccination candidates, this research involves molecular docking and dynamics modeling experiments.The proposed vaccine is a peptide generated from RdRp DENV-2 that will serve as an epitope-based vaccine, as predicted by the Immune Epitope Database Anal-ysis Resource (IEDB-AR).As DCs play a significant role in a variety of diseases, they have become excellent targets for the development of new strategies to prevent and treat these diseases [29].In this study, DC is utilized in the development of a dengue vaccine for DENV-2.

DENV-2 Sequences Alignment and Reconstruction of Phylogenetic Tree
The DENV-2 sequences were retrieved from NCBI (https://www.ncbi.nlm.nih.gov,accessed on 2 April 2022).Using the FASTA format, the NS5 region encoding RNAdependent RNA polymerase was then aligned.The phylogenetic tree was reconstructed using MEGAX with the neighbor joining method and 1000 replicates of the bootstrap.

B-Cell Linear Epitope Analysis
The DENV-2 sequence from Sumatra (Indonesia) with accession number BAD42415 was selected for linear B cell analysis using the Immune Epitope Database (IEDB) (http: //tools.iedb.org/,accessed on 4 April 2022).The results of the linear epitope analysis for B cells are located in the B cell region or the prediction region.

Receptor Preparation
The three-dimensional protein structure of the dendritic cell that acts as a receptor was downloaded from the protein data bank with PDB ID 3WBP.The Biovia Discovery Studio software was used to prepare the receptor.This preparation aims to obtain a structure that is devoid of water molecules, making molecular docking easier to accomplish.Polar hydrogen atoms were added to the receptor using Autodock Tools [33].Dendritic cell protein structure quality was analyzed on the PROCHECK web server (https://saves.mbi.ucla.edu/,accessed on 11 April 2022) [34] and displayed in a Ramachandran plot [35].

Peptide Modeling and Molecular Docking
The selected epitope obtained based on the results of the analysis was then modeled using the PEP-FOLD web server (https://mobyle.rpbs.univ-paris-diderot.fr/cgi-bin/ portal.py#forms::PEP-FOLD,accessed on 12 April 2022).The peptide epitope was subsequently docked with a prepared dendritic cell receptor.The process of molecular docking is performed on ClusPro web server designed specifically for protein-protein docking [36].The results of molecular docking were visualized using the PyMOL software incorporated into the PDBsum web server (http://www.ebi.ac.uk/thornton-srv/databases/pdbsum/Generate.html,accessed on 15 April 2022) [37].

Molecular Dynamics Simulation Study
Simulations on the molecular dynamics (MD) of peptide and receptor interactions were carried out using a web-based CHARMM-GUI interface system [38][39][40] with a CHARMM36 force field [41].All simulations were performed using the NAMD 2.13 package [42].Utilizing the explicit solvation model TIP3P [43], periodic boundary conditions with dimensions of 106, 106, and 106 in x, y, and z, respectively, were defined.The system was then neutralized with 5Cl − ions.The MD protocol involves minimization, equilibration, and production.The 2 femtoseconds time integration step was selected for all MD simulations.Equilibration was carried out in a canonical ensemble (NVT), while production took place in an isothermal-isobaric ensemble (NPT).Through the production of an MD of 100 ns, the pressure was set at 1 atm using a Langevin Nosé-Hoover piston barostat [44,45] with a Langevin piston decay of 0.05 ps and 0.1 ps period.The temperature was set at 298.15 K using a Langevin thermostat [46].A distance limit of 12.0 Å was applied to short-range non-bonded interactions with a pair list distance of 16.0 Å, and the Lennard-Jones interaction smoothly truncated at 8.0 Å. Long-distance electrostatic interactions were treated using the particle mesh Ewald (PME) method [47,48], where a 1.0 Å grid spacing was used for all simulated cells.All covalent bonds involving hydrogen atoms are restricted using the SHAKE algorithm [49].

MM/GBSA Binding-Free Energy Calculation
Molecular mechanics/generalized born surface area (MM/GBSA) [50,51] is one method implemented in the MOLAICAL code [52] for calculating the relative binding-free energy that occurs when a peptide forms a complex with a protein receptor.The following formula was employed: which can be represented by different interaction contributions, where, the change in gas phase molecular mechanics (∆E MM ), Gibbs solvation energy (∆G Sol ), and conformational entropy (T∆S) is determined as follows: ∆E MM is the sum of the changes in electrostatic energy (∆E ele ), van der Waals energy (∆E VaW ), and internal energy (∆E inter )(bound interaction).∆G Sol is the total of both polar solvation (calculated using the born generalized model) and nonpolar solvation (calculated using the solvent accessible surface area).T∆S was calculated by normal mode analysis.However, this section was neglected because this study focused only on relative binding energies.The dielectric constant (ε) of the solvent was 78.5 and the surface tension constant of 0.03012 kJ mol −1 Å 2 was used for the calculation of MM/GBSA.

DENV-2 Sequences Collection
The first step in this study was to collect the NS5 region of RdRp of DENV-2 sequences from various countries.The collection of DENV sequences from diverse nations yielded 39 DENV-2 and 2 DENV-4 sequences for further study (Table 1).The DENV-4 sequences were served as outgroups.Some researchers suggested that RdRp can serve as the foundation for vaccine development [21,[53][54][55].

Multiple Sequence Alignment
Multiple sequence alignment (MSA) was performed to determine the degree of similarity between each amino acids' sequences.This technique is typically employed to examine the evolution of sequences from a common ancestor in order to identify sequence database commonalities [56].The discrepancy in the alignment findings shows a mutation process in a sequence, whereas the gap suggests the presence of an insertion or deletion.Similarity between protein sequences can facilitate the identification of protein structure and function [57].The MSA results (Supplementary File S1) reveal several variations that result in mismatches due to amino acid alterations.In addition, there is a gap characterized by amino acid sequence deletions.Mutations in the gene producing the NS5 protein generate the amino acids variations.A single mutation can be fatal to the virus.One finding revealed that the presence of a single mutation to alanine at NS5 position R888, which is a residue conserved in all Flavivirus, causes the virus to become completely non-viable [58].This MSA was carried out to determine the conserved region of DENV-2, which will then be used to design peptide-based vaccine epitopes.

Phylogenetic Tree Analysis
Phylogenetic tree analysis using molecular data, such as DNA or proteins, aims to construct precise relationships between genes or taxa.Figure 1 displays the findings of the reconstruction of the phylogenetic tree for the NS5 region.The tree consists of four clades for DENV-2 and one clade for DENV-4.The results of this reconstruction show that Indonesian specimens, particularly Sumatra and Jakarta, are in the same clade as Singapore.This indicates that these specimens are closely related.Sequences that are closely related can be identified by occupying neighboring branches in the tree [59].The genotypes of the specimens are presented in Table 2.According to previous research, each DENV-2 genotype is neutralized differently by vaccine antibodies [60].Therefore, it is assumed that the peptide-based vaccine epitope derived from the DENV-2 conserved region is a viable option.

B-Cell Epitope Analysis
The BepiPred [61] was used for the analysis of B-cell epitopes because the application has a high prediction accuracy on the test data collection, and it has been demonstrated to perform significantly better than other methods.Antibodies in the immune system recognize a molecule by its B-cell epitopes [62].The DENV-2 sequence from Sumatra (BAD42415) was used to predict B-cell epitopes.Figure 2 depicts the outcomes of the prediction.Residues with a score greater than 0.35 are displayed in yellow and are candidates for predicted epitopes [63].There are several predicted peptides that seem able to bind antibodies, as shown in Table 3.A peptide vaccine consists of 8 to 30 amino acid residues and is considered adequate for eliciting the proper cellular and human immune response, while eliminating allergenic and/or reactogenic responses.These vaccines can also be used to protect against multiple strains or serological variants of a given pathogen [64].Using this criteria, seven peptides with lengths ranging from 15 to 25 amino acids (shown in bold) were identified as prospective DENV-2 vaccine candidates.These peptides subsequently underwent predictive immunological testing.

Analysis of Epitope Immunological Properties
The selected candidate peptides for vaccines were next evaluated for immunological features, such as antigenicity, allergenicity, toxicity, and homology with human peptides.Vaccines must have antigenic potential, be non-allergenic, non-toxic, conserved, and not be identical to the human peptide/protein [65,66].Table 4 displays the results of the predicted immunological properties of the vaccine candidates.Almost all of the peptides are antigenic, non-toxic, non-allergenic, and non-homologous to humans, and fit the criteria for epitope vaccine candidates except for GKVRKDIQQWEPSRGWNDWTQ.The peptide MGKREKKLGEFGKAKG, having an antigenicity score of 1.3969, was chosen as the epitope vaccine candidate.Additionally, this peptide is located in a conserved catalytic core domain of RdRp from the positive-sense single-stranded RNA [(+)ssRNA] viruses and closely related viruses.Its antigenicity score was the highest among the peptides studied, indicating that its ability as an antigen to bind to specific antibodies is highly strong.

Target Receptor Preparation
The receptor used is a dendritic cell-associated C-type lectin 2 (PDB ID: 3WBP), which is an innate receptor produced on antigen-presenting cells that recognizes glycosylated pathogens and self-glycoprotein [67].The resolution value of this protein is 1.8 Å and it belongs to Homo sapiens.One of the criteria utilized in receptor selection is the resolution value, where a value less than 3 Å can impair the stability of the receptor during the molecular docking process [68].A low-resolution score is indicative of a more stable receptor.Table 5 displays the outcomes of the validation of the quality of the protein structure using the Ramachandran Plot.Ramachandran plots are the "gold standard" for evaluating new crystal structures [34].
Table 5.The results of the validation of the quality of the 3WBP structure using a PROCHECK program that relies on the Ramachandran Plot.

Ramachandran Criteria Percentage
Most Favored Regions 90.1% Additional allowed regions 9.6% Generously allowed regions 0.1% Disallowed regions 0.1% The results of the analysis show that 90.1% of non-glycine residues are in the most favored regions, 9.6% are in additional allowed regions, 0.1% are in generously allowed regions, and 0.1% are in disallowed regions.A good model, in addition to having very few residues in the disallowed regions, should also have a large number of residues in the most favored regions [35].Thus, the outcomes of this analysis are, therefore, deemed satisfactory.Based on these factors, the protein structure of dendritic cells has a higher quality; thus, it can be stated that the quality of the produced structural model is also good.

Protein Modeling and Molecular Docking
The MGKREKKLGEFGKAKG peptide sequence was modeled using the PEP-FOLD webserver.In addition, the modeling result with the optimal conformation was chosen based on the minimum energy indicated by the optimized potential for efficient structure prediction (sOPEP).This sOPEP energy specifies the modeling structure of the peptide that is close to its original state so that it can interact with the target receptor and give high affinity [69].The peptide was subsequently docked with protein of dendritic cells.The binding-free energy (BFE) value was used to assess the strength of the binding affinity, which indicates how strong a reversible bond can be formed between two or more molecules [70].This is because the vaccine candidate must be able to interact well with the receptor in order to activate the immune response [71].The interaction resulted in a BFE value of −474.4 kcal/mol.BFE is released whenever a ligand forms an association with a receptor, as this results in a reduction in the overall energy of the complex.The release of BFE also compensates for any transformation of the ligand from its minimum energy to its bound conformation with the receptor.This transformation can occur at any point during the process.As a result, the affinity of a ligand for a particular receptor increases in proportion to the amount of energy that is liberated as a result of the ligand's binding to the receptor.The interaction results with the lowest BFE value are optimal.More negativity is better.With this BFE value, it is expected that the epitope vaccine candidate peptide has a very strong interaction with the receptor; hence, possessing a high immunogenic potential.
Figure 3 depicts a molecular docking visualization.This visualization is required to ensure that the epitope can interact with the paratope at the dendritic cell receptor's active site.According to the results of two-dimensional (2D) imaging, there are amino acid residues that establish hydrogen bonds (H-bonds) and hydrophobic interactions with the active site of the receptor.These amino acid residues form H-bonds with interaction distances ranging from 2.51 to 3.02 Å.There are seven H-bonds that bind to the active site region of the paratope at residues Arg151, Gln154, Asp157, Thr159, and Glu172.H-bonds contribute significantly to ligand-receptor interactions and contribute to the stability of the complex [72].The angle at which the bond is formed is another important factor to consider.The strength of the hydrogen bond is directly proportional to the degree to which its geometry is accurate [73].Additionally, the epitope interacts hydrophobically with the paratope at residues Trp153, Val156, Trp155, Tyr161, Asn162, Thr166, Val165, Trp168, Glu172, Pro173, and Asn174.There is also the formation of salt bridges between Lys15 and Asp157, as well as Lys7 and Glu 172.In many cases, the structural driving forces that increase the interaction stability are related to salt bridges [74].The more bonds formed between the peptide and the active site of the receptor, the more stable the interaction, and the more negative the energy, the stronger the potential activity of the epitope vaccine.

Molecular Dynamics Simulations Study
Molecular dynamics simulations (MDS) are an advanced simulation used to examine the movement and interactions of molecules.On the basis of a general physics model that governs interatomic interactions, these simulations make predictions regarding the movement of each atom that makes up a protein or other molecular system over the course of time [75].Moreover, MDS can be used to determine the conformational changes of the docked molecules [76].The simulations are able to analyze the stability and the interaction mechanism formed between the complexes through the analysis of root-meansquare deviation (RMSD) and radius of gyration (Rg) during the simulation based on a predetermined time.RMSD analysis was carried out to compare the conformation of the structure at a certain time to the conformation of the initial structure.The RMSD for the protein-peptide complex was calculated based on the backbone atoms using the visual molecular dynamics (VMD) program, which resulted in an average value of 2.397 ± 0.449 Å (Figure 4).The RMSD graph for the protein backbone reveals that the structure undergoes an increase in RMSD for up to 20 ns, after which the structure remains stable with some fluctuations in the ~1 range, as is typical for globular proteins.The increase in RMSD indicates that the protein structure is beginning to unfold [77], at which point the peptide will search for the appropriate binding site or coordinates.A stable structure denotes that the conformation of the protein bound to the peptide has begun to be attained, allowing the protein to maintain its position.The average RMSD value for the peptide was 4.762 ± 0.311 Å.During the simulation, the RMSD of the peptide remained quite stable, indicating that it remained bound to the protein.
During the simulation, the flexibility of the complex was analyzed using root-meansquare fluctuation (RMSF) (Figure 5), which allowed for the acquisition of additional information regarding the fluctuating protein regions.RMSF is the square root of the average fluctuation between the positions of atoms and the reference structure.Protein flexibility is calculated for each amino acid residue to see how each residue fluctuates or moves during the simulation.During the simulation, the fluctuation occurred around 0.59 nm.This is demonstrated by the presence of multiple peaks on the graph.
The compactness or density of protein molecules can be determined by analyzing Rg.The reduction in the number of fluctuations that occurred while running the simulation is evidence of the increased cohesiveness of the system [78].The Rg for the protein-peptide complex is calculated based on the C-alpha atom with an average value of 24.815 ± 0.372 Å (Figure 6).There is a slight fluctuation with the Rg value of 1 Å during the simulation time, indicating a slight opening and closing of the N-and C-terminal domains.Hydrogen bonding in the complex is crucial for structure stability [33].The total number of H-bonds formed between the protein and peptide during the 100 ns simulation is shown in Figure 7A.All ligands display a fluctuating number of H-bonds with the protein and maintain their bound state.The fluctuations indicated that the peptide altered the protein pocket's conformation.The H-bonds occupancy was calculated as the conformational fraction of 1000 conformations of the protein-ligand complex in which a given residue participates in hydrogen bonding.The 1000 conformations of each complex are derived from the corresponding 100 ns molecular dynamics trajectory.Meanwhile, the average center of mass distance between the ligand and protein over 100 ns of simulation time is shown in Figure 7B.The mean distance is 13.00 ± 0.9766 Å. Contact frequency (CF) analysis was performed to further evaluate the binding between the protein-ligands tested.Interactions between proteins and ligands are vitally important to almost every process that takes place in living organisms [79].The analysis was performed using the contactFreq.tclmodule in VMD with a cutoff of 4 Å.Residues with a higher percentage of CF (>90%) are Arg150, Pro147, HSD152 (prototropic tautomer of histidine), and Glu163 (Figure 8).This indicates that the interaction is within the cutoff, which is typical of van der Waals interactions.The potential energy, pressure, and temperature of the system over a 100 ns MD simulation, as obtained from the NAMD (nanoscale molecular dynamics) log file, is shown in Figure 9.The graph demonstrates that the potential energy, pressure, and temperature remained constant throughout the 100 ns simulation.

Binding-Free Energy
The molecular mechanics-generalized born surface area (MM-GBSA) method was chosen for complex rescoring because it is the fastest force field-based method that calculates binding-free energy compared to other computational free energy methods, such as free energy perturbation (FEP) or thermodynamic integration (TI) methods.Comparative studies also show that MM-GBSA outperform molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) [80].MM-GBSA calculations were conducted using the drug design software MolAICal, which combines artificial intelligence and classical programming [52].The software provides a method for calculating MM-GBSA based on the output results of MDS performed by NAMD [42].

Figure 1 .Table 2 .
Figure 1.Evolutionary relationships of DENV-2 RNA-Dependent RNA Polymerase Protein.The evolutionary history was inferred using the UPGMA method.The evolutionary distances were computed using the Jukes-Cantor Method.Evolutionary analyses were conducted in Geneious.The DENV-4 sequences, which serve as outgroup, are from Brazil.Table2.The genotypes of DENV-2 generated from their evolutionary relationship.ClusterOrigin of Specimen I (American Genotype) Brazil, Puerto Rico, Peru

Figure 2 .
Figure 2. Graph displaying the results of B-cell epitope prediction using BepiPred Linear Epitope Prediction 2.0.The graph illustrates two prediction areas: the area with a yellow peak over the red line (threshold), indicating the area of the B-cell epitope, and the area with a green peak, which is not the prediction area for the B-cell epitope.

Figure 3 .
Figure 3.The molecular docking between the epitope MGKREKKLGEFGKAKG and paratope in the dendritic cells: (a) three-dimensional (3D) visualization showing that epitope binds in the pocket site of paratope; (b) two-dimensional (2D) visualization of the interaction of the complex using LigPlot+ V.2.2; (c) 2D visualization of the interaction generated from PDBSum webserver.The epitope is indicated in green color and the paratope is indicated in blue color.

Figure 4 .
Figure 4. RMSD graph for protein backbone (in blue), peptide (in red), and the protein-peptide complex (in black) during 100 ns simulation.

Figure 5 .
Figure 5. RMSF graph of the protein complex based on C-alpha atoms, which shows the fluctuation of each residue below 0.59 nm during simulation.

Figure 6 .
Figure 6.Rg graph during the simulation of 100 ns.

Figure 7 .
Figure 7. (A) Hydrogen bonding of protein-ligand complex; (B) average distance between peptide and protein during 100 ns simulation.

Figure 8 .
Figure 8. Histogram of contact frequency analysis with a cutoff of 4 Å.

Figure 9 .
Figure 9. Graph of analysis of temperature (K) (in red), pressure (bar) (in green), and potential energy (kcal/mol) (in purple) throughout the 100 ns simulation time.

Table 1 .
The sequences of the NS5 region of DENV-2 RdRp from several countries accessed from NCBI.

Table 3 .
The predicted B-cell epitopes derived from the amino acid sequence of Sumatran specimen.

Table 4 .
Predicted results of the immunological properties of the epitopes as vaccine candidates.The values are enclosed in brackets.
Table 6 displays the BFE results, which indicate that the ligand remained bound throughout the entire simulation time.