Development of a Digital Twin for Enzymatic Hydrolysis Processes

Abstract

Enzymatic hydrolysis processes can be used to produce organic nutrient media from renewable raw materials. However, many of these processes are not optimally designed, so expensive enzymes and substrates are wasted. Mathematical models and Digital Twins (DTs) are powerful tools, which can be used to optimize bioprocesses and, thus, increase the yield of the desired products. Individual enzymatic hydrolysis processes have already been modeled, but models for the combined starch hydrolysis and proteolysis, or DTs, are not available yet. Therefore, an easily adaptable, dynamic, and mechanistic mathematical model representing the kinetics of the enzymatic hydrolysis process of the combined starch hydrolysis and proteolysis was developed and parameterized using experimental data. The model can simulate the starch hydrolysis process with an agreement of over 90% and the proteolysis process with an agreement of over 85%. Subsequently, this model was implemented into an existing DT of a 20 L stirred tank reactor (STR). Since the DT cannot only map the kinetics of the enzymatic process, but also the STR with the associated periphery (pumps, heating jacket, etc.), it is ideally suited for future process control strategy development and thus for the optimization of enzymatic hydrolysis processes.

1. Introduction

Enzymatic hydrolysis processes, such as starch hydrolysis and proteolysis, play a key role in the production of organic nutrient media, as no toxic chemicals or other additives are needed in these processes. Renewable raw materials containing starch or proteins are split into their main components (e.g., glucose, free amino acids (fAA)) by the enzymes and can then be used by various microorganisms as the basis for their growth and product formation. However, there is still room for optimization of many enzymatic hydrolysis processes due to high enzyme consumption or low substrate conversion. Seasonal fluctuations in the quality of the organic raw materials may also make it necessary to regularly adjust the process conditions and control to achieve the optimal result.

Mathematical models can be used to support the optimization of bioprocesses. Usually, these models represent the micro-kinetics of the reactions under consideration. However, it is beneficial to integrate these kinetic models into a Digital Twin (DT). In many application areas, DTs are becoming more important. Especially in “Industry 4.0”, DTs are becoming increasingly valuable and their application is being studied frequently [1,2,3]. DTs also may represent the reactor macro-kinetics and the dynamics of the system. They are thus ideally suited for systematic optimization of dynamic processes [4,5].

DTs are digital representatives of material or immaterial objects or processes. It is irrelevant whether the counterpart already exists in the real world or will exist in the future. They consist of models of the represented object or process and can also contain algorithms and services that describe or influence the properties or behavior of the represented object or process. The coupling between the real process and the associated DT can take place directly on the running process (process accompanying) or cyclically after the execution of the real process [6,7,8,9,10]. DTs for bioprocesses may be designed, using a shell structure [11,12] (Figure 1), including a biological, a physicochemical, a reactor, a plant and periphery as well as a control and automation sub-model. To facilitate the operation of DTs, they can be equipped with graphical user interfaces (GUIs). These GUIs are designed similarly to the user interfaces of the reactors’ process control systems.

Within the scope of this work, a DT for enzymatic hydrolysis processes was developed. To achieve this, a generic, dynamic, and mechanistic mathematical model for the combined enzymatic starch hydrolysis and proteolysis was developed and integrated with an existing stirred tank reactor (STR) model. During the enzymatic starch hydrolysis, α-amylase splits starch into oligosaccharides, which are subsequently converted into glucose by glucoamylases. In the first step of the enzymatic proteolysis, proteins are converted into peptides by endopeptidases; these peptides are subsequently split into fAA using exopeptidases. Models already exist for the individual hydrolysis processes [13,14,15], but a generic model for the combined process has not yet been developed. The model was designed in such a way that it can be easily adapted to changes in enzyme or substrate properties.

Subsequently, the developed mathematical model for the combined starch hydrolysis and proteolysis was integrated into an existing DT developed by our group [3,16,17,18,19,20,21,22]. This DT has been shown to simulate the cultivation of S. cerevisiae and the whole-cell biocatalysis of ethyl-3-hydroxy-butyrate (E3HB) from the substrate ethyl acetoacetate in a 20 L STR under aerobic and anaerobic conditions [23]. Instead of developing a completely new DT, this DT was used as a basis, because of its ability to represent the characteristics of the STR in which the enzymatic hydrolysis processes are carried out and its associated equipment. Furthermore, once the new process model has been integrated, the entire production cycle, from nutrient media production and cultivation (S. cerevisiae) to the generation of the products ethanol and E3HB can be simulated.

The new DT is intended to provide a tool for the systematic optimization of enzymatic hydrolysis processes.

The materials and methods used for modeling the DT for enzymatic hydrolysis processes, as well as for the experimental realization and the analysis of the enzymatic starch hydrolysis and proteolysis are presented in Section 2. The developed dynamic, mechanistic, mathematical model for the enzymatic hydrolysis processes (combined enzymatic starch hydrolysis and proteolysis) is presented in Section 3 and it is shown how it was combined with an already existing DT. In addition, the comparison between simulation results with the new DT and real experiments is made. Furthermore, it is shown how control strategies can be developed with the new DT for enzymatic hydrolysis processes. Section 4 concludes with a discussion of the results and an outlook.

2. Materials and Methods

This section describes how the dynamic and mechanistic mathematical model of the combined enzymatic starch hydrolysis and proteolysis was created. This model was then implemented into an existing DT representing a 20 L STR. For the generation and validation of the model, experiments of the enzymatic starch hydrolysis and proteolysis were carried out at a small scale (5–10 mL) as well as in STRs (6 L and 20 L).

2.1. Development of the Enzymatic Hydrolysis Process Model

The dynamic and mechanistic mathematical model for the enzymatic hydrolysis processes was developed, applying the modeling cycle [12] shown in Figure 2.

