Towards a Novel Computer-Aided Optimization of Microreactors: Techno-Economic Evaluation of an Immobilized Enzyme System

Immobilized multi-enzyme cascades are increasingly used in microfluidic devices. In particular, their application in continuous flow reactors shows great potential, utilizing the benefits of reusability and control of the reaction conditions. However, capitalizing on this potential is challenging and requires detailed knowledge of the investigated system. Here, we show the application of computational methods for optimization with multi-level reactor design (MLRD) methodology based on the underlying physical and chemical processes. We optimize a stereoselective reduction of a diketone catalyzed by ketoreductase (Gre2) and Nicotinamidadenindinukleotidphosphat (NADPH) cofactor regeneration with glucose dehydrogenase (GDH). Both enzymes are separately immobilized on magnetic beads forming a packed bed within the microreactor. We derive optimal reactor feed concentrations and enzyme ratios for enhanced performance and a basic economic model in order to maximize the techno-economic performance (TEP) for the first reduction of 5-nitrononane-2,8-dione.


Introduction
The application of miniaturized continuous reactors has increased rapidly over the past few decades. Microreactors in particular possess great potential for synthesizing various compounds, including substances of pharmaceutical interest [1,2]. Flow reactors [3], and especially micro flow reactors, benefit from increased control of reaction conditions, enhanced mass transport, increased safety, the possibility of automation and easier scaleup [3][4][5][6], as compared to batch reactors. Simultaneously, they enable sustainable and green processes, due to potential recycling of substrates and less chemical waste [7]. Many enzymes require an expensive cofactor, which is consumed [8,9]. In order to create a more cost-efficient process, additional reactions are utilized to regenerate the cofactor [10]. Ideally, this is done in close proximity, e.g., in the same reactor compartment, provided the required reaction conditions such as temperature and pH value for both enzymatic reactions are similar. This can be achieved by enzyme immobilization on macroscopic supports, which can be placed at a desired position in the reactor [11]. To this end, enzymes can be (non-)covalently attached to matrices such as beads, enabling a potentially increased enzyme stability and recycling of the biocatalyst [12,13]. Reactors with immobilized multienzyme systems enabled by different immobilization techniques are increasingly gaining attention due to their versatility and application in flow chemistry [14][15][16][17][18]. Applications include packed beds with micrometer-and nanometer-scale particles, wall-coated systems and matrix-entrapped enzymes. We extend this list with an investigation of an overflown packed bed reactor with micrometer-scale particles.
Understanding of such complex multi-enzyme systems is vital to utilize them as efficiently as possible. Thereby material usage, productivity and production cost are optimized Figure 1. Reaction scheme of chiral reduction of NDK to the corresponding hydroxyketone (HK) and the diol using immobilized Gre2. The cofactor NADPH is regenerated during oxidation of glucose with glucose 1-dehydrogenase (GDH) [36].
In a first approach, we immobilized Gre2, which is genetically fused to a Halo-Tag variant, on Streptavidin-coated beads (Dynabeads M280), using a HaloTag-PEG-Biotin linker and kept GDH in the aqueous solution. The magnetic beads (MBs) were held inside the reactor with magnets placed below the reactor, forming a packed bed. Burgahn et al. clearly indicated a need for higher regeneration rates inside the packed bed [36]. To achieve a sufficient concentration of NADPH within the packed bed layer, GDH is immobilized on separate beads bearing reactive epoxy groups (Dynabeads M270 Epoxy) by employing a His-SpyCatcher linker. The tetrameric GDH is equipped with SpyTag fusion proteins and reacts with the linker in a covalent manner, thereby efficiently binding GDH to the beads.
The reactor encompasses a single microchannel with 200 µm height, 2 mm width and a particle bed with a length of 48 mm. The MBs are held within the reactor with eight magnets and form a packed bed at the bottom of the channel, schematically shown in Figure 2. This system was previously investigated by Burgahn et al., additionally applying a mathematical model to create a deeper understanding of the transport and reaction phenomena within the reactor. This model is also used here and is addressed as a basic Matlab model in the following. It can be regarded as a one-dimensional model, as it divides the reactor into sequential segments in the flow direction. Each segment consists of a packed bed and an ideal mixed free volume above the packed bed. Reactions and mass transport transversal to the flow direction in the packed bed are described by a set of ordinary differential equations, which are solved numerically assuming a finite mass transfer between the packed bed and the free volume above according to a Sherwood number correlation. Mass transport within the packed bed is calculated using effective diffusivities based on porosity and tortuosity of a compact packing of uniform, spherical particles. The particle bed is assumed to be a homogeneous package with constant height. Due to the high magnetic forces and the low flow rates, the packed bed is stable; in particular, there are no magnetic beads dragged out. The actual packing height depends on the amount of loaded particles. Inhomogeneity within the particle bed was not considered. The reactions are described with multi-substrate reaction mechanisms [38]. The model approximated the three-dimensional rectangular microchannel by a one-dimensional model in Matlab. A more detailed description of the model can be found in the work of Burgahn et al. [36]. Reaction scheme of chiral reduction of NDK to the corresponding hydroxyketone (HK) and the diol using immobilized Gre2. The cofactor NADPH is regenerated during oxidation of glucose with glucose 1-dehydrogenase (GDH) [36].
In a first approach, we immobilized Gre2, which is genetically fused to a Halo-Tag variant, on Streptavidin-coated beads (Dynabeads M280), using a HaloTag-PEG-Biotin linker and kept GDH in the aqueous solution. The magnetic beads (MBs) were held inside the reactor with magnets placed below the reactor, forming a packed bed. Burgahn et al. clearly indicated a need for higher regeneration rates inside the packed bed [36]. To achieve a sufficient concentration of NADPH within the packed bed layer, GDH is immobilized on separate beads bearing reactive epoxy groups (Dynabeads M270 Epoxy) by employing a His-SpyCatcher linker. The tetrameric GDH is equipped with SpyTag fusion proteins and reacts with the linker in a covalent manner, thereby efficiently binding GDH to the beads.
The reactor encompasses a single microchannel with 200 µm height, 2 mm width and a particle bed with a length of 48 mm. The MBs are held within the reactor with eight magnets and form a packed bed at the bottom of the channel, schematically shown in Figure 2. This system was previously investigated by Burgahn et al., additionally applying a mathematical model to create a deeper understanding of the transport and reaction phenomena within the reactor. This model is also used here and is addressed as a basic Matlab model in the following. It can be regarded as a one-dimensional model, as it divides the reactor into sequential segments in the flow direction. Each segment consists of a packed bed and an ideal mixed free volume above the packed bed. Reactions and mass transport transversal to the flow direction in the packed bed are described by a set of ordinary differential equations, which are solved numerically assuming a finite mass transfer between the packed bed and the free volume above according to a Sherwood number correlation. Mass transport within the packed bed is calculated using effective diffusivities based on porosity and tortuosity of a compact packing of uniform, spherical particles. The particle bed is assumed to be a homogeneous package with constant height. Due to the high magnetic forces and the low flow rates, the packed bed is stable; in particular, there are no magnetic beads dragged out. The actual packing height depends on the amount of loaded particles. Inhomogeneity within the particle bed was not considered. The reactions are described with multi-substrate reaction mechanisms [38]. The model approximated the three-dimensional rectangular microchannel by a one-dimensional model in Matlab. A more detailed description of the model can be found in the work of Burgahn et al. [36].
The kinetic parameters were fitted for NDK concentrations ranging from 0.1 to 2 mmol/L for Gre2 and 1 to 100 mmol/L glucose for GDH reactions, covering the concentrations in this work. NADP/H ranged for both reactions between 0.1 and 1 mmol/L. For increased NADP/H concentration, an extrapolation was applied. The model was found to show good agreement with experimental data, which justifies its use for the optimization intended in this work. For bed heights of 88 µm, good agreement with experimental data could be shown. For flow rates from 0.5 to 2 µL/min, correct tendencies were predicted by the model; however, conversion was slightly overestimated for higher flow rates. This was attributed to remaining uncertainties of the mass transport parameters such as tortuosity and porosity. Reactor behavior for process parameters such as bed heights above 88 µm Symmetry 2021, 13, 524 4 of 20 (more than 4.5 mg MBs) and high NADPH concentrations are not in the validated range of parameters, hence are extrapolations. This needs to be done in future. Furthermore, the results hint at sufficiently and insufficiently used parts of the reactor, resulting in our desire to enhance the system. The kinetic parameters were fitted for NDK concentrations ranging from 0.1 to 2 mmol/L for Gre2 and 1 to 100 mmol/L glucose for GDH reactions, covering the concentrations in this work. NADP/H ranged for both reactions between 0.1 and 1 mmol/L. For increased NADP/H concentration, an extrapolation was applied. The model was found to show good agreement with experimental data, which justifies its use for the optimization intended in this work. For bed heights of 88 µm, good agreement with experimental data could be shown. For flow rates from 0.5 to 2 µL/min, correct tendencies were predicted by the model; however, conversion was slightly overestimated for higher flow rates. This was attributed to remaining uncertainties of the mass transport parameters such as tortuosity and porosity. Reactor behavior for process parameters such as bed heights above 88 µm (more than 4.5 mg MBs) and high NADPH concentrations are not in the validated range of parameters, hence are extrapolations. This needs to be done in future. Furthermore, the results hint at sufficiently and insufficiently used parts of the reactor, resulting in our desire to enhance the system.
To enable a sufficiently accurate model for optimization of the present reactor, models ranging from a one-dimensional to a three-dimensional representation of the microchannel were investigated and compared with each other regarding their applicability for the MLRD methodology. For this reason, the basic Matlab model was translated to ANSYS Fluent.

