Epitope-Based Potential Vaccine Candidate for Humoral and Cell-Mediated Immunity to Combat Severe Acute Respiratory Syndrome Coronavirus 2 Pandemic.

The emergence of severe acute respiratory syndrome from novel Coronavirus (SARS-CoV-2) has put an immense pressure worldwide where vaccination is believed to be an efficient way for developing hard immunity. Herein, we employ immunoinformatic tools to identify B-cell, T-cell epitopes associated with the spike protein of SARS-CoV-2, which is important for genome release. The results showed that the highly immunogenic epitopes located at the stalk part are mostly conserved compared to the receptor binding domain (RDB). Further, two vaccine candidates were computationally modeled from the linear B-cell, T-cell epitopes. Molecular docking reveals the crucial interactions of the vaccines with immune-receptors, and their stability is assessed by MD simulation studies. The chimeric vaccines showed remarkable binding affinity toward the immune cell receptors computed by the MM/PBSA method. van der Waals and electrostatic interactions are found to be the dominant factors for the stability of the complexes. The molecular-level interaction obtained from this study may provide deeper insight into the process of vaccine development against the pandemic of COVID-19.

The sudden outbreak of febrile respiratory syndrome caused by a novel β-coronavirus (2019-nCoV) has created a global catastrophe where death rate is increasing every day.[1−3] The rapid propagation of novel-coronavirus or SARS-CoV-2 has outcompeted the past epidemics caused by Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV) and Middle East Respiratory Syndrome Coronavirus (MERS-CoV),[4] and therefore, it has been declared as a global pandemic by the WHO. As of now, more than 29 million confirmed cases and 924 000 deaths are caused by SARS-CoV-2. The genome of SARS-CoV-2 is a single-stranded positive-sense RNA and reported to be largest viral genome until date.[5,6] The genome entry in the host cell is guided by crown-shaped glycosylated spike protein (S-protein) located at the envelope surface. The S-protein is a trimeric class I fusion protein, which exists in a metastable prefusion conformation and undergoes a dramatic structural rearrangement to fuse into the host cell membrane.[7−9] The ectodomain of S protein includes the receptor binding S1-sub unit and the trimeric membrane fusion S2 or stalk domain.[10−12] The S1 subunit receptor-binding domain (RBD) specifically interacts with the host receptor known as angiotensin-converting enzyme 2 or ACE 2.[13] At the moment of binding of the S1 domain to a host-cell receptor, the perfusion trimer is destabilized, which results in shedding of the S1 subunit, and transition of the S2 subunit to a stable post fusion conformation.[7,10,14] The study of crucial activity of the S-protein can therefore provide a breakthrough in vaccine design and development compared with other structural proteins of SARS-CoV-2 such as nucleocapsid phosphoprotein (N-protein) or envelop protein (E-protein).[15]

Many efforts have been made for the discovery of antiviral drugs against SARS-CoV-2,[16−19] but no such licensed therapeutic are available in the market until date. Therefore, the development of an effective treatment strategy for the pandemic is a research priority. Moreover, the design of a novel vaccine against viruses using kits and related antibodies is time-consuming and expensive.[20,21] Previously, numerous approaches including the whole virus, viral-DNA, subunit, and virus-like particles have been employed in developing epitope-based vaccines against SARS and MERS virus.[22−27] These epitopes can be prepared by chemical synthesis techniques and are easier in quality control.[28,29] There are evidence which supports that in silco predictions are helpful in successful production of commercially important vaccines.[30] However, the structural modifications, delivery systems, and adjuvants are the additional requirements in the formulation because of low immunogenicity caused by their structural complexity and low molecular weight.[31] Recently, a set of B and T cell epitopes from the highly conserved region in SARS-CoV-2 were identified, which may help in developing vaccine candidates.[32−34] However, very little is known about the dynamic stability and affinity of the predicted epitope toward the interacting domain of antibody and T-cell receptors (TCRs), which is crucial for validating and improving the efficacy of predicted vaccines. In this respect, we apply a combination of immuno-informatic[10,35] approach to identify potent epitopes to design the vaccine candidates followed by computational chemistry analysis to check their effectiveness. With the help of molecular docking, MD simulations and free-energy calculations, an analysis of all the important interactions necessary to give stability to the immuno-receptor complexes have been performed.

