A computational approach to design a multiepitope vaccine against H5N1 virus

Abstract

Since 1997, highly pathogenic avian influenza viruses, such as H5N1, have been recognized as a possible pandemic hazard to men and the poultry business. The rapid rate of mutation of H5N1 viruses makes the whole process of designing vaccines extremely challenging. Here, we used an in silico approach to design a multi-epitope vaccine against H5N1 influenza A virus using hemagglutinin (HA) and neuraminidase (NA) antigens. B-cell epitopes, Cytotoxic T lymphocyte (CTL) and Helper T lymphocyte (HTL) were predicted via IEDB, NetMHC-4 and NetMHCII-2.3 respectively. Two adjuvants consisting of Human β-defensin-3 (HβD-3) along with pan HLA DR-binding epitope (PADRE) have been chosen to induce more immune response. Linkers including KK, AAY, HEYGAEALERAG, GPGPGPG and double EAAAK were utilized to link epitopes and adjuvants. This construct encodes a protein having 350 amino acids and 38.46 kDa molecular weight. Antigenicity of ~ 1, the allergenicity of non-allergen, toxicity of negative and solubility of appropriate were confirmed through Vaxigen, AllerTOP, ToxDL and DeepSoluE, respectively. The 3D structure of H5N1 was refined and validated with a Z-Score of − 0.87 and an overall Ramachandran of 99.7%. Docking analysis showed H5N1 could interact with TLR7 (docking score of − 374.08 and by 4 hydrogen bonds) and TLR8 (docking score of − 414.39 and by 3 hydrogen bonds). Molecular dynamics simulations results showed RMSD and RMSF of 0.25 nm and 0.2 for H5N1-TLR7 as well as RMSD and RMSF of 0.45 nm and 0.4 for H5N1-TLR8 complexes, respectively. Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) confirmed stability and continuity of interaction between H5N1-TLR7 with the total binding energy of − 29.97 kJ/mol and H5N1-TLR8 with the total binding energy of − 23.9 kJ/mol. Investigating immune response simulation predicted evidence of the ability to stimulate T and B cells of the immunity system that shows the merits of this H5N1 vaccine proposed candidate for clinical trials.

Introduction

A member of the Orthomyxoviridae family is influenza virus A [1]. This virus is responsible for causing acute respiratory disease and consistently changes while circulating among various animal hosts, such as wild birds, pigs, poultry, and horses as well as humans. The differentiation of these viruses is primarily according to the surface glycoproteins, specifically Hemagglutinin (HA) along with Neuraminidase (NA) [2]. There exist eighteen distinct variations of HA as well as eleven different NAs, each of which can be distinguished serologically. Notably, antibodies produced against a virus genotype don’t exhibit a reaction through other genotypes [3].

H5 influenza viruses have been typically found in avian populations and are naturally hosted by birds. However, they can infect a range of mammal species as well. There are nine recognized subtypes of H5 viruses, which include H5N1, H5N3, H5N2, H5N5, H5N4, H5N7, H5N6, H5N9 and H5N8. Among these subtypes, Worldwide, the viruses known as Low Pathogenic Avian Influenza (LPAI) are the most often seen in poultry and wild birds. Nevertheless, Highly Pathogenic Avian Influenza (HPAI) viruses are sporadically identified as well. The very first known instance of a highly pathogenic H5N1 virus infection in humans was reported by Yuen et al. 1998 [4]. These H5N1 viruses have significantly spread across regions in Asia, and the world. They have undergone rapid evolution, giving rise to ten distinct clades, labeled as clades 0 through 9. Since 2007, H5N1 viruses have been identified to be a significant hazard to men, wild birds, as well as poultry. Specifically, clade 2.2.2 was prevalent from 2007 to 2011, and from 2011 onwards, clade 2.3.2.1a has been in circulation [5].

The genotypes of H5N1 are determined by eight single-stranded RNA segments that encode a total of 11 proteins, among them surface glycoproteins known as neuraminidase (NA) and hemagglutinin (HA) [6]. Both HA and NA play essential roles that revolve around their contact with sialic acid. A terminal structure seen on glycoproteins or glycolipids that are expressed on the surface of cells is sialic acid [7]. The mechanism known as clathrin-mediated endocytosis, which is facilitated by HA's attachment to sialic acids on cellular receptors, allows the virus to enter cells. It's worth noting that the virus may also employ alternative endocytic pathways, like micropinocytosis, for cell entry [8].

NA is important in the final phases of infection. In particular, it functions to eliminate sialic acids from both cellular receptors and the freshly produced NA and HA present upon developing virions. The glycosylation activities of the host cell include the addition of these sialic acids [9, 10]. Moreover, NA's cleavage of sialic acids inhibits virions from clumping together and from attaching themselves to dying host cells via HA. This, in turn, facilitates the efficient release of virions and their spread to new target cells [9]. Given their crucial roles in the virus's pathogenicity, In the realm of vaccine research and development, HA and NA have become the main targets to neutralize avian influenza viruses (AIV) [10,11,12].

At present, there isn’t a vaccine strain on the market that can successfully produce protective immunity against the various clades and subclades of highly pathogenic avian influenza viruses (HPAIV) in the H5 lineage [13]. Nonetheless, Previously, research has indicated that vaccination against seasonal influenza A strains like H1N1 and H3N2 may improve cross-type cellular and humoral protection against the highly pathogenic avian influenza virus subtype H5N1 [14]. Both our research and the previously mentioned findings underscore the significant public health threat posed by H5N1 viruses, even though the number of human cases may be relatively small.

Consequently, the development of vaccines holds promise in helping to mitigate the impact of disease outbreaks. Reverse vaccinology is a recent breakthrough that, in addition to the conventional and laborious vaccine creation approaches, integrates immunogenomics and bioinformatics to quickly generate unique multi-epitope-based subunit vaccines. This approach offers the potential to expedite the vaccine development process [15, 16]. Today, there are extensive studies based on the use of bioinformatics and immunoinformatics tools to design new-generation vaccines [17,18,19,20,21,22,23].

