The present study identified the potential cytotoxic T-lymphocyte immune signature from the infectious sequence of the SARS-CoV-2 strain, which could be used as antiviral therapeutics to evoke the specific immune response to fight against the coronavirus. Putative CTLs have been designed with advanced immunoinformatics approaches to be used as a vaccine candidate.
3.1 Coronavirus infectious sequence retrieval and assessment
The infectious sequence of coronavirus was availed from the reports and molecular evidence of coronavirus from the protein database available on NCBI (National Center for Biotechnology Information). The updated version of the infectious sequence of spike glycoprotein of coronavirus disease-2019 (COVID-19) from Wuhan, China, was retrieved with GenBank ID QHR63290.2. The spike glycoprotein of coronavirus is reported to play a vital role in infection and progression of virus to the host. It was said to isolate from the early-stage corona infected patients. Furthermore, we have analyzed the retrieved infectious sequence through the protein-protein BLAST analysis. The outcome results showed the spike glycoprotein sequence is substantially conserved and has high similarity with surface glycoprotein of severe acute respiratory syndrome coronavirus and Bat coronavirus. These shreds of evidence signified the spike glycoprotein infectious sequence, a potential candidate for cytotoxic T-lymphocytes epitopes derivation, to evoke immune response against the coronavirus.
3.2 Putative Cytotoxic T-lymphocytes epitopes prediction
The Cytotoxic T-lymphocytes (CTL) epitopes can be used as a potential vaccine candidate to induce active immunity to fight against the coronavirus. Further, they provide the immunological memory to prepare the immune system against viral re-infection. To identify the most precise and specific CTLs, five different immunological algorithms based servers were employed. The outcome CTLs from the five servers (Rankpep, BIMAS, NetMHC 4.0, CTLPred, and NetCTL-1.2) are enlisted in table-1. Among the resulting CTLs, top-scored ten epitopes were shortlisted; FVSGNCDVV, FGGFNFSQI, TLDSKTQSL, SIIAYTMSL, MFVFLVLLP, FVFLVLLPL, FRVQPTESI, NYNYLYRLF, LTDEMIAQY, and WTAGAAAYY. Shortlisted top-scored CTLs were further investigated for various immunological parameters (Antigenicity, immunogenicity, C-terminal proteasomal affinity, and TAP affinity) essential for an efficacious vaccine.
Table 1: Putative Cytotoxic T-lymphocytes epitopes for the infectious sequence of coronavirus identified by employing the five different algorithms. Top ten ranked predicted epitopes from each server are depicted here.
Epitope Rank
|
Rankpep
|
BIMAS
|
NetMHC 4.0
|
CTLpred
|
NetCTL-1.2
|
1
|
FVSGNCDVV
|
TLDSKTQSL
|
MFVFLVLLP
|
FRVQPTESI
|
LTDEMIAQY
|
2
|
FGGFNFSQI
|
SIIAYTMSL
|
FVFLVLLPL
|
NYNYLYRLF
|
WTAGAAAYY
|
3
|
FQFCNDPFL
|
GLTVLPPLL
|
VFLVLLPLV
|
AQVKQIYKT
|
TSNQVAVLY
|
4
|
KVGGNYNYL
|
VLNDILSRL
|
FLVLLPLVS
|
YTNSFTRGV
|
CVADYSVLY
|
5
|
CVNFNFNGL
|
SALEPLVDL
|
LVLLPLVSS
|
TRFQTLLAL
|
KTSVDCTMY
|
6
|
FCGKGYHLM
|
ALNTLVKQL
|
VLLPLVSSQ
|
NENGTITDA
|
STECSNLLL
|
7
|
VVNQNAQAL
|
MFVFLVLLP
|
LLPLVSSQC
|
YQTSNFRVQ
|
GAEHVNNSY
|
8
|
VLYENQKLI
|
GLTVLPPLL
|
LPLVSSQCV
|
KCYGVSPTK
|
NIDGYFKIY
|
9
|
YHKNNKSWM
|
WTFGAGAAL
|
PLVSSQCVN
|
GKIADYNYK
|
YSSANNCTF
|
10
|
VSGNCDVVI
|
SLSSTASAL
|
LVSSQCVNL
|
KSTNLVKNK
|
WMESEFRVY
|
3.3 Optimization of CTL epitopes with immunological parameters
Top scored CTL epitopes were investigated for MHC binding, Transporter associated with antigen processing (TAP) affinity, C-terminal cleavage efficiency, antigenicity, and non-toxic nature characterization. The MHC binding affinity, TAP transport capability, and peptide C-terminal cleavage efficiency were analyzed using the netCTLpan server. The obtained results depicted the higher cleavage efficiency for antigen processing with a weighted score for proteasomal C-terminal affinity with a cut-off value of 0.225; the epitopes with higher scores are susceptible to proteasomal C-terminal sympathy. It will help in antigen processing into the short peptides (8-10 amino acids) for antigen presentation to Major Histocompatibility Complex -1 (MHC-1) T-cells to generate an immune response. The CTL epitopes were also evaluated for TAP affinity, which demonstrated the high TAP affinity of epitopes with a cut-off value of 0.025, shown in the table-2.
Moreover, epitopes were further analyzed for MHC binder affinity and antigenicity profiles through the VaxiJen server. Obtained results showed the high antigenic characteristics of the CTL epitope and MHC binding affinity with a threshold value of 0.4. After that nature of epitopes was examined which, depicted the non-toxic nature of all shortlisted epitopes using the Toxinpred server enlisted in table-2. These significant results demonstrated the high ranked epitopes (FGGFNFSQI, TLDSKTQSL, FVFLVLLPL, FRVQPTESI, and WTAGAAAYY) with high immunological parameters. The high-rank epitopes were further examined for their specific binding with HLA-A*0201 receptor and virus infection host membrane specific toll-like receptor-2 (TLR2).
Table 2: Optimization of CTL epitopes with immunological parameters for MHC binding, Transporter associated with antigen processing (TAP) affinity, C-terminal cleavage efficiency, combinatorial rank, antigenicity, and non-toxic nature characterization.
|
Peptide
|
MHC
|
TAP
|
Cle
|
Comb
|
%Rank
|
Antigenicity
|
Nature
|
1.
|
FVSGNCDVV
|
0.49600
|
0.17200
|
0.73303
|
0.66523
|
3.00
|
Non-antigenic
|
Non-Toxin
|
2.
|
FGGFNFSQI
|
0.16600
|
0.25400
|
0.65415
|
0.30683
|
16.00
|
Highly-antigenic (1.27)
|
Non-Toxin
|
3.
|
TLDSKTQSL
|
0.48400
|
0.58000
|
0.97582
|
0.71806
|
2.00
|
Highly-antigenic (1.06)
|
Non-Toxin
|
4.
|
SIIAYTMSL
|
0.77400
|
1.27400
|
0.97302
|
1.02478
|
0.15
|
0.5234
|
Non-Toxin
|
5.
|
MFVFLVLLP
|
0.05000
|
0.39800
|
0.13377
|
0.09005
|
50.00
|
0.4166
|
Non-Toxin
|
6.
|
FVFLVLLPL
|
0.68300
|
1.04400
|
0.97043
|
0.92745
|
0.80
|
Highly-antigenic (0.86)
|
Non-Toxin
|
7.
|
FRVQPTESI
|
0.04400
|
0.75000
|
0.93979
|
0.27420
|
32.00
|
Highly-antigenic (0.93)
|
Non-Toxin
|
8.
|
NYNYLYRLF
|
0.02500
|
2.73800
|
0.43657
|
0.19168
|
50.00
|
Non-antigenic
|
Non-Toxin
|
9.
|
LTDEMIAQY
|
0.08000
|
2.65200
|
0.97316
|
0.36526
|
9.00
|
Non-antigenic
|
Non-Toxin
|
10.
|
WTAGAAAYY
|
0.09900
|
2.83000
|
0.94703
|
0.38283
|
9.00
|
Highly-antigenic (0.630)
|
Non-Toxin
|
3.4 Binding analysis of selected epitopes with HLA-A*0201
Selected CTL epitopes were investigated for binding with the HLA-A*0201 receptor to assess their specificity for humans. HLA-A*0201 is a major type of human major histocompatibility complex class-1 cell surface receptor. The three-dimensional structure of HLA-A*0201 was availed from the RCSB protein data bank with PDB ID 1I4F. Retrieved HLA-A*0201 structure available in co-crystal form was processed and prepared molecular docking analyses. HLA-A*0201 consisted of two chains; Chain-A of size 275 amino acids and Beta-2 microglobulin chain of size 100 amino acids. The processed HLA-A*0201 structure was analyzed for stereochemical parameters through the Ramachandran plot. The Ramachandran plot illustrated that 93.7 % of residues of HLA-A*0201 lie in the favorable region, and 6.3% residues lie in the allowed region with no residues in the outlier region of the plot (Figure 1).
Moreover, the Verfiy3D server also showed the high quality of structure with more than 97% residues having average 3D-1D scores >=0.2. The optimized HLA-A*0201 structure was employed for binding analysis of selected epitopes using the HPEPDock server. The 3D structures of selected epitopes were modulated with the PepFold server. Designed 3D structures of CTL epitopes were energy minimized and optimized for molecular docking analysis with HLA-A*0201. HPEPDock server worked on the bases of the hierarchical algorithm for molecular docking with more than 72% success rate. The Top-scored high docking energy minimized models were analyzed for their molecular binding with the HLA-A*0201 receptor (Table 3). The obtained results depicted that CTL-1 and CTL-5 epitopes have the highest binding energy score of -251.779 kJ/mol and -246.834 kJ/mol, respectively. Among them, most proximal binding epitopes (CTL-1 and CTL-5) with the highest docking score were selected and further validated with virus infection specific host membrane receptor TLR-2 in humans.
Table 3: Molecular binding analyses of shortlisted CTL epitopes with HLA-A*0201 receptor.
S. No.
|
Docked CTL Models
|
CTL
Epitopes
|
Molecular Docking Scores with HLA-A*0201 receptor
|
1
|
CTL-1
|
FGGFNFSQI
|
-251.779 kJ/mol
|
2
|
CTL -2
|
TLDSKTQSL
|
-172.727 kJ/mol
|
3
|
CTL -3
|
FVFLVLLPL
|
-240.909 kJ/mol
|
4
|
CTL -4
|
FRVQPTESI
|
-205.960 kJ/mol
|
5
|
CTL -5
|
WTAGAAAYY
|
-246.834 kJ/mol
|
3.5 Validation of lead epitopes with coronavirus infection specific human membrane TLR2 receptor
Toll-like receptors play a main regulatory role in viral proteins reception and further derivation of active immunity to diminish viruses. TLR2 interaction with the ligand evokes the critical regulators of the immune system and stimulates the interferons and interleukin to fight the infections. The lead epitopes were examined for their binding with the human TLR2 receptor. The 3D crystal stricter of TLR2 was retrieved from the protein data bank with PDB ID 1E3G and energy minimized and prepared for molecular docking assays. TLR2 structure was also assessed for the stereochemical properties. Ramachandran plot analysis depicted 68.9% residue of TLR2 structure lie in the favorable region; 31.1 % residues lie in the allowed region and none residues in the outlier region. After that, lead epitopes CTL-1 and CTL5 were docked with the TLR2 receptor. Notably to strengthen our results and to define highly specific binder epitope, we employed three different molecular docking programs; HPEPDock, Hex 8.0, and Cluspro suite. HPEPDock worked based on the hierarchy algorithm, Hex 8.0 by fast Fourier transformation, and Cluspro worked based on piper rigid body methods. Outcome results depicted CTL-5 (WTAGAAAYY) interacted proximally to the binding groove of the TLR2 receptor from all three molecular docking methods.
Moreover, we analyzed and compared the top-scored CTL-1 and CTL-5 though Ligplot and Chimera molecular modeling suite. Obtained results showed the strong interaction of CLT-5 than CTL-1 with the involvement of large hydrogen bonds and hydrophobic interactions. CTL-5 formed the three hydrogen bonds with TLR2 receptor; first Thr2 with Leu392 residue of TLR2 of bond length 2.55Å, Secoddn with Thr2 with Asp419 residue of TLR2 of bond length 2.90Å and Tyr9 with Asp463 residue of TLR2 of bond length 2.48Å. Also, CTL-5 showed the involvement of a large number of Hydrophobic interactions (Figure 2). Besides, the CTL-1 epitope only showed the formation of two hydrogen bonds and lesser hydrophobic interactions shown in the table-4. These outcomes significantly elucidated the potential CTL binder with host-specific membrane receptors.
Table 4: Molecular binding analyses of shortlisted lead CTL epitopes with Human TLR2 receptor with different docking programs.
TLR2 Receptor
|
CTL Epitope
|
HPEPDock Score
|
HEX 8.0 Score
|
Cluspro Score
|
Hydrogen Bonds
|
Hydrophobic interaction
|
CTL Model-1
|
FGGFNFSQI
|
-187.772 kJ/mol
|
-311.25 kJ/mol
|
-848.6 kcal/mol
|
Gly2-Thr416 of bond length 2.32Å
|
Arg340, Gln363, Tyr364, Ile393, Ile461, Lys439, Tyr440 and Gln481
|
Gly3-Gln390 of bond length 2.22Å
|
CTL Model-5
|
WTAGAAAYY
|
-210.148 kJ/mol
|
-360.37 kJ/mol
|
-911.5 kcal/mol
|
Thr2-Leu392 of bond length 2.55Å
|
Tyr364, Thr391, Leu392, Ile393, Arg395, Asn417, Ile418, Tyr440, Ile 461, Glu481 and Tyr483
|
Thr2-Asp419 of bond length 2.90 Å
|
Tyr9-Asp463 of bond length 2.48Å
|
3.6 Molecular dynamics simulation of TLR2 and TLR2-epitope complex
The MD simulations for TLR2 and TLR2-epitope of 20 ns were assessed through trajectory analysis. For the course of 20 ns MD simulation, the stable trajectory was observed and the representative structures were obtained. The deviation of the backbone atoms for simulated structures relative to the starting structures used as a reference was evaluated through RMSD (Figure 3A). Steep magnitude RMSD variation during the entire simulation can be an implication of a malleable and free instinctive protein or the alteration of the force field. Depending upon the outcomes of the RMSD evaluation, Figure 3A represents that the RMSD fluctuation leads to stabilize at about 8 ns and 5 ns MD simulations for TLR2 and TLR2-epitope, respectively and the simulation time was acceptable. In the time from 8-20 ns and 5-20 ns, RMSD for TLR2 and TLR2-epitope have approximate values about 1.1-1.3 nm and 0.8-0.9 nm, respectively.
The integral extent of this criterion is achieved by assessing the variations arising from shifts of each of the protein residues that majorly feature the most flexible chain frames. Hence, we validated the residual fluctuations by calculating the mean fluctuation for stable trajectories of each simulation. The RMSF evaluation of all protein residues were achieved in order to check the residues that may have tend to an enhancement in the RMSD results (Figure 3A). Compelling fluctuations existed in terminal residues and in various loops of the β linker segments of the TLR2 (residues 285-348) to about 0.45 nm and TLR2-epitope (residues 105-156 and 239-254) to about 0.4 nm. Also, we noticed that residues 390 to 480 showed comparatively reduced discrepancy in RMSF values, which were the corresponding binding segment of the TLR2-epitope complex after MD simulation. Interestingly, the binding region for TLR2-epitope before MD simulation was also within the residues 390 to 480, which remains the same after simulation, suggesting as a strong binding interaction of the epitope towards TLR2 and provides stability to the TLR2-epitope (Figure 3B).
The radius of gyration analysis was performed to determine the change in compactness of TLR2 and TLR2-epitope throughout the MD simulations. The Rg plots for TLR2 and TLR2-epitope show slight fluctuations at the initial frame and attains stability after 8 ns and 5 ns with Rg score of 3.4 nm and 3.3 nm, respectively which remained approximately similar till 20 ns run (Figure 3C). When compared with TLR2, Rg value for TLR2-epitope is reduced, suggesting more compactness due to the strong binding interaction of the epitope. Similar observations were determined through SASA analysis representing the solvent defined protein surface and its orientation through folding, making the alterations in the exposed and buried regions of the surface area of proteins. SASA values for both the simulation systems were about 280 nm/S2/N and 275 nm/S2/N after 8 ns and 5 ns, respectively (Figure 3D). Also, it remains similar till 20 ns in the case of TLR2-epitope, whereas it reduces slightly about to 272 nm/S2/N after 15 ns in the case of TLR2. Here, TLR2-epitope solvation profile shows a convincing SASA value suggesting a stable structure and strong binding interaction with the epitope.
Further, cluster analysis having a RMSD based cut-off value of 0.25 nm demonstrated the development of five and ten distinctive clusters for TLR2 and TLR2-epitope simulation systems. The most dominant cluster attained after 20 ns of MD simulation for TLR2 and TLR2-epitope are shown in Figure 4. Also, secondary structure analyses of the stable trajectory for both the simulation systems were performed using the DSSP tool of GROMACS. Both, TLR2 and TLR2-epitope were formed mainly of conserved β-sheet along with various loops of the β linker segments and small coil regions as secondary structure elements infused with various small segments of bend, α-helix, turn, and β-bridge (Figure 4). Both cluster analysis and secondary structure analysis reveals the conformational changes before and after simulations for TLR2 and TLR2-epitope structures. A noticeable observation through these analyses supports the rationale that protein turns into somewhat unfolded and implements further surface to bind with the epitope.
Inaddition, the hydrogen bond landscape was assessed, which revealed the dynamic equilibration of the complex trajectories with a high number of hydrogen bonds, as shown in Figure-5. The consistent high numbers of hydrogen bonds were observed which contributed significantly to the proximal binding of the CTL epitope with the TLR2 receptor. Further, these results were strengthening by the vital contribution by the complex binding energies throughout the simulation run. These calculations with consistent high binding energies and large hydrogen bonds involvement demonstrated the stable binding of epitope with TLR2.
3.7 In Silico immunization analysis of vaccine construct
The vaccine construct was analyzed for immune response capability through the in silico immune simulation approach for 100 simulation steps. The immune simulation method was computed through machine learning with a mesoscopic scale simulator to assess the generation of key regulators of the immune system. Obtained results depicted the substantial elicitation of cytotoxic T-cell (Tc) population in response to vaccine construct. Major regulator immune systems (T- cells) were found to evoke above 1100 cells per mm3 and maintained for the long span to 30 days (Figure 6). Moreover, through the mesoscopic analyses, it was found that T- cytotoxic cells count was perpetuated in resting state (elevated level) till 30 days span and contributed to the memory cells on the immune system, which also prepare the immune system to fight off with the re-infection of specific antigenic determinants. These outcomes demonstrated the potency of immune CTL epitope to be used as a vaccine candidate against COVID-19.
3.8 3D Structural conformational and physicochemical analysis of lead CTL
The 3D structure of lead CTL was designed by the Pepfold server and examined through the Ramachandran plot. It showed the existence of structural coordinates in the favorable region of the Ramachandran plot (Figure 7). The local structure profiling depicted the stable composition of secondary structural coordinates through the probabilities of each residue at each position of the sequence. Moreover, the CTL-5 was found non-allergenic by the Allertop server. The physicochemical properties were investigated by using the Protparam server. The obtained results showed that the lead CTL-5 of size nine amino acids has a molecular weight of 905.03 Daltons and theoretical pH 5.52. The molar extinction coefficient was computed to 28480 M-1 Cm-1 with an absorption matrix of 0.1% (1g/l) 8.715 through the cysteine coupling at 280 nm in solvent condition. Lead CTL-5 was found to be globular nature with an aliphatic index score of 44.44 and classified to stable with an instability index score of -3.53, the output instability index score less than 40 reported to steady and stable coordinates. Moreover, it was estimated using the computational assays, which showed it has a long span (>2.8 hrs) half-life in mammalian reticulocytes by in vitro data statistics, >2 minutes for Escherichia coli and Yeast by in vivo data statics. The GRAVY score (grand average of hydropathy) was computed to be 0.289 and indicated the hydrophilic nature with optimal cellular localization.