Modeling and simulation‐based design of electroenzymatic batch processes catalyzed by unspecific peroxygenase from A. aegerita

Unspecific peroxygenases have attracted interest due to their ability to catalyze the oxygenation of various types of C–H bonds using only hydrogen peroxide as a cosubstrate. Due to the instability of these enzymes at even low hydrogen peroxide concentrations, careful fed‐batch addition of the cosubstrate or ideally in situ production is required. While various approaches for hydrogen peroxide addition have been qualitatively assessed, only limited kinetic data concerning enzyme inactivation and peroxide accumulation has been reported so far. To obtain quantitative insights into the kinetics of such a process, a detailed data set for a peroxygenase‐catalyzed benzylic hydroxylation coupled with electrochemical hydrogen peroxide production is presented. Based on this data set, we set out to model such an electroenzymatic process. For this, initial velocity data for the benzylic hydroxylation is collected and an extended Ping‐Pong‐Bi‐Bi type rate equation is established, which sufficiently describes the enzyme kinetic. Moreover, we propose an empirical inactivation term based on the collected data set. Finally, we show that the full model does not only describe the process with sufficient accuracy, but can also be used predictively to control hydrogen peroxide feeding rates To limit the concentration of this critical cosubstrate in the system.

UPOs, not unlike nonbacterial P450s, have been proven to be notoriously difficult to express at significant titers in the usual laboratory hosts (Bormann, Gomez Baraibar, Ni, Holtmann, & Hollmann, 2015). This is why the UPO from Agrocybe aegerita has become the de facto model UPO when Molina-Espeja et al. (2014) evolved an enzyme variant that can be recombinantly expressed in Saccharomyces cerevisiae (AaeUPO-PaDa-I, herein referred to as rAaeUPO) and further increased expression levels of the variant by transferring it to Pichia pastoris (Molina-Espeja, Ma, Mate, Ludwig, & Alcalde, 2015). While the wide substrate scope of the enzyme had already been somewhat explored with the native enzyme (Kinne et al., 2009;Kluge, Ullrich, Scheibner, & Hofrichter, 2012;Peter, Kinne, Ullrich, Kayser, & Hofrichter, 2013;Peter et al., 2011), the improved availability allowed for more detailed investigations into the use of this enzyme for synthetic chemistry purposes (Gomez de Santos et al., 2019;Ni et al., 2016;Perz et al., 2020). These efforts have focused on two main goals: achieving preparative product concentrations  and reducing the inactivating effects of hydrogen peroxide on the biocatalyst to increase biocatalyst efficiency (Burek, Bormann, Hollmann, Bloh, & Holtmann, 2019). While unsurpassed as an oxidant concerning its economic and "green chemistry" properties, hydrogen peroxide is a deleterious reactant as it quickly deactivates UPOs, which is why great effort has been put into the development of H 2 O 2 delivery system that compensate this downside and maximize the total turnover number (TTN) of the biocatalyst (Hernandez, Berenguer-Murcia, Rodrigues, & Fernandez-Lafuente, 2012).
To date, the highest amount of turnovers for rAaeUPO has been achieved using an enzymatic cascade that produces H 2 O 2 from methanol, achieving a TTN of 470,000 for a benzylic hydroxylation . The use of an electrochemical gas diffusion electrodebased in situ hydrogen peroxide production (GDE) system has been established as a viable alternative that can achieve up to 85% of the highest TTN of the enzymatic system (Horst et al., 2016). Potential advantages of the electrochemical system are the higher turnover frequencies achieved so far, the simpler process design and a higher flexibility concerning the hydrogen peroxide production rate. While great efforts have been made to elucidate the mechanism of UPO inactivation (Karich, Scheibner, Ullrich, & Hofrichter, 2016) and to identify an ideal means of providing hydrogen peroxide in situ, less data concerning the product formation kinetics and enzyme inactivation kinetics has been reported to date.
In this study, we sought to investigate an UPO-catalyzed hydroxylation reaction concerning reaction kinetics and enzyme inactivation kinetics while collecting a data set amenable to modeling with the goal of simulating such a process. Due to the flexibility and simplicity of the system, we chose to concentrate on the electrochemical in situ production of H 2 O 2 . This combination results in an electroenzymatic cascade reaction. Understanding, analysis, and control of enzymatic cascade reactions are challenging problems, for example, due to their complexities, nonlinearities, and delays. The investigation of a novel model-based control of biosystems is of crucial significance for the improvement of operational stability and processes efficiency. Given the industrial interest in this type of reaction, we focused our efforts on developing a model that could combine a reasonable level of accuracy with the requirement for a low amount of easy to generate experimental data.

