Kinetic modeling of phosphorylase-catalyzed iterative β-1,4-glycosylation for degree of polymerization-controlled synthesis of soluble cello-oligosaccharides

Background Cellodextrin phosphorylase (CdP; EC 2.4.1.49) catalyzes the iterative β-1,4-glycosylation of cellobiose using α-d-glucose 1-phosphate as the donor substrate. Cello-oligosaccharides (COS) with a degree of polymerization (DP) of up to 6 are soluble while those of larger DP self-assemble into solid cellulose material. The soluble COS have attracted considerable attention for their use as dietary fibers that offer a selective prebiotic function. An efficient synthesis of soluble COS requires good control over the DP of the products formed. A mathematical model of the iterative enzymatic glycosylation would be important to facilitate target-oriented process development. Results A detailed time-course analysis of the formation of COS products from cellobiose (25 mM, 50 mM) and α-d-glucose 1-phosphate (10–100 mM) was performed using the CdP from Clostridium cellulosi. A mechanism-based, Michaelis–Menten type mathematical model was developed to describe the kinetics of the iterative enzymatic glycosylation of cellobiose. The mechanistic model was combined with an empirical description of the DP-dependent self-assembly of the COS into insoluble cellulose. The hybrid model thus obtained was used for kinetic parameter determination from time-course fits performed with constraints derived from initial rate data. The fitted hybrid model provided excellent description of the experimental dynamics of the COS in the DP range 3–6 and also accounted for the insoluble product formation. The hybrid model was suitable to disentangle the complex relationship between the process conditions used (i.e., substrate concentration, donor/acceptor ratio, reaction time) and the reaction output obtained (i.e., yield and composition of soluble COS). Model application to a window-of-operation analysis for the synthesis of soluble COS was demonstrated on the example of a COS mixture enriched in DP 4. Conclusions The hybrid model of CdP-catalyzed iterative glycosylation is an important engineering tool to study and optimize the biocatalytic synthesis of soluble COS. The kinetic modeling approach used here can be of a general interest to be applied to other iteratively catalyzed enzymatic reactions of synthetic importance. Supplementary Information The online version contains supplementary material available at 10.1186/s13068-021-01982-2.


