A Customized Bayesian Algorithm to Optimize Enzyme-Catalyzed Reactions

Design of experiments (DoE) plays an important role in optimizing the catalytic performance of chemical reactions. The most commonly used DoE relies on the response surface methodology (RSM) to model the variable space of experimental conditions with the fewest number of experiments. However, the RSM leads to an exponential increase in the number of required experiments as the number of variables increases. Herein we describe a Bayesian optimization algorithm (BOA) to optimize the continuous parameters (e.g., temperature, reaction time, reactant and enzyme concentrations, etc.) of enzyme-catalyzed reactions with the aim of maximizing performance. Compared to existing Bayesian optimization methods, we propose an improved algorithm that leads to better results under limited resources and time for experiments. To validate the versatility of the BOA, we benchmarked its performance with biocatalytic C–C bond formation and amination for the optimization of the turnover number. Gratifyingly, up to 80% improvement compared to RSM and up to 360% improvement vs previous Bayesian optimization algorithms were obtained. Importantly, this strategy enabled simultaneous optimization of both the enzyme’s activity and selectivity for cross-benzoin condensation.

The optimization of experimental reaction conditions to maximize the figures of merit of a catalytic reaction (FoM, i.e. yield, turnover number, selectivity, rate etc.) plays an essential role both in academic and industrial settings. 1The outcome of a chemical reaction is governed by a complex network of interactions between reactants, catalysts, solvent, and other ingredients, as well as temperature, pressure, pH, etc.To maximize a given FoM, multiple continuous variables need to be optimized.Since most variables influence each other, it is challenging to determine the global optimum by optimizing individual variables one by one (i.e.one-factor-at-a-time, OFAT, Figure 1A). 2,3When considering interactions among variables, the number of experiments required for the optimization increases exponentially as the number of parameters increases.Therefore, in the majority of catalyst optimization campaigns, it is impossible to comprehensively evaluate all possible combinations of experimental variables.
][8] One of the most commonly-used methods in DoE is the response surface methodology (RSM). 2,9,10This method computes a response surface model (relying on an approximate equation) between the experimental conditions and FoMs, then improves the accuracy of this model through a limited number of experiments.Various strategies have been proposed to reduce the number of experiments, such as Box-Behnken design, 11,12 central composite design, 13,14 etc. 15 Using a software package, the experimenter can automatically select the best of these strategies and apply DoE without any need for in-depth prior knowledge of the mechanism. 2owever, RSM remains challenging as the number of experiments grows exponentially with a linearly-increasing number of variables.If the amount of data is insufficient for the complexity of the model, the optimized model will be inaccurate, and over-fitted only on existing data.Therefore, when optimizing multiple variables simultaneously, there is limited accuracy in modeling the complex relationships among these variables while relying on a reasonable number of experiments.
One of the reasons for the massive increase in experimental cost is that the sampling points are selected nearly homogeneously from the entire parameter hyperspace.Indeed, RSM is a single-iteration optimization process, whereby all experiments are performed simultaneously to determine the optimal parameters.As the process relies on a single iteration, all regions of the parameter hyperspace are sampled without bias.This procedure thus leads to sampling portions of the parameter hyperspace that are inherently not worth exploring (Figure 1B). Figure 1 Comparison of the sampling process.(A) In OFAT, one variable of the reaction conditions is selected and optimized while the other parameters are held constant.By iterating this process, the variables are changed one by one.(B) In the RSM, the entire hyperspace of reaction conditions is sampled homogeneously.Then, analysis of the single-iteration leads to the prediction of the optimized reaction conditions for a given FoM.Accordingly, regions of hyperspace leading to sub-optimal FoM are sampled unnecessarily.(C) In Bayesian optimization, experimental sampling and prediction are iteratively repeated.For each iteration, the reaction conditions are suggested for which the FoM is predicted to be best.Accordingly, the better the resultant hyperspace, the denser the sampling density is applied.
7][18][19] In this method, the Gaussian process regression 20,21 first predicts the FoM value and its uncertainty under untested reaction conditions.From these data, and with the aim improving the FoM, the algorithm computes which conditions should be experimentally tested.For this purpose, an evaluation formula, coined acquisition function (AF), is applied.In the evaluation of the FoM, two types of strategies are applied: exploration and exploitation.(1) exploration of areas that have not yet been tested and (2) exploitation of good conditions that have already been obtained.By iteratively repeating the process of prediction and experiments, regions with greater potential for an improved FoM can be preferentially explored, enabling a more efficient optimization (Figure 1C).In addition, in a Gaussian process regression, a virtually infinite dimension of functions are assumed as a model, and the likelihood of the model functions are computed based on Bayes' theorem, according to the existing data. 21Compared to the conventional RSM whose fitting function is limited to a few dimensions, it is expected that Bayesian optimization can identify complex relationships among variables under experimental conditions.Due to these advantages, the Gaussian process regression and Bayesian optimization are expected to have potential applications in a variety of fields such as flow chemistry, 22,23 automated chemical design, 24 adjustment of machine learning systems 25,26 or experimental device, 27 materials engineering, 28 biological engineering, 29,30 etc.2][33] In this study, we compared the performance of RSM and Bayesian optimization to optimize an FoM for two enzyme-catalyzed reactions.In the process, we identified limitations of the acquisition function used in existing Bayesian optimization 17 and adjusted it.The resulting Bayesian Optimization Algorithm (BOA) led to the rapid identification of significantly-improved FoMs for enzyme-catalyzed reactions.The BOA workflow we applied in this study is outlined in Figure 2. Initially, the FoM and continuous experimental variables tobe-optimized are identified, and initial experimental data are collected, either from experiments or from the literature (A).In a Gaussian-process regression, the fitting curve is not fixed via a single type of equation, but various functions are assumed as stochastically valid.Based on the experimental data, the Bayesian inference computes predicted values for the FoM and their uncertainties for untested experimental conditions (B).An AF is computed to identify the values of the experimental conditions, which should be tested in the next iteration (C).The data point that possesses the highest value of the AF is selected and the resulting experimental data resulting from these conditions are added to the data set (return to A).Then, the Gaussian process-regression and the evaluation of the AF are performed again (B, C).This iterative process is repeated until a targeted FoM is reached.As more data are collected, the model is refined and becomes more accurate: the probability of identifying improved experimental conditions increases accordingly.In an ideal Bayesian optimization, only the single point with the largest AF for each iteration is evaluated, thus significantly slowing the iterative optimization process.To expedite the procedure, we adopted a batch optimization process to evaluate a handful of reaction conditions during each iteration.The Kriging believer algorithm 34 was selected as batch optimization method.In this method, when identifying the second and subsequent experimental conditions to be tested, Gaussian-process regression is performed including conditions identified in the previous iterations (Figure S1).
9][40][41][42] The optimization was performed with the aim of maximizing the total turnover number (TON), used as the FoM.Five parameters were selected as continuous variables for both the PAL-and BFD-catalyzed reactions (Scheme 1).In enzymatic reactions, enzyme and substrate concentrations are directly related to the TON, and the pH often has a significant effect on enzyme activity.

