## 3 Results
### 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 (Figures 4D and 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.