Dynamic Regulation of Phenylalanine Hydroxylase by Simulated Redox Manipulation

Recent clinical studies revealed increased phenylalanine levels and phenylalanine to tyrosine ratios in patients suffering from infection, inflammation and general immune activity. These data implicated down-regulation of activity of phenylalanine hydroxylase by oxidative stress upon in vivo immune activation. Though the structural damage of oxidative stress is expected to be comparably small, a structural rationale for this experimental finding was lacking. Hence, we investigated the impact of side chain oxidation at two vicinal cysteine residues on local conformational flexibility in the protein by comparative molecular dynamics simulations. Analysis of backbone dynamics revealed a highly flexible loop region (Tyr138-loop) in proximity to the active center of phenylalanine hydroxylase. We observed elevated loop dynamics in connection with a loop movement towards the active site in the oxidized state, thereby partially blocking access for the substrate phenylalanine. These findings were confirmed by extensive replica exchange molecular dynamics simulations and serve as a first structural explanation for decreased enzyme turnover in situations of oxidative stress.


Introduction
The mononuclear non-heme iron containing monooxygenase phenylalanine hydroxylase (PAH, EC 1.14.16.1) is a member of the aromatic amino acid hydroxylase family as tyrosine hydroxylase or tryptophan hydroxylases 1, 2. It catalyzes the oxidation of phenylalanine (Phe) to tyrosine (Tyr). This oxidation of the substrate amino acid is achieved by molecular oxygen and the reductive co-factor 5,6,7,8-tetrahydrobiopterin (BH 4 ) [1]. Hydroxylation of Phe to Tyr is the committed step in Phe catabolism and therefore requires strict regulation to ensure homeostasis of the essential amino acid Phe. Dysfunction of PAH causes phenylketonuria (PKU) [2,3], a common and well-examined genetic disease [4] leading to mental retardation upon accumulation of Phe [5].
The catalytic mechanism of PAH involves the reduction from ferric (III) to ferrous (II) form by the co-factor BH 4 . Ferrous iron is subsequently oxidized by molecular oxygen to a Fe(IV)O intermediate, which in turn hydroxylates Phe to Tyr [6][7][8][9]. Eukaryotic PAH is found as a homotetramer, where each subunit contains three distinct domains: an N-terminal autoregulatory domain (residues 1-142), the catalytic domain (residues 143-410) and a tetramerization domain stabilizing the quaternary structure (residues 411-452) [10,11]. BH 4 acts as a negative regulator of PAH activity by stabilizing an inactive form [12], whereas the enzyme is activated by phosphorylation of Ser16 [13] as well as binding of Phe [14]. Activity of PAH increases upon preincubation with Phe up to 100-fold [15]. Allosteric effects are held responsible for the activating effect of Phe by conformational changes in the protein [16]. The positive cooperativity of Phe binding allows for fast responses to increased Phe levels in order to avoid damages to the brain [17]. Recently, a potential binding site of Phe distal to the active site was discovered by a crystallographic study of PAH from Chromobacterium violaceum [18].
Expression of the truncated PAH catalytic domain does not result in a loss in activity suggesting the autoregulatory domain to be a key factor in allosteric regulation [19]. The structural background of PAH activation remains still unclear, as crystal structures are not yet available for the full sequence protein. The only available crystal structures containing the autoregulatory domain of a dimeric rat PAH lack parts of the tetramerization domain and structural insights into residues 1-18 including the phosphorylation site at Ser16 [20,21]. Nevertheless, the position of the additional residues of the regulatory domain suggest a complete physical blockage of the active site in the inactive enzyme [22]. However, it remains unclear, if activation is due to a localized movement of the regulatory sequence or a global conformational reorganization in the protein.
Current studies found immune activation and inflammation to be associated with an increase in the ratio of Phe to Tyr in blood of patients suffering from trauma, sepsis, carcinoma and infection [23][24][25] suggesting an inhibition of PAH [26]. Recently, treatment with immune modulating interferon-a of patients with hepatitis C virus, resulted in an increase of Phe and the Phe to Tyr ratio [27] consistent with earlier findings concerning patients with multiple melanoma treated with interferon-a [28]. As immune activation of macrophages is paralleled by the release of toxic reactive oxygen species (ROS) and neopterin [29,30], oxidative stress is discussed as chemical background for PAH deactivation consistent with studies correlating PKU with oxidative stress [31]. This hypothesis is further strengthened, as phenylalanine concentrations in blood samples of cancer patients were shown to correlate to the oxidative stress marker isoprostane-8 [24]. Immune activation marker neopterin was shown to further increase oxidizing capacity of ROS [32] and therefore raise concentrations of protein oxidation products [33] whilst decreasing antioxidant levels [34]. Modulation of PAH activity has also been reported after reaction with various disulfide reagents [35,36], pointing particular attention to cysteine residues in PAH.
To investigate the impact of oxidative stress on PAH structure and dynamics at atomic level, molecular dynamics (MD) simulations of the catalytic domain were performed. Inspection of 23 available X-ray structures of the catalytic domain of PAH in the PDB [37,38] revealed one single apparent oxidative site at two free cysteine residues (Cys203, Cys334) in close proximity to each other (4 Å ) but distant from the active center of the enzyme (15 Å ) (see Figure 1). An intrinsic connection between these cysteine residues and PAH activity in vivo is documented by mutation data from genetic studies of PKU. Point mutations of both residues are associated with mild forms of PKU [39][40][41].
Comparative MD simulations were performed to analyze conformational differences between native PAH (PAH nat ) and PAH in oxidized state (PAH ox ) containing an artificially introduced disulfide bond between these residues. Findings from these simulations were confirmed by follow-up replica exchange molecular dynamics (REMD) simulations increasing sampling of the observed conformational transitions in the protein.

