Modelling the elusive conformational activation in kynurenine 3-monooxygenase

Abstract

Kynurenine 3-monooxygenase (KMO) regulates the levels of important physiological intermediates in the kynurenine pathway [Guillemin, et al., Journal of Neuroscience, 2007, 27, 12884], which is the major route for L-tryptophan catabolism. Its catalytic activity (hydroxylation) is dependent on the formation of a short-lived intermediate that forms after the reduction of the coenzyme FAD. The reduction takes place fast when the substrate binds to KMO. Crystal structures of the apo form and in complex with an effector inhibitor, which prevents the hydroxylation of the substrate but also stimulates KMO like the substrate, and a competitive inhibitor, which suppresses the substrate hydroxylation, are available for the resting in conformation only. The active out conformational state that enables the reduction of FAD at an exposed location of KMO after its stimulation by an effector, however, was implicated but not resolved experimentally and has remained elusive so far. Molecular dynamics simulations of apo KMO and the inhibitor–KMO complexes are carried out using extensive multi-dimensional umbrella sampling to explore the free-energy surface of the coenzyme FAD's conformational conversion from the in state (buried within the active site) to the out state. This allows a discussion and comparison with the experimental results, which showed a significant increase in the rate of reduction of FAD in the presence of an effector inhibitor and absence of enzymatic function in the presence of a competitive inhibitor [Kim, et al., Cell Chemical Biology, 2018, 25, 426]. The free-energy barriers associated with those conformational changes and structural models for the active out conformation are obtained. The interactions during the conformational changes are determined to identify the influence of the effector.

---

1 Introduction

Tryptophan is catabolized in the kynurenine pathway,¹ in which L-kynurenine (L-Kyn) is the branching point for several pathways with contrasting biochemical effects. Kynurenine 3-monooxygenase (KMO)² hydroxylates L-Kyn and releases the product 3-hydroxykynurenine (3-HK); this itself and its subsequent metabolites (3-hydroxyanthranilate (3-HanA) and quinolinate (Quin)) are responsible for neurodegenerative diseases.^1,3–5L-Kyn is also a substrate to kynurenine aminotransferase (KAT), which converts it to the neuroprotective compound kynurenic acid (KynA).^6–11 Decreasing the Quin/KynA ratio by inhibiting KMO is regarded as a possible intervention point to treat neurodegeneration. This point of intervention allows the targeted treatment of central-nervous-system damage.

Several successful inhibitors were developed¹² over the years, but it was soon realized^12,13 that most of the inhibitors themselves cause the release of neurodegenerative hydrogen peroxide (H₂O₂). Although inhibiting KMO to prevent neurodegeneration and neuroinflammation, which cause diseases like Huntington's,³ Alzheimer's,⁴ and Parkinson's,⁵ seems to be a straightforward strategy, formation of H₂O₂ complicates this therapeutic route. Therefore, understanding the role of an effector on the activation of KMO is crucial in the development of optimal inhibitors.

The requirement of substrate binding for the fast reduction with nicotinamide adenine dinucleotide phosphate (NADPH) is a common feature of hydroxylases like KMO of the class A monooxygenase family.²L-Kyn binding to the enzyme initiates the reduction of FAD's isoalloxazine ring system by NADPH^13,16 (Scheme 1). According to kinetic studies,¹³ the enzymatic rate constant of the reduction of L-Kyn–KMO is 2500-fold in comparison to that of free KMO (0.0020 s⁻¹). ‘Non-substrate effector’ ligands that inhibit the hydroxylation were also shown^13,16 to stimulate the reduction of the isoalloxazine ring system.¹⁶ In the presence of such a non-substrate effector, the highly reactive hydroperoxide intermediate (Scheme 1) decays to the oxidized FAD by releasing harmful H₂O₂, since the oxygen from the reactive hydroperoxide intermediate is not transferred to the substrate.^13,16

Scheme 1 Reaction mechanism¹³ of KMO. L-Kyn binding to KMO initiates FAD_ox reduction by NADPH, then molecular oxygen binds to FAD's 4a position, resulting in the formation of the reactive hydroperoxide intermediate that is transformed into FAD_ox upon release of H₂O and the hydroxylation of L-Kyn, which yields the product 3-HK. The oxygen transfer step has been computationally investigated^14,15 in detail. 3-HK is released in the rate-determining step of the enzymatic cycle.

In p-hydroxybenzoate hydroxylase (from Pseudomonas aeruginosa), which is also a class A monooxygenase, this initiation was suggested to take place via a substrate-induced conformational change.¹⁷ According to the suggested mechanism,¹⁸ upon substrate binding, the enzyme ‘senses’ the presence of the substrate via an H-bonding network, then a conformational change takes place, involving the movement of FAD's isoalloxazine ring system from its buried position (the in conformational state) to another position (the out conformational state),^17,19 where it is closer to the surface of the enzyme and is reduced by the bulky NADPH (Fig. 1). KMO, however, does not have such an H-bonding network and was never captured in the out conformational state.^13,16 This implies that it utilizes a different activation mechanism in comparison to similar flavoproteins and its in conformational state appears to be more thermodynamically stable. Thus, initiation of KMO reduction must take place via a different substrate-induced activation mechanism that also involves the in to out conformational transition, since the possibility of the reaction between the two large molecules, FAD and NADPH, within the narrow active site is not feasible.

Fig. 1 Structural details of the two co-crystallized conformational states of the FAD of p-hydroxybenzoate hydroxylase from Pseudomonas aeruginosa (pdb id: 1DOB). The isoalloxazine is buried within the active site in the in conformation (red), while it is closer to the surface of the enzyme in the out conformation (blue), where it can react easier with NADPH. In the crystal structure, both conformational states were captured simultaneously with a ratio in : out = 24 : 76.

