Network of Conformational Transitions Revealed by Molecular Dynamics Simulations of the Carbonic Anhydrase II Apo-Enzyme.

Human carbonic anhydrase II (HCA II) is an enzyme that catalyzes the reversible hydration of CO2 into bicarbonate (HCO3-) and a proton (H+) as well as other reactions at an extremely high rate. This enzyme plays fundamental roles in human physiology/pathology, such as controlling the pH level in cells and so on. However, the binding mechanism between apo-HCA II and CO2 or other ligands as well as related conformational changes remains poorly understood, and atomic investigation into it could promote our understanding of related internal physiological/pathological mechanisms. In this study, long-time atomic molecular dynamics simulations as well as the clustering and free-energy analysis were performed to reveal the dynamics of apo-HCA II as well as the mechanism upon ligand binding. Our simulations indicate that the crystallographic B-factors considerably underestimate the loop dynamics: multiple conformations can be adopted by loops 1 and 2, especially for loop 1 because loop 1 is one side of the binding pocket, and its left-to-right movement can compress or extend the binding pocket, leading to one inactive (closed) state, three intermediate (semiopen) states, and one active (open) state; CO2 cannot get into the binding pocket of the inactive state but can get into those of intermediate and active states. The coexistence of multiple conformational states proposes a possible conformational selection model for the binding mechanism between apo-HCA II and CO2 or other ligands, revising our previous view of its functional mechanism of conformational change upon ligand binding and offering valuable structural insights into the workings of HCA II.

Introduction

Enzymes, or proteins in general, are not static and can adopt multiple conformations, as demonstrated in three-dimensional atomic detail by the observation of structures of the same enzyme in different liganded states through X-ray crystallography and nuclear magnetic resonance (NMR) or even by the multiple conformational equilibria of an enzyme in the same liganded state revealed by NMR.[1,2] It is well-known that these multiple conformational transitions are essential for ligand binding because of the conflicting structural requirements: an unbound enzyme must adopt an open conformation to allow its substrate to enter its active site to form a Michaelis complex; subsequently, to maximize transition state stabilization and prevent side reactions, the enzyme typically assumes a closed conformation, such that the substrate is buried in its active site; after the chemical reaction, the enzyme must again open up its active site to allow the product to exit.[3] However, how a protein could transit from an unbound conformation to a bound conformation in a complex with a ligand, or specifically, whether an unbound enzyme can undergo large conformational changes and assume bound as well as other multiple conformations remains unknown.[4−7] Insights into it will promote our understanding of protein–ligand binding, which is critical to molecular recognition as well as the drug design.[8]

Human carbonic anhydrase II (HCA II) has been long known as an enzyme that catalyzes the reversible hydration of CO2 into bicarbonate (HCO3–) and a proton (H+) at an extremely high rate,[9−13] but its catalytic activity is not limited to the hydration of CO2: it is reported that HCA II can also catalyze the hydrolysis of 1-fluoro-2,4-dinitrobenzene[14] and sulfonyl chlorides[15] as well as the hydration of cyanamide to urea.[16] This enzyme plays fundamental roles in human physiology/pathology: it is essential in keeping the adequate balance between carbon dioxide and bicarbonate, thus controlling the pH level in cells, and is vital to a variety of biological processes, such as vision,[17] development, and function of bone,[18] calcification,[18] and so forth. As a cytosolic and monomeric protein with a small molecular weight of 30 kDa but multifunctions,[9] it is often used as a prototypical model system for biophysical studies[9,13] and medicinal chemistry applications.[13,19] The 15 Å deep conical binding pocket[20] (illustrated in Figure ) is formed with 10-stranded β-sheets on one side,[21] and on the other side, with the active site loop of residues 197–206 (loop 1), which has been identified to be important for ligand binding,[22] especially Leu198, Thr199, and His200, which have direct contact and form hydrogen bonds with ligands;[22−25] at the bottom of the pocket, a Zn2+ ion is tetrahedrally coordinated with three histidine residues (His94, 96, and 119) and a bound water/hydroxyl (H2O/OH–).[26] Another loop related to ligand binding is the surface loop of residues 230–240 (loop 2), which can control the size of the entrance of the binding pocket.

Figure 1

