Prediction of Substrate Affinity by Molecular Dynamics Simulations and WaterSwap

(S)-selective pose with a hydride attack angle of 90.5° and distance of 3.68 Å is observed, which has a lower score than the (R)-selective poses.

In the induced-fit docking of C25G/I67T with 17a, all found poses are in the flipped orientation which will give rise to (S)-17b. This result is in quite good agreement with the experimental value of 76% (S)-selectivity. Notably, rear-rangements of H172 and H175 are seen, by which hydrogen bonds of both resi-dues to the carbonyl can occur.

The docking results show 17a in C25G with the flipped pose and C25D with the normal pose and therefore give an explanation for the facial selectivity of these two variants. Overall, docking revealed that the larger hydrophobic pocket opened through C25G and I69T exchange allows sufficient space and stabilizing van-der Waals interactions to host smaller substituents like methyl, ethyl or methoxy groups, allowing flipped binding poses similar to the one of 17a in the wt to occur. When C25D is present, a polar pocket is formed, disfavouring such interactions. Substrates are forced to flip into normal binding poses, thereby presenting the other face of the carbon-carbon bond to the hydride. Neverthe-less, a bit more space is generated in the C25D/I67T variant, which allows the normal poses to approach closer to the hydride in an optimal angle (Figure 37B).

3.5.2 Prediction of Substrate Affinity by Molecular Dynamics

understanding of the molecular interactions between the TsER proteins and substrates, enabling analysis of the effects of mutation on protein dynamics and substrate binding.

The results of the experimental studies should be examined with the help of fundamental techniques in theoretical chemistry. Theoretical studies, such as those employing molecular dynamics (MDs) simulations, can provide im-portant insights into mechanistic details that may not be possible via experi-mental means. Relatively few theoretical studies have been performed on the OYE family.^[130–132] An explanation for the different reactivities of TsER wt and the double variants C25D/I67T and C25G/I67T towards especially 17a is ex-pected from the MD simulations.

3.5.2.2 Molecular Dynamics Simulations with TsER Wildtype

The IFD pose of the flipped 3-methylcyclohex-2-en-1-one (17a, MCH) in the TsER wt structure based on pdb 3HGJ with FMNH² described in chapter 3.5.1 was used as basis for simulations (see Figure 37A). PROPKA^[276] was used to predict the protonation state of titratable residues in the protein at working pH 7.4. This method predicted both histidine residues in positions 172 and 175 to be neutral, which was considered in the MD simulations.

During the early stages of MD simulations some small restraints were applied to all C^α atoms, see chapter 6.10 for details. Analysis of the trajectories showed that the hydrogen bonds between the substrate and H172 and H175 were bro-ken after 2 ns. The distance between the carbonyl oxygen and H172 stayed in range of 2-4 Å between 5-20 ns before moving to much larger distances (Figure 38E). The substrate remained above the flavin in the flipped orientation, which is illustrated by the dominant cluster from the 100 ns MD simulations (see Fig-ure 38A).

Figure 38. MD simulations of TsER wildtype with reduced FMNH2 and 3-methyl cyclohexenone (17a). A) Active site of TsER wt in the dominant cluster from MD simulations with 17a (green lines) in complex. Water molecules within 4 Å are shown in stick form. B) Root mean square deviation (RMSD) of all protein over 100 ns. C) RMSD of 17a in TsER wt over 100 ns. D) The distance between N5H of FMNH2 and C3 of 17a over the course of 100 ns MD simulations. E) The distance between NH of H175 (black) and NH of H172 (red) to O of 17a over the course of

The high fluctuation of 17a can be seen in the root-mean-square deviation (RMSD) of atomic positions (Figure 38C). The distance between the substrate C3 with the flavin N5 hydride was monitored over the course of the simulations, giving an average value of 4.55 Å (Figure 38D) for this distance.

LONSDALE et al. discussed the role of Y169 as proton donor in OYE YqjM.^[132]

They observed that the proton of Y169 points away from the substrate during the majority of the simulation, even though they suggest that Y169 is more like-ly to protonate the substrate directlike-ly, rather than via a bridging water mole-cule.^[132] In TsER the corresponding residue is Y177. In the clustered pose the OH group of Y177 points towards the C2 atom of 17a, with a distance of 3.4 Å, indicating that protonation would be possible. Additionally, the OH group of Y177 interacts with nearby water molecules which could lead to protonation over a water network, in agreement with the study of YqjM.^[132]