Background
Cellodextrin phosphorylase (CdP; EC 2.4.1.49) catalyzes the consecutive (non-processive) depolymerization of cello-oligosaccharides (COS) in the presence of phosphate, forming α-d-glucose 1-phosphate (αGlc1-P) as the cleavage product [1,2]. The COS are β-1,4-linked gluco-oligosaccharides well-known for the fact that they are released during the hydrolysis of cellulose [3]. COS are soluble in water up to a degree of polymerization (DP) of about 6 [4,5]. COS of higher DP self-assemble spontaneously in solution and thus precipitate as solid cellulose [6][7][8]. The CdP reaction is freely reversible, with the forward direction referred to as phosphorolysis and the reverse direction referred to as synthesis [1,2,9].
Due to mass action of phosphate present in excess over αGlc1-P in the cell, phosphorolysis is the preferred direction of reaction in vivo [10]. However, the CdP reaction in reverse direction under ex vivo conditions can provide an interesting route for the bottom-up synthesis of COS. Iterative β-1,4-glycosylation of cellobiose from αGlc1-P was shown in several studies [6-8, 11, 12] and it was recently demonstrated for the enzymatic production of COS at ~ 100 g/L final concentration [13]. The COS isolated from the process were shown to stimulate the growth of certain probiotic bacteria (e.g., Clostridium butyricum; Lactococcus lactis subsp. lactis), suggesting that they could be interesting dietary fibers providing a selective prebiotic effect [13]. Based on the emerging evidence on possible applications of COS for food and feed use [14], there is considerable interest in the intensification of the CdP-catalyzed conversion for the development of an efficient biocatalytic production.
The bottom-up synthesis of COS is kinetically complex, not only for the iterative glycosylation process that it involves, but also because precipitation-prone oligosaccharide products of DP ≥ 6 can be released in its course [6][7][8]. An efficient synthesis of soluble COS must therefore include control over the DP distribution of the oligosaccharide products formed, so as to avoid loss into insoluble material. It is intuitive, and has been demonstrated in previous studies, that the average DP of the COS depends on the molar ratio of cellobiose acceptor and αGlc1-P donor used in the reaction [6,8,15]. However, additional process parameters must be considered, in particular the substrate concentration and the reaction time, and there exists a complex relationship between these process parameters and the reaction output obtained (i.e., product yield, composition of the COS, insoluble portion).
We considered that a mathematical model describing the CdP-catalyzed build-up of the COS from cellobiose linked to the formation of insoluble material from the aggregation-prone portion of the COS would be fundamentally important to promote the enzymatic synthesis. Besides kinetic-mechanistic insight leading to advanced process understanding, the model could represent an important engineering tool for in silico window-ofoperation analysis and reaction optimization. This would enable rational development of a biocatalytic process designed to have tailored performance characteristics. The approach worked out here for the CdP-catalyzed synthesis of COS could be of general importance for enzymatic iterative glycosylation applied to oligosaccharide synthesis [16][17][18][19][20]. Besides CdP, other glycoside phosphorylases [16][17][18][19] as well as trans-glycosidases [21] and sugar nucleotide-dependent glycosyltransferases [22,23] are able to catalyze polymerization reactions and their corresponding products have drawn substantial interest across scientific disciplines and industrial sectors.
Development of a mechanism-based kinetic model for the CdP reaction was built upon well-established relationships between structure and function for the enzyme [24]. CdP is classified by sequence similarity as member of the glycoside hydrolase family GH94 [2,9]. CdP uses a ternary complex kinetic mechanism in which both substrates must bind to the enzyme before the catalytic reaction happens [25][26][27]. The glycosyl transfer is a single-step catalytic process involving a stereo-chemically inverting, nucleophilic substitution at the anomeric carbon of the transferred glucosyl residue [2,28]. The selectivity of CdP is strictly β-1,4 [29,30]. A crystal structure of the CdP from Clostridium thermocellum in complex with cellotetraose reveals three subsites + 1 to + 3 for binding of the glucosyl/glucose residues of the acceptor substrate [27]. The catalytic subsite -1 accommodates the glucosyl residue of αGlc1-P. Substrate binding is ordered whereby αGlc1-P binds before the disaccharide or oligosaccharide acceptor [25][26][27]. It is relevant for the kinetic model of iterative β-1,4-glycosylation from cellobiose that enzyme subsites + 1 and + 2 show much stronger interaction with the bound sugar residues than subsite + 3 does [27]. Results of previous kinetic studies are in agreement with the structural evidence, showing that in terms of k cat (turnover number) and K M (apparent Keywords: Cellodextrin phosphorylase (EC 2.4.1.49), Cello-oligosaccharides, Cellulose, Kinetic analysis and modeling, Iterative glycosylation, Bottom-up oligosaccharide synthesis Michaelis constant), cellobiose is comparable an acceptor substrate for β-1,4-glycosylation from αGlc1-P as are soluble COS of DP 3 to 5 [25,27,31].
Here, we present a detailed time-course analysis of COS synthesis from cellobiose using the CdP from C. cellulosi (CcCdP) [6] and develop a new Michaelis-Menten type mathematical model for the enzymatic conversion. We expand this model with an empirical description of the kinetics of the DP-dependent self-assembly of the COS into insoluble cellulose. A novel type of hybrid (empirical-mechanistic) model (for general case, see reference [32,33]) is thus obtained for the enzymatic polymerization process. Based on kinetic parameters determined from time-course fits, this hybrid model was shown to give an excellent description of the experimental dynamics of the COS in the DP range 3-6 and also accounts well for the insoluble product formation. We demonstrate application of the model to establish conditions for COS production that maximize the soluble product concentration at minimal loss of the product to insoluble cellulose.