Therefore, the current research aimed to produce a chimeric vaccine versus avian influenza A virus subtype H5N1 that is non-allergenic and immunogenic. This was achieved by utilizing a vaccinia’s approach, and It is advised that model animals be used in a laboratory environment to carry out the experimental validation of the suggested vaccination candidate.

Material and method

Antigen retrieval

The Influenza NCBI Research Database included the HA and NA protein sequences from several influenza H5N1 virus strains as described by Squires et al. 2012 [24].

Liner B cell epitopes and discontinues B cell epitopes prediction

To identify B cell antigenicity, Researchers used the Immune Epitope Database (IEDB) and BepiPred-2.0, a tool available at http://www.cbs.dtu.dk/services/BepiPred/, to detect linear B cell epitopes within the sequences of HA and NA proteins to predict B cell epitopes with 75% accuracy, 0.49 sensitivity, and 0.75 specificity, as detailed in the studies by Saadi et al. (2017) and Saha and Raghava [25](2006a). These tools incorporated several algorithms such as Emini surface accessibility prediction, and Bepipred linear epitope prediction 2.0 as well as Kolaskar and Tongaonkar antigenicity scale, which can be accessed at http://tools.iedb.org/bcell/. The semi-empirical Kolaskar and Tongaonkar technique takes into account the physicochemical characteristics of amino acid residues and their frequency in experimentally established segmental epitopes to predict antigenic determinants on proteins [26].

We also used the ElliPro server to forecast the discontinuous B cell epitopes' three-dimensional shape, which can be accessed at http://tools.iedb.org/ellipro/. ElliPro is a program that takes the three-dimensional structure of the protein into account to help anticipate conformational or discontinuous B cell epitopes. This is crucial for understanding how these epitopes are identified via the immunity system as well as for designing effective vaccines or therapeutic antibodies.

MHC-I epitope prediction

The T cell epitopes’ prediction, specifically cytotoxic T cell epitopes, has been carried out utilizing an enhanced neural network methodology offered via the NetMHC-4.0 server, accessible at https://services.healthtech.dtu.dk/service.php?NetMHC-4.0. In this analysis, all epitopes have been regarded to be 9-mers in length, as recommended by the server. These 9-mer epitopes are recognized by HLA Class I molecules, including HLA-A, HLA-C as well as HLA-B.

MHC-II epitope prediction

The epitope binding to MHC-II molecules has been predicted using the NetMHCII 2.3 server. This service predicts epitope affinity for a variety of HLA-II molecules, such as HLA-DR, HLA-DQ, and HLA-DP, using an artificial neural network. By using this prediction tool, Researchers can evaluate the probability that epitopes will attach to MHC-II molecules, which is critical for understanding and designing immune responses, particularly those involving helper T cells.

Population coverage determination

The Population Coverage instrument in the IEDB server, accessible at http://tools.iedb.org/population, was used to determine the population's coverage rate for particular epitopes. Determining the population coverage is essential to forecast the right epitopes for various HLA-binding interactions. Epitopes that exhibit strong binding to a variety of HLA alleles can provide broader coverage across different ethnicities, reducing the constraints of MHC limitations on T cell reactions. In such context, the evaluation of MHC-I HLA binding included a range of HLA alleles such as HLA-B5701, HLA-B5801. HLA-B5802, HLA-B2705, HLA-A0203, HLA-A2603, HLA-B8301, HLA-B4002, HLA-A3002, HLA-A2902, HLA-B4402, HLA-C0401, and HLA-C0303 for the NA and HA proteins. Additionally, MHC-II HLA binding assessments were performed for DRB1_0101, DPA10201-DPB10501, DPA10103-DPB10301, DQA10501-DQB10201, DRB1_0403, DPA10103-DPB10401, DRB1_0404, DRB1_0403, DRB1_0402, DQA10102-DQB10501, DRB1_0401, DRB1_1302, DQA10301-DQB10302 and DQA10401-DQB10402 for the NA and HA proteins. This comprehensive analysis helps in identifying epitopes that have the potential to generate immunity reactions across a broad range of human populations.

Analysis of epitopes conservation

The degree of conservation of particular epitopes was evaluated using the IEDB Conservancy Analysis tool as a monitoring tool about other sequences that were comparable. Epitopes can exhibit varying levels of conservation, ranging from a minimum of 0 percent to a maximum of 100 percent, depending on their specific length and sequence. This analysis is important in understanding how well-conserved an epitope is within a given protein or among related sequences. Highly conserved epitopes are more likely to be effective in stimulating immune responses because they are present in a broader range of strains or related proteins, making them potentially useful for vaccine design or therapeutic development. On the other hand, less conserved epitopes may be strain-specific and may have limited utility in generating immune responses across different variants or strains.

Physicochemical, antigenicity, allergenicity, toxicity properties of the vaccine

In the design of the vaccine construct, various linkers were strategically employed to connect different components and epitopes. These linkers serve as connectors between different functional elements and help optimize the vaccine’s overall design. AAY linkers, specific HEYGAEALERAG linkerGPGPG linkers, KK linkers and EAAAK linkers are used to bond Bcell, CLT and HLT epitopes. Also,HβD-3 and PADRE are used as adjuvants in the N-terminal of the construct.

The ProtParam instrument, available at https://web.expasy.org/protparam/, has been utilized for predicting the physicochemical characteristics of the final protein construct. It provides information about various protein parameters, such as molecular weight, and isoelectric point (pI), as well as amino acid composition.

The VaxiJen v2.0 server, accessible at http://www.ddg-pharmfac.net/vaxijen/VaxiJen.html, has been utilized for predicting the antigenic properties of the protein construct. It's considered an independent alignment method and is employed for antigen prediction. The server's accuracy typically ranges from 70 to 89%. AllerTOP v. 2.0, found at https://www.ddg-pharmfac.net/AllerTOP/, has been utilized for calculating the total protein sensitivity. This server uses an auto cross-covariance (ACC) as well as a k-nearest neighbor algorithm (kNN, k = 1) for predicting the proteins’ allergenicity. It considers various factors, including hydrophobicity, molecular weight, secondary structure characteristics, as well as the amino acids’ relative abundance. The ToxDL server, available at http://www.csbio.sjtu.edu.cn/bioinf/ToxDL/index.html, has been used to assess the protein's toxicity. This server employs an interpretable deep learning-based methodology for categorizing proteins into 2 categories: toxic and non-toxic. It uses a multi-modal approach like Convolutional Neural Networks (CNNs), and the InterProscan database. These tools and servers collectively provide valuable insights into physicochemical characteristics, antigenicity, allergenicity, along toxicity of protein construct, which are essential considerations in vaccine development and protein engineering.