Experimentally, TsER wt is found to reduce 17a with just a 1% yield; however, these simulations indicate that reactive conformations are accessible. Structures from these MD simulations were used to initiate binding energy calculations using the WaterSwap method, (described in chapter 3.5.2.5), to predict absolute binding free energies of 17a to TsER wt.

As mentioned in chapter 3.5.1, why the conversion of ketoisophorone 22a by most OYEs is possible but not isophorone 44a is still unresolved. The RBD stud-ies gave an initial explanation by indicating different binding poses for both compounds (Figure 36). In the pose for 22a, the methyl group is in α-position for hydrogenation, whereas for 44a it is in β-position. These docking poses were taken as a starting point for MD simulations.

Figure 39 shows the dominant cluster from 70 ns MD simulations with 22a (KIP). The enzyme substrate complex stayed stable over the whole simulation, with the methyl group in α-position for hydrogenation. According to the analy-sis of the distances between C2 of 22a and the N5 hydride of FMN as well as

between the carbonyl oxygens of 22a and H172/H175, the substrates stayed in the active side for at least 60 ns.

Figure 39. MD simulations of TsER wildtype with reduced FMNH2 and ketoisophorone (KIP, 22a). A) Active site of TsER wt in the dominant cluster from MD simulations with 22a (green lines) in complex. Water molecules within 4 Å are shown in stick form. B) Root mean square deviation (RMSD) of protein over 70 ns. C) The distance between N5H of FMNH² and C3 of 22a over the course of 70 ns MD simulations. D) The distance between NH of H175 (black) and NH of H172 (red) to O1 of 22a over the course of 70 ns MD simulations. E) The distance between NH of H175 (black) and NH of H172 (red) to O2 of 22a over the course of 70 ns MD simulations.

In the starting conformation from docking, 22a is orientated in the flipped mode, whereas in the dominant cluster from MD simulations it is oriented in the normal pose. Experimentally TsER wt produces 90% (R)-22b, which means that the hydride must attack C2 from the bottom (distance 4.7 Å) and C3 gets protonated from the top in the normal binding pose. In Figure 39A, it can be

seen that Y177 is too far from the substrate to interact directly and the protona-tion can only occur via a nearby water molecule.

Figure 40 shows the dominant cluster from 70 ns MD simulations with 44a (IPN). The distance between the substrate C3 with the flavin N5 hydride was monitored over the course of the simulations and is most often between 4-6 Å.

The initial docking pose showed 44a in the flipped orientation with C3 in β-positon for hydride transfer, whereas in the dominant cluster from MD simula-tions it is oriented in the normal pose and not stacking above FMN. The posi-tioning of 44a makes it impossible for hydride transfer with a distance of 6.5 Å and a hydride attack angle of 51.7°.

Figure 40. MD simulations of TsER wildtype with reduced FMNH² and isophorone (44a). A) Active site of TsER wt in the dominant cluster from MD simulations with 44a (green lines) in complex. Water molecules within 4 Å are shown in stick form. B) Root mean square deviation (RMSD) of protein over 70 ns. C) The distance between N5H of FMNH² and C3 of 44a over the course of 70 ns MD simulations.

Comparison of the dominant structures from simulations of both 22a and 44a gives an explanation of the experimental results. 22a is consistently in a reactive conformation, allowing quantitative conversion by TsER wt, whereas 44a can-not adopt a reactive conformation and is can-not converted at all.

3.5.2.3 Molecular Dynamics Simulations with TsER C25G/I67T

The IFD of 3-methylcyclohex-2-en-1-one (17a, MCH) in the prepared TsER C25G/I67T structure based on the wildtype crystal structure (pdb 3HGJ) with FMNH2 was used as basis for simulations.

Figure 41 shows the dominant cluster from 70 ns MD simulations with 17a. The RMSD of the whole protein over 70 ns MD simulations is shown in Figure 41B.

The average value for the RMSD is 3 Å, indicating that there is quite a signifi-cant difference in atomic positions to the crystal structure. According to the analysis of the distances between C3 of 17a and the N5 hydride of FMN (Figure 41D) as well as between the carbonyl oxygens of 17a and H172/H175 (Figure 41E) the substrate leaves the active site after 20 ns.