Time-course analysis of the enzymatic COS synthesis
The molar substrate ratio of αGlc1-P and cellobiose determines the DP distribution in the COS products released from the CdP reaction [6,8,15]. The larger this ratio, the greater is the abundance of high-DP products (DP ≥ 5) and so the portion of total product going into insoluble material. Synthesis reactions were performed at different αGlc1-P/cellobiose ratios (0.2-4.0) to represent situations of a variable extent of COS precipitation in our experimental data set. Using HPLC method [6], the COS of DP 2-6 were well separated for quantitative analysis, as illustrated in Fig. 1. The phosphate release was additionally measured. The experimental time courses are shown in Fig. 2. Overall, there was close mass balance between substrates consumed and products released in soluble and insoluble form. In the following, we identify cellobiose as G2 and the individual COS as Gn where n indicates the DP.
Using the lowest αGlc1-P/G2 ratio of 0.2 (= 10 mM/50 mM), only G3 and G4 were formed ( Fig. 2A). Consistent with the iterative nature of the COS synthesis, the G3 was formed slightly faster than the G4. The G2 consumption paralleled the G3 formation. Increase in αGlc1-P/G2 ratio to 1.0, 2.0 and 4.0 (Fig. 2, panels B-D in that order) resulted in the formation of G5 and G6 in concentrations that increased dependent on the substrate ratio used. The G2 consumption and the phosphate release also increased, as expected. The time courses of G3 and G4 both passed through a maximum only to decrease later in the reaction. At the highest αGlc1-P/G2 ratio of 4.0, the decrease in concentration was observed even for G5. The rate and the extent of the decrease in G3 and G4 were dependent on the αGlc1-P/ G2 ratio. We analyzed the soluble and insoluble COS pools, both expressed as αGlc1-P equivalents incorporated. Comparison of the two pools revealed that the portion of insoluble product increased from effectively zero at αGlc1-P/G2 ratios of 0.2 and 1.0 to 16.1 ± 5.7 mol.% and 43 ± 7 mol.% at the higher ratios of 2.0 and 4.0, respectively (Fig. 3A).

Additional analysis of initial rates
To support the kinetic analysis of time courses with relevant boundary conditions as explained later, we determined a select set of kinetic parameters for the enzyme reaction. In particular, reaction was studied with G2 or G3 when αGlc1-P (25 mM) was constant and saturating. Initial rate data and the corresponding fitting results are shown in the Additional file 1: Figure  S1. The apparent K M for G2 was 6.0 ± 0.2 mM, that for G3 was 21 ± 3 mM. The apparent k cat was 15.9 ± 0.6 s −1 and 25.5 ± 1.3 s −1 for G2 and G3, respectively. These  x-axis and y-axis are represented by time in minutes and concentration in mM, respectively parameters of CcCdP are comparable to reported parameters for other CdP enzymes (e.g., K G2 13.2 mM and k cat_G2 47.1 s −1 for the CdP from Ruminococcus albus [25]; K G2 0.89 mM and k cat_G2 10.1 s −1 for the CdP from C. thermocellum [34]). In addition, using G2 at a saturating (100 mM) or non-saturating concentration (5 mM), we obtained kinetic parameters for αGlc1-P. Based on initial rate data shown in Additional file 1: Figure S2, the K M was determined as 0.47 ± 0.24 mM at 100 mM G2 and 0.35 ± 0.15 mM at 5 mM G2.

Mechanism-based kinetic model of the iterative CdP reaction
An ordered kinetic mechanism with αGlc1-P as the first binding substrate was assumed for each reaction of the iterative polymerization of G2 to G12 [25,27,28]. The G12 was chosen as the largest COS considering that, as shown previously, the insoluble cellulose from the CcCdP reaction has an average DP of 7-9 depending on the conditions used [6]. G6 is the largest of the soluble COS present in noticeable amount. We concluded, therefore, that the abundance of products of DP larger than 12 would have to be irrelevantly small. The CcCdP used was not active on the insoluble cellulose that precipitates from its reaction with αGlc1-P and G2, consistent with previous study stating that the precipitated oligo-cellulose is hardly accessible for further elongation by the CdP [8].
The basic reaction is shown in Eq. 1 for the conversion of G2 to G3 and is analogously used for each further step of the enzymatic polymerization.
The rate equation for the enzymatic reaction in the forward direction is shown in Eq. 2.
V is the initial reaction rate, V max is the maximum rate, K iGlc1P is the constant for binding of αGlc1-P to the free enzyme and K Glc1P and K G2 are the Michaelis constants for αGlc1-P and G2, respectively. However, Eq. 2 only would be the correct rate law to describe the conversion of G2 into G3 if there was no further conversion of the G3 into G4, and so on. To account for the iterative nature of the CdP reaction, it is necessary to develop a rate law describing the enzymatic glycosylation from αGlc1-P under conditions when multiple acceptor substrates are present at the same time. We considered that an alternative substrate (e.g., G4) affects the kinetics of an enzymatic reaction (e.g., conversion of G2 to G3) in the same way as a competitive inhibitor would affect it, except that the relevant constant describing the inhibition is not an inhibitor binding constant (K i ), but effectively the Michaelis constant (e.g., K G4 ) for the enzymatic reaction of the alternative substrate. We therefore expanded Eq. 2, by including the proper inhibition terms for the competition from all other COS that could be glycosylated from αGlc1-P. In Eq. 3, eff K G2 is an effective Michaelis constant that includes the effect from all alternative substrates present in the reaction. For the sake of simplicity, we show in Eq. 3 only the terms for the mainly formed COS. For each reaction (G2-G11), the effective Michaelis constant is defined analogously.
Each COS substrate involves its own reaction rate for glycosylation from αGlc1-P (V G2 , V G3 , V G4 …) and the rate of αGlc1-P consumption (V Glc1P ) is the sum of the individual COS rates. The rate of phosphate release (V Pi ) equals -V Glc1P . From the rate equations for the glycosylation of each COS, a set of coupled ordinary differential equations is established based on mass balance.

