The 3D structure of the SARS-CoV-2 spike ectodomain in an open state, represented by PDB ID 6VYB (Walls et al., 2020), was chosen for the monomeric spike protein. This structure was obtained through electron microscopy, with a resolution of 3.20 Å, and was released on 2020–03-11. The selection of 6VYB was based on its representation of the spike protein in a conformation relevant to its interaction with host cells, which is crucial for eliciting an immune response against the virus antigen. The crystal structure of human L-ferritin, represented by PDB ID 2FG8 (Wang et al., 2006), was selected for the ferritin protein. This structure was determined using X-ray diffraction, with a resolution of 2.50 Å, and was released on 2006–07-04. PDB ID 2FG8 was chosen for its well-defined atomic details, which is crucial for accurate structural analysis of ferritin in the designed vaccine. The lectin protein structure was obtained from PDB ID 5B1W (Nagae et al., 2016). This structure was determined using X-ray diffraction, with a resolution of 3.05 Å, and was released on 2016–05-11.

The selected structure of the monomeric spike was used to dock to one unit of ferritin with the HADDOCK server (van Zundert et al., 2016), where the HADDOCK server performed rigid body energy minimization in the initial step of the docking process. Residues at the C-terminal of monomeric spike (residues: 1,136–1,147) were chosen as a binding interface area according to previous studies (Powell et al., 2021; Wuertz et al., 2021; Gordon et al., 2022) and several docking scores (Supplementary Figure S1). One exposed loop on the ferritin structure (residues:73–91) was selected as the most plausible binding site (Masoomi Nomandan et al., 2022), which was confirmed by comparing the energetic terms and clustering of the resulting docked poses. The docked pose with the minimum HADDOCK score and lowest RMSD values was considered for model development. The poses were validated via further molecular docking and processed by the ClusPro web server (Supplementary Table S1) (Kozakov et al., 2017). Details outlining the atomistic interactions of the modeled construct are provided (Supplementary Table S2; Supplementary Figure S2). Biomolecular visualization of the complex was carried out by PyMOL (LLC, 2015).

Lectins are a unique class of proteins or glycoproteins known for their ability to selectively recognize carbohydrate structures and form reversible linkages upon interacting with glycoconjugate complexes (Raposo et al., 2021). The molecular docking analysis examined the interaction between the monomeric spike-ferritin complex and lectin molecules, aiming to explore the potential interference of lectin with vaccine efficacy. HADDOCK server (van Zundert et al., 2016) was employed to dock full-length human lectin and monomeric spike-ferritin nanoparticle vaccine together (Supplementary Figure S3). It was observed that seven residues involved in the ferritin-lectin binding interface were also found to participate in the ferritin-monomeric spike interface. However, the most stable cluster exhibited a Root Mean Squared Deviation (RMSD) value of 4.7±0.0 Å and a HADDOCK score of 163.2±41.3. This suggests that the connection between lectin and the constructed vaccine is weak and unstable, and might not significantly impact the vaccine efficacy. These observations were further validated through ferritin and lectin docking (Supplementary Figure S4). The most stable cluster displayed an RMSD value of 20.1±0.2 Å and a HADDOCK score of 188.3±23.4. Based on these results, it is unlikely that ferritin or the monomeric spike would be bound to lectin proteins in a manner that would significantly interfere with their intended roles in the vaccine’s function.