The structure of the binding site from the dominant cluster of the simulations with TsER C25G/I67T and 17a confirmed that the substrate leaves the active site (Figure 41A). That a substrate is able to leave the active site of TsER is quite im-portant for the catalytic reaction, because otherwise it would not be possible to follow the bi-bi ping pong mechanism with NADPH. In contrast to an inhibitor, a converted substrate should not have an enormous binding affinity, otherwise the reaction will stop. FITZPATRICK et al. experimentally observed weak binding affinities of substrate cyclohexenone to YqjM.^[277]

Figure 41. MD simulations of TsER C25G/I67T with reduced FMNH² and 3-methyl cyclohex-enone (17a). A) Active site of TsER C25G/I67T in the dominant cluster from MD simulations with 17a (green lines) in complex. Water molecules within 4 Å of FMN are shown in stick form.

B) Root mean square deviation (RMSD) of protein over 70 ns. C) RMSD of 17a in TsER C25G/I67T over 70 ns. D) The distance between N5H of FMNH² and C3 of 17a over the course of 70 ns MD simulations. E) The distance between NH of H175 (black) and NH of H171 (red) to O of 17a over the course of 70 ns MD simulations.

3.5.2.4 Molecular Dynamics Simulations with TsER C25D/I67T

The IFD pose of 3-methylcyclohex-2-en-1-one (17a) in the prepared TsER C25D/I67T structure based on the wildtype crystal structure (pdb 3HGJ) with FMNH2 was used as basis for simulations. The variant model based on the wildtype structure was used due to the higher resolution compared to the crys-tal structure of C25D/I67T (pdb 5NUX).

The first round of MD simulations was performed with the same parameters as described for the wildtype and C25G/I67T variant with neutral H172 and H175.

In these simulations the protein-substrate complex was not stable and the sub-strate 17a left the active site within the first nanosecond of simulation. Figure 42C illustrates the dominate cluster from the 100 ns MD simulations.

For this reason, alternative protonation states of H172 and H175 were tested, as in the study of LONSDALE et al. for YqjM.^[132] The MD simulations were repeated with protonated H172 and H175. The dominant cluster from the 100 ns MD simulations is shown in Figure 42D, which illustrates a stable complex of FMN with 17a, but shows a flip of the cofactor-substrate complex out of the binding pocket. Further inspection of the non-covalently bound cofactor FMN reveals that the hydrogen bond between the phosphate moiety and R319 is broken, al-lowing this flip to occur.

Hence, the MD simulations were repeated a third time with harmonic restraints of 20 kcal/mol · Ų from FMN to the backbone of R319 (Figure 42 grey, Figure 43). The first 10 ns of simulations showed a stable substrate-enzyme complex with the flavin correctly positioned in the binding site. After 10 ns, substrate 17a moved from its position above the flavin. For this reason the simulations were stopped after 30 ns.

Figure 42A shows the structural alignment of the enzyme backbone over the three different MD simulations, which did not indicate any extreme changes, like unfolding in the dimeric structure. Whereas Figure 42B clearly illustrates

the different binding modes of the flavin in the active site of TsER C25D/I67T depending on the protonation states of H172/H175.

Figure 42. Comparison of MD simulations of three different TsER C25D/I67T structures with reduced FMNH² and 3-methyl cyclohexenone (17a). A) Structural alignment of MD simulations of TsER C25D/I67T with neutral H172/H175 (purple), with H172/H175 protonated (cyan) and with H172/H175 protonated and distance restraint of FMN to backbone (grey). B) Illustration of FMN and H172/H175/Y177 in alignment. C) Active site with 17a and neutral H172/H175, D) with 17a and protonated H172/H175, E) with 17a and protonated H172/H175 and backbone restraints to FMN.

Figure 43. MD simulations of TsER C25D/I67T with reduced FMNH2, restraints of FMN to backbone, protonated H172/H175 and 3-methyl cyclohexenone (17a). A) Active site of TsER C25D/I67T after 10 ns MD simulations with 17a (green lines) in complex. Water molecules with-in 4 Å are shown with-in stick form. B) Root mean square deviation (RMSD) of protewith-in over 10 ns. C) RMSD of 17a in TsER C25D/I67T over 10 ns. D) The distance between N5H of FMNH2 and C3 of 17a over the course of 10 ns MD simulations. E) The distance between NH of H175 (black) and NH of H171 (red) to O of 17a over the course of 10 ns MD simulations.