Systematic Design Methodology
The reactor design as described above holds potential for optimization [36]. While some reactor parameters are fixed due to the conditions, there are several others, such as feed concentration, MB enzyme loading and flow rate, which can be optimized. In the following sections, we give a short overview of the systematic of the multi-level reactor design (MLRD) approach [39] to enhance the reactor. We applied and adjusted this methodology to fit the present reactor system.
As mentioned before, there are a variety of different optimization methodologies. With the concept of elementary process functions (EPF), Freund and Sundmacher [29] took a very basic approach to the task of reactor optimization. They split chemical systems into functional modules regarding physical and chemical processes, which leads to an abstract module-specific system analysis. Based on the idea of EPF, Peschel [28] and To enable a sufficiently accurate model for optimization of the present reactor, models ranging from a one-dimensional to a three-dimensional representation of the microchannel were investigated and compared with each other regarding their applicability for the MLRD methodology. For this reason, the basic Matlab model was translated to ANSYS Fluent.

Systematic Design Methodology
The reactor design as described above holds potential for optimization [36]. While some reactor parameters are fixed due to the conditions, there are several others, such as feed concentration, MB enzyme loading and flow rate, which can be optimized. In the following sections, we give a short overview of the systematic of the multi-level reactor design (MLRD) approach [39] to enhance the reactor. We applied and adjusted this methodology to fit the present reactor system.
As mentioned before, there are a variety of different optimization methodologies. With the concept of elementary process functions (EPF), Freund and Sundmacher [29] took a very basic approach to the task of reactor optimization. They split chemical systems into functional modules regarding physical and chemical processes, which leads to an abstract module-specific system analysis. Based on the idea of EPF, Peschel [28] and Freund [39] developed a multi-level reactor design (MLRD), which has been used successfully for chemical reactors with catalyst deactivation [31], flexible feedstock [30] and catalytic particles [32].
The MLRD methodology contains three levels, shown in Figure 3-level 1: optimal reaction conditions, level 2: optimal reactor concept, and level 3: technical reactor. For level 1, generally, no limitations for heat or mass transport are considered, and thus the catalyst can work under optimal conditions at a given position. It is important to mention that these conditions are the utopian optimum for the reaction and reactor.
In level 2, transport processes are included, aiming to optimize driving forces. This includes the specific dimensions or shapes within a reactor. Thereby, optimal reaction volume designs are evaluated and eventually chosen.
Finally, in level 3 technical approximations are evaluated based on the previously defined optimal reactor conditions. The final set of reactor parameters will be a tradeoff between optimal reactor conditions, technical approximation and productivity of the investigated reactor system.
The MLRD methodology contains three levels, shown in Figure 3-level 1: optimal reaction conditions, level 2: optimal reactor concept, and level 3: technical reactor. For level 1, generally, no limitations for heat or mass transport are considered, and thus the catalyst can work under optimal conditions at a given position. It is important to mention that these conditions are the utopian optimum for the reaction and reactor.  [40] and further adaptation to suit a multi-enzyme cascade reactor with mass transport limitations with a fixed channel length.
In level 2, transport processes are included, aiming to optimize driving forces. This includes the specific dimensions or shapes within a reactor. Thereby, optimal reaction volume designs are evaluated and eventually chosen.
Finally, in level 3 technical approximations are evaluated based on the previously defined optimal reactor conditions. The final set of reactor parameters will be a tradeoff between optimal reactor conditions, technical approximation and productivity of the investigated reactor system. This method was already successfully applied for different chemical reactor systems, for example, synthesis of methanol [32], of ethylene oxide [28] and of biopharmaceutical manufacturing in Pichia pastoris [33]. The abstract approach of the EPF and MLRD methodology, based on fundamental processes, is a tool with great potential in the field of computer-aided design. It enables a systematic design of enhanced chemical reactors based on physical phenomena.
Considering our specific EPF, we adapted the levels of the MLRD methodology. Optimizing a process requires an understanding of the major contributing sources to enable the best possible reaction conditions. The present multi-enzyme system is realized with a heterogeneous catalyst bed, whereby the prevailing fluxes of materials must be balanced in the optimization. Regarding the reduction reaction, important fluxes include the cofactor regeneration and the mass transport. Within the first level, these fluxes were set to "unlimited" separately, resulting in three separate cases: This allows a clear distinction between the limitation by regeneration and mass transport. The second level will address all limitations due to fluxes and hence will pro- This method was already successfully applied for different chemical reactor systems, for example, synthesis of methanol [32], of ethylene oxide [28] and of biopharmaceutical manufacturing in Pichia pastoris [33]. The abstract approach of the EPF and MLRD methodology, based on fundamental processes, is a tool with great potential in the field of computer-aided design. It enables a systematic design of enhanced chemical reactors based on physical phenomena.
Considering our specific EPF, we adapted the levels of the MLRD methodology. Optimizing a process requires an understanding of the major contributing sources to enable the best possible reaction conditions. The present multi-enzyme system is realized with a heterogeneous catalyst bed, whereby the prevailing fluxes of materials must be balanced in the optimization. Regarding the reduction reaction, important fluxes include the cofactor regeneration and the mass transport. Within the first level, these fluxes were set to "unlimited" separately, resulting in three separate cases: • Case 1.1-unlimited cofactor regeneration and mass transport • Case 1.2-unlimited cofactor regeneration (mass transport limitation) • Case 1.3-unlimited mass transport (cofactor regeneration limitation) This allows a clear distinction between the limitation by regeneration and mass transport. The second level will address all limitations due to fluxes and hence will provide reasonably optimized reactor parameters. We formulated the condition for the optimization so as to utilize the reactor PMMA chip at its fullest, resulting in a technical solution for the investigated system. By this approach, we combine level 2 and level 3 into one case. Further technical measures to improve the reactor would be installed upstream or downstream of the reactor. These measures were not investigated within the scope of this work.

