A computational study to disclose potential drugs and vaccine ensemble for COVID-19 conundrum

2.1. SARS-COV-2 N protein preparation

The study commenced with the retrieval of the crystal structure of N protein from the PDB database available with a PDB tag of 6M3M [17]. The structure was determined by X-ray diffraction up to a resolution of 2.70 Å. The 3D structure was then treated in UCSF Chimera [24] for structure editing to optimize receptor energy. During the process, the protein first underwent processing using a steepest descent algorithm keeping the step size to 0.02 Å and number of cycles to 100. Afterwards, the conjugate gradient algorithm was applied on the structure for the default ten steps. Both algorithms tended to clean the structure by improving local interactions, fine-tuning widespread structural errors and moving the structure towards local minima. Next, missing hydrogen atoms were added to the receptor and charges to standard and non-standard residues were assigned by mean of AMBER ff14SB [25] and AM1-BCC [26], respectively. Both initial crystal N protein and post-treated minimized structures were examined using a PDBSum PROCHECK analysis [27] to assist in the selection of starting receptors in a high throughput anti-viral scaffold screening process.

2.2. Ligands preparation

The Asinex anti-viral library delivers a meaningful starting point by arranging chemical entities of potent anti-viral activities and good safety profile for the discovery of new powerful leads. The compounds are also easy to access and purchasable for testing in experimental assays. The anti-viral library was retrieved in SDF format and subsequently filtered in PyRx software [28] to select compounds that fulfil the Lipinski rule of five [29]. As per this rule, only compounds with molecular weight ≤ 500 Da, MlogP ≤4.15, N or O ≤ 10, and NH or OH ≤ 5 are selected. The primary library had 6827 compounds that, after deletion of molecules violating Lipinski rule of five, reduced to 4860. This new list of compounds was minimized for energy minima and converted to .pdbqt format to be ready for docking study.

2.3. Uncovering potential epitopes for SARS-COV-2 N protein

The amino acid sequence of the SARS-COV-2 N protein was scanned for potential B and T cell epitopes capable of evoking strong but protective immunological responses. To accomplish this objective, we employed the Immune Epitope Database (IEDB) [30] where first linear B cell epitopes were predicted using Bepipred Linear Epitope Prediction 2.0 setting cutoff score of >0.5. In parallel, T cell epitopes were also predicted starting with MHC-I alleles using IEDB recommended 2020.04 (NetMHCpanEL 4.0) considering a reference set of alleles (S-Table 1 ). The epitopes were sorted on IC50 score basis, and those with score < 500 nM were selected. The MHC-II alleles prediction was accomplished by selecting IEDB 2.22 method and full HLA reference set (S-Table 1). Each of the predicted epitopes was then evaluated for their potential of evoking immunological response by scanning the epitopes in VaxiJen (cut off antigenic score ≥ 0.4) [31]. Next, the shortlisted immunodominent candidates were checked for allergenicity via AllerTOP version 2.0 [32]. Non-allergens were then toxicity filtered using ToxinPred [33] to opt for non-toxic epitopes. Additionally, the non-homologic and virulent nature of the epitopes was verified by using the epitopes in Blastp search against human (Homosapien: 9606) and VirulentPred [34], respectively. Lastly, IFN-gamma inducing epitopes were filtered via IFNepitope web server [35].

2.4. Epitope conservation and population coverage

All entries of the SARS-COV-2 N protein available in the NCBI COVID-19 datahub were retrieved and used to examine predicted IFN epitopes conservation through an online IEDB conservancy analysis tool [36]. It was important to include population coverage of the epitopes in the study as it gave clear directions about the percentage of a specific population likely to respond to the epitopes. This was achieved by employing the IEDB epitope conservation analysis tool [37].

2.5. Construction of vaccine ensemble

