key: cord-0994277-smb8ptn2
authors: Lamichhane, Tika Ram; Ghimire, Madhav Prasad
title: Evaluation of SARS-CoV-2 main protease and inhibitor interactions using dihedral angle distributions and radial distribution function
date: 2021-10-19
journal: Heliyon
DOI: 10.1016/j.heliyon.2021.e08220
sha: 3bcffefecc8bc27193bdeee3952f089ccf43afc2
doc_id: 994277
cord_uid: smb8ptn2

In order to evaluate the interactions between a potential drug candidate like inhibitor N3 and the residues in substrate binding site of SARS-CoV-2 main protease ([Formula: see text]), we used molecular docking and dynamics simulations. The structural features describing the degrees of folding states of [Formula: see text] formed by beta-barrels and alpha-helices were analyzed by means of root mean square deviation, root mean square fluctuation, radius of gyration, residue velocity, H-bonding, dihedral angle distributions and radial distribution function. All of the residues forming ligand binding domain (LBD) of [Formula: see text] lie within the allowed region of the dihedral angle distributions as observed from the equilibrating best pose of [Formula: see text]-N3 system. Sharp peaks of radial distribution function (RDF) for H-bonding atom pairs (about 2 Å radial distance apart) describe the strong interactions between inhibitor and SARS-CoV-2 [Formula: see text]. During MD simulations, HSE163 has the lowest residue speed offering a sharp RDF peak whereas GLN192 has the highest residue speed resulting a flat RDF peak for the H-bonding atom pairs of [Formula: see text]-N3 system. Along with negative values of coulombic and Lenard-Jones energies, MM/PBSA free energy of binding contributed by the non-covalent interactions between [Formula: see text] and N3 has been obtained to be -19.45 ± 3.6 kcal/mol. These physical parameters demonstrate the binding nature of an inhibitor in [Formula: see text]-LBD. This study will be helpful in evaluating the drug candidates which are expected to inhibit the SARS-CoV-2 structural proteins.

From the origin of COVID-19 until now, many researches have sought to identify the potential inhibitors of SARS-CoV-2 structural proteins. The SARS-CoV-2 proteins have ligand binding domains (LBD) which are targets for the potential inhibitors [1, 2, 3, 4, 5] . The researchers are working clinically as well as computationally on the different types of potential antiviral drugs in order to search for the effective control of the ongoing pandemic. The best binding drugs like inhibitor N3 [1, 2] , -ketoamide derivative [6] , cannabisin A and darunavir [7] execute stable interactions to inhibit the SARS-CoV-2 main protease (M pro ). These types of drugs prevent the viral transcription and replication by interfering with the polypeptide transnational activities [3] . Early research focused on the study of covalent inhibitors targeting the catalytic HIS41-CYS145 dyad which may result in off-target side effects and toxicity [8] . Noncovalent inhibitors are beneficial for the easy reaction pathway of ligand binding and dissociation. Some vaccines have been brought into practice along with the clinical trials to * Corresponding authors.

E-mail addresses: tika.lamichhane@cdp.tu.edu.np (T.R. Lamichhane), madhav.ghimire@cdp.tu.edu.np (M.P. Ghimire). prevent coronavirus transmission. Understanding drug impacts on the molecular level through the analysis of protease-inhibitor interactions is a good way to evaluate effective drugs.

The patients with SARS-CoV-2 are susceptible to pneumoniainduced acute respiratory distress syndrome and multi-organ failure. Human genetics, host immune response, age and comorbidities influence the cytokine release syndrome [9] . The coronavirus consists of a positive-sense, single-stranded RNA genome with about 30 thousand nucleotides which encodes two overlapping polyproteins pp1a and pp1ab. The M pro , also known as 3CL pro , is a 33.8 kDa cysteine protease that conserves functional polypeptides released by proteolytic processing of pp1a and pp1ab [2] . The key enzyme, SARS-CoV-2 M pro plays an important role in gene transcription activity. The M pro -LBD has the catalytic HIS41-CYS145 dyad in between I and II domains comprising of residues 8-101 and 102-184, respectively, and including an anti-parallel beta-barrel [2, 4] . Flexible alpha-helices and strong beta-barrel structures forming three domains of the M pro are shown in Fig. 1 drug candidates with the LBD residues. The peptidomimetic inhibitor N3 with well designed side chains fits finely in the substrate binding pocket. The M pro -LBD residues are: MET49, PHE140, ASN142, GLY143, CYS145, HSE163-164, MET165, GLU166, LEU167, PRO168, HSE172, GLN189, ALA191 and GLN192 as identified from the known structure, 6LU7.PDB [4] .

