key: cord-320490-3jmo35jc
authors: Ismail, Saba; Ahmad, Sajjad; Azam, Syed Sikander
title: Immuno-informatics Characterization SARS-CoV-2 Spike Glycoprotein for Prioritization of Epitope based Multivalent Peptide Vaccine
date: 2020-04-12
journal: bioRxiv
DOI: 10.1101/2020.04.05.026005
sha: 
doc_id: 320490
cord_uid: 3jmo35jc

The COVID-19 pandemic caused by SARS-CoV-2 is a public-health emergency of international concern and thus calling for the development of safe and effective therapeutics and prophylactics particularly a vaccine to protect against the infection. SARS-CoV-2 spike glycoprotein is an attractive candidate for vaccine, antibodies and inhibitor development because of many roles it plays in attachment, fusion and entry into the host cell. In this study, we characterized the SARS-CoV-2 spike glycoprotein by immune-informatics techniques to put forward potential B and T cell epitopes, followed by the use of epitopes in construction of a multi-epitope peptide vaccine construct (MEPVC). The MEPVC revealed robust host immune system simulation with high production of immunoglobulins, cytokines and interleukins. Stable conformation of the MEPVC with a representative innate immune TLR3 receptor was observed involving strong hydrophobic and hydrophilic chemical interactions, along with enhanced contribution from salt-bridges towards inter-molecular stability. Molecular dynamics simulation in solution aided further in interpreting strong affinity of the MEPVC for TLR3. This stability is the attribute of several vital residues from both TLR3 and MEPVC as shown by radial distribution function (RDF) and a novel analytical tool axial frequency distribution (AFD). Comprehensive binding free energies estimation was provided at the end that concluded major domination by electrostatic and minor from van der Waals. Summing all, the designed MEPVC has tremendous potential of providing protective immunity against COVID-19 and thus has the potential to be considered in experimental studies.