The final set of conserved epitopes that provided the greatest population coverage were fused to each other, followed by addition of an adjuvant molecule to the epitope peptide. The 3D structure of the full ensemble was created ab initio using 3Dpro [38]. Loops of the model structure were refined to strengthen structural stability using GalaxyRefine version 2 [39]. Quality assessment of the ensemble 3D model was made employing free available tools: Ramachandran plot of PDBsum [27], ERRAT [40], VERIFY-3D score [41], PROSA Z-score [42]. Then, the structure was minimized to the local energy minima and passed through downward analysis of vaccine design. Different physicochemical parameters of the vaccine were predicted using ProtParam [43]. Solubility and aggregation-prone regions of the vaccine were predicted using Protein-Sol [44] and Aggrescan3D 2.0 [45], respectively.

2.6. MHC clustering analysis and immune simulation of the vaccine ensemble

During the vaccine design process, identification of candidates showing an optimized affinity for a wide range of MHC HLA alleles is fundamental. The MHC clustering analysis was performed using MHCcluster v2.0 [46]. Further, the immune response profile of the vaccine construct was understood using an agent based immune simulator - the C-ImmSim server [47]. The server employed a position specific scoring matrix to spot immune dominant epitopes and machine learning methods to elucidate immune interactions. Most of the simulation parameters were treated as default: adjuvant = 100, number of antigen injection = 1000, random seed chosen = 12,345, and the vaccine was injected with non LPS. The 1000 units of antigen are considered suitable to induce an appropriate immune response to the viral antigen [48]. The simulation steps allowed are 1100 and volume to 110. C-ImmSim server has been successfully applied in several studies to understand host immune system dynamics in response to an antigen [[49], [50], [51], [52], [53]].

2.7. Molecular docking of antiviral ligands and vaccine ensemble

A blind docking approach was applied for docking both anti-viral ligands with the N protein and vaccine ensemble with TLR3 (PDB tag: 1ZIW). In the case of anti-viral ligands, docking was performed with AutoDock4 [54] allowing the central XYZ dimensions search to be restricted to 9.9921 Å on the X-axis, −3.8536 Å on the Y-axis and − 12.7487 Å on the Z-axis. This gives dimensions on the XYZ plane as 31.6231 Å, 46.5257 Å, and 43.2407 Å, respectively, and as net, the whole surface of the receptor molecule was covered to allow ligand molecules to bind freely to the hotspot points of the N protein. Each ligand molecule was docked 10 times to the receptor and the best binding pose was selected by looking for the one with best binding affinity score in kcal/mol (more negative score indicates good binding affinity). The vaccine ensemble docking with innate immune receptor (as a test case here we used TLR3) was assessed in Patchdock [55], the generated complexes were refined with FireDock [56], and the best complex with minimum global energy was considered for visualization and assay [[57], [58], [59], [60]]. Both AutoDock4 and FireDock are best for performing blind docking and generating intermolecular poses that bind best to each other, hence achieving highly stable complexes. Both sets of software are the most citable forms of docking software, widely used and are freely available [56,[61], [62], [63], [64], [65]]. A blind docking approach was used in the present work by providing the complete surface of receptors (N protein in case of drug molecules identification and TLR3 in case of vaccine ensemble docking). This overcomes the limitations of the specific docking that, in majority cases gives false-positive results [66]. Additionally, we proved the docking procedure by docking the co-crystalized ligands to the receptors through the same procedure used in virtual screening of Asinex antiviral library and vaccine ensemble docking to the TLR3 [67]. Both docking protocols results revealed coherent results, thus validating the docking protocol. Further, we employed widely accepted and more accurate molecular dynamics simulation and binding free energies methods to validate the good affinity of the drug molecules and vaccine ensemble [[68], [69], [70]].

2.8. Molecular dynamics simulation assay

