Multiple protonation states in ligand-free SARS-CoV-2 main protease revealed by large-scale quantum molecular dynamics simulations.

The main protease (Mpro) in severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) catalyzes the cleavage of polyproteins for viral replication. Here, large-scale quantum molecular dynamics and metadynamics simulations for ligand-free Mpro were performed, where all the atoms were treated quantum-mechanically, focusing on elucidation of the controversial active-site protonation state. The simulations clarified that the interconverting multiple protonation states exist in unliganded Mpro, and the catalytically relevant ion-pair state is more stable than the neutral state, which is consistent with neutron crystallography. The results highlight the importance of the ion-pair state for repurposing or discovering antiviral drugs that target Mpro.

Introduction

In 2019, a novel coronavirus, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), was observed in Wuhan, Hubei province, China, and rapidly spread worldwide, resulting in a worldwide pandemic of coronavirus disease 2019 (COVID-19) [1], [2]. Until January 2022, more than 300 million cases and more than 5 million deaths had been reported. The development of effective vaccines and/or drugs, particularly oral antiviral drugs, is still urged [3], [4], [5].

The main protease (Mpro) or 3CLpro, which is one of the target proteins in SARS-CoV-2 against COVID-19, has an essential role in viral replication via proteolytic processing [6], [7], [8], [9]. Mpro in the dimeric form (Fig. S1) catalyzes the hydrolysis reaction for the cleavage of polyproteins generated after SARS-CoV-2 infection. As SARS-CoV-2 Mpro has 96% identity with SARS-CoV Mpro in the amino-acid sequence [10], [11], therapeutics targeting Mpro are expected to be effective even for mutated SARS-CoV-2. As a result, covalent drugs such as N3 [12], GC373 [13], YH-53 [14], and PF-07321332 [15], and non-covalent drugs such as ebselen [16] and S-217622 [17] have been proposed as effective drug candidates for Mpro.

Recently, docking and/or classical molecular dynamics (MD) simulations for Mpro [18], [19], [20], [21], [22], [23], [24], [25], starting from the crystal structures captured from X-ray crystallography, have been widely conducted to search for effective inhibitors by evaluating the affinity for the active site of Mpro, mainly consisting of the catalytic dyad, i.e., Cys145 and His41 (Fig. 1 a). As hydrogen/proton is not resolved in the usual X-ray crystallography due to the spatial resolution, the protonation states should be presumed in these simulations. The neutral state is the protonation state of the catalytic dyad usually employed in these simulations, except for the limited case [22], [24], [25], to the best of our knowledge. The neutral state is consistent with the findings from previous theoretical studies based on the quantum-mechanics/molecular-mechanics (QM/MM) method, in which Mpro initiates the catalytic cycle via proton transfer from Cys145 to His41 [26], [27], indicating that the ligand-free Mpro remains in the neutral form as the resting state before binding the substrate or inhibitors. In contrast, a recent experimental study based on neutron crystallography at room temperature with the resolution of 2.50 Å, the PDB ID of 7JUN, which can determine the positions of hydrogen/deuterium in principle, revealed that the active site of Mpro forms an unusual ion-pair (IP) state instead of a neutral state in the absence of ligands (Fig. 1 b) [28].

Fig. 1

(a) Close view of the active site in one of the protomers (protomer A) in the crystal structure of the Mpro dimer captured with XFEL at room temperature with the resolution of 1.90 Å (PDB ID: 7CWB). (b) Corresponding active site from neutron crystallography at room temperature with the resolution of 2.50 Å (PDB ID: 7JUN), where the deuteron in the catalytic dyad is highlighted with the light green circle. Note that the numbers indicate the distances in Å. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

In the present study, large-scale quantum MD (QMD) and metadynamics (MetaD) simulations for ligand-free Mpro were performed focusing on the protonation states of the active site. The QMD/MetaD simulations for large-scale systems including more than 10,000 atoms can be achieved using the divide-and-conquer density-functional tight-binding (DC-DFTB) method [29], [30]. Thus far, DC-DFTB-MD/MetaD simulations have been successfully applied to solutions [31], [32], batteries [33], [34], perovskites [35], catalysts [36], [37], and biomolecular systems [38], [39]. Here, the DC-DFTB-MD/MetaD and free energy analysis revealed that the functionally relevant IP state consistent with neutron crystallography [28] pre-exists even without ligands, and that this IP state is the most stable protonation state.

Computational method

The crystal structure of apo-Mpro resolved from serial femtosecond crystallography with X-ray free electron laser (XFEL) at room temperature with the resolution of 1.90 Å, the PDB ID of 7CWB [9], was adopted as the initial structure. The standard protonated states at pH 7, including the neutral catalytic dyad suggested in previous QM/MM studies [26], [27], were employed as the initial protonation states. The missing hydrogen atoms of Mpro were added, and the apo-Mpro dimer was solvated with 23,521 water molecules and 8 Na+ ions for neutralization, using the LEaP module in Amber14 software [40]. As shown in Fig. S2, the total number of atoms was 79,935 with the periodic boundary condition, transcending the previous QM/MM studies [26], [27] by two orders of magnitude in the number of atoms in the QM region.

