Free-Energy Landscape Analysis of Protein-Ligand Binding: The Case of Human Glutathione Transferase A1

Abstract

Glutathione transferases (GSTs) are a superfamily of enzymes which have in common the ability to catalyze the nucleophilic addition of the thiol group of reduced glutathione (GSH) onto electrophilic and hydrophobic substrates. This conjugation reaction, which occurs spontaneously but is dramatically accelerated by the enzyme, protects cells against damages caused by harmful molecules. With some exceptions, GSTs are catalytically active as homodimers, with monomers generally constituted of 200 to 250 residues organized into two subdomains. The first is the N-terminal subdomain, which contains an active site named G site, where GSH is hosted in catalytic conformation and which is generally highly conserved among GSTs. The second subdomain, hydrophobic, which binds the substrate counterpart (H site), can vary from one GST to another, resulting in structures able to recognize different substrates. In the present work, we performed all-atom molecular dynamics simulations in explicit solvent of human GSTA1 in its APO form, bound to GSH ligand and bound to GS-conjugated ligand. From MD, two probes were analyzed to (i) decipher the local conformational changes induced by the presence of the ligand and (ii) map the communication pathways involved in the ligand-binding process. These two local probes are, first, coarse-grained angles $(\theta$,$\gamma )$, representing the local conformation of the protein main chain and, second, dihedral angles $\chi$ representing the local conformation of the amino-acid side chains. From the local probes time series, effective free-energy landscapes along the amino-acid sequence were analyzed and compared between the three different forms of GSTA1. This methodology allowed us to extract a network of 33 key residues, some of them being located in the experimentally well-known binding sites G and H of GSTA1 and others being located as far as 30Å from the original binding sites. Finally, the collective motions associated with the network of key residues were established, showing a strong dynamical coupling between residues Gly14-Arg15 and Gln54-Val55, both in the same binding site (intrasite) but also between binding sites of each monomer (intersites).

1. Introduction

Proteins are biological macromolecules that perform a large variety of functions in living cells comprising biochemical (enzymes), structural, mechanical, and signaling functions. To perform their functions, proteins interact with small molecules referred to as ligands, which are able to bind to a protein with high affinity and specificity [1]. These protein/ligand interactions are crucial in biology, particularly in the context of drug design [2]. Since proteins interact with a broad range of drugs, it is of particular interest to study the mechanisms of binding of ligands to proteins and its impact on the structural dynamics to gain insights into (i) phenomena involved in the biological process and related to diseases [3] (misfolding, aggregation), and (ii) discovery, design, and development of new drugs [4]. The experimental structural data (e.g., X-ray crystallography, NMR, or cryo-EM) provide key structural information of the ligand-bound and ligand-unbound (APO) proteins [5]. Nevertheless, the static information is not always sufficient for understanding protein–ligand binding mechanisms, especially when pockets are highly flexible and contain several binding sites. Therefore, molecular dynamics (MD) is a powerful tool that provides a description of the dynamics and structures of protein–ligand systems with a high spatial and temporal resolution.

Glutathione transferases (GSTs) belong to a ubiquitous superfamily of enzymes that metabolize a broad range of reactive toxic compounds by catalyzing the conjugation of reduced tripeptide glutathione ($\gamma$-Glu-Cys-Gly; named GSH) to the electrophilic center of a second substrate [6,7,8], the reactivity of GSH being due to the thiol group SH of the cysteine residue. The conjugation reaction occurs spontaneously but GST accelerates it dramatically. This process of detoxification protects cells against damages caused by both exogenous and endogenous molecules. GSTs were first discovered in liver cells [9,10], and since then, they have been found to exhibit ligand-binding properties for a large variety of compounds, which are not always their enzymatic substrates [11]. Therefore, GSTs participate in diverse biological processes, making them multifunctional proteins. Moreover, GSTs are classified into three families according to their location in the cell: cytosolic, mitochondrial, and microsomal, which is not evolutively related to the two other classes [12]. First-discovered and most-abundant cytosolic GSTs are divided into 13 classes based on homology of their sequences. Members of the same cytosolic class have at least 40% of sequence identity, while members of different classes must have at most 25% of sequence identity. Even if they present a low homology with the cytosolic GST, mitochondrial GSTs can be considered as a particular class of GSTs (Kappa). Humans possess GST members in seven different classes [13], particularly the Alpha-class (cytosolic GSTs), in which the GSTA1 protein, which is the protein of interest in the present work, belongs. An alternative classification is possible on the basis of the residue located in the G site and which favors the activation of GSH (deprotonation of GSH plus reduction of pKa): the Cys-GSTs, whose structures are very similar to the ancestral precursor of all GSTs, the Ser-GSTs (Delta, Theta, Zeta, and Phi classes, including also Nu-GST activated by a threonine), and the Tyr-GSTs (Alpha, Pi, Mu, and Sigma classes); this latter subfamily comprises the more recently evolved GSTs.

