Kinetic Parameter Determination for Depolymerization of Biomass by Inverse Modeling and Metaheuristics

Abstract

A computational methodology based on inverse modeling and metaheuristics is presented for determining the best parameters of kinetic models aimed to predict the behavior of biomass depolymerization processes during size scaling up. The Univariate Marginal Distribution algorithm, particle swarm optimization, and Interior-Point algorithm were applied to obtain the values of the kinetic parameters (K_M and V_max) of four mathematical models based on the Michaelis–Menten equation: (i) Traditional Michaelis–Menten, (ii) non-competitive inhibition, (iii) competitive inhibition, and (iv) substrate inhibition. The kinetic data were obtained from our own experimentation in micro-scale. The parameters obtained from an optimized micro-scale experiment were compared with a bench scale experiment (0.5 L). Regarding the metaheuristic optimizers, it is concluded that the Interior-Point algorithm is effective in solving inverse modeling problems and has the best prediction power. According to the results, the Traditional model adequately describes the micro-scale experiments. It was found that the Traditional model with optimized parameters was able to predict the behavior of the depolymerization process during size scaling up. The methodology followed in this study can be adopted as a starting point for the solution of future inverse modeling problems.

1. Introduction

Vegetal biomass (VB) is an alternative for biofuel production. VB is mainly composed of cellulose, hemicellulose, and lignin and can be transformed into different products through four fundamental stages: (i) conditioning, (ii) depolymerization or enzymatic hydrolysis, (iii) fermentation, and (iv) the downstream process [1,2]. In this context, the depolymerization stage is the most expensive, and predicting the reaction behavior is essential to achieve a higher reaction yield. One of the most widely used models to describe the behavior of this phenomenon is the Michaelis–Menten (M–M) equation. The M–M equation consists of two parameters: the Michaelis constant (K_M) and the maximum reaction velocity (V_max), which describe the speed of the enzymatic reactions related to the reaction rate (r_s) and the concentration of a substrate (S), as described in Equation (1) [3].

$$-r_s=-\frac{dS}{dt}=\frac{V_{max}\times S}{K_M+S}$$

(1)

The determination of these parameters (V_max and K_M) allows one to predict the behavior of the depolymerization reaction. This equation has been used to follow various enzymatic reactions, e.g., hydrolysis of proteins, biomass, and even the elimination or interaction of drugs in the human body [2,3,4]. A methodology that can be followed to obtain these parameters is inverse modeling (IM) supported by computational intelligence (CI) methods. On the one hand, IM consists of estimating the parameters in a mathematical model from measured observations [5]. This type of estimation begins with comprehensive data at the level of intact systems. It leads to inferences of parameter values at lower levels of organization—namely, features of individual processes. Thus, it is called “top-down” or “inverse” [6].