The X-ray structure of human KMO (hKMO) was recently published.¹⁶ As predicted,²⁰ in the absence of its membrane-targeting domain, which binds to the mitochondria, the enzyme lost its functionality. KMO structures from Pseudomonas fluorescens (pfKMO)^21,22 and Saccharomyces cerevisiae (scKMO)^16,23 were also published. The release of H₂O₂ was detected in the presence of the substrate L-Kyn or the inhibitors BA or UPF-648 with hKMO, scKMO or pfKMO.¹⁶ Reported²³ scKMO structures enabled a comparison between the apo (4J33) and holo (4J36) forms of the enzyme (Fig. 2). scKMO is a superior model system¹⁶ for hKMO because neither hKMO nor scKMO cause H₂O₂ formation while pfKMO does when hydroxylation is inhibited by Ro 61-8048. Since truncated hKMO is inactive, scKMO must be used in structure-based drug-discovery studies.¹⁶

Fig. 2 Superposition of scKMO crystal structures, in the presence of the effector inhibitor (Inh_E) UPF-648 (pdb id: 4J36) and in its absence (pdb id: 4J33). The positions of the Phe322 (orange) and Tyr323 (green) residues that are part of a loop above the isoalloxazine (cyan) in UPF-648–KMO are different in comparison to their original positions in apo-KMO, where they are shown in red and blue, respectively. This change is induced by the steric clashes due to the presence of the dichlorobenzyl group of the inhibitor.² The C-terminal α helix (in red in apo-KMO) was not resolved in UPF-648–KMO, but it must be shifted due to its steric clash with Tyr323. These conformational changes were related to the H₂O₂ production observed in the presence of UPF-648–KMO. The 2D structure of UPF-648 is shown on the right.

In the loop above the isoalloxazine, residues Phe312 and Phe313 (hKMO numbering), which have aromatic side chains, correspond to Phe322 and Tyr323 in scKMO, and Phe319 and His320 in pfKMO.¹⁶ In all KMOs mentioned, the reduction was completely abolished¹⁶ when these two residues were replaced by Ala. As shown in Fig. 2, this relation is apparent^2,23 from analysis of the X-ray structures of an scKMO complex that contains a non-substrate effector UPF-648 inhibitor (Inh_E) and apo-KMO. In contrast, the competitive inhibitor (Inh_C) Ro 61-8048 does not cause such structural changes in scKMO.¹⁶

In this work, non-equilibrium molecular dynamics (MD) simulations with extensive and multi-dimensional umbrella sampling^24,25 (US) are used to model the in to out conformational activation step in KMO. The investigation of the thermodynamic driving forces and transition barriers of this conformational change is critical for the entire KMO mechanism. These simulations showed that conformational activation has a high transition barrier height in the apo form, is facile in the complex with a non-substrate effector and is abolished in the complex with a competitive inhibitor.

2 Results and discussion

2.1 The in to out FAD conformational transition in apo-KMO

In the text, the resulting computational systems associated with the structures apo-scKMO, UPF-648–scKMO and Ro 61-8048–scKMO are concisely referred to as apo-KMO, Inh_E–KMO and Inh_C–KMO, respectively, while the atomic positions of the isoalloxazine ring system are referred to as shown in Scheme 1. Please see the “Details of non-equilibrium MD simulations” subsection for further details and definitions of the collective variables X and Y. As described in detail therein, the 2D X–Y space is spanned by using US simulations at two levels of accuracy: ‘fast’ US simulations refer to initial scans by simulations in which the windows are at an interval of 5° and each window is run for 1 ns; in the more accurate ‘slow’ US simulations, they are at an interval of 1° and each one of them is run for 7 ns after 1 ns of equilibration.

According to the analyses performed using the Weighted Histogram Analysis Method (WHAM, in Fig. 3) for the ‘fast’ US(X) simulations, the system can evolve through a saddle point only at values of Y = 145°, 150°, 155° and 180°. For other values of Y, the free energy of the system kept increasing continuously. Dihedral angles of Y = 145° and 180° were selected for the subsequent more accurate ‘slow’ US(X) simulations. The simulations at Y = 150° and 155° were not directly utilized as these were close to the simulations at Y = 145°, and higher in energy. However, during the generation of the two-dimensional local free-energy surface, the points near to the maxima in these ‘fast’ simulations were scanned (Fig. S14†). According to Fig. 4, the run at Y = 180° (‘slow’ US(X, 180)) does not reach a maximum as it keeps increasing, while the run at Y = 145° (‘slow’ US(X, 145)) displays a maximum at X = −115° with a relative free energy of +11.3 kcal mol⁻¹. However, the free-energy profile of the ‘slow’ US(−175, Y) run on the right-hand side of Fig. 4 reveals an important point: in this run, the free energy at (X, Y) = (−175°, 145°) is 6.4 kcal mol⁻¹ higher compared to the point at (−175°, 180°). Therefore, calculation of the free-energy barrier requires the combination of several data.

Fig. 3 Changes in free energy from X = −175° to X = −85° for the ‘fast’ US(X) scan of apo-KMO, obtained for different values of dihedral angle Y at an increment of 5°. The windows of the collective variable of US (dihedral angle X) are also at an interval of 5° for each run (19 = 19 × 1 ns). The change in free energy is given relative to the minimum obtained for each run.

Fig. 4 Changes in free energy from X = −179° to X = −85° in the ‘slow’ US(X, 145) and US(X, 180) runs (left-hand side), and from Y = 145° to Y = 194° in the ‘slow’ US(−175, Y) run (right-hand side) for apo-KMO. Windows are at an interval of 1° for each run. The total production runtime is 665 = 95 × 7 ns for each US(X), and 350 = 50 × 7 ns for the US(Y).