Molecular dynamics simulations for both top complexes of anti-viral ligand with N protein and vaccine ensemble with TLR3 were performed using AMBER18 package [71]. The N protein parameters were generated using a ff14SB force field [25] whereas an Amber force field (GAFF) [72] was chosen for anti-viral ligands. To record topology files of both complexes, the leap module [73] was employed. The systems were neutralized by adding an appropriate number of counterions. Next, both systems were submerged in a TIP3 water box, allowing a padding distance of 12 Å. The waterbox with submerged N protein-drug complex and TLR3-vaccine complex are depicted in Fig. 2 . Systems minimization was achieved by running 1500 rounds of steepest descent and conjugate gradient to clean the complexes for unfavorable structural clashes. Heating of systems was done for 100 ps with a gradual increase from 0 K to 300 K, applying pressure of 1 atm. Afterheat, systems were equilibrated for a time period of 100 ps at a constant temperature of 300 K. A production run for each system was completed at time scale of 100 ns. In the process, SHAKE algorithm [74] was applied to constrain all covalently bonded hydrogen atoms of the systems. Periodic boundary conditions were used in the solvation box by the canonical ensemble. Temperature was kept constant at 300 K using a Langevin thermostat [75] and the non-bounded interaction threshold was treated as 8.0 Å. Ewald simulations were utilized for long range interactions. Structural dynamics of the complexes were elucidated through several statistical parameters analyzed through the CPPTRAJ module [76] of the AMBER. Visualization of the snapshots was done by means of UCSF Chimera [24] and VMD software [77].

2.9. Binding free energy estimation

Estimation of binding free energies for biomolecular complexes is a good way of validating intermolecular strength of interactions and shedding light on the dominancy of a particular chemical energy contributing to overall stability. The MMPBSA method of AMBER is an easy and relatively straight forward approach for quantifying binding free energies of ligand(s) docked to a receptor, though this method does not account for entropy contribution [78]. The MMPBSA binding free energy was estimated using equation given below,

where,

ΔG is the net binding free energy for a given system by subtracting the sum of receptor and ligand combine binding energy from the complex. ΔG_MM is the gas phase energy change estimated by molecular mechanics and consists of van der Waals and electrostatic energy. ΔG_solv is the solvation free energy change and is the product of polar and non-polar energy. The latter term is calculated via the solvent accessible surface area (SASA).

2.10. Waterswap absolute protein-ligand binding free energies

The Waterswap method implemented in the Sire package permits estimation of absolute protein-ligand binding free energies [79]. The reaction coordinate is constructed during the process which swaps the protein bound ligand to an equivalent volume of water present in the binding pocket. This method uses all available processing cores of a computer node, slow in-process and converges the free energy averages that took at least five days. The average binding free energy is estimated simultaneously through thermodynamic integration (TI), free energy perturbation (FEP) and Bennetts Acceptance Ratio (BAR) methods [70]. Waterswap was run on simulation trajectories of last 10 ns to estimate the absolute binding free energy of the systems.

3.1. SARS-COV-2 N protein and anti-viral library minimization

A primary phase of protein minimization was applied to SARS-COV-2 N protein to lower its overall potential energy. This was necessary to make protein conformation as close as possible to natural biological systems that are dynamic and low in potential energy to ease spontaneous interactions. However, a minimization event may introduce bad contacts in the structure and disturb the conformation of the molecule, which in turn might affect the compound ranking in the docking-based virtual screening process. The minimization process revealed that minimized N protein is a bit better compared to the pre-minimized original N protein structure. Ramachandran plot investigation showed that both original and energy minimized structures have the same percentage of residue distribution in all four quadrants as can be seen in S-Fig. 1. Statistically, both structures secured 89.4% and 10.6% of residues in the most favoured and additionally allowed regions whereas no residue was plotted in generously allowed and disallowed regions. According to the G-factor assay, the minimized protein had a better overall G-factor score of 0.06, reflecting no unusual features in the structure opposed to the original structure (G-factor score of −0.14). Secondly, the PROSA Z-score was calculated which indicated overall good quality and a non-erroneous structure of the minimized N protein. The Z-score for pre-minimized N protein was −5.06 while it was −5.16 for the minimized form. Based on this evidence, we used the minimized N protein structure in the downward framework.

3.2. Virtual screening