The dynamic behavior of the monomeric spike SARS-CoV-2-ferritin complex was investigated through molecular dynamics simulations using GROMACS 2021.4 software (Hess et al., 2008). OPLS-AA/M force field (Jorgensen et al., 1996) was employed to generate parameters and topology files. The system was solvated in a cubic box of water molecules with a volume of 3,140.21 nm³, ensuring a 10 Å boundary around the proteins. TIP3P (Jorgensen et al., 1983) was used to simulate water molecules and counter ions were added to neutralize the system (13 Na + ions). A two-step approach involving 50,000 steps of steepest descent and 50,000 steps of conjugate gradient energy minimization was utilized to relieve the system from steric clashes. In the first step, the steepest descent algorithm was used with a maximum force tolerance of 500 kJ/mol/nm and an energy step size of 0.01 nm. Subsequently, the conjugate gradient algorithm was applied with a maximum force tolerance of 100 kJ/mol/nm and an energy step size of 0.01 nm, until convergence was achieved. The LINCS algorithm (Hess et al., 1997) was employed to constrain the bond lengths. Nonbonded interactions (electrostatic and VDW) were calculated with a cutoff value of 10 Å. The Particle Mesh Ewald method (Darden et al., 1993) with 0.16 nm Fourier grid spacing (Saini et al., 2019) was used for long-range electrostatic interactions. In the initial equilibration phase, the NVT ensemble was utilized to maintain a constant temperature of 310 K throughout a 100 ps simulation with a time step of 1 fs. The system was heated up from 0 K to 310 K over four steps using an annealing method. Following this, the temperature was maintained at 310 K using the V-rescale temperature coupling algorithm, a modified version of the Berendsen thermostat. To ensure structural stability, position restraints were applied to the proteins to prevent large fluctuations. Following the NVT equilibration, the system underwent further equilibration through a 250 ps simulation under an NPT ensemble at 310 K and 1 bar pressure, with a time step of 1 fs. Temperature and pressure coupling were achieved employing V-rescale and Berendsen methods, respectively. Finally, under NPT, with the V-rescale thermostat maintaining 310 K and the Parrinello-Rahman pressure coupling algorithm maintaining 1 bar, three replicates of 300 ns with a time step of 2 fs were conducted.

Bio3D package (Grant et al., 2006) in R was used for calculating the Root Mean Squared Deviation (RMSDs), Root Mean Squared Fluctuations (RMSFs), Principal Component Analysis (PCA), and Cross-Correlation of Cα atoms (DCCMs) for each replicate simulation during 300 ns of simulation time. Distances between centers of masses of spike and ferritin were calculated using VMD (Visual Molecular Dynamics) software (Humphrey et al., 1996). Secondary structure timeline analyses were calculated over 300 ns of simulation time with VMD STRIDE plugin. VMD was also utilized for calculating Solvent Accessible Surface Area (SASA). The visualization of structures was performed using PyMOL (LLC, 2015) and Figure 2 was created using BioRender.com.

The immunological simulation of the monomeric spike-ferritin vaccine was performed using the C-IMMSIM server (Rapin et al., 2010; Castiglione et al., 2021) to describe the immune response profile and immunogenicity of the antigenic peptide. The experiment was carried out by adding FASTA sequences of vaccine constructions into the simulation, which lasted approximately 12 months (a time step is about 8 h). The random seed and simulation volume was set to 12,345 and 50. Three in silico injections with no LPS were administered at time steps 1, 84, and 168, with a minimum time interval of 30 days between two injections. Additionally, to ensure the reliability of the simulation results, two additional simulations were carried out, with the random seed set to different numbers: one with the seed 3,452 and the other with 15,342 (Supplementary Figures S5, 6).

The peptides from monomeric spike proteins of the designed vaccine that are naturally processed by Major Histocompatibility complex (MHC) class I molecules were predicted using the Immune Epitope Database (IEDB) (Vita et al., 2019). This analysis employed the IEDB MHC prediction tool, MHC-NP (Giguère et al., 2013), which utilizes machine learning algorithms to predict peptides with potential binding affinity to MHC class I molecules. The ectodomain sequence of the monomeric spike protein of SARS-CoV-2 was utilized as the input sequence for the analysis. To ensure relevance to human immune responses, the chosen MHC source species was human. Additionally, to ensure comprehensive coverage of potential T-cell epitopes, five common Human Leukocyte Antigen (HLA)-B variants were selected as the MHC alleles. Peptide lengths ranging from 8 to 11 amino acids (all lengths) were assigned to capture a diverse range of potential epitopes. Finally, 3,620 peptides were identified for each of the examined variants, each demonstrating different affinities for MHC class I binding. Furthermore, MHC-II binding predictions were conducted to identify potential T-cell epitopes, utilizing the NetMHCIIpan EL prediction method (Reynisson et al., 2020). This analysis spanned across the 19 most common HLA-DR alleles, ensuring a broad representation of human genetic diversity. The default peptide length of 15 amino acids was employed for these predictions. This approach aimed to identify epitopes with high binding affinity to MHC-II molecules, thereby enabling investigation into vaccine efficacy across genetic variations in the host immune system.

