Fragment molecular orbital calculations containing Mg2+ ions: PPlase domain of Cyclophilin G

Masayasu Fujii; Chiduru Watanabe; Kaori Fukuzawa; Shigenori Tanaka

doi:10.1273/cbij.22.55

Abstract

We investigated fragmentation methods around metals when performing fragment molecular orbital (FMO) calculations for metal-containing proteins, as well as appropriate structural preprocessing. The protein structure data was employed from the peptidyl prolyl cis/trans isomerase (PPlase) domain of human Cyclophilin G, which contains Mg²⁺ ions (PDBID: 2WFI). The results of the present study revealed three issues: First, it was better to contain Mg²⁺ ions in the same fragment as that for water molecules in the first hydration shell, which was revealed thorough PIEDA (pair interaction energy decomposition analysis) and atomic charge analysis obtained by FMO calculations. Second, while there were two conformers of Phe72 in the PPlase domain of Cyclophilin, we could determine the more appropriate conformer for computational analysis by comparing each energy component of PIEDA. Finally, we derived the optimal constraint conditions for the structural optimization of this molecular system in which the exchange repulsion (EX) component of PIEDA took relevant values without deforming the initial structure too much. These findings could be applied to FMO calculations for other proteins as well.

1. Molecular modeling and FMO calculation setup containing Mg²⁺ ions

The peptidyl prolyl cis/trans isomerase (PPlase) domain of human Cyclophilin G containing Mg²⁺ ions (PDBID: 2WFI) [1] (Figure 1) was used for FMO calculation. This data was obtained by high-resolution X-ray crystallography with the resolution of 0.75 Å. Therefore, since the coordinate data for this structure was of very high-resolution and reliable, we examined the conditions for proper structure optimization without deforming it too much from the initial structure. The structures were prepared using the molecular operating environment (MOE, v2020.09; Chemical Computing Group Inc., Montreal, QC, Canada) with the following procedure. First, five water molecules with a temperature factor greater than 59.56, which was the maximum temperature factor of this protein, were removed because they were supposed to be structurally inaccurate (Figure S1). Then, there remained 427 water molecules for the FMO calculations. Since the closest interatomic distance between the five removed water molecules and the two Mg²⁺ ions was 13.57 Å, the following study on the Mg²⁺ ion was unaffected by the removal of the five water molecules. If there were multiple conformers in a PDB structure, the conformer with the highest occupancy was used first. The missing residues were formally complemented by the “Structure Preparation” function (built-in functions in MOE) and hydrogen atoms were added using the “Protonate3D” function at pH 7.0, while this PDB structure had no missing residues. The residues at the N- and C-termini were treated as electrically neutral with NH₂ and COOH, respectively. Subsequently, only the positions of the complementary atoms and water molecules were energetically optimized using the Amber10:EHT force field.

All hydrogen atoms and the oxygen atoms of water molecules to be optimized were constrained to the starting coordinates by the harmonic potential with the force constant ( 3 2 ) k T σ 2 , where the heavy atoms of the other parts were fixed at the PDB coordinates. Here, σ is a value specified by Deviation (in Angstrom). Subsequently, all heavy atoms were fixed, and only the hydrogen atoms were structurally optimized without constraint. A force constant of zero σ means that the atom is fixed, and as the value increases, the force constant becomes smaller. Since the user could adjust the force constant by changing the value of σ , we performed fragment molecular orbital (FMO) calculations for each optimized structure by changing σ from 0.5 to 1.0 Å in increments of 0.1 Å. Due to the high resolution of the X-ray data, the results of these calculations were then compared to determine the value of σ at which the structure could be properly optimized without deforming it too much from the initial structure.

Figure 1. (a) Structure of Cyclophilin G (PDB ID: 2WFI). The green sphere refers to Mg²⁺ ion. (b) and (c) Two views of molecular structure around Mg²⁺ ion