With some exceptions, GSTs are catalytically active as dimers. The GSTA1 dimer is stabilized by a “lock and key” motif consisting of two key residues (Met51 and Phe52) fitting into a hydrophobic cavity of the other dimer [14]. Human GSTA1 dimer revealed negative cooperativity properties depending on the substrates. It was proposed that this negative cooperativity allows the self-preservation of their functions [15,16]. Additionally, there is no clear evidence that monomeric forms of GSTA1 are active as well as folded [17]. GST monomers are, in general, made of 200 to 250 residues with a molecular weight generally comprising between 25 and 30 kDa [18,19] and are organized into two subdomains (Figure 1A): the typically 80-residue-long N-terminal subdomain (I) has the typical fold of thioredoxin. It contains a first active site where GSH is hosted in catalytic conformation, named G site (Figure 1B). The thioredoxin fold is composed of a characteristic sequence of $\alpha$-helices and $\beta$-strands encountered in the thioredoxin protein family, i.e., ${\beta }_1-{\alpha }_1-{\beta }_2-{\alpha }_2-{\beta }_3-{\beta }_4-{\alpha }_3$, which is characteristic of enzymes dealing with gluthatione, such as glutaredoxin or glutathione peroxydase [20]. In addition, it has been shown that the region ${\beta }_3-{\beta }_4-{\alpha }_3$ is well conserved among GSTs and enables GSH recognition by the enzyme [7]. Residues forming the G site are generally conserved among GSTs [19]. A noticeable exception is the residue of ${\alpha }_1$-helix which interacts with the sulfur atom of GSH and which can be a cysteine, a serine, or a tyrosin (as is the case for human GSTA1). Depending on the nature of this residue, GSTs develop slightly different catalytic functions and target a different range of substrates [21]. The second subdomain (II) is all-helical and contains the H site which binds the substrate. The number of helices varies between four and seven among GSTs. Subdomain (II) is hydrophobic, and therefore attractive for hydrophobic molecules. Together with the G site, the combined architecture of GST monomers is adapted to bind GSH to hydrophobic substrates [22]. Contrary to the G site, the residues of the H site strongly vary from one GST to the others, resulting in H sites of different natures and the ability to recognize different substrates [23].

Hereafter, we focus our interest on the structure of human GSTA1 (hGSTA1), which is a homodimer with each monomer made of 222 residues [24]. Among the 29 experimental structures of hGSTA1 available in the Protein Data Bank, 21 of them contain at least one ligand, which is GSH, substrates (R), or GS-conjugates (GS-R). For each structure, a sequence analysis is provided, featuring the residues of hGSTA1 which are in contact with the ligand. Collecting these data (Figure S1), we determined the list of residues bound to GSH and therefore involved in the G site of hGSTA1: Tyr9, Arg15, Arg45, Gln54, Val55, Pro56, Gln67, Thr68, Asp101, Arg131, and Phe220 (Asp101 and Arg131 belong to the opposite monomer; see Figure 1B). We performed the same analysis to identify residues bound to GS-conjugates and therefore involved in the H site of hGSTA1: Phe10, Gly14, Ser18, Arg69, Leu72, Ile96, Glu97, Ala100, Ile106, Leu107, Leu108, Val111, His159, Met208, Leu213, and Phe222. From one experimental X-ray resolved structure, we performed all-atom classical molecular dynamics (MD) in explicit solvent of hGSTA1 enzyme (see Supplementary Materials for details about the MD protocol) in three different ligand-binding forms of the conjugation reaction cycle shown in Figure 1C: (i) its APO form, when there is no compound bound to GST, (ii) its GSH form, when the glutathione is bound to the G-site, and (iii) its GS-R form, when the glutathione-S-conjugated substrate is bound to the G and H sites, the substrate considered here being the 1-chloro-2,4-dinitrobenzene (CDNB). Another existing form of the hGSTA1 enzyme is not considered in the present work, i.e., the R form (substrate-bound). Indeed, the binding of the gluthatione ligand and of the substrate is not sequential. The substrate can bind first to GST and then the GSH, or the opposite way, as presented in Figure 1C.

Overall, the present study aims at deciphering the local conformational changes and shifts of populations of different local minima of the main chain and side chains of hGSTA1 enzyme in its three different ligand-binding forms. Predicting ligand-induced local free-energy changes is relevant both for understanding sequence–structure–function relationships in enzymes and also for structure-based drug design. From MD, we explore local conformational changes using internal coordinates, i.e., coarse-grained angles of the main chain and dihedral angles of the side chains. Dissimilarities between effective free-energy landscapes of the internal coordinates were established in order to identify the network of residues involved in the ligand-binding process of hGSTA1. Finally, the coupling between internal coordinates was revealed by analyzing the correlation of their motions using the principal component analysis, allowing the definition of collective motions to which these degrees of freedom contribute the most.

2. Materials and Methods

2.1. Internal Angles as Local Probes of the Protein Main and Side-Chain Conformational Changes

