Radical Stabilization Energies for Enzyme Engineering – Tackling the Substrate Scope of the Radical Enzyme QueE

: Experimental assessment of the reaction mechanisms and profiles of radical enzymes can be severely challenging due to the reactive nature of the intermediates, and sensitivity of cofactors such as iron sulfur clusters. Here we present an enzyme-directed computational methodology for the assessment of thermodynamic reaction profiles and screening for radical stabilization energies (RSEs) for the assessment of catalytic turnovers in radical enzymes. We have applied this new screening method to the radical SAM enzyme CPH 4 synthase (QueE), following a detailed molecular dynamics (MD) analysis that clarifies the role of both specific enzyme residues and bound Mg 2+ , Ca 2+ or Na + . The MD simulations provided the basis for a statistical approach to sample different conformational outcomes. RSE calculation at the M06-2X/6-31+G* level of theory provided the most computationally cost-effective assessment of enzyme-based energies, facilitated by an initial triage using semi-empirical methods. The impact of intermolecular interactions on RSE was clearly established and application to the assessment of potential alternative substrates (focusing on radical clock type rearrangements) proposes a selection of carbon-substituted analogues that would react to afford cyclopropylcarbinyl radical intermediates, as candidates for catalytic turnover by QueE.


INTRODUCTION
Radical intermediates are extremely versatile for chemical functionalization and transformation reactions. Due to their high reactivity, radicals can facilitate these reactions with otherwise non-activated, unreactive substrates. However, this advantage comes at the cost that these highly reactive intermediates are particularly hard to control. This can lead to a multitude of possible unwanted side reactions and is one of the reasons why radical chemistry is predominantly found in industry for processes where either such side reactions are desirable or side reactions can be controlled or eliminated, for example in the downstream processing of crude oil (cracking) or in polymerization chemistry.
In nature, radical reactions play an important role in enzyme catalysis. The radical SAM enzyme superfamily is one group of enzymes that is capable of exploiting the potential of radical reactions in a very controlled way. These enzymes are able to initiate radical formation and direct their reaction by both preventing side reactions and facilitating the desired reaction simultaneously. This result of millions of years of evolution harnesses key similarities in catalytic mechanism, and yet results in a broad chemical reaction space (see reviews by Broderick et al. 1 , Dowling et al. 2 , and Jaeger and Croft 3 for a critical summary). Radical SAM enzymes catalyze reactions that include C-C bond formations, 4-6 decarboxylation reactions, 7 functional group migrations (1,2-shifts), [8][9] sulfur insertions, [10][11][12][13] methylations , [14][15] and more complex radical rearrangement mechanisms, [16][17][18][19] with these radical-mediated transformations involved in a multitude of biochemical synthesis routes that include compounds with antibiotic and antiviral activity. As such, it would be highly beneficial to gain access to and adapt these biotransformations for their use in industrial biotechnological applications, facilitating sustainable routes towards fine chemicals, pharmaceuticals, or bulk chemicals that would be highly challenging to synthesize by alternative methods. 3 The key commonality for the catalysis of radical SAM enzymes is that they use S-adenosylmethione (SAM) either as cofactor or co-substrate. SAM is bound to a central Fe4S4 iron sulfur cluster responsible for initiating the redox reaction. The cluster is embedded via binding to cysteine residues of a conserved CX3CXφC motif (with φ representing a conserved aromatic) and transfers an electron upon reduction to the SAM molecule which subsequently cleaves to afford the 5'-adenosyl radical (Ado•) and methionine (Met), which remains bound to the cluster. The Ado• radical represents the first reactive intermediate that then abstracts a hydrogen from a bound substrate to initiate catalysis.
Generally, it is believed that the control of radical intermediates in radical SAM enzymes is based on perfect positioning of the substrate towards the cluster bound SAM and the stabilization of the radicals. The argument is that more stable radicals are less likely to undergo unwanted side reactions, for example with the enzyme itself. However we recently showed with the example of the radical rearrangement in 7-carboxy-7-deazaguanine (CDG) synthase (QueE), 20 that this assumption is not always true. In the case of QueE, the unmodified substrate radical is so stable that the energy barrier for the subsequent radical rearrangement is too high for efficient catalysis, unless the radical is structurally or electronically perturbed. 21 The reaction catalyzed by QueE represents a central step in queosine synthesis and facilitates the formation of the 7deazapurine scaffold through a radical-mediated ring contraction. 22 Scheme 1 shows the radical rearrangement during the catalysis of QueE propagating through an azacycloproylcarbinyl intermediate followed by NH2 elimination. This mechanism has been first described as one of two potential pathways by Drennan and coworkers, 20 and was confirmed computationally by Zhu and Liu 23 and us. 21 The rearrangement of CPH4 proceeds after initial hydrogen abstraction from the C6 position through a cyclic azacyclopropylcarbinyl intermediate (3), before hydrogen reabstraction from AdoH and deamination to form the final product (6). Intermediate 3 and its analogues also represent a class of structures known as radical clocks, 24 which undergo quick unimolecular rearrangements and have been described extensively by experiment and theory. [25][26][27][28][29][30][31][32] We recently demonstrated that it is necessary to hold the substrate radical 2 in an energetically unfavorable configuration to overcome the rearrangement barrier for the ring conversion. This conformation is achieved by binding the substrate already in this conformation (which represents for the substrate a local energy minimum only slightly higher in energy than the unbound conformer) and hold this conformation after hydrogen abstraction through electrostatic fixation by a Mg 2+ ion in the active site of the enzyme. Without this constraint, the radical would fall into a very stable, planar conformation and the energy barrier for the rearrangement would be too high for efficient reaction turnover. In other words, the substrate radical needs to be destabilized in order to facilitate the reaction.
This observation was also demonstrated by evaluating the stabilization of the radicals represented by radical stabilization energies (RSEs). By calculating the energy of a formal hydrogen transfer reaction (as described in the Experimental Section), RSEs inform on the reactivity of radicals in comparison to a given reference species. Thus, the higher the radical stabilization energy, the less energized and reactive the corresponding radicals are, also diminishing their role in unwanted side reactions. The concept of calculating RSEs has been used extensively for evaluating accurate absolute bond dissociation energies (BDEs) for a broad range of radicals including many systems related to peptides and enzyme catalysis. [33][34][35][36][37][38][39][40] More recently Hioe and Zipse also demonstrated the usefulness of calculating radical stabilization energies for radical SAM enzymes. 40 They have calculated the thermodynamic reaction profiles for selected radical SAM enzymes on the basis of high-level RSE calculations in gas phase and demonstrated that enzymes using SAM as cofactor generally combine an initial exothermic hydrogen abstraction step with a subsequent endothermic step, while enzymes using SAM as co-substrate perform significantly exothermic Habstraction reactions. What was not incorporated in this study was the effect on the radical stabilities and on the reaction profiles of binding to the enzyme and additional cofactors. As shown above, QueE is one such example where the intermolecular environment has a significant influence on these observables.