Methods
Two MD simulations were set up in parallel using the aforementioned starting structures PAH nat and PAH ox . In simulation PAH nat of native PAH we treated Cys203 and Cys334 as individual cysteine residues, as reported in all available X-ray structures. In simulation of PAH ox these cysteine residues were connected via an introduced disulfide bond to represent an oxidized state of PAH. This introduced disulfide bond between Cys203 and Cys334 is the only difference between the two simulated states PAH nat and PAH ox .
As a starting point for simulations we chose the X-ray structure of human PAH with bound co-factor and highest available resolution (PDB: 1J8U [42], resolution:1.50 Å ), although the bound co-factor and the iron atom were both found in reduced form. All atom root mean square deviation (RMSD) values to high resolution X-ray structures containing the oxidized species (PDB: 1TG2 [43], 1LTZ [44]) were found to be reasonably low (RMSD ,1.13 Å ) to set up simulations using template structure 1J8U with both oxidized co-factor 7,8-dihydrobiopterin BH 2 and Fe(III).
Both all-atom MD simulations of the PAH catalytic domain containing a Fe(III)-ion and the co-factor BH 2 were carried out using the AMBER10 package [45] using the ff99SB parameter set [46] for the protein and GAFF parameters [47] for the co-factor. Partial charges for the co-factor were derived by RESP fitting [48] at HF-6/31G*-level using Gaussian03 [49]. Van der Waals parameters for Fe(III) were derived by distance scanning against a helium atom at HF-6/31G*-level and found to be r = 1.20 Å and e = 0.0140 kcal/(mol*cm) in agreement with other transition metals in ff99SB.
Glutamate as well as aspartate residues were deprotonated, lysine and arginine residues protonated. Histidine (His) protonation was carried out manually according to the respective chemical environment: His200, His201, His208 and His271 at the protein surface were protonated at both imidazole nitrogens, His146, His170, His264 were assigned e-hydrogen atoms. His285 and His290 coordinating the Fe(III)-ion at the active center were assigned d-hydrogen atoms. The total charge of the system was found to be +3, which was neutralized by a uniform neutralizing plasma for Particle Mesh Ewald simulations [50].
The solvent was treated explicitly in the simulations including all resolved water molecules from the X-ray structure. The system was soaked with a wall distance of 10 Å by superposing an octahedral box of 7554 TIP3P [51] water molecules, resulting in a box edge of 70.8 Å after NpT equilibration.
After relaxation with harmonic restraints on heavy protein atoms and subsequent minimization, the systems were heated gradually from 100 K to 300 K over 200 ps in NVT ensemble. After heating and 1 ns of pressure and density equilibration, a NpT MD production run of 200 ns was performed to allow for larger conformational rearrangement in the enzyme. Simulations were carried out at 300 K using a Langevin thermostat [52] at 1.0 bar pressure with 8.0 Å nonbonded cutoff. SHAKE algorithm [53] on hydrogen atoms allowed a 2 fs time step. Snapshots were saved to trajectory every 500 steps or equivalent 1 ps for further analysis.
REMD simulations increase conformational sampling compared to standard MD simulations by exchanging replica of the simulated system at different temperatures, if an energy criterion is met. Hence, energy barriers are more likely to be overcome extending the sampling of conformational space [54]. In this study, equilibrated systems of standard MD simulations were used as input coordinates for two REMD simulations (PAH nat , PAH ox ). Sixteen parallel simulations were performed for both systems in a temperature spacing of 2 K as suggested by Patriksson et al [55] in a temperature range from 300 to 330 K. Exchange probability was set to 0.25 whilst attempting 10 exchanges per nanosecond. Besides a change to NVT ensemble, which is justified after volume adaption during pressure equilibration, identical parameters as for standard MD simulations were applied for REMD simulations. Production runs were performed for the two systems at 16 temperatures for 100 ns, resulting in total REMD trajectory length of 1.6 ms each.
Trajectories were analyzed using ptraj (version 4/2010) from AmberTools [45]. Positional fluctuations were calculated as RMSD of Ca-atoms to assess stability of the simulations. Bfactors (temperature factors) were calculated residue-wise as a measure of individual local flexibility. Structures were visualized using pymol [56]. This software package was also used to calculate the accessible surface area of the catalytic Fe(III)-ion over the simulation time by the command 'get_area'. Binding site volumes were estimated using POVME [57] by creating a 6 Å sphere around the center Fe(III)-ion filling the active site cavity with 1 Å spaced grid points. Subsequently, grid points near protein atoms were removed yielding an estimate of the actual binding site by summing up volume elements of residual grid points.