Optimization Methodology
Usually the performance of a catalytic reactor is characterized by yield, selectivity, conversion, space time yield (STY) and catalyst productivity. As the literature shows, these measures proved to be appropriate to optimize biochemical processes [27]. The values reached should be as high as possible. These scalars are defined by Equations (1)- (5). For our optimizations, the intermediate (HK) is the desired product, and therefore STY, yield and biocatalytic productivity Prod cat . STY is calculated with the concentration of the desired product (HK) c HK , its molecular mass M and the residence time τ. The yield relates the product concentration c HK and the feed concentration of the starting material NDK. The biocatalytic productivity is related to the mass of Gre2 m E , which is the enzyme catalyzing the reduction and the amount of HK m HK after 24 h.
The reactor is divided into N segments as outlined in Section 2.1 and in more detail in Burgahn et al. [36]. For each segment i, a mean enzyme utilization (EU) η j within the packed bed is defined. The segment's packed bed contains M cells of the height ∆y i with a respective reaction rate R based on the present concentration c i,j in the cell I and segment j. The sum of all the cells' reaction rates in relation to the volume averaged reaction rate R j and the bed thickness d bed results in the mean EU η j (Equation (4)). This provides an assessment of how well the enzymes in a given segment are utilized. The value of the mean EU ranges from 0 to 1, with 1 showing no transport and cofactor regeneration limitation. Values approaching 0 represent an increased transport and cofactor regeneration limitation. Combining all segments' mean EU, with respect to the length of the segment, results in a single scalar for the reactor's EU η R of length L (Equation (5)). η R is calculated for the reduction reaction of NDK to HK, since this is the desired reaction and ranges from 0 to 1.

Measure of Process Economics
Evaluating the economic performance of a biocatalytic process can be split into capital expenditure (CapEx) and operating expenditure (OpEx) [41]. CapEx covers the cost for installed equipment and materials in €. This includes in this particular case the magnetic beads (MBs) inside this reactor, since they are assumed to be reusable and stationary during the operation of the process. To enable immobilization, the enzymes have to be modified. This adds 90% of the pure enzyme costs, due to the procedure, equipment and labor, excluding the carrier [41]. MB costs are differentiated between streptavidin (Gre2) and epoxy (GDH) immobilization techniques, since streptavidin MBs are almost five times more expensive [42]. This leads to CapEx calculations as described in Equation (6), with the amount of carriers for Gre2 and GDH and the cost for the modified enzymes.
Operating costs (OpEx) consist of direct, indirect and fixed costs in € per time. The latter two can be calculated from capital investment as converted annual cost. Process parameters, conversion, etc. result in certain amounts of raw materials consumed in the process as given by the material balances. Raw material costs are listed in the Supplementary Materials. Additionally, energy consumption, maintenance, requirements for heating, sterilization and labor costs are included in the operation costs.
A complete economic analysis is a major effort, so we apply here a simplified cost estimation for the present reactor system. Assuming that all investigated cases are using the same equipment, we only investigated the case-specific CapEx, namely the MBs and the related immobilized enzymes. OpEx is calculated for pumped media and their components, excluding buffer and water. This results in a sum of all raw materials fed with the adapted flow rate of . V , molecular weight M i and their respective costs as shown in Equation (7). Deriving an overall scalar to measure the economic performance, we calculate product cost with CapEx and OpEx as given by Equation (8). This is achieved by relating costs and the amount of HK produced after one week of operation time (OT).
Due to the simplifications, the cost derived in this work only provides an economic trend for a cost-effective reactor design. For a detailed and realistic cost estimation, a complete economic study has to be performed, which is, however, beyond the scope of this paper.

Boundary Conditions
Due to the variety of adjustable parameters, we define boundary conditions with regard to the reactions and the reactor chip used in the experiments [35,36]. The desired yield for the intermediate should be achieved at the end of the reactor. The flow rate . V and hence the residence time τ will be adjusted to achieve the desired yield. To keep the case studies consistent and comparable, the feed concentrations of NDK and glucose are kept constant. For all cases, fixed values of temperature of 30 • C and pH value of 7.6 are assumed. These values had been identified as an appropriate choice in the experimental screening phase. The pH value is constant within the reactor due to the addition of a T-TEMg buffer to the feed. A steady temperature was achieved with a thermostat for heating and cooling the PMMA chip. Of course, temperature and pH value have a great impact on the reaction kinetics, and their influences have to be built into the kinetic model if variations of temperature and pH value are to be considered. In our case we excluded both as variable process parameters. Considering this in the optimization would require a kinetic model describing their influence on the reaction rate, which has not yet been established. Furthermore, all results depict stationary reactor conditions.
Enzyme immobilization provides a general system parameter, whether GDH is fixed in the reactor or dissolved in the feed. The latter benefits from cofactor regeneration in the liquid volume with an equilibrated NADP/NADPH concentration in the feed. However, in that case expensive GDH flows through the reactor, which increases the product cost. Although the immobilization demands further costs in preparation, we obtain an enzyme fixed within the reactor. Compared with feeding GDH, the average enzyme concentration within the layer is higher due to the loading capacity of the beads.
Kinetic parameters from Burgahn et al. are used in all calculations in this work [36]. As reported in the literature, the immobilization of GDH has a negligible impact on its activity [35]. Therefore, kinetic parameters for immobilized GDH are assumed to be equal to those of free GDH.
The packed bed within the reactor is assumed to be homogeneous with a fixed height. The bed height is proportional to the amount of beads loaded in the channel. The basic case contains 4.5 mg of beads in the reactor, resulting in a bed height of 88 µm [36]. Due to different methods of immobilization used for Gre2 and GDH, they cannot be provided co-immobilized on the same magnetic beads. The total amount of beads held in the reactor is comprised of both types of loaded MBs (m total = m Gre2 + m GDH ). Hence, adding immobilized GDH will reduce the amount of Gre2 in the reactor for the same amount of magnetic beads or bed height. For both immobilization methods, a maximum enzyme loading is assumed (Gre2: 24 pmol/mg; GDH: 55 pmol/mg).
These conditions result in a set of parameters for the investigated reactor system. As shown in Figure 4, the parameters are classified into reactor and flow parameters. The latter describe the concentrations in the feed stream. NADP in particular is an interesting parameter, since NADPH is three times more expensive then NADP [43,44]. For the different cases of an immobilized or dissolved regeneration enzyme GDH, its concentration in the feed, the enzyme loading on the beads and the bead ratio are included in the reaction and reactor parameters. Reactor parameters are bound to fixed costs for the investigated case. Bed thickness is proportional to the amount of loaded beads. For the case of immobilized GDH, the ratio of Gre2 versus GDH beads is an essential reactor parameter.
Flow rate . V and residence time τ are adapted due to the condition of maximal HK yield at the end of the reactor. bilized GDH will reduce the amount of Gre2 in the reactor for the same amount of magnetic beads or bed height. For both immobilization methods, a maximum enzyme loading is assumed (Gre2: 24 pmol/mg; GDH: 55 pmol/mg).
These conditions result in a set of parameters for the investigated reactor system. As shown in Figure 4, the parameters are classified into reactor and flow parameters. The latter describe the concentrations in the feed stream. NADP in particular is an interesting parameter, since NADPH is three times more expensive then NADP [43,44]. For the different cases of an immobilized or dissolved regeneration enzyme GDH, its concentration in the feed, the enzyme loading on the beads and the bead ratio are included in the reaction and reactor parameters. Reactor parameters are bound to fixed costs for the investigated case. Bed thickness is proportional to the amount of loaded beads. For the case of immobilized GDH, the ratio of Gre2 versus GDH beads is an essential reactor parameter. Flow rate and residence time are adapted due to the condition of maximal HK yield at the end of the reactor.