T-cells are known to recognize and activate defense responses against viral infection; B-cells on the other hand can have antibody reactions which help in recovering extreme respiratory infection. Therefore, we have done a detailed analysis of the viral antigens to predict B-cell, T-cell linear epitopes located at the S-protein of SARS-CoV-2, evaluated their immunogenicity, and designed chimeric vaccines. The conservation of all B- and T- cell epitopes were assessed across most of the isolates and coronavirus species from different parts of the globe. Furthermore, we carried out in silico cloning of the linear vaccine construct to design a recombinant plasmid that can help in expressing the vaccines in E. coli expression system.

The methodology of the exploration of potential epitopes from the S-protein, vaccine design and their validation are discussed in Supporting Information 1 (SI-1). All protein molecules are protonated at the biological pH of 7.0. Then the chimeric vaccine and peptide epitopes are docked at the antigen binding domain of respective immune-receptors by Hex software. We carried out atomistic MD simulations of all the systems using GROMACS-2016.5.[36−38] Amber 99-SB[39] force-field was employed due to its better balance in β-sheet and helicity propensity compared to other force fields.[40−42] SPC/E[43] water molecules were used to solvate the receptor-peptide systems because of their compatibility with AMBER force field. All the systems underwent a 50 000 step energy minimization by steepest-descent algorithm[44,45] to remove the steric clash. Leapfrog integrator algorithm was used to integrate the equation of motions with a time step of 2 fs. LINCS[46] algorithm was applied to constrain all the bonds in the peptide molecule, and SETTLE algorithm was employed to constrain the geometry of water molecules. The systems were equilibrated in canonical ensemble (NVT) followed by isothermal–isobaric ensemble (NPT) for 5 and 10 ns respectively by restraining the solute heavy atoms. Next, the restraints were removed, and the protein molecules were allowed to move freely during the production run for 100 ns. The temperature and pressure of the system was maintained employing Velocity rescale[47] (τt = 0.1 ps) and Parrinello–Rahman coupling algorithm[48] (τp = 0.2 ps). The cutoff for short-range electrostatic and van der Waals interactions was assigned to 1.2 nm for the particle-mesh Ewald method.[49,50]

The amino-acid sequences of the spike-proteins for twenty-six different coronavirus species (SI-1) were considered for multiple sequence analysis (MSA) to figure out the conserved amino-acid region and the variable region which differentiate the SARS-CoV-2 from other classes of coronavirus. MSA reveals that the C-terminal end is highly conserved compared to the N-terminal end of the amino acid sequence of spike protein retrieved from different coronavirus species (Figure S1). This C-terminal perfusion or shaft domain (S2-domian) is thus common in all coronavirus spices for genome transfer while the receptor binding domain (RBD or S1-domain) is unique in SARS-CoVid2 (Figure A). The phylogenic tree constructed from MSA indicates that the spike protein of β-coronavirus species (SARS-CoVid, MERS-CoVid) has maximum structural similarity with the γ-coronavirus (Figure S-2). The antibody-mediated defense responses or humoral immunity against viral infection is believed to be guided by B-cell epitopes.[51] The linear B-cell epitopes of S-glycoprotein sequence (accession no: MT-415323) were predicted by a set of physicochemical parameters, such as the exposed surface propensity by Emini et al.,[52] hydrophilicity by Parker et al.,[53] flexibility by Karplus et al.,[54] antigenic propensity by Kolaskar et al.,[55] β-turn propensity by Chou et al.,[56] and so on, implemented in the immune epitope database and analysis resource (IEDB) tool[57] and are documented in Tables S1–S9 and in S10. The predicted B-cell and T-cell epitopes with high antigenicity scores are shown in Table .

Figure 1

(A) The structure of spike glycoprotein in trimeric form. (B) The location of the linear B-cell epitope on the chain B of spike glycoprotein. The yellow, red, pink, green, orange, and firebrick colored regions are predicted by Emini surface accessibility, Parker hydrophilicity, Karplus and Schulz flexibility, Kolaskar and Tongaonkar antigenicity and Chou–Fasman method, respectively. (C) The location of the conformational B-cell epitopes on the chain B of spike glycoprotein in their open state and closed state. The residues of conformational B-cell epitope are indicating in yellow color. (D) The R1, R2, and R3 are indicating the three different conformational epitope regions in the single chain of spike protein. The location of T-cell peptide epitopes corresponds to (E) MHC-I and (F) MHC-II. The brown, pink, olive, green, red, and orange color corresponds to the peptide associated with A*02:01, A*24:02, B*40:01, B*58:01, DRB1*04:01 and DRB1*07:01 respectively.