The computational methods such as molecular docking and dynamics simulations are widely used to identify the binding efficacy of the potential drug and its stability into the target viral receptors. Here, after the co-crystallized docking of inhibitor N3 into the M pro -LBD ( Fig. 1-b) , structural features of the SARS-CoV-2 M pro -LBD have been evaluated by means of physical parameters such as root mean square deviation [10, 11] , root mean square fluctuation [11] , radius of gyration [11, 12] , radial distribution function [13, 14] , dihedral angle distributions [15, 16] and solvent accessibility [17, 18] . Furthermore, the M pro -N3 interactions are evaluated by means of time-based changes of H-bonds, residue speed, coulombic and Lenard-Jones energies, and binding energy calculations during MD simulations [19] . The free energy of binding of the M pro -inhibitor system is obtained by an end-point method, i.e. Molecular Mechanics / Poisson Boltzmann Surface Area (MM/PBSA) calculations [20, 21] . This type of study is important to explore the nature of interactions between the potential drug candidate and the viral structural protein. As a biophysical perspective, the concept of clinical data modeling for the trendline of coronavirus cases and deaths, and SARS-CoV-2 structural features are shortly reviewed by Lamichhane and Ghimire, 2020 [22] .

The primary structure of SARS-CoV-2 M pro -N3 system was taken from the crystallographic data, 6LU7.PDB [4] . The protein was prepared with PyMOL deleting everything except amino acid residues and it was protonated by adding hydrogen atoms into it [23] . Thus prepared M pro structure was uploaded to Swiss PDB Viewer and it was optimized selecting all residues using gromos96 43B1 force field [24] . The structure of inhibitor N3 obtained from 6LU7.PDB was cross-checked with PyMOL, protonated and then uploaded to ACD/ChemSketch 12.0 for its geometrical optimization [25] . The electrically neutral N3 with 99 atoms has the mass of 682.81 amu. The ligand moiety, N3 was docked non-covalently in the co-crystallized form of its original position in 6LU7.PB with the resulting RMSD 1.66 Å as determined from the Dock-RMSD tool ( Fig. 1-b) [26] . The docking software used was AutoDock 4.2 [27] . The AutoDock Tool performs grid-based energy evaluation using Lamarckian genetic algorithm [28] . During the docking simulations, side chains and torsional bonds of N3 were rotated freely by keeping the receptor M pro system rigid. Gasteiger charges [29] were computed by adding polar hydrogens to M pro . Fifty runs of docking simulations were performed with 300 population size, 2,500,000 evaluations and 27,000 generations. To generate the docking box of actual mapping, grid point spacing of 0.375 Å, coordinates of central grid points (x = -9.204 Å, y = 14.323 Å, z = 67.590 Å) and grid (x, y, z) points of (64, 90, 70) were employed in the receptor M pro -LBD. To obtain the best docking pose of N3 in the M pro system, the process was revised five times. The best docked conformer of main protease-inhibitor (M pro -N3) was then subjected to nanoscale molecular dynamics (NAMD) simulations for the purpose of evaluating the stable interactions by means of equilibration techniques.

A protein structure file (PSF) of M pro -N3 system was generated using topology files and NAMD simulations were performed using parameter files of CHARMM 36 force fields [30, 31, 32] . The PSF contains an information of atom types, atomic masses and partial atomic charges. The missing H-atoms in protein data bank (PDB) file of M pro -N3 system were assigned automatically during the PSF generation. The parameter file contains the values of numerical constants which are used to calculate bonded and non-bonded energies of CHARMM potential function [19] . The N3 was obtained from the best pose conformer of M pro -N3 system after the molecular docking simulations. The missing hydrogen atoms were added to N3 and it was subjected to SwissParam server [33] to generate the related parameter and topology files. The related topology files were used to generate PSF of both M pro and inhibitor N3, and both structures were merged to form the full PSF and PDB files of M pro -N3 complex with the help of MergeStructs plugin-1.1 of visual molecular dynamics (VMD) [34] . This method preserves all of the original bonds and existing coordinates in both M pro and N3 molecules so that the original pose of N3 remains unchanged while liganded in the M pro -LBD. Thus prepared M pro -N3 structure assigning with all H-atoms was fully solvated in a periodic box of 41040 TIP3P-water molecules. The system was neutralized with 0.15 mol/L concentration of Na + and Cl − ions in order to provide a cellular environment. The software packages used were VMD-1.9.3 for structure generation and analysis, and NAMD-2.12 for geometrical optimization, equilibration run and production run of MD simulations.