Scheme 1.
Proposed ring contraction mechanism in the radical SAM enzyme QueE.
In this study, a combination of molecular dynamics (MD) simulations and RSE calculations is presented in order to investigate the effect of intermolecular binding on RSE estimation. It is demonstrated that a quick and affordable strategy can highlight if the enzyme environment facilitates significant changes to the reaction profile in comparison to the reaction assessed in gas phase or solution. Further, this information can subsequently be used to screen alternative substrate candidates for their potential for conversion, by comparing their reaction profiles to the native reaction. Potential enzyme mutation candidates can also be screened that may stabilize radical intermediates differently to the wild type enzyme. Thus, the combination of MD sampling and RSE calculations has the potential be used as an initial tool for in silico enzyme design of radical SAM enzymes, bringing these enzymes one step closer to efficient biotechnological applications.

RESULTS AND DISCUSSION
MD SIMULATIONS As previously indicated, QueE is found to facilitate catalysis only effectively in presence of Mg 2+ . 20,22 To investigate the effects on substrate binding and conformation in the presence of different ions in the enzyme, a series of molecular dynamics simulations have been carried out. Starting from the crystal structures of QueE 20 in complex with SAM, the substrate CPH4 (1) and Na + (pdb code 4NJH) or Mg 2+ (pdb code 4NJI) in the active site, the corresponding ions were either left unaltered, switched between the two structures, or substituted by Ca 2+ . All simulations were run for at least 100 ns and maximal ~1 µs resulting in a total simulation time of approximately 6 μs. During the timeframe of all simulations, binding of SAM to the iron-sulfur cluster and the overall behavior of the dimeric enzyme showed high stability and low fluctuations. In a small number of simulations, conformational changes for SAM could be observed after several hundred nanoseconds, which was followed by the movement of a loop involved in SAM binding (see Figure S20 in the Supporting Information). Substrate unbinding, on the other hand, was observed in all cases after varying simulation time.
All simulations could reproducibly demonstrate significant differences for different ions in the active site. Based on the results of our previous study, which highlighted the need to bring the substrate into a bent conformation for optimal catalysis, the initial simulation analysis focused on analyzing this conformation and the positioning of the substrate for the initial hydrogen abstraction step.
Only simulations with Mg 2+ in the active site kept the substrate in the bent conformation for significant simulation times, with the average retention time in the preferred conformation of 150 ns. The range covered binding of only 10 ns in one case and over 226 ns in other. Such long retention times were rarely seen for simulations containing Ca 2+ and even less frequently for those containing Na + , demonstrated by more flexible binding and much shorter retention times of the correct binding mode. The observed average retention times for these two ions were 20 ns for Ca 2+ (0-57 ns) and 10 ns (0-26 ns) for Na + . It should be noted that these averages should not be used as a truly quantitative measure but as a trend. To present a statistically rigorous binding time, a much larger number of simulations and starting configurations would be needed. Figure 1 depicts the essential binding analysis of substrate binding for three selected simulations over 200 ns each (more detailed analysis for all simulations can be found in Section 1.4 of the Supporting Information). All three simulations started from the crystal structure 4NJI with the substrate bound in the reactive conformation and the analysis started after careful and extensive restrained equilibration for 9 ns. The RMSD of all simulations remains low (Figure 1 left) and the SAM molecule remains bound to the iron sulfur cluster interacting with the unique iron through one carboxylate oxygen and the nitrogen of the methionine moiety. The carbon-carbon distance, essential for the hydrogen transfer between the C5' carbon of SAM and the C6 carbon of CPH4, already depicts significant variation and differences. Without the need to calculate the transition state structure and energy for the hydrogen abstraction and on basis of structural sampling with SAM and not the cleaved Ado• radical, this C-C distance can already provide information on the relative likeliness for this initial hydrogen abstraction.
While for the Mg 2+ simulation the substrate stays in adequate position for almost the entire 200 ns in both active sites of the dimeric enzyme, both substrates leave the correct position in the Ca 2+ simulation within a few nanoseconds and even earlier during the simulation with Na + . In most cases, leaving the initial position is directly coupled to changing into the more planar conformer of CPH4 (also depicted in Figure 1). Only in some cases is the substrate seen leaving the pocket together with the tightly bound Mg 2+ before the substrate changes its conformation. In contrast, in simulations with Ca 2+ and Na + ions left the binding pocket without the substrate.
A more detailed look at substrate binding with different ions reveals further key differences ( Figure 2). In all simulations, the substrate is fixed in the active site by hydrogen bonds between its carboxylate oxygens and the side chains of Arg27 and Thr90. Additional frequently found hydrogen bonds, common for all simulations, involve the side chains of residues Glu13 and Gln15 and the backbone of Gly14. The main variations between the simulations arise from the ion coordination within the active site.
Mg 2+ coordinates throughout the whole simulation to one carboxylate and carbonyl oxygen, fixing the substrate in the bent conformation ( Figure 2a). The ion is further coordinated by 4 water molecules but with no direct interaction observed with residues in the active site. Coordination to Asp50 and Thr51 is mediated through first shell water molecules. This is in slight disagreement with the crystal structure 4NJI, which suggests a direct coordination between Mg 2+ and Thr51, but with a slightly long interatomic distance of 3.1 Å. 20 Further, the two monomers in the crystal structure dimer also depict slight differences in their electron density in their active sites, attributed to differing positions of water. Interestingly, an additional simulation with the proposed intermediate 5 and Mg 2+ in the active site (Figure 2d) resulted in a coordination chemistry resembling the crystal structure, with direct coordination between the ion and Thr51, 3 water molecules and the substrate. Further, the hydrogen bonding interaction between the carboxylate of the C-terminal Pro209 residue and the substrate, described in the crystal structure, is indicated as being much stronger in this intermediate-containing simulation. This contrasts with simulations containing bound substrate, where this interaction is loose and the C-terminus shows flexibility during the simulations. Distinct from Mg 2+ , Ca 2+ coordinates the carbonyl oxygen of CPH4 and Asp50 directly, while the contact to the carboxylate of the substrate is water-mediated (Figure 2b). The frequency of hydrogen bonding to Glu13 and Gln15 is also significantly reduced. This suggests additionally that, during the period where the substrate adopts the bent conformation, it is fixed less tightly in this conformation. This is certainly also the case in the very short time of correct binding observed with Na + . The ion either leaves the pocket quite rapidly or is coordinated to Asp50 and the carboxylate of the substrate. Thus, the substrate flips rapidly into the planar conformation and also migrates from the optimal position for hydrogen abstraction.
In summary, clear differences for the binding of the substrate and the different ions in the active site are seen by MD simulation analysis. Mg 2+ seems to place the substrate much better position for abstraction and in the correct conformation, relative to what is seen for the other ions. However, the question as to what this effectively means for the radical rearrangement remains open.