| Chemicals
If not stated otherwise, chemicals were purchased from Sigma-Aldrich or Carl Roth in >98% purity. 4-(1-Hydroxyethyl)benzoic acid was purchased from BLD Pharm in >97% purity. Recombinant UPO from A. aegerita produced with P. pastoris was graciously provided by Prof. F. Hollmann (TU Delft) as crude enzyme .

| Electroenzymatic process
Electroenzymatic reactions were carried out in a magnetically stirred (120 rpm) glass reactor (250 ml working volume), which served as a reservoir for the reaction mixture. Reaction mixture was circulated through a flow-through GDE with a dead volume of 8 ml as described by Krieg, Hüttmann, Mangold, Schrader, and Holtmann (2011) at 40 ml min −1 using a peristaltic pump. The two electrode setup consisted of a pressed carbon black gas diffusion working electrode (Gaskatel) with a projected surface area of 5.5 cm 2 and a platinum counter electrode with an equal surface area. Currents were applied to the electrode using a Gamry Reference 600 potentiostat. From a bypass reaction medium was continuously pumped through a flowthrough sample chamber from which samples for hydrogen peroxide determination were withdrawn approximately every 120 s (vide infra). Reaction medium consisted of 10 mM 4-ethylbenzoic acid (EBA) in 100 mM potassium phosphate buffer pH 7. The pH was adjusted after the solubilization of EBA. The reaction volume was 100 ml. Before starting the current, rAaeUPO was added to a concentration of 12.5 nM and a sample was withdrawn for determination of enzyme concentration using the 2,2′-azino-bis(3-ethylbenzothiazoline-6-sulfonic acid) (ABTS) assay as described below ("100% active enzyme"). All activity measurements are given relative to this sample, as absolute concentrations could not be determined due to the simultaneous presence of EBA and ABTS in the assay. Experiments were initiated by applying a constant current to the GDE, except for simulation experiments, were a time-dependent current was applied. Samples for offline substrate and product analysis were immediately mixed with sodium azide (1 mM), diluted appropriately, and analyzed by high-performance liquid chromatography (HPLC; vide infra).

| Enzymatic assays
Enzyme concentrations of stock solutions were determined using the ABTS assay (Molina-Espeja et al., 2014). The assay consisted of 0.3 mM ABTS and 2 mM H 2 O 2 in 90 mM citric acid/Na-phosphate buffer pH 4.4 (100 μl sample volume). Samples were appropriately diluted to achieve an absorption increase of 0.1-1 cm −1 min −1 at 420 nm (ε ABTS = 36,000 cm −1 M −1 ). Enzyme concentrations were determined using the single substrate Michaelis-Menten equation with K M and k cat values as published before for AaeUPO-PaDa-I produced in P. pastoris (Molina-Espeja et al., 2015).
For the determination of initial velocity data concerning the hydroxylation of EBA, H 2 O 2 and 9.5 nM UPO in 100 mM potassium phosphate buffer pH 7 were mixed and incubated at 25°C for 1 min.
Every 20 s a sample was withdrawn and mixed with sodium azide to a final concentration of 1 mM. Samples were appropriately diluted and the 4-(1-hydroxyethyl)benzoic acid (HEBA) concentration determined as described below.

