key: cord-0781843-b2gboumv
authors: Farhadian, Sadegh; Heidari-Soureshjani, Ehsan; Hashemi-Shahraki, Fatemeh; Hasanpour-Dehkordi, Ali; Uversky, Vladimir N.; Shirani, Majid; Shareghi, Behzad; Sadeghi, Mehraban; Pirali, Esmaeil; Hadi-Alijanvand, Saeid
title: Identification of the SARS-CoV-2 surface therapeutic targets and drugs using molecular modeling methods for inhibition the virus entry
date: 2022-01-28
journal: J Mol Struct
DOI: 10.1016/j.molstruc.2022.132488
sha: b26350abd2fd7e029500c189bc3a5037dad08b59
doc_id: 781843
cord_uid: b2gboumv

Although COVID-19 emerged as a major concern to public health around the world, no licensed medication has been found as of yet to efficiently stop the virus spread and treat the infection. The SARS-CoV-2 entry into the host cell is driven by the direct interaction of the S1 domain with the ACE-2 receptor followed by conformational changes in the S2 domain, as a result of which fusion peptide is inserted into the target cell membrane, and the fusion process is mediated by the specific interactions between the heptad repeats 1 and 2 (HR1 and HR2) that form the six-helical bundle. Since blocking this interaction between HRs stops virus fusion and prevents its subsequent replication, the HRs inhibitors can be used as anti-COVID drugs. The initial drug selection is based on existing molecular databases to screen for molecules that may have a therapeutic effect on coronavirus. Based on these premises, we chose two approved drugs to investigate their interactions with the HRs (based on docking methods). To this end, molecular dynamics simulations and molecular docking were carried out to investigate the changes in the structure of the SARS-CoV-2 spike protein. Our results revealed, cefpiramide has the highest affinity to S protein, thereby revealing its potential to become an anti-COVID-19 clinical medicine. Therefore, this study offers new ways to re-use existing drugs to combat SARS-CoV-2 infection.

Coronaviruses (CoVs) are important zoonotic pathogens, which are well-known for their capability to infect many mammalian species [1, 2] . CoVs often have effective strategies for transmission and immune evasion, particularly if epidemic and pandemic occur within dense populations of humans [3, 4] . CoVs are enveloped viruses with RNA genomes [5, 6] . They are classified into four genera, α-, β-, γ-, and δ-coronaviruses [7] . Among them, β-CoVs are responsible for previous Severe Acute Respiratory Syndrome (SARS) and Middle Eastern Respiratory Syndrome (MERS) epidemics [8] that were localized and have similar characteristics. The newly discovered 2019-nCoV (also known as SARS-CoV-2) is a β-CoV [9] is responsible for the newest CoV pandemic, which is mysterious, widespread, and dangerous. It has started in December 2019 in Wuhan, China and spread rapidly in more than two hundred countries, infecting almost 35 million people and causing more than a million deaths.

The entry of all CoVs into the human host cells is mediated by the homotrimeric spike proteins (S proteins or SPs) that protrude from the viral envelope forming a crown-like halo that surrounds the viral particle and that defines the name of these viruses. Sequence of the SARS-CoV-2 S protein contains 1273 residues and includes a signal peptide (residues 1-13), the S1 subunit (residues 14-685 residues), and the SP2 subunit (residues 686-1273). The S1 and S2 subunits, which are responsible for receptor binding and membrane fusion, respectively [10] , are further subdivided into several functional domains and regions. In the S1 subunit, one can find an N-terminal domain (residues 14-305) and a receptor-binding domain (RBD, residues 319-541); whereas the SP2 subunit includes the fusion peptide (FP) (residues 788-806), heptapeptide repeat sequence 1 (HR1, residues 912-984), HR2 (residues 1163-1213), single-pass transmembrane (TM) domain (residues 1213-1237), and a short cytoplasm domain (residues 1237-1273). An important feature of the S2 subunit is the presence of two repeat domains, HR1 and HR2 that include a repetitive heptapeptide: HPPHCPC, where H, P, and C correspond to a hydrophobic, polar and charge residue respectively [11] , and which interact with each other, forming a six-helical bundle (6-HB), causing virus fusion to host cell [12] .