RSE ASSESMENT
The calculation of radical stabilization energies (RSEs) relating to different substrates and intermediates can provide a clearer thermodynamic picture of radical reactions in radical SAM enzymes. These calculations alone, though, lack information about the influence of enzyme binding on the thermodynamic reaction profile. From the MD simulations and our previous work on a DFT model system21 we already know that intermolecular interactions seem to have an important influence on the conformational space available for the reaction and thus on the stability of specific intermediates.
This study uses MD simulations as a sampling method for the system incorporating the bound substrate to then subsequently calculate the radical stabilization energies from multiple snapshots of these simulations. Firstly, this is applied without any further optimization of the structures from the simulations in order to get a rapid insight into changes in radical stabilization energies upon enzyme binding. Since this lacks the structural relaxation of the substrate and, more importantly, the radical, absolute RSE values cannot be compared to accurate RSE values from high level DFT gas phase optimizations. However, when compared to the similar sampling in the gas phase it can provide direct information regarding how conformational restrictions in the active site might influence the stability of the intermediate generated during formation. In other words, it can signpost significant influences brought about by enzyme-induced binding and interactions.
In a second step, the protein environment is added to the calculations to look for further effects influencing radical stability. Additional residues can either explicitly be taken into account in the calculations or be represented by partial point charges. In that way, the impact of electrostatic effects in the active site in further influencing radical stability can be observed.
In a final step, the substrate and its radical form are optimized within a field of point charges representing the enzyme. This computationally more expensive QM/MM-type step confirms if the radical is able to substantially change its conformation once formed, and if this would have an effect of its stability and the energetic profile of the radical reaction.
An initial simulation where the substrate is bound correctly with Mg 2+ was investigated. From this simulationfive thousand snapshots along a simulation time of 500 ns were extracted. For comparison, the same number of snapshots have been generated from a 500 ns simulation of the substrate in the gas phase. Subsequently, RSEs have been calculated for the C6 hydrogen abstracted radical (2) of the substrate for all snapshots on basis of the closed shell molecular mechanics geometries with a number of methods and basis sets.
Three exchange functionals (and standard Hartree-Fock calculations) and six basis sets have been considered and the results are shown in Table 1. Double and triple-ζ basis with polarisation and diffuse functions have been cross-analysed with functionals of varying degrees of sophistication. As expected, the Hartree-Fock (HF) calculations performed purely as standard implementations of HF are known to overestimate energies in protein systems and fail to describe dispersion correctly, making it an inappropriate method for thermochemical calculations. The Minnesota suite of functionals have been well received as suitable candidates for studying kinetics, thermochemistry and non-covalent interactions, 41 hence why two (M06-2X and M11) have been included in the preliminary calculations. Here the M06-2X functional was designed for these type of calculations, and the M11 has been presented as an improvement 42 ). M06-2X was used for further calculations as it was a more dependable solution than M11, which had difficulty converging for this particular set of systems. All further calculations have also been corrected for dispersion by use of Grimme's D3 dispersion correction with correction values taken from literature. 43 As the focus of this study depends more on relative shifts rather than absolute values, the choice of functional is not as significant as the choice of basis set. As long as a sensible choice is made, relative energy shifts are similar (see also Table S26 in the Supporting Information). In summary, calculations at the M06-2X/6-31+G* level of theory presented best-balanced cost-accuracy relation of all DFT methods tested and all following DFT results are presented using this approach. RSE values have also been tested with semi-empirical (SE) methods to test their applicability for a quick RSE screening from MD data since the use of SE is significantly faster, with a speed up of around 1200-fold. A positive linear correlation between the calculated RSE values from both DFT (M06-2X/6-31+G*) and SE (PM6-D3) (r2 = 0.86Å) could be observed for the given example (see Figure S25 of the Supporting Information). This correlation allows the subsequent construction of a Gaussian distribution of the SE RSE's and then use of QM to calculate the RSEs of higher accuracy on a subset of selected frames near the Gaussian peak position (see Table S28 of the Supporting Information). This delivers RSE values of high accuracy at reduced calculation time.
The RSE values have then been calculated for the Mg 2+bound substrate with and without optimization in the electrostatic field of the protein and have been compared to the initially calculated vertical RSE values. As can be seen from the data presented in Table 2, the vertical radical stabilization of the substrate radical bound correctly in the active site of QueE together with Mg 2+ drops significantly by 34.0 kJ mol -1 in comparison to gas phase sampling only considering the substrate/radical itself in the conformation retrieved from the MD simulations directly. Thus, it can clearly be seen that the main contribution to the change in stability originates from the conformational change. Further, the standard deviation of the energies appears to be higher (SD ±20.2) for the bound conformations. This is also expected, as these radical conformations are not close to a stable conformational minimum and thus slight geometric changes represent larger energetic changes on the potential energy surface. When the same procedure is repeated optimizing both the substrate and the radical in the electrostatic field of the protein, represented by atomic partial charges, the RSE values drop to more negative values. The shift in RSE upon enzyme binding even increases to -70.7 kJ mol -1 and while the standard deviation for the unconstrained gas phase system drops to a very low value of ±2.3, it remains high (±23.3) for the bound system. This once more indicates a radical in an uncomfortable conformation far from the preferred optimum and thus with large energy changes upon small structural changes.
When including the electrostatic field without further optimization (and comparison to the gas phase system) it can be seen that the electrostatic field seems to only have a small effect on radical stabilization for this example. Comparing the results after QM/MM optimization between DFT and SE calculations moreover shows that optimizations at the SE level are not reliable for calculating relative shifts of RSE values for constrained molecules (see Table S27 in the Supporting Information). Thus, the semi-empirical sampling can only be suggested for an initial sampling of MD trajectories as described before.
When comparing simulations of the substrate bound in the active site together with different cations, the consequences of the structural and dynamic differences seen in the MD simulations on the radical stability of the substrate in the enzyme are demonstrated very clearly. The average stabilization in cases where the substrate is bound in the correct conformation is shown to be significantly lower compared to unreactive binding of the substrate. As can be seen from Table 3 in a simulation where unbinding can be observed during the simulation with Ca 2+ in the binding pocket, the vRSE value drops by over 25 kJ mol -1 from -33.9 (bound) to -59.1 kJ mol -1 (unbound). The effect of unbinding on the radical stability of the substrate (and thus on the thermodynamic reaction profile) can also be monitored in form of QM post-processing of the underlying MD simulation. In this way, changes of this central feature can be monitored in quasi-real (simulation) time without other analysis needed. Figure 3 demonstrates this for the Ca 2+ simulation. The stability of the substrate radical shifts significantly after ~56 ns. Additional structural analyses confirm that this shift is correlated with movement of the substrate slightly out of the pocket, which changes the complexation to the cation and results in conformational change of the substrate. While the substrate is still anchored in the active site by the strong salt bridge between Arg27 and the substrate's carboxylate, it is not binding in a reactive fashion anymore. As demonstrated in the MD discussion above, this behavior is much more likely when Mg 2+ is substituted by other cations in the active site.