Techno-Economic Performance (TEP)
With the measures described in Section 2.3.1, during the optimization, a tradeoff between a high yield and economic viability with regard to the defined boundary conditions has to be met. To estimate STY, biocatalytic productivity and enzyme utilization (EU) are calculated for maximum yield , obtained with the respective value of ′.
Summarizing the measures to fully evaluate the reactor with a focus on performance and economics, the techno-economic performance (TEP) in Equation (9) was introduced and is used as an objective function to be maximized for our system. The case from the investigations by Burgahn

Techno-Economic Performance (TEP)
With the measures described in Section 2.3.1, during the optimization, a tradeoff between a high yield and economic viability with regard to the defined boundary conditions has to be met. To estimate STY, biocatalytic productivity and enzyme utilization (EU) are calculated for maximum yield Y , obtained with the respective value of τ .
Summarizing the measures to fully evaluate the reactor with a focus on performance and economics, the techno-economic performance (TEP) in Equation (9) was introduced and is used as an objective function to be maximized for our system. The case from the investigations by Burgahn et al. was set as a reference for the TEP [36]. All measures refer to the adapted space-time τ for achieving a maximum intermediate yield Y . Since the productivity is related to the amount of enzyme within the system, the EU is covered in the productivity.

Reactor Model
With an increase in dimensions, mathematical models gain accuracy in modeling fluid dynamics at a tradeoff to the computational effort. For the investigation and especially the optimization of a reactor setup, an accurate and sufficiently fast model is needed, since a lot of calculations have to be performed. Translating the basic 1D Matlab model briefly explained in Section 2.1 into two and three-dimensional computational fluid dynamics (CFD) models in ANSYS Fluent (in the following named model F2D and F3D) enables more detailed results on the concentration patterns in the liquid volume above the packed bed. In addition to a more realistic description of the flow behavior of the fluid, the CFD-based models also account for diffusional transport in all directions, which allows for a more accurate quantification of the mass transport from the bulk liquid flow along the channel to the packed bed.
The kinetics of the enzymatic reactions were included via User-Defined-Functions. Transient calculations were performed until steady values were obtained. To keep calcu-lation time reasonable only the entrance region of the reactor was simulated, where the highest reaction rates and therefore also concentration gradients are present.
As a result of the channel geometry and the assumed constant bed height throughout the reactor, the investigated reactor exhibits axial symmetry. Thus, a two-dimensional representation was applied first. The principle of reducing the dimension from a realistic three-dimensional model to a two-dimensional model is commonly used [20,45].
As the Matlab model basically chains axial segments of the free channel volume sequentially, the calculations result in a two-dimensional depiction of the reactor, which can be compared with the model F2D. Figure 5 shows the concentration distribution as calculated by the Matlab and F2D models for the first 2 mm of the reactor. The Matlab model assumes a constant flow velocity over the whole cross-section of the free channel volume, whereas the F2D model yields a parabolic flow velocity profile typical of laminar channel flow. However, the concentration distributions in the packed bed show only small differences. This is even more obvious from Figure 6, which compares the concentration profiles in the y-direction at the beginning of the reactor. The only notable difference is that the Matlab model slightly overestimates the mass transport resistance from the bulk liquid flow to the surface of the packed bed. This is obviously due to the correlation used to derive the Sherwood number. Thereby, a more pronounced drop in the concentration is predicted, and hence small deviations at the top of the catalytic bed are present. However, the concentration profiles within the catalytic bed are almost identical.  Concentration distribution for feed NDK, product HK and diol for the first 2 mm of the reactor. The lower part (brown) depicts the particle bed, and the volume above is the free volume with convective flow. An additional calculation for 25 mm reactor length was performed to further investigate differences between the models. The results are shown in Figure 7, which reports the average concentration in the volume above the bed. As the curves for F2D and the basic Matlab model are very similar, the most pronounced differences occur in the entrance region of the reactor. The effects are either diminished, maintained or slowly accu-    Concentration distribution for feed NDK, product HK and diol for the first 2 mm of the reactor. The lower part (brown) depicts the particle bed, and the volume above is the free volume with convective flow. An additional calculation for 25 mm reactor length was performed to further investigate differences between the models. The results are shown in Figure 7, which reports the average concentration in the volume above the bed. As the curves for F2D and the basic Matlab model are very similar, the most pronounced differences occur in the en-   An additional calculation for 25 mm reactor length was performed to further investigate differences between the models. The results are shown in Figure 7, which reports the average concentration in the volume above the bed. As the curves for F2D and the basic Matlab model are very similar, the most pronounced differences occur in the entrance region of the reactor. The effects are either diminished, maintained or slowly accumulated along the channel. Since the reduction rate of NDK gets smaller with increasing x-position in the reactor, less NADPH is needed and hence the cofactor regeneration system reaches its equilibrium. Thereby the differences between the Matlab and Fluent models reduce with increasing x-position. Altogether, however, the deviations are small. At 25 mm (half) of the reactor, a maximum deviation of 3.4% of the HK concentration between the F2D model and the basic Matlab model is present for the case shown in Figure 7. Aiming to select a suitable model to be used for the intended optimization, the effort in computing time for the whole reactor has to be considered. A comparison is given in Table 1 for the three different models. With respect to hundreds of reactor simulations to be performed in the course of the optimization, the decision for the basic Matlab model is obvious. It cuts the computational time immensely and provides the best compromise between accuracy and computational time for our investigated systems.

GDH Immobilization
The immobilization of GDH influences the economic efficiency and the reactor behavior. To compare the effects of immobilization, conditions leading to the same STY are investigated. For both cases 4.5 mg of total MBs was investigated, since this is the easiest realizable amount of beads in the experiments. With a lower amount of beads, the EU and hence the catalytic productivity for aqueous GDH will increase as shown by Burgahn et al. [36]. For all cases, the feed consists of the raw materials NDK, glucose and NADP. In the case of aqueous GDH, NADPH is already produced before the feed is applied to the We further investigated the dimensional influence applied to this reactor. The F3D model expands the advantages of the F2D model in comparison to the basic Matlab model, with a consideration of influences by the corners of the channel and possible diffusioninduced inhomogeneity across the width of the channel. Assuming no influences by the solute species, the flow regime is found to be fully laminar (Re = 0.5). Due to the increase in complexity, the simulated length of the reactor had to be reduced to 1 mm to keep manageable computing times. Differences between models F2D and F3D are negligible with a deviation of 0.1% as shown in Figure S2.
Aiming to select a suitable model to be used for the intended optimization, the effort in computing time for the whole reactor has to be considered. A comparison is given in Table 1 for the three different models. With respect to hundreds of reactor simulations to be performed in the course of the optimization, the decision for the basic Matlab model is obvious. It cuts the computational time immensely and provides the best compromise between accuracy and computational time for our investigated systems.

