H172Y mutation perturbs the S1 pocket and nirmatrelvir binding of SARS-CoV-2 main protease through a nonnative hydrogen bond

Nirmatrelvir is an orally available inhibitor of SARS-CoV-2 main protease (Mpro) and the main ingredient of PAXLOVID, a drug approved by FDA for high-risk COVID-19 patients. Although the prevalent Mpro mutants in the SARS-CoV-2 Variants of Concern (e.g., Omicron) are still susceptible to nirmatrelvir, a rare natural mutation, H172Y, was found to significantly reduce nirmatrelvir’s inhibitory activity. As the selective pressure of antiviral therapy may favor resistance mutations, there is an urgent need to understand the effect of the H172Y mutation on Mpro’s structure, function, and drug resistance. Here we report the molecular dynamics (MD) simulations as well as the measurements of stability, enzyme kinetics of H172Y Mpro, and IC50 value of nirmatrelvir. Simulations showed that mutation disrupts the interactions between the S1 pocket and N terminus of the opposite protomer. Intriguingly, a native hydrogen bond (H-bond) between Phe140 and the N terminus is replaced by a transient H-bond between Phe140 and Tyr172. In the ligand-free simulations, strengthening of this nonnative H-bond is correlated with disruption of the conserved aromatic stacking between Phe140 and His163, leading to a partial collapse of the oxyanion loop. In the nirmatrelvir-bound simulations, the nonnative H-bond is correlated with the loss of an important H-bond between Glu166 and nirmatrelvir’s lactam nitrogen at P1 position. These results are consistent with the newly reported X-ray structures of H172Y Mpro and suggest a mechanism by which the H172Y substitution perturbs the S1 pocket, leading to the decreased structural stability and binding affinity, which in turn explains the drastic reduction in catalytic activity and antiviral susceptibility.


Introduction
The COVID-19 pandemic is still ongoing and remains a major global health threat. At the end of 2021, the U.S. Food and Drug Administration (FDA) issued an Emergency Use Authorization for Pfizer's PAXLOVID to treat mild-to-moderate COVID-19 cases. 1,2 In a recent clinical trial for high-risk non-hospitalized adults with COVID-19, 3 PAXLOVID reduced the risk of progression to severe disease by 89% as compared to placebo. This antiviral drug is a ritonavir-boosted formulation of nirmatrelvir (PF-07321332), an orally available inhibitor of the SARS-CoV-2 main protease (Mpro). Mpro, which is also known as 3CLpro or Nsp5, is a cysteine protease essential to the viral replication process as it cleaves the majority of the polyproteins pp1a and pp1ab. 4 Nirmatrelvir is a reversible covalent peptidomimetic inhibitor, which binds to the active site of Mpro and inhibits its proteolytic activity. 5 Although Mpro is one of the most conserved proteins among coronaviruses, 4 the rapid and constant evolution of the viral genome raises great concern of potential emergence of antiviral resistance. Several biochemical studies, however, showed that the prevalent Mpro mutants in the Variants of Concern or Variants of Interest declared by the World Health Organization (WHO), such as G15S (Lambda), K90R (Beta), and P132H (Omicron), are still susceptible to nirmatrelvir, with IC50 values and catalytic efficiencies similar to the wild type (WT) Mpro. [6][7][8] Nevertheless, biochemical assays of several infrequent natural substitutions, e.g., H164N, H172Y, and Q189K, are associated with reduced activities of nirmatrelvir, among which H172Y caused the largest reduction in the inhibitory activity, with a 233-fold increase in the K i value of nirmatrelvir according to a disclosure by Pfizer. 1 Although H172Y is a rare mutation (found in only a few entries of the database GISAID 9 ), it may become favored in the future under the selection pressure of nirmatrelvir therapy. Thus, understanding the antiviral resistance mechanism is important and urgently needed.
Motivated by the aforementioned need, we investigated the H172Y mutation effect on the Mpro's conformation and binding interactions with nirmatrelvir using all-atom molecular dynamics (MD) simulations. As an experimental structure of H172Y Mpro was unavailable at the start of the study (see later discussion), the simulations were initiated from a mutant model. The simulations showed that the H172Y substitution disrupts the interactions between the S1 pocket and N terminus of the opposite protomer. Importantly, the formation of a nonnative hydrogen bond (H-bond) between Tyr172 and the conserved Phe140 is correlated with a partial collapse of the S1 pocket in the free enzyme and the loss of an important H-bond with nirmatrelvir. These results offer a mechanism for the reduced inhibitory activity of nirmatrelvir and also suggest that the H172Y mutation may decrease the catalytic efficiency and dimer stability of the Mpro. Following the simulation study, we have obtained the enzyme kinetic data of H172Y Mpro, which confirmed the reduced catalytic efficiency. The simulation findings regarding the loss of N-finger interactions and formation of a nonnative hydrogen bond corroborate with a recent bioRxiv paper 10 which reports a free and ligand-bond X-ray structures of H172Y Mpro.