POTENTIAL SUBSTRATE SCOPE OF QueE
In the view of the possibility to apply the method to calculate vRSE values for radical enzymatic reactions in the context of protein engineering, a set of alternative substrates have been tested for their potential to be converted through a similar reaction mechanism to the natural substrate QueE. A set of structures were selected that are also able to react via an analogous radical rearrangement through either an azacyclopropylcarbinyl radical or their carbon-and sulfursubstituted analogues. The structures were chosen to represent different substitutions next to the radical center that either stabilize or destabilize the radical by electron pushing or pulling effects, or might influence substrate binding through the presence or absence of functional groups necessary for hydrogen bonding. The selected structures are listed in Figure 4. The gas phase radical stabilization energies of the alternative substrates were evaluated in an analogous fashion to our previous paper 21 at the M06-2x/6-311++G(3df,3p)//B3LYP/6-31+G(d) level of theory and have been compared to the corresponding radical rearrangement barriers at the same level of theory (and to even more accurate G3B3 data, where affordable). The results presented in Table 4 and Figure 5 clearly show two remarkable trends. Firstly, the radical rearrangement barriers for azacyclopropylcarbinyl appear to be significantly higher compared to their cyclopropylcarbinyl counterparts. Secondly, the trend that higher radical stabilization correlates with higher rearrangement barriers is confirmed for the selected examples. The trend is more prominent for rearrangements through the azacycloproylcarbinyl and outliers can mainly be attributed to heterocyclic structures that either contribute with additional spin delocalization within the ring or represent structural constraints due to the ring structure, which additionally hinder the rearrangement.
Subsequently, the alternative substrates examined were docked into the crystal structure of QueE using the docking program GLIDE from the Schroedinger suite of programs. For the resulting enzyme substrate complexes, the radical stabilization energies were then calculated as single point energies of the docked substrate conformations and included optimizing the substrate and substrate radical structures in the electrostatic field of the enzyme. The results were evaluated taking into account the following selection criteria for potential alternative substrates of QueE: 1) the correct positioning of the substrate carbon involved in the necessary first hydrogen abstraction step between the Ado• radical and the potential substrate to initiate the radical rearrangement; 2) the corresponding docking score as an indicator for substrate affinity to the pocket (This score is also compared to the lowest docking score of the corresponding substrate to see whether alternative binding conformations are more likely); 3) the radical stabilization energy of the substrate in order to evaluate if a high rearrangement barrier is to be expected; 4) the RSE after optimization in the electrostatic field. This last assessment adds information about the likelihood of a substrate to be bound in a preferred conformation for catalysis, but where the radical would undergo quick relaxation to form a stable inactive radical intermediate.
All these four criteria do not only give valuable information about whether a potential substrate might be a good candidate for catalysis, but also indicate potential ways to improve the candidacy by signposting additional mutations within the enzyme's active site. The full docking results are presented in the Table S29-30 in the Supporting Information and the most interesting results are briefly discussed below.
From the docking calculations, the three conformers of each molecule with the best docking score have been taken for further RSE analysis. Additional docking studies without including the Mg 2+ ion in the active site were also performed to investigate if some substrates might be suitable for transformation without the support of the ion. Applying this protocol to the natural substrate CPH4 (N8) (see Table 5) delivered reasonable docking poses with and without Mg 2+ in the active site. Subsequent RSE calculations including QM/MM geometry optimisations could confirm an increase in radical stabilization but the structures did not relax to a very stable planar conformation, making the rearrangement in principle possible.
Across the complete dataset, substrates docked in preferential position and conformation often drop into a more preferable and less reactive radical conformation with high radical stabilization upon optimization. This is indicates that these substrates are unlikely to undergo a catalytic rearrangement and that additional mutations of active site residues might be necessary to further maintain reactive conformations. Further, the observation that the single point calculations including the point charge (PC) field, but neglecting further optimizations, differ much more significantly from single point calculations without the PC field (compared to the results from vRSE single point calculations based on MD simulations) indicates that equilibration after docking is necessary for an adequate interpretation and reliable results, either through MD or geometry optimizing the structures. Figure 5. Calculated radical stabilization energies and rearrangement barriers (Ea) for QueE. The trend for the correlation between the two values is highlighted in yellow. The green area represents rearrangement barriers approximately suitable for catalysis. The red circles highlight the values for the natural substrate CPH4 in gas phase and in the model system from our previous study. 21 Unsurprisingly, substrates without carbonyl or carboxylic functional groups cannot be attached strongly to residues in the active site (in particular Arg27) and thus result in very low docking scores and variable docking orientations (Table S30 in the Supporting Information). Thus, substrates C1, C2, C4, and C5, for example, show relative low RSE values, but do not bind strongly enough or in the correct position for catalysis and are very unlikely to act as alternative substrates. For structures abele to facilitate supporting anchoring by Arg27, on the other hand, catalytic turnover might be possible (see ligands C3, C6, C7, and C8 in Table 5) with improved positioning and even without Mg 2+ support in some cases. Figure 6b depicts ligand C3 docked and subsequently optimized in its radical form in the pocket of QueE without Mg 2+ demonstrating adequate positioning and radical stabilization for enzymatic turnover. In general, the cyclopropylcarbinyl precursors, that do not need to overcome a similar high transition barrier to the aza-cyclopropylcarbinyl analogues, might function better without Mg 2+ due to better positioning. Also, ligand N3 (which includes the anchoring carbonyl but lacks the ring structure of the natural substrate) might undergo turnover without Mg 2+ as it's radical does not need to be fixed in a tightly bended conformation.
On the other hand, structures that are placed well for initial hydrogen abstraction but where significant stabilization occurs upon optimization effectively might act as inhibitors. Although they show the principle potential for rearrangement, they are trapped in a low energy minimum after hydrogen abstraction and thus are inhibited from further turnover. Examples showing this behavior for all conformers analyzed are ligands C4 and N4 (see Table S29-30 in the Supporting Information for details). Table 5. Selected RSE values with and without QM/MM geometry optimization in the point charge (PC) environment, based on docked structures of alternative substrates from Scheme 2. XPscore denotes the extra precision Glide docking score, and the C6-C5' distance indicates the crucial distance for the initial hydrogen abstraction reaction in the docked structure.
This first combined docking and RSE assessment gives new insights into how the radical enzyme QueE controls the central catalytic radical rearrangement. The protocol developed here offers a new and suitable approach to assess the thermodynamic reaction profile of natural and alternative substrates in radical SAM enzymes. Excitingly, it has good potential to serve as a pre-screening tool for alternative substrates in radical enzymes, as a first step towards enzyme engineering in this underexploited domain.