The M pro -N3 water box had the cell basis vectors of 71.38 Å, 86.70 Å and 79.88 Å with center coordinates (-25.83, 12.78, 58.05 Å). Particle mesh Ewald (PME) was active with Ewald coefficient 0.25, grid spacing 1.0 Å and grid dimensions (72, 88, 80 Å). The NAMD simulations were performed in NPT ensemble using Langevin piston temperature, pressure and damping coefficient as 300 K, 1.0 bar and 1 ps −1 , respectively. The MD simulations of M pro -N3 system were performed using a time step of 2 fs and constraining the H-atoms with rigid bonds. The coulombic and van der Waals interactions were evaluated using 1-4 scaling and distance related parameters: 10 Å, 12 Å and 14 Å for switching, cut-off and pair-list distances, respectively. The number of steps implemented for energy minimization or geometrical optimization of M pro -N3 system were 3000 CG steps. The system was subjected to 50 ns equilibrating simulation at the constant temperature 300 K. Verlet algorithm [35] was used to calculate the atomic trajectories. With the help of NAMD-Energy Plugin-1.4, the interaction energies (electrostatic and van der Waals) were evaluated. The residue speed or displacement per frame were observed by means of VMD Timeline-2.3 [34] . The timebased changes of the physical parameters (root mean square deviation, root mean square fluctuation, residue velocity, dihedral angle distributions, radial distribution function, solvent accessible surface area and H-bonding) describe the structural properties of protein-ligand systems as studied in the previous researches [36, 37, 38, 39] . These tools were also used to evaluate the structural changes of the M pro -inhibitor system during the MD simulations. All the interaction visualization analysis were performed by PyMOL molecular visualization tools [23] , VMD tools [34] , BIOVIA discovery studio visualizer [40] and LIGPLOT + [41] .

The radial distribution function (RDF) denoted by ( ) gives the probability of finding the positions of an atom-pair within a specified range of radial distance during MD simulations. The normalized RDF for the atom pairs , within a spherical shell having radius and thickness is defined by Equation (1) [13, 14, 36] .

where , ( ) is the number of occurrences of atom pairs between and + . Likewise, solvent accessible surface area (SASA) is also observed to know the structural properties and nature of interactions of a biomolecule. The SASA is an area described by the solvent molecule when it rolls around the molecular complex approaching to van der Waals contacts of the atoms. It is calculated by using Equation (2) [17, 18, 36] .

where is obtained by adding radius of solvent molecule and van der Waals radius of the atom, is the distance between sphere center and section , is the distance between two consecutive ℎ and ( + 1) ℎ sections, and is arc length computed on the section .

The binding free energy of M pro -N3 complex was calculated by using MM/PBSA method in CaFE1.0 plugin [20] . The CaFE uses NAMD trajectories supporting CHARMM force field parameters. The end-point energy method, MM/PBSA consists of three terms: molecular mechanics (MM) for gas-phase energy, Poisson-Boltzmann (PB) for polar solvation energy and surface area (SA) for nonpolar solvation energy [21] . The free energy of binding, Δ can be estimated by using Equation (3) [20] .

where the MM term, i.e. gas-phase energy (Δ ) was obtained from non-bonded (coulombic and Lennard-Jones) interactions between inhibitor and protein. The PB term, i.e. Δ was calculated by setting interior dielectric constant 2 and exterior dielectric constant 80 in APBS program [42] interfaced to CaFE. The SA term, i.e. Δ was calculated by setting an offset 0.92 kcal/mol and surface tension 0.00542 kcal/mol/Å 2 . Here, the SA term has an approximate linear relation with SASA. The entropy term Δ is ignorable for a larger ligand [43] and also for a group of structurally similar compounds [20] . Due to inconsistency in the results by roughly entropic approximation, CaFE excludes the entropy term and sums the results from MM, PB and SA terms averaged over the saved snapshots of simulated system.

The molecular docking based analysis explains that the antiviral drugs against SARS-CoV-2 reveal strong interactions with higher binding affinity or smaller inhibition constant. From our observations, the inhibitor N3 has binding energy of -9.53 kcal/mol, inhibition constant of 103.5 nM and the best pose RMSD of 1.66 Å with respect to its cocrystallized form in the M pro -LBD (6LU7.PDB). To evaluate the binding stability of the docked complexes, MD simulations are performed and the equilibrated structures are analyzed using the different tools as explained in the methodology section.

