Structural Stability and Binding Ability of SARS-CoV-2 Main Protease with GC376: A Stereoisomeric Covalent Ligand Analysis by FMO calculation

Yuya Seki; Chiduru Watanabe; Norihiko Tani; Kikuko Kamisaka; Tatsuya Ohyama; Daisuke Takaya; Teruki Honma

doi:10.1273/cbij.24.13

Abstract

In modeling structures with multiple conformers, the one with the higher occupancy is typically used under implicit understanding. Multiple conformers have rarely been studied to determine whether a structure with a higher ligand occupancy conformer is more stable and has stronger protein-ligand binding than the one with the lower occupancy. We performed fragment molecular orbital (FMO) calculations on complexes with high and low ligand occupancies, comparing their energies, such as total energy and ligand binding energy, to investigate this question. The structures for FMO calculation were the SARS-CoV-2 main protease (Mpro) and the covalent inhibitor, GC376. The initial X-ray crystal structure of the complex (PDB ID: 7CB7) possessed stereo isomers of GC376, namely B1S (R form) and K36 (S form), with different occupancies. Our findings showed that for all structural optimization conditions, the total energy of the high ligand occupancy B1S complex was lower than that of the low ligand occupancy K36 complex. In addition, stronger interfragment interaction energy (IFIE) of B1S than of K36 was confirmed for most structural optimization conditions. These differences between B1S and K36 are caused by the presence and absence of a hydrogen bond between the ligand and His41, as revealed by IFIE and structural geometry analyses between the ligand and each amino acid residue of Mpro. The present study, under several structural optimization conditions, would have the potential to model energetically stable structures by efficiently selecting high occupancy conformers from structures with different ones.

1. Introduction

GC376, a broad-spectrum antiviral drug, has been evaluated for its efficacy in laboratory cats with feline infectious peritonitis [1]. GC376 binds covalently to the amino acid residue of Cys145 (Figure 1A), converting to an aldehyde group by releasing the bisulfite group in vivo (Figure 1B) [2]. GC376 exists in both R and S forms (PDB ligand names: B1S and K36, respectively), each with different stereo configuration of hydroxy groups (Figures 1A and 1B).

Figure 1. M^pro and GC376 complex. Superposition of B1S and K36 for the A-chain in the X-ray crystal structure (PDB ID: 7CB7) in (A). Mechanism of covalent binding between M^pro and GC376 in (B)

