Modeling and simulation of the enzymatic degradation of 2,4,6-trichlorophenol using soybean peroxidase

The enzymatic degradation of organic pollutants is a promising and ecological method for the remediation of industrial effluents. 2,4,6-Trichlorophenol is a major pollutant in many residual waters, and its consumption has been linked to lymphomas, leukemia, and liver cancer. The goal of the present work is to comprehend the enzymatic degradation of 2,4,6-trichlorophenol using soybean peroxidase. Different assumptions for the kinetic model were evaluated, and the simulations were compared to experimental data, which was obtained in a microreactor. The literature pointed out that the bi-bi ping-pong model represents well the kinetics of soybean peroxidase degradation. Since it is a complex model, some reactions can be considered or not. Six different possibilities for the model were considered, regarding different combinations of the generated enzyme forms that depend on the hypotheses for simplifying the model. The adjustment of the models was compared based on different metrics, such as the value of the objective function, coefficient of determination and root-mean-square error. The process modeling was obtained by the mass balance of all the reaction components, and all the simulations were performed in MATLAB ﻿R2015a. Reaction parameters were estimated based on the weighted least squares between the experimental data set and the values predicted by the model. The results showed that the data were better adjusted by the model that considers all the enzyme forms, including enzyme inactivation. Therefore, a better comprehension of the reaction mechanism was achieved, which allows a more precise reactor project and process simulation.