Table 1

Highly Immunogenic Linear B-Cell and T-Cell Epitope from the Spike Glycoprotein of SARS-CoVid2

linear B-cell epitope
methods	peptides	start_position	end_position	hotspot residues	score	length
Emini surface accessibility	DPSKPSKRSF	808	817	P, S	5.67	10
Parker hydrophilicity	DSTECSN	745	751	E	6.34	7
	GTNTSNQ	601	607	T	6.09	7
Karplus and Schulz flexibility	EGKQGNF	180	186	K, Q	1.1	7
	PGQTGKI	412	418	T	1.1	7
BepiPred linear epitope	ILPDPSKPSKRS	805	816	P, S, K, P, S	1.67, 2.06, 2.29, 2.16,	12
	QTQTNSPRRRARSV	675	687	T, N, S	1.51, 1.69, 1.54	13
Kolaskar and Tongaonkar antigenicity	FLVLLPLVSSQCVNL	4	18	L, L, P, L, V	1.243, 1.239, 1.22, 1.261, 1.227	15
	PHGVVFLHVTYVPA	1057	1070	F, L	1.21, 1.21	14
Chou–Fasman beta-turn	PGTNTSN	600	606	N	1.36	7
	KGCCSCG	1245	1251	C	1.3	7

The size of the predicted linear B-cell epitopes is found to be varied between 7mer- 15mer which is believed to be of optimum size. The hot-spot residues which are considered as the key to the antigenicity of the epitopes, are found to be mostly polar in nature which may contribute in stable hydrogen bonds to the residues present in the antigen binding region of the antibodies. In Figure B, we showed the location of linear B-cell epitopes on the single chain of SARS-CoV-2 spike protein. It is evident from Figure B that most of the epitopes are located at the RDB domain and the shaft region of the S-protein. The peptide 808DPSKPSKRSF817 and 600PGTNTSNQ607 are found to be common in Emini-surface accessibility, Bepipred linear epitope and Chou–Fasman methods indicating their greater importance as linear B-cell epitopes.

The 3-D structure of the spike protein was generated by homology modeling implemented in Swiss Model[58] using PDB ID: 6VIB and 6VXX as template. The conformational epitopes of the spike protein of SARS-CoV-2 in their closed and open state were determined by Disctope[59,60] (V-1.1) server and are documented in Tables S11, S12. It can be seen that the conformational B-cell epitopes of spike protein are located almost in the same position (RDB region) in their both closed and open state (Figure S3A,B). The locations of the conformational epitopes on the tertiary structure of S-protein are mainly found at 405–427 and 439–505 residue stretches and shown in Figure C. Based on the location of the epitopes, the domains are marked as R1, R2, and R3. It can be noted that in both closed and open state, all antibody binding regions (R1, R2, and R3) have extended antiparallel β-sheet or β-barrel structure which is an important factor for antigenicity (Figure -D).

The identification of CD8+ and CD4+ T-cell epitopes are known to be important for eliciting cell-mediated immunity or generation of memory B-cell against viral infections.[61] Such peptide epitopes are generally presented by the Major Histocompatibility complex (MHC) -class I and MHC-class II molecules expressed on the surface of helper T-cell (TH).[61,62] We employed IEDB server to identify and validate the T-cell peptide epitopes of the S-protein of SARS-CoV-2 which can bind with MHC-I and MHC-II receptors. All MHC-I and based epitope-peptides obtained from IEDB T-cell epitope prediction tools are documented in Tables S13–S16 and that of MHC-II are documented in Tables S17 and S18. The peptides with highest antigenicity and high affinity are listed in Table , and their location at the S-protein of SARS-CoV-2 is depicted in Figure E,F. The size of the peptides for MHC binding are found to vary from 8mer to 15 mer. It is evident from Figure that peptide epitope 781VFAQVKQIY789 corresponding to MHC-I and 960NTLVKQLSSNFGA972 associated with MHC-II are located at the stem part of the S-protein. Unlikely, the epitopes 567RDIADTTDAV576, 76TKRFDNPVLPF86, 747VRFPNITNL754, (related to MHC-I), and 116SLLIVNNATNVVIK129 (related to MHC-II) are found to be located at the surface (i.e., the receptor binding domain (RBD) of S-protein). The location of this peptides shares the binding domain of Angiotensin-Converting Enzyme 2 (ACE2) receptor which make them more interesting. These selected peptides are used to design multiepitope vaccine for cell mediated immunity.