The points of the slow US(X, 145) run were used as starting conformations in the search for the saddle point, since the US(X, 180) run does not reach a maximum. The details about how the contour plot were constructed can be found in Fig. S14.† In Fig. S14,† the saddle point is clearly visible between X = [−119°, −118°] and Y = [180°, 181°]. In addition, three degenerate shallow local minima are located in this contour map. The one that is between X = [−114°, −113°] and Y = [183°, 184°] is furthest in the X-angle direction. This shallow minimum is a more likely candidate for the out conformational state, since the X-angle is the leading coordinate in the conformational change. Further exploration of the free-energy surface was initiated from (X, Y) = (−111°, −179°) in Fig. S14† using ‘slow’ US(X). As shown in Fig. 5, a local minimum can be obtained at (X, Y) = (−84°, −179°). However, a detailed contour analysis (Fig. S15†) in the vicinity of this region showed that this minimum can only correspond to a shallow out conformational state, and the free energy sharply increases beyond this point. In addition, this local minimum between X = [−87°, −86°] and Y = [180°, 181°] is higher in free energy compared to the shallow local minimum between X = [−114°, −113°] and Y = [183°, 184°] in Fig. S14.† Therefore, these correspond to two shallow out conformational states, and the conformer in the latter minimum is more stable.

Fig. 5 The free-energy profile of US(X, −179) from X = −111° (as obtained from the contour map shown in Fig. S14†) to X = −70°. The free-energy minimum is found at (X, Y) = (−84°, −179°). The free energy increases sharply beyond this point.

The in conformational state of apo-KMO (Fig. 6) was obtained under the physiological conditions of the H₂O₂ detection experiments,¹⁶ starting from selected coordinate points from the slow US(X, 180). The most stable in conformation was found between X = [−165°, −164°] and Y = [181°, 182°].

Fig. 6 Free-energy contour map for the in conformational state of apo-KMO. The in conformational state is located within the area between X = [−165°, −164°] and Y = [181°, 182°]. Both in increasing and decreasing directions in free energy, ten ‘slow’ US(Y) simulations (770 = 10 × 11 × 7 ns) were carried out for X = [−169, −160] using the output conformers of slow US(X, 180).

Obtaining the free-energy barrier and free-energy change accompanying the conformational change (Fig. 7) can only be achieved indirectly by combining several scan data (‘Details of Relative Free Energy Calculations of KMO Conformational States’ in the ESI†).

Fig. 7 The reaction coordinate of FAD's conformational change in apo-KMO. The vertical axis is the relative free-energy change (kcal mol⁻¹). The representative structures of the most populated clusters displayed were obtained from cluster analysis using Amber's CPPTRAJ. The four trajectory files analyzed for each structure are the ones that span the extremum points defined in each caption of the contour map (Fig. 6 and S14†). The residues that are <3.5 Å from any of the selected atoms of FAD (N1, O2, H3, O4, N5, H6, 7a Hs, 8a Hs and H9) during at least half of the time of the respective four trajectory files are shown, while nonpolar hydrogen atoms are not displayed.

Polar sites of the isoalloxazine are stabilized by five H-bonds (dashed lines between the atoms in Fig. 7) when the system is at its global minimum, the in conformational state (Fig. 7). In addition, the numbers of water molecules surrounding each polar site of the isoalloxazine are listed in Table S1.† O2 is not in contact with water at all since it forms two H-bonds with Asn328, while O4 and N5 are surrounded by several water molecules, indicating the presence of an H-bonding network in their surroundings. The side chain of Phe322 of the loop above the isoalloxazine bends towards the space in front of N5. This space will be occupied by the isoalloxazine ring system during the conformational change. On the other hand, the side chain of Phe246 bends towards the same space and forms a stacking interaction with Phe322's side chain. The in to out conformational change of the isoalloxazine is blocked through this interaction. In addition, the movement of the isoalloxazine in the vertical direction is restricted by another residue in the loop above the isoalloxazine, Pro321, which is seen above the ring with methyl substituents in the isoalloxazine.

In the TS (Fig. 7), H-bonding interactions are retained to a great extent and the isoalloxazine is still confined within the active site, as can be seen from the number of water molecules surrounding its polar sites (Table S1†). These observations are in line with the fact that the high-energy TS is early along the reaction coordinate. The side chain of Tyr195 blocks the movement of the isoalloxazine ring system. This side chain cannot occupy a position to the left because it would clash with Phe246. Pro321 of the loop is on top of the ring with methyl substituents in the isoalloxazine. There, it restricts its alignment in the vertical direction (Fig. S20A and B†).

In the out conformational state (Fig. 7), the number of H-bonds of the isoalloxazine is significantly reduced to one, while the number of water molecules surrounding it is increased (Table S1†). As exemplified in Fig. S21,† the out conformational state lies at a high and shallow relative free-energy level due to the repelling interactions that the isoalloxazine must form to maintain it. As expected from the free-energy contour map in Fig. S14,† the repelling intramolecular interactions start to emerge as the isoalloxazine further rotates in the reaction coordinate (Fig. S21†).

Overall, in the apo form of the enzyme, the side chain of Phe322 from the loop was shown to interact with Phe246's side chain with stacking interactions, blocking the way in which FAD's isoalloxazine ring system was supposed to move. As shown in Fig. S22,† this interaction is considerably weakened during the transition state as the closest distance between the two rings increases from 3.6 to 5.4 Å. However, this separation is only enabled by a high-energy TS that is obtained at an early stage of the conformational change in the X-angle direction. Tyr195 is a blocking residue that obstructs the conformational change, because the loop above the isoalloxazine that includes the Pro321 residue prevents the movement of the isoalloxazine in the vertical direction (Fig. S20A and B†) and makes the escape from the repelling interaction with Tyr195 difficult. In the TS, there is a repelling H–H interaction at an H-bond distance, 2.6 Å, due to this confinement (Fig. S20C†). As a result, the preservation of the position of the loop above the isoalloxazine is critical in the high free-energy barrier.