High throughput structure-based virtual screening was performed using SARS-COV-2 N protein as a receptor and drug-like molecules from the Asinex anti-viral library as ligands in the process. As a blind docking strategy was employed, the entire surface of the protein was exposed for binding of the drug molecules. The molecules bound to three different binding pockets of loop region 1, β-sheet core, and loop region 2 of the enzyme as shown in Fig. 3 . The majority of the molecules tended to interact with the pockets of the loop 2 region compared to the loop 1 region and β-sheet core pockets. Top 5 complexes were picked based on the binding affinity of the molecules for the pockets and analyzed. The 2D structures of the compounds are provided in S-Fig. 2. The high affinity top 1 compound with binding energy of −8.9 kcal/mol demonstrated binding to the front pocket of a loop 2 region. This pocket is the result of the long loop of β6 connecting to β7, and β2 ending at the N-terminal site through β1. The bulk of the compound structure (2,5,5,7,9-pentamethyl-2,3,4,4a,5,10b-hexahydropyrano[3,2-c]chromene-3-carboxylic acid) aligned horizontally inside the long loop originating from β6 of the β-sheet core, covering 3₁₀ helix at the acidic wrist of the protein and ending at the carboxyl site. This chemical moiety formed weak van der Waals, pi-pi stacked, and pi-alkyl interactions with the loop of β6. The 4-isopropylphenol moiety was deep in the cavity, making both hydrogen and hydrophobic contacts with the loops. Compound 3 (binding energy: −8.6 kcal/mol) and 5 (binding energy: −8.2 kcal/mol) followed the same pattern of binding with minimal interaction from strong hydrogen bonding and maximum interaction from weak hydrophobic binding. Compound 2 (binding energy: −8.7 kcal/mol) binding mode was stable at loop region 1 and produced a rich network of hydrogen and van der Waals contacts. The 1-(3-(1,2,4-triazolidin-4-yl)phenyl)ethanone part of the molecule mainly drove the interactions of strong hydrogen bonding with Thr92 and Arg94 β5 of β2 and Asp129 of β6 loop. Compound 4 (binding energy: −8.4 kcal/mol) favoured binding to the back pocket of loop region 2, collectively formed by loops of β2, β5 and β6. The 4-formyl-6-methyl-3-methylene-2,3-dihydroisoxazolo[5,4-b]pyridin-7-ium chemical moiety made hydrogen bonds with Pro74 and Asn76 of β2. The rest of the chemical portion reacted with all three loops through van der Waals bonding. The chemical interactions of the compounds are illustrated in Fig. 4 .

3.3. SARS-COV-2 N protein epitopes prediction

A complete protein sequence was first analyzed for continuous B cell epitopes that predicted 6 epitopes ranging in length from 7mer to 28mer of score greater than default cut off (0.5). Such B cell epitopes are the recognition and binding sites of the adaptive immunity B cells. Once activated, the B cell transformed into a mature form, differentiate and produce soluble antibodies. Afterward, the antibodies bind to the epitopes and activate humoral adaptive immunity by triggering formation of neutralizing toxins, labeling the pathogen for destruction via T cell immunity. Therefore, B cell epitopes mapping holds a central role in vaccine designing [80]. The predicted B cell epitopes were examined for T cell epitopes. T cell vaccine induces protective cellular immunity and are significant in targeting mutating viruses like SARS-COV-2. Furthermore, T cell vaccines are regarded more effective than conventional B cell epitopes [81,82]. In the present investigation, both B and T cell epitopes were used to design highly efficacious vaccines. Antigens that are displayed on the antigen presenting cells bound to the major histocompatibility complex molecules are recognized by the T cell receptors present on the surface of T cells. T cell epitopes are surface epitopes via two classes: MHC-I and MHC-II, identified by two distinct T cytotoxic and T helper cells. Several T cell epitopes are reported for each B cell epitopes with IC₅₀ values <50 nM depicting high affinity for the reference set of MHC alleles. Next, the T cell epitopes were subjected to MHCphred assay to filter only those epitopes that associated with great affinity with the DRB*0101 allele, a prominent allele in the human population [83]. All T cell epitopes showed great binding with the DRB*0101 allele with an IC50 value of <100 nM. T cell epitopes then analyzed for their ability of provoking the host immune system. This was achieved by passing the T cell epitopes through the VaxiJen server and selecting those with predicted score higher than the default 0.4. Antigenic T cell epitopes were subsequently filtered through allergenic and toxicity check, and only non-allergenic and non-toxic epitopes were picked. The final set of potential epitopes selected for vaccine ensemble design is listed in Table 1 and their exo-membrane topology is shown in Fig. 5 .