The knowledge regarding the linear-epitope-conservation is necessary to design vaccines capable of inducing adaptive immune response for all coronavirus species. It is found that all the epitopes considered here contain at least two conserved residues (Figure S1). Epitopes with three identical or similar residues are considered to be highly conserved. Among all the selected B-cell epitopes in Table , 412PGQTGKI418, 745DSTECSN751, 805 ILPDPSKPSKRS 816, 808 DPSKPSKRSF 817, 1057PHGVVFLHVTYVPA1070, 1245KGCCSCG1251 are found to be highly conserved (Figure S1). In case of the T-cell epitope, 116SLLIVNNATNVVIK129, 781VFAQVKQIY785, and 960NTLVKQLSSNFGA972 are highly conserved. It is evident from Figure B that most of the highly conserved B-cell epitopes are located at the stalk or shaft region of the S-protein except 412PGQTGKI418 and 1057PHGVVFLHVTYVPA1070. Similarly, in the case of T-cell epitopes, highly conserved 116SLLIVNNATNVVIK129, 781VFAQVKQIY785 epitopes are located at the shaft of the S-protein. The rest of the epitopes are less conserved and mostly located at the RBD region of S-protein. This indicates that amino-acid composition of RDB of SARS-CoV-2 S-protein is different from other coronavirus species.

Further, we have designed two separate vaccines from the epitopes predicted in Table , and named as Vac-COVID-B, Vac-COVID-T respectively, which can be used as a combination for both humoral and cell-mediated immunity. The linear B-cell epitopes and T-cell epitopes are linked by the GPGPG linker to avoid the formation of junctional epitope[34] and Cholera Toxin B (CTB) adjuvant linked by EAAAK linker for immune regulation. The total number of amino-acids of the vaccines are 133 and 111 amino acid (aa) long with molecular weight of 13.11 and 11.4 kDa, respectively. The 3D-models of the chimeric vaccine shown in Figure are predicted to be nonallergenic by the AllerTOP server.[63]

Figure 2

Tertiary structure of chimeric vaccines made of (A) B-cell linear epitopes (B) T-cell linear epitopes. The amino acid sequence of the corresponding vaccines is shown below. The adjuvant is shown in red; the adjuvant linker is shown in green, and the epitope linkers are indicated in blue.

The structural quality of the vaccines is assessed by Z-score and Ramachandran plot (Figure S4). The Z-score of the vaccines made of B-cell and T-cell epitopes are −3.43 and −4.85, respectively, which confirm the reliability of our model. Additionally, 95.42% and 92.7% amino acid residues are found to be in the favored region in the Ramachandran plot. Further, the antigenicity of the two vaccines are found to be 0.73 and 0.58 (score >0.4 is considered to be antigenic).

Next, it will be interesting to study the interaction of the multiepitope vaccines designed with the immune-cell receptor for eliciting stable immune response. We considered the structure of the humanized antibody (7BZ5) as the immune receptor for B-cell epitope, where the antigen binding region of 7BZ5 is formed by a variable region of light chain (VL) and heavy chain (VH) (Figure S6A). The antigen binding pocket of the MHC-I molecule (2GTZ, 5WWJ, 6IEX, 5IM7) is formed by the interaction of α1 and α2 domain of α-chain (Figure S6B), whereas the antigen presenting domain of MHC-II (2SEB, 6BIY) molecules is formed by the association of α1 and β1 domains (Figure S6C). However, the peptide presenting platform is similar for both the MHC class-I and class-II molecules. The detailed interaction of the chimeric vaccine with the variable region of 7BZ5 is depicted in Figure A by molecular docking. The linear vaccine from the B-cell epitope exhibited 12 hydrogen bond interactions with the variable region of the antibody. The residues Gly38, Gly 50, Glu 41, Gly76, Gly 54, Thr56, Tyr32, Ser31 of Vac-COVID-B exhibited hydrogen bond with Ser56, Tyr58, Tyr52, Ser53, Tyr33, Tyr94 of VH and Tyr97, Ala99, Ile59 and Lys58 of VL respectively. The interaction of conformational epitopes with the 7BZ5 at the molecular level is depicted in Figure B–D. The conformational epitope R1 (Arg 102, Ser98, Phe140, and Leu242) was found to make hydrogen bond with Tyr100, Tyr32, Tyr 33, and Tyr94 of the antibody molecule (Figure B). Similarly, Gly28, Tyr32, Asn92, Tyr58, Tyr52, Tyr33 of 7BZ5 formed hydrogen bonds with Tyr449, Ser494, Tyr453, Arg403, Asp405, and Lys417 of R2 (Figure C).