Scheme 1 Selected enzyme-catalyzed reactions used to benchmark and improve the performance of the DoE algorithms.
For both reactions, DMSO was used as a co-solvent to increase the solubility of the substrate in the aqueous medium.For BFD, the concentration of the TPP cofactor was selected as a variable.Indeed, TPP is abundant in E. coli cells 43 , thus maximizing the formation of the holoenzyme TPP (TPP displays rather low affinity for its cognate apoenzyme).Accordingly, TPP concentration has been shown to play a critical role in maximizing the TON with purified enzymes. 36,44or the PAL-catalyzed reaction, the NH4 + concentration was selected as a variable as it acts as a co-substrate for the hydroamination of trans-cinnamic acid.All other factors, including buffer concentration, temperature, reaction volume, and reaction time, were kept constant throughout the DoE.Based on the experimental table generated by MODDE ® , 29 experiments were performed (Table S1, S2).The contour plot for the BFD reaction at a TPP concentration of 1 mM is displayed in Figure S2.As can be appreciated, the response surface increases monotonically with increasing substrate concentration and decreasing enzyme concentration.The activity of the enzyme tends to decrease with higher concentrations of DMSO, and the BFD performs best at pH 8.With the help of the program, optimal experimental conditions were proposed.Gratifyingly, when the BFD reaction was performed using the predicted best conditions, a TON = 3289 was obtained, significantly higher than the TON = 2776 predicted by MODDE ® for these experimental conditions (Table 1, S3).For the PAL reaction, experimental tables and predicted optimal conditions were also prepared (Table S2, S4).The response surfaces in the presence of 5% DMSO are displayed in Figure S3.As expected, the enzyme's performance is poor at low NH4 + concentrations.At pH 8, low enzyme concentration combined with high substrate concentration result in high TON.However, the PAL reaction under the predicted optimal conditions affords a TON = 1544, which is significantly lower than the predicted TON = 2764 (Table S4).Next, a Bayesian optimization was performed, starting with five initial reaction conditions randomly selected.Unexpectedly, the optimization process led to the identification of optimized reaction conditions resulting in significantly lower TON compared to the RSM described above: TON = 1573 (for BFD) and TON = 705 (for PAL).The five reaction conditions proposed for each iteration of optimization turned out to be similar to each other and concentrated around the reaction conditions previously explored (Table S5, S6).Inspection of the AF, eq. 1, reveals the following trends of the algorithm: i) it leans strongly toward the exploitation of experimental conditions resulting in high FoM within the existing dataset and ii) neglects sampling the unexplored regions of the experimental hyperspace (Figure 3A).We hypothesized that this AF might be challenged to identify the global optimum conditions.To address this issue, we considered modifying the AF and the corresponding sampling algorithm used in the above Bayesian optimization.As can be appreciated from eq. 1, the AF includes an evaluation term for i) the exploration of untested regions and ii) the exploitation of the best result in the existing data.With a rapid optimization process in mind, it is indispensable to maintain a balance between these two factors.Vari-ous formulas have been used to address this delicate balance including: Upper/Lower Confidence Bound, 45 Probability of Improvement, 45,46 and Expected Improvement. 46,47nspired by the previous report 17 , we selected the Expected Improvement function (EI), eq. 1.Where f(x + ) is the current best data, μ(x) is the predicted mean value, σ(x) is the predicted standard deviation, Φ and ϕ represent the cumulative distribution function and the probability density function of standard normal distribution, and ξ is a parameter to determine the balance between exploration and exploitation, which is set to 0.01.The Expected Improvement function computes how much the FoM will improve compared to the best FoM resulting from the experimental data.Shields et al. reported that this function enables the efficient optimization of chemical reactions. 17o overcome the challenge of over-emphasizing exploitation vs. exploration 48 observed with both BDF-and PALcatalyzed reactions (Table S5, S6), we sought to implement a mechanism that favors a more active exploration of the untested experimental regions.With this goal in mind, we evaluated two complementary strategies (i.e.BOA-1 and BOA-2 respectively) to balance the exploration/exploitation components of the AF.