The optimization condition is given by the constraint with 𝛔 = 0.5 Å using conformer B of Phe72.

In order to confirm whether the structural optimization was properly performed, we investigated fragment pairs with the exchange repulsion (EX) component being more significant than 30 kcal/mol and larger than the electrostatic (ES) component in absolute value as unstable substructures by the interaction energy analysis using Pair Interaction Energy Decomposition Analysis (PIEDA) [2]. Since the EX component represents the exchange repulsion term, its too large value indicates that the atoms are closer together than the allowed range. As a rule of thumb, if a dimer pair has an EX value of more significant than 30 kcal/mol and the absolute value of EX is greater than that for ES, the dimer pair may be causing steric hindrance. Here, the FMODB data [3] were used to reconsider the validity of this rule. In Figure S2, we focus on the interaction energy distribution of fragment pairs holding a single hydrogen bond as a reference interaction. From Figure S2(a), it is considered that the IFIE begins to become unstable at distances closer than 2.45 Å between the heavy atoms forming the hydrogen bond. Moreover, we presume that the maximum value of the box-plot of IFIE (MP2) turns into a repulsion, and structures containing steric hindrance begin to appear. Therefore, the median of box-plot on EX data at 2.45 Å, as a value of 30 kcal/mol (Figure S2(b)), is recommended as a threshold value for the structural check because the EX term mainly causes repulsive interaction.