The crystal structure of apo-HCA II (PDB entry 3KS3). The zinc ions are colored orange, the loops 1 and 2 are highlighted in red, whereas residues 68, 202 and 235 are marked in green. His94, Val143, Leu198, Thr199, and Val207 are displayed as the licorice style.

The flexibility of the two loops in the apo form of the enzyme is evidenced by their higher B-factors in the crystal structure; however, in the crystal structure, the flexibility is always greatly hindered by crystal contacts between adjacent copies and lower cryogenic temperature;[27−29] more importantly, the detail of their conformational changes cannot be revealed either. The conformational changes of the two loops could impact the ligand binding of the enzyme, especially for loop 1, although its B-factor values are not as large as those of loop 2, even a tiny conformational change of this loop could change the volume of the binding pocket as well as the hydrogen-bond network, which is critical to ligand binding and catalytic reactions. Therefore, in our work, conventional molecular dynamics (MD) simulations totaling 1 μs at room temperature (300 K) were applied to research the conformational transitions of apo-HCA II. Our simulations revealed multiple conformations of loops 1 and 2, including the open, semiopen, and closed conformations, thereby revising our previous view of the functional mechanism of conformational change upon CO2 or other ligand binding.

Results

Stability of the MD Trajectories

Because a stable protein structure was essential to our analysis, the root-mean-square deviation (rmsd) values were evaluated for all MD trajectories. The average rmsd value from the initial structure calculated using all Cα atoms is ∼2 Å (Figure a), indicating that the protein was stable during each of the 200 ns MD simulations. As measured by the B-factor values calculated for the MD trajectories (Figure b), loop 2 is one of the most flexible regions with its MD B-factors larger than 100 Å2; loop 1 also has a peak of MD B-factor values, revealing some large conformational change on it which, although higher than experimental values as expected (crystal contacts between adjacent copies and lower cryogenic temperature as mentioned in the Introduction section), is still well-consistent. In the simulations, the residues with the highest B-factor values in loops 1 and 2 were Pro202 and Gly235, respectively; correspondingly, Val68 is located at the edge of the protein cavity and has little mobility. Therefore, the distance between Pro202 and Val68 can be used to represent the motion of loop 1 and its opening or closing; likewise, the distance between Gly235 and Val68 can reflect the motion of loop 2 (see Figure ).

Figure 2

(a) Time evolution of Cα rmsds during the MD simulations and (b) the B-factor distributions derived from the MD simulations and the crystal structure.

Clustering and Conformational Network Analysis for the Trajectories

To map out the conformational states and transitions of loops 1 and 2, clustering and conformational network analyses were performed on the conformations of the two loops in the simulations, which is depicted in Figure . Five clusters are extracted from the MD simulations, and in them, multiple conformations of loops 1 and 2 are observed (Figure a); especially for loop 1, even though its B-factors are not as large as loop 2, loop 1 is one side of the binding pocket and so its movement can directly compress or extend the binding pocket, representing closed (clusters 2), semiopen (clusters 1, 3, and 4), as well as open conformations (cluster 5) of the binding pocket. For loop 2, although its B-factors are larger than loop 1, its conformations are open (clusters 1, 2, 3, and 5) or even superopen (cluster 4), keeping the entrance of the binding pocket always large enough for ligands.

Figure 3

Distance root mean square (DRMS)-based clustering and network analysis of the MD simulations: (a) center structures of the clusters, (b) time evolution of the clusters, and (c) distribution and interconversion network of the clusters. In panel b, the different clusters are represented by different colors. In panel c, the relative populations of the clusters are indicated by percentage values, and the thickness of each connecting line is proportional to the corresponding frequency of interconversion.

Figure b shows the time series of the cluster assignments for the five trajectories. Although all trajectories start with the same initial structure belonging to cluster 1, the observed cluster types are different in different trajectories, which means that different trajectories have different conformational spaces. Clusters 1 and 2 are all included in the five trajectories, representing the broad observation of the transitions between closed and semiopen conformations of loop 1, as well as that between open and superopen conformations of loop 2; except trajectory 2, cluster 3 is present in the other four trajectories; the most dramatic conformational transitions of the two loops occur in trajectories 3 and 5, in which all the five clusters are present, which leads to multiple transitions between closed, semiopen, and open conformations of loop 1.