Introduction
The environmental impact caused by industries occurs mainly due to several pollutants present in industrial effluents, such as 2,4,6-trichlorophenol (TCP). This substance is a common chemical intermediate and is usually generated as a by-product of the chlorination of water and combustion processes. It can also appear in the water supply as a result of the hydrolysis of chlorinated benzene contaminants (Gowdal et al. 1985). TCP has already been found in large quantities in aquatic environments in many countries. The exposure to humans can thus occur by the ingestion of such waters (Zhang et al. 2018).
TCP is an organochlorine compound used in the manufacturing of antiseptics, pesticides in general (bactericides, fungicides, herbicides, insecticides, germicides) and wood and glue preservatives. The consumption of TCP has been linked to lymphomas, leukemia, and liver cancer; it is classified as Group B2 (probable human carcinogen) by the United States Environmental Protection Agency (Gowdal et al. 1985). Therefore, there is considerable interest in studying alternatives for the decomposition of this pollutant.
The removal of undesirable components from water may be classified under physical, chemical, and biological treatments. The physical treatments include unit operations such as solid-liquid separations, screening, mixing, and adsorption. The chemical ones encompass flocculants, chemical precipitation, oxidation, and advanced oxidation processes. The biological ones, in turn, involve biological oxidation, denitrification, anaerobic oxidation and activated sludge processes (Tchobanoglous et al. 2014). However, these treatments may not be suitable for some recalcitrant effluents, such as TCP and other complex organic molecules. Enzymatic reactions should be used in these special cases (Husain 2010). There are several recent reviews about the use of enzyme-based technologies for bioremediation, which cover all the treatments of soil and water polluted with hazardous environmental pollutants using enzymes (Sharma et al. 2018).
Enzymes are organic substances that act as biological catalysts, decreasing the activation energy of a chemical reaction under mild environmental conditions and with high efficiency. They can be obtained from renewable resources, which is an advantage over traditional chemical catalysts. They act specifically on a substance or a group of substrates, so they can be used to modify particular components. Many authors state that soybean peroxidase (SBP) shows good potential for the degradation of organic pollutants (Nissum et al. 2001;Al-Ansari et al. 2011;Steevensz et al. 2014).
SBP is a Fe(III)-heme enzyme that belongs to the oxidoreductase family, which oxidize organic substrates using hydrogen peroxide as an electron acceptor. The heme-containing peroxidases in general can be divided into 2 groups: one group found only in animals and other group found in fungi, bacteria, and plants. The second one is further divided into 3 classes (Sharma et al. 2018): • Class I intracellular enzymes, including cytochrome-C peroxidase produced by yeast, ascorbate peroxidase (APX) produced by some species of plants and bacterial catalase peroxidases. • Class II secreted fungal enzymes, including lignin peroxidase (LiP) and manganese peroxidase (MnP). • Class III plant secreted peroxidases such as SBP from soybean seed hulls and horseradish peroxidases (HRP) from horseradish plants.
Although peroxidases have been proposed for applications in several fields, there are few industrial processes that utilize peroxidases. The commercial applications of these enzymes are reduced to diagnosis and research. Unfortunately, the economic feasibility of peroxidase-based processes has not always been assessed (Torres and Ayala 2010). In this context, in recent years SBP has been the focus of research, due to its great availability and many other factors. The advantages of SBP compared to other enzymes are: it is cheaper, since it is a byproduct of the soybean processing industry; it is easier to obtain, since a complex purification process is not necessary; it can operate at relatively high temperatures and over a wide pH range; it is less susceptible to irreversible inactivation by hydrogen peroxide and chemical denaturation; it is the most effective peroxidase (Al-Ansari et al. 2011). SBP removal from soybean seed hulls can be processed in a sustainable and low-cost way (Steevensz et al. 2014;Calza et al. 2016), so it can be the basis for a very promising and environmentally friendly industrial effluent remediation method.
There are several studies of SBP uses for industrial wastewater treatment in the literature, as well as SBP extraction and purification processes (Calza et al. 2016;Silva et al. 2013;Graham et al. 1991;Gacche et al. 2003;Bassi and Gijzen 2004;Ghaemmaghami et al. 2010;Steevensz et al. 2013;Tolardo et al. 2019). Over the years some of the limitations and advancements in the field have been presented, considering aspects such as cost, use of additives for increased enzyme economy, enzyme recycling and studies already completed on industrial wastewaters. There are many more studies on synthetic wastewater compared with a limited number of real effluents, and the feasibility to treat various industrial wastewaters has not yet been proven (Steevensz et al. 2014).
The ability of SBP in the enzymatic degradation of TCP has been previously studied (Calza et al. 2014). Hypothetical schemes for the catalytic mechanism of TCP oxidative dechlorination by the peroxidase/H 2 O 2 system have been presented, in which the final metabolite is obtained by two one-electron oxidations (Torres and Ayala 2010;Ferrari et al. 1999). Furthermore, the degradation products can be identified directly (Laurenti et al. 2002(Laurenti et al. , 2003, since the products usually show different toxicity than the reagents (Silva et al. 2013). The most recent study of phenolic compound removal using SBP was in the treatment of aqueous solutions containing pentachlorophenol (Tolardo et al. 2019), based on previous studies of TCP degradation (Calza et al. 2014). Besides, it was observed that the removal of pentachlorophenol was not complete, maybe because of the interference of intermediates or reaction products. In this context, TCP degradation is still one of the most consistent in this area, and further research is needed in order to determine other pollutants for which this approach can be applicable (Tolardo et al. 2019).
Enzymes of the peroxidase class such as SBP and HRP catalyze the degradation of compounds such as phenols and amines following a mechanism called modified bi-bi pingpong, involving a sequence of reactions with intermediate oxidized components (compound I, II and III) (Al-Ansari et al. 2011), as will be detailed in the next section. The representation of this enzymatic kinetic model, as other chemical reaction models, can be made employing mathematical systems composed of differential and/or algebraic equations of the process variables, such as the concentrations, so that predictions based on the model response can be obtained. These predictions allow the simulation, design, and optimization of the process. Nicell (1994) developed steady-state, fully transient and pseudo-steady-state models for the polymerization and precipitation of 4-chlorophenol catalyzed by HRP in the presence of hydrogen peroxide. The author implemented several simplifying assumptions, finding satisfactory results for the fully transient and pseudo-steady-state models incorporating inactivation mechanisms. Besides, the computation time to solve the model was excessive for the fully transient model. Nicell and Wright (1997) studied the effect of hydrogen peroxide concentration with an assay based on phenol as substrate and developed a steady-state kinetic model based on experimental data obtained for reaction with both SBP and HRP. The model indicated that SBP tends to form more compound III and is catalytically slower than HRP during phenol oxidation. Buchanan and Nicell (1997) developed a pseudo steady-state kinetic model for the system HRP/H 2 O 2 /phenol, improving significantly its predictive ability by incorporating enzyme inactivation mechanisms. The kinetic constants were adjusted using a series of experimental data sets. The rate of enzyme inactivation due to end-product polymer was modeled as being proportional to the rate of polymer formation. They also confirmed the apparent rate of compound III decomposition as being first order with respect to compound III concentration.
However, currently there is a lack of more work in the literature with an approach focused on the mathematical modeling of these enzymatic degradation processes, and this is one of the main objectives of this work. In addition, this work aims to contribute to the literature in order to conduct the enzymatic degradation reactions of TCP in a microreactor, which is a suitable equipment for fast reactions such as enzymatic ones. The main advantages of using microreactors are the strict control of reaction conditions, safety in operation, improvement in mass and heat transfer and energy efficiency (Laurenti and Vianna-Jr 2016). Recently, a study on the degradation of TCP in a microreactor using the enzyme SBP was published by Costa et al. (2020). The present work aims to focus more deeply on the mathematical modeling of the ping-pong model, aiming to complement Costa et al.'s work. For that, different scenarios of kinetic models were evaluated, considering the bi-bi ping-pong model, and the simulations were compared to experimental data, obtained in a microreactor.

Fundamentals
The peroxidase reactions follow a bi-bi ping-pong mechanism, as shown in Fig. 1 (Costa et al. 2020): This mechanism involves two-electron transfer to the first substrate ( H 2 O 2 ) that creates a first enzyme oxidized intermediate called 'compound I' (CpI) (Eq. 1), and two successive one-electron oxidations of the second substrate (TCP) to return the enzyme to its reduced state via a second enzyme intermediate called 'compound II' (CpII) (Eqs. 2, 3) (Nicell 1994).
In Fig. 1, it can also be observed that peroxidase can be deactivated due to an excess of H 2 O 2 , generating a third enzyme oxidized intermediate called 'compound III' (CpIII) (Eq. 4), which is catalytically inactive but can be decomposed back to the enzyme native form (Eq. 5). CpIII can also be reduced to CpI (Eq. 6), which in turn can be oxidized back to the enzyme native form (Eq. 7) (Nicell 1994).
It is important to emphasize that all forms of the enzyme can interact with the substrates forming enzyme-substrate complexes (Fersht 1999), which was not considered in reactions 1-7. So, the kinetic constants k 1 , k 2 , k 3 , k app , k a , k b and k c are the apparent constants of the reaction. However, there is no evidence of enzyme-substrate complex formation for Peroxidase mechanism and suicide pathways reactions with phenolic compounds, that is, only the interaction of H 2 O 2 with the enzyme oxidized intermediates might form enzyme-substrate complexes (Nicell 1994). Besides, only the existence of the complex CpI⋅H 2 O 2 has been confirmed, and in this case the interaction of CpI with H 2 O 2 (Eq. 7) can be replaced by Eq. (8) (Nicell 1994): As illustrated in Fig. 1, the CpI⋅H 2 O 2 complex can also generate an enzyme inactive form called verdohemoprotein (P-670). Besides, all the enzyme inactive forms will be represented in the work by Eq. (9), which also involves other factors such as free-radical inactivation and adsorption or entrapment on the end-product polymers. Once the reaction has gone to completion, the number of phenolic compounds precipitated per molecule of the enzyme in solution is constant. So, SBP can be treated as a pseudo-substrate of the reaction, as represented below (Nicell 1994): K s is the number of TCP molecules removed per enzyme molecule after complete reaction, also known as the turnover number of an enzyme, and is calculated by: Then the amount of permanently inactivated enzyme can be calculated by: The kinetic constants of the reactions listed in Table 1 will be used as a reference in the model simulation. All of these rate constants were obtained for HRP enzyme at 25 • C and pH 7. In this case, although the SBP and HRP enzymes are very similar, a different behavior in these aspects cannot be excluded in the absence of experimental evidences.

Experiments
To develop an appropriate model, experimental data reported by Costa et al. (2020) were used as reference. The reactions were carried out in a Syrris 250 μ L microreactor in the laboratory of the Chemical Engineering Department of the Polytechnic School of the University of São Paulo, at pH 5.4 and temperature of 25 • C. 3 experimental data sets were obtained, considering the SBP concentrations of 0.001 mg/mL, 0.002 mg/mL and 0.0005 mg/mL. All measurements were performed in triplicate. The concentrations of TCP and H 2 O 2 were 0.25 mM and 0.6 mM, respectively. For the inhibition of the reaction process, 37% hydrochloric acid was used. At the end of the reaction process, a 0.1 mL sample was collected at the microreactor outlet, which was analyzed off-line by high performance liquid chromatography (Shimadzu LC-20AD Prominence). The HPLC system operated at: • Column: Phenomenex Luna C18; • Eluent: 60% acetonitrile + 40% demineralized water with 0.2% acetic acid; • Flow: 0.6 mL/min; • Temperature: 40 • C; • Estimated analysis time: 10 min; • Detector, = 220 nm.

Process modeling
The process modeling was obtained from the mass balance of all the species considering six scenarios with different hypotheses for simplifying the bi-bi ping-pong model, as below: • Model 1 Formation of CpI and CpII (reactions of Eqs. 1, 2 and 3), in order to estimate parameters k 1 , k 2 , and k 3 ; • Model 2 Reactions considered in Model 1 including the formation of CpIII (reactions from Eqs. 4-7), in order to estimate parameters k 1 , k 2 , k 3 , k app , k a , k b , and k c ; • Model 3 Reactions considered in Model 2 including the formation of CpI⋅H 2 O 2 (reaction of Eq. 8 instead of Eq. 7), in order to estimate parameters k 1 , k 2 , k 3 , k app , k a , k b , k 4 , k −4 , and k 5 ; • Model 4 Reactions considered in Model 1 including enzyme inactivation (Eq. 11), in order to estimate parameters k 1 , k 2 , and k 3 ; • Model 5 Reactions considered in Model 2 including enzyme inactivation (Eq. 11), in order to estimate parameters k 1 , k 2 , k 3 , k app , k a , k b , and k c ;  Arnao et al. (1990) • Model 6 Reactions considered in Model 3 including enzyme inactivation (Eq. 11), in order to estimate parameters k 1 , k 2 , k 3 , k app , k a , k b , k 4 , k −4 , and k 5 .
The system of differential equations generated by Model 1 is presented below: Steady state can also be assumed for all the enzyme forms (they quickly reach this condition): From Eqs. (13), (14) and (18), CpI and CpII concentrations can be obtained: The total (initial) amount of enzyme is given by:  In Models 4, 5, and 6, the equations are almost the same as in Models 1, 2, and 3, respectively, with only the term [SBP] inactive being considered in Eqs. (22), (29) and (32)

Model simulation
Based on the equations described in the previous section, the procedure for estimating kinetic constants was performed, which is basically an optimization procedure with restrictions, with the values reported in Table 1 being considered as the initial estimates. The following restrictions for the parameters were considered (Nicell 1994): All computational calculations were performed with the software MATLAB version R2015a. A parameter estimation procedure was implemented with a main routine that reads all the experimental reaction data, in order to carry out the optimization problem using fmincon considering the interior point method, tolerances of 10 −10 and all other settings as default. The parameter estimation procedure also takes into account an objective function (OF) and the respective mathematical model consisting of a DAE, which were implemented in two respective subroutines.
The OF was minimized by the weighted least squares between the calculated values predicted by the model ( û i ) and their respective experimental data set ( u i ) for all the number of points (n) and for all the experiments (m), as shown in Eq. (35). This equation also takes into account the standard deviation of the experimental data ( ), which in this case was calculated for each set of experiments performed according to Eq. (36), where u i is the average experimental data.
The model solution was calculated by the backward differentiation formula (BDF) method (ode15s MATLAB's self-implemented solver). Relative and absolute tolerances of 10 −8 and 10 −10 were considered, respectively, and other configurations as default.
The model fit was evaluated from the plots of model prediction and experimental data, and the calculated values of root-mean-square error (RMSE) and the adjusted coefficient of determination (R 2 adjusted ), Eqs. (37) and (38), respectively, where p is the number of parameters.
The smaller the numerical value of RMSE and the higher the numerical value of the adjusted R 2 , the better the model fit. Adjusted R 2 was used instead of R 2 because it is more appropriate for evaluating model fit and comparing alternative models due to the inclusion of a penalty for additional parameters that do not contribute to model fit. However, this criterion is not always reliable for non-linear models, so it was used in conjunction with other criteria.
The model fit could also be evaluated by the value of the OF obtained from Eq. (35) in each case, because a lower value indicates a better minimization of the OF and consequently a greater success in the optimization procedure.

Results and discussion
As informed in Sect. 2.3, six possibilities were carried out in the model simulation, in which different hypotheses were made to simplify the bi-bi ping-pong model. It should be noted that several possibilities were tested before obtaining the model presented, including the Michaelis-Menten model, considering only TCP as a substrate, as well as other models considering different combinations of intermediate reaction complexes, whether or not reverse reactions, etc. A very large set of parameters was considered, but this did not allow obtaining a single set of parameters in the estimation procedure. The strategy then was to start from the simplest model and gradually add the assumptions.
The parameter estimation procedure was performed as described in the previous section. Table 2 presents the estimated parameters and Fig. 2 presents the model prediction and experimental data for the models. It is important to emphasize that the results of Model 5 have already been published by Costa et al. (2020), and the present work is a general analysis involving the entire process of formulating the model.
In Fig. 2, all measurements of the experimental data were performed in triplicate, with average errors of ± 0.00033 mM, ± 0.00024 mM, and ± 0.00027 mM for the experiments with enzyme concentrations of 0.0005 mg/mL, 0.001 mg/mL, and 0.002 mg/mL, respectively. These values were so small that it was not possible to display the error bars on the experimental points. These experimental results show the good repeatability and reproducibility of the experiments in the microreactor, which is an adequate equipment for studies of the kinetics of fast reactions such as enzymatic ones (Costa et al. 2020). Still in Fig. 2, in general a coherent fit of the experimental data for the bi-bi ping pong model can be observed. It can be observed that Models 1, 2, and 3 presented almost equal adjustments, as well as Models 4, 5, and 6, for the three sets of experiments with different enzyme concentrations. However, the adjustment of these two groups of models was considerably different, which is explained by the addition of the enzyme inactivation in the latter. Models 4, 5, and 6 visually present curves closer to the experimental points, except for experiments with an enzyme concentration of 0.0005 mg/ mL at longer reaction times. Anyway, all the curves respect the expected behavior of a greater degradation of TCP at higher enzyme concentrations.
The results presented in Fig. 2 can be explained quantitatively by the values of OF, RMSE, and adjusted R 2 presented in Table 2, where the adjusted R 2 values refer to the experiments at the three enzyme concentrations. In general, it can be observed that the models that consider the enzyme inactivation (Models 4, 5, and 6) presented a considerably better fit than the models that do not consider the enzyme inactivation (Models 1, 2, and 3), because they presented 2.0 × 10 7 2.0 × 10 7 2.0 × 10 7 2.0 × 10 7 2.0 × 10 7 2.0 × 10 7 M −1 s −1 k 2 3.4 × 10 6 2.52 × 10 6 3.36 × 10 6 6.38 × 10 6 6.79 × 10 6 8.19 × 10 6 M −1 s −1 k 3 8.8 × × 10 5 9.72 × 10 5 8.85 × 10 5 2.47 × 10 6 2.75 × 10 6 3.01 × 10 6 M −1 s −1 k app -7.92 × 10 −1 8.99 × 10 −2 -6.6 × 10 0 1.26 × 10 −1 M −1 s −1 k a -4.28 × 10 −2 6.48 × 10 −2 -1.09 × 10 −2 9.56 × 10 −3 higher values of adjusted R 2 and lower values of OF and RMSE. Therefore, it can be inferred that considering the inactivation of the enzyme resulted in a great gain for the adjustment of the model. Comparing the adjustment for the three sets of experiments using the different values of adjusted R 2 in each model, in cases where the enzyme inactivation was not considered (Models 1, 2, and 3), a better fit for the first set of experiments can be observed when there is a lower concentration of enzyme, while in the models considering enzyme inactivation (Models 4, 5, and 6) a better adjustment is observed for the highest concentrations of enzyme. This could be explained by the simple fact that, in the first set of experiments, a smaller amount of initial enzyme would result in a smaller amount of inactivated enzyme, so this consideration would not be so relevant in this specific case.
Comparing the effect of considering the formation of CpIII (Models 2 and 4) and CpI⋅H 2 O 2 (Models 3 and 6), in turn, it can be observed that the consideration of a more complex model results in lower values of adjusted R 2 , which would mean a worsening of the model's fit a priori. However, it can be seen that the values of OF and RMSE reflect an improvement in the adjustment. It is known that adjusted R 2 includes a correction for the additional parameters added to the model, but perhaps this penalty was greater than it should be, so that the number of experimental points was not sufficient to compensate. That is why this criterion is not always reliable for evaluating nonlinear model adjustments. So, regarding only the OF and RMSE criteria, in the models without considering enzyme inactivation, the adjustment worsens from Model 1 to Model 2 and improves from Model 2 to Model 3, while in models considering enzyme inactivation the adjustment improves sequentially from Model 4 to Model 6. Because the models with enzyme inactivation showed an improvement in the model in general and the OF and RMSE criteria were more reliable than adjusted R 2 in this case, the global analysis can infer that Model 6 presented more consistent results among the considered hypotheses. This result corroborates the considerations made regarding the formation of compound CpIII and CpI⋅H 2 O 2 complex. However, considering that few data were used and the experimental errors involved, it cannot be concluded with certainty that these assumptions were completely adequate, and it is necessary to verify further by performing new experiments.
Simulations of all the enzyme forms and the enzymesubstrate complexes were performed, as can be seen in Figs. 3, 4, 5, 6, 7 and 8. Although the enzyme can be regenerated to its initial state, all the models indicate that its concentration always decreases throughout the reaction, since it can generate its oxidized forms or be inactivated. In Models 1, 2, and 3, all other forms of the enzyme grow slightly until they reach a steady state, according to model considerations. In Models 4, 5, and 6, all forms of the enzyme decrease throughout the reaction until reaching steady state, except the inactive form, which substantially increases. The models also indicate that the most significant enzyme form in the reaction is CpII, followed by CpI, CpIII and CpI⋅H 2 O 2 , respectively, when applicable.
Besides, it is worth mentioning that these results are based exclusively on the model considered and the set of experimental data obtained. An identification of the intermediate forms of the enzyme should be made to clearly state the mechanism of the reaction.

Conclusions
This study presents a significant example on the use of the SBP enzyme in the degradation of chlorinated phenolic compounds such as TCP. Mathematical modeling and simulation were developed in MATLAB considering six different hypotheses based on the bi-bi ping pong model in order to match the experimental data. These assumptions consist of considering or not the enzyme inactivation, and the formation of the additional intermediate compound CpIII and the CpI⋅H 2 O 2 complex. The experiments were carried out in a microreactor with very good repeatability and reproducibility of the experimental results.
A parameter estimation procedure was implemented in order to obtain the kinetic constants of the reactions using the MATLAB's self-implemented solver fmincon, which performs the optimization procedure considering the interior point method, and ode15s, which solves the systems of algebraic-differential equations by the backward differentiation formula (BDF) method. The objective function was the weighted least squares between the experimental data set and the values predicted by the respective model. All curves obtained by the models showed the expected behavior of a greater degradation of TCP at higher enzyme concentrations. The results of the model fit evaluation were compared based on the criteria of the OF values, RMSE, and adjusted R 2 . However, the latter was not very suitable for evaluating the model adjustment, because it adds a penalty for the parameters included in the model, so that the low number of experimental points was probably not sufficient to compensate.
Finally, a global analysis concluded that considering all the enzyme forms (CpI, CpII, CpIII, and CpI⋅H 2 O 2 ), as well as the inactivation model, provides more consistent results. The models also indicate that the most significant enzyme form in the reaction is CpII, followed by CpI, CpIII, and CpI⋅H 2 O 2 , respectively. However, more experiments must be carried out in order to identify these intermediate forms of the enzyme and, consequently, establish the reaction mechanism more clearly.
Thus, it can be concluded that the results presented in this work can serve as the basis for a better understanding of the reaction mechanism, which allows a more accurate reactor project and process simulation.