Three-level hybrid modeling for systematic optimization of biocatalytic synthesis: α-glucosyl glycerol production by enzymatic trans-glycosylation from sucrose.

Mechanism-based kinetic models are rigorous tools to analyze enzymatic reactions, but their extension to actual conditions of the biocatalytic synthesis can be difficult. Here, we demonstrate (mechanistic-empirical) hybrid modeling for systematic optimization of the sucrose phosphorylase-catalyzed glycosylation of glycerol from sucrose, to synthesize the cosmetic ingredient α-glucosyl glycerol (GG). The empirical model part was developed to capture nonspecific effects of high sucrose concentrations (up to 1.5 M) on microscopic steps of the enzymatic trans-glycosylation mechanism. Based on verified predictions of the enzyme performance under initial rate conditions (Level 1), the hybrid model was expanded by microscopic terms of the reverse reaction to account for the full-time course of GG synthesis (Level 2). Lastly (Level 3), the application of the hybrid model for comprehensive window-of-operation analysis and constrained optimization of the GG production (~250 g/L) was demonstrated. Using two candidate sucrose phosphorylases (from Leuconostoc mesenteroides and Bifidobacterium adolescentis), we reveal the hybrid model as a powerful tool of "process decision making" to guide rational selection of the best-suited enzyme catalyst. Our study exemplifies a closing of the gap between enzyme kinetic models considered for mechanistic research and applicable in technologically relevant reaction conditions; and it highlights the important benefit thus realizable for biocatalytic process development.

SucP from Bifidobacterium adolescentis

fructose

α‐d‐glucose 1‐phosphate

α‐d‐glucosyl glycerol

glucose

glycerol

SucP from Leuconostoc mesenteroides.

phosphate

sucrose

sucrose phosphorylase

INTRODUCTION

Enzyme kinetic studies of approximately 50 years ago have provided theory, together with the corresponding algebra, to describe the initial rate behavior of enzymatic transformations from the underlying microscopic reaction steps (Cleland, 1990; Cornish‐Bowden, 2012; Johnson, 2013; Segel, 1993). The rate laws for the full‐fledged kinetic mechanism are rigorous tools of reaction analysis and optimization (e.g., Berendsen et al., 2006; Buchholz et al., 2019; Kasche et al., 1987; Ohs et al., 2019; Rios‐Solis et al., 2015; Straathof and Heijnen, 1996; Sun et al., 2015; Willeman et al., 2000; Youshko et al., 2002; Zavrel et al., 2008). Often derived for mechanistic inquiries, however, these rate laws are not widely applied, with exceptions noted above, to model‐assisted development of biocatalytic processes (Gernaey et al., 2010; Vasić‐Rački et al., 2011). Mechanistic‐kinetic models reveal their practical importance when apparent (lumped) Michaelis–Menten parameters no longer capture the essential complexity of the overall transformation (e.g., Kasche, 1986; Ohs et al., 2019; Youshko et al., 2002). Chemical group transfer reactions catalyzed by hydrolases represent this situation characteristically and have considerable importance in industry (Adlercreutz, 2017; Giordano et al., 2006; Vera et al., 2020). The enzymatic production reflects the actual synthesis reaction overlapped with the hydrolysis of substrate, product or both, as illustrated in Figure 1 on the example of α‐glucosyl glycerol (GG) produced from sucrose (Suc) and glycerol (GOH) (Klimacek et al., 2020). Parameters of reaction efficiency and selectivity can exhibit complex dependence on key process variables, such as the concentrations of donor and acceptor substrate or the molar ratio of the two (Kasche et al., 1987; Klimacek et al., 2020; Youshko et al., 2002). Lumped kinetic parameters (e.g., the donor substrate ) explain these dependencies only to the extent that they are understood as composites of the individual reaction steps (Buchholz et al., 2019; Klimacek et al., 2020; Ohs et al., 2019; Rios‐Solis et al., 2015), hence the need for a mechanistic‐kinetic model.

Figure 1

Kinetic mechanism of SucP‐catalyzed glycosylation of GOH from Suc (1) for GG (2) synthesis. Microscopic rate constants affected by high [Suc] are indicated with asterisks. α,α‐Trehalose (3) is structurally similar to sucrose but unreactive in the enzymatic reaction. It was used to mimic the general solute effect of high [Suc]

Limitation on the applicability of mechanistic‐kinetic models can arise when their extension to the actual conditions of the biocatalytic synthesis proves difficult (Gonçalves et al., 2002; Grosch et al., 2017; Wandrey et al., 1979; Zhang et al., 2019). Although bulk conditions (e.g., pH, temperature, ionic strength, and enzyme preparation) may be invariant in the course of development, the high substrate concentrations of the synthetic process are typically avoided in studies of the kinetic mechanism, bearing in mind the possible convolution of specific and nonspecific effects of substrate used at high concentration. To close the gap in the use of mechanistic‐kinetic models in fundamental enzyme studies and in enzyme applications in bio‐catalysis, we here propose an approach by hybrid modeling. The hybrid model expands the mechanistic‐kinetic model by an empirical description of the effect of the real process conditions on individual steps of the kinetic mechanism. Hybrid models have attracted considerable interest across different fields of biochemical engineering (Narayanan et al., 2019; Smiatek et al., 2020; Zhang et al., 2019), but their use in applied bio‐catalysis research has, to the best of our knowledge, not been shown (Bulik et al., 2009).

