H5N1 influenza multi-epitope vaccine design through immunoinformatics approaches

1 Introduction

H5N1 is an avian influenza virus primarily affecting birds. While rare in humans, infections are severe and often fatal. Preventing its spread is crucial to protect poultry industries and, more importantly, to minimize the risk of the virus mutating to gain the ability to sustain human-to-human transmission, which could trigger a devastating global pandemic [1, 2]. Furthermore, its increasing detection in other mammals highlights its evolving threat. Robust surveillance and control measures are vital to prevent further spread of H5N1 and safeguard both animal and human health [3].

Today, immunoinformatics provides powerful computational approaches and a range of accessible web services and resources for predicting crucial properties related to the immunogenic potential of a protein, starting from its primary amino acid sequence and/or the physicochemical properties. These methods leverage our understanding of the complex interactions within the immune system to identify potential immunogenic regions of the viral constituent molecules that can constitute viable components for vaccines. They enable the prediction of key immunogenicity aspects, such as B-cell epitopes and binding to the Major Histocompatibility Complex (MHC). Furthermore, these approaches facilitate the evaluation of crucial factors related to the safety and biological feasibility of vaccine components, including potential toxicity and allergenicity. Additionally, the structural stability of the designed vaccine construct can also be assessed through suitable computational analyses.

B-cell epitopes need to be identified as targets of immune responses. B cells recognize conformational or linear epitopes on the surface of intact antigens. Immunoinformatics tools like ABCpred [4], BepiPred [5], and CBTOPE [6] predict B cell epitopes based on various sequence-derived features such as hydrophilicity, surface accessibility, flexibility, and the presence of specific motifs. Some tools also attempt to predict discontinuous epitopes if the 3D structure of the protein is available or can be predicted.

MHC binding of T-cell peptides is another crucial step in both humoral and cellular responses. A fundamental step is to identify which peptide fragments derived from the proteins of the antigen are likely to bind to MHC (or Human Leukocyte Antigen, HLA in humans) molecules. MHC class I molecules present peptides to cytotoxic T cells (CD8+ T cells), while MHC class II molecules present peptides to helper T cells (CD4+ T cells). NetMHC (and its variants like NetMHCpan, NetMHCII, NetMHCIIpan [7]) developed by the the Technical University of Denmark, this is one of the most widely used and highly regarded suites of tools. NetMHCpan predicts peptide binding to any MHC Class I molecule (including non-human) using artificial neural networks trained on large datasets. NetMHCIIpan predicts peptide binding to any MHC Class II molecule using similar pan-specific approaches. Other variants exist for specific species or types of MHC molecules. The IEDB Analysis Resource is a comprehensive online platform that hosts many prediction tools, besides NetMHC tools, along with tools for predicting antigen processing (proteasomal cleavage, TAP transport) and other aspects of epitope presentation. SYFPEITHI is an older but still useful database and prediction tool based on position-specific scoring matrices derived from known MHC binding peptides [8]. RANKPEP predicts MHC binding using position-specific scoring matrices and is also available through the IEDB [9]. MHCflurry is a newer, open-source tool for MHC Class I binding prediction based on deep learning, known for its speed and accuracy.

While MHC binding is a critical step, the peptide must also be generated through antigen processing and transported to where it can bind MHC [10]. Some tools attempt to incorporate these steps are NetChop which predicts proteasomal cleavage sites, that is a key step in generating peptides for MHC Class I presentation [11]. TAPpred which predicts the likelihood of a peptide being transported into the endoplasmic reticulum by the Transporter associated with Antigen Processing (TAP), another step in MHC Class I presentation [12]. These tools, utilize machine learning algorithms trained on experimental data to predict peptide-MHC binding affinity based on the amino acid sequence of the peptide and the specific MHC allele.

As for toxicity and allergenicity predictions to ensure the safety of potential vaccine candidates, immunoinformatics tools can predict if a peptide sequence has the potential to be toxic or allergenic based on comparisons with known toxic and allergenic peptides. AllergenOnline is widely used database and sequence similarity-based tool for comparing novel sequences against known allergens [13]. AllergyPred is a web server for allergen prediction [14]. AlgPred is another web server that uses various approaches, including amino acid composition, motif searching, and machine learning, to predict allergenicity [15]. Other examples of web server for predicting protein allergenicity potential from protein sequences are AllerCatPro [16] and AllerTOP [17]. Concerning the prediction of toxicity, ToxinPred is a tool specifically designed for predicting toxic peptides using machine learning techniques [18]. ToxPi (Toxicity Prioritization Index), while not a peptide-specific predictor, is a framework and tool used for integrating various sources of toxicity data for prioritization and assessment, which can be relevant in a broader safety context [19]. VISH-pred is an ensemble of fine-tuned ESM models for protein toxicity prediction which achieved to date the best accuracy and area under receiver operating curve scores on independent test sets [20].

Overall, immunoinformatics offers a powerful suite of computational tools for analyzing the primary sequence of target antigens to predict various properties crucial for immunogenicity. These In silico methods significantly accelerate vaccine design and the identification of promising immunotherapeutic targets by prioritizing candidates for experimental validation. The extensive use of these tools in recent scientific literature demonstrates their value in identifying potential vaccine candidates against infectious diseases and active immunotherapy targets for cancer [21–24].

In this work, we outline a step-by-step process for creating a potential multi-epitope vaccine for H5N1. The design focuses on eliciting an immune response against the most conserved molecules of the pathogen.