The use of NetSurfP-2.0 online software, available at https://services.healthtech.dtu.dk/service.php?NetSurfP-2.0, has been implemented to predict the secondary structure of the protein. NetSurfP-2.0 is a sequence-based instrument that offers predictions regarding the secondary structure of amino acids within a protein.

3D modeling, refinement and validation of the H5N1 construct

The RoseTTAFold online software, accessible at https://robetta.bakerlab.org/, has been employed for constructing the three-dimensional (3D) structure of the protein. RoseTTAFold is considered a sophisticated “three-track” neural network, which means it regards multiple aspects simultaneously for predicting protein structures that have high accuracy. RoseTTAFold considers patterns in protein sequences, allowing it to understand the amino acid sequences that make up the protein.

To guarantee the precision and excellence of the protein structure simulation, PROCHECK and SAVES online servers, available at http://services.mbi.ucla.edu/PROCHECK and http://services.mbi.ucla.edu/SAVES, respectively, were used for structural validation.

The PROSA server, accessible at https://prosa.services.came.sbg.ac.at/prosa.php, has been utilized to identify the protein energy and Z-Score balance. This helps assess the quality and stability of the 3D structure of the protein.

The Ramachandran plot, which depicts the phi as well as psi angles of amino acids in protein structure, was generated using the online website http://molprobity.biochem.duke.edu/. It provides insight into the quality of the protein's backbone torsion angles.

The optimum spatial resolution for the study's optimal energy was found using the PyMol program. It is likely used for visualizing and analyzing the 3D protein structure. The website http://galaxy.seoklab.org/cgi-bin/submit.cgi?type=REFINE has been employed for protein structure refinement.

Docking analysis

The 3D structures of TLR7 (PDB code: 7CYN) and TLR8 (PDB code: 3w3g), were obtained from the RCSB server. The TLRs’ 3D structures have been adjusted for energy using the PyMOL v2.3.4 program. The vaccination and TLR PDB structures were cleared of water molecules. The prepared 3D structures were further processed within the Chimera V 1.13.1 program. H5N1 Influenza construct, TLR7 and TLR8 have been then submitted to the HDOCK server to identify interaction regions. Taking into account the complete molecular interface between H5N1-TLRs, the complexes with the lowest intermolecular binding energy have been ranked best. The choice was made with the lowest mean Root Mean Square Deviation (RMSD). H-bond formation has been evaluated using LigPlot⁺. These steps describe the process of retrieving, preparing, and analyzing the 3D structures of the components and their interactions to comprehend the binding and molecular contacts among the H5N1 influenza construct and TLRs.

Molecular dynamics and MM/PBSA analysis

Molecular dynamics simulations are a powerful tool for understanding the dynamic behavior of biomolecular complexes and are commonly used to investigate the interactions and stability of such complexes over time. The use of Gromacs and the CHARMM36 force field in this context helps provide valuable insights into the function of vaccine complexes in a simulated environment. Molecular dynamics calculations for the vaccine-m826 and vaccine-TLR7/8 complexes have been carried out using Gromacs v4.6.5 and the CHARMM36 all-atom force field. The protein–protein complexes have been solvated by adding water molecules, as well as the entire system has been then neutralized to ensure that it had an overall neutral charge. An energy minimization step was performed to optimize the geometry of the system and remove any steric clashes or unfavorable interactions. This helps to stabilize the system before proceeding to simulations. Optimization of the system was done for both NPT (constant number of particles, temperature and pressure) and NVT (constant number of particles, temperature and volume). These optimizations were carried out at 300 K (Kelvin), a 1 bar pressure, and with restraint forces set to 1000 kJ/mol. Following the optimization steps, molecular dynamics simulations were conducted for a total of 40,000 picoseconds (40 ns). During this simulation, the positions of all atoms in the system were updated over time to study the dynamics and behavior of the H5N1-TLR7 and H5N1-TLR8 complexes. To maintain stability during the simulation, all bonds have been constrained utilizing the LINCS (LINear Constraint Solver) algorithm, ensuring that bond lengths between atoms remain fixed.

Simulations of H5N1 immune response

The C-ImmSim server was used to model and forecast the immunological response to recombinant influenza structures. To anticipate immunological interactions, C-ImmSim combines an agent-based model with machine learning methods and immune epitope prediction. It does this by using a position-specific scoring matrix (PSSM). C-ImmSim replicates three different mammalian anatomical areas at the same time: the bone marrow (representing hematopoietic stem cells, lymphoid cells, and myeloid cells), the thymus (mimicking native T cell selection), as well as a third lymph node resembling typical lymph nodes. The Simulation Steps parameter has been set to 1050 to account for the vaccine's 4-week gap between doses one and two. Time-steps of 1, 336, and 672, each corresponding to eight hours in real life, have been considered for the three injections. Different parameters used in this study were applied as default settings.

Insilico cloning and vaccine optimization

The GeneScript program was used to optimize codons for the production of a recombinant construct in the pET-26b(+) vector. Because the expression host used in the investigation was different from the natural host for influenza, this optimization was required. To adapt the codon usage for the new host, based on the recommended codon use of influenza for production in Escherichia coli K12, the codon optimization has been modified. This ensures that the gene sequence is better suited for efficient expression in the new host. Two important metrics for evaluating codon optimization are the Codon Adaptation Index (CAI) score as well as the GC content. An ideal CAI score is 1.0, indicating a high level of adaptation to the host's codon usage, while a GC content of 30–70% is typically considered optimal for predicting protein expression levels in the host. To facilitate the cloning process, the locations of the NdeI as well as XhoI restriction enzyme sites have been recognized in the optimized nucleotide sequences. The insertion of the recombinant influenza gene sequence into the pET-26b(+) vector and effective target protein production in the selected host are made possible by these restriction sites. SnapGene has been utilized as a tool to aid in this identification and cloning process.