Here, we demonstrate hybrid modeling at three levels coordinated in content for systematic step‐by‐step development of a thoroughly optimized enzymatic production of GG. Glycoin® (50% solution of GG; bitop AG) is a commercialized cosmetic ingredient marketed for its excellent skin‐moisturizing properties (Breitenbach et al., 2006; Iki et al., 2007; Novejarque, 2012). The GG is produced at multiton/year industrial scale via a trans‐glycosylation between sucrose and glycerol catalyzed by sucrose phosphorylase (SucP) (Goedl et al., 2008). At the first level of hybrid modeling, the mechanistic‐kinetic model of the SucP reaction (Klimacek et al., 2020) was expanded with an empirical description of the kinetic effect of high sucrose concentrations (≥100 g/L) used for GG production. Using two representative SucP enzymes from Leuconostoc mesenteroides (LmSucP) and Bifidobacterium adolescentis (BaSucP) (Franceus & Desmet, 2020; Schwaiger et al., 2021), we show that the empirical model part required enzyme‐specific parametrization, but was in each case indispensable to account for the enzyme kinetic behavior under real process conditions. The second level of hybrid modeling moved from the analysis of initial rates to the integrated description of full reaction time courses up to high degrees of substrate conversion with large concentrations of product accumulating (for the general case, see Arcos et al., 2001; Bauer et al., 1999; Buchholz et al., 2019; Chen et al., 2008; Flores & Halling, 2002; Johnson, 2009; Rakels et al., 1994; Straathof, 2001; Straathof & Heijnen, 1996). We show that evidence from both levels of the study was important to inform the rational selection of the best‐suited candidate enzyme. The third level of development involved the application of the validated hybrid model for computational reaction optimization, combined with experimental verification, in a window of operation analysis based on predefined processing tasks (for the general case, see Lima‐Ramos et al., 2014). The hybrid model was thus shown as a powerful tool of “process decision making” to efficiently customize the GG synthesis to the requirements of the envisaged production. Our study can be generally important in demonstrating a structured approach of hybrid modeling with broad applicability to enzymatic transformations. It furthermore highlights the important benefit thus achievable for biocatalytic process development.

MATERIALS AND METHODS

Enzymes and materials

Purified preparations of LmSucP and BaSucP were obtained as reported in Klimacek et al. (2020). Unless noted, the materials used were also from that article. α,α‐Trehalose (≥99%) was from Sigma‐Aldrich. PEG 35,000 was from Sigma‐Aldrich. Protein was determined with bovine serum albumin as reference.

All experiments were performed at 30°C in 50 mM MES buffer (pH 7.0), with incubations in Eppendorf tubes (1500 µl working volume) on an Eppendorf ThermoMixer C with agitation at 750 rpm.

Initial rate studies

Experiments

Previous data sets (Klimacek et al., 2020) were extended to high [Suc] (1200 mM), with reactions done in duplicate with GOH present (trans‐glycosylation) or absent (hydrolysis exclusively), as shown in Figure S1a. Enzyme was applied at 1–3 μg/ml. [Suc] or [GOH] was varied to the relevant extent of saturation, with the other substrate concentration kept constant. Samples (200 µl) were withdrawn at five times up to 40 min. Sample processing to inactivate the enzyme and analysis of the glucose (Glc) and fructose (Fru) released were done as in Klimacek et al. (2020).

Reaction in reverse direction from GG (50–400 mM, Figure S1b) was performed in the presence of phosphate (Pi; 50 mM, saturating), with incubation done for 48 h. Samples were analyzed for release of α‐d‐glucose 1‐phosphate (G1P) using a reported enzyme‐coupled assay (Wildberger et al., 2011).

Volumetric rates (μmol·L−1·s−1) determined from linear time courses of product release were converted into turnover rates (s−1) using molar enzyme concentrations calculated from the measured protein concentration and mole mass of enzyme (LmSucP: 55.750 kDa; BaSucP: 56.20 kDa) (Schwaiger et al., 2021).

Kinetic parameter determination

Nonlinear least‐squares fitting was done with OriginPro 2019b. The phenomenological Equations (1)–(4) are defined analogously to Klimacek et al. (2020). They were fitted to single sets of initial rate data. In these equations, (s−1), (s−1), and (s−1) are the release rates of Fru, GG, and Glc, respectively. Note that due to mass balance,  =  + . is the G1P release rate. We use superscripts in kinetic parameters to denote their apparent nature (maximum rate, app V; Michaelis constant, app K) and the constant substrate concentration (constant [Suc], , constant phosphate [Pi], ) used in their determination. Subscript is used to indicate the product measured in V (e.g., app V G1 P) and the varied substrate in K (e.g., app K Suc). The transfer coefficient [M−1] is the rate ratio dependent on [GOH].

Effect of nonreactive solutes