Effect of the reaction equilibrium
To describe the full time-course of the enzymatic conversion in solution, it was necessary to consider the chemical equilibrium for the reversible glycosylation reactions [35,36]. The equilibrium constant (K eq ) for the reaction in , with reactant concentrations being at equilibrium, and was defined analogously for all other reactions. Using the online tool eQuilibrator [37] which computes equilibria based on chemical group contribution, we obtained K eq estimates for reactions of G2, G3, G4 and G5 in a similar range (2.63-3.58). However, the K eq estimates were afflicted with errors of up to 570%, making the calculated numbers inappropriate for direct use. For further analysis including data fitting, we set a value of 3.6 as an upper boundary for the K eq .
For modeling of the reaction time courses, we chose a simple mass action term to describe the effect of the reverse reaction and so the approach to chemical equilibrium. Equation (4) was used where net V G2 is the net rate and V G2 the initial rate of G2 consumption, and Γ is the mass action ratio, expressed for the reaction in Eq. 1 as Γ = ([phosphate] t × [G3] t )/([G2] t × [αGlc1-P] t ) and for other reactions in the same way. Reactant concentrations in Γ refer to a certain time during conversion.
Combining Eqs. 2, 3 and 4, we obtained a full rate equation for each step of the iterative glycosylation. The coupled set of rate equations under the constraint of mass balance describes the polymerization process as a whole and can be used to analyze time courses in Fig. 2 that do not involve product precipitation.

Constraints in kinetic analysis
A complete summary of the constraints used is given in Table 1. Equations 5 and 6 provide relationships between kinetic parameters (K iGlc1P , K Glc1P , K G2 , K G3 ) and apparent K M values obtained from initial rate experiments described above. (3) Equation 6 was applied analogously to G3. A further constraint, dictated by the ordered mechanism of CdP [25][26][27], was that K iGlc1P and V max /K Glc1P (the catalytic efficiency for αGlc1-P) are the same for all COS substrates (G2-G11). The constraint on V max /K Glc1P was implemented by way of scaling factors (k 1 -k 5 ) that normalized the V max and K Glc1P for each COS substrate (G3-G6; k 2 -k 5 ) on the V max and K Glc1P for G2 (k 1 = 1.0). For COS substrates larger than G6, the scaling factor was assumed as k 5 . Lastly, since CdP exhibits only 3 sugar-binding subsites (+ 1-+ 3) to position the COS substrate for glycosylation from αGlc1-P [27], we assumed that the Michaelis constant for COS of DP ≥ 5 was the same as K G4 (Table 1). Preliminary fitting analysis (data not shown) revealed that the kinetic model developed to this point was able to describe well the experimental time courses of G3 and G4 for reaction at the lowest substrate ratio (αGlc1-P/ G2 = 0.2) that effectively prohibits the insoluble product formation ( Fig. 2A). We thus were encouraged to consider kinetic modeling of the COS precipitation.