2.2 An effector inhibitor lowers the free-energy barrier of FAD activation

According to ‘fast’ US(X) simulations (Fig. 8), the system can evolve towards a saddle point at Y = −170, −175, 180, 175, 170, 165, 155, 150 and 145° in the presence of the effector inhibitor molecule Inh_E.

Fig. 8 Changes in free energy in kcal mol⁻¹ between X = −175° and X = −85° for the ‘fast’ US(X) scan (19 = 19 × 1 ns) of the Inh_E–KMO complex obtained from WHAM, at various values of Y.

Values of Y = 155° and 170° were selected for further ‘slow’ US(X). A saddle point could not be localized (Fig. S16†) utilizing the coordinate points of US(X, 155) in Fig. 9. The US(X, 170) in Fig. 9 presents a TS candidate that is followed by a deeper well, and starting from the coordinate points obtained from US(X, 170), a saddle point is clearly found between X = [−113°, −112°] and Y = [170°, 171°] (Fig. 10). The global minimum that corresponds to the in conformational state was found between X = [−169°, −168°] and Y = [163°, 164°] (Fig. 10) while the product state that corresponds to the out conformational state was found between X = [−80°, −79°] and Y = [166°, 167°] (Fig. 10).

Fig. 9 Changes in free energy in kcal mol⁻¹ between X = −179° and X = −85° in the ‘slow’ US(X, 155) run (665 = 95 × 7 ns), between X = −179° and X = −60° in the ‘slow’ US(X, 170) run (840 = 120 × 7 ns), and between Y = 136° and Y = 180° in the ‘slow’ US(−175, Y) run (315 = 45 × 7 ns) for the Inh_E–KMO complex.

Fig. 10 Free-energy contour maps of the in state (left), TS (middle) and out state (right) of Inh_E–KMO. Each contour map was obtained using the points of the ‘slow’ US(X, 170). The in state that is between X = [−169°, −168°] and Y = [163°, 164°] was obtained via nineteen slow US(Y) runs (1862 = 19 × 14 × 7 ns) for X = [−176, −158]. The saddle point that is between X = [−113°, −112°] and Y = [170°, 171°] was obtained via ten ‘slow’ US(Y) runs (770 = 10 × 11 × 7 ns) for X = [−120, −111]. The out conformational state that is between X = [−80°, −79°] and Y = [166°, 167°] was obtained via ten slow US(Y) runs (770 = 10 × 11 × 7 ns) for X = [−83, −74].

WHAM gives the free-energy barrier and the relative free energy of the out conformational state (Fig. 11) as explained in the ESI (‘Details of Relative Free Energy Calculations of KMO Conformational States’†). As expected,^13,16 the conformational change is accompanied by an endergonic free-energy change.

Fig. 11 The reaction coordinate of FAD's conformational change in Inh_E–KMO. The selection scheme applied to visualize the surrounding residues of the isoalloxazine is explained in the caption of Fig. 7.

In the in conformational state, the isoalloxazine is stabilized by six H-bonds with the surrounding residues (Fig. 11) and it is more confined from water in comparison to apo-KMO's in conformational state, since only N5 is closely interacting with a water molecule (Table S1†). Pro321 is not in the range of the selection scheme applied in Fig. 11. A stacking interaction is formed between the rings of the inhibitor and the side chain of Phe246 like the interaction between Phe322 and Phe246 in apo-KMO. Through this interaction, the movement that should take place for the in to out conformational change is blocked, but with one important difference: with Inh_E being a small molecule and loosely bound, it can move easier in comparison to Phe322, whose position is mostly dictated by the protein loop. The change in the position of the inhibitor can be visually inspected via the superposition of the in conformational state and TS shown in Fig. S22.†

In the TS of Inh_E–KMO (Fig. 11), the isoalloxazine forms significantly fewer (only two) H-bonds and is surrounded by more water molecules in comparison to the TS of apo-KMO. These observations are in line with the later location of the TS in the reaction coordinate in comparison to that of apo-KMO. Interestingly, interactions with water molecules increased for every polar site but N5 (Table S1†) in comparison to the in conformational state. As seen for the TS of apo-KMO, Tyr195 does not form an H-bond with the isoalloxazine ring system and remains perpendicular at its methyl-substituent end. The transition barrier energy is mostly due to this blocking effect of Tyr195. However, unlike apo-KMO, once it rotates further along the reaction coordinate, it relaxes to a deep local minimum, as seen in Fig. 10.

The TS barrier in Inh_E–KMO is reduced in comparison to that of apo-KMO because the space that is occupied by the loop above the isoalloxazine becomes available for the isoalloxazine and it is efficiently occupied (Fig. S20A and B†) and thereby the blocking interaction due to Tyr195, which is also seen in the TS of apo-KMO, is less severe. In apo-KMO, the side chain of the loop residue Pro321 prevents the movement of the isoalloxazine in the vertical direction, and therefore, it cannot escape from the blocking of Tyr195 as easily as seen in Inh_E–KMO. In the TS of Inh_E–KMO, the repelling H–H interaction is absent due to the availability of space above the isoalloxazine, which bends into it (Fig. S20C and D†). As shown in Fig. S20A,† if the Pro321 residue was not displaced as in apo-KMO, the distorted isoalloxazine's H and Pro321's H would clash. In addition, the distance between Inh_E and Phe246 increases from 3.3 to 5.2 Å (Fig. S22†) with relative ease, due to the mobility of Inh_E (Fig. S22†).