Local conformational changes of the main chain of human GSTA1 were analyzed using coarse-grained angles (CGAs) $(\theta ,\gamma$), which can be represented by a unit vector in spherical coordinates, where $\gamma$ is the azimuth angle and $\theta$ the polar angle [25]. For a residue i along the amino-acid sequence (see Figure 2A), ${\theta }_i$ is the bond angle formed by the virtual bonds joining three successive C${}^{\alpha }$ atoms ($i-1$, i and $i+1$) and ${\gamma }_i$ is the dihedral angle formed by the virtual bonds joining four successive C${}^{\alpha }$ atoms ($i-1$, i, $i+1$ and $i+2$). The first pair of CGA $(\theta ,\gamma )$ along the sequence is $({\theta }_2,{\gamma }_2)$ and the last one is $({\theta }_{N-2},{\gamma }_{N-2})$, where N is the total number of residues. The convention for $\gamma$ angles is the following: each angle varies between −180${}^{\circ }$ and +180${}^{\circ }$, with $\gamma =0^{\circ }$ being chosen when C${}_{i-1}^{\alpha }$ is cis to C${}_{i+2}^{\alpha }$ and the clockwise rotation of C${}_{i+1}^{\alpha }-C_{i+2}^{\alpha }$ is positive when looking from C${}_i^{\alpha }$ to C${}_{i+1}^{\alpha }$. Because the length of the $C^{\alpha }-C^{\alpha }$ virtual bond between two consecutive residues is nearly constant, the main-chain conformation is entirely described by the main-chain bond angles $\theta$ and the main-chain torsional angles $\gamma$ (Figure 2A). These CGAs $(\theta ,\gamma )$, which represent the torsion and curvature of the protein main chain and form a complete set of order parameters for protein folding [26], are part of coarse-grained protein models [27] and were used to analyze large conformational changes of proteins [28], protein folding, and dynamics in all-atom simulations [25,29,30] and conformational ensemble of intrinsically disordered proteins [31].

Moreover, local conformational changes of the side chains of human GSTA1 were analyzed using side-chain dihedral angles (SCAs) ${\chi }^k$. SCAs capture the rotation around $C-C$ or $C-N$ bonds of the side chain from its $C^{\beta }$ to its extremity. Each SCA is built from the coordinates of four successive atoms along the side chain of an amino acid. First of all, dihedral angle ${\chi }^1$ is made of $N-C^{\alpha }-C^{\beta }-X$, where X depends on the amino acid. X can be $S^{\gamma }$ (Cys), $O^{\gamma }$ (Ser, Thr), or $C^{\gamma }$ for all the other amino acids except Gly and Ala, for which ${\chi }^1$ dihedral angle is not defined. Second, ${\chi }^2$ corresponds to the rotation around the bond $C^{\beta }-C^{\gamma }$ and is not defined for amino acids Gly, Ala, Cys, Ser, Thr, and Val. Third, ${\chi }^3$, corresponding to the rotation around the bond $C^{\gamma }-C^{\delta }$ or $C^{\gamma }-S^{\delta }$ (Met), is only defined for amino acids Gln, Glu, Met, Arg, and Lys. Finally, ${\chi }^4$ are only defined for Arg and Lys amino acids and correspond to the rotation around the bond $C^{\delta }-C^{ϵ}$ (Lys) or $C^{\delta }-N^{ϵ}$ (Arg), and ${\chi }^5$ is only defined for Arg and corresponds to the rotation of the side chain around the bond $N^{ϵ}-C^{\zeta }$ (Figure 2A).

2.2. Free-Energy Surface, Free-Energy Profile, and Similarity Index

Effective 2D Free-Energy Surfaces (FESs) $V({\theta }_i,{\gamma }_i)$ were computed for each pair of CGA $({\theta }_i,{\gamma }_i)$ by using

$$V({\theta }_i,{\gamma }_i)=-k_{B}T\log P({\theta }_i,{\gamma }_i)$$

(1)

where $k_B$ is the Boltzmann constant, T is the temperature, and $P({\theta }_i,{\gamma }_i)$ is the probability density function (PDF) of CGA pair $({\theta }_i,{\gamma }_i)$. Two-dimensional PDFs were computed on a sphere, using CGA time series (Figure 2B) extracted from concatenated MD trajectories for each hGSTA1 binding-state: APO (2 runs of 900 ns), GSH (2 runs of 900 ns), and GS-R (1 run of 1200 ns). An example of free-energy surface for CGA pair $({\theta }_{84},{\gamma }_{84})$ (Gly83-Lys84-Asp85-Ile86) is shown in Figure 2C.

Then, to quantify the modifications between two FESs of identical CGA pair $({\theta }_i,{\gamma }_i)$ due to different ligand-binding form of hGSTA1, i.e., APO vs. GSH, GSH vs. GS-R, or GS-R vs. APO, as shown in Figure 1C, we computed the similarity index H between their associated 2D PDFs using

$$H_{1|2}({\theta }_i,{\gamma }_i)={\displaystyle \frac{{\displaystyle 2{\int }_0^{\pi }\sin {\theta }_{i}d{\theta }_i{\int }_{-\pi }^{+\pi }{P}_1({\theta }_i,{\gamma }_i)P_2({\theta }_i,{\gamma }_i)d{\gamma }_i}}{{\displaystyle {\int }_0^{\pi }\sin {\theta }_{i}d{\theta }_i{\int }_{-\pi }^{+\pi }{P}_1{({\theta }_i,{\gamma }_i)}^{2}d{\gamma }_i+{\int }_0^{\pi }\sin {\theta }_{i}d{\theta }_i{\int }_{-\pi }^{+\pi }{P}_2{({\theta }_i,{\gamma }_i)}^{2}d{\gamma }_i}}}$$