On the other hand, CI can be defined as a set of nature-inspired computational models that comprehend the subfields of fuzzy logic, artificial neural networks, and evolutionary algorithms. These models are designed to address real-world problems of prediction, control, classification, optimization, and other tasks [7]. Optimization methods of CI, also known and hereafter referred to as metaheuristics, have been successfully applied to identify a vast number of nonlinear dynamic systems [8] and parameter estimations and facilitate the structure identification of biochemical systems [6]. To address the IM, not only metaheuristics (such as genetic algorithms, simulated annealing particle swarm optimization (PSO), or Estimation of Distribution algorithms (EDAs) but also other algorithms for specific biochemical systems (such as alternating regression, eigenvector optimization, and dynamic flux estimation) have been applied [6]. Among the available methods, EDAs and PSO have been recently studied for biochemical parameter estimation and have shown to be optimal algorithms for some hyper-parameter tuning applications [9,10].

In the biofuel area, depolymerization experiments are developed at the flask level, evaluating multiple parameters resulting in large assay experiments with intensive requirements. To reduce these intensive requirements and enable the evaluation of numerous parameters, micro-scale experiments have been proposed as an alternative tool to optimize some aspects of bioprocesses (e.g., biomass depolymerization) [11]. In addition, few reports have been published on an assay of the depolymerization process at a micro-scale [11,12,13]. This may be because the methodology has not been validated as a tool that can predict the behavior of depolymerization at larger scales.

Given the high number of experiments and parameters to evaluate, the implementation of a numerical analysis is necessary to study and improve the depolymerization process. The main objective of this research is to obtain the best parameters for kinetic models based on the Michaelis–Menten equation to identify the model with the highest predictive power at micro-scale data so that it can be used to predict the depolymerization process during size scaling up. To achieve this objective, a computational methodology based on IM and metaheuristics was applied to 14 micro-scale (1 mL) depolymerization experiments previously reported in [12].

The next section describes how the experimental observations were obtained, the formulation of IM as an optimization problem, the metaheuristics employed for kinetic parameter determination, and the procedures followed to apply the computational methodology.

2. Materials and Methods

2.1. Experimental Data Collection

The experimental procedures for biomass pretreatment and enzymatic hydrolysis were previously reported by Molina et al. [12]. Nevertheless, the trajectories of the experimental kinetics were not reported, and the kinetic parameters were not determined in that study. Briefly, the substrate used was wheat straw (WS) delignified with the alkaline peroxide method. The enzyme used was Accellerase 1500^® gently donated by Genencor^® (Danisco, Rochester, NY, USA). The experimental observations were obtained from micro-scale experiments (1 mL of work volume). These experiments were subject to an experimental optimization design (Central Composite Design (CCD)) while varying the operating parameters (potential of hydrogen (pH), temperature (°C), and enzyme–substrate ratio (μL/mg)) as shown in

Table 1. Experimental conditions for enzyme depolymerization at micro-scale. pH stands for potential of hydrogen.

Experiment	pH	°C	µL/mg	Experiment	pH	°C	µL/mg
1	4.3	48	1.4656	8	4.0	55	0.9000
2	4.3	62	0.3343	9	6.0	55	0.9000
3	4.3	62	1.4656	10	5.0	45	0.9000
4	5.7	48	0.3343	11	5.0	65	0.9000
5	5.7	48	1.4656	12	5.0	55	0.1000
6	5.7	62	0.3343	13	5.0	55	1.8000
7	5.7	62	1.4656	14	5.0	55	0.9000

.

The CCD allowed us to obtain a statistical model that showed the optimal operating parameters of the depolymerization process. A micro-scale experiment was developed using the optimized operative parameters, and the trajectory was followed. Finally, a bench scale experiment (0.5 L) was performed, thereby reproducing the operative parameters obtained in the optimized micro-scale experiment. The main difference between the micro-scale and the bench scale experiments was the stirring. The micro-scale experiments used high orbital shaking (700 rpm), while the bench scale experiment used propeller shaking (150 rpm). At both scales, the orbital and propeller shaking were high, allowing us to maintain spatial homogeneity, which was previously suggested for the study of reactors as ideal systems [14]. However, during scaling, homogeneity may be reduced and cause interference in the mass transfer. The details are described in [12].

2.2. Inverse Modeling

The ability of four mathematical models based on the Michaelis–Menten equation (Equations (2)–(5)) to fit the experimental observations of biomass depolymerization was evaluated. For these models, the parameters V_max, K_M, I, and K stand for the maximal velocity, substrate concentration (S) at half-maximal velocity, inhibitor concentration, and characteristic constant, respectively. The K parameter is an offset that enables the kinetic models to adapt to the initial conditions.

Michaelis–Menten:

$$-r_s=-\frac{dS}{dt}=\frac{V_{max}\times S}{K_M+S}+K$$

(2)

Non-Competitive Inhibition:

$$-r_s=-\frac{dS}{dt}=\frac{V_{max}\times S}{(K_M+S)\times (1+\frac{I}{K_I})}+K$$

(3)

Competitive Inhibition:

$$-r_s=-\frac{dS}{dt}=\frac{V_{max}\times S}{K_M\times (1+\frac{I}{K_I})+S}+K$$

(4)

Substrate Inhibition:

$$-r_s=-\frac{dS}{dt}=\frac{V_{max}\times S}{\left(K_M+S+\frac{S^2}{K_I}\right)}+K$$

(5)

The depolymerization data were calculated indirectly according to the concentration of the reducing sugar (RS) for each time step using Equation (6):

$${\mathrm{C}}_{ws}^0=X_{ws}{C}_{ws}; {\mathrm{C}}_{\mathrm{ws}}=C_{ws}^0-\frac{M{W}_{cellulose}}{M{W}_{RS}}{C}_{RS}$$

(6)

where $C_{ws}^0$ is the initial concentration of $WS$, $X_{ws}$ is the fraction of cellulose + hemicellulose in the dried substrate (0.825 mg/mg), $C_{ws}$ is the dried $WS$ concentration in the enzymatic hydrolysis media (10 mg/mL), $MW$ is the molecular weight of the cellulose, $C_{ws}$ is the concentration of the substrate, and $C_{RS}$ is the concentration of reducing sugars reported as glucose equivalents [12].

In this work, the aim of the inverse modeling is to find the Michaelis–Menten model that best fits the substrate concentration $S={\left(S_1,\dots S_n\right)}^{\mathrm{T}}$ to the experimental data $C_{ws}={\left(C_{ws}^1,\dots ,C_{ws}^n\right)}^{\mathrm{T}}$ obtained from $n$ instants during the depolymerization process. To measure the quality of the fit achieved by each model, the weighted least square cost function ${\chi }^2={\displaystyle \sum }_{i=1}^n{w}_i·{\left(S_i-C_{ws}^i\right)}^2$ is implemented. The weight factor is defined by $w={\left(w_1,\dots ,w_n\right)}^{\mathrm{T}}={\left(1/{\sigma }_1^2,\dots ,1/{\sigma }_n^2\right)}^{\mathrm{T}}$, where ${\sigma }_n$ is the standard deviation [15]. The four models require a set of parameters ($\mathsf{\theta }\in {\mathrm{R}}^4$) to be specified, such that $S$($\mathsf{\theta }$) can be obtained by solving the associated differential equation. Thus, the measure ${\chi }^2$ can be defined to quantify how well an enzymatic model reproduces $C_{ws}$ for a given $\mathsf{\theta }$. Therefore, the aim of the inverse modeling can be achieved through a parameter identification process formulated as the following optimization problem:

$${\mathsf{\theta }}^{*}=\underset{\mathsf{\theta }}{\min }{\chi }^2\left(\mathsf{\theta }\right);\mathrm{subject}\mathrm{to}0<{\theta }_{\mathrm{i}},\mathrm{i}=1,\dots ,4$$

(7)

where ${\mathsf{\theta }}^{*}={\left(V_{max}^{*},K_m^{*},I^{*},K^{*}\right)}^{\mathrm{T}}$ is the optimal kinetic parameter vector obtained from the four-dimensional space of the positive real values excluding the origin. Since it is unknown if two different sets of parameters ${\mathsf{\theta }}_1$ and ${\mathsf{\theta }}_2$ could lead to the same measure ${\chi }^2$, we cannot assume a functional relationship between $\mathsf{\theta }$ and ${\chi }^2$. Furthermore, the relevant characteristics of this measure, such as the convexity or roughness of its corresponding response surface, are also unknown. Consequently, sophisticated optimizers need to be applied to solve Equation (7).

2.3. Metaheuristics

Metaheuristics are approximate optimization techniques that provide acceptable solutions for solving complex problems [16]. These techniques are frequently used to deal with optimization problems where the solution cannot be found with exact methods or where the computational time required to draw the solution is so high that the problem becomes intractable.

For the inverse modeling problem stated in Equation (7), preliminary results were reported with the Nelder–Mead algorithm in [17]. This study showed that Nelder–Mead is effective at minimizing ${\chi }^2$; however, the best kinetic parameters found (${\mathsf{\theta }}^{*}$) with this method lacked physical meaning. The reason for this result is that Nelder–Mead is an unconstrained method; therefore, the search for the optimal set of parameters is prone to reaching infeasible regions.

The problem presented by the Nelder–Mead algorithm can be solved via constrained optimization methods that retain the candidate solutions inside feasible regions. Three generic constrained optimization methods were employed: (i) the Univariate Marginal Distribution algorithm ($UMD{A}_c^G$, hereafter referred to as UMDA) [18]; (ii) particle swarm optimization [19]; (iii) the Interior-Point algorithm [20]. These three methods differ in their exploration of search spaces and are selected to understand the different aspects of the quality measure ${\chi }^2$. The next lines describe the implementation of these methods.

UMDA and PSO are metaheuristics from different families. UMDA is an example of the Estimation of Distribution algorithms. EDAs explore the solution space by iteratively building and sampling explicit probabilistic models of candidate solutions that guide their search. The generic procedure, presented in Algorithm 1, is as follows [21]: The method starts by measuring the quality of a set of random candidate solutions sampled from a uniform distribution. Then, a fraction of the best solutions is selected to construct an estimate of their probability distribution (the mean and standard deviation for a normal distribution in the case of the UMDA). By sampling this distribution, new solutions are generated to update the model. Thus, the estimated probability distribution is improved. A bias to the best solutions is the fundamental idea behind EDAs; the aim is that after some iterations, the random candidate solutions will converge to an optimal solution (assuming that this solution exists and is unique, i.e., assuming that the solution space is convex). The proof of the convergence of the algorithm implemented in our work (UMDA), as well as that of other EDAs, is provided in [22]. This process is repeated until a termination criterion is met (usually, when a maximum number of iterations is reached). Depending on the specific probabilistic model built, an EDA takes a name. The UMDA then builds a Gaussian model $P$ from a set of solutions $\mathsf{\Theta }={\left\{{\mathsf{\theta }}_j\right\}}_{j=1}^M$, with the parameters given by (Equation (8)):

$${\mu }_{\kappa }=\frac{1}{M}{\displaystyle \sum }_{j=1}^M{\theta }_{j,\kappa }; {\sigma }_k={\left(\frac{1}{M-1}{\displaystyle \sum }_{j=1}^M{\left({\theta }_{j,\kappa }-{\mu }_{\kappa }\right)}^2\right)}^{1/2}$$

(8)

where $M$ is the number of solutions considered to compute these parameters, and ${\theta }_{j,\kappa }$ represents the $\kappa$-th component of a solution ${\mathsf{\theta }}_j$, which is a $D$-dimensional vector in ${\mathrm{R}}^D$. The UMDA was implemented in this work to observe if the measure ${\chi }^2$ can be regarded as convex. A drawback of the UMDA is that the algorithm is sensitive to local minima. Thus, if different runs of the UMDA lead to different answers, the response surface of the quality measure may contain more than one minimum, i.e., the quality measure may not be convex.

Algorithm 1 General Estimation of Distribution Algorithm
1:	Generate an initial population of $N$ candidate solutions: ${\mathsf{\Theta }}^{\left(0\right)}={\left\{{\mathsf{\theta }}_j\right\}}_{j=1}^N$ at iteration $t=0$, uniformly
2:	WHILE (stopping criteria are not met) // $N$ = 100 in this work
3:	Compute the quality measure of each candidate solution $g\left({\mathsf{\theta }}^{\left(t\right)}\right)$ // ${\chi }^2$ in this work
4:	Select a subset of the best solutions, ${\mathsf{\theta }}_M^{\left(t\right)}$ from ${\mathsf{\theta }}^{\left(t\right)}$ // $M=0.5 N$ in this work
5:	Build a probabilistic model $P^{\left(t\right)}$ based on ${\mathsf{\theta }}_M^{\left(t\right)}$
6:	Sample $P^{\left(t\right)}$ to generate new solutions ${\mathsf{\theta }\prime }^{\left(t\right)}$
7:	Substitute/Incorporate ${\mathsf{\theta }\prime }^{\left(t\right)}$ into ${\mathsf{\theta }}^{\left(t\right)}$
8:	Update $t=t+1$
9:	END WHILE
10:	Output the best solution found, ${\mathsf{\theta }}^{*}$, as a result

PSO is a member of the bio-inspired metaheuristics family that emulates the collective behavior of swarms (such as a flock of birds or a school of fish) to explore continuous search spaces looking for the optimal value. In this simulation, each particle of the swarm corresponds to a candidate solution ($\mathsf{\theta }$). An iterative process like the UMDA is then followed, except for the model-updating mechanism. For each iteration ($t$) of the PSO, the particle with the best quality measure (${\mathsf{\theta }}^{*\left(\mathrm{t}\right)}$) is identified. Then, the rest of the particles (${\left\{{\mathsf{\theta }}_j\right\}}_{j=1}^{N-1}$) update their positions according to the directional vector (${\mathsf{\theta }}_j-{\mathsf{\theta }}^{*\left(\mathrm{t}\right)}$) plus a random walk, described as Equation (9) [23]:

$${\mathrm{v}}_j^{t+1}={\mathsf{\omega }\mathrm{v}}_j^t+\mathrm{U}\left(0,{\varphi }_1\right)\otimes \left({\mathsf{\theta }}_j^b-{\mathsf{\theta }}_j\right)+\mathrm{U}\left(0,{\varphi }_2\right)\otimes \left({\mathsf{\theta }}^{*\left(\mathrm{t}\right)}-{\mathsf{\theta }}_j\right); {\mathsf{\theta }}_j^{t+1}={\mathsf{\theta }}_j^t+{\mathsf{\theta }}_j^{t+1}$$

(9)

where $\mathrm{U}\left(0,\varphi \right)$ represents a vector of random numbers, uniformly distributed in $\left[0,\varphi \right]$ and generated at each iteration for each particle, and the operator $\otimes$ represents a component-wise multiplication. The parameters ${\varphi }_1$ and ${\varphi }_2$ are known as acceleration coefficients. They determine the magnitude of the random forces in the direction of the best position in the trajectory of a given particle (${\mathsf{\theta }}_j^b$) and the neighborhood best ${\mathsf{\theta }}^{*\left(\mathrm{t}\right)}$. The values of ${\varphi }_1={\varphi }_2=2.0$ are generally adopted. $\omega$ is termed the inertia weight, which can be interpreted as the fluidity of the medium in which a particle moves; this value is commonly initialized to a high value (e.g., 0.9) and gradually reduced to a low value (e.g., 0.4). The PSO is considered in this study since it has been proven to solve highly complex optimization problems. Importantly, it has the advantage of avoiding local minima. Thus, if the quality measure of the inverse modeling problems is rough, the PSO can offer better answers than the UMDA. The main drawback of PSO is its stochastic nature, which means that no convergence warranties can be offered.

Interior-Point algorithms (IPAs) are not metaheuristics in a strict sense but exact methods that improve the performance of the optimization process in terms of both computational time and the ability to handle degenerate problems (by employing one of several heuristics) [24]. IPAs emerged as an extension of the simplex method to solve complex linear programming problems (involving thousands of constraints and variables). In contrast to the simplex algorithm, where the movement of candidate solutions ($\mathsf{\theta }$) occurs along the boundary of the feasible region, the points generated by IPAs lie in the interior of the inequality constraints, hence the name “Interior-Point”. IPAs provide a natural way to incorporate inequality constraints into the objective function and transform a constrained optimization problem into an unconstrained one. The idea is to add barrier functions to the original objective function like in the Lagrange multiplier method [25]. A common formulation of the barrier method is the following:

$$\underset{\theta ,s}{\min }f\left(\theta \right)-\mu {\displaystyle \sum }_{i=1}^m\ln s_i; \mathrm{subject}\mathrm{to}h\left(\theta \right)=0;g\left(\theta \right)+s=0$$

(10)

where $\mu >0$ is the barrier parameter, and the slack variable $s$ is assumed to be positive. From among the several existing Interior-Point algorithms, we implemented the version in [20] (which starts from Equation (10)) because it includes heuristics in the trust region to 1) avoid approaching the boundary of the feasible regions too soon and 2) allows the algorithm to treat convex and non-convex problems uniformly by facilitating the use of second derivative information when the problem is non-convex. In this work, the Interior-Point algorithm was implemented to evaluate if the measure ${\chi }^2$ could be treated as a function of the parameter set $\left(\mathsf{\theta }\right)$. Thus, for this method, it is assumed that $f\left(\mathsf{\theta }\right)={\chi }^2$. The main drawback of IPAs is that they are sensitive to an initial estimate of the solution required to perform the optimization process, suggesting that a different initial point can lead to a different solution [26].

