Designing AbhiSCoVac - A single potential vaccine for all 'corona culprits': Immunoinformatics and immune simulation approaches.

The coronaviridae family has generated highly virulent viruses, including the ones responsible for three major pandemics in last two decades with SARS in 2002, MERS outbreak in 2012 and the current nCOVID19 crisis that has turned the world breadthless. Future outbreaks are also a plausible threat to mankind. As computational biologists, we are committed to address the need for a universal vaccine that can deter all these pathogenic viruses in a single shot. Notably, the spike proteins present in all these viruses function as credible PAMPs that are majorly sensed by human TLR4 receptors. Our study aims to recognize the amino acid sequence(s) of the viral spike proteins that are precisely responsible for interaction with human TLR4 and to screen the immunogenic epitopes present in them to develop a multi-epitope multi-target chimeric vaccine against the coronaviruses. Molecular design of the constructed vaccine peptide is qualified in silico; additionally, molecular docking and molecular dynamics simulation studies collectively reveal strong and stable interactions of the vaccine construct with TLRs and MHC receptors. In silico cloning is performed for proficient expression in bacterial systems. In silico immune simulation of the vaccine indicates highly immunogenic nature of the vaccine construct without any allergic response. The present biocomputational study hereby innovates a vaccine candidate - AbhiSCoVac hypothesized as a potent remedy to combat all the virulent forms of coronaviruses.

Introduction

The COVID-19 pandemic that started in Wuhan, China has now conquered the world leaving 194 million infected cases and 4.16 million fatalities [1], generating an even more substantial impact on the world economy. The causative virus, the SARS-CoV-2 is a 29,909 bp positive-sense single-stranded RNA virus belonging to the Coronaviridae family. However, SARS-CoV-2 is not the first outbreak of coronavirus in history. In 2002, SARS-CoV-1 threatened the world with severe acute respiratory syndrome (SARS) outbreak with 10% fatality rate [2] followed by the Middle East respiratory syndrome (MERS) outbreaks in 2012 caused by the MERS-CoV killing 34% of its infected populace [3]. Moreover, human coronavirus hCoV-HKU1 is associated with respiratory disorders ranging from common cold to pneumonia and bronchitis. Coronaviruses of hCoV-229E and hCoV-OC43 are responsible for 10–15% of the common cold cases worldwide in a seasonal pattern [4]. hCoV-NL63, however, is the least virulent species of human coronaviruses and is impaired to elicit any pathogenesis to healthy adults [5], thus, evading research scope. Nonetheless, as witnessed with spawning of more than 5700 geographical variants of SARS-CoV-2 [6], other coronaviruses are also known to infect numerous animals ranging from minks and pigs to bats and pangolins. Most of these viruses have been found to possesses remarkable ability to mutate, recombine and replicate in the organisms leading to the generation of an outstanding degree of viral diversity resulting in a vast number of new strains and even newer species in the near future. Such increasing number of coronaviruses ensemble a triple-front threat right in front of us, just waiting for the clock to tick.

The presence of spike ‘S’ proteins on the surface is a unique characteristic feature of the viruses of Coronaviridae [7]. Spike proteins are highly glycosylated type I transmembrane fusion protein composed of 1,160 to 1,400 amino acid residues. Recent studies establish human Toll-like receptor 4 (TLR4) as the native immune sensor for recognizing these spike proteins [8], [9], [10], [11]. Upon binding with the spike protein, TLR4 activates various signaling cascades via the intracellular adaptor proteins like MyD88, TRAF6, TIRAP, etc. All these signaling cascades finally lead to the activation of the NF-κB transcription factor that initiates the expression of the pro-inflammatory cytokines such as IL-1β, TNF- α, and IL-6. These cytokines in turn recruit other white blood cells and generate an abrupt as well as cumulative surge of cytokines secretion namely the “cytokine storm”. Essentially, this uncontrolled and incessant hyper inflammation combined with the direct virus-induced damage features a multi-organ level clinical presentation including an acute respiratory distress syndrome (ARDS), thromboembolism, sepsis, acute cardiac injury, cardiovascular hypertension, nephric cellular damage, hematuria, hypoxemic injury to brain tissue and ultimately death [12], [13].

A clear understanding of the role of the innate immune receptors including TLRs in augmenting protection against several viral infections has provided a new dimension to the researchers. Therefore, TLRs are now considered a major target for developing vaccines against viral infections of human [14], [15], [16]. In this context, reverse vaccinology-based approaches following experimental validation have been the popular choices behind designing most of the vaccines against SARS-CoV-2 [17], [18], [19]. Interestingly, the major immunocompetent antigen of SARS-CoV-2, i.e., spike glycoprotein has been established as a ligand for TLR4 [8], [9], [10], [11], [20]. Moreover, the interaction between spike protein and TLR4 is common across all the coronaviruses that infect human [21]. Thus, TLR4-spike glycoprotein interaction can be exploited and coupled with reverse vaccinology technique to design a multi-epitope chimeric vaccine that could efficiently generate prophylaxis against all the six coronaviruses depicted in the earlier section. Interestingly, TLR signaling pathway has long been considered as a potent therapeutic target to develop advanced vaccination strategies owing to the ability to generate a robust immune response possessing a binary nature of cytokine-mediate innate as well as immunoglobin-mediated adaptive reflexes [22], [23], [24], [25], [26].

