extfieldsofstudy: Biology Medicine

s2fieldsofstudy: Chemistry Biology

text:
Ligand-induced Conformational Changes in the Capping Subdomain of a Bacterial Old Yellow Enzyme Homologue and Conserved Sequence Fingerprints Provide New Insights into Substrate Binding*

We have recently reported that Shewanella oneidensis, a Gram-negative γ-proteobacterium with a rich arsenal of redox proteins, possesses four old yellow enzyme (OYE) homologues. Here, we report a series of high resolution crystal structures for one of these OYEs, Shewanella yellow enzyme 1 (SYE1), in its oxidized form at 1.4Å resolution, which binds a molecule of PEG 400 in the active site, and in its NADH-reduced and p-hydroxybenzaldehyde- and p-hydroxyacetophenone-bound forms at 1.7Å resolution. Although the overall structure of SYE1 reveals a monomeric enzyme based on the α8β8 barrel scaffold observed for other OYEs, the active site exhibits a unique combination of features: a strongly butterfly-bent FMN cofactor both in the oxidized and NADH-reduced forms, a collapsed and narrow active site tunnel, and a novel combination of conserved residues involved in the binding of phenolic ligands. Furthermore, we identify a second p-hydroxybenzaldehyde-binding site in a hydrophobic cleft next to the entry of the active site tunnel in the capping subdomain, formed by a restructuring of Loop 3 to an “open” conformation. This constitutes the first evidence to date for the entire family of OYEs that Loop 3 may indeed play a dynamic role in ligand binding and thus provides insights into the elusive NADH complex and into substrate binding in general. Structure-based sequence alignments indicate that the novelties we observe in SYE1 are supported by conserved residues in a number of structurally uncharacterized OYEs from the β- and γ-proteobacteria, suggesting that SYE1 represents a new subfamily of bacterial OYEs.

The OYE family has been studied since the 1930s when Warburg and Christian isolated a yellow enzyme from brewers' bottom yeast (Saccharomyces carlsbergensis) (1). Theorell purified the enzyme in 1935 and showed that it was composed of a colorless apoprotein and a yellow dye. The dye turned out to be a vitamin B 2 -derived molecule, FMN, and to be essential for catalysis. When a second yellow enzyme was isolated from yeast, the first enzyme was named old yellow enzyme (OYE), 3 a name that has persisted till this day. Since then, OYE has become synonymous with the study of flavoproteins, but despite decades of research culminating with a dissection of the catalytic mechanism (2)(3)(4), the physiological role of OYE has only recently begun to emerge (5)(6)(7).
The OYE family has grown steadily in recent years and now features the bacterial nitro-ester reductases PETN reductase (8), glycerol trinitrate reductase (9), XenA/B reductase (10), the bacterial morphinone reductase (11), the YqjM from Bacillus subtilis (41), the plant oxophytodienoic acid reductases (12,13), several yeast OYEs (14 -17), and an enzyme involved in prostaglandin synthesis in Trypanosoma cruzi (18). Although these enzymes originate from different organisms and catalyze different reactions, they share several common functional characteristics. They are able to reduce simple and complex unsaturated aldehydes and ketones (11,19), nitro-esters (8 -10, 20), and nitro-aromatic substrates (21)(22)(23), and yet clear preferences are observed. Furthermore, OYEs can form long wavelength charge transfer interactions with phenolic compounds (24), which in turn exhibit a distinct binding mode in the active site involving hydrogen bonding of the phenolate oxygen to a strictly conserved histidine/asparagine or histidine/histidine pair (25). The reaction mechanism of OYE proceeds by a ping-pong mechanism consisting of an oxidative and reductive half-reaction. Consistent with the plethora of functional similarities among members of the OYE family, the first crystal structures reported revealed a common ␣ 8 ␤ 8 barrel fold with significant variations in the capping subdomain (25)(26)(27). Nonetheless, interesting differences have begun to emerge in recent studies. For example, morphinone reductase lacks the general acid (tyrosine) that protonates the substrate (28). In another case, enzymes possessing two histidines (23) in the active site for substrate binding can catalyze the reduction of TNT to the H-Meisenheimer complex, whereas the enzymes with only one histidine cannot. Finally, YqjM has a flavin that is butterfly-bent in both the reduced and the oxidized state (29).
We have recently probed the genome of Shewanella oneidensis, a Gram-negative ␥-proteobacterium with a remarkable respiratory versatility (30) (www.tigr.org), using yeast OYE1. We identified four OYE homologues (NP718044, NP718043, NP719682, and NP718946) that we named SYE1-4 for Shewanella yellow enzyme (31). SYE1 (NP718044) is the Shewanella OYE homologue with the greatest sequence identity to OYE1. Here, we report structural studies of SYE1 at high resolution and show that SYE1 exhibits a number of novel features that can be mapped to strongly conserved sequence fingerprints in several structurally uncharacterized OYE homologues from the ␤and ␥-proteobacteria. This suggests that SYE1 may be the prototype of a new subfamily of bacterial OYEs.