In the first step a dynamic and mechanistic mathematical model [12] representing the kinetics of the enzymatic hydrolysis processes (combined starch hydrolysis and proteolysis) derived from appropriate process descriptions and the knowledge gained from small scale characterization experiments, was developed in the C++-based modeling environment C-eStIM [24]. Subsequently, the model was first parameterized using the experimental data from the characterization experiments.

The model was then used to predict further batch experiments in STRs. After conducting and evaluating these batch experiments for starch hydrolysis and proteolysis, the measured data were compared to the data simulated by the model. Finally, the model was integrated into the existing DT for a 20 L STR with peripheral equipment, mentioned before.

To determine the accuracy of the fit of the model after its parametrization, the coefficient of determination (R²) was calculated by dividing the difference between experimental y_i and simulated data y_s,i with the difference between the experimental data and their mean value ӯ [25] (Equation (1)).

$$R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-y_{s,i}\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(1)

R² can take on values between minus infinity and one. An R² close to one signifies a good fit of the model to the experimental data. If R² is less than zero, the mean of the measured data points is closer to the mean value than the simulated results.

Initially, the mathematical model was written in C-eStIM, and its parameters were adjusted manually to approximate the measured data. For further parameterization, the model parameters were adapted by minimizing the weighted mean-square deviation (wMSD) using the Nelder–Mead algorithm [25,26]. The Nelder–Mead algorithm was selected because the method is effective and computationally efficient [26]. In addition, the programming environments C-eStIM [24] and R [27] used for this work, offer special Nelder–Mead packages. The wMSD was calculated from the squared difference between the measured value y_m and the simulated value y_s, divided by the number of data points n in the data set and multiplied by a factor for weighting individual data points k_weighting [25] (Equation (2)).

$$wMSD={\displaystyle \sum }_{i=1}^n\frac{{\left(y_{s,i}-y_{m,i}\right)}^2}{n}·k_{weighting}$$

(2)

2.2. Integration of the Mathematical Model for the Combined Starch Hydrolysis and Proteolysis into an Existing DT of a 20 L STR

The DT, which was presented in the introduction, consists of four sub-models representing the biological and physical-chemical processes in the reactor and the interactions with associated equipment (plant and periphery). The sub-models were written in the C++-based modeling and simulation environment C-eStIM [24]. The C++ models were implemented in the modular process control and automation system WinErs [28] using C-eStIM-DLL interfaces. In addition, the DT is equipped with a control and automation sub-model as well as the GUI, also created using WinErs.

In the biological sub-model of the DT, the kinetics for the growth of S. cerevisiae, as well as the kinetics of the whole-cell biocatalysis of E3HB are described. The developed mathematical model of the enzymatic hydrolysis processes was embedded in the source code of the biological sub-model. Subsequently, the adapted C-eStIM [24] model was implemented in WinErs [28] using a C-eStIM-DLL interface.

2.3. Enzymatic Hydrolysis Processes (Starch Hydrolysis and Proteolysis)

The experiments of the enzymatic hydrolysis processes (starch hydrolysis and proteolysis) were carried out in small-scale test tubes (V = 1–5 mL) as well as in 6–20 L STRs.

2.3.1. Starch Hydrolysis

