2024 年 47 巻 5 号 p. 967-977
Ensitrelvir is a noncovalent inhibitor of the main protease (Mpro) of severe acute respiratory syndrome coronavirus 2. Acquisition of drug resistance in virus-derived proteins is a serious therapeutic concern, and drug resistance occurs due to amino acid mutations. In this study, we computationally constructed 24 mutants, in which one residue around the active site was replaced with alanine and performed molecular dynamics simulations to the complex of Mpro and ensitrelvir to predict the residues involved in drug resistance. We evaluated the changes in the entire protein structure and ligand configuration in each of these mutants and estimated which residues were involved in ensitrelvir recognition. This method is called a virtual alanine scan. In nine mutants (S1A, T26A, H41A, M49A, L141A, H163A, E166A, V186A, and R188A), although the entire protein structure and catalytic dyad (cysteine (Cys)145 and histidine (His)41) were not significantly moved, the ensitrelvir configuration changed. Thus, it is considered that these mutants did not recognize ensitrelvir while maintaining Mpro enzymatic activities, and Ser1, Thr26, His41, Met49, Leu141, His163, Glu166, Val186, and Arg188 may be related to ensitrelvir resistance. The ligand shift noted in M49A was similar to that observed in M49I, which has been shown to be experimentally ensitrelvir resistant. These findings suggest that our research approach can predict mutations that incite drug resistance.
On November 22, 2022, ensitrelvir (Xocova), developed by Shionogi, was urgently approved as an antiviral agent for severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) that caused the global coronavirus disease 2019 (COVID-19) pandemic. Ensitrelvir impedes viral growth by inhibiting the SARS-CoV-2 main protease (Mpro).1)
SARS-CoV-2 binds its spike glycoprotein to angiotensin converting enzyme 2, which allows SARS-CoV-2 to invade host cells and release its positive-sense single-stranded RNA.1) The single-stranded RNA is then translated into two large polyproteins, pp1a and pp1ab. Pp1a and pp1ab are cleaved by Mpro and papain-like protease into 16 nonstructural proteins (nsp1–16). SARS-CoV-2 Mpro is a cysteine protease that cleaves 11 distinct sites on pp1a and pp1ab, including nsp5 (Mpro).1,2) Therefore, pp1a and pp1ab undergo initial cleavage of their own N- and C-termini via a mechanism called autoprocessing.3–5) Mpro cleaves pp1a and pp1ab to generate 11 nonstructural proteins that eventually become essential proteins for virus replication, including nsp12, that is, RNA-dependent RNA polymerase (RdRp).2,3,6) Thus, Mpro inhibitors were anticipated to be anti-SARS-CoV-2 drugs that hamper viral replication, as well as neutralizing antibodies7) and RdRp inhibitors.8) Mpro inhibitors have also been noted from the point of view of drug repositioning because Mpro inhibitors have already been confirmed to be effective as antiviral drugs for treating human immunodeficiency virus (HIV) and hepatitis C virus (HCV).9,10) Additionally, Mpro has a lower likelihood of mutation compared to spike proteins and has different substrate specificity from human protease.1,11) Therefore, Mpro inhibitors are expected to be effective for treating COVID-19, and researches have been conducted to explore these inhibitors. Consequently, Pfizer has successfully developed the first orally administrable covalent Mpro inhibitor, nirmatrelvir.12) SARS-CoV-2 Mpro is a homodimer consisting of two subunits, chains A and B. Each subunit contains 306 residues and comprises three domains (domains I–III). Domains I and II, composed of residues 8–101 and 102–184, respectively, are composed of antiparallel β-barrels with similarities to trypsin-like serine proteases. Domain III (residues 201–306) is characterized by a cluster of five α-helices and connected to domain II via a loop region (residues 185–200). The substrate-binding sites, including histidine (His)41–cysteine (Cys)145 catalytic dyads, are located at the cleft between domains I and II and play important roles in proteolytic activities. N-terminus residues 1–7 are located between domains II and III and play roles in the formation of the substrate-binding site and the joint of each subunit.13–15) The substrate-binding site consists mainly of four pockets, S1′, S1, S2, and S4. In the apo structure, the S1 pocket in each subunit requires glutamine residues. The S1 pocket comprises Phe140, Leu141, Asn142, His163, Glu166, and His172. Its adjacent S1′ pocket, which mainly recognizes serine (Ser)/alanine (Ala)/glycine (Gly), is lined with residues Thr24 and Thr26. The S2 pocket is a highly hydrophobic subpocket surrounded by His41, Met49, Tyr54, Met165, and Asp187 and has a preference for leucine (Leu)/Met/phenylalanine (Phe)/Val. The S4 pocket is a semienclosed subsite composed of Leu165, Phe185, and Gln189.6,15–18)
The His41–Cys145 catalytic dyad plays a central role in enzymatic function. The side-chain thiol sulfur of Cys145 nucleophilically attacks the amide carbon of the substrate, and a tetrahedral intermediate is formed. In this process, the side-chain imidazole group of His41 plays a role in increasing the nucleophilicity of Cys145. Subsequently, proton donation from His41 to the tetrahedral intermediate results in the cleavage of the peptide bond of the substrate, and the side chain of Cys145 is acylated. Finally, the acylated Cys145 undergoes nucleophilic attack by a water molecule, and the catalytic dyad is regenerated via a thiohemiketal intermediate.14,15)
Although ensitrelvir is an Mpro inhibitor, it is a noncovalent inhibitor contrary to nirmatrelvir. Ensitrelvir satisfies the drug metabolism and pharmacokinetics profile and dose-dependent antiviral activity in vivo and proven efficacious and administered orally once per day. For the development of ensitrelvir, the structure-based drug design strategy was adopted. Figure 1A shows the chemical structure of ensitrelvir. Ensitrelvir combines a triazinane-2,4-dione with 1-methyl-1H-1,2,4-triazole (P1), 6-chloro-2-methy-2H-lindazole (P1′), and a 2,4,5-trifluorobenzyliz (P2) substituent to optimize the binding mode. In the X-ray cocrystal structure of Mpro complexed with ensitrelvir (PDB ID: 7vu6),1) the P1 unit, which forms a hydrogen bond with His163, fits into the S1 pocket. The P2 moiety that interacts with His41 occupies the S2 pocket. The P1′ moiety has a hydrogen bond with Thr26 and hydrophobic contact with Met49. Figure 1B shows the three dimensional (3D) structures of the significant interacting residues and the binding poses of ensitrelvir.
(B) The 3D structures of the significant interacting residues and the binding poses of ensitrelvir. The amino acid residues are represented by salmon pink. The ligand is represented by faded green.
Although Mpro is considered relatively difficult to mutate,1) F305L has already been found mainly in an England outbreak.19) Furthermore, P132H and K90R have been considered circulating strains of SARS-CoV-2, which have been reported to occur in 2505897 and 99665, respectively.20,21) Because amino acid mutations can result in drug resistance, Mpro mutations can reduce or nullify the efficacy of ensitrelvir. Although 11 major Mpro mutations (G15S, T21I, T24I, K88R, L89F, K90R, P108S, P132H, A193V, H246Y, and A255V) are also inhibited by ensitrelvir, M49I/L/T/V and E166A/V mutations have experimentally confirmed resistance to ensitrelvir.21–28)
If the Mpro mutants retain the protein structure but fail to recognize ensitrelvir, the mutants may exhibit drug resistance to ensitrelvir. Thus, in this study, we computationally constructed artificial mutants, in which each residue around the ligands was replaced with alanine and estimated which residues important for ensitrelvir recognition by performing molecular dynamics (MD) simulations (referred to as virtual alanine scan methods). Some computational approaches combine alanine scan and MM/PBSA.29,30) Because our study aimed at a rapid estimation of candidates of drug resistance mutants, we did not perform free energy calculations that require longer MD simulations. The comparison between our results and the results of MM/GBSA indicates that our method without free energy calculations can obtain similar results to MM/GBSA. In a previous study,31) we performed virtual alanine scan with a complex of Mpro and NCL-00024905, a noncovalent Mpro inhibitor. The ligand configuration was altered in H163A and E166A mutants, although the structures of the entire protein and catalytic-dyad did not significantly change. Therefore, mutations in His163 and Glu166 could lead to drug resistance. Additionally, the E166A mutation affected the positions of Ser1 and His172 and changed the shape of the pocket, resulting in a reduction in the interaction between Ala166 and the neighboring residues. Hence, Glu166 is important for the retention of pocket shape. Virtual alanine scan is a useful tool not only for predicting drug resistance but also for revealing the Mpro structural features. We expect that this study can predict the possible drug resistance to ensitrelvir.
The initial structure of the complex consisting of Mpro and ensitrelvir was constructed from an experimentally determined structure registered in the Protein Data Bank (PDB ID: 8dz0). In 8dz0, Mpro is a symmetric homodimer consisting of chains A and B, with the ligands located in both active sites of the dimer. The atomic coordinates of the four C-terminal residues in chain A were not experimentally determined, and these coordinates were complemented using those of chain B before the simulation. The protonation states of the titratable groups at pH 7.0 were determined based on pKa values calculated using H++ web server.32) For the calculations of pKa values, the internal and external dielectric constants of protein were 10 and 80, respectively. In 8dz0, the side chains of the two His residues (His41 and His163) were in close proximity to the ensitrelvir. Since the chlorine atom of the ensitrelvir was located near the δ-nitrogen atom of His41, the side-chain imidazole group of His41 was considered to be protonated at δ-position. Since the triazole group of the ensitrelvir was located near the ε-nitrogen atom of His163, the side-chain imidazole group of His163 were considered to be protonated at ε-position. The amino group of Ser1 is surrounded by the side chain carbonyl group of Glu166 and the main chain carbonyl group of Phe140, so it is possible that the N-terminus of Ser1 is protonated state (NH3+). MD simulations were conducted on the wild-type Mpro and ensitrelvir complex according to the procedure described below. We performed virtual alanine scan of the structure after MD simulations. Comparison of the computational structures of the wild-type and alanine-scanned mutants contributes to estimating the role of premutated amino acid residues in Mpro. Figure 2 illustrates the residues within a radius of 4 Å from the ligands in Mpro after the simulations. Each of these residues was replaced with alanine by erasing the side chain coordinates, rewriting the residue name to “ALA,” and complementing the methyl group of the side chain with tleap module of AmberTools. Table 1 summarized the residues replaced in this study and their major interactions. We generated 24 artificial-mutant Mpro–ensitrelvir complexes, and conducted the energy minimization and MD simulations for each complex.
Alanine-scanned residues are illustrated in gray. Ensitrelvir is illustrated in white.
| Amino acid residue within a radius of 4 Å from the ligands | Interaction |
|---|---|
| Ser1 | Hydrogen bond with Phe140 and Glu166 |
| Thr25 | — |
| Thr26 | Hydrogen bond with Gln19, Val20, Thr21, and ligand P1′ moiety |
| Leu27 | — |
| His41 | Recognition of ligand P2 moiety |
| Thr45 | — |
| Ser46 | — |
| Met49 | Hydrophobic and/or CH–π interaction with ligand and P2 moiety |
| Phe140 | π–π Interaction with His163 hydrogen bond with Ser1 |
| Leu141 | Hydrogen bond with Ser144 |
| Asn142 | — |
| Gly143 | Hydrogen bond with a triazinane-2,4-dione of ligand |
| Ser144 | Hydrogen bond with Leu141 |
| Cys145 | Hydrogen bond with a triazinane-2,4-dione of ligand |
| His163 | π–π interaction with Phe140 hydrogen bond with Ser144 and P1 |
| His164 | — |
| Met165 | Hydrophobic and/or CH–π interaction with ligand P2 moiety |
| Glu166 | Hydrogen bond with Ser1 and a triazinane-2,4-dione of ligand |
| His172 | — |
| Val186 | Hydrogen bond with Arg188 |
| Asp187 | — |
| Arg188 | — |
| Gln189 | Hydrogen bond with ligand P1′ |
| Thr190 | — |
AMBER ff14SB force field33) and general AMBER force field 2 (GAFF2)34) were employed to the protein and ligands, respectively. The ligand atoms were assigned atom types using the antechamber module of AmberTools20, and their atomic charges were determined using the AM1-BCC methods.35,36) The systems were solvated with TIP3P water molecules under periodic boundary conditions that were spaced with 10 Å margins from the protein surface. To neutralize the system charges, chloride (Cl−) ions were placed. First, 1000-cycle energy minimizations were performed to optimize the coordinates of water molecules and Cl− ions with strong positional restraints of 500 kcal mol−1 Å−2 applied to the proteins and ligands, and then 2500-cycle minimization for the whole system was performed. After minimizations, temperature-increasing MD simulations were performed for 20 ps from 0 to 300 K with weak positional restraints of 10.0 kcal mol−1 Å−2 employed to the proteins and ligands, followed by equilibration of MD simulations for the wild type and mutants under constant pressure of 1 bar and temperature of 300 K for 1000 and 100 ns, respectively. For all MD simulations, the time step was 2 fs. The SHAKE algorithm37) was used to constrain bond stretching involving hydrogen atoms. Amber20 software38) was used for energy minimization and MD simulations.
To analyze the calculated structures, the root mean square deviations (RMSDs) for the main-chain atoms, catalytic dyads (His41 and Cys145), and ligands were calculated along the MD trajectories using the initial structure of the equilibration MD simulations as the reference structure. RMSDs were used as an indicator of the change in the structure of proteins and ligands. The RMSDs of the ligand and catalytic dyad in chains A and B were separately determined. To investigate the structural flexibility of the protein, root mean square fluctuations (RMSFs) for all Cα atoms were calculated. The hydrogen bond frequencies were analyzed on the criteria for the hydrogen bonds with the cutoff distances of 3.5 Å and the angles of greater than 120°. Moreover, π–π stacking interactions were analyzed based on the distance from the center of mass of one aromatic ring to that of the other. CH–π interactions were analyzed based on the distance from the center of mass of the aromatic ring to the hydrogen atom. The radius of gyration (Rg) was calculated for the Cα atoms. Rg values are indicators of the protein compactness and size.39) The RMSDs, RMSFs, hydrogen bond frequencies, and distances were analyzed using the cpptraj module of AmberTools 20 for the final 10-ns MD trajectories.
The binding free energies of complexes of wild-type/mutated Mpro and ensitrelvir were calculated using MM/GBSA approaches.40) Final 10-ns trajectories of the complexes were used for MM/GBSA calculations. In this study, only complexes were simulated, and the trajectories of free Mpro and ensitrelvir were obtained by removing either the latter and former from the trajectories of complex, respectively. The modified generalized Born (GB) model (igb = 5)41) was used as the continuum solvent model, and the salt concentration was set to 0.15 M. MM/GBSA calculations were performed using MMPBSA.py program of AmberTools.42)
Figure 3 shows the three-dimensional structures of wild-type SARS-CoV-2 Mpro before and after the simulation. Figure 4 shows the RMSDs of the main chain, ligand, and catalytic dyad throughout the MD simulations. As demonstrated in Fig. 4, the entire protein structure, ligand configurations, and catalytic-dyad structures did not remarkably change during the simulations, and the RMSDs of protein backbones and ligands converged. It is considered to maintain the appropriate structure for the wild-type Mpro, in which its enzyme activity was inhibited by ensitrelvir.
(B) Ensitrelvir in chain A. (C) Ensitrelvir in chain B. Before and after MD simulations are reflected as cyan and faded yellow stick models, respectively. Cyan and faded yellow ribbons indicate the overall complex structure before MD simulations (at 0 ns) and after MD simulations (at 1000 ns), respectively.
RMSDs for chains A and B are shown in red and blue lines, respectively.
Table 2 presents the average RMSDs of the main chain, ligands, and catalytic dyad, and Supplementary Figs. S1–S3 show the RMSDs of the main chains, ligands, and catalytic dyads, respectively, throughout the MD simulations for each mutant. As in our previous study,31) we regarded the RMSDs ≥ 1.5 Å as changing compared with the wild type. Although the criteria to identify incorrect ligand pose is set to 2.0 Å in many cases, we determined that this criterion was inappropriate for our study, which involved only short MD simulations. In this study, we adopted the criteria that keep the number of identified drug resistance candidates to 10 or fewer. Thus, it was possible to adjust the number of results obtained in this study by changing the criteria. Drug resistance becomes a problem only when a protein does not recognize an inhibitor retaining its functions. When mutations affect the entire protein structure, it is not important in terms of drug resistance because the enzyme will no longer function. That is, the candidates of drug-resistant Mpro are the mutants in which the RMSDs of ligands show a large value, although the RMSDs of the main chains and catalytic dyads do not. The large ligand RMSDs suggest that the recognition of ensitrelvir is weakened in the mutant Mpro compared with the wild type. In this study, the structural features of mutants were discussed if large RMSDs were obtained from at least one of the subunits because SARS-CoV-2 Mpro is a homodimer with structural symmetry. Additionally, for the fast evaluation of drug resistance, 100-ns simulations were conducted, similar to our previous study.31) Regarding the simulation time, the calculations for each mutant were set to be completed within 24 h of real time to enable rapid execution of numerous computations.
 |