EXPERIMENTAL PROCEDURES
Materials-The Q-Sepharose FF and Source 30Q matrices and the MonoQ and Superdex 200 columns were purchased from Amersham Biosciences, and the pACYC-Duet-1 vector was from Novagen. Structure screens 1 and 2 were from Molecular Dimensions Ltd., and the additive screen was from Hampton Research. Trizma base, ammonium sulfate, FMN, MES, Na 2 HPO 4 , and NaH 2 PO 4 were purchased from Sigma; HEPES and PEG 400 were from Fluka; NaCl and NADH were from Merck; t-hexenal, t-decenal, and t-dodecenal were from Acros; and isopropyl ␤-D-thiogalactopyranoside was from Duchefa.
Cloning, Expression, and Purification-The sye1 gene was excised from a pGEX-SYE1 construct using EcoRI and NotI and cloned in the first multiple cloning site of pACYC-Duet-1, generating pACYC-SYE1. This cloning strategy puts the sye1 gene out of frame with the His tag. The pACYC-SYE1 vector was introduced into Escherichia coli BL21(DE3) cells. The cultures were grown at 28°C under constant shaking, and SYE1 expression was induced at 0.5 A 600 nm using 0.5 mM isopropyl ␤-Dthiogalactopyranoside. The cultures were left to grow for another 7 h, and the cells were harvested by centrifugation and resuspended using 30 ml of 50 mM Tris, pH 8.0, per liter of culture grown. The resuspended cells were lysed by sonication, and the soluble fraction was clarified by centrifugation. SYE1 was purified in four steps. Lysate containing SYE1 was manually loaded onto a Q-Sepharose FF column equilibrated with 50 mM Tris, pH 8.0. The column was washed with 50 mM Tris, pH 8.0, 50 mM NaCl, and SYE1 was eluted with 50 mM Tris, pH 8.0, 300 mM NaCl. The eluate was dialyzed to remove the salt and was loaded onto a Source 30Q column equilibrated with 50 mM Tris, pH 8.0. SYE1 eluted at 80 -100 mM NaCl. The pooled fractions were dialyzed against 20 mM MES, pH 6.5, and loaded onto a MonoQ column from which SYE1 eluted at 40 -60 mM NaCl. As a final step the SYE1 fractions were pooled, concentrated, and dialyzed against 50 mM Tris, pH 8.0, 100 mM NaCl and loaded onto a Superdex 200 column. The elution profile of SYE1 was consistent with a 40-kDa protein, indicating that SYE1 was a monomer. Final purity was confirmed by silver staining of an SDS-PAGE gel. The pure fractions were pooled and dialyzed against 35 mM Tris, pH 8.0, 70 mM NaCl. The protein was subsequently concentrated to 12 mg/ml and stored in 100-l aliquots at Ϫ80°C.
Crystallization of SYE1-Crystallization trials were set up by the hanging drop vapor diffusion method, using 200 l of reservoir solutions (Structure Screens 1 and 2; Molecular Dimensions Ltd.) and mixing equal volumes of protein solution and reservoir (1 l ϩ 1 l) to make the droplets. Within 2 days needle-like crystals were observed in condition 30 (2% PEG 400, 2 M ammonium sulfate and 0.1 M HEPES, pH 7.5). Optimization of the crystals by varying the pH and the concentrations of PEG 400 and ammonium sulfate was unsuccessful. We therefore screened a series of additives (Additive screen; Molecular Dimensions Ltd.) against the condition with the best needles (100 mM Tris, pH 8.5, 1.8 M ammonium sulfate, 2% PEG 400). Substantially improved crystals were observed in the presence of 0.5% ␤-octyl glucoside. Upon further optimization, diffraction quality crystals were grown in droplets containing 100 mM Tris, pH 8.0 -8.5, 1.75 M ammonium sulfate, 2% PEG 400, and 0.25% ␤-octyl glucoside. This condition gave rise to two distinct crystal morphologies, rectangular block-shaped (Form A) versus hexagonal-shaped crystals (Form B), which diffracted to 1.3 Å versus 2.4 Å resolution, respectively, using synchrotron radiation. We thus chose to focus on Form A crystals. They belong to space group P2 1 2 1 2 1 with a ϭ 48.45 Å, b ϭ 83.60 Å, and c ϭ 87.53 Å and contain one molecule per asymmetric unit with a V m value of 2.24 Å 3 /Da. Data Collection, Structure Solution, and Refinement-X-ray diffraction data were collected at 100 K at the EMBL beamlines X13 and BW7A in DESY (Hamburg, Germany) using single crystals that were previously cryo-cooled in liquid nitrogen following a brief incubation (ϳ20 s) in a cryoprotectant solution (100 mM Tris, pH 8.0, 2 M ammonium sulfate, 2% PEG 400, 20% (v/v) glycerol). All of the diffraction data were processed using the programs Denzo and Scalepack (32). Data statistics are presented in Table 1. The native SYE1 structure was solved by molecular replacement using maximum-likelihood methods implemented in the program PHASER (33). The search model was generated from the coordinates of morphinone reductase (Protein Data Bank entry 1GWJ) which exhibited a 50% sequence identity with SYE1. Nonconserved residues were replaced by alanine (or glycine), and all of the insertions were deleted, as well as all ligands, including water molecules and the FMN cofactor. The correctness of the structure solution was assessed from the quality of the electron density for missing side chains and other unique structural features. Several rounds of model building were subsequently alternated with further crystallographic refinement employing simulated annealing, conjugate gradient minimization, and individual B-factor refinement using a maximum-likelihood target function, as implemented in CNS (34) (see Table 1). A random 5% set of the possible reduced data were selected and used to flag reflections for the calculation of R free values. The same set of reflections was used for the calculation of the R free for all data sets. The FMN cofactor was initially modeled with a planar isoalloxazine ring, but difference electron density maps after the first round of refinement showed clear evidence for a butterfly-bent ring about the imaginary N-5-N-10 axis. The planarity weights were therefore relaxed during subsequent rounds of refinement to facilitate proper refinement of the FMN ring atoms. Ensuing electron density maps indeed contained no residual positive or negative difference electron density for the FMN substructure. A PEG 400 molecule, a sulfate ion, a ␤-octyl glucoside molecule, and water molecules were modeled in the final stages of refinement. The electron density maps allowed us to model the entire PEG 400 molecule, which is consistent with the refined temperature factors of the PEG 400 molecule (see Table 1). All of the ligands were refined using atomic parameters extracted from the HICUP server. Ligand p-hydroxyacetophenone was composed by combining two molecules (phenol and an acetyl group). For all three soaking experiments difference density maps were initially calculated using coefficients F obs,soak Ϫ F obs,native and calculated phases from the high resolution native structure. All of the residues were found in the most favored, additional, or generously allowed regions of the Ramachandran plot (program Procheck (35)) except for Trp 274 , an aromatic amino acid that lines the active site tunnel. The r.m.s.d. values for the different models were calculated by the program SEQUOIA (36), and the figures were drawn with Pymol (www.pymol.org).
Soaking Experiments-Form A crystals of SYE1 were washed two times with 100 mM Tris, pH 8.0, 2 M ammonium sulfate and were subsequently incubated in 100 mM Tris, pH 8.0, 2 M ammonium sulfate, supplemented with 50 mM of the desired ligand. Complex formation in the crystals was monitored via the rapid change of color from yellow to lime-green, in the case of the phenolic ligands and from yellow to colorless in the case of NADH. Incubation times for the phenolic ligands were typically 2-5 min depending on the size of the crystal used and in the case of NADH 30 min to several hours. After incubation, the crystals were flash frozen in liquid nitrogen.
Kinetic Studies-All of the kinetic parameters were determined under strict anaerobic conditions, and all of the assays were performed in triplicate. The reaction mixtures consisting of 0.5ϫ phosphate-buffered saline buffer, pH 7.3 (5.75 g of Na 2 HPO 4 , 1.49 g of NaH 2 PO 4 , and 2.92 g of NaCl/liter), 60 -120 M NADH, and varying concentrations of substrate were prepared within a glove box (Coy Laboratories, Grass Lake, MI) and transferred to sealable quartz cuvettes (Hellma). The concentration of NADH was determined spectrophotometrically using a molar extinction coefficient of 6220 M Ϫ1 cm Ϫ1 at 340 nm. Enzyme solutions were prepared in sealed flasks and made anaerobic by bubbling with N 2 gas. The reactions were initiated by the addition of enzyme to the reaction mixtures using a Hamilton needle. The enzyme concentration was kept constant at 100 nM. SYE activity was assayed by following the oxidation of NADH using a double beam Uvikon 943 spectrophotometer. Apparent steady state kinetic constants were determined by least squares fitting procedures to the standard Michaelis-Menten equation.
Protein Data Bank Accession Codes-The coordinates and structure factors for the four structures of SYE1 reported here are available via the Protein Data Bank (www.rcsb.org) with accession codes 2GOU, 2GQ8, 2GQ9, and 2GQA.