The capability of virus binding to host cell receptors is determined by the spike protein. It is well documented that SARS-CoV-2 binds to the angiotensin-converting enzyme-2 (ACE-2) receptor located on the surface of some human cells [13] . Blocking the spike protein binding to human cell receptors and inhibiting the subsequent virus fusion are among the effective therapeutic strategies for fighting this deadly virus. These are also the antiviral effects of the drugs proposed in this study. Since there are no specific drugs or strategies for the SARS-CoV-2 infection treatment, finding means to efficiently control COVID-19 represents an important scientific task that drives efforts of many scientists around the globe. For example, Canrong Wu and et al. examined the effects of multiple drugs on the intracellular phase of the virus and on the interactions of SARS-CoV-2 proteins with the intracellular proteins and constructed a database of anti-viral drugs for SARS-CoV-2 [14] . In the same line, we investigated the effects of drugs on the extracellular phase of the viral life and on the SARS-CoV-2 S2. Functions of many proteins encoded by the CoV genome and associated with SARS-CoV-2 pathology remain a mystery. We are presenting here a computational analysis of one of these mysterious proteins, spike protein, which in addition to playing a role in the host cell entry might function as a potential modulator of host immunity to delay or attenuate the immune response against the viruses.

In this study, we investigated the structure of the spike protein of SARS-CoV-2 after binding to various drugs (Table S1 ) and identified drugs that potentially interact with spike protein with high-affinity. We further discussed potential candidates for the treatment of SARS-CoV-2 infection and COVID-19, the changes they induce in the spike structure, and strategies to inhibit virus entry. In this study, we investigated the structural changes induced in the fusion core of the S protein by all the drugs listed in Table S1 , and the best drug is suggested for the potential utilization in the treatment of COVID-19. Table S1 2. Materials and Methods