Figure 3

Nonbonded interaction of the vaccines and conformational epitopes with the immune receptor. (A) The residues of the linear vaccine involved in the formation of hydrogen bond with 7BZ5. The vaccine is shown in blue color. The interaction of (B) R1 (red), (C) R2 (green), and (D) R3 (brown) with the variable region of 7BZ5. The T-cell epitope peptide presented at the peptide presenting groove of (E) HLA-A*02:01, (F) HLA-A*24:01, (G) HLA-B*40:01, (H) HLA-B*58:01, (I) DRB1*04:01, (J) DRB1*07:01. The interacting residues of the receptor are shown in violet color, whereas the residues of the vaccine or the discontinuous epitope are marked in pink color.

In addition, Tyr 94 of the antibody and Tyr 505 of R2 showed π–π interaction which is crucial for the antibody-epitope stabilization. In case of R3, no nonbonded interaction was found (Figure D) except one π–π stacking interaction between Tyr 32 and Phe562. The T-cell peptide epitopes were also found to be stabilized by the hydrogen bond formed with the peptide presenting groove of the MHC-receptors. The TCRs are known to recognize the antigens in pieces that are presented by MHC molecules.[61] Therefore, the interaction of full-length Vac-COVID-T with MHC molecules is scientifically not required. The epitope 567RDIADTTDAV576 formed hydrogen bonds with Thr80, Lys146, Trp147, Tyr116, Thr73, Arg97, Tyr7, and Glu63 of HLA-A*02:01 (Figure E). 76TKRFDNPVLPF86 was found to form hydrogen bonds with Arg83, Asn77, Thr73, Ala69, Thr163, Lys146, and Trp147 of the HLA-A*24:02 peptide presenting pocket (Figure F), and in the cases of HLA*40:01, Ala150, Tyr116, Asn114, Arg62, Glu58, and Trp167 were involved in hydrogen bond formation with the 327VRFPNITNL335 peptide (Figure G). The epitope 781VFAQVKQIY789 exhibited hydrogen bonds with Ala139, Thr 143, Trp147, Tyr74, Tyr159 of HLA-B*58:01 (Figure H). Further, the MHC-II based T-cell epitope 960NTLVKQLSSNFGA972, 116SLLIVNNATNVVIK129 were found to be involved in hydrogen bond formation with Phe48, His81, Gln70, Trp62 of DRB1*-04:01 (Figure I) and Arg71, Ser53 of DRB1*07:01, respectively (Figure J).

