key: cord-281528-xy8j5jiv authors: Di Paola, Luisa; Hadi-Alijanvand, Hamid; Song, Xingyu; Hu, Guang; Giuliani, Alessandro title: The Discovery of a Putative Allosteric Site in the SARS-CoV-2 Spike Protein Using an Integrated Structural/Dynamic Approach date: 2020-06-17 journal: J Proteome Res DOI: 10.1021/acs.jproteome.0c00273 sha: doc_id: 281528 cord_uid: xy8j5jiv [Image: see text] SARS-CoV-2 has caused the largest pandemic of the twenty-first century (COVID-19), threatening the life and economy of all countries in the world. The identification of novel therapies and vaccines that can mitigate or control this global health threat is among the most important challenges facing biomedical sciences. To construct a long-term strategy to fight both SARS-CoV-2 and other possible future threats from coronaviruses, it is critical to understand the molecular mechanisms underlying the virus action. The viral entry and associated infectivity stems from the formation of the SARS-CoV-2 spike protein complex with angiotensin-converting enzyme 2 (ACE2). The detection of putative allosteric sites on the viral spike protein molecule can be used to elucidate the molecular pathways that can be targeted with allosteric drugs to weaken the spike-ACE2 interaction and, thus, reduce viral infectivity. In this study, we present the results of the application of different computational methods aimed at detecting allosteric sites on the SARS-CoV-2 spike protein. The adopted tools consisted of the protein contact networks (PCNs), SEPAS (Affinity by Flexibility), and perturbation response scanning (PRS) based on elastic network modes. All of these methods were applied to the ACE2 complex with both the SARS-CoV2 and SARS-CoV spike proteins. All of the adopted analyses converged toward a specific region (allosteric modulation region [AMR]), present in both complexes and predicted to act as an allosteric site modulating the binding of the spike protein with ACE2. Preliminary results on hepcidin (a molecule with strong structural and sequence with AMR) indicated an inhibitory effect on the binding affinity of the spike protein toward the ACE2 protein. The recent pandemic caused by SARS-CoV-2 (COVID-19) has posed the most challenging global health and sanitation concern of the twenty-first century, with 6,663,304 confirmed cases of infected and 392,802 confirmed fatalities (World Health Organization, updated on June 6th, 2020 1 ). Moreover, this pandemic has caused a widespread lockdown, affecting more than 1 in 2 individuals in the world. The virus responsible for COVID-19, SARS-CoV-2, emerged in late 2019, in the region of Wuhan, China. 2 SARS-CoV-2 belongs to the family Coronaviridae, which is characterized by large, enveloped, positive-stranded RNA viruses that can potentially cause gastrointestinal, nervous system, and respiratory distress. 3 This family also includes SARS-CoV and MERS-CoV, 4 which, despite several similarities, present some interesting differences from SARS-CoV2. 4, 5 SARS-CoV-2, MERS-CoV, and SARS-CoV use an envelope protein, termed the spike (S) protein, to gain access to host cells. 4, 6 The S protein displays the highest degree of genetic variability in the virus genome. 6 Additionally, the SARS-CoV-2 S protein exhibits a greater affinity for the ACE2 receptor with respect to analogs in SARS-CoV and MERS-CoV. 7, 8 The SARS-CoV-2 S protein shares a 96.2% sequence similarity with the bat SARS-like RaTG13 S protein, indicating the natural origin of the virus. 5, 9 Due to its crucial role in the very first stage of viral infection, several studies 7,10−14 point to the S protein as an elective target to develop coronaviruses therapies and vaccines (the S protein is also the main target of neutralizing antibodies 15−17 ) . The S protein belongs to the viral fusion class I proteins, 7 has a trimeric, crown-like shape, protrudes from the virus envelope, and targets diverse host cell receptors in different species (see Table 1 in ref 4) . Both SARS and SARS-like (MERS and SARS-CoV-2) coronaviruses target the angiotensin-convertingenzyme 2 (ACE2) present in pneumocytes and enterocytes in the respiratory system and other susceptible cells (e.g., intestinal mucosal, renal tubular epithelial, cerebral neurons, and immune cells). 18 The S-ACE2 interaction occurs at the very first stage of viral entry to initiate viral fusion with the host cell. 19, 20 The hallmark of SARS-CoV-2 is that its S protein is optimized for ACE2 binding. 9 The S protein comprises two subunits, termed S1 and S2, which are cleaved during the first stage of viral infection. The S1 subunit is responsible for the interaction of the S protein with the host receptor (ACE2), while S2 promotes viral fusion. 21 The receptor binding domain (RBD) resides between residues 318 and 513 22 in SARS-CoV-2 and 318 and 510 in SARS-CoV. 13 The development of therapeutics targeting the S-ACE2 complex requires a comprehensive study of this complex interface and its associated dynamics. 23 Thus, investigating the allosteric nature of the S-ACE2 complex will help to understand and quantify the extent by which the S protein binds to ACE2, and therefore the molecular basis of virus infectivity. The identification of allosteric modulation regions (AMRs) allows for the design and testing of allosteric drugs (i.e., drugs targeting a different site from the complex interface) with enhanced bioactivity with respect to orthosteric drugs targeting the protein−protein interaction interface, 24 which are often tightly packed and scarcely accessible. The field of allosteric drugs is a rapidly developing field in drug discovery and may play a crucial role in the therapeutic strategy for SARS-CoV-2. 25−28 In this study, we adopted a multifaceted computational approach to identify promising "druggable" allosteric sites of viral S protein. We first applied the method of protein contact networks (PCNs) to the S-ACE2 complex in SARS-CoV-2 and SARS-CoV. In previous studies, we adopted this method to investigate both protein−protein interactions 29−31 and the allosteric nature of binding. 32, 33 In this work, the application of PCNs clearly highlighted a region in the SARS-CoV and SARS-CoV-2 S protein, which acts as the putative locus for allosteric modulation (AMR). AMR was further investigated in terms of its relevance in the stability of the complex and ability to transmit modulating signals to the binding site by means of independent computational approaches. Specifically, we adopted a monomer-based approach 34 which can predict the affinity of the S protein to its ligand. In this step, the ensembles for the ground state of the structures have been produced using the anisotropic network model (ANM) approach. 35 Moreover, we adopted an unsupervised blind procedure to predict possible interaction sites on S proteins and their affinities for tentative partners using a softness-based prediction of intersubunits affinity (SEPAS). In addition, the dynamical difference in RBD between the S proteins of the two virus strains suggested that they may have different binding and allosteric properties. We also provided biophysical evidence based on the elastic network modeling (ENM) approach, combined with perturbation-response scanning (PRS) 36 that AMRs in both viruses acted as a mediator of intermolecular allostery between the S protein and ACE2. "The three methods converged in allosteric character for residues in AMR in both complexes. This allowed us to state that residues in AMR for SARS-CoV-2 AMR are more susceptible to allosteric drug targeting than for SARS-CoV. A recent study 37 suggested that the many residues in the AMR may be one of the most efficient epitopes for antibody recognition. Preliminary docking studies individuated some molecules that bind to AMR (hepcidin), which were independently shown to exhibit strong sequence similarities with the AMR. 38 The analyzed structures consist of the following: complex SARSspike glycoprotein-human ACE2 complex (stabilized variant, all ACE2-bound particles, PDB code 6CS2, 39 termed the SARS-CoV S/ACE2 complex) and the SARS-CoV-2 analogous spike glycoprotein-human ACE2 complex 40 . Purposed software was used to transform the full structural information in the PDB files into a protein contact network (PCN): the network nodes are the amino acid residues represented by α-carbons. Links between nodes (residues) exist if the mutual distance of the residues (centered on αcarbons) were in the range between 4 and 8 Å, thus including only significant (<8 Å) noncovalent (>4 Å) bonds, while discarding "obliged" contacts due to proximity along the sequence. The adjacency matrix A mathematically represented the PCN in terms of undirected, unweighted network and is defined as The topological role of nodes (residues) addressed the functional role at the corresponding residues, based on the value of network descriptors. The method is widely discussed elsewhere. 41 The basic network descriptor was the node degree k i , defined as To detect allosteric sites and, more generally, the functional regions activating upon binding, we adopted a method based on network spectral clustering. 42, 43 Spectral clustering was based on the Laplacian matrix spectral decomposition, where the Laplacian matrix L is defined as where D is the degree matrix (i.e., a diagonal matrix whose diagonal is the degree vector). The spectral clustering was based on the eigenvalue decomposition of Laplacian L. The sign of the Fiedler vector (the eigenvector corresponding to the second minor eigenvalue) was used for binary partition (each network was parted into two clusters, which in turn could be parted into four, and so on), as outlined in refs 43 and 44. We previously demonstrated that network clusters (group of residues) correspond to functional regions in the protein. 43 Once network nodes (again, protein residues) were parted into the given number of clusters (an independent user-defined parameter for the clustering software), it was possible to define the participation coefficient P where k si is the degree of the ith residue computed in the cluster it belongs (i.e., accounting only for links with nodes pertaining to the same cluster). This network descriptor has been demonstrated to represent the functional role of residues in binding and stability. 33, 41, 43, 45 Residues with a high P value are responsible for the communication between clusters (functional regions) and are thus addressed via the role in allosteric communication. 46 In this study, the participation coefficient maps were projected onto ribbon protein structures (as heat maps), to highlight activating hotspots in the spike-ACE2 complex. Participation coefficient maps are visualized by PyMol (https://pymol.org/2/). We characterized the spike protein/ACE2 interface by means of descriptors stemming from the structural information and the network analysis described above: • the interchain degree k i IC was defined as the node degree, but was only computed over contacts between nodes (residues) belonging to different chains The interface roughness Q/R 29 was defined for each chain participating in an interface, Q was the number of chain residues in the interface and R, the sequence range of chain residues in the interface. • the interface amino acid range 29 IAR = R/N, where N was the total number of residues in the chain • the interface energy matrix, E, defined as Elastic network models (ENMs) provide efficient methods for investigating the intrinsic dynamics and allosteric communication pathways in proteins. We adopted two elastic network models: (1) anisotropic network model (ANM); and (2) Gaussian network model (GNM). In ANM, 47 the interaction potential for a protein of N residues was: In ANM, the 3N × 3N Hessian matrix, H, determined the systems dynamics. where X i , Y i , and Z i represented the Cartesian components of the ith residue, V was the potential energy of the system. We selected r c = 13 Å. ANMs provide information about the amplitudes, as well as the direction of residue fluctuations. The similarity between two ANM modes, u k and v l , evaluated for proteins with two different conformations, can be quantified in terms of the inner product of their eigenvectors: GNM 36 was based on the description of the protein structure as a network of C α connected by springs of uniform force constant γ if they were located within cutoff distance r c . In GNM, the interaction potential for a protein of N residues was where R ij 0 and R ij were the equilibrium and instantaneous distance between residues i and j, and Γ was the N × N Kirchhoff matrix: Thus, the square fluctuations were ( 3 ) The normal modes were extracted by eigenvalue decomposition of the matrix Γ = UΛU T , with U being the orthogonal matrix whose kth column u k was the kth mode eigenvector. Λ was the diagonal matrix of eigenvalues λ k . The perturbation response scanning (PRS) approach was based on the linear response theory and allowed evaluation of residue displacements in response to external forces. 48 Our study performed a PRS analysis based on GNM by constructing the inverse Laplacian matrix, Γ −1 . The N-dimensional vector ΔR of node displacements in response to the application of a perturbation (a N-dimensional force vector F) obeys Hooke's law (F = H • ΔR). The purpose of PRS was to exert a force of a given magnitude on the network, one residue at a time, and F for other residues not being perturbed, equal to zero. Thus, the force exerted on residue i was expressed as 12) and the resulting response was ΔR (i) was an N-dimensional vector that described the deformation of all the residues in response to F (i) . A metric for the response of residue k was the magnitude ⟨∥ΔR k (i) ∥ 2 ⟩ of the kth block of ΔR (i) averaged over multiple F (i) , expressed as the ikth element of the N × N PRS matrix, S PRS . The elements of S PRS referred to the unit (or uniform) perturbing force. The response to unit deformation at each perturbation site was obtained by dividing each row by its diagonal value i k j j j j j j j j j j j y { z z z z z z z z z z z The average effect of the perturbed effector site i on all other residues was computed by averaging over all sensors (receivers) residues j and can be expressed as ⟨(ΔR i ) 2 ⟩ sensor . The effector profile ⟨(ΔR i ) 2 ⟩ effector described the average effect that the local perturbation in the effector site i had on all other residues. The maxima along the effector and sensor profiles corresponded to the functional mobile residues that underwent allosteric structural change. Because the crystallized structure presented one snapshot of the protein dynamic life at the bottom of a potential well, we created a trajectory of protein dynamics in its ground state by setting the RMSD cutoff to 5 Å using the ANM 36,47 approach implemented in Prody. 49 To predict the most possible protein binding patches (PBPs) on spike proteins, we utilized ISPRED to predict the interface residues. 50 The ISPRED-predicted interface residues acted as a seed to create a possible PBP based on the procedure as previously described. 34 The predictions of monomer-based affinities based on the mechanical softness of PBPs were performed by utilizing SEPAS ver. 1. 34 We introduced ALA mutations into the Covid-spike using FoldX suite. 51 We computed the general interface properties of the spike/ ACE2 complex using PISA software, 52 Figure 1B ). The chains are also shown to highlight which chains participate in the two clusters ( Figure 1A ). The interface between the fusion peptide and the ACE2 ectodomain was so strong (in terms of number of contacts between viral and host proteins) that the clustering partition algorithm recognized the fusion+ACE2 region as a single cluster, even though it comprised sequences belonging to different chains. The map of the participation coefficient projected onto the ribbon structure of the SARS-CoV/ACE2 complex ( Figure 1C ) shows an active region (P > 0) in the junction between the fusion peptide and the trimeric bulk phase of the spike protein. Active residues are shown more in greater detail in red in Figure 1D . As previously demonstrated, the participation coefficient describes the attitude of residues toward participating in the intercommunication between clusters; thus, it was a putative score of the involvement of the residues in allosteric communication. Therefore, this region was a good candidate to intervene in the allosteric regulation of complex formation and there was value in investigating allosteric drugs that target this portion of the protein. Figure 2 refers to the analogous analysis for the SARS-CoV2 S/ACE2 complex. A similar active region appeared at the junction of the fusion peptide to the body of the spike protein; however, the active region in the SARS CoV S/ACE2 was more compact than that of SARS-CoV2 S/ACE2. Active regions in both complexes were characterized by two β-sheets and unfolded traits, and could thus be targeted to peptides due to flexibility of the region and absence of a pocket to bind small organic molecules. 53 Table 1 reports the values of the participation coefficient P in residues with P > 0, all of which were located in the chain carrying the fusion peptides of the SARS-CoV S/ACE2 (chain B, in cyan in Figure 1A ) and SARS-CoV2 S/ACE2 (chain C, in magenta in Figure 1A) complexes. The SARS CoV S/ACE2 and SARS-CoV2 S/ACE2 complexes accounted for 21 and 22 active residues, respectively. The average value was similar for the two distributions (0.47 for SARS CoV S/ACE2, 0.48 for SARS-CoV2 S/ACE2). However, the absolute values of higher P residues (>0.6) were higher for COVID19 S/ACE2, indicative of a corresponding larger reactivity of the active region. Characterization of the S-ACE2 Interface Table 2 reports the properties of interface S/ACE2 in the two complexes. The interface area and the absolute value of energy of the SARS CoV S/ACE2 complex were higher than the corresponding area for the SARS-CoV-2 S/ACE2 complex; however, the specific values for the unit area and single residues were higher for SARS-CoV-2 S/ACE2, indicating a more efficient yet less stable interface. This implied an optimized complex formation with the ACE2 of SARS-CoV-2, as observed in ref 9 . This helped to explain the more rapid kinetics of infection, mediated by a faster complex formation with main receptor, ACE2. The RBDs for the two proteins are reported as a wide peptide, which comprised approximately 200 residues for both structures; however, a more limited core region, termed the receptor binding motif (RBM), was claimed to actively participate in the complex formation of the spike protein with its receptor, ACE2. We independently identified the residues in the spike protein for SARS-CoV and SARS-CoV2 using the method of PCNs. The position and the value of the interchain degree are reported in Table 3 for both complexes. The two residues in the two structures with the highest value of interface degree were 488G in SARS-CoV and 502G in SARS-CoV-2, respectively. These two residues could be considered the hotspots of the PPI interface in the S/ACE2 complex; both are glycine residues, likely due to the high steric adaptability of the small glycine residues in the most tightly packed region of the S-ACE2 interface. This confirmed the high density of the interface and its extremely scarce accessibility to small molecules targeting the region, which also suggested the need to interfere with the S/ ACE2 interface through allosteric regulation. In the previous sections, the AMRs of the S/ACE2 complexes have been identified for SARS-CoV-2 and SARS-CoV using the PCNs approach ( Table 1 ). The AMR is accessible for binding, which is at odds with the PBD. One critical property of the protein complexes was the intersubunit affinity. Therefore, besides the S/ACE2 stability, we predicted the affinity between the AMR and incomingdesigned molecule (peptide molecules), which will be designed to affect the function of AMR. This was a crucial step that will provide an independent proofof-concept of both the allosteric features of AMR (by definition, an allosteric site must "sense" the microenvironment) and the druggability of the predicted AMR. We reviewed the need for complete 3D knowledge of the S/ ACE2 complex by means of a recently introduced method to predict the intersubunit affinity of protein complexes using protein binding patch (PBP) on a single subunit by computing the PBP mechanical softness. 34 To alleviate the static presentation of the protein structure in a crystallized state, we created ensembles of protein complexes in their ground state using the ANM approach. In the first step, we predicted the affinity of the SARS-CoV-2 and SARS-CoV S proteins for ACE2. We compared the predicted affinity between the S protein and ACE2 when the S protein was in its monomeric or trimeric state to obtain some hints about the effect of adjacent subunits on the affinity between spike and its ligand (Figure 3 ). Subunits (chains) are denoted as SA, SB, and SC. To manage the potential states of the spike protein, we propose different forms of PBPs on the S protein: • The SA subunit is alone and waiting for its partner (left column, Figure 3 ). We suppose that SA may interact with ACE2 using PBPs similar to those in the ensemble of S/ ACE2 in its dimeric state (the first row, Figure 3 ) or by using PBPs similar to ones SA displays in the trimer spike and ACE2 ensemble (second row, Figure 3 ) • Subunit SA finds ACE2 and binds it directly (right column, Figure 3 ). SA binds to ACE2 via the PBPs that it displays, either when no SB subunit is available (the first row, Figure 3 ) or in association with SB monomers (second row, Figure 3 ). Our computations suggested that SARS-CoV-2 SA (dark blue curves, Figure 3 ) was more sensitive to the presence of the SB subunit than SARS-CoV SA (red curves, Figure 3 ). This result suggested approaches that potentially disrupted the SA-SB association may provide an efficient therapeutic route. We observed a higher affinity between S and ACE2 in SARS-CoV-2 in SARS-CoV. This observation may explain the higher infectivity of SARS-CoV-2. We are interested in predicting the affinities between possible PBPs on spike proteins and tentative protein partners. In this case, we did not know the PBPs on spike proteins because we did not have the resolved 3D structures of all possible protein complexes between the spike protein and other natural or synthesized partners. We predicted the seed interface residues on the spike surface using the ISPRED approach. 50 By expanding the selection radius for defining the PBP region, 34 many possible PBPs were generated for both SARS-CoV-2 and SARS-CoV spike proteins. Using SEPAS, the affinities between the generated PBPs and unknown partners were predicted. The affinities averaged over the whole are presented in Figure 4a and b (for SARS-CoV-2 and SARS-CoV spikes, respectively); the average affinities mapped on the surface of the spike proteins are reported in Supporting Information Figures 1 and 2 . Some of the predicted regions with an affinity < −5 kcal/mol resided in the interface region of the spike subunits. There were some residues in the proposed PBPs that belonged to AMR (indicated as red triangles in Figure 4 ). This meant that AMR could act as a distinct PBP for some partners, suggesting that targeting AMR may be a good choice for drug design. Thus, in considering AMR as distinct PBP, we were able to 1. Predict the affinities between AMR and unknown druglike protein partners, simply by using a 3D structure of AMR on the spike subunit in different ensembles (monomer, case "C" in Figure 5a ) 2. Account for the effect of ACE2 binding on AMR affinity for its possible partners. We assumed the spike monomer in association with ACE2 then predicted the affinity between AMR and its possible partners (dimer, case "DC" in Figure 5a ) 3. Address the effect of the SB spike chains on the affinity of AMR in the SA spike chain for its partners, feeding SEPAS with the trimeric state of spikes (trimer, case "ABC" in Figure 5a ) 4. Analyze the condition of a trimer spike and ACE2 in terms of a holo-complex (S/ACE2 complex, case "ABCD" in Figure 5a ) These results indicated that the AMR of SARS-CoV-2 spike exhibited a higher affinity for generic peptides than the AMR of SARS-CoV spike. Interestingly, the affinity of AMR for their partners was increased upon binding to ACE2, which suggested an interplay between the two distant regions, the ACE2-binding site of spike and AMR. This last result could be considered to be another proof-of-concept of the allosteric character of the AMR. . Predicted affinities of spikes to ACE2 are presented. The affinity of SARS-CoV 2 spike (blue curve) and SARS-CoV spike (red curve) to ACE2 is predicted using SARS-CoV 2/SARS_X_C/B_Y summarized the state of spike proteins in affinity prediction: X denotes number of chain in predictions. X = 1 means S1 protein in monomer form is considered, X = 2 means S1 + ACE2 complex, and X = 4 means whole spike complex + ACE2. Y notes the considered PBP in computations: Y = 2 means the PBP's residues of S1 in S1-ACE2 complex is used to predict the affinity and Y = 4 means the PBP's residues of S1 in whole spike proteins-ACE2 complex is used to predict the affinity between S1 and ACE2. Horizontal axis presents predicted affinity (kcal/mol). Normalized density denotes in Y-axis. Journal of Proteome Research pubs.acs.org/jpr Article To further investigate the effect of the induced structural perturbation in AMR on the binding affinity of the spike protein to ACE2, we introduced two sets of ALA mutations in AMR. In one set, the AMR residues that play a major role in allosteric modulation (P > 0.5 in Table 1 ) were mutated to ALA, and in the other set, all residues of AMR on SARS-CoV-2 SA were mutated to ALA. To simulate an SA perturbation, we used the ensemble of structures generated using the ANM approach after setting the RMSD cutoff to 8 Å for selecting normal mode states. Such global perturbation only distorts the ACE2-binding site of SA to less than 1.5 Å. In Figure 5b , we reported that wild-type SA exhibits a shoulder in its affinity density curve at a high affinity region. Our observations suggested that the structural perturbation of AMR due to the ALA mutation of native residues decreased the affinity of SA for ACE2. This was further confirmation of the allosteric properties of the AMR. In order to explore AMR druggability, we started by recognizing the strict structural similarity of the AMR region with hepcidin, a crucial protein for iron regulation. Moreover, a recent preprint reported a distant (albeit relevant) sequence similarity between hepcidin and the SARS-CoV-2 spike protein. 38 It has been established that there was some proof of the similarity between the hydrophobicity spectra of ligand− receptor pairs, which fostered the observed sequence similarity between hepcidin and the SARS-CoV-2 S-protein. 55, 56 Moreover, hepcidin is strictly involved in interleukin 6 (IL-6) regulation, 57 given that IL-6 is the main actor associated with the "cytokine storm" that is linked to fatalities provoked by SARS-CoV-2 infection. 58 A preliminary docking analysis showed a perfect fit between hepcidin and the AMR region; however, this now must only be considered as a first step with which to prioritize different peptides as allosteric effectors. In addition to the characterization of the S/ACE2 interface, a comparison of RBD dynamics between SARS-CoV and SARS-CoV-2 was also important to understand the functional mechanism of the S/ACE2 interaction. To this aim, the overlapping values between the slowest ANM modes (eigenvectors) were calculated to compare the intrinsic dynamics of RBDs in both S proteins (eq 8 in Materials and Methods). In Figure 6a , the dot products (overlap values) were plotted in the map for the RBDs in SARS-CoV versus SARS-CoV-2 for the first 10 ANM modes. In this map, the upper limit of 1 indicated a perfect overlap (red) correspondent to a perfectly shared dynamics, while 0 indicated no overlap (blue), Figure 4 . Affinity of possible PBPs on the S1 subunit of the spike to tentative proteins are predicted. Possible PBPs are defined for the S1 subunit. The affinity of the PBPs on SARS-CoV2 (a) and SARS-CoV (b) S1 protein for their possible partners is presented. Left vertical axis points to the averaged predicted affinity (kcal/mol) of the proposed PBPs. Right vertical axis denotes the P metric (Table 1 ) of the AMR participat residues. X-axis notes the residue of the PBP's central residue. Red triangles represent the position of AMR residues on the S1 protein sequence. Figure 5 . Affinity of AMR for tentative protein partners and the effect of AMR perturbation on ACE2 affinity are presented. (a) Affinity (kcal/ mol) of AMR to possible partners are predicted when S1 is a monomer (state C), S1 is associated with ACE2 (stade DC), S1 is accompanied by two S2 subunits (state ABC), and when the whole spike binds to ACE2 (state ABCD). (b) Affinity of S1 for ACE2 is predicted when AMR is intact (WT), residues of AMR with P > 0.5 are mutated to ALA (P > 0.5), or all AMR residues are mutated to ALA (P > 0). The affected region of affinity density curve is marked with vertical gray arrow. X-axis notes predicted affinity (kcal/mol) and vertical axis presents normalized density of affinities. indicating that they showed independent dynamics. Although RBDs in the two complexes showed very high structural similarity with an RMSD of 0.68 Å, the overlap between the first 10 ANM modes of two RBDs exhibited some differences. Overlap values with ∼0.6 for the first three modes indicated a weaker dynamical overlap, even at the lowest frequencies, which were those endowed with most important functional meaning. To further understand the relative motions and dynamics of the interface structures, the ANM motions of the first mode for two RBDs are shown in Figure 6b and c. In ANM 1, most of the interfacial residues were very stable (blue), while the peripheral loop exhibited the highest fluctuation (red). The most relevant dynamical difference between the two systems referred to this loop motion. For SARS-CoV RBD, the highest fluctuations point to residues 463D and 465 K, while the SARS-CoV-2 RBD displayed higher flexibilities at 476G to 479P (orange beads). This loop was found to be located at the C terminal of the α1 helix of the RBD, interacting with ACE2 through van der Waals forces, hereinafter referred to as an "interfacial C loop". We suppose that the difference in the dynamics of such an interfacial C loop may be caused by the amino acid mutations in their respective interfaces with ACE2 (i.e., 442Y to 455L, 443L to F, 460F to 473Y, and 479N to 493Q) from SARS-CoV-RBD to SARS-CoV-2-RBD. 59 Thus, we further mapped these interfacial mutations onto the structures (yellow beads). Thus, the lower fluctuations of these mutations in SARS-CoV-2-RBD may account for the higher affinity of such interface. Accordingly, despite the overall structural similarity, RBDs in the two S proteins showed different dynamics, especially for the S-ACE2 interfaces. These differences in dynamics may also indicate that these two S-ACE2 interfaces achieved their activities, including substrate binding and allosteric regulation. We performed perturbation response scanning (PRS) based on GNM to quantify the allosteric properties of AMRs and identify the key residues for which perturbation provokes structural dynamical changes at a long distance. First, the PRS map provided information regarding the sensitivity and effectiveness of a given residue in transmitting signals. The column of the matrix corresponding to a single residue (e.g., for r313V) described how well-connected that residue was within the AMR region and its likely relevance in transmitting allosteric signals to distal parts of the S-ACE2 complex. Here, we adapted this method to probe the allosteric response of the AMR region upon perturbation on residues with P > 0.5 from the PCN prediction (Table 1 ). For the S/ACE2 complex in SARS-CoV, the dynamical changes of the other residues in AMR (red arrows in Figure 7a ) have been caught due to their close distance in space. In addition, all AMR residues have allosteric effects on the interfacial C loop with high fluctuation. Because the effectiveness of the three residues with P > 0.5 was almost identical, we only reported in Figure 7a for the case of the perturbation of 568S that had the highest P value (0.64). It is appreciated that three residues (313V, 529N, and 568S) displayed long-range allosteric effects on some particular regions in ACE2 through the S/ACE interface. The residues relative to the AMR in SARS-CoV-2 could be classified into three types according to their allosteric character. The first class included most residues, including 320Q, 321P, 322T, 323E, 539N, 548T, and 578P. Upon perturbation, these residues exerted a moderate dynamical effect on other AMR residues, and a relatively weak but sensible intramolecular allosteric effect on residues of the S/ACE2 interface, especially on the highly flexible C loop. The perturbation of these AMR residues largely affected the dynamics of ACE2 through stronger long-range allosteric communications. The perturbation of the second type of (531N and 532L) residues provokes a large intermolecular allosteric effect in the RBD, including some peaks at the interfacial C loop (Figure 7c ). In addition, the perturbation weakened the intermolecular communications from AMR to ACE2. The last type only included 580T, which exhibited similar dynamical effects on three regions (AMR, RBD, and ACE2) via long-range communications. Accordingly, our PRS analysis provided important insight into AMR allostery. Among other findings, 3 residues in SARS-CoV AMR and 10 residues in SARS-CoV-2 AMR were predicted to have wide-ranging allosteric effects on the ACE2 protein. This result in conjunction with the statistically significant higher effectiveness of the SARS-CoV-2 AMR perturbation (SI Figure 3) , suggested that its allosteric character was stronger than that of AMR in SARS-CoV. Notably, in both complexes, the interfacial C loop formed a bridge to transmit the signal from AMR to ACE2. The reported RBD-ACE2 interfacial mutations from SARS-CoV to SARS-CoV-2 were close to the interfacial C loop. Therefore, we propose that the different allosteric properties of AMR residues were due to mutations in the S-ACE2 interface. In Figure S3 , the box plot shows that AMR in SARS-CoV-2 had a statistically higher effectiveness than that of SARS-CoV (under the p-value of 7.4 × 10 −5 by Wilcoxon signed-ranked test). The application of PCNs, SEPAS (Affinity by Flexibility) and PRS based on Elastic Network modes to two complexes of the spike protein with its receptor, ACE2, generated highly consistent results. The integrated computational approach disclosed active zones in the SARS-CoV and SARS-CoV2 complexes with the ACE2 ectodomain. These active zones, as demonstrated in previous studies 41, 60, 61 and in the present research by the convergence among different approaches, have an extremely high probability of acting as allosteric sites, which can modulate spike/ACE2 protein binding. These regions are worthy of further investigation from the perspective of drug or vaccine development as indicated by a preliminary docking simulation. From a methodological perspective, our results highlighted the need to identify convergence of different simulation methods arising from independent theoretical premises, in order to obtain a reliable picture of the associated biomolecular systems. The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jproteome.0c00273. Figure S1 : Affinity for the SARS-CoV S/ACE2 complex. Figure S2 : Affinity for the SARS-CoV2 S/ACE2 complex. Figure S3 : Wilcox test for effectiveness (PDF) Coronavirus disease COVID-2019 (In Russ.). Safety and Risk of Pharmacotherapy A Novel Coronavirus from Patients with Pneumonia in China Genome Structure, Replication, and Pathogenesis Mechanisms of Coronavirus Cell Entry Mediated by the Viral Spike Protein The 2019-new Coronavirus Epidemic: Evidence for Virus Evolution Structure, Function, and Evolution of Coronavirus Spike Proteins Characterization of the Receptor-Binding Domain (RBD) of 2019 Novel Coronavirus: Implication for Development of RBD Protein as a Viral Attachment Inhibitor and Vaccine Cryo-EM Structure of the 2019-NCoV Spike in the Prefusion Conformation. Science (80-.) 2020 The Proximal Origin of SARS-CoV-2 MERS-CoV Spike Protein: Targets for Vaccines and Therapeutics Peptide-Protein Interaction Studies of Antimicrobial Peptides Targeting Middle East Respiratory Syndrome Coronavirus Spike Protein: An In Silico Approach Advances in MERS-CoV Vaccines and Therapeutics Based on the Receptor-Binding Domain The Spike Protein of SARS-CoV -A Target for Vaccine and Therapeutic Development Receptor-Binding Domains of Spike Proteins of Emerging or Re-Emerging Viruses as Targets for Development of Antiviral Vaccines Receptor-Binding Domain of SARS-CoV Spike Protein Induces Highly Potent Neutralizing Antibodies: Implication for Developing Subunit Vaccine Preliminary Identification of Potential Vaccine Targets for the COVID-19 SARS-CoV-2) Based on SARS-CoV Immunological Studies A Sequence Homology and Bioinformatic Approach Can Predict Candidate Targets for Immune Responses to SARS-CoV-2 Pathology and Pathogenesis of Severe Acute Respiratory Syndrome Structure, Function, and Evolution of Coronavirus Spike Proteins Angiotensin-Converting Enzyme 2 Is a Functional Receptor for the SARS Coronavirus Functional Assessment of Cell Entry and Receptor Usage for SARS-CoV-2 and Other Lineage B Betacoronaviruses Composition and Divergence of Coronavirus Spike Proteins and Host ACE2 Receptors Predict Potential Intermediate Hosts of SARS-CoV-2 Cryo-EM Structures of MERS-CoV and SARS-CoV Spike Glycoproteins Reveal the Dynamic Receptor Binding Domains Emerging Computational Methods for the Rational Discovery of Allosteric Drugs Exploring and Exploiting Allostery: Models, Evolution, and Drug Targeting From Allosteric Drugs to Allo-Network Drugs: State of the Art and Trends of Design The Different Ways through Which Specificity Works in Orthosteric and Allosteric Drugs Integrated Computational Approaches and Tools for Allosteric Drug Discovery Exploring the Stability of Dimers through Protein Structure Topology Non-Symmetrical Structural Behavior of a Symmetric Protein: The Case of Homo-Trimeric TRAF2 (Tumor Necrosis Factor-Receptor Associated Factor 2) Molecular Features of Interaction between VEGFA and Anti-Angiogenic Drugs Used in Retinal Diseases: A Computational Approach GH32 Family Activity: A Topological Approach through Protein Contact Networks Protein Contact Network Topology: A Natural Language for Allostery Complex Stability Is Encoded in Binding Patch Softness: A Monomer-Based Approach to Predict Inter-Subunit Affinity of Protein Dimers Global Motions Exhibited by Proteins in Micro-to Milliseconds Simulations Concur with Anisotropic Network Model Predictions The Gaussian Network Model. In Normal Mode Analysis Cross-Neutralization of SARS-CoV-2 by a Human Monoclonal SARS-CoV Antibody Distant Sequence Similarity between Hepcidin and the Novel Coronavirus Spike Glycoprotein: A Potential Hint at the Possibility of Local Iron Dysregulation in COVID-19 A Molecular Docking Model of SARS-CoV S1 Protein in Complex with Its Receptor, Human ACE2 Shedding Light on Protein-Ligand Binding by Graph Theory: The Topological Nature of Allostery GH32 Family Activity: A Topological Approach through Protein Contact Networks Modules Identification in Protein Structures: The Topological and Geometrical Solutions Protein Structure: Insights from Graph Theory Proteins as Networks: A Mesoscopic Approach Using Haemoglobin Molecule as Case Study Modules Identification in Protein Structures: The Topological and Geometrical Solutions Anisotropy of Fluctuation Dynamics of Proteins with an Elastic Network Model Manipulation of Conformational Change in Proteins by Single-Residue Perturbations ProDy: Protein Dynamics Inferred from Theory and Experiments ISPRED4: Interaction Sites PREDiction in Protein Structures with a Refining Grammar Model The FoldX Web Server: An Online Force Field Inference of Macromolecular Assemblies from Crystalline State Protein-Protein Interactions as Druggable Targets: Recent Technological Advances Soft Regions of Protein Surface Are Potent for Stable Dimer Formation Mode Matches in Hydrophobic Free Energy Eigenfunctions Predict Peptide-Protein Interactions. Biopolymers Nonlinear Signal Analysis Methods in the Elucidation of Protein Sequence-Structure Relationships Interleukin-6 Induces Hepcidin Expression through STAT3 COVID-19: Consider Cytokine Storm Syndromes and Immunosuppression Structural Basis for the Recognition of the SARS-CoV-2 by Full-Length Human ACE2 Characterization of Protein-Protein Interfaces through a Protein Contact Network Approach The authors declare no competing financial interest. The authors thank the encouraging support and the fruitful discussion from Prof. C. Serangeli. The project was supported by National Natural Science Foundation of China (31872723) and the Priority Academic Program Development (PAPD) of Jiangsu Higher Education Institutions.