3.4. Population coverage of T cell epitopes

The MHC molecules are extremely polymeric and thousands of MHC molecules are known that are spread across the world population and in different ethnicities. Thus, peptide-based vaccine design that covers the majority of these alleles could provide a wonderful approach for the design of a broad spectrum and highly effective vaccine that does not show any ethnically biased population coverage. The selected B cell derived T cell epitopes are showing average population coverage of 95.04%. This can be divided up as: Central Africa (94.79%), Central America (53.8%), East Africa (97.08%), East Asia (97.08%), Europe (99.96%), North Africa (99.01%), North America (99.89%), Northeast Asia (97.88%), Oceania (97.88%), South Africa (95.27%), South America (95.15%), South Asia (98.7%), Southeast Asia (97.66%), Southwest Asia (95.79%), West Africa (98.43%), and West Indies (99.69%). The number of combined HLA alleles recognized and cumulative percent of population coverage of the epitopes are given in S-Fig. 3.

3.5. Construction of a chimeric multi-epitope peptide vaccine

Peptide-based vaccines are weakly immunogenic, therefore, the final set of epitopes were fused with the help of rigid AYY linkers to produce a multi-epitope peptide [84,85] (Fig. 6A). AYY linkers aid in keeping the epitopes separated and aid in easy recognition and processing of the epitopes by the host immune system. Also, an adjuvant in the form of β-defensin was added to the N-terminal site of the epitope peptide for further strengthening of the immune provoking ability of the epitope's peptide. β-defensin activates lymphokines production which in turn results in Ig production specific to the antigen along with activation of cellular immunity [86]. Ligation of the adjuvant with the multi-epitope peptide was done through a rigid EAAAK linker. This linker is also rigid and promotes easy immune recognition and processing. The complete construct structure was modelled ab initio as no suitable template was available for homology-based modelling (Fig. 6B). The structure has the majority of its residues (80%) plotted in the most favoured region, 12.90% in the additional allowed region, 5.70% in the generously allowed region, and 1.40% in the disallowed regions (Fig. 6C). From secondary structure point of view, the structure maximum of alpha helix (53.0%), 9.6% of 3–10 helix, and 37.3% of beta turns, and gamma turns (Fig. 6D). The overall vaccine ensemble is antigenic with a score of 0.5060, soluble with a probability score of 0.940 and does not harbour any transmembrane helices; thus it is a good candidate for experimental follow up. The weight of the construct is 8.9 KDa, theoretical PI of 9.36, stability index of 39.9, GRAVY score of −0.602 and aliphatic index of 53.01. The vaccine ensemble is also predicted to be highly soluble with predicted scaled solubility value of 0.775 (threshold, 0.45) and does not contain any aggregation-prone regions. The average aggregation score of the vaccine is −1.11, which is far less than the cutoff score of 0.0.

3.6. MHC restricted HLA alleles cluster analysis

The vaccine ensemble molecule interacted with a large numbers of HLA alleles in MHC I class compared to MHC II class. This infers the increased potential of the designed vaccine ensemble to be recognized by the vast majority of polymorphic HLA alleles, maximizing many fold epitopes presentation the immune cells leading to providing strong immune protection capabilities. The MHC I and MHC II HLA clusters are tagged as S-Fig. 4 and S-Fig. 5, respectively.