The ferritin nanoparticle vaccine was modeled with the monomeric spike antigen by molecular docking. The docking score obtained for the selected binding pose and the RMSD value were −88 and 1.1 (Å), respectively (Figure 3). The ferritin full-cage structure was assembled with 24 ferritin single units and has shown twenty monomeric spikes on ferritin structure, with four ferritin units remaining without adjuvants due to steric hindrance (Figure 4).

By examining the mutations that occurred in the spike protein of variants of concern, including Alpha, Beta, Gamma, Delta, and Omicron (Supplementary Figure S7), it was noticed that not only do none of the mutations exist in the ferritin binding interface but the center of occurrence of mutations is also far from the binding interface of ferritin in the monomeric spike vaccine model. While the mutations D1118H and H1101D, associated with the Alpha and Beta variants respectively, are in close proximity to the binding interface, it is crucial to note that they are not directly situated within it. The binding interface is specifically defined by residues within the range of 1,136–1,147, where neither D1118H nor H1101D fall within this range.

Thus, the monomeric spike-ferritin vaccine model can still be effective despite the known mutations.

Molecular dynamics simulations were utilized to thoroughly examine the structural stability and flexibility of the monomeric spike-ferritin complex and its constituent parts. The RMSD analysis was employed to assess the global dynamics of the monomeric spike-ferritin structure. The general trend of time-dependent RMSD plots indicated an initial increase, followed by stabilization after 200 ns (Figure 5A). The average calculated RMSD values across three replicate simulations were 16.5 Å, 27.5 Å, and 17.4 Å, respectively. Notably, the observed high RMSD values are primarily attributed to the substantial size of the monomeric spike-ferritin nanoparticle structure (1,200 a. a) (Figure 2). These observations aligned with the expectation of reduced system stability due to the size of the spike monomer, which naturally results in a higher degree of conformational changes. The distribution patterns of calculated RMSD values derived from the RMSD histogram aligned with the findings of the time-dependent RMSD plots (Figure 5B). All replicates of the overall system revealed stable RMSD values with peaks at 16 Å. That is expected due to multiple domains within its structure (Figures 1, 2).

The flexibility of the complex was further investigated through RMSF analysis, focusing on the Cα atoms per residue (Figure 6A). The average calculated RMSF values were 6.3 Å, 8.9 Å, and 8.4 Å for the three replicate simulations. These findings corroborated the insights gained from RMSD analysis and indicated that the monomeric spike-ferritin complex is highly stable. The RMSF plots revealed multiple peaks, which indicated the presence of distinct flexible loops and their dynamic motions within the structure. Specifically, the RBD region exhibited prominent peaks at residues (476–490) and (498–507) (Figure 6B), highlighting the dynamic motions within these specific loops. Similarly, fluctuations were observed in the ferritin structure, particularly at residues (155–158), suggesting the dynamic nature of the loops within its structure (Figure 6B). However, it is important to note that while these loops were flexible, they did not directly interfere in the spike-ferritin interface and thus the stability of the designed vaccine may not be affected. Loops one and 2, which were located in the interface of RBD with ACE2, were far from the monomeric spike interface with ferritin, and, therefore, would not affect the stability of the modeled vaccine. Loop three in ferritin was also located far from the interface. Thus, despite the flexibility indicated by the peaks in the RMSF plots, these loops did not influence the stability of the designed vaccine.