Solid underline: >2.0 Å and dotted underline: >1.5 Å.
After the MD simulations, the mutants in the RMSDs of the main chain and the catalytic dyad did not exceed 1.5 Å; however, those of the ligands that exceeded were S1A, T26A, H41A, M49A, L141A, H163A, E166A, V186A, and R188A. These nine types of mutants were candidates for drug resistance.
Hydrogen Bond Frequencies and Distance AnalysisHis41, Met49, Val186, and Arg188 are residues that comprise the highly hydrophobic S2 pocket.5,6,43,44) Among the nine types of mutants, the residues constituting the S2 pocket were His41, Met49, Val186, and Arg188. M49I/L/T/V and R188K/S mutations have experimentally demonstrated ensitrelvir resistance,22,24–28) and M49I mutations have been detected from 1,062 samples until July 2023 in GISAID.20) Thus, MD simulations can predict drug resistance of Mpro mutants, and structural elucidation of mutants is useful for investigating the cause of drug resistance.
In the wild-type chain B obtained after simulations, the side chain of Met49 contacts the side chain of His41, which recognizes a P2 substituent of ensitrelvir and forms hydrophobic interactions and/or CH–π interactions with P1′. Furthermore, Gln189 in close proximity to Met49 formed a hydrogen bond at 76.5% with the ligand (Fig. 5A, Table 3). Gln189 was the only residues around S2 pocket that formed highly-frequent hydrogen bond with the ligand. The remarkable differences in hydrogen bond frequencies were observed at these residues between the wild type and mutant, and were considered to affect the changes of protein structures and ligand positions. Only residues mentioned in the main text as being involved in hydrogen bonds were included in Fig. 5, and hydrogen bonds were indicated by dotted lines. The differences in hydrogen bond frequencies between the wild type and mutant were summarized in Table 3.
(A) Superimposition of the wild type and M49A. The wild-type and M49A atoms and chains are presented in white and red, respectively. Ensitrelvir is represented by the ball-and-stick model. His41, Met49, Met165, and Gln189 are represented by the thin tube model. (B) Superimposition of the wild type and V186A. Wild-type and V186A atoms and chains are presented in white and yellow–green, respectively. Ensitrelvir is represented by the ball-and-stick model. His41, Met49, Cys145, Met165, Val186, Arg188, and Gln189 are represented by the thin tube model. (C) Superimposition of the wild type and R188A. Wild-type and R188A atoms and chains are represented in white and faded azure, respectively. Ensitrelvir is represented by the ball-and-stick model. Arg188 and Gln189 are represented by the thin tube model.
| Hydrogen bond | Protein name | Frequency/% |
|---|---|---|
| Ligand@N12-Gln189@HE21 | Wild type | 76.5 |
| M49A | 71.5 | |
| V186A | 0 | |
| R188A | 0.04 | |
| L141A | 59.8 | |
| Val186@O-Arg188@NH | Wild type | 92.2 |
| V186A | 0 | |
| Ligand@O09-Cys145@NH | Wild type | 84.6 |
| V186A | 67.9 | |
| T26A | 2.32 | |
| L141A | 93.2 | |
| Met49@O-Arg188@HE | V186A | 89.3 |
| Met49@O-Arg188@HH21 | V186A | 61.2 |
| Thr26@OG1-Gln19@He22 | Wild type | 94.7 |
| T26A | 0 | |
| Thr26@HG1-Val20@O | Wild type | 70.1 |
| T26A | 0 | |
| Thr26@HG1-Thr21@OG1 | Wild type | 54.8 |
| T26A | 0 | |
| Thr26@O-Asn119@HD22 | Wild type | 0 |
| T26A | 62.5 | |
| Leu141@O-Ser144@NH | Wild type | 66.1 |
| L141A | 29.0 | |
| Ser144@OG-His163@HE2 | Wild type | 86.1 |
| L141A | 78.7 | |
| Phe140@NH-Ser1@O | Wild type | 83.0 |
| L141A | 62.5 |
While the average distance in the final 10 ns trajectories between the ring of His41 and P2 was 5.13 Å and between the ring of P2 and the side-chain CH of Met165 was 4.61 Å in the wild type chain B, the corresponding distances were 4.49 and 3.94 Å, respectively, in the M49A mutant chain B (Figs. 5A, 6A, B, Table 4). In the B chain of M49A, the thioether of the side chain was replaced with a methyl group, and the frequency of hydrogen bonds between the ligand and Gln189 near Met49 slightly decreased from 76.5% (wild type) to 71.5% (M49A) (Fig. 5A, Table 3). Although the recognition around P1′ weakened and the affinity between P2 of the ligand and the S2 pocket increased in M49A, the difference in interactions between the wild-type and M49A mutant was small. Fragment molecular orbital calculations revealed that the ability of Met49 to recognize P1′ was weak, and its dispersion interaction was −1.3 kcal mol−1.18) Thus, the interaction with Met49 may be not essential for ligand recognition. In another study, M49I demonstrated ensitrelvir resistance, and its crystal structure (PDB ID: 8dz1) showed that ensitrelvir shifted toward the S2 pocket compared with the wild type (PDB ID: 8dz0) because of the hydrophobicity of the isoleucine side chain.27) As the ligand shift toward the S2 pocket was noted in this study, this ligand shift may have a significant effect on ensitrelvir resistance. In the B chain of the wild type, the average distance from Met49 to P1′ was 4.42 Å, while it was 8.70 Å in the wild-type A chain (Figs. 6A, B, Table 4). In the wild-type and M49A A chain, P1′ interacted with His41. The interactions between the P2 and S2 pocket weakened in the A chain of the wild type, suggesting that the ring of P2 significantly shifted in the A chain of M49A during MD simulations. Shifts of His41 into the space vacated by Ala substitution were observed in the A and B chains in M49A. Thus far,1,18,27,45) the calculated structure of the B chain is considered to be a more appropriate structure in the wild type and M49A. When the simulation time was extended to 3000 ns, Met49 approached P1′ in the active site in the A chain of the wild type, and the active-site structure of the A chain changed to a structure similar to that of the B chain.
The distance was illustrated by black dashed line. The amino acid residues are represented by the thin tube model. Ensitrelvir is represented by the ball-and-stick model. (A) The wild type B chain (B) M49A B chain (C) The wild type A chain (D) V186A B chain (E) T26A B chain.
| Distance | Protein name | Å |
|---|---|---|
| Ligand P2-His41 | Wild type | 5.13 |
| M49A | 4.49 | |
| V186A | 3.81 | |
| T26A | 4.64 | |
| Ligand P2-Met165@HG3 | Wild type | 4.61 |
| M49A | 3.93 | |
| V186A | 4.04 | |
| Ligand P1′-Met49 in the B chain | Wild type | 4.42 |
| Ligand P1′-Met49 in the A chain | Wild type | 8.70 |
| Val186@HG12-Arg188@HG3 | Wild type | 2.12 |
| Ala186@HB3-Arg188@HB3 | V186A | 4.30 |
Ligand P2, His41, and P1′ are assumed to be the distance from the center of mass of 2,4,5-trifluorobenzyliz, imidazole, and 6-chloro-2-methy-2H-lindazole. Met49 is assumed to be the distance from the center of mass of CE and SD.
As mentioned above, the distance between ensitrelvir and His41 or Met165 can be used as an indicator of the affinity between the P2 and S2 pocket. Previous studies have suggested that the side-chain imidazole of His41 and Met165 recognizes the ensitrelvir P2 moiety.1,18,27,45) As demonstrated in Table 2, for H41A, although the RMSDs of the entire protein and the catalytic dyad were not large values, the ligand RMSDs of the A and B chains exceeded 1.5 Å. Both RMSDs of the ligands in both A and B chains were larger than 1.5 Å, only for H41A and M49A out of the candidates of drug resistance. The changes in ligand configurations may be attributed to the loss of the interaction that was formed in the wild-type Mpro-ligand complex. Thus, evaluating the interactions between His41 and the ligand can be an indicator of ensitrelvir affinity. The S2 pocket in M165A was largely transformed, and the helix constructed by residues 49–51 located on the opposite side of Ala165 across the ligand was unfolded. The ligand RMSDs in M165A continued to increase during the simulation and exceeded 2.4 Å. Thus, Met165, which constitutes the S2 pocket, is considered to be important in recognizing ensitrelvir.
In the wild-type structure obtained after the simulation, Val186 and Arg188 were included in the loop located at the bottom of the S2 pocket, and their side chains were exposed to the solvent (Fig. 5B). Val186 was mutated residue, and Met49 and Arg188 were residues around the mutation site. The average distance in the wild type B chain between these side chains were 2.12 Å (Fig. 6A, Table 4), and the hydrogen bond frequency between the main chains of Val186 and Arg188 was 92.2% (Table 3). In V186A B chain, the hydrogen bond between the main chains of Ala186 and Arg188 disappeared (Table 3). The average distance between the side chains of Val186 and Arg188 increased to 4.30 Å, and Arg188 approached Met49 rather than Ala186 (Figs. 5B, 6D, Table 4). The Arg188 side chain formed two hydrogen bonds with frequencies of 89.3 and 61.2% with the Met49 main chain (Table 3). Consequently, the helix including Met49 and the loop of residues 186–188 were moved (Fig. 5B), and the hydrogen bond between Gln189 and ensitrelvir disappeared (Table 3). Conversely, the average distances between the ring of His41 and P2 decreased from 5.13 Å (the wild type) to 3.81 Å (V186A), and the average distance between the ring of P2 and the side-chain CH of Met165 decreased from 4.61 Å (the wild type) to 4.04 Å (V186A) (Figs. 6A, D, Table 4), suggesting that the affinity between the S2 pocket and P2 increased similar to M49A. In R188A B chain, P2 rotated to match the shape of the S2 pocket, and the hydrogen bond between Gln189 and ensitrelvir disappeared (Table 3). The hydrogen bond cleavage was considered to be attributed to the conformational change of the long side chain of Gln189, and the hydrophilic side chain was excluded from the S2 pocket. Thus, the S2 pocket becoming more hydrophobic is a possibility. Cys145 was the residue consisting of catalytic dyad and formed highly-frequent hydrogen bond with the ligand, suggesting the importance in ligand recognition. In V186A B chain, the hydrogen bond between Cys145 and ensitrelvir decreased from 84.6% (the wild type) to 67.9% (V186A) (Fig. 5B, Table 3), which was the only decrease in the top five hydrogen bond frequencies with the ligand. Thus, the interaction between the ligand and Cys145 may be important for ensitrelvir stability.
Although Thr26 is not a residue constructing the S2 pocket, in the B chain of T26A, ensitrelvir demonstrated a significant decrease in the hydrogen bond frequency with Cys145 from 84.6% (wild type) to 2.32% (T26A) (Fig. 7, Table 3). Thr26 was mutated residue, and Gln19, Val20, Thr21, and Asn119 were residues around the mutation site. The remarkable differences in hydrogen bond frequencies were observed at these residues between the wild type and mutant, and these residues are considered to affect the protein structure and ligand position. Only residues mentioned in the main text as being involved in hydrogen bonds were included in the Fig. 7, and hydrogen bonds were indicated by dotted lines. One of the residues that construct antiparallel β-barrels is Thr26.13,15) In the wild type B chain, the Thr26 side chain formed hydrogen bonds with Gln19, Val20, and Thr21 in a neighboring β-strand at frequencies of 94.7, 70.1, and 54.8%, respectively (Table 3). However, in T26A B chain, the corresponding hydrogen bonds disappeared and a new hydrogen bond whose frequency was 62.5% was formed with Asn119 (Table 3). Hence, P1′ moved to the created space owing to the shifting of the β-strand of Ala29 and P2 approached His41 from 5.13 Å (the wild type) to 4.64 Å (T26A) (Figs. 6A, E, Table 4).
Superimposition of the wild type and T26A. The wild-type and T26A atoms and chains are presented in white and faded violet, respectively. Ensitrelvir, Gln19, Val20, Thr21, Thr26, His41, Asn119, and Cys145 are represented by the stick model.
Leu141 formed a loop with Cys145 (Fig. 8A). In B chain, there were differences in the hydrogen bond frequencies in the loop, decreasing from 66.1% (wild type) to 29.0% (L141A) between residue 141 and Ser144 and from 86.1% (wild type) to 78.7% (L141A) between Ser144 and His163, which constructed a β-strand along with Met165 and Glu166 (Fig. 8A, Table 3). Furthermore, His163 is believed to have a π–π interaction with Phe140 in L141A, and the hydrogen bond frequency between Phe140 and Ser1 decreased from 83.0% (the wild type) to 62.5% (L141A) (Table 3). Ser1 formed a hydrogen bond with Glu166 (Fig. 8B). The loop of residues 140–145 formed a hydrogen bond network with the β-strand including His163, Met165, and Glu166. The hydrogen bond network caused a positional change in Met165, and the shape of the entire S2 pocket was altered (Fig. 8C). The hydrogen bond between Gln189 and ensitrelvir decreased from 76.5% (wild type) to 59.8% (L141A) (Table 3). The ligand shifted in a direction opposite to the S2 pocket and the ring of P2 revolved, leading to an increase in the hydrogen bond frequency between Cys145 and ensitrelvir to 84.6% (the wild type) and 93.2% (L141A), respectively (Table 3). Although the L141A mutation affected the residues constructing the S2 pocket, such as Gln189, His41, and Met165, the helix including Met49 did not significantly move. Thus, the surface area of the hydrophobic region almost changed and the hydrophobic effect may not have increased. In H163A and E166A of the A chain, the shape of the S2 pocket was slightly changed and the ligand shifted in the direction opposite to the S2 pocket. In a previous experimental study, passage of the wild-type virus in the presence of ensitrelvir revealed that E166A was associated with a reduced sensitivity to ensitrelvir.22,26) If the crystal structure of E166A is analyzed, a shift may be observed.
Superimposition of the wild type and L141A. Wild-type and L141A atoms and chains are presented in white and orange, respectively. Distinctive residues are represented by the stick model. Interactions are indicated by dashed lines. (A) Overall illustration of the loop, including the Leu141 and S2 pocket. (B) Loop including L141A and the β-strand including Glu166. (C) S2 pocket.
No significant changes were observed in the loop of residues 140–145 in the B chains of the wild type and L141A (Fig. 8A). RMSFs of these residues were approximately 0.46 Å for both the wild type and L141A B chain, indicating that this loop has minimal flexibility.
Affinity AnalysisCalculated binding free energies were shown in Table 5. In the wild type, the binding free energies to two ligands were comparable. Although the ligand RMSD(s) in the chain A of S1A and in the chains A and B of M49A exceeded 1.5 Å, the binding free energies were lower than those of the wild type. Distance analysis suggested the change of ligand configuration toward S2 pocket and affinity increase in M49A, which supported by MM/GBSA calculations. In the remaining seven mutants, the binding free energies of all ligands with RMSD exceeding 1.5 Å were higher than that of the wild type. In chain B of V186A and chain B of R188A, the binding free energies of ligands with RMSD exceeding 2.0 Å were −42.88 and −44.61 kcal mol−1 and were the second and third highest among the nine mutants. For more detailed study, we plan to perform longer simulations for a limited number of mutants in future study.
| Mutants | Mean | S.D. | |
|---|---|---|---|
| Wild type | Chain A | −53.22 | 5.37 |
| Chain B | −53.15 | 4.06 | |
| S1A | Chain A | −58.17 | 3.78 |
| Chain B | −52.91 | 4.09 | |
| T26A | Chain A | −57.40 | 3.72 |
| Chain B | −50.00 | 4.06 | |
| H41A | Chain A | −48.26 | 6.01 |
| Chain B | −49.18 | 5.20 | |
| M49A | Chain A | −57.17 | 3.93 |
| Chain B | −58.87 | 4.38 | |
| L141A | Chain A | −55.57 | 3.82 |
| Chain B | −50.66 | 3.97 | |
| H163A | Chain A | −40.49 | 3.51 |
| Chain B | −44.65 | 4.28 | |
| E166A | Chain A | −52.11 | 4.00 |
| Chain B | −50.16 | 4.13 | |
| V186A | Chain A | −58.31 | 3.83 |
| Chain B | −42.88 | 3.53 | |
| R188A | Chain A | −49.16 | 3.59 |
| Chain B | −44.61 | 4.10 |
S.D.: standard deviation.
Figure 9 shows the Rg for Cα atoms of His41, Met49, Tyr54, Met165, and Asp187 composed of the S2 pocket. The Rg values of M49A, V186A, and R188A were larger than those of the wild type. The ensitrelvir shift toward S2 pocket may be due to structural change.
In this study, we performed MD simulations and virtual alanine scans to estimate the residues involved in the emergence of resistance to the SARS-CoV-2 Mpro inhibitor ensitrelvir. In the nine mutants (S1A, T26A, H41A, M49A, L141A, H163A, E166A, V186A, and R188A), the entire structure of the Mpro and catalytic dyad did not remarkably change; however, the ensitrelvir configuration changed. Thus, mutations in Ser1, Thr26, His41, Met49, Leu141, His163, Glu166, Val186, and Arg188 can weaken the affinity for the inhibitor while maintaining the Mpro enzymatic function. Our study revealed the positional change of the ensitrelvir to the S2 pocket due to hydrophobicity and spatial influence in M49A, V186A, and R188A. This ligand shift noted in M49A was also reported in the resistance-acquired mutant M49I. These results demonstrated that the virtual alanine scan is a useful tool for predicting drug resistance and is effective in estimating the function of the residues for recognizing inhibitors.
This work was supported by Grants-in-Aid for Scientific Research (17K08257) and (19J23595) from the Japan Society for the Promotion of Science.
The authors declare no conflict of interest.
This article contains supplementary materials.