RESULTS AND DISCUSSION
Overall Structure of SYE1-The structure of oxidized SYE1 was determined by molecular replacement methods using morphinone reductase as the search model and was refined at 1.4 Å resolution (Table 1). SYE1, like all OYE homologues characterized to date (26,27,29,37), folds into an ␣ 8 ␤ 8 barrel containing one noncovalently bound FMN molecule located at the top of is the intensity of the i th measurement of reflection h , and ͗I(h)͘ is the average value over multiple measurements.
the ␤ barrel (Fig. 1, A and B). Like morphinone reductase (37) and PETN reductase (27), SYE1 has four extra barrel elements: (a) the N-terminal ␤-hairpin that closes the bottom of the bar-rel (Fig. 1A, bottom arrow); (b) a capping subdomain formed by secondary structure elements between ␤3 and ␣3; (c) an ␣-helix linking ␤4 and ␣4, which in OPR1 (26) is thought to aid in substrate recognition; and (d) the ␣-helix between ␤8 and ␣8 that contributes to the binding of FMN.
In the course of refinement difference electron density maps showed clear additional density for three nonprotein elements: a PEG 400 molecule, a sulfate ion, and a ␤-octyl glucoside molecule (Fig.  1A). PEG 400 was an essential component of the crystallization condition and could be traced in the active site of SYE1 (Fig. 1C). OPR1 also binds PEG in the active site, but in that case the PEG molecule is doubly bent and resembles the substrate 12-oxophytodienoic acid. In SYE1 the PEG 400 molecule is extended, and only the first six atoms (corresponding to two ethylene glycol units) enter the active site tunnel, whereas the rest of the molecule bends around the hydrophobic active site entrance and extends outwards like the tail of a pig (Fig.  1C). Interactions between the PEG molecule and SYE1 are largely based on van der Waals' interactions with the hydrophobic lining of the active site. The sulfate ion and the ␤-octyl glucoside molecule participate in intermolecular contacts between two adjacent SYE1 molecules in the crystal lattice. The ␤-octyl glucoside molecule lies on one side parallel to His 39 and is coordinated on the other side by Lys 261 from a symmetry-related molecule. In contrast to the sugar ring, the tail of the ␤-octyl glucoside molecule is not well defined in the electron density maps. The sulfate ion is coordinated by His 337 and Arg 336 and by His 344 , Arg 301 , and Arg 323 from a symmetry-related SYE1. The last two amino acids are involved in the coordination of the flavin phosphate oxygens.
The  r.m.s.d. over 293 C␣ atoms), and YqjM (1.7 Å r.m.s.d. over 217 C␣ atoms). The largest differences can be found in the capping subdomain and the loop between ␤6 and ␣6. The ␤-hairpin that connects the two ␤-strands of the capping subdomain is orientated in SYE1 toward the inside of the barrel.
The Active Site Structure-The active site of SYE1 is composed of a narrow and shallow tunnel about 10 -15 Å long and is lined by aromatic and hydrophobic residues. At its bottom, the tunnel opens up toward the flavin. Interestingly, Loop 6 4 folds toward the center of the barrel and contributes to a significant part of the lining of the tunnel. This makes the tunnel markedly narrower than any other OYE characterized thus far, which is surprising given the high overall structural similarity with the bacterial homologues PETN reductase and morphinone reductase. We propose that this feature may be of functional importance.
FMN is tightly bound in the protein and exhibits the lowest B-factors in the structure. The amino acids that contact the riboflavin moiety of FMN are conserved among the different known OYE structures (26,27,29,37): the flavin N-5 interacts with the main chain amide of Thr 26 , O-4 is hydrogen-bonded to the O␥ of Thr 26 and the main chain amide of Gly 58 , and N-3 and O-2 interact with the side chain of Asn 100 . O-2 makes an additional hydrogen bond with the amide of His 181 and the guanidino group of Arg 233 ( Fig. 2A).
Differences can be seen in the amino acids binding the ribityl hydroxyl groups and the phosphate oxygen atoms. The first two hydroxyl groups are each hydrogen-bonded to the guanidino group of Arg 233 , and to the main chain carbonyl oxygen of Pro 24 and a water molecule, respectively. The third OH group and the phosphate oxygens are involved in an extensive hydrogen bonding network involving the side chains of Arg 323 and Arg 301 , the main chain amides of these two arginines and of Gly 300 , and three water molecules. In contrast to the si-face of the flavin, the ribityl side chain, and the dimethylbenzene ring, which are all solvent-accessible, the re-face of the flavin is buried in the structure and is surrounded by residues from Loop 1.
The tricyclic ring of the FMN in SYE1 is strongly butterflybent around the hypothetical N-5-N-10 axis in both the reduced and oxidized states (Fig. 2B). Structural analysis of NADH-reduced SYE1 shows that reduction of the enzyme does not result in a further bending of the flavin (see below). A similar feature is reported for the flavin group in Bacillus YqjM (29), yet only 10% of all flavoproteins are thought to have a bent FMN group (38). Both SYE1 and YqjM have a methionine at position 25 just like many other structurally uncharacterized OYE homologues from the ␤-proteobacteria and the ␥-proteobacteria we have identified via a BLAST search using SYE1 4 Loops are numbered after the ␤-strand they follow. Thus Loop 6 lies between ␤-strand 6 and ␣-helix 6. as a probe (Fig. 2C). Interestingly, all of the other structurally characterized OYEs to date have a leucine instead of methionine at position 25, and they all exhibit planar FMNs in the oxidized state. It has been proposed that a bent flavin would destabilize the formation of semiquinones and thus one-electron reactions of the flavin (39). It is interesting to note that semiquinone species were observed during the photoreduction of OPR1 (13) and OYE (40), but not of YqjM (41) and SYE1 (31). In this respect, the behavior of YqjM and SYE1 is similar to other bacterial OYE homologues such as PETN reductase (42), morphinone reductase (43), XenA reductase, and XenB reductase (10), all of which have planar flavin groups. Clearly then, a mechanism other than the bending of the flavin is responsible for the lack of semiquinone formation in these enzymes. Although the importance of a butterfly-bent FMN is not yet known, our analysis suggests that a methionine residue at position 25 may be a general prerequisite for a butterfly-bent flavin group across all OYEs.
Structures of SYE1 in Complex with p-HBA and p-ACE-Our choice of phenolic ligands was limited by the fact that previous experiments had shown that SYE1 is not capable of binding phenolic ligands with high pK a values (pK a Ͼ 9), such as p-cresol (31). We therefore decided to pursue complexes with p-hydroxybenzaldehyde (p-HBA) and p-hydroxyacetophenone (p-ACE). Upon incubation with the phenolic compounds, the crystals turned from yellow to lime-green. Structures of the complexes were determined by difference Fourier analysis and were refined with full occupancy of the ligands in the active site to 1.7 Å resolution ( Table 1).
The two phenolic compounds were well defined in the difference electron density maps, and the orientation of the molecules could be readily determined in the active site. p-HBA and p-ACE are orientated in the same fashion (Fig. 3A). The phenolic oxygen is placed above the C-2 atom of the flavin, forming a hydrogen bond with His 181 and Asn 184 . Although the aldehyde carbonyl of p-hydroxybenzaldehyde in OYE (25) is hydrogen-bonded to Tyr 375 , no bond of this kind is possible because this residue is replaced by Phe 350 in SYE1 (Fig. 3B). Instead, the aldehyde oxygen is hydrogenbonded to Trp 274 from Loop 6, which is forced to adopt an otherwise unfavorable conformation to serve this role. This novel interaction compensates for the loss of Tyr 375 (OYE numbering) because of a mutation to phenylalanine in SYE1. Morphinone reductase has a Trp 279 corresponding to Trp 274 , but it is located toward the outside of the barrel. Two evolutionary scenarios may explain the role of Loop 6; the restructuring of Loop 6 in SYE1 arose when the enzyme had to substitute for the lost Tyr 375 interaction, or Loop 6 appeared first in the active site. Because Trp 274 may substitute effectively for the Tyr 375 interaction, a Tyr to Trp mutation was not harmful for catalysis.
It has been shown for morphinone reductase that a Y356F mutation (equivalent of Tyr 375 ) increases the K m for NADH by a factor of ϳ7. For SYE1 the K m for NADH binding is 9.3 M, which is consistent with the values for other OYE homologues ( Table 2). Thus Trp 274 appears to compensate well for the absence of Tyr 375 both in terms of phenolic ligand binding and NADH binding. Interestingly, several OYE homologues exhibit the Trp 274 and Phe 350 substitutions. We identified the consensus sequence (A/S)E(A/V)DWDDA flanking the Trp 274 residue, whereas the Phe 350 is always preceded by a leucine residue   (Fig. 3B). We propose that these sequence fingerprints, combined with the structural features in the active sites of SYE1, indicate that OYE homologues possessing these features may constitute a new subfamily of OYEs in ␤and ␥-proteobacteria. NADH-reduced Structure-Upon incubation of SYE1 crystals with NADH, the yellow color of the crystals fades away, indicating a reduction of the flavin. However, the occupancy of NADH in the active site was too low to allow its modeling in the active site. This appears to be a recurring theme in the OYE superfamily, because attempts to obtain such complexes in morphinone reductase (37), OPR1 (26), and OYE (25) also failed. It has been noted for other homologues that reduced flavin harbors ordered water molecules above the ring that replace bound ions (Cl Ϫ ) (25,27). In native SYE1, there are no such water molecules because of the presence of PEG 400 in the active site, but upon reduction with NADH, PEG 400 no longer occupies the active site, and instead three ordered water molecules move into its place.
The "Open" Conformation of Loop3 Provides Insights into NADH Binding-The NADH-reduced and the p-ACE bound structures differ little from the PEG bound SYE1 structure, the most notable feature being the disappearance of the ␤-octyl glucoside molecule in the NADH-reduced structure, which brings about a 90°f lipping of His 39 around its hinge. The structure of the p-HBA complex, on the other hand, was not straightforward to interpret. Incubation with p-HBA increased the internal disorder of the crystals, resulting in the disappearance of the electron density of many side chains. As in the NADH-reduced structure, the ␤-octyl glucoside molecule was washed away in the soak, whereas the sulfate ion became highly disordered.
At first glance, a p-HBA molecule could be unambiguously modeled in the active site of SYE1 stacked against the flavin ring in an orientation similar to what we observed for p-ACE. Quite unexpectedly, a second p-HBA molecule (p-HBA-B) was found sandwiched between two Phe residues (Phe 132 and Phe 350 ) at the entrance to the active site at the capping subdomain (Fig. 4, A and  B). The occupancy of this second site has been estimated to be 0.6 based on a comparison of the refined B-factors of p-HBA-B at that occupancy, with the surrounding residues of SYE1. To accommodate for this molecule, Phe 132 was flipped by 90°a round its 1, and the hairpin loop of the capping subdomain was twisted by about 20°around Val 133 and 80°around Gly 141 toward the outside of the barrel. This open conformation of the capping subdomain creates a hydrophobic cleft next to the active site entrance defined by Tyr 68 , Trp 70 , Phe 132 , Phe 142 , and Phe 350 (Fig. 4C). Modeling of a p-ACE molecule at the site of the second p-HBA explains why a similar interaction was not observed in the p-ACE complex. In that case, the methyl group of p-ACE would approach only 2.4 Å from Trp 70 and would cause a steric clash. Trp 70 is part of Loop 2, which is structurally conserved in SYE1, morphinone reductase, and PETN reductase.
The fortuitous discovery of a second binding site for p-HBA at the capping subdomain provides the first direct evidence to date for any OYE homologue that Loop 3 of the capping subdomain can indeed participate in ligand binding. This provides important new insights into the elusive NADH complex and into ligand binding in general. It is tempting to propose that the two binding platforms for p-HBA reflect the binding sites for the nicotinamide and adenine moieties of NADH, respectively. However, careful modeling of NADH, with the nicotinamide ring occupying the site of p-HBA-A at the flavin ring, shows that the aromatic platform presented by Phe 132 and Phe 350 can only serve as the binding site for the ribose group preceding the adenine moiety of NADH. Such a stacking interaction is frequently observed in protein carbohydrate interactions. The adenine moiety of NADH would therefore have to extend toward the exit of the active site channel, but in the absence of an obvious nearby binding cleft, it is difficult to speculate whether the adenine will bind to SYE1 at all. Interestingly, we find that all structurally uncharacterized SYE1 homologues we have identified in our sequence alignments exhibiting a tryptophan residue and a phenylalanine or tyrosine at positions 274 and 350, respectively (Fig. 3B), also have a phenylalanine or a tryptophan at position 132 (Fig. 4D). The strong conservation of such sequence features suggests that this group of OYE homologues likely binds NADH and other ligands in a similar fashion and that SYE1 can serve as a structural prototype for this subfamily of OYEs.
Our observation of a second binding site for p-HBA is consistent with the recent proposal that the distal part of the NADPH moiety interacts with Loop 3 in PETN (27). It has also been proposed that the preference for NADPH in PETN reductase is brought about by the presence of two arginine residues (Arg 142 and Arg 130 ) in Loop 3 (27,44). In morphinone reductase (which prefers NADH), these arginines are replaced by a glutamate residue (Glu 135 ), whereas in SYE1 we observe a lysine (Lys 130 ) and a