The PSO and the Interior-Point algorithms were adapted to solve the inverse modeling problem by using the particle swarm and fmincon MATLAB^® procedures (with default hyper-parameter values), respectively. As the UMDA is not a MATLAB procedure available in the optimization toolbox, it was implemented by following Algorithm 1. The differential equations of the models described in Equations (2)–(5) were solved using the ODE45 procedure as an implicit step in the quality measure computations. The initial values of $V_{max}$ and $K_M$ were calculated by the Eaddie–Hofstee diagram [27].

Although the quality measure ${\chi }^2$ is effective to measure the fit level accuracy of kinetic models to experimental data, this measure is not informative. To ease the interpretation of the results provided by the optimization methods, the proportions of the variation between the fitted models and experimental data were obtained by calculating the determination coefficients $\left(R^2\right)$. The closer to the unit this coefficient is, the more qualitative the model. A zero value of $R^2$ means that the kinetic model behaves as a constant function, and a negative value indicates that more variables need to be integrated into the model to explain the experimental data. As suggested in [28], the performance of predictive models should be measured by following a multidimensional approach. In this work, the $R^2$ was considered to evaluate the correspondence between predictions and experimental data, the mean square error (MSE) was used to evaluate the prediction accuracy, and Akaike’s information criterion (AIC) was used to quantify the model complexity. These three metrics are described as Equations (11)–(13):