Results
After equilibration, for both simulated systems stable molecular dynamics trajectories over 200 ns with mean backbone Ca RMSD values below 2 Å could be generated (see Figure 2). After 20 ns simulation PAH ox shows a sharp increase in backbone dynamics as a loop region around residue 140 distant from the artificially linked cysteine residues starts to rearrange. Despite this structural transition in the so-called Tyr138-loop [58][59][60], simulation PAH ox yields a stable trajectory for the whole simulated system resulting in an RMSD plateau below 2 Å at comparable level to simulation PAH nat .
Calculated B-factors also highlight the Tyr138-loop as region of particular interest with B-factors reaching beyond 1000 Å 2 in PAH ox compared to 300 Å 2 in PAH nat . Other regions show less pronounced flexibility and similar behavior in both MD simulations. This finding also holds true for the vicinities of Cys203 and Cys334, where the artificial oxidation was introduced in PAH ox , which are found to be rigid in both simulations.
As indicated by the distance plot in Figure 3, the flexible Tyr138-loop performs a movement towards the active site of the protein after an initial phase of elevated dynamics in PAH ox (15-70 ns) leading to a contraction of the binding site. This movement is not present in the simulation of the native state PAH nat . This movement in PAH ox is paralleled by a decrease in accessibility of the catalytic iron center indicated by lowered accessible surface area of below 10 Å 2 in PAH ox compared to over 15 Å 2 in PAH nat (see Figure 3 and Table 1). Furthermore, the total binding site volume is reduced significantly in the course of this loop rearrangement in PAH ox : An average active site volume of over 150 Å 3 in PAH nat is reduced to less than half in PAH ox by the introduced perturbation.
REMD simulations yielded comparably stable trajectories (see Figure S1), confirming trends observed from for standard MD simulations. B-factors in the Tyr138-loop are generally elevated in PAH ox compared to PAH nat (see Figure S2). This observation is true for every single direct comparison of 16 temperatureequivalent replicas. Likewise, PAH ox shows a lowered average accessibility of the catalytic iron of 13.7 Å 2 compared to 15.4 Å 2 in PAH nat . The volume of the active site is decreased similarly in the oxidized state from 116.2 Å 3 in PAH nat to 104.8 Å 3 in PAH ox .