3.7. Loop modelling and refinement of the vaccine ensemble

Loop regions of the vaccine ensemble were modelled/refined to present its closest possible structure for the downward analysis. The structure was refined in consecutive rounds covering regions: Arg14- Ser34, Asn53- Ser57, Gly63- Ser70, Tyr74-Ala83, and Ala47-Ala48. The output structure was further refined for global and local conformations. An improved vaccine ensemble structure was revealed in which the percent of residues in the Rama favoured regions was increased to 96.3% compared to 93.8% in the original model. The clash error score in the original structure was 29.4 and was found only 2.2 in the refined model. The MolProbity score is relatively low (1.241) illustrating the good quality of the refined model. Overall, the galaxy energy of the refined vaccine ensemble is −1744.48 which is significantly less than the initial structure (132.01), a confirmation of good stability.

3.8. Immune responses simulation of the vaccine ensemble

An in silico immune response model of the vaccine ensemble was generated that illustrated primary, secondary and tertiary immune responses are credible in clearing the pathogen (Fig. 7 ). The IgM + IgG is critical in showing an immune protection response to the antigen whereas IgG2 is the least productive. Similarly, significant interleukin and cytokine reactions are witnessed, IFN-gamma being the most key immune protective factor (Fig. 7). Formation of different isotypes of the B cell immunity to the antigen demonstrates a fundamental role of humoral immunity against the pathogen and subsequent creation of antigen memory. Increased production of cytotoxic and helper T cell populations with corresponding memory formation further affirm the role of T cell immunity complementing B cell immunity in protecting the host from the pathogen. A high population of dendritic cells, as well as that of macrophages, makes it evident that the set of epitopes used in the construction of multi-epitope vaccine ensembles is attractive for activation of all components of the host immune system aiding in combating the pathogen.

3.9. Vaccine ensemble docked conformation with TLR3

The vaccine ensemble was docked with TLR3, an immune receptor from innate immunity that prompts viral recognition and induces type I interferon production. Blinded protein-peptide docking was performed to predict the predominant orientation of the vaccine ensemble with respect to the TLR3 receptor. In total, 100 conformations of the vaccine ensemble were generated, followed by refinement of each complex to opt for the highly stable conformation. Different energy values for the top 10 refined solutions of the docked vaccine ensemble and TLR3 complexes are presented in Table 2 . Solution 9 was ranked as the top solution due to global energy of −7.93 kJ/mol (attractive van der Waals energy: −33.10 kJ/mol, repulsive van der Waals energy: 14.31 kJ/mol, atomic contact energy: 15.51 kJ/mol, and hydrogen bonding energy: −3.66 kJ/mol). Visual interpretation of solution 9 unveils deep central docking of the vaccine ensemble at the TLR3 cavity. The vaccine generates a network of hydrophobic and hydrophilic chemical interactions with TLR3 residues. Ser387, Glu460, Tyr462, Tyr465, and Asn662 are hydrogen bonding contributing residues from the TLR3 with Ser56, Thr68, Asp59, Tyr62, and Arg12 of the vaccine, respectively. Van der Waals interactions from the TLR3 that played a significant role in binding and anchoring of the vaccine are: Lys330, Glu363, Tyr383, Ile411, Thr415, Lys416, Val435, Phe459, Arg484, Arg488, Ile534, His565, Ile566, Phe657, and Asn659. From the vaccine, residues Val13, Arg14, Lys50, Asn55, Ser57, Pro58, Ile66, Ala73, Tyr74, Ile77, and Asn81 participated in van der Waals interactions. The docked conformation of the vaccine ensemble with TLR3 can be seen in Fig. 8 .

3.10. Disulfide engineering and codon optimization of the vaccine ensemble