The docking predicted interaction of the vaccine–receptor complex in motion were assessed by all atom MD-simulation studies.[64−66] In order to check the reproducibility of our MD results, we carried out another set of simulation and provided in Supporting Information 2 (Figure S20, S21) and Supporting Information 3 (Table S21). The time evolution of the RMSD, RMSF, radius of gyration (rGyr) of B-cell linear vaccine and epitopes located at R1 R2, R3 region are illustrated in Figure . The RMSD of the linear vaccine (Vac-COVID-B) was found to increased approximately 0.7 nm up to initial 8 ns and remain stabilize until 40 ns. The RMSD fluctuation and rGyr profile of Vac-COVID-B near 60–90 ns fluctuation indicate minor secondary structure modification of Vac-COVID-B during the course of simulation (Figure A,B). The RMSF profile of the designed vaccine is shown in Figure S7A. The residue stretches 39–48, 57–73 were found to highly fluctuate during the course of simulation because of the interaction with the antibody molecule. The RMSF value of the residues located at antigen binding domain of 7BZ5 is found to be less than 0.3, which helped the chimeric vaccine to get stabilized at the immune complex (Figure S8). The average number of hydrogen bonds between the antibody and the vaccine is calculated to be 9 and remain intact throughout the simulation trajectory (Figure E, Figure S9). The docking predicted residue pairs Gly54-Tyr97, Gly54-Tyr52, Thr56-Tyr33 have hydrogen bond occupancy values of 59.4%, 13.7%, 11.3%. The highest hydrogen bond occupancy of 69.7% is found between Tyr58 of the linear vaccine and Pro51 of 7BZ5, which newly evolved during the course of simulation (Table S20A). The structural stability of the conformational epitopes located at R1, R2, R3 is assessed in Figure C,D. The RMSD and rGyr profile of the conformational epitope located at R1 and R3 is found to be stable, whereas that of R2 is highly unstable which indicates the conformational change of the epitope during the course of the simulation. The C-terminal and the N-terminal end of all the epitopes have greater fluctuation because of solvent exposure (Figure S7B–D). The average number of hydrogen bonds between the epitope located at R1, R2, and R3 are 2, 7, and 6, respectively. The docking predicted residues pairs were not found to make hydrogen bonds during the course of simulation; rather, the new hydrogen bonds evolved during the course of simulation with a lesser hydrogen bond occupancy percentage (Table S20B). In the case of R2, the docking predicted residue pair Asn92-Tyr453 showed the highest hydrogen bond occupancy of 99.9% (Table S20C). The epitope located at R3 was found to make stable hydrogen bond with Asn32 and Asn92 with H-bond occupancy of 55.9% and 51.6%, respectively (Table S20D). The number of hydrogen bond profile of R1 is fluctuating compared with other conformational epitopes, whereas the number of hydrogen bond is found to be increasing between the antibody and R3 with respect to simulation time to stabilize the immune complex (Figure S9B–D). The conformational free-energy landscape of the antibody and antigen binding is depicted in Figure S10 to confirm the adequate sampling of the immune complexes. The solvent accessible surface and RMSD of the chimeric vaccine, B-cell conformational epitopes are considered as reaction coordinate. It is evident from FEL graphs that R1 is trapped in a deep minimum, which indicates the single binding conformation throughout the simulation. In the case of other antigens, two prominent binding conformations along with others are present, which evolved during the course of simulation.

Figure 4

(A) The RMSD profile of the chimeric vaccine (Vac-COVID-B) (B) The gyration radius (rGyr) profile of Vac-COVID-B throughout the simulation. (C) The time evolution of RMSD of epitopes located at R1, R2 and R3. (D) The gyration radius (rGyr) of the conformational epitope located at R1, R2, R3. (E) The average number of hydrogen bond between the antibody and vaccine molecules.

We have done MD simulation of individual peptide epitopes with their corresponding HLA compounds. The dynamic nature of the T-cell epitopes with their respective HLA-complex are assessed in the Figure . The RMSD of the peptides correspond to MHC-I increased up ∼20 ns and stabilized with an average value of 0.45 nm for the rest of the simulation. The peptide associated with MHC-II is found to be stable for the last 30 ns of simulation. Further, we calculated the time evolution of the distance between the peptide and the surface of the binding pocket of MHC molecules (as shown in Figure S11). The distance between the surface of the peptide epitope and the platform of antigen presenting domain was found to be stable throughout the simulation, indicating that the peptides did not diffuse from the peptide presenting grooves. The solvent assessible surface area of the peptide epitopes are found to be stabilized with their initial value, which confirms the constant solvent exposure at the peptide presenting groove (Figure S12).

Figure 5

RMSD profile of the T-cell epitopic peptides corresponding to (A) A*-02:01, A*-24:02, B*-40:01 and B*-58:01. (B) DRB1*04:01, DRB1*07:01. The RMSF profile of (C) MHC-I and (D) MHC-II based peptide epitopes. (E) The average number of hydrogen bond formed between the peptide epitope and MHC molecules during the simulation.

The peptide associated with HLA-B*40:01 and DRB1*04:01 is found to have less fluctuation in MHC-I and MHC-II based epitopes, respectively. MHC-I based peptides mostly fluctuate at their N-terminal end whereas; MHC-II based peptides have higher fluctuation at the C-terminal end (Figure C,D). The RMSF profile of the binding domain of MHC-I and MHC-II molecules are below 0.3 which helps to stabilize the peptide epitopes inside the peptide presenting pocket (Figure S13).