While several immunoinformatics studies have proposed multi-epitope constructs against H5N1 or related influenza strains, the present work differs in the sequence dataset which is drawn exclusively from recent human infections, capturing the most current circulating variants including the 2024 Clade 2.3.4.4b spillover events. Moreover, unlike most comparable studies that focus on surface glycoproteins alone, epitope mining here spans both surface and internal structural proteins, broadening the potential for cross-protective immunity. HLA alleles were chosen to cover the regions most affected by H5N1 outbreaks and human cases. The selected epitopes were assembled into a construct specifically designed for recombinant production and subjected to a rigorous multi-step safety filtering pipeline. Finally, immune efficacy was evaluated using an agent-based simulator that captures a richer stochastic dynamics than aggregated differential equations models.

Section 2 reports the starting data set and the tools used in the analysis and Section 3 describes the results obtained including the vaccine construct and its validation through structural modeling and computer simulations.

2 Materials and methods

The H5N1 virus, like other influenza A viruses, is composed of several proteins. We have downloaded data from the curated collection of virus sequences at NCBI Virus [25] by querying for VirusLineage = Alphainfluenzavirus influenzae, taxid:2955291, Serotype = H5N1, and ProtNames = hemagglutinin, neuraminidase, matrix protein 1, and matrix protein 2.

Specifically, we have downloaded all curated structural proteins of H5N1 virus, namely Hemagglutinin (HA), Neuraminidase (NA), Matrix protein 1 (M1), and Matrix protein 2 (M2), of H5N1 that have been found in human and avian infections in the period from Jan 1, 2023 to April 17, 2025. This particular initial period of 2024 is significant in the context of H5N1 due to a couple of key developments: (i) in late March 2024, H5N1 avian influenza was detected in dairy cows in the United States (Texas). This was notable because H5N1 was primarily known to affect birds, and its emergence in cattle marked a change in the behavior of the virus; (ii) following the detection in dairy cows, the first human cases of H5N1 infection linked to exposure to infected dairy cows were reported in early April 2024. This raised concerns about potential transmission pathways and the risk to individuals working with these animals.

As of April 2025 we obtained 1118, 987, 989, and 1,090 sequences, for HA, NA, M1, and M2 respectively. Sequences were retrieved from a curated repository in which submissions undergo automated and manual quality review; the dataset includes isolates from both human and avian infections, broadening alignment depth and strengthening the identification of conserved regions across host-adapted variants. Any incomplete or low-quality sequences not filtered at source contribute gap-rich columns and are therefore implicitly excluded during conserved-region extraction.

The selection of the four key structural proteins (Hemagglutinin, Neuraminidase, Matrix protein 1, and Matrix protein 2) is based on their vital roles in the viral life cycle and their potential as vaccine targets. Hemagglutinin (HA), represented by the “H5” in H5N1, is a surface glycoprotein that initiates infection by binding the virus to sialic acid receptors on host cells. For H5N1, the HA protein often prefers alpha-2,3 sialic acid receptors, which are more common in the lower human respiratory tract. This preference may contribute to the severe pneumonia seen in human infections. Critically, HA is the main antigen that triggers neutralizing antibodies, making it the primary target for vaccine development. Neuraminidase (NA), the “N1” in H5N1, is another crucial surface glycoprotein. It facilitates the release of new virus particles from infected cells by cleaving sialic acid and also aids the virus in moving through mucus. Like HA, NA is an important target for both antiviral drugs and antibodies. While HA and NA are the primary focus for current H5N1 vaccine development, researchers are also exploring other proteins, specifically the more conserved internal proteins M1 and M2, to enhance protection and broaden the immune response. Matrix protein 1 (M1) forms a layer directly beneath the viral envelope, providing essential structural support and assisting in virus assembly. It is also involved in transporting the viral ribonucleoprotein complex into the host cell's nucleus. Matrix protein 2 (M2) is an ion channel protein embedded within the viral envelope. This protein plays a critical role in the virus's uncoating process once it enters the host cell. Targeting these internal proteins, M1 and M2, offers the potential for broader protection because they are more conserved across different influenza strains and subtypes. This concept is central to the development of a universal influenza vaccine. Some vaccine candidates are already designed to elicit T-cell responses against M1 and M2. This approach can provide cross-protective immunity by enabling the immune system to recognize and target infected cells, regardless of the specific HA or NA subtype present on the viral surface.

For vaccine design, it is important to identify epitopes that are conserved across different strains of a pathogen and that can elicit an immune response in a broad range of individuals with diverse MHC alleles. To this purpose, for each of the four molecules HA, NA, M1 and M2, we have performed a multiple sequence alignment (MSA) using Clustal Omega [26, 27]. From the alignment we have selected the conserved regions (with a minimal length of 9 amino acids) as subsequences with potential to be targeted. In the alignment we have considered exact amino acid match, strong amino acids substitutions (i.e., substitutions with a strongly similar amino acid), as well as weak conservation (i.e., substitutions with a weakly similar amino acid) for an extended coverage of the molecular diversity observed in infectious variants of the virus. To maximize the robustness of the conservation analysis, the alignment included sequences from both human and avian isolates. Conservation was assessed globally in a single per-protein alignment without within-clade stratification. Clustal Omega annotates each column with ^* (all residues identical), : (all residues within a strongly similar substitution group), or . (all residues within a weakly similar group), imposing an effective 100% per-column conservation requirement. Gap-containing columns receive a blank annotation and are implicitly excluded. Conserved blocks were extracted with a custom Python script retaining contiguous annotated runs of length ≥9 amino acids.

