J. Phys. Chem. B 2009, 113, 10097–10103
10097
Umbrella Sampling Simulations of Biotin Carboxylase: Is a Structure with an Open ATP Grasp Domain Stable in Solution? Brian R. Novak,*,† Dorel Moldovan,*,† Grover L. Waldrop,*,‡ and Marcio S. de Queiroz*,† Departments of Mechanical Engineering and Biological Sciences, Louisiana State UniVersity, Baton Rouge, Louisiana, 70803, USA ReceiVed: December 3, 2008; ReVised Manuscript ReceiVed: May 20, 2009
Biotin carboxylase is a homodimer that utilizes ATP to carboxylate biotin. Studies of the enzyme using X-ray crystallography revealed a prominent conformational change upon binding ATP. To determine the importance of this closing motion, the potential of mean force with the closure angle as a reaction coordinate was calculated using molecular dynamics simulations and umbrella sampling for a monomer of Escherichia coli biotin carboxylase in water with restraints to simulate attachment to a surface. The result suggests that the most stable state for the enzyme is a closed state different from both the ATP-bound and open state X-ray crystallography structures. There is also a significant motion of a region near the dimer interface not predicted by considering only open and closed configurations, which may have implications for the dynamics and activity of the dimer. Introduction Biotin carboxylase catalyzes the ATP dependent carboxylation of biotin. Biotin carboxylation comprises the first halfreaction of the multifunctional enzyme acetyl-CoA carboxylase, which catalyzes the first committed step in fatty acid synthesis. In eukaryotes, acetyl-CoA carboxylase is a single polypeptide chain, whereas in prokaryotes each function is on a separate protein. Biotin carboxylase from E. coli occurs as a homodimer with each monomer having a complete active site. However, mutation of the dimer interface can lead to monomer formation without significant loss of activity.1 The monomer consists of three domains. Two domains, the A and C, interact to form a large globular section of the protein. Connected to the A and C domains is the smaller, flexible, B domain (also called the ATP grasp domain). See Figure 1. There are X-ray crystallography structures with the B domain in open and closed states. In the crystal structure2 of unliganded E. coli biotin carboxylase, the B domain is in an open state (PDB code 1DV1). In contrast, in the structure2 of an inactive mutant form (PDB code 1DV2) of biotin carboxylase with the substrate ATP or the structure of wild-type enzyme with the ATP analog AMPPNP3 bound, the B domain is in the closed down position. Even though an open structure forms in the solid phase, this structure may not be stable in solution. Computational studies provide a method to determine the enzyme’s stability relative to open or closed configurations. This information is very useful for confirming or complementing X-ray crystallography studies. For example, it can resolve the issue of whether the open structure of the B domain is stable or not. There have been few computational studies of biotin carboxylase to date. Lill et al.4 used molecular dynamics simulations to study the physical bonding of Mg2ATP in the active site and performed a covariance analysis to determine hinge regions in the enzyme. Bordelon et al.5 studied physical * To whom correspondence should be addressed. E-mail: bnovak1@ lsu.edu (B.R.N.),
[email protected] (D.M.),
[email protected] (M.S.d.Q),
[email protected] (G.L.W.). † Department of Mechanical Engineering. ‡ Department of Biological Sciences.
Figure 1. Domains, hinge regions, restrained atoms, and angle for umbrella sampling. The domains from FIRST are shown as blue and orange spheres. The backbone atoms of the hinge residues (Cys130, Val131, Glu205, and Asp206) from DynDom are shown as green spheres. The backbone of the domains from DynDom are shown in cyan and yellow, and the hinge axis is shown in purple. DynDom domain 1 corresponds to domains A and C, and DynDom domain 2 corresponds to domain B. The restrained atoms are shown as red spheres. The atoms used to define the angle are shown as black spheres connected by dashed lines and labeled as atoms 1-3. The angle θ is represented by a pink arc.
bonding of Mg2ATP in the active site when Gly165 or Gly166 were mutated using molecular dynamics simulations. Sueda et al.6 studied binding of ATP starting from the open X-ray structure using molecular dynamics simulations. They also docked biotin to the closed, ATP bound state. Ito et al.7 studied the reaction mechanism using density functional theory. The free energy difference between states of a system indicates their relative stability. The state corresponding to the global minimum in free energy is the most stable, and a system will be in this state with the highest probability unless it becomes trapped in a local minimum due to high energy barriers. At a given temperature, the free energy is a function of all of the positions and momenta in the system. Calculation of the free energy as a function of all of these variables is practically very
10.1021/jp810650q CCC: $40.75 2009 American Chemical Society Published on Web 07/08/2009
10098
J. Phys. Chem. B, Vol. 113, No. 30, 2009
TABLE 1: Values of θ0 Used for Biasing (eq 4), Equilibration Time, and Force Constants for Each Umbrella Sampling Window θ0 (°)
kθ (kJ/mol)
equilibration time (ns)
67.0 72.0 77.0 82.0 85.0 90.0 93.0 94.0 99.2 104.4 109.6 114.8 120.0 124.0 129.0 134.0
1500 1500 1500 1500 1500 1000 3000 1000 1000 1000 1000 1000 1000 1500 1500 1500
23.0 11.0 10.0 8.0 12.0 16.0 18.0 8.0 10.0 8.0 15.0 17.0 17.0 12.0 12.0 8.0
difficult. To simplify things, everything can be integrated out except one or a few internal coordinates (reaction coordinates). The goal is to choose coordinates that can approximately represent the behavior of the system between initial and final states and capture the major features of the energy landscape. Integrating everything except the reaction coordinates leaves a free energy profile that is also called the potential of mean force (PMF) introduced by Kirkwood.8 The PMF difference between two values of a reaction coordinate is related to the relative probability of observing those values of the reaction coordinate.9 This requires sampling of all values of the reaction coordinate between those two values. Since, for our system, the only experimental data are the X-ray crystallography structures for open and closed states, simulations can be used to sample states in between them. Direct calculation of the PMF difference between open and closed states can be performed in a molecular dynamics or Monte Carlo simulation only if the enzyme can frequently open and close on a time scale of nanoseconds. If there are significant barriers to opening or closing, there will be poor or no sampling near those barrier regions (transition states). There are several ways to indirectly calculate the PMF so that sampling near transition states is sufficient. These include thermodynamic integration using constrained molecular dynamics,10,11 umbrella sampling,12 Jarzynski’s equality13,14 and steered molecular dynamics,15 and adaptive force biasing.16 In this work, umbrella sampling is used. In umbrella sampling, bias potentials are imposed on the system
Novak et al. to keep it near various values of the reaction coordinate(s) (windows) to provide more uniform sampling over the desired range of the reaction coordinate(s). This yields a biased probability distribution for each bias potential. These are combined to obtain the unbiased probability and the PMF. The goal of this work is to use computational methods to elucidate some of the mechanistic properties of biotin carboxylase. First, we use umbrella sampling and molecular dynamics simulations to determine the PMF of a biotin carboxylase monomer in solution with a closure angle as the reaction coordinate. This PMF indicates the relative stability of open and closed states of the enzyme and the size of barriers to opening or closing. The results show that a closed state of the monomer is more stable than an open configuration corresponding to the dimer crystal structure, but the probability of the monomer being open is only about 16 times less than the probability of being closed. From the umbrella sampling simulations, we also calculate the Euler angles for rotation of the plane defined by the points used to determine the closure angle. The Euler angles are a measure of the motion of the B domain, since only it undergoes significant conformational change. Besides the closure motion (roll angle), the motion consists of a small, nearly monotonic change in the twist angle, with no significant tilt motion. Finally, there is significant motion of the A and C domains relative to each other as a function of closure angle during the umbrella sampling simulations that may be significant for the dimer form of biotin carboxylase. Methods Simulation Details. A monomer of E. coli biotin carboxylase was simulated using GROMACS 3.3.217 with the GROMOS 43a1 force field. The starting structure was the mutant form closed down on an ATP molecule with the mutated residue changed back to the wild type4 and the ATP removed. Three extra residues at the N-terminus were removed, and two missing ones at the C-terminus were added. GROMACS position restraints with force constants of 1000 kJ mol-1 nm-2 were used on 11 atoms of 3 residues (Met1: Cβ, Cγ, Sδ, Cε; Lys27: Cγ, Cδ, Cε, Nζ; Glu71: Cδ, Oε1, Oε2) that could be modified for bonding to a nickel surface. They will be referred to as the nickel restraints. Extended system equations of motion were used to maintain constant temperature or pressure. A constant average temperature of 300 K was maintained using Nose-Hoover thermostats.18,19 The enzyme and solvent were thermostatted separately with a time constant of 0.1 ps for each. For constant pressure
Figure 2. The plot on the left is the PMF difference divided by RT (2.49 kJ/mol) calculated using umbrella integration with θ divided into 250 bins. The reference is the global minimum at around θ ) 71°. The solid line is the average of three, 6 ns segments of data taken every 100 time steps. The dashed lines are two times the standard deviation of the mean away from the solid curve. The plot on the right is the relative probability density corresponding to the average PMF calculated using eq 1. The inset is a blow up of θ > 82°.
Umbrella Sampling Simulations of Biotin Carboxylase
J. Phys. Chem. B, Vol. 113, No. 30, 2009 10099
Figure 3. (a) A comparison of closed states from umbrella sampling at θ ) 70° with the backbone shown as a cyan ribbon and 1DV2 with the backbone shown as a red ribbon and bound ATP shown in green. The structures were superimposed by least-squares fitting of the positions of the backbone atoms in the A domains. (c) With bound ATP, the flexible BC P-loop is attracted to the phosphate groups. Residues 164-166 are closely associated with the phosphates and, in this configuration, Gly166 is hydrogen bonded to one of the phosphate oxygens. (b) In the absence of ATP, it is attracted to an aspartic acid residue (Asp115) at the junction between two helices. Residues Gly164, Gly163, and Ser161 are hydrogen bonded to Asp115. This shifts the configuration of the entire B domain. In the regions showing atomic details, carbon is green, oxygen is red, nitrogen is blue, hydrogen is white, and phosphorus is tan.
Figure 4. The roll, tilt, and twist angles. The upper right figure shows the definition of the axes relative to the three atoms used to define θ. Atom 1 is in the yz-plane. The upper left figure defines the rotations. The lower figures show 20° rotations of the B domain starting from a configuration with θ ) 105° made by applying the appropriate rotation matrix to the B domain coordinates. They correspond to rotation around the x-axis, which is the vector normal to the plane defined by the same atoms used to define θ (roll), the y-axis (tilt), and the z-axis, which is the vector from atom 2 to atom 3 in Figure 1 (twist). The A and C domains are shown in cyan, the original B domain is shown in green, the rotated B domain is shown in red, and the hinges are not shown.
simulation, the Parinello-Rahman barostat20,21 with isotropic scaling was used. The time constant and the compressibility were 1.0 ps and 4.5 × 10-5 bar-1, respectively. The set point pressure was 1.0 bar.
All bonds were constrained; LINCS22 was used on the enzyme bonds, and SETTLE23 was used with the SPC24 water model. The particle mesh Ewald algorithm25,26 was used for evaluation of the long-range electrostatic interactions with a 0.12 nm grid spacing. GROMACS uses a twin range cutoff for pair interactions. The short-range van der Waals and Coulomb cutoff was 1.0 nm, and the long-range van der Waals cutoff was 1.4 nm. Interactions for pairs with distances between the short and longrange cutoffs were evaluated every 5 time steps. System Preparation. Preparation consisted of adding internal water molecules using DOWSER,27 steepest descent energy minimization, addition of water molecules and a sodium ion to neutralize the negative charge on the enzyme, and additional steepest descent energy minimization. Two iterations of the DOWSER program were run, and 24 waters in addition to those in the X-ray crystallography structure were added. Some of these filled the spot left by removing the ATP. The first round of energy minimization was on the enzyme plus waters from the X-ray crystallography structure plus those added by DOWSER. A 0.01 nm maximum step size and a 10 kJ mol-1 nm-1 force tolerance were used. The enzyme backbone was restrained using position restraints with 1000 kJ mol-1 nm-2 force constants, and the water bonds and angles were made flexible. Then, 14 137 waters and the sodium ion were added in a rectangular box large enough to give at least 0.8 nm of water between the enzyme and the box edge on any side or 1.6 nm between the enzyme and its images. This was followed by the same energy minimization procedure with the additional water.
10100
J. Phys. Chem. B, Vol. 113, No. 30, 2009
Novak et al.
Figure 5. Roll, tilt, and twist angles as a function of the biasing angle (θ0) for each window. The bars are two standard deviations from the mean to show the range of values. Error bars, which are related to the standard deviations of the mean values, are not shown and are much smaller.
Figure 7. DynDom analysis on the A and C domains for θ ) 80° (green and blue) and 110° (yellow and red) from umbrella sampling simulations. The view is down the rotation axis shown as a black X. The blue region rotates 18° and the C-terminus changes location when θ goes from 80 to 110°.
Figure 6. Snapshots from the umbrella sampling simulations showing the enzyme backbone with the B domain colored cyan at various angles. Note the opening and closing of the A and C domains, which was not originally predicted by DynDom.
Initial equilibration consisted of several steps. Three 100 ps simulations were ran with various 1000 kJ mol-1 nm-2 position restraints under isobaric and isothermal conditions. The first was with the backbone atoms restrained. Then, only the N and C terminal atoms were restrained. Finally, only the nickel restraints were used. The volume was then averaged over a 25 ps run at
constant temperature and pressure. Finally, a 6 ns simulation was run at constant volume and temperature. Domain Analysis. Since the closure angle was used as the reaction coordinate, the GROMACS angle restraint potential was used as the bias potential for umbrella sampling. This potential applies to only atoms, not groups of atoms. A good choice of the three atoms to use for definition of the angle was one of the restrained atoms which are in the A domain, one in the region that acts as a hinge for closure, and one in a relatively rigid region of the B domain. The restrained atom and the hinge atom did not move very much and the third atom was in the middle of a rigid region, so the bias potential would not distort the structure. Determining the rigid atom and the hinge atom to use required that two types of domain analysis be performed. The floppy inclusions and rigid substructure topography (FIRST)28 program was used to determine the rigid regions of the enzyme. This program uses covalent bonding, hydrogen bonding, and hydrophobic tethers along with an algorithm
Umbrella Sampling Simulations of Biotin Carboxylase
J. Phys. Chem. B, Vol. 113, No. 30, 2009 10101
known as the pebble game to determine rigidity. The cutoff energy for hydrogen bonds can be varied. Since there are X-ray crystallography structures for both an open (1DV1) and a closed configuration (1DV2), the DynDom program29,30 was used. This program gives the angle of rotation, hinge residues, and hinge axes as well as the domains involved in the motion. PMF and Umbrella Sampling. In the canonical ensemble, the PMF (A(ξ)) along a reaction coordinate ξ is related to the probability (P) of observing particular values of ξ by
[
A(ξ) ) -kBT ln
∫ δ[ξ - ξ(r)] exp[-H(r)/kBT] dr ∫ exp[-H(r)/kBT] dr
]
+
F ) -kBT ln[P(ξ)] + F (1)
In this equation, δ is the Dirac delta function, H is the total energy of the system, r represents all coordinates for the system, F is an undetermined constant, kB is Boltzmann’s constant, and T is temperature. The delta function is 1 when ξ(r) ) ξ and is 0 otherwise. All coordinates except ξ are integrated out. Equation 1 does not include a bias potential, so there will be poor or no sampling near transition states. The expression corresponding to eq 1 when a bias potential (Ubias) is applied along the reaction coordinate is
A(ξ) ) -kBT ln Pbias(ξ) - Ubias(ξ) + F
(2)
where Pbias is the biased probability.12 Equation 2 applies to all windows with a different Ubias and calculated Pbias for each. There are different ways to combine the probability distributions, including the weighted histogram analysis method31 and the umbrella integration method.32 In this study, umbrella integration was used. This method is based on the assumption that the biased probability distribution in each window is a normal distribution. Since Pbias(ξ) for each window is then an analytical function of ξ, the PMF is also an analytical function of ξ. The derivative of A with respect to ξ for each window is then obtained analytically,
∂Ubias ∂A ξ - ξ¯ ) kBT 2 ∂ξ ∂ξ σ
(3)
where ξj is the mean value of ξ for the window, and σ is the standard deviation of ξ for the window. The whole range of the reaction coordinate is then divided into bins, and the derivative for each bin is calculated from a weighted average from all windows, since the windows can overlap and multiple windows can contribute to a bin. Physically, the result is the equivalent of the force along the ξ direction required to maintain the system at a certain value of ξ. The PMF was obtained by numerical integration of this force using the trapezoidal rule. Results and Discussion Domain Analysis. Since only three atoms were used to define a closure angle and one of those needed to be in the B domain, it was desirable to choose an atom in the middle of a relatively rigid region of the B domain to reduce the potential for artificial distortion of the enzyme structure by the force due to the bias potential. To find the rigid regions, the FIRST program was
run with a hydrogen bond cutoff energy of -1.6736 kJ/mol and the default rule for hydrophobic tethers on 50 equilibrated closed configurations taken every 10 ps from a 500 ps trajectory. The rigid atoms that were conserved in all 50 runs yielded 3 domains. Two of these were single helices in the B domain that are shown in Figure 1, and the third one was in the A and C domains and is not indicated in the figure. The hinge regions of the enzyme were found using the DynDom program using open and closed configurations in order to choose an atom from those regions for use in defining the closure angle. DynDom also gives the angle(s) of rotation and hinge axes. The open X-ray crystallography structure of biotin carboxylase has electron density missing in the B domain due to flexibility. Therefore, an open configuration was generated from the closed one by pulling the closed structure open. To determine approximately how to pull the structure open, an initial hinge axis was calculated using DynDom with the equilibrated closed configurations and the incomplete open X-ray crystallography structure (1DV1). Then, the closed configuration was opened by pulling around this hinge axis at a rate of 0.01 nm/ps using steered molecular dynamics.15 A small linear correction was also needed to get domain 2 in Figure 1 to approximately line up with the one in the open structure. This new open configuration was then equilibrated for 6 ns. DynDom was then run again using the equilibrated open and equilibrated closed configurations. The result was two domains with an angle of rotation of 46.5°. The rotation axis and hinge residues are shown in Figure 1. The hinge is the same as one of those found by Lill et al.4 using HingeMaster.33 Note that changing the two configurations used in DynDom yields slightly different hinge residues. Instead of just 2 pairs, there are probably 8-10 residues involved in the hinge motion. The motion was classified as 92.8% closure motion and 7.2% twist motion. Closure motion is rotation around the axes perpendicular to the line connecting the centers of mass of the two domains, and twist motion is rotation around the axis parallel to this line. Umbrella Sampling. The results from FIRST and DynDom were used to define an angle (θ) to use as the reaction coordinate for umbrella sampling. The angle was defined between the C (backbone carbonyl carbon) atom in Met184, which is in domain 2, the C atom in Val131 which is in the hinge region, and the Cδ atom in Glu71, which is a restrained atom. Figure 1 shows the angle and atoms. The GROMACS angle restraint potential that was used for the bias potential is
Ubias(θ) ) kθ(1 - cos(θ - θ0))
(4)
Sixteen windows were used, and the initial configuration for a window was taken from an adjacent one, except for the first one at 120°. Each initial configuration was then equilibrated with the new bias potential before sampling. The bias potential parameters and equilibration times are shown in Table 1. The Mann-Kendall test for trend34 for the mean and standard deviations of segments of data was used to determine if a window was equilibrated. This requires that uncorrelated data be used. The von Neumann test for serial correlation34 indicated that the correlation length was 150-250 ps, and the angle autocorrelation function indicated that it was at least 250300 ps. The segment length for the Mann-Kendall test was 300 ps. The PMF is shown in Figure 2. A closed state at θ ) 71° is the global minimum of the PMF. At most, there are only shallow minima for an open state around 4.0 RT. However the relatively large error in the PMF around these states is too large for them
10102
J. Phys. Chem. B, Vol. 113, No. 30, 2009
to be characterized. They are not very significant, since the barriers are likely less than the thermal energy. This means that, before binding, the enzyme is found most often in a closed state. The global minimum is significantly different from the closed configuration with bound ATP (1DV2) despite having similar values of θ (75.9° in the 1DV2 X-ray structure). This is due to the flexible BC P-loop3 that is attracted to the phosphate groups in ATP being instead attracted to the Asp115 residue located at the junction of two helices. See Figure 3. The fact that a configuration that is similar to the closed one with bound ATP was not sampled means that it is either at a higher energy than the observed configurations or there is an energy barrier between the two. In either case, another reaction coordinate(s) would be needed. The relative probability density of the enzyme having a configuration corresponding to a particular value of θ is also shown in Figure 2. If open is defined as θ g 82° and closed as 60° e θ < 82°, then comparing the areas under the probability density curve shows that the open state is about 16 times less likely to occur than the closed state. Structural Changes Along the Reaction Coordinate. Although only one reaction coordinate could practically be studied using umbrella sampling, other changes in the enzyme were observed as it closed. To study these structural changes, first, the Euler angles of the B domain were calculated as a function of θ. Next, motion of the A and C domains was studied since they form the interface between monomers in the homodimer, which could have implications for the dynamics of the dimer. Conformational Changes of the B Domain. The DynDom program gave an estimate of the deviation of the motion of the B domain from simple rotation around an axis, 7.2% twist. However, this did not give any information about the motion of the B domain as it closed, since it only compared the open and closed states. To quantify the motion of the B domain, Euler angles were calculated as a function of closure angle. These angles were calculated for rotation of the plane defined by the atoms used to define θ (see Figures 1 and 4) for each window to give an idea of how the closure occurs. There are many possible conventions for calculation of the Euler angles. The same convention35 that is used in the 3DNA program36 was followed. It is illustrated in Figure 4. The x-axis was defined as the unit vector normal to the plane defined by the positions of the atoms that define the angle θ. The z-axis was defined as the normalized vector from atom 2 to atom 3, and the y-axis was orthogonal to both of them. To define the angles, a set of reference axes (xref, yref, zref) must be chosen. The roll angle is rotation about the xref-axis, the tilt angle around the yref-axis, and the twist angle around the zref-axis. The reference axes were defined relative to the first data point at θ0 ) 67°. The closure angle (θ) is strongly correlated with the roll angle, so a large change in roll was expected. The amount of tilt and twist and how they changed as θ changed were of most interest. Figure 5 shows the average roll, tilt, and twist angles for each window. The roll angle increases as expected, and changes by about 35° as θ0 changes from 104.4 to 67°. Beyond 104.4° the PMF increases rapidly. The tilt angle shows essentially no change over this range of interest; it only changes for θ0 > 104.4°, and some of this change is due to rotation of the whole enzyme, which was not taken into account since it is only significant for tilt due to the restraints. The twist angle changes by about 20° as the enzyme closes. Another way to gauge the amount of motion that occurs outside the original plane is to calculate the percentage of the surface of a unit sphere that an axis samples. An upper bound
Novak et al. on this can be found by integrating over the range of the angles covered in spherical coordinates. This percentage was only 7.2% for the z-axis, and would be smaller if overall rotation was removed. In summary, the motion of the B domain as a whole from 104.4 to 67.0° is fairly simple. It involves almost no tilt and a small amount of twist motion compared to roll motion. The twist motion is nearly monotonic except for the dip at a θ0 of 99.2°. The motion is qualitatively similar to the result of the DynDom calculation, which predicted mostly closure motion (primarily roll in this case) with a small amount of twist motion. If two reaction coordinates were used for the motion, roll and twist would be reasonable choices. Changes in the A and C Domains. Figure 6 shows configurations at several different values of θ from the umbrella sampling simulations. Notice that the A and C domains (DynDom domain 1 defined in Figure 1) open and then close again as the B domain closes. This motion was not predicted by DynDom since these domains have roughly the same structure in the open and closed configurations that were used. The motion is possible because the A and C domains are not covalently bonded. An analysis of the A and C domain motion was performed using DynDom between the configurations shown in Figure 3 at θ ) 80° and θ ) 110°. The result is shown in Figure 7. A small part of the C domain rotates by 18° and the C-terminus changes from pointing toward the A domain to pointing away from it. The motion of the A and C domains has implications for the function of the homodimer form of the enzyme since the contact between monomers occurs in this domain. A motion similar to this may be required for the closure of one side of the dimer as well, although it is likely inhibited by the presence of the other monomer. It is thought that a reaction only occurs in one-half of the dimer at a time,37 and it has been observed that mutation of the active site of only one monomer can substantially reduce the overall activity.38 If a large motion of the connecting region is required for binding the substrates in the first monomer, then correct binding by the second monomer may be prevented by a conformational change due to binding by the first monomer. Conclusions The results of these simulations of the unliganded biotin carboxylase monomer in solution suggest that the open state observed in the X-ray crystallography structure of the dimer is a consequence of the crystal contacts that atoms in the B domain form with symmetry related molecules or that the monomer behaves differently than the dimer due to the lack of specific monomer-monomer interactions. Hinge analysis programs only give an estimate of the overall motion between states, and not what happens in between. Euler angles for rotation of the plane defined by the points used to determine the closure angle were calculated from the umbrella sampling windows as a function of closure angle to further quantify the motion of the B domain. The motion is relatively simple. Besides the closure motion (roll angle), the motion consists of a small, nearly monotonic change in the twist angle, with no significant tilt motion. A study of adenylate kinase using various experimental and computational techniques39 showed that even in the absence of ligands, the enzyme sampled a minor state similar to the catalytically competent ligand bound state. Although a state
Umbrella Sampling Simulations of Biotin Carboxylase similar to the 1DV2 X-ray structure was not sampled in the umbrella sampling simulations, this does not mean that it would not be sampled if an additional reaction coordinate(s) were considered. For the value of the closure angle corresponding to the 1DV2 structure there is simply another configuration at a lower energy. A structure similar to 1DV2 might be rarely sampled, but is not seen on molecular dynamics time scales without adding a bias along another coordinate. The lack of a significant energy minimum for an open state raises two questions. Does the B domain open for substrates to enter in the way suggested by the X-ray crystallography structures, or can substrates come in from the side of the closed conformation accompanied by a different conformational change? If the B domain opens, by how much? Answering these questions will be the subject of future investigation. The motion in the A and C domains may be important for the dynamics and function of the dimer since it affects the subunit interface in the homodimer. Future work will simulate the dimer to help elucidate the communication mechanism between the monomers. The probability of being open or closed could be measured experimentally to obtain a qualitative comparison with the potential of mean force results calculated here, using single molecule fluorescence resonance energy transfer (FRET). In FRET, two different types of dye molecules are attached to two different locations in a molecule. After being excited, the donor molecule transfers energy to the acceptor molecule through induced-dipole-induced-dipole interactions. The magnitude of this transfer is a function of distance between the donor and acceptor. The excited donor and acceptor then relax by fluorescence. The relative intensity of fluorescence of the dye molecules is therefore a function of their distance (primarily due to the donor at large distances and mostly due to the acceptor at short distances), so the probability as a function of dye molecule distance can be calculated. For biotin carboxylase, the dye molecules would be attached to the B domain and either the A or C domain. In fact, FRET studies of adenylate kinase exhibited similar results to the computational results found here with biotin carboxylase.40 Adenylate kinase undergoes large conformational changes upon binding ATP and AMP. In the absence of substrates, the probability of being in the closed conformation was greater than being in the open conformation with little or no barrier to closure. Acknowledgment. This work was supported in part by the LSU Innovation in Engineering Research Fund and the LSU Department of Mechanical Engineering. Computational resources were provided by HPC @ LSU and the Louisiana Optical Network Iniative. References and Notes (1) Shen, Y.; Chou, C. Y.; Chang, G. G.; Tong, L. Mol. Cell 2006, 22, 807–818.
J. Phys. Chem. B, Vol. 113, No. 30, 2009 10103 (2) Thoden, J. B.; Blanchard, C. Z.; Holden, H. M.; Waldrop, G. L. J. Biol. Chem. 2000, 275, 16183–16190. (3) Mochalkin, I.; Mill, J. R.; Evdokimov, A.; Lightle, S.; Yan, C.; Stover, C. K.; Waldrop, G. L. Protein Sci. 2008, 17, 1706–1718. (4) Lill, S. O. N.; Gao, J. L.; Waldrop, G. L. J. Phys. Chem. B 2008, 112, 3149–3156. (5) Bordelon, T.; Lill, S. O. N.; Waldrop, G. L. Proteins 2009, 74, 808–819. (6) Sueda, S.; Ito, Y.; Yoshida, T.; Kondo, H. Seikagaku 2004, 76, 777. (7) Ito, Y.; Kondo, H.; Shiota, Y.; Yoshizawa, K. J. Chem. Theory Comput. 2008, 4, 366–374. (8) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300–313. (9) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: Oxford, 1987. (10) den Otter, W. K.; Briels, W. J. J. Chem. Phys. 1998, 109, 4139– 4146. (11) Sprik, M.; Ciccotti, G. J. Chem. Phys. 1998, 109, 7737–7744. (12) Torrie, G. M.; Valleau, J. P. Chem. Phys. Lett. 1974, 28, 578–581. (13) Jarzynski, C. Phys. ReV. Lett. 1997, 78, 2690–2693. (14) Jarzynski, C. Phys. ReV. E 1997, 56, 5018–5035. (15) Isralewitz, B.; Baudry, J.; Gullingsrud, J.; Kosztin, D.; Schulten, K. J. Mol. Graph. Model. 2001, 19, 13–25. (16) Darve, E.; Pohorille, A. J. Chem. Phys. 2001, 115, 9169–9183. (17) Van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701–1718. (18) Nose´, S. J. Chem. Phys. 1984, 81, 511–519. (19) Hoover, W. G. Phys. ReV. A 1985, 31, 1695–1697. (20) Nose´, S.; Klein, M. L. Mol. Phys. 1983, 50, 1055–1076. (21) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182–7190. (22) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. J. Comput. Chem. 1997, 18, 1463–1472. (23) Miyamoto, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 952–962. (24) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. In Intermolecular Forces; Pullman, B., Ed.; D. Reidel Publishing Company: Dordrecht, 1981; pp 331-342. (25) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089– 10092. (26) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577–8593. (27) Zhang, L.; Hermans, J. Proteins 1996, 24, 433–438. (28) Jacobs, D. J.; Rader, A. J.; Kuhn, L. A.; Thorpe, M. F. Proteins 2001, 44, 150–165. (29) Hayward, S.; Kitao, A.; Berendsen, H. J. C. Proteins 1997, 27, 425–437. (30) Hayward, S.; Berendsen, H. J. C. Proteins 1998, 30, 144–154. (31) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. J. Comput. Chem. 1992, 13, 1011–1021. (32) Ka¨stner, J.; Thiel, W. J. Chem. Phys. 2005, 123, 144104. (33) Flores, S.; Echols, N.; Milburn, D.; Hespenheide, B.; Keating, K.; Lu, J.; Wells, S.; Yu, E. Z.; Thorpe, M.; Gerstein, M. Nucleic Acids Res. 2006, 34, D296-D301. (34) Schiferl, S. K.; Wallace, D. C. J. Chem. Phys. 1985, 83, 5203– 5209. (35) Gonzalez, O.; Maddocks, J. H. Theor. Chem. Acc. 2001, 106, 76–82. (36) Lu, X. J.; Olson, W. K. Nucleic Acids Res. 2003, 31, 5108–5121. (37) de Queiroz, M. S.; Waldrop, G. L. J. Theor. Biol. 2007, 246, 167– 175. (38) Janiyani, K.; Bordelon, T.; Waldrop, G. L.; Cronan, J. E. J. Biol. Chem. 2001, 276, 29864–29870. (39) Henzler-Wildman, K. A.; Thai, V.; Lei, M.; Ott, M.; Wolf-Watz, M.; Fenn, T.; Pozharski, E.; Wilson, M. A.; Petsko, G. A.; Karplus, M.; Hubner, C. G.; Kern, D. Nature 2007, 450, 838–844. (40) Hanson, J. A.; Duclerstacit, K.; Watkins, L. P.; Bhattacharyya, S.; Brokaw, J.; Chu, J. W.; Yang, H. P. Natl. Acad. Sci. USA 2007, 104, 18055–18060.
JP810650Q