𝐸𝐼(𝒙
For BOA-1, the algorithm used to determine the next round of experimental conditions was modified while keeping the AF unchanged.In the batch optimization process, the second and subsequent proposed experimental conditions are selected among points that are sufficiently distant from the previously-sampled conditions (Figure 3B).This strategy avoids accumulating experimental conditions located in the same area during the same iteration.In this way, it is expected that a globally optimal solution may be reached more rapidly.Figure 3 Modification of the algorithm and AF in BOA-1 and BOA-2.In the conventional method (A), the prediction leans toward the exploitation of experimental conditions resulting in high FoM within the existing dataset (red cross) and neglects sampling the unexplored regions of the experimental hyperspace (gray area).In BOA-1 (B), when determining the next experimental conditions to test, points close to the previously evaluated conditions (in blue disk) are excluded.This method avoids concentrating suggested points in the same area during one iteration.In BOA-2 (C), the AF is modified.By squaring the standard deviation of the predictions (σ 2 (x)), the function is expected to give more emphasis to regions far away from the existing data, which have high uncertainty.It enables fully exploring the entire hyperspace spanned by the experimental conditions.
For BOA-2, the algorithm used to determine the next round of experimental conditions was kept, and the Expected Improvement function was modified, as illustrated in Figure 3C.By squaring the standard deviation of the predictions, σ 2 (x), the experimental values in areas of higher uncertainty are given higher priority.This favors a selection of experimental conditions to be tested that are farther away from existing data points.In applying this function, FoM values in the existing dataset were scaled so that the maximum value is 10 to prevent the predicted uncertainty from becoming too large.
The results of BOA-1 and BOA-2 for both BDF-and PALcatalyzed reactions are compiled in Table 1 (entries 3,4) and Table S7-S10.For both enzymatic reactions, BOA-1 and BOA-2 afforded higher TONs than either RSM and the conventional Bayesian optimization.Moreover, the number of required experiments tended to be lower for BOA.In the case of BOA-1, the highest FoM from the RSM was matched in the second iteration already, ultimately leading to the identification of significantly-improved TONs after the five iterations systematically applied in this study (Table S7, S8).In the latter half of the optimization, when fine-tuning of parameters was required, the search efficiency was reduced by the limitation that closely-related experimental conditions were not proposed by the BOA during this iteration (See Figure 3B).For example, in the BFD reaction, the optimal solution for RSM was reached in the second iteration, but subsequent fine-tuning of the pH required three additional iterations to ultimately achieve the highest TON.In contrast, the BOA-2 achieved high efficiency without limiting the search range (Table S9, S10).In the first iteration, experimental conditions were selected mainly at the boundaries of the search range for each variable.This enables examining trends in the entire variable hyperspace.In the latter part of the optimization, variables were fine-tuned while remaining in the high-FoM region.It reveals that the exploration of untested region and exploitation of existing best results is well balanced.Based on these observations, the BOA-1 tends to be effective during the early optimization stage and the BOA-2 in the late optimization stage.Having identified promising BOAs, we set out to test these towards the optimization of two FoMs simultaneously.0][51][52] With this goal in mind, we selected the benzaldehyde lyase carboligation between benzaldehyde and isobutyraldehyde (cross-BAL). 53wo α-hydroxyketones may result from cross-carboligation.While the enantioselectivity of both products mostly exceeds 92 % ee, the ratio of the yield of these products is 1:1, leaving plenty of room for the improvement of chemoselectivity. 53We thus set to optimize simultaneously both chemoselectivity (eq.2) and TON (Scheme 2).

