Submitted:
02 November 2026
Posted:
03 November 2026
You are already at the latest version
Abstract
Flagellin activates innate, cellular and humoral immunity through Toll-like receptor 5 (TLR5), and this interaction can be exploited to improve vaccines against oncogenic human papillomavirus (HPV). The binding mechanism between flagellin-fabricated peptide vaccines (FFPVs) and TLR5, however, remains unclear. Here, three FFPV candidates (C58, C52 and C16) were designed from overlapped epitope segments in the conserved N-terminal region of the L2 minor capsid protein of HPV58, HPV52 and HPV16. These segments were predicted to induce cytotoxic T-lymphocyte, T-helper, IFNγ and B-cell responses. Each construct also carried two TLR agonist adjuvants (flagellin and RS09), two universal T-helper epitopes and defined peptide linkers. Physicochemical and immunological properties were evaluated, and three-dimensional models of the FFPVs and TLR5 were generated, refined and validated. Protein–protein docking with HADDOCK and PRODIGY ranked the TLR5-C58 complex highest, ahead of TLR5-C16 and TLR5-C52. To test the docking result, the TLR5-C58 complex was simulated for [100 ns] in explicit solvent, and the root mean square deviation, total energy, hydrogen bonds and key-residue distances were analyzed. The simulation identified a binding mode with ten key residues in the leucine-rich repeat (LRR) region of TLR5, and the binding free energies agreed with experimental reports. These findings offer a structural basis for the rational design of L2-based vaccines against cervical cancer.
Keywords:
homology modelling
; HPV
; molecular dynamics simulation
; minor capsid protein L2
; overlapped epitope segments
; prophylactic vaccine
; Rosetta
; TLR5
1. Introduction
Cervical cancer is the fourth most common cancer in women worldwide, and persistent infection with high-risk HPV (hrHPV) is its principal cause [1]. Beyond HPV16 and HPV18, other hrHPV genotypes make a large contribution to disease in particular regions. In Chinese women, HPV52 and HPV58 are among the most frequently detected hrHPV types after HPV16, and surveys in several provinces report HPV52 infection rates that have risen or remained high over recent years [2,3]. Three licensed prophylactic vaccines (Cervarix®, Gardasil® and Gardasil 9®) are made from recombinant L1 virus-like particles (VLPs). The nine-valent vaccine covers HPV6, 11, 16, 18, 31, 33, 45, 52 and 58 and is highly efficacious against infection and high-grade lesions caused by these types [4]. These vaccines nevertheless protect mainly against the genotypes they contain, have no therapeutic effect in people already infected, and remain costly to manufacture and deliver, which limits access [5,6]. The minor capsid protein L2 is an attractive alternative antigen. About 72 copies of L2 (~470 amino acids) sit beneath the L1 capsomeres of each virion, and its N-terminal region is highly conserved across low-risk and high-risk HPV types, in contrast to the type-restricted epitopes of L1 [7,8]. Vaccination with full-length L2 or with N-terminal L2 peptides elicits cross-neutralizing antibodies in animal models, and the 17–36 region carries a broadly protective epitope [7,9]. The main obstacle is the weak immunogenicity of monomeric L2, which does not multimerize or assemble into VLPs on its own [6,10,11]. A rational remedy is to present L2 epitopes together with an adjuvant. Flagellin–L2 fusion vaccines protected mice and rabbits against challenge with diverse papillomavirus genotypes [12,13].
TLRs are pattern recognition receptors that sense invading pathogens and trigger pro-inflammatory signalling [14]. TLR5 recognizes bacterial flagellin through the horseshoe-shaped leucine-rich repeat (LRR) ectodomain, and the structure of the TLR5–flagellin complex has been solved [15]. TLR5 expression increases during the progression of cervical neoplasia, which supports its relevance in the cervical setting [16]. Flagellin has four domains (D0–D3). The highly conserved N-terminal and C-terminal D0/D1 regions are required for TLR5 activation [17], whereas the hypervariable D2/D3 domains differ markedly among bacteria and are dispensable for it. Replacing D2/D3 with foreign antigens promotes mucosal IgA production and attenuates flagellin-induced inflammation after intranasal immunization [18]. We therefore retained only the D0/D1 regions of flagellin in our vaccine design. Immuno-informatics allows candidate epitopes to be screened rapidly and at low cost before experimental work [19]. We previously designed chimeric L2 constructs from overlapped CTL, helper T-cell and B-cell epitope segments of HPV58 (C58) [20,21], and an in silico designed HPV16 L2 epitope vaccine that also uses flagellin (C16) has been reported by others [22]. How these flagellin-fabricated peptide vaccines engage TLR5 at the atomic level, and which of them binds best, has not been examined. In this study we describe a third candidate (C52) derived from the HPV52 L2 N-terminus and compare it with C58 and C16. We evaluated the physicochemical and immunological properties of the three constructs and modelled, refined and validated their three-dimensional structures. We then compared their binding to TLR5 using protein–protein docking with two independent scoring schemes and analysed the intermolecular hydrogen-bond networks. The most favorable complex was examined by molecular dynamics (MD) simulation. The resulting binding mode is a starting point for designing FFPV candidates that engage both innate and adaptive immunity against cervix papilloma.
2. Materials and Methods
2.1. Sequence Retrieval
2.2. Epitope Prediction and Selection for HPV52
Epitopes were predicted for the conserved N-terminal region of HPV52 L2 (residues 12–114). MHC class I-binding epitopes were predicted with the IEDB consensus method [24] and NetMHC 4.0 [25] for eleven HLA class I alleles (HLA-A*01, A*02:01, A*02:07, A*11:01:01, A*24:01, HLA-B*40, B*46:01, B*58:01, HLA-C*03:01, C*07:02 and C*12:03). MHC class II-binding epitopes were predicted with the IEDB consensus approach [26] for ten alleles (HLA-DQB1*03:01, DQB1*03:03, DQB1*06:01, DRB1*07, DRB1*09, DRB1*14:01:01, DRB1*15:01, DRB1*15:07:01, DPB1*05:01 and DPB1*05:02:01). The alleles were chosen because they are common in Chinese populations [27]. CTL epitopes were cross-checked with CTLPred [28] and PAComplex [29]. IFN-γ-inducing epitopes were predicted with IFNepitope [30], linear B-cell epitopes with ABCpred [31], allergenicity with AllerTOP v2 [32] and AllergenFP [33], and antigenicity with ANTIGENpro [34]. Cross-reactivity with the human proteome (Taxonomy ID 9606) was assessed by peptide matching, and conservancy was calculated across the L2 sequences of 15 hrHPV genotypes (HPV16, 18, 31, 33, 35, 39, 45, 51, 56, 58, 59, 68, 69, 73 and 82). Only overlapped segments that were predicted to be MHC-I and MHC-II binders, IFN-γ inducers, non-allergenic and antigenic, with no similarity to human proteins, were retained.
2.3. Vaccine Construct Design and Evaluation
Each construct comprised two overlapped epitope segments, two TLR agonist adjuvants (the N-terminal 5–143 and C-terminal 419–504 regions of S. dublin flagellin as the TLR5 agonist, and the TLR4 agonist peptide RS09 [35]), and two universal T-helper epitopes (PADRE [36] and TpD [37]). The seven parts were joined with GPGPG, AAY and EAAAK linkers [38]. Molecular weight, theoretical pI, charged residues, extinction coefficient, half-life, instability, aliphatic and GRAVY indices were computed with ProtParam [39]. Solubility upon overexpression in E. coli was predicted with SOLpro [40], and protective antigenicity with VaxiJen at a threshold of 0.4 [41] and with ANTIGENpro [34].
2.4. Structure Modelling, Refinement and Validation
Three-dimensional models of C58, C52, C16 and TLR5 were built with the Robetta server, which combines ab initio and comparative modelling of protein domains [42]. Five models of each protein were refined with GalaxyWEB [43] and 3Drefine [44]. Model quality was assessed with ProSA-web (Z-score) [45], Ramachandran plot analysis with RAMPAGE [46] and ERRAT [47]. Structures were visualized with UCSF Chimera [48].
2.5. Protein–Protein Docking and Binding-Affinity Estimation
HADDOCK. Each FFPV was docked to TLR5 with the HADDOCK 2.2 web server [49]. The HADDOCK score is a weighted sum of van der Waals, electrostatic, desolvation and restraint-violation energies and buried surface area (BSA); lower scores indicate more favourable binding.
PRODIGY. Binding affinity (ΔG) and the dissociation constant (Kd, at 25 °C) of each complex were predicted from intermolecular contacts with PRODIGY [50]. Interface analysis. Intermolecular hydrogen bonds and non-bonded contacts were identified and plotted with LIGPLOT v4.5.3 [51]. The three complexes are referred to as complex I (TLR5–C58), complex II (TLR5–C52) and complex III (TLR5–C16).
2.6. Molecular Dynamics Simulation
The TLR5–C58 complex, which had the best docking scores, was simulated with GROMACS 5.1.2 [52] using the CHARMM force field [53]. The complex was placed in a cubic box of TIP3P water [54], and counter-ions were added to neutralize the system. Energy minimization used the steepest-descent algorithm for up to 50,000 steps with a force tolerance of 1000 kJ mol⁻¹ nm⁻¹. Short-range van der Waals and Coulomb interactions were truncated at 1.0 nm, and long-range electrostatics were treated with the particle mesh Ewald method [55]. Bonds and water geometry were constrained with LINCS [56] and SETTLE [57]. The system was equilibrated in the NVT and NPT ensembles at 300 K and 1 bar, with temperature controlled by the V-rescale thermostat [58] and pressure by the Parrinello–Rahman barostat [59]. Production simulation was run for [100 ns] with a 2 fs time step, and coordinates were saved every 2 ps. Trajectories were analysed for the root mean square deviation (RMSD), total energy, number of intermolecular hydrogen bonds and distances between key interface residues. Binding free energy (BFE) was calculated from the trajectory.
3. Results
3.1. Selection of Overlapped Epitope Segments
Overlapped epitope segments from the N-terminal region of HPV52 L2 were identified by combining MHC class I, CTL, MHC class II and IFN-γ predictions (Supplementary Table S1 and Table S2). Three segments (residues 10–23, 23–36 and 29–42) were predicted by all tools. Of these, residues 23–36 and 29–42 were non-allergenic by both AllergenFP and AllerTOP, antigenic by ANTIGENpro and IFN-γ-inducing (scores 0.389 and 0.505), and showed no similarity to human proteins (Table 1). Three linear B-cell epitopes overlapped these segments: VIPKVEGTTIADQLLK (residues 20–35), YQTCKASGTCPPDVIP (7–22) and GTTIADQLLKYGSLGV (26–41). The 23–36 and 29–42 segments therefore have the potential to induce both cellular and humoral immunity, and they were conserved across the 15 hrHPV genotypes examined (Supplementary Table S3). The corresponding segments of HPV58 (C58) and the 10–36 and 65–89 segments of HPV16 (C16) were included for comparison (Table 1).
3.2. Design and Properties of the Vaccine Constructs
The C52 and C58 constructs each contain 318 amino acids and C16 contains 341. All three comprise the epitope segments, the TLR adjuvants (flagellin and RS09), the T-helper epitopes (PADRE and TpD) and the GPGPG, AAY and EAAAK linkers, with the N-terminal (5–143) and C-terminal (419–504) flagellin regions at the head and tail (Table 1 and Table 2).
All three constructs were predicted to be stable (instability indices 35.88, 32.38 and 35.47 for C58, C52 and C16; Table 2). C16 had the highest molecular weight (35.57 kDa), the highest pI (9.28) and the most charged residues (52). VaxiJen scores (0.4301, 0.4468 and 0.4870 for C58, C52 and C16) exceeded the 0.4 threshold, and ANTIGENpro scores were 0.9438, 0.9490 and 0.9482. The predicted solubility upon overexpression was 0.621, 0.587 and 0.701 (Table 2).
3.3. Structure Modelling, Refinement and Validation
Robetta produced five models each for TLR5 and the vaccine constructs. Models 1 and 4 were de novo models and models 2, 3 and 5 were ab initio models. The Ginzu domain-prediction confidence scores were 0.9375 for TLR5 and 0.6502 for C58. All models were refined with GalaxyWEB and 3Drefine, and the refined structures were validated by ProSA-web Z-score, ERRAT overall quality factor and RAMPAGE favoured, allowed and outlier regions (Table 3; Supplementary Tables S4 and S5). On the basis of the combined results, Robetta model 3 of C58 and Robetta model 5 of TLR5 were selected for docking (Supplementary Figure S2).
HADDOCK. Complex I (TLR5–C58) had the lowest HADDOCK score (204.0 ± 27.6), followed by complex III (218.2 ± 14.1) and complex II (269.4 ± 12.4) (Table 4). Complex I also had the largest buried surface area (4542.46 Ų), compared with 4367.7 Ų for complex II and 4064.31 Ų for complex III, and the lowest restraint-violation energy (5026.5 ± 167.8, compared with 5378.9 ± 256.5 and 5368.8 ± 112.1). Its van der Waals energy (−157.4 ± 14.1 kcal mol⁻¹) was the most favourable, although the difference from complex III (−153.4 ± 14.0 kcal mol⁻¹) lies within the standard deviation. Electrostatic energies were similar among the three complexes (−572.5 ± 102.2, −578.7 ± 50.2 and −568.5 ± 102.3 kcal mol⁻¹). Desolvation energy was most favourable for complex III (−66.8), followed by complex I (−30.0) and complex II (−2.4). Taken together, C58 had the best overall HADDOCK profile and engaged TLR5 over the largest interface.
PRODIGY. PRODIGY gave the same overall ranking (Table 5). The TLR5–C58 complex had ΔG = −21.0 kcal mol⁻¹ (Kd 4.0 × 10⁻¹⁶ M), the TLR5–C16 complex −20.3 kcal mol⁻¹ (1.3 × 10⁻¹⁵ M) and the TLR5–C52 complex −17.4 kcal mol⁻¹ (1.7 × 10⁻¹³ M). The two docking schemes agree that C58 forms the strongest complex with TLR5 and that C52 forms the weakest.
Hydrogen bonds. LIGPLOT analysis identified 27 intermolecular hydrogen bonds in complex I and 23 each in complexes II and III (Table 6, Figure 1 and Figure 2). Several residues recurred in the hydrogen-bond networks of all three complexes, including Asp133, Lys137, Thr85, Asp87 and Ser159 in the first residue set of Table 6 and Arg62, Gln72, Thr73, Asn65 and Asn129 in the second. This points to a shared interface hotspot. Complex I also showed a dense cluster of bonds (rows 17–24 of Table 6) that includes Asp215, Lys218, Ser213 and Asn245, residues that do not form hydrogen bonds in the other two complexes. The larger number of hydrogen bonds is consistent with the better docking scores of C58. The docking results motivated us to test the stability of the TLR5–C58 complex on the nanosecond time scale by MD simulation.
4. Discussion
L2 is a promising basis for second-generation HPV vaccines because its N-terminal region is conserved and carries cross-neutralizing epitopes, but its low immunogenicity has so far prevented clinical use [6,11]. Combining L2 epitopes with a TLR agonist in one construct is one way to solve this problem, and a flagellin–L2 multimer protected animals against diverse papillomavirus genotypes [12,13]. In our constructs, the highly conserved D0/D1 regions of flagellin supply TLR5 agonism, RS09 adds TLR4 agonism [35], and PADRE and TpD supply universal T-helper help that is independent of HLA polymorphism [36,37]. Flexible and rigid linkers separate the functional units so that each can fold and be recognized independently [38].
The three FFPVs differ only in their L2 epitope segments, yet they bind TLR5 differently. Both docking schemes ranked the TLR5–C58 complex first, C16 second and C52 third, and the order was reflected in buried surface area and hydrogen-bond count. The agreement between HADDOCK and PRODIGY is reassuring because the two methods score the same docked models differently: HADDOCK uses an energy-based function, whereas PRODIGY uses counts of interfacial contacts. Differences between individual HADDOCK terms were, however, often within the standard deviation, and the PRODIGY ΔG values for C58 and C16 differ by only 0.7 kcal mol⁻¹. We therefore regard the ranking of C58 above C52 as clear and the ranking of C58 above C16 as suggestive, and we would not read much into the absolute values of the scores.
Hydrogen bonds are a major source of specificity in biological recognition, and the 27 hydrogen bonds in the TLR5–C58 interface (versus 23 in the other complexes) indicate a more extensive polar network. That several residues recur in all three interfaces suggests that the flagellin-derived portion of the vaccines engages TLR5 through a common surface and that the L2 epitope segments modulate the strength of the interaction. Flagellin binds the lateral face of the TLR5 LRR ectodomain [15], and the LRR region is also where the ten key residues identified by MD simulation are located. Together, these observations support the view that the constructs bind in a flagellin-like mode. Given that TLR5 expression increases during cervical neoplasia [16], a TLR5-targeted vaccine could also engage the receptor at the site of disease.
Several limitations should be noted. The vaccine and TLR5 structures are predicted models and not experimental structures, and the docking scores depend on model quality (Table 3). Docking and PRODIGY scores are ranking tools that are not substitutes for measured affinities. The MD simulation covers a single trajectory of one complex. Computational prediction of allergenicity and antigenicity should also be confirmed experimentally, since safety evaluation remains a central issue for vaccine components [60]. The next steps are to express the constructs, to measure binding to TLR5 (for example by surface plasmon resonance), to test TLR5 and NF-κB activation in reporter cells, and to assess immunogenicity and neutralizing activity in animals. The combination of computational prediction and experimental validation is increasingly the standard in data-driven vaccine design [61]. MD-based validation of docked complexes has been used in our earlier work on ligand–protein systems [62,63,64], and here we apply it to a protein–protein complex.
5. Conclusions
We designed three flagellin-fabricated chimeric L2 peptide vaccines (C58, C52 and C16), predicted them to be stable, soluble and antigenic, and compared their binding to TLR5. Two independent docking schemes and the hydrogen-bond analysis ranked the TLR5–C58 complex highest, and MD simulation of this complex identified a binding mode with ten key residues in the TLR5 LRR region. The study provides a structural hypothesis for how flagellin-fabricated L2 vaccines engage TLR5 and a basis for experimental testing of the C58 candidate against cervix papilloma.
Acknowledgments
The authors are grateful for the Key Research Area Grant 2016YFA0501703 from the Ministry of Science and Technology of China, Henan Natural Science Grant 162300410060, grants from the State Key Laboratory of Microbial Metabolism, and Joint Research Funds for Medical and Engineering and Scientific Research at Shanghai Jiao Tong University, supported by D.Q.W. G.S. and S.K. are grateful for Postdoctoral Research Grants 21450003 and 21450004 from Henan University of Technology, Henan Province, China. K.G. is supported by research grant 2013BAB11B02 from the Ministry of Science and Technology of China.
Competing interests
The authors declare that they have no competing interests.
References
- Bray, F.; Laversanne, M.; Sung, H.; Ferlay, J.; Siegel, R.L.; Soerjomataram, I.; et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin. 2024, 74(3), 229–263. [Google Scholar] [CrossRef] [PubMed]
- Zhou, H.L.; Zhang, W.; Zhang, C.J.; Wang, S.M.; Duan, Y.C.; Wang, J.X.; et al. Prevalence and distribution of human papillomavirus genotypes in Chinese women between 1991 and 2016: a systematic review. J. Infect. 2018, 76(6), 522–528. [Google Scholar] [CrossRef] [PubMed]
- Zheng, L.; Chen, S.; Yang, F.; Wang, W.; Xu, C.; Zheng, L. High-risk HPV prevalence and genotype distribution among women in Liaocheng, Shandong Province, China from 2016 to 2022. Front Public Health 2023, 11, 1145396. [Google Scholar] [CrossRef] [PubMed]
- Joura, E.A.; Giuliano, A.R.; Iversen, O.E.; Bouchard, C.; Mao, C.; Mehlsen, J.; et al. A 9-valent HPV vaccine against infection and intraepithelial neoplasia in women. N Engl. J. Med. 2015, 372(8), 711–723. [Google Scholar] [CrossRef] [PubMed]
- Monie, A.; Hung, C.F.; Roden, R.; Wu, T.C. Cervarix: a vaccine for the prevention of HPV 16, 18-associated cervical cancer. Biologics 2008, 2(1), 97–105. [Google Scholar] [CrossRef]
- Tyler, M.; Tumban, E.; Chackerian, B. Second-generation prophylactic HPV vaccines: successes and challenges. Expert Rev. Vaccines 2014, 13(2), 247–255. [Google Scholar] [CrossRef] [PubMed]
- Karanam, B.; Jagu, S.; Huh, W.K.; Roden, R.B. Developing vaccines against minor capsid antigen L2 to prevent papillomavirus infection. Immunol. Cell Biol. 2009, 87(4), 287–299. [Google Scholar] [CrossRef] [PubMed]
- Wang, J.W.; Roden, R.B. L2, the minor capsid protein of papillomavirus. Virology 2013, 445(1–2), 175–186. [Google Scholar] [CrossRef] [PubMed]
- Gambhira, R.; Karanam, B.; Jagu, S.; Roberts, J.N.; Buck, C.B.; Bossis, I.; et al. A protective and broadly cross-neutralizing epitope of human papillomavirus L2. J. Virol. 2007, 81(24), 13927–13931. [Google Scholar] [CrossRef] [PubMed]
- Schiller, J.T.; Müller, M. Next generation prophylactic human papillomavirus vaccines. Lancet Oncol. 2015, 16(5), e217–e225. [Google Scholar] [CrossRef] [PubMed]
- Jiang, R.T.; Schellenbacher, C.; Chackerian, B.; Roden, R.B. Progress and prospects for L2-based human papillomavirus vaccines. Expert Rev. Vaccines 2016, 15(7), 853–862. [Google Scholar] [CrossRef] [PubMed]
- Kalnin, K.; Tibbitts, T.; Yan, Y.; Stegalkina, S.; Shen, L.; Costa, V.; et al. Low doses of flagellin-L2 multimer vaccines protect against challenge with diverse papillomavirus genotypes. Vaccine 2014, 32(28), 3540–3547. [Google Scholar] [CrossRef] [PubMed]
- Kalnin, K.; Chivukula, S.; Tibbitts, T.; Yan, Y.; Stegalkina, S.; Shen, L.; et al. Incorporation of RG1 epitope concatemers into a self-adjuvanting flagellin-L2 vaccine broaden durable protection against cutaneous challenge with diverse human papillomavirus genotypes. Vaccine 2017, 35(37), 4942–4951. [Google Scholar] [CrossRef] [PubMed]
- Takeuchi, O.; Akira, S. Pattern recognition receptors and inflammation. Cell 2010, 140(6), 805–820. [Google Scholar] [CrossRef] [PubMed]
- Yoon, S.I.; Kurnasov, O.; Natarajan, V.; Hong, M.; Gudkov, A.V.; Osterman, A.L.; et al. Structural basis of TLR5-flagellin recognition and signaling. Science 2012, 335(6070), 859–864. [Google Scholar] [CrossRef] [PubMed]
- Kim, W.Y.; Lee, J.W.; Choi, J.J.; Choi, C.H.; Kim, T.J.; Kim, B.G.; et al. Increased expression of Toll-like receptor 5 during progression of cervical neoplasia. Int. J. Gynecol. Cancer 2008, 18(2), 300–305. [Google Scholar] [CrossRef] [PubMed]
- Forstnerič, V.; Ivičak-Kocjan, K.; Plaper, T.; Jerala, R.; Benčina, M. The role of the C-terminal D0 domain of flagellin in activation of Toll like receptor 5. PLoS Pathog. 2017, 13(8), e1006574. [Google Scholar] [CrossRef] [PubMed]
- Yang, J.; Zhong, M.; Zhang, Y.; Zhang, E.; Sun, Y.; Cao, Y.; et al. Antigen replacement of domains D2 and D3 in flagellin promotes mucosal IgA production and attenuates flagellin-induced inflammatory response after intranasal immunization. Hum. Vaccin Immunother. 2013, 9(5), 1084–1092. [Google Scholar] [CrossRef] [PubMed]
- Kaliamurthi, S.; Selvaraj, G.; Junaid, M.; Khan, A.; Gu, K.; Wei, D.Q. Cancer immunoinformatics: a promising era in the development of peptide vaccines for human papillomavirus-induced cervical cancer. Curr. Pharm. Des. 2018, 24(32), 3791–3817. [Google Scholar] [CrossRef] [PubMed]
- Kaliamurthi, S.; Selvaraj, G.; Chinnasamy, S.; Wang, Q.; Nangraj, A.S.; Cho, W.C.; et al. Exploring the papillomaviral proteome to identify potential candidates for a chimeric vaccine against cervix papilloma using immunomics and computational structural vaccinology. Viruses 2019, 11(1), 63. [Google Scholar] [CrossRef] [PubMed]
- Kaliamurthi, S.; Selvaraj, G.; Kaushik, A.C.; Gu, K.; Wei, D.Q. Designing of CD8+ and CD8+-overlapped CD4+ epitope vaccine by targeting late and early proteins of human papillomavirus. Biologics 2018, 12, 107–125. [Google Scholar] [CrossRef] [PubMed]
- Negahdaripour, M.; Eslami, M.; Nezafat, N.; Hajighahramani, N.; Ghoshoon, M.B.; Shoolian, E.; et al. A novel HPV prophylactic peptide vaccine, designed by immunoinformatics and structural vaccinology approaches. Infect. Genet Evol. 2017, 54, 402–416. [Google Scholar] [CrossRef] [PubMed]
- UniProt Consortium. UniProt: the universal protein knowledgebase. Nucleic Acids Res. 2017, 45(D1), D158–D169. [Google Scholar] [CrossRef] [PubMed]
- Moutaftsi, M.; Peters, B.; Pasquetto, V.; Tscharke, D.C.; Sidney, J.; Bui, H.H.; et al. A consensus epitope prediction approach identifies the breadth of murine T(CD8+)-cell responses to vaccinia virus. Nat. Biotechnol. 2006, 24(7), 817–819. [Google Scholar] [CrossRef] [PubMed]
- Andreatta, M.; Nielsen, M. Gapped sequence alignment using artificial neural networks: application to the MHC class I system. Bioinformatics 2016, 32(4), 511–517. [Google Scholar] [CrossRef] [PubMed]
- Wang, P.; Sidney, J.; Dow, C.; Mothé, B.; Sette, A.; Peters, B. A systematic assessment of MHC class II peptide binding predictions and evaluation of a consensus approach. PLoS Comput Biol. 2008, 4(4), e1000048. [Google Scholar] [CrossRef] [PubMed]
- Gonzalez-Galarza, F.F.; Takeshita, L.Y.; Santos, E.J.; Kempson, F.; Maia, M.H.; da Silva, A.L.; et al. Allele frequency net 2015 update: new features for HLA epitopes, KIR and disease and HLA adverse drug reaction associations. Nucleic Acids Res. 2015, 43(D1), D784–D788. [Google Scholar] [CrossRef] [PubMed]
- Bhasin, M.; Raghava, G.P.S. Prediction of CTL epitopes using QM, SVM and ANN techniques. Vaccine 2004, 22(23–24), 3195–3204. [Google Scholar] [CrossRef] [PubMed]
- Liu, I.H.; Lo, Y.S.; Yang, J.M. PAComplex: a web server to infer peptide antigen families and binding models from TCR–pMHC complexes. Nucleic Acids Res. 2011, 39, W254–W260, (Web Server issue). [Google Scholar] [CrossRef] [PubMed]
- Dhanda, S.K.; Vir, P.; Raghava, G.P.S. Designing of interferon-gamma inducing MHC class-II binders. Biol. Direct 2013, 8, 30. [Google Scholar] [CrossRef] [PubMed]
- Saha, S.; Raghava, G.P.S. Prediction of continuous B-cell epitopes in an antigen using recurrent neural network. Proteins 2006, 65(1), 40–48. [Google Scholar] [CrossRef] [PubMed]
- Dimitrov, I.; Bangov, I.; Flower, D.R.; Doytchinova, I. AllerTOP v.2—a server for in silico prediction of allergens. J. Mol. Model. 2014, 20(6), 2278. [Google Scholar] [CrossRef] [PubMed]
- Dimitrov, I.; Naneva, L.; Doytchinova, I.; Bangov, I. AllergenFP: allergenicity prediction by descriptor fingerprints. Bioinformatics 2014, 30(6), 846–851. [Google Scholar] [CrossRef] [PubMed]
- Magnan, C.N.; Zeller, M.; Kayala, M.A.; Vigil, A.; Randall, A.; Felgner, P.L.; et al. High-throughput prediction of protein antigenicity using protein microarray data. Bioinformatics 2010, 26(23), 2936–2943. [Google Scholar] [CrossRef] [PubMed]
- Shanmugam, A.; Rajoria, S.; George, A.L.; Mittelman, A.; Suriano, R.; Tiwari, R.K. Synthetic Toll like receptor-4 (TLR-4) agonist peptides as a novel class of adjuvants. PLoS ONE 2012, 7(2), e30839. [Google Scholar] [CrossRef] [PubMed]
- Alexander, J.; Sidney, J.; Southwood, S.; Ruppert, J.; Oseroff, C.; Maewal, A.; et al. Development of high potency universal DR-restricted helper epitopes by modification of high affinity DR-blocking peptides. Immunity 1994, 1(9), 751–761. [Google Scholar] [CrossRef] [PubMed]
- Fraser, C.C.; Altreuter, D.H.; Ilyinskii, P.; Pittet, L.; LaMothe, R.A.; Keegan, M.; et al. Generation of a universal CD4 memory T cell recall peptide effective in humans, mice and non-human primates. Vaccine 2014, 32(24), 2896–2903. [Google Scholar] [CrossRef] [PubMed]
- Chen, X.; Zaro, J.L.; Shen, W.C. Fusion protein linkers: property, design and functionality. Adv. Drug Deliv. Rev. 2013, 65(10), 1357–1369. [Google Scholar] [CrossRef] [PubMed]
- Gasteiger, E.; Hoogland, C.; Gattiker, A.; Duvaud, S.; Wilkins, M.R.; Appel, R.D.; et al. Protein identification and analysis tools on the ExPASy server. In The Proteomics Protocols Handbook; Walker, J.M., Ed.; Humana Press: Totowa (NJ), 2005; pp. 571–607. [Google Scholar]
- Magnan, C.N.; Randall, A.; Baldi, P. SOLpro: accurate sequence-based prediction of protein solubility. Bioinformatics 2009, 25(17), 2200–2207. [Google Scholar] [CrossRef] [PubMed]
- Doytchinova, I.A.; Flower, D.R. VaxiJen: a server for prediction of protective antigens, tumour antigens and subunit vaccines. BMC Bioinform. 2007, 8, 4. [Google Scholar] [CrossRef] [PubMed]
- Kim, D.E.; Chivian, D.; Baker, D. Protein structure prediction and analysis using the Robetta server. Nucleic Acids Res. 2004, 32, W526–W531, (Web Server issue). [Google Scholar] [CrossRef] [PubMed]
- Ko, J.; Park, H.; Heo, L.; Seok, C. GalaxyWEB server for protein structure prediction and refinement. Nucleic Acids Res. 2012, 40, W294–W297, (Web Server issue). [Google Scholar] [CrossRef] [PubMed]
- Bhattacharya, D.; Cheng, J. i3Drefine software for protein 3D structure refinement and its assessment in CASP10. PLoS ONE 2013, 8(7), e69648. [Google Scholar] [CrossRef] [PubMed]
- Wiederstein, M.; Sippl, M.J. ProSA-web: interactive web service for the recognition of errors in three-dimensional structures of proteins. Nucleic Acids Res. 2007, 35, W407–W410, (Web Server issue). [Google Scholar] [CrossRef] [PubMed]
- Lovell, S.C.; Davis, I.W.; Arendall, W.B., 3rd; de Bakker, P.I.W.; Word, J.M.; Prisant, M.G.; et al. Structure validation by Cα geometry: φ, ψ and Cβ deviation. Proteins 2003, 50(3), 437–450. [Google Scholar] [CrossRef] [PubMed]
- Colovos, C.; Yeates, T.O. Verification of protein structures: patterns of nonbonded atomic interactions. Protein Sci. 1993, 2(9), 1511–1519. [Google Scholar] [CrossRef] [PubMed]
- Pettersen, E.F.; Goddard, T.D.; Huang, C.C.; Couch, G.S.; Greenblatt, D.M.; Meng, E.C.; et al. UCSF Chimera—a visualization system for exploratory research and analysis. J. Comput Chem. 2004, 25(13), 1605–1612. [Google Scholar] [CrossRef] [PubMed]
- van Zundert, G.C.P.; Rodrigues, J.P.G.L.M.; Trellet, M.; Schmitz, C.; Kastritis, P.L.; Karaca, E.; et al. The HADDOCK2.2 web server: user-friendly integrative modeling of biomolecular complexes. J. Mol. Biol. 2016, 428(4), 720–725. [Google Scholar] [CrossRef] [PubMed]
- 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 2016, 32(23), 3676–3678. [Google Scholar] [CrossRef] [PubMed]
- Wallace, A.C.; Laskowski, R.A.; Thornton, J.M. LIGPLOT: a program to generate schematic diagrams of protein–ligand interactions. Protein Eng. 1995, 8(2), 127–134. [Google Scholar] [CrossRef] [PubMed]
- Abraham, M.J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J.C.; Hess, B.; et al. GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1–2, 19–25. [Google Scholar] [CrossRef]
- Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; de Groot, B.L.; et al. CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nat. Methods 2017, 14(1), 71–73. [Google Scholar] [CrossRef] [PubMed]
- Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79(2), 926–935. [Google Scholar] [CrossRef]
- Essmann, U.; Perera, L.; Berkowitz, M.L.; Darden, T.; Lee, H.; Pedersen, L.G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103(19), 8577–8593. [Google Scholar] [CrossRef]
- Hess, B.; Bekker, H.; Berendsen, H.J.C.; Fraaije, J.G.E.M. LINCS: a linear constraint solver for molecular simulations. J. Comput Chem. 1997, 18(12), 1463–1472. [Google Scholar]
- Miyamoto, S.; Kollman, P.A. SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput Chem. 1992, 13(8), 952–962. [Google Scholar] [CrossRef]
- Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126(1), 014101. [Google Scholar] [CrossRef] [PubMed]
- Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: a new molecular dynamics method. J. Appl. Phys. 1981, 52(12), 7182–7190. [Google Scholar] [CrossRef]
- Selvaraj, G.; Kaliamurthi, S.; Peslherbe, G.H.; Wei, D.Q. Are the allergic reactions of COVID-19 vaccines caused by mRNA constructs or nanocarriers? An immunological insight. Interdiscip. Sci. 2021, 13(2), 344–347. [Google Scholar] [CrossRef] [PubMed]
- Selvaraj, G.; Kaliamurthi, S.; Wei, D.Q. Editorial: Data-driven vaccine design for microbial-associated diseases. Front Immunol. 2025, 16, 1749184. [Google Scholar] [CrossRef] [PubMed]
- Kaliamurthi, S.; Selvaraj, G.; Selvaraj, C.; Singh, S.K.; Wei, D.Q.; Peslherbe, G.H. Structure-based virtual screening reveals ibrutinib and zanubrutinib as potential repurposed drugs against COVID-19. Int. J. Mol. Sci. 2021, 22(13), 7071. [Google Scholar] [CrossRef] [PubMed]
- Selvaraj, C.; Selvaraj, G.; Ismail, R.M.; Khan, J.; Baazeem, A.; Wei, D.Q.; et al. Interrogation of Bacillus anthracis SrtA active site loop forming open/close lid conformations through extensive MD simulations for understanding binding selectivity of SrtA inhibitors. Saudi J. Biol. Sci. 2021, 28(7), 3650–3659. [Google Scholar] [CrossRef] [PubMed]
- Chinnasamy, S.; Selvaraj, G.; Kaushik, A.C.; Kaliamurthi, S.; Chandrabose, S.; Singh, S.K.; et al. Molecular docking and molecular dynamics simulation studies to identify potent AURKA inhibitors: assessing the performance of density functional theory, MM-GBSA and mass action kinetics calculations. J. Biomol. Struct. Dyn. 2020, 38(14), 4325–4335. [Google Scholar] [CrossRef] [PubMed]
Table 1.
Overlapped epitope segments from the N-terminal region of the L2 protein of HPV58 (C58), HPV52 (C52) and HPV16 (C16) retained for vaccine construction.
Table 1.
Overlapped epitope segments from the N-terminal region of the L2 protein of HPV58 (C58), HPV52 (C52) and HPV16 (C16) retained for vaccine construction.
| Candidate | Start | End | Overlapped epitope segment | IFN-γ-inducing epitope (score) | AllergenFP | AllerTOP | ANTIGENpro |
|---|---|---|---|---|---|---|---|
| C58 [20] | 23 | 36 | KVEGTTIADQILRY | 0.528 | Non-allergen | Non-allergen | Antigen |
| C58 | 29 | 42 | IADQILRYGSLGVF | 1.000 | Non-allergen | Non-allergen | Antigen |
| C52 | 23 | 36 | KVEGTTIADQLLKY | 0.389 | Non-allergen | Non-allergen | Antigen |
| C52 | 29 | 42 | IADQLLKYGSLGVF | 0.505 | Non-allergen | Non-allergen | Antigen |
| C16 [22] | 10 | 36 | TKRASATQLYKTCKQAGTCPPDIIPKV | 3.21 | Non-allergen | Non-allergen | Antigen |
| C16 | 65 | 89 | TGGRTGYIPLGTRPPTATDTLAPV | 2.73 | Non-allergen | Non-allergen | Antigen |
IFN-γ production is indicated by positive scores. Cross-reactivity of the selected epitopes with the human proteome was zero.
Table 2.
Physicochemical and immunological properties of the designed vaccine constructs.
| Property | C58 | C52 | C16 |
|---|---|---|---|
| Number of amino acids | 318 | 318 | 341 |
| Molecular weight (Da) | 33394.15 | 33338.13 | 35567.69 |
| Theoretical pI | 8.00 | 7.98 | 9.28 |
| Negatively charged residues (Asp + Glu) | 24 | 24 | 23 |
| Positively charged residues (Arg + Lys) | 25 | 25 | 29 |
| Extinction coefficient (M⁻¹ cm⁻¹) | 12950 | 12950 | 12950 |
| Estimated half-life | 20 h (mammals), 30 min (yeast), >10 h (E. coli) | 20 h (mammals), 30 min (yeast), >10 h (E. coli) | 20 h (mammals), 30 min (yeast), >10 h (E. coli) |
| Instability index | 35.88 (stable) | 32.38 (stable) | 35.47 (stable) |
| Aliphatic index | 94.62 | 94.62 | 87.98 |
| GRAVY | −0.190 | −0.190 | −0.257 |
| Solubility probability (SOLpro) | 0.621 | 0.587 | 0.701 |
| Antigenicity (VaxiJen) | 0.4301 (antigen) | 0.4468 (antigen) | 0.4870 (antigen) |
| Antigenicity (ANTIGENpro) | 0.943820 (antigen) | 0.949012 (antigen) | 0.948248 (antigen) |
GRAVY, grand average of hydropathicity.
Table 3.
Structure validation of the refined models of the vaccine constructs (ProSA-web, RAMPAGE and additional stereochemical checks). Only the values for C52 model 1 were available in the files supplied; "–" marks values to be completed from the validation output.
Table 3.
Structure validation of the refined models of the vaccine constructs (ProSA-web, RAMPAGE and additional stereochemical checks). Only the values for C52 model 1 were available in the files supplied; "–" marks values to be completed from the validation output.
| Model | ProSA Z-score | Favoured (%) | Allowed (%) | Outlier (%) | Chi1–chi2 better | Chi1–chi2 inside | Chi1–chi2 worse | Bond length/angle | Bad contacts | Dihedrals | Covalent | Planar within limits (%) | Planar highlighted (%) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| C52 model 1 | – | 90.5 | 8.1 | 1.1 | 5 | 0 | 0 | 4.8 | 0 | 0.43 | 0.40 | 100 | 0 |
| C52 model 2 | – | – | – | – | – | – | – | – | – | – | – | – | – |
| C52 model 3 | – | – | – | – | – | – | – | – | – | – | – | – | – |
| C52 model 4 | – | – | – | – | – | – | – | – | – | – | – | – | – |
| C52 model 5 | – | – | – | – | – | – | – | – | – | – | – | – | – |
| C16 model 1 | – | – | – | – | – | – | – | – | – | – | – | – | – |
| C16 model 2 | – | – | – | – | – | – | – | – | – | – | – | – | – |
| C16 model 3 | – | – | – | – | – | – | – | – | – | – | – | – | – |
| C16 model 4 | – | – | – | – | – | – | – | – | – | – | – | – | – |
| C16 model 5 | – | – | – | – | – | – | – | – | – | – | – | – | – |
Table 4.
HADDOCK results for the docking of the vaccine constructs to TLR5.
| Complex | HADDOCK score | Binding energy (kcal mol⁻¹) | van der Waals energy (kcal mol⁻¹) | Electrostatic energy (kcal mol⁻¹) | Desolvation energy (kcal mol⁻¹) | Restraint-violation energy (kcal mol⁻¹) | Buried surface area (Ų) | Hydrogen bonds |
|---|---|---|---|---|---|---|---|---|
| TLR5–C58 | 204.0 ± 27.6 | −36341.2 | −157.4 ± 14.1 | −572.5 ± 102.2 | −30.0145 | 5026.5 ± 167.8 | 4542.46 | 27 |
| TLR5–C52 | 269.4 ± 12.4 | −33712.6 | −142.4 ± 4.2 | −578.7 ± 50.2 | −2.36732 | 5378.9 ± 256.5 | 4367.7 | 23 |
| TLR5–C16 | 218.2 ± 14.1 | −33604.5 | −153.4 ± 14.0 | −568.5 ± 102.3 | −66.8167 | 5368.8 ± 112.1 | 4064.31 | 23 |
Table 5.
Binding affinity of the TLR5–vaccine complexes predicted by PRODIGY from intermolecular contacts.
Table 5.
Binding affinity of the TLR5–vaccine complexes predicted by PRODIGY from intermolecular contacts.
| Complex | ΔG (kcal mol⁻¹) | Kd (M) at 25 °C |
|---|---|---|
| TLR5–C58 | −21.0 | 4.0 × 10⁻¹⁶ |
| TLR5–C52 | −17.4 | 1.7 × 10⁻¹³ |
| TLR5–C16 | −20.3 | 1.3 × 10⁻¹⁵ |
ΔG, binding affinity; Kd, dissociation constant.
Table 6.
Intermolecular hydrogen bonds, interacting residue pairs and distances (Å) in the three complexes. Each pair lists the two residues (and atoms) forming the bond.
Table 6.
Intermolecular hydrogen bonds, interacting residue pairs and distances (Å) in the three complexes. Each pair lists the two residues (and atoms) forming the bond.
| No. | TLR5–C58: residue pair | Distance | TLR5–C52: residue pair | Distance | TLR5–C16: residue pair | Distance |
|---|---|---|---|---|---|---|
| 1 | Glu191:OE1 – Arg62:NH2 | 3.25 | Lys275:NZ – Tyr283:O | 2.88 | Asp133:OD1 – Gln59:HE2 | 3.29 |
| 2 | His188:O – Arg62:NH2 | 2.88 | Asp241:OD1 – Ser276:OG | 3.07 | Asp133:OD2 – Arg62:NH1 | 2.77 |
| 3 | Asp133:OD1 – Thr73:OG1 | 2.69 | Arg59:NH1 – Asn83:ND2 | 3.23 | Ser132:OG – Arg62:NH1 | 3.11 |
| 4 | Ser132:O – Gln72:HE21 | 2.85 | Thr60:OG1 – Asn83:ND2 | 2.97 | Ile107:O – Arg62:NH2 | 3.33 |
| 5 | Lys137:O – Asn129:ND2 | 2.68 | Thr62:OG1 – Glu80:OE1 | 2.89 | Thr85:OG1 – Lys55:NZ | 2.80 |
| 6 | Tyr108:OH – Glu75:OE2 | 2.65 | Thr85:OG1 – Gly76:O | 2.93 | His111:NE1 – Asn295:OD1 | 2.80 |
| 7 | Thr62:OG1 – Glu118:OE2 | 2.93 | Lys137:O – Asn129:ND2 | 2.87 | Lys88:NZ – Glu302:OE2 | 2.69 |
| 8 | Thr85:OG1 – Gly76:O | 2.55 | Lys137:NZ – Asp66:OD2 | 2.67 | Asp87:OD1 – Arg298:NH1 | 2.94 |
| 9 | Thr85:OG1 – Glu80:N | 2.64 | Val135:N – Ser69:OG | 2.97 | Asp87:OD2 – Arg298:NH2 | 2.58 |
| 10 | Val61:N – Glu80:OE1 | 2.65 | Asp133:OD1 – Thr73:OG1 | 2.77 | Ser159:OG – Asp66:OD2 | 2.71 |
| 11 | His111:NE2 – Gln125:O | 2.82 | Gly130:O – Gln72:NE2 | 3.22 | Leu136:NZ – Asp66:OD1 | 2.79 |
| 12 | Asp87:OD1 – Gln125:HE22 | 2.70 | Asp133:OD2 – Gln72:NH | 2.77 | Lys137:O – Asn65:ND2 | 2.74 |
| 13 | Arg287:NE – Ser29:O | 2.82 | Ile107:O – Gln72:OE1 | 2.96 | Lys137:NZ – Ser69:CB | 2.64 |
| 14 | Ser159:NH2 – Asn65:O | 2.63 | Lys106:NZ – Ser69:OG | 2.99 | Arg158:NE – Asn129:OD1 | 2.79 |
| 15 | Ser159:NH2 – Ser30:OG | 2.65 | Ala134:N – Arg121:NH2 | 2.74 | Arg158:NH2 – Asn129:OD1 | 3.09 |
| 16 | Asn245:ND2 – Asp180:OD1 | 3.31 | Asp87:OD2 – Arg121:NH2 | 2.60 | Glu189:OE1 – Thr73:OG1 | 2.57 |
| 17 | Trp216:NE1 – Asp280:OD2 | 2.77 | Asp87:OD1 – Gln125:NE1 | 2.58 | Cys186:SG – Thr73:O | 2.98 |
| 18 | Trp216:O – Arg275:NH2 | 2.98 | Glu89:OE1 – Asp280:NE2 | 3.10 | Val214:O – Arg121:NH2 | 2.87 |
| 19 | Val214:O – Arg275:NH2 | 3.02 | Gln156:O – Asn65:OD1 | 3.20 | Glu187:OE1 – Asn79:ND2 | 3.11 |
| 20 | Val214:O – Arg275:NH1 | 2.85 | Arg158:N – Asn65:OD1 | 2.85 | His188:O – Gln72:NE2 | 2.81 |
| 21 | Asp215:OD1 – Thr51:OG1 | 2.60 | Ser159:OG – Asn65:OD2 | 2.78 | His188:ND1 – Glu75:OE1 | 3.13 |
| 22 | Lys218:NL – Thr51:OG7 | 2.68 | Arg158:NH2 – Asp258:OD1 | 2.67 | Glu191:OE2 – Gln72:HE21 | 2.69 |
| 23 | Lys218:NL – Glu279:OE2 | 2.63 | Arg158:NH1 – Arg62:NE | 2.65 | Trp216:N – Glu80:OE1 | 2.92 |
| 24 | Lys218:NL – Ala47:O | 3.26 | NA | NA | NA | NA |
| 25 | Glu187:OE2 – Lys55:NZ | 2.63 | NA | NA | NA | NA |
| 26 | Cys186:SG – Lys55:NZ | 3.23 | NA | NA | NA | NA |
| 27 | Ser213:N – Asn272:OD1 | 2.85 | NA | NA | NA | NA |
NA, not applicable.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. |
© 2026 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.