Design of multi-epitope vaccine candidate against SARS-CoV-2: a in-silico study.

The best therapeutic strategy to find an effective vaccine against SARS-CoV-2 is to explore the target structural protein. In the present study, a novel multi-epitope vaccine is designed using in silico tools that potentially trigger both CD4 and CD8 T-cell immune responses against the novel Coronavirus. The vaccine candidate was designed using B and T-cell epitopes that can act as an immunogen and elicits immune response in the host system. NCBI was used for the retrieval of surface spike glycoprotein, of novel corona virus (SARS-CoV-2) strains. VaxiJen server screens the most important immunogen of all the proteins and IEDB server gives the prediction and analysis of B and T cell epitopes. Final vaccine construct was designed in silico composed of 425 amino acids including the 50S ribosomal protein adjuvant and the construct was computationally validated in terms of antigenicity, allergenicity and stability on considering all critical parameters into consideration. The results subjected to the modeling and docking studies of vaccine were validated. Molecular docking study revealed the protein-protein binding interactions between the vaccine construct and TLR-3 immune receptor. The MD simulations confirmed stability of the binding pose. The immune simulation results showed significant response for immune cells. The findings of the study confirmed that the final vaccine construct of chimeric peptide could able to enhance the immune response against nCoV-19.

Introduction

COVID-19 pandemic is result of the infection caused by Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) and it attacks the vital organs of body and targets pneumocytes in lungs which leading to fatal respiratory distress (Elfiky, 2020; Galante et al., 2016; Joshi et al., 2020; Tse et al., 2004; Yi et al., 2020). SARS, MERS, and SARS-CoV-2 caused diseases are characterized by a lower respiratory ailment like bronchitis, bronchiolitis, and pneumonia (Bogoch et al., 2020; Elfiky & Azzam, 2020). Coronaviruses are enveloped, single stranded RNA viruses which bears club shaped glycoproteins on their surface (Prajapat et al., 2020; Sarma et al., 2020). Corona virus is contagious and spreads through inhalation, ingestion of viral droplets resulting from coughing and sneezing. The coronavirus genome is comprised of ∼30000 nucleotides and it encodes with four structural proteins, Nucleocapsid (N) protein, Membrane (M) protein, Spike (S) protein and Envelop (E) protein and several non-structural proteins (nsp) (Boopathi et al., 2020; Gupta et al., 2020; Hasan et al., 2020; Khan et al., 2020). The non-structural protein 16 (nsp16) or 2′-OMTase is the crucial protein responsible for viral replication and expression in host cells (Khan et al., 2020). The SARS-CoV-2 has a surface spike like glycoprotein (S-Spike), it binds specifically to angiotensin-converting enzyme 2 (ACE2) to access host type-2 pneumocyte cells that could bind to human cells (Hoffmann et al., 2020; Verdecchia et al., 2020). The enveloped virus covered with spike surface glycoprotein and after binding the host cell is engulfed in to the cell and releases positive sense single stranded RNA in to type-2 pneumocyte. RNA dependent RNA polymerases and proteinases makeup the components of viral proteins like nucleocapsid spike proteins and enzymes necessary for viral replication and damage to host cells. Type-2 pneumocytes that produces a surfactant molecule and decreases surface tension in the alveoli and reduce the collapsing pressure (Belouzard et al., 2012; Weiss & Navas-Martin, 2005). However, at present researchers are aiming to uncover this spike protein processing proteases for early drug development using various approved inhibitor drugs/other compounds (Aanouz et al., 2020; Elmezayen et al., 2020; Enmozhi et al., 2020; Islam et al., 2020; Muralidharan et al., 2020; Pant et al., 2020; Sinha et al., 2020; Wahedi et al., 2020). The viral proteins encounters immune cells and release specific cytokines leads to vasodilation, capillary permeability, aleveolar edema and finally increase the collapsing pressure to burst out the pneumocyte (van de Veerdonk et al., 2020). SARS-CoV-2 surface spike (S) protein contains two subunits S1 and S2, of which S1 is sole responsible for host cell receptor and S2 harbors the membrane fusion machinery. Spike protein from SARS-CoV-2 shares the high structural similarity with SARS-CoV spike (Weiss & Navas-Martin, 2005). The vaccine development against SARS-CoV-2 spike protein is an important approach and hence, the present study is focused on epitope prediction analysis for construction of vaccine candidate by computational methods.