Results

Linear and conformational B cell epitope prediction

B cell epitopes have a crucial part in the H5N1 vaccine due to their ability to stimulate humoral immunity. To identify these epitopes, within the final design, both discontinuous and linear B cell epitopes have been predicted. Linear B cell epitopes have been determined utilizing BepiPred 2.0, while discontinuous epitopes were predicted with the IEDB servers. In total, 14 and 10 linear B cell epitopes have been detected for HA and NA, respectively, and detailed in Table 1. Additionally, Additional file 1: Tables S1 and S2 provide a summary of all B cell epitopes found within the NA and HA proteins.

Table 1 Liner B-cell epitopes predicted from NA and HA antigens

Selection of CTL epitopes prediction

To predict CTL (Cytotoxic T Lymphocyte) epitopes, the NetMHC4 server has been employed specifically for the NA and HA proteins. The selection of MHC class I epitopes was based on their high levels of immunogenicity and antigenicity. Additionally, epitopes with a strong binding score, as determined by the NetMHC4 server, were considered. Select MHC-I epitopes were merged when they overlapped with one another. Ultimately, a total of 21 MHC class I epitopes were chosen, comprising 9 from HA and 12 from NA, which were integrated into the influenza vaccine construct (as detailed in Table 2). There are all predicted CTL epitopes, those that bind to MHC-I, within Additional file 1: Tables S3 and S4.

Table 2 MHC-I epitopes prediction. HA and NA epitopes have been chosen according to MHC-I HLAs binding affinity (nM), %rank and antigenicity rate

Selection of HTL epitopes prediction

For the NA and HA antigens, the method used to predict how peptides will attach to MHC class II alleles has been the NetMHCII-2.3 server. In this process, all MHC class II epitopes have been thoroughly examined, and the epitopes have been chosen according to their strong immunogenicity, and antigenicity, as well as the highest binding score. Additionally, independent MHC class II epitopes were considered when they exhibited overlapping regions. Ultimately, an overall of 9 MHC class II epitopes have been selected for the construction of the H5N1 Vaccine, consisting of 6 HA epitopes and 3 NA epitopes, as outlined in Table 3. Thus, all the predicted HTL (Helper T Lymphocyte) epitopes, which bind to MHC class II, can be found in Additional file 1: Tables S5 and S6.

Table 3 MHC-II epitopes prediction from HA and NA proteins

Population coverage

The IEDB population coverage instrument, accessible at http://tools.iedb.org/population, has been employed to assess the extent of population coverage for each chosen epitope. As illustrated in Table 4, the greatest population coverage was found when both Class I and Class II epitopes were taken into account. This is noteworthy because NA and HA antigens linked to the [Influenza A virus (A/chicken/EastJava/UT551/2010(H5N1))] antigen originate in the Asian region. A comprehensive breakdown of population coverage in terms of Class I, and Class II, as well as the combined Class I and II epitopes associated with various continents has been presented in Table 4 and Fig. 1.

Table 4 Population coverage (%) for specific epitopes containing HLA-binding alleles of MHC class I, class II, as well as combined class I and class II

Fig. 1

Population coverage (%) for specific epitopes containing HLA-binding alleles. For each chosen epitope, a global and average population coverage percentage is taken into account. The colors of MHC class I, and class II, as well as combined class I and class II are displayed in purple, green, and olive-green, correspondingly

Conservancy selected epitopes

The IEDB conservancy analysis instrument has been used to be a monitoring instrument for examining the level of conservation in chosen epitopes relative to other similar sequences. As shown in Table 5 epitopes exhibit varying lengths of minimum 0 and maximum 100 percent conservation.

Table 5 Conservancy selected epitopes from HA and NA epitopes

Final epitopes selection, construction design and physicochemical properties

The nucleotide sequence of H5N1 can be found in Additional file 1: Table S7. Final selected epitopes and adjuvants are presented in Table 6.

Table 6 Final selected epitopes. All overlapping sequences and epitopes are labeled as underlined and in bold form

The final selection of epitopes was interconnected using specific linkers. AAY linkers were chosen to link the MHC-I epitopes, and GPGPGPG linkers have been used for linking the MHC-II epitopes. Adjuvants such as HβD-3 and the PADRE sequence have been connected using EAAAK linkers at the C-terminal site. To bridge the final MHC-I epitope to the initial MHC-II epitope, the HEYGAEALERAG linker was employed. KK linkers were utilized for connecting B-cell epitopes. To complete the structure, a "Histidine Tag" has been included at the C-terminal and attached to the last construct as depicted in Fig. 2A. Also, the secondary structure of the H5N1 construct is shown in Fig. 2B.

Fig. 2

A Amino acid sequence of H5N1 construct. Amino acids in gray are KK linkers, Amino acids in purple are AAY linkers, amino acids in blue are GPGPGPG linkers, amino acids in green are HEYGAEALERAG and amino acids in yellow are EAAAK linkers. B Secondary structure of H5N1 construct. Helixes are shown as orange color, Beta strands are shown as purple color and coils are shown as pink color and disorders are shown as grey color

The outcomes obtained from the ProtParam server revealed that the H5N1 vaccine amino acid sequence, which consists of 350 amino acids, is expected to produce a protein with a molecular weight of ~ 38.46 kDa. The theoretical isoelectric point (pI) and instability index (II) have been determined at 9.32 and 34.44, respectively, indicating that the influenza vaccine falls within the category of stable proteins. Moreover, this vaccine's half-life in a variety of prokaryotic and eukaryotic hosts is predicted to be between 10 and 30 h. The aliphatic index is 84.37 and the grand average of hydropathicity (GRAVY) to be 0.018.

Discontinues B cell from final H5N1 construct

Table 7 provides an overview of the conformational B cell epitopes from the final construct of the H5N1 vaccine. Each of the eight predicted scores for conformational B cell epitopes is shown in Fig. 3.

Table 7 Predicted discontinuous epitope of final H5N1 construct

Fig. 3