Classical MD simulations were performed for the pre-equilibration using the Sander module of Amber 14 [40] with the restraint for the non-hydrogen atoms of Mpro and the internal water molecules with a harmonic force constant of 100.0 kcal/mol/Å2. Amber ff99SB-ILDN [41] for standard amino acids and TIP3P [42] for internal water molecules and water solvent were used as the empirical force fields. The NPT ensemble simulations were performed with the restraint described above for 10 ns under a pressure of 1.0 atm and temperature of 300.0 K. The pre-equilibrated structures were used as the initial structures for the subsequent DC-DFTB-MD simulations, where the restraint was completely eliminated.

DC-DFTB-MD simulations were conducted for equilibration using the Dcdftbmd program [30]. The DFTB3 method was employed with the 3OB [43] and 3OBw [44] parameter sets. Grimme’s DFT-D3 dispersion correction [45] was applied with the Becke-Johnson damping scheme. The accuracy of the DFTB parameters for describing the proton transfer between the catalytic dyad was assessed for the small model system, i.e., imidazole and methanethiol (CH3SH) in the gas phase (15 atoms in total), through a comparison with the corresponding DFT and MP2 calculations. The details are described in the Supplementary Material (Figs. S3 and S4).

The criteria of the self-consistent charge (SCC) convergence were 10−9 and 10−6 au for the total energy and Mulliken charge, respectively. The modified Broyden method was used with a mixing parameter of 0.3. The DC method with 3 × 3 × 3 Å3 cubic grids and a buffer radius of 5.5 Å was employed to reduce the computational costs. The resulting number of subsystems was 21,700. The pre-equilibrated structures were thermalized for 30 ps under the NVT ensemble at 298.15 K, controlled by the Andersen thermostat with a time step of 0.5 fs. Note that the proton transfer reaction in the active site was not observed during the unbiased 30 ps DC-DFTB-MD simulation due to the limitation of the timescale (Fig. S5), indicating the necessity of the enhanced sampling method such as MetaD for free energy analysis of the protonation state. The structural analysis such as the probability function was performed on the basis of the last 20 ps trajectories. The other analysis of the distances between the counter ions (Na+) and the catalytic dyad was provided in the Supplementary Material (Fig. S6).

DC-DFTB-MetaD simulations were conducted using a well-tempered scheme [46]. The effective coordination number for the S atom on Cys145 with respect to the proton was selected as a collective variable (CV). The distance between the S and H atoms of Cys145 was used as a variable of coordination number expressed in rational functional form with a cutoff distance of 1.6 Å. The initial height of the Gaussian function for the biasing potential was set to 0.63 kcal/mol. The width of the Gaussian function was set to 0.04, which was estimated from the standard deviation of the CV in unbiased DC-DFTB-MD simulations. The Gaussian biasing potential was deposited at a time interval of 20 fs. The bias factor of the well-tempered MetaD technique was set to 15, corresponding to a bias temperature of 4,200 K. All the biased histograms of MetaD trajectories were merged using the weighted histogram analysis method, which effectively estimates the converged unbiased histogram [47]. Furthermore, the reweighting technique was employed to obtain the free energy profiles with respect to the reaction coordinates different from the CV [47].

Results and discussion

Fig. 2 shows the two-dimensional (2D) free energy surface (FES) obtained from the reweighting method merged from four 70-ps DC-DFTB-MetaD trajectories. The horizontal axis denotes the difference between the Sγ(Cys145) − Hγ(Cys145) distance and the Nε(His41) − Hγ(Cys145) distance, associated with the proton transfer between Cys145 and His41 in protomer A. A negative value corresponds to the neutral state, whereas a positive value is related to the IP state. The vertical axis indicates the Sγ(Cys145) − Nε(His41) distance, which is related to the hydrogen bond (or salt bridge) distance in the neutral (or IP) state. Both the horizontal and vertical axes are given in Å. The 2D FES shown in Fig. 2 is characterized by the most stable state, IP1, and three metastable states, N1, N2, and IP2. The global minimum located at (1.8, 3.6) of the horizontal and vertical axes, corresponding to IP1, is in reasonable agreement with the neutron crystal structure located at (2.0, 3.8) in Fig. 2. The overall conformation of Mpro in the IP1 state is also consistent with the neutron structure (Fig. S7). These results support the experimental finding that the IP state is the most stable in ligand-free Mpro [28].

Fig. 2