This study reports quantum chemical calculations using the fragment molecular orbital (FMO) method [3] on the X-ray crystal structure of the complex between the main protease (M^pro) of SARS-CoV-2 and its inhibitor GC376 which is a reported therapeutic target against COVID-19 (PDB ID: 7CB7) [4]. These results have been available in FMODB (URL: https://drugdesign.riken.jp/FMODB/) [5, 6]. We previously reported FMO calculations for the complex between M^pro and K36 (PDB ID: 6WTT, FMODB ID: 4NVVN) [7, 8], where the structure with the highest occupancy conformer was adopted as a calculation target. In the new X-ray crystal structure used in this study (PDB ID: 7CB7), B1S and K36 were identified as stereo isomers of GC376 in the same binding site, although their occupancies in chains A and B of M^pro differed by 0.53 and 0.47, and 0.30 and 0.70, respectively. In the previous study [8], we adopted the structure with the highest occupancy conformer as a calculation target. In modeling structures with multiple conformers, the conformer with the higher occupancy has generally been employed because of computational cost. This conventional approach is based on the implicit understanding that highly occupied conformations, which are observed in large numbers in a crystal structure, are energetically stable structures. Is it guaranteed that a structure with high occupancy conformers based on conventional modeling is structurally more stable and has stronger protein–ligand binding ability than the one with lower ligand occupancy? There are few cases where multiple conformers have been used and examined to answer this question. Therefore, this paper compares the energies of M^pro complexes with high-occupancy B1S and the low-occupancy K36 using FMO calculations to determine if structures modeled with high occupancy isomers are energetically stable.

2. Results and Discussion

2.1 Structural preparation and FMO calculation

Based on the X-ray crystal structures with multiple conformers (PDB ID: 7CB7), we prepared complexes of the SARS-CoV-2 M^pro with the stereo isomers B1S and K36 as ligands. We confirmed that no amino acid residues with multiple conformers showed different occupancies within a 4.5 Å range around the ligand. Thus, for amino acid residues distant from the ligand in chains A and B that did not affect the ligand-binding mode, the highest occupied conformations were adopted. In case of the A chain, Figure 2 shows conformers of the Arg60 and Arg298 that correspond to this condition and the conformation of the ligand. The Auto-FMO protocol was used for the structural preparation of the complexes [5]. Missing heavy- and hydrogen atoms were added using the “Structure Preparation” and “Protonate 3D” functions in Molecular Operating Environment (MOE) [9]. Structural optimizations were performed while gradually expanding the range of atomic constraints. The five different optimization conditions considered in the study are listed in Table 1. The structural optimization was performed for five conditions: hydrogen atoms only (OptH); hydrogen atoms and solvent molecules (OptHSol); hydrogen atoms, solvent molecules, and ligand molecules (OptHLSol); hydrogen, solvent, ligand molecules, and amino acid residue side chains (OptHSideLSol); and all atoms (OptAll). During structural optimization, the hydrogen atoms were unfixed, while the heavy atoms to be optimized were constrained with tether weight (1.0 kcal Å⁻¹), and the remaining heavy atoms were fixed to the PDB coordinates. Finally, twenty modeling structures for the complex of M^pro and GC376 were obtained. Five optimization ranges were applied to each chain (A and B) and conformer of the ligands (B1S and K36). FMO calculations at the MP2/6-31G* level of these complexes were performed using the ABINIT-MP program [10]. These structures were divided into small fragments, each consisting of an amino acid residue, ligand, and a water molecule. As the ligand was covalently bound to Cys145, fragmentation was performed between them. (Figure 3). The structural stability of the entire system was evaluated using total energy, and protein–ligand interactions were evaluated with interfragment interaction energy (IFIE) and pair interaction energy decomposition analysis (PIEDA) [11]. IFIE comprises electrostatic (ES), exchange-repulsion (EX), charge-transfer with higher order mixed terms (CT+mix), and dispersion (DI) contributions as analyzed with PIEDA. The results of these FMO calculations were registered in FMODB [6]. The experimental and calculation data, such as occupancy and temperature factor, structurally optimized conditions, and FMODB ID, are summarized in Table 1. The temperature factors in the X-ray crystal structure, taken from Table 1, were examined to determine whether the A or B chain is a reasonable structure to use for detailed total energy and IFIE analyses. The complex between M^pro and GC376 of the A chain has lower maximum, minimum, and average temperature factors than that of the B chain; the complex of the A chain has less atomic fluctuations as compared with that of the B chain. Therefore, using total energy and IFIE, we analyzed the FMO calculations for the A chain to reveal the detailed structural stability of the molecular system and molecular recognition mechanism between M^pro and GC376.

Figure 2. Ligand and multiple conformer positions and occupancies in the A chain of the complex between M^pro and GC376 (PDB ID: 7CB7)

Ligands K36 and B1S are illustrated with light blue and pink CPK models, respectively. Arg60 and Arg298 conformers, A and B, are illustrated with orange and green CPK models, respectively. An entire M^pro structure was represented by ribbon model.

Figure 3. Fragmentation between Cys145 and ligand in the case of B1S (A) and K36 (B)

The SG sulfur atom of Cys145 was split as bond attached atom (BAA) and the C21 carbon atom of the ligand as bond detached atom (BDA).

Table 1. List of occupancies, temperature factors, and FMODB IDs for several complexes regarding chain IDs, ligand conformers, and conditions of structural optimizationsThe temperature factors show the maximum, minimum, and average of the heavy atoms on the ligand. The results of FMO calculations were registered in FMODB regarding FMODB IDs.

	Occupancies of ligands
Chain	^aA			^aB
Ligand	^bB1S	^bK36		^bB1S		^bK36
Occupancy	0.53	0.47		0.30		0.70
	Temperature factors of ligands
	^aA				^aB
Statistics of temperature factor	^bB1S		^bK36		^bB1S		^bK36
Max	28.73		28.78		35.03		35.36
Min	7.70		8.53		8.74		6.71
Average (mean)	16.46		16.55		18.89		18.89
	FMODB IDs
	^aA				^aB
Optimization condition	^bB1S		^bK36		^bB1S		^bK36
^cOptH	GN571		43JKN		66YRZ		LJN69
^dOptHSol	14NVZ		K3GZ3		22N5R		3QJ4L
^eOptHLSol	VRQZ1		Q1V4Y		7G6ZK		J3QR9
^fOptHSideLSol	YYV72		R58K8		M3N2Z		N1N7Q
^gOptAll	53JQZ		ZYVGN		9G822		8296Y

^aChain IDs, ^bLigand names (which identify the conformer), ^cStructure optimization of hydrogen atoms only, ^dStructure optimization of hydrogen atoms and solvent molecules, ^eStructure optimization of hydrogen atoms, solvent molecules, and ligand molecule, ^fStructure optimization of hydrogen atoms, solvent molecules, ligand molecule, and side chain of amino acid residues, ^gStructure optimization of all atoms.

2.2 Details for total energy and IFIE

We compared the total energy of the entire system and ligand binding energy evaluated by the sum of IFIEs between M^pro and the ligands (IFIE-sum) for each condition of structural optimization (Tables 2 and 3). Since the ligand and Cys145 are covalently bound, IFIE between Cys145 and the ligand was excluded from the IFIE-sum in the framework of conventional IFIE analysis by ABINIT-MP [10,12]. Table 2 shows that the total energy and its difference between B1S and K36, Δ(B1S–K36), changed toward convergence to an energetically stable state as the optimization range expanded. As an exception, the total energies of OptHSideLSol structures were more energetically stabilized than those of OptAll structures in FMO-based evaluation. The B1S complex exhibited higher structural stability of the entire system than the K36 complex. Since the ligand occupancy of B1S is higher than that of K36, the correlation between the structural stability of the molecular system and the ligand occupancy is consistent.

Table 2. Total energy of molecular system and their differential energies for each condition of structural optimization of B1S complex and K36 complex

Optimization condition	Total energy /hartree		Δ(B1S–K36) /hartree	Δ (B1S–K36) /kcal mol⁻¹
Optimization condition	B1S	K36	Δ(B1S–K36) /hartree	Δ (B1S–K36) /kcal mol⁻¹
^aOptH	−144716.986136	−144716.888094	−0.098042	−61.5
^bOptHSol	−144718.642741	−144718.564790	−0.077951	−48.9
^cOptHLSol	−144718.851410	−144718.810831	−0.040578	−25.5
^dOptHSideLSol	−144719.463004	−144719.403286	−0.059718	−37.5
^eOptAll	−144719.323931	−144719.303072	−0.020859	−13.1

^aStructure optimization of hydrogen atoms only, ^bStructure optimization of hydrogen atoms and solvent molecules, ^cStructure optimization of hydrogen atoms, solvent molecules, and ligand molecule, ^dStructure optimization of hydrogen atoms, solvent molecules, ligand molecule, and side chain of amino acid residues, ^eStructure optimization of all atoms.

Next, comparing the IFIE-sum of B1S and K36, we revealed that the IFIE difference, Δ(B1S–K36), became more negative with the expansion as the range of atoms for the optimization region was expanded. Here, the interaction energy components of IFIE-sum were also investigated by PIEDA and are presented in Table 3. This result indicated that compared with that of K36, the protein binding ability of the B1S complex gradually increased. In addition, Δ(B1S–K36) for optimization range including ligand (OptHLSol, OptHSideLSol, and OptAll) became more negative even for the ES and CT+mix components contributed to hydrogen bonds [6]. These results suggested that expanding the optimization range increases the difference in the protein binding strength between B1S and K36 with different occupancies. This difference may be attributed to the stability of hydrogen bonds.

Table 3. Sum of IFIE (IFIE-sum) and its energy decomposition by PIEDA and their difference energies between B1S and K36 for each condition of structural optimization of complexes

Optimization Condition	Ligand	IFIE-sum (Ser1–Ser144, Gly146–Gly302) /kcal mol⁻¹
Optimization Condition	Ligand	IFIE	ES	EX	CT+mix	DI
^aOptH	B1S	−116.4	−91.6	72.0	−31.1	−65.7
	K36	−117.1	−93.9	81.7	−35.6	−69.3
	Δ(B1S–K36)	0.6	2.3	−9.7	4.5	3.5
^bOptHSol	B1S	−114.9	−88.9	70.5	−30.9	−65.6
	K36	−109.7	−86.5	79.5	−34.3	−68.3
	Δ(B1S–K36)	−5.2	−2.4	−9.0	3.5	2.7
^cOptHLSol	B1S	−128.1	−97.5	65.4	−30.9	−65.1
	K36	−119.3	−85.7	63.2	−30.4	−66.3
	Δ(B1S–K36)	−8.8	−11.7	2.1	−0.5	1.2
^dOptHSideLSol	B1S	−126.4	−97.9	69.7	−31.6	−66.6
	K36	−115.5	−87.7	69.0	−29.4	−67.4
	Δ(B1S–K36)	−10.9	−10.2	0.8	−2.3	0.8
^eOptAll	B1S	−127.7	−103.7	80.0	−33.4	−70.6
	K36	−113.5	−90.2	76.6	−29.3	−70.6
	Δ(B1S–K36)	−14.3	−13.5	3.4	−4.1	0.0

2.3 Interactions between amino acid residues and ligands

To determine the key amino acid residues for ligand binding in the M^pro–GC376 complex, we analyzed IFIE and PIEDA for amino acid residues in M^pro whose IFIE with B1S or K36 is below −5.0 kcal mol⁻¹ (Figure 4A or 4B). Here, we investigated the complex structure data of B1S and K36 optimized under the OptAll condition, using FMODB entries 53JQZ and ZYVGN, respectively. As a result, His163, Met165, Glu166, Leu167, and Gln189 were identified as the amino acid residues with IFIEs less than −5 kcal mol^-1 with both B1S and K36 (Figure 4A). PIEDA results showed that the main components ES and CT+mix energies [6] stabilize both B1S and K36 through hydrogen bonds with these residues (Figure 4B). The IFIE with His41 showed a significant difference of −15.3 kcal mol⁻¹ between B1S (strong attraction) and K36 (weak attraction). IFIE with PIEDA between His41 and B1S showed ES energy (−16.9 kcal mol⁻¹) as the main component, while IFIE with PIEDA between His41 and K36 showed that ES energy had a positive value (1.2 kcal mol⁻¹), indicating that weak repulsion was occurred between them. In addition, CT+mix components of B1S and K36 were −3.9 kcal mol⁻¹ and −1.9 kcal mol⁻¹, respectively. This indicates that B1S has a stronger attraction towards M^pro than K36. The difference of His41 between B1S and K36 in ES and CT+mix components suggests that only B1S can form a hydrogen bond with His41. DI energy of His41 was close for both B1S and K36 (B1S-His41: −7.0 kcal mol⁻¹, K36-His41: −5.1 kcal mol⁻¹). DI energy was identified as the main component in the interaction between K36 and His41. This suggests that there is a dispersion interaction between B1S and His41 as well as between K36 and His41 to a comparable extent. Tables S1-S5 show the interaction energy between ligand and each amino acid residue of M^pro by IFIE and PIEDA for the OptH, OptHSol, OptHLSol, OptHSideLSol, and OptAll conditions; IFIE and ES components were stabilized by including of ligand in the optimization range.

Figure 4. Interaction energies between each amino acid residue of M^pro and GC376 by IFIE (A) and PIEDA (B)

FMO data for the complex between M^pro and B1S or K36 under OptAll condition were used (FMODB IDs: 53JQZ and ZYVGN). Fragment pairs between ligand and amino acid residues listed on the horizontal axis indicate with IFIE of −5.0 kcal mol^-1 or less.

Figure 5. Diagram of the interaction between the ligand and the surrounding amino acid residues using complexes under OptAll condition. 3D structure (A) and schematic diagram of ligand interaction (B) of the complex between Mpro and B1S (FMODB ID: 53JQZ); 3D structure (C) and schematic diagram of ligand interaction (D) of the complex between Mpro and K36 (FMODB ID: ZYVGN)

The red dotted lines and values represent hydrogen bonds and their distances (between heavy atoms and hydrogen); the green dotted lines and values represent CH/π interactions and their distances (between the center of gravity of the imidazole ring of His41 and the hydrogen of the ligand).

The interaction modes between these amino acid residues and the ligands are shown in Figure 5. The ligands form hydrogen bonds with heteroatoms and amides in the side and main chains of amino acid residues, respectively. CH/π interactions (2.9 Å) between a hydrogen atom of the alkyl chain in ligands and a π-orbital in the imidazole ring of His41 were also observed. These hydrogen bonds and the CH/π interaction were also confirmed through FMO calculations in the complex of M^pro with K36 in the previous paper [8]. In addition, the difference in interaction mode between B1S and K36 is because of the hydrogen bond (1.8 Å) between the oxygen atom of the hydroxy group in B1S and the NH of the imidazole ring in His41 (Figures 5A and 5B). These results indicated that B1S and His41 are stabilized by both the CH/π interaction and hydrogen bond. In contrast, the hydroxy group of K36 that identifies the R and S form was oriented opposite to His41 around the covalent bond between Cys145 and the ligand (Figure 5C). Then, hydrogen bonds between K36 and His41 were not observed. The orientation of the hydroxy group suggests that K36 forms a hydrogen bond with the NH in the main chain of Cys145. Thus, it would be more desirable to add an evaluation of the interaction between K36/B1S and Cys145 to discuss protein−ligand binding strength. However, as mentioned earlier, direct evaluation of IFIE between fragments covalently bound via BDA is not recommended within the framework of the conventional FMO calculation by ABINIT-MP [10,12]. Although not a desirable analysis, as a preliminary study, we confirmed that bare IFIE between the Cys145 and B1S/K36 (ca. −9.6×10³ kcal/mol) and its containing IFIE-sum between M^pro and ligand, Table S6 summarizes the difference between B1S and K36 of the followed values: the IFIE between Cys145 and ligand, the IFIE-sum between M^pro (Ser1–Gly302) and ligand. For all optimized structures, as shown in Table S6, the CT+mix components derived from the hydrogen bonding contribution [6] were a few kcal/mol stronger for K36 than for B1S, indicating an attractive interaction; However, the IFIE-sum of B1S was attractively stronger than that of K36 by several to ca. 20 kcal/mol among all optimized structures. This preliminary study suggests that even if we consider the IFIE between Cys145 and the ligand containing the hydrogen bond between Cys145 and K36, the conclusion that the protein–ligand binding energy strength is more stable for B1S than for K36 would not change. Here, note that the value of the bare IFIE between Cys145 and the ligand is much larger than the IFIE-sum of the ligand without the contribution of Cys145 as the ligand binding energy and the difference values of Table S6 may be in the error range. To address this issue, several advanced methods [11,13,14] can analyze the IFIE between covalently bound fragment pairs. Therefore, future studies will use these advanced methods to precisely assess covalent ligand binding ability.

3. Conclusion

In this study, we used FMO calculations on the complex of SARS-CoV-2 M^pro with its inhibitor, GC376, to investigate the effect of conformational differences in the stereo isomer of GC376 on the structural stability and the protein−ligand binding ability. A comparison of the total energy representing the structural stability between B1S and K36 reveals that the complex with higher ligand occupancy B1S is consistently lower than that with lower ligand occupancy K36 across structural optimization conditions. In addition, a comparison of IFIE-sum of the ligand and amino acid residues, representing the protein−ligand binding strength, between B1S and K36, revealed that the interaction tends to be stronger for the high-occupancy B1S than the low-occupancy K36 for most structural optimization conditions. IFIE-sum indicated that the difference between B1S and K36 in the protein binding strength increases as the optimization range expanded. The IFIE and PIEDA results for the optimization of all atoms indicated that both B1S and K36 bound to His163, Met165, Glu166, Leu167, and Gln189 by hydrogen bonds and His41 by CH/π interaction. The similarity of these interactions arises because the structures of B1S and K36 differ only at one stereocenter. On the other hand, only B1S can form a hydrogen bond with His41 via a hydroxy group of the stereo isomer. Since B1S has a higher protein binding strength than K36, a hydrogen bond with His41 contributes to the difference in the protein binding strength between B1S and K36. The results of this study, along with analysis of the interaction patterns between M^pro and various ligands registered in FMODB [8, 15], can be useful for structure-based drug design of the inhibitors of M^pro.

To our knowledge, this is the first study that investigates the structural stability and the protein–ligand binding strength of complexes containing ligands of different occupancies, using FMO calculations. In the all atoms optimization of the M^pro–GC376 of this study, structures modeled with the high occupancy conformers are energetically more stable. However, there were modeling structures by the hydrogen atom-only optimization, a common optimization condition, for which this relationship did not apply. Therefore, we could not clearly state a significant relationship between occupancy variance and structural stability or protein–ligand binding strength. In addition, selecting the conformer with the higher occupancy in modeling structures is not always right. For example, the modeled structure may crash when the conformer with the higher occupancy is selected [16]. It is preferable to confirm, by combining all conformers, whether a structure with varied occupancies yield high structural stability or strong protein–ligand binding. At least three M^pro complexes (PDB IDs: 6WTT, 7CB7, and 7T3Y) contain multiple ligands with different ligand occupancies in PDB. In addition, some protein structures of X-ray crystal structures registered in PDB have multiple conformers at amino acid residues. Moreover, as demonstrated by this study, interaction energy yielded through FMO calculations of protein–ligand binding strength could be affected by the choice of conformers. FMO-based analysis of more complexes containing multiple conformers may pave the way for future research examining a relationship between occupancies and structural stability or the protein–ligand binding ability.

Acknowledgements

This research was performed as part of the activities of the FMO drug design consortium (FMODD). The authors acknowledge Prof. Kaori Fukuzawa of Osaka University for her useful support of COVID-19 data collection. We also thank Prof. Yuji Mochizuki of Rikkyo University, Dr. Tatsuya Nakano of the National Institute of Health Sciences, and Dr. Yoshio Okiyama of Kobe University for providing ABINIT-MP and general discussions related to FMO calculations. This study was partially supported by Research Support Project for Life Science and Drug Discovery (Basis for Supporting Innovative Drug Discovery and Life Science Research (BINDS)) from AMED under Grant Number JP23ama121030. We also thank Editage ( www.editage.jp) for English language editing.

References

Corresponding author

Correction information

Funder information

1.Fund name: Research Support Project for Life Science and Drug Discovery (Basis for Supporting Innovative Drug Discovery and Life Science Research (BINDS)) from AMED

Register with J-STAGE for free!