$$R^2=1-\frac{S{S}_{res}}{S{S}_{tot}}=\frac{{{\displaystyle \sum }}_i{\left(C_{ws}^i-S_i\right)}^2}{{{\displaystyle \sum }}_i{\left(C_{ws}^i-{\overline{C}}_{ws}\right)}^2}$$

(11)

where $i$ refers to the i-th datum for each experiment. Further, the mean square errors were calculated by:

$$MSE=\frac{S{S}_{res}}{n}$$

(12)

where n is the number of data in each experiment.

The AIC combines the quality of fit and the parameter number $\left(prm\right)$ of models, allowing to compare models of varying complexity [29]. The AIC is defined by the Equation (13):

$$\mathrm{AIC}=n·ln\left(\frac{{SS}_{res}}{n}\right)+2\left(prm+1\right)$$

(13)

2.4. Experimental Methodology

The procedure followed for the evaluation of kinetic models is summarized in Algorithm 2. Each experimental dataset, with its conditions reported in Table 1, was partitioned to conduct a 4-fold cross validation scheme. For the sake of a fair comparison between models, each model employed the same data for training along with the same test data. The final performance was computed as an average of the four cross-validations, and these results are reported in

Table 2. Best-fitting models and parameters found by optimizer: Interior-Point Algorithm (IPA), Particle Swarm Optimization (PSO) and Univariate Marginal Distribution Algorithm (UMDA). V_max (mg/mL × s), K_M (mg/mL), I (mg/mL). MSE: mean square error.