| Analytics
EBA and HEBA concentrations were determined using a NEXERA X2 line HPLC (Shimadzu). The instrument was equipped with a Luna C18 RP column (Phenomenex). The analysis was carried out at 35°C using a binary gradient of 0.1% vol/vol formic acid in water and acetonitrile (AcN) with the following gradient: 0 min, 35% AcN; to 7 min, 80% AcN; to 10 min, 35% AcN. Analytes were detected using a SPDM20A UV/vis detector (Shimadzu) at 237 nm by comparison to authentic standards. The sample volume was 3 μl.
Hydrogen peroxide concentrations were determined electrochemically using a Yellow Springs Instruments Select 2700 equipped with blank membranes and calibrated with a 30 mg L −1 H 2 O 2 solution in 20 mM phosphate buffer pH 7. Typical calibration samples yielded a current of 10 nA. The sample volume was 25 μl.

| Kinetic modeling
The model identification process was based on the systematic approach presented by Ohs, Leipnitz, Schöpping, and Spiess (2018). For the parameter estimation and simulation studies, the ordinary differential equations describing the reaction system were implemented in the MATLAB computing environment. The rough outline of the individual steps is given below. standing of the kinetic system akin to a Michaelis-Menten system, the microkinetic parameters were substituted by macrokinetic parameters using Cleland's notation (Cleland, 1963).