In the out conformational state of Inh_E–KMO, in which the isoalloxazine is in an open position (Fig. 12), there is only one H-bond involving the isoalloxazine: between O4 and Tyr195's side-chain polar hydrogen (Fig. 11). Tyr195's side chain occupies the space in front of the isoalloxazine ring system and is to the left side of N5. In the out conformational state, the number of direct interactions with water molecules increases for O2 and remains constant for all other atoms compared to the TS. N5's relatively low number of interactions with water molecules is interesting, considering its immediate solvent exposure. This is the pre-reactive orientation for the FAD reduction. As the thermodynamics in this work shows, the out conformation is higher in free energy and probably short-lived and expected to reconvert back to the more stable in conformational state. However, its shielding from direct interactions with solvent molecules can stabilize the critical hydroperoxide intermediate that might form after the reduction within its short lifetime.

Fig. 12 The out conformational state of Inh_E–KMO in which the isoalloxazine is clearly in an open position. The utilized frame is the same as in Fig. 11 but with a different selection scheme; only the residues that are within 5 Å of N5 of FAD are shown explicitly.

2.3 A competitive inhibitor suppresses the conformational activation of KMO

In the presence of the competitive inhibitor Inh_C, the ‘fast’ US simulations (Fig. 13) cover a wider range of dihedral angles X in comparison to the previous two complexes. The runs at Y = −175, −150 and 175 were selected for further ‘slow’ US(X) refinement simulations (Fig. 14) since those were indicative of a saddle point with a reasonable barrier height. According to the corresponding contour maps (Fig. S17–S19†), there is neither a saddle point nor a local minimum in these ranges. This agrees with experiments,¹⁶ which showed that FAD reduction is not taking place in the presence of a competitive inhibitor molecule. Nonetheless, an area from the contour of Fig. S19† will be investigated to fully rationalize the inhibitor's mode of action. As displayed in Fig. 15, the in conformational state is located between the points X = [−178°, −177°] and Y = [−169°, −168°].

Fig. 13 Changes in free energy in kcal mol⁻¹ for the scan from X = −180° to X = −85° in the ‘fast’ US(X) run (20 = 20 × 1 ns) for the Inh_C–KMO complex at various values of Y.

Fig. 14 Changes in free energy in kcal mol⁻¹ for scans from X = −180° to X = −85° (672 = 96 × 7 ns) in the ‘slow’ US(X, −150), US(X, −175) and US(X, 175) runs, and from Y = 141° to Y = 225° in the ‘slow’ US(180, Y) run (595 = 85 × 7 ns) for the Inh_C–KMO complex.

Fig. 15 Free-energy contour map for the in conformational state of the Inh_C–KMO complex. The in conformational state is localized between X = [−178°, −177°] and Y = [−169°, −168°]. Using the slow US(X, −175) run, both in the decreasing and the increasing directions in free energy, ten slow US(Y) simulations (980 = 10 × 14 × 7 ns) were carried out for X = [−178°, −169°].

In the in conformational state (Fig. 16), the isoalloxazine forms six H-bonds with KMO residues. These are similar to the ones in Inh_E–KMO, and these polar sites do not directly interact with the solvent molecules as in the in conformation for Inh_E–KMO. The ring of the inhibitor Inh_C aligns horizontally relative to the isoalloxazine ring system allowing the formation of a ‘parallel offset’ stacking interaction.²⁶ The position of the inhibitor is enhanced by additional stacking interaction with the side chain of Phe322 of the loop above the isoalloxazine (Fig. S23†).

Fig. 16 A. The in conformational state of FAD of Inh_C–KMO. The selection scheme applied to visualize the surroundings of the isoalloxazine is explained in the caption of Fig. 7. B. The same as in A, but from a different point of view.

In order to understand how Inh_C prevents the conformational transition of FAD, a high-energy region (located at X = [−112°, −111°] and Y = [174°, 175°]) on the free-energy contour map in Fig. S19† was analyzed (Fig. 17). The most significant result from Fig. 17 is the highly restrained conformation of the isoalloxazine in comparison to the in conformational state in Fig. 16. This is due to the unfavorable ‘parallel face-centered’ stacking interaction²⁶ between the rings of the inhibitor and the isoalloxazine. In contrast, the interaction between Phe322 of the loop and the inhibitor is conserved (Fig. S23†), therefore, the positioning of the inhibitor is dictated by the loop region above the isoalloxazine. This loop stabilizes the position of the inhibitor at where it is in the in conformational state while the putative rotational transition of the isoalloxazine mode is destabilized. Comparing this high-energy state to the in conformational state shows that the side chain of Tyr195 (Fig. S24†) is occupying a position close to the isoalloxazine. The side chain of Tyr195 is partly responsible for the restrained isoalloxazine binding mode in Fig. 17. The presence of the inhibitor between the isoalloxazine and the KMO loop further enhances the restriction of the rotation of the isoalloxazine in the vertical direction. Therefore, Tyr195 must be displaced to create space for the rotation of the isoalloxazine. In addition, Tyr195's side chain cannot move towards the left-hand side (Fig. 17) as it would clash with the inhibitor, so it moves forward in the positive X-angle direction of the conformational change and further clogs it (Fig. S24†).

Fig. 17 A. A structural view of a high-energy contour region of the Inh_C–KMO complex. The four runs cover dihedral angles ranging between X = [−112°, −111°] and Y = [174°, 175°] on the free-energy contour map seen in Fig. S19.† The selection scheme applied to visualize the surroundings of the isoalloxazine is explained in the caption of Fig. 7. B. The same as in A but from a different point of view.

3 Conclusions