Materials and methods

CTL epitope prediction

The complete amino acid sequence of spike glycoprotein (NCBI Accession id = “QHD43416.1) of SARS-CoV-2 strain was retrieved from NCBI. Cytotoxic T-cell lymphocyte (CTL) epitopes were predicted for spike glycoprotein using artificial neural network algorithm based online server NetCTL 1.2 which predicts the MHC-class I binding, and then followed by submitting predicted NetCTL generated results to VaxiJen v2.0 (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html), ToxinPred (http://crdd.osdd.net/raghava/toxinpred/) servers to predict protective nontoxic antigens (Larsen et al., 2007). After the screening of epitopes with VaxiJen v2.0 and ToxinPred servers, the resultant epitopes were subjected to immunogenicity prediction using IEDB server (https://www.iedb.org/).

HTL and linear epitope prediction

IEDB server was used to predict the MHC-II restricted epitopes as it uses special patterns for HLA-DRB1*01:01, HLA-DPA1*01/DPB1*04:01, HLA-DQA1*03:01/DQB1*03:02 alleles, further, T-helper 1-type immune response activation and IFN-γ production was predicted using IFNepitope server (http://crdd.osdd.net/raghava/ifnepitope/). Toxicity was predicted using ToxinPred server (http://crdd.osdd.net/raghava/toxinpred/) and BCPRED 2.0 online server was used to predict the linear B-cell epitopes of spike protein (http://crdd.osdd.net/raghava/bcepred/). The results of linear B-cell epitopes and HTL epitopes of overlapping regions were assembled and considered as final predicted epitopes.

Construction of multi-epitope vaccine sequence

The potential non toxic and probable antigenic vaccine was constructed using selected epitopes. The linear B-cell and HTL epitopes were joined with GPGPG linker peptides and CTL epitopes were joined using AAY linker. The N- terminal position of the vaccine construct was linked with the sequence of 125 amino acid residue 50S ribosomal L7/L12 peptide which acts as an adjuvant and C-terminal portion was linked with HHHHHH (6HIS) linker. The basic property of allergenicity was assessed using an online server, AllerTOP v. 2.0 (https://www.ddg-pharmfac.net/AllerTOP/) and then submitted to VaxiJen server to find out whether the construct vaccine could be a probable antigen to elicit immune response.

Homology modeling and molecular docking

PSIPRED web tool was employed for finding the secondary structure analysis (http://bioinf.cs.ucl.ac.uk/psipred/). Homology modeling of the vaccine protein tertiary structure was performed using fully automated protein modeling I-TASSER server (https://zhanglab.ccmb.med.umich.edu/I-TASSER/) and the best model was selected and then optimized using SPDB viewer (https://spdbv.vital-it.ch/). Loop refinement was done using ModLoop (https://modbase.compbio.ucsf.edu/modloop/). Ramachandran plot and ERRAT server (https://servicesn.mbi.ucla.edu/ERRAT/) analysis were performed for further validation study. The designed vaccine candidate was subjected to docking using GRAMM-X Simulation web server (http://vakser.compbio.ku.edu/resources/gramm/grammx/) for docking vaccine protein model with TLR-3 (PDB ID: 1ziw) and interaction were visualized using LIGPLUS 1.2 software (https://www.ebi.ac.uk/thornton-srv/software/LIGPLOT/).

Molecular dynamic simulations

The protein & receptor complex was subjected to MD simulations. The MD simulations were done by GROMACS 2018 package to carry out 20 ns simulations using OPLS force field. The TIP3P water model was selected for solvating complexes followed by the addition of ions to neutralize. Periodic boundary conditions were used and Equilibration of the system was done using NVT and NPT ensemble for 100 ps. The trajectory was set to be generated every 2 fs and save every 2 ps. The protein-protein complex result was then analyzed (Enayatkhani et al., 2020).

Immune simulation

The immune response profile of vaccine construct was recorded by in silco method C-ImmSim, online simulation server (http://150.146.2.1/C-IMMSIM/index.php). The C-ImmSim model describes both humoral and cellular response of a mammalian immune system against vaccine construct. The target product profile of a prophylactic vaccine, three injections were given at different intervals of four weeks. The time step of simulation corresponds few hours of real life and simulation was performed with default parameters. The sequence of injections is Ag1, Ag2, Ag3 were administered four weeks apart. The simulation volume and simulation steps were set at 1000, (random seed = 12345 with an injection of vaccine containing no LPS.

Analysis of vaccine construct and in silico cloning

The physiochemical properties of the vaccine construct was assessed using online web tool the ProtParam (https://web.expasy.org/protparam/), where as the solubility, allergenicity and probable antigenic prediction were performed using protsol, AllerTOP v. 2.0 and VaxiJen v2.0 online web servers. J-CAT tool (http://www.jcat.de/) is used for codon optimization of vaccine construct and E.coli (K12) strain is selected as source organism. SnapGene software (https://www.snapgene.com/try-snapgene/) was used for in silico cloning of vaccine construct in to pET-28a vector.

Results and discussion

The amino acid sequence was used to predict the possible probable antigenic epitopes of linear B-cell, HTL and CTL epitopes for designing the multi-epitope vaccine. The vaccine construct consisted of 425 amino acid residues derived from different peptide sequences. CTL epitopes of 9-mer lengths were predicted using NetCTL1.2 (Table 1). Based on high binding affinity score, the results were submitted to VaxiJen v2.0 and predicted the 16 protective probable antigens. The non antigenic epitopes were removed and subjected to predict the toxicity using ToxinPred and after removing two toxin epitopes 14 non-allergenic epitopes were selected using toxinpred and, the IEDB immunogenicity server produced the results of seven epitopes and were given in the Table 2. The predicted probable antigenic HTL epitopes were selected for further screening of toxigenicity prediction using vaxigen 2.0 server and 13 HTL epitopes were selected and further classified as non-toxins using Toxinpred server (Tables 3 and 4). The final HTL epitopes were selected as a result of IFN-γ inducing epitopes (Table 5). The linear B-cell epitopes were used in vaccine construct as overlapping B-cell and T-cell epitopes.

Table 1.

List of predicted T-cell epitopes on the basis of C-terminal cleavage and TAP scores.

Residue No	Peptide Sequence	MHC binding affinity	Rescale binding affinity	c-terminal cleavage affinity	Transport affinity	Prediction score
865	LTDEMIAQY	0.7953	3.3768	0.9723	2.7790	3.6616
258	WTAGAAAYY	0.6735	2.8596	0.7339	2.8630	3.1128
604	TSNQVAVLY	0.6559	2.7847	0.9440	2.9910	3.0758
361	CVADYSVLY	0.5348	2.2705	0.9764	3.1800	2.5759
733	KTSVDCTMY	0.4908	2.0840	0.9649	3.0160	2.3795
746	STECSNLLL	0.5136	2.1808	0.8879	0.7030	2.3492
652	GAEHVNNSY	0.4042	1.7163	0.9769	2.6630	1.9960
196	NIDGYFKIY	0.3921	1.6649	0.9664	3.0150	1.9606
160	YSSANNCTF	0.3975	1.6878	0.9032	2.5980	1.9531
152	WMESEFRVY	0.3902	1.6569	0.7993	2.9290	1.9232
162	SANNCTFEY	0.3737	1.5865	0.9196	2.9900	1.8739
687	VASQSIIAY	0.3529	1.4986	0.9656	3.0890	1.7978
30	NSFTRGVYY	0.3389	1.4389	0.6421	3.1240	1.6915
136	CNDPFLGVY	0.2613	1.1095	0.6900	2.4500	1.3355
392	FTNVYADSF	0.2704	1.1480	0.3800	2.3170	1.3208
261	GAAAYYVGY	0.2253	0.9568	0.7608	2.9690	1.2194
357	RISNCVADY	0.2106	0.8941	0.9292	3.3940	1.2032
465	ERDISTEIY	0.2097	0.8903	0.9744	2.6460	1.1687
285	ITDAVDCAL	0.2350	0.9979	0.8708	0.7900	1.1680
1039	RVDFCGKGY	0.2036	0.8644	0.7618	3.2320	1.1403
343	NATRFASVY	0.1955	0.8300	0.9342	2.8730	1.1138
1237	MTSCCSCLK	0.2260	0.9595	0.7525	0.4790	1.0963
50	STQDLFLPF	0.1974	0.8383	0.5530	2.5110	1.0468
1096	VSNGTHWFV	0.2012	0.8544	0.6143	0.2180	0.9574
880	GTITSGWTF	0.1656	0.7031	0.7489	2.5570	0.9433
815	RSFIEDLLF	0.1421	0.6035	0.5938	3.0320	0.8441
1264	VLKGVKLHY	0.1262	0.5356	0.9783	2.8590	0.8253
748	ECSNLLLQY	0.1413	0.6000	0.5316	2.7470	0.8171
370	NSASFSTFK	0.1671	0.7093	0.5456	0.5070	0.8165
372	ASFSTFKCY	0.1180	0.5010	0.9587	3.2750	0.8085
628	QLTPTWRVY	0.1189	0.5047	0.9661	2.7820	0.7887
296	LSETKCTLK	0.1515	0.6432	0.8919	0.2200	0.7879
192	FVFKNIDGY	0.1358	0.5767	0.4093	2.9130	0.7837
445	VGGNYNYLY	0.1164	0.4941	0.9518	2.6580	0.7698
83	VLPFNDGVY	0.1130	0.4797	0.9703	2.8460	0.7675
1095	FVSNGTHWF	0.1232	0.5231	0.7203	2.6210	0.7622
612	YQDVNCTEV	0.1531	0.6501	0.5870	0.2420	0.7502

Table 2.

List of seven predicted T-cell epitopes predicted after IEDB immunogenicity prediction server screen.

S.No.	Peptide	Length	Score
1.	QLTPTWRVY	9	0.31555
2.	VLPFNDGVY	9	0.1815
3.	WTAGAAAYY	9	0.15259
4.	CNDPFLGVY	9	0.15232
5.	GAAAYYVGY	9	0.09963
6.	ITDAVDCAL	9	0.08501
7.	STQDLFLPF	9	0.06828

Table 3.

Helper T-Cell (HTL) epitopes for spike glycoprotein using IEDB MHC-II module.

Allele	Start	End	Length	Peptide	VaxiJen probability (threshold 0.4)
HLA-DRB1*01:01	513	527	15	LSFELLHAPATVCGP	Antigen
HLA-DRB1*01:01	512	526	15	VLSFELLHAPATVCG	Antigen
HLA-DRB1*01:01	511	525	15	VVLSFELLHAPATVC	Antigen
HLA-DRB1*01:01	514	528	15	SFELLHAPATVCGPK	Non-antigen
HLA-DRB1*01:01	510	524	15	VVVLSFELLHAPATV	Antigen
HLA-DRB1*01:01	509	523	15	RVVVLSFELLHAPAT	Antigen
HLA-DPA1*01/DPB1*04:01	504	518	15	GYQPYRVVVLSFELL	Antigen
HLA-DPA1*01/DPB1*04:01	507	521	15	PYRVVVLSFELLHAP	Antigen
HLA-DPA1*01/DPB1*04:01	506	520	15	QPYRVVVLSFELLHA	Antigen
HLA-DPA1*01/DPB1*04:01	505	519	15	YQPYRVVVLSFELLH	Antigen
HLA-DPA1*01/DPB1*04:01	508	522	15	YRVVVLSFELLHAPA	Antigen
HLA-QA1*03:01/DQB1*03:02	761	775	15	TQLNRALTGIAVEQD	Antigen
HLA-QA1*03:01/DQB1*03:02	762	776	15	QLNRALTGIAVEQDK	Antigen

Table 4.

Helper T-Cell epitopes for spike glycoprotein using Toxinpred server.

Peptide sequence	SVM score	Prediction	Hydrophobicity	Hydropathicity	Hydrophilicity	Charge	Mol.wt
LSFELLHAPATVCGP	−1.10	Non-Toxin	0.12	0.85	−0.60	−0.50	1555.04
VLSFELLHAPATVCG	−1.31	Non-Toxin	0.16	1.23	−0.70	−0.50	1557.06
VVLSFELLHAPATVC	−1.39	Non-Toxin	0.19	1.54	−0.80	−0.50	1599.14
VVVLSFELLHAPATV	−1.49	Non-Toxin	0.22	1.65	−0.83	−0.50	1595.14
FVFLVLLPLVSSQCV	−1.06	Non-Toxin	0.28	2.23	−1.23	0.00	1664.31
RVVVLSFELLHAPAT	−1.58	Non-Toxin	0.07	1.07	−0.53	0.50	1652.19
GYQPYRVVVLSFELL	−1.40	Non-Toxin	0.04	0.66	−0.70	0.00	1783.33
PYRVVVLSFELLHAP	−1.50	Non-Toxin	0.06	0.81	−0.63	0.50	1740.30
QPYRVVVLSFELLHA	−1.65	Non-Toxin	0.02	0.68	−0.61	0.50	1771.32
YQPYRVVVLSFELLH	−1.53	Non-Toxin	0.00	0.47	−0.73	0.50	1863.42
YRVVVLSFELLHAPA	−1.63	Non-Toxin	0.08	1.03	−0.66	0.50	1714.26
TQLNRALTGIAVEQD	−1.87	Non-Toxin	−0.17	−0.26	0.06	−1.00	1629.03
QLNRALTGIAVEQDK	−1.56	Non-Toxin	−0.23	−0.47	0.29	0.00	1656.10

Table 5.

Final selection of HTL epitopes after IFN-epitope analysis.

S.no	Peptide	Method	Result	Score
1	VVLSFELLHAPATVC	SVM	POSITIVE	0.47729732
2	VVVLSFELLHAPATV	SVM	POSITIVE	0.48135984
3	RVVVLSFELLHAPAT	SVM	POSITIVE	0.5092653
4	GYQPYRVVVLSFELL	SVM	POSITIVE	0.6533982
5	PYRVVVLSFELLHAP	SVM	POSITIVE	0.56872818
6	QPYRVVVLSFELLHA	SVM	POSITIVE	0.60855322
7	YRVVVLSFELLHAPA	SVM	POSITIVE	0.77626818
8	TQLNRALTGIAVEQD	SVM	POSITIVE	0.43169637
9	QLNRALTGIAVEQDK	SVM	POSITIVE	0.44573276
10	WTAGAAAYY	SVM	POSITIVE	0.57688791
11	GAAAYYVGY	SVM	POSITIVE	1.0042404

Homology modeling, tertiary structure prediction and validation

The final multi-epitope subunit vaccine model was generated using I-TASSER server. The top finest threading templates for building the protein models were selected (1rqu, 6f0k, 1dd3, 3j4a, 2ftc) and based on high c-score value -0.65 and an estimated TM score of 0.63, model protein was selected. Energy minimization and refinement of the modeled structure was carried out with SPBD viewer and loop regions were identified and refined. Final structure was checked with Ramachandran plot and showed 99% of the residues are in favorable region (Figure 1a) and ERRAT server showed 81% quality score (not shown). The results of ProSA-web obtained for vaccine construct was provided the z-score value of −0.79, as the major parts of the energy plot with N-terminal region and C-terminal region showed highly positive energy values (Figure 1b).

Figure 1.

(a) Ramachandran plot is showing 99% of the residues are in favorable region; (b) ProSA-web z-score determined by X-ray crystallography (light blue) and NMR spectroscopy (dark blue) according to the length.

Analysis of vaccine construct and insilico cloning

PSIPRED produced the secondary structural information of the vaccine construct (Figure 2a). The Prosol server provided the solubility prediction calculations and average of all residues produced the value of 0.46 and indicated good solubility of the vaccine construct (Figure 2b).

Figure 2.

(a) Secondary structure predictions of vaccine construct using PSIPRED and (b) Solubility analysis vaccine construct using ProtSol.

The physicochemical features of the vaccine were analyzed by ProtParam tool and were given the Molecular weight of 44 kda, therotical pI caluculated value of 5.16. The instability index (II) is computed to be 16.88 and protein classified as stable and the aliphatic index is 94.92%. Grand average of hydropathicity (GRAVY): 0.262. The estimated half-life is: 30 h (mammalian reticulocytes, in vitro). The overall prediction of the vaccine construct is found to be probable antigen with a score of 0.5848 generated by vaxigen 2.0 and AllerTOP 2.0 has classified the construct to be non-allergen and predicted that the nearest protein to be Scarecrow 1 in Oryza sativa (Figure 3). Optimized codons sequence length was found to be 1175 nucleotides and average GC content was 59.1%. The pET28a (+) vector was used to clone the vaccine construct DNA sequence using SnapGene software (Figure 4).

Figure 3.

AllerTOP 2.0 predicted the vaccine construct is a non allergen.

Figure 4.

Cloning of the final vaccine construct into the pET28a (+) expression vector and purple represents the gene coding vaccine DNA sequence.

Molecular docking analysis and MD simulation

The adjuvant which is a 50S ribosomal protein has the ability to stimulate TLR3, Molecular docking using the GRAMXX server (http://vakser.compbio.ku.edu/resources/gramm/grammx/) produced the best protein-protein docking model complex with binding score for molecular docking produced the best structure with global energy -35.98 and interactions and attractive vanderwaals (-32.28) was selected and DIMPLOT of LIGPLUS version 1.2 visualized the interaction between chain A- TLR-3 protein and chain B – Vaccine construct (Figure 5a–c). The RMSD values of protein-ligand complexes were recorded from 0 to 20 ns. The RMSD values steadily increased from 0 to 5 ns and reached a stable state throughout the simulation. The average RMSD values of the complex were found to be 0.27 nm (Figure 6).

Figure 5.

(a) Homology modeling of vaccine construct generated by I-TASSER; (b) Docking complex of TLR-3 (PDB ID: 1ziw) and vaccine construct; (c) Docking complex of TLR-3 and vaccine construct interactions visualized using DIMPLOT (Ligplus 1.2 version) (chain A- TLR-3 protein; chain B – Vaccine construct).

Figure 6.

The RMSD values of simulated complex of the TLR-3 and Vaccine construct.

C-ImmSim considers the successive and successful immune responses of the state of the cell and model the memory of immune cells by a mechanism that increases their half-life. The result of process is that few cells increase their half-life considerably and live longer than other cells. ImmSim server immune simulation results confirmed consistency with actual immune responses. High levels of IgM indicated the primary response. Furthermore, an increase in the B-cell population was characterized by an increase in the expression of immunoglobulins which resulted in a decrease in the concentration of the antigen. Also, there is a consistent rise in Th (helper) cell population with memory development (Figure 7a–c). It was also observed that the production of IFN-γ was stimulated after immunization (Figure 7d). The results clearly explained the T cell population was highly responsive as the memory developed and all other immune cell population shown to be consistent.

Figure 7.

C-ImmSim server prediction results of immune response after administering vaccine construct; (a) Antigen and immunoglobulins; (b) B-lymphocytes cell population; (c) CD4+ helper T cells population per state; (d) Induced levels of the cytokine and Simpson index, D.

Conclusion

The vaccine candidate against spike viral surface glycoprotein of SARS-CoV-2 was designed by in silico methods. The epitopes predicted with different web servers and adjuvant linkers were used to construct a potent antigenic, non – allergenic vaccine that could elicit strong immune response against SARS-CoV-2. Docking analysis provided the validation in the form of affinity between two molecules (TLR-3 and vaccine) and stability of complex was supported by MD simulations. The in silico immune simulation confirmed immune cell response against antigen clearance rate. The computational cloning by SnapGene confirmed the strong expression of proteins. However, the experimental validation could be essential to ensure to vaccine construct efficacy against COVID-19.