The average numbers of hydrogen bonds formed between MHC-I based TCRs T-cell epitopes are 4, 8, 6, and 7, respectively, for the four alleles during a 100 ns simulation (Figure E). It is evident from Figure S14 that the number of hydrogen bonds was almost constant throughout the simulation time. Hydrogen bond occupancy of 41.2% is found between Ala5 and Thr73 of HLA-A*02:01- RDIADTTDAV complex (Table S20E). The residue pair Leu10-Trp147, Phe12-Tyr84 of HLA-A*24:02- TKRFDNPVLPF complex and Pro5-Tyr159, Arg3-Glu63 of HLA-B*40:01- VRFPNITNL has high hydrogen bond occupancy which stabilizes the peptide-ligand at the peptide presenting domain (Table S20F, G). The higher hydrogen bond occupancy of Gln8-Asn77 (87.8%), Tyr10-Tyr80 (89.9%) Tyr10-Thr143 (83.0%) residue pair stabilized the peptide VFAQVKQIY at the binding domain of HLA-B*58:01 (Table S20H).MHC-II based peptides epitope were found to have higher number of hydrogen bond due to their higher size (Figure E). The residue pair Lys6-Asn82, Asn2-Ala52 showed hydrogen bond occupancy of 95.0%, 87.7% for the peptide associated with DRB1*04:01, whereas Leu4-Asn82, Ser2-His81, Leu3-Ser53 of DRB1*07:01-peptide complex showed high hydrogen bond occupancy (Table S20I, J). This above-mentioned hydrogen bond between the T-cell epitope peptide and HLA-molecules stabilized the complex to transduce stable immune response against SRAS-CoV-2.

The free-energy landscape of peptide epitope binding at the peptide presenting pocket of MHC molecules is shown in Figure S15. It is evident that all the peptides crossed a higher free-energy barrier and stabilized at deep minima at the higher RMSD values.

Lastly, the dynamic stability of the designed chimeric vaccines at aqueous solvent were assessed by MD simulation studies (Figure S16). The time evolution of RMSD, radius of gyration, and number of native hydrogen bonds were calculated to check their conformation in motion. The RMSD of both the vaccines were found to be converged with an average value of 0.6 nm. It is evident from the gyration radius profile (rGyr) of Vac-COVID-B that compactness of the protein is decreased and stabilized with an average value of 1.5 nm. The compactness of Vac-COVID-T is found to be constant throughout the simulation. Further, the number of native hydrogen bond is found to be constant throughout the simulation trajectory which indicates the conformation stability of the vaccines in polar solvent.

The free energy of binding (ΔGbind) is believed to be an important thermodynamic quantity to assess the favorable protein–protein interaction as well as their affinity for accurate modeling of biological systems. In this regard, we calculated the binding free energy of the epitopes at the binding domain of immune cell receptors in implicit solvent model by the end-point free-energy method such as MM/PBSA[67] and summarized in Table . The designed linear vaccine showed remarkably high binding affinity of −453.59 kJ/mol with the variable region of 7BZ5 which confirms the thermodynamic stability of the complex. The ΔGelec term (−469.00 kJ/mol) between the residues of vaccine and the humanized antibody was found to have important contribution toward the stability of the complex (Table). The affinity of the epitope located at R1 is the lowest compared with other discontinuous epitopes because of favorable van der-Waals energy (−198.54 kJ/mol). The lower affinity of R2 and R3 is due to higher ΔGsol penalty and positive electrostatic interactions respectively. It is found that both the MHC-II based peptides 960NTLVKQLSSNFGA972 and 116SLLIVNNATNVVIK129 exhibit the lowest binding energy of −493.66 kJ/mol and −538.71 kJ/mol, respectively. Among the MHC-I peptide epitopes, 76TKRFDNPVLPF86 associated with HLA-A*24:02 has the lowest binding energy of −430.3 kJ/mol. All the peptide epitopes related to MHC-molecules have considerably low Gibbs free energy, indicating the stable epitope-TCR complex. It is evident that the electrostatic and van der Waals free-energy terms of the T-cell peptide epitopes are crucial for their stability. The hydrogen bond observed during our MD study is mainly responsible for favorable electrostatic energy contribution, which is dominating the positive solvation free-energy term for the stable interaction of the epitopes at the immune cell receptors.

Table 2

Binding Affinities (kJ/mol) of the Vaccines toward the Immune Cell Receptors by MM/PBSA Methoda

