{'fraction_non_alphanumeric': 0.029504698592474045, 'fraction_numerical': 0.02285690836718782, 'mean_word_length': 4.890742083635485, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computational Analysis of the SARS-CoV-2 RBD–ACE2-Binding Process Based on MD and the 3D-RISM Theory

The binding process of angiotensin-converting enzyme 2 (ACE2) to the receptor-binding domain (RBD) of the severe acute respiratory syndrome-like coronavirus 2 spike protein was investigated using molecular dynamics simulation and the three-dimensional reference interaction-site model theory. The results suggested that the protein-binding process consists of a protein–protein approaching step, followed by a local structural rearrangement step. In the approaching step, the interprotein interaction energy decreased as the proteins approached each other, whereas the solvation free energy increased. As the proteins approached, the glycan of ACE2 first established a hydrogen bond with the RBD. Thereafter, the number of interprotein hydrogen bonds increased rapidly. The solvation free energy increased because of the desolvation of the protein as it approached its partner. The spatial distribution function of the solvent revealed the presence of hydrogen bonds bridged by water molecules on the RBD–ACE2 interface. Finally, principal component analysis revealed that ACE2 showed a pronounced conformational change, whereas there was no significant change in RBD.


INTRODUCTION
The global spread of the novel severe acute respiratory syndrome-like coronavirus 2 (SARS-CoV-2) has continued to cause serious damage to human life, health, and economy. Many experimental and computational studies have been devoted to the investigation of the mechanism underlying SARS-CoV-2 infection and spread. Coronaviruses initiate cell entry via the binding of a viral spike protein to a receptor protein, angiotensin-converting enzyme 2 (ACE2), on the surface of cells. 1,2 The spike protein is located on the surface of the viral particle and is formed by a trimer of proteins consisting of two subunits, S1 and S2. The domain that binds to ACE2, which is called the receptor-binding domain (RBD), is located at the C-terminus of the S1 subunit. The RBD of the spike protein selectively recognizes ACE2. In turn, the binding of ACE2 to the RBD triggers various processes of invasion and proliferation. Therefore, the binding of spike proteins to ACE2 is an essential initial process of SARS-CoV-2 infection.
Because of its importance, the binding process of ACE2 to RBD has been studied with special attention. Molecular dynamics (MD) simulation is a promising tool to investigate the process at the atomic level, and many studies on the ACE2 and RBD interactions have been reported. 3−14 A comprehensive study of the ACE2−RBD-binding process was conducted by Barros et al. 10 They performed an all-atom MD simulation of the full-length glycosylated ACE2 homodimer including membrane and explicit solvent and revealed the importance of the flexible motion in the head-transmembrane domain linker region and helix mobility in the membrane of ACE2. In addition, it was suggested that specific glycans contribute to the inhibition of RBD binding. MD simulation studies of the full-length glycosylated spike protein were also conducted. 11,12,14 These studies revealed the role of the glycans on open−close or up−down conformational change of the spike protein. Differences in the interaction and structure between SARS-CoV-2 and SARS-CoV have also been intensively studied. 4,5,7,9 Those studies revealed the importance of the mutation of specific residues of RBD. Interactions with other coronavirus spike proteins have also been studied, and specific interactions of SARS-CoV-2 have been identified. 15 These interaction analyses have pointed out the importance of hydrogen bonding between amino acids on the RBD−ACE2binding interface, as well as the importance of the hydration structure changes and hydrogen-bond bridging afforded by water molecules. 5,16 The design of RBD−ACE2-binding inhibitors is also being actively studied. 17−20 Chemical modification of amino acids in the spike protein RBD and the use of nanobodies to inhibit binding have also been examined. 21,22 MD simulations are useful for analyzing the binding affinity of mutants, and new analyses are being rapidly performed with each successive appearance of SARS-CoV-2 variants, including the most recent Omicron variant. 4,23−25 As mentioned above, RBD−ACE2 binding has been investigated in various ways, 5,16 and a large amount of knowledge has been accumulating regarding the protein−protein interaction and the hydration structure after binding; however, studies focusing on the solvent effect on the entire binding process are scarce.
In the present study, the binding process of the SARS-CoV-2 spike protein RBD to ACE2 was investigated using MD and the three-dimensional reference interaction-site model (3D-RISM) theory. The 3D-RISM theory is a statistical mechanics theory that is applicable to the solvation of biomolecules. 26−29 This theory has been applied to various biological and chemical processes of biomolecules in solution, including the RBD−ACE2-binding process. 5,16,30,31 In a previous study, we applied this theory to investigate the role of solvation in the binding process and reported the solvation structure and the thermodynamics of the RBD−ACE2 structure, which was optimized based on the experimentally observed structure. In the present study, we employed MD simulation to examine the structural fluctuation of the proteins during the binding process. In general, the determination of protein-binding pathways is computationally exceedingly expensive and, thus, requires state-of-the-art sampling methods and huge computational resources. For the RBD−ACE2-binding process, D. E. Shaw Research provided several representative pathways of binding that were evaluated using long-time-scale MD simulations and accelerated weighted ensemble techniques performed on the Anton 2 supercomputer. 32 In the present study, we employed one of the trajectories provided by D. E. Shaw Research as a reference for the computational efficiency. According to the original trajectories, we assumed that the RBD takes already "up" conformation and only monomer RBD is considered. Although previous studies have suggested that both interprotomer and intraprotomer interactions affect the conformational changes of spike proteins, 11,12,14 which in turn affect the RBD−ACE2 binding, such interactions are ignored in this study. Therefore, the focus here is on the interaction between the RBD monomer and ACE2. Umbrella sampling was performed with an explicit solvent based on the reference trajectory, and fluctuating protein structures along the pathway were obtained. For the obtained protein structures, the 3D-RISM theory was applied to evaluate the solvation structure and thermodynamics. Based on these results, we clarified the correlation between protein structural changes and solvation thermodynamics during the RBD−ACE2-binding process.

METHOD
In the present study, an energetical analysis of the RBD− ACE2-binding process was conducted using the following procedure. First, the structural sampling of RBD−ACE2 was performed by MD simulation using an umbrella sampling algorithm. The structures of RBD−ACE2 were extracted from the obtained trajectory, and the solvation properties, such as solvation free energy (SFE), partial molar volume (PMV), and solvation structure, were determined using the 3D-RISM theory. We also analyzed the structural changes of RBD− ACE2 via a principal component analysis (PCA). 33−35 The details of the computational conditions are described below.
2.1. Structure Sampling. We started the sampling simulation from the structures proposed by D. E. Shaw  32 The simulation employed the structure of PDB entry 6VW1 as a target system. 36 According to the trajectory data, ASN53, 90, 103, 322, and 546 of ACE2 and ASN343 of RBD were glycosylated. In the present study, only the glycans included in the PDB data were introduced. Four possible binding paths from the free energy surface were proposed, one of which (trajectory number 000713) was referred to as the initial structure for the present calculation. Based on the 89 structures that were contained in the binding path, we performed umbrella sampling with constraints on the distance between the centers of mass of RBD and ACE2. Sampling was performed at 10 ns for each window. Prior to sampling, a 1 ns equilibration simulation was performed. The simulations were conducted at 300 K and 1 bar in the NPT ensemble under the MC barostat and Langevin thermostat. Time step was set to 0.002 ps, and the SHAKE method was applied for hydrogen constrain. Amber ff14SB, GLYCAM 06j, and TIP3P force fields were employed for the proteins, glycans, and water, respectively. 37−39 A cubic water box with the periodic boundary condition was used, with box lengths of about 180 Å. Na + and Cl − ions were added to a final concentration of 0.2 M.

3D-RISM Analysis.
We employed the 3D-RISM theory coupled with the Kovalenko-Hirata closure to evaluate the correlation functions and the SFE. 26−28 The SFE was calculated for 20 snapshots that were extracted every 500 ps from each trajectory. Here, all of the explicit water molecules and ions, with the exception of Zn 2+ , were stripped from the snapshots. The same potential parameters as those used in the MD simulation were employed for all species. The temperature was 300 K, and a 0.2 M NaCl aqueous solution under ambient condition was assumed. The number of grid points in the 3D-RISM-KH calculations was 512, 3 with a spacing of 0.5 Å. All calculations were performed using the in-house 3D-RISM codes. 40,41 3. RESULTS AND DISCUSSION 3.1. Binding Process. In Figure 1a, the distance between the center of mass of ACE2 and RBD was plotted against the umbrella sampling windows. Here, the averaged distances of 20 snapshots for each window were plotted. At the unbind state, window = 1, the interprotein distance was about 85 Å, and the distance decreased as the binding process progressed. At the binding state, window = 89, the interprotein distance reached a value of about 48 Å. At around window = 50, the interprotein distance seemed to have reached a plateau, with a distance around 50 Å. This behavior continued up to about window = 75, and the distance decreased to about 48 Å. Based on the distance change, the binding process can be considered to consist of two steps. The approaching step, in which the distance between the proteins decreases drastically, that is, windows = 1−49, is followed by a moderately changing step, that is, windows = 50−89. In Figure 1b, the snapshots of the protein structure taken from windows = 1, 50, and 89 are shown. At window = 50, the distance between the proteins was sufficiently close that the structure was similar to the bound state at window = 89. Therefore, the process after window = 50 would correspond to the arrangement of the local protein structure. In Figure 1c, the free energy of the system is plotted, which was defined by where E MM and Δμ denote the protein structure energy and the SFE, respectively. The plotted values were relative to the average value of window = 1. The change in the total free energy was rather small compared with its components. The total free energy change at window = 89 from window = 1 was about −10 kcal/mol. Because of computational cost limitations, only 20 snapshots were calculated for each window; however, to evaluate the average value of the free energy with high accuracy, the initial and final states were calculated for 100 structures, and the free energy change was still about −10 kcal/mol. The components of the free energy, namely, the protein structure energy and the SFE, exhibited dramatic changes through the binding process. As the proteins approached each other, the conformational energy decreased, whereas the SFE increased. This was attributed to the penalty of free energy caused by desolvation that occurred as the distance between the proteins decreased. Because the first noticeable change in this process was observed at window = 34, it is likely that partial protein−protein contact occurred at this distance, which stabilized a large protein−protein interaction and caused a desolvation penalty. Subsequently, both the structure energy and the SFE changed gradually.
To examine the changes in interprotein interactions that occur during the binding process, the structural energy and number of hydrogen bonds were plotted (Figure 2). The structure energy of the protein can be split into the   Journal of Chemical Information and Modeling pubs.acs.org/jcim Article particular, significant stabilization was observed around window = 34, which was attributed to the formation of interprotein hydrogen bonds, as shown in Figure 2b. The number of protein−protein hydrogen bonds increased rapidly around window = 34, followed by a further increase after window = 70. As observed in Figure 1a, after window = 50, the distance between the centers of mass changed only slightly, whereas the stabilization energy changed significantly. This may be attributed to the contribution of local structural changes to hydrogen-bond formation. Figure 2c shows a clear correlation between hydrogen-bond formation and interprotein interaction energy. However, it is known that the stabilization energy of one hydrogen bond is about 10 kcal/ mol, and the stabilization of the interprotein interaction (about −400 kcal/mol) was too large for the number of hydrogen bonds (up to about 9) in this case. Therefore, the interprotein interaction might be stabilized by the electrostatic interaction and van der Waals (vdW) interaction, in addition to hydrogen bonds.
To address the local structure of the hydrogen bonds between ACE2 and RBD, lists of the hydrogen-bond pairs with high probability of formation are listed in Tables 1 and S1. At window = 50, the glycan of ACE2 formed hydrogen bonds with RBD residues. In addition, the frequency of hydrogenbond formation between RBD/ASN501 and ACE2/LYS353 and between RBD/ASN487 and ACE2/SER19 or ACE2/ GLN24 was also high. However, the hydrogen-bond fraction detected at window = 50 was less than 20%, even for the largest fraction one, suggesting that a strong interprotein interaction is not yet established at this stage. At window = 89, the hydrogen bond fraction increased, indicating a strong interprotein interaction. The residues forming the interprotein hydrogen bond between ACE2 and RBD at window = 89 in higher fractions are illustrated in Figure 3. ACE2/LYS353 and RBD/GLY502 were located almost at the center of the protein−protein interface, which has the highest fraction of hydrogen bonds. The hydrogen bonds of ACE2/TYR83-RBD/ ASN487 and ACE2/SER19-RBD/ALA475 are located right at the edge of the protein−protein interface. These residues correspond to those in "hotspots" suggested by previous studies. 3 In this snapshot, ACE2/GLN42, RBD/THR446, and RBD/TYR449 appeared to be too far away to form hydrogen bonds, in contrast with their high hydrogen-bond-formation fraction. This suggests that the loop structure including RBD/ THR446 and RBD/TYR449 is highly fluctuating. Most of the residues were located at the N-terminal helix of ACE2, and the loop consisted of ASN487 to TYR505 of RBD; however, interestingly, the glycan connected to ACE2/ASN90 and THR415 of RBD formed a hydrogen bond with high probability. This hydrogen bond showed highest probability at window = 50 and even earlier windows. In fact, the hydrogen bond formed at window = 34, as shown in Figure 2b, is between this glycan connected to ACE2/ASN90 and THR415 of RBD. This suggests that the hydrogen-bond formation between the flexible glycans of ACE2 and RBD amino acids plays an important role in the early process of RBD−ACE2 binding. Although a previous study observed ASN53 contacts in addition to ASN90, ASN53 contacts were weak in the present study. 10 This is considered to be due to the short glycans used in the present model compared with that study. For further investigation of the role of glycans, a more realistic model for glycans needs to be used.
3.2. Solvation. In this subsection, let us consider the solvation structure and its role in protein binding. In Figure 4a, the SFE and its components are plotted against the window steps. The SFE can be written as in the 3D-RISM framework. The summation on the right-hand side of the equation is taken for all the solvent species; thus, the SFE can be split into components corresponding to water and ions in the present case. The total SFE increased with the progression of the binding process, and the components of SFE also increased. At about window = 34, the SFE showed a drastic change, which was in accordance with the structural change of the protein, as mentioned above. The SFE of water showed a larger fluctuation compared with that of ions. It is indicated that the ions may coordinate to the relatively rigid structure of the protein. The increase in SFE was caused by the desolvation associated with protein binding. To investigate desolvation in terms of protein conformational changes, the solvent-accessible surface area (SASA) and the PMV are plotted in Figure 4b,c, respectively. The SASA was calculated using the "surf" protocol in cpptraj of the amber suite. 42 The PMV was calculated from the distribution function evaluated by the 3D-RISM theory. 43,44 The SASA decreased with the progression of the protein-binding process. Conversely, the PMV increased with protein binding. This was probably caused by the decrease in the solvent distribution in the first solvation shell because of the decrease in the SASA. The decrease in the interprotein distance resulted in an increase in interprotein interactions, rather than in protein−solvent interactions. Figure 4d shows the snapshots and vdW surfaces for windows = 1, 50, and 89. No contact between the vdW surfaces of the two proteins was shown in window = 1, whereas windows = 50 and 89 showed ∑ μ To visualize the interprotein interaction and the solvation structure of the protein, the snapshots and solvent distribution for windows = 1, 50, and 89 are shown in Figures 5−7. Figure  5 shows the structure near the RBD−ACE2 interface in the initial state (window = 1), displaying the atom model of the hydrogen-bond donor−acceptor residues listed in Table 1. At window = 1, the distance between ACE2 and RBD was sufficiently large for them to be fully solvated independently. The 3D distribution functions (3D-DFs) of solvent species were also considered. (See also Figure S1.) The water molecules were widely distributed on the protein surface of both ACE2 and RBD. The sodium ion peaks on the RBD side correspond to the coordination to RBD/GLU484 and RBD/ GLU471. For chloride ions, there was a prominent peak at  Journal of Chemical Information and Modeling pubs.acs.org/jcim Article ACE2/LYS31 near the center of the interface. The left peak corresponded to the coordination to ACE2/LYS68, and the RBD showed coordination to RBD/LYS452 and 458, which is not involved in hydrogen bonding. At window = 50 (Figures 6 and S2), the RBD−ACE2 interface was quite close, and hydrogen-bond formation was evident. Interestingly, there was sufficient distribution of water at the interface, and ion coordination to ACE2/GLU35 and ACE2/ASP38 was also maintained. Conversely, there was no distribution of chloride ions to ACE2/LYS31 at this threshold of the isovalue surface plot, indicating that the ion coordination was weakened.
At window = 89 (Figures 7 and S3), the RBD−ACE2 interface was well contacted and interprotein hydrogen bonds had formed. Some water distributions were maintained at the interface, although they were weaker compared with the former window. In particular, there were strong water peaks at ACE2/LYS31 and RBD/GLN493, 498. A broad sodium ion peak was observed at ACE2/GLU35, whereas the peak at ACE2/ASP38 was reduced compared with window = 50. Chloride ion peaks were rarely detected on the interface at this stage.
In Figure 8, the Placement algorithm was applied to investigate the position of explicit water molecules at the RBD−ACE2 interface. 45 Water molecules coordinated to ACE2/LYS31, GLU35, and LYS353 and RBD/TYR489 and GLN498 were observed. The water molecule on the left side was coordinated to three residues at the same time, namely, ACE2/ASP38, LYS353, and RBD/GLN498. Sodium ions were also examined in the same way. The sodium ion corresponding to the peak near ACE2/ASP38 in Figure 7 was observed, and coordination to RBD/TYR449 was also observed.  The cases in which RBD is the hydrogen-bond receptor and ACE2 is the donor and vice versa are listed. The evaluation was performed on snapshots of 100 frames, and the fractions in which hydrogen bonding was observed are shown. The atomic labels are according to the convention of the AMBER force field. 0MB618 denotes the β-Lmannose at the terminal of the glycan connected to ACE2/ASN90. The results of the first PC (PC1) and the second PC (PC2) for ACE2 and PC1 for RBD are shown in Figure 9a−c, respectively. ACE2 underwent a characteristic conformational change upon protein binding, whereas the RBD showed no specific change. PC1 of ACE2 corresponded to the open−close structure fluctuation of the hydrophobic pocket. The eigen vectors on the first and second helices of the C terminus were pointing downward, whereas those of the third and fourth helices were pointing upward, as shown in Figure 9a. As seen in the lower panel of Figure 9a, the PC seemed to change at around window = 50. After window = 50, PC took a negative value, indicating that the hydrophobic pocket opened slightly. The change corresponded to the interprotein hydrogen-bond formation discussed above, suggesting that there is a correlation between the open−close structural flexibility of the hydrophobic pocket of ACE2 and the binding process. The PC2 of ACE2 depicted in Figure 9b shows the sliding motion, namely, the eigen vectors of first and second helices of the C terminus were pointing to the right, whereas those of the third and fourth helices were pointing to the left. Conversely, the PC1 of RBD was rather small. The PC1 of RBD corresponded to the local structural fluctuation in the loop structures. The magnitude of the RBD structural fluctuation remained almost unchanged after complex formation.

SUMMARY
The binding process of ACE2 to the RBD of the spike protein of SARS-CoV2 was investigated using MD and the 3D-RISM theory.
The binding process can be regarded as consisting of two distinct steps, namely, a protein−protein approaching step and a local structural rearrangement step. In the approaching step, the interprotein interaction energy decreased as the proteins approached each other, whereas the SFE increased. As the distance between the proteins decreased to about 55 Å, interprotein hydrogen bonds began to form. At that time, it was revealed that the flexible glycan chain of ACE2 established hydrogen bonds with RBD. The importance of glycans has been highlighted in previous studies. [10][11][12]14 In the present analysis, only the glycans included in the PDB were considered, but more sophisticated models need to be used to examine the function of glycans in greater detail. 13 Thereafter, the number of interprotein hydrogen bonds increased rapidly, and eventually about 6−10 hydrogen bonds were established. The SFE increased because of the desolvation of the protein as it approached its partner. The SASA and PMV analyses suggest that the dehydration of hydrophilic residues on the protein surface contributed significantly to this increase in SFE. The spatial distribution function of the solvent revealed the presence of hydrogen bonds bridged by water molecules on the RBD−ACE2-binding interface. Furthermore, a PCA of the protein structures showed structural changes associated with complex formation. ACE2 exhibited a pronounced conformational change, such as the closing of the hydrophobic pocket, upon binding to RBD, whereas there was no significant change in RBD.
It is expected that the findings reported in this study will contribute to a deeper understanding of viral infection and replication, which is required to overcome the current humanscale crisis caused by coronaviruses.
Additional hydrogen bonding list and details of solvent distribution (PDF).