Disulfide engineering is a directed approach of introducing disulfide bonds in a protein, and is a logical approach of emulating molecular interactions stabiilty required in many industrial and biomedical applications [87]. The chimeric vaccine sequence was subjected to Design 2.0 server which highlighted 5 pairs of residues with high bond energy ranging from 1.26 kcal/mo to 6.75 kcal/mol. The pairs of residues are: Ile2-Ala61, Arg17-Cys40, Val20-Cys41, Cys21-Lys26, and Pro51-Pro62 were mutated to enhanced vacine stability (S-Table 2). The original and mutated vaccine structures are shown in Fig. 6E. Furthermore, the vaccine sequence was reverse translated to improve the codon usage of the sequence as per Escherichia coli expression system to get maximum expression of the vaccine during wet lab experimentation. The codon adaptation index of the vaccine is 1 which is an ideal value for enhance expression. Lastly, the vaccine sequence was cloned into pET28a(+) expression vector (Fig. 6F).

3.11. Evaluation of conformational stability

The conformational stability of the SARS-COV-2 N protein in complex with the top drug molecule and vaccine ensemble was deciphered by running 50-ns of MD simulations. Structural stability of the systems was monitored first by calculating carbon alpha distance of superimposed 50,000 snapshots of the MD simulation (Fig. 9 (top left)). For the drug complex, a small number of minor structure fluctuations were noticed in the protein suggesting the structure has endured less conformational changes during the course of simulation. The average rmsd estimated for the drug complex is 1.38 Å with maximum observed is 2.53 Å. The variations were investigated by visualizing snapshots at regular intervals (0-ns, 10-ns, 20-ns, 30-ns, 40-ns, and 50-ns) and superimposed in UCSF Chimera [88] (Fig. 10A). The 3D alignment revealed an rmsd of 1.3 Å that clearly showed the drug complex stability. In the process, the compound was seen in different conformations where 4-isopropylphenol region of the molecule is showing stability at the initial docked site but the 2,5,5,7,9-pentamethyl-2,3,4,4a,5,10b-hexahydropyrano[3,2-c]chromene-3-carboxylic acid is stretching along the length allowing the molecule entry deep in the cavity with formation of extra van der Waals contacts as depicted by the binding free energy estimation in section. Further, these movements are the result of flexible loop regions of the N protein that propel the molecule to alter its conformation by adjusting the molecule along the channel of the receptor. This depicts binding free energy methods, with the compound enjoying binding at the cavity site and is an attempt at shwing stable binding as simulations progress. In the same way, rmsd analysis of the vaccine ensemble with TLR3 was carried out (Fig. 9 (top left)). The average rmsd of the TLR3-vaccine ensemble is 3.34 Å. Variations in the rmsd are mainly because of the ligand structure move in an attempt to get a highly stable and immune cell recognizable pose, resulting in force applied on the flexible regions of the TLR3. This was affirmed by taking regular snapshots at interval of 10-ns and superimposition in UCSF Chimera that revealed an rmsd of 1.124 Å (Fig. 10B). The receptor rmsd was followed by ligand rmsd to affirm the ligand binding stability with the receptors (Fig. 9 (top right)). An average rmsd of 1.30 Å and maximum of 2.41 Å were reported for the drug compound. Again minor ups and downs in ligand rmsd were spotted, reflecting conformational moves as stated earlier, though these adjustments are in favour of the increased system stability. The vaccine ensemble rmsd is very low (0.15 Å) demonstrating very stable behavior of the vaccine with TLR3. Next, compactness of the N protein structure was tested using the radius of gyration tool (Fig. 9 (bottom left)). The average gyration value of the N protein in the presence of the drug along the simulation is 14.56 Å and maximum value observed is 14.87 Å. This clearly shows that the secondary elements are highly compact in the N protein 3D structure and show stability in the drug molecule presence. The same compact complex system was noted for the TLR3 and vaccine ensemble with a mean radius of gyration value of 33.72 Å. The number of hydrogen bonds formed between the drug/vaccine with the receptors is shown in Fig. 9 (bottom right). On average, the drug molecule is in contact with the N protein docked site via 1 hydrogen bond through the simulation period whereas the vaccine produced a robust network of hydrogen bonds with TLR3 that, on average, is higher in each simulation frame.