The activation of KMO requires an enzymatic conformational transition towards the reactive out conformation. Extensive exploration of the free-energy surface by umbrella sampling revealed the free-energy differences between the inactive in and the reactive out conformational states of FAD in KMO and the transition-state barriers. The studies in this work confirmed the experimental results that showed^16,23 a relation between the positioning of the loop region above the isoalloxazine of FAD and the conformational changes.

In the apo form of the enzyme, this loop above the isoalloxazine is maintained at its original position of the in state, and the loop orientation is restricting the transition towards the out conformation in the vertical direction of the isoalloxazine via the interactions with a surrounding residue, Phe246.

The competitive inhibitor (Inh_C) fully inhibits the conformational change of FAD and a transition state connecting the in and out KMO states could not be localized. Furthermore, the structural restraints due to the presence of the inhibitor do not allow the KMO enzyme to occupy an out conformational state. The interactions between the inhibitor and the isoalloxazine, and the inhibitor and the loop residues, are critical for this inhibition.

The non-substrate effector (Inh_E) significantly lowers the free-energy barrier for the conformational activation of FAD. This lower transition-state barrier is enabled by the displacement of the KMO loop region above the isoalloxazine in the presence of the effector molecule. This enables a faster orientation of FAD towards the reactive conformation.

For the first time, non-equilibrium MD simulations and the free-energy changes show that the out conformational state can be localized and corresponds to a local minimum during KMO activation. It corresponds to a near surface position of the isoalloxazine of FAD in the enzyme, as was already implied by experiments but could not be verified.^13,16 Residues Tyr195, Pro321 and Phe322 are critical for the conformational transition regardless of the type of inhibitor or whether an inhibitor is present or not. However, in the presence of the effector inhibitor UPF-648, the positions of residues Pro321 and Phe322 are modified, relieving the repulsive interactions between the isoalloxazine and Tyr195 or Pro321 during the transition and preventing the stacking interaction between Phe322 and Phe246 (which blocks the transition).

The displacement of the KMO loop region by an effector molecule is essential for the conformational change to occur with reasonable rates. This displacement is induced by the vertical positioning of the ring of the inhibitor compound, relative to the isoalloxazine ring system of FAD. Human KMO responds to inhibitor binding in the same manner¹⁶ as Saccharomyces cerevisiae KMO. Therefore, the results of this study for scKMO are also transferable to human KMO. To avoid H₂O₂ formation during inhibition of human KMO, the inhibitors must be designed in such a way that their structural clash with the loop of KMO is avoided, and instead a constructive interaction between them is formed to maintain the binding mode of the isoalloxazine ring system.

4 Computational methods

4.1 Preparation of protein complexes

X-ray structures of scKMO were obtained from the protein databank. They include 4J33.pdb (resolution = 1.82 Å) in which the protein is in the apo form, 4J36.pdb (resolution = 2.13 Å) in which the inhibitor UPF-648 is present within the active site and 5X6R.pdb (resolution = 1.91 Å) in which the inhibitor Ro 61-8048 is in the active site. Although both inhibitors successfully inhibit the oxygen transfer step of the hydroxylation, UPF-648 is a non-substrate effector while Ro 61-8048 is a competitive inhibitor. Therefore, UPF-648 stimulates the reduction of FAD's isoalloxazine ring system, which further causes H₂O₂ formation, while Ro 61-8048 is only inhibiting the hydroxylation.¹⁶

Missing residues were re-constituted by homology modeling using SWISS-MODEL.²⁷ The missing C-terminal α-helix in 4J36 was modeled using the Builder module of PyMOL²⁸ and Amber20²⁹ sequential minimization. According to the Ramachandran plots, which were generated using Maestro,³⁰ these models are structurally reliable representations of each protein structure (Fig. S1–S3†). For each protein structure, SWISS-MODEL reported a “Ramachandran Favoured” value above 95%. The protonation states of His residues were manually assigned depending on their environments, while those of other residues were assigned automatically and all of them were checked using the H++ server.³¹

The topology and the coordinate files were generated using the tLEAP module of Amber20, and ff19SB³² was used as the force field. Parametrizations of ligands (FAD, UPF-648 and Ro 61-8048) were carried out using Amber's antechamber program suit.³³ The RESP model³⁴ was used to derive the atomic charges. The quantum mechanical calculations were performed with Gaussian 16 Revision A.03³⁵ for the generation of the electrostatic potentials. In these calculations, a B3LYP^36,37/6-31+G(d,p) functional/basis set combination was used in conjunction with the MK method.^38,39 During the optimizations, the heavy atoms were fixed to their X-ray coordinates. The parameters of the ligands were generated using GAFF2.⁴⁰

The protein was solvated with the TIP3P model⁴¹ of water in an octahedral simulation box. The minimum distance from the protein structure to an edge of the box was set to 15 Å. The protein complexes were sequentially minimized with a non-bonded cut-off distance of 12 Å (same throughout the simulations as well): (1) the water molecules were minimized with 10 000 steepest descent steps before switching to the conjugate gradient algorithm with a convergence criterion of 0.01 kcal mol⁻¹ Å⁻¹, while the protein was constrained with a weight of 250 kcal mol⁻¹ Å⁻². (2) The constraints were removed, and the entire system was optimized with the same minimization parameters. This procedure was applied to each computational system.

4.2 Details of MD simulations

All the MD simulations were performed using the cuda version of Amber20's pmemd program. Utilizing the Langevin thermostat, the equilibration was carried out in three steps in which the time step was set to 2.0 fs, the SHAKE algorithm was turned on for H-atom–heavy-atom bonds, and the random collision frequency was 2 ps: (1) starting from 0 K, the system was heated to 310 K within 0.1 ns with weak constraints (10.0 kcal mol⁻¹ Å⁻²) on the protein complex in the NVT ensemble. (2) The simulation was run with the same simulation parameters as in step 1, but without the constraints. (3) To bring the system to its physiological condition, the density of the system was equilibrated using the Langevin thermostat in an NPT ensemble for 100 ns. According to Fig. S4–S6,† each protein complex was well equilibrated.