𝐶ℎ𝑒𝑚𝑜𝑠𝑒𝑙𝑒𝑐𝑡𝑖𝑣𝑖𝑡𝑦 = |𝑌𝑖𝑒𝑙𝑑 𝐴−𝑌𝑖𝑒𝑙𝑑 𝐵|
+  × 100% eq. 2 Considering the example of the Expected Improvement modification for optimization under constraints with a second FoM value, 54 the AF was modified to identify experimental conditions whereby the chemoselectivity (FoM 2) is expected to be greater than a set threshold, eq. 3.
Where EI(x) is the function used in BOA-1 or 2, μ'(x) and σ'(x) are the predicted mean value and predicted standard deviation for the FoM 2. The threshold is the target value of FoM 2 (i.e., chemoselectivity).This value should be set by weighing the relative importance of the two FoMs: the higher the threshold value for FoM 2, the less accurate the prediction of the main FoM 1 will be.For this validation, we set a higher emphasis on the TON (i.e., FoM 1) aimed at reaching a > 70% chemoselectivity, corresponding to a > 85:15 ratio of isomers.
The number of variables in the experimental conditions was increased from the two previous examples from five to six.This resulted in a larger experimental table for optimization by RSM and required 47 experiments (Table S11).In general, the higher the TON, the lower the chemoselectivity.There were however a few experimental conditions for which both FoMs were improved simultaneously.The best score was TON = 1426 and 94% chemoselectivity.However, the optimal reaction conditions predicted by the software based on these measured data did not lead to improved FoM (Table S12).The response surfaces generated for each of the TON and chemoselectivity reveal very different trends for the six parameters to achieve their respective maxima (Figure S4, S5).This highlights the difficulty of optimizing both FoMs simultaneously using the RSM method.In contrast, the results obtained using BOA-1 and BOA-2 outperformed RSM while requiring only 10-25 experiments -vs.47 for the RSM.-The maximum TON = 5441 with a chemoselectivity = 90.5 %. (Table 2, 3, S13).
For both BOA-1 and BOA-2, the algorithms succeeded in identifying the optimal conditions for RSM within a smaller number of experiments, and then reached even better conditions through further rounds of iterations.From these re-sults, it can be concluded that Bayesian optimization is potentially more efficient than RSM in enzyme reaction optimization.Moreover, the modifications that we introduced in the AF enables addressing the exploration-exploitation bias encountered in conventional Bayesian optimization. 17OA-1 favors more extensive exploration within a single iteration.Thus, it is more likely to yield valuable data in the early stages of optimization.BOA-2 provides effective suggestions, by balancing the exploration of unknown regions with the exploitation of existing results.The python scripts used can be downloaded from GitHub (see Supporting Information).
As highlighted in the last example whereby two FoMs are optimized simultaneously, the method based on Bayesian optimization can be adapted to the optimization of complex and highly interdependent experimental conditions.By tailoring the AF and algorithm which determines the condition to be tested, each FoM can be optimized using the required evaluation method.This method can be applied not only to reactions with purified enzymes, but also to any chemical reactions, including reactions in complex systems such as cells.We anticipate that straightforward optimization algorithm of the experimental conditions based on BOA will reduce the cost of optimizing experimental conditions, which can be the rate-limiting step in the development of chemical reactions and the engineering of biological systems.by the Swiss National Science Foundation.TRW further acknowledges support from an advanced ERC grant (DreAM: the Directed Evolution of Artificial Metalloenzymes, grant agreement 694424).Novartis University Basel Excellence Scholarships for Life Sciences (grant number 3CH1046).