The electrically neutral M pro -N3 solvated water box has the mass of 283297 amu consisting of 45903 atoms, 32196 bonds, 22385 angles, 12738 dihedrals, 890 impropers, 304 cross terms, 16179 H-groups and 94305 degrees of freedom. These parameters collectively describe the conformational and binding properties of the receptor-ligand system [19] . The protein structure and stability depends on the number of H-bonding residues in the protein backbone and side chains. As a binary two dimensional matrix, Fig. 2 -a shows the contact map of H-bonding amino acid residues obtained by machine learning technique [38] after 50 ns equilibrations of SARS-CoV-2 M pro system. The contact map is invariant of rotations and translations which distinguishes nature of H-bonding interactions between main chains and/or side chains as indicated by dark, blue and red points in Fig. 2-a. Here, alpha-helices, antiparallel beta-sheet and parallel beta-sheet appear adjacent to diagonal, cross-diagonal and parallel to diagonal, respectively. In this way, root mean square deviation (RMSD), root mean square fluctuation (RMSF) and radius of gyration (RG) describe the conformational stability of the protein backbone forming alpha-and beta-strands and side chain residues. Their time-based changes explain the degree of folding states and compactness of the system by acting as the reaction coordinates. The average RMSD of the viral receptor M pro is 2.40 ± 0.34 Å and that of the inhibitor N3 in LBD is 1.41 ± 0.32 Å provided with standard deviation (SD) as the error value. The small fluctuations seen in the RMSD plots ( Fig. 2-b) are due to the minutely changing number of H-bonds among the protein residues during the equilibration run. The mean RG is 22.70 ± 0.20 Å for M pro and 6.61 ± 0.25 Å for N3. The RMSF of each residue has been shown in Fig. 2 -c whose value lies between 0.3 Å and 1.5 Å. The RMSF of terminal residues SER1 and GLN306, and the intermediate residues TYR154 and ASP155 has larger values than the other residues of showing higher flexibility in the M pro system. Thus, the average measures of RMSD, RMSF and RG resulting in small fluctuations indicate the conformational stability of protein-drug systems during MD simulations.

The M pro residue speed (displacement/frame) ranges from 0.2 Å/frame for CYS160 to 2.51 Å/frame for TYR154 as demonstrated in Fig. 3 -a. In case of H-bonding residues, the speed seems to be slightly fluctuating over the time frame of simulation (20 ps/frame) as shown in Fig. 3 -b. As demonstrated by Ligplot (Fig. 4) , the stable H-bonding residues interacting with N3 at the donor-acceptor distance of about 3 Å are PHE140, ASN142, HSE163, HSE164, GLU166 and GLN192 whose speeds are measured to be 0.87 ± 0.38, 0.95 ± 0.41, 0.70 ± 0.31, 0.73 ± 0.31, 0.80 ± 0.35 and 1.00 ± 0.45 Å/frame, respectively. Here, each of these average values does not exceed 1 Å/frame that signifies the stability of N3 in the viral M pro -LBD.

The hydrogen bonds existing among the amino acid residues preserve the protein conformation. Due to its strict rules governed by dipole-dipole interactions, even breaking of a strong H-bond between M pro helices influences the activities of viral gene transcription. The strong interactions between an inhibitor and M pro residues are ex- pected to restrict the replication phenomena. The H-bonds distribution between polar N and O atoms by the M pro -N3 interactions is 6 ± 1 throughout the simulation time ( Fig. 5-b) . The number of H-bonds was estimated by assuming donor-acceptor distance 3.5 Å and cut-off angle 30 o . Even the energy due to thermal fluctuation makes the Hbonds distribution change. From in-silico docking followed by 50 ns equilibration runs of M pro -N3 best pose, we have found that the M pro -LBD residues forming stable H-bonds with N3 are PHE140, ASN142, HSE163, HSE164, GLU166 and GLN192 (Fig. 4) . According to Jin et al. (2020) [2] , there are seven stable H-bonds while binding the computer designed drug N3 by the SARS-CoV-2 M pro -LBD residues.