4.3 Details of non-equilibrium MD simulations

As can be seen for the FAD of p-hydroxybenzoate hydroxylase in Fig. 1, the in state to the out state conversion involves two collective variables with a torsion around the dihedral angles. These dihedral angle changes take place around the blue (X) and red (Y) axes shown in Fig. 18, where the specifications related to the applied forces are explained.

Fig. 18 Definitions of dihedral angles in the umbrella sampling of the conformational changes of FAD in scKMO. The collective variables involve two dihedral angle rotations around the axes X (blue) and Y (red). The yellow atoms were constrained (500 kcal mol⁻¹ Å⁻²) throughout all US simulations. The rotation around the X-axis was performed by applying a harmonic force (force constant = 0.25 kcal mol⁻¹ deg⁻²) to the dihedral angle formed by atoms a–b–c–d, while the rotation around the Y-axis was obtained by applying a harmonic force with the same constant to the dihedral angle formed by atoms b–c–d–e.

In the text, for simplicity, these collective variables are going to be referred to as X and Y, respectively. In addition, twelve atoms (shown in yellow in Fig. 18) were constrained during the US simulations to prevent structural relaxation to the starting configuration and drive the isoalloxazine ring system from the in state to the out state.

Since the generation of FAD in the out conformation involves a two-dimensional dynamical change, three US schemes with different levels of accuracy were used consecutively. First, approximate US simulations were carried out for a fast initial exploration around the X-axis in steps of 5° windows. For each window, MD simulations were carried out for 1 ns with a 1 fs time-step. In this one-dimensional screening, the dihedral angle around the X-axis was scanned while the dihedral angle around the Y-axis was fixed with the harmonic force and treated as a parameter. These simulations were repeated at various fixed values of the dihedral angle Y in steps of 5°. These initial simulations will be referred to as the ‘fast’ US(X) simulations. They were carried out in an NPT ensemble in which the temperature was increased from 0 K to 310 K using the Langevin thermostat. Those fast US(X) simulations were only used to obtain a qualitative understanding of the free-energy landscape related to the conformational transformation of the in state to the out state, and to give initial coordinates for more accurate scanning of the free-energy surface.

Second, selected ‘fast’ US(X) simulations were repeated using their final coordinates and with a refined increment of 1° (‘slow’ US(X)). For each window, first an equilibration run was carried out for 1 ns with a 1 fs time-step (with the same simulation parameters as described for the fast US(X)), then a production run was carried out for 7 ns with a 1 fs time-step (with the same simulation parameters as described for the fast US(X), but at constant temperature, 310 K).

Finally, in the third US scheme (‘slow’ US(Y)), which involved rotating around the dihedral angle Y and the dihedral angle X as the parameter, the coordinates obtained from the selected consecutive ‘slow’ US(X) simulations were used as the initial coordinates. In this scheme, the windows were at an interval of 1°; first, an equilibration run was carried out for 1 ns with a 1 fs time-step, followed by a production run for 7 ns with a 1 fs time-step, as above. Each ‘slow’ US(Y) was repeated for incremental changes of the parameter X by 1°. Since the windows of ‘slow’ US(Y) are very close in the two-dimensional X–Y conformational space, the obtained data were combined, and the local free-energy surface related to this X–Y space was analyzed in two dimensions. This procedure allowed the identification and localization of both minima and the saddle point of the in to out transition. Although some overlap of windows of initial ‘fast’ simulations was already obtained, these simulations did not allow extraction of quantitative information about the free-energy surface. Therefore, they merely served as qualitative starting points for consecutive ‘slow’ simulations. On the other hand, the overlap of windows in the ‘slow’ simulations was excellent throughout. Due to the large number of windows, we are discussing only representative ones here. For these selected windows (which refer to minima and saddle points), the overlapping data are displayed for each protein complex (in Fig. S7–S13†). Free-energy analyses were carried out using the Weighted Histogram Analysis Method²⁵ (WHAM). One- and two-dimensional analyses were performed using the wham and wham-2d programs,⁴² respectively.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by The Scientific and Technological Research Council of TÜRKİYE (TÜBİTAK application number: 1059B192202527). We thank the Max Planck Society for the Advancement of Science for financial support. This work is part of the Action CA21162 ‘Computational Redesign of Enzymes (COZYME)’ of COST (European Cooperation in Science and Technology). The numerical calculations reported in this paper were partially performed at TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources).

References