CONCLUSION
Radical stabilization plays an important role in the enzymatic catalysis of radical SAM enzymes. Often the way intermediate radicals are stabilized or destabilized within the active site determines the rate determining steps of the catalysis of these enzymes. In cases where the reaction mechanism is known, evaluating the thermodynamic reaction profiles for these enzymes offers an affordable gateway into computational predictions of substrate scope and initial steps into reengineering substrate scope, turnover, or promiscuity of these enzymes.
In this study on QueE we show that the combination of MD simulations and the evaluation of radical stabilities delivers deeper understanding in how the enzyme controls enzymatic turnover and what the crucial role of the central Mg 2+ ion is, namely that tight control of the reactive radical conformation Applying a combination of substrate docking and radical stabilization assessment on a set of alternative substrates delivered detailed information about potential alternative substrates and inhibitors. We could confirm that only ligands that show significant lower stabilization in their preferred conformation, or substrates that can be stabilized in a less preferred conformation in the enzyme, are likely to act as alternative substrates. Cyclopropylcarbinyl precursors were shown to be more likely to undergo turnover than heteroatom substituted analogues, based on their smaller radical stabilization and lower rearrangement barriers. Also, anchoring to Arg27 by a functional group of the ligands (preferably by carboxylates) is necessary to bring the ligands in optimal position (for initial hydrogen abstraction) and stabilize reactive conformations for potential turnover.
These results show that there is significant potential in the presented methodology to be used as a screening approach for enzyme engineering of radical SAM enzymes. Screening thermodynamic reaction profiles is also easier and quicker to perform than more costly transition state searches for multistep reactions, like the reaction presented here. The next step is using this methodology to screen for mutations in the active site further, either supporting turnover for alternative substrates or altering turnover for the natural substrate CPH4. This will serve as the entry point for the computational design of radical SAM enzymes, facilitating reactions with non-natural substrates for the generation of novel, bioactive compounds.