GDH Immobilization
The immobilization of GDH influences the economic efficiency and the reactor behavior. To compare the effects of immobilization, conditions leading to the same STY are investigated. For both cases 4.5 mg of total MBs was investigated, since this is the easiest realizable amount of beads in the experiments. With a lower amount of beads, the EU and hence the catalytic productivity for aqueous GDH will increase as shown by Burgahn et al. [36]. For all cases, the feed consists of the raw materials NDK, glucose and NADP. In the case of aqueous GDH, NADPH is already produced before the feed is applied to the reactor, hence enabling an immediate onset of reduction reactions, compared with generating NADPH only within the bed in the case of immobilized GDH. NADPH is regenerated in situ in both cases. However, the mean GDH concentration in the bed is more than 15-fold higher when immobilized, compared with aqueous GDH. Aiming to select a suitable model to be used for the intended optimization, the effort in computing time for the whole reactor has to be considered. A comparison is given in Table 1 for the three different models. With respect to hundreds of reactor simulations to be performed in the course of the optimization, the decision for the basic Matlab model is obvious. It cuts the computational time immensely and provides the best compromise between accuracy and computational time for our investigated systems.

GDH Immobilization
The immobilization of GDH influences the economic efficiency and the reactor behavior. To compare the effects of immobilization, conditions leading to the same STY are investigated. For both cases 4.5 mg of total MBs was investigated, since this is the easiest realizable amount of beads in the experiments. With a lower amount of beads, the EU and hence the catalytic productivity for aqueous GDH will increase as shown by Burgahn et al. [36]. For all cases, the feed consists of the raw materials NDK, glucose and NADP. In the case of aqueous GDH, NADPH is already produced before the feed is applied to the reactor, hence enabling an immediate onset of reduction reactions, compared with generating NADPH only within the bed in the case of immobilized GDH. NADPH is regenerated in situ in both cases. However, the mean GDH concentration in the bed is more than 15-fold higher when immobilized, compared with aqueous GDH.
Since applying immobilized GDH reduces the amount of Gre2 MBs in the reactor when the total amount of MBs is fixed, an optimal relation had to be found. For the same amount of total beads in the reactor, the same STY was iteratively derived for a bead ratio Aiming to select a suitable model to be used for the intended optimization, the effort in computing time for the whole reactor has to be considered. A comparison is given in Table 1 for the three different models. With respect to hundreds of reactor simulations to be performed in the course of the optimization, the decision for the basic Matlab model is obvious. It cuts the computational time immensely and provides the best compromise between accuracy and computational time for our investigated systems.

GDH Immobilization
The immobilization of GDH influences the economic efficiency and the reactor behavior. To compare the effects of immobilization, conditions leading to the same STY are investigated. For both cases 4.5 mg of total MBs was investigated, since this is the easiest realizable amount of beads in the experiments. With a lower amount of beads, the EU and hence the catalytic productivity for aqueous GDH will increase as shown by Burgahn et al. [36]. For all cases, the feed consists of the raw materials NDK, glucose and NADP. In the case of aqueous GDH, NADPH is already produced before the feed is applied to the reactor, hence enabling an immediate onset of reduction reactions, compared with generating NADPH only within the bed in the case of immobilized GDH. NADPH is regenerated in situ in both cases. However, the mean GDH concentration in the bed is more than 15-fold higher when immobilized, compared with aqueous GDH.
Since applying immobilized GDH reduces the amount of Gre2 MBs in the reactor when the total amount of MBs is fixed, an optimal relation had to be found. For the same amount of total beads in the reactor, the same STY was iteratively derived for a bead ratio Aiming to select a suitable model to be used for the intended optimization, the effort in computing time for the whole reactor has to be considered. A comparison is given in Table 1 for the three different models. With respect to hundreds of reactor simulations to be performed in the course of the optimization, the decision for the basic Matlab model is obvious. It cuts the computational time immensely and provides the best compromise between accuracy and computational time for our investigated systems.

GDH Immobilization
The immobilization of GDH influences the economic efficiency and the reactor behavior. To compare the effects of immobilization, conditions leading to the same STY are investigated. For both cases 4.5 mg of total MBs was investigated, since this is the easiest realizable amount of beads in the experiments. With a lower amount of beads, the EU and hence the catalytic productivity for aqueous GDH will increase as shown by Burgahn et al. [36]. For all cases, the feed consists of the raw materials NDK, glucose and NADP. In the case of aqueous GDH, NADPH is already produced before the feed is applied to the reactor, hence enabling an immediate onset of reduction reactions, compared with generating NADPH only within the bed in the case of immobilized GDH. NADPH is regenerated in situ in both cases. However, the mean GDH concentration in the bed is more than 15-fold higher when immobilized, compared with aqueous GDH.
Since applying immobilized GDH reduces the amount of Gre2 MBs in the reactor when the total amount of MBs is fixed, an optimal relation had to be found. For the same amount of total beads in the reactor, the same STY was iteratively derived for a bead ratio  Since applying immobilized GDH reduces the amount of Gre2 MBs in the reactor when the total amount of MBs is fixed, an optimal relation had to be found. For the same amount of total beads in the reactor, the same STY was iteratively derived for a bead ratio of 0.6 m Gre2 m total . The resulting reactor characteristics are shown in Table 2. Immobilization of GDH results in benefits in all aspects. The same STY means an increase in productivity from 155.9 to 263.5 gHK gGre2 after 24 h operating time for immobilized GDH, since less Gre2 is applied. A batch reactor with the same reaction system with aqueous GDH was experimentally investigated previously [36]. It achieved a maximum HK yield after 2 h, resulting in a productivity of 63.2 gHK gGre2 and an STY of 11.6 g Lday . To uphold this productivity and STY would demand an ideal situation with 100% catalyst recyclability and zero downtime between two consecutive batches. The productivity of the flow system increases with OT as the amount of enzyme is constant. To achieve the idealized batch productivity, the setup with aqueous and immobilized GDH would need to run 5 h for aqueous and 3 h for immobilized GDH. For immobilized GDH, the bed has a higher η R , which results from an increased GDH concentration in the bed as compared to aqueous GDH. Immobilization is reported to reduce the production costs due to the fixation within the reactor [10,41]. While the operating cost decreases by more than 50%, the capital cost is reduced as well, since GDH beads are cheaper than Gre2 beads.
As indicated in Figure 8, immobilization enables higher NADPH concentration in the bed, due to ongoing regeneration. This is a consequence of high enzyme loadings on the magnetic beads. The reactor behavior with aqueous GDH is set as the reference case and has the TEP of 1. The TEP for immobilized GDH has been found to be 2.28, which shows the overall benefits of GDH immobilization. Due to this we studied only the reactor system with immobilized GDH. Figure 8, immobilization enables higher NADPH concentration in the bed, due to ongoing regeneration. This is a consequence of high enzyme loadings on the magnetic beads. The reactor behavior with aqueous GDH is set as the reference case and has the TEP of 1. The TEP for immobilized GDH has been found to be 2.28, which shows the overall benefits of GDH immobilization. Due to this we studied only the reactor system with immobilized GDH.