Experiment	Optimizer	Best-Fitting Kinetic Model	R2 Train	AIC Mean ± Std	${\mathit{\chi }}^2$	VMax*	Km*	K*	${\mathit{I}}^{*}$	MSE Train	MSE Test
1	IPA	Traditional	0.98 ± 0.03	−15.29 ± 11.52	2.77 ± 5.0	1.43	0.58	1.25		0.02 ± 0.03	0.01 ± 0.02
	PSO	Traditional	0.88 ± 0.10	−2.01 ± 5.80	15.76 ± 3.59	0.83	8.25	0.32		0.19 ± 0.19	1.05 ± 1.38
	UMDA	Non-Competitive I.	0.81 ± 0.12	2.88 ± 5.11	23.55 ± 14.6	2.19	3.29	0.80	4.33	0.31 ± 0.23	0.76 ± 0.57
2	IPA	Traditional	0.86 ± 0.21	−11.08 ± 5.45	12.72 ± 19.33	1.05	8.25	0.44		0.06 ± 0.06	0.08 ± 0.07
	PSO	Non-Competitive I.	0.82 ± 0.25	−8.92 ± 6.32	31.17 ± 55.71	3.21	8.25	0.73	4.02	0.06 ± 0.05	0.09 ± 0.06
	UMDA	Non-Competitive I.	0.82 ± 0.14	−6.33 ± 5.81	16.19 ± 19.77	3.47	4.21	1.04	6.66	0.09 ± 0.07	0.1 ± 0.09
3	IPA	Traditional	0.90 ± 0.16	−4.63 ± 8.43	6.32 ± 8.78	0.55	8.25	0.17		0.17 ± 0.22	0.39 ± 0.02
	PSO	Competitive I.	0.90 ± 0.14	−2.51 ± 6.85	7.3 ± 9.95	1.15	7.84	0.15	8.11	0.16 ± 0.15	0.18 ± 0.11
	UMDA	Substrate I.	0.81 ± 0.10	1.63 ± 4.06	22.12 ± 14.81	4.96	7.74	1.23		0.39 ± 0.17	0.81 ± 0.63
4	IPA	Traditional	0.98 ± 0.01	−7.79 ± 3.66	43.98 ± 58.02	0.17	8.25	0.04		0.07 ± 0.05	0.73 ± 1.14
	PSO	Non-Competitive I.	0.80 ± 0.33	6.67 ± 9.17	144.2 ± 188.17	8.04	0.10	2.62	8.25	0.57 ± 0.87	0.85 ± 1.25
	UMDA	Substrate I.	0.86 ± 0.07	3.73 ± 3.13	227.62 ± 121.3	4.80	6.40	1.26		0.49 ± 0.21	1.12 ± 0.91
5	IPA	Non-Competitive I.	0.87 ± 0.06	3.97 ± 2.1	37.01 ± 65.49	0.38	8.25	0.12	0.001	0.37 ± 0.13	0.65 ± 0.62
	PSO	Non-Competitive I.	0.84 ± 0.08	5.06 ± 2.55	37.19 ± 65.92	2.16	2.31	0.6	7.97	0.44 ± 0.20	1.03 ± 1.01
	UMDA	Non-Competitive I.	0.78 ± 0.14	6.48 ± 1.85	38.14 ± 66.29	1.52	5.78	0.34	7.94	0.56 ± 0.16	0.75 ± 0.60
6	IPA	Non-Competitive I.	0.72 ± 0.15	−20.73 ± 6.98	12.86 ± 4.66	2.06	8.25	0.99	5.58	0.01 ± 0.01	0.01 ± 0.00
	PSO	Traditional	0.81 ± 0.09	−20.64 ± 4.58	8.89 ± 3.07	5.37	0.66	4.90		0.01 ± 0.01	0.01 ± 0.01
	UMDA	Traditional	0.64 ± 0.04	−21.18 ± 1.99	15.48 ± 1.56	6.36	6.35	3.50		0.01 ± 0.00	0.02 ± 0.01
7	IPA	Competitive I.	0.76 ± 0.33	−22.2 ± 3.87	47.33 ± 83.02	2.13	1.50	1.44	5.58	0.01 ± 0.00	0.02 ± 0.01
	PSO	Traditional	0.77 ± 0.16	−22.35 ± 2.68	50.41 ± 78.26	3.81	8.25	1.80		0.01 ± 0.00	0.02 ± 0.02
	UMDA	Non-Competitive I.	0.82 ± 0.17	−23.27 ± 3.49	48.61 ± 78.27	6.12	5.15	2.06	5.22	0.01 ± 0.00	0.02 ± 0.01