In recent past, a number of strategies have been adopted to design and/or develop vaccine candidates against SARS-CoV-2 and a few of them have been approved for clinical application. For example, mRNA-coupled lipid nanoparticle (LNP) delivery systems (Pfizer and Moderna), DNA embedded within non-replicating recombinant adenoviral vector (AstraZeneca, Johnson and Johnson and Sputnik V) are in use for combating pandemic currently [27]. Intriguingly, spike glycoprotein has been the antigen of choice for developing vaccine against SARS-coV-2 till date. Experimental studies have suggested induction of humoral- and cell- mediated immunity against SARS-CoV-2 in murine and human model after immunization with spike protein as well as vaccine designed by triggering the spike glycoprotein [27], [28], [29], [30]. However, the efficacy of the aforementioned vaccines varies with the immunogenetic makeup of the recipient. Moreover, emergence new SARS-CoV-2 strains like delta and delta-plus having the ability to escape the neutralizing effects of vaccines need special attention. In this context, a new vaccine targeting the major route of the pathogenesis i.e., spike protein-TLR4 interaction could provide an efficacious mean to combat not only SARS-CoV-2 but also other members of the coronaviridae family that infect human. The present study hopes to engineer a vaccine peptide having potency for the same (Fig. 1, Fig. 2, Fig. 3, Fig. 4, Fig. 5, Fig. 6 ).

Fig. 1

Process flow in designing of AbhiSCoVac.

Fig. 2

A. showing the modelled tertiary structure of the vaccine beneath its van der waals surface, B. shows the antigenicity plot of the vaccine construct derived using Kolaskar and Tongaonkar method. C. shows the secondary structural map. D. is the Ramachandran plot computed using Structural Assessment tool from Swiss ExPasy and E. denotes the Z-score of the model.

Fig. 3

The vaccine peptide (green) is bound to A. TLR2 (blue) and B. TLR4-MD2 complex (purple), being shown in zoom out configuration as well as zoomed in interaction pocket.

Fig. 4

The vaccine peptide (green) is bound to A. MHC I (blue) and B. MHC II receptor (brown), being shown in zoom out configuration as well as zoomed in interaction pocket.

Fig. 5

Different stabilization parameters extracted from MD simulation studies of the complexes between the vaccine candidate and sensor proteins: A. (I) RMSD (in nm) and (II) Rg (in nm), B. (I) SASA (in nm2) and (II) number of hydrogen bonds, and C. (I) RMSF (in nm) and (II) the contributions of the van der Waals, electrostatic, polar solvation and SASA to the total binding energy (in kJ/mol). The MHC-I-vaccine, MHC-II-vaccine, TLR2-vaccine and TLR4-vaccine complexes are shown in black, red, blue and green colors, respectively.

Fig. 6

Cloned map of the vaccine cDNA sequence cloned into a pET-28a(+) plasmid for further expression in bacterial systems.

Materials and methods

Sequence retrieval and phylogenetic analysis

Phylogenetic relationship amongst various pathogenic coronaviruses was studied using their respective amino acid sequences through multiple sequence alignment and phylogenetic tree using MEGA 11 (see supplementary material for detail).

Selection of the antigenic peptides for developing the vaccine

TLR4 serves as the principal sensor in recognizing the spike glycoproteins of the virions that are considered as the principal pathogen-associated molecular patterns (PAMPs) present in coronaviruses [8], [9], [10], [11], [20], [21]. To identify the TLR4-interacting PAMPs, native TLR4-MD-2 complex (PDB ID: 3FXI) was docked with the spike proteins of different pathogenic coronaviruses (PDB IDs: 6U7H, 5I08, 6Q05, 6NZK, 5XLR, and 6VYB) employing ClusPro 2.0 () and binding affinity (ΔG) was determined by PRODIGY () [31]. Interacting peptides/amino acid sequences present in the TLR4-spike interface were extracted and analyzed using PyMol. The refined antigens were used for further studies (see supplementary material for more).

The presence of B-cell and T-cell epitopes are the critical components of an ideal vaccine candidate. Herein, the presence of B-cell cell epitopes in the refined spike protein antigen was determined using immune epitope database (IEDB) () and BepiPred 2.0 () while T-cell epitopes including both MHC-I and MHC-II epitopes were mapped using IEDB T-cell epitope prediction tool () and ProPred server (). Furthermore, antigenicity of the predicted epitopes was evaluated using VaxiJen 2.0 (). Experimental detail on the same is given in the supplementary material.

Designing and validation of the vaccine candidate