In December 2019, a new strain of coronavirus emerged in Wuhan city of Hubei province in China and has since spread globally. The virus belongs to clade B of family Coronaviridae in the order Nidovirales, and genera Betacoronavirus and caused pulmonary disease outbreak [1, 2] . It is positive-sense RNA, enveloped and non-segmented virus and named as SARS-CoV-2 as it share 82% sequence identity with SARS coronavirus (SARS-CoV) [3, 4] . SARS-CoV-2 caused coronavirus disease-19 and evidence suggest a zoonotic origin of this disease [5] . Though the zoonotic transmission is not completely understood but facts provide the ground that it proliferates from the seafood market Huanan in Wuhan and human-to-human transmission resultant into the exponential increase in number of cases [6, 7] . As of March 24, 386,332 cases are reported worldwide with 16,747 deaths and 102,333 recovered patients. Among the active cases, 267,252 are currently infected, 255,166 (95%) are in mild conditions and 12,086 (5%) are seriously ill. serisouly illed. Among the 119,080 closed cases, 102,333 (86%) are recovered whereas 16,747 (14%) die. On March 11, the World Health Organization (WHO) affirmed COVID-19 as a pandemic (https://www.worldometers.info/coronavirus/). SARS-CoV-2 utilizes a highly glycosylated, homotrimeric class I viral fusion spike protein to enter into host cells [8] . This protein is found in a metastable pre-fusion state which undergoes a structural rearrangement facilitating viral membrane fusion with the host cell [9] [10] [11] . The binding of S1 subunit to a host-angiotensin-converting enzyme initiates this process and disrupts the prefusion trimeric structure resulting into S1 subunit dispersion and stabilizes the S2 subunit to a post-fusion conformation [12] . The receptor-binding domain (RBD) of S1 goes through a hingelike conformational change that temporarily hides or exposes the determinants of receptor binding in order to occupy a host-cell receptor [11] . Down and up conformation states are recognized where former is related to the receptor-inaccessible state and the later one explains receptor-accessible state and considered as less stable [13] [14] [15] [16] . This critical role of the spike protein makes it an important target for antibody-mediated neutralization, and detailed study of the pre-fusion S structure would provide information at atomic-level helping in the design and development of a vaccine [17] [18] [19] [20] [21] . Current data indicates that SARS-CoV-2 spike and SARS-CoV spike both share the same functional receptor (host cell) -angiotensin-converting enzyme 2 (ACE2) [22, 23] . Interestingly, ACE2 binds to SARS-CoV-2 spike ectodomain with ∼15 nM affinity, about 10-20 folds higher than ACE2 binding to SARS-CoV spike [24] . One possible reason for SARS-CoV-2 capability of spreading infection from human-to-human is SARS-CoV-2 spike's high affinity for human ACE2 [25] . Series of cellular immune and humoral responses can be triggered by SARS-CoV-2 infections [26] . Immunoglobulin G (IgG) and IgM were noticeable after the 2 weeks of onset of infection which are specific antibodies to SARS-CoV-2. High titers of neutralizing antibodies and SARS-CoV-2 specific cytotoxic T lymphocyte responses have been identified in the patients who have improved from SARS-CoV-2. This phenomenon clearly suggest that both cellular and humoral immune reactions are vital to clearing the SARS-CoV-2 infection [26] [27] [28] [29] [30] .

The study presented, herein, is an attempt to get insights about antigenic determinants of SARS-CoV-2 spike glycoprotein and highlight all antigenic epitopes [31] of the spike that can be used specifically for the design of a multi-epitope peptide vaccine construct (MEPVC) [32] to counter COVID-19 infections. The epitopes predicted by immunoinformatics techniques were fused together as well as to β-defensin adjuvant [33, 34] to boost the antibody production and long-

The MEPVC affinity for an appropriate immune receptor as an agonist was checked in the step of molecular docking [60] . TLR3 available under PDB id of 1ZIW was retrieved and used as a receptor molecule. TLR3 also named as CD283 is a transmembrane protein belongs to the family of pattern recognition receptor [61] . It detects viral infection-associated dsRNA and evoke the activation of interferon regulatory transcription factor (IRF3) and (Nuclear Factor kappa-lightchain-enhancer) NF-kB [62] . Unlike other TLRs, TLR3 uses TIR-domain-containing adapterinducing interferon-β (TRIF) as a primary adapter [63] . IRF3 eventually induces the development of type I interferons leading to the activation of innate immune system and eventually to long lasting adaptive immunity [64] . The TLR3 receptor and MEPVC were used in a blind docking approach through an online PATCHDOCK server interface [65] . The interacting molecules were docked according to the shape complementarity principle. The clustering RMSD is allowed to default 4.0 Å. The output docked solutions were immediately refined with Fast Interaction Refinement in Molecular Docking (FireDock) server [66] which provides an efficient framework for refining PATCHDOCK complexes. The refined complexes were examined and one with lowest global energy was considered as top ranked. The opted complex was subjected to in-depth MEPVC conformation with respect to the TLR3 using UCSF Chimera 1.13.1 [67] .

MD simulation was applied on the selected top complex for 50-ns to understand complex dynamics and stability for practical applications. This assay was categorized into three phases: (i) parameters file preparation (ii) pre-processing, and (iii) simulation production [68] . In first phase using an antechamber module of AMBER16 [69] , complex libraries and set of parameters for TLR3 and MEPVC were generated. The complex system was solvated into 12 Å TIP3P solvation box achieved through Leap module of AMBER. The intermolecular and intramolecular interactions of the system were determined by ff14SB force [70] field. Counter ions in the form of Na + were added to the system for charge neutralization. In the system pre-processing stage, complex energy was optimized through several rounds: minimization of complete set of 6 hydrogen atoms for 500 steps, minimization of system solvation box energy for 1000 steps with restraint of 200 kcal/mol -Å 2 on the remaining system, minimization of complete set of system atoms again for 1000 steps with applied restraint of 5 kcal/mol -Å 2 applied on system carbon alpha atoms, and 300 steps of minimization on system non-heavy atoms with restraint of 100 kcal/mol -Å 2 on other system components. The complex system then underwent a heating step where the complex was heated gradually to 300K through NVT ensemble, maintained through Langevin dynamics [71] and SHAKE algorithm [72] to restrain hydrogen bonds. Complex equilibration was achieved for 100-ps. Pressure on the system was maintained using NPT ensemble allowing restraint on Cα atoms of 5 kcal/mol -Å 2 . In the simulation production, trajectories of 50-ns were produced on time scale of 2-fs. Non-bounded interactions were differentiated by describing cut-off distance of 8.0 Å.CPPTRAJ module [73] was lastly used for statistical computation of different structure parameters to probe complex stability. The MD simulation trajectories were visualized and analyzed in Visual Molecular Dynamics (VMD) 1.9.3 [74] .

Axial frequency distribution or simply AFD [75] is a novel analytical technique run on simulation trajectories to access ligand 3D conformation with respect to reference receptor atom. Such local structural movements are not captured by any other available technique. AFD can be mathematically presented by Eq.I,

where, i and j are ligand atom coordinates on X and Y axis with cut-off value k and l, respectively. The mi,j sums interactions frequency that fall in the coordinate (i,j).

The interaction energy and solvation free energy for TLR3 receptor, MEPVC, TLR3-MEPVC complex were calculated utilizing the MMPBSA.py module [76] of AMBER16. Also, an average of the above was estimated as a net binding free energy of the system. The binding free energy was computed through MM-PBSA method and its counterpart MM-GBSA of AMBER with objective to derive the difference between bound and unbound states of solvated conformations of the same molecule [77] . 1 . Computational approach adopted for the design of a SARS-CoV-2 spike protein based MEPVC.

The SARS-CoV-2 spike protein was targeted for MEPVC designing because of many filters it fulfilled required for a potential vaccine candidate. First, it does not share any significant homology to the human host and as such chances of autoimmune responses are negligible [78] . Second, the protein is also not found to have any sequence identity to the mouse proteome and thus accurate immunological findings can be deciphered from in vivo mice experimentations [79] . This spike protein only harbored one transmembrane helix ensuring the wet lab protein cloning and expression for antigen analysis easy [80] . Antigenicity is another factor that make this candidate highly suitable for vaccine designing as this allows efficient binding to the products of host immune system [81] . Further, this protein is strongly adhesive which makes it an excellent target for creation of adhesion based vaccine [82] . Lastly, all the sequences of SARS-CoV-2 spike protein are highly conserved thus a vaccine based on its sequence will be highly likely to have broad spectrum immunological implications [83] . Prioritization of potential epitopes for the SARS-CoV-2 spike protein commenced with the mapping of B-cell epitopes that predicted total of 34 epitopes of vary length ranging from one to 62 (S Table 1 ). The average score predicted for these B cell epitopes is 0.470, maximum (max) of 0.696 and minimum (min) of 0.188 (Fig.2) . Each Bcell epitope was then analyzed in MHC-1 alleles binding regions prediction [84] . The predicted epitopes were then screened and stringent criteria of lowest percentile score was used to choose the excellent binders. Afterward, the B cell epitopes were simultaneously run in MHC-II alleles binding [85] . Likewise MHC-I, reference set of MHC-II binding were: HLA-DRB4*01:01, HLA-DRB1*04:01, HLA-DRB1*04:05, HLA-DRB1*07:01, HLA-DRB1*09:01, HLA-DRB1*11:01, HLA-DRB1*03:01, HLA-DRB1*13:02, HLA-DRB1*15:01, HLA-DRB3*01:01, HLA-DRB1*12:01, HLA-DRB3*02:02, HLA-DRB1*08:02, HLA-DRB1*01:01, and HLA-DRB5*01:01. The MHC-II predicted epitopes were also filtered on basis of percentile score and then cross checked with the selected MHC-I allele and those common in both classes were considered only which were 50 in numbers. The shortlisted common MHC-I and MHC-II epitopes then subjected to antigenicity check. In this check, ability of the filtered B-cell derived T-cell epitopes ability to evoke and bind to products of adaptive immunity. This yielded 38 epitopes all of which have strong ability to bind to the most prevalent DRB*0101 with average IC50 score of 35.6552, max of 98 and min of 0.89.The antigenic epitopes then underwent allergenicity check to discard allergic peptides that may cause allergic reactions [86] . This resulted into 31 epitopes. Non-toxic epitopes were 7 whereas 6 were IFN-gamma producer (Fig.3) . The set of epitopes obtained at different stages of epitope mapping phase is tabulated in Table 1 . These epitopes appear to provide coverage to 98% of the world population (Fig.4) 

and green peaks are those predicted as epitopes and non-epitopes, respectively. 

Prioritized T cell epitopes derived from B cells were fused together tandemly by AAY linkers to make a multi-epitope peptide (MEP). AAY linker avoid formation of junctional epitopes and as such enhance epitope presentation [87] . To the N-terminal of MEP, EAAAK linker was added to attach β-defensin as an adjuvant leading to the design of a MEPVC. The MEPVC is schematically shown in Fig.5A . MEPVC offers many advantages compared to a separate antigenic peptide. Such vaccines induce both CD4+ and CD8+ responses and the antigens optimization are optimal. EAAAK is a rigid spacer and allow separation of the attached domain and promoting efficient immune processing of the epitopes [88] . β-defensins are potent immune adjuvants as they are capable of significantly enhancing production of lymphokines resulting into antigen-specific Ig production and T cell-dependent cellular immunity.

The sequence of MEPVC is: GIINTLCKYYCRVRGGRCCVCSCCPKEEQIGKCSTRGRKCCRRKKECAAKYAWNRK CISACYIAPGQTGKICCYPRRARSVCSACYTVYDPCQPCAAYVYDPLCPELCAYCKN HTSCDV.

Physicochemical properties of the MEPVC were evaluated in order to assist experimentalists in the field to setup experiments accordingly in vitro and in vivo. The length of MEPVC is spanned across 110 amino acids and has molecular weight of 13.30 kDa. Vaccine construct with weight less than 110 kDa is generally believed to effective vaccine target because of its easier purification. Theoretical pI of MEPVC is 9.8 and aids in locating MEPVC on 2D gel. MEPVC aliphatic index is 69.08 projecting the vaccine thermostable at different temperatures. The total number of negatively charged and positively charges residues are 8 and 22, respectively. The grand average of hydropathicity (GRAVY) score computed for the MEPVC is -0.545, illustrating hydrophilic nature of the protein and is likely to interact with water molecules. The estimated half-life of MEPVC in mammalian reticulocytes, in vitro is 30 hours, yeast, in vivo is greater than 20 hours, and Escherichia coli, in vivo higher than 10 hours. The antigenicity of the MEPVC was cross-checked and predicted highly antigenic with value of 0.69. Total entropy of the protein is 17.0170 which is considered ideal and also the vaccine has no transmembrane helices (alpha helical transmembrane protein, 0.0474783 and beta barrel transmembrane protein, 0.0060384) hence no difficulties can be anticipated in cloning and expression analysis. The predicted solubility upon overexpression of MEPVC is 0.965751 reflecting higher solubility of MEPVC.

The 3D model of the MEPVC was constructed using ab initio 3Dpro predictor as no appropriate template was available for homology modeling and threading methods. The 3D structure of MEPVC is shown in Fig.5B . The structure secured 85.4% of residues in the Ramachandran favored, 12.6%, 1.9% and 0% residues in additionally allowed, generously allowed and disallowed regions, respectively. As the predicted MEPVC unit has number of loop regions that need to be modeled proper before moving forward. In total, five sets of residues: Alas7-Lys32, Ile63-Gly69, Cys73-Arg77, Thr87-Pro102, and Asn113-Val119 were loop modeled. The loop modeled structure increased the Ramachandran favored residues percentage to 92.3%, residues in allowed region reduced to 6.7%, residues of generously allowed region to 1.0% and disallowed remained to 0%. The structure was subjected to structure perturbations and relaxations to obtain a refined model. Among the generated structures (S- Table 2 ), the first model was selected as it has improved Rama favored score, lowest stable galaxy energy of 0.96, improved clash score of 23.1 and good MolProbity value. Similarly, the structure lacks poor rotamers in contrast to the original structure. The Ramachandran statistics for the refined 3D structure are in following order: Ramachandran favored residues (93.2%), additionally allowed region (5.8%), generously allowed region (1.0%) and disallowed region (0 %). The Z-score of the refined MEPVC is -4.3 and within the score range of same size protein in structure data bases (S- Fig.1).  Fig.5 . A. Schematic depiction of the MEPVC. B. The original predicted 3D MEPVC structure and refined along with respective Ramachandran plots. AAY linkers are shown in red while epitopes are in coal and yellow is for EAAAK linker. Cyan color represents the β-defensin adjuvant. In the Ramachandran plot, the torsion angles are shown by black squares dispersed across the core secondary structures (colored as red). The allowed regions can be understand by yellow, generously allowed by pale yellow and disallowed by white region. The top right, top left, bottom right and bottom left represent quadrants for left handed alpha helices, beta sheets, right handed alpha helix, and no elements, respectively.

Further, disulfide engineering of the MEPVC was performed in order to optimize molecular interactions and confer considerable stability by attaining precise geometric conformation [89, 90] . Eight pairs of residues were selected to be replaced with cysteine amino acid. These pairs are: Gln7-Ala19 (χ3 angle,+118, energy value, 4.20 kcal/mol), Cys18-Leu21 (χ3 angle,+84.35, energy value, 3.69 kcal/mol), Lys44-Ala47 (χ3 angle,+74.17, energy value,5.59 kcal/mol), Arg57-Ala61 (χ3 angle,+122.71 , energy value,6.14 kcal/mol), Ala72-Ala85 (χ3 angle,-62.92, energy value,4.40 kcal/mol), Ala73-Ala82 (χ3 angle, , energy value, kcal/mol), Leu92-Glu95 (χ3 angle, -102.37 , energy value, 3.86 kcal/mol), and Phe111-Pro117 (χ3 angle,-96.20 , energy value,4.14 kcal/mol). These residues have either higher energy level i.e. > 2 kcal/mol and χ3 angle out of range (< −87 and > + 97) were selected on purpose to make them stable. The original and disulfide mutant MEPVC structures are shown in Fig.6 . The primary purpose of in silico cloning of the MEPVC was guide molecular biologist and genetic engineers about the possible cloning sites and predicted level of expression in a specific expression system for instance here in this study we used E. coli K12 system. Prior to cloning, reverse translation of the MEPVC sequence was conducted to have an optimized codon usage as per E. coli K12 to yield its max expression. The CAI value of the improved MEPVC sequence is 1 indicating ideal expression of the vaccine [91] . The GC content whereas is 53.2 % nearly to the E. coli K12 and range within the optimum ranged between 30 % and 70%. The cloned MEPVC is shown in Fig.7 . 

Both primary and secondary immune responses seem to play a significant contribution against the pathogen and may be compatible to the actual immune response. The in silico host immune system response to the antigen is shown in Fig.8 . High concentration of IgG +IgG and IgM was characterized at the primary response, followed by IgM, IgG1+ IgG2 and IgG1 at both primary and secondary stages with concomitant of antigen reduction. Additionally, robust response of interleukins and cytokines were observed. All this suggest the efficient immune response and clearance of the pathogen upon subsequent encounters. Elevated B cell population including memory cells and different isotypes in response to the antigen, points to the long lasting formation of memory and isotype switching. The T helper cell population additionally with the cytotoxic T cell and their respective memory development are in strong agreement of strong response to the antigen. 

Bioinformatic modelling driven molecular docking of the desingned MEPVC to one representative innate immune response receptor TLR3 was carried out in order to decode MEPVC potential of binding to the innate immune receptors. This was fundamental to understand as TLR3 is significant in recognition of virus associated molecular patterns and of activaiton of type I interferons and NF-kappa B. The docking assesment predicted top 20 complexes sorted mainlny on scoring functions along with interacting molecules area size, desolvation energy, and complexes actual rigid transformation. Following, the complexes were subjected to FireDock web server for refinement assay. This ease in discarding flexibility errors of the docking procedure and provide a deep refinement of the predictions thus limiting the chances of false positive docking calculations. According to the global energy, solution 8 was considered as a best complex with net global energy of -20.78 kJ/mol ( Table 2 ). This energy is the output of -16.88 kJ/mol attractive van der Waals (vdW), 3.81 kJ/mol repulsive vdW, 8.14 kJ/mol atomic contact energy (ACE), and -0.93 kJ/mol hydrogen bond energy. The docked conformation and chemical interacting residues of the MEPVC with TLR3 is shown in Fig.9 . Visual inspection of the complex leads to observation of deep binding of the MEPVC at the center of TLR3 and favor rigorously rigoursly hydrogen and weak van dar Waals interactions with various residues of TLR3. Within 3 Å, the MEPVC was noticed to formed interactions with His39,Val55,Asn57,Asp81,Phe84,Val103,Asn105,Gln107,His108,Thr126,Glu127,Ser132,Hi s129,Thr151,His156,Gln174,GLu175,Lys200,Lys201,Glu203,Asn229,Ser256,Asp280,Ser28 2,Tyr283,Tyr302,Phe304,Tyr307,Lys330,Tyr383,Tyr326,Asn328,His359,Asn361, and Glu363. 

The stability of MEPVC with TLR3 was further investigated through MD simulations. The trajectories of MD simulations were used in vital statistical analysis to decode backbone stability and residual flexibility. Root mean square deviation (RMSD) [92, 93] was performed first that compute average distance of backbone carbon alpha atoms of superimposed frames (Fig.10A) . The average RMSD for the system is 3.23 Å with max of 4.4 Å at 24-ns. An initial sudden change in RMSD can be seen up to 2.7-ns that may be due to adjustments adopted by the complex when exposed to dynamics forces and milieu. The second minor RMSD shift can be noticed between 22-ns to 26-ns. Afterward, the system is quite stable with not global and local conformational changes specified. Next, root mean square fluctuations (RMSF) [94] was applied on the system trajectories (Fig.10B) . RMSF is the average residual mobility of complex residues from its mean position. Mean RMSF for the MEPVC-TLR3 complex calculated is 1.60 Å with max of 8.6 Å pointed at the N-terminal of the MEPVC. Most of the interacting residues of the MEPVC with TLR3 are subject to minor fluctuations, a fact in analogy to complex high stability. The thermal residual deviation was assessed afterward by beta-factor (β-factor) [95] , the outcomes of which is strongly correlated to the RMSF and hence further affirming system stability (Fig.10C) . The average β-factor of the system analyzed is 88.64 Ų with max of 1956.23 Ų. Lastly, we evaluated the compactness of the system by means of radius of gyration (Rg) [96] analysis (Fig.10D) . High Rg and low Rg illustrate the magnitude of system compactness and system less tight packing. It further tell us the whether the system of interest in order or not. Highly compact system is an indication of system stability and vice versa. The mean Rg for our system is 55.8411 Å with max score of 74.884 reflecting higher ordered and compact nature of the system. 

Hydrogen bonds are dipole-dipole attractive forces and formed when a hydrogen atom bounded to a highly electronegative atom such as F, N, and O is attracted by another electronegative atom [97] . The strength of a hydrogen bond vary from 4 kJ to 50 kJ per mole. Hydrogen bonds are deemed vital in molecular recognition and provide rigidity in achieving stable conformation [98] . The frequency of hydrogen bonds in each frame of the MD simulation trajectories can be visualized in Fig.11A . These hydrogen bonds are extracted by mean of VMD HBonds plugin and are 104 in number as tabulated in Table 3 . The cut-off distance set is 3.0 Å and cut-off angle 20 degrees. Each residue pair may for one, two or more each of which is counted separately. The min, mean and max number of hydrogen bonds between MEPVC and TLR3 are 1, 5, and 12, respectively.

The distribution and bonding pattern of intermolecular interactions of the MEPVC residues atom(s) with respect to the TLR3 were studied through radial distribution function (RDF) (Abbasi et al., 2016; Donohue, 1954; Kouetcha et al., 2017) . RDF mainly describes distance 'r' between two entities and is represented by g(r). The factor 'r' is extracted from simulation trajectories and range from o to ∞ [75] . The hydrogen bonds predicted by VMD were utilized in RDF that shown only 8 interactions between MEPVC and TLR3 with good affinity for each other. In these interactions, TLR3 residues (atoms) are: Asp52:OD1, Gln78:HE21, Glu98:OE1, Asp124:OD2, Lys171:HZ3, Asp251:OD2, Glu323:OE2, and Glu328:OE1 are found to have strong radii distribution to their counterpart MEPVC residues (atoms): Arg705:HH11, Arg705: O, Lys679:HZ3, Arg706:HH12, Glu675:OE1, Asn633:HD22, Arg643:HH22, and Tyr638: HH, respectively. The RDF plots for the above said interactions are illustrated in Fig.11B . The interaction between Asp124-OD2 and Arg706-HH12 has a refined distribution pattern and highest density distribution among all. The max g(r) value for this interaction is 3.51 observed at distance range of 1 Å. This is followed by Glu328-OE1-Tyr638-HH with max g(r) value of 3.26 mostly interaction at distance range of 0.6 Å. The Glu323-OE2-Arg643-HH22 is also much refined having g(r) value of 1.98 and mostly interacts within distance range of 0.6 Å. The remaining interactions density distribution is not confined and vary considerably but important from MEPVC and TLR3 interaction point of view.

Salt bridges are non-covalent in nature and the outcome of interactions between two ionized states [101] . These interactions comprised two parts: an electrostatic interaction and a hydrogen bond. In salt bridges, lysine or arginine typically behave as base where glutamine or aspartate as acid and the bridge is created when carboxylic acid group allows a proton migration to guanidine and amine group in arginine. Salt bridges are the strongest among all non-covalent interactions and contribute to a major extent in biomolecular stability [102] [103] [104] . In total, 17 salt bridges were identified between TLR3 and MEPVC within the cut-off distance of 3.2 Å as can be depicted from Fig.11C . The higher numbers of salt bridges were recorded for TLR3-Glu628 and MEPVC-Lys685. The mean number of salt bridges for this interaction is 18, max, 35 and min, 3. The count for other salt bridges from TLR3 to MEPVC is in following order: Asp124-Arg706 (mean, 3, max,7 and min,3), Glu98-Lys679 (mean, 5, max,10 and min,4), Asp402-Arg641 (mean, 12, max, 18 and min, 4), Glu146-Arg706 (mean,7 max,11 and min,3), Glu146-Lys679 (mean,5 max,13 and min,2), Glu272-Lys655 Gln722-Side 0.02%

The vital hydrogen bond interactions involved between TLR3 receptor and MEPVC shortlisted by VMD were subjected to a novel AFD analysis to elucidate 3D movements of MEPVC atoms with respect to a reference TLR3 residues atom in simulation time. To this objective, interactions mentioned in the RDF were used in AFD. Preliminary investigation suggested that only three interactions: TLR3-Asp52-MEPVC-Arg705, TLR3-Glu328-MEPVC-Tyr638, and TLR3-Glu323-MEPVC-Arg643 are mainly represented frequently and found in most of the simulation frames. The TLR3-Asp52-MEPVC-Arg705 is uncovered in 4997 frames, TLR3-Glu328-MEPVC-Tyr638 in 4988, and TLR3-Glu323-MEPVC-Arg643 in 4985 making these interactions ideal for interpreting density distribution of the interactions on XYZ planes and also appropriate for gaining ideas about conformational changes of the interacting atoms with respect to each other. As the local structure movements and rotations are responsible for functional shifts, their understanding in our system is important to be unveiled. For TLR3-Asp52-MEPVC-Arg705 (Fig.12) , the density distribution is not uniform, dispersed and behave flexibility in affinity on all three axis for the receptor atom. Parallel, the strength of interaction is also observed affected due to these minor structural movements of the MEPVC residue atom. Though, the mentioned interaction depicts MEPVC is still within the vicinity of the TLR3 reference residue and enjoys this interaction flexibility with the said MEPVC residue during simulation. TLR3-Glu328-MEPVC-Tyr638 interaction (Fig.13 ) has less distribution area and has much higher intensity illustrating strong affinity of the interacting atoms for each other. It also gives an idea of the lesser movements of the atoms with respect to each other, an indication of a correct system conformation. The distribution area TLR3-Glu323-MEPVC-Arg643 is much dispersed though high intensity of the interaction can be seen in close vicinity (Fig.14) . 

The net free energy of binding (ΔTOTAL) in both GB and PB models are revealed favorable MEPVC-TLR3 complex in pure water. The net GB and PB energy for the MEPVC-TLR3 complex is -53.81 kcal/mol and -89.02 kcal/mol, respectively. To this energy, high contribution was noticed from gas phase energy (ΔG gas) compared to highly insignificant contributions from solvation energy (ΔG solv). In GB model, the ΔG gas energy for the system is -1889. 76 

The net free energy of the simulated system was subjected to per residues and pairwise residues decomposition to point residues that contribute majorly in system stabilization and lower energy. Molecular docking simulation studies demonstrated 64 residues from the TLR3 receptor that are in direct contact with the MEPVC but per residue decomposition assay illustrated that among the residues only Hie31, Phe55, Glu98, Hie100, Met102, Ile121, Thr122, Asp124, Glu146, Glu146, Glu174, Ser176, Phe198, Asn200, Ser225, Met249, Asp251, Tyr254, Tyr273, Phe275, Tyr278, Tyr297, Glu323, Hie324, and Tyr348 are hotspot as they contribute rigoursly in binding interaction with MEPVC at the docked side. The side chain of Hie100, Thr122, Asp124, Glu174, Ser176, Tyr254, and Glu323 contribute significantly in chemical interactions and have energy value in following order: -2.86602 kcal/mol, -3.71782 kcal/mol, -3.77019 kcal/mol, -3.80724 kcal/mol, -2.71475 kcal/mol, -2.40187 kcal/mol, and -3.54158 kcal/mol, respectively. To these TLR3 hotspot residues, the MEPVC interacting residues were also observed in quite lower energies illustrating high affinity for the receptor residues for chemical interactions. 

The binding free energy of the TLR3-MEPVC complex, TLR3 receptor, MEPVC and the net system energy is further decomposed into 100 frames extracted from simulation trajectories (S- Fig.3 ). This information deemed vital in predicting the simulation time where higher intermolecular affinity was observed and can guide about the most suitable docked conformation. In general the complex, receptor and construct energies are higher in PB compared to GB but for the total energy, the PB energies are quite lower for frames in contrast to GB. For the complex, the min, max and average binding energy reported are -67264.7 kcal/mol, -66901.5 kcal/mol, and -67069.5 kcal/mol, respectively in GB. The PB max frame energy is -66751.4 kcal/mol, min is -67120.2 kcal/mol and average is -66921.6 kcal/mol. The GB receptor max is -57111.5 kcal/mol whereas the min is -57381. 8 

Pair-wise energy contribution to the net energy of the system was accomplished in order understand pair residues role from both TLR3 and MEPVC in complex stability. We found that the Thr122 and Asp124 (-4.56 kcal/mol in GB and -5.45 kcal/mol in PB), Glu174 and Ser176 (-3.45 kcal/mol in GB and -3.77 kcal/mol in PB), Glu323 and Hie324 (-2.86 kcal/mol in GB and -3.99 kcal/mol in PB) of TLR3 receptor have high combine contribution to the net energy. In case of MEPVC, Asn633 and Thr634 (-3.21 kcal/mol in both GB and PB), Val642 and Arg643 (-5,87 kcal/mol in GB and -3.27 kcal/mol in PB) and Arg705 and Arg708 (-2.74 kcal/mol and 2.04 kcal/mol).

Taken together, we characterized SARS-CoV-2 spike glycoprotein for antigenic peptides and proposed a MEPVC by means of several computational immunological methods and biophysical calculations. The outcomes of this study could save time and associated cost that go into experimental epitope targets study. The MEPVC is capable of activating all components of the host immune system, have suitable structural and physicochemical properties. Also, it seems to have very stable dynamics with TLR3 innate immune receptor and thus has higher chances of presentation to the host immune system. However, additional in vivo and in vitro experimentations are needed to disclose its potential in fight against COVID-19.

No conflict of interest was reported by the authors.

Authors are highly grateful to the Pakistan-United States Science and Technology Cooperation Program (Grant No. Pak-US/2017/360) for granting the financial assistance. Fig.1 . PROSA-Z energy plot for the MEPVC. Fig.2 . Residue wise decomposition of net binding energy into TLR3 receptor and MEPVC interacting residues. Top (GB) and bottom (PB). Fig.3 . Binding energy decomposition per frame for TLR3-MEPVC complex (A), TLR3 receptor (B), MEPVC (C) and net energy (D).

S- Table 1 . B cell epitopes predicted for the SARS-CoV-2 spike glycoprotein.

S- Table 2 . Top 5 refined model of the MEPVC. The input MEPVC is also provided.

Glu323-Arg646 (mean,5 max,10 and min,3), Glu323-Lys655 (mean,13 max,25 and min,3), Glu328-Arg641 (mean,10 max,17 and min,4), Glu628-Lys699 (mean,9 max,28 and min,2), Glu675-Lys171 (mean,4 max,13 and min,2), Glu675-Lys172 (mean,10 max,16 and min,2), Glu98-Arg705 (mean,9 max,14 and min, 5), Glu98-Arg706 (mean

COVID-19, an emerging coronavirus infection: advances and prospects in designing and developing vaccines, immunotherapeutics, and therapeutics

A review of the 2019 Novel Coronavirus (COVID-19) based on current evidence

Emerging coronaviruses: genome structure, replication, and pathogenesis

The deadly coronaviruses: The 2003 SARS pandemic and the 2020 novel coronavirus epidemic in China

Zoonotic origins of human coronaviruses

Estimating clinical severity of COVID-19 from the transmission dynamics in Wuhan, China

Novel, others, The epidemiological characteristics of an outbreak of 2019 novel coronavirus diseases (COVID-19) in China, Zhonghua Liu Xing Bing Xue Za Zhi= Zhonghua Liuxingbingxue Zazhi

Structure, function, and antigenicity of the SARS-CoV-2 spike glycoprotein

Structure, function, and evolution of coronavirus spike proteins

The coronavirus spike protein is a class I virus fusion protein: structural and functional characterization of the fusion core complex

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

Tectonic conformational changes of a coronavirus spike glycoprotein promote membrane fusion

Cryo-electron microscopy structures of the SARS-CoV spike glycoprotein reveal a prerequisite conformational state for receptor binding

Shi, others, Immunogenicity and structures of a rationally designed prefusion MERS-CoV spike antigen

Lanzavecchia, others, Unexpected receptor functional mimicry elucidates activation of coronavirus fusion

Cryo-EM structures of MERS-CoV and SARS-CoV spike glycoproteins reveal the dynamic receptor binding domains

others, Characterization of novel monoclonal antibodies against MERScoronavirus spike protein

Wu, others, Identification of the immunodominant neutralizing regions in the spike glycoprotein of porcine deltacoronavirus

Development of epitope-based peptide vaccine against novel coronavirus 2019 (SARS-COV-2): Immunoinformatics approach

Therapeutic options for the 2019 novel coronavirus (2019-nCoV)

Receptor recognition by the novel coronavirus from Wuhan: an analysis based on decade-long structural studies of SARS coronavirus

Greenough, others, Angiotensin-converting enzyme 2 is a functional receptor for the SARS coronavirus

The novel coronavirus 2019 (2019-nCoV) uses the SARS-coronavirus receptor ACE2 and the cellular protease TMPRSS2 for entry into target cells

Stabilized coronavirus spikes are resistant to conformational changes induced by receptor recognition or proteolysis

others, A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-to-person transmission: a study of a family cluster

Perspectives on therapeutic neutralizing antibodies against the Novel Coronavirus SARS-CoV-2

Longitudinal profile of immunoglobulin G (IgG), IgM, and IgA antibodies against the severe acute respiratory syndrome (SARS) coronavirus nucleocapsid protein in patients with pneumonia due to the SARS coronavirus

Disappearance of antibodies to SARS-associated coronavirus after recovery

An efficient method to make human monoclonal antibodies from memory B cells: potent neutralization of SARS coronavirus

Anti--spike IgG causes severe acute lung injury by skewing macrophage responses during acute SARS-CoV infection

Peptide vaccine: progress and challenges

Multi-epitope vaccines: a promising strategy against tumors and viral infections

Avian antimicrobial peptides: the defense role of $β$-defensins

Exploring the Zika Genome to Design a Potential Multiepitope Vaccine Using an Immunoinformatics Approach

Exploring NS3/4A, NS5A and NS5B proteins to design conserved subunit multi-epitope vaccine against HCV utilizing immunoinformatics approaches

Molecular dynamics simulations of biomolecules

The MM/PBSA and MM/GBSA methods to estimate ligandbinding affinities

The immune epitope database (IEDB): 2018 update

BepiPred-2.0: improving sequence-based B-cell epitope prediction using conformational epitopes

Gapped sequence alignment using artificial neural networks: application to the MHC class I system

Peptide binding predictions for HLA DR, DP and DQ molecules

MHCPred: a server for quantitative prediction of peptide--MHC binding

Pangenome and immuno-proteomics analysis of Acinetobacter baumannii strains revealed the core peptide vaccine targets

VirulentPred: a SVM based prediction method for virulent proteins in bacterial pathogens

VaxiJen: a server for prediction of protective antigens, tumour antigens and subunit vaccines

AllerTOP-a server for in silico prediction of allergens

Peptide toxicity prediction

Designing of interferon-gamma inducing MHC class-II binders

Development of an epitope conservancy analysis tool to facilitate the design of epitope-based diagnostics and vaccines

Exploring dengue genome to construct a multi-epitope based subunit vaccine by utilizing immunoinformatics approach to battle against dengue infection

SCRATCH: a protein structure and structural feature prediction server

others, Galaxy: a platform for interactive large-scale genome analysis

GalaxyRefine: protein structure refinement driven by sidechain repacking

Disulfide by Design 2.0: a web-based tool for disulfide engineering in proteins

Codon usage: nature's roadmap to expression and folding of proteins

JCat: a novel tool to adapt codon usage of a target gene to its potential expression host

Innate immune pattern recognition: a cell biological perspective

Molecular docking, in: Mol. Model. Proteins

A family of human receptors structurally related to Drosophila Toll

Signaling to NF-$κ$B by Toll-like receptors

Medical Microbiology E-Book: A Guide to Microbial Infections: Pathogenesis, Immunity, Laboratory Diagnosis and Control. With STUDENT CONSULT Online Access

Recognition of doublestranded RNA and activation of NF-$κ$B by Toll-like receptor 3

PatchDock and SymmDock: servers for rigid and symmetric docking

FireDock: fast interaction refinement in molecular docking

UCSF Chimera-a visualization system for exploratory research and analysis

Combating tigecycline resistant Acinetobacter baumannii: A leap forward towards multi-epitope based vaccine discovery

The FF14SB force field

Langevin stabilization of molecular dynamics

A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations

CPPTRAJ: software for processing and analysis of molecular dynamics trajectory data

VMD: visual molecular dynamics

AFD: an application for bi-molecular interaction using axial frequency distribution

MMPBSA.py: An efficient program for end-state free energy calculations

Assessing the Performance of the MM_PBSA and MM_GBSA Methods. 1. The Accuracy.pdf

Identification of putative vaccine candidates against Helicobacter pylori exploiting exoproteome and secretome: a reverse vaccinology based approach

Vaxign: the first web-based vaccine design program for reverse vaccinology and applications for vaccine development

PanRV: Pangenome-reverse vaccinology approach for identifications of potential vaccine candidates in microbial pangenome

Prioritization of potential vaccine targets using comparative proteomics and designing of the chimeric multi-epitope vaccine against Pseudomonas aeruginosa

Adhesins as targets for vaccine development

Exoproteome and secretome derived broad spectrum novel drug and vaccine candidates in Vibrio cholerae targeted by Piper betel derived compounds

The MHC class I antigen presentation pathway: strategies for viral immune evasion

MHC class II proteins and disease: a structural perspective

Vaccine allergies

Predefined spacers between epitopes on a recombinant epitope-peptide impacted epitope-specific antibody response

In silico design of multimeric HN-F antigen as a highly immunogenic peptide vaccine against Newcastle disease virus

Disulphide bonds and protein stability

Protein disulfide engineering

Novel immunoinformatics approaches to design multi-epitope subunit vaccine for malaria by investigating anopheles salivary protein

Molecular dynamics simulation studies of novel $β$-lactamase inhibitor

Significance of root-mean-square deviation in comparing three-dimensional structures of globular proteins

Interaction mechanisms of a melatonergic inhibitor in the melatonin synthesis pathway

Binding mode analysis, dynamic simulation and binding free energy calculations of the MurF ligase from Acinetobacter baumannii

Radius of gyration as an indicator of protein structure compactness

The hydrogen bond in molecular recognition

Hydrogen bonds in proteins: role and strength

Ultrafast scalable parallel algorithm for the radial distribution function histogramming using MPI maps

Radial Distribution Functions of Some Structures of the Polypeptide Chain

Defining the role of salt bridges in protein stability

Do salt bridges stabilize proteins? A continuum electrostatic analysis

Protein stabilization by salt bridges: concepts, experimental approaches and clarification of some misunderstandings

Contribution of surface salt bridges to protein stability: guidelines for protein engineering