(2)

where $P_1({\theta }_i,{\gamma }_i)$ is the PDF of CGA pair $({\theta }_i,{\gamma }_i)$ in binding form 1 and $P_2({\theta }_i,{\gamma }_i)$ is the PDF of CGA pair $({\theta }_i,{\gamma }_i)$ in binding form 2. The similarity index H varies between 0 (dissimilar) and 1 (identical). We consider similarity between 2 FESs to be (i) large if $H>0.70$, (ii) moderate if $0.30\le H<0.70$, and (iii) low if $H<0.3$, as per a previous work [28].

Similarly, effective 1D free-energy profiles (FEPs) $V\left({\chi }_i^k\right)$ were computed for each SCA ${\chi }_i^k$ ($k=1,...,5$) by using

$$V\left({\chi }_i^k\right)=-k_{B}T\log P\left({\chi }_i^k\right)$$

(3)

where $P\left({\chi }_i^k\right)$ is the PDF of the dihedral angle ${\chi }_i^k$. One-dimensional PDFs were computed on a circle from time series extracted from MD (Figure 2B) using concatenated MD trajectories, as performed for 2D PDFs (see above). An example of free-energy profile for side-chain dihedral angle ${\chi }_{84}^1$ (Lys) is shown in Figure 2D. Moreover, modifications between two FEPs of identical side-chain dihedral angle ${\chi }_i^k$ due to different ligand-binding forms of GST were quantified using the similarity index H between their associated 1D PDFs [28]:

$$H_{1|2}\left({\chi }_i^k\right)={\displaystyle \frac{{\displaystyle 2{\int }_{-\pi }^{+\pi }{P}_1\left({\chi }_i^k\right)P_2\left({\chi }_i^k\right)d{\chi }_i^k}}{{\displaystyle {\int }_{-\pi }^{+\pi }{P}_1{\left({\chi }_i^k\right)}^{2}d{\chi }_i^k+{\int }_{-\pi }^{+\pi }{P}_2{\left({\chi }_i^k\right)}^{2}d{\chi }_i^k}}}$$

(4)

where $P_1\left({\chi }_i^k\right)$ is the PDF of SCA $({\theta }_i,{\gamma }_i)$ in binding form 1, and $P_2\left({\chi }_i^k\right)$ is the PDF of SCA $\left({\chi }_i^k\right)$ in binding form 2. The same scale of similarity, as described above for FESs, is used to quantify similarities (large, moderate, and low).

3. Results

3.1. Structural Flexibility of Human GSTA1 in Its APO, GSH, and GS-R Form

Figure 3A presents thermal B-factors computed during MD simulations of hGSTA1 in its APO, GSH, and GS-R form. First, the correlation between B-factors of monomers A and B in each of the three forms independently is very high (more than 80%). Therefore, we compare average B-factors computed over the two monomers for each form. Overall, four regions of the hGSTA1 enzyme were found to exhibit a large flexibility (Figure 3B): the ${\alpha }_2$ region (residues 38–50), which contains residue Arg45 that belongs to the G site; the loop between ${\alpha }_4$–${\alpha }_5$ helices (residues 107–119), i.e., L${}_{4,5}$, which contains residues Leu107, Leu108, and Val111 that belong to the H site, the loop between ${\alpha }_8$–${\alpha }_9$ helices, i.e., L${}_{8,9}$ (residues 199–206), and, finally, the C-terminal region (residues 210–222) containing the ${\alpha }_9$ helix with residues Phe220 and Phe222 that belong to the G and H sites, respectively. By computing the difference of flexibility of hGSTA1 between the three forms shown in Figure 1C, we found that the APO and GSH forms are more flexible than the GS-R form in the ${\alpha }_2$ region, which contains residue Arg45. This flexibility in the absence of the substrate (APO and GSH forms) may promote its binding, the same region being constrained by the presence of the conjugated substrate (GS-R form). On the opposite, the loop and extremities of ${\alpha }_4$ and ${\alpha }_5$ helices are more flexible in the GS-R and GSH forms of hGSTA1 than in its APO form. It means that the presence of the ligand in the G and H sites of GST and, particularly, the interactions with residues Ala100, Asp101, Ile106, Leu107, Leu108, Val111, (in ${\alpha }_4$) and Arg131 (in ${\alpha }_5$) modify the flexibility of this region (Figure 1B). Moreover, the presence of GSH in the G site strongly increases the flexibility of the loop $L_{8,9}$ located a few residues away for the C-terminal helix ${\alpha }_9$ of hGSTA1. This loop is also slightly flexible in the APO form and is not at all in the GS-R form. Last, but not least, the APO form of hGSTA1 shows a much larger flexibility in the C-terminal region (residue 210–222) than both the GSH and the GS-R forms. The sum of the B-factors in this region for each of the forms is 39, 21, and 12 nm${}^2$ for APO, GSH, and GS-R, respectively. This result is in agreement with experimental observations showing that the protein flexibility of hGSTA1, including the dynamics of C-terminal ${\alpha }_9$ helix on nanosecond-millisecond timescale and the protruding extremities of ${\alpha }_4$–${\alpha }_5$ helices, contributes remarkably to the catalytic and noncatalytic ligand-binding functions of GSTs [32,33].