High [Suc] affects fluid micro‐viscosity, and it can impact the enzyme activity by ways unrelated to the immediate catalytic process. To mimic these effects, we used PEG‐35000 and α,α‐trehalose, respectively. Referring to literature values with suitable interpolation (Bechekh & Ghaouar, 2014; Telis et al., 2007), we reproduced the expected fluid micro‐viscosity of 800 mM Suc (2.14 Ns/m2) using PEG‐35000 at 24.44 g/L. The α,α‐trehalose was used at 780 mM. The Suc substrate was added at 20 mM in both reactions. GOH was present as indicated. Initial rates (, and ) were obtained as described above. Using reaction with Suc in the absence of GOH, we ruled out by the criterion of = that PEG‐35000 functioned as acceptor for glycosylation.

To use α,α‐trehalose in the way proposed, purity of the commercial reagent (0.83% glucose) had to be taken into account. The glucose interferes with determination of and it serves as glucosyl acceptor in the SucP reaction (Goedl & Nidetzky, 2009; Verhaeghe et al., 2016). To remove glucose, we incubated α,α‐trehalose (150 g/L; in water) with commercial baker's yeast (5 g/L) for 2 h (ambient temperature, pH 7.0, 300 rpm agitation in water). After cell removal (3650 rcf, 20 min at ambient temperature, Centrifuge 5804R, Eppendorf), the supernatant was sterile‐filtered (Minisart 0.2 µm, Sartorius), heat treated to remove protein (20 min, boiling water bath) and cleared with centrifugation (3650 rcf, 30 min at ambient temperature). The α,α‐trehalose was crystallized (Naoyuki, 1994) and dried at 60°C. The product did not contain glucose above the detection limit (0.014%, by weight) of the enzymatic assay used (K‐SUFRG; Megazyme). Using the K‐GCROL assay (Megazyme), GOH was detected in low amounts (0.15%, by weight) compatible with the intended use of the α,α‐trehalose in kinetic assays. HPLC analysis (Aminex HPX‐87C column; (Kruschitz & Nidetzky, 2020b) showed the α,α‐trehalose preparation to be at least 99% pure.

Reaction time course studies

Industrially relevant concentrations of Suc (≥300 mM) and GOH (≥1.00 M) were used (Goedl et al., 2008). Reactions were performed identically as for the initial rate measurements. A reaction scheme can be found in Figure S1c. Unless indicated, LmSucP and BaSucP were used at 3.7 U/ml. Enzyme activity was determined under initial rate conditions applying 800 mM Suc and 2 M GOH. Samples (200 µl) were taken at suitable times, processed as above and analyzed by HPLC. A Shimadzu LC‐20AD system equipped with refractive index detection was operated with a YMC‐Pack Polyamine II/S‐5um/12 nm column (YMC) at ambient temperature and a flow rate of 1 ml/min (Kruschitz & Nidetzky, 2020b). Suc, Fru, Glc, GG, and GOH were determined. The GG released was an isomeric mixture of α‐glycosylated glycerol. The 2‐O‐regioisomer was the major compound and the 1‐O‐regioisomer was formed at ~10% (LmSucP) and ~30% (BaSucP) of the total GG product (Schwaiger et al., 2021). The GG concentrations shown later are the sum of the 2‐O‐ and the 1‐O‐regioisomer. The isomeric composition of GG was independent of the conditions used and did not change dependent on the degree of substrate conversion. In all analyses performed, therefore, GG was treated as a single product. Close mass balance for substrates consumed and products formed was ensured for the data presented.

Modeling methods

Derivation of mathematical models