Level 1: Infinite Fluxes
Within the reactor, a tradeoff between optimal reaction conditions and limiting physical phenomena has to be met. To investigate the influences of cofactor regeneration and mass transport fluxes, they were artificially set to unlimited in three different cases. As described in Section 2.2, this level shows the optimal reaction conditions with a view to bead ratio and NADP concentration in the feed. The NADP feed concentration sets the maximal possible equilibrium NADPH concentration within the system. In regard to this, we address this correlation with "NADP/H" concentration. In cases 1.1 and 1.2 infinitely fast cofactor regeneration enables a constant, maximum equilibrium NADPH concentration and hence will minimize the cofactor influences on the reduction rates. Raising the NADP/H concentration increases conversion and yield while increasing the operating costs. For case 1.2, the optimum NADP feed concentration is evaluated. Enabling a limited cofactor regeneration (case 1.3) aims at optimizing the bead ratios. This is crucial for the reactor behavior, since this parameter balances the enzyme ratios within the packed bed. Besides the importance of the enzyme ratio, its optimization is rarely applied for enzymatic systems [27]. A genetic algorithm [46] was used to derive optimum reaction conditions, resulting in a feed concentration of 8.4 mmol/L NADP and a bead ratio of 0.743 m Gre m total , respectively. Details on the results of applying the generic algorithm can be found in the Supplementary Materials.
The respective results include the EU η R , the STY compared with the productivity, the production cost, and the TEP, which are all shown in Figure 9. Distinct behaviors resulting from the investigated physical phenomena are observable. As case 1.1 is not limited, the EU is 1, and the lowest cost and highest STY and productivity are achieved. Notable is a constant productivity over all amounts of beads due to the same reaction conditions for the enzymes. This is reflected in constant production cost ( Figure 9C) and consequently an almost linearly rising TEP ( Figure 9D). If we compare cases 1.1 and 1.2 in their STY-productivity trends ( Figure 9B), a connection to the case of unlimited mass transport is evident, underlined by the value of 1 for the EU. Reduced productivity and STY are both due to changing equilibria of the regeneration reaction, caused by a rising concentration of gluconolacton within the reactor. This lowers the concentration of NADPH and hence reduces the reduction rate. To achieve the same yield, flow rates have to be reduced, resulting in higher residence time and lower STYs and productivities and thus higher costs and lower TEPs.

Levels 2 and 3: (Technical) Reactor Concept
A genetic algorithm [46] was employed to derive optimal reactor parameters with regard to performance and economics. Hence, the TEP described in Section 2.3.3 was maximized with the given boundary conditions from Section 2.3.2. The bead ratio is bound from 0 to 1, ranging from solely GDH beads to solely Gre2 beads. The NADP concentration was capped at 20 mmol/L in the feed stream. The results for immobilized GDH enzyme from Section 3.1 were used as the starting point. A TEP value of 1 is related to the result for aqueous GDH (base case).
The genetic algorithm derived optimal parameter sets for the different MB amounts as shown in Figure 10. The algorithm provides a lot of data points, which show a considerable variation. It can easily be seen that for each bead loading a limiting curve is established. These curves obtain successively higher maxima up to about 4.5 mg of MBs. For even higher amounts of beads, no positive effect is achieved anymore. It is evident that both functional relationships exhibit a maximum and neither curve is symmetrical.
Regarding the concentration of NADP plus NADPH ( Figure 10A), it is observed that up to a value of about 10 the TEP rises fast, at first almost linear, then a maximum is achieved. This is caused by the availability of hydrogen for the desired enzyme reaction. Nevertheless, this effect reaches an economic limit, once the supply of hydrogen is not a crucial parameter anymore. We can conclude from these results that an amount of 4.5 mg of MBs is optimal. At that point, the concentration for the hydrogen carrier should be 10 mmol/L.
The curve for the bead ratio ( Figure 10B) is a bit more complex. On both extremes a very steep rise is observed, as expected, since both enzymes are needed to obtain high conversion. The optimal bead ratio is found at 0.74 . The optimum bead ratio and NADP feed concentration differ from the optimum values in Section 3.3. The reactor behavior is a combination of the limitations by mass transport and cofactor regeneration and cannot be represented by either one of those. Thereby new optima for bead ratio and NADP feed concentration appear. By enabling mass transport limitation in case 1.3, an additional consequence is a reduced EU for higher MB loads. Due to an increased bead loading, the bed height and thus the diffusion path length increase. For small bead quantities, case 1.3 approaches case 1.1. The EU explains lower STY and productivity in Figure 9B, as different reaction conditions are present within the bed. Thereby Gre2 enzymes deeper within the bed are less productive. For smaller amounts of beads, this effect of productivity inhomogeneity is reduced, approaching the outcomes of case 1.1. The TEP for case 1.3 shows an approaching plateau, which foreshadows a possible saturation for the investigated system. It also shows that cases 1.1 and 1.2 are probably too strongly simplified.
In Figure 9B-D, for a bead loading of 2 mg an interception between cases 1.2 and 1.3 is visible. At this loading, the limitations induced by mass transport and cofactor regeneration are equal. For higher bead loadings, mass transport limitations are of higher importance than cofactor regeneration limitations.