EXPERIMENTAL SECTION
Additional information and more detailed methods are provided in the Supporting Information. Simulation input files, important output files, analysis and scripts used for analysis are also provided on on the figshare repository (DOI: https://doi.org/10.6084/ m9.figshare.c.4290332). Computational workflows developed for the calculation of vRSE values from MD data on GitHub (https://github.com/ChrisSuess/RSE-Calc).

Molecular Dynamics Simulations
All molecular-dynamics simulations were performed using the GPU implementation [44][45][46] of the Amber16 47 molecular dynamics package. The force field parameters for SAM are based on electrostatic reparametrized force field parameters from Saez and Vöhringer-Martinez 48 as described and tested previously. 49 The parameters for the 4Fe4S cluster are based on a recent parametrisation of biological relevant iron-sulfur clusters by Carvalho and Swart. 50 These parameters showed good structural identity of the clusters as shown in the MD analysis in Section 1.4 of the Supporting Information. The interaction of SAM and the cluster was treated by electrostatic interactions only. No restraints have been applied to the cluster, the substrate, or SAM apart from the simulation equilibration.
All simulations have been performed in explicit solvent, using the SPC/E 51 water model. The simulations were conducted at a temperature of 300 K using periodic boundary conditions. Thorough restrained equilibration as described in Section 1.2 of the Supporting Information was followed by multiple simulations in between 100 and over 1 µs.

DFT and RSE calculations
All accurate RSE calculations for alternative substrates were performed analogue to our previous publication on QueE 21 and as previously outlined for reactions in radical SAM enzymes by Hioe and Zipse. 40 The determined RSE values were corrected with unscaled zero-point energies on the level of their geometry optimization. RSE energies were then calculated applying thermal corrections to enthalpies at 298.15 K at the level of their geometry optimization. All stationary points have been characterized by frequency calculations.
For the calculation of the vRSE values snapshots from the trajectories of the dynamic simulations were analyzed with single point energy calculations, at both a semi-empirical (SE) and density functional level of theory (DFT). All semi-empirical calculations reported were carried out using MOPAC 52 a semi-empirical quantum chemistry package based on Dewar and Thiel's NDDO approximation 53 with the PM6-D3 method which uses Grimme's D3 dispersion corrections for correlation. 43 All DFT calculations use the quantum chemistry package Q-Chem 54 at an M062X/6-31+G(d) level of theory. Grimme's D3 dispersion corrections are applied where: s6 = 1.0, sr,6 = 1.619 and s8 = 0.0. All energies have been corrected for dispersion by use of Grimme's D3 dispersion correction with correction values taken from literature.

Docking calculations
The prepared set of substrates was docked into the receptor QueE using a combined standard precision (SP) extra precision (XP) protocol with Glide 55-56 as described in detail in Section 3.1 of the Supporting Information. Following an exhaustive sampling search to predict orientation, conformation and binding position of a structure inside the rigid receptor pocket by Glide SP, the best conformers of each substrate were selected and docked with Glide XP to retain more accurate results. The OPLS3 force field 57 was used for the docking calculations and no constraints were applied.