Abstract
Nirmatrelvir is a severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) main protease (Mpro) inhibitor that exerts its antiviral activity by covalently binding to the catalytic cysteine (Cys145) of Mpro; however, the emergence of drug-resistant variants remains an obstacle to successful antiviral therapy. In this study, we established 23 artificial SARS-CoV-2 Mpro mutants by substituting each active-site residue with alanine in the Mpro–nirmatrelvir complex using a computational approach referred to as a virtual alanine scan. Although the methods were primarily used for non-covalent inhibitor complexes, we conducted a virtual alanine scan for a protein–covalent inhibitor complex. Mutants, in which the structural changes of the main chain and the catalytic dyad were minimal, while the ligand configuration was significantly shifted, were considered to potentially confer drug resistance. The analysis revealed 13 residues: Ser1, His41, Tyr54, Phe140, Leu141, Gly143, Ser144, Cys145, Met165, Glu166, Pro168, Gln189, and Gln192 that are important for the recognition of nirmatrelvir by Mpro, and mutations at these residues may result in drug resistance. The ligand shifts observed in the experimentally reported resistant mutant G143S and the artificially mutant G143A were very similar. These results also indicate that a virtual alanine scan can be applied to covalent inhibitors.
1. Introduction
Coronavirus Disease 2019 (COVID-19), caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), was first reported in China in December 2019 and subsequently caused a global pandemic [1]. SARS-CoV-2 is a positive-sense, single-stranded RNA virus with a lipid bilayer envelope. Its genome encodes four structural proteins (spike, envelope, membrane, and nucleocapsid proteins) and two open reading frames (ORF1a and ORF1b). The structural proteins form the viral particle, which is approximately 100 nm in diameter. ORF1a and ORF1b encode nonstructural proteins (nsp 1–16), from which polyproteins pp1a (containing nsp 1–11) and pp1ab (containing nsp 1–10 and 12–16) are translated [2,3]. These polyproteins are cleaved at 11 distinct sites by the SARS-CoV-2 main protease (Mpro) and papain-like protease to yield mature nonstructural proteins that are required for viral replication [4,5]. SARS-CoV-2 Mpro is a cysteine protease with a homodimeric structure consisting of three domains. The active site is located between domains I and II and contains a catalytic dyad consisting of His41 and Cys145. For the catalytic reaction, His41 deprotonates the thiol side chain of Cys145, which facilitates nucleophilic attack on the peptide substrate [6-9]. The deprotonated Cys145 attacks the carbonyl carbon of the polyprotein backbone to form a tetrahedral thiohemiacetal intermediate, which also contributes to the formation of an oxyanion within the substrate [10-12]. This oxyanion is stabilized by the oxyanion hole, which consists of hydrogen bond-donating amide groups from His41, Gly143, and Ser144. Cleavage of the thiohemiacetal into a thioester breaks the peptide bond, thereby releasing the C-terminal portion of the polyprotein substrate [13]. Because of the indispensable role of SARS-CoV-2 Mpro in viral replication, the inhibition of its function suppresses the proliferation of SARS-CoV-2. SARS-CoV-2 Mpro cleaves the peptide bond on the C-terminal side of glutamine residues in its substrates [14-18]. This substrate specificity is not observed in any known human protease [19], which suggests that Mpro is a target for the selective anti-SARS-CoV-2 drug. Thus, SARS-CoV-2 Mpro inhibitors have been proposed and marketed as therapeutic agents for the treatment of COVID-19.
Paxlovid (nirmatrelvir/ritonavir), which was developed by Pfizer and launched in February 2022, is an orally administered SARS-CoV-2 Mpro inhibitor [20]. It is a covalent inhibitor that forms a reversible covalent thioimidate ester with the catalytic dyad cysteine residue (Cys145) of Mpro through its nitrile group [20]. The chemical structure of nirmatrelvir is shown in Fig. 1a, and its position in the Mpro active site is shown in Fig. 1b. Nirmatrelvir is a tripeptidomimetic compound with a nitrile group at its C-terminus and a trifluoroacetamide moiety at its N-terminus. It consists of four groups: P1 is β-(S-2-oxopyrrolidin-3-yl)-methyl side chain, P2 is dimethylcyclopropylproline residue, P3 is tert-butylglycine residue, and P4 is trifluoroacetamide moiety. The P1 group mimics the glutamine residue of substrates and acts as a recognition motif for the S1 pocket, which consists of Phe140, Leu141, Asn142, Glu166, and His172. The P2 group fits deeply into the S2 pocket, which is composed of His41, Met49, Met165, His172, and Gln189. The P4 moiety interacts with the S3 and S4 subpockets and contributes to the favorable pharmacokinetic profile [14,21-24]. Ritonavir has a pharmacokinetic boosting role by irreversibly inhibiting CYP3A4, which is the metabolic enzyme responsible for degrading nirmatrelvir, thus prolonging its plasma half-life and enhancing its antiviral activity [25-27].
In this study, we performed molecular dynamics (MD) simulations on the complex between SARS-CoV-2 Mpro and nirmatrelvir to elucidate its three-dimensional (3D) conformation. The SARS-CoV-2 genome continuously mutates during replication, and even a single amino acid substitution in Mpro can confer drug resistance [28]. To identify the amino acid residues involved in drug resistance, we analyzed the 3D structures of artificial mutants, in which each active site residue of SARS-CoV-2 Mpro was substituted with alanine using the Mpro–nirmatrelvir complex. This is known as a virtual alanine scan, in which each amino acid residue of a protein is individually replaced with alanine to determine the functional or structural significance of the original residue. A mutation that causes the protein to lose its ability to recognize the inhibitor while retaining its enzymatic function may indicate drug resistance. For the calculations, the covalent bond between nirmatrelvir and Mpro was removed from the input files for the MD simulations. Thus, we hypothesized that nirmatrelvir is first recognized by Mpro and subsequently forms a covalent bond. Previously, we conducted a virtual alanine scan on SARS-CoV-2 Mpro complexed with the non-covalent inhibitors, NCL-00024905 and ensitrelvir [29,30]. Ensitrelvir, which is a non-covalent inhibitor marketed by Shionogi, was shown in our virtual alanine scan to be influenced by several mutations experimentally associated with resistance. In the present study, we performed a virtual alanine scan for the covalent inhibitor and identified residues that may confer nirmatrelvir resistance.
2. Materials and Methods
The complex structure of SARS-CoV-2 Mpro and nirmatrelvir was generated based on the crystal structure available in the Protein Data Bank (PDB ID: 7RFS) [20]. The crystal structure of Mpro is a homodimer, with nirmatrelvir positioned as a ligand in the active site of both subunits. Although nirmatrelvir was covalently bound to Mpro to generate the initial wild-type structure of the virtual alanine scan, a pre-covalent binding conformation of the SARS-CoV-2 Mpro–nirmatrelvir complex was used, in which the covalent bond between Mpro and nirmatrelvir was eliminated. Therefore, the atomic coordinates of nirmatrelvir were extracted from the crystal structure, and the only nitrile moiety in nirmatrelvir undergoing nucleophilic attack by the Cys145 in Mpro was optimized using the B3LYP density functional method with 6-31G(d) basis set. Then, nirmatrelvir was redocked in the Mpro active site. The protonation states of the protein at pH 7.4 were calculated based on the predicted pKa values using the H++ web server [31]. Sodium ions (Na⁺) were added to neutralize the net charge of the system. Mutant Mpro–nirmatrelvir complexes were generated by introducing amino acid substitutions into the wild-type structure following MD simulations. The target residues for the mutations were the amino acids located within 4 Å of the ligand at the active site. Each was individually substituted with alanine to construct an artificial mutant. Active site residues, which were analyzed by a virtual alanine scan, are shown in Fig. 2, whereas the generated mutants are listed in Table 1. Structural optimization and MD simulations were conducted using the AMBER20 software package [32]. The AMBER ff14SB force field [33] was used for the proteins, and GAFF2 [34] was applied to the small-molecule inhibitors. The TIP3P water model was used as the solvent. Energy minimization was done in two steps: 1,000 steps for water molecules and counterions, followed by 2,500 steps for the entire system. The system was heated from 0 K to 300 K over 20 ps using heating MD simulations. Equilibration was conducted at 300 K, and the simulation times were 1,000 ns for the wild-type complex and 100 ns for each mutant. The simulation time step was set to 2 fs. Snapshots were saved every 2 ps throughout the simulations. The SHAKE algorithm was applied to constrain bonds involving hydrogen atoms [35].
Table 1. Artificial mutants established by a virtual alanine scan
| S1A |
H41A |
M49A |
Y54A |
F140A |
L141A |
| N142A |
G143A |
S144A |
C145A |
H163A |
H164A |
| M165A |
E166A |
L167A |
P168A |
H172A |
V186A |
| D187A |
R188A |
Q189A |
T190A |
Q192A |
|
Based on the results of MD simulations, we identified amino acid residues potentially involved in drug resistance. Mutants, in which the main chain and catalytic dyad structure were largely unshifted compared with the wild type, while the ligand position exhibited significant displacement, were considered to retain the enzymatic function of Mpro, while being capable of excluding nirmatrelvir. For the wild-type and mutant complexes, the structure obtained after the heated MD simulation served as a reference structure, and the root-mean-square deviation (RMSD) values of the main chain, ligand, and catalytic dyad were calculated. The RMSDs for the main chain served as an indicator of the overall structural stability of Mpro, whereas the RMSDs of the ligands were used to assess whether nirmatrelvir remained properly positioned within Mpro. The RMSDs of the catalytic dyad were used to estimate the potential disruption of enzymatic activity. High RMSD values for the main chain and the catalytic dyad suggest that Mpro may not maintain its proper structure and could have lost enzymatic function. Average RMSD values were calculated using the final 10 ns of each MD trajectory.
To validate the method and determine the threshold value for identifying the drug-resistant mutants, a comparative analysis of G143S, a mutation that has been experimentally reported to confer resistance to nirmatrelvir [36], and G143A was performed. Similar structural features and average RMSD values for G143S and G143A, particularly the RMSDs of the ligands, were similar, which suggests that G143A may also confer drug resistance. The G143S mutant structure was generated under the same conditions as the other mutants, by introducing the mutation into the wild-type Mpro–nirmatrelvir complex and conducting MD simulations. Next, we evaluated the similarity between the MD-derived G143S structure and the corresponding experimentally determined crystal structure (PDB ID: 8DZ9) [36], and compared the computational structures of G143S and G143A to investigate whether alanine scan can be useful for predicting the structural features of proteins with mutations other than Ala.
For further analyses of the MD trajectories, the frequency of hydrogen bonds between the protein and ligand in the resulting complex structures was calculated. Hydrogen bonds were defined using a distance cutoff of 3.5 Å and an angle cutoff of 120°. Analyses of the hydrogen bond frequencies were performed for the final 10 ns of the MD trajectories. The calculations for the RMSD and hydrogen bond frequencies were done using the cpptraj module of AmberTools [32]. Binding free energy calculations for the wild-type and mutant Mpro–nirmatrelvir complexes were performed using the MM/GBSA method [37]. The MM/GBSA calculations were performed for the final 10 ns of each MD trajectory. Only the protein–ligand complex system was explicitly simulated. The trajectories for ligand-free Mpro and nirmatrelvir were established by stripping the ligand or protein, respectively, from the complex trajectory. A modified Generalized Born (GB) model (igb=5) was used as the implicit solvent model [38], with the salt concentration set to 0.15 M. MM/GBSA calculations were performed using the MMPBSA.py script from AmberTools [39].
3. Results and Discussion
3.1 RMSD
The 3D structures of the wild-type SARS-CoV-2 Mpro complex before and after the MD simulations are presented in Fig. 3. The RMSDs for the main chain, ligand, and catalytic dyad throughout the MD simulations of the wild-type Mpro complex as shown in Fig. 4. After a 1000-ns MD simulation, the wild-type SARS-CoV-2 Mpro and nirmatrelvir complex retained a structure that was very similar to its initial conformation. The RMSD values for the main chain and catalytic dyad converged, which indicated structural stability. The RMSD values of the ligand exhibited convergence, and the configuration did not significantly shift during the simulation, indicating that the system had reached an equilibrium state. This supports the conclusion that the wild-type Mpro complex represents a plausible structure, in which nirmatrelvir is recognized by Mpro and inhibits its enzymatic activity. In the present study, the covalent bond between Mpro and nirmatrelvir was virtually cleaved to enable MD simulations of the pre-covalent binding state. The ligand RMSD indicates that the ligand remained stable during the simulation, and it consistently stayed within the active site, even in the absence of a covalent bond. This suggests that nirmatrelvir is recognized by Mpro before the formation of a covalent bond. Table 2 lists the average RMSD values for the main chain, ligands (chains A and B), and the catalytic dyads (chains A and B), calculated from the final 10 ns of the MD trajectories for each mutant Mpro complex. The plots for the ligand and the catalytic dyad RMSDs for each mutant are shown in Figs. S1 and S2, respectively. A small average RMSD for the main chain and catalytic dyad indicates that Mpro maintains its overall structure and retains its enzymatic activity. In contrast, a large ligand RMSD suggests that nirmatrelvir does not remain within the active site of the mutant Mpro. Based on the results of the 100-ns MD simulations, all mutants had average RMSD values of less than 2 Å for the main chain, indicating that the overall structure of Mpro remained. However, in the N142A, L167A, V186A, and T190A mutants, the average RMSDs for the catalytic dyad in chain A exceeded 2 Å. In H163A and H164A, the dyad RMSDs in chain B were greater than 2 Å, whereas in D187A and R188A, the catalytic dyads in both chains exhibited average RMSD values greater than 2 Å. These results suggest that the structures of the catalytic dyad are altered in the mutants compared with the wild-type, potentially leading to reduced or lost enzymatic activity. Mutants N142A, L167A, V186A, D187A, R188A, and T190A exhibited large fluctuations in catalytic dyad RMSDs during MD simulation (Fig. S1). This further supports the possibility that these mutations disrupt the enzymatic function of Mpro. Drug resistance typically occurs when the protein retains its function, but fails to recognize the inhibitor. Therefore, mutations that preserve the overall protein structure and the catalytic dyad structure, but exhibit high ligand RMSD values, indicating that the ligand is no longer in the binding pocket, are considered candidates for drug resistance. Thus, the mutants with low average RMSD values for the main chain and catalytic dyad, and high average RMSD values for the ligand, are associated with resistance to nirmatrelvir.
Table 2. Average RMSDs/Å for the last 10 ns trajectories of the artificial mutants
|
Main Chain
|
Chain A
|
Chain B
|
| Ligand
|
Catalytic dyad
|
Ligand
|
Catalytic dyad
|
| S1A |
1.518
|
1.855
|
1.247
|
1.740
|
1.304
|
| H41A |
1.386
|
1.759
|
1.478
|
1.772
|
1.435
|
| M49A |
1.629 |
1.917
|
1.130
|
1.265 |
1.154
|
| Y54A |
1.373
|
1.918
|
1.263
|
1.671
|
0.937
|
| F140A |
1.290
|
1.640
|
1.014
|
1.365 |
1.160
|
| L141A |
1.369
|
1.826
|
1.114
|
1.599
|
1.107
|
| N142A |
1.257
|
1.530 |
2.721 |
1.221 |
0.966
|
| G143A |
1.401
|
1.374 |
1.154
|
1.552
|
0.907
|
| S144A |
1.272
|
1.885
|
0.990
|
1.566
|
1.195
|
| C145A |
1.339
|
1.607
|
1.027
|
1.228 |
0.738
|
| H163A |
1.382
|
2.743
|
1.204
|
4.273
|
2.077 |
| H164A |
1.506
|
1.477 |
1.628 |
1.679
|
2.058 |
| M165A |
1.341
|
2.005
|
0.961
|
1.481 |
1.207
|
| E166A |
1.551
|
2.211
|
1.221
|
1.295 |
1.138
|
| L167A |
1.755 |
1.750
|
2.380 |
1.601
|
1.124
|
| P168A |
1.537
|
1.642
|
1.088
|
1.912
|
1.256
|
| H172A |
1.337
|
1.477 |
0.861
|
1.365 |
1.130
|
| V186A |
1.323
|
1.925
|
2.396 |
1.294 |
1.504
|
| D187A |
1.604 |
1.679
|
2.472 |
1.640
|
2.117 |
| R188A |
1.646 |
1.842
|
2.866 |
1.381 |
2.284 |
| Q189A |
1.325
|
1.686
|
1.193
|
1.741
|
1.207
|
| T190A |
1.315
|
1.756
|
2.554 |
1.624
|
1.199
|
| Q192A |
1.293
|
1.823
|
1.130
|
1.409 |
1.271
|
Underline: average RMSDs/Å of main chain and catalytic dyad < 1.552
Bold: average RMSDs/Å ≥ 1.552
A comparison between the G143S crystal structure, the G143S computational structure, and the G143A computational structure is shown in Fig. 5. The structures and ligand conformations for G143S (crystal structure), G143S (computational structure), and G143A (computational structure) were very similar (Fig. 5 and Table 3). The RMSD between the G143S (crystal) and G143S (computational) ligands was 1.190 Å, whereas the RMSD between G143S (computational) and G143A (computational) was 0.490 Å, indicating that ligand poses for these structures were similar. Thus, the experimental crystal structure of the mutant G143S can be reproduced using the MD simulation with our settings. The importance of G143 in ligand recognition may be assessed by the MD simulation of the G143A mutant generated by the virtual alanine scan protocol used in the present study. Based on these results, we defined an RMSD value of 1.552 Å, the average ligand RMSD of G143A, as the threshold to define “high” RMSD values or “significant” configuration shifts. Mutants with average RMSD values for the main chain and catalytic dyad below 1.552 Å, and ligand average RMSD values ≥1.552 Å, were considered drug resistance candidates. Because the catalytic dyad is comparable in size to the ligand, the same RMSD value (1.552 Å) was used as the threshold. For the main chain RMSDs, the values obtained from the simulations of the wild type (Fig. 4) were not markedly different between the main chain (Fig. 4a) and the ligands (Fig. 4b). Therefore, we judged that using the same threshold was acceptable. A more detailed and comprehensive examination of the RMSD threshold is currently underway. Based on this criterion, 13 mutants were identified (S1A, H41A, Y54A, F140A, L141A, G143A, S144A, C145A, M165A, E166A, P168A, Q189A, and Q192A), which were considered likely to exhibit nirmatrelvir resistance. Although the C145A mutant is presumably not drug resistant because Cys145 is the catalytic residue, it was included among the candidates.
Table 3. Average RMSDs/Å for the final 10 ns trajectories of MD simulations in mutant G143S and G143A
|
Main Chain
|
Chain A
|
Chain B
|
| Ligand
|
Catalytic dyad
|
Ligand
|
Catalytic dyad
|
| G143S |
1.410 |
1.503 |
1.094 |
1.565 |
1.321 |
| G143A |
1.401 |
1.374 |
1.154 |
1.552 |
0.907 |
3.2 Hydrogen Bond Frequency Analysis
The hydrogen bond frequencies between the ligand and surrounding residues for the wild-type and 13 candidate nirmatrelvir-resistant mutants are listed in Table S1. Fig. 6 shows the residues that form hydrogen bonds with the ligand in the wild-type Mpro complex. During the final 10 ns of MD simulations, the wild-type ligands formed four hydrogen bonds with a frequency >80% in chains A and B. In contrast, for most of the 13 mutants identified as drug resistance candidates, at least one hydrogen bond between the ligands and Mpro exhibited a significant decrease in formation frequency or was completely lost. Of these interactions, hydrogen bonding with Glu166 was previously shown to be important for the binding affinity between nirmatrelvir and Mpro [14]. In the wild-type, Glu166 formed three hydrogen bonds with the ligand. Of these, the frequency for the hydrogen bond between the ligand 2-oxopyrrolidine moiety and the Glu166 side chain showed a significant reduction in several mutants compared with the other two Glu166–ligand hydrogen bonds. For mutants S1A, G143A, S144A, C145A, M165A, and E166A, the hydrogen bonds between Glu166@OE1 and the ligands@NH9 were either lost or decreased in chains A and B. For mutants H41A, Y54A, P168A, and Q192A, the hydrogen bond frequencies between Glu166@OE1 and the ligand@NH9 were significantly decreased in one of the two chains.
The 3D structures obtained from the MD simulations of several mutants are shown in Fig. 7. The interaction between Ser1 and Glu166 is involved in Mpro dimerization [24], and contributes to the formation of the S1 pocket. In the S1A mutant, the structure around the mutation sites is shown in Fig. 7a. In the S1A mutant, the decreased frequency of the hydrogen bond between the Glu166 side chain and the ligand results in a change in ligand configuration. Because structural changes were observed around the mutation site, the conformational change of the S1 pocket may affect the ligand position. In the E166A mutant (Fig. 7b), the side chain of Glu166 was lost due to the mutation, which resulted in the disappearance of the hydrogen bonds in the Glu166 side chain–ligand interaction (Table S1). Moreover, Cys145 and His164 exhibited a decreased frequency of hydrogen bond formation with the ligand, which likely affected the ligand configuration. The 3D structures of G143A, S144A, and C145A are shown in Fig. 7c–e. Previous studies indicate that alanine substitutions of residues within the loop formed by Gly143, Ser144, and Cys145 alter the hydrogen bond network, thereby changing the shape of the binding pocket [30]. In the present study, a marked decrease in the Glu166–ligand hydrogen bond frequency was observed in the G143A, S144A, and C145A mutants. The alanine substitutions affected the surrounding interaction network and altered the shape of the binding pocket. Gly143 plays an important role in stabilizing the oxyanion hole [13], with its backbone amide group forming a hydrogen bond with Glu166. Substitution with Ala143 may destabilize the oxyanion hole, leading to structural changes in the active site (Fig. 7c). For the S144A mutant (Fig. 7d), the hydrogen bond frequency between His164 and the ligand was also significantly decreased. Although the main chain and the catalytic dyad in chain A of the S144A mutant showed low average RMSD values (1.272 Å and 0.990 Å, respectively) and were stable during the simulation, the ligand in chain A showed a high average RMSD of 1.885 Å, suggesting that the hydrogen bond between His164 and the ligand may contribute to ligand stabilization within the binding pocket. The structure of the C145A mutants is shown in Fig. 7e. The C145A mutant has been shown to reduce the binding affinity of nirmatrelvir [40], which is consistent with our results. For the M165A mutant (Fig. 7f), the hydrogen bond between the ligand and the binding pocket was significantly decreased (Table S1). In the wild-type structure, the Met165 side chain forms hydrophobic interactions with the dimethylcyclopropylproline moiety (P2) of nirmatrelvir and interacts directly with the N-terminal trifluoroacetamide (P4) of the ligand [22,36]. These interactions were considered to be important for the proper recognition of nirmatrelvir, and their loss in the M165A mutant may contribute to decreased ligand binding.
3.3 Binding Free Energy Calculation
Binding free energy calculations for nirmatrelvir were conducted for the wild-type and 13 mutant nirmatrelvir-resistant candidates (Table 4). The bold and underlined ΔG values indicate that the ligands were largely moved (ligand RMSDs >1.552 Å). As listed in Table 4, when the ligand RMSD value exceeded the threshold (1.552 Å), the corresponding binding free energy was higher (lower affinity) compared with that of the wild-type. The binding free energies support the validity of our virtual alanine scan, although the candidate drug-resistant mutants were identified based only on RMSD threshold criteria without binding free energy calculations. Except for the L141A chain B (ΔΔG = +0.11 kcal mol−1) and S144A chain B (ΔΔG = +0.97 kcal mol−1), all boldfaced ligands in Table 4 indicate that the ΔΔG values for nirmatrelvir are ≥ +2.0 kcal mol-1. The largest ΔΔG for nirmatrelvir was observed in E166A chain A (+8.13 kcal mol−1), followed by S1A chain A (+6.71 kcal mol−1). As described previously, Glu166 forms hydrogen bonds with the ligand and interacts with Ser1, whereas MM/GBSA analysis supports the importance of these interactions. Met165 and Pro168 mutations, which are located near Glu166, resulted in notable changes in hydrogen bond frequency between Glu166 and the ligand, but did not increase the ΔΔG as high as that of E166A. Specifically, the binding free energy between M165A chain A and the ligand was ΔΔG = +2.00 kcal mol−1, whereas P168A exhibited +3.02 kcal mol−1 (chain A) and +3.53 kcal mol−1 (chain B). The results suggest that mutations at Met165 or Pro168 may have less of an impact on ligand recognition compared with mutations at Glu166 and Ser1. The F140A chain A showed a ΔΔG of +5.75 kcal mol−1. Because Phe140 interacts with Ser1, this further supports the importance of the interaction between Glu166 and Ser1.
Table 4. Binding free energy for nirmatrelvir to SARS-CoV-2 Mpro calculated by MM/GBSA method
| Protein name
|
Chain
|
Mean
|
SD
|
ΔΔG
|
| Wild-type |
A |
−51.37 |
3.57 |
— |
| B |
−54.41 |
3.51 |
— |
| S1A |
A |
−44.66
|
3.38 |
6.71 |
| B |
−49.55
|
3.44 |
4.86 |
| H41A |
A |
−46.65
|
2.92 |
4.72 |
| B |
−49.15
|
3.32 |
5.26 |
| Y54A |
A |
−47.58
|
3.56 |
3.79 |
| B |
−50.10
|
3.10 |
4.31 |
| F140A |
A |
−45.62
|
3.32 |
5.75 |
| B |
−49.27 |
3.09 |
5.14 |
| L141A |
A |
−47.19
|
3.22 |
4.18 |
| B |
−54.30
|
3.16 |
0.11 |
| G143A |
A |
−49.84 |
3.33 |
1.53 |
| B |
−53.52
|
3.46 |
0.89 |
| S144A |
A |
−45.27
|
3.15 |
6.10 |
| B |
−53.44
|
3.16 |
0.97 |
| C145A |
A |
−49.02
|
2.99 |
2.35 |
| B |
−58.02 |
3.30 |
−3.61 |
| M165A |
A |
−49.37
|
3.37 |
2.00 |
| B |
−48.50 |
3.08 |
5.91 |
| E166A |
A |
−43.24
|
3.17 |
8.31 |
| B |
−47.84 |
3.42 |
6.57 |
| P168A |
A |
−48.35
|
3.45 |
3.02 |
| B |
−50.88
|
3.29 |
3.53 |
| Q189A |
A |
−46.06
|
3.05 |
5.31 |
| B |
−52.10
|
3.07 |
2.31 |
| Q192A |
A |
−45.94
|
3.32 |
5.43 |
| B |
−51.71 |
3.20 |
2.70 |
Binding free energies are expressed in kcal mol−1. Standard deviations are expressed as SD. ΔΔG was the difference of ΔG, that is, ΔG of the wild-type minus ΔG of the mutant. Binding free energies of the ligands whose RMSD values exceeded the threshold (1.552 Å) are shown in bold and underlined.
4. Conclusion
In the present study, MD simulations were conducted on the complex structure between SARS-CoV-2 Mpro and nirmatrelvir to analyze its 3D structure, along with a virtual alanine scan to predict amino acid residues that potentially confer drug resistance. We constructed 23 artificial mutants of the Mpro–nitmatrelvir complexes. The 13 mutants (S1A, H41A, Y54A, F140A, L141A, G143A, S144A, C145A, M165A, E166A, P168A, Q189A, and Q192A) may reduce nirmatrelvir recognition, while maintaining the overall structure of the protein and the structure of the catalytic dyad. These 13 mutants are potential candidates for nirmatrelvir resistance. The results suggest that our approach can apply not only to non-covalent inhibitors, but also to covalent inhibitors, and is useful for predicting drug-resistant mutations. The MD simulations of the wild-type Mpro–nirmatrelvir complex indicate that nirmatrelvir is recognized by Mpro before covalent bond formation. Based on our results, a virtual alanine scan and MD simulations are expected to apply to the pre-covalent state of covalent inhibitor–protein complexes.
The inhibitory activities of nirmatrelvir for the 13 mutants other than G143S, which was simulated in this study, have not yet been experimentally evaluated. In the future, we plan to collaborate with experimental researchers to verify whether the 13 mutants truly confer resistance to nirmatrelvir. In addition, if new nirmatrelvir-resistant Mpro mutants are experimentally identified, we will reassess the validity of this approach based on the corresponding experimental findings. This study aims to rapidly identify amino acid residues that contribute to drug resistance in emerging infectious diseases. Thus, we did not perform simulations of Mpro-substrate complex because of the time cost and predicted resistance using only simulations of the Mpro-inhibitor complex. Therefore, we estimated the enzymatic activity based on the RMSD values of the catalytic dyad and the protein main chain. The crystal structure of the substrate-bound SARS-CoV-2 Mpro is available in the Protein Data Bank (PDB ID: 7MGS), and we are currently conducting MD simulations using this natural substrate complex. We plan to further refine our approach for evaluating enzymatic activity based on the results of these ongoing simulations.
Supplemental information
The ligand RMSDs of each virtual-alanine-scanned mutant are shown in Figure S1, and the RMSDs of the catalytic dyad are shown in Figure S2. The hydrogen bond frequencies with the ligand during the final 10 ns for each virtual-alanine-scanned mutant are summarized in Table S1.
Acknowledgments
This work was supported by Nagai Memorial Research Scholarship from the Pharmaceutical Society of Japan and a Grant-in-Aid for Scientific Research (25H01706) from the Japan Society for the Promotion of Science.
References
- [1] Andersen, K. G.; Rambaut, A.; Lipkin, W. I.; Holmes, E. C.; Garry, R. F. The proximal origin of SARS-CoV-2. Nat. Med. 2020, 26, 450–452. DOI: 10.1038/s41591-020-0820-9
- [2] Wu, C. R.; Yin, W. C.; Jiang, Y.; Xu, H. E. Structure genomics of SARS-CoV-2 and its Omicron variant: drug design templates for COVID-19. Acta Pharmacol. Sin. 2022, 43 (12), 3021–3033. DOI: 10.1038/s41401-021-00851-w
- [3] Wu, F.; Zhao, S.; Yu, B.; Chen, Y.-M.; Wang, W.; et al. New Coronavirus Associated with Human Respiratory Disease in China. Nature 2020, 579, 265–269. DOI: 10.1038/s41586-020-2008-3
- [4] Pillaiyar, T.; Manickam, M.; Namasivayam, V.; Hayashi, Y.; Jung, S. H. An Overview of Severe Acute Respiratory Syndrome-Coronavirus (SARS-CoV) 3CL Protease Inhibitors: Peptidomimetics and Small Molecule Chemotherapy. J. Med. Chem. 2016, 59, 6595–6628. DOI: 10.1021/acs.jmedchem.5b01461
- [5] Jin, Z.; Du, X.; Xu, Y.; Deng, Y.; Liu, M.; et al. Structure of Mpro from SARS-CoV-2 and discovery of its inhibitors. Nature 2020, 582, 289–293. DOI: 10.1038/s41586-020-2223-y
- [6] Huang, C.; Wei, P.; Fan, K.; Liu, Y.; Lai, L. 3C-like Proteinase from SARS coronavirus catalyzes substrate hydrolysis by a general base mechanism. Biochemistry 2004, 43, 4568–4574, DOI: 10.1021/bi036022q
- [7] Ratia, K.; Saikatendu, K. S.; Santarsiero, B. D.; Barretto, N.; Baker, S. C.; et al. Severe acute respiratory syndrome coronavirus papain-like protease: structure of a viral deubiquitinating enzyme. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 5717– 5722. DOI: 10.1073/pnas.0510851103
- [8] Li, C.; Qi, Y.; Teng, X.; Yang, Z.; Wei, P.; et al. Maturation mechanism of severe acute respiratory syndrome (SARS) coronavirus 3C-like proteinase. J. Biol. Chem. 2010, 285, 28134– 28140. DOI: 10.1074/jbc.M109.095851
- [9] Zhang, L.; Lin, D.; Sun, X.; Curth, U.; Drosten, C.; et al. Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved α-ketoamide inhibitors. Science 2020, 368, 409–412. DOI: 10.1126/science.abb3405
- [10] Anand, K.; Palm, G. J.; Mesters, J. R.; Siddell, S. G.; Ziebuhr, J.; et al. Structure of coronavirus main proteinase reveals combination of a chymotrypsin fold with an extra alpha-helical domain. EMBO J. 2002, 21, 3213–3224. DOI: 10.1093/emboj/cdf327
- [11] Chuck, C.-P.; Chong, L.-T.; Chen, C.; Chow, H.-F.; Wan, D. C.-C.; et al. Profiling of substrate specificity of SARS-CoV 3CL. PLoS One 2010, 5 (10), e13197. DOI: 10.1371/journal.pone.0013197
- [12] Kumar, V.; Shin, J. S.; Shie, J.-J.; Ku, K. B.; Kim, C.; et al. Identification and evaluation of potent Middle East respiratory syndrome coronavirus (MERS-CoV) 3CLPro inhibitors. Antiviral Res. 2017, 141, 101–106. DOI: 10.1016/j.antiviral.2017.02.007
- [13] Ferreira, J. C.; Fadl, S.; Villanueva, A. J.; Rabeh, W. M. Catalytic dyad residues His41 and Cys145 impact the catalytic activity and overall conformational fold of the main SARS-CoV-2 protease 3-chymotrypsin-like protease. Front. Chem. 2021, 9, 692168. DOI: 10.3389/fchem.2021.692168
- [14] Naderi Beni, R.; Elyasi-Ebli, P.; Gharaghani, S.; Seyedarabi, A. In silico studies of anti-oxidative and hot temperament-based phytochemicals as natural inhibitors of SARS-CoV-2 Mpro, PLoS One 2023, 18, e0295014. DOI: 10.1371/journal.pone.0295014
- [15] Ullrich, S.; Nitsche, C. The SARS-CoV-2 main protease as drug target. Bioorg. Med. Chem. Lett. 2020, 30, 127377. DOI: 10.1016/j.bmcl.2020.127377
- [16] Rut, W.; Groborz, K.; Zhang, L.; Sun, X.; Zmudzinski, M.; et al. SARS-CoV-2 Mpro inhibitors and activity-based probes for patient-sample imaging. Nat. Chem. Biol. 2021, 17 (2), 222– 228. DOI: 10.1038/s41589-020-00689-z
- [17] Anand, K.; Ziebuhr, J.; Wadhwani, P.; Mesters, J. R.; Hilgenfeld, R. Coronavirus main proteinase (3CLpro) structure: Basis for design of anti-SARS drugs. Science. 2003, 300, 1763– 1767. DOI: 10.1126/science.1085658
- [18] Zhang, L.; Lin, D.; Kusov, Y.; Nian, Y.; Ma, Q.; et al. α-Ketoamides as broad-spectrum inhibitors of coronavirus and enterovirus replication: Structure-based design, synthesis, and activity assessment. J. Med. Chem. 2020, 63 (9), 4562– 4578. DOI: 10.1021/acs.jmedchem.9b01828
- [19] Unoh, Y.; Uehara, S.; Nakahara, K.; Nobori, H.; Yamatsu, Y.; et al. Discovery of S-217622, a noncovalent oral SARS-CoV-2 3CL protease inhibitor clinical candidate for treating COVID-19. J. Med. Chem. 2022, 65, 6499–6512. DOI: 10.1021/acs.jmedchem.2c00117
- [20] Owen, D. R.; Allerton, C. M. N.; Anderson, A. S.; Aschenbrenner, L.; Avery, M.; et al. An oral SARS-CoV-2 Mpro inhibitor clinical candidate for the treatment of COVID-19. Science 2021, 374, 1586–1593. DOI: 10.1126/science.abl4784
- [21] Citarella, A.; Dimasi, A.; Moi, D.; Passarella, D.; Scala, A.; et al. Recent advances in SARS-CoV-2 main protease inhibitors: From Nirmatrelvir to future perspectives. Biomolecules 2023, 13, 1339. DOI: 10.3390/biom13091339
- [22] Yang, K. S.; Leeuwon, S. Z.; Xu, S.; Liu, W. R. Evolutionary and Structural Insights about Potential SARS-CoV-2 Evasion of Nirmatrelvir. J. Med. Chem. 2022, 65, 8686–8698. DOI: 10.1021/acs.jmedchem.2c00404
- [23] Mengist, H. M.; Dilnessa, T.; Jin, T. Structural basis of potential inhibitors targeting SARS-CoV-2 main protease. Front. Chem. 2021, 9, 622898. DOI: 10.3389/fchem.2021.622898
- [24] Frecer, V.; Miertus, S. Antiviral agents against COVID-19: Structure-based design of specific peptidomimetic inhibitors of SARS-CoV-2 main protease. RSC Adv. 2020, 10, 40244–40263. DOI: 10.1039/d0ra08304f
- [25] Hammond, J.; Leister-Tebbe, H.; Gardner, A.; Abreu, P.; Bao, W.; et al. Oral Nirmatrelvir for High-Risk, Nonhospitalized Adults with Covid-19. N. Engl. J. Med. 2022, 386, 1397–1408. DOI: 10.1056/NEJMoa2118542
- [26] JA OS. EUA 105 Pfizer Paxlovid LOA (12222021). In: Administration USFaD, editor, 2021.
- [27] Hashemian, S. M. R.; Sheida, A.; Taghizadieh, M.; Memar, M. Y.; Hamblin, M. R.; et al. Paxlovid (Nirmatrelvir/Ritonavir): A new approach to COVID-19 therapy?. Biomed. Pharmacother. 2023, 162, 114367. DOI: 10.1016/j.biopha.2023.114367
- [28] Sheik Amamuddy, O.; Verkhivker, G. M.; Tastan Bishop, O. Impact of early pandemic stage mutations on molecular dynamics of SARS-CoV-2 Mpro. J. Chem. Inf Model 2020, 60, 5080– 5102. DOI: 10.1021/acs.jcim.0c00634
- [29] Nakayoshi, T.; Kato, K.; Kurimoto, E.; Oda, A. Virtual alanine scan of the main protease active site in severe acute respiratory syndrome coronavirus 2. Int. J. Mol. Sci. 2021, 22, 9837. DOI: 10.3390/ijms22189837
- [30] Mizuno, A.; Nakayoshi, T.; Kato, K.; Kurimoto, E.; Oda, A. Computational estimation of residues involving resistance to the SARS-CoV-2 main protease inhibitor ensitrelvir based on virtual alanine scan of the active site. Bio. Pharm. Bul. 2024, 47(5), 967–977. DOI: 10.1248/bpb.b24-00031
- [31] Gordon, J. C.; Myers, J. B.; Folta, T.; Shoja, V.; Heath, L. S.; et al. H++: A server for estimating pKas and adding missing hydrogens to macromolecules. Nucleic Acids Res. 2005, 33 (Web. Server), W368–W371. DOI: 10.1093/nar/gki464
- [32] Case, D. A.; Belfon, K.; Ben-Shalom, I. Y.; et al. Amber 2020. University of California. San Francisco 2020.
- [33] Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; et al. ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696–3713. DOI: 10.1021/acs.jctc.5b00255
- [34] Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 1157– 1174. DOI: 10.1002/jcc.20035
- [35] Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. Numerical integration of the Cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes. J. Comput. Phys. 1997, 23, 327–341. DOI: 10.1016/0021-9991(77)90098-5
- [36] Noske, G. D.; de Souza Silva, E.; de Godoy, M. O.; Dolci, I.; Fernandes, R. S.; et al. Structural basis of nirmatrelvir and ensitrelvir activity against naturally occurring polymorphisms of the SARS-CoV-2 main protease. J. Biol. Chem. 2023, 299, 103004. DOI: 10.1016/j.jbc.2023.103004
- [37] Genheden, S.; Ryde, U. The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities, Expert Opin. Drug Discov. 2015, 10, 449–461. DOI: 10.1517/17460441.2015.1032936
- [38] Onufriev, A.; Bashford, D.; Case, D. A. Exploring protein native states and large-scale conformational changes with a modified generalized born model. Proteins 2004, 55, 383–394. DOI: 10.1002/prot.20033
- [39] Miller 3rd, B. R.; McGee Jr, T. D.; Swails, J. M.; Homeyer, N.; Gohlke, H.; et al. MMPBSA.py: An efficient program for end-state free energy calculations. J. Chem. Theory Comput. 2012, 8, 3314–3321. DOI: 10.1021/ct300418h
- [40] Kovalevsky, A.; Aniana, A.; Coates, L.; Bonnesen, P. V.; Nashed, N. T.; et al. Contribution of the catalytic dyad of SARS-CoV-2 main protease to binding covalent and noncovalent inhibitors, J. Biol. Chem. 2023, 299 (7), 104886. DOI: 10.1016/j.jbc.2023.104886