| Structural model identifiability
The differential equations describing the concentration changes over time were implemented in the DAISY-tool (Bellu, Saccomani, Audoly, & D'Angiò, 2007) to investigate whether the parameters were globally or only locally identifiable, or structurally nonidentifiable.

| Parameter estimation
Parameters were estimated using the Levenberg-Marquardt algorithm of the MATLAB-internal nonlinear least-square solver isqnonlin by minimizing the residual sum of squares (RSS) between simulated and experimental data. Initial parameter guesses for the solver were obtained by first generating 100,000 random parameter sets in the range of 0.00001 ≤ k i ≤ 100,000 in the parameter specific units and selecting the set with the smallest residual sum of squares.

| Model evaluation
The results of the parameter estimation step were used as the basis for statistical analysis of the model to evaluate its predictive quality for subsequent simulations. The relative standard deviation (relSD) of each parameter k i was used as a measure of confidence in the value and was calculated according to the following equation: with variance v i as the ith diagonal element of the variance-covariance matrix V(θ), which was approximated by the following equation: Here, J is the Jacobian matrix as one output from the lsqnonlin solver and describes the parameter sensitivities, and s² is the estimated measurement variance (Asprey & Naka, 1999) defined by the following equation: The variance is calculated with the number of measured data points n and the number of parameters p.
Additionally, the Akaike information criterion (AIC) was used as a means for evaluating the goodness of fit of individual model candidates. The AIC was calculated via the following equation: The number of parameters in a model serves as a penalty in the calculation of the AIC to avoid reaching a good fit by overfitting.
A lower value therefore suggests a model that describes the experimental data well with a minimal number of different parameters.

| Process design
To gain better insights into the kinetics of UPO process performance and enzyme inactivation, we opted to carry out the benzylic hydroxylation of EBA to HEBA coupled with electrochemical in situ hydrogen peroxide production using the GDE system. This reaction was chosen since it is similar to the model reaction of ethylbenzene to 1-phenyl ethanol that has routinely been used for the evaluation of various UPO process strategies so far. In contrast to ethylbenzene, EBA is readily soluble in buffered aqueous systems and not prone to evaporation or adsorption, thus affording the possibility for a closed mass balance over the course of the whole reaction, which has proven to be difficult when ethylbenzene was used as a substrate. In contrast to previous work, we also employed automated at line hydrogen peroxide analysis to achieve better data coverage. Sampling points for offline substrate and product analysis and enzyme activity were chosen so that the same number of data points were collected for all experiments, independent of process duration, to avoid biases when fitting the model. The reactor setup consisted of a magnetically stirred reaction vessel, which served as a reservoir from which the reaction medium was continuously pumped through a flow-through GDE. Samples for automated H 2 O 2 analysis were taken from the reaction vessel through a second bypass.
Before experiments involving enzyme, the reactor setup was characterized abiotically in terms of electrochemical hydrogen peroxide production. As it has been shown before, current and hydrogen peroxide production rates were mostly linearly correlated ( Figure S1; Horst et al., 2016;Krieg et al., 2011).

| Process data
With the described setup, UPO-catalyzed EBA hydroxylation reactions were carried out at various hydrogen peroxide production rates by applying currents ranging from −10 to −80 mA to the GDE, which corresponds to approximately 1.3-10.5 mM h −1 H 2 O 2 ( Figures   S2-S9). Exemplarily, processes carried out at −30 and −60 mA are shown in Figure 1.
In the beginning of these reactions, the hydrogen peroxide concentration quickly increased to a pseudo-steady-state concentration, presumably the concentration at which enzyme velocity and hydrogen peroxide production rate were equal. Following this plateau, hydrogen peroxide continued to accumulate, albeit at a slower rate, which was most probably caused by enzyme inactivation. Product formation rates were near-linear for the ma-

| Modeling
Modeling of enzyme-based processes is a key tool in process development to support optimization and large scale operation. We We specifically chose not to record oxygen production rates (i.e. the product of the catalase reaction) due to the inherent difficulty in obtaining such data in a reliable way, given our goal of establishing an UPO model that can be created with a limited set of experiments and a setup that is readily available in a regular biochemical laboratory.
F I G U R E 2 Summary of the process performance of rAaeUPO catalyzed EBA at various currents. TOFs ( ) were determined from the part of experiments were product formation rates were approximately linear, that is, as initial frequencies.  -Bowden, 1979). This led us to include this substrate inhibition into the reaction mechanism (Figure 3, yellow). We therefore chose to use the full catalase mechanism including inhibition by the organic substrate (Figure 3, green, blue, and yellow) for initial data fitting.
While the catalase malfunction pathway (Figure 3, red) would have been a superior choice as it can mechanistically describe both product formation and enzyme inactivation, it was disregarded in this study, since it does not conform to the steady-state assumption that is required for the derivation of kinetic equations using methods such as proposed by King and Altman (1956). Moreover, it would require data difficult to obtain from routine initial velocity experiments, that is, the simultaneous determination of active enzyme concentrations and initial velocity data, which is unfeasible given the fast sampling times of such experiments.
Rate equations of the full catalase mechanism including substrate inhibition were derived using the method of King and Altman (1956). Estimating the differential equation system's kinetic parameters with the obtained initial velocity data, it became obvious that the lack of oxygen production rates resulted in very high parameter uncertainty of k 6 (catalase model in Having established a kinetic equation that should sufficiently describe product formation and hydrogen peroxide consumption rate, we sought to establish an enzyme inactivation kinetic. The inactivation parameters were estimated by fitting time-dependent process data with the inhibCPDI model as the base for reaction kinetics. Modeling the enzyme inactivation using a single inactivation rate based on a first or higher-order hydrogen peroxide concentration did not result in a satisfactory fit of the obtained experimental data as the enzyme inactivation was usually underestimated at lower hydrogen peroxide production rates. It has been proposed before that P450-like hydroxylation reactions exhibit an inherent "maximum" number of turnovers, independent of extrinsic inactivating factors  Brummund et al., 2016). The existence of a similar hydrogen peroxide concentration-independent limit for UPOs is supported by the data obtained by Horst et al. (2016). Therein, they showed that in an electroenzymatic system similar to that employed in this study, a reduction in hydrogen peroxide reduction rate below a certain threshold could not further increase the maximum TTN (400,000).
The work of Ni et al. (2016), reported a relatively similar TTN (470,000) for the same reaction in a completely different reaction setup at significantly lower hydrogen peroxide production rates, which gives rise to the notion that an infinitely low hydrogen peroxide production rate will not result in ever-higher TTN for this type of catalyst.
Based on this hypothesis and given the near linear correlation between hydrogen peroxide production rate and TTN observed in our experimental dataset (Figure 2), we introduced an empirical term describing a maximum TTN that would be achieved at an infinitely low H 2 O 2 production (i.e., the intercept of the linear regression of the TTN at 0 mA, TTN max = 400,000). This TTN based term was combined with an H 2 O 2 -based inactivation rate as well as the previously established reaction kinetics. This resulted in the ordinary differential equation system (ODE) depicted in Equations (5)-(8). Therein, the substrate consumption rate is presumed to be equivalent to the product formation rate, given that the carbon mass balances were closed and no product overoxidation to the corresponding ketone was observed. The H 2 O 2 production rate is based on the current (I) dependent formation rate (K H O 2 2 ) and the H 2 O 2 consumption rate, which is equal to the product formation rate. The product formation rate is described by the kinetic equation of the simplified double inhibition mechanism as discussed above, where K M depicts the Michaelis-Menten constants of the respective substrates given in the subscripts and K i depicts their inhibition constants. Enzyme inactivation is described by a thirdpower inactivation constant (k inact ), the hydrogen peroxide concentration in the third power, based on the catalase malfunction mechanism ( Figure 3) and the empirical term describing the maximum turnover number interpolated from the experimental data set (vide supra).
Using this set of ODEs, it was possible to reflect our experimental data set with very high accuracy. However, we noticed minor discrepancies between the total amount of experimentally measured H 2 O 2 in the system (sum of free compound and HEBA) and the projected H 2 O 2 produced by the GDE, which indicates that the real currents applied during the experiments slightly deviated from the set values. By fitting the projected H 2 O 2 production to the experimental values with different currents, the discrepancies could be eliminated. Figure 5 shows the experimental data and the respective simulations with adjusted H 2 O 2 production rates. The value of the H 2 O 2 dependent inactivation constant k inact was estimated at 2.12 × 10 −2 mM −3 min −1 with a low relative standard deviation of 5%.
Note that the parameters for the reaction kinetics were estimated with a completely different set of initial velocity data but still predict process data with such high precision. This further strengthens our assumption that the inhibCPDI model includes all relevant mechanisms in a relatively simple and easy to use form. Next, we investigated if the model can be used for process control under novel conditions.
The previous experimental data set had been produced by applying a constant current to the GDE, that is, at constant hydrogen peroxide production rates. As is apparent from the collected data, this leads to ever-increasing H 2 O 2 concentrations due to enzyme inactivation, which in turn results in increasing enzyme inactivation rates. Therefore, we chose to use the model to predict a variable current that would result in an increase in hydrogen peroxide concentration only up to a certain level, that is, a H 2 O 2 -stat mode. This operational mode has before been shown to be superior for peroxygenase catalyzed processes (Seelbach, van Deurzen, van Rantwijk, Sheldon, & Kragl, 1997).
For this, two cases of progress curve data were simulated with an initial current of −80 mA for fast H 2 O 2 production in the beginning, and the constraint to keep the upper H 2 O 2 concentration limited to arbitrarily chosen values 0.5 or 1.2 mM for a proof-of-concept.
The time step for current adjustments was set to 1 s. The simulations were experimentally validated and are shown in Figure 6. In both cases, the current stays at the initial value of −80 mA, resulting in the maximal production rate of H 2 O 2 , until the respective concentration  (Figure 6a). The deviation between simulated and experimental data in Figure 6b can be explained by an unsteady H 2 O 2 production in the GDE during this experiment, for example, caused by the bubble formation, changes in the ambient temperature or aging of the electrode. The high agreement between experimental and predicted data shows that the established model can not only be used descriptively but also predictively and thus lays