Kinetic constants for the binding of t-2-ketanals to SYE1
Apparent kinetic constants were derived from steady state kinetic analyses using the software Graphpad Prism4. Assays were performed under strict anaerobic conditions in 1 ml of 0.5ϫ phosphate-buffered saline buffer, pH 7.3, and 60 -100 M NADH. Enzyme concentration was 100 nM. Also given, is the structure of each t-2-ketanal and of PEG 400 for comparison purposes. glutamate (Glu 139 ). Because SYE1 exhibits a preference for NADH, this theory cannot be extended to all OYE homologues. When we compare the utilization of reducing power among OYE homologues, we find that the ones preferring NADPH accepted both NADPH and NADH, but with a preference for the former (18,31,37,41,45). On the other hand, OYEs that use NADH cannot accept NADPH. We therefore conclude that one should not search for a mechanism that discriminates between NADH and NADPH, but for one that excludes NADPH.
Further Insights into Substrate Selection and Binding-Thus far, activity of SYE1 has only been shown with relatively small molecules: N-ethylmaleimide, acrolein, and 2-cyclohexenone (31), whereas bulkier molecules, such as nitroglycerine and androstadiene/progesterone, exhibit no activity with SYE1. The shallow and narrow landscape of the active site of SYE1 inevitably imposes size restraints on candidate substrates. The identification of a PEG 400 molecule in the active site tunnel of SYE1 prompted us to explore whether SYE1 accepted long chain unbranched ␣,␤-unsaturated aldehydes as a substrate. Indeed, a steady state kinetic analysis showed that SYE1 showed reactivity toward trans-2-hexenal, trans-2-decenal, and trans-2-dodecenal, three naturally occurring components of plant oils with bactericidal properties against Salmonella (46,47) (Table 3). The k cat /K m values of SYE1 for these substrates are rather similar, differing only by a factor of two, in contrast to tomato OPR1, for which the k cat /K m values for t-2-hexenal and t-2-dodecenal differ by a factor of ϳ20 (13). This is probably due to the much narrower active site of SYE1, which likely can only accommodate substrates with short chains. OPR1, on the other hand, binds a complete PEG 400 molecule in its active site (26) and can thus accommodate larger substrates. Furthermore, SYE1 only binds the first six atoms of PEG 400 in the active site tunnel with the rest of the molecule extending outwards. A preference can be seen for the shorter chain of t-2-hexenal. A similar observation was made for OPR1 (13), which has the highest catalytic efficiency with the small molecule t-2-hexenal.
Considering the fact that PEG 400 was absolutely necessary for protein crystallization and that it binds in the active site, one can argue for a physiological substrate similar to t-2-alkenal. trans-2-Alkenals are breakdown products in lipid peroxidation reactions (46), produced from polyunsaturated fatty acids. Substances such as t-2-hexenal and t-2-dodecenal are volatile compounds of plant oils, produced upon bacterial injury to the plant. S. oneidensis, however, is not a plant pathogen, and polyunsaturated fatty acids are only present in psychrophilic bacteria (48,49). An involvement of SYE1 in detoxification of these compounds is thus unlikely. Nonetheless, the identification of a second ligand-binding site at a hydrophobic cleft near the mouth of the active site entrance, created by a restructuring of Loop 3 in the capping subdomain, indicates that the physiological substrate of SYE1 may be a molecule consisting of a bulky (possibly aromatic) head connected to a thin tail. The "tail" of the substrate could then intercalate to the active site, whereas the bulky head could be stabilized at the entrance of the active site.


extfieldsofstudy:
Medicine
Physics
Computer Science

text:
Mechanochemical Coupling in the Myosin Motor Domain. I. Insights from Equilibrium Active-Site Simulations

Although the major structural transitions in molecular motors are often argued to couple to the binding of Adenosine triphosphate (ATP), the recovery stroke in the conventional myosin has been shown to be dependent on the hydrolysis of ATP. To obtain a clearer mechanistic picture for such “mechanochemical coupling” in myosin, equilibrium active-site simulations with explicit solvent have been carried out to probe the behavior of the motor domain as functions of the nucleotide chemical state and conformation of the converter/relay helix. In conjunction with previous studies of ATP hydrolysis with different active-site conformations and normal mode analysis of structural flexibility, the results help establish an energetics-based framework for understanding the mechanochemical coupling. It is proposed that the activation of hydrolysis does not require the rotation of the lever arm per se, but the two processes are tightly coordinated because both strongly couple to the open/close transition of the active site. The underlying picture involves shifts in the dominant population of different structural motifs as a consequence of changes elsewhere in the motor domain. The contribution of this work and the accompanying paper [36] is to propose the actual mechanism behind these “population shifts” and residues that play important roles in the process. It is suggested that structural flexibilities at both the small and large scales inherent to the motor domain make it possible to implement tight couplings between different structural motifs while maintaining small free-energy drops for processes that occur in the detached states, which is likely a feature shared among many molecular motors. The significantly different flexibility of the active site in different X-ray structures with variable level arm orientations supports the notation that external force sensed by the lever arm may transmit into the active site and influence the chemical steps (nucleotide hydrolysis and/or binding).