The software ABINIT-MP [4] was used to perform FMO calculations with the prepared structures, where the MP2/6-31G* level was employed to treat the electron-correlation effects. Because there were water molecules around Mg²⁺ ions (Figure 1(b) and 1(c)), we examined the calculation results for two fragmentation schemes, treating Mg²⁺ ion as a single fragment or merging it with the water molecules in the first hydration shell, to investigate whether the difference in the fragmentation scheme affected the calculated interaction energies [5]. To this end, the following two fragmentation methods were attempted: No merging: amino acid residues, Mg²⁺ ions and water molecules were treated as one fragment; Merging: one fragment of amino acid residues, Mg²⁺ ions and water molecules within 2.2 Å of Mg²⁺ ions were merged into one fragment and the other water molecules were treated as one fragment each. The FMO calculation results were registered in FMO database (https://drugdesign.riken.jp/FMODB/) [3, 6], and their entry IDs (FMODB IDs) are listed in Table S1.

2. Results and Discussion

First, to examine the validation of fragmentation around Mg²⁺ ion, we extracted fragment pairs with EX component larger than 30 kcal/mol, when treating Mg²⁺ ion as a single fragment or when merging it with the water molecules in the first hydration shell (Table 1). Here, the constraint condition for the structure used for verification was σ = 0.5 Å. The EX and ES components of several fragment pairs between Mg²⁺ ion and the water molecules in the first hydration shell were large in the former case. The fragments between Mg²⁺ ion and selected water molecules are found to be close together at a distance of ca. 1.92 Å, while the van der Waals radii of oxygen and magnesium atoms are ca. 1.40 and 1.60 Å, respectively. As the IFIE and its energy components alone were not sufficient to determine the validity of the magnesium ion treatment in the FMO context, an atomic charge analysis of magnesium ion was also performed and shown in Table 2. In the absence of Mg²⁺ fragment merging, the atomic charges of MG1001 and MG1002 are 0.949 and 0.169, respectively, by the Mulliken population analysis; 1.613 and 1.207, respectively, by the natural population analysis. In the case of merging Mg²⁺ fragment with water molecules, the atomic charges of MG1001 and MG1002 are 1.454 and 1.359, respectively, by the Mulliken population analysis; 1.871 and 1.765, respectively, by the natural population analysis. Without merging the Mg²⁺ fragment with water molecules, the atomic charges of magnesium ion are less than 1.0 for Mulliken population analysis, suggesting that the electronic state is not treated relevantly. By merging the Mg²⁺ fragment with water molecules, on the other hand, the results show that the atomic charges of magnesium ion became closer to +2 for both Mulliken and natural population analyses (ca. 1.4 and 1.8). It is noted here that the charge transfer from water molecule to magnesium ion is greater than 0.15 (Table 1). Thus, it was better to merge Mg²⁺ fragment with water molecules in the first hydration shell in the FMO calculation. Therefore, we adopted the method of “merging water” for fragmentation in the subsequent calculations. For the other fragment pairs in Table 1, it can also be seen that there are several pairs such as PHE72 and HOH2090 pair where the interaction energy components satisfy the conditions of EX > 30 kcal/mol and |EX| > |ES|.

Table 1. Fragment pairs (I, J) with EX larger than 30 kcal/mol under different merging conditions in PIEDA

The optimization condition is given by the constraint with σ = 0.5 Å using conformer B of Phe72. All the PIEDA energies are shown in units of kcal/mol. “Dist” and “q(I=>J)” mean the distance between nearest neighboring fragments and the charge transfer from fragment I to fragment J, respectively.

merge	fragment pair (I, J)	Dist	total	ES	EX	CT+mix	DI	q(I=>J)
No merge	SER53, GLU88	1.763	−33.03	−48.28	33.58	−11.06	−7.27	−0.088
	ASN18, HOH2043	1.546	−15.21	−31.44	30.80	−11.56	−3.00	0.068
	PHE72, HOH2090	1.674	−0.54	−29.13	51.45	−15.61	−7.25	0.006
	HOH2142, HOH2144	1.514	−7.16	−27.42	31.05	−9.01	−1.78	0.057
	GLU81, HOH2246	1.524	−30.23	−44.53	30.67	−12.59	−3.78	0.079
	GLU96, HOH2266	1.593	−23.79	−41.64	32.06	−10.51	−3.71	0.079
	GLU96, HOH2267	1.530	−13.58	−30.86	30.30	−9.98	−3.05	0.060
	MG1001, HOH2268	1.914	−118.03	−138.68	33.42	−8.68	−4.10	−0.154
	MG1001, HOH2322	1.918	−118.58	−138.66	33.11	−8.83	−4.21	−0.163
	MG1001, HOH2323	1.920	−116.41	−136.47	32.76	−8.62	−4.09	−0.158
	GLU155, HOH2379	1.525	−14.86	−33.44	34.76	−12.17	−4.01	0.057
	HOH2423, HOH2432	1.551	−3.65	−23.62	32.83	−11.08	−1.78	−0.062
Merge water	SER53, GLU88	1.763	−33.02	−48.27	33.58	−11.06	−7.27	−0.088
	ASN18, HOH2043	1.546	−15.21	−31.44	30.80	−11.56	−3.00	0.068
	PHE72, HOH2090	1.674	−0.54	−29.13	51.45	−15.60	−7.25	0.006
	HOH2142, HOH2144	1.514	−7.16	−27.42	31.05	−9.01	−1.78	0.057
	GLU81, HOH2246	1.524	−30.24	−44.53	30.67	−12.59	−3.78	0.079
	GLU96, HOH2266	1.593	−23.55	−41.53	32.11	−10.44	−3.68	0.078
	GLU96, HOH2267	1.530	−13.43	−30.74	30.32	−9.97	−3.04	0.059
	GLU155, HOH2379	1.525	−14.86	−33.43	34.76	−12.17	−4.01	0.057
	HOH2423, HOH2432	1.551	−3.65	−23.62	32.83	−11.08	−1.78	−0.062

Table 2. Atomic charge and population of Mg²⁺ ion evaluated by Mulliken and natural population analyses

The optimization condition is given by the constraint with σ = 0.5 Å using conformer B of Phe72. The atomic charge and population are shown in units of e.

		Mulliken atomic population		Natural atomic population
merge	Atom	Charge	Population	Charge	Population
No merge	MG1001	0.949	11.051	1.613	10.387
No merge	MG1002	0.169	11.831	1.207	10.793
Merge water	MG1001	1.454	10.546	1.871	10.129
Merge water	MG1002	1.359	10.641	1.765	10.235

Table 3. Fragment pairs (I, J) with EX larger than 30 kcal/mol and larger than ES in absolute value under different constraint conditions in PIEDA

The optimization conditions are given by the several constraints using conformer B of Phe72. All the PIEDA energies are shown in units of kcal/mol. “Dist” and “q(I=>J)” mean the distance between nearest neighboring fragments and the charge transfer from fragment I to fragment J, respectively.

𝛔 (Å)	fragment pair (I, J)	Dist	total	ES	EX	CT+mix	DI	q(I=>J)
0.5	PHE72, HOH2090	1.674	−0.54	−29.13	51.45	−15.60	−7.25	0.006
	HOH2142, HOH2144	1.514	−7.16	−27.42	31.05	−9.01	−1.78	0.057
	GLU155, HOH2379	1.525	−14.86	−33.43	34.76	−12.17	−4.01	0.057
	HOH2423,HOH2432	1.551	−3.65	−23.62	32.83	−11.08	−1.78	−0.062
0.6	PHE72, HOH2090	1.725	−0.92	−21.57	37.78	−10.59	−6.56	0.012
0.6	GLU155, HOH2379	1.517	−14.61	−33.85	35.50	−12.35	−3.97	0.057
0.7	PHE72, HOH2090	1.761	−1.71	−18.59	31.72	−8.72	−6.14	0.016
0.7	GLU155, HOH2379	1.512	−14.53	−34.18	36.01	−12.49	−3.94	0.058
0.8	GLU155, HOH2379	1.508	−14.46	−34.36	36.33	−12.57	−3.92	0.058
0.9	GLU155, HOH2379	1.505	−14.52	−34.65	36.65	−12.67	−3.91	0.058
1.0	GLU155, HOH2379	1.501	−14.40	−34.85	37.06	−12.73	−3.88	0.059

Next, Table 3 showed the extraction of fragment pairs with the EX component being larger than 30 kcal/mol and larger than the ES component in absolute value for each constraint condition. The fragmentation of Mg²⁺ ion and the water molecules in the first hydration shell was performed based on the optimized structural coordinates. It was found that the larger the value of σ , the fewer residues were extracted. However, for the pairs of HOH2090 and Phe72, and HOH2379 and Glu155, changing the constraint conditions did not eliminate the tendency for EX to be larger than ES.

The structure of Phe72 was investigated, and it was found that there are two conformers of Phe72 in the PDB structure 2WFI (Figure 2). The structures examined in Tables 1 and 3 were taken from the conformer structure with the highest occupancy rate in the presence of multiple conformers, and the conformer B structure (occupancy rate 0.52) was used for Phe72. Therefore, we performed the structural modeling again using the conformer A structure (occupancy 0.48) for Phe72 and investigated the EX components of Phe72 and HOH2090 under different constraint conditions. The results showed that the EX values of Phe72 and HOH2090 were less than 30 kcal/mol even when the constraint condition was σ = 0.5 Å (Table 4). The data show that the conformer with the highest occupancy is not always the optimal one. These results indicated that the conformer A structure of Phe72 was preferable for 2WFI, and the optimal constraint was σ = 0.5 Å.

On the other hand, we think the structure of Glu155 had no problem (Figure 3). Because the ES component could almost cancel out the EX component and the IFIE value (ca. −14 kcal/mol) was similar to the standard hydrogen-bond energy value seen in FMODB (−13.7 kcal/mol) [3], no adjustments were made for the fragment pair of Glu155 and HOH2379.

Figure 2.Figure 2. Location of Phe72 and HOH2090 for each conformer. Green color for conformer A (occupancy 0.4 Location of Phe72 and HOH2090 for each conformer. Green color for conformer A (occupancy 0.48) and red color for conformer B (occupancy 0.52)

The optimization condition is given by the constraint with σ = 0.5 Å.

Table 4. Each energy component of PIEDA between Phe72 (conformer A) and HOH2090 under different constraint conditions

All the PIEDA energies are shown in units of kcal/mol. “Dist” and “q(I=>J)” mean the distance between nearest neighboring fragments and the charge transfer from Phe72 to HOH2090, respectively.

𝛔 (Å)	Dist	total	ES	EX	CT+mix	DI	q(I=>J)
0.5	2.397	−2.96	−0.90	1.00	−1.54	−1.50	−0.015
0.6	2.401	−2.75	−0.85	0.93	−1.41	−1.41	−0.013
0.7	2.466	−3.21	−1.15	1.04	−1.57	−1.51	−0.015
0.8	2.586	−1.69	−0.33	1.05	−1.02	−1.37	−0.013
0.9	2.597	−1.73	−0.35	1.00	−1.01	−1.36	−0.013
1.0	2.603	−1.72	−0.33	0.96	−1.00	−1.33	−0.012

Figure 3. Location of Glu155 and HOH2379 The optimization condition is given by the constraint with σ = 0.5 Å using conformer B of Phe72.

3. Conclusion

We investigated the fragmentation method around the metal and the appropriate structural preprocessing scheme when performing FMO calculations for the PPlase domain of human Cyclophilin G containing Mg²⁺ ions (PDB ID: 2WFI). For the fragmentation method around the metal, the effects of the different fragmentation methods on the interaction energies and atomic charges were investigated in one case where Mg²⁺ ion was treated as a single fragment and in another case where Mg²⁺ ion and the water molecules in the first hydration shell were merged. In the former case, the EX component of fragment pairs between Mg²⁺ ion and the water molecules in the first hydration shell was larger, so the latter seems to be better. However, the IFIE and its energy components alone were insufficient to determine the validity of the treatment for magnesium ion in fragmentation. In comparison, atomic charges of Mg²⁺ ion were more reasonable for both Mulliken and natural population analyses in merging the Mg²⁺ fragment with water molecules than without merging fragments between Mg²⁺ ion and water molecules. Therefore, the magnesium fragment should be treated in this system by merging it with the water molecules in the first hydration shell. For the appropriate structural preprocessing, by comparing the EX components of PIEDA obtained by FMO calculations, the conformer A structure was preferable for Phe72, which had two conformers in the PDB structure. Then, by comparing the calculated results of the optimized structures with different constraints, we identified the optimized structure so that the EX component of PIEDA was appropriate without deforming the initial structure too much. As a result, the optimal constraint for this molecular system was found to be σ = 0.5 Å. These findings can be applied to FMO calculations for other proteins as well.

Acknowledgements

This research was performed in the activities of the FMO drug design consortium (FMODD). The authors thank Prof. Yuji Mochizuki of Rikkyo University, and Dr. Tatsuya Nakano and Dr. Yoshio Okiyama of the National Institute of Health Sciences for providing the Fugaku version of ABINIT-MP and general discussions related to FMO calculations. We also acknowledge Dr. Teruki Honma, Dr. Daisuke Takaya and Dr. Kikuko Kamisaka of RIKEN for useful comments and technical support of FMODB registration. The FMO calculations were performed using Fugaku of RIKEN (project ID: hp210130). This study was partially supported by the Platform Project for Supporting Drug Discovery and Life Science Research (Basis for Supporting Innovative Drug Discovery and Life Science Research (BINDS) from the Japan Agency for Medical Research and Development (AMED) (grant number JP21am0101113) as well as JSPS Kakenhi (JP17H06353, JP18K03825, and JP21K06098) and MEXT Quantum Leap Flagship Program (MEXT QLEAP) Grant Number JPMXS0120330644. CW acknowledges JST PRESTO grant (JPMJPR18GD).

References

責任著者(Corresponding author)

J-STAGEへの登録はこちら（無料）