1. G. J. Guillemin, K. M. Cullen, C. K. Lim, G. A. Smythe, B. Garner, V. Kapoor, O. Takikawa and B. J. Brew, J. Neurosci., 2007, 27, 12884–12892.
2. J. R. Smith, J. F. Jamie and G. J. Guillemin, Drug Discovery Today, 2016, 21, 315–324.
3. P. Guidetti, R. E. Luthi-Carter, S. J. Augood and R. Schwarcz, Neurobiol. Dis., 2004, 17, 455–461.
4. D. Zwilling, S.-Y. Huang, K. V. Sathyasaikumar, F. M. Notarangelo, P. Guidetti, H.-Q. Wu, J. Lee, J. Truong, Y. Andrews-Zwilling and E. W. Hsieh, et al. , Cell, 2011, 145, 863–874.
5. T. Ogawa, W. Matson, M. Beal, R. Myers, E. Bird, P. Milbury and S. Saso, Neurology, 1992, 42, 1702–1702.
6. M. Perkins and T. Stone, Brain Res., 1982, 247, 184–187.
7. T. W. Stone, Trends Pharmacol. Sci., 2000, 21, 149–154.
8. C. Hilmas, E. F. Pereira, M. Alkondon, A. Rassoulpour, R. Schwarcz and E. X. Albuquerque, J. Neurosci., 2001, 21, 7463–7473.
9. C. Prescott, A. M. Weeks, K. J. Staley and K. M. Partin, Neurosci. Lett., 2006, 402, 108–112.
10. A. C. Foster, J. A. Kemp, P. D. Leeson, S. Grimwood, A. E. Donald, G. R. Marshall, T. Priestley, J. Smith and R. Carling, Mol. Pharmacol., 1992, 41, 914–922.
11. E. X. Albuquerque and R. Schwarcz, Biochem. Pharmacol., 2013, 85, 1027–1032.
12. R. S. Phillips, E. C. Iradukunda, T. Hughes and J. P. Bowen, Front. Mol. Biosci., 2019, 6, 3.
13. K. R. Crozier-Reabe, R. S. Phillips and G. R. Moran, Biochemistry, 2008, 47, 12420–12433.
14. Y. Özkılıç and N. Ş. Tüzün, J. Phys. Chem. A, 2019, 123, 3149–3159.
15. Y. Özkılıç and N. Ş. Tüzün, J. Phys. Chem. A, 2021, 125, 9459–9477.
16. H. T. Kim, B. K. Na, J. Chung, S. Kim, S. K. Kwon, H. Cha, J. Son, J. M. Cho and K. Y. Hwang, Cell Chem. Biol., 2018, 25, 426–438.
17. D. L. Gatti, B. A. Palfey, M. S. Lah, B. Entsch, V. Massey, D. P. Ballou and M. L. Ludwig, Science, 1994, 266, 110–114.
18. B. A. Palfey, G. R. Moran, B. Entsch, D. P. Ballou and V. Massey, Biochemistry, 1999, 38, 1153–1158.
19. H. A. Schreuder, A. Mattevi, G. Obmolova, K. H. Kalk, W. G. Hol, F. J. van der Bolt and W. J. van Berkel, Biochemistry, 1994, 33, 10161–10170.
20. K. Hirai, H. Kuroyanagi, Y. Tatebayashi, Y. Hayashi, K. Hirabayashi-Takahashi, K. Saito, S. Haga, T. Uemura and S. Izumi, J. Biochem., 2010, 148, 639–650.
21. J. P. Hutchinson, P. Rowland, M. R. Taylor, E. M. Christodoulou, C. Haslam, C. I. Hobbs, D. S. Holmes, P. Homes, J. Liddle and D. J. Mole, et al. , Nat. Commun., 2017, 8, 15827.
22. S. Zhang, M. Sakuma, G. S. Deora, C. W. Levy, A. Klausing, C. Breda, K. D. Read, C. D. Edlin, B. P. Ross and M. Wright Muelas, et al. , Commun. Biol., 2019, 2, 271.
23. M. Amaral, C. Levy, D. J. Heyes, P. Lafite, T. F. Outeiro, F. Giorgini, D. Leys and N. S. Scrutton, Nature, 2013, 496, 382–385.
24. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187–199.
25. S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem., 1992, 13, 1011–1021.
26. C. R. Martinez and B. L. Iverson, Chem. Sci., 2012, 3, 2191–2201.
27. A. Waterhouse, M. Bertoni, S. Bienert, G. Studer, G. Tauriello, R. Gumienny, F. T. Heer, T. A. P. de Beer, C. Rempfer and L. Bordoli, et al. , Nucleic Acids Res., 2018, 46, W296–W303.
28. L. Schrödinger and W. DeLano, PyMOL, 2020, https://www.pymol.org/pymol.
29. D. A. Case, K. Belfon, I. Y. Ben-Shalom, S. R. Brozell, D. S. Cerutti, T. E. Cheatham III, V. W. D. Cruzeiro, T. A. Darden, R. E. Duke, G. Giambasu, et al., AMBER 2020, University of California, San Francisco, 2020.
30. Schrödinger Release 2020-2: Maestro, Schrödinger, LLC, New York, NY, 2020.
31. J. C. Gordon, J. B. Myers, T. Folta, V. Shoja, L. S. Heath and A. Onufriev, Nucleic Acids Res., 2005, 33, W368–W371.
32. C. Tian, K. Kasavajhala, K. A. Belfon, L. Raguette, H. Huang, A. N. Migues, J. Bickel, Y. Wang, J. Pincay and Q. Wu, et al. , J. Chem. Theory Comput., 2020, 16, 528–552.
33. J. Wang, W. Wang, P. A. Kollman and D. A. Case, J. Mol. Graphics Modell., 2006, 25, 247–260.
34. C. I. Bayly, P. Cieplak, W. Cornell and P. A. Kollman, J. Phys. Chem., 1993, 97, 10269–10280.
35. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, et al., Gaussian 16, Revision A.03, Gaussian, Inc., Wallingford, CT, 2016.
36. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
37. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
38. U. C. Singh and P. A. Kollman, J. Comput. Chem., 1984, 5, 129–145.
39. B. H. Besler, K. M. Merz Jr. and P. A. Kollman, J. Comput. Chem., 1990, 11, 431–439.
40. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
41. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
42. A. Grossfield, An implementation of WHAM: the Weighted Histogram Analysis Method Version 2.0. 11, https://membrane.urmc.rochester.edu/wordpress/?page_id=126.

---

Footnote

† Electronic supplementary information (ESI) available: A pdf file containing Ramachandran plots, equilibration plots, US overlapping examples, additional contour maps, details of relative free-energy calculations, structural details, and analysis of the number of water molecules around the isoalloxazine, and a pdb file for the out conformational state of the effector-inhibitor bound complex. See DOI: https://doi.org/10.1039/d4ob00862f

---