With comprehensive studies, the most efficient compounds were selected as virtual ligands. The approved drugs were selected. Various software was used to investigate the changes in S2 conformation. The chemical structure of selected drugs were downloaded from the chem-spider (http://www.chemspider.com) site. The three-dimensional structure of the 6-HB post-fusion core containing HR1 and HR2 domains in the SARS-CoV-2 S protein S2 subunit (PDB ID 6LXT) was obtained from Protein Data Bank (PDB) [15] .

The molecular docking assays for the Spike protein in the presence of all drugs was carried out using Auto dock 4.2.6 software. The structure of the S2 and all drugs were obtained from PDB (RCSB), chem-spider, and drug bank site, respectively. They were saved in PDB format after optimizing by chimera 1.13.1 software [16] . For determining the best drugs that bind to Spike protein, for molecular dynamic simulation study, the initial docking was performed. All receptor and ligand files were prepared in PDBQT format using Auto-Dock Tools. All possible binding sites were explored through blind docking as there was no assumption regarding the relative affinity of the receptors various binding sites. Blind docking was done with a grid box that covering all dimensions of the receptor. A standard docking protocol with the addition of Gasteiger atomic charges and assignment of default atom-types was implemented. In the molecular docking assays, the best conformation for all drugs were selected. After computing the binding free energy and binding mode between all drugs and spike protein-2 using Auto-Dock 4.2.1 suites, the two lowest docking energy of all drugs were analyzed in MD simulation. The grid box size for the both Conivaptan and Cefpiramide-spike protein-2 complex were set 80 Å × 80 Å × 80 Å. The required hydrogen atoms were attached to the spike protein-2. The grid spacing was 1 Å. For determining the type of interaction, R-studio discovery (version V 16.1) and chimera software were applied. The choice of the best binding mode in each binding site was based on the higher affinity of the obtained structures, and then the molecular dynamic simulation was performed.

The GROMACS 2018.1 package conducted 30-ns MD simulation on the elected complex of protein-drugs obtained from docking studies [17] . Molecular dynamics simulations analysis of spike protein-2 and both conivaptan and cefpiramide were performed. The G43a1 force-field parameters were utilized for spike protein, counter ions, and water. The G43a1force field, with the intermolecular potential denoted as a sum of Lennard-Jones (LJ) force and pairwise Coulomb interaction, was applied to provide the topology and interaction parameters. The prodrug server was employed as an automated interface to allocate the parameters [18] . The two complexes were immersed in a defined triclinic box of the SPC water model with 10 Å padding from the walls [19] . The box of complex (2.4 × 2.4 × 6.82 Å 3 ) was filled with the appropriate amount of water molecules and five Na + for neutralization. The 9054 and 5765 water molecules with the SPC model were added to the protein-conivaptan and protein-cefpiramide boxes, respectively.

The system was minimized energy with a maximum of 100,000 steps by the steepest descent method to remove inappropriate contacts. Succeeding, the position restrained simulation with confining the position of heavy atoms of the system, at 300 K has been performed for 500 ps (NVT condition). Then, another 10000-ps simulation has been performed to equilibrate the system at 1 bar pressure while keeping the temperature of the system constant at 300 K (NPT condition). Also, the leap-frog algorithm with a time step of 1 fs was utilized to support the periodic boundary condition and integration of motion equations. The Berendsen thermostat and the Parrinello-Rahman barostat were applied under NVT and NPT conditions, respectively.

Then were used to analyze the GROMACS simulation results [20] .

In this study, we utilzied in silico methods, such as molecular docking and Molecular Dynamic (MD) simulation to work with SARS-CoV-2 spike protein, which is located on the surface of the virus and is responsible for the two processes (receptor binding fusion to the host cell) assocate with the viral cell entry. It is likely to be cleaved into S1 and S2 subuits after interaction with the host receptor. Molecular docking represents a broadly used approah to determine the binding site of the selected drugs and was successfully used to characterize the interaction between drug candidates and the spike protein of the coronavirus [21] . In this study, we selected two drugs, conivaptan and cefpiramide, from the list of the 80 approved drugs in Table S1 Among the 80 drugs anayzed in our study, cefpiramide was ranked second (after conivaptan) based on the significant binding affinity (-6.04 kcal/mol) to S2 and the presece of several binding modes. Many S2 residues can be involved in cefpiramide-S2 interactions and several pockets capable of interaction with the cefpiramide were found in the HR1 domain.

Based on the consensus binding affinity and meaningful interactions with the significant S2 residues, cefpiramide is being recommended as another possible binding partner of the S2 protein. The chosen bound conformation of cefpiramide with S2 is represented in Figures 3 and   4 , and in Table 1 

The study of the structural properties, flexibility, compaction, and conformational changes of the SARS-CoV-2 S2 subunit in the absence and presence of the approved drugs, conivaptan and cefpiramide, were surveyed using Molecular Dynamic (MD) simulations. These target drugs were selected from the set of 80 approved drugs based on the results of their molecular docking to the SARS-CoV-2 protein S2, and the best-docked complexes were taken for the MD simulations.

In the MD simulations, evaluation of the Root Mean Square Deviations (RMSDs) is the most accurate method for investigation of the average distances between the backbone atoms of proteins and a ligand or other proteins [22] . Also, the RMSD can be used for monitoring the equilibrium processes of the system and the stability of protein structure upon the binding of a small ligand [23] . In general, RMSD in MD simulation represents a means to assess conformational deviation from the initial protein structure [24] . Furthermore, the RMSD could be used to compare the stability of protein structure considering the absence and presence of a When a biomolecule, such as protein, undergoes changes in compaction, these changes can be estimated. To gain the corresponding information on the structural compactness of protein conformation during MD simulation, the protein gyration radius (R G ) was investigated. In general, the radius of gyration of a molecule is the root-mean-square distance of all atoms from its center of gravity [26] . Therefore, radius of gyration provides a means to estimate the correlations between protein dimension, protein compaction, and protein folding [27] . Figure 6 shows that the time course of changes in the R G value of the SARS-CoV-2 S2 protein and its complexes with conivaptan and cefpiramide show a certain downward trend. It is also seen that binding of cefpiramide caused larger changes in compaction than binding of conivaptan.

obviously, the RG changes in S2 complexed with conivaptan and cefpiramide were easily distinguishable from the RG changes in S2 alone. The changes in cefpiramide are larger than conivaptan. This might be because of the differences between the conivaptan and cefpiramide binding to S2, where conivaptan binds to the interior of the HR2 helix and fusion core region, whereas cefpiramide interacts with the HR1 head region. There is a strong interaction between HR1 and HR2 domains within the fusion core region. conivaptan binds to the fusion core region disrupting its structure. cefpiramide binds to HR1 domain and affects the interaction between the two HR domains. 

When a dynamic system fluctuates in some portions, the extent of fluctuations can be computed.

One of the widely used methods to study the protein fluctuations is the roots mean square fluctuation (RMSF), which is a key measure of the flexibility of the residue over the whole simulation time [28] . The greater the amino acid RMSF value, the more flexibility it shows in the process of binding. Therefore, the regions, where conivaptan and cefpiramide have largest impacts on protein residues, were evaluated by RMSF analysis (Figure 7) . This analysis revealed that binding of both conivaptan and cefpiramide caused the increase in the average RMSF values of S2 reflecting the increase in the local flexibility of the protein.

Based on these results, we conclude that the structural dynamics and flexibility of S2 protein increases in the presence of these drugs. Furthermore, in the presence of cefpiramide, S2

protein was more flexible, which could be related to the mode of the drug coordination. The higher backbone RMSF values with higher standard deviation (Figure 6 and Table 3 ) represent higher thermal instability, a conclusion supported by subsequent thermodynamics analysis (see below).

The RMSF values of unliganded SARS-CoV-2 spike protein-2 and SARS-CoV-2 spike protein-2-complex were plotted against residue numbers. 

During the simulation process, the conformational changes in the SARS-CoV-2 spike protein-2 can be further clarified by calculating the secondary structure content in this protein alone and in it complexes with conivaptan and cefpiramide. As shown in Figure 8 and Table 4 , the structure of the SARS-CoV-2 protein S2 was characterized by the presence of 62.6%, 23.4%, 10.12%, and 0% of the α-helical, random coil, β-turn, and β-sheet structure, respectively. Therefore, the dominant type of secondary structure in the SARS-CoV-2 spike protein-2 is α-helix. This analysis also showed that the α-helical content of the S2 protein increased to 74.93% and 71.58%

as a result of protein interaction with conivaptan and Cefpiramide, respectively. The random coil content decreased slightly from 23.01% and 23.08 %, respectively, whereas the β-turn and βsheet contents were stable at around 0%, which indicated that conivaptan and cefpiramide binding influenced the secondary structure of the SARS-CoV-2 spike protein-2. Hence, the drug binding-induced increase in the α-helical structure of S2 was accompanied by some decrease in its random coil and β-turn content. However, these changes in the secondary structure of SARS-CoV-2 spike protein-2 during its binding of conivaptan and cefpiramide were not considerable, which is consistent with the results of R G , RMSD, and H-bond analyses. On the other hand, although the drug-induced changes in the secondary structure content of the protein were rather minor, they do not exclude a possibility that the conformational changes induced by cefpiramide and conivaptan in the S2 protein can cause alterations in the SARS-CoV-2 binding to the human cells. 

Next, changes in the Solvent Accessible Surface Area (SASA) upon binding of drugs to the SARS-CoV-2 spike protein-2 were analyzed, and the corresponding results are presented in Figure 9 . The unbound form of the S2 protein was characterized by the SASA of 66.59 nm 2 , whereas SASA decreased slightly after the protein complexation with drugs. Combined with the R G data, these observations suggested some drug-induced compaction of S2 protein. We speculate that the attraction between the hydrophilic residues of the S2 protein and drugs resulted in the observed SASA changes.

Since hydrophobic interactions and hydrogen bonds (H-bonds) are important for protein structure, the H-bonding of three structures were analyzed, and the corresponding results are reported in Table 5 . The binding of cefpiramide and conivaptan resulted in some decrease in the number of the internal H-bonds in S2 protein ( Table 5 ), indicating that the conformation of SARS-CoV-2 spike protein-2 changed after drugs binding. 

The accessible surface area (ASA) is the surface area of a biomolecule that is accessible to a solvent [29] that can be applied to forecast the secondary structure of protein [30, 31] . According to the results of the analysis of ASA of all the S2 amino acids (see Table S2 ), binding of conivaptan and cefpiramide drugs causes considerable changes in ASA of several residues due to the hydrogen bond with the drugs. Also, computed ASA of some residue indicated that the ASA of conivaptan and cefpiramide become smaller after interaction S2 with conivaptan and cefpiramide. Table S2 3.

Proteins are stable due to cooperative intramolecular forces that maintain the native state.

Transitions to non-native states are prevented by free energy barriers that typically cannot be overcome by the natural conformational fluctuations that occur in free proteins. Interaction of a protein interacts with the ligand can increase these structural fluctuations. In the cases where the magnitude of the fluctuations is sufficiently high, the extra kinetic energy imparted on the molecule becomes great enough for free energy barriers to be overcome, allowing the transitions to non-native conformations. At this point, the protein experience conformational alterations and starts to unfold. Base on the information shown in Table 6 , cefpiramide has strong interaction with spike protein-2 and therefore can be considered as a potential medicine for COVID-19.

Approximately 300 K temperature controlled systems. The temperature of the box was kept constant by the thermostat Nose-Hoover. 

The Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA), developed by Kollman and Case, is a method for precisely determining the extent of ligand-protein interaction [33, 34] .To the calculation of binding free energy of docked system (drugs-spike protein complex) by 