3.2. Identification of Key Residues Involved in Ligand Binding to hGSTA1 from Free-Energy Landscape Analysis of CGAs and SCAs

The network of residues influenced by ligand binding into hGSTA1 in the three different forms presented in Figure 1C is determined by computing similarity indices H between (i) 2D free-energy surfaces of CGAs $(\theta ,\gamma )$ and (ii) 1D free-energy profiles of SCAs $\left(\chi \right)$ along the amino-acid sequence. First, only FESs that show large or moderate similarities ($H>0.3$) between monomers A and B of the homodimer are considered as the DATA. This ensures that non-converged FESs and FEPs between the two monomers are not taken into account. In fact, even though monomers A and B can have an asymmetrical dynamical behavior from MD due to the fact they are flexible monomers (sampling issue), no information can be exploited from these non-converged characteristics from a thermodynamical point of view. They might or might not be relevant to the conformational fluctuations of the protein. In total, 11 CGAs were excluded ($i=50$ for APO, $i=141,145,220$ for GSH, and $i=144,146,147,149,217,218,219,220$ for GS-R form). These CGAs are mainly located in the loop between helices ${\alpha }_5$ and ${\alpha }_6$ and in the C-terminal region, and they represent only 5% of the total number of CGAs ($N=218$).

Second, we computed similarity indices H between FESs of CGAs (from DATA) between the different forms of hGSTA1, (i) APO vs. GSH, (ii) GSH vs. GS-R, and (iii) GS-R vs. APO. Figure 4 shows dissimilarity $1-H$ along the amino-acid sequence for each comparison for both CGAs $(\theta ,\gamma )$ (panel A) and SCAs ${\chi }^1$ (panel B). On one hand, from APO vs. GSH comparison, only one CGA ${(\theta ,\gamma )}_i$ presents moderate similarity of its FES between the two forms. This result is not very surprising since at the microsecond timescale, the association/dissociation of GSH ligands from G sites of monomer A and B of hGSTA1 were observed during MD simulations (see Figure S3). It concerns CGA ${(\theta ,\gamma )}_{143}$ and involves residues Ser142, His143, Gly144, and Gln145, which are not directly involved in the ligand-binding sites of hGSTA1. However, these residues are located at the end of ${\alpha }_5$-helix, which contains residue Arg131 of the G site. On the other hand, comparisons involving the GS-R form of hGSTA1 present many more structural modifications. From GSH vs. GS-R comparison, seven CGAs ${(\theta ,\gamma )}_i$ present moderate or low similarities of their FESs between these two forms, for $i=13,14,15,16,17,18,51$. These seven CGAs involve 13 different amino acids with Gly14, Arg15, and Ser18, which are directly involved in the ligand-binding site of hGSTA1 (Figure 4C). Among the ten other residues detected and which are not located in the binding sites of hGSTA1, Arg20 is of great interest since this amino acid is highly conserved through GST classes. Moreover, by looking at FES of CGA ${(\theta ,\gamma )}_{15}$ in the three ligand-binding forms of hGSTA1 (Figure 4D), we show that both $\theta$ and $\gamma$ angles are modified by the presence of gluathione-S-conjugated ligand compared to the APO and GSH forms. For instance, the global minimum of the FES $V\left({\theta }_{15}\right)$ is significantly shifted to lower values of $\theta$ and there is the creation of a new global minimum for $V\left({\gamma }_{15}\right)$ in the GS-R form compared to the GSH form. It is also the case for residues 13 and 17 (see Figure S5). Therefore, the presence of GS-R ligand modifies the structure of hGSTA1 for some residues, which corresponds to an induced fit model of protein–ligand binding. On the opposite, for the comparison between the APO and GSH forms, the population shift model of protein–ligand binding is more pronounced since the similarity between FESs and FEPs is larger than comparisons with GS-R forms (Figure 4A,B).

Finally, 10 FESs of CGAs were identified as moderately or largely modified due to ligand binding for the comparison between GS-R and APO forms of hGSTA1. In detail, CGAs ${(\theta ,\gamma )}_i$ for $i=10,12,13,14,15,16,17,18,53$, and 158 were identified. Among these 10 CGAs, six were already detected from GSH vs. GS-R comparison (Table S1). It means that these six CGAs are GS-R specific, i.e., they are similar in the APO and GSH forms compared to the GS-R form. Particularly, CGA ${(\theta ,\gamma )}_{53}$, which involves residues Gln54 and Val55, is of special interest since both residues belong to the G site. In addition, CGA ${(\theta ,\gamma )}_{158}$, which involves residue His159, is also of special interest since this residue belongs to the H site. In total, from the three different comparisons of ligand-binding forms of hGSTA1 shown in Figure 4A, 12 unique CGAs were identified to be influenced by the absence/presence of a ligand; some of them being detected in two different comparisons (Table S1). It corresponds to a total of 26 different residues with four residues which belong to the G site (Tyr9, Arg15, Gln54, and Val55) and four residues which belong to the H site (Phe10, Gly14, Ser18, and His159). All 2D FESs of the 12 CGAs identified to be involved in the ligand binding process in hGSTA1 are presented in Figure S5.