The occurrence percentages of these clusters and the conformational space network among them are shown in Figure c. Clusters 1 and 2 are dominant, with the percentages up to 45.2 and 53.4%, respectively; the other three clusters totally occupy only 1.4%. From the conformational space network, one can observe that cluster 4 is a key intermediate state: first, it provides another way from cluster 1 to 2 or vice versa; more importantly, cluster 5 connects only to cluster 4; therefore, only through cluster 4, cluster 5 can be reached by clusters 1 and 2. Because loop 1 is closed in cluster 2 while open in cluster 5, with loop 1 semiopen, cluster 4 is the intermediate state between closed and open conformations of loop 1. Another interesting cluster is cluster 3: it only has direct transitions from or to cluster 2.

To precisely determine which conformations in our simulations are accessible to CO2 and other ligands, the experimental structure of CO2-bound HCA II (PDB entry 5DSO) and the multiple conformations of apo-HCA II from our MD simulations are aligned. As shown in Figure , in the closed conformations of loop 1 (Figure a), the entrance to the binding pocket is partially blocked, and the binding pocket is also too small to accommodate even the smallest ligand–CO2, and detailed investigation revealed that the binding pocket is occupied by His94, Val143, Leu198, Thr199, and Val207 (See Figure ); in the semiopen conformations of loop 1, the entrance and the binding pocket are enough for CO2 entering and binding; in the open conformations, the entrance and the binding pocket are even larger, provide a possibility to accommodate and catalyze larger ligands, such as 1-fluoro-2,4-dinitrobenzene, sulfonyl chlorides, and cyanamide. Therefore, the five clusters can be identified as one inactive state (cluster 2), three intermediate states (clusters 1, 3, and 4), and one active state (cluster 5) with their percentages 53.4, 46.1, and 0.5% (see Figure c), which means that about one-half of the conformations from our simulations are accessible to CO2, whereas the other half are not.

Figure 4

The alignment of the experimental structure of CO2-bound HCA II (PDB entry 5DSO) and the multiple conformations of apo-HCA II from our MD simulations: (a) closed loop 1; (b) semiopen loop 1; (c) open loop 1. With the experimental structure of the protein hidden, the graphs show the relative positions of bounded CO2 and the multiple conformations from our MD simulations. The CO2 molecule is displayed as three connected spheres with the C atom in green and the two O atoms in red, and the conformations from our simulations are in “surface” view.

The Free-Energy Landscape of Loop 1

As loop 1 is critical to ligand binding to HCA II, its two-dimensional free energy was estimated from our MD simulations and is shown in Figure . The Cα distance between Pro202 and Val68 was chosen as the x-coordinate, whereas the rmsd of loop 1 to the initial structure (PDB entry 3KS3) was chosen as the y-coordinate. The reason for such a choice is that the Cα distance between Pro202 and Val68 can reasonably discriminate different conformations (closed, semiopen, open) of loop 1, whereas the rmsd of loop 1 represents its conformational changes. Points representative of the central conformations of the five clusters are also shown in the graphs.

Figure 5

Two-dimensional free energy for loop 1 from the MD simulations. The x axis represents the Cα distance between Pro202 and Val68, whereas the y axis represents the rmsd of loop 2 to the initial structure (PDB entry 3KS3). Points representative of the central conformations of the five clusters are also shown in the graphs.

In Figure , clusters 1 and 2 are located at the area with lower free energies, which are consistent with those of Figure c; they have the largest populations. Cluster 4 is just between cluster 5 and cluster 1 or 2, thus agreeing well also with that of Figure c. Cluster 4 is a key intermediate state between closed and open conformations of loop 1 because the representative point of cluster 1 is closer to that of cluster 4 than cluster 2; just as shown in Figure c the transitions between cluster 1 and cluster 4 are more frequent than those between cluster 2 and cluster 4. Estimated from this figure, the free-energy difference between closed and open conformations of loop 1 is about 6.75–3.75 = 3 kcal/mol.

Principal Component Analysis

In our principal component analysis (PCA), all Cα atoms are used to define the backbone conformation of HCA II, resulting in 771 principal components from the 257 Cα atoms in the form of 771 eigenvectors and their associated eigenvalues. The validation of the PCA is certified in the Supporting Information (Table S1). Figure S2 shows that, consistent with earlier research,[30] 30 out of the 771 PCs (principal components) capture approximately 80% of the protein’s motion. The first three PCs, in particular, have much larger eigenvalues than the rest, and in every trajectory, they capture about 50% of the motion (see Figure S2). Hence, the first three PCs are sufficient for the analysis.