The fusion process between the SARS-CoV-2 and host cells is mediated by the heptad repeats (HRs) located within the S2 subunit of the SARS-CoV-2 spike protein. Altering HR interaction can stop virus fusion and prevent virus replication. Since it has been reported that the HR inhibitors can be used as anti-COVID drugs, we selected two approved drugs to investigate their binding to HRs. The results of this study provide a suitable opportunity for a clinical assessment of these drugs. If the clinical results are consistent with the results of this article, these two drugs may be involved in the treatment of the pandemic severe acute respiratory syndrome coronavirus 2. In this study, molecular docking and molecular dynamics (MD) simulations were carried out to investigate the changes in the SARS-CoV-2 protein S structure. Through these MD simulations and molecular docking studies, we explain how the structure and conformation of the SARS-CoV-2 Spike protein-2 were affected by conivaptan and cefpiramide. The conivaptan binding site was located at the top of S2 crowns at the position of HR2 according to the molecular docking results. This drug is anticipated to bind to the Gln920 and Asp1192 amino acids of the spike protein, which play a key role in ACE-2 attachment. The cefipiramid docking results revealed that the binding site was located on the Ser937, Lys933, Ala930, Ser929, Gln926, and Phe927 on the HR1 domain. Drug binding caused a conformational hindrance.

Furthermore, the RMSD of the S2 protein in the presence of conivaptan and cefpiramide was increased and protein structure was destabilized. In addition to the fact that cefpiramide had larger number of hydrogen bonds with S2, its binding energy to S2 compared to conivaptan was slightly lower. Therefore, cefpiramide is expected to have higher binding affinity to the S2 and 

☒ 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.

Emerging Microbes & Infections

Emerging Microbes & Infections