The Principal Component Analysis (PCA) and its results were useful in revealing the structural changes in the monomeric spike-ferritin complex (Supplementary Figure S8). PCA, which effectively reduces the dimensionality, was performed to identify the most significant fluctuation modes of the complex. Notably, the percentage of variance captured by PC1 and PC2 of the system amounted to 73.23% and 13.88%, respectively. These findings directly demonstrated that structural fluctuations are noticeable in the systems.

To gain deeper insights into the interactions within the RBD of the monomeric spike and ferritin, separate RMSD and RMSF analyses were conducted for these components (Figure 7). Notably, three replicate simulations displayed rather consistent structural stability for RBD, with average RMSD values of 4.2 Å, 3.0 Å, and 3.6 Å (Figure 7A). Suggesting that when RBD is present within the full ectodomain spike monomer, it remains structurally stable, as evidenced by the RMSF analysis as well. The average RMSF values obtained for RBD were 1.3 Å, 1.3 Å, and 1.5 Å, respectively (Figure 7B).

Similarly, the examination of ferritin across three replicates yielded RMSD values averaging at 3.0 Å, 2.4 Å, and 2.5 Å (Figure 7C). These RMSD plots indicated the significant stability of ferritin when it creates a complex with the full ectodomain spike monomer.

To gain insights into crucial changes in protein structure during the molecular dynamics simulations, the secondary structures were calculated for the trajectories (Figure 8). The secondary structure analysis of the monomeric spike-ferritin complex did not show any significant changes in its dynamics, indicating overall stability of the complex.

In summary, it was demonstrated that the RBD and ferritin components within the monomeric spike-ferritin complex maintain strong stability. These findings have significant implications, suggesting that the full ectodomain spike monomer, despite its substantial size, holds promise for vaccine development against SARS-CoV-2.

In the analysis of spike-ferritin interactions, the cross-correlation matrices of the Cα atoms were calculated in all simulations (Figure 9). Positive and negative values were used to denote correlated and anti-correlated motions, respectively. Anti-correlated motions occurring between spike and ferritin supported the stable interaction between the two proteins. It was also observed that when the spike becomes compact, the anti-correlation between RBD and ferritin grows as they come closer to each other (Figure 9B, C; Figure 10E, F).

Furthermore, visualizations of the dynamics indicated that the observed anti-correlated motions would not disturb the ‘‘up’’ active conformation of RBD (Figure 10). This is crucial because RBD needs to remain in an upward orientation to be accessible to ACE2.

The distance between the center of masses of ferritin and specific regions of the spike protein, including RBD, S1, and S2, was calculated over three replicates of 300ns MD simulations (Figure 11). In general, a descending trend was evident in the calculated distances during the simulation. Among calculated distances, S2-ferritin exhibited the least variation, while RBD-ferritin showed the greatest change in the overall trend. The observed decline in the distance was supported by the visualization of dynamics, indicating that the spike becomes more compact, with the S1 and S2 subunits of the spike approaching each other. However, the RBD maintains its upward conformation consistently throughout the entire simulation time (Figure 10).

Solvent Accessible Surface Area (SASA) calculations were performed on the RBD to track changes in its exposure and accessibility over time (Figure 12). The stable trend observed in the RBD SASA values suggested that despite structural changes in the spike protein, the RBD maintains its accessibility to the ACE2 receptor. This indicates that the RBD’s binding ability to the ACE2 receptor remains relatively constant, even during the compaction that happens to the spike protein.