Levels 2 and 3: (Technical) Reactor Concept
A genetic algorithm [46] was employed to derive optimal reactor parameters with regard to performance and economics. Hence, the TEP described in Section 2.3.3 was maximized with the given boundary conditions from Section 2.3.2. The bead ratio is bound from 0 to 1, ranging from solely GDH beads to solely Gre2 beads. The NADP concentration was capped at 20 mmol/L in the feed stream. The results for immobilized GDH enzyme from Section 3.1 were used as the starting point. A TEP value of 1 is related to the result for aqueous GDH (base case).
The genetic algorithm derived optimal parameter sets for the different MB amounts as shown in Figure 10. The algorithm provides a lot of data points, which show a considerable variation. It can easily be seen that for each bead loading a limiting curve is established. These curves obtain successively higher maxima up to about 4.5 mg of MBs. For even higher amounts of beads, no positive effect is achieved anymore. It is evident that both functional relationships exhibit a maximum and neither curve is symmetrical.
Regarding the concentration of NADP plus NADPH ( Figure 10A), it is observed that up to a value of about 10 the TEP rises fast, at first almost linear, then a maximum is achieved. This is caused by the availability of hydrogen for the desired enzyme reaction. Nevertheless, this effect reaches an economic limit, once the supply of hydrogen is not a crucial parameter anymore. We can conclude from these results that an amount of 4.5 mg of MBs is optimal. At that point, the concentration for the hydrogen carrier should be 10 mmol/L. As a result, cost evaluations allow for higher product cost, due to lower flow rates for maximum HK yield. As the parameters are not influenced by the amount of beads, the results establish a simple relation between MB load and optimal conditions. Higher amounts of MBs are advantageous at these parameters, as all enzymes are used sufficiently to compensate for an increase in costs. As an increased NADP feed concentration increases the regeneration rate within the packed bed, their limitation is lowered. As implied in Figure 10, an increased TEP is shown for higher quantities of beads. Figure 11 shows all cases calculated by the genetic algorithm in the relationship between STY and productivity. For improvement in visualization, grey dashes show the limits for different amounts of MBs. The symbols are color-coded according to the TEP from red (low) via yellow to green (high). The initial case was 263 g/g productivity and 96.3 g/L day for STY and is marked with a star. An optimum in relation to TEP is indicated by the black dashed line. The curve for the bead ratio ( Figure 10B) is a bit more complex. On both extremes a very steep rise is observed, as expected, since both enzymes are needed to obtain high conversion. The optimal bead ratio is found at 0.74 m Gre2 m total . The optimum bead ratio and NADP feed concentration differ from the optimum values in Section 3.3. The reactor behavior is a combination of the limitations by mass transport and cofactor regeneration and cannot be represented by either one of those. Thereby new optima for bead ratio and NADP feed concentration appear.
As a result, cost evaluations allow for higher product cost, due to lower flow rates for maximum HK yield. As the parameters are not influenced by the amount of beads, the results establish a simple relation between MB load and optimal conditions. Higher amounts of MBs are advantageous at these parameters, as all enzymes are used sufficiently to compensate for an increase in costs. As an increased NADP feed concentration increases the regeneration rate within the packed bed, their limitation is lowered. As implied in Figure 10, an increased TEP is shown for higher quantities of beads. Figure 11 shows all cases calculated by the genetic algorithm in the relationship between STY and productivity. For improvement in visualization, grey dashes show the limits for different amounts of MBs. The symbols are color-coded according to the TEP from red (low) via yellow to green (high). The initial case was 263 g/g productivity and 96.3 g/L day for STY and is marked with a star. An optimum in relation to TEP is indicated by the black dashed line.
The shape of the limiting grey dashes for different amounts of beads overlaps for very high productivities. With increasing bead loadings, more enzymes are in the reactor enabling faster NDK reduction rates. As the overlapping dashes indicate, cofactor regeneration limitations are very small, as the curvature is similar to the results of case 1.3 in Figure 9. The same can be concluded for the optimum line in Figure 11. The colors indicate an improved TEP for higher amounts of beads. Although specific conditions allow for higher possible productivity, economic considerations shift the optimal operating points to slightly lower STY and productivity.
Symmetry 2021, 13, 524 16 of 21 Figure 11. Optimization of a microchannel enzymatic reactor system with an overflown particle bed and immobilized enzyme cascades. Results generated by a genetic algorithm with characterizing Pareto for different amounts of beads, indicated by grey dashed lines. An increase in TEP is depicted by symbols' color, ranging from red (low TEP) to green (high TEP).
The shape of the limiting grey dashes for different amounts of beads overlaps for very high productivities. With increasing bead loadings, more enzymes are in the reactor enabling faster NDK reduction rates. As the overlapping dashes indicate, cofactor regeneration limitations are very small, as the curvature is similar to the results of case 1.3 in Figure 9. The same can be concluded for the optimum line in Figure 11. The colors indicate an improved TEP for higher amounts of beads. Although specific conditions allow for higher possible productivity, economic considerations shift the optimal operating points to slightly lower STY and productivity.
Since NADPH replenishment is secured, the selectivity to HK for higher bead loadings slightly decreases as shown in Figure 12A. This includes a higher amount of wasted NDK. For thicker particle beds, diffusion path length increases, increasing the time for all species in the bed, and hence giving more time for the consecutive reaction decreasing HK. The EU decreases with an increased amount of beads indicating mass transport limitation, as can be seen from Figure 12B.
OpEx results are similar for all bead loadings, due to the same optimal NADP/H concentration and HK yield. CapEx in turn is determined by the amount of MBs and hence increases ( Figure 12C). As explained above, all CapEx values are calculated for one week of operation, thus the effect of increasing CapEx on specific product cost is very high. After a realistic estimation of catalyst deactivation, which would give us a reasonable time on stream estimation for the reaction system, we could redo the calculation using reasonable assumptions. This could lead to different results, depending on the achievable time on stream.
The correlation of TEP with MBs indicates a maximum of 13.25 at a value of about 5 mg of beads ( Figure 12D). This shows a balance between higher product cost and decreased selectivity for an increased quantity of beads on the one hand and the profit of higher flow rates on the other hand. Differences in TEP between 4 and 6 mg MBs are only small. Increasing the operation time would increase the amount of product and thereby decrease the contribution of CapEx to the production cost. It should be kept in mind that this maximum is valid for one week of operation. A different operation time impacts the CapEx influence within the cost calculations. For shorter operation times, CapEx makes a stronger contribution to the cost, and hence lower bead amounts are beneficial, as the m↑ TEP Figure 11. Optimization of a microchannel enzymatic reactor system with an overflown particle bed and immobilized enzyme cascades. Results generated by a genetic algorithm with characterizing Pareto for different amounts of beads, indicated by grey dashed lines. An increase in TEP is depicted by symbols' color, ranging from red (low TEP) to green (high TEP).
Since NADPH replenishment is secured, the selectivity to HK for higher bead loadings slightly decreases as shown in Figure 12A. This includes a higher amount of wasted NDK. For thicker particle beds, diffusion path length increases, increasing the time for all species in the bed, and hence giving more time for the consecutive reaction decreasing HK. The EU η R decreases with an increased amount of beads indicating mass transport limitation, as can be seen from Figure 12B.  Figure 13A,B shift. TEP refers to the corresponding TEP of the basic case and operation time. With increased time of operation, higher bead amounts become viable as shown in Figure 13D. Additionally, we can note that increasing the operation times reduces the maximum achievable adapted TEP, but it is still a major improvement compared with the base case. An amount of beads of 4.5 mg shows great TEPs for all operation times with minor drawbacks for higher operation times. Although these evaluations are promising and show great enhancement in performance and economics, experiments should be performed to validate the calculations and trends. OpEx results are similar for all bead loadings, due to the same optimal NADP/H concentration and HK yield. CapEx in turn is determined by the amount of MBs and hence increases ( Figure 12C). As explained above, all CapEx values are calculated for one week of operation, thus the effect of increasing CapEx on specific product cost is very high. After a realistic estimation of catalyst deactivation, which would give us a reasonable time on stream estimation for the reaction system, we could redo the calculation using reasonable assumptions. This could lead to different results, depending on the achievable time on stream.
The correlation of TEP with MBs indicates a maximum of 13.25 at a value of about 5 mg of beads ( Figure 12D). This shows a balance between higher product cost and decreased selectivity for an increased quantity of beads on the one hand and the profit of higher flow rates on the other hand. Differences in TEP between 4 and 6 mg MBs are only small. Increasing the operation time would increase the amount of product and thereby decrease the contribution of CapEx to the production cost. It should be kept in mind that this maximum is valid for one week of operation. A different operation time impacts the CapEx influence within the cost calculations. For shorter operation times, CapEx makes a stronger contribution to the cost, and hence lower bead amounts are beneficial, as the maxima in Figure 13A,B shift. TEP refers to the corresponding TEP of the basic case and operation time. With increased time of operation, higher bead amounts become viable as shown in Figure 13D. Additionally, we can note that increasing the operation times reduces the maximum achievable adapted TEP, but it is still a major improvement compared with the base case. An amount of beads of 4.5 mg shows great TEPs for all operation times with minor drawbacks for higher operation times. Although these evaluations are promising and show great enhancement in performance and economics, experiments should be performed to validate the calculations and trends. maxima in Figure 13A,B shift. TEP refers to the corresponding TEP of the basic case and operation time. With increased time of operation, higher bead amounts become viable as shown in Figure 13D. Additionally, we can note that increasing the operation times reduces the maximum achievable adapted TEP, but it is still a major improvement compared with the base case. An amount of beads of 4.5 mg shows great TEPs for all operation times with minor drawbacks for higher operation times. Although these evaluations are promising and show great enhancement in performance and economics, experiments should be performed to validate the calculations and trends.   As mentioned previously, the boundary conditions make the results from level 2 a viable technical solution. Additional measures that are not related to the reactor can influence the TEP. Since OpEx is a major contributing factor to production cost, measures to reduce this would enhance the performance with regard to TEP. NADP/H is the most expensive content in the feed. By recycling NADP/H OpEx will decrease, while CapEx will be raised due to the separation method. In the long term, this would allow higher NADP/H con-centration, boosting reduction rates significantly. The potential of an NADP/H recycling system was reported by Baumer et al. with a two-phase system and a flow liquid-liquid extraction (FLLEX) [47]. Another option for single-phase recycling is micro simulated moving bed (µSMB) systems [48]. These recycling methods can be implemented in the model of this work by adding a recycling factor for the OpEx for NADP and NADPH, alongside additional CapEx for the separation unit and the separation efficiency. However, as we could not retrieve quantitative information from the literature, we did not include this topic in the present study.