Hybrid model including DP-dependent precipitation of the COS
COS solubility dependent on the DP was taken from literature (G3-G5) and extrapolated from the available data to higher DPs (G5-G12) [5]. As shown in Additional file 1: Figure S3, the solubility followed a linear trend (R 2 = 0.93) to decrease as the DP increases. For modeling, the solubility of G5-G12 was allowed to vary between boundaries of 0.72 and 1.28 times the value obtained by extrapolation ( Table 1). Precipitation of the COS was modeled as a first-order reaction, shown in Eq. 7 for G5 and used for other COS of DP ≥ 6 analogously.
The precipitation was modeled empirically as a kinetic process controlled by the precipitation rate constant. The rate constant varied between zero and its maximum value ( max k G5 ) dependent on the COS concentration in relation to the COS solubility limit, as shown in Eq. 8 for the example of k G5 In Eq. 8, [G5] SL is the G5 solubility limit and [G5] t is the soluble G5 concentration at a certain time. The parameter a relative to the solubility limit determines the steepness of the increase from zero to max k G5 . We set the a (G5) to have a value of 10 -4 , to ensure a sharp increase in the precipitation rate once 99.99% of the solubility limit was reached. k G5 is 50% of the max k G5 at the solubility limit and approaches 99.995% of max k G5 at 100.005% of the solubility limit (Additional file 1: Figure S4). To ensure the same precipitation dynamics for each COS, the parameters a and max k G had to be adjusted based on relative solubility limit. Therefore, the parameter a (Gn) was scaled from a (G5) with the ratio [Gn] SL /[G5] SL and max k G (Gn) was scaled from max k G5 with the inverse ratio [G5] SL /[Gn] SL . The max k G5 was assumed to be fast and its value set to 10 6 1/min.  The mechanistic model of CcCdP-catalyzed polymerization of COS chains in solution was thus expanded into a hybrid model that included the empirical description of COS precipitation.

Data fitting and parameter estimation
The mechanistic model of the CdP reaction involved 9 unique kinetic parameters and 1 thermodynamic parameter (K eq ) for estimation. Table 1 shows these parameters together with the associated upper and lower boundaries applied in data fitting. Initial conditions for the substrate concentrations and the CcCdP activity (V max for G2) were allowed to vary by experimental error estimated to be ± 0.5 mM and ± 0.05 mM/min, respectively ( Table 1). Data fitting was done using simultaneously all 22 time courses from Fig. 2 that represent four different synthesis experiments. Parameter estimation was done in 3 independent fittings (models PE18-PE20) and consistent results were obtained for all parameters ( Table 2). The fitted time courses are shown together with the experimental data in Fig. 2. To evaluate the overall quality of the models we created and analyzed the correlation plots (Additional file 1: Figure S5) and calculated bias and accuracy factors (Additional file 1: Table S1) [38]. The overall fit quality of model PE18-PE20 was hardly distinguishable with respect to correlation parameters (R 2 (= 0.991); slope value (= 0.991) and root mean square error (RMSE) (= 1.323-1.326)), bias factor (= 1.01) and accuracy factor (= 1.14). Obtained values of model quality parameters imply that the three models are comparable; and they represent experimental data with low bias (1.0%) and at high accuracy (predicted data are only by a factor of 1.14 different from the observed value). The correlation coefficient R 2 for the fit of the respective time course (G2-G6, phosphate) in each experiment is summarized for each model in Table 3. Again, the fit quality of model PE18-PE20 was almost identical and R 2 differed by only ≤ 0.51%. With the exception of G6 in experiments carried out at an initial αGlc1-P/G2 ratio of 2.0 (R 2 = 0.6833 ± 0.0035) and 4.0 (R 2 = 0.7556 ± 0.0014), the R 2 exceeded a value of 0.857 and generally was 0.90 or higher (Table 3). This indicated the model fit to have been very good and to cover the dynamics of each individual COS appropriately. Therefore, the models predicted equally well the soluble and insoluble COS pools accumulated in total (Fig. 3A) and over time, as can be seen in panels B and C of Fig. 3. The predicted composition of the insoluble product (G6: 44%, G7: 30%, G8: 22%, G9: 4%, G10: 1%; by mol), with an average DP of 7.0 (1145.7 Da), was in good accordance with experimental findings.
Parameter estimates are summarized in Table 2. The correlation matrix associated with these estimates is shown in Additional file 1: Table S2. The maximum reaction rate with G3 was about twofold that with G2. It was further increased about 2.7-fold and 4.6-fold when G4 and G5 was the substrate, respectively. For COS with a DP ≥ 6, the V max was assumed to be constant. The fitting results showed this V max to be enhanced not at all, or just slightly (~ 1.55-fold), as compared to the V max for G3. The K iGlc1P (αGlc1-P binding constant) and the K Glc1P (αGlc1-P Michaelis constant) were estimated at the respective upper boundaries of the fit. The K G2 was well determined whereas the K G3 and the K G4 (= K G5 -K G11 ) were all at their respective upper boundary, with K G4 in model PE20 as the sole exception ( Table 2). The reason for parameter estimates at upper boundary was identified from Additional file 1: Table S2 which shows strong statistical correlation between the estimate of K G2 and the estimates of K G3 and K G4 as well as the estimated scaling factors k 1 -k 5 . These results emphasize the importance of additional data from complementary initial rate studies (Additional file 1: Figures S1 and S2) to provide unique constraints for the time-course fitting. It is interesting that the K G2 was considerably smaller than both the K G3 and the K G4 . Therefore, occupancy of subsite + 3 (G3) and subsites + 3/ + 4 (G4) did not lead to stronger apparent binding of the acceptor oligosaccharide compared to G2. The K eq was well determined from the model fit. Its estimated value was 1.037 ± 0.004.