Results and Discussion
Simulations of free H172Y Mpro revealed perturbed S1 pocket and N terminus. To build a structure model of H172Y Mpro for initiating MD simulations, we used the X-ray crystal structure of WT Mpro in complex with nirmatrelvir (PDB id 7vh8, resolution 1.58 Å) 5 as a template in the homology modelling software Modeller. 11 Our previous simulation work 12 showed that His172 is predominantly neutral at physiological pH, and changing to the doubly protonated charged state results in a partial collapse of the S1 pocket. Thus, we rationed that substitution of a neutral histidine with tyrosine would not cause a significant change in the Mpro structure. Indeed, the modeled H172Y Mpro structure is nearly superimposable with the WT, except for a slight displacement of the backbone of Phe140, resulting in a 0.3 Å larger distance between the backbone carbonyl oxygen of Phe140 a b I II III E166 H163 F140 S1* H172Y Figure 1: Overlay of the WT and modeled H172Y Mpro structures. a. Overlay of the modeled structure of H172Y mutant (blue) on the X-ray structure of WT SARS-CoV-2 Mpro dimer (cyan, PDB id 7vh8). 5 I, II, and III represent the three domains that compose the Mpro monomer. A surface rendering is shown for protomer B, and inhibitor is hidden. b. A zoomed-in view of the S1 pocket and the region of the H172Y mutation. Dashed lines represent conserved hydrogen bonds or electrostatic interactions involving Phe140, Glu166, and Ser1 * . and the amino nitrogen of Ser1 * (asterisk indicates the opposite protomer). Starting from this modeled structure, we carried out MD simulations of the ligand-free as well as the nirmatrelvir-bound H172Y Mpros using the Amber20 program. 13 Each simulation lasted 2 µs and was repeated once (Table S1). The protonation states were the same as in our previous study of the WT Mpro. 12 As a control and reference, the ligand-free and nirmatrelvir bound WT Mpros were also simulated using the same settings.
The Mpro's S1 pocket conformation is upheld by several key interactions, one of which is the aromatic stacking between the absolutely conserved Phe140 and His163 (Fig. 1).
This interaction is stable in the present WT simulation and the previous 12 WT simulations started from a different crystal structure (PDB id 6Y2G). 4 The center-of-mass (COM) distance is about 4 Å between the aromatic rings of Phe140 and His163 (Fig. 2a, magenta line). However, while the stacking interaction was stable in simulation run 2 of the H172Y mutant, it got lost after about 1 µs in simulation run 1, with the F140-H163 stack-ing distance increased to above 7 Å (Fig. 2a). The breakage of the aromatic stacking is concurrent with a sudden ∼1-Å increase in the heavy-atom root-mean-square deviation (RMSD) of the oxyanion loop (L1, residues 139-145) which harbors Phe140 (Fig. S3), and a sudden ∼2-Å decrease in the COM distance between L1 and Glu166 sidechain.
The latter is reminiscent of the L1 collapse observed in the simulations of the H172 protonated WT Mpro 12 as well as an X-ray structure of SARS-CoV Mpro determined at pH 6 (PDB id 1uj1). 14 The S1 pocket conformation as well as the stability of the Mpro dimer are also supported by the interactions with the so-called N-finger (residues 1-9) of the opposite protomer. 14 In particular, Ser1 * plays an important role, as its backbone and sidechain maintain H-bond or salt bridge with three absolutely conserved residues in the S1 pocket, His172, Glu166, and Phe140, according to the X-ray structures and the previous 12 as well as the current WT Mpro simulations. We first consider the H-bond between the hydroxyl group of Ser1 * and the imidazole group of His172, which was stable in the WT simulation ( Fig. 2d, magenta line). Interestingly, to compensate for the missing imidazole group, the hydroxyl group of Tyr172 initially interacted with the amine nitrogen of Ser1 * in both runs of H172Y Mpro (Fig. 2d, Fig. S5 and S6). However, the Tyr172-Ser1 * interaction was completely lost after about 0.75 µs in run 1 and 1.8 µs in run 2. The distance between the hydroxyl oxygen and amine nitrogen increased to 5-10 Å in both protomers of run 1 (Fig. 2d) and 14-18 Å in protomer B of run 2 (Fig. S6).
Another prominent N-terminus interaction is with Glu166. Glu166 and the terminal amine form either a H-bond/salt bridge or loose electrostatic interaction in the WT Mpro simulation, with the minimum charge-center distance of 2.5 or 5 Å, respectively (  between the COM of the carboxylate oxygens of Glu166 and the Cα atoms of L1; c. between the backbone amide nitrogen of Phe140 and the hydroxyl oxygen of Tyr172; d. between the hydroxyl oxygen of Tyr172 and the N-terminus nitrogen * ; e. between the nearest carboxylate oxygen of Glu166 and the N-terminus nitrogen * ; f. between the backbone carbonyl oxygen of Phe140 and the N-terminus nitrogen * . The magenta dashed lines represent the corresponding average distances sampled by the simulations of the free WT Mpro. In c. and d., the magenta lines refer to the distances involving the imidazole nitrogen of H172. g. A representative structure of H172Y Mpro taken from the clustering analysis of the last 1-µs trajectory. h. Probability density as a function of the F140-H163 and F140-Y172 distances from the last 1-µs trajectory. Analysis here uses the data of protomer A from the simulation run 1. Analysis of protomer B of simulation run 1 and simulation run 2 is given in Fig. S5 and S6, respectively.

Formation of a nonnative H-bond disrupts the S1 pocket and nirmatrelvir binding.
A curious question is why the H172Y mutation perturbs the S1 pocket conformation and its interactions with the N-finger. In comparing the H172Y and WT trajectories, we noticed that the hydroxyl group of Tyr172 can occasionally accept a H-bond from the backbone amide nitrogen of Phe140 (Fig. 2c, Fig. 2g, Fig S5, and Fig. S6), whereas the analogous H-bond between the imidazole of His172 and the carbonyl of Phe140 is not possible ( Fig. 2c and Fig. 2g). Interestingly, around the same time as the Phe140-His163 stacking got lost in the simulation run 1 of H172Y Mpro, the distance between the hydroxyl oxygen of Phe140 and the amide nitrogen of Phe140 suddenly decreased (Fig. 2c), which resulted in a significant increase of the H-bond occupancy from about 10% to about 50% ( Fig. S7). A representative structure (cluster centroid) obtained from clustering analysis of the second half of the trajectory confirms a perturbed S1 pocket, whereby the Phe140-His163 stacking is abolished and Ser1 * is moved away from the S1 pocket; however, Tyr172 is in a tight H-bond with the backbone of Phe140. These data suggest that the nonnative H-bond between Tyr172 and Phe140 plays a key role in perturbing the conformational dynamics of the S1 pocket and its interactions with the N-finger. The interaction between Phe140 and Tyr172 appears to compete with the Phe140-Ser1 * H-bond and facilitates its breakage.
We further hypothesized that a strong Tyr172-Phe140 H-bond would disrupt the aromatic stacking between Phe140 and His163. To test this hypothesis, we calculated the two-dimensional probability densities for the distances between Phe140 and His163 and between Tyr172 and Phe140. The density map shows a maximum located around the Phe140-His163 and Tyr172-Phe140 distances of 7.5 Å and 3.0 Å, respectively (Fig. 2h), representing a perturbed state in which the Phe140-His163 stacking is disrupted but a stable H-bond between Tyr172-Phe140 is formed. The density map also shows a local density maximum located at the Phe140-His163 and Tyr172-Phe140 distances of 4 Å and 3.1-3.6, respectively (Fig. 2h), representing a state in which the Phe140-His163

stacking interaction is intact and an occasional H-bond is formed between Tyr172 and
Phe140. This analysis demonstrates that the formation of the Tyr172-Phe140 H-bond is correlated with the loss of the aromatic stacking between Phe140 and His163, which may be responsible for the partial collapse of the oxyanion loop.
The mutation induced perturbation of the S1 pocket may explain the 233-fold decrease in the nirmatrelvir activity. 2 To test this hypothesis, we performed two replicas of 2-µs simulation of the WT and H172Y Mpros in complex with nirmatrelvir (Table S1). To facilitate possible conformational changes, the noncovalent binding mode captured by the X-ray structure (PDB id 7vh8) 5 was used to initiate the simulations. Although the aromatic stacking between Phe140 and His163 was intact in these simulations, the interactions between Ser1 * and Glu166 or Phe140 were lost in both simulation runs ( Fig. S8 and S9), and the RMSD of the oxyanion loop was unstable ( Fig. S3c and S3d). These data are consistent with the simulations of the free H172Y Mpro, and suggest that the S1 pocket is perturbed by the mutation.
In both simulation runs of the inhibitor-bound H172Y Mpro, nirmatrelvir remained bound to H172Y Mpro, although in run 2, the RMSD briefly increased by as much as 1 Å in protomer B (Fig. S10d). The distribution of the RMSD of nirmatrelvir relative to the X-ray structure is broader and the peak is slightly shifted to larger value in the H172Y than the WT simulations, which indicates that nirmatrelvir is more flexible and has a small conformational change when complexed with the mutant Mpro (Fig. S11a). Examination of the protein-ligand interactions showed that the change mainly affects the γ-lactam ring in the P1 position, of which the amide nitrogen is in a H-bond with the carboxylate oxygen of Glu166 in the X-ray structure (PDB id 7vh8). 5 This H-bond was stable in the WT simulation, with an occupancy over 60% (Fig. S11b), but it was significantly weakened in the simulation of the H172Y mutant, with an occupancy about 20% (Fig. S11b). By contrast, the H-bond between the lactam nitrogen and the backbone carboxyl oxygen of Phe140 was stabilized in the mutant simulations, with an occupancy of about 70% as compared to 40% in the WT simulation (Fig. S11b). The representative structures of the nirmatrelvir bound WT and H172Y Mpros confirmed the mutation-induced change of the P1 site, i.e., the H-bond between the lactam ring and Glu166 is absent (Fig. 3). Importantly, similar to the ligand-free simulations, they also display the nonnative H-bond between Tyr172 and Phe140 as well as the departure of Ser1 * from the S1 pocket (Fig. 3). This consistency suggests that the perturbation of the S1 pocket by the H172Y mutation is responsible for the change in the P1 site binding, which we speculate may contribute to the decreased affinity for nirmatrelvir.
H172Y mutation reduces Mpro's stability, catalytic activity, and susceptibility to nirmatrelvir. Following the simulation study, we determined the thermal stability and enzyme kinetics of WT and H172Y Mpros as well as the IC 50 values of nirmatrelvir (Fig. 4).
Thermal-shift assays were used to determine the unfolding temperatures (Tm) of the WT Mpro and its H172Y mutant (Fig. 4a). The Tm for the WT was found to be 57.9 • C, whereas that of the H172Y mutant was lower by 4.2 • C (Fig. 4d). This is consistent with the MD results of the MD simulations, which reveal considerable disturbance of the structure.
Kinetic analysis revealed a significant decline in enzymatic activity of H172Y compared to WT Mpro (Fig. 4a). The K m value obtained for the H172Y Mpro was found to be 802.7±273.5 µM, which is 69% larger than the value for the WT Mpro, indicating that the H172Y mutation significantly reduces substrate binding affinity. The k cat /K m value obtained for the WT enzyme was 5355.3±1447.8 M −1 s −1 , while that for H172Y was only 863.3±380.3 M −1 s −1 , i.e. only 16% enzyme activity remains in the FRET assay, compared to the WT. The IC 50 of nirmatrelvir against the H172Y mutant is 344.2±89.0 nM, around 24 times higher than that for the WT protein. In conclusion, even though the H172Y mutant conveys substantial resistance to nirmatrelvir, its catalytic activity is so disturbed that it is unlikely to become a dominant nirmatrelvir resistance mutation.
As we were preparing the manuscript, a bioRxiv paper 10 was published that reports the X-ray structures of the free and inhibitor-bound H172Y Mpros. In the ligand-free X-ray structure (PDB id: 8d4j), 10 the salt bridge between Glu166 and the N-terminus * is lost in one protomer, and the H-bond between Phe140 and the N-terminus * is lost in both protomers (Table S2). These data corroborate the simulation finding of the abolished interactions between Phe140/Glu166 and the N-terminus * (Fig. 2e and f). Note, in the inhibitor bound structure (PDB id: 8d4k), 10 the position of Ser1 is not resolved. As to the interaction between Phe140 and Tyr172, the X-ray structures (PDB ids: 8d4j and 8d4k) 10 display the Phe140:O-Tyr172:OH distance of 3.5/3.6 in the free enzyme and 3.3/3.3 Å in the ligand-bound form (Table S2). These distances are similar to the average distances sampled by the simulation trajectories when the S1 pocket is stable (Fig. 2, Fig. S6, S8, and S9). In this work, 10 the k cat /K m value of H172Y Mpro was reported as 13.9-fold lower than the WT (compared to the 6.2-fold decrease determined by us), and the K i value of nirmatrelvir was reported as >113.7 fold higher than the WT. Taken together, our data offer a mechanism in which the H172Y mutation perturbs the N-finger and S1 pocket, leading to the decreased structural stability and binding affinity of the H172Y Mpro. The latter is responsible for the significant reduction in the catalytic activity and susceptibility to nirmatrelvir. Our finding of the altered P1 site interaction has implication for inhibitor optimization to restore antiviral potency. Although the drastic decline in catalytic activity suggests that H172Y may not become a dominant nirmatrelvir resistance mutation, it is conceivable that double or triple mutations in combination with H172Y may emerge to rescue function while maintaining antiviral resistance (e.g., as demonstrated by a recent work 10 ). On that note, our work demonstrates the synergy between predictive MD simulations and biochemical experiments in the active surveillance of continually emerging SARS-CoV-2 mutations.

Computational methods and protocol
The Modeller software 11 was used to generate an initial structural model of the H172Y mutant of SARS-CoV-2 Mpro, with the X-ray crystal structure of the wild-type (WT) Mpro in complex with nirmatrelvir (PDB id 7vh8) 5 as a template. Next, the WT and H172Y Mpro were prepared for MD simulations using the LEAP utility of Amber, 13 with the termini left free. The protein was represented by the Amber ff14SB force field 15 and water molecules by the TIP3P model. 16 The protonation states were determined using the GPUaccelerated GBNeck2-CpHMD titration 17 with asynchronous pH replica exchange 18 as in the previous work. 12 The octahedral water box was used to solvate the protein, with a distance of at least 11 Å between the protein heavy atoms and the edges of the box. Sodium and chloride ions were added to neutralize the system and create an ionic strength of 150 mM. For the nirmatrelvir-bound Mpro simulations, the reversible bound model in the X-ray structure was used. The ligand parameters were generated using the general Amber force field (GAFF2) 19 and the AM1 BCC method was used to derive the charges 20 All simulations were carried out using Amber20. 13 First, energy minimization with a harmonic restraint of 100 kcal/mol/Å 2 on the protein heavy atoms was performed for 10000 steps using the steepest descent algorithm followed by 10000 steps using the conjugate gradient algorithm. Next, the system was heated from 100 K to 300 K using the same harmonic restraint in the canonical ensemble by 1 ns. Five equilibration stages using harmonic forces of 10, 5, 2, 1, and 0.1 kcal/mol/Å 2 were then performed for 50 ns in the NPT ensemble. The pressure was maintained at 1 atm using the Berendsen barostat with a relaxation time of 0.1 ps, and the temperature was maintained at 300 K using the Langevin thermostat with a collision frequency of 1.0 ps −1 . 13 The particle-mesh Ewald 21 method was used to treat the long-range electrostatics with a grid spacing of 1 Å. A cutoff of 8 Å was used for van der Waals interactions. SHAKE was used to increase the time step to 2 fs. Finally, the production simulations were performed for 2 µs for both the ligand-free and nirmatrelvir-bound WT and H172Y Mpros. Each simulation was repeated once.

Protein production and characterization
Cloning of SARS-CoV-2 Mpro H172Y. The H172Y mutation was inserted by overlap extension-PCR reaction. A pair of special primers, H172Y_forward (ACTGGTGTATATGCCGGGACGGACT; the underlined sequence corresponds to the mutated H172Y codon) and H172Y_reverse (AGTCCGTCCCGGCATATACACCAGT) were designed. The first PCR reaction was performed to generate two splice fragments containing a 5 ′ overhang. The WT Mpro coding gene with BamHI and XhoI sites was amplified from the Mpro construct as described previously, 22 and was used as template. The second PCR joined these two spliced fragments to generate the PCR product encoding the Determination of protein stability of SARS-CoV-2 Mpro WT and H172Y by nano differential scanning fluorimetry (nanoDSF). Thermal-shift assays of SARS-CoV-2 Mpro and its H172Y mutant were carried out using the nanoDSF method as implemented in the Prometheus NT.48 (NanoTemper Technologies). The nanoDSF method is based on the autofluorescence of tryptophan (and tyrosine) residues to monitor protein unfolding. As the temperature increases, the protein will unfold and the hydrophobic residues of the protein get exposed, the ratio of autofluorescence at wavelengths 350 nm and 330 nm will change. The first derivative of 350/330 nm can be used to determine the melting Inner-filter effect corrections were applied for the kinetic measurements according to Liu et al. 22 The fluorescence of the substrate (in RFU) dissolved in 100 µL final volume of reaction buffer at the corresponding concentrations used for the kinetic assay was measured and defined as f (substrate). Afterwards, 1 µL free Edans was added (final concentration: 1 µM) to each well, and the fluorescence reading was taken as f (substrate