3  3. Results and Discussions

3.1  Retrieval of protein sequences
The NCBI most recently dedicated a coronavirus disease data hub containing all nucleotide and protein sequences information published from across the world. In this study, we aimed at in silico prioritization potential vaccine candidates and designing a chimeric peptide vaccine for COVID-19 based on all available protein sequences in the data hub. Several bioinformatic and immunoinformatics techniques are employed with the aim to assist experimentalists in vaccine development against the virus.

3.2  Identification of Potential Vaccine Candidates
Prioritization of potential vaccine candidates could help in minimizing time, labor cost and resources for developing and optimizing the success of getting an effective vaccine against the pathogen. In total, 193 protein entries were retrieved (S-Table 1 ) and analyzed first for sequence homology with the human host proteome. This was significant to evaluate as homology between virus protein (s) to be used in vaccine designing and the host is likely to cause strong autoimmune reactions in the host [58]. This check identified two proteins: (orf1a polyprotein (Accession id, YP_009725295) and nsp3 (Accession id, YP_009725299) as host homologous thus discarded from further evaluations. Experimental studies of the human vaccines are usually done in mice because of the many practical advantages they provide compared to vaccine research in higher animal models [22]. The appropriateness of mice as an animal model for vaccine research is defined by its ability to reproduce relevant human physiology [59]. Considering this, a homology check was devoted to the pipeline and applied to the filtered human non-similar proteins (191 in number) to ensure the selection of non-similar mice proteins. This assessment resultant into shortlisting of Orf1ab polyprotein (Accession ids, QHQ82463, QHD43415, QHR63279, QHR63259, QHQ71962, QHO62106, QHQ71972, QHR63249, QHO60603, QHO62111, QHN73794, QHR84448, QHR63269, QHN73809, QHR63289, and QHO62876) and nsp3 (Accession id, YP_009724389) as mice similar proteins. The mice non-similar proteins will aid in discouraging false positive results during in vivo experimentations and accurate interpretation of immune protective efficiency of the prioritized vaccine candidates against the virus [22,60]. Next, enumeration of transmembrane helices in the pooled 174 mice non-similar proteins was accomplished. This transmembrane topology characterization is deemed vital in relatedness to the afterward experimental expression studies of the proteins. Protein with transmembrane helices less than 2 in numbers are often considered as best vaccine candidates as multiple helices make recombinant proteins purification and expression difficult in vaccine development [20]. This result 33 proteins spanning across five types, and include: Orf3a polyprotein (Accession id, QHQ71974, QHR84450, QHD43417, QHQ71964, QHQ82465, QHO60595, YP_009724391, QHN73811, QHN73796, QHO62878), nsp4 (Accession id, YP_009725300), membrane glycoprotein (Accession id, QHD43419, QHQ71966, QHQ71976, QHQ82467, YP_009724393, QHO60597, QHN73813, QHN73798, QHO62880, QHR84452), membrane protein (Accession id, QHR63283, QHR63263, QHR63273, QHR63293, QHR63253), matrix protein (Accession id, QHO62109, QHO62114), and nonstructural protein NS3 (Accession id, QHR63251, QHR63281, QHR63261, QHR63271, QHR63291) containing multiple helices therefore not proceeded further. The creation of adhesin-based vaccines is considered an attractive and effective strategy and is being explored as a solution to number of infectious pathogens [61]. The idea behind exploiting adhesin for a vaccine is based on the promising preclinical findings. The aim is to confer protective immunity via two main mechanisms: (i) opsonization driven by opsonising antibodies that is capable of binding the target antigen as an immunological tag leading to activation of other components of the host immune system for enhance recognition of the pathogen and subsequent complement system activity and virus killing by phagocytosis, (ii) neutralization driven by adhesin-specific antibodies that block virus binding ability to host tissues. The adhesion probability computation revealed 26 protein to have adhesion probability value greater than a threshold as tabulated in S-Table 2 and can be ideal putative vaccine candidates against COVID-19. The adhesin probability of protein ranges from 0.593 to 0.796 (mean, 0.645).Antigenicity of proteins was predicted to reflect their ability of binding to products of adaptive immunity: antibodies or T-cell receptors. In total, 7 proteins: nsp8, nsp9, nsp10, 3C-like proteinase, spike glycoprotein, surface glycoprotein, and ORF1ab polyprotein were recognized as antigenic and scored higher than the threshold. Coronavirus nsp8 suggested having diverse activities, including template-dependent RNA polymerase activities, canonical RNA-dependent RNA polymerases, cofactor function of nsp8 for nsp12-mediated RNA-dependent RNA polymerase activity, and metal ion-dependent RNA 3′ polyadenylation activities [62]. Nsp9 is a non-structural protein 9, key to coronavirus replication and is a single-stranded RNA-binding protein [63]. Nsp10 is a critical cofactor that switches on multiple enzymes in replication cycle. It is known to interact with nsp14 and nsp16 subunits activating their respective 3′-5′ exoribonuclease and 2′-O-methyltransferase functions [64]. The 3C-like proteinase is main cysteine protease and is nonstructural protein number 5 (nsp5) and essential in mediating cleavage of nsp4 to nsp16 [65]. The trimeric transmembrane coronavirus spike glycoprotein initiates infectious cycle by binding to a specific receptor on the host membrane followed by viral fusion [66]. The surface glycoprotein was analyzed to be different in a sequence patch (Leu3-Phe11) at the start of the spike glycoprotein. The ORF1ab is replicase polyprotein cleaved by papain-like protease and 3C-like protease at specific cleavage sites to yield 15 to 16 non-structural proteins (nsps) [67].  The final numbers of potential vaccine candidates obtained in this step by step subtraction phase are presented in Figure 2 .
Table 1 The final set of selected B-cell derived T-cell epitopes for the potential three vaccine candidates against COVID-19.
Protein Common B and T-cell Epitopes Antigenicity (cut off score, 0.4) MHCphred (IC50 score, 100 nM) Allergenicity Virulentpred (cut off score, 0.5)
Nsp8 DRDAAMQRK 0.8641 35.81 non-allergen 1.0606
QARSEDKRA 0.5770 24.15 non-allergen 1.0606
Proteinase EDMLNPNYEDL 1.0913 21.68 non-allergen 1.0600
EFTPFDVVR 1.6049 5.07 non-allergen 1.0604
Spike Glycoprotein VNNSYECDIPI 1.0996 23.93 non-allergen 1.0593
Table 2 Top 10 refined models of the MEPVC along with the initial input structure.
Model RMSD MolProbity Clash score Poor rotamers Rama favored GALAXY energy
Initial 0.000 3.312 95.9 3.8 93.6 20995.70
MODEL 1 1.734 1.134 1.2 0.0 95.7 -3784.31
MODEL 2 1.920 1.444 2.4 0.0 93.6 -3750.72
MODEL 3 1.766 0.997 0.6 0.0 95.7 -3743.81
MODEL 4 1.667 0.997 0.6 0.0 95.7 -3739.43
MODEL 5 2.611 1.187 0.9 1.3 95.2 -3731.44
MODEL 6 1.782 1.187 0.9 1.3 95.2 -3724.93
MODEL 7 2.907 1.144 1.5 0.0 96.3 -3722.41
MODEL 8 1.999 1.459 3.0 0.6 94.7 -3718.59
MODEL 9 3.194 1.192 1.8 0.0 96.3 -3715.99
MODEL 10 2.235 1.310 2.7 0.6 96.3 -3713.14
Fig. 2 The final set of potential vaccine candidates filtered through several in silico checks presented as interactive Venn. The total number of proteins (TPs), filtered to 191 host non-homologous proteins (HNHPs), 174 mouse non-similar proteins (MNSPs), 26 adhesive proteins (APs), 7 antigenic proteins and finally to 3 potential vaccine candidates (PVCs).

3.3  B and T-cell Epitopes Mapping
Identification of epitopes in given antigens is vital for a number of practical reasons, including understanding etiology of a disease, monitoring of immune system, development of diagnostic assays, and epitope-based vaccines designing [68]. From vaccine designing point of view, host adaptive immunity is highly specific and is able to recognize and destroy the invading pathogen [69]. Additionally, adaptive immunity is able to remember the pathogens, creating long-lasting pathogen-specific protective memory enabling stronger attacks against the pathogen reencountered on successive times [70]. This arm of host immune system is driven by lymphocytes of two types: B and T-cells responsible for the humoral and cell-mediated immunity, respectively [71]. Both cells recognize pathogen molecular components called antigens. The antigens interact with specific receptors present on the surface of B and T-cells. The activation of both these cells required antigen recognition by these receptors, in addition, to the second activation signals from the innate system. The vaccine candidates prioritized in the first phase were deeply investigated for B-cell epitopes. Different lengths of linear B-cell epitopes were predicted for each protein and only recurrent epitopes simultaneously predicted by different servers were selected for chimeric vaccine designing. The B-cell epitopes predicted for the vaccine candidates were in the following order: nine for Nsp8 and 3C-like proteinase, five for Nsp9, eight for Nsp10, 34 for spike glycoprotein and surface glycoprotein, and four for ORF1ab polyprotein| partial. These B-cells epitopes are recognized as solvent-exposed antigens through B-cell receptors (BCR) and upon activation, B-cells secrete antibodies. Antibodies have different functions including neutralizing pathogens, toxins and labeling pathogens for destruction [72]. The proteins were also analyzed for T-cell epitopes through very stringent criteria of p-value less than 0.005. The epitopes were of different lengths and interact with several different alleles of MHC-I and MHC-II. For Nsp8, 94 epitopes predicted whereas Nsp9, Nsp10, 3C-like proteinase, spike glycoprotein, surface glycoprotein and ORF1ab polyprotein| partial protein were mapped for 61, 76, 153, 657, 651, and 44, respectively T-cell epitopes. These epitopes are presented on the surface through specific receptors known as T-cell receptor (TCR) allowing recognition of these antigens when displayed by antigen-presenting cells bound to MHC molecules [73]. Epitopes presented by MHC I are recognized by CD8 (cytotoxic T lymphocytes) [74] whereas those presented by MHC II are recognized by CD4 T-cells [75]. The CD4 T-cells later become helper T-cells that amplify the immune responses against the pathogen. Comparative analysis of the predicted epitopes was further carried out to select epitopes that are common to B-cell, CD4 T-cell and CD8 T-cell alleles in order to design a specific, effective, and strong vaccine. On this basis, three epitopes from Nsp8 (DRDAAMQRK, QARSEDKRA, EQAVANGDSEV), none for Nsp9 and ORF1ab polyprotein| partial protein, four for Nsp10 (GCSCDQLREP, YLASGGQPIT, YLASGGQPI, TVTPEANMDQESFG), three for 3C-like proteinase (EDMLNPNYEDL, KYNYEPLTQDHV, EFTPFDVVR), five for spike glycoprotein (RVYSTGSNVFQ, VNNSYECDIPI, LADAGFIKQYGDCLG, GQSKRVDFC, RNFYEPQIITTD) and surface glycoprotein (RVYSTGSNVFQ, VNNSYECDIPI, LADAGFIKQYGDCLG, GQSKRVDFC, RNFYEPQIITTD). The epitopes were then reevaluated in antigenicity check to make sure their binding potential of binding to immune cells. For Nsp8, DRDAAMQRK and QARSEDKRA were found antigenic with score higher than the default threshold of 0.4 whereas the third epitope EQAVANGDSEV was found non-antigenic hence removed. All the four epitopes of Nsp9 revealed non-antigenic therefore not processed further. In case of 3C-like proteinase, KYNYEPLTQDHV was found non-antigenic whereas EDMLNPNYEDL and EFTPFDVVR were antigenic therefore considered in afterward analysis. The spike glycoprotein contains epitopes VNNSYECDIPI, and RNFYEPQIITTD as antigenic whereas none of the surface glycoprotein provided epitopes were antigenic. Following, the affinity of the filtered antigenic epitopes for the most prevalent DRB*0101 allele in humans was evaluated through IC50 value and those with value < 100 nM were classified as high affinity binders. All the pooled antigenic epitopes were found to have great ability of binding to the mentioned allele. Similarly, these epitopes were evaluated in allergenicity and virulent potential check and only virulent and non-allergen were selected. Virulent check was significant in ensuring selection of epitopes mediating infectious pathways in the host. The final selected epitopes that cleared all these checks are tabulated in Table 1.

3.4  Construction of MEPVC
A MEP was constructed first comprising epitopes finalized in the previous phase. MEP based vaccines are considered an ideal approach to prevent and treat viral infections [24]. The epitopes shown in Table 1 were linked to each other through flexible AAY linkers as such it allows efficient separation required for the effective working of each epitope. Once the MEP was designed, to its N-terminus an adjuvant of Cholera toxin subunit B (CTB) was added [76]. The schematic representation of the MEPVC is shown in Figure 3 . CTB is nontoxic part of cholera toxin and is considered an accelerator in protective immunity and a break in auto-immunity. It shows high affinity for monosialotetrahexosylganglioside displayed on variety of cell types, including gut epithelial cells, antigen-presenting cells (dendritic and macrophages) and B-cells [77]. CTB is a preferred choice as an adjuvant because of its ability of self-expression in variety of organisms and can be coupled to antigens through several approaches involving chemical manipulation and genetic fusion resulting in strong immunological responses against the antigens to which it is attached. The vaccine construct was then used in a comparative 3D structure prediction to ensure confidentiality in selection of the most suitable model for the construct with minimum structural errors.
Fig. 3 Schematic representation of the designed MEPVC.

3.5  Evaluation of Physicochemical Properties
Several different physicochemical properties of the MEPVC were deduced from its sequence. The vaccine construct is 189 amino acids long with total number of 2984 atoms. The molecular weight of the construct is ideal i.e. 21.36 kD as small size construct is easy to handle and purify during experimental evaluation. The construct has instability index value of 35.43, signifying its high stability. The aliphatic index computed for the construct is 79.68, reflecting high thermostability. The estimated half-life in mammals, yeast, and Escherichia coli is 30 hours, > 20 hours, and > 10 hours, respectively. The Grand average of hydropathicity (GRAVY) score is -0.315 which highlights hydrophilic nature of the construct. The theoretical pI is 6.10 pointing to construct slightly acidic nature.

3.6  MEPVC Secondary and Tertiary Structure
The secondary structure elements of the vaccine construct can be divided into the following order: alpha helix (56.08 %), 310 helix (0 %), Pi helix (0 %), beta bridge (0 %), extended strand (16.40 %), beta-turn (7.94 %), bend region (0 %), and random coil (19.58 %). Compared to the I-tasser, phyre2, and Swiss-model, the 3Dpro predicted structure was determined as the most suitable structure based on the complete modeling of the given length of the amino acid sequence. Loop modeling was done at Leu29-Gln37, Ser51-Gln70, Glu72-Gln77, Glu57-Ser76, Val103-Thr113, and Glu167-Asp186. The model was refined to minimum RMSD of 1.734 Å and molprobity score of 1.134 that is quite low compared to the original structure score of 3.312, reflecting good quality of the modelled structure. Similarly, the clash score in contrast to the original structure is 94.7 times lower demonstrating steric clash free structure. The galaxy energy of the structure is very stable (-3784.31) and Ramachandran favored distribution increased from 93.6 to 95.7 %. The top 10 refined models of the MEPVC are tabulated in Table 2. The 3D models of the vaccine construct after loop modelling and refinement is presented in Figure 4.A. The overall Z-score of the modelled structure is -4 and the score is within the range of same size proteins in the pdb illustrating good quality as depicted in Figure 4.B. Refinement of the structure Ramachandran plot demonstrated the construct to contain 93.2 % of its residues in the most favored regions, while 5.7%, 0.0%, 1.1% residues are in additional allowed region, generously allowed region, and disallowed regions, respectively (Figure 4.C). The overall average G-factor of the construct is 0.05.
Fig. 4 A. The tertiary structure of the designed MEPVC. The red, green, purple, and yellow colors represent AAY linkers, epitopes, adjuvant, and EAAK linkers, respectively. B. The z-score plot indicating the overall good quality of the MEPVC. The z-score of the input MEPVC is shown by a black dot, validating the score within the range revealed for similar size proteins. C. Ramachandran plot for the MEPVC illustrating distribution of torsion angles (blue squares) comparative to the core (shown in red) and allowed (shown in brown) regions. Residues in the generously allowed region are shown in dark yellow whereas disallowed regions are depicted in pale yellow.

3.7  Disulfide Engineering of MEPVC
Enhancing protein stability is important in many biomedical applications and is an appealing approach to emulate nature stabilizing molecular interactions [78]. The covalent disulfide bonds provide substantial stability to target proteins and disulfide engineering had achieved considerable success in broad range of applications [33]. In total, 11 pairs of residues were selected for the purpose of disulfide engineering. These include Met1-Ser81, Lys5-Ala59, Val12-Asp28, Ser16-Asp28, Thr40-Ser47, Val73-Gln77, Ala119-Ala133, Met122-Ala133, Ala126-Asp129, Leu156-Leu163, and Asn159-Asp162. The average Chi3 and energy value for the pairs is 12.54 (max, 108.65 and min, -110.64) and 3.12 (max, 4.39 and min, 1.14). The disulfide engineered MEPVC structure is presented in Figure 5 .
Fig. 5 Disulfide engineering of the MEPVC. The yellow spheres represent mutated residues.

3.8  Codon Optimization and In Silico Cloning
In the follow up experimental studies, the maximum expression of MEPVC is highly desirable [35]. One requirement for that is the codon usage of MEPVC that must be adapted according to the expression system, for instance, here we used E. coli K12 as a MEPVC expression system. The codon adaptation index (CAI) and GC content revealed for the improved sequence are highly satisfactory with value of 0.96 and 48.85, respectively strongly indicating high MEPVC expression. The MEPVC then enclosed on both sites by 6x histidine tag to ease its purification process and inserted at appropriate sites of pET28a(+) vector as shown in Figure 6 .
Fig. 6 In silico restriction cloning of the MEPVC into pET28a(+) vector. The insert is shown in red.

3.9  MEPVC Interactions with Immune Receptors
Molecular interactions and binding conformation of the designed MEPVC with TLR3 and TLR4 innate immune receptors were deciphered via a protein-peptide docking approach. Both TLR3 and TLR4 belong to toll-like receptor family of pattern recognition receptor and function to activate intracellular signaling NF-κB pathway and production of inflammatory cytokines responsible for the development of effective innate immunity [79,80]. These receptors recognize viral associated molecular patterns and induce the production of interferon leading to activation of strong host defense responses. Also, the specific adaptive immunity takes time to establish against antigens therefore it's important to evaluate MEPVC affinity for the innate immune receptors. In case of MRPVC-TLR3 complex, the patch dock predicted 10 best solutions sorted based on the docking geometric shape complementarity score (S-Table 3). A high score implies enhanced affinity of the interacting molecules and best docked conformations of the molecules with respect to each other. Solution 3 was visualized for docked conformations and intermolecular forces responsible for such high affinity of the molecules. The selection was based on the FireDock analysis (S-Table 4) which is an efficient package for refinement and reassigning procedure of docking scores to rigid body docking solutions. The global binding energy of solution 3 is better i.e. -6.39 kJ/mol compared to the rest of predicted solutions. The contribution to the total score form attractive van der Waals (VdW) energy is -12.12 kJ/mol, repulsive (VdW) energy (4.79 kJ/mol), hydrogen bond (HB) energy (-1.63 kJ/mol), and atomic contact energy (ACE) (5.95 kJ/mol). The MEPVC, within 3 Å, was noticed to posed right in the center of the TLR3 receptor (Figure 7 ) interacting via hydrogen and hydrophobic bonds with His156, Asp180, Lys201, Glu203, Ser206, Phe227, Asp229, Asp230, Ser254, Ser256, Asp257, Ser282, Tyr283, Asp284, Asp285, Glu301, Tyr302, Phe304, Glu306, Tyr307, Arg325, Glu358, His359, Lys382, Tyr383, P408, His410, Iso411, Gly431, His432, Glu434, Pro408, Asp457, Phe459, Gln483, and Glu533. Similarly, among 10 predicted MEPVC-TLR4 complexes (S-Table 5), solution 7 (S-Table 6) was affirmed as best with total global energy of -10.39 kJ/mol whereas the rest of complexes were noticed as highly unstable with score in positive. The attractive VdW, repulsive VdW, HB and ACE contribution to the global energy is -35.75 kJ/mol, 22.04 kJ/mol, -5.51 kJ/mol, and 12.61 kJ/mol, respectively. Visual analysis of the complex revealed binding of the MEPVC at the interface of chains B and D (Figure 8 ). The MEPVC is surrounded by chain B residues: Pro23, Glu24, Ser25, Asp44, Lys47, Asp50, Asp51, Arg67, Arg87, Glu89, Pro113, Gln115, Asp137, His159, Asp160, Ser184, and Lys186 and chain D residues: Lys20, Phe64, and Asp114.
Fig. 7 Binding conformation of MEPVC with respect to TLR3 innate immune receptor. The yellow labeled region pointing to the residues involved in both hydrophobic and hydrophilic interactions with MEPVC.
Fig. 8 Binding conformation of MEPVC with respect to TLR4 innate immune receptor. The yellow labeled region pointing to the residues involved in both hydrophobic and hydrophilic interactions with MEPVC.

3.10  Computational Immune Simulation
The dynamic simulations of the human immune system in response to the designed vaccine construct were deciphered through C-immsim server [40]. The vaccine construct upon administration revealed to generate robust primary immune responses. As can be seen in Figure 9A that combine IgM and IgG antibodies has a titer scale close to 10,000/ml followed by IgM antibody (> 6000 antibody titer per ml). The combined IgG1 and IgG2 and IgG1 were seen to generate high titer scale of around 6300/ml, 3100/ml, and respectively. The IgG2 antibody response revealed to be low throughout post vaccine administration period. The dimerized soluble cytokine IFN-g produced against the antigen is > 400000 ng/ml (Figure 9B).
Fig. 9 In silico simulation of the host immune system using MEPVS as an antigen. A. Antibodies titer (A) and cytokines and interleukins (B) in response to the antigen.

3.11  Molecular Dynamics Simulation
The stability and dynamics of the designed vaccine construct ensemble docked to innate immune receptors were disclosed through 100-ns of MD production run and interpreted through the root mean square deviation (RMSD) [81], root mean square fluctuation (RMSF) 82, the radius of gyration (Rg) [83] and beta factor (β-factor) [84] assays as depicted in Figure 10 .
Fig. 10 MD trajectories based calculation of Cα RMSD (top left), RMSD (top right), Rg (bottom left), and β-factor (bottom right).
The average Cα atomic distance over 10,000 frames of TLR3-MEPVC and TLR4-MEPVC was decoded through RMSD assay. An average RMSD of 3.44 Å (maximum, 6.30 Å) and 3.36 Å (maximum, 5.22 Å) were estimated for TLR3-MEPVC and TLR4-MEPVC, respectively. The TLR3 and TLR4 receptors are observed more compact than the MEPVC and as a result, continuous movements of the vaccine construct through its length are noticed at the exposed regions though rooted stable at the docked position. This seems to be responsible for bringing small structural deviation as probed by RMSD. Visual inspection of the trajectories illustrated no major global and local secondary structure conformation changes in the receptor TLR3 and TLR4 structure. The MEPVC movements with respect to the TLR3 are shown at different snapshots as illustrated in Figure 11 . The MEPVC seems responsible for the deviations in the receptor TLR4 at positions Glu1-Ser79, Hie540-Asn614, Gln1144-Gln1209, Gln1347-Pro1371, and Lys1489-Lys1491. Continuous flexibility of the loops of MEPVC exposed region is observed and is highly flexible responsible for the movements of the MEPVC at the docked site though rooted stably (Figure 12 ).
Fig. 11 Binding mode of MEPVC at different snapshots at the TLR3 docked side.
Fig. 12 Binding mode of MEPVC at different snapshots at the TLR4 docked side.
The second statistical parameter computed for both complexes was the root mean square fluctuations (RMSF) that demonstrate average dynamical residue fluctuations over a specific length of time. The mean RMSF concluded for TLR3-MEPVC system is 2.1 Å whereas, an average RMSF of 2.03 Å for TLR4-MEPVC was demonstrated. Majority of the receptors residues are showing less variability and are satisfactory stable with mean RMSF < 3 Å. The residues range mentioned above were reported to have higher RMSF values and revealed fluctuating throughout the simulation time. These fluctuations are in all likelihood as an outcome of the moving adjuvant of the MEPVC. This is supported by the fact that the mean RMSF of the MEPVC is comparatively very high than that of receptors. The RMSF of MEPVC in TLR3-MEPVC and TLR4-MEPVC is 1.86 Å (maximum, 10.16 Å) and 7.45 Ǻ (maximum, 13.04 Ǻ). The receptors stability in both complexes is reaffirmed by radius of gyration (Rg) analysis that depicted stable plot with mean value of 33.97 Å (maximum, 35.62 Å) and 41.01 Ǻ (maximum, 41.65 Ǻ) for TLR4-MEPVC. The Rg findings are coherent with that of RMSD in interpreting systems stable behavior and compact nature of the receptors. Lastly, thermal beta factor (β-factor) indicated the same pattern of residues dynamics as shown by RMSF. An average β-factor derived from TLR3-MEPVC and TLR4-MEPVC is 121.34 Ų and 123.47 Ų, respectively.

3.12  Hydrogen bond analysis
Hydrogen bonding results when a hydrogen atom attached to a highly electronegative atom is attracted by another electronegative atom [85]. In a biological system, these hydrogen bonds are vital in determining specificity and directionality fundamental in molecular recognition [86]. The patterns of hydrogen bonds for both complexes were illustrated in each frame within 3 Å in order to probe the strength of intermolecular association across the simulation period. The maximum number of hydrogen bonds of MEPVC with TLR3 and TLR4 are 11 and 12, respectively demonstrating the high strength of interactions. The number of hydrogen bonds for both complexes are shown in Figure 13 .
Fig. 13 The number of hydrogen bonds formed during simulation between MEPVC and TLR3 and TLR4 receptors.

3.13  Analysis of Salt bridges
In a protein molecule, salt bridges are formed between charged side chains of amino acids at neutral pH. The residues mainly involved in these interactions include negative full electron charge glutamine and aspartate and opposite positive full electron charge arginine and lysine [87]. The presence of salt bridges between the interacting molecules is a clear sign of strengthening interaction stability. For TLR3-MEPVC complex high numbers of salt bridges were estimated within 3.2 Å between receptor Glu8, Glu276, Arg306, Glu333 with MEPVC Lys3, Lys8, Glu83, Arg16, respectively. In case of TLR4-MEPVC complex, receptor residues Arg646, Asp268, Arg1396 are involved in salt bridging with Asp146, Arg175 and Glu184 of the MEPVC, respectively.

3.14  Estimation of Binding Free Energies
The binding free energies of both TLR3-MEPVC and TLR4-MEPVC complexes were computed using continuum solvation MMPBSA and its complementary MMGBSA. The binding free energies of the vaccine construct for the receptors were estimated considering molecular mechanics energies as well as solvation energies. The net binding free energy for both complexes revealed to be very high illustrating the high interacting affinity of the molecules. For TLR3-MEPVC, the total free energy of binding is -41.4273 kcal/mol in MMGBSA while in MMPBSA it is -84.4908 kcal/mol (Table 3 ). In case of TLR4-MEPVC, the net binding energy is -13.9690 kcal/mol and -21.2289 kcal/mol in MMGBSA and MMPBSA, respectively. The gas phase energy in TLR3-MEPVC complex is highly dominating and contributes significantly to the overall binding energy. The total gas phase energy of TLR3-MEPVC is -530.5792 kcal/mol in MMGBSA and -530.5792 kcal/mol in MMPBSA. The electrostatic contribution to the net gas phase energy is high i.e. -393.2724 kcal/mol in both MMGBSA and MMPBSA compared to van der Waals energy (-137.3068 kcal/mol). The electrostatic energy contribution to the solvation energy is 507.5596 kcal/mol in MMGBSA and 462.9381 kcal/mol in MMPBSA and is highly non-favorable to the net solvation energy. For TLR4-MEPVC, the net MMPBSA is the comparatively high (-21.2289 kcal/mol) than the MMGBSA (-13.9690 kcal/mol) with solvation energy noticed to dominate the overall interaction energies (Table 4 ). The solvation energy in MMGBSA and MMPBSA is -245.2265 kcal/mol (electrostatic contribution, -240.2983 kcal/mol) and -252.4865 kcal/mol (electrostatic contribution, -248.2779 kcal/mol). The net gas phase energy is 231.2576 kcal/ mol in both MMGBSA and MMPBSA and major favorable contribution is from van der Waals is -35.8526 kcal/mol opposed to the non-favorable contributions is 267.1101 kcal/mol. Residue wise, most of the interacting residues of TLR3 reported in docking assay have very low binding energy reflecting the highly stable nature of these residues in interaction and the key role they are playing in holding the MEPVC at the docked region. The binding energy of the hotspot residues are as follow: His156 (-2.14 kcal/mol), Asp180 (-6.21 kcal/mol), Lys201 (-1.23 kcal/mol), Glu20 (-2.13 kcal/mol), Ser206 (-0.43 kcal/mol), Phe227 (-3.45 kcal/mol), Asp229 (-1.51 kcal/mol), Asp230 (-1.88 kcal/mol), Ser254 (-2.98 kcal/mol), Ser256 (-1.52 kcal/mol), Ser282 (-1.21 kcal/mol), Asp284 (-1.99 kcal/mol), Asp285 (-0.23 kcal/mol), Glu301 (-3.13 kcal/mol), Tyr302 (-1.48 kcal/mol), Phe304 (-3.52 kcal/mol), Glu306 (-3.43 kcal/mol), Tyr307 (-2.89 kcal/mol), Arg325 (-8.12 kcal/mol), His359 (-0.17 kcal/mol), Lys382 (-0.43 kcal/mol), Tyr383 (-3.27 kcal/mol), Pro408 (-2.22 kcal/mol), His410 (-1.16 kcal/mol), Iso411 (-1.81 kcal/mol), His432 (-1.76 kcal/mol), Glu434 (-2.77 kcal/mol), Pro408 (-3.11 kcal/mol), Asp457 (-1.15 kcal/mol), Phe459 (-.5.12 kcal/mol), Gln483 (-4.11 kcal/mol), and Glu533 (-1.23 kcal/mol). In case of TLR4-MEPVC, docking residues Pro23 (-1.54 kcal/mol), Glu24 (-2.31 kcal/mol), Ser25 (-3.33 kcal/mol), Asp44 (-4.43 kcal/mol), Lys47 (-5.23 kcal/mol), Asp50 (-2.21 kcal/mol), Asp51 (-1.99 kcal/mol), Arg67 (-5.65 kcal/mol), Arg87 (-1.21 kcal/mol), Glu89 (-2.11 kcal/mol), Pro113 (-0.12 kcal/mol), Gln115 (-1.21 kcal/mol), Asp137 (-0.23 kcal/mol), His159 (-3.22 kcal/mol), and Asp160 (-2.43 kcal/mol) from chain B and Lys20 (-2,89 kcal/mol), Phe64 (3.21 kcal/mol), and Asp114 (-2.11 kcal/mol) from chain D are hotspot residues.
Table 3 Differences of binding free energy (Complex - Receptor - Ligand) in MMGB/PBSA for the TLR3-MEPVC complex.
Method Net Energy Component Average Std. Dev. Std. Err. of Mean
MMGBSA DWAALS -137.3068 7.0875 0.7088
EEL -393.2724 43.8913 4.3891
EGB 507.5596 41.1568 4.1157
ESURF -18.4077 0.8700 0.0870
DELTA G gas -530.5792 41.8517 4.1852
DELTA G solv 489.1519 41.4134 4.1413
DELTA TOTAL -41.4273 7.7453 0.7745
MMPBSA VDWAALS -137.3068 7.0875 0.7088
EEL -393.2724 43.8913 4.3891
EPB 462.9381 39.2924 3.9292
ENPOLAR -16.8497 0.6609 0.0661
EDISPER 0.0000 0.0000 0.0000
DELTA G gas -530.5792 41.8517 4.1852
DELTA G solv 446.0884 39.5288 3.9529
DELTA TOTAL -84.4908 10.3941 1.0394
Table 4 Differences of binding free energy (Complex - Receptor - Ligand) in MMGB/PBSA for the TLR4-MEPVC complex.
Method Net Energy Component Average Std. Dev. Std. Err. of Mean
MMGBSA DWAALS -35.8526 3.0453 0.3045
EEL 267.1101 29.5270 2.9527
EGB -240.2983 27.5172 2.7517
ESURF -4.9282 0.3849 0.0385
DELTA G gas 231.2576 29.5487 2.9549
DELTA G solv -245.2265 27.4463 2.7446
DELTA TOTAL -13.9690 4.2201 0.4220
MMPBSA VDWAALS -35.8526 3.0453 0.3045
EEL 267.1101 29.5270 2.9527
EPB -248.2779 27.2761 2.7276
ENPOLAR -4.2085 0.2122 0.0212
EDISPER 0.0000 0.0000 0.0000
DELTA G gas 231.2576 29.5487 2.9549
DELTA G solv -252.4865 27.2706 2.7271
DELTA TOTAL -21.2289 5.3805 0.5380

3.15  Probable Prevention Strategy
After performing the immune system simulation to MEPVC and a thorough peptide dynamics analysis, a novel MEPVC is proposed as a probable solution to the widely spread viral coronaviruses outbreak. Meanwhile, this seems necessary to have a deep down lesson as COVID-19 viral outbreak is propagating to the human species on this part of the universe called Earth. As we are already living under a serious threat of antibiotic and antimicrobial resistance [88] which if ignored can cause havoc and the current spread of COVID-19 is a clear example of even a viral resistance coming in line. Referring to the introduction section where the probable cause of this viral attack is discussed, this is mandatory to mention that there is retaliation from animals and animal-associated viruses or bacteria observed against Homo sapiens quite frequently in the recent past [89]. This not only alarms the way of consuming animals in our food in the form of livestock but at the same time explains the advancing defense system of various species around. Either it is horizontal gene transfer among bacteria or transient nature of viral transmission throughout the world as one way or the other this is a form of advanced defense of organisms. It is meant to state that even before and after the Cognitive Revolution, Homo sapiens cannot be exempted from biological laws. While considering under the realm of these biological laws this is much expected from all other species to have the right of keeping and exercising a defense mechanism [90]. Instead of spending on wars among humans there is a dire need to be equipped for all wars to come against human species under the disguise of either environmental hazards or antimicrobial resistance. Even viruses are finding new ways of infecting hosts. Under this perspective the role of avian has been observed as much dominating in the recent past [91]. After having insights from comparative advancing mechanism we mention here and thereafter the preliminary existence of a “Theory of Retaliation” on the basis of currently available facts and observations. This theory of retaliation will play its role in coming future. The continuous attacks in the form of outbreaks by them have raised many challenges and threats to human beings. A recent mechanism study by Bazel et al reported H5N1 influenza viruses is posing threat to human and animal health and is currently unclear what restricts these interspecies jumps on the host side. It is further signified that PB1-F2 (a short viral peptide) assists H5N1 bird influenza viruses to overcome a human restriction factor of the viral polymerase complex HAX-1. It is also evidenced that a functional PB1-F2 aids in direct transmission of viruses from birds to humans [92]. Additionally, there is already a mechanism existing for the antimicrobial defense of avian eggs which illustrates their efficacy in defense mechanisms [93]. Furthermore, discovery of avian antimicrobial peptides, classified as B-defensins, present in chicken and turkey found active against bacteria, fungi, and yeast. This is another example of advance mechanism showed by avian and the possibility of a common ancestral gene between avian and other mammalian peptides seems obvious [94]. WHO alarmed quite often that influenza viruses with a vast silent reservoir in aquatic birds are impossible to eradicate while avian influenza is proved to be a threat to human health [95]. Even having the phylogenetic identity between SARS-CoV-2 and SARS-CoV, some clinical characteristics differentiate SARS-CoV-2 from SARS-CoV, MERS-CoV, and seasonal influenza infections [2,96]. This leads to the fact that viruses are attacking humans with different clinical features every time. The complex relationships between the human and animal species never faced a halt in evolving giving rise to numerous infectious pathogens. Whatsoever is the case, the dramatic impact of infectious diseases affecting the modern human population worldwide is evident and unexpected rise of coronavirus infections raised to 100000 while writing this research and getting beyond the normal control (https://www.nature.com/articles/d41586-020-00154-w). After all we belong to the kingdom Animalia with even ten on ten embarrassing similarities to chimpanzees [97]. Probable prevention in this context is to somehow avoid the excruciating utilization and consumption of other mammals as edibles from Homo sapiens and stop becoming a reason for the extinction of certain species. Avenging from these species from time to time by utilizing the microbes associated with them against human has become obvious and vaccine design and discovery is of utmost importance. This is crucial to strategize a preventive control not only for symptomatic relief but in general taking steps to prevent hunting and strategy required to maintain a balanced ecosystem where specifically avian will not dwell under threat by sapiens. The advent of various new mechanisms of viral survival has remained sapiens perplexed and a broader strategy is required to circumvent this problem.