2D FES in kcal/mol with respect to the distance difference between Sγ(Cys145) − Hγ(Cys145) and Nε(His41) − Hγ(Cys145) and Sγ(Cys145) − Nε(His41) distance projected with the reweighting method, merged from four 70-ps DC-DFTB-MetaD trajectories. The white cross indicates the neutron crystal structure. Representative snapshots of N1, N2, IP1, and IP2 states are also shown, in which the proton in the catalytic dyad is highlighted with the light green circle. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The metastable IP state, IP2, located at (2.8, 4.5) in Fig. 2, can be accessed from IP1 across the free energy barrier of 4.2 kcal/mol by cleaving the salt bridge between the charged Cys145 and His41. This structure is stabilized by the solvent water molecule bridging the charged catalytic dyad, as shown in Fig. 2. The elongation of the Sγ(Cys145) − Nε(His41) distance and the cleavage of the salt bridge resulted in the flip of the Nε(His41) − Hγ(Cys145) bond from Sγ(Cys145) to the P1′ position, which is occupied by the solvent water molecules in the absence of ligands. This flip motion is relevant to the acylation process in the catalytic cycle of holo-Mpro, where the proton transfer occurs from Nε(His41) to N(P1′) during the cleavage of the C(P1) − N(P1′) bond and the formation of the substrate-enzyme complex [26], [27]. The existence of IP2 in apo-Mpro implies that the functionally relevant structure of His41 pre-exists without the substrate.

The N1 state located at (−0.6, 3.2) in Fig. 2 can be accessed from the IP1 state across the free energy barrier of 11.5 kcal/mol, which is the rate-limiting step in the interconversion between the neutral and IP states in unligated Mpro. Although the N1 state was adopted as the protonation state of the active site in the previous docking and classical MD simulations for Mpro, the relative free energy compared with the most stable IP1 was found to be 6.4 kcal/mol in the present free energy analysis, indicating that the N1 state is not the resting state, at least in the ligand-free situation. The elongation of the hydrogen bond between Cys145 and His41 in the N1 state resulted in the transition to the N2 state with the free energy barrier of 2.8 kcal/mol, in which Cys145 exchanged the hydrogen-bond partner from Nε(His41) to O(His165), as shown in Fig. 2 , via the side-chain flip of Cys145. This hydrogen bond was apparently irrelevant for the catalytic proton transfer from Cys145 to His41 in Mpro. Therefore, the development of an inhibitor that enhances this inactive hydrogen-bond structure via steric and/or electrostatic interactions is essential for prohibiting the enzymatic function of Mpro.

The free energy profile for proton transfer in the catalytic dyad was also calculated using QM/MM simulations for SARS-CoV-2 Mpro with and without ligands [48]. The previous study concluded that the free energy cost of the reaction from the N1 to IP1 states in Mpro without any ligands (2.9 kcal/mol) is lower than that with the peptide substrate or inhibitor because the charged residues in the IP1 state can be stabilized by solvent water molecules in the apo-form [48]. In the present study, the difference in the solvent environment surrounding the catalytic dyad between the N1 and IP1 states was analyzed using radial distribution functions (RDFs). Fig. 3 a and b show the RDFs for the O and H atoms, respectively, of the solvent water molecules around the Sγ atom of Cys145 obtained from the unbiased DC-DFTB-MD simulations for ligand-free SARS-CoV-2 Mpro in both the N1 and IP1 states. In the Sγ(Cys145) − H RDF, the first broad peak was located at 3.00 Å in the N1 state, whereas the first sharp peak was located at 1.75 Å followed by the second peak at 2.25 Å in the IP1 state, indicating that the hydrogen-bond network of water molecules around the catalytic dyad in the IP1 state was considerably more structured (Fig. 3 b). A similar tendency was also observed in the Sγ(Cys145) − O RDF (Fig. 3 a). Therefore, the stability of the IP1 state found in the present simulations can be associated with the enhancement of the hydrogen-bond network by the solvent water molecules.

Fig. 3

RDFs for (a) O and (b) H of the solvent water molecules around the Sγ atom of Cys145, which were obtained from the unbiased DC-DFTB-MD simulations for ligand-free SARS-CoV-2 Mpro in both the N1 and IP1 states.

In addition to the catalytic dyad, the functional role of the conserved internal water molecule, referred to as the catalytic water, Wcat, was investigated in a previous study using X-ray crystallography and kinetic experiments together with classical MD simulations for the wild-type and mutants of MERS-CoV and SARS-CoV Mpro [49]. The catalytic water, Wcat, located in the cavity surrounded by His41, Gln167, and Asp190 in MERS-CoV Mpro (corresponding to His41, His164, and Asp187 in SARS-CoV-2 Mpro, respectively), was suggested to mediate the remote interaction between the catalytic dyad and a partial negative charge cluster composed of Arg-Tyr-Asp, which is essential for catalytic activity [49]. In particular, changing the protonation state of the catalytic dyad from the neutral form to the catalytic ion pair was clarified to induce the approach of Wcat to Met168 via the weakening of the hydrogen-bond interaction between Wcat and Gln167/Asp190 in MERS-CoV Mpro [49].