The immunogenic profile of the monomeric spike nanoparticle vaccine candidate was obtained from the C-IMMSIM server. The vaccine was injected in three doses over 2 months with 30-day intervals (Figure 13). It was shown that the modeled vaccine was able to elicit both humoral and cellular mediated immune responses. Following each vaccine dose, there was an increase in B cell population, indicating that the vaccine effectively stimulated the immune system and led to a more robust antibody response (Figure 13A). Additionally, the CD4⁺ T cell population peaked after the third injection (Figure 13B). Investigation of the CD8⁺ T cell population revealed a higher population of active cells compared to anergic cells (Figure 13C). Furthermore, the plasma cell population revealed a substantial increase (>80 cells per mm) following the third dose, suggesting robust antibody production (Figure 13D). The antibody response to the monomeric nanoparticle vaccine model (Figure 13E) exhibited a slight increase in the graphs after the first two doses. From the second dose onwards, there was a rise in antibody secretions attributed to memory T lymphocytes. An upward trend of IgM and IgG antibody titer was observed after the third injection, alongside a decline in antigen levels. Additionally, elevated levels of cytokines such as IFN-γ and IL-2 were observed (Figure 13F). Clearance of antigenic molecules was demonstrated after three doses of vaccination. There was also a significant increase in the populations of both B and T memory cells, which remained elevated for several months.

Upon analysis of two additional simulations conducted with different random seeds (Supplementary Figures S5, S6), it became evident that the immune responses observed across all three simulations were comparable. Minor variations were noted, with indications of improved immune responses observed in some cases. Importantly, the pattern of immune system stimulation remained consistent across all three simulations.

The adaptive immune response to SARS-CoV-2 relies on the activity of T cells, wherein MHC molecules play a crucial role in processing and presenting antigenic peptides to T cells, initiating immune responses (Lin et al., 2023). Human MHC molecules, also referred to as Human Leukocyte Antigen (HLA), are divided into two classes. MHC class I (HLA-A, HLA-B, and HLA-C) present peptides from endogenous sources to CD8⁺ T cells, while MHC class II molecules (HLA-DR, HLA-DP, and HLA-DQ) predominantly process exogenous peptides, presenting them to CD4⁺ T cells (Traherne, 2008). The diversity in HLA alleles contributes to variations in peptide processing and presentation, thereby influencing immune responses (Wolday et al., 2023).

In order to ensure the effectiveness of the designed vaccine against different variants of HLA, the prediction of T cell epitopes on the monomeric spike protein of the vaccine was performed. However, due to the high diversity in HLA alleles, assessing all variants was not feasible. Instead, the study focused on investigating the two most polymorphic MHC genes. The HLA-B gene, with 1,077 reported alleles, is recognized as the most polymorphic among MHC class I genes, while the HLA-DRB1 gene, with 669 alleles, stands out as the most polymorphic among MHC class II genes (Shiina et al., 2009; Aliseychik et al., 2018). The analysis of monomeric spike peptides processed by HLA-B variants revealed notable findings. For each of the examined variants, 3,620 peptides were predicted. These peptides were then assessed based on their affinity for binding to HLA-B molecules. The prob score distribution of monomeric spike peptides processed by five common HLA-B variants revealed distinct patterns (Supplementary Figure S9). Overall, variant HLA-B*57:01 exhibited the highest prob score distribution, with three peptides scoring as high as 1, indicating a strong probability of binding to HLA-B molecules. Following closely were variants HLA-B*53:01 and HLA-B*44:03, which also showed relatively high prob score distributions. In contrast, the low prob score distribution and median for variants HLA-B*35:01 and HLA-B*07:02 suggested a lower probability of binding processed peptides to HLA-B molecules. Remarkably, despite the diversity in HLA-B alleles, peptides with a high probability of binding to HLA-B molecules (>70%) were identified for all of five most common HLA-B variants. Further exploration involved predicting peptides that bind to MHC class II variants. A thorough examination assessed the binding affinity of predicted peptides across 13 alleles of HLA-DRB1 and its paralogues: HLA-DRB3, HLA-DRB4, and HLA-DRB5 (Supplementary Figure S10). Remarkably, a consistent pattern was observed across all alleles. While the majority of predicted peptides exhibited a binding affinity distribution ranging from 0 to 0.1, it is noteworthy that all alleles (except for HLA-DRB4*01:01) presented peptides with a high probability (>60%) of binding to MHC class II molecules. These observations underscore the robustness of the immune system’s response to the monomeric nanoparticle vaccine across both MHC class I and class II pathways, despite variations in MHC alleles.