3D visualization of discontinuous B cell epitopes in H5N1 construct. Yellow regions are related to discontinuous amino acids in the final H5N1 construct. A Amino acids of F169, G170, P171, G172, P173, G174, P175, G176, P177, T178, T179, Y180, I181, S182, I183, G184, T185, S186, T187, L188, N189, Q190 and R191. B Amino acids of K315, K316, E317, A318, A319, A320, K321, A322, K323, F324, V325, A326, A327, W328, T329, L330, K331, A332, A333, A334, E335, A336, A337, A338, K339, E340, A341, A342, A343, K344, H346 and H347. C Amino acids of R209, T210, L211, D212, F213, H214, G215, P216, G217, P218, G219, P220, G221, N222, L223, Y224, D225, K226, V227, R228, L229, Q230, L231 and R232. DAmino acids of D233, N234, A235, K236, E237, E239, and Y240. E Amino acids of P78, L79, T80, I81, A82 and I85. F Amino acids of M250, R251, and Y254. G Amino acids of G241, A242, E243, A244, E246, and R247. H G161, L162, D163, C164, I165, and R166. (I)Amin acids of V291, L292, S293, C294, L295, P296, K297, S305, T306, R307, G308, and R309

Modeling of 3D construct, refinement and validation

The 3D structure of H5N1 has been modeled using the RoseTTAFold server, as illustrated in Fig. 4. Subsequently, the obtained H5N1 3D structure was subjected to the GalaxyRefine server for further refinement. Model validation was carried out using online software PROCHECK and PROSA. After refining, the findings showed a Z-Score of -0.87, which is used to evaluate the quality of the modeled protein based on NMR or X-ray methodologies (Fig. 5A). For proteins with less than 200 amino acids, the NMR methodology is usually utilized, however for proteins with more than 200 amino acids, the X-ray methodology is applied. The (blue-blue) NMR and (pale-blue) X-ray zones on the PROSA graph denote great simulation accuracy, the lowest error rate, as well as the strongest confidence in the simulated model, respectively, if the Z-Score dot falls inside them. Local model quality associated with the protein's structure is depicted in Fig. 5B.

Fig. 4

3D structure of H5N1. Coils are shown as brown, α-Helixes are shown as blue and β-Sheets are shown as purple

Fig. 5

A Z-Score validation of H5N1 3D homology modeling. A value of − 0.87 shows that the 3D H5N1 structure is located in X-ray crystallographic proteins. B The local model quality of H5N1 confirms that the energy of this structure is close to a stable manner. C Ramachandran plot associated with vaccine 3D structure. 98% of residues are modeled in favored regions and more than 99.8% of residues are modeled in allowed regions

Based on the Ramachandran graph analysis of the anticipated model, 336 out of 350 residues, or 96% of all residues, were found inside the preferred areas (with a 98% probability). In total, approximately 99.7% of all residues, or 349 out of 350, were within the allowed regions (with a likelihood exceeding 99.8%). The Ramachandran graph organizes amino acids based on the angles of phi and psi, as shown in Fig. 5C. This data reflects the overall quality and reliability of the model.

Protein–protein docking

The structures of H5N1 including TLR7 and TLR8 proteins have been submitted to the HDOCK server to determine the contact areas by protein–protein docking. The highest-ranking results for every complex have been chosen based on the least intermolecular binding energy among entire H5N1-Proteins complexes with the lowest docking score as shown in Tables 8 and 9.

Table 8 An overview of the top ten H5N1-TLR7 models

Table 9 An overview of the top ten H5N1-TLR8

The docking results indicate that H5N1 is capable of forming interactions with the 3D structures of TLR7 and TLR8. Within the final construct, two sequences, PADRE and HβD-3, were selected and included as adjuvants at the C-terminal. As displayed in Figs. 6A and 7A, the initial amino acid regions of H5N1 were observed to be engaged with TLR7 and TLR8. Moreover, H5N1 directly established interactions with several amino acids within the binding sites of TLR7 and TLR8. This demonstrates the potential for a strong interaction between H5N1 and these immune-related receptors.

Fig. 6

A Docking HADDOCK outcomes for the H5N1-TLR7 complex are visualized, B and the residues' interaction with the complex has been magnified. The H5N1 structure is shown as a green stick-filled carton. The blue cartons and sticks that make up the TLR7 structure are displayed. The dashed line in yellow represents hydrogen bonds

Fig. 7

A Docking HADDOCK outcomes for the H5N1-TLR8 complex are visualized, B and the residues' contact with the complex has been magnified. An illustration of the H5N1 structure is a carton filled with green sticks. The cooper carton and sticks that make up the TLR8 structure are displayed. The dashed line in yellow represents hydrogen bonds

As illustrated in Figs. 6B and 7B, it was observed that specific amino acids in H5N1, including SER846, SER850, SER860 and PHE864 interacted with amino acids in the TLR7, such as ARG232, TYR224, PHE214 and LEU211. Also, GLN58, SER59 and ASN 75 amino acids from the H5N1 structure interacted with GLU438 and ARG628 amino acids of the TLR8 structure. This interaction demonstrates the potential binding between H5N1 and the TLR7 and TLR8.

The H5N1 exhibited interaction with TLR7 through 4 hydrogen bonds. Similarly, in the H5N1-TLR8 complex, interactions were facilitated by 3 hydrogen bonds to the interactions in this complex. These findings demonstrate the significance of hydrogen bond contacts in forming these complexes.

The 2D interactions are summarized and visualized, providing a detailed depiction of these interactions (H5N1-TLR7 in Fig. 8A and H5N1-TLR8 in Fig. 8B).

Fig. 8

A Presentation of amino acids’ 2D interaction in H5N1-TLR7 complex. Green residues are presented as H5N1 structure and blue residues are presented as TLR7 structure. B Presentation of 2D interaction of amino acids in H5N1-TLR8 complex. Green residues are presented as H5N1 structure and blue residues are presented as TLR8 structure

Molecular dynamic and MM/PBSA analysis

The complexes underwent molecular dynamics (MD) simulations for a duration of 100 ns. The MD simulation results indicated the hydrogen bonds’ creation in H5N1-TLR7 and H5N1-TLR8 complexes with root mean square deviations (RMSD) of 2 and 3.5 nm, respectively, as shown in Figs. 9A and 10A. These RMSD values provide insights into the stability and behavior of the simulated complexes during the MD simulations.