Fig. 4a and b show the 2D FESs for the hydrogen-bond distances between Wcat and His164/Asp187 in the N1 and IP1 states, respectively, obtained from the corresponding potential of mean force calculated from the unbiased DC-DFTB-MD simulations for ligand-free SARS-CoV-2 Mpro. Here, the horizontal (vertical) axis denotes the distance between the O atom of Wcat and the Nδ atom of His164 (the Oδ atom of Asp187). In the N1 state, the shape of the 2D FES was relatively narrow, indicating that the motion of Wcat was restricted by the stable hydrogen-bond interaction between Wcat and His164/Asp187 (Figs. 4a and c). In contrast, the basin of the 2D FES in the IP1 state was relatively broad, indicating that the hydrogen bond between Wcat and His164/Asp187 in this state was more fragile than that in the N1 state (Figs. 4b and d). The flexibility of the hydrogen bonds for Wcat arising from the transition from the N1 to IP1 states was consistent with the previous finding for MERS-CoV and SARS-CoV Mpro determined from the integration of experimental and theoretical studies described above [49].

Fig. 4

2D FESs in kcal/mol for O(Wcat) − Nδ(His164) and O(Wcat) − Oδ(Aps187) distances obtained from the corresponding potential of mean force calculated from the unbiased DC-DFTB-MD simulations for ligand-free SARS-CoV-2 Mpro in (a) N1 and (b) IP1 states. The representative snapshots of the hydrogen-bond structure around Wcat are also shown in (c) N1 and (d) IP1 states.

Finally, we compare the present results with the previous theoretical and experimental studies focusing on the relationship between the protonation state and binding affinity of ligands in SARS-CoV-2 Mpro. Gumbart and co-workers performed the classical MD simulations of SARS-CoV-2 Mpro in both the apo form and the ligand-bound form with N3 and α-ketoamide for several different protonation states of the catalytic dyad and histidine residues in its close proximity [24]. The previous study clarified that the relevant protonation state for the structural stability of the bound inhibitor depends on the ligands, illustrating the importance of evaluating the binding affinity for the appropriate protonation state. Procacci and co-workers carried out the classical MD simulations of SARS-CoV-2 Mpro in the ligand-bound form with PF-07321332 for the IP and neutral states using two different kinds of force fields [25]. They found that the active-site protonation state and force fields have implications for the inhibitor-bound structure, highlighting the necessity of the QMD simulation in addition to the classical MD simulation based on empirical force fields. Moreover, as in the apo form [28], Kovalevsky and co-workers determined the crystal structure of SARS-CoV-2 Mpro in the ligand-bound form with α-ketoamide inhibitor telaprevir using neutron crystallography [50]. A comparison with the corresponding apo form [28] reveals that the binding of the inhibitor might give rise to the different protonation state in the active site, indicating that the ligand binding is coupled with the alteration of the protonation state in the active site. Considering these previous results and the present finding that the IP state is the most stable in the apo form, in silico drug discovery targeting inhibitors with high binding affinity for the IP state rather than the neutral state could be beneficial to developing potent inhibitors against COVID-19. Furthermore, investigating the protonation-coupled binding mechanism using the QMD simulation could be helpful for gaining atomic-level insight into the reactive inhibition process in SARS-CoV-2 Mpro as one of the future directions for solving the COVID-19 pandemic.

Conclusion

In the present study, large-scale DC-DFTB-MD/MetaD simulations for ligand-free Mpro dimer in aqueous solution and free energy analysis of the active-site protonation states were conducted. The results indicate that the functionally relevant IP state, IP1, was more stable than the conventionally accepted neutral state N1. This behavior is consistent with recent neutron crystallography. In addition, the cleavage of the salt bridge in the IP1 state resulted in a pre-existing structure relevant for proton transfer from His41 to N(P1′) in the catalytic cycle, whereas the breaking of the hydrogen bond in the N1 state yielded an inactive hydrogen-bond formation between Cys145 and His164. The present results highlight the importance of the IP1 state for estimating the binding affinity of covalent/non-covalent inhibitors targeting Mpro.

CRediT authorship contribution statement

Junichi Ono: Conceptualization, Methodology, Investigation, Validation, Data curation, Formal analysis, Writing – original draft, Funding acquisition. Uika Koshimizu: Investigation, Validation. Yoshifumi Fukunishi: Methodology, Validation, Writing – review & editing. Hiromi Nakai: Validation, Resources, Data curation, Writing – review & editing, Supervision, Project administration, Funding acquisition.

Declaration of Competing Interest

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.