Computational identification, characterization and validation of potential antigenic peptide vaccines from hrHPVs E6 proteins using immunoinformatics and computational systems biology approaches

High-risk human papillomaviruses (hrHPVs) are the most prevalent viruses in human diseases including cervical cancers. Expression of E6 protein has already been reported in cervical cancer cases, excluding normal tissues. Continuous expression of E6 protein is making it ideal to develop therapeutic vaccines against hrHPVs infection and cervical cancer. Therefore, we carried out a meta-analysis of multiple hrHPVs to predict the most potential prophylactic peptide vaccines. In this study, immunoinformatics approach was employed to predict antigenic epitopes of hrHPVs E6 proteins restricted to 12 Human HLAs to aid the development of peptide vaccines against hrHPVs. Conformational B-cell and CTL epitopes were predicted for hrHPVs E6 proteins using ElliPro and NetCTL. The potential of the predicted peptides were tested and validated by using systems biology approach considering experimental concentration. We also investigated the binding interactions of the antigenic CTL epitopes by using docking. The stability of the resulting peptide-MHC I complexes was further studied by molecular dynamics simulations. The simulation results highlighted the regions from 46–62 and 65–76 that could be the first choice for the development of prophylactic peptide vaccines against hrHPVs. To overcome the worldwide distribution, the predicted epitopes restricted to different HLAs could cover most of the vaccination and would help to explore the possibility of these epitopes for adaptive immunotherapy against HPVs infections.


Introduction
Human papillomaviruses (HPVs), cervical cancer causing agents, are known to be involved in both morbidity and mortality.Annual epidemics of HPV is approximately 0.5 million while the death rate is about 0.25 million worldwide.Many other disorders such as genital, respiratory, warts and hyper proliferative abrasions are associated with these small DNA viruses [1,2].More than 200 different genotypes of HPVs are characterized.The phylogenetic reconstruction of these genotypes, classified them as Alpha, Beta, Gamma, Mu and Nu.Alpha genus of papillomaviruses are known to be involved in human diseases [3].Among the characterized species of genus Alpha papillomavirus, most of them are associated with the infection of genital tracts [4,5].Sexual intercourse is one of the common ways in the transmission of these viruses.However, fomite transmission as a non-major route of transmission has also been reported [6].
High-risk HPVs (hrHPVs) and low-risk HPVs are the two broad categories of HPV Viruses.Out of the total, 99% of cervical cancers are associated with High-risk HPVs (hrHPVs) species (HPV 16,18,26,31,33,34,35,39,45,51,52,56,58,59, 66, 68 and 70) [7][8][9][10][11].Among the hrHPVs, HPV16 and 18 are responsible for approximately 75% of the total cases.However, low-risk HPV species (i.e., HPV 6,7,11,32,42,43,44,54,61,and 71) are not widely associated with cervical cancer but lead to infection like non-proliferative warts [5,12].Despite the diversity in pathogenicity, all HPVs shares common genome organization.Core and accessory proteins are the two types of genes products in papillomaviruses.Core proteins, E1 and E2, are reported to be directly involved in the viral replication while L1 and L2 are involved in structural assembly.E4, E5, E6 and E7 are considered as accessory proteins, which show variability in both functional aspects and in expression control.The accessory proteins are reported to be involved in virus replication inside infected cell.E6 and E7, two important oncoproteins, are found to be expressed in all positive cases of cervical cancer and are responsible for viral entry, cellular alteration and tumor induction [13][14][15][16].Experimental results proved that the expression of E6 and E7 proteins is the primary cause of the immortalization of primary human keratinocytes in a genome wide study [17][18][19].Beside the carcinogenesis, the deactivation of the tumor suppressor proteins, such as p53 and the retinoblastoma (pRb), is due to the continuous expression of E6 proteins in the cellular environment.Interaction of E6 proteins with E6AP alter the substrate specificity substrate specificity and polyubiquitylates p53, leading to the in degradation of p53 aided by 26S proteasome [20,21].Therefore, the important role of E6 in causing and developing cervical cancer is important and clear.On the other hand, E7 protein perform the function of degradation of pRb and p130 which is a proteasome-dependent process [12,22].
Prediction and development of novel vaccine candidates against the complex diseases has sophisticatedly provoked the desired response and has greatly aided the work of molecular and chemical biologists to expose safe and effective vaccines [23].Immunological mechanisms of surface presentation of antigen along with MHC protein directs the activation of cytotoxic Tlymphocyte (CTL), being an effector, to kill the infected cell.Upon the interaction of CTLs on the infected cell, self-destruction or apoptosis is observed typically.Usually the peptide fragment of the pathogen confers this signaling process and thus provoke the immunity.The underlying mechanism is the attachment of peptide fragment, usually a virulent factor, binds to the MHC molecule and is presented on the surface of infected cells.This process rely on the proteasomal cleavage and transportation to the endoplasmic reticulum (ER) along with MHC molecule.The antigen processing channels (TAP) are required to present the peptide-MHC complex on the surface of the cell for immune response.Therefore, considering the c-terminal cleavage activity and TAP efficiency greatly help in the selection of effective vaccine candidates [24][25][26][27].
The purpose of our present study is to promote the designing of a vaccine against hrHPVs species using in silico methods, taking the most important protein E6 into consideration.The important role of this protein in the cervical cancer carcinogenesis is clear.Therefore, we designed an epitope-based peptide vaccine against hrHPVs that are important but neglected by most of the researchers.To date only HPV 16 & 18 are studied for immunological studies.
To fill this gap and study other hrHPVs this meta-analysis was carried out to predict antigenic potential peptide vaccines using immunoinformatics approach.