Fig. 9

MD simulation outcomes of H5N1 complex. A RMSD of H5N1-TLR7 in Pico seconds, B RMSD of H5N1-TLR7 in residues, C Gyration of H5N1-TLR7 in Pico seconds, D H-bond of H5N1-TLR7 in Pico seconds and E SASA graphs of H5N1-TLR7 complex

Fig. 10

MD simulation outcomes of H5N1-TLR8 complex. A RMSD of H5N1-TLR8 in Pico seconds, B RMSD of H5N1-TLR8 in residues, C Gyration of H5N1-TLR8, Pico seconds (D) H-bond of H5N1-TLR8 in residues and E SASA graphs of H5N1-TLR8 complex

The root mean square fluctuation (RMSF) diagrams revealed that residues in the H5N1-TLR7 and H5N1-TLR8 complexes during MD simulation exhibited very weak volatility, as depicted in Figs. 9B and 10B.

These results suggest that the H5N1 structure is maintained in the binding interactions with TLR7 and TLR8 complexes. Additionally, it's worth noting that the H5N1-TLR7 complex returned to equilibrium after approximately 2 ns, while the and H5N1-TLR8 complex reached equilibrium around 3 ns after the start of the MD simulations. These findings indicate the stability and equilibration of the simulated complexes throughout the simulations.

The gyration radius data associated with H5N1-TLR7 and H5N1-TLR8 complexes can be found in Figs. 9C and 10C, providing insights into the compactness and overall structure of these complexes. Additionally, the H-bond and SASA overall results of H5N1-TLR7 and H5N1-TLR8 complexes are included in Figs. 9D and 10D as well as Figs. 9E and 10E, respectively, offering further details about the interactions within these complexes.

The H5N1-TLR7 and H5N1-TLR8 complexes' binding affinities were verified using the use of Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) modeling. Tables 10 and 11 provide an overview of the findings of these calculations, which identify the energy connected to the bonds that are created in each complex between the vaccine and the receptor. In the H5N1-TLR7 complex, the total binding energy was calculated to be − 29.97 kJ/mol, indicating a strong affinity between H5N1 and TLR7. On the other hand, in the H5N1-TLR8 complex, the total binding energy was found to be − 23.9 kJ/mol. These calculations provide insights into the thermodynamic stability and binding affinities of these complexes, with a notably high affinity observed in the case of H5N1-TLR7.

Table 10 MMPB/SA calculation from H5N1-TLR7 complex

Table 11 MMPB/SA calculation from H5N1-TLR8 complex

Additional investigations have shown that apart from total binding energy, van der Waals and electrostatic energies also have a substantial part in influencing complexes. In the H5N1-TLR7 complex, according to Table 10 and Fig. 11, the electrostatic energy has been determined to be − 100.18 kJ/mol, and the van der Waals energy to be − 43.49 kJ/mol. Conversely, in the H5N1-TLR8 complex, the van der Waals energy contribution was − 38.06 kJ/mol and the electrostatic energy was − 132.02 kJ/mol (Table 11 and Fig. 12).

Fig. 11

Graphical presentation of energy contribution of H5N1-TLR7 complex. Vander Waals (VDWAALS) energies are shown as pink, Electro static (EEL) energies are shown as brown, effective polarizable bonds (EPB) are as green, effective non-polarizable bonds (ENPOLAR) are shown as pale-green, gas phases (GGAS) are shown as blue, solvent phases (GSOLV) are shown as purple and total energies contribution (TOTAL) are shown as hot-pink

Fig. 12

Graphical presentation of energy contribution of H5N1-TLR8 complex. Vander Waals (VDWAALS) energies are shown as pink, Electro static (EEL) energies are shown as brown, effective polarizable bonds (EPB) are as green, effective non-polarizable bonds (ENPOLAR) are shown as pale-green, gas phases (GGAS) are shown as blue, solvent phases (GSOLV) are shown as purple and total energies contribution (TOTAL) are shown as hot-pink

These energy components further characterize the nature of the contacts in complexes, with both van der Waals and electrostatic forces playing a role in influencing the binding affinities and stability of these complexes. The balance of van der Waals as well as electrostatic energies contributes to the overall stability of these interactions.

Immune simulation

The C-ImmSim web server's results showed a rise in secondary immune response generation that is consistent with the true immunological response. Notably, As the immunological response grew more intense, there was a rise in total IgM + IgG1 and IgG2 antibody levels, as depicted in Fig. 13A. Moreover, the quantity of B memory (y2) increased as the levels of IgG1 and IgG2 isotypes decreased, as depicted in Fig. 13B. Additionally, the outcomes demonstrated that an increase in TH memory (y2), as depicted in Fig. 13C, correlated with an increase in both TH and TC cell populations, as illustrated in Fig. 13D. Furthermore, there was an elevation in the focus of produced IFN-γ and an increase in TH cell population, as seen in Fig. 13E, F. Together, these results show that the C-ImmSim server's simulation of the immunological response to influenza H5N1 vaccination closely resembles the real immune response, underscoring the simulation's accuracy in simulating the activity of the immune system.

Fig. 13

Immune response simulation against the H5N1 strain. A Immunoglobulin levels upon H5N1 construct. B Total B isotype (IgM, IgG1, IgG2) as the immunological response, and amount of B memory (y2). C CD4 T-helper cells as a result of the development of T-helpers. D Displaying both memory and total T-helper cells. E Cell populations of CD8 T-cytotoxic lymphocytes in each condition after being injected with antigen. F Generated IL-4, IL-6, and IL-12 concentrations together with IFN-γ

Insilico cloning and vaccine optimization

The nucleotide sequences of the finished construct were optimized using GenScript. The organism chosen for host expression was E. Coli K12. To prevent bacterial ribosome binding sites, cleavage sites for restriction enzymes, and rho-independent transcription terminators, certain settings were made in GenScript. GC-Content for the improved Influenza construct, which consists of 1032 nucleotides (excluding the His-Tag), was determined to be 54.46%. These values affirm the likelihood of successful protein expression. To facilitate intracellular expression, to exclude the PelB sequence from the final structure, Nde I and Xho I restriction enzyme sites have been introduced to the N- and C-terminals of the nucleotide sequence, respectively. Additionally, a stop codon has been included following the His-Tag sequence. Ultimately, SnapGene was used to clone the design into the pET-26b(+) plasmid, as shown in Fig. 14.