We applied the exact same procedure to compare the one-dimensional FEPs of SCAs ${\chi }^k$. First, three amino acids were excluded from the DATA since they do not exhibit converged FEPs of SCAs between monomers A and B of hGSTA1. It concerns residues Glu146 (${\chi }^{1,2}$), His159(${\chi }^{1,2}$), and Arg221 (${\chi }^{1,2,3,4,5}$). It represents only 1.5% of the total number of residues characterized by SCAs ($N=192$). Then, from dissimilarity indices, $1-H$, computed along the amino-acid sequence for the three different comparisons APO vs. GSH, GSH vs. GS-R, and GS-R vs. APO, as shown in Figure 4B for SCAs ${\chi }^1$ and in Figure S3 for the other SCAs ${\chi }^{2,3,4,5}$, a total of 19 SCAs were identified to be influenced by the absence/presence of a ligand in the binding sites of hGSTA1 homodimer. It concerns 15 different residues, i.e., Phe10, Arg13, Arg15, Met16, Glu17, Thr19, Arg45, Gln54, Leu102, Phe136, His143, Arg155, Leu163, Phe220, and Phe222 (Table S1). Among these 15 residues, four are located in the G site (Arg15, Arg45, Gln54, and Phe220) and two in the H site (Phe10 and Phe222). In detail, the comparison between the APO and GSH forms of hGSTA1 highlights three SCAs, ${\chi }_{143}^{1,2}$ and ${\chi }_{222}^1$ (C-term). Particularly, residue His143 was also found to be crucial (the only one) for ligand binding from the analysis of the CGA (${\theta }_i,{\gamma }_i$) FESs. It means that residue His143 is an unexpected probe of ligand binding as both its local main chain conformation and its side chain are modified by the absence (APO) or presence (GSH) of the ligand, even though His143 is not directly located in the G or H binding site. In fact, this residue is at the surface of the hGSTA1 enzyme, at a distance around 30 Å from the binding sites.

Furthermore, 11 SCAs are found to show free-energy dissimilarities between their GSH and GS-R forms, i.e., ${\chi }_{10}^{1,2},{\chi }_{15}^3,{\chi }_{16}^3,{\chi }_{19}^1,{\chi }_{45}^3,{\chi }_{54}^2,{\chi }_{136}^2,{\chi }_{155}^{1,3}$, and ${\chi }_{220}^1$. The six residues mentioned above for the comparison of APO vs. GSH forms and which belong to the G and H sites were also identified above from the comparison of CGA FESs for APO vs. GSH. In particular, SCA ${\chi }_{15}^3$ is associated with Arg15 which is conserved in GST class Alpha and located at the interface between G and H sites. For this SCA, as shown in Figure 4E, the presence of GS-R ligand modifies the position of the global minimum of the FEP, from −180${}^{\circ }$ in both APO and GSH forms to −60${}^{\circ }$. This local conformation exists in APO and GSH forms but as a local minimum (population shift). Finally, 12 SCAs were detected for the comparison GS-R vs. APO. Among these 12 SCAs, seven were already identified from the comparison GSH vs. GS-R, i.e., residues Phe10, Arg15, Met16, Thr19, Phe136, and Phe220 are, therefore, GS-R specific. The same kind of behavior was already observed for CGAs of the main chain with six coordinates GS-R specific. Among these six coordinates, residue Arg15 was identified, which support its case as one of the most important residues involved in ligand binding of hGSTA1. The five supplementary SCAs identified are ${\chi }_{13}^3,{\chi }_{17}^3,{\chi }_{102}^1$, and ${\chi }_{163}^{1,2}$, which correspond to residues Arg13, Glu17, Leu102, and Leu163. Particularly, Leu102 is located very close to the G site, Asp101 belonging to the G site and located in the opposite monomer of the dimer. All 1D FEPs of the 19 SCAs identified to be involved in the ligand binding process in hGSTA1 are presented in Figure S6.

3.3. Collective Motions of the Network of Residues Involved in Ligand Binding in hGSTA1 Revealed by PCA

As described in detail above, the free-energy landscape (FEL) comparison of CGAs $(\theta ,\gamma )$ reveals that the absence/presence of ligands in hGSTA1 modifies at least 12 internal coordinates, located in both binding sites G and H or even further in the sequence (Figure 4C). In other words, the modification of the (unknown) free-energy landscape of hGSTA1 upon ligand binding can be represented here quantitatively by its projection over 12 dimensions. Obviously, these 12 coordinates are not independent and the local motions of the CGAs are coupled to each other. Hereafter, we establish the relation between these 12 local conformational changes by applying principal component analysis (see Supplementary Materials for details). From PCA for CGAs $(\theta ,\gamma )$ in the three different forms of hGSTA1, we found that more than 30% of the total contribution to the MSF of the 12 CGAs is captured by the first two collective modes over the 72 existing ones (Figure S7). Then, we computed time series of the projections of vectors $q_i$ along the eigenvectors of collective modes 1 and 2, which defined the collective coordinates PC${}_1$ and PC${}_2$. From the two-dimensional probability density functions of PC${}_1$ and PC${}_2$ time series computed from MD trajectories, the effective 2D FEL $V_{1,2}$, shown in Figure 5, was computed for the three forms APO, GSH, and GS-R. The largest motions defined by the first and second PCs can be interpreted in terms of specific features of the two-dimensional FEL $V_{1,2}$. The 2D FEL of the APO form of hGSTA1 presents two minima defined from the contour lines of free-energy (<2 k${}_B$T), whereas both GSH and GS-R forms present three minima. In addition, for GSH and GS-R forms, the different minima are well separated by energy barriers from around 3 to 4 k${}_B$T, whereas in the APO form, the two minima are separated by a barrier around 1 k${}_B$T. These 2D FELs are, in fact, in agreement with the observation that FEPs of $\theta$ and $\gamma$ angles (Figure 4D and Figure S5) are more spread out and therefore contain more wells for GSH and GS-R forms in addition to the one(s) of the APO form.