Prediction of linear B-cell epitopes
Interaction of B-lymphocytes with antigen B-cell epitope directs the differentiation of B-lymphocytes into memory cells and antibody secreting-plasma [28].The characteristics properties such as accessibility for flexible region and hydrophilic nature are important for B-cell epitopes [29].Different in silico peptide development approaches such Parker hydrophilicity prediction [30], Emini prediction of surface accessibility [31], Kolaskar and Tongaonkar's antigenicity [32], Karplus and Schulz Flexibility Prediction were used using an online analysis resource at IEDB (http://www.iedb.org/).ElliPro [33] (http://tools.immuneepitope.org/toolsElliPro/) is an integrated tool in an online IEDB server which can predict B-cell epitopes using both structural information or protein sequences.ElliPro employ three different algorithms including Protrusion Index (PI) of residues, protein shape approximation and the final neighboring residues clustering which rely on PI.

Prediction of potential cytotoxic T-lymphocyte (CTL) epitopes
NetCTL.1.2[34] (http://www.cbs.dtu.dk/services/NetCTL/) is an online most widely using server for the prediction of CTL epitopes.The Proteomic data from all the hrHPVs were screened to predict potential T-CD8+ (MHC class I binding epitopes) epitopes by using algorithms NetCTL and NetMHC [35,36].NetCTL accept FASTA sequence as an input that perform different analysis such as prediction of MHC class I binding affinity, TAP transport efficiency and C-terminal Cleavage activity.Concerning the MHC alleles, the predictions were restricted to 12 human alleles HLA-A Ã 0101, HLA-A Ã 0201, HLA-A Ã 2402, HLA-A Ã 2601, HLA-B Ã 0801, HLA-B Ã 2705, HLA-B Ã 3901, HLA-B Ã 4001, HLA-B Ã 501, HLA-B Ã 1501, HLA-C Ã 0801 and HLA-C Ã 0202.The weight matrix and artificial neural network was used for the prediction of MHC-I binding and proteasome C-terminal cleavage.

Allergenic and antigenic profiling of B & T-cell predicted epitopes
In order to validate the non-allergenic potential of the predicted B-cell and T-cell epitopes an online web tool AlgPred [37] (http://crdd.osdd.net/raghava/algpred/)was utilized by using multiple algorithms (SVMc, IgEepitope, ARPs BLAST and MAST) to predict the allergenic peptides with an accuracy of 85% by combining these methods.Therefore, we used the primary amino acid sequences to test the allergenic potential of all the selected E6 proteins.On the other hand, to map the antigenic index of our predicted epitopes ANTIGENpro [38] (http://scratch.proteomics.ics.uci.edu/) was used.This server access five different machinelearning algorithms and multiple representation of primary sequences to pile up the antigenicity results by protein microarray data analysis.

Peptides libraries construction and molecular docking
The 3D coordinates of all the selected peptides were predicted by using an online effective web server PEP-FOLD3.For sampling the conformations of predicted peptides simulation runs was set 200 [39].sOPEP energy function integrated in PEP-FOLD3 was applied to cluster the diverse conformational models [40].Selection of specific epitope from all the species was based on low percentile rank and high C-terminal cleavage activity with good TAP score.Sharing of amino acids between the B-cell epitopes and T-cell epitopes were also selected as a docking criteria.Afterward, the best peptide coordinates were docked to the class I MHC molecules HLA-A Ã 0101 (PDB ID 4NQV), HLA-B Ã 1501 (PDB ID 1XR9), HLA-B Ã 5801 (PDB ID 5IM7) and HLA-C Ã 0801 (PDB ID 4NT6) using the PatchDock rigid-body docking server based on the defined threshold [41].PatchDock uses a geometry based docking algorithm to find docking transformations with good molecular shape complementarity [42].Scoring and refining of the docked complexes produced by fast rigid-body docking was performed by employing FireDock server [43,44].Complexes with high global docking energy, Attractive Vander Waal Energy and Hydrogen Bonding energy were subjected to molecular dynamics simulation [45].

Kinetics simulation for the validation of predicted epitopes
A computational systems biology workbench [46] was used to design and execute an in silico biochemical pathway to confirm the antigenic potential of the peptides.Kinetics simulation or pharmacokinetics simulation [47] is useful tool to describe the sufficient dose of a testing drug.The pharmacokinetics of most drugs is first order at therapeutic doses.This non-linear kinetics scheme follow simple Mass kinetics equation or Typically, Michaelis-Menten equation [48] as shown below; This equation can be transformed to Here C represents the steady state concentration, C max the theoretical maximum for C, A the amount absorbed, A max the theoretical maximum for A, and D the dose.
The literature survey was done to collect the necessary information for the hrHPVs.A pharmacokinetics pathway was established to validate the peptides for their potential action.Nodes in the pathway represent the entities, and edges represent the connectivity of one node to another node, which is closely related to each other.In order to carry out pharmacokinetic, concentration doses (0.2 μm) were assigned from available research [49].

Epitope cluster analysis
Clustering of epitopes into groups based on identity among the selected proteins sequences was carried out with the aid of an online Epitope Cluster Analysis tool (http://tools.immuneepitope.org/main/index.html).In the current study, a cluster is a group of sequences sharing a minimum of 80% of the sequence identity is known to be a cluster.

Molecular dynamics simulations
MD simulations of all the selected complexes were carried out by using AMBER 14 molecular dynamics package [50].To neutralize the systems counter Na+ ions and hydrogens were added.The tleap package of Amber was utilized to perform this process.A TIP3P water box of 8.0 A˚was used.A two stages energy minimization, each of 6000 steps, of the complexes using the SANDER module of AMBER 14, was performed to remove the constraints all atoms in the systems except those from the water molecules.PMEMD.cuda[51] unit of AMBER 14 was used to accomplish MD simulations of the minimized complexes.For long-term interactions, the SHAKE algorithm and Particle-Mesh Ewald (PME) method was used and a non-bond contacts cutoff radius of 10A˚was kept.Using the Langevin temperature 310K and constant pressure (1atm) with isotropic molecule-based scaling was considered for equilibration of 10,000 ps time, followed by a total of 20ns simulation was carried.Sampling of MD trajectories was carried out after every 2.0 ps time scale.Analysis such as RMSD and Hydrogen bonding analysis was carried out by using an integrated programs CPPTRAJ and PYTRAJ [52] in AMBER 14.The following equation was solved to calculate the stability of the complexes after 20ns.

RMSD
Where N is the number of atoms, m i is the mass of atom i, x i is the coordinate vector for target atom i, Y i is the coordinate vector for reference atom i, and M is the total mass.If the RMSD is not mass-weighted, all m i = 1 and M = N.

Hydrogen bonding analysis
Hydrogen bonds are an important non-covalent structural force in molecular systems.They are formed when a single hydrogen atom is effectively shared between the heavy atom it is covalently bonded to (the hydrogen bond donor) and another heavy atom (the hydrogen bond acceptor).Here, we analyzed the hydrogen bonds between all the selected complexes.Hydrogen Bonds were analyzed at three different stages.The bonds were checked before the simulation and after the simulation.After the minimization and production in the PDB coordinated of the complexes were saved from the .rstfiles and were analyzed by using an online server PDBePISA [53] UCSF Chimera [54] and PyMOL [55] visualization software.

Antigenic B-cell epitopes prediction
The antigenic epitopes were determined by using Tongaonkar's method [32] using the physiochemical properties of amino acids.Experimental precision for this method is observed to be 75% [32].Four antigenic peptides in each E6 protein of HPV31, HPV35, HPV45, and HPV68 were predicted that range from 7-15 amino acids.Moreover, HPV39, HPV51, HPV52 and HPV58 possess five antigenic epitopes, ranging from 6-14 amino acids respectively.HPV56 contain three antigenic epitopes while HPV33 possess six antigenic epitopes.The range of HPV56 epitopes is 9-14 amino acids while HPV33 peptides range from 6-14 amino acids.The B-cell epitopes predicted are shown in the Table 1.
To predict the maximum residual score for each amino acid in the E6 of hrHPVs species Kolaskar and Tongaonkar's was used.Proteins with residual score >1 were quantified.Among the total selected proteins large number of residues with score greater than 1 were found, which is showing the antigenic potential of E6 protein

Surface accessibility for E6 proteins
The surface probability of each residue was predicted using a threshold 1.0.Amino acids with the score greater than 1 has the highest probability to be found on the surface [31].The minimum surface probability score 0.05 for HPV31 from amino acid position (ICDLLI 96-101 ), 0.032 for HPV33 from amino acid position (IRCIIC 101-106 ), 0.034 for HPV35 from amino acid position (ICLNCV 26-31 ), 0.035 for HPV39 from amino acid position (IRCMCC 103-108 ), 0.032

Surface flexibility for E6 proteins
Temperature or B factor is used to demonstrate the back and forth motion of atoms within a protein coordinates.To calculate the motion of atoms Karplus and Schulz's flexibility method was implemented.Atoms with profoundly systematized structure appeared to have low B-factor while the distorted appeared higher [56]

Parker hydrophilicity prediction for E6 protein
Hydrophilicity of the predicted peptide was calculated based on retention times of a peptide during HPLC using reversed phase column.Here, we used Parker hydrophilicity prediction method to predict the water loving potential of the predicted antigenic peptides.Immunological studies reported the direct association of hydrophilic region with the antigenic sites [30].S4 Fig is showing the graphical illustration of predicted Parker Hydrophilicity of E6 Proteins of hrHPVs on the basis of the x-axis is showing the position of the amino acids and y-axis is plotting the hydrophilicity.Among the selected species, the lowest hydrophilicity was calculated as -5.086 from all the E6 protein of HPV35 and HPV58 from amino acid position LCHLLIR 96-102 and ILIRCII 99-105 .These regions were predicted to act as active T-cell epitopes.
On the other hand the maximum hydrophilicity score 6.371 was calculated for HPV35 E6 protein for the peptide sequence QDTEEKP 3-9 .Moreover, the maximum and minimum hydrophilicity score for all the other HPV species are shown in the S3 Table.

Antigenic T-cell epitopes prediction
Cytotoxic T-lymphocyte (CTL) epitopes were explored from the E6 proteins of the selected hrHPVs.NetCTL 1.2 server [34] was utilized to predict the CTL epitopes.MHC binding affinity, proteasomal C-terminal cleavage, TAP transport affinity and potential MHC ligands were recognized by following > 0.75000 threshold as criteria.A total of 38 peptide sequences from E6 Proteins of all the selected hrHPVs were predicted as CTL epitopes whose prediction score were > 0.75000 (Table 2).

Molecular docking of E6 proteins with HLA-A Ã 0101
Among the total 38 epitopes, only 9 epitopes, were docked to MHC class I HLA-A Ã 0101, 2 epitopes against the HLA-B Ã 1501, 1 epitope against HLA-B Ã 5801 and 3 epitopes against HLA-C Ã 0801 from all the selected HPV species as shown in the Table 3.Initially, epitopes predicted against the 12 human alleles with good percentile rank, TAP transport efficiency and C-terminal Cleavage activity were selected to be analyzed for the binding efficiency [57] approach.A 2% percentile was used in the present study.Because it has been reported that using a defined threshold of percentile rank and MHC binding affinity, most of the predicted epitopes provoked the immune response in experimental condition.The global and attractive van der Waals energy (vdW) were computed ranging from -22.25 to -54.31 kcal/mol and -17.52 to -31.80 to determine the binding efficiency of each epitope respectively.The docking scores along with the bonding pattern is presented in the Table 4.Among the essential features Asn77, Tyr99, Arg114 and Arg156 residues from the MHC-I groove were most abundantly involved in bonding with different predicted peptides.Among the total 10 epitopes docked against HLA-A Ã 0101 (ETEVLDFAF, RSEVYDFAF and KTLQRSEVY) showed the highest binding affinities and highest number of hydrogen bonding within 3Å.However, the other epitopes also showed good affinities.This analysis established a good interaction of the modeled antigenic peptides with the MHC-I molecules.Nevertheless, residues shared by MHC-I in bonding with different peptides are also reported by Mirza, Rafique et al. [58] in the same computational study while predicting antigenic epitopes.Furthermore, residues from different epitopes such as (Arg5, Glu3 and Thr4) were predicted to act as antigenic.The graphical representation of these docked complexes are given in the Fig 2 .The hydrogen bonds with length less than 3Å were most frequently found in all complexes.Overall stability of the docked complex seems to be well preserved by the formation of hydrogen bonds.

Molecular docking of E6 proteins with HLA-B Ã 1501 and HLA-B Ã 5801
Epitope ATLERTEVY was found to be good against HLA-B Ã 1501 with the percentile rank 1.9.While epitope KTLQRSEVY with percentile rank 1.3 and 0.9 restricted to both HLA-B Ã 1501 and HLA-B Ã 5801 showed good affinity.The global energies of each of these epitopes restricted  specific allele were -30.47, -24.89 and -21.19 kcal/mol while the attractive van der Waals energy (vdW) were ranging from -18.84 to -20.29 kcal/mol respectively.Asp61, Tyr59, Asn66, Arg97 and Ser73 from these three MHC-I molecules were frequently in vigorous interaction with the predicted peptides.Complexes after the docking were subjected to analyze the hydrogen bonding within 3Å.Variable hydrogen bonds (2 to 4) within 3Å were found.Molecular interactions of the specific peptides and their respective docked poses are tabulated in Table 5 is tabulating the atomic feature involved in interactions while the graphical illustrations are given in the Fig 3.

Molecular docking of E6 proteins with HLA-C Ã 0801
Among the total epitopes, only three epitopes restricted to HLA-C Ã 0801 were found to have the percentile rank below 2%.Docking of these peptides against MHC class I HLA-C Ã 0801 revealed good atomic interactions.The interacting global energies and attractive van der Waals energy (vdW) were measured ranging from -15.97 to -30.36 kcal/mol and -22.99 to -24.88 respectively.These peptides were selected as specie specific.Arg97, Ser77, Asn114, Gln155 from the MHC-I groove were most abundantly involved in bonding with different predicted peptides.Post-docking analysis confirmed the stability of the complexes by observing  6, while the structural complexes before and after simulation are given in the S5 Fig.

Pharmacokinetics simulation based validation
Biological network representations have been used to explain interaction as well as mechanisms between entities, and the biochemical mathematical models of the network have been used to study the pharmacokinetic mechanism in the specific systems.Pharmacokinetics simulation, as shown in the Fig 5, reported that during High Risk HVPs Infection, pRb activates the E2F which infect E7, and combined E2F and E7 activates the E2F, which make basal, and parabasal layer and convert into transactivation which up regulate the genes necessary for S-phase progression.Up regulation of E7 and down regulation of p21 form E7+p21 complex, where p21 inactivated and Cyclin E/cdk present in activated form at low levels.While down regulation of E7 and up regulation of p21 form p21+E7+Cyclin E complex, where Cyclin E/cdk inactivated and present in activated form at high levels.Peptide inactivates the MDM through p14 ARF, which degrade p53.This mechanism of degrading the tumor suppressor genes MDM and p53 is already reported by the experimental studies and thus increasing the potential of our study [21,22].Furthermore, we have used an experimentally reported concentration to obtain the similar results.The binding affinities of all the docked epitopes are given in the graph Fig 6(A).

Epitopes cluster analysis
All of the predicted epitopes for hrHPVs E6 were grouped into 16 clusters as shown in the (S4 Table ).Among the 16 hrHPVs' E6 clusters, cluster 3 and 5 contain the most 7 epitopes out of the total selected 10.On the other, hrHPVs E6 clusters, cluster 4 contained the most epitopes (4 epitopes), while the rest of the clusters possess two and one epitopes respectively.

Discussion
Vaccination is one of the significant approaches to improve the standard of community health, provide the best and safe way to control the increasing infections leading to complex diseases.Modern therapeutics developments greatly rely on immunoinformatics resources for the synthesis of antigen-specific epitopic therapeutic vaccines against the wide diversity of pathogenic infections [34,59].Using these in silico methods, highest accuracy is reported by different applied groups.Epitopic vaccine against HIV, malaria and tuberculosis provided promising results and supported the defensive and therapeutic uses of these vaccines [60].Genomic and proteomic information of Human Papillomaviruses delineate that E6 and E7 proteins has significant therapeutic value that could be targeted to prevent the progression of HPVs persistent infection.Targeting E6 protein is important for the immunotherapy of cancer as these are reported to be important in the survival and maintenance of oncoproteins in the malignant cell [61][62][63].The present study predicted epitopes for twelve different MHC class I HLAs to provide wide spectrum of peptides from HPVs E6 proteins that could strongly provoke the immune response.Detail analysis of the predicted epitopes such as percentile rank, MHC binding affinity, TAP, C-terminal cleavage activity, antigenicity and allergenic profiling was carried out to select the most promising epitopes as this criteria is experimentally validated for immune response potential [64][65][66].The epitopes predicted in this study could become a clinical candidate sooner or later for the treatment of HPVs infection and cervical cancer, as epitopes such as KLPQLCTEL 18-26 and FAFRDLCIV 52-60 of E6 proteins, have been tested in a transgenic mice for IC 50 value which resulted in immune response in experimental conditions [67].Previous studies already verified that peptide FAFRDLCIVYR 52-62 possess antitumor effect and is reported to be processed by T-cell endogenously [68,69].The current study predicted epitopes for the E6 proteins, which also contain regions from epitope, such as FAFRDLCIVYR 52-62 , thus increasing the accuracy and validating our study to be taken to the experimental room.Furthermore, we tested our predicted epitopes RREVYDFAF 46-54 (Cluster 3) for the optimal peptide properties, using PepCalc, which resulted in promising water solubility and other good experimental properties.Our results showed that combining cluster 3 and cluster 4 will form a continues epitope that will elicit robust immune response against multiple species of HPVs thus acting as more active prophylactic vaccine candidate, as present in seven species, for the clinical exploration.Previous study also reported that the combined epitopes PYAVCDKCLKF 66-76 of the E6 protein presented by HLA-A Ã 1101 and HLA-A Ã 2402 with a worldwide population frequency ranging 29.5-30.5% possess positive behavior.As a result, the epitope combinations we predicted (Cluster 5 65-75 ) might be feasible in therapeutic vaccines of hrHPVs [70].The purpose of our study was to predict and select T-cell epitopes because they are more promising and evoke a long-lasting immune response, and because with antigenic drift, an antigen can easily escape the memory response of antibody.However, in vaccine development, allergenicity is a prominent hurdle [71].Today, most vaccines stimulate the immune system into an 'allergic' reaction, through induction of type 2 T helper cells and immunoglobulin E. Therefore, we also confirmed the non-allergenic potential of our validated epitopes.The allergenicity testing results of these epitopes showed that our epitopes are safe and will not show any allergenic reaction possibly.Furthermore, Docking, MD simulations and PKPD modeling approach also verified that our peptides has experimental potential.
The current state of the study strongly support the development of peptide based vaccines against the HPVs species that are active and has largest number of infections but neglected.To date only one or two species of HPVs are studied but no such large meta-analysis integrated with dynamics is reported.This meta-analysis of multiple therapeutically important species has increased the scope and accuracy of this study.Our results specified the regions from 46 to 54, 65 to 76 and a combined epitope from 46 to 62 could be used for development of candidate CTL epitopes.To overcome the worldwide distribution of the major epidemic of hrHPVs, the predicted epitopes restricted to different HLAs could cover most of the vaccination and desire outcome.This study aided the development of prophylactic vaccines that will provide cross immunity against multiple cervical causing species.

Conclusion
In the current study, we utilized advanced immunoinformatics and molecular dynamics simulations approaches to predict and verify potential T-cell epitopes to act as vaccine candidates against the hrHPVs.Many previous studies conveyed different information only about the major HPV16, 18 and HPV45 but none of the study focused the other important performers of cervical cancer.We have disclosed a wide range of information and predicted potential vaccine agents that will act against multiple species of cervical cancer causing agents.Prediction, docking, post docking and simulation analysis will aid the development of prophylactic vaccines against hrHPVs.The advantage of this study is that it covers many aspects along with multiple risk species of cervical cancer.However, further experimental insight will be fruitful in taking the predicted potential peptides into the clinical room and market.

Fig 1 .
Fig 1.Schematic of the workflow for the prediction and validation of peptide vaccines from hrHPVs E6 proteins.The above figure is showing the whole methodology including the resources and results obtained from these analyses.https://doi.org/10.1371/journal.pone.0196484.g001 . The graphical illustration, given in S1 Fig, of predicted antigenic propensity, maximum and minimum residual score and number of residues with residual score >1 are given in the S2 Table.

Fig 2 .Fig 3 .
Fig 2. Molecular interaction analysis of predicted HPVs E6 peptides docked to MHC-I HLA-A Ã 0101.(A) ETEVLDFAF, (B) RSEVYDFAF, (C) FQDPAERPY, (D) QTEVYEFAF, (E) ATLERTEVY, (F) EATIKKSLY, (G) FTDLRIVYR, (H) LCDLLIRCY, (I) KTLQRSEVY, (J) ETITNTKLY.https://doi.org/10.1371/journal.pone.0196484.g002 Stability of the complexes was examined by analyzing the trajectories from simulation.Root mean square deviation (RMSD) of the Cα atoms was calculated to analyze the stability using CPPTRAJ.The RMSD as shown in the Fig6(B) of all the simulated systems are reaching the equilibrium is shorter and the averaged RMSD values of all the bound complexes range from 0.1 to 0.25 (nm).These results imply that the bound peptides restrain the motions of MHC-I and favor the stability of the complex structure.

Fig 5 .Fig 6 .
Fig 5. (A) Pharmacokinetics simulation showing the pathway of degrading the p53 genes in HPV infections, while in (B) X-axis represents the transition time of entities and Y-axis represent the concentration of the peptides.https://doi.org/10.1371/journal.pone.0196484.g005 Table.The schematic flow of this work is given in the Fig 1.

Table 2 . Predicted CTL epitopes from (HPV31 to HPV68) E6 protein. Prediction score threshold was set at > 0.75000.
Bold indicates amino acids that were also predicted as antigenic sites.

Table 3 . Binding affinity against each allele is determined in terms of percentile rank by IEDB for MHC class I.
Based on the low percentile rank (2%) these epitopes were docked against the specific allele.

Table 4 . HPV E6 peptides-HLA-A Ã 0101 interaction.
FireDock energy for the best ranked complex initial distance between the H-bond donor and the acceptor; measured with the Find H.Bond tool in Chimera (H-Bond constraints were relaxed by 1 Å and 20.0 degrees) distance between the H-bond donor and the acceptor after molecular dynamics simulation (MD); measured in PyMOL, nd = no detected H-bond.

Table 6 . HPV E6 peptides-HLA-C Ã 0801 interaction.
FireDock energy for the best ranked complex initial distance between the H-bond donor and the acceptor; measured with the Find H.Bond tool in Chimera (H-Bond constraints were relaxed by 1 Å and 20.0 degrees) distance between the H-bond donor and the acceptor after molecular dynamics simulation (MD); measured in PyMOL, nd = no detected H-bond.