The vaccine is chimeric in nature and aims to function against six different corona viruses, therefore, all the high-ranking epitopes so screened from the spike sequences of the viruses were connected to each other using appropriate utilization of linker/spacer sequences. The MHC I epitopes were linked by glycine-proline rich GPGPG linkers while the MHC II epitopes were connected by AAY linkers. Additionally, Cholera Toxin B (CTB) protein was linked by EAAAK linker to the N-terminal of the vaccine. CTB operates as an adjuvant and boosts the immunogenicity of the vaccine construct [32]. The secondary protein structure was obtained from PolyView () and finally the tertiary structure model was constructed from the Phyre2 server (). It uses an intuitive and advanced remote homology detection method for building coordinate models by predicting ligand binding sites and analyzing the effects of the amino acid variants. Further, the validation of the construct 3D model was done by different methods. Primarily, the vaccine protein was subjected to the ERRAT tool offered by the SAVES 6.0 server (), followed by assessment with Pro-SA web server (https://prosa.services.came.sbg.ac.at) that evaluates the modelling quality based upon the Z-score it computes. Additionally, the Ramachandran plots were graphically plotted by using the SWISS-MODEL Structure Assessment tool from the ExPASy server ().

Physicochemical analysis of the vaccine construct

Various biophysical and biochemical properties of the vaccine construct were assessed in silico using VaxiJen 2.0, Ellipro, AllerTop 2.0, AllergenFP 1.0 and Protparam tool (see supplementary material for more)

Docking studies of vaccine with TLR4, TLR2, MHC I and MHC II

For the development of a brisk and consistent immune response, it is necessary for the vaccine molecule to interact with the target immune cell receptors. To analyze the key interactions, molecular docking studies were performed. In essence, the interaction of the construct with TLR4 natively complexed with MD2 protein, and also the interaction with native TLR2 dimer (PDB ID: 2Z7X) were studied as they reside on the cell surface when they are activated by the vaccine molecule. Further, the peptide molecule was also docked with MHC class I and MHC class II receptors and their crystal structures were retrieved from PDB ID 1I1Y and 1 KG0, respectively.

The docking was performed using ClusPro 2.0 server and after that the various interactions were studied. Firstly, overall binding affinity (ΔG) among the receptor proteins and the vaccine molecule and also the dissociation constant (Kd) at 25.0 ℃. were studied using PRODIGY server. The interaction zone surfaces of the proteins were then keenly analyzed using PyMol for studying the residues involved and the presence of bonding interactions.

Molecular dynamics (MD) and binding free energy calculation of the complexes

MD simulations were performed for the docked complexes of vaccine peptide with TLR4, TLR2, MHC I and MHC II. As depicted in our earlier reports [33], [34], [35], [36], the Amber99sb [37] force field and GROMACS 5.1.4 were employed for the simulation studies. In a cubic box of simple point charge (SPC) water molecules, the complexes were immersed. By adding a suitable number of Na+ and Cl- ions, every system was neutralized. The energy minimization was carried out using the steepest descent method for 50,000 steps to reduce short range bad contacts. The systems were then equilibrated for 5 ns followed by 50 ns MD production runs at 298 K temperature and 1 bar pressure with a 2 fs time step. To check the conformational changes and stability of vaccine bound to TLR4, TLR2, MHC I and MHC II, the root-mean-square deviation (RMSD) and fluctuation (RMSF), the radius of gyration (Rg), hydrogen bond and solvent accessible surface area (SASA) were calculated after the MD runs. The Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) [38], [39] method was used to calculate the binding free energy (see supplementary material for further detail). The binding free energies of the complex were calculated based on 5000 snapshots taken at an equal interval of time from 50 ns MD simulations. More information about the methods used in this study can be found in our previous works [40], [41].

Normal modal analyses of the complexes

In addition of atomistic molecular dynamics simulation, the iMODS server () was employed to accomplish the modal dynamic simulation of the docked complexes. The combined motion of molecule simulated in the normal mode analysis (NMA) considers dihedral coordinates of Cα atoms. The disarticulations of atoms from its structural equilibria were revealed through B-factors. Similarly, deformability plot marks the non-rigid portion of the simulated structure that represents coil or domain linkers usually offer flexibility to the protein. Also, NMA generates a key parameter of eigenvalue, an indicator of stability. Higher eigenvalue is related inversely to the B-value and indicative of a stable structure.

Reverse translation and in silico cloning