Due to this observation constant pH MD simulations^[278] for C25D/I67T with 17a in complex were performed. With this method pK^a values for the titratable side chains of residues D25, H172, H175 and Y177 were calculated, to gain more in-sight into their different protonation states at the working pH of 7.4.

Table 7. Predicted pK^a values for titratable residues D25, H172, H175 and Y177 from constant pH MD simulations at pH 7.4.

residue Offset^(a) predicted pK^a

Fraction protonated^(b)

Transitions^(c)

D25 2.479 9.879 0.997 1

H172 -4.108 3.292 0.000 1

H175 -2.754 4.646 0.002 23

Y177 4.848 12.248 1.000 2

a) is the difference between the predicted pKa and the system pH, b)is the fraction of time the residue spends protonat-ed, c) gives the number of accepted protonation state transitions over 100 fs.

The results of populations of every state for every titratable residue show that D25 is deprotonated for 0.33% of the time and only syn-protonated (>99% of the time) on each oxygen. The H172 residue is deprotonated for 0.0078% of the time, syn-protonated for 0.24% and anti-protonated for 99.7% on each nitrogen.

Residue H175 is deprotonated for 0.17% of the time, syn-protonated for 74%

and anti-protonated for 25.7% on each nitrogen. Y177 stayed always protonated during the simulations.

These results indicated the different protonation states of at least three residues, which are important for the catalytic activity. To achieve stable molecular dy-namics structures for TsER C25D/I67T all of these protonation states have to be considered. This would be a highly computational intensive calculation which needs further investigation.

3.5.2.5 Calculations of Binding Free Energies with WaterSwap and MM-PBSA

The WaterSwap method works by constructing a reaction coordinate that swaps the ligand bound to the protein with an equivalent volume of water mol-ecules. The affect is to move the substrate from being bound to the protein, to

being free in solution, while simultaneously transferring an equivalent volume of water from being free in solution to being bound to the protein.^[279] The ad-vantage of WaterSwap is that it allows calculation of absolute binding free en-ergies.^[280] This simplifies investigation of the effect of protein mutation on bind-ing compared to relative bindbind-ing free energy methods, as differences in abso-lute binding free energies capture changes in binding mode of the ligand upon protein mutation.^[280]

The input files for a WaterSwap calculation are the AMBER format coordinate and topology files from MD simulations that represent the solvated protein-ligand complex.^[279] WaterSwap^[279] absolute binding free energy calculations of the ligand 3-methyl cyclohex-2-en-1-one (17a) bound to TsER wt followed the MD simulations, using the structure after 100 ns from the dynamics simula-tions.

MM-PBSA is another post-processing end-state method to calculate free ener-gies of molecules in solution from MD simulations.^[281] This end-state free ener-gy method reduces computational cost by eliminating the need for simulating intermediate states. Modelling the solvent implicitly further reduces the compu-tational cost by eliminating the noise caused by explicit solvent molecules.^[282]

To experimentally determine binding free energies of substrates to OYEs is a challenging task. Most often steady state kinetics under anaerobic conditions with purified enzyme are performed.^[95,283] Using WaterSwap as in-silico pre-screening for substrate libraries in OYEs would be a great breakthrough.

The WaterSwap calculation for the dominant cluster of TsER wt with 17a (Figure 38) resulted in a calculated binding free energy of 12 ± 2.266 kcal/mol.

Unfortunatly, there is no comparable value from an experiment, therefore the second in-silico binding free energy method MM-PBSA was used. With this method it was not possible to get any analysable data. Table 8 shows that all

other performed calculations with the variants C25D/I67T and C25G/I67T and also with the substrates 22a and 44a failed.

Table 8. Performed calculations of binding free energies with WaterSwap and MM-PBSA.

enzyme substrate WaterSwap MM-PBSA wildtype 17a 12 kcal/mol non analysable C25D/I67T 17a non analysable non analysable C25G/I67T 17a non analysable non analysable wildtype 22a non analysable non analysable C25D/I67T 22a non analysable non analysable C25G/I67T 22a non analysable non analysable wildtype 44a non analysable non analysable C25D/I67T 44a non analysable non analysable C25G/I67T 44a non analysable non analysable

A reason for this might be the unstable MD simulations with the variants and all tested substrates and also the presence of a non-covalently bound cofactor.

These calculation methods are not yet optimized for such complex systems and need further investigations.