Figure shows the vector-field representations of the PCs obtained from the combined data of all trajectories. The figure clearly shows that loops 1 and 2 have one of the largest motions of the protein, which is well-consistent with the B-factor graphs. From the figure, one can observe that, in the first three PCs, both loops 1 and 2 always move from left to right or vice versa; however, the consequences are quite different because loop 1 is one side the binding pocket, and such a movement can compress or extend the binding pocket; but for loop 2, it is a loop on the surface, and such a movement has little impact on the binding pocket. This is why loop 1 has closed, semiopen, and open conformations, but loop 2 has only open or superopen conformations.

Figure 6

Vector-field representations and corresponding eigenvalues of the first three PCs obtained from the combined data of the MD trajectories.

Discussion

The Dynamics of Loops 1 and 2 of Apo-HCA II Revealed by the MD Simulations

HCA II catalyzes the reversible hydration of CO2 into HCO3– and a H+ at an extremely high rate, thus playing very fundamental roles in human physiology/pathology, such as controlling the pH level in cells. Besides that, HCA II can also catalyze the hydrolysis of 1-fluoro-2,4-dinitrobenzene and sulfonyl chlorides as well as the hydration of cyanamide to urea. Atomic investigation into this enzyme will promote our understanding of related internal physiological/pathological mechanisms.

In this study, long time (1 μs) atomic MD simulations of apo-HCA II as well as the clustering and free-energy analyses were reported, which has provided novel insights into the dynamics of apo-HCA II as well as the mechanism upon ligand binding. Multiple conformations of loops 1 and 2 can be observed in our MD simulations especially for loop 1 because loop 1 is one side of the binding pocket, and its left-to-right movement can compress or extend the binding pocket, leading to one inactive state (cluster 2 with loop 1 closed), three intermediate states semiopen (clusters 1, 3, and 4, with loop 1 semiopen), and one active state (cluster 5 with loop 1 open); both the intermediate and active states are accessible to CO2, the active state can accommodate large ligands. Loop 2 is a surface loop, so although it has larger B-factors than loop 1, its left-to-right movement can represent only open or even superopen conformations, never blocking the entrance to the binding pocket.

“Induced Fit” versus “Conformational Selection”

At present, two different models, which are generally referred to as “induced fit (IF)” and “conformational selection (CS)”, have been proposed to describe the conformational changes that occur upon ligand binding.[31−33] The IF model[34] presumes that the binding of the ligand may induce a conformational change in the receipt or enzyme, resulting in an optimized geometry that exists only in the complex state. In this scenario, the protein energy landscape changes upon the binding of the ligand. By contrast, the CS model presumes that more than one structural conformation pre-exists in conformational equilibrium before ligand binding.[35−37] This model involves the stabilization of an accessible conformation, and the ligand binding causes a shift in the pre-existing conformational equilibrium. It is important to note that these two models are not mutually exclusive; a survey of the recent literature indicates that many processes can simultaneously include certain elements of both the IF and CS models.[32,33,38−41]

In our simulation, inactive, intermediate, as well as active states coexist, and CO2 can easily get into the binding pocket of the active state but cannot get into the pocket of the inactive states; the ratio between CO2 accessible and inaccessible conformations is about 1:1. The longest continuous periods of CO2 accessible and CO2 inaccessible conformations are all about 20% of total simulation time (see Figure b). Therefore, for the apo form of HCA II, the lasting time of CO2 accessible conformations is long enough for its entry, but the lasting time of CO2 accessible conformations is also not short, and CO2 must “select” a right time to enter HCA II and bind. In other words, the conformational changes of apo-HCA II upon ligand binding favor the CS model as the conformation of the loop 1 is variable, and its open conformations can exist prior to antibiotics binding; when the loop 1 is semiopen or open, antibiotics enter and bind to HCA II; when loop 1 is closed, such an entry is forbidden. The situation for larger ligands is similar, but the lasting time of their accessible conformations is shorter. However, in this process, the entry of CO2 can also induce minor a conformational change of HCA II, so there will still be some minor IF elements. Once CO2 enters the binding pocket, loop 1 will close and compress the pocket so that CO2 is buried into the active site, then the side chains of His119, Leu198, Thr199, and Trp209 are attracted by CO2, rotated toward and bound to CO2.