The alignment resulted in conserved subsequences for each of the protein under consideration. We got 7, 7, 1, and 2 respectively for HA, NA, M1, and M2. These are shown in Table 1. Worth to note that, since the number of contiguous subsequences of length l within an amino acid sequence of length L is calculated as L−l+1, the number of potential class I peptides (L = 9, 10, 11, 12) for the proteins studied are: 1,794 (HA), 1,418 (NA), 970 (M1), and 232 (M2), summing up to 4,414. As for potential class II peptides (with L = 9, 10, 11, 12, 13, 14, and 15), the corresponding counts are: 3,066 (HA), 2,408 (NA), 1,687 (M1), and 385 (M2), for a total of 6,546.

Table 1. Conserved amino acid sequences of the proteins Hemagglutinin, Neuraminidase, and Matrix protein 1 and 2 of the H5N1 virus after multiple-sequence alignment of curated sequences of human and avian strains from Jan 1, 2023 to April 17, 2025.

The selection of vaccine candidates from the generated pool of 10960 peptides is dependent on the HLA haplotype chosen. Specifically, we identified peptides predicted to exhibit strong binding to the relevant class I and II molecules of that haplotype, as here described.

To maximize the immunogenicity and global applicability of our multi-epitope vaccine, we selected HLAs that are common in a large fraction of the world's population with a focus on regions where one or more outbreaks in birds have been reported, as well as regions where human cases of highly pathogenic H5N1 infection have been observed [28]. These strategies ensured that the vaccine candidate could be effective in diverse populations with varying distributions of HLA alleles.

In order to capture the diversity of multiple regions in the world, we focused on high-frequencies class I and class II alleles [29]. Specifically, our class I selection comprised HLA-A02:01 [this allele is very common in all EUR populations and about 75% of the populations of North America and South America. It is also frequent in South-West, Central and South Asia, and in Oceania (Pacific, New-Guinea, and Australia)], HLA-A24:02 (frequent in the majority of populations in all regions but Sub-Saharan Africa), and HLA-C04:01 [also frequent in most regions of Sub-Saharan Africa, South-West, Central and South Asia, North and South America, and Oceania (Pacific, New-Guinea and Australia)]. As for class II we considered DQA10501 [common in large proportions of Sub-Saharan Africa, Europe, North-East Asia, South-East Asia and Oceania (Pacific, New-Guinea and Australia) populations], DQA10101 (common in most populations of all regions), and DPA10103 (this allele is very frequent in almost all populations of all regions). For NetMHCIIpan predictions, complete alpha–beta heterodimer pairs were specified as DQA1^*05:01/DQB1^*02:01, DQA1^*01:01/DQB1^*05:01, and DPA1^*01:03/DPB1^*04:01, corresponding to the most globally frequent pairing for each alpha chain. The same pairings are used internally by C-ImmSim, which does not allow independent specification of the beta chain for DQ and DP molecules.

The vaccine design process followed a sequence of steps aimed at filtering the large initial pool of potential epitopes and peptides to select the most promising components, as illustrated in Figure 1.

The computational tools employed for this process are briefly described here, with the detailed results presented in Section 3.

For predicting B-cell epitopes, we used the BepiPred-3.0 tool, a sequence-based epitope prediction tool that exploits protein language model embeddings and residue annotation for predicting linear and conformational epitopes [5].

To predict class I peptides or cytotoxic T lymphocyte (CTL) peptides, we utilized NetMHCIpan-4.1, a cutting-edge tool designed to predict peptide binding to MHC class I molecules [7]. NetMHCIpan-4.1 employs machine learning algorithms trained on experimental binding data and can predict peptide-MHC interactions for any known or unknown MHC class I allele. The tool provided binding affinity scores, which were used to classify peptides into strong binders (high affinity) and weak binders (low affinity). Peptides with high affinity were preferred for their potential to induce significant immunological responses.

For predicting class II petides or Helper T lymphocyte (HTL) peptides, we used the NetMHCIIpan-3.2 tool, a well-established method for predicting peptide binding to MHC class II molecules [30]. This tool leverages a pan-specific model that predicts binding affinities for any MHC class II allele, even those without experimental data. The tool provided binding affinity scores indicating how well peptides bind to the selected MHC-II molecules. Again, peptides were classified based on their binding strength, with high-affinity peptides being selected for their potential to elicit strong immunological responses.

We used ToxinPred3, an enhanced web server and standalone tool for predicting peptide toxicity [18]. The tool was used to screen all predicted epitopes and peptides for toxicity. The input consisted of the peptide sequences, and the output provided insights into the potential cytotoxicity of each peptide, which were crucial for ensuring the safety of the vaccine candidates.

AllerTOP is an online server used for In silico allergen prediction [17]. It utilizes a machine learning approach based on amino acid E-descriptors (i.e., key physicochemical properties), auto-cross covariance transformation, and dataset filtering. The tool was used to calculate the allergenicity scores for each of the predicted B-epitopes, CTL peptides, and HTL peptides. Peptides flagged as allergenic were filtered out to ensure the safety of the vaccine candidates.