From the 2D FEL shown in Figure 5A, the coupling between the 12 CGAs of the network identified as important for ligand binding in hGSTA1 was estimated by computing the influences of each coordinate into the largest amplitude collective motions, which are of particular interest for the biological function of proteins [34]. First of all, as shown in Figure 5B for each form, the dynamics of the 12 coordinates are highly asymmetrical between monomers A and B of hGSTA1. For collective mode 1 of APO, GSH, and GS-R forms, the contribution of monomer A is 66, 24, and 13%, respectively (with 34, 76, and 87% for monomer, B). For collective mode 2, the contribution of monomer A is even larger. Moreover, for the collective mode 1 of GSH form, the Pearson correlation of influences along the 12 CGAs is large, i.e., ${\rho }_{A|B}=0.84$, whereas for the collective mode 2, the correlation is low, i.e., ${\rho }_{A|B}=0.23$. However, PCA reveals two major dynamical couplings between residues. First, there is a coupling between CGAs 13–14–15 (including residues Gly14 of the H site and Arg15 of the G site) and CGAs 51, 53 (including residues Gln54 and Val55 of the G site) in the same monomer A or B. This coupling can be more or less pronounced between monomers and is modified by the presence of GS-R ligand, which modifies the location of CGA fluctuations in collective mode 1 (Figure 5B). Second, there is a coupling between the two binding sites of the two monomers, as shown particularly for CGAs 51, 53 in collective mode 1 for both APO and GSH form. This coupling is also present but is less pronounced in collective mode 2 of GS-R form. Therefore, dynamical coupling of intra-binding sites and inter-binding sites of hGSTA1 along the main chain of the enzyme was found and appears to be crucial for ligand binding in hGSTA1.

We applied the same procedure to the 19 SCAs detected from FEPs similarities. First, the two collective modes 1 and 2 contribute to 17, 18, and 36% for APO, GSH, and GS-R forms of hGSTA1, respectively (Figure S6). As shown in Figure 5C, 2D FELs $V_{1,2}$ from SCAs present much more minima than for CGAs. It comes from the fact that SCAs show unrestricted FEPs with values between −180 and +180${}^{\circ }$ compared to CGAs. Moreover, the shape of the FEL for the GS-R form of hGSTA1 is much more different compared to the APO and GSH forms. The GS-R form of hGSTA1 explores the free-energy landscape less and is more restricted. However, the free-energy barriers between minima are much larger than the APO and GSH forms (up to 6 k${}_B$T), for which the multiple minima are more easily accessible to the system. This behavior is related to the fact that the binding of GS-R ligand, which is bigger than GSH ligand, involves more constraints on the side chains of the residues and, therefore, larger free-energy barriers to overcome. In addition, compared to CGAs, collective modes 1 and 2 are more spread out over the two monomers, particularly in the GSH form (Figure 5D). The analysis of the influences of each SCA in the collective motions of mode 1 and 2 reveals a strong coupling Leu163 and residue Glu17 and its neighbors in the G site in the APO form. In the GS-R form, the coupling is more pronounced with residue Arg155.

4. Discussion

From the three different comparisons of ligand-binding forms of human GSTA1 shown in Figure 1C, 12 CGAs and 19 SCAs were identified to be influenced by the absence/presence of ligands, most of them being gluthatione-S-conjugated-specific (Table S1). It corresponds to a total of 33 different amino acids (15% of the total sequence), eight of them being identified from both CGAs and SCAs analyses (

Table 1. Summary of amino acids involved in the ligand binding to hGSTA1 and extracted from similarity of FESs of CGAs and FEP of SCAs. Amino acids, which were detected from both analyses, are underlined. Amino acids, which are highly conserved among GSTs, are denoted in bold. Secondary structure location of amino acids are indicated in parentheses.

Location	List of Amino Acids
G site	Tyr9 (${\beta }_1$), Arg15 (L${}_{{\beta }_1,{\alpha }_1}$), Arg45 (${\alpha }_2$), Gln54, Val55 (L${}_{{\alpha }_2,{\beta }_3}$), Phe220 (${\alpha }_9$)
H site	Phe10, Gly14 (L${}_{{\beta }_1,{\alpha }_1}$), Ser18 (${\alpha }_1$), His159 (${\alpha }_6$), Phe222 (${\alpha }_9$)
Others	Asn11, Ala12, Arg13 (L${}_{{\beta }_1,{\alpha }_1}$), Met16, Glu17, Thr19, Arg20 (${\alpha }_1$), Leu50, Met51, Phe52, Gln53 (L${}_{{\alpha }_2,{\beta }_3}$), Leu102 (${\alpha }_4$), Phe136, Ser142, His143, Gly144, Gln145 (${\alpha }_5$), Arg155, Asp157, Ile158, Leu160, Leu163 (${\alpha }_6$)

). Among these 33 residues, 11 of them (32%) are directly located to the two binding sites of hGSTA1, with 6 which belong to the G site, i.e., Tyr9, Arg15, Arg45, Gln54, Val55, Phe220, and 5 which belong to the H site, i.e., Phe10, Gly14, Ser18, His159, and Phe222 (Figure S8). Furthermore, a collective motion involving residues Arg15, Gln54, and Val55 of the G site and residue Gly14 of the H site was identified using principal component analysis, revealing a strong coupling between these residues. Particularly, residue Arg15, which is found here to be involved in both local conformational changes of the main chain and of the side chains of hGSTA1, was demonstrated to have crucial catalytic activity both theoretically [35] and experimentally [36,37]. Finally, residues Phe220 (G site) and Phe222 (H site) belong to the C-terminal part of hGSTA1 (${\alpha }_9$-helix; see Figure S8) and it has been demonstrated that this secondary structure contributes remarkably to the catalytic and noncatalytic ligand-binding functions of GSTs [32,33].