For efficient expression of the vaccine, the sequence must be propagated to a bacterial system. In this connection, Java Codon Adaptation Tool (JCat) (http://www.jcat.de) was first used for codon optimization and consequent reverse translation to generate a cDNA sequence of the construct, so optimized for expression in E. coli K-12 strain. The results comprise of GC content and codon adaptation index (CAI) score were determined to predict protein expression levels. Moreover, the optimized vaccine cDNA sequence was inserted into a pET-28a(+) plasmid vector using SnapGene tool 5.2.1.

Immune simulation studies

C-IMMSIM server () was employed for performing the immune simulation of the vaccine to characterize the immune response profile and immunogenicity of the chimeric vaccine. C-IMMSIM is an agent-based model that uses position-specific scoring matrices (PSSM) for peptide prediction derived from machine learning techniques for predicting immune interactions. The simulation was set for 1000 days or 3000 time steps wherein one time step being 8 h. The minimum recommended gap between two doses for most vaccines is 4 weeks. 12 consequent injections of each vaccine comprising 1000 peptide molecules were administered and the doses were kept at time steps of 10, 94, 178, 262, 346, 430, 514, 598, 682, 766, 850, and 934 i.e., 4 weeks apart from each other. On 1100th step, a live self-replicating virus was administered. Simultaneously, in a separate control experiment, only a virus was injected on the 366th day without any prior vaccination. This was done to simulate and compare immune capacity of the host with and without vaccination.

Results

Selection of the antigenic peptides from spike protein for vaccine development

The present study seeks to develop a vaccine that could be efficacious against all the pathogenic coronaviruses pathogenic for human. Studies conducted on SARS-CoV-2 and other coronaviruses, collectively resemble spike protein as the key mediator of virulence, inflammation, cytokine storm and death of the human [12]. In the context of coronavirus-induced immunopathological consequences, TLR4 play is known to play a major role [8]. Therefore, targeting the interaction between the spike protein and TLR4 through developing a vaccine could provide us a useful option for combating the pandemic. Therefore, we compared the spike protein sequence of each pathogenic coronavirus reported to date and also studied the comparative interaction between each spike protein subtype and human TLR4 (Supplementary Figs. 2 and 3, and Supplementary Table 1). The spike protein patches across the different coronaviruses that are involved in spike protein-TLR4 interaction were selected for further studies. Presence of B-cell and T-cell (MHC-I and MHC-II) epitopes were determined and listed in Supplementary Tables 2 and 3.

Construction and validation of the multi-epitope multi-target chimeric vaccine, AbhiSCoVac

Presence of different epitopes against different receptors were predicted and the antigenic segments were combined using specific linker sequences to create a single vaccine construct. Each epitope present in the vaccine was checked for its capacity of developing an effective immune response. The linker sequences were used to separate the different domains of the peptides without disrupting their functionality. Linkers are also important to maintain the stability and flexibility to the protein complex as well as conferring immunogenicity to the vaccine of interest. The MHC-I epitopes were joined by “GPGPG” linker sequences while MHC-II epitopes were linked by “AAY”. To further boost the long-standing effect of the vaccine, a Cholera Toxin B (CTB) adjuvant was further added to the N-terminal of the protein using a “EAAAK” linker. The final construct was found to contain 434 residues with a molecular weight of 46.74 kDa. The coordinate model of the vaccine was constructed using Phyre2 server. Stereochemical quality of the model was validated with the help of Ramachandran plot, Z-score and ERRAT score. The model has a Z-score of −6.62 that lies within the range of scores of proteins of comparable sizes and this indicates towards the reliability of the modelling output [32]. Further, Ramachandran plots were constructed upon the model using Structure Assessment tool of the ExPASy server which revealed that 91.67% of the residues are present in the favoured region of the plot. A model is considered as a reliable one if > 90% of its residue are present in the favoured region. In addition to Ramachandran plot, ERRAT analysis was also performed. Any score above 50 in ERRAT output suggests a good modelling quality. Herein, our model showed a score of 78.534 that enhances our confidence on our modelled structure. Furthermore, various important biophysical, biochemical and immunological characteristics of the vaccine construct were also studied and are presented in Supplementary Table 4.

Docking of vaccine construct with TLR2 and TLR4

Toll-like receptors (TLRs) play a key role in innate immunity as they recognize distinct pathogen-associated molecular pattern (PAMP) present in different pathogens including viruses, and lead to the activation of an adaptive immune response. TLR2 [42], [43] and TLR4 have been recognized of being involved in the detection of SARS-CoV-2 VAMPs resulting in profound release of proinflammatory cytokines. Therefore, an ideal vaccine against SARS-CoV-2 needs to have TLR-binding efficacy. In this study, we investigated the binding efficacy of the vaccine construct against TLR2 and TLR4 using molecular docking. The docking was performed using ClusPro server. While docking the construct with TLR2, the highest-ranking cluster consisted of 49 output models having the lowest energy weighted score of −985.3. The binding was orchestrated mainly by the residues of Tyr389, Asn392, Phe426 and Ser329 on the vaccine protein and the Arg321, Tyr326, Lys373, and Tyr377 on TLR2 via hydrogen bonding. Similarly, the docking results of TLR4 with the vaccine exhibited that the strongest cluster having the score of −1064.7 consisted of 55 models and the complex stability was mediated by hydrogen bonding involving Ser419, Ser333, Leu335, Leu52, and Ser76 on vaccine protein, and Gln39, Lys47, Arg290, Ala291, and Tyr317 on TLR4. PRODIGY server was further used to calculate the binding free energy (ΔG) involved in the formation of both complexes. It was revealed that TLR2 and TLR4 bound to the vaccine molecule possess ΔG of −15.7 kcal/mol and −16.2 kcal/mol respectively. However, the TLR2 complex seemed to have a tiny larger dissociation constant (Kd) of 2.9 × 10-12 M than 1.3 × 10-12 M of TLR4 complex.

Docking of vaccine construct with MHC-I and MCH-II receptors

The vaccine construct comprises of MHC-I and MCH-II epitopes that are likely to bind with their cognate receptors to form immune-complexes that can lead to the activation of cytotoxic T-lymphocytes and helper T-cells to generate an adequate immune response. To analyze the efficacy of the vaccine to encounter MHC-related responses, molecular docking of the vaccine construct with MHC-I and II was conducted. MHC-I receptor bound to the vaccine with a score of −1004.8 in a cluster having 62 members. The bonding was conveyed mainly by Asp362, Arg312, Tyr317, Tyr232, and Thr382 residues on vaccine chain and Arg65, Asn86, Glu89, Lys146, and Gln155 on MHC-I receptor protein. However, the MHC-II protein bound to the vaccine using His420, Thr418, Thr416, Tyr336, and Glu298 on the construct hydrogen-bonded to His16, Tyr78, Arg93, Tyr123, and Trp153 on itself. It reflected a lowest energy weighted score of −1040.2 in a 73 membered top-ranking cluster. MHC-I and MHC-II produced feasible complexes with ΔG of −13.8 kcal/mol and −9.9 kcal/mol. Nevertheless, their Kd differed by a small amount, from being 7.1 × 10–11 for MHC-I complex to 5.3 × 10–8 in MHC-II complex (Table 1 ).

Table 1

Recognition of vaccine peptide by innate sensor proteins of MHC I, MHC II, TLR2 and TLR4, quantified in terms of docking score, binding affinity (ΔG) and dissociation constant (Kd) at 25.0 ℃.

Sensor protein	Lowest Energy Weighted Score	Binding Affinity ΔG (kcal/mol)	Dissociation Constant Kd (M) at 25.0 ℃	Interacting residues
				Receptor protein	Vaccine peptide
MHC I	−1004.8	−13.8	7.1 × 10-11	Glu58, Glu58, Arg65, Asn86, Glu89, Lys146, Gln155, Asn86, Asn86, Thr142, Glu89, Thr73, Lys146, Ala149, Arg65, Glu58, Glu154	Arg384, Lys403, Asp362, Thr382, Tyr232, Arg312, Tyr317, Leu335, Lys261, Gly262, Arg290, Arg312, Ser333, Ser334, Tyr341, Arg343, Asp362, Arg384, Lys403, Ser422
MHC II	−1040.2	−9.9	5.3 × 10-8	His16, Arg93, Tyr123, Trp153, Arg93, Tyr78, Asn82, Asp152, Asp57, Phe13, Leu11	His420, Glu298, Thr418, Thr416, Glu298, Tyr336, Arg343, Arg346, Asn392, Ala424, Phe426
TLR 2	−985.3	−15.7	2.9 × 10-12	Arg321, Tyr326, Asp327, Ser329, Asp288, Glu321, Cys343, Ser345, His340, Asn367, His370, Asp288, Gly313, Gly285, Gln286, Arg290, Cys343	Tyr389, Asn392, Phe426, Lys373, Tyr377, Ala387, Ala388, Ser391, Asn392, Gly394, Ala395, Tyr401, Ser422, Tyr425, Phe426, His420, Ser422, Tyr380, Ser391
TLR 4	−1064.7	−16.2	1.3 × 10-12	Gln39, Lys47, Leu52, Ser76, Asn51, Asp50, Glu27, Asn44, Pro28, Met40, Gln39, Glu144, Met41, Glu42, Met40	Ser419, Ser333, Leu335, Arg290, Ala291, Tyr317, Gln344, Tyr377, Tyr380, Tyr414, Thr417, Ser419, His420, Ser422, Tyr425, Thr417

MD simulation of the docked complexes

The importance of MD simulation in vaccine design is undeniable as it provides information that can direct experimental validation. It has great significance to analyze the internal motions, conformational changes, stability, etc., of the complexes [40], [41]. Using the MD trajectories, RMSD, RMSF, Rg, hydrogen bond, SASA and binding free-energy of the complexes were computed to investigate their structural stabilities and binding strengths.

RMSD is the measure of stability where a larger RMSD value is indicative of the lower stability of a protein complex. The RMSD plots against simulation time for the complexes of vaccine with TLR4, TLR2, MHC I and MHC II are shown in Fig. 7 A(I). The lower RMSD values with less fluctuation suggest stability of all of the complexes. The average RMSD as shown in Table 2 for the vaccine-TLR4, vaccine-MHC II, vaccine-TLR2 and vaccine-MHC I complexes are, respectively, in the following order: 1.7 nm > 1.1 nm > 1.09 nm > 0.97 nm. The least average RMSD as well as fluctuation of vaccine-MHC I [Fig. 7A (I), black curve] imply its largest stability amongst all the complexes. Apart from slight fluctuation of the vaccine-TLR4 complex, other complexes show minimal fluctuations during the whole MD simulation period.

Fig. 7

Shows the comparative picture among A. Immunogenecity profile of the vaccine in addition to administration of live virus on 366th day. And B. Control experiment involving administration of live virus at 366th day without any prior vaccination. In both cases, (I) shows antigen count per ml, (II) Antibody titers (III) B-cell population per mm3(IV) Population of helper-T cells per mm3(V) Macrophage population per state per mm3 and (VI) shows the population of cytotoxic T-lymphocytes per state per mm3.

Table 2

The contributions of the van der Waals, electrostatic, polar solvation and SASA to the total binding energy (in kJ/mol) of the MHC-I-vaccine, MHC-II-vaccine, TLR2-vaccine and TLR4-vaccine complexes.

Stabilization parameters	Complexes of the vaccine peptide and different sensor proteins
	MHC-I	MHC-II	TLR2	TLR4
van der Waals energy	−213.6 ± 5.0	−134.6 ± 5.1	−251.9 ± 10.1	−244.1 ± 9.1
Electrostatic energy	−326.7 ± 16.2	−55.5 ± 5.0	−249.3 ± 16.3	−327.6 ± 14.1
Polar solvation energy	435.0 ± 14.0	152.5 ± 6.0	333.6 ± 25.9	387.6 ± 20.7
SASA energy	−35.6 ± 0.6	−17.8 ± 0.5	−31.5 ± 1.1	−30.3 ± 0.6
Total binding free energy (kJ/mol)	−141.4 ± 11.7	−55.2 ± 8.0	−199.4 ± 18.9	−213.8 ± 19.2
Average RMSD (nm)	0.97	1.1	1.09	1.7
Average Rg (nm)	4.5	5.8	5.1	5.1
Average SASA (nm2)	653.3	541.3	793.1	995.1
Average H-bond	6	9	5	6
Average RMSF (nm)	0.57	0.85	1.1	1.7

Rg describes the level of compaction of a protein. Similar to RMSD, the least fluctuation of Rg can be observed in the vaccine-MHC I complex [Fig. 7A (II), black curve], indicating its most compact behavior due to the strongest interaction. Although the Rg values of vaccine-MHC II are higher, the complex remains stable throughout the simulation period. The average Rg value of the vaccine-MHC II, vaccine-TLR4, vaccine-TLR2, and vaccine-MHC I complexes are, respectively, in the order: 5.8 nm > 5.1 nm ∼ 5.1 nm > 4.5 nm.

SASA accounts for the surface area of a complex which directly interacts with the solvent molecules. The increase in SASA denotes expansion of the protein backbone. The average values of SASA of the vaccine-TLR4, vaccine-TLR2, vaccine-MHC I and vaccine-MHC II complexes [Fig. 7B (I), Table 2] are, respectively, in the order: 995.1 nm2 > 793.1 nm2 > 653.3 nm2 > 541.3 nm2, suggesting more solvent can access the vaccine-TLR4 complex.

Hydrogen bonds in protein–protein interaction are very crucial to check the stability of the complexes. From Fig. 7B (II, red curve), it can be observed that the vaccine-MHC II complex possesses the highest number of hydrogen bonds that retained throughout the MD simulation. Likewise, the hydrogen bonds of other three complexes were also retained throughout the MD simulation, indicating their stability. The average numbers of hydrogen bonds of the vaccine-MHC II, vaccine-MHC I, vaccine-TLR4 and vaccine-TLR2 complexes [Fig. 7B (II), Table 2] are 9 > 6 ∼ 6 > 5, respectively.

RMSF is a dynamic parameter that measures the average main chain flexibility at each residue position. To know flexibility, fluctuations of each amino acid residue presented by the RMSF plots as a function of residue number are shown in Fig. 7C (I). A large RMSF value indicates flexibility of the protein structure, loose bonding, or the presence of loops; contrarily, a small value indicates stability and also the presence of secondary structures such as sheets and helices. While analyzing the vaccine-bound complexes, we observed RMSF to be in a similar trend of RMSD. The RMSF values infer that the vaccine-MHC I complex is the most stable and the vaccine-TLR4 complex to be the least stable and more fluctuating among all the complexes. The average RMSF as shown in Table 2 for vaccine-TLR4, vaccine-TLR2, vaccine-MHC II and vaccine-MHC I follows the trend 1.7 nm > 1.1 nm > 0.85 nm > 0.57 nm, respectively.

The binding free energy (ΔG) has been computed with MM/PBSA [39], which is considered as a very efficient and reliable approach to study crucial molecular recognition processes [39], [40], [41]. Currently, this method has been widely applied in various studies including stability, protein–ligand and protein–protein interactions, mutation, drug design and re-scoring [40], [41].

Table 2 presents the contributions of different interactions to the overall binding free energy of the complexes appearing either positive or negative. Decomposition into separate energy terms manifests that the polar solvation energy opposes binding of the complexes significantly, and thereby reducing the total binding energy, Fig. 7C (II). ΔG provided in Table 2 and Fig. 7C (II) for the vaccine-TLR4, vaccine-TLR2, vaccine-MHC I, and vaccine-MHC II complexes shows the following order −213.8 kJ/mol > -199.4 kJ/mol > -141.4 kJ/mol > -55.2 kJ/mol, respectively. Some inconsistency of the MD parameters in the case of vaccine-TLR4 complex may be due to the large size of TLR4 in which the fluctuating regions are generally consists of the non-interacting residues of less importance.

Overall, the results of MD simulations and binding free energies strongly suggest the stability of the complexes of vaccine with TLR4, TLR2, MHC I and MHC II.

Normal modal analyses of complex stabilization

NMA is a form of the analyses whereby the internal coordinates of the vaccine-receptor complex reveal insights about protein mobility and stability. We conducted NMA analysis to verify the stability of TLR4-vaccine complex using iMODS server and the results are depicted in Suppl. Fig. 4. Eigenvalue is a key indicator of stability in a given structure and higher the value the better [44]. For the first mode, the four vaccine-receptor complexes analyzed exhibited values ranging from 2.043587 × 10-6 as in vaccine-TLR2 complex to 6.657389 × 10-6 as in vaccine-MHC I complex, moreover, the values increase with later modal vibrations which is a clear indicator of high stability of the complex. Nevertheless, the fluctuation depicted in the deformability plots indicates the residues forming coiled structures that impart flexibility to the complexes. Additionally, the covariance matrices described the occurrence of anti-correlated, uncorrelated and correlated atomic motion signified with blue, white and red colour respectively [45].

In silico cloning of AbhiSCoVac

In silico cloning was performed to analyze the expression of the vaccine protein in bacterial system like E. coli. Thus, it was essential to have a codon optimized copy of the reverse translated cDNA of the vaccine. JCat tool was used for the same purpose and a 1302 long cDNA sequence of the vaccine construct was generated for E. coli K-12 strain [46]. For effective expression of a protein within E. coli, the codon adaptation index (CAI) value should be > 0.8 while the GC content must be around 30–70% [47]. Our construct shows a CAI of 1.0 and GC content of 53.76%, this reflects proficient protein expression in the bacterial cell system. Finally, a recombinant plasmid was designed by inserting the codon adapted cDNA sequences into pET-28a(+) plasmid [48] using SnapGene tool ver. 5.2.1.

Immune simulation of the vaccine

C-ImmSim server was used to simulate the immunogenic profile of the designed vaccine, AbhiSCoVac. It was straightaway apparent that secondary and tertiary responses developed in the simulation were much higher as compared to the primary response. The immune response patterns also revealed a decrease in the concentration of antigens against the normal high immunoglobin activity i.e., IgG1 + IgG2, IgM, and IgG + IgM antibodies. Multiple long-lasting B-cell isotypes were also detected which indicated towards the isotype switching potentials as well as memory formation [49]. Similarly, helper T-cells and cytotoxic T-cell populations reflected a higher response with pre-activation of T-cells following vaccination. Dendritic cells and natural killer cells activity have been found to be consistent with the higher macrophage activity levels. Additionally, high levels of IFN-γ and IL-2 after administration of the vaccine suggest a strong immune response. Once the vaccination course ended, a simulated live virus was injected at the 366th day, however, there was no antigenic increment detected. This suggests an efficacious immune response exerted by the specific antibodies developed due to vaccination. In contrast, a control simulation set was studied by injecting the virus on 1100th step without any prior vaccination. The results were quite different as the host was absolutely incapable to contain the viral infection cycle that further resulted in a mass surge of viral load.

Discussion

Immunoinformatics approaches have been found extremely useful in developing novel vaccines to combat numerous fatal diseases. In recent times, in silico approaches have been exploited to develop multi-epitope vaccine peptides against several viral pathogens including dengue virus, ebola virus and most importantly the SARS-CoV-2 virus [50], [51], [52], [53]. Interestingly, a number of vaccine candidates against coronaviruses designed thorough immunoinformatics approaches have shown promising efficiency in ex vivo conditions [54], [55], [56], [57]. However, they are limited to offer protection against a single coronavirus type at a given time. Recurrent emergence of the viral family has been the consistent and pernicious threat to mankind in last two decades. Particularly, millions of lives were shaded due to the pathogenic outbreaks ranging from SARS to MERS and COVID-19 [58]. In this connection, our study has attempted to perform a feat never made before - to design a novel vaccine candidate that can enable armor against all the pathogenic human coronaviruses except hCoV-NL63 species which is avirulent to healthy adults.

Earlier studies have already shown that spike ‘S’ protein of the coronaviruses acts as a credible signatory PAMP for recognition of the viral infection by the surface TLRs of human, especially TLR4 [8], [9]. Taking clue from the studies, we began with the molecular docking of spike proteins of all the six virulent coronaviruses with human TLR4. After examining the interacting complexes comprising viral spike protein and TLR4, energetically most stable and naturally favorite conformations were chosen for further study. The amino acid sequences of spike protein responsible for binding to TLR4 in the conformations were thereby extracted and used for determining the presence of epitopes. B-Cell Epitope Prediction tools from the Immune Epitope Database (IEDB) and BepiPred 2.0 were used to identify the B-cell epitopes, while IEDB T-cell Epitope Prediction tools for recognizing the MHC I epitopes and ProPred server was used to detect MHC II epitopes. Further, different immunological filters were used to screen the predicted epitopes for their antigenicity, immunogenicity, and binding efficacy. Based on the epitopes, the vaccine construct was designed.

The chimeric vaccine peptide was constructed using glycine-proline rich GPGPG and alanine-rich AAY linkers to connect the epitopes at respective positions. Further, a EAAAK linker was employed to connect the N-terminal of the peptide to a Cholera Toxin B (CTB) adjuvant to enhance the recognition efficacy and overall immunogenicity of the vaccine. Intriguingly, the novel vaccine was found to be a non-allergen by AllerTOP v.2.0 server and the result was further verified by AllergenFP v.1.0. The physicochemical properties were analyzed using the ProtParam tool offered by the SWISS ExPaSy server. The molecular weight of the peptide was 46.748 kDa and the Guruprasad’s instability index was evaluated to be 26.94 which classify the vaccine to be highly stable. Generally, a protein whose instability index is lesser than 40 is predicted to be stable and values above that predict the protein as unstable. The theoretical pI of the vaccine was calculated to be 9.58. The Grand average of hydropathicity (GRAVY) index of the vaccine was − 0.534 which indicated towards the polar nature of the vaccine and its effective interaction with water resulting in high solubility. The aliphatic index of 56.98 of the vaccine candidate indicated the protein to be thermostable. Also, the predicted half-life of the vaccine was evaluated to be 30 h (mammalian reticulocytes, in vitro), > 20 h (yeast, in vivo) and >10 h (Escherichia coli, in vivo) which collectively denotes the time taken by the protein to reach 50% of its concentration after its synthesis in the cell. Further, the coordinate structure of the vaccine construct was validated using Ramachandran plot exhibited 91.67% amino acid residues within the favorable zone. The ERRAT analysis had a good score of 78.534 that further validated the overall quality of the designed vaccine. In addition, Z-score assessment by the ProSA web server also revealed a score of −6.62 indicating that the theoretically modelled protein structure is at par with the crystal structures determined by NMR and X-ray.

The major aim of designing the vaccine candidate was to target TLR of human. TLR4 is the principal sensor of SARS-CoV-2 spike protein in humans and it is expressed in various immune cell types such as macrophages, granulocytes, monocytes and immature dendritic cells. In order to trigger TLR4, we incorporated CTB peptide in the vaccine construct. Direct interaction between TLR4 and CTB stimulates TLR4 activity to induce NF-κB activation [59]. On the other hand, the dimer of TLR1/TLR2 has been found to recognize the viral envelop glycoproteins [60]. In TLR-mediated signaling pathways, MyD88 acts as the primary adaptor to guide NF-κB and MAPK activation leading to the generation of proinflammatory cytokines. All these information were indeed useful during the design of the vaccine. The interaction map of the vaccine peptide with TLR4 and TLR1/TLR2 as well as MHC I and MHC II receptors were analysed by molecular docking studies. The docking study of TLR4 with the construct exhibited strongest binding affinity of −16.2 kcal/mol via interacting residues of Ser333, Leu335, Leu52, and Ser76 on the vaccine protein, and residues of Gln39, Lys47, Arg290, Ala291, and Tyr317 on TLR4. However, strong binding propensities were also followed by TLR2, MHC I and MHC II. Additionally, the molecular dynamics simulation and binding free energy calculation demonstrated the stability of all the vaccine-target complexes.

To ensure the efficient expression of the vaccine in E. coli host, codon optimization of the vaccine peptide was performed. The linear vaccine construct was reverse translated into its complementary DNA sequence wherein the GC-content was found to be 53.76% reflecting a high possibility of efficient expression of the vaccine candidate in E. coli host. Further, in silico cloning methods were used to insert the vaccine sequence into the expression vector of pET-28a(+) for perpetuation and expression in the bacterial system.

In addition, the present study also involves computational simulation techniques that use agent-based immunoinformatics model to simulate the deterministic immunological behavior, avoiding the mean field approximation and the effect of the spatial distribution. The simulation model is based upon only published epitopes as well as their corresponding binding strengths and it is an excellent method to screen different vaccine candidates differentiated by their respective immunogenicity [61]. Similarly, the molecular dynamics simulation studies tend to emulate the absolute atomistic level molecular behavior in near-cytosolic conditions by employing the relevant force-fields and water model, nevertheless, it awaits physiological validation by consequent in vitro and ex vivo studies to be made in near future [62].

Conclusions

Hitherto, the coronaviridae family constitutes several virulent viruses that have repeatedly threatened mankind in the past two decades with pandemics like SARS, MERS as well as the COVID-19 outbreaks. Thousands of different virulent strains are positioned across Animalia and waiting to jump species to humans. A universal vaccine is the need of the hour. In this study we employed immunoinformatic approaches to design AbhiSCoVac - a multi-epitope multi-target chimeric peptide that would be able to generate protective immunity against all six virulent members of the family hCoV-229E, hCoV-HKU1, hCoV-OC43, SARS-CoV, MERS-CoV as well as SARS-CoV-2. The designed vaccine was found to be highly stable, antigenic and immunogenic. Moreover, AbhiSCoVac being a precision-directed vaccine, it employs the epitopes present in conserved structural domains of the spike proteins of each virus species. Consequently, the paratope diversity generated by it is hoped to stochastically retard virtually any attempt of infection made by any present or future member of the family including diverse number of strains evolved in the course of pandemic, essentially, protecting against the breakthrough infections.

CRediT authorship contribution statement

Abhigyan Choudhury: Data curation, Formal analysis, Investigation, Methodology, Writing-original draft. Parth Sarthi Sen Gupta: Data curation, Formal analysis, Investigation, Methodology, Writing-original draft. Saroj Kumar Panda: Data curation, Formal analysis, Investigation, Methodology. Malay Kumar Rana: Formal analysis, Investigation, Methodology, Writing-review & editing. Suprabhat Mukherjee: Coneptualization, Formal analysis, Investigation, Methodology, Project Administration, Resources, Supervision, Writing-review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.