The final set of epitopes and peptides was assembled into a multi-epitope vaccine construct. It incorporates β-defensin-3 as an adjuvant, leveraging its immunomodulatory properties, activating and recruiting immune cells, thereby enhancing the response [31, 32]. The EAAAK linker then connects the adjuvant to the B epitopes to provide a rigid, alpha-helical structure, crucial for maintaining optimal spatial presentation and reducing steric hindrance for effective immune recognition [33]. Subsequently, the KK linker, a short and positively charged sequence, is used to facilitate specific molecular interactions and provide minimal, defined spacing between B epitopes. The AAY linker is used to link CTL peptides, potentially influencing local structure and stability [34]. Lastly, the linker GPGPG between HTL peptides is utilized for its inherent flexibility, ensuring individual peptides can fold correctly and remain accessible for optimal immune receptor engagement [35].

The tertiary structure of our multi-epitope vaccine construct was computationally determined using a hierarchical modeling approach. This method systematically builds the structure, moving from simpler, local predictions to more complex, global arrangements. It starts from the secondary structure prediction, which identified the fundamental local folding patterns of the polypeptide chain, distinguishing between alpha-helices, beta-sheets, and disordered coils (used as crucial initial constraints. Next, threading (or fold recognition) attempts to fit the vaccine's sequence onto a library of known 3D protein structures, leveraging existing knowledge of protein folds, and using the secondary structure predictions to guide optimal alignments. For regions without suitable templates or for novel segments, ab initio modeling was employed, predicting structures from scratch solely based on the amino acid sequence and physical principles, often incorporating the initial secondary structure predictions as guides. Finally, the resulting preliminary 3D model underwent rigorous atomic-level refinement techniques. This critical last step involved computationally adjusting the atomic positions, optimizing bond lengths and angles, and minimizing the overall energy of the structure to correct any minor inaccuracies, resolve steric clashes, and ensure the predicted tertiary structure was physically realistic and stable.

We have used PSIPRED [36] served as the first tool to predict secondary structural elements which revealed the presence of α-helices β-sheets and coil regions. The construct's structural motifs received further analysis through PDBsum [37] to compare them with existing structural templates. The three-dimensional structure prediction involved running i-TASSER [38–40] alongside Phyre2 [41]. The structural modeling process of i-TASSER depends on threading and ab initio methods to create reliable models through multiple threading alignments and iterative assembly simulations yet Phyre2 uses profile-profile matching and hidden Markov models to effectively model distant homologs. The final structure underwent energy minimization and structural refinement through ModRefiner [42] and GalaxyRefine [43] to optimize side-chain conformations and backbone geometry into physically realistic and energetically favorable models. The stereochemical quality and stability assessment of the final structure involved Ramachandran plot analysis [44] for backbone dihedral angle evaluation, ProSA for the recognition of errors in the three-dimensional structure [45, 46], and the Structure Validation Server SAVESv6.1 which includes ERRAT for non-bonded atomic interaction analysis [47], MolProbity for local stereochemical correctness and atom clashes (Ramachandran plots) [48], Verify3D for the assessment of the models with three-dimensional profiles [49, 50], ProCheck to check the stereochemical quality of a protein structure by analyzing its overall and residue-by-residue geometry [51, 52], and WhatCheck which is a suite of programs to assess the stereochemical quality of a given protein structure [53, 54].

3 Results

Starting from the list of conserved molecular sequences reported in Table 1, we have applied the prediction tools mentioned above to predict B-epitopes and T-peptides.

3.1 Prediction of epitopes and peptides

Linear B-cell epitopes, have been predicted with BepiPred-3.0 [5]. This resulted in 21 epitopes of varying length (Table 2).

Table 2. B epitope sequences of conserved targeted proteins of H5N1.

Chosen epitopes are marked with (+).

For class I peptides we obtain 40 strong binders (i.e., 22 for HA, 8 for NA, 8 for M1, and 2 for M2). These are shown in Table 3.

Table 3. CTL peptide sequences of conserved targeted proteins of H5N1.

Chosen peptides are marked with (+).

S_*EL* (Score_EL) and R_*EL* (%Rank_EL) are output values that help predict the likelihood of a peptide being presented via the Eluted Ligand (EL) prediction method. R_*EL* (Percentile Rank - Eluted Ligand) is a normalized score that expresses the S_*EL* as a percentage compared to the scores of a large set of random natural peptides (often around 10⁶) of the same length. It indicates the percentile of the predicted binding affinity (inferred from eluted ligand data) of the query peptide relative to this background set. Lower R_*EL* rank means the peptide is predicted to be a better ligand (more likely to be eluted) compared to most random peptides. Common thresholds for identifying potential strong and weak binders based on R_*EL* are often used: R_*EL* < 1% for Strong binders, R_*EL*∈(1, 2] for Weak binders (R_*EL*>2 means no-binders).

For class II peptides we obtain 58 strong binders (i.e., 15 for HA, 23 for NA, 13 for M1, and 2 for M2). These are shown in Tables 4, 5.

Table 4. HTL peptide sequences of conserved targeted proteins of H5N1 (continue on Table 5).

Chosen peptides are marked with (+).

Table 5. HTL peptide sequences of conserved targeted proteins of H5N1.

M1 and M2 ONLY (continue from Table 4). Chosen peptides are denoted with a (+).

Besides the B-epitopes that are not allergenic and not toxic and are marked with a (+) in Table 2, only peptides determined not to be allergenic or toxic were considered for selection. From this subset, to maximize coverage, a maximum of two peptides per allele exhibiting the highest S_*EL* and lowest R_*EL* were chosen, as are indicated by (+) in Tables 3–5.

Population coverage analysis was performed by submitting the full set of epitope–HLA pairs to the IEDB Population Coverage tool. Across the combined Southeast Asia, Africa, Europe, and North America populations, the 90th percentile population coverage (PC90) reached 3.48 for CTL epitopes and 10.56 for HTL epitopes. When extended to the global population, PC90 values were 3.93 and 10.45, respectively (global population coverage reached 74% for CTL epitopes and 84% for HTL epitopes when evaluated against their respective HLA alleles). Both epitope sets exceeded the threshold for excellent coverage (PC90 ≥ 2.0) in all evaluated populations, with HTL epitopes showing particularly high redundancy across HLA class II alleles. The near-identical PC90 values between the combined regional and global analyses indicate that the selected epitope–HLA pairs capture global HLA diversity without substantial regional bias.

To assess the risk of autoimmune cross-reactivity, all B cell, CTL, and HTL epitopes were screened against the human proteome via BLASTp. The reference database comprised 197,870 human proteins retrieved from NCBI RefSeq (filtered by Organism: Homo sapiens, taxid:9606). B cell epitope searches were run without length restrictions, whereas CTL and HTL peptide searches were constrained to matches of length 9, consistent with our restrictions. The PAM30 and PAM70 substitution matrices were selected for their sensitivity to short, highly similar alignments, with compositional adjustment disabled. Alignment results for T cell epitopes were subsequently filtered to retain only length-9 matches; no exact hits to human proteins were identified, indicating a low theoretical risk of autoimmune cross-reactivity for the candidate epitopes.

3.2 Assembling the multi-epitope vaccine

The assembled polypeptide chain vaccine construct comprises the adjuvant (i.e., β-defensin-3) linked to the epitopes via the EAAAK linker, followed by the 13 B-epitopes joined by the linker KK, followed by the (14+15) 29 T-peptides linked with AAY and GPGPG linkers (Figure 2). The construct is preceded at the N-terminus by the sequences HHHHHH and LVPRGS, serving specific functional purposes in recombinant DNA technology and protein engineering for the development and production of subunit vaccines [55, 56]. In particular, HHHHHH is a Histidine tag (His-tag) whose primary purpose is for affinity purification of the recombinant vaccine protein. His-tags bind strongly and specifically to immobilized metal ions (like nickel or cobalt) on chromatography columns, allowing the tagged protein to be easily separated from other host cell proteins. The tag is placed at the N-terminus of the protein to ensure it is accessible for binding to the purification resin. LVPRGS is a protease cleavage site, specifically a Thrombin cleavage site. Thrombin is an enzyme that recognizes and cleaves polypeptide chains at specific amino acid sequences. The sequence Leu-Val-Pro-Arg•Gly-Ser is a recognition and cleavage site for Thrombin, where the cleavage occurs after the Arginine residue.

The 3D structure of the folded multi-epitope vaccine is presented in Figure 3. The leftmost panel details the secondary and tertiary structure of the vaccine construct. Progressing across the figure from left to right, the remaining panels sequentially locate the positions of the adjuvant, the B epitopes, and the T peptides relative to the complete architecture.

3.2.1 Physicochemical characterization of the vaccine construct

The physicochemical properties of the designed multi-epitope vaccine were analyzed using the ProtParam tool (ExPASy [57]). The vaccine construct is composed of 698 amino acids, with a predicted molecular weight of 76.05 kDa and a theoretical isoelectric point of 9.39 (Table 6). An instability index of 36.91 indicates that the protein is stable, while the aliphatic index (68.22) suggests good thermostability. The negative GRAVY value (–0.593) reflects the hydrophilic nature of the construct. In addition, the estimated half-life of the vaccine protein was approximately 30 h in mammalian cells, over 20 h in yeast, and more than 10 h in Escherichia coli, supporting its suitability for recombinant expression and downstream experimental studies.

Table 6. Physicochemical properties of the multi-epitope vaccine construct predicted using ProtParam (ExPASy).

3.3 Modeling of the 3D construct, refinement and validation

The removal of the protease cleavage site (LVPRGS) and His-tag (HHHHHH) from the polypeptide chain led to a final construct of 699 amino acids. This multi-epitope vaccine construct has then been inspected for structural features to verify stability before proceeding with downstream simulation studies.

The H5N1 vaccine construct was subjected to 3D structure modeling using the i-TASSER and Phyre2 servers, integrating motif information from PDBsum and secondary structure data from PSIPRED. ModRefiner and GalaxyRefine were used to enhance bond angles, together with side-chain conformations and steric interactions within the initial 3D model.

The model validation process utilized PROCHECK and ProSA as well as Ramachandran analysis, together with the SAVESv6.1 suite. The ProSA Z-score result of −3.67 falls within the range typical of computationally predicted models and is consistent with structural plausibility for models of this size and complexity, although the graphical output remains inaccessible. The X-ray methodology proves more reliable for proteins that exceed 200 amino acids.

The Ramachandran plot of the vaccine model appears in panel A of Figure 4. The stereochemical acceptability of the structure is confirmed by the distribution of most residues in the favored and additionally allowed regions. The geometry of the structure remains stable despite the presence of a few expected outliers in the large recombinant construct.

The local model quality evaluation in panel B of Figure 4 shows both averaged and raw scores for each amino acid residue. Adequate local structural quality is suggested by the majority of residues falling within the favorable range.

The constructed model received a score of 69.19% on the ERRAT quality factor assessment, as shown in panel C of Figure 4. Despite falling short of ideal values from high-resolution crystal structures, the score is consistent with expectations for In silico models and indicates that the model is a suitable starting point for further refinement. The measured ProSA Z-score and ERRAT values are within the expected range for In silico-constructed multi-epitope vaccines, though further optimisation is ongoing.

3.4 Structural validation of the vaccine construct

The predicted three dimensional structure of the multi-epitope vaccine construct H5N1-MEVC was assessed using multiple well-established structural validation tools. The ProSA analysis produced a Z-score of –3.67, suggesting that the overall quality of the predicted model is consistent with experimentally solved protein structures of similar size. ERRAT analysis yielded an overall quality factor of 81.07%, indicating favorable non-bonded atomic interactions across most regions of the model. Verify3D evaluation further supported the structural reliability, with 80.95% of residues showing acceptable compatibility with their respective local 3D environments, meeting the recommended threshold. In addition, MolProbity analysis resulted in an overall score of 2.13 with a clash score of 10.40, reflecting reasonable stereochemical quality. Ramachandran plot analysis revealed that 88.22% of residues were located in favored regions, while only 1.44% were identified as outliers, Overall, these validation metrics indicate that the predicted structure is of moderate quality. While the model provides a reasonable basis for downstream in-silico analyses, further experimental validation will be required to confirm its structural accuracy and biological relevance

3.5 Binding interaction analysis of the H5N1 multi-epitope vaccine construct with Toll-like receptors

Toll-like receptors play a central role in sensing viral antigens and initiating innate immune signaling, which subsequently supports adaptive immune responses. Therefore, the binding behavior of the designed vaccine construct H5N1-MEVC was evaluated against TLR2 and TLR4 through molecular docking. The docked complexes showed that the H5N1-MEVC was positioned at the receptor interface of both TLR2 and TLR4, indicating favorable structural accommodation of the vaccine construct within the receptor-binding region. The residue-level interaction maps further supported the docking results. In the TLR2-H5N1-MEVC complex, the vaccine construct formed several close contacts, hydrogen bonds, hydrophobic interactions, and salt-bridge interactions with receptor residues, mainly involving chain B of TLR2 and chain V of the vaccine construct. Similarly, the TLR4-H5N1-MEVC complex showed multiple receptor-vaccine contacts involving TLR4 receptor chains and the vaccine construct. These interactions included hydrogen-bonding, electrostatic, hydrophobic, and close-contact interactions, suggesting that both complexes maintained stable molecular interfaces (Figure 5).

Overall, the docking and residue-interaction profiles indicate that the designed H5N1-MEVC can interact effectively with both TLR2 and TLR4. These findings support the potential immunostimulatory relevance of the vaccine construct and provide a structural basis for further stability assessment using molecular dynamics simulation.

3.6 Molecular dynamics simulation of H5N1-MEVC complexes with TLR2 and TLR4

To examine whether the docked vaccine-receptor complexes remained stable under dynamic conditions, 50-ns molecular dynamics simulations were carried out for the H5N1 multi-epitope vaccine construct (H5N1-MEVC) in complex with TLR2 and TLR4. The trajectories were analyzed using RMSD, radius of gyration, RMSF, and hydrogen-bond profiles. The RMSD plots showed that both complexes underwent an initial adjustment phase and then remained within a stable range during the simulation (Figure 6). The TLR2-H5N1-MEVC complex showed an average RMSD of 1.30 ± 0.39 nm, with a final RMSD of 1.48 nm. The TLR4-H5N1-MEVC complex showed an average RMSD of 1.41 ± 0.32 nm, with a final RMSD of 1.11 nm. These results suggest that both docked complexes retained their overall structural arrangement during the 50-ns simulation and did not show continuous conformational drift.

The radius of gyration analysis was used to assess the compactness of the complexes. The TLR2-H5N1-MEVC complex showed an average Rg of 5.89 ± 0.07 nm, while the TLR4-H5N1-MEVC complex showed an average Rg of 6.88 ± 0.04 nm. The small variation in Rg values indicates that both complexes maintained stable compactness throughout the simulation. The higher Rg value of the TLR4 complex is consistent with the larger size of the TLR4 receptor complex and does not indicate instability.

RMSF analysis further supported the stability of the simulated complexes. Most residues showed moderate fluctuations, while higher flexibility was limited to specific regions, likely corresponding to terminal or loop regions. The average RMSF values were 0.87 ± 0.52 nm for the TLR2 complex and 0.77 ± 0.41 nm for the TLR4 complex. This indicates that the overall receptor-vaccine frameworks remained stable, with flexibility restricted to localized regions.

Hydrogen-bond analysis showed that both complexes maintained persistent interaction networks during the simulation. The TLR2-H5N1-MEVC complex showed an average of 1,188 ± 26 hydrogen bonds, whereas the TLR4-H5N1-MEVC complex showed an average of 1,467 ± 20 hydrogen bonds. The final hydrogen-bond counts were 1,220 for TLR2-H5N1-MEVC and 1,468 for TLR4-H5N1-MEVC. The higher number of hydrogen bonds in the TLR4 complex suggests a more extensive stabilizing interaction network; however, this should be interpreted as interaction persistence rather than direct binding affinity, since MM-PBSA/MM-GBSA free-energy calculations were not performed.

Overall, the MD simulation results indicate that the H5N1-MEVC remained structurally stable in complex with both TLR2 and TLR4 (Figure 7). The stable RMSD behavior, consistent compactness, localized residue fluctuations, and persistent hydrogen-bond profiles support the structural reliability of the docked receptor-vaccine complexes.

3.7 IFN-g induction potential

To evaluate the interferon-gamma (IFN-g) induction potential, both B-epitopes and T-peptides were screened using the IFNepitope server (http://crdd.osdd.net/raghava/ifnepitope/) [58]. This platform predicts the likelihood of peptide-mediated IFN-g secretion by leveraging a combination of motif-based mapping and Support Vector Machine (SVM) models. Based on these predictions, the potential B-cell epitopes (QNALNGNGDPNNMDRAVKL and ESMREEYRQEQQSAVD), the CTL peptide (GILGFVFTL), and the HTL binders (RQEQQSAVD and QQSAVDVDD) were all classified as positive IFN-g inducers.

3.8 In silico cloning and expression optimization

To evaluate the expression feasibility of the vaccine construct, codon optimization was performed using the Java Codon Adaptation Tool (JCat, [59]) with Escherichia coli K12 selected as the host expression system. During optimization, options were selected to avoid rho-independent transcription terminators, prokaryotic ribosome binding sites, and undesirable restriction enzyme cleavage sites (Figure 8).

The optimized nucleotide sequence consisted of approximately 2,102 bp, corresponding to the 698-amino-acid vaccine construct. The optimized sequence exhibited a Codon Adaptation Index (CAI) value of 1.0 and a GC content of 36.9%, indicating favorable compatibility with the E. coli expression machinery and suggesting a high probability of efficient recombinant expression.

To facilitate cloning, NcoI and XhoI restriction sites were introduced at the 5′ and 3′ termini of the optimized coding sequence, respectively. The codon-optimized insert was subsequently cloned In silico into the pET-28a(+) expression vector using SnapGene software (SnapGene^® software by Dotmatics available at: https://www.snapgene.com). Restriction cloning generated a recombinant plasmid of approximately 7,333 bp containing the vaccine coding sequence under the control of the T7 promoter and compatible with expression in E. coli BL21(DE3). The resulting construct provides a theoretical framework for downstream recombinant protein expression and purification using a His-tag-based system.

3.9 Simulations of the response to the vaccine

Finally, to check the vaccine for efficacy, we have used C-ImmSim¹ to simulate the immune response. This is an agent-based computational simulator rooted in the Celada-Seiden model [60], designed to mimic the intricate dynamics of the human immune response to antigens. The tool allows researchers to perform In silico immune simulations, modeling the complex interactions between various immune cells, such as B lymphocytes, T helper cells (including their subsets), cytotoxic T lymphocytes, and antigen-presenting cells, as well as the molecular recognition events involving HLA presentation [61].

C-ImmSim can integrate immunoinformatics predictions to simulate the expected humoral and cellular immune responses elicited by candidate vaccines, including multi-epitope constructs designed through computational approaches [62]. It is commonly used in research to computationally evaluate the potential efficacy of these designed vaccines in eliciting a robust immune response, providing valuable predictive data before costly and time-consuming in vitro and in vivo experiments are conducted

We have simulated three vaccination protocols, corresponding to (i) single-dose, (ii) prime-boost, and (iii) prime-boost-challenge to a virtual individual with the same HLA haplotype used above, namely A0201, A2402, C0401 for class I, and DPA1_0103, DQA1_0101, DQA1_0501 for class II. In all three protocols the first injection (priming) has been performed at time 0, corresponding to a virtual individual that has never been exposed to the epitopes and peptides contained in the vaccine (Figure 2). Subsequently, the boost injection is given at day 30 and the challenge injection at day 60.

The full simulation settings were as follows. Each run covered 547 time steps of 480 min each (time step dt= 8h), corresponding to approximately 182 days of simulated time. The simulated compartment representing muscle tissue at the injection site had a volume of 10 μL on a 3D spatial grid. Each injection delivered 50,000 antigen entities with an antigen replication rate of 0.75. Injections were administered at simulation steps 0, 90, and 180.

When the vaccine is injected in the simulated volume representing the muscle tissue of an individual, it elicits an immune response through a series of coordinated events involving both the innate and adaptive immune systems. When administered, these particles, which carry the vaccine antigens (such as the B-epitopes and the T-peptides), are recognized by the host as foreign. The presence of an adjuvant significantly influences and enhances the phases of immune recognition. Components of the vaccine, including the multi-epitope peptides and the adjuvant, are immediately encountered by resident APCs of the innate immune system at the injection site, such as macrophages (MA) and dendritic cells (DCs). Adjuvants triggers the activation of these innate cells, leading to the release of pro-inflammatory cytokines (such as IFN-g, TNF-a and a generic danger signal, see Figure 9D). This creates a local inflammatory environment that is essential for attracting more immune cells to the site and initiating the adaptive response.

Activated APCs (including B-cells that recognize the epitopes in Table 2 specifically), and particularly dendritic cells, efficiently engulf and process the multi-epitope peptides and any associated vaccine components into smaller fragments. These fragments are then loaded onto the HLA molecules (Figure 9E). Recognition and presentation are modeled through a stochastic selection process where the probability is determined by the normalized S_*EL* scores presented in Tables 3–5. Peptides designed to stimulate T helper cells are presented on class II molecules, while epitopes designed for T cytotoxic cells are cross-presented on HLA class I molecules.

Upon activation and antigen uptake, the now activated APCs undergo maturation and migrate from the injection site to the nearest draining lymph nodes. In C-ImmSim, these two compartments are represented by different volumes (set of voxels) where cells encounter and interact. Migration from one volume representing the tissue (e.g., the muscle, in this case) and the lymph node, is modeled through stochastic diffusion with a drift provided by the chemotaxis.

In the lymph node compartment, T helper cells recognize peptides presented on class II HLA molecules. Upon successful recognition and co-stimulatory signaling from the APC, naive CD4+ T cells become activated and duplicate (Figure 9C).

Naive B cells can encounter intact multi-epitope peptides or vaccine particles directly in the lymph node. They receive a secondary signal (T-dependent response) from activate T helper cells and start proliferating. This leads to the differentiation of B cells into antibody-producing plasma cells and memory B cells (Figure 9B). Plasma cells secrete large quantities of antibodies. Initially IgM, followed by IgG (Figure 9A). The antibody levels rise over time, reaching a peak, generally 2–6 weeks after injection. Inset plot of Figure 9A shows IgG titers exceeding IgM approximately 9 days after first inoculum.

While an multi-epitope vaccine does not cause intracellular pathogens infection, the cross-presentation of peptides can prime CTLs for future encounters (not shown).

A critical outcome of the primary response is the generation of long-lived memory T and B cells. This is shown in panel B and C of Figure 9. Upon subsequent exposure to the actual pathogen displaying these epitopes, these memory cells are rapidly activated, leading to a much faster and stronger secondary immune response.

These simulations show how the multi-epitope vaccine leads to the activation of mainly the humoral arms of the adaptive immune system and the establishment of immunological memory.

The kinetics and intensity of the immune response, including the antibody peak and cellular immunity, are highly dependent on the vaccine components including the adjuvant. Although a single dose primes the immune system, a prime-boost regimen is typically required for optimal maturation and durable protection. We extended the investigation simulating also a three dose protocol (prime-boost-challenge). Our simulations compared these three protocols, revealing increased antibody peaks about 4 weeks after each injection (Figure 10). The third curve corresponds to the prime-boost-challenge protocols which achieves an even higher antibody counts at day 82.5 (i.e., about 3 weeks after the 3rd dose, hence earlier compared to the other two). The Figure 10 also highlights the rapid and significant increase in antibody production, boosted by approximately 140% between single and triple dose.

4 Conclusions

This study demonstrates the significant potential of immunoinformatics as a powerful suite of computational tools to streamline vaccine design and the identification of promising immunotherapeutic targets. By analyzing the primary sequence of target antigens and predicting crucial immunogenic properties, In silico approaches can substantially accelerate the process of prioritizing candidates for experimental validation.

In this work, we applied a systematic immunoinformatics pipeline to design a multi-epitope vaccine candidate against the H5N1 avian influenza virus. Our approach integrated the identification of conserved regions across key viral proteins (HA, NA, M1, and M2), the prediction of B-cell, CTL, and HTL epitopes, and rigorous filtering for potential toxicity and allergenicity. This comprehensive In silico analysis successfully narrowed down a large initial pool of potential vaccine components (~11, 000) to a targeted set of promising candidates (42). We then assembled these selected components into a multi-epitope polypeptide chain using suitable linkers and performed structural analysis and computer simulations to evaluate its potential to induce a lasting immune response.

Our findings highlight that this computational strategy, leveraging existing immunoinformatics tools, not only accelerates the initial stages of vaccine development but also provides a robust foundation for subsequent experimental investigation and validation. This capability is particularly valuable for enabling a more rapid response to emerging infectious threats like novel influenza strains.

It is important to acknowledge that immunoinformatics relies on experimental data for its development and validation, and In silico predictions serve as powerful guides rather than replacements for traditional experimental research methods. Crucially, certain aspects are not yet fully accounted for in the development of epitope/peptide prediction tools, nor their integration into the proposed pipeline. Machine-learning-based B epitope prediction, for instance, only indirectly incorporates wet-lab findings by learning from positive samples in training datasets. This can lead to overlooking nuances such as the well-known fact that B cells and antibodies targeting conserved epitopes are often non-neutralizing. Nevertheless, these antibodies can induce significant indirect protective efficacy through antibody-dependent cellular cytotoxicity (ADCC) mediated by NK cells and antibody-dependent cellular phagocytosis (ADCP) driven by macrophages. Therefore, these Fc-mediated effector functions remain crucial for conferring protection against infection [64].

Several limitations of the present work should be acknowledged. The adjuvant choice (β-defensin-3) and epitope ordering follow conventions established in the multi-epitope vaccine literature and were not subject to comparative or combinatorial optimisation; the optimal adjuvant, formulation, and delivery system remain open questions requiring dedicated experimental investigation. The predicted 3D structure is currently undergoing further refinement; molecular dynamics simulation and receptor docking, together with experimental structure determination, would ultimately be required before drawing conclusions about structural stability or innate immune activation potential. Finally, the C-ImmSim outputs are theoretical predictions of immunogenic potential and should not be interpreted as evidence of actual vaccine efficacy; in vitro antigen presentation assays, animal immunization studies, and viral challenge experiments remain necessary for any translational claim.

To conclude, the core strength of an immunoinformatics approach lies in its ability to significantly reduce the number of vaccine candidates requiring extensive laboratory testing. This computational efficiency directly alleviates many funding challenges inherent in traditional vaccine design, thereby allowing for the reallocation of critical resources toward the experimental validation of promising In silico-predicted targets and constructs.

Generative AI statement

The author(s) declared that Generative AI was used in the creation of this manuscript. Exclusively for English language proofreading and grammatical correction.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Gennady Bocharov, Marchuk Institute of Numerical Mathematics (RAS), Russia

Reviewed by: Otun Saha, Washington State University, United States

Jitendra Malviya, Barkatullah University Bhopal, India

Footnotes

http://c-immsim.iac.cnr.it