immune receptor	complexes	ΔGvdw	ΔGelec	ΔGsol	ΔGSASA	ΔGbind
antibody	Ab-VAC-COVID B	–145.07	–469.00	345.5	–158.02	–453.59
	Ab-R1	–198.54	–80.89	73.3	–12.27	–218.4
	Ab-R2	–111.79	–41.15	103.12	–28.95	–78.77
	Ab-R3	–75.258	27.69	37.84	–8.76	–18.52
TCR (MHC-I)	HLA-A*-02:01-RDIADTTDAV	–199.89	–365.33	367.24	–27.98	–225.95
	HLA-A*-24:02-TKRFDNPVLPF	–340.64	–826.49	779.80	–43.00	–430.3
	HLA-B*-40:01- VRFPNITNL	–328.14	–528.54	602.63	–40.70	–294.76
	HLA-B*-58:01-VFAQVKQIY	–176.52	–384.67	368.43	–20.62	–213. 39
TCR-(MHC-II)	HLA-DRB1*04:01- NTLVKQLSSNFGA	–272.99	–652.81	463.34	–31.19	–493.66
	HLA-DRB1*07:01- SLLIVNNATNVVIK	–339.73	–608.54	450.9	–41.33	–538.71

The ΔGelec, ΔGvdw, ΔGsol, ΔGSASA and are indicating the electrostatic, van der Waals, polar solvation, solvent accessible surface energies respectively.

Further, we calculated the free-energy contribution of residues located at the antigen binding domain of immune cell receptors (Figures S17 and S18). The docking predicted residues and residues with higher hydrogen bond occupancy percentage was found to contribute lower energies which stabilize the immune complexes. In the case of chimeric vaccine (Vac-COVID-B), Ile 51 of VH and Pro 51 of VL has maximum free-energy contribution. The amino acid residues with benzene ring (Phe 27, Tyr 32, Tyr 33, Tyr 58, Tyr 94, Tyr 100) have higher energy contribution because of stacking interactions (Figure S17). In the case of T-cell epitopes, the hydrophilic or polar amino acids located at the peptide presenting groove have higher contribution to stabilize the peptide epitopes (Figure S18)

In order to design a recombinant plasmid, we back translated the protein sequences of the vaccines and optimized the codons[68] in the E. coli system for successful expression of linear B-cell (Vac-COVID-B) and T-cell (Vac-COVID-T) vaccines proposed in our immunoinformatic study. The size of c-DNA sequences of the vaccines made from linear B-cell, T-cell epitopes are 399 base pair (bp) and 333 bp long, respectively. The codon adaptation index (CAI) value for both the vaccines were computed to be 1.0 and the percentage of GC content for Vac-COVID-B, Vac-COVID-T are 60.65%, 56.75% respectively which are in the permissible range[69] and hence confirmed their proficient expression in the E-coli K-12 strain. Finally, the c-DNA sequences of the vaccines were inserted computationally at the multiple cloning site (MCS) of the pUC-19 plasmid vector. The restriction map of the recombinant vector is shown in Figure S19.

To conclude, in the present article we tried to design a vaccine based on B-cell and T-cell epitopes present in spike glycoprotein of highly infective SARS-CoV-2. With the help of immunoinformatic studies, we identified the most promising epitopes which are found to be scattered on the RBD and shaft region of the spike protein. The epitopes located in the shaft region are highly conserved. The hotspot residues that are considered as key to the antigenicity are mostly found to be polar in nature, which contributes to stable electrostatic interaction with respective immune receptors. Docking calculations showed the major interaction between the immune–receptor complexes were hydrogen bond and π–π stacking interactions which were found to contribute maximum in the stability, as evident from free-energy decomposition studies. The hydrogen bond evolved during MD simulation studies between the antigen and receptors was found to be the main contributor of the electrostatic energy. Vaccines designed from linear B-cell epitopes were found to exhibit a higher number of hydrogen bonds at the binding domain of the antibody compared with the conformational epitopes. All peptide epitopes corresponding to MHC-I and MHC-II showed remarkable stability because of the van der Waals and the electrostatic energy terms. The present article therefore provides deeper biophysical insights toward the stabilization of predicted vaccine candidates with immune cell receptors which will be helpful in further experimental design of potential vaccine against SARS-CoV-2.