Discussion
As a structural explanation of the impairment of PAH caused by oxidative stress has been lacking by now, we propose a computational model of an oxidized PAH specie. This protein structure shows an altered behavior in MD simulations in the Tyr138-loop in close proximity to the active site [61] leading to a stable structure showing a contracted binding site. As the loop rearrangement occurs as a single event in the 200 ns MD trajectory, the statistical significance of this observation cannot be inferred from standard MD simulations of that time scale [62]. This intramolecular movement hence appears on a time scale beyond sampling of standard state-of-the-art MD simulations.
To statistically substantiate this proposed structural rationale for PAH down-regulation, REMD simulations providing improved conformational sampling were performed. These computationally demanding simulations agreed with trends seen in standard MD simulations, namely elevated loop dynamics paralleled with a decrease in accessibility of the catalytic iron center and binding site contraction. This contraction of the binding pocket would obviously lead to reduced substrate binding and hence enzyme down-regulation or even dysfunction. The enhanced sampling by REMD simulations decreased the observed differences between the simulations of PAH nat and PAH ox , but still supported findings from standard MD simulations, as 16 additional parallel simulations confirmed the trends observed from a single standard MD run.
Findings from MD simulations linking dynamics of the Tyr138loop region with catalytic activity are in best agreement with experimental data on allosteric regulation of PAH. Li and coworkers showed by an increase in exchange kinetics of hydrogen/ deuterium exchange mass spectrometry, that region 129-143 is indeed a solvent-accessible and highly flexible region in PAH that is influenced by presence of Phe [63]. This observation is further confirmed, as this region is not resolved in several crystal structures of PAH in absence of a substrate stabilizing the active site. Furthermore, modulation of the dynamics of region 129-143 by incubation with phenylalanine, an activator of PAH, indicates a mechanistic role in allosteric regulation [63].
The critical regulatory role of surface loops has been discussed in a recent review on allosteric regulation of PAH by Fitzpatrick [60]. He describes the Tyr138-loop as undergoing the most drastic change upon ligand binding in accordance with crystal structure analyses of Andersen et al [58,59]. Tyr138 included in that loop closes down on the active site leading to a major contraction restricting further access similar to the movement observed in our simulations (see Figure 4). A similar movement of a surface loop has been proposed by Sura et al. for the homologous enzyme tyrosine hydroxylase [64].
This dynamic behavior of the loop region could be seen as an analogy to the behaviour of the autoregulatory sequence of PAH. This domain of the protein physically blocks entrance to the active site of PAH in the inactive form. A smaller but similar impact of the surface loop can be inferred from presented molecular dynamics simulations directly explaining the observed PAH impairment upon oxidative stress.
The importance of conformational flexibility for PAH activity has been revealed in several studies: Cerreto et al and Gersting et al showed correlations between PAH conformational stability and PAH down-regulation by PKU-causing mutations [65,66]. Changes in stability were shown by MD simulations upon phosphorylation of Ser-16 [67]. This computational technique was also used to infer factors for specificity and affinity of PAH ligand binding [68,69].
As presented MD simulations allow to keep track of alterations in molecular motions within the protein, a pathway for the inferred allosteric signaling between Cys203/Cys334 and the active site could be anticipated. An analysis of hydrogen bonding within the protein did not highlight a particular alteration in molecular interactions within the protein but rather suggested a concerted flow of dynamic behavior throughout the protein backbone of PAH (see Figure 5). The differences in molecular flexibility between PAH ox and PAH nat REMD simulations are highlighted by a color-coding from blue for regions more flexible in PAH nat over white to red for regions more flexible in PAH ox . The dominating red coloring of PAH ox starting from the artificial oxidation site propagates over the protein backbone towards the surface loops in close proximity to the active site, whereas only minor blue regions are found in a small part of the binding site next to the biopterin co-factor. This highlights a putative in- Figure 3. Dynamics of the binding site region of native and oxidized phenylalanine hydroxylase. A: Distance between the catalytic iron center and the Ca-atom of Glu141 included in the flexible Tyr138-loop near the active site. Whereas this distance remains mostly stable in the native state PAH nat over the simulation time (blue), simulation of PAH ox shows an initial phase of increased loop flexibility that is followed by decrease of that distance indicating a loop movement towards the active site. This finding is confirmed by calculation of the accessible surface area of the iron center (B). After an initial adaption to the perturbation, PAH ox shows a constantly reduced accessibility of the catalytic iron center compared to PAH nat . The overall active site volume (C) is similarly reduced over simulation time in PAH ox compared to PAH nat , suggesting a reduced enzyme turnover. doi:10.1371/journal.pone.0053005.g003 teraction path for the allosteric effect of disulfide formation on the active center via the central helices directly to the Tyr138-loop near the active site. A loop movement near the active site in turn restricts accessibility of the catalytic site, hence explaining reduced enzyme turnover. Our study offers a new example in the emerging field of redox signaling in enzyme activity, that is facilitated by novel findings concerning molecular mechanisms of redox switches in proteins [70]. Schroeder et al for example recently highlighted how the reduction of disulfide bonds activates an antimicrobial protein [71]. Tanner et al. recently reviewed the redox regulation of protein tyrosine phosphatases, where an active site cysteine takes several oxidation states and is inactivated by redox signaling [72]. A regulatory mechanism similar to the proposed signal propagation in PAH has been found by Wang et al. in C-reactive protein, where a single disulfide bond serves as a switch for the function of the protein [73].
An obvious limitation of the presented study is the underlying simplistic model of PAH ox . Although reduced Cys203 and Cys334 in close proximity are a clear weak-point for oxidation in the protein, this model might not cover all effects of oxidative stress. Cysteine residues might undergo further oxidation processes besides simple disulfide formation. On the other hand, also the active site of PAH is susceptible to oxidative species, as the enzyme co-factor is readily oxidizable under oxidizing conditions. These latter two effects could not be covered in the presented simulations and hence might be starting points for further investigations from theoretical and experimental side.
In conclusion, our study provides a structural basis for the down-regulation of PAH activity upon oxidative stress. MD simulations highlight a flexible loop region around residue 140 (Tyr138-loop) as susceptible to side chain oxidation in PAH. A complex movement over the whole protein leads to a contraction of the binding site leading to the experimentally observed reduced enzyme activity. . The pronounced loop movement towards the catalytic center (iron ion as brown sphere, co-factor as sticks) was measured as distance of the Ca-atom of Glu141 (red and blue sphere respectively) to the catalytic iron (brown). This rearrangement is paralleled by a reorientation of Tyr138 (shown as sticks, left), in turn concertedly blocking access to the active site of PAH ox . doi:10.1371/journal.pone.0053005.g004