Methods

MD Simulations

The MD simulations were performed using GROMACS 5.1.4[42,43] and the AMBER99sb-ildn force field.[44] The initial coordinates for the MD simulations were taken from the crystal structure of apo-HCA II with one Zn2+ ion and 9 crystal water molecules (PDB entry 3KS3; Figure ).[23] The protein molecule was placed in a cubic TIP3P water box,[45] with the edges of the water box located at least 10 Å from the protein atoms. The system consisted of 4055 protein atoms and 45 900 solvent atoms, neutralized by 2 Cl– ions under periodic boundary conditions. To remove potential bad contacts, the initial model was subjected to two rounds of energy minimization: first, 1000 steps of steepest descent followed by 1000 steps of conjugate gradient, with the protein atoms constrained by a harmonic potential with a force constant of 100 kcal/mol, and then, 2000 steps of steepest descent and 3000 steps of conjugate gradient without any constraints. Five 200 ns MD trajectories at 300 K with different initial random velocities were obtained, resulting in a total of 1 μs of simulation. In each trajectory, the system was heated at constant volume from 0 to 300 K over 60 ps, and then, a 5 ns equilibration MD simulation and a 200 ns production MD simulations were run at constant temperature and pressure, with the temperature maintained at the temperatures set by the velocity-rescale thermostat[46] with a collision time of 0.2 ps and the pressure maintained at 1.0 atm using the Parrinello–Rahman barostat[47] with the coupling time of 2.0 ps. The long-range electrostatic interaction is calculated by the particle-mesh Ewald method[48] with a real space cutoff of 10 Å; the cutoff distance of the van der Waals interaction is also 10 Å as well. All bonds involving hydrogen atoms were fixed to their equilibrium values using the LINCS algorithm[49] to allow a 2.0 fs time step. The MD trajectories were recorded at 4 ps intervals.

Multiple Trajectories in the Simulations

In our simulations, five 200 ns MD trajectories rather than a single 1000 ns trajectory were generated. The benefits of this approach are as follows: first, five shorter simulations with different initial velocities can more efficiently cover the relevant basin in configuration space and minimize force-field-induced artifacts;[50] second, the five trajectories can be run simultaneously, thereby saving time.

Conformational Clustering and Network Analysis

The MD snapshots were clustered based on the DRMS using the leader algorithm implemented in WORDOM.[51,52] The DRMS between two structures or snapshots (a and b) is defined as follows.where d is the distance between atoms i and j and N is the total number of pairs of atoms. The DRMS was chosen for the conformational clustering because it is sensitive to not only the conformations of individual loops but also the conformational coupling of the loops. The backbone atoms Cα, C, and N of the loops 1 and 2 (residues 197–206 in loop 1 and residues 230–240 in loop 2) were used for the calculations. The DRMS cutoff for the conformational clustering was 1.5 Å. For the MD trajectories, each cluster was treated as a node of a conformational network, with the node size proportional to the number of snapshots in the cluster and the edge thickness proportional to the number of direct transitions between two connected clusters. The visual representations of the conformational networks were created in Pajek.[53]

PCA was performed on the Cα atoms using the “gmx covar” command in GROMACS 5.1.4. The coordinates of each trajectory were superposed by means of a least-squares fit using the first snapshot as a reference. A covariance matrix of atomic fluctuations, C, was then calculated, with each element, cij, for a given pair of atoms defined as followswhere r1, r2, ..., r3 are the Cartesian coordinates of the Cα atoms and the angular brackets denote an ensemble average calculated from all snapshots of the MD trajectory.[54−57] Diagonalization of the covariance matrix yielded the eigenvalues λ, that is, the diagonal elements of the diagonalized matrix, and the eigenvectors L, that is, the columns of the orthonormal transformation matrix. These eigenvectors are the principal modes of motion, and the eigenvalues are the mean-square fluctuations along the eigenvector coordinates. PCA is not restricted to harmonic motions; it can also describe collective transitions between structures that differ greatly.

The percentage v of the total variance contained in a given eigenvector L is given bywhere v measures the contribution of a given PC to the overall protein fluctuations. The cumulative contribution from a set of PCs is

Typically, the first few PCs capture most of the intramolecular fluctuation of a protein.[30,57−63]