Fig. 14

Cloned nucleotide sequence of H5N1 into pET-26b(+) plasmid. The H5N1 sequence is labeled as a green ribbon. NdeI and XhoI enzymes are located at N-terminal and C-terminal of H5N1 construct. His-Tag (6xHis) sequence is located at the C-terminal of H5N1 construct

Discussion

Influenza H5 viruses are responsible for respiratory infections in humans, and the severity of these infections can vary from no symptoms to being life-threatening. In the realm of infectious diseases, the influenza virus stands out as among the most substantial hazards to both human and avian populations worldwide. Virologists have issued a warning over the possibility of a brand-new, extremely damaging influenza pandemic, highlighting the pressing necessity of being ready and vigilant [27]. In the two decades following its initial emergence in China in 1996, approximately 60 nations have documented ongoing cases of H5 viruses affecting domestic poultry, wild bird populations, and human beings. This underscores the widespread and continuous nature of H5 virus incidents around the world [28, 29].

Recent progress in the fields of immunoinformatics, bioinformatics, and structural vaccinomics has led to a significant transformation in how antigens are identified. These advancements have also facilitated the development of the novel vaccines creation approach, enabling the vaccines targeting a broad spectrum of bacterial and viral infections. Our prior research, has demonstrated that multi-epitope peptide vaccines can use more diverse and higher effective factors against infectious agents compared to vaccines based on single epitopes [16, 30]. In line with this approach, the present study concentrated on multiepitope subunit vaccines. These vaccines are performed using multiple immunogenic factors of pathogens and possess the capability to direct the humoral immunity reaction towards particular antigenic epitopes, ultimately leading to a safer and more efficient immune response.

In this current work, we employed the NCBI database to obtain the sequences of two viral HA/NA proteins from AIV (H5N1) strain, following a comprehensive review of the existing literature. Two viral proteins, namely HA (GenBank: BAL61222.1) and NA (GenBank: BAL61230.1), serve as surface glycoproteins responsible for binding to host cells; are two of the proteins among the most abundant within the influenza A virus and play crucial functional and structural parts throughout life cycle of viruses [31].

Vaccines stimulate B cells to generate antibodies that, in turn, carry out effector functions by specifically interacting with pathogens or toxins. Due to their role in preserving memory cells, providing longer-term immunity, and defending versus reinfection, the humoral immunity reaction plays a vital part in safeguarding versus infections through vaccination. Numerous antigens and vaccinations also cause T-cell responses in addition to B-cell responses [32,33,34].

Well-known CD4⁺ T cell is another crucial type of immune cell capable of adopting Th1 or Th2 phenotypes as well as orchestrating immune reactions [35]. The Th1 response is responsible for activating cytotoxic CD8⁺ T lymphocytes natural killer cells, as well as macrophages. In contrast, the Th2 response has a crucial part in activating B cells, facilitating isotype switching, promoting affinity maturation, as well as generating antibodies that can target and neutralize external pathogens [36]. T cell epitope-based vaccination is therefore a unique strategy for generating a strong immune response against pathogenic pathogens [37]. Here, to found possible CTL and HTL immunogenic epitopes inside HA and NA protein, we used topological screening, the MHC-I and MHC-II binding predictions, and the VaxiJen server. These epitopes have been selected according to their ability to link to a wide range of HLA-A and HLA-B alleles having high binding affinity.

In this study, we utilized tools to evaluate the allergenicity of chosen epitopes. During the immune stimulation procedure, allergenicity, which poses a significant challenge in vaccine development, is detected in many vaccine candidates, potentially leading to an "allergic" response. According to a criterion for predicting allergenicity, a sequence is deemed possibly allergenic if, when compared to known allergens, it shares the identity of at least six consecutive amino acids within an 80-amino acid window [38]. In addition to that, one important factor in the reverse vaccinology technique is population coverage. As per the results, all the predicted T cell epitopes from HA and NA are capable of providing coverage for populations across the majority of geographic regions worldwide, with rates exceeding 90%. Moreover, since MHC superfamilies are important in the creation of vaccines and pharmaceuticals, it has become possible to ascertain the functional link among MHC variants by MHC cluster analysis.

In our study, we used a total of 18 B cell epitopes, CTL and HTL. We used two algorithms of Linear B-cell epitope prediction methodology, including IEDB server as well as BepiPred 2.0 for predicting amino acid likely B-cell epitopes. While over the past few decades, numerous B-cell epitope prediction methods have been developed, but they have mostly failed, the majority of experimental epitope determination efforts have concentrated on finding a linear B-cell epitope prediction technique [39]. Furthermore, accurate prediction and a quicker, less expensive vaccine design process can both be aided by linear B-cell epitope prediction [40].

Based on the Bepipred linear epitope prediction 2.0 for the AIV-A (H5N1) strain, the most effective B cell epitopes for each of the three proteins were selected as possible vaccination candidates. Linear epitopes from B cells, epitopes recognized by cytotoxic T cells (CTL) and helper T cells (HTL), along with Beta-defensin 2, as well as the PADRE peptide sequence, have been employed in the formulation of the ultimate vaccine proteins. In the past, adjuvants have been used as immunomodulators to increase the effectiveness of certain vaccination candidates [41,42,43]. The effectiveness of the Beta-defensin adjuvant as an immune booster has been demonstrated in numerous experiments against various organisms [44,45,46]. The VaxiJen server has been utilized to evaluate all of the acquired protein sequences, determining the most potential antigenic protein with the capacity to induce immunity.