Model verification by simulation
To verify the hybrid model of the CcCdP reaction as obtained by fitting (Table 1, Fig. 2), we performed new synthesis experiments and used all three models to simulate the time courses of G2 consumption and product formation. The experimental conditions involved new αGlc1-P/G2 substrate ratios that lead to primarily soluble   Figure 3D, E compares the simulated solution space for the time-dependent formation of the pools of soluble (panel D) and insoluble COS (panel E) with the corresponding experimental data obtained from reaction at αGlc1-P/G2 ratio of 3.0. The soluble COS were predicted well. For the insoluble COS, the overall trend was captured correctly. It is worth noting that the simulation solution space resulted solely from the allowed variation in the initial conditions for substrate concentrations and enzyme loading. It was not due to difference in the models PE18-PE20.

Model-based analysis and optimization
The hybrid model PE18 with parameters from Table 2 was used for computational window-of-operation analysis to facilitate optimization of the enzymatic synthesis of soluble COS. As models PE18-PE20 had yielded identical results in simulations (substrate consumption, formation of soluble and insoluble COS), window-of-operation analysis with only one of the three models was considered to be sufficient. In the example chosen, the concentrations of αGlc1-P and G2 were varied each at 5-mM intervals between 10 and 100 mM. The enzyme activity was 1 U/mL and the reaction time was scanned at 2-min intervals up to 500 min. The processing objective was the maximum soluble COS concentration ([COS] max ; g/L) for each condition used. Results are shown in a condensed form in Fig. 5. They are presented fully in the Excel file in the Additional file 2. We show in Fig. 5A that the [COS] max increased with increasing ratio of αGlc1-P/G2; and that it also increased with increasing αGlc1-P concentration at constant αGlc1-P/G2 ratio.   Fig. 5B, the productivity exhibited complex dependence on the αGlc1-P/ G2 ratio. Interestingly, the productivity showed a pronounced minimum at αGlc1-P/G2 of ~ 2.0, only to increase at lower and higher substrate ratios. The productivity increases at αGlc1-P/G2 ratios greater than 2.0 involved an additional promoting effect of high αGlc1-P concentration (Fig. 5B). Figure 5 furthermore shows that the G2 conversion efficiency (panel C) increased with increasing αGlc1-P/G2 ratio whereas the αGlc1-P conversion efficiency (panel D) exhibited the opposite trend dependence. Additional observation of interest concerned the release of insoluble COS. Figure 5A shows that the [COS] max were accompanied by insoluble COS when the αGlc1-P/G2 ratio exceeded Fig. 6 Relative enrichment of G3-G6 in COS max fractions (A-D) and window-of-operation to produce G4 enriched COS max fractions (E, F) depending on the αGlc1-P/G2 ratio. Full circles and empty circles indicate soluble COS max fraction without and with co-production of insoluble COS, respectively. E, F Dashed lines correspond to data obtained at identical initial αGlc1-P concentrations. Variation of αGlc1-P and G2 concentrations within the substrate range screened (10-100 mM, with 5 mM intervals for both substrates) is indicated by arrows a value of ~ 2.0. The plots of productivity (Fig. 5B) and substrate conversion efficiency (Fig. 5C/D) also emphasize the transition to formation of insoluble product dependent on the αGlc1-P/G2 ratio. The product released at [COS] max was composed of G3 (56.3%), G4 (32.7%), G5 (9.2%) and G6 (1.7%) by weight. At maximum productivity, the COS product was composed of G3 (90.3%), G4 (9.3%) and G5 (0.4%) by weight. In summary, the window-of-operation analysis here exemplified can quickly identify reaction conditions (substrate concentration, substrate ratio, reaction time) aligned to the immediate task of the biocatalytic synthesis (e.g., avoidance of insoluble product). The maximum [COS] for each condition can be combined with the usage efficiency of αGlc1-P and G2 to identify point(s) of practical operation.
We then moved on to analyze the COS composition. Figure 6 shows for each oligosaccharide of DP3 to DP6 on the basis of mass, how the relative portion in total COS changes dependent on the αGlc1-P/G2 ratio. Whereas the G3 portion decreased continuously with increasing αGlc1-P/G2 (Fig. 6A), the G4 portion passed through a distinct maximum of ~ 0.34 at a αGlc1-P/G2 ratio of 1.2 (Fig. 6B). The G5 portion was not significant up to a αGlc1-P/G2 ratio of ~ 0.2 (Fig. 6C). It increased sharply at higher αGlc1-P/G2 ratios, reaching a portion of ~ 0.2 at αGlc1-P/G2 of 1.9-2.0. Due to the onset of insoluble product formation at αGlc1-P/G2 greater than 2.0, increase in the relative portion of G5 in total COS product was attenuated. The G6 formation started to be significant at a αGlc1-P/G2 ratio of ~ 0.5. The portion G6 increased to value of ~ 0.13 before insoluble COS are coproduced (Fig. 6D).
In a next step, we were interested in a more detailed analysis of the G4 and show the calculated concentration (Fig. 6E) and productivity (Fig. 6F) of G4 dependent on the αGlc1-P/G2 ratio. The αGlc1-P/G2 ratio giving maximum [G4] was ~ 1.0, that giving maximum productivity was ~ 0.84. Comparable to the trend for [COS] max (Fig. 5B), the productivity of G4 displayed a minimum at a αGlc1-P/G2 ratio of 2.0. From these results, conditions can be selected that tune the concentration of G4 product or the relative portion of it in the total COS released. Importantly, predicted relationships identified through window-of-operation analysis would not have been accessible without the model.