Conclusions
Immobilized enzyme cascades in a flow system show great benefits over traditional batch systems. The application of complex reaction systems exhibits different demands on suitable flow systems. Combined with the possible limitations due to cofactor regeneration and mass transport, optimization of such systems is challenging. Optimization in combination with mathematical modeling allows us to predict and derive enhanced reactor systems, while reducing the experimental effort. We showed that a simple Matlab model, which describes a diffusion-driven processes in combination with a complex reaction network, is feasible, especially for optimization purposes. Mathematical modeling of such reactor systems enables insights into the processes controlling the reactor performance, indicating insufficiently used regions, major contributing limitations and bottlenecks. Maintaining low production cost while improving the performance is an important aspect in scale-up and competitiveness of continuous flow reactors. Within this work, we utilize a modified multi-level reactor design (MLRD) approach to optimize selective reduction with cofactor regeneration catalyzed by immobilized enzymes. Optimization aims to maximize the techno-economic performance (TEP).
In a first step, we showed that immobilization of the regeneration enzyme GDH alongside the ketoreductase Gre2 is beneficial in all aspects of performance, economics and lastly TEP. To achieve the same space time yield (STY), catalytic productivity increases while the production cost decreases, as depicted in Figure 14. NADP feed concentration and enzyme ratios are crucial parameters of the system. High NADPH demand of the reduction reactions can be compensated by an increased NADPH concentration in the feed and additional immobilized GDH in the packed bed. Optimal reactor conditions were derived for those two parameters and for different amounts of beads with the aim of maximizing hydroxyketone (HK) yield. For all cases investigated in level 2, an NADP feed concentration of 10 mmol/L and a bead ratio of 0.74 m Gre2 m total gave the best TEP independent of the amount of beads inside the reactor. An increase in bead loading results in higher STY while reducing catalyst productivity, as the EU is decreasing. While production cost for HK rises with the utilization of more beads, the overall TEP shows the highest improvement with 5 mg of loaded magnetic beads for an operation time of 1 week. Figure 14D summarizes the evolution of TEP in this work from the base case with aqueous GDH to the derived optimum. While immobilization improves the TEP by a factor of 2.3, an optimal NADP feed concentration and bead ratio increases the TEP further to 13.25 after evaluating more than 1300 reactor configurations. We predict an increase of the STY from 100 to 260 gHK L day and of the specific productivity from 150 to 430 gHK gGre2 , while reducing the production cost by 50% for an operation time of 1 week. Further improvements in the economics by including NADP/H recycling would further increase the TEP. It has been shown conclusively how the MLRD approach led to a deeper understanding of the process and enabled its optimization. Such continuous-flow microreactor systems could replace conventional batch reactors. We look at the system as a possible new platform technology for continuous flow biochemical production. Upscaling of the production would be possible by internal and external numbering-up. Internal numbering-up brings about the challenge of uniform distribution of the flow to a multitude of channels while avoiding the need for duplication of peripheral components (see, for example, [49]). For single-phase flow, the internal flow distribution can be solved with the aid of computational fluid dynamics simulations and precise fabrication methods. Well-engineered continuous-flow systems can deliver productivities and yields comparable to or better than conventional batch systems over a wide range of production volumes relevant for specialty products.
place conventional batch reactors. We look at the system as a possible new platform technology for continuous flow biochemical production. Upscaling of the production would be possible by internal and external numbering-up. Internal numbering-up brings about the challenge of uniform distribution of the flow to a multitude of channels while avoiding the need for duplication of peripheral components (see, for example, [49]). For singlephase flow, the internal flow distribution can be solved with the aid of computational fluid dynamics simulations and precise fabrication methods. Well-engineered continuous-flow systems can deliver productivities and yields comparable to or better than conventional batch systems over a wide range of production volumes relevant for specialty products. (B), economics (C) and overall TEP (D) in relation to the base case with aqueous GDH (black), comparable immobilization (red) and optimized system parameters (blue).

Supplementary Materials:
The following are available online at www.mdpi.com/2073-8994/13/3/524/s1. Figure S1: Development of the cross-section averaged concentration in the free channel volume in the flow direction for the basic Matlab model and F2D and F3D models, Figure  S2: Development of the cross section averaged concentration in the free channel volume in flow direction for the basic Matlab model , model F2D and model F3D, Figure S3: Genetic algorithm results for NADP/H concentration for case 1.3 and bead ratio for case 1.2, Table S1: Kinetic parameters for NDK reduction and cofactor regeneration reactions, Table S2: Material prices for CapEx estimations in €/g, Table S3: Cost of species for OpEx estimation in €/g. Author Contributions: Conceptualization, P.P. and R.D.; methodology, P.P. and M.K.; validation, P.P.; formal analysis, P.P. and M.K.; investigation, M.K.; writing-original draft preparation, P.P.; writing-review and editing, P.P., M.K. and R.D.; visualization, P.P.; supervision, M.K. and R.D.; project administration, R.D.; funding acquisition, R.D. All authors have read and agreed to the published version of the manuscript.   Figure S1: Development of the cross-section averaged concentration in the free channel volume in the flow direction for the basic Matlab model and F2D and F3D models, Figure S2: Development of the cross section averaged concentration in the free channel volume in flow direction for the basic Matlab model, model F2D and model F3D, Figure S3: Genetic algorithm results for NADP/H concentration for case 1.3 and bead ratio for case 1.2, Table S1: Kinetic parameters for NDK reduction and cofactor regeneration reactions, Table S2: Material prices for CapEx estimations in €/g, Table S3: Cost of species for OpEx estimation in €/g. Author Contributions: Conceptualization, P.P. and R.D.; methodology, P.P. and M.K.; validation, P.P.; formal analysis, P.P. and M.K.; investigation, M.K.; writing-original draft preparation, P.P.; writing-review and editing, P.P., M.K. and R.D.; visualization, P.P.; supervision, M.K. and R.D.; project administration, R.D.; funding acquisition, R.D. All authors have read and agreed to the published version of the manuscript.