For the characterization of the enzymes used in the starch hydrolysis process, 5 mL buffer (pH 2–4 (phosphate citrate buffer), pH 5–8 (phosphate buffer, pH 9–10 (Tris-HCl buffer)), containing 20 g L⁻¹ soluble potato starch (substrate, Carl Roth, Karlsruhe, Germany), 0.05 mg L⁻¹ Termamyl SC (α-amylase, Novozymes, Bagsværd, Denmark) or 0.01 mg L⁻¹ Spirizyme Ultra (glucoamylase, Novozymes, Bagsværd, Denmark) and 0.02 g L⁻¹ calcium chloride dihydrate (CaCl₂·2H₂O, Merck, Darmstadt, Germany) was added to a 15 mL test tube (VWR, Darmstadt, Germany), mixed and placed in a water bath (T = 0–100 °C). After 10 min, 90 µL of sample was drawn and the reaction was stopped by adding 600 µL of 5 molar hydrochloric acid and heating the sample for 5 min at 80 °C in a heating block (VWR). For the investigation of the pH dependency, the temperature was set to 60 °C. For investigation of the temperature dependency, the pH was set to 5.

The batch starch hydrolysis experiment for the parameterization of the model was performed in a 20 L STR (Biostat C, B. Braun, Melsungen, Germany). Therefore, 13.5 L phosphate citrate buffer (pH = 5), with 160 g L⁻¹ soluble potato starch (substrate, Carl Roth, Karlsruhe, Germany), 0.006 g L⁻¹ Termamyl SC (α-amylase, Novozymes, Bagsværd, Denmark), 0.4 g L⁻¹ Spirizyme Ultra (glucoamylase, Novozymes, Bagsværd, Denmark), and 0.02 g L⁻¹ calcium chloride dihydrate (CaCl₂·2H₂O, Merck, Darmstadt, Germany) was added to the reactor. The temperature was set to 60 °C and the stirrer speed to 300 rpm. Every 5–10 min, 90 µL of sample were drawn from the reactor. The reaction was stopped by adding 600 µL 5 molar hydrochloric acid, and heating the sample for 5 min at 80 °C in a heating block (VWR, Darmstadt, Germany).

Starch concentrations were determined using the potassium iodide method [29]. A potassium iodide solution was prepared, containing 20 g L⁻¹ potassium iodide (Carl Roth, Karlsruhe, Germany) and 2 g L⁻¹ iodine (Riedel-de Haën, Honeywell, Seelze, Germany), dissolved in H₂O. Moreover, 50 µL of potassium iodide solution was added to each sample and hydrochloric acid mixture (690 µL). After that, the absorbance of each sample was measured in a UV–visible spectrophotometer (UV mini 1240, Shimadzu, Nakagyo-ku, Kyoto, Japan) at 623 nm. A calibration line was determined to calculate the concentration of starch in the samples, using standards with a defined amount of soluble potato starch (Carl Roth, Karlsruhe, Germany).

The concentrations of glucose (product), maltose, and sucrose (intermediate products) were determined using an enzymatic assay kit (Maltose/Sucrose/D-Glucose, UV-method, R-Biopharm, Darmstadt, Germany).

2.3.2. Proteolysis

For the characterization of the enzymes used in the proteolysis process, 1 mL buffer (pH 4 (phosphate citrate buffer), pH 5–8 (phosphate buffer, pH 9 (Tris-HCl buffer)) with 2 g L⁻¹ casein powder (substrate, Carl Roth, Karlsruhe, Germany), 0.2 mg L⁻¹ EnerZyme^® P7 (endopeptidase, Erbslöh, Geisenheim, Germany) or 0.2 mg L⁻¹ Flavourzyme™ (exopeptidase, Novozymes, Bagsværd, Denmark) was added to a 2 mL test tube (VWR, Darmstadt, Germany), mixed and placed in a water bath (T = 25, 35, 45, 55, 65, 70 °C). After 30 min, 100 µL of sample was drawn and the reaction was stopped by heating the sample for 5 min at 95 °C in a heating block (VWR, Darmstadt, Germany). For the investigation of the pH dependency, the temperature was set to 50 °C. For investigation of the temperature dependency, the pH was set to 7.

The batch proteolysis experiment for the parameterization of the model was performed in a 6 L STR (BioFlo 3000, New Brunswick Scientific, Eppendorf, Hamburg, Germany). Therefore, 3 L phosphate buffer (pH 7.5) with 40 g L⁻¹ organic sunflower seed meal (substrate, All Organic Treasures, Wiggensbach, Germany) containing 14 g L⁻¹ protein, 0.2 g L⁻¹ EnerZyme^® P7 (endopeptidase, Erbslöh, Geisenheim, Germany) and 0.2 g L⁻¹ Flavourzyme™ (exopeptidase, Novozymes, Bagsværd, Denmark) were added to the STR. The temperature was set to 50 °C and the stirrer speed to 300 rpm. Every 15 min, samples (2 mL) were drawn from the reactor. The reaction was stopped by heating the sample for 5 min at 95 °C in a heating block (VWR, Darmstadt, Germany).

The fAA (product) concentration was determined using the ninhydrin method [30]. The ninhydrin solution, prepared, was a 1:1 mixture of a solution containing 1.6 g L⁻¹ tin chloride dihydrate (Carl Roth, Karlsruhe, Germany) in citrate-phosphate buffer (pH 5.2) and a solution containing 40 g L⁻¹ ninhydrin (Carl Roth, Karlsruhe, Germany) in ethylene glycol (Carl Roth, Karlsruhe, Germany). 1 mL of ninhydrin solution was added to 100 µL of sample solution and incubated at 80 °C for 20 min in a water bath (WBT 12, Carl Roth, Karlsruhe, Germany). After the addition of 5 mL of isopropyl alcohol (Carl Roth, Karlsruhe, Germany), the absorbance of the mixture was measured in a UV–visible spectrophotometer (UV mini 1240, Shimadzu, Nakagyo-ku, Kyoto, Japan) at 570 nm. Using standards with defined L-glutamic acid (Fluka, Honeywell, Seelze, Germany) concentrations, a calibration curve was determined to calculate the concentration of fAA in the samples.

3. Results

In the subsequent section, brief process descriptions of the enzymatic hydrolysis processes (starch hydrolysis and proteolysis), from which the corresponding mathematical model was derived are presented. In addition, data from experimental studies concerning the reaction rate dependency on temperature and pH value are shown, serving as a basis for the calibration of the corresponding part of the kinetic model. For the parameterization of the model, reactor experiments were pre-simulated and then carried out. After successful reparameterization of the model, the model was embedded into the existing DT. Finally, a simulation study on the interaction between temperature and pro-duct concentrations of starch hydrolysis in a continuously operated STR using the DT is shown.

3.1. Process Description of the Enzymatic Hydrolysis Processes (Starch Hydrolysis and Proteolysis)

In starch hydrolysis, starch is converted into glucose molecules using α- and glucoamylases. α-amylase (α-1,4 glucan glucanohydrolase, EC. 3.2.1.1) catalyzes the hydrolysis of 1,4-α-D-glucosidic linkages in polysaccharides, containing three or more 1,4-α-linked D-glucose units randomly. Starch is degraded to low molecular weight dextrin, limited by α-1,6 bonds in a way that a-amylase is unable to hydrolyze. Glucoamylase (1,4-α-D-glucan glucohydrolase, EC 3.2.1.3) is a multi-domain exo-glycosidase that catalyzes the hydrolysis of α-1,4 and α-1,6 glucosidic linkages of starch and related polysaccharides to release β-D-glucose from the non-reducing ends [31,32,33].

During proteolysis, proteins are split into fAA and di- and tripeptides using endo- and exopeptidases. Peptidases split proteins by hydrolysis of the peptide bonds. Endopeptidases enzymatically split peptide bonds within the protein. In contrast to the endopeptidases, the exopeptidases only split terminal peptide bonds of the amino acid sequence [34].

For starch hydrolysis and proteolysis, reaction rates can be calculated using Michaelis–Menten kinetics, with maximal reaction rates (r_max), half-saturation constants (K_S) and the substrate concentration (c_substrate) (Equation (3)).

$$r=\frac{r_{max}·c_{substrate}}{K_S+c_{substrate}}$$

(3)

Michaelis–Menten kinetics was chosen because it has already been successfully applied by other research groups for the modeling of the enzymatic hydrolysis processes under investigation [13,14,15,35]. In addition, the relatively simple structure of the kinetics offers rapid adaptability, which supports the generic approach of the model.

3.2. Mathematical Modeling of the Enzymatic Hydrolysis Process Combining Starch Hydrolysis and Proteolysis

The structure of the mathematical model is illustrated in Figure 3.

It is assumed that the substrates for starch hydrolysis and proteolysis are composed of hydrolysable and non-hydrolysable components. Depending on the substrate used, the ratio must be set before the simulation experiment. The yield coefficients (Y_IPS,E1, Y_PS,E1, Y_PIP,E2) used in the model were based on assumptions derived from literature data [31,32,33,34].

During starch hydrolysis, 95% of the hydrolysable components of the substrate are converted by α-amylase into the intermediate product and 5% directly into glucose. The intermediate product formed is then converted by glucoamylase to 100% into glucose. Similarly, in proteolysis, 95% of the hydrolysable components of the substrate are converted by endopeptidase into the intermediate product and 5% directly into fAA. The intermediate product formed is then converted by exopeptidase to 100% into fAA.

The model equations for starch hydrolysis and proteolysis have a similar structure and are indicated by suffixes P1 (starch hydrolysis) and P2 (proteolysis). In the interest of simplicity, the general equations are presented below without the suffixes (Equations (4)–(20). The abbreviations used in the model are shown in

Table 1. Abbreviations used in the enzymatic hydrolysis processes model.

Abbreviations	Description
rSIP,E1	Degradation rate of substrate to intermediate product.
rSP,E1	Degradation rate of substrate to product.
rIPP,E2	Degradation rate of intermediate product to product.
rmax,SIP,E1	Maximum degradation rate of substrate to intermediate product.
rmax,SP,E1	Maximum degradation rate of substrate to product.
rmax,IPP,E2	Maximum degradation rate of intermediate product to product.
KM,SIP,E1	Half-saturation constant for degradation of substrate to intermediate product.
KM,SP,E1	Half-saturation constant for degradation of substrate to product.
KM,IPP,E2	Half-saturation constant for degradation of intermediate product to product.
fT,act,E1	Factor for temperature-dependent activity of enzyme 1.
fT,act,E2	Factor for temperature-dependent activity of enzyme 2.
fpH,act,E1	Factor for pH-dependent activity of enzyme 1.
fpH,act,E2	Factor for pH-dependent activity of enzyme 2.
rIPS,E1	Conversation rate of intermediate product from substrate.
rPS,E1	Conversation rate of product from substrate.
rPIP,E2	Conversation rate of product from intermediate product.
YIPS,E1	Yield coefficient for conversation of intermediate product from substrate.
YPS,E1	Yield coefficient for conversation of product from substrate.
YPIP,E2	Yield coefficient for conversation of product from intermediate product.
rden,E1	Denaturation rate of enzyme 1.
rden,E2	Denaturation rate of enzyme 2.
rmax,den,E1	Maximum denaturation rate of enzyme 1.
rmax,den,E2	Maximum denaturation rate of enzyme 2.
fT,sta,E1	Factor for temperature-dependent stability of enzyme 1.
fT,sta,E2	Factor for temperature-dependent stability of enzyme 2.
fpH,sta,E1	Factor for pH-dependent stability of enzyme 1.
fpH,sta,E2	Factor for pH-dependent stability of enzyme 2.
VB	Volume of broth in the reactor.
S1	Concentration of starch in the reactor.
SH1	Concentration of the hydrolysable components of the substrate in the reactor.
SNH1	Concentration of the non-hydrolysable components of the substrate in the reactor.
E11	Concentration of enzyme 1 in the reactor.
E21	Concentration of enzyme 2 in the reactor.
IP1	Concentration of the intermediate product in the reactor.
P1	Concentration of the product in the reactor.
FS,0	Inflow of substrate solution to the reactor.
FE1,0	Inflow of enzyme 1 to the reactor.
FE2,0	Inflow of enzyme 2 to the reactor.
SH0	Concentration of the hydrolysable components of the substrate in FS,0.
SNH0	Concentration of the non-hydrolysable components of the substrate in FS,0.
E10	Concentration of enzyme 1 in FE1,0.
E20	Concentration of enzyme 2 in FE2,0.
IP0	Concentration of the intermediate product in FS,0.
P0	Concentration of the product in FS,0.
Fin	Absolute inflow into the reactor.
Fout	Absolute outflow from the reactor.

.

Michaelis–Menten kinetics were implemented in the model for the calculation of the reaction rates, where the maximum reaction rates were modulated with functions describing the influence of temperature and pH (Equations (4)–(6)). Equations (7)–(9) representing the conversation rates. The denaturation rates of the enzymes are calculated using Equations (10) and (11).

$$r_{SIP,E1}=\frac{r_{max,SIP,E1}·(f_{T,act,E1 }·f_{pH,act,E1 })· S{H}_1}{K{s}_{SIP,E1}+S{H}_1}$$

(4)

$$r_{SP,E1}=\frac{r_{max,SP,E1}·(f_{T,act,E1 }·f_{pH,act,E1 })· S{H}_1}{K{s}_{SP,E1}+S{H}_1}$$

(5)

$$r_{IPP,E2}=\frac{r_{max,IPP,E2}·\left(f_{T,act,E2 }·f_{pH,act,E2 }\right)· I{P}_1}{K{s}_{IPP,E2}+I{P}_1}$$

(6)

$$r_{IPS,E1}=Y_{IPS,E1}·r_{SIP,E1}$$

(7)

$$r_{PS,E1}=Y_{PS,E1}·r_{SP,E1}$$

(8)

$$r_{PIP,E2}=Y_{PIP,E2}·r_{IPP,E2}$$

(9)

$$r_{den,E1}=r_{max, den,E1}·f_{T,sta,E1}·f_{pH,sta,E1}$$

(10)

$$r_{den,E2}=r_{max, den,E2}·f_{T,sta,E2}·f_{pH,sta,E2}$$

(11)

The differential equations in the model were created using the general dynamic mass balance (Equation (12)).

$$accumulation=input-output+generation –consumption$$

(12)

The reactor model of an ideally mixed STR consists of seven differential equations (Equations (13)–(19)). The influence of temperature and pressure on the differential equations is neglected in the model. In addition, the overall substrate concentration is given by Equation (20).

$$\frac{d{V}_B}{dt}=F_{in}-F_{out}$$

(13)

$$\frac{dS{H}_{1 }}{dt}=\frac{F_{S,0}·S{H}_0}{V_B}-\frac{S{H}_1·\left(F_{in}-F_{out}\right)}{V_B}-r_{SIP,E1}·E{1}_1-r_{SP,E1}·E{1}_1$$

(14)

$$\frac{dE{1}_1}{dt}=\frac{F_{E1,0}·E{1}_0}{V_B}-\frac{E{1}_1·\left(F_{in}-F_{out}\right)}{V_B}-r_{den,E1}·E{1}_1$$

(15)

$$\frac{dE{2}_1}{dt}=\frac{F_{E1,0}·E{2}_0}{V_B}-\frac{E{2}_1·\left(F_{in}-F_{out}\right)}{V_B}-r_{den,E2}·E{2}_1$$

(16)

$$\frac{dI{P}_{1 }}{dt}=\frac{F_{S,0}·I{P}_0}{V_B}-\frac{I{P}_1·\left(F_{in}-F_{out}\right)}{V_B}+r_{IPS,E1}·E{1}_1-r_{IPP,E2}·E{2}_1$$

(17)

$$\frac{d{P}_{1 }}{dt}=\frac{F_{S,0}·P_0}{V_B}-\frac{P_1·\left(F_{in}-F_{out}\right)}{V_B}+r_{PIP,E2}·E{2}_1+r_{PS,E1}·E{1}_1$$

(18)

$$\frac{dSN{H}_{1 }}{dt}=\frac{F_{S,0}·SN{H}_0}{V_B}-\frac{SN{H}_1·\left(F_{in}-F_{out}\right)}{V_B}$$

(19)

$$S_1=S{H}_1+SN{H}_1$$

(20)

3.3. Temperature- and pH-Dependency of the Enzymes Used in Starch Hydrolysis and Proteolysis

The results of the enzyme characterization experiments were normalized to values between 0 and 1 for the individual temperatures or pH values. A value of 1 corresponds to the temperature or pH value at which the highest concentration of product was formed. Furthermore, double sigmoidal (DSig) equations (Equation (21)) [36,37] were fitted to the experimental data.

$$f_{DSig}\left(x\right)=\left(Y_{LS}+\frac{Y_{mid}-Y_{LS}}{1+e^{-K_{LS}·\left(x-r_{max,low}\right)}}\right)·\left(1+\frac{\left(\frac{Y_{RS}}{Y_{mid}}-1\right)}{1+e^{-K_{RS}·(x-r_{max,high)}}}\right)$$

(21)

The value of a state variable is described by x. Y_LS is the value at low x, Y_RS is the value at high x, Y_mid is the value between r_max,low and r_max,high, which are location parameters of the low/high side of the function, K_LS determines the slope on the low side and K_RS determines the slope on the high side of the function.

These DSig equations are used in the mathematical model of the enzymatic hydrolysis processes to fit the temperature and pH dependency of the enzymes. In the model, the factors (f_DSig(x)) resulting from the double sigmoidal equations are multiplied by the maximum reaction rates (r_max).

Figure 4 shows the experimentally determined temperature and pH dependency of the α-amylase and glucoamylase used in the starch hydrolysis process, as well as the DSig equations adapted to it (temperature (a) and pH (b)).

For the α-amylase and the glucoamylase, almost no activity could be determined at temperatures below 20 °C. The activity of the α-amylase increased at the highest rate between 40 and 70 °C. The maximal activity was reached at 100 °C. Similarly, the highest increase in activity of the glucoamylase was found between 40 and 60 °C. The maximum activity was reached at 70 °C and above. The results correspond to the manufacturer’s specifications for these thermostable enzymes [38,39]. At temperatures between 40 and 80 °C, varying temperatures lead to clear changes of enzymatic activity. Here, good temperature control is required to achieve a stable process performance.

The α-amylase shows almost no activity at pH-values below 3 and above 8. The maximum activity was reached at a pH of around 5. For the glucoamylase, almost no activity could be determined at pH-values below 3 and above 10. The highest activity was found in a pH range between 5 and 7. At pH values between 3 and 8 (10), varying pH values lead to clear changes of enzymatic activity. Thus, in this range, good pH control is required to achieve a stable process performance.

The DSig functions could be fitted to the experimentally determined temperature and pH dependency of the glucoamylase and α-amylase investigated, by setting the parameters of the DSig equations as shown in

Table 2. Parameters of the DSig equations for mapping the temperature (T) and pH dependency of the α-amylase (E1) and glucoamylase (E2) used in starch hydrolysis (P1).

Parameter	Description	Value
YLS,T,act,E1,P1	DSig value at low temperature.	0.0
Ymid,T,act,E1,P1	DSig value between rmax,Tlow,E1,P1 and rmaxT,high,E1,P1.	1.0
YRS,T,act,E1,P1	DSig value at high temperature.	1.0
KLS,T,act,E1,P1	Slope on the low side of the DSig function.	0.1
KRS,T,act,E1,P1	Slope on the high side of the DSig function.	0.1
rmax,Tlow,act,E1,P1	Location parameter of the low side of the DSig function.	56.0 °C
rmax,Thigh,act,E1,P1	Location parameter of the high side of the DSig function.	56.0 °C
YLS,pH,act,E1,P1	DSig value at low pH.	0.0
Ymid,pH,act,E1,P1	DSig value between rmax,pHlow,act,E1,P1 and rmax,pHhigh,act,E1,P1.	1.4
YRS,pH,act,E1,P1	DSig value at high pH.	0.0
KLS,pH,act,E1,P1	Slope on the low side of the DSig function.	3.5
KRS,pH,act,E1,P1	Slope on the high side of the DSig function.	1.6
rmax,pHlow,act,E1,P1	Location parameter of the low side of the DSig function.	4.5
rmax,pHhigh,act,E1,P1	Location parameter of the high side of the DSig function.	6.0
YLS,T,act,E2,P1	DSig value at low temperature.	0.0
Ymid,T,act,E2,P1	DSig value between rmax,Tlow,E2,P1 and rmaxT,high,E2,P1.	1.0
YRS,T,act,E2,P1	DSig value at high temperature.	1.0
KLS,T,act,E2,P1	Slope on the low side of the DSig function.	0.1
KRS,T,act,E2,P1	Slope on the high side of the DSig function.	0.1
rmax,Tlow,act,E2,P1	Location parameter of the low side of the DSig function.	45.0 °C
rmax,Thigh,act,E2,P1	Location parameter of the high side of the DSig function.	45.0 °C
YLS,pH,act,E2,P1	DSig value at low pH.	0.0
Ymid,pH,act,E2,P1	DSig value between rmax,pHlow,act,E2,P1 and rmax,pHhigh,act,E2,P1.	1.0
YRS,pH,act,E2,P1	DSig value at high pH.	0.0
KLS,pH,act,E2,P1	Slope on the low side of the DSig function.	3.0
KRS,pH,act,E2,P1	Slope on the high side of the DSig function.	3.0
rmax,pHlow,act,E2,P1	Location parameter of the low side of the DSig function.	4.0
rmax,pHhigh,act,E2,P1	Location parameter of the high side of the DSig function.	8.5

.

Figure 5 shows the experimentally determined temperature and pH dependency of the relative activity of the endo- and exopeptidase used in the proteolysis process, as well as the DSig equations adapted to it (temperature (a) and pH (b)).

For the endopeptidase and exopeptidase, only low activities could be determined at temperatures below 30 °C and above 70 °C and at pH values below 5 and above 10. The endopeptidase showed the highest activity at a temperature of around 55 °C and a pH of about 7. Similarly, the exopeptidase showed the highest activity at a temperature of about 50 °C and a pH value of about 7. At temperatures between 40 and 80 °C and pH values between 5 and 9, varying temperatures and pH values lead to clear changes of enzymatic activity. Here, good temperature and pH control is required to achieve a stable process performance.

The DSig functions could be fitted to the experimentally determined temperature and pH dependency of the endo- and exopeptidase investigated, by setting the parameters of the DSig equations as shown in

Table 3. Parameters of the DSig equations for mapping the temperature (T) and pH dependency of the endo-(E1) and exopeptidase (E2) used in proteolysis (P2).

Parameter	Description	Value
YLS,T,act,E1,P2	DSig value at low temperature.	0.0
Ymid,T,act,E1,P2	DSig value between rmax,Tlow,E1,P2 and rmaxT,high,E1,P2.	1.5
YRS,T,act,E1,P2	DSig value at high temperature.	0.0
KLS,T,act,E1,P2	Slope on the low side of the DSig function.	0.1
KRS,T,act,E1,P2	Slope on the high side of the DSig function.	0.2
rmax,Tlow,act,E1,P2	Location parameter of the low side of the DSig function.	43.0 °C
rmax,Thigh,act,E1,P2	Location parameter of the high side of the DSig function.	65.0 °C
YLS,pH,act,E1,P2	DSig value at low pH.	0.0
Ymid,pH,act,E1,P2	DSig value between rmax,pHlow,act,E1,P2 and rmax,pHhigh,act,E1,P2.	1.1
YRS,pH,act,E1,P2	DSig value at high pH.	0.0
KLS,pH,act,E1,P2	Slope on the low side of the DSig function.	2.5
KRS,pH,act,E1,P2	Slope on the high side of the DSig function.	2.0
rmax,pHlow,act,E1,P2	Location parameter of the low side of the DSig function.	5.8
rmax,pHhigh,act,E1,P2	Location parameter of the high side of the DSig function.	8.6
YLS,T,act,E2,P2	DSig value at low temperature.	0.0
Ymid,T,act,E2,P2	DSig value between rmax,Tlow,E2,P2 and rmaxT,high,E2,P2.	1.1
YRS,T,act,E2,P2	DSig value at high temperature.	0.0
KLS,T,act,E2,P2	Slope on the low side of the DSig function.	0.1
KRS,T,act,E2,P2	Slope on the high side of the DSig function.	0.2
rmax,Tlow,act,E2,P2	Location parameter of the low side of the DSig function.	28.0 °C
rmax,Thigh,act,E2,P2	Location parameter of the high side of the DSig function.	64.0 °C
YLS,pH,act,E2,P2	DSig value at low pH.	0.0
Ymid,pH,act,E2,P2	DSig value between rmax,pHlow,act,E2,P2 and rmax,pHhigh,act,E2,P2.	1.8
YRS,pH,act,E2,P2	DSig value at high pH.	0.0
KLS,pH,act,E2,P2	Slope on the low side of the DSig function.	1.2
KRS,pH,act,E2,P2	Slope on the high side of the DSig function.	1.0
rmax,pHlow,act,E2,P2	Location parameter of the low side of the DSig function.	6.4
rmax,pHhigh,act,E2,P2	Location parameter of the high side of the DSig function.	8.4

.

Due to the detailed representation of the temperature and pH dependency of the enzymes, the model can be used to find the temperature and pH values at which the two enzymatic processes run best in combination. It is of major importance that the subsequent DT can reproduce the properties of the enzymes, to predict the optimal process settings.

3.4. Designing Experiments Using the Model of the Enzymatic Hydrolysis Processes

To ensure the usability of the subsequent DT, the accuracy of the developed model for enzymatic hydrolysis processes was investigated. Therefore, batch experiments for starch hydrolysis and proteolysis in the STR were planned and simulated with the developed model. The preliminary parametrization was derived from the characterization experiments described before. Subsequently, the experiments were carried out and the model parameterization was adjusted if necessary.

3.4.1. Starch Hydrolysis

Figure 6 compares the results of the starch hydrolysis experiments in the STR with the results simulated in advance with the model.

About 120 g L⁻¹ of the product (glucose) was formed from 125 g L⁻¹ hydrolysable substrate (starch) in a processing time of 140 min. The results simulated in advance with the model, deviate from the experimentally determined results particularly for the product formed. In the simulation, the product concentration (P_sim) reaches a maximum value of around 60 g L⁻¹ and is, thus, only half as large as the experimentally determined product concentration of around 120 g L⁻¹. In addition to the substrate (S_sim), intermediate product (IP_sim) and product (P_sim) concentrations, the model also calculates the concentrations of non-hydrolysable components of the substrate (SNH_sim) and the two enzymes (E1_sim, E2_sim) used. E1_sim and E2_sim show slight decreases due to enzyme denaturation. Since the experiment shown was carried out close to the temperature and pH optimum of the enzymes, this denaturation is very low. To investigate the critical process areas with the model or the DT, it is necessary to be able to simulate the concentrations of the enzymes. In the batch experiment shown, SNH_sim is constant. However, if the model or the DT is used to simulate fed-batch experiments, it is important to know the concentration of the non-hydrolysable components of the substrate to avoid accumulation in the reactor.

For higher accuracy of the model, it was reparametrized using the new measurement data. For reparameterization, the maximum reaction rates (r_max), as well as the half-saturation constants (K_M), were adjusted according to

Table 4. Adjustment of starch hydrolysis model parameterization, with α-amylase (E1) and glucoamylase (E2).

Parameter	Description	Initial	Adjusted
rmax,SIP,E1,P1	Maximum degradation rate of substrate to intermediate product.	1.10 s−1	2.00 s−1
rmax,SP,E1,P1	Maximum degradation rate of substrate to product.	1.30 s−1	0.10 s−1
rmax,IPP,E2,P1	Maximum degradation rate of intermediate product to product.	1.80 s−1	1.99 s−1
KM,SIP,E1,P1	Half-saturation constant for degradation of substrate to intermediate product.	1.40 g L−1	15.04 g L−1
KM,SP,E1,P1	Half-saturation constant for degradation of substrate to product.	9.00 g L−1	20.00 g L−1
KM,IPP,E2,P1	Half-saturation constant for degradation of intermediate product to product.	9.00 g L−1	8.32 g L−1

.

After reparameterization, r_{max,SIP,E1,P1} was increased from 1.10 s⁻¹ to 2.00 s⁻¹, r_max,SP,E1,P1 was lowered from 1.30 s⁻¹ to 0.10 s⁻¹, r_{max,IPP,E2,P1} was increased from 1.80 s⁻¹ to 1.99 s⁻¹, K_M,SIP,E1,P was increased from 1.40 g L⁻¹ to 15.04 g L⁻¹, K_M,SP,E1,P1 was increased from 9.00 g L⁻¹ to 20.00 g L⁻¹ and K_M,IPP,E2,P1 was lowered from 9.00 g L⁻¹ to 8.32 g L⁻¹. Figure 7 shows the comparison between the simulated and the experimental results after reparameterization of the model.

After adjustment of the model parameters, the simulated trajectories fit the experimental data much better than before. Reparameterization increased the coefficient of determination (R²) for the substrate from 0.95 to 0.96, and for the product from 0.20 to 0.88.

3.4.2. Proteolysis

Figure 8 compares the results of the proteolysis experiment in the STR with the results simulated by the model.

About 5–6 g L⁻¹ of product (fAA) were formed from, approximately, an estimated 14 g L⁻¹ hydrolysable substrate (protein) in a processing time of 300 min. The results simulated in advance with the model, deviate from the experimentally determined results. In the simulation, the product concentration (P_sim) reaches a maximum value of around 15 g L⁻¹ and is thus approximately 3 g L⁻¹ higher than the experimentally determined product concentration of around 12 g L⁻¹. In addition to the substrate (S_sim), intermediate product (IP_sim) and product (P_sim) concentrations, the model also calculates the concentrations of non-hydrolysable components of the substrate (SNH_sim) and the two enzymes (E1_sim, E2_sim) used. E1_sim and E2_sim show slight decreases due to enzyme denaturation. Since the experiment shown was carried out close to the temperature and pH optimum of the enzymes, this denaturation is very low. Moreover, in proteolysis, it is necessary to be able to simulate the concentrations of the enzymes and the non-hydrolysable components of the substrate. The reasons for this correspond to those described in the previous section.

For higher accuracy of the model, it was reparametrized using the new measurement data. Reparameterization led to new values of the maximum reaction rates (r_max), as well as the half-saturation constants (K_M) that were adjusted according to

Table 5. Adjustment of proteolysis model (P2) parameterization with endo- (E1) and exopeptidase (E2).

Parameter	Description	Initial	Adjusted
rmax,SIP,E1,P2	Maximum degradation rate of substrate to intermediate product.	0.70 s−1	0.10 s−1
rmax,SP,E1,P2	Maximum degradation rate of substrate to product.	0.16 s−1	0.16 s−1
rmax,IPP,E2,P2	Maximum degradation rate of intermediate product to product.	0.15 s−1	0.15 s−1
KM,SIP,E1,P2	Half-saturation constant for degradation of substrate to intermediate product.	6.94 g L−1	5.94 g L−1
KM,SP,E1,P2	Half-saturation constant for degradation of substrate to product.	7.80 g L−1	9.80 g L−1
KM,IPP,E2,P2	Half-saturation constant for degradation of intermediate product to product.	2.10 g L−1	0.10 g L−1

.

During reparameterization, r_{max,SIP,E1,P2} was lowered from 0.70 s⁻¹ to 0.10 s⁻¹, K_M,SIP,E1,P2 was lowered from 6.94 g L⁻¹ to 5.94 g L⁻¹, K_M,SP,E1,P2 was increased from 7.80 g L⁻¹ to 9.80 g L⁻¹, and K_M,IPP,E2,P2 was lowered from 2.10 g L⁻¹ to 0.10 g L⁻¹. No changes were made to any other model parameters. Figure 9 shows the comparison between the simulated and the experimental results after reparameterization of the model.

After adjustment of the model parameters, the simulated trajectories for the product concentration show a much better fit to the experimental data than before. No experimental data were available for the substrate concentration. Thus, the initial value of the simulation was set to 40 g L⁻¹, as in the experiment. The course of the simulated substrate concentration is plausible, as it is strongly related to the increase of the product. Reparameterization of the model increased the coefficient of determination (R²) for the product concentration from −6.76 to 0.76.

3.5. Application of the DT for Enzymatic Hydrolysis Processes

Once the mathematical model of the enzymatic hydrolysis processes was able to satisfactorily represent the enzymatic processes of starch hydrolysis and proteolysis, the model was implemented into the existing DT described before. Using the DT for the enzymatic hydrolysis processes, it is possible to accelerate the simulation of the enzymatic hydrolysis processes in a 20 L STR up to 100-fold.

For the process of the enzymatic starch hydrolysis in the DT, PID controls for temperature and product (glucose) concentration were implemented using WinErs. The temperature in the reactor is controlled via the temperature of the heating fluid in the inflow to the heating jacket, while the product concentration is controlled via the inflows of the two enzymes (α-amylase and glucoamylase) into the reactor. For simplification, it was assumed in the study that continuous glucose measurement was available in the DT.

In the first step, the parameters of the PID controllers were determined using the methods developed by Ziegler and Nichols [40]. Then, the controller parameters were adjusted manually and the control result was examined with the help of simulations in WinErs. During manual PID tuning, special care was taken to ensure that the state variable does not leave the specified tolerance range and that the settling time is kept as short as possible. A tolerance range of ±2 g L⁻¹ was specified for product concentration control and ±0.5 °C for temperature control. For initial investigations, the manual PID tuning procedure was applied. In the future, it is planned to determine the controller parameters using automated optimization algorithms, such as, e.g., Nelder–Mead [26]. Figure 10 shows the simulation results of the developed control strategies.

Initially, the process was in a steady-state condition. At a temperature of 60 °C, the product concentration could be maintained at the setpoint of 100 g/L. After a processing time of approximately 1.5 h, the temperature was lowered to 55 °C. This led to an increase in product concentration to approximately 105 g/L. After a processing time of about 6 h, the inflow of enzymes could be adjusted by the controller so that the product concentration approached the setpoint of 100 g/L again. After a processing time of 10 h, the temperature in the reactor was increased from 55 °C to 65 °C, which resulted in the product concentration dropping to near 90 g/L. The PID controller thereafter increased the inflow of enzymes from 0.04 mL/min to over 0.6 mL/min to get the product concentration back to 100 g/L.

The process of starch hydrolysis is strongly dependent on the temperature of the broth in the reactor. Changes of only 5 °C led to strong changes in the enzyme characteristics, to which the control must respond. Figure 10b shows how the developed DT for enzymatic hydrolysis processes can simulate the temperature behavior of the represented STR. When the setpoint increases from 55 °C to 65 °C, the temperature of the broth in the reactor does not rise abruptly to 65 °C. First, the heating jacket is heated by the heating fluid, which then heats the reactor wall, which in turn heats the broth in the reactor. In all these heat transfer steps, as in reality, there are heat losses to the environment (T = 25 °C).

For the development of efficient controls for enzymatic hydrolysis processes, it is therefore necessary that not only the reaction kinetics (micro-kinetics), but also the dynamic behavior (macro-kinetics) of the STR and its periphery can be simulated in detail. For this reason, the use of the developed DT is beneficial for process development.

4. Discussion and Outlook

Within the scope of this work, a dynamic and mechanistic mathematical model for enzymatic hydrolysis processes, which is easily adaptable (generic), has been developed. It could be shown that the model can be adapted to experimental data by adjusting the parameters of the Michaelis–Menten kinetics and the DSig equations used. As a result, the model can simulate the starch hydrolysis and proteolysis process with an agreement of over 90% and over 85%, respectively. Since the model is adaptable, it is possible to respond to changes in enzyme and substrate properties. Due to the model structure used, it is also possible to replace the Michaelis–Menten kinetics with other enzyme kinetics.

The model could successfully be integrated into a DT of a 20 L STR. Through this integration, not only the reaction kinetics of the enzymatic hydrolysis processes (micro-kinetics) can be simulated, but also the macro-kinetics of the STR and its periphery (e.g., pumps, heating jacket, etc.). It could be shown that the temperature behavior of the STR can affect the reaction kinetics (starch hydrolysis). Thus, the developed DT is beneficial for the design of control and automation strategies since it can realistically reproduce the interaction between micro- and macro-kinetics.

Process control strategies developed with the new DT for enzymatic hydrolysis processes in the STR can be transferred to the real process quickly and with a high probability of success. This has already been demonstrated for the cultivation process of S. cerevisiae [4].

In the future, it is planned to further develop the coupling of the real process and the DT. Currently, the DT is being adapted to the real process after the experiments have been carried out. In the future, it is planned that the DT will continuously adapt to the running process.

By combining the newly developed DT and model-based tools for process optimization [25,41], the enzymatic hydrolysis process to produce organic nutrient media can be systematically improved. Since both the model of the enzymatic hydrolysis processes and the DT were built generically, the DT for the enzymatic hydrolysis processes can be quickly and easily adapted to other enzymatic processes as well as other reactor configurations.