The sharp and flat RDF peaks, as shown in Fig. 5 -a, are obtained by the radial distributions of H-bonding atom pairs between M pro and inhibitor N3 during MD simulations. The sharper RDF peaks suggest more potent interactions. A strong H-bonding atom pair: HE2(HSE163)-O7(N3) results a sharp RDF peak in about 2.10 ± 0.55 Å radial distance and a weaker pair: O(GLN192)-H2(N3) results a flat RDF in about 4.20 ± 2.65 Å radial distance. Here, the strongly interacting residue GLU166 forms two stable H-bonds with the inhibitor N3 resulting the sharp RDF peaks. Among the H-bonding M pro -LBD residues, the sharpest RDF peak for HSE163 is also associated with the smallest speed of 0.70 ± 0.31 Å/frame during the equilibrating simulation. Fig. 6 -a shows a plot for dihedral angle distributions which is also known as Ramachandran plot. In this plot, there are 15 amino acid residues of the M pro -LBD whose dihedral angles ( , ) distribute mainly in three quadrants. The characteristic features of the Ramachandran plot are such that the right twisted antiparallel beta-barrels lie in the first quadrant (+ , + ), right-handed alpha-helices fall in the second quadrant (-, + ) and the left-handed alpha-helices fall in the third quadrant (-, -). Because the residue GLY has no side chains, it performs no steric hindrance with the neighboring residues and offers a larger crude velocity. Due to this fact that it is allowed in the white region of the Ramachandran plot [14, 15] . The constraints imposed by backbone H-bonding induce the changes in dihedral angle distributions of alpha-helix and beta-forms. The torsion angles of M pro residues obtained after in-silico docking and 50 ns equilibrating run are plotted in Fig. 6 -b and these values for M pro -LBD residues are listed in Table 1 . The dihedral angles of the amino acid residues fluctuate during the MD simulation as given by the ranges of and in Table 1 . However, all H-bonding residues and most of the M pro -LBD residues lie in the allowed region of the Ramachandran plot.

Coulombic interactions in the biomolecules arise from the partial charge distributions in their atoms and these are mainly related to the electrostatic dipole-dipole interactions. The M pro -inhibitor binding efficacy depends on electrostatic and van der Waals contributions of the LBD residues. It is known that the positively charged amine groups are histidine, lysine and arginine, and negatively charged groups are phosphate, carboxyl and glutamate side chains. The negative term in Lenard-Jones potential comes from the induced dipole moments by fluctuating charge density and positive term originates from repulsive effects of the overlapped electron cloud of nearby atoms. NAMD-Energy Plugin calculates the coulombic and van der Waals interaction energies which are observed as -49.74 ± 6.21 kcal/mol and -58.24 ± 4.94 kcal/mol, respectively for the M pro -N3 system. These values indicate the higher binding efficacy of such protein-ligand complex. Fluctuations on the interaction energies as explained by the inclusion of standard deviation in their mean values are also reflected by the time-based plots shown in Fig. 6 -c. The interaction energies depend on the type of drug candidates and the nature of point mutation occurred in the viral structural proteins.

SASA measures the surface area of amino acid residues that are accessible to the solvent molecules. The SASA of the LBD residues highlights nonpolar solute-solvent interactions in the drug binding site. As plotted in Fig. 7-a, SASA ranges from 0 to 300 Å 2 [17] . The M pro residues have SASA ranging from 174.53 Å 2 for GLY138 to 360.48 Å 2 for TRP31 and such variation can be seen in Fig. 7 -b. Here, the probe radius used to find SASA is 1.4 Å which is equivalent to the radius of a solvent molecule.

Since the drug like N3 relies on the shallow binding cleft of M pro , it has larger surface accessibility to the solvent molecules and executes strong electrostatic interactions. The inhibitor also gets hydrophobic interactions with SARS-CoV-2 M pro forming H-bonds with the LBD residues [1] . Fig. 8 shows the equilibrated M pro -LBD residues interacting with N3 where the residues and the bond with green color represent van der Waals interactions and conventional H-bonds, respectively. In this way, C-H bonds, pi-Sulpher, alkyl and pi-alkyl residues are distinguished with separate colors. The characteristic features of M pro -LBD such as interpolated charge of receptor atoms, H-bondig with donor and acceptor atoms, hydrophobicity and basic as well as acidic ionizability can be understood with color codes as shown in Fig. 9 . The colored surfaces of M pro -LBD in Fig. 9 -c indicate the hydrophobicity of amino acid residues representing hydrophilic (blue) and hydrophobic (brown) interactions.