Notes
The authors declare no competing financial interest.

ACKNOWLEDGMENT
This research was supported in part by a grant from the Naito Foundation (to RT).This publication was created as part of NCCR Catalysis (Grant 180544), a National Center of Competence in Research funded by the Swiss National Science Foundation.TRW further acknowledges support from an advanced ERC grant (DreAM: the Directed Evolution of Artificial Metalloenzymes, grant agreement 694424).SB thanks Novartis University Basel Excellence Scholarships for Life Sciences (grant number 3CH1046).ABBREVIATIONS AF, acquisition function; BFD, benzoylformate decarboxylase; BOA, Bayesian optimization algorithm; cross-BAL, benzaldehyde lyase; DMSO, dimethyl sulfoxide; DoE, design of experiments; EI, expected Improvement; FoM, figures of merit; OFAT, one-factor-at-a-time; PAL, phenylalanine ammonia lyase; RSM, response surface methodology; TON, turnover number; TPP, thiamine pyrophosphate.

Figure 2
Figure2Overview of the Bayesian Optimization Algorithm (BOA).This iterative process is repeated until the figure of merit (FoM) reaches a target value.(A) An initial dataset (black crosses) is experimentally acquired, enabling the algorithm to build a model.As fitting curve, various types of functions are assumed (dotted lines) and their likelihood is computed by Bayesian inference.(B) A Gaussian-process regression computes the predicted FoM value (solid red line) and its uncertainty (blue area) for untested experimental conditions.(C) By applying the acquisition function (AF, orange line), the experimental conditions predicted to improve FoM most (white cross) is computed.The FoM for these experimental conditions are i) computed (white cross), ii) validated experimentally, and iii) added to the existing experimental data set.The process is iterated until the FoM reaches a target value.

Table 1
Comparison of the optimized reaction conditions determined for the BFD and PAL-catalyzed reactions using DoE algorithms (*MODDE ® , **Bayesian optimization with conventional AF and BOA).

Scheme 2 Table 3
The benzaldehyde lyase as a testbed for the simultaneous optimization of two Figures of Merit: total turnover number and chemoselectivity.Table 2 Summary of the iterative optimization process for the cross-BAL-catalyzed reaction between benzaldehyde and isobutyraldehyde using BOA-1.Comparison of the optimized reaction conditions resulting from DoE based on the RSM, BOA-1 and BOA-2 for the cross-BAL-catalyzed reaction.