At large distance from the ligand-binding sites of hGSTA1, 12 residues (from Leu102 to Leu163; see Table 1) were identified as sensitive to the ligand disturbance. For example, His143 is involved in the main-chain and side-chain conformational changes of hGSTA1 and is located around 30 Å away from the G and H sites Figure S8). Moreover, GST sequence analysis indicates that some residues form an highly conserved local sequence GXXh(T/S)XXDh (h: hydrophobic), constituted by the ${\alpha }_6$-helix and its preceding loop [38,39,40]. This motif belongs to a substructure named N-capping box and hydrophobic staple motif of GST and has been shown to be critical for the protein folding and stability [41]. In the present work, six residues were identified from the ${\alpha }_6$-helix, particularly Arg155, Asp157 (highly conserved), and Ile158 (hydrophobic) which belongs to the corresponding motif GNKLSRADI. According to experimental structures, Ile158 can form hydrophobic interactions with the ${\alpha }_1$-helix [41], which is an important structural element supporting the active binding site of GSTs.

Finally, it has been reported that the sequence alignment of all known GST structures (more than 100) shows that only 6–7 residues, i.e., less than 5% of the entire polypeptide chain, are strictly conserved [42]. For instance, Pro56 of the GST Alpha class, which belongs to the G site, is conserved through GST classes Mu, Pi, Sigma, Phi, Tau, Theta, Zeta, Omega, and Beta. In the present work, four amino acids, i.e., Arg13, Arg20, Val55, and Asp157, were detected and are highly conserved among GST classes. Particularly, Arg20 is strictly conserved through GST classes Alpha, Mu, Pi, Sigma, Tau, Zeta, and Omega; Arg13 is strictly conserved through GST classes Alpha, Mu, and Pi [43].

5. Conclusions

In the present work, we performed all-atom classical MD simulations in explicit solvent of human GSTA1 enzyme in the three different following forms of the protein: (i) the APO form, when no ligand is bound to hGSTA1, (ii) the GSH form, when the gluthatione ligand is bound to the G site of hGSTA1, and (iii) the GS-R form, when the gluthatione-S-conjugate is bound to the G and H sites of hGSTA1. From MD runs, we performed a free-energy landscape analysis of internal coordinates along the amino-acid sequence of the protein. The analysis was conducted using coarse-grained angles $(\theta ,\gamma )$ of the main-chain (CGAs) and side-chain angles $\chi$ (SCAs), as local probe to track conformational changes. The comparison between each form of hGSTA1 reveals 12 CGAs and 19 SCAs influenced by the ligand binding into the G and H sites. It corresponds to a total of 33 residues of hGSTA1 which are involved in the protein–ligand binding process. Among the 33 residues identified, 11 of them belong to the two binding sites of hGSTA-1 identified experimentally from XRD structures, which confirms that the method developed and applied here is very robust. Finally, the dynamics of coarse-grained and side-chain angles were studied using principal component analysis. It shows an asymmetrical behavior between the two monomers of hGSTA1, particularly from CGAs. Moreover, the analysis of the first collective motions reveals a strong dynamical coupling between residues Arg15-Gln54-Val55 of the G site and Gly14 of the H site. These residues are coupled in the same binding site (intrasite) but also between the two binding sites of each monomer of hGSTA1 (intersites), which is an important result of the present work. Finally, residue Glu17 shows important coupling with residues Arg155 and Leu163, depending on the binding form of hGSTA1 (APO and GS-R, respectively). This free-energy landscape methodology, already applied successfully to chaperone proteins [28] and now to protein–ligand binding of human GSTA1 enzyme, is very powerful to identify network of residues and their dynamical coupling involved in the communication between domains or binding sites of biological systems. An important outcome of the present approach is its ability to decipher long-range effects, i.e., influence of residues far from the binding site on the ligand–protein association. This might guide engineered mutagenesis to tailor specific ligand–protein interactions.