Conclusions
In a novel approach of kinetic modeling, we herein developed a hybrid kinetic model for the iterative β-1,4glycosylation of G2 from αGlc1-P catalyzed by CcCdP. The hybrid model combined a detailed mechanism-based description of kinetic and thermodynamic characteristics of the enzymatic reaction in multiple steps with an empirical description of the spontaneous (uncatalyzed) self-assembly of COS into insoluble material. The hybrid model involved 9 parameters that were estimated from reaction time-course fits performed with well-defined constraints derived from accompanying initial rates studies. Simulations showed the hybrid model to predict the time courses of further experiments, not previously used in fitting, with excellent quality. Model simulations were also used in a window-of-operation analysis to study and optimize the biocatalytic synthesis of COS. Key process parameters (i.e., yield, COS concentration, composition of soluble COS) were predicted in dependence of the main process variables (substrate and enzyme concentrations, αGlc1-P/G2 ratio, reaction time). It was shown that a αGlc1-P/G2 ratio of ~ 2.0 marked the transition in the reaction to insoluble product formation. This ratio was also the point of lowest productivity of the soluble COS.
The mechanistic part of the hybrid model can be relevant to other enzymatic reactions used to synthesize product oligomers via iterative polymerization. Various phosphorylases (e.g., α-glucan phosphorylase [19], laminarin phosphorylase [16,17], β-glucan phosphorylase [18]) catalyze the polymerization of a primer substrate from αGlc1-P. Sugar nucleotide-dependent glycosyltransferases [22,23] catalyze the polymerization of neutral and acidic oligosaccharides. These oligosaccharides have considerable interest for technological application in different industrial sectors, including food and feed, cosmetics and health-care. The current study provides an important resource of methodology for the kinetic analysis of the enzymatic polymerization; and it shows use of the model for reaction analysis and optimization.

Enzyme
N-terminally His-tagged CcCdP (GenBank identifier CDZ24361.1) expressed and purified according to reported method [6] was used in all experiments. Enzyme stock solution (~ 20 mg protein/mL) was made in 50 mM MES buffer (pH 7.0). It was stored at − 20 °C for several weeks, and could be thawed and frozen repeatedly, without significant loss of activity (data not shown). Enzyme purity was verified by SDS PAGE and determination of specific activity (13.3 ± 0.4 U/mg purified protein).