AAY linkers have been utilized to bond CTL (Cytotoxic T Lymphocyte) epitopes to each other. These linkers help separate and align CTL epitopes effectively. Specific HEYGAEALERAG linker serves as the interface among CTL epitopes as well as HTL (Helper T Lymphocyte) epitopes, ensuring a proper connection and interaction between these two critical components of the vaccine. GPGPG linkers act as connectors between different HTL epitopes. They help maintain structural integrity and spacing within the HTL epitope region. KK linkers have been used to connect B-cell epitopes, as well as have been utilized to link HβD-3 and PADRE adjuvants in the N-terminal of the construct. These linkers ensure that B-cell epitopes are properly connected, enhancing their accessibility for immune recognition. EAAAK was used to bond the vaccine construct to the HisTag sequence at the C-terminal. EAAAK Linkers help maintain flexibility and proper spacing between these components. The strategic use of these linkers is essential in vaccine design to ensure that different epitopes and components work together efficiently, promoting a robust immune response when administered.

Additionally, The ProtParam server has been utilized to examine the final construct's physicochemical characteristics. The molecular weight of the H5N1 vaccine construct is 38.46 kDa. This is deemed a suitable vaccine candidate, as proteins having a molecular weight of less than 110 kDa are generally regarded as more favorable for vaccine development because of their ease of purification [47]. After expression, the construct's heat stability was evaluated using the aliphatic index. The H5N1 protein is classified as stable by the estimated instability index since its value is less than 40. An instability index of less than 40 indicates stability, whereas a value above 40 suggests potential instability. The calculated GRAVY index of H5N1 serves as an indicator of its polar nature and its interaction with water. A greater solubility for this structure is indicated by a lower GRAVY score. Furthermore, the selected vaccine candidate exhibited a higher proportion of coil structure in its secondary structure. This suggests that random coils are significant in conferring greater protein flexibility and might contribute to an enhanced antibody binding capability [48].

Incorporating adjuvants into the subunit vaccine could enhance the immunity reaction it elicits. Adjuvants serve to be activators for Toll-like receptors (TLRs) and could be effectors to broaden antibody recognition [49]. It was observed that TLR7 and TLR8 are capable of recognizing Influenza Vaccine Antigen (IVA) due to their ability to couple with single-stranded RNAs (ssRNAs) [50]. TLR-7 was demonstrated to trigger both a humoral and long-lasting memory response in reaction to ssRNA-mediated infections as well as vaccinations, such as those involving Influenza and HIV [51, 52]. Additionally, It has been discovered that synthetic TLR7 and TLR8 agonists enhance both innate and adaptive immune responses [53].

The interaction between H5N1 3D structure and TLR7 and TLR8 structures has been examined through protein–protein docking. In the case of H5N1-TLR7 as well as H5N1-TLR8 complexes, it was observed that these interactions involved 4 and 6 hydrogen bonds, respectively.

After docking and determining how to connect complexes H5N1-TLR7 and H5N1-TLR8 we performed the MD test. The MD technique was performed to check the stability and behavior of the complexes per unit of time. Since it is possible that over time, the proteins in a complex move away from each other and the complex becomes unstable, MD can show the binding quality of two proteins resulting from their docking (within a certain time).

So in our study to do MD:

• (1) First, we ran each complex separately for 100 ps and analyzed the results through RMSD, RMSF, Gyration, H-bond and SASA graphs.

• (2) Then we went one step further and analyzed the MD results using MMPBSA to estimate the contribution of each of the energies involved in the complex such as electrostatics, van der Waals, etc.

The complexes of H5N1 with TLR7 and TLR8 demonstrated interactions between the proteins via the creation of hydrogen bonds. Analysis of RMSD diagrams for the H5N1-TLR7 and H5N1-TLR8 complexes revealed that the simulations were relatively stable. However, it's noteworthy that Compared to the H5N-TLR8 complex, which had a mean RMSD of 0.3 nm, the H5N1-TLR7 complex had a mean RMSD of 0.6 nm, indicating that its intermolecular energies are higher. The stability of the complex is confirmed by the gyration graph, which shows that the amino acids stay inside and don't go beyond their radial axis range during the simulation. Furthermore, the RMSF plot aligns with these outcomes, indicating that the amino acids within binding sites of H5N1-TLR7 and H5N1-TLR8 complexes keep their stability during simulation.

Overall, the results from the docking and molecular dynamics (MD) simulations suggest that these complexes are capable of maintaining their stability during the simulation. The outcomes from such a study indicate that recombinant construct may effectively stimulate TLR7 and TLR8. The fact that TLR7 and TLR8 attach to the residues in the vaccine structure suggests that the vaccination suggested in this study would probably cause the development of antibodies against it. This expectation stems from the importance of having the accurate 3D structure of antigen for it to be identified effectively via the humoral immunity mechanism. Consequently, all selected epitopes within the 3D-vaccine structure should closely mimic the correct fold found in the native structure of the antigen to elicit a comparable immunity response [54].

Additionally, the MM/PBSA outcomes revealed that, apart from the overall energy, the complexes are influenced by the energies arising from electrostatic and covalent bonds. Notably, the positive binding energy observed in the H5N1-TLR8 complex could be influenced by presence of salt bridges in contact between the TLR8 and vaccine structure [55].

Influenza viruses are ssRNA viruses that require the activation of T cell-dependent immune responses plus humoral immunity induced via B cells. This activation is crucial for preventing the proliferation and intracellular survival of the virus. According to predictions from the C-ImmSim Online server, the growth of both B cells and T cells leads to the establishment of long-term immune memory. Thus, IgG2 and IgG1 which are indicative of Th2 and Th1 responses to influenza antigens, respectively, play important roles in protecting against influenza infection [56, 57]. Past research has indicated that both avian and human hosts experience a reduction in their T-cell populations in reaction to avian influenza infection [58]. A decrease in T lymphocytes is frequently linked to proinflammatory cytokines’ upregulation, particularly interferons (IFNs), notably IFN-γ, and interleukins (ILs) like IL-6. These changes can, in the end, result in a condition called hypercytokinemia and trigger the cell apoptosis’ activation [59]. Increased levels of IL-12 and IFN-γ are indicative of a robust cellular immune response. While previous research has demonstrated that T epitopes from H5N1 strains designed to mimic human sequences may lead to weakened IFN responses in humans, the results from the ICM server indicate that this vaccine has an excellent capacity to stimulate IFN-γ, suggesting its potential to trigger a robust immune response [60,61,62,63,64]]. Nonetheless, the findings from the ICM server suggest that a rise in the count of Th1 cells is linked to the elevation in the number of cytotoxic T cells.