Introduction
Coupling between chemical events and conformational transitions is crucial for the function of many biomolecules [1]. Notable examples include molecular motors [2][3][4], signaling proteins [5], and membrane transporters [6], in which large-scale conformational transitions are often coupled to the binding and hydrolysis of Adenosine triphospate (ATP) or phosphorylation. Despite many previous efforts from both experimental and theoretical studies [2][3][4]7,8], the complex nature of the phenomena still makes it an outstanding challenge to establish a concrete understanding of ''mechanochemical coupling'' in those ''molecular machines'' with structural and energetic details.
One particularly interesting question in this context is whether the ATP hydrolysis reaction itself is capable of driving any major structural transitions. Although the answer appears to be positive according to most biochemistry textbooks [1], the popular modern view in the biophysical community [2,9,10] is that the large-scale structural transitions implicated in most ''molecular machines'' are directly coupled to the binding of ATP, and the hydrolysis of ATP mainly functions as a switch to allow the system to move forward into the next functional cycle. Indeed, considering the highly charged nature of ATP (or ATPÁMg 2þ ), it is physically reasonable to expect the binding process to involve substantial local structural adjustments that can then propagate into larger structural transitions. The hydrolysis of ATP, by contrast, involves only subtle structural changes and therefore may not be particularly capable of driving major structural transitions.
Along this line, conventional myosin (which will be referred to as myosin-II or myosin in the subsequent discussions) is interesting to study because it is a prototypical molecular motor [11] that utilizes ATP; more importantly, there is solid experimental evidence [12] that supports the explicit role of ATP hydrolysis in driving the key structural transition (see below).
The fundamental functional cycle of myosin-II has been well-established via kinetic studies and is summarized in Figure 1 [13][14][15][16][17]. The ATP hydrolysis step implicates the detached kinetic states, for which relevant high-resolution Xray structures are available for the truncated motor domain. As shown in Figure 2A, two representative X-ray structures were solved with ATP (PDB code 1FMW [18]) and ADPÁVO 4 À (PDB code 1VOM [19]) as the ligands, respectively, and they were postulated to reflect the conformations before and after (or during) ATP hydrolysis, correspondingly. Fitted into the kinetic scheme of the functional cycle (Figure 1), they are referred to as the ''post-rigor'' and ''pre-powerstroke'' states, respectively. The most visible structural change occurs in the C-terminal converter region, which undergoes a major rotation between the two structures and is believed to propagate into the lever-arm rotation; with the two structures best-fitted based on the first 650 residues, the backbone atom root mean square deviation (RMSD) in the converter is as large as 18.7 Å while it is only 2.4 Å for the rest. In addition to this striking rearrangement, there are also more subtle structural changes in the nucleotide binding site ( Figure 2B) and the ''relay helix'' ( Figure 2C) that connects the binding site and converter.
There has been a large body of experimental [12,[17][18][19][20][21][22][23][24][25][26][27][28] and computational [29][30][31][32] studies aimed at understanding the mechanism of the ''recovery stroke'' between the postrigor and pre-powerstroke states and its connection to the hydrolysis of ATP. The FRET study of Shih et al. [12] has firmly illustrated that the rotation of the converter, which is the major part of the recovery stroke, is dependent on and tightly coupled to the hydrolysis of ATP. This is consistent with the functional cycle of myosin ( Figure 1): (i) the binding of ATP is used to detach the motor domain from the actin and therefore the recovery stroke needs an alternative driving force; (ii) without a tight coupling between the hydrolysis and recovery stroke, the motor may rebind to the actin with a lever-arm configuration that is incapable of performing the powerstroke, which leads to futile ATP hydrolysis.  [13][14][15][16][17] Only a myosin monomer is shown for simplicity. The binding of ATP to the actin-myosin complex (the ''rigor state'') leads to rapid dissociation of myosin from actin without immediate hydrolysis of ATP. Coupled with a major structural change in the orientation of the lever arm (recovery stroke), ATP hydrolysis proceeds, after which the motor domain rebinds to actin weakly. Following the release of P i , the motor domain undergoes ''powerstroke'' during which the orientation of the lever arm changes back to that in the ''rigor state'' and the motor domain becomes strongly bound to actin. Dissociation of ADP leads the system back to the ''rigor state.'' Note that the sum of free-energy drops in an entire cycle is equal to the hydrolysis free energy of ATP in solution, which is the ultimate thermodynamic driving force of the system. doi: 10

Author Summary
The hydrolysis of Adenosine triphosphate (ATP) provides the energy for most life processes, including the motility of molecular motors. How the chemical energy of hydrolysis is converted into mechanical work in these fascinating ''nanomachines'' is a central question that has only been answered in an outline form for almost all molecular motors. The fundamental challenge is that the working cycle of molecular motors involves processes of different physicochemical natures and scales, including ATP chemistry and protein structural transitions of diverse magnitudes (from a few Angstroms to a few nanometers), which makes mechanistic analysis using experiments alone difficult. Combined with previous computational studies from this lab, molecular dynamics simulations help identify energetic and structural properties of myosin, a prototypical molecular motor, that are essential to its energy conversion function. In addition to the role of flexibilities at the domain scale, which has been emphasized in previous studies of similar systems, the current results highlight the comparable significance of local flexibilities and how these flexibilities at different scales are modulated by ATP chemistry.
An interesting point further arises when one considers the energetics of the functional cycle (see Discussion). Since ATP hydrolysis occurs in the detached kinetic states, which are not capable of exerting force on actin and performing work, the reaction is expected to have a small free-energy drop [33], which is consistent with experimentally measured equilibrium constants [14,28]. Naively, the nearly thermoneutral nature of the hydrolysis reaction seems inconsistent with its observed role [12] in driving the major structural transitions in the recovery stroke.
Therefore, a question of critical functional importance is: what are the properties of the myosin motor domain that ensure a tight coupling between ATP hydrolysis and the converter rotation, despite not only the long distance separating the nucleotide binding site and the converter but also the nearly thermoneutral nature of the hydrolysis reaction? To answer this question, it is important to address not only the structural features of the active site but also the energetics of various processes involved in the recovery stroke. In this work, we meet this challenge with a combination of equilibrium simulations of the myosin motor domain. Specifically, potential of mean force (PMF) simulations [34,35] are used to quantitatively characterize the open/closed transition of the nucleotide binding site and analyze how this transition is coupled, thermodynamically, to the structure of the relay helix and orientation of the converter domain. Equilibrium simulations, using both molecular mechanics (MM) and combined quantumn mechanics/molecular mechanics (QM/MM), are also carried out for different nucleotide states (ATP versus ADPÁP i ) to probe the effect of ATP hydrolysis on the active-site behavior, which The difference between the post-rigor (1FMW [18]) and pre-powerstroke (1VOM [19]) states. The structures are aligned based on backbone atoms in the first 650 residues; red and blue correspond to 1FMW, green and yellow correspond to 1VOM. Trp501 is shown in the van der Waals form. (B) Superposition of the active-site region in 1FMW and 1VOM (same color coding as in (A); the nucleotide is shown in the van der Waals form). (C) Superposition of the relay helix and relay loop in 1FMW and 1VOM with the same color coding as in (A). doi:10.1371/journal.pcbi.0030021.g002 aims to explore how ATP hydrolysis might induce the subsequent rebinding of the motor domain to actin (Figure 1).
In the accompanying article [36], several alternative computational approaches, which include targeted molecular dynamics (MD) [37,38], normal mode [39,40] hinge analysis, and statistical coupling analysis [41,42] are used to explore transition pathways for the recovery stroke and to identify residues important to the coupling between distant functional sites.
Compared with previous simulation studies on the mechanochemical coupling in myosin [29][30][31][32]43] and other molecular motors [8,[44][45][46], the present set of studies represents the first serious attempt to examine energetic properties concerning the conformational transitions in these complex systems and a systematic analysis of residues critical to mechanochemical coupling. With all the results taken together, a clearer mechanistic picture emerges for the coupling between ATP hydrolysis and various essential structural motifs, which includes some features that might be shared by different classes of motors.

Results
To facilitate discussion, we summarize the nomenclature and conformational characteristics for key structural motifs in Table 1 following the discussions in previous structural studies [17,47]. Specifically, we refer the active site with the critical salt bridge (Arg238-Glu459) formed as the closed configuration; otherwise as the open configuration. In the 1FMW structure, the converter domain is in a down position with the relay helix straight, while in the 1VOM structure, the converter domain is in an up position with a kink in the relay helix (see Table 1 and Figure 2A).

Energetics for the Open/Closed Transition of the Nucleotide Binding Site in Different Conformational States
Although the converter domain and the active site were seen to adopt different conformations in the two X-ray structures (1FMW and 1VOM) [18,19], to what degree these distant regions are coupled thermodynamically is not well-established; whether the nucleotide chemical state (e.g., ATP versus ADPÁP i ) contributes directly to the coupling is also not clear. These issues are probed with PMF calculations for the open/closed transition in the active site with different X-ray structures of the motor domain, which mainly differ in the lever-arm orientations (down versus up) and relay helix conformations.

Thermodynamics Based on PMF Calculations
The PMFs for the open/closed transition of the active site have been calculated in both one dimension (DRMSD) and two dimensions (RMSD FMW , RMSD VOM ); see Methods. The general trend is very similar; thus, we will only refer to the 1-D PMFs in the discussion, although both sets are shown in Figure 3 for comparison.
With the lever arm down (1FMW:ATP), the PMF has the double-well feature, with the open configuration slightly preferred (;À0.6 kcal/mol) and a low barrier of ;5 kcal/mol ( Figure 3A and 3B). Experimentally it was found that Switch II closure is a rapid (;1000 s À1 ) event with an equilibrium constant close to 1 (DG 0 ; 0 kcal/mol) [14,28]. The drastic difference between the measured open/closed rate (;1000 s À1 ) and the calculated barrier (;5 kcal/mol, which corresponds to a rate of 1.4 3 10 9 s À1 using a pre-factor of k B T/h in a transition state rate expression, ignoring the effect of damping) might be due to the approximate nature of the reaction coordinate used in the PMF calculations; an alternative explanation is that the reporting Trp501 fluorescence might be more directly related to the position of the lever arm rather than the open/closed of the Switch II motif [48] (also see [36]). Nevertheless, the trend that the open and closed configurations of the active site have similar energetics in post-rigor is found in both experiments [27] and the present PMF calculations.
When the lever arm is up (1VOM:ATP, 1VOM:ADPÁP i ), the free-energy well corresponding to the open active site (DRMSD ; À1.5 Å ) disappears and the closed state becomes the much preferred one ( Figure 3C-3F). The PMF profile is qualitatively the same with ATP and ADPÁP i bound, although the open configuration (DRMSD ; À1.5 Å ) appears to be even  higher in free energy in the latter case; this is in qualitative agreement with the expectation based on QM/MM studies of ATP hydrolysis in different active-site conformations that ADPÁP i favors the closed configuration compared with ATP [29,49]. Due to the high energy cost, the open configuration of the active site is essentially inaccessible when the lever arm is up, regardless of the nucleotide chemical state. In other words, the lever-arm orientation is clearly tightly coupled to the active-site conformation, presumably through residues in the relay helix (see below).

Behavior of Individual Residues
It is worth noting that the free-energy minimum corresponding to the closed active site has an RMSD ; 1.5 Å ( Figure 3G and 3H) relative to the 1VOM X-ray structure in all three sets of simulations shown in Figure 3, which suggests that the DRMSD reaction coordinate offers an unbiased description for the open/closed transition in different X-ray structures. The effectiveness of the simulation protocol is also demonstrated by the water structure in the active site. As emphasized in previous simulations and experimental studies [22,29], there are two water molecules in the closed active site and they form stable hydrogen bonds with a number of residues and the c-phosphate of ATP; this hydrogen-bonding network, which does not exist in the open configuration of the active site, helps to orient the lytic water for a favorable nucleophilic attack. Encouragingly, the correct hydrogenbonding network is formed spontaneously in 1FMW:ATP simulations when the active site is closed with umbrella sampling simulations based on DRMSD ( Figure 3G).
Although the active site and the converter are separated by more than 40 Å , the different active-site flexibilities reflected by the open/closed PMFs indicate that variation in the leverarm orientation produces changes that propagate to the active site, which in turn affects the relative stability of the open and closed configurations of the active site. To identify residues that make major contributions to the open/closed PMF and their dependence on the lever-arm orientation, free-energy maps for local torsions of active-site residues are obtained using Equations 5 and 6. The (v 2 , v 3 ) conformational free-energy profile for Arg238 and (u, w) conformational free-energy profiles for residues Asp454, Ile455, Ser456, Gly457, Phe458, and Glu459 are shown in Figure 4. As references, the dihedral angles describing the conformation of those residues in the crystal structures are summarized in Table 2.
The main chain conformational spaces accessed in the PMF simulations are within allowed regions in the Ramachandran plot [50,51], which reassures us that the sampled conformations do not involve unphysical ones. The most striking difference between the 1FMW:ATP and 1VOM:ATP simulations is the Ramachandran plot for the side-chain conformation of Arg238, which involves flipping around v 3 during the open/closed transition. When the lever arm is down, the side chain can adopt two different conformations about (À808, 908) (A) and (À508, À608) (B). Region A is preferred over region B with a free-energy difference of ;3 kcal/mol. When the lever arm is up (1VOM), by contrast, region A is no longer accessible. All other residues within Switch II involve a population shift between different conformations, most notably Gly457 and Glu459. This is consistent with the proposal that Gly457 is involved in a rotation to allow the formation of the Arg238-Glu459 salt bridge and to make an amide hydrogen bond to the c-P of ATP during the Switch II closure [15]; a mutation of Gly457 to even an Ala affects the flexibility in this region and therefore reduces the stabilization of the c-P, which in turn causes the loss of the basal ATPase activity [52].
Although only Switch I/II residues are subject to bias in the umbrella sampling simulations, substantial changes in other regions of the motor domain are observed. As shown in Figure 5, opening the active site in the 1VOM:ATP simulations induces significant changes in the actin-binding interface (regions C,D), the loop between Switch II, and the N-terminal of the relay helix (region B), as well as the Cterminal of the second central ß-sheet (region A); the displacements in those regions are substantially larger than their equilibrium fluctuations (unpublished data). As illustrated by the Ramachandran (u, w) conformational freeenergy profiles ( Figure 6) for the residues between Switch II and the N-terminal of relay helix (region B as defined in Figure 5), the accessible conformational space is more restricted in the lever-up (1VOM:ATP) simulations compared with the lever-down case (1FMW:ATP), which partly explains why the closed state is preferentially stabilized when the lever arm is up. For example, the (u, w) of Ile460 in the X-ray structures of 1FMW and 1VOM are (À1308, 858) and (À1338, 1278), respectively. With the lever arm down, Ile460 can have a widely accessible conformation space surrounding its conformations in both X-ray structures; with the lever arm up, Ile460 preferentially adopts its conformation in the closed state.

Structural Flexibility of the Closed Active Site Before and After ATP Hydrolysis
To further explore how ATP hydrolysis initiates structural changes that ultimately induce the rebinding of the motor domain to actin, substantially longer equilibrium simulations (up to 12 ns) of the closed active site are performed for the ATP and ADPÁP i chemical states with both MM and QM/MM potentials. The results of MM and QM/MM simulations are largely consistent; thus, only MM values are discussed below, although both sets of results are shown in Figure 7 for comparison. We focus on the interaction between Mg 2þnucleotide and the Switch I region because it has been speculated that changes in Switch I are coupled to transitions in actin-binding motifs [17].
In the ATP state ( Figure 7A), the two distances between Mg 2þ and the nucleotide oxygen atoms are significantly shorter than those between Mg 2þ and Thr186/Ser237. The average distances for O 1ß -Mg 2þ and O 1c -Mg 2þ are 1.83 6 0.04 Å and 1.78 6 0.03 Å , respectively; the corresponding values for O Ser237 -Mg 2þ and O Thr186 -Mg 2þ are 2.14 6 0.11 Å and 2.09 6 0.08 Å , respectively. This trend is expected considering that the oxygen atoms in the nucleotide are more negatively charged (;À0.85e) than those in polar amino acid side chains (;À0.66e) [53].
In the ADPÁP i simulations ( Figure 7B), the most striking change is observed for the O Ser237 -Mg 2þ distance, for which the average is 2.48 Å , 0.34 Å longer than that in the ATP simulation. The fluctuation in the distance also increases more than twice that in the ATP state, to 0.25 Å . The largest distance sampled in ;12 ns reaches almost 4 Å , as compared with 2.75 Å in the ATP simulations. In the QM/MM simulations ( Figure 7C and 7D), the Mg 2þ -O distances are consistently shorter, but the trend of a weaker interaction between Mg 2þ and Ser237 in the ADPÁP i simulation is also evident. Therefore, both MM and QM/MM simulations show that following ATP hydrolysis the interaction between the Mg 2þ -nucleotide and Switch I becomes substantially weak-ened; the fact that MM simulations produce similar results as QM/MM simulations suggests that the cause is largely electrostatic in nature and due to the difference in the charge distribution in the nucleotide before and after hydrolysis. A careful examination of the structures in the simulations suggests that the increase in the distance has   contributions from both the structural relaxation of Mg 2þnucleotide following hydrolysis and the slight displacement of the Ser237 side chain ( Figure 7E); the average backbone configuration of the Switch I motif, however, is not perturbed significantly during the nanosecond scale simulations (unpublished data). Nevertheless, the weakened interaction between Mg 2þ -nucleotide and Ser237, as reflected by the larger fluctuations of Ser237 ( Figure 7F), is expected to make structural transition of Switch I more facile. Analysis of sidechain internal fluctuations also indicates that Ser237 is the only residue that is affected significantly by hydrolysis (see Protocol S1).

Defining Mechanistic Questions Regarding Mechanochemical Coupling Based on Functional Constraints
Before discussing the results from the current work in detail, it is useful to recapitulate the fundamental questions regarding mechanochemical coupling in myosin, and other molecular motors, that ultimately drive our study. As pointed out elegantly by Hill and others in their pioneering analysis of coupled biological processes [33,54,55], biological systems need to meet specific thermodynamic and kinetic constraints to acquire the corresponding characteristics necessitated by their biological functions. Specifically for myosin, the constraints based on functional considerations are [55]: (i) processes that occur in the detached states have small freeenergy drops; (ii) in particular, the ATP hydrolysis reaction itself is close to being thermoneutral; (iii) the ATP hydrolysis step is tightly coupled to the reorientation of the lever arm. The first two constraints arise because for a motor to achieve a high thermodynamic efficiency, only the steps capable of performing useful work should be associated with large freeenergy drops, since no work can be done in the detached states (no force can be exerted on the actin), and, considering that the sum of free-energy changes through a functional cycle is strictly a constant (the hydrolysis free energy of ATP in solution), the detached states should have small free-energy drops. In principle, large free-energy increases and drops associated with the detached states may cancel out to still leave a large free-energy drop for the powerstroke, but the large free-energy increases for certain step(s) may compromise the speed of the motor; in other words, meeting the first two constraints is likely correlated with a fast motor without any significant kinetic bottleneck. In this context, we emphasize that one has to distinguish the free-energy drop for ATP hydrolysis in the motor domain, which is small in magnitude (as shown experimentally for not only myosin [14] but also a number of other motors such as F 1 -ATPase [56]), from that for the entire functional cycle; the latter is strictly equivalent to the hydrolysis free energy of ATP in solution and forms the ultimate thermodynamic driving force for the vectorial motion of the motor. The third constraint minimizes the amount of futile ATP hydrolysis, because otherwise a loose coupling may lead to rebinding of the motor to actin with a lever-arm orientation incapable of powerstroke.
Although these ''idealized'' constraints are derived by considering a motor with the optimal efficiency, and most biological motors do not need to have perfect thermodynamic efficiency for their biological functions [7], they provide an interesting set of guidelines to investigate the properties of molecular motors in the context of mechanochemical coupling. Indeed, it is of interest to investigate to what degree these constraints are satisfied in realistic systems and, if so, how they are implemented in structural and energetic terms [49,57].

Constraints I and II in the Myosin Motor Domain
For constraint I, although our study here explicitly deals with only a few steps in the functional cycle, a number of relevant observations can be made. First, normal mode analyses of multiple conformational states clearly demonstrate that the structural transitions are highly correlated with the low-frequency modes and, therefore, intrinsic structural flexibility of the motor domain [36,43]. Although no explicit free-energy cost can be inferred from such normal mode analyses, the significant degree of observed correlation strongly suggests that conformational transition among those states is associated with a modest free-energy change. A series of hydrophobic and polar interactions is observed in both minimum energy path [30] and biased MD simulations [36] to facilitate the large structural rearrangements in the relay helix and converter domain and ensures the corresponding low energy cost. Second, for the open/close of the active site, which has been shown to activate ATP hydrolysis [22,29,49], both current PMF calculations and recent fluorescence experiments [27] suggest that the transition is nearly thermoneutral with the lever arm down (the expected orientation following detachment from actin). The PMF calculations are worthwhile because the interpretation of the fluorescence is not straightforward and the signal due to Trp501 might be better correlated with the lever-arm position rather than the active-site configuration per se [48]. Third, for processes other than the powerstroke, although not explicitly studied here, it is natural to envision how they can be nearly thermoneutral. The detachment of the motor domain from actin is compensated for by the binding of ATP; the rebinding of the motor domain to actin following ATP hydrolysis is known to be only a weak association; the weak ! strong actin-binding transition is compensated for with dissociation of hydrolysis product(s) and may partially couple to the powerstroke.
For constraint II, which is consistent with the measured equilibrium constant [14,27], both our [29] and others' [31] QM/MM simulations suggest that the small free-energy drop is likely controlled by well-orchestrated electrostatic environment of the active site. Both favorable and unfavorable contributions from residues/water in the closed active site have been identified, which together result in a small hydrolysis-free energy. Importantly, as argued below, the small free-energy drop for the hydrolysis step itself does not compromise the coupling between ATP hydrolysis and leverarm orientation as far as both couple tightly to the open/close transition of the active site.

Constraint III in the Myosin Motor Domain
As the major focus of this work, understanding how constraint III is implemented is fascinating because of its allosteric nature: ATP hydrolysis in the active site is strongly (E) Overlay of the active site in two snapshots from MM simulations with ATP and ADPÁP i , respectively. The one with color coding is for the ATP state while the one in grey is for the ADPÁP i state. The Mg 2þ -O Ser237 distance is much longer in the ADPÁP i -state snapshot, due to displacements in both the Ser237 sidechain and the hydrolyzed nucleotide. (F) Root mean square fluctuation (RMSF) for Ca atoms in the three critical active-site motifs based on equilibrium MM simulations. Ser237 has notably larger fluctuations in the ADPÁP i state, although the rest residues do not show distinct differences. doi:10.1371/journal.pcbi.0030021.g007 coupled to the rotation of the converter domain (thus, lever arm) at 40 Å away. It is physically reasonable to argue that the two processes are coupled indirectly by both being directly coupled to variations in the active-site configuration. Indeed, both structural studies [19,22] and QM/MM simulations [29] support the idea that ATP hydrolysis is only possible with the closed active site and is prohibited in the open active site. In addition, PMF calculations ( Figure 3) show explicitly that the active site open/close transition is sensitive to the conformation of the lever arm and relay helix. In the post-rigor state (lever arm down, relay helix straight), the open/close transition is a rapid and thermoneutral process, while with the pre-powerstroke state (lever arm up, relay helix kinked) the active site is strongly biased (by more than 15 kcal/mol) toward the closed configuration.
Taken together, at the thermodynamic level, the simulation results show explicitly how ATP hydrolysis is coupled to the lever-arm rotation through a series of population shifts in local conformations (Figure 8). Starting with the post-rigor conformation, the active site (more precisely, Switch II) rapidly explores the open and closed configurations; with the closed configuration, ATP hydrolysis is a relatively rapid and reversible process. Based on the thermodynamic linkage between the hydrolysis and active-site closure (see the O-C-DPÁP i -TP plane in Figure 8), the dominant population of the active site is shifted toward the closed configuration once ATP is hydrolyzed (with ADPÁP i bound). In turn, based on the thermodynamic linkage between the active-site open/close transition and the converter rotation (see the D-U-O-C plane in Figure 8, assuming that the converter and relay helix are tightly coupled), closure of the active site drives the dominant population of the lever arm from the post-rigor to the prepowerstroke orientation.
Although this description is purely thermodynamic in nature and does not specify the sequence of events (see discussions in the accompanying paper [36]), the clear energetic relations make it possible to refine models concerning the functional cycle, which is difficult to do with a structure-only analysis [30]. For example, most kinetic schemes in the literature assume that activating the hydrolysis depends on (at least partial) rotation of the lever arm [12], because it is difficult in experiments to simultaneously monitor the active-site motion and lever-arm rotation due to their significantly different length scales [26,27,48]. Our computational analyses taken together, however, suggest that the activation of ATP hydrolysis does not require the rotation of the lever arm per se, because closure of the active site, which directly regulates ATP hydrolysis, is a facile process even with the lever arm in the post-rigor orientation. This is consistent with the phenotype of many ''uncoupling mutants'' [23,24], in which the ATPase activity is similar to the wild-type but motility (which presumably is correlated with lever-arm rotation) is significantly compromised or abolished.
As mentioned in the Introduction, the tight coupling between ATP hydrolysis and converter rotation is even more striking considering the near-thermoneutral nature of processes in the detached states. With the energetics-based analysis, it is worth noting that a nearly thermoneutral ATP hydrolysis (in the motor) does not imply a weak coupling to lever-arm orientation. Provided that the lever arm is strongly correlated to the active-site open/close transition through the relay helix, the only requirement for ATP hydrolysis (for the sake of avoiding futile hydrolysis) is that it is tightly coupled to the same active-site transition; the magnitude of this coupling depends only on the relative hydrolysis free energy in the open and closed active-site configurations. As a nearthermoneutral reaction in the closed state but unfavorable reaction in the open state, the ATP hydrolysis is strongly coupled to the open/close transition and, simultaneously, satisfies the thermodynamic constraint II imposed on the motor.
Another point clearly underlined by the energetics-based analysis is that the coupling between ATP hydrolysis and the lever-arm rotation relies on the intrinsic flexibility of the motor domain at both the small (active-site) and large (converter-rotation) scales [29]. In this sense, the current description is reminiscent of the ''rectified Brownian diffusion'' framework [58] proposed for a number of motors including myosin. The current analysis, however, is built on an atomistic description of the relevant processes with concrete structural and energetic features; therefore, we are able to provide mechanistic details regarding how ''Brownian motions'' of the motor domain are ''rectified'' by the hydrolysis of ATP, or, in the language of physical chemistry, how the dominant population of different structural motifs is shifted by events that occur elsewhere in the motor domain. Moreover, we note that although many previous analyses of ''molecular machines'' based on normal modes [29,44,[59][60][61][62][63] emphasized the role of intrinsic structural flexibility at the domain scale, the current study highlights the comparable significance of local flexibility (e.g., active site) and, more importantly, how these flexibilities at different scales are modulated by chemistry. Conversely, the significantly different flexibility of the active site, as captured by the PMF calculations, in different X-ray structures of the motor domain, supports the notion that external force sensed by the lever arm may transmit into the active site and perturb processes such as ATP hydrolysis or product release, which ultimately leads to the change in the motor speed [64].

Communication between Active Site and the Actin-Binding Interface
In addition to the coupling between hydrolysis and the converter, an issue of major importance concerns the coupling between the active-site activities and the actin interface. In the LymnÀTaylor scheme [13][14][15], hydrolysis of ATP leads to the reattachment of the motor domain to actin, although the corresponding mechanism is not well-understood, due largely to the lack of a high-resolution structure for the myosinÀactin complex. The current set of equilibrium simulations of the active site, although limited by the nanosecond time scale and the boundary condition, have provided useful hints.
First, significant distortion of the actin-binding interface and the central ß sheet has been observed to be induced by the opening of the active site in the umbrella sampling simulations ( Figure 5), strongly suggesting the coupling between the open/close of the active site and motion in the actin-binding interface. This hypothesis is consistent with the structural analysis of Geeves and Holmes [17], who proposed that during the powerstroke the opening of Switch I induces the twisting of the central ß sheet as well as the closure of the actin-binding cleft.
In addition, an interesting observation from simulations with different nucleotide chemical states found that the interaction between the Mg 2þ and Switch I (in particular, Ser 237) is substantially reduced once ATP is hydrolyzed; this is observed with both MM and QM/MM potential functions, which suggests that the cause is largely electrostatic and structural in nature rather than due to any subtle quantum mechanical effects. This supports the hypothesis of Reubold et al. [65] and Coureux et al. [66], who analyzed the X-ray structures of nucleotide-free myosin II and myosin V and proposed that hydrolysis of ATP leads to the destabilization of the closed configuration of Switch I. This might prepare for the opening of the Switch I loop, which has been postulated to occur upon binding of actin and is necessary for the release of the hydrolysis products [17]. We note that as this work was in progress, a similar observation was reported in a QM/MM minimization study of myosin [31] and an MD simulation using MM potential [32]. Our unique contribution here is the confirmation of this important effect with both MM and QM/MM potential functions with extensive sampling. An interesting observation in [32] is that Asn475 breaks its hydrogen bond with Switch II upon ATP hydrolysis. This, however, was not observed in either the MM or QM/MM simulations here; to what extent this is due to the spherical boundary condition of these simulations is being investigated (also see Protocol S1). Interestingly, we note that Koppole et al. [32] used a mixture of the force-field parameters in their model, which involves polar-hydrogen parameters for non-aromatic residues and all-atom parameters for the aromatic residues. It is not clear what justifies such mixing of forcefield parameters and what is the impact of this on the simulation result.

Conclusions
Vectorial biological systems such as molecular motors and pumps are fascinating to study [57] because their unique functional characteristics may impose specific thermodynamic and kinetic constraints [33,54,55] that distinguish them from regular enzymes. With remarkable progress in structural and dynamical characterizations of complex biomolecules in both experimental and theoretical arenas, the challenge is to establish to what degree these ''idealized'' constraints are satisfied and to reveal the physical and chemical principles behind their implementation and regulation in structural and energetic terms.
Specifically for myosin II, a characteristic molecular motor, we argue that the key constraints the motor domain is likely to satisfy include small free-energy drops for all steps in the detached kinetic states, and that events at the active site are tightly coupled to distant structural motifs, which are consistent with a broad range of experimental observations. Molecular simulations are then carried out to reveal how these constraints, or ''mechanochemical coupling,'' are reflected in the structural, energetic, and dynamical features of the motor domain. Specifically in this work, equilibrium active-site simulations with explicit solvent molecules have been carried out to probe the behavior of the motor domain as functions of the nucleotide chemical state and conformation of the converter/relay helix. In conjunction with our previous QM/MM studies of ATP hydrolysis in different active-site conformations [29] and normal mode analysis of the motor domain [43], the results here help establish an energetics-based framework for understanding the mechanochemical coupling, in which the active-site closure tightly regulates both the hydrolysis activity and orientation of the converter (thus, lever arm). The flexibilities at both the small (active-site open/close) and large (converter-rotation) scales inherent to the motor domain make it possible to implement these couplings while maintaining small free-energy changes for all key changes occurring in the detached states; in this sense, the current framework is reminiscent of the ''rectified Brownian motion'' model of molecule motors [58]. The contribution of this work and the accompanying paper [36] is to propose the actual mechanism for ''rectifying'' these ''Brownian motions'' and residues that play important roles in the process. In addition, the significant coupling between the converter orientation/relay helix and the active-site configuration as identified by the PMF calculations provides a concrete energetics basis for how active-site activities can be affected by the external force sensed by the lever arm. It is anticipated that many motors, not only those in the myosin family, share many similar features in their mechanochemical coupling cycles.
In the near future, the results reported here together with additional simulations of the lever-arm transition will facilitate the development of an effective phenomenological model [67,68] that can simulate the entire functional cycle of myosin(s), which ultimately provide the missing link between microscopic structural and energetic details with the behavior of motors at the macroscopic length and time scales. Before such a model can be developed, however, key questions regarding myosin-actin interactions need to be answered [17,69], which requires contributions from both experimental and theoretical efforts.

Materials and Methods
PMF simulations with explicit solvent. To estimate the PMF for the active-site open/close transition, umbrella sampling calculations with a stochastic boundary condition and explicit solvent molecules are carried out using the CHARMM program (c32a2 version) [70]. Three systems are simulated: post-rigor with ATP bound (denoted as 1FMW:ATP), pre-powerstroke with ATP bound (denoted as 1VOM: ATP), and pre-powerstroke with ADP and P i bound (denoted as 1VOM:ADPÁP i ). As starting coordinates, the X-ray structures solved by Rayment and co-workers for the Dictyostelium discoideum myosin II motor domain are used; the lever arm is down and the active site is open in 1FMW [18] (resolution 2.15 Å ), while the lever arm is up and the active site is closed in 1VOM [19] (resolution 1.90 Å ). The missing residues (1, 203-208, 500-508, 622-626 in 1FMW, and 1, 205-208, 711, 716-719, 724-730 in 1VOM) are modeled by using the program ModLoop [71,72]. The protein atoms and the Mg 2þ ion in the active site are described with the all-atom CHARMM force field for proteins [53], and the water molecules are described with the TIP3P model [73]. The hydrogen atoms are added with HBUILD [74]. In all cases the system is partitioned into a 32-Å inner region centered at the geometric center of ATP, P-loop, Switch I, and Switch II (about 12,800 atoms), with the remaining portion of the system in the outer region. Newtonian equations-of-motion are solved for the MD region (within 28 Å ), and Langevin equations-of-motion are solved for the buffer region (28-32 Å ) with a temperature bath of 300 K [75]. All the atoms in the inner region are subject to a weak GEO type of restraining potential to keep them inside the inner sphere with the MMFP module of CHARMM; the effect of the restraint on most inner-region atoms is negligible. All protein atoms in the buffer region are harmonically restrained with force constants determined directly from the B-factors in the PDB file [75]. Langevin atoms are updated heuristically during the simulation to consistently treat protein groups and water molecules that may switch regions during the simulation. All bonds involving hydrogen are constrained using the SHAKE algorithm [76], with a relative geometric tolerance of 10 À10 , and the time step is set to 1 fs. Nonbonded interactions within the inner sphere are treated with an extended electrostatics model, in which groups beyond 12 Å interact as multipoles [77]. The entire system is heated gradually to 300 K and equilibrated for 160 ps prior to the production simulations.
To account for the electrostatics between the inner and outer region atoms and the effect of solvation, the generalized solvent boundary potential approach developed by Roux and co-workers is used [78,79]. In the generalized solvent boundary potential, the contribution from distant protein charges (screened by the bulk solvent) in the out region DW where q a is the partial charge on atom a. The dielectric effect on the interactions among inner-region atoms is represented through a reaction field term where Q m/n is the generalized multipole moment and M mn is the element of the reaction field matrix, M. The static field due to the outer region / ðioÞ s is evaluated with the linear Poisson-Boltzmann equation using a focusing scheme that place a 70-Å cube of fine grid (0.4 Å ) into a larger 120-Å cube of coarse grid (1.2 Å ). The reaction field matrix, M, is evaluated using 400 spherical harmonics. In the Poisson-Boltzmann calculations, the protein dielectric constant e p ¼ 1, the water dielectric constant e w ¼ 80, and 0.0 M salt concentration is used. The optimized radii of Roux and Nina [80,81], based on experimental solvation energies of small molecules as well as the calculated interaction energy with explicit waters, are adopted to define the solvent-solute dielectric boundary.
In the umbrella sampling simulations, intermediate structures along the transition pathway between the open and closed active site are generated by using a 1-D RMSD restraint on all atoms in the active site (Switch I: 233-238 and Switch II: 454-459). Two harmonic restraining potentials are applied to rmsd (X t À X open ) and rmsd (X t À X closed ), respectively, with the condition rmsd (X t À X open ) þ rmsd (X t À X closed ) ¼ 2.4 Å ; X t is the instantaneous structure, and X open and X closed are the two reference configurations of the active site based on the 1FMW and 1VOM X-ray structures. The intermediates are generated at 0.1-Å RMSD intervals, which produced 25 intermediate structures spanning the range of ;2.4 Å RMSD separating the two reference structures. These intermediate structures are then used as the starting points for umbrella sampling simulations [34,35,82] with the restraining potential w j on DRMSD as a reaction coordinate [83,84] where DD min specifies the position of the minimum for the harmonic umbrella potential in a specific window. The DRMSD is the difference between the all-atom RMSD values for the Switch I and Switch II regions of the instantaneous structure (X t ) from the two reference structures (X open and X closed ), DRMSD ¼ rmsdðX t À X open Þ À rmsdðX t À X closed Þ ð 4Þ The force constant K RMSD is gradually reduced from 500 kcal/(molÁÅ 2 ) to 25 kcal/(molÁÅ 2 ) in the 50-ps equilibration simulation for each window and kept fixed for another 100 ps of production simulation. Additional windows are added or additional simulations are performed for some windows by inspecting the 2-D distribution of RMSD, q bias (rmsd (X t À X open ), rmsd (X t À X closed )). The Weighted Histogram Analysis Method (WHAM) [82,85,86] is used to obtain the PMF along the 1-D reaction coordinate (DRMSD) from the instantaneous values of DRMSD saved every step during the MD simulations. The convergence of the PMF profiles is checked by calculations using substantially longer trajectories, up to 450 ps per window and in total about 11 ns (see Protocol S1). The MD simulations with the umbrella potential applied along DRMSD as the reaction coordinate contain information about many individual degrees of freedom that contribute to the overall transition between the open and closed states. It is possible to estimate the free-energy profile along other additional degrees of freedom using trajectories from these umbrella sampling MD simulations [83,84]. For example, one can obtain the PMF surface as a function of the 2-D reaction coordinates defined by (rmsd (X t À X open ) and (X t À X closed ). The procedure to obtain such PMF begins by considering the biased probability distribution q bias (n 1 , n 2 , n 3 ) where n 1 ¼ rmsd (X t À X open ), n 2 ¼ rmsd (X t À X closed ), and n 3 ¼ DRMSD. The unbiased probability distribution q(n 1 , n 2 , n 3 ) can be obtained using the standard Weighted Histogram Analysis Method approach (by setting the biased potential in n 1 , n 2 , direction to zero). The probability q(n 1 , n 2 ) can then be obtained by integrating q(n 1 , n 2 , n 3 ) along n 3 : qðn 1 ; n 2 Þ ¼ Z qðn 1 ; n 2 ; n 3 Þdn 3 ð5Þ and the corresponding PMF surface can be then obtained from Wðn 1 ; n 2 Þ ¼ Àk B Tlnqðn 1 ; n 2 Þ ð 6Þ In a similar way, the PMF along the local (u, w) space for some key residues are calculated from the trajectories saved every 100 steps during the MD simulations.
Equilibrium MD simulations of the closed state with different nucleotides using MM and QM/MM potentials. To probe the effect of ATP hydrolysis on the active-site structure, especially the Switch I region that has been proposed to communicate with the distant actin-binding interface, we carry out both MM and QM/MM [87,88] MD simulations for the pre-powerstroke state (1VOM [19]) with different nucleotides (ATP or ADPÁP i ). QM/MM simulations are also carried out to investigate if the qualitative behavior is sensitive to the treatment of charge transfer and polarization associated with the nucleotide.
The simulation setup is essentially the same as that in the PMF simulations described above, except that a smaller inner region of 20-Å radius is used to allow much longer simulations of the active site (;12 ns for MM and ;4 ns for QM/MM simulations). In the QM/MM simulations, 48 atoms are described quantum mechanically, including tri-phosphate and part of the ribose group of ATP, the lytic water (or the corresponding atoms in ADPÁP i ), the Mg 2þ ion and all its ligands (the side chains of Thr186, Ser237, and two water molecules), as well as the side chain of the conserved Ser236. These atoms are described at the SCC-DFTB level [89,90] with the recent parameterization for phosphate systems (Yang et al., unpublished data). Four link atoms are introduced to saturate the valence of the QM atoms at the QM/MM interface [88], which interact with all MM atoms electrostatically except with the link host (Ca atoms in Thr186, Ser236, Ser237, and C59 in the ribose); no van der Waals interactions are considered for the link atoms. This link atom treatment has been shown to work well when the nearby MM charges are small [91][92][93], which is the case here. The generalized solvent boundary potential approach is used for treating electrostatics in both MM [78,79] and QM/MM [94] simulations.