Reaction schemes were translated into mathematical models using the King–Altman procedure (http://www.biokin.com/tools/king-altman/index.html). Differential equations of the general form, , were thus obtained. Reaction rates were expressed as functions of the reactant concentrations (c) as well as of microscopic rate constants for the forward (+) and reverse (−) reaction ( and , with i referring to the respective reaction step). Coefficients based on microscopic rate constants in the numerator and denominator terms were grouped into kinetic parameters generally used in enzyme kinetics (Segel, 1993). Algebraic expressions describing overall reaction comprised of glucosyl transfer to GOH and hydrolysis were obtained from the derived rate equations (Tables S1 and S2).

Initial rate modeling methods

Full sets of initial rate data for reactions at varied [Suc] in the absence and presence of GOH at varied concentration were fitted with the relevant kinetic model as described under Section 3. Specifically, the dependencies of on [Suc] and [GOH], on [Suc] and on [GOH] were fitted simultaneously. The dependence of on [GOH] was excluded from the fit. Its simulation with the fitted model was used to verify the quality and internal consistency of the fit obtained. Microsoft Excel Solver add‐in was used with GRG Nonlinear Solver as the solving method, with the following set values: constraint precision of 10−6, integer optimality of 0, and convergence criterion of 10−4. The sum of relative errors squared was minimized. Constraints for the fitting were implemented as shown under Section 3. From their algebraic relationship with the microscopic rate constants estimated during fitting (Tables S1), lumped kinetic parameters from the experiment (Table S3) were used as constraints for the fit. Multiple fittings were performed to ensure that a unique solution was obtained independent of the start values for the fitted parameters. Indicated R 2 was calculated including all data shown in the respective Graph.

Modeling of time courses

Fits of integral data (concentration‐time courses) were done using the parameter estimation tool of COPASI 4.28 (Build 226) (Hoops et al., 2006). The hybrid model shown under Section 3 was used (Table S2). Reversible reaction from the products released was included. The relevant set of coupled ordinary differential equations was formulated based on mass balance (Table S4). The parameter estimation tool uses an evolutionary strategy with stochastic ranking (SERS) optimization method set to a total of 4000 generations and a population size of 60, other parameters were kept at standard settings. Parameters from initial rate analysis were used as start values for the fitting. They were additionally used to define upper and lower boundaries for parameter estimation. Constraints on rate constants from lumped kinetic parameters were applied identically as for fitting the initial rate data.

Window of operation analysis

This was performed based on simulations done with COPASI, applying the program's parameter scan tool. Simulations varied [Suc] and [GOH] in 100 mM intervals up to 2.00 M and 4.00 M, respectively. Using an enzyme loading of 11.1 U/ml, the maximum reaction time was 48 h. The simulated time courses were evaluated for substrate conversion , yield , and final [GG]. Based on processing tasks described under Section 3, an operational window was thus defined from the data.

RESULTS AND DISCUSSION

We initially performed GG synthesis under conditions of the original paper ((Goedl et al., 2008); 0.8 M Suc, 2.0 M GOH) using LmSucP or BaSucP. The idea was to apply the mechanistic‐kinetic model from Figure 1 for simulation/fitting of the reaction time courses. It must be emphasized that the model was established on initial rate data (Klimacek et al., 2020) and that it excluded microscopic rate constants of reversible product binding, that is, and (Figure 1). The model substantially deviated from the experiment, underestimating the conversion in the case of LmSucP (Figure S2a) and overestimating it in the case of BaSucP (Figure S2b). Additionally, the Fru/Glc product ratio was measured in the reaction to indicate the portion transferred of the total Suc converted. For both enzymes, as shown in Figure S2c,d, the model underestimated the product ratio by far. Moreover, in the case of BaSucP (Figure S2d), the Fru/Glc ratio showed a dependence on [Suc] conversion that was considerably larger than explicable by the decrease in [GOH] due to reaction giving GG. In an effort at reconciliation, we considered that the mechanistic‐kinetic model was developed on data acquired at relatively low [Suc] (≤ 20 mM) (Klimacek et al., 2020). We, therefore, sought to extend the model's scope to the high substrate concentrations used industrially (≥ 0.3 M).

Enzymatic reaction at high suc concentration

Reactions were done with [Suc] varied from 1.0 mM to 1.20 M with GOH absent or present at 2.00 M. Initial rates for release of Fru () and Glc (i.e., hydrolysis; ) were measured. Note, the requirement of the kinetic mechanism (Figure 1), that = when [GOH] = 0, was verified experimentally for the whole range of [Suc]. Results in Figure 2 revealed to exhibit complex “multi‐phased” dependence on [Suc]. A common feature to LmSucP and BaSucP reactions was that the did not level out at [Suc] (≥15 mM) expected to be saturating according to the apparent Suc , but continued to increase into the region of high [Suc]. The increase was largely linear, except for the BaSucP reaction at [GOH] = 2.0 M in which it was dome‐shaped. Effect of high [Suc] on was more pronounced for LmSucP (∼2.0‐fold increase) than BaSucP (∼1.5‐fold increase). The initial hydrolysis rate (), however, increased linearly with [Suc] in the absence of GOH (Figure 2). It remained constant, or even decreased, when [GOH] was increased (Figures 3b and 4b). Figures 3c and 4c reveal that for conditions of [GOH] = 2.0 M, there was a relatively larger / (LmSucP: ∼1.7‐fold; BaSucP: ∼1.2‐fold) at high (800 mM) compared with low [Suc] (20 mM). This, therefore, explains the change in the Fru/Glc ratio when during the enzymatic conversion the [Suc] decreased (Figure S2d).

Figure 2

Results of initial rate analysis extended to high [Suc]. The Fru release rates (, circles) and the hydrolysis rates (, triangles) are shown together with their fits by the first‐level hybrid model (solid lines). Panel (a) shows the results for LmSucP, panel (b) the results for BaSucP, with goodness of fit (R 2) indicated. Panel (c) shows effect of α,α‐trehalose (Tre) or PEG‐35000 on ([Suc20] = 20 mM; [GOH] = 2 M; white bars) and ([Suc20] = 20 mM; shaded bars) in comparison to [Suc800] = 800 mM; experiments were performed using LmSucP. Data are averages (N = 2) and error bars show the corresponding standard deviation [Color figure can be viewed at wileyonlinelibrary.com]

Figure 3

Results of initial rate analysis for reaction by LmSucP. Symbols show the data (circles, [Suc] = 20 mM; triangles, [Suc] = 800 mM). Solid lines show fit of the data with the first‐level hybrid model, with goodness of fit (R 2) indicated. The dependence of v H on [GOH] was not fitted and shows results of a simulation. Panel (d) compares apparent kinetic parameters determined directly (white bars) and calculated from the microscopic constants of first‐ (gray bars) and second‐level (shaded bars) hybrid model. and were determined at [Suc] of 20 mM and 800 mM, respectively. The symbol * indicates that the parameter was applied as restriction in second‐level hybrid model fitting. Data are averages (N = 2) and error bars show the corresponding standard deviation [Color figure can be viewed at wileyonlinelibrary.com]

Figure 4

Results of initial rate analysis for reaction by BaSucP. Symbols show the data (circles, [Suc] = 20 mM; triangles, [Suc] = 800 mM). Solid lines show fit of the data with the first‐level hybrid model, with goodness of fit (R 2) indicated. The dependence of on [GOH] was not fitted and shows results of simulation. Panel (d) compares apparent kinetic parameters determined directly (white bars) and calculated from the microscopic constants of first‐ (gray bars) and second‐level (shaded bars) hybrid model. and were determined at [Suc] of 20 mM and 800 mM, respectively. The symbol * indicates that the parameter was applied as restriction in second level hybrid model fitting. Data are averages (N = 2) and error bars show the corresponding standard deviation [Color figure can be viewed at wileyonlinelibrary.com]

Increase in or at high [Suc] was not accounted for by the proposed enzymatic mechanism. We speculated that the effect might be “allosterically” caused and involve Suc in a role different from that of enzyme substrate. Previous studies of Suc‐induced protein stabilization suggest that preferential exclusion of the sugar from the protein surface leads to increase in the protein's chemical potential (Kendrick et al., 1997; Kim et al., 2018; Lee & Timasheff, 1981). This in turn can favor a compactly folded protein that exposes the smallest surface area. To mimic the effect of high [Suc], we used the disaccharide α,α‐trehalose (Figure 1) which is structurally similar to Suc but not an enzyme substrate. Like Suc, α,α‐trehalose interacts with proteins via preferential exclusion at high solute concentration (Jain & Roy, 2009; Xie & Timasheff, 1997). We performed experiments with the LmSucP for the relatively larger effect of high [Suc] on this enzyme compared to BaSucP. We show in Figure 2c that α,α‐trehalose affected the LmSucP similarly overall as Suc. The ([GOH] = 2.0 M) was increased (∼1.5‐fold) compared to the control while the (no GOH added) was hardly changed. We considered that high [Suc] also affects the micro‐rheological properties of the bulk solution. Using PEG‐35,000 to modulate viscosity to approximately that of 0.8 M Suc, we show that was not increased. Evidence that was less affected by high [Suc] than indicated that decrease in water activity at high [Suc] was not a relevant factor for the phosphorylase activity. Collectively, therefore, these results suggested that an increase in at high [Suc] was caused by a general solute effect of the Suc, most likely preferential exclusion resulting in increased protein compactness. Effect on larger for the monomeric LmSucP than for the homo‐dimeric BaSucP was consistent with the idea that in enzymes of comparable subunit size (∼55 kDa), the exposed surface of the protein monomer is decreased upon oligomerization (Gunasekaran et al., 2004; Sprogøe et al., 2004). Protein structural compaction due to solute preferential exclusion is generally below detection by global spectroscopic probes of the protein conformation, as shown in several studies (Das et al., 2017; Kendrick et al., 1997; Kim et al., 2018). Development of dedicated methodology for its demonstration in the sucrose phosphorylases used here was beyond the scope of the current study.

Hybrid model for the enzymatic reaction

The hybrid model here envisaged involved empirical description of the effect of high [Suc] on microscopic rate constants of the kinetic mechanism. This “microscopic” level of integration was fundamental due to the compounded nature of the observable rate parameters. According to the proposed mechanism (Figure 1, Table S1; Klimacek et al., 2020), increases in and the transfer coefficient must arise from corresponding increase(s) in the rate constants , and . The empirical model assumed constancy of the respective rate constant at [Suc] ≤ 20 mM as well as increase linearly dependent on [Suc] at higher concentrations, as shown in Equation (5). Mathematical description in the form of a jump function, clearly, is a simplification. However, it was chosen to represent the phenomenological dependence of microscopic rate constant on [Suc] in a parameter‐economic fashion. Additionally, the meaning of the lumped kinetic parameters was maintained.

Fitting was done on the data in Figure 2 plus initial rate data acquired at [Suc] = 20 mM and 800 mM, both with [GOH] varied between 0 mM and 2.0 M. The hybrid model used is shown in Table S1. Relationships between microscopic rate constant and steady‐state parameters are also shown in Table S1. Together with steady‐state parameters known from the literature (Cerdobbel et al., 2011; Klimacek et al., 2020; Mueller & Nidetzky, 2007), these relationships were implemented as constraints for the fitting (Tables S5 and S6). The fitting results for LmSucP are shown in Figures 2a and 3; those for BaSucP in Figures 2b and 4. The estimated rate constants are summarized in Tables S5 (LmSucP) and S6 (BaSucP), with model‐calculated effects of [Suc] on the (apparent) kinetic parameters shown in Figures S3 (LmSucP) and S4 (BaSucP). For both enzymes, the slope factor () from Equation (5) that describes linear change of rate constant dependent on [Suc] was positive (increase) for k +3 and k +5 while it was negative (decrease) for k +4. Considering the increase in [Suc] from 20 mM to 800 mM, the k and the kor LmSucP increased 2.3‐ and 1.8‐fold, respectively, while the k +4 decreased 1.2‐fold. For BaSucP in the same range of [Suc], the k +3 and the k +5 increased 2.0‐ and 1.5‐fold, respectively, while the k +4 decreased 2.3‐fold. Mechanistic interpretation of these rate constant changes would require speculation not sufficiently supported by direct evidence from the current study. However, the proposed compaction of protein structure due to preferential exclusion effect of Suc at high concentration could arguably benefit the immediate catalytic steps of glycosylation of glycerol (k +3) and hydrolysis (k +5). Similarly, it would not be unreasonable for it to cause a slowing down of the GG release (k +4).

Overall, the hybrid model described the rates at low and high [Suc] quite well (Figure 2a,b). Effect of [GOH] on (Figures 3a and 4a) and (Figures 3c and 4c) was captured properly by the model, with some deviation in the dependence of on [GOH] noted. Additional evidence of the hybrid model's validity was obtained from a comparison of model‐derived and experimentally determined apparent kinetic parameters. The experimental parameters were from fits of the phenomenological Equations (1)–(3) to the data. Their algebraic equivalents (Tables S2 and S3) were calculated from estimated rate constants (Tables S5 and S6). For both enzymes, the lumped kinetic parameters (Table S3) agreed excellently with the corresponding parameters determined experimentally (Figures 3d and 4d). The simulated dependence of on [GOH] was in good accordance with the experimental results, as shown in Figure 3b for LmSucP and in Figure 4b for BaSucP. Therefore, hybrid modeling to extend the mechanistic‐kinetic model of enzymatic glycosylation of glycerol (Klimacek et al., 2020) to conditions of high [Suc] was considered to have been successful. Aggregate evidence from the initial rate analysis suggested BaSucP as the enzyme of choice for GG synthesis. In particular, about twofold higher for BaSucP than LmSucP seemed to be a difference of decisive importance.

From initial rates to full‐time courses

Denoting the analysis of initial rates as the first level of hybrid modeling, we proceeded to the evaluation of full reaction time courses as the second level. Experimental data were from reactions performed at industrially relevant substrate concentrations, using [Suc] at 800 mM ([GOH] = 2.0 M), 550 mM ([GOH] = 1.0 M) and 300 mM ([GOH] = 2.0 M). Results are shown in Figure 5 (LmSucP) and Figure 6 (BaSucP). The mechanistic‐kinetic model in Table S1 was extended by relevant terms of the reverse reaction from GG (rate constants and ). Of note, the overall reaction from GG not only proceeds to synthesize sucrose, but it also involves hydrolysis (). Release rates for the individual products as well as their rate ratios are thus affected. Equations (6) and (7) (Table 1) were derived to show effects of the reverse reaction on and respectively. Additionally, mathematical expression was derived to describe the rate ratio (Equation 8, Table 1) over the full course of conversion. Table S2 shows the relevant algebra relating the microscopic rate constants to the compounded kinetic parameters. Validity of the rate Equations (6)–(8) was confirmed numerically for the reaction conditions investigated in the experiment. Note that, due to canonical steady‐state assumptions made in their derivation (Segel, 1993), Equations (6)–(8) only hold up to the point in the conversion when the remaining substrate concentration approaches the enzyme concentration within the same order of magnitude.

Figure 5

Results of reaction time course analysis for LmSucP. [Suc]/[GOH]: 800 mM/2.00 M (a); 300 mM/2.00 M (b); and 550 mM/1.00 M (c). Symbols show the data (Suc, diamonds; Fru, squares; GG, circles; Glc, triangles), dashed lines the fit by the second‐level hybrid model and solid lines corresponding simulations, with goodness of fit (R 2) indicated. Panel (d) shows the initial rate of phosphorolysis of GG determined at 50 mM phosphate. The solid line is the fit of the data with Equation (4), with goodness of fit (R 2) indicated. Results of panel (d) are used as constraints for the fits in panel (a) [Color figure can be viewed at wileyonlinelibrary.com]

Figure 6

Results of reaction time course analysis for BaSucP. [Suc]/[GOH]: 800 mM/2.00 M (a); 300 mM/2.00 M (b); and 550 mM/1.00 M (c). Symbols show the data (Suc, diamonds; Fru, squares; GG, circles; Glc, triangles), dashed lines the fit by the second‐level hybrid model and solid lines corresponding simulations, with goodness of fit (R 2) indicated. Panel (d) shows the initial rate of phosphorolysis of GG determined at 50 mM phosphate. The solid line is the fit of the data with Equation (4), with goodness of fit (R 2) indicated. Results of panel (d) are used as constraints for the fits in panel (a) [Color figure can be viewed at wileyonlinelibrary.com]

Table 1

Equations of product release (Equations 6 and 7) and partitioning (Equation 8) between glucosyl transfer and hydrolysis derived for second‐level hybrid model

Rate equations
v_Fru=V_rV_fSucGOH−FruGGK_eq+V_rV_SucHK_GOHSucV_rK_GOHSuc+V_rK_SucGOH+K_GGV_fK_eqFru+K_FruV_fK_eqGG+V_rSucGOH+K_GGV_fK_iSucK_eqSucFru+V_fK_eqFruGG+K_SucV_rK_iGGGOHGG+K_SucHK_GOHV_r	(6)
v_H=V_rV_fHK_GOHSuc+V_fV_GGHK_FruGGK_eqV_rK_GOHSuc+V_rK_SucGOH+K_GGV_fK_eqFru+K_FruV_fK_eqGG+V_rSucGOH+K_GGV_fK_iSucK_eqSucFru+V_fK_eqFruGG+K_SucV_rK_iGGGOHGG+K_SucHK_GOHV_r	(7)
v_GGv_H=V_rV_fSucGOH−Fru[GG]K_eq−V_fV_GGHK_FruK_eq[GG]V_rV_SucHK_GOHSuc+V_fV_GGHK_FruK_eq[GG]	(8)
Kinetic parameters indicated: v, reaction rate of respective component; V_f and V_r, maximal velocity in the forward and reverse direction, index H is for hydrolysis (k+5*); K, Michaelis constants; K_i, dissociation constant, K_eq, equilibrium constant.

Note: A summary of kinetic parameters expressed with microscopic rate constants can be found in Table S2.

Based on the derived model, a set of constraints on the reverse reaction of GG was obtained from apparent kinetic parameters for GG conversion in the presence of phosphate, yielding G1P as the product. Results in Figures 5d and 6d show the relatively low reactivity of GG, with (= ∼) of 1.97 (±0.14) ⋅ 10−3 s−1 for LmSucP, in agreement with earlier work (Goedl et al., 2008), and 3.68 (±0.12) ⋅ 10−3 s−1 for BaSucP. The (Michaelis constant for GG) was 133 (±23) mM for LmSucP and 49 (±6) mM for BaSucP. Expression of in terms of microscopic rate constants is shown in Table S2.

With constraints set as indicated in Tables S7 and S8, the second‐level hybrid model was fitted to the time course data in Figures 5a (LmSucP) and 6a (BaSucP). The estimated rate constants are summarized in Tables S7 (LmSucP) and S8 (BaSucP). The fits are shown as dashed lines in the figures. Excellent description of the overall dynamics of the enzymatic conversions was obtained, as indicated by visual inspection but also from the calculated R 2 value. The time course of each reactant analyzed was captured very well (Figures 5a and 6a). Model simulations reproduced experimental results from different reaction conditions precisely (Figures 5b,c and 6b,c), with R 2 ≥ 0.80 and more typically around 0.95. To clearly see the important advance made with the second‐level hybrid model, it should be useful to compare the model predictions in Figures 5 and 6 with those in Figures S2a and S2b that were obtained with the initial rate model parametrized on data for low [Suc] (Klimacek et al., 2020).

We further showed that, in the case of LmSucP, rate constants and kinetic parameters from the time course fits (Tables S7 and S9) agreed very well with those obtainable from fits of the initial rates (Tables S5 and S9). In contrast, there was substantial variation among the corresponding rate constants for BaSucP (initial rates: Table S6; time courses: Table S7). However, kinetic parameters calculated from these different rate constants were excellently matched one with another (Table S10). This suggested that despite the differences noted, each set of rate constants for BaSucP was internally consistent. One can explain these results by fit constraints which based on known kinetic parameters were much better definable for LmSucP (Klimacek et al., 2020; Mueller & Nidetzky, 2007) than BaSucP (Cerdobbel et al., 2011). Rate constants not substantial for the fit quality (e.g., binding/release of Suc and Fru) can thus show large fluctuation, with smaller associated variation in constants partially correlated to them. Furthermore, these insensitivities may cause variations of kinetic parameters while not affecting the quality of fit, as observed in case of (the Suc release rate). Model‐calculated values of varied between 29.7 nM and 2.97 mM when altering between 10−3 s−1 and 100 s−1, although the simulated time courses did not change (Figure S5). A of 5.57 s−1 was calculated from the equilibrium constant estimated for LmSucP, which must be equal for both enzymes. These findings, therefore, strongly underline the importance of a basic set of steady‐state kinetic parameters for model parametrization via fit to reaction time course data. Note that our use of steady‐state parameters as constraints of the progress curve fit differs fundamentally from the reported idea (Chen et al., 2008) of using these parameters as mere initial guesses of the fit parameters. Importantly, as pointed out by Straathof and Heijnen (1996), the problem of parameter identifiability must be considered. Although a full‐fledged sensitivity analysis of each microscopic rate constant was beyond the scope of this study, the results (Tables S5 and S7 for LmSucP; Tables S6 and S8 for BaSucP) provide an excellent basis for the assessment of parameter sensitivity.

Enzyme comparison revealed that under all reaction conditions used (Figures 5 and 6), BaSucP released a smaller amount of glucose than LmSucP, consistent with its higher as discussed above. However, in terms of Suc conversion and GG concentration formed, LmSucP was approximately 1.5‐fold more efficient than BaSucP. For GG production optimized for conversion efficiency, therefore, LmSucP rather than BaSucP should be chosen. Interestingly, the superior performance of LmSucP was not explained from the basic kinetic mechanism but arose from the relative stronger effect of high [Suc] on LmSucP than BaSucP, as shown in Figure 2. Initial rate data and their associated parameters for both enzymes were excellently reproduced from the “reversible” hybrid model as parametrized from fits to the time course data.

Model‐based window of operation analysis for efficient optimization

Third level of the herein developed modeling approach was to use the hybrid model for targeted optimization. Our aim was to demonstrate model‐based window of operation analysis for GG production based on set reaction conditions and processing tasks. We defined a volumetric enzyme concentration of 11.1 U/ml and a reaction time of 48 h. Substrate conditions in the range [Suc] ≤ 2.0 M and [GOH] ≤ 4.0 M were screened computationally in steps of 100 mM. The conversion targets were defined as [GG] ≥ 983 mM (250 g/L), conversion of Suc () ≥ 0.98, and GG yield () ≥ 0.90. Requirements of the GG downstream processing, in particular the need to have Suc largely removed for its difficult separation from the GG, were taken into account (Kruschitz & Nidetzky, 2020a, 2020b). Results are shown in a two‐dimensional map (Figure 7), showing regions according to degree of fulfillment of the processing targets. The region in which all targets were fulfilled was larger for LmSucP than BaSucP. Point at the lower‐left corner of the region indicated the reaction conditions for properly targeted GG synthesis at minimum substrate input. This “corner point” would plausibly represent the operational optimum for production. For LmSucP, we thus find [Suc] and [GOH] at 1.1 M and 2.3 M, respectively. For BaSucP, the required [GOH] was higher (2.7 M) while [Suc] was the same (1.1 M). Predictions from the window of operation analysis were examined experimentally and the results (measured data and model predictions) are displayed in Figure 8. The shown evidence was in excellent agreement with prediction and thus verified the hybrid model for optimization purpose.

Figure 7

Results of window of operation analysis for LmSucP (a) and BaSucP (b). Marked areas correspond to reaction conditions yielding  ≥ 0.98 (I–III),  ≥ 0.90 (II–III), and [GG] ≥ 250 g/L (III). Data shown is based on time courses simulated with second‐level kinetic model. Marks show the determined operational optimum. Gray area corresponds to initial conditions the model is not defined at (negative ) [Color figure can be viewed at wileyonlinelibrary.com]

Figure 8

Experimental verification of the model‐predicted operational optimum for LmSucP (a) and BaSucP (b). The [Suc]/[GOH] was 1.10 M/2.30 M (a) and 1.10 M/2.70 M (b). Symbols show the data (Suc, diamonds; Fru, squares; GG, circles; and Glc, triangles) and lines the simulated reaction time courses by the second‐level hybrid model, with goodness of fit (R 2) indicated [Color figure can be viewed at wileyonlinelibrary.com]

CONCLUSIONS

Hybrid mathematical model for biocatalytic production of GG at industrially relevant substrate concentrations was developed and its flexible use for computational process optimization (e.g., enzyme catalyst selection; best use of substrates) was demonstrated. The hybrid model involves microscopic‐level extension of the mechanism‐based kinetic model for glycosylation of glycerol from sucrose by an empirical description of the general solute effect, likely preferential exclusion, of sucrose at high concentration. Solute‐induced compaction of the enzyme structure might explain the enhanced catalytic rates of glycosylation and hydrolysis as well as the decreased rate of GG dissociation observed at high [Suc] (≥20 mM). As demonstrated here by simulation and experimental verification, the hybrid model has immediate practical importance as an advanced tool of process engineering to customize the GG synthesis to specific requirements of industrial production. The computationally optimized production of 250 g/L GG at full conversion of Suc (1.1 M) and minimum usage of excess GOH (2.3 M) is shown.

The hybrid modeling approach from this study can be generally relevant in the field of applied bio‐catalysis (Illanes, 2008; Liese et al., 2006; Vasić‐Rački et al., 2011) in two distinct ways relating to methodology. Firstly, the evidence shown underlines the fundamental importance of mechanistic‐kinetic models for detailed engineering analysis of enzymatic transformations, those of hydrolases in synthesis in particular (Adlercreutz, 2017; Giordano et al., 2006; Vera et al., 2020). By way of its integrative structure, the hybrid model for GG production exemplifies a closing of the gap between the “microscopic” models widely used for mechanistic enzyme research and the empirical models for application in technologically relevant reaction conditions. The important benefit thus realizable for biocatalytic process development is highlighted from the current study, in particular with the window‐of‐operation analysis performed. Secondly, a step‐by‐step systematic approach to model parametrization is demonstrated based on fitting analysis with an evolving set of constraints, done on initial rate and time‐course data. In pointing out the essential role of fitting constraints for parameter identification, we take a stand on a debate about the use of initial rates or time courses (progress curves) for determining the parameters of enzyme kinetic models (Eisenkolb et al., 2020; Ohs et al., 2019; Rakels et al., 1994; Straathof, 2001; Sun et al., 2015). The question is not one of either/or. A suitable combination of both will usually be required.

Supporting information.

Click here for additional data file.