3.12. Binding free energy calculation

The MMPBSA is now commonly applied on biological systems to model molecular recognition and is a central focus in molecular simulations [89]. This method uses reasonable computational cost and is in routine use to estimate binding free energies of small ligand molecules bound to a large biomolecule receptor [89]. The MMGBSA analysis, which is complementary to the MMPBSA, revealed binding with considerable affinity to the N protein as depicted by the net total energy of −30.4504 kcal/mol. This net energy can be split into components as: complex (−10,452.14 kcal/mol), receptor (−10,414.58 kcal/mol), and ligand (−7.10 kcal/mol). The gas phase energy and solvation energy contribution to the net energy of the system are significant; the latter dominates by polar energy as the favourable (−15.27 kcal/mol) opposed to minor support from non-polar energy (−4.27 kcal/mol). The gas phase energy is the attribute of van der Waals energy (−37.49 kcal/mol) in contrast to the non-favourable part from electrostatic energy (26.60 kcal/mol). The detailed N protein and drug molecule MMGBSA binding energies of the complex, receptor, ligand and the net energies are provided in S-Table 3 . In the MMPBSA method, the system net total binding energy is −23.59 kcal/mol (delta gas phase (10.89 kcal/mol) + delta solvation energy (−12.70 kcal/mol)). Together, both gas and solvation contributions play a balanced role in the binding affinity of the compound for the receptor. As in MMGBSA, the solvation energy is the output of a good polar solvation role (−9.43 kcal/mol) compared to less non-polar energy (−3.27 kcal/mol) in MMPBSA. The gas phase energy van der Waals component as reported in the MMGBSA is the prime factor in the compound interaction (37.49 kcal/mol) with the protein whereas, due to the lack of formation of ionic moieties that can give the electrostatic nature of interactions, the columbic interaction is insignificant (26.60 kcal/mol). The complete data of the MMPBSA analysis for N protein and drug molecule are presented in S-Table 4. Likewise, the TLR3-vaccine ensemble was highly stable with MMPBSA net binding energy of −47.96 kcal/mol. The stabilizing factor was revealed to be the gas phase electrostatic energy (−1359.83 kcal/mol) in the system lower energy state along with −67.79 kcal/mol of van der Waals energy. The polar solvation energy is the non-favourable contributor (1390.67 kcal/mol) to the net solvation energy (1379.67 kcal/mol), in contrast to little favourable energy from non-polar energy (−10.99 kcal/mol). The full data of MMGBSA and MMPBSA for TLR3 and vaccine ensemble are given as S-Table 5 and S-Table 6, respectively.

3.13. Energy contribution of key receptor residues and ligands

To underline residues that are essential to ligand binding and are favourable in complex stability, per residue decomposition was performed. Common residues of MMGBSA and MMPBSA with negative average binding energy were categorized as essential amino acids and were vital in interactions with the ligand (Table 3). For instance, residues such Trp5, Thr29, Asn30, Ile94, and Ile105 were demonstrated as hot spot residues due to their profound contribution in complex stabilization. For the TLR3-vaccine, Ser387, Glu460, Tyr462, Tyr465, and Asn662 are hotspots in anchoring the vaccine ensemble at the docked position.

3.14. WaterSwap absolute binding free energy

To estimate the absolute binding free energy of the system, WaterSwap was used. Unlike MMPBSA, this uses an explicit water model thus does not suffers from the shortcomings of the usually applied MMPBSA method. The binding free energy is calculated by means of three highly sophisticated approaches including TI, FEP and BAR. The resultant value is considered good when the mentioned techniques are converged with difference of ≤1 kcal/mol. The drug complex produced a TI value of −22.314 kcal/mol, FEP of −22.145 kcal/mol, and BAR of −23.87 kcal/mol. Similarly, the vaccine ensemble in complex TLR3 was reported to have TI, FEP and BAR values of −47.45 kcal/mol, −48.74 kcal/mol, and − 48.91 kcal/mol, respectively.