The calculations based on MM/PBSA method of binding free energy assist to obtain the important biophysical information of drug binding mechanism in the protein-enzyme systems. The MM/PBSA free energy of M pro -N3 binding has been obtained to be -19.45 ± 3.6 kcal/mol for which the MM contributions from coulombic and Lennard-Jones interactions are -47.06 ± 5.01 kcal/mol and -61.16 ± 3.82 kcal/mol, respectively. Here, the PB term is 95.54 ± 4.77 kcal/mol as polar solvation energy and the SA term is -6.77 ± 0.24 kcal/mol as nonpolar solvation energy. Due to the changing nature of dielectric constant ( ) with the protein residues, the MM/PBSA results may be more accurate for slightly higher whose default value is taken to be unity in MM electrostatic calculation [44] . Though there are some evidences of poor correlation between MM/PBSA results and experimental data for the systems like SARS-CoV-2 M pro having shallow binding cleft [45, 46] , our result is in a close agreement with the MM/PBSA free energy values for the non-covalent ligands reported by Muhtar et al. (2021) [47] . According to Gupta et al. (2020), the MM/PBSA free energy for di-hydroergotamine is -17.9 kcal/mol and that for midostaurin is -16.2 kcal/mol. Due to the higher binding affinity, these were reported as potential drug candidates in SARS-CoV-2 M pro [43] .

The N3 like peptidomimetic inhibitors form the most extensive and stable H-bonds and hydrophobic interactions with the LBD residues resulting the decreased flexibility of SARS-CoV-2 M pro . The evaluation based on the results obtained from hydrogen bonding analysis, radial distribution of donor-acceptor pairs, dihedral angle distributions of the LBD residues and free energy calculations suggests that the SARS-CoV-2 M pro stands as one of the most promising targets for the potential drug candidates because no human protease has been recognized as such sort of analogous cleavage procedures.

The SARS-CoV-2 M pro organized by stable beta-barrels and flexible alpha-helices strongly adapts N3 like drugs in its binding domain. The small fluctuations seen in RMSD and RG values explain the conformational stability of M pro -inhibitor system during the molecular docking and dynamics simulations. The RMSD of both receptor M pro and inhibitor N3 lies below 3 Å and the RG values are 22.70 ± 0.20 Å and 6.61 ± 0.25 Å, respectively. The RG with small standard deviation measures the compactness of the protein-ligand system. Also, the RMSF values of M pro residues lie below 1.5 Å. Residue speed or displacement/frame is the least (0.2 Å/frame) for CYS160 and the greatest (2.51 Å/frame) for TYR154. The sharpness of RDF curve about small radial distance of H-bonding atom pairs reflects the nature of strong interactions between inhibitor and SARS-CoV-2 M pro . A way of evaluating the binding mechanism of an inhibitor in the M pro -LBD is to study the nature of RDF curves. In the M pro -N3 binding, HSE163 has a sharp RDF for the strong H-bonding atom pairs within 2.10 ± 0.55 Å radial distance, whereas GLN192 has a flat RDF for the weak H-bonding atom pairs within 4.20 ± 2.65 Å radial distance. Being adapted in the M pro -LBD, N3 acquires hydrophobic interactions along with coulombic energy -49.74 ± 6.21 kcal/mol and van der Waals energy -58.24 ± 4.94 kcal/mol. The MM/PBSA free energy of M pro -N3 binding has been calculated to be - 19.45 ± 3.6 kcal/mol as the contribution of non-covalent interactions.

In conclusion, the physical parameters explained in this research will be useful tools for evaluating the binding mechanism of potential drug candidates that inhibit SARS-CoV-2 structural proteins.

Tika Ram Lamichhane: Conceived and designed the experiments; Analyzed and interpreted the data; Wrote the paper.

Madhav Prasad Ghimire: Contributed reagents, materials, analysis tools or data; Wrote the paper.

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Data will be made available on request.

The authors declare no conflict of interest.

No additional information is available for this paper.

Protein DataBank

The PyMOL Molecular Graphics System, Version 1

Olson

Thermodynamics-Interaction Studies-Solids, Liquids and Gases, In-techOpen

Discovery Studio Visualizer, Version 20.1.0.19295, Dassault Systemes

The authors thank Prof. Binil Aryal for providing the necessary facilities to complete this work and Dr. Bhanu Neupane for a careful reading of this manuscript. The computations were performed on the supercomputer machine of Central Department of Physics, Tribhuvan University (TU) supported by the Alexander von Humboldt Foundation, Germany under the equipment grants program and HERP-TU, Nepal.