8	IPA	Traditional	0.89 ± 0.16	−8.22 ± 6.91	396.46 ± 784.32	0.97	8.25	0.38		0.09 ± 0.11	0.06 ± 0.05
	PSO	Traditional	0.89 ± 0.14	−7.5 ± 6.02	398.22 ± 782.66	0.88	3.82	0.5		0.10 ± 0.10	0.10 ± 0.07
	UMDA	Non-Competitive I.	0.75 ± 0.12	1.12 ± 1.54	428.95 ± 764.32	2.57	3.61	1.40	0.97	0.27 ± 0.07	0.38 ± 0.71
9	IPA	Traditional	0.87 ± 0.15	−6.75 ± 4.9	237.31 ± 444.35	0.58	8.25	0.22		0.11 ± 0.10	0.25 ± 0.36
	PSO	Non-Competitive I.	0.86 ± 0.12	−3.38 ± 2.71	243.96 ± 441.48	7.21	0.18	6.86	0.0001	0.14 ± 0.07	0.11 ± 0.12
	UMDA	Non-Competitive I.	0.8 ± 0.09	−0.3 ± 3.58	259.85 ± 419.39	4.58	4.53	1.30	7.01	0.22 ± 0.10	0.23 ± 0.13
10	IPA	Traditional	0.88 ± 0.2	−9.58 ± 6.29	57.97 ± 83.54	0.83	8.25	0.33		0.07 ± 0.09	0.08 ± 0.03
	PSO	Competitive I.	0.9 ± 0.07	−7.79 ± 6.27	88.4 ± 132.63	1.28	2.47	0.85	0.23	0.07 ± 0.05	0.11 ± 0.08
	UMDA	Non-Competitive I.	0.86 ± 0.06	−3.85 ± 5.28	91.28 ± 70.49	5.05	4.26	1.39	7.60	0.13 ± 0.09	0.13 ± 0.13
11	IPA	Traditional	0.86 ± 0.24	−18.82 ± 7.1	57.15 ± 102.38	3.83	7.26	1.83		0.01 ± 0.00	0.07 ± 0.09
	PSO	Traditional	0.86 ± 0.23	−16.58 ± 2.91	63.04 ± 103.31	3.92	8.20	1.78		0.01 ± 0.00	0.01 ± 0.00
	UMDA	Traditional	0.86 ± 0.22	−16.31 ± 2.64	62.35 ± 98.77	5.32	5.06	3.04		0.01 ± 0.00	0.02 ± 0.01
12	IPA	Traditional	0.89 ± 0.07	1.44 ± 7.89	28.75 ± 27.83	0.19	8.25	0.05		0.38 ± 0.24	0.73 ± 0.51
	PSO	Substrate I.	0.87 ± 0.08	2.5 ± 3.77	68.16 ± 41.92	0.53	8.25	0.09		0.43 ± 0.22	1.09 ± 1.02
	UMDA	Substrate I.	0.76 ± 0.18	5.99 ± 2.27	104.04 ± 69.11	5.15	8.21	1.27		0.76 ± 0.28	0.75 ± 0.29
13	IPA	Traditional	0.93 ± 0.07	−3.06 ± 6.19	15.03 ± 23.47	0.55	8.25	0.16		0.20 ± 0.19	0.24 ± 0.15
	PSO	Competitive I.	0.93 ± 0.1	−1.24 ± 6.89	15.11 ± 25.02	0.42	3.99	0.14	1.91	0.19 ± 0.23	0.14 ± 0.08
	UMDA	Substrate I.	0.86 ± 0.13	1.02 ± 4.87	29.65 ± 35.81	3.64	7.61	0.84		0.35 ± 0.23	0.73 ± 0.70
14	IPA	Traditional	0.95 ± 0.02	−7.72 ± 1.46	164.65 ± 309.48	0.78	8.25	0.24		0.10 ± 0.03	0.15 ± 0.10
	PSO	Traditional	0.85 ± 0.19	−4.17 ± 4.47	138.93 ± 249.87	1.07	7.15	0.38		0.17 ± 0.09	0.18 ± 0.11
	UMDA	Non-Competitive I.	0.81 ± 0.23	−0.38 ± 2.95	144.33 ± 247.4	3.94	2.83	1.14	6.80	0.22 ± 0.08	0.26 ± 0.23

. The criterion to select the best model was Akaike’s information since it achieves a trade-off between model complexity and data fitting.

Algorithm 2 Experimental Methodology for Kinetic Model Evaluation
INPUT: A set of experimental data files (D), a set of kinetic models (K), a set of optimizers (O), user-defined search ranges for parameters $V_{max},K_m,I,$ and $K$. Number of folds for data partition and performance metrics OUTPUT: Optimal optimizer, kinetic models, kinetic parameters ${\mathsf{\theta }}^{*}={\left(V_{max}^{*},K_m^{*},I^{*},K^{*}\right)}^{\mathrm{T}}$, and Performance Indexes on train ($P{I}_{train}=R^2,MSE, and AIC$) and test ($P{I}_{test}=MSE$) data.
1:	FOR EACH experimental data file d in D:					// 14 data files
2:		Partition $d$ into $F={\{Trai{n}_j,Tes{t}_j\}}_{j=1}^4$				// 4-fold cross-validation
3:		FOR EACH pair of train and test data folds (f_tr, f_te) in F:				// f_tr = 75% of experimental data
4:			FOR EACH kinetic model mod in K:			// Described in Equations (2)–(5)
5:				FOR EACH optimizer opt in O:		// IPA, PSO and UMDA
6:					$\left({\mathsf{\theta }}^{*}, P{I}_{train}\right)\leftarrow opt\left(mod,f_{tr},{\chi }^2\right)$	// Calibrate and evaluate kinetic models
7:					$P{I}_{test}\leftarrow \mathrm{eval}\left(mod,{\mathsf{\theta }}^{*},f_{te}\right)$	// Compute prediction error (MSE_test)
8:				END FOR
9:			END FOR
10:			Save_file(${\mathsf{\theta }}^{*}, P{I}_{train},P{I}_{test},mod,opt,f\_tr,f\_te)$			// Store partial results
11:		END FOR
12:		Compute Average and Standard deviation on PI				//Compute Statistics and plots
13:	END FOR
14:	Output the best optimizer, kinetic model, and corresponding optimal parameters as a result

3. Results

3.1. The Best Kinetic Model

Table 2 shows the best kinetic model found by each optimizer and its corresponding optimal parameter values ${\mathsf{\theta }}^{*}=\left(V_{max}^{*},K_M^{*},I^{*},K^{*}\right)$. It also presents the average and standard deviation of the quality measures R², AIC, ${\chi }^2$, and MSE. As can be observed from Table 2, the least complex and best-fitting kinetic model was the Michaelis–Menten Traditional model for most of the experiments.

For the optimizers, the UMDA presented the worst performance in terms of both error and quality measures. In addition to its bad performance, the UMDA showed premature convergence on different runs. The quality measure of the UMDA was insignificantly improved as the number of iterations and the population size increased. Because of this behavior, we concluded that the UMDA is easily trapped in local minima; consequently, we hypothesize that the quality measure ${\chi }^2$ is not convex.

Surprisingly, the Interior-Point algorithm dominated in determining the best-fitting model for most of the experiments. Furthermore, the MSEs in training and testing were similar. However, as can be observed in the results of Experiments 08 and 11 (Table 2), different parameter sets lead to the same quality measures. From this result, it can be concluded that ${\chi }^2$ is a non-injective function of $\mathsf{\theta }$. Despite this result, the Interior-Point algorithm was robust enough to carry out the optimization process. The good results of the PSO in conjunction with the bad performance of the UMDA provide evidence that the quality measure ${\chi }^2$ is not smooth but presents a certain degree of roughness. Therefore, optimizers that can avoid local minima like the PSO or the Interior-Point should be preferred.

Experiments 04 and 06 are two representative cases of the numerical results reported in Table 2. Figure 1 illustrates a scenario where the Nelder–Mead, PSO, and Interior-Point optimizers achieve similarly high performance, while the UMDA leads to an unsatisfactory fit, and Nelder–Mead leaves the feasible regions in the parameter search space. Figure 2 depicts one example of the two cases (Experiments 06 and 07) with the lowest performance. The graphical comparisons for the rest of the experiments are provided at the Supplementary Materials section.

3.2. Computation Time for Kinetic Parameter Optimization

A comparison of the computation time required for each optimizer to find the best kinetic parameters is depicted in Figure 3. Here, the computation time not only depends on the optimizer but also on the kinetic model being optimized. For the first three models (Michaelis–Menten, Non-competitive inhibition, and competitive inhibition), the UMDA and PSO required between 10 and 100 times the average amount of time required by the Nelder–Mead and Interior-Point algorithms to converge to their best solutions. Notably, for the Substrate Inhibition model, the UMDA and PSO (and sometimes the Nelder–Mead) optimizers presented convergence problems. These three algorithms required around 100–1000 times the amount of time compared to the time required by the Interior-Point method to converge. These problems could be due to two reasons: (i) All methods except for the Interior-Point are completely stochastic; thus, some of the random candidate solutions can be stuck in valleys of the response surface, thereby hindering the convergence of the optimizer; (ii) The Interior-Point was provided with an initial candidate solution that favors its convergence.

The optimizers were run on a computer with an Intel^® (Intel Corporation, Santa Clara, CA, USA) Core™ i7 CPU @ 2.60 GHz processor, 16GB RAM, and a 256 GB hard disk. The total computation times to calibrate the first three models with the available experimental data were 0.19, 1.05, 13.06, and 33.37 min for the Nelder–Mead, Interior-Point, PSO, and UMDA algorithms, respectively. For the Substrate Inhibition model, the total computation time was 0.13, 27.54, 592.41, and 832.15 min for the Interior-Point, Nelder–Mead, UMDA, and PSO, respectively. In conclusion, the Interior-Point algorithm is the fastest method among the four optimizers analyzed in this work for kinetic parameter optimization. Nelder–Mead is also fast, but it is not recommended because it converges to unfeasible regions and is sensitive (like the UMDA and PSO) to the typical randomized initialization of metaheuristics for complex models such as Substrate Inhibition.

3.3. Comparison between Micro-Reaction and Bench Scale Reactor

As a final experiment, the values of the kinetic parameters for the optimized experiment in the micro-reaction were obtained. The values of the parameters V_max, K_m, and K were 1.06 mg/Ml × s, 8.24 mg/mL, and 0.27, respectively. The data were obtained using the IPA method and the traditional kinetic model according to the conclusions derived from Table 2. Figure 4 depicts the trajectory of the optimized micro-reactor compared to the three bench scale reactors. According to the results presented in [12], a micro-scale experiment was developed using the optimal parameters obtained by the central composite design, and these parameters were also used to test the experiment on a scale of 0.5 L.

Notably, the difference in size between the two reported experiments is 500 times higher. Furthermore, the reported micro-scale experiments were tested by orbital shaking at a speed of 700 rpm and a very stable temperature control with variations of ±1 °C. The difference with the 0.5 L reactor was the use of propeller agitation (150 rpm), ensuring that all the reactor contents were well mixed to ensure spatial homogeneity. Further studies should be done to ensure spatial homogeneity, but this goal is beyond the scope of this work.

According to the results, the Traditional model with its optimized parameters was able to predict the behavior in the micro-reaction, showing $R^2=0.99$ and $MSE=0.06$. In addition, the same model and parameters predicted the behavior of the reaction on a bench scale for three different reactors, achieving performance of ($R^2=0.96$; $MSE=0.30)$, ($R^2=0.93; MSE=0.54)$, and $\left(R^2=0.97;MSE=0.22\right)$ for reactors 1, 2, and 3, respectively. Therefore, the behavior of depolymerization on a bench scale can be predicted with acceptable errors using micro-reaction experiments.

4. Discussion

The results in Table 2 present the best-fit model found by each optimizer. To clarify the differences between the best-fitted model and the rest of the models, Figure 5 provides a graphical comparison of the four M–M models calibrated using the four optimization strategies for Experiment No. 13. Figure 5a–c depict the fit obtained by both IPA and PSO, which reached the same R² of 0.93, but the AIC value for the IPA method was lower (−3.06). Consequently, we observed a better fit in Figure 5a than in Figure 5b. The UMDA optimizer, shown in Figure 5c, presents an R² of 0.86 and an AIC of 1.02. UMDA presents a poor fit compared to the three previous methods, but it also adequately recapitulates the trajectory shown by the experimental data.

The parameters V_max and K_M characterize the enzymatic reactions that are described with the kinetics of M–M. V_max depends on the total concentration of the enzyme, but K_M does not [14]. K_M is equal to the concentration of the substrate at half V_max. It can also be associated with the affinity of the enzyme for the substrate [3]. As reported by Cho et al. [3], a low K_M value indicates a sizeable binding affinity, as the reaction approaches V_max more rapidly. Furthermore, a high K_M indicates that the enzyme does not bind as efficiently with the substrate, and V_max will only be reached if the substrate concentration is high enough to saturate the enzyme [3]. The $V_{max}$ and $K_m$ parameters were searched in (0, 8.25) and (0, 8.25), respectively. Parameter I (the inhibitor) was calculated using a range between 0 and 8.25. This range was employed because it is well known that there is an intermediate product called cellobiose during the production of glucose that exerts an inhibitory effect over biomass depolymerization reactions [2]. The maximum amount of glucose that can be ideally produced is 8.25 mg.

As previously explained, Table 2 shows the results obtained for the parameters V_max, K_M, I, K, ${\chi }^2$, and R², as well as the error value. To explain the kinetic parameters (V_max and K_M), Experiment 11 and 12 were used (Figure 6). Experiments 11 and 12 have the same characteristics with the exception of temperature: Experiment 11 @ pH = 5; T = 65 °C; enzyme–substrate ratio = 0.9 µL of enzyme/g of substrate; Experiment 12 @ pH = 5; T = 45 °C; enzyme–substrate ratio = 0.9 µL of enzyme/g of substrate. The V_max and K_M values for Experiment 11 were 3.83 mg/mL·s and 7.26 mg/mL, indicating that the V_max was quickly reached due to the affinity of the enzyme with the substrate, probably associated with the effect of the temperature. It is well known that a high temperature favors enzymatic reactions. The results agree with the low K_M value obtained, which is associated with a high affinity between the enzyme and the substrate.

On the other hand, Experiment 12 shows a low V_max (0.19 mg/mL·s) and a high K_M value (8.25 mg/mL). The high K_M value indicates that the enzyme–substrate affinity is low; however, it allows the experiment to attain a higher depolymerization value. Both experiments show a high inhibitor value, indicating that the reaction is being inhibited. This result could be associated with a high accumulation of cellobiose during the reaction favored by temperature [2]. A graphic compilation of all the methods obtained is shown in the complementary information.

5. Conclusions

For the metaheuristic optimizers, we conclude that the Interior-Point algorithm is effective for solving the inverse modeling problem of kinetic parameters for the depolymerization of biomass. If a priori information, such as the initial candidate solution, is available, the Interior-Point should be the first option; if not, the PSO may be considered. The default values were shown to be sufficient for all the optimizers to provide acceptable solutions. However, pure metaheuristic algorithms (like the UMDA and PSO) might further improve convergence efficiency by tuning their hyper-parameters.

The quality measure employed to evaluate the best-fitting kinetic model is likely non-convex and rough. Although many other quality measures could be used, the Weighted Least Squares metric ${\chi }^2$ was enough to guide the optimizer exploration and correctly identify the best model for experimental data, at both the micro and macro scale.

The methodology applied in this study can be adopted as a starting point for solving forthcoming inverse modeling problems. The use of different metaheuristic and hybrid optimizers can reveal the nature of the quality measures, thereby easing the implicit optimization problems.

In this study, we found an effective computational procedure to determine a better calibration for mathematical models of kinetic polymerization process. However, new experiments should be considered in the future to study micro-reactions, compare these with larger scales and test multiple operating parameters. In addition, important considerations such as heat and mass transfer should be studied but are beyond the scope of this investigation.

Finally, according to our observations, the models with the best fit are not necessarily consistent with a fixed set of kinetic parameters, K_M and V_max. Therefore, to inform the best decision for the process it is necessary to discern which K_M and V_max values best describe the behavior of the evaluated reaction.