Model-Based Design of Biochemical Microreactors

Mathematical modeling of biochemical pathways is an important resource in Synthetic Biology, as the predictive power of simulating synthetic pathways represents an important step in the design of synthetic metabolons. In this paper, we are concerned with the mathematical modeling, simulation, and optimization of metabolic processes in biochemical microreactors able to carry out enzymatic reactions and to exchange metabolites with their surrounding medium. The results of the reported modeling approach are incorporated in the design of the first microreactor prototypes that are under construction. These microreactors consist of compartments separated by membranes carrying specific transporters for the input of substrates and export of products. Inside the compartments of the reactor multienzyme complexes assembled on nano-beads by peptide adapters are used to carry out metabolic reactions. The spatially resolved mathematical model describing the ongoing processes consists of a system of diffusion equations together with boundary and initial conditions. The boundary conditions model the exchange of metabolites with the neighboring compartments and the reactions at the surface of the nano-beads carrying the multienzyme complexes. Efficient and accurate approaches for numerical simulation of the mathematical model and for optimal design of the microreactor are developed. As a proof-of-concept scenario, a synthetic pathway for the conversion of sucrose to glucose-6-phosphate (G6P) was chosen. In this context, the mathematical model is employed to compute the spatio-temporal distributions of the metabolite concentrations, as well as application relevant quantities like the outflow rate of G6P. These computations are performed for different scenarios, where the number of beads as well as their loading capacity are varied. The computed metabolite distributions show spatial patterns, which differ for different experimental arrangements. Furthermore, the total output of G6P increases for scenarios where microcompartimentation of enzymes occurs. These results show that spatially resolved models are needed in the description of the conversion processes. Finally, the enzyme stoichiometry on the nano-beads is determined, which maximizes the production of glucose-6-phosphate.

1. Introduction

One of the greatest challenges in biology is to understand the fundamental principles on how evolution has selected networks to fulfill specific functional needs in the control of metabolism or transcription. Synthetic biology approaches may help to shed light on such principles by identifying modular functional units of a network and uncover how units can be linked together to yield new function (Bashor et al., 2010). There has already been early success in engineering simple regulatory circuits that recapitulate some of the behaviors of natural regulatory circuits, like, e.g., circuits that regulate gene expression oscillations, bistable switches, or circuits that perform combinatorial logic operations, see Boyle and Silver (2009) and Purnick and Weiss (2009) for reviews. However, synthetic biology approaches are not limited to this field. An other important area of synthetic biology is in the development of synthetic organelles, which host metabolic processes, and which are able to communicate with the outside environment via transport processes over semipermeable membranes, which can either be built from natural constituents or from synthetic polymers. Such robust bioreactors can be, e.g., applied in biotechnology or drug delivery for the production of bioactive ingredients (Roodbeen and van Hest, 2009; Marguet et al., 2013).

A fundamental principle in the development of bioreactors is compartmentation. Hereby, the strategy is to mimic cellular organization, see, e.g., Roodbeen and van Hest (2009). The presence of compartments (organelles) inside living cells allows for a better regulatory control over the biological processes that occur inside these compartments, e.g., pathways competing for intermediates can occur in spatially separated compartments and can be regulated differently. Furthermore, microcompartmentation by means of metabolic channeling prevents the loss of intermediates and minimizes competing cross-reactions. However, clearcut experimental evidence demonstrating the importance of metabolic channeling for metabolic flux in vivo is lacking. Mimicking the natural situation, nanoreactors can be built by encapsulating enzymes in vesicular compartments [as summarized in Peters et al. (2012)], and in addition microcompartmentation of enzymes can be achieved with the help of synthetic protein scaffolds (Chen and Silver, 2012). Microreactors based on microfluidic devices carrying out enzymatic processes, reviewed in Nomura et al. (2004) and Asanomi et al. (2011), have also been into the focus of research in the past years. Several methods used to immobilize enzymes are available, like, e.g., immobilization of enzymes on magnetic microparticles. An approach in the development of microreactors for biosynthesis is the creation of fluidic assay-based microreactors with membrane-bounded subcompartments, carrying out biochemical conversion, and allowing exchange of substrates and products between individual subcompartments.

In our paper, we are focusing on the model-based design of biochemical microreactors able to carry out enzymatic reactions and to exchange metabolites between individual subcompartments. The bioreactor is built according to reported model predictions and is based on a microfluidic system consisting of chambers separated by porous walls, which are interspersed with biomembranes (Figure 1). The membranes carry specific transporters for input of substrates and export of products. Inside the chambers magnetic nano-beads are present that are positioned by an outer magnetic field.

FIGURE 1

Figure 1. A biochemical microreactor consisting of an array of compartments separated by membranes carrying specific transporters for input of substrates and export of products. This microreactor is currently under construction and the modeling results reported here are being utilized to influence the conceptual design of the microreactor.

Multienzyme complexes assembled on these nano-beads by peptide adapters are used to carry out metabolic reactions. The nano-beads allow a maximal enzyme concentration, which depends on the beads’ surface. Furthermore, it is possible to provide beads with a given enzyme stoichiometry. We exploit the microcompartmentation offered by the immobilized enzymes on the bead surface to address the question, if the spatial proximity of the individual enzymes and their stoichiometry has an impact on the productivity of the microreactor. In Figure 1, a microreactor consisting of an array of compartments is illustrated.

Our goal is to describe the spatio-temporal dynamics of metabolite concentrations involved in the biosynthetic conversions, and to optimize the microreactor, in order to increase the accumulation of the final product. To achieve this goal, we have developed a mathematical model describing the ongoing metabolic processes in spatial resolution. By applying a model with spatial resolution, we can assess, if the spatial arrangement of the individual enzymes and their stoichiometry has an impact on the productivity of the microreactor. Hence, our model will be able to predict, if metabolic channeling plays a role in the described in vitro assembly. Our model consists of a system of diffusion equations together with boundary and initial conditions. The boundary conditions model the exchange of metabolites with the neighboring chambers and the reactions at the surface of the nano-beads carrying the multienzyme complexes. For the numerical simulation of the mathematical model, efficient and accurate approaches are developed. These allow the computation of the spatio-temporal distribution of the metabolite concentrations inside the compartments and the flux of products through the export boundaries. Here, different computational scenarios can be considered including different numbers and loading capacities of nano-beads, different enzyme stoichiometries on the surface of the beads, as well as the distribution of enzymes both inside the fluid and on the beads’ surfaces. Comparing the metabolic flux through the system for those scenarios allow to test the hypothesis that microcompartmentation of enzymes increases the efficiency of the metabolic pathway. Based on the mathematical model and the simulation approach, the optimal design of the microreactor is performed.

As a proof-of-concept scenario, we choose a synthetic pathway for the conversion of sucrose to glucose-6-phosphate. This metabolic pathway is localized in one chamber of microreactor allowing the import of sucrose and the export of glucose-6-phosphate. Based on the mathematical model, several investigations are performed. First, the spatio-temporal distributions of the metabolite concentrations, as well as application relevant quantities like the production rate of the metabolites and outflow rate of G6P are computed. These computations are performed for different scenarios, where the number of beads as well as their loading capacity are varied. Furthermore, for the calculations, we consider two different hexokinases, namely HsHK2 and ScHK2, and we suppose sucrose and ATP to be present in surplus quantities. Finally, for the mentioned scenarios, we determine the stoichiometry of enzyme concentrations on the nano-beads which maximizes the production of glucose-6-phosphate.

2. Materials and Methods

2.1. The Mathematical Model

In this section, we set up a mathematical model describing the conversion of sucrose (S) to glucose-6-phosphate (G6P). This metabolic pathway is carried out in a microreactor consisting of a chamber separated from the surrounding medium by membranes carrying transport proteins for the input of sucrose and export of glucose-6-phosphate. Nano-beads loaded with multienzyme complexes are distributed inside the chamber. The enzymatic reactions constituting the metabolic pathway, as well as the corresponding enzymes and metabolites are listed in Table 1.

TABLE 1

Table 1. Metabolic reactions and the corresponding enzymes and metabolites.

The layout of the microreactor is given in Figure 2. We denote the reactor chamber by Ω_c ⊂ ℛⁿ (n = 2 or n = 3). The space inside the chamber occupied by nano-beads, which are balls with diameter d, is denoted by Ω_b, whereas the remainder, representing the domain occupied by the bulk solution, is denoted by $\mathrm{\Omega }:={\mathrm{\Omega }}_c\setminus \overline{{\mathrm{\Omega }}_b}$. We assume that Ω_b is strictly included in Ω_c, this means that the beads do not touch the walls of the chamber. The number of beads in the reactor is denoted by n_b. The boundary ∂Ω of the domain Ω is decomposed into ∂Ω_c, the boundary of the chamber, and Γ_b: = ∂Ω_b, the reactive surface of the nano-beads. Furthermore, the boundary ∂Ω_c of the chamber consists of Γ_i, Γ_e, and Γ₀, i.e.,

$\partial {\mathrm{\Omega }}_c=\overline{{\mathrm{\Gamma }}_i}\cup \overline{{\mathrm{\Gamma }}_e}\cup \overline{{\mathrm{\Gamma }}_0},$

FIGURE 2

Figure 2. Geometric representation of the modeled microreactor: reaction chamber Ω_c with boundary ∂Ω_c consisting of import boundary Γ_i, export boundary Γ_e, and impermeable boundary Γ₀, and ensemble of nano-beads Ω_b with reactive boundary Γ_b.

and the three sets are pairwise disjoint. The sets Γ_i, respectively, Γ_e represent the boundary parts where metabolites are transported into, respectively, out of the chamber Ω_c. These import/export boundaries have a complex geometric structure. They consist of fenestrations lined with lipid-membranes carrying transporters for the exchange of metabolites. On the boundary Γ_i the sucrose/proton cotransporters are distributed, whereas Γ_e contains the transporters for the exchange of glucose-6-phosphate and inorganic phosphate. In our model, the microscopic structure of the boundaries is taken into account in an averaged (homogenized) way, which makes the model amenable for numerical calculations. More precisely, we assign to each of the boundaries Γ_i, Γ_e an effective permeability for the transported metabolites denoted by θⁱ, respectively, θ^e (θⁱ, θ^e ∈ [0,1]). These permeabilities can be calculated by an averaging approach (homogenization), assuming that the pores of the transporters are very small compared to the dimension of the reactor chamber and occur in a large number, and that the transporters are uniformly distributed within the lipid-membrane. A sketch illustrating the idea behind the homogenization approach is given in Figure 3.

FIGURE 3

Figure 3. For pores of size ε > 0, the boundary part ${\mathrm{\Gamma }}_{\ast }^{\mathrm{ϵ}}$, with * ∈ {i,e} contains the union of all pores in the inflow, respectively, outflow boundary. When performing the homogenization approach, we suppose that the number of the pores gets larger, while the size of the pores gets smaller, in such a way that the ratio between the surface occupied by pores and the surface of the fenestration remains constant. The value of this constant is the permeability θ^*, * ∈ {i,e}, in the boundary conditions [equations (1c)–(1d)]. (A). Inflow/outflow boundary for pores with size ε₁. (B). Inflow/outflow boundary for pores with size ε₂ < ε₁. (C). Via homogenization, i.e., ε → 0, we obtain a homogenized model formulated on inflow/outflow boundaries without complex pore structure.

The spatio-temporal dynamics of the concentrations y_j, j = 1, …, m of metabolites involved in the metabolic pathway is governed by a system of reaction–diffusion equations of the form:

${\partial }_t{y}_j\left(t,x\right)-D_j\mathrm{\Delta }{y}_j\left(t,x\right)=R_j^{\mathrm{\Omega }}\left(y\left(t,x\right),\mathrm{\lambda }\right)\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}for}\left(t,x\right)\in \left(0,T\right)\times \mathrm{\Omega },$ (1a)

together with the boundary conditions

$\begin{array}{ll}\hfill & -D_j\nabla y_j\left(t,x\right)\cdot \nu \left(x\right)=-R_j^b\left(y\left(t,x\right),\mathrm{\lambda }\right)\text{\hspace{0.17em}for}\left(t,x\right)\in \left(0,T\right)\times {\mathrm{\Gamma }}_b,\end{array}$ (1b)

$\begin{array}{ll}\hfill & -D_j\nabla y_j\left(t,x\right)\cdot \nu \left(x\right)=-{\theta }^e{R}_j^e\left(y\left(t,x\right)\right)\text{for}\left(t,x\right)\in \left(0,T\right)\times {\mathrm{\Gamma }}_e,\end{array}$ (1c)

$\begin{array}{ll}\hfill & -D_j\nabla y_j\left(t,x\right)\cdot \nu \left(x\right)=-{\theta }^i{R}_j^i\left(y\left(t,x\right)\right)\text{for}\left(t,x\right)\in \left(0,T\right)\times {\mathrm{\Gamma }}_i,\end{array}$ (1d)

$\begin{array}{ll}\hfill & -D_j\nabla y_j\left(t,x\right)\cdot \nu \left(x\right)=0\text{for}\left(t,x\right)\in \left(0,T\right)\times {\mathrm{\Gamma }}_0,\end{array}$ (1e)

and the initial condition

$y_j\left(0,x\right)=y_j^0\left(x\right)\text{for}x\in \mathrm{\Omega }.$ (1f)

Here, ν denotes the outer unit normal on ∂Ω with respect to Ω. Equation (1a) describes diffusion with diffusivity D_j and enzymatic reactions with kinetics $R_j^{\mathrm{\Omega }}$ for the metabolite number j. In case that enzymes are localized on beads (and thus no reactions are carried out in the bulk domain), the terms $R_j^{\mathrm{\Omega }}$ are equal to zero. The boundary conditions involve the quantity −D_j▽y_j ⋅ ν which describes the normal component of the diffusive flux of the j-th metabolite at the boundary of the domain Ω. Thus, condition (1e) models an impermeable boundary, where the normal flux is equal to zero. Conditions (1c)–(1d) model the flux of metabolites through the import and export boundary, respectively. This flux is proportional to the permeability of the boundary, and the kinetics of the corresponding transport protein. Finally, condition (1b) describes the normal flux of metabolites at the boundary of the beads, generated by enzymatic reactions at the beads’ surface. We emphasize that the reaction rates $R_j^b$ and $R_j^{\mathrm{\Omega }}$ depend on an additional parameter $\mathrm{\lambda }=\left({\mathrm{\lambda }}_1,\dots ,{\mathrm{\lambda }}_{n_e}\right)\in {ℛ}^{n_e}$. This parameter describes the stoichiometry of enzymes on the beads and in the bulk, and thus for i = 1, …, n_e, where n_e denotes the number of enzyme species involved in the reactions, holds

${\mathrm{\lambda }}_i\ge 0\text{and}{\displaystyle \sum _{i=1}^{n_e}{\text{λ}}_i=1.}$ (2)

(The case λ_i* = 1 for some i* ∈ {1, … n_e} means that all binding sites on the bead surface are occupied by the enzyme i*, whereas the case ${\mathrm{\lambda }}_i=\frac{1}{n_e}$ for all i ∈ {1, … n_e} means that all enzymes occupy the same amount of binding sites).

The existence and uniqueness of positive solutions for the model (1) can be shown by arguments similar to those in Gahn et al. (under review)¹. We also emphasize that the model (1), valid for one conversion chamber, can easily be extended to model a multicompartment microreactor. This is done by adding transmission conditions at the interfaces separating the compartments, which model the exchange of metabolites between neighboring compartments.

In this paper, we want to investigate the interplay between the metabolic processes at the beads and the export of the product glucose-6-phosphate. We assume that during the conversion processes sucrose and ATP concentrations can be regarded as constant, since these two metabolite species are present in excess. As a consequence, we drop reaction (0). The equation for ADP can be neglected, since ADP is just a product of irreversible reactions and therefore has no direct effect on the reaction rates in equation (1). Finally, we assume that the concentrations of G6P_e and Pi_e (glucose-6-phosphate and inorganic phosphate outside the chamber) are constant. This is motivated by the fact that G6P_e may be consumed in a following reaction, and Pi_e may be delivered at a desirable rate in the space outside the chamber. Hence, we can drop the equations for G6P_e and Pi_e. The constant values of the concentrations mentioned above are given by the corresponding initial concentrations (Table 3). The scenario taking into account these assumptions shall be referred to as the sucrose excess scenario (SE-scenario).

The reactions relevant for the SE-scenario are reactions (1)–(5) (Table 1). These are in general multisubstrate enzymatic reactions. Their reaction mechanisms and the corresponding reaction kinetics are displayed in Table 2A. Note that in Table 2, the concentrations of metabolites are denoted with brackets (e.g., for sucrose, we use [S] instead of y_S). This is chosen in order to keep the notation clear.

TABLE 2

Table 2. (A) Metabolic reactions relevant for the sucrose excess scenario and the corresponding reaction kinetics; reaction mechanisms: (1) irreversible Michaelis–Menten; (2), (3) irreversible bi–bi ordered; (4) reversible Michaelis–Menten; (5) bi–bi ping pong. See, e.g., Segel (1975) for an overview on reaction mechanisms for multisubstrate enzymatic reactions; (B) reaction terms in equation (1) corresponding to the concentrations y_j, j = G, …, G6P, for the SE-scenario with enzymes distributed on beads, i.e., ${\mathit{\text{R}}}_{\mathit{\text{j}}}^{\mathrm{\Omega }}=0$.

The unknowns of the model in the SE-scenario are the following metabolite concentrations: y_G, y_F, y_F6P, y_Pi, y_G6P. For each unknown, a reaction–diffusion equation of type (1a) complemented by boundary conditions of type (1b)–(1e), and initial conditions holds. The reaction terms occurring in the equations are denoted by $R_{\text{G}}^{\mathrm{\Omega }},\dots ,R_{\text{G}6\text{P}}^{\mathrm{\Omega }}$, whereas the fluxes at the boundaries Γ_b, Γ_e, Γ_i are given by reaction terms denoted by $R_{\text{G}}^b,\dots ,R_{\text{G}6\text{P}}^b$, $R_{\text{G}}^e,\dots ,R_{\text{G}6\text{P}}^e$, and $R_{\text{G}}^i,\dots ,R_{\text{G}6\text{P}}^i$. The precise form of these reaction terms, in case when enzymes are distributed on beads, i.e., $R_j^{\mathrm{\Omega }}=0$, is given in Table 2B. We emphasize that, when the metabolite j, j = G, …, G6P participates in different reactions, the reaction term $R_j^b$ is given as a sum of all relevant reaction kinetics. If reactions take place also in the bulk, the reaction terms $R_j^{\mathrm{\Omega }}$, for j = G, …, G6P, have the same structure as $R_j^b$ with potentially different parameters. Finally, we mention that in the SE-scenario the inflow boundary Γ_i is impermeable, thus the reaction terms $R_{\text{G}}^i,\dots ,R_{\text{G}6\text{P}}^i$ are set to zero.

The values $v_{\text{max}}^i$ for an enzyme i, i ∈ {inv, hk, pgi} in Table 2, can be calculated by $v_{\text{max}}^i=k_{\text{cat}}^i{\left[E\right]}_0$, with a turnover number $k_{\text{cat}}^i$ given in Table 3, and the concentration [E]₀ of occupied binding sites on a bead. The enzymatic activity for the enzyme i,i ∈ {inv, hk, pgi} on each bead is then given by ${\mathrm{\lambda }}_i{v}_{\text{max}}^i$, where the vector λ = (λ_inv, λ_hk, λ_pgi) describes the enzymes stoichiometry on the bead. For the simulations, we choose $\mathrm{\lambda }={\mathrm{\lambda }}^E:=\left({\mathrm{\lambda }}_{\text{inv}}^E,{\mathrm{\lambda }}_{\text{hk}}^E,{\mathrm{\lambda }}_{\text{pgi}}^E\right)$, where λ^E satisfies

$\begin{array}{l}\hfill \begin{array}{rl}& k_{\text{cat}}^{\text{inv}}{\mathrm{\lambda }}_{\text{inv}}^E=k_{\text{cat}}^{\text{hk}}{\mathrm{\lambda }}_{\text{hk}}^E=k_{\text{cat}}^{\text{pgi}}{\mathrm{\lambda }}_{\text{pgi}}^E,\\ & {\mathrm{\lambda }}_{\text{inv}}^E+{\mathrm{\lambda }}_{\text{hk}}^E+{\mathrm{\lambda }}_{\text{pgi}}^E=1.\end{array}\end{array}$ (3)

TABLE 3

Table 3. Parameter values for the SE-scenario, corresponding to two different hexokinases HsHK2 and ScHK2.

This enzyme stoichiometry leads to equal enzymatic activity for all enzyme species. The values of λ^E are given in Table 3. Please note that in Section 2.3 optimal values for the parameter λ are computed.

2.2. Numerical Methods

We use a finite element method in order to find an approximation to the solution of equation (1) on a fixed time interval [0, T] and with a given initial state.

For the numerical accuracy and the computational complexity of an implementation, also in view of the optimization, it is crucial to choose a suitable discretization. For our problem, we use lowest order Raviart–Thomas elements (Raviart and Thomas, 1977). On a first glance, linear finite elements seem superior to Raviart–Thomas elements due to a lower number of unknowns and a higher order of convergence. However, for our problem, it turns out that the use of linear finite elements leads to solutions with negative concentrations when the triangulation𝒯_h is not extremely fine. More precisely, for n_b = 1 we investigated two scenarios, and compared the number of unknowns that are needed for each discretization, to obtain realistic results. For Ω_c = (0.50 μm)², we need 3200 degrees of freedom per species in order to obtain realistic results for linear finite elements, whereas Raviart-Thomas elements only need 993 degrees of freedom. For Ω_c = (0.300 μm)² we need 122,544 degrees of freedom per species in order to obtain realistic results for linear finite elements, whereas Raviart–Thomas elements only need 38,256 degrees of freedom. The computational complexity for Raviart–Thomas elements can be further reduced by hybridization, in this case only 607 degrees of freedom for the smaller domain and 23,016 degrees of freedom for the larger domain are needed.

The resulting space-discrete system is stiff, therefore, we use the implicit Euler method for the time integration. The resulting finite dimensional non-linear problems are solved with Newton’s method. Whenever the non-linear solver can not find a non-negative solution, the time step is rejected and a smaller time step size is used (until the non-linear solver can find a non-negative solution or the time step size is less than a predefined minimal time step size). We emphasize that decreasing the time step size was not necessary on the meshes used for spatial discretization (see Supplementary Material) with Raviart–Thomas elements and a time step size of 0.25 s.

In order to decrease the computational complexity of problem (1), arising from the number of species, we use the scheme presented in Kräutle and Knabner (2005) with minor modifications: instead of finding the basis of the image of the stoichiometric matrix S on the right hand side of the partial differential equation, we find the basis of the image of the matrix (S_b|S_e), where S_x refers to the stoichiometric matrix of the reactions on Γ_x, x ∈ {b,e}. The decoupled problem consists of 6 coupled species, whereas the original problem consists of 11 coupled species. In the SE-scenario, the decoupling scheme does not decrease the computational complexity, since the resulting stoichiometric matrix has full rank.

The numerical scheme has been implemented using the software package M++ that is based on the data structure in Wieners (2005). The algorithm uses MPI parallelization and can therefore use an arbitrary number of processors for computation. This is particularly useful for problems involving many species.

2.3. Optimization Methods

Besides setting up a mathematical model for the microreactor to describe the temporal and spatial distribution of enzymes and reactants, our goal also was to determine optimal parameters for enhanced performance of the microreactor based on this model. We demonstrate here how to succeed in a systematical way using derivative-based non-linear optimization techniques. Exemplary, we focus on how to determine the real parameters λ modeling the stoichiometry of the enzymes loaded on the beads entering in the non-linear terms in equation (1). We then discuss extensions of such methods to other parameters such as the discrete number of beads n_b, the combinatorial problem of which specific enzymes should be used for the metabolic pathway and continuous geometry parameters such as the bead diameters or locations, the total activity of the transporter proteins, the shape and size of the chamber itself, etc.

The criterion for the performance of the reactor will be the total outflow of a desired product, say $y_{j_{\ast }}$(λ) on Γ_e, over a fixed production time horizon [0,T],

$\begin{array}{llll}\hfill J\left(y\right)& ={\int }_0^T{\int }_{{\mathrm{\Gamma }}_e}-D_{j_{\ast }}\nabla y_{j_{\ast }}\left(t,x\right)\cdot \nu \left(x\right)d\mathrm{\sigma }\left(x\right)\mathit{dt}& \hfill & \\ \hfill & ={\int }_0^T{\int }_{{\mathrm{\Gamma }}_e}-{\theta }^e{R}_{j_{\ast }}^e\left(y\left(t,x\right)\right)d\mathrm{\sigma }\left(x\right)\mathit{dt},\end{array}$ (4)

where we have used equation (1c) for $j=j_{\ast }$ in equation (4). For our sucrolytic chamber, we have $j_{\ast }=j_{\text{G6P}}$. By optimal parameters, we mean that any feasible choice close to the optimal one does not lead to a higher total outflow predicted by the model.

While derivative-based non-linear optimization techniques are conceptually well-known, we emphasize that due to the large-scale, non-linear system of equations to be solved for each direct simulation of the model, the main challenge is an efficient computation of the derivative of the reduced cost function

$Ĵ\left(\mathrm{\lambda }\right)=\left(y\left(\mathrm{\lambda }\right)\right),$ (5)

where y(λ) denotes the numerical solution of the system [equations (1a)–(1f)] as a function of the parameter λ. We will exploit that the dimension of the parameter space n_p is typically small compared to the dimension obtained from discretization of the infinite-dimensional state space for y in the model and will therefore follow a sensitivity-based approach to compute the derivative [see, e.g., Hinze et al. (2009)]. Letting e_i, i = 1, …, n_p, being a bases of the parameter space this means that we can compute the derivative $Ĵ\prime \left(\mathrm{\lambda }\right)$ as

$Ĵ\prime \left(\mathrm{\lambda }\right)={\displaystyle \sum _{i=1}^l}{J}_y\left(y\left(\mathrm{\lambda }\right)\right){\delta }_{e_i}y,$ (6)

where the sensitivity ${\delta }_{e_i}y=y\prime \left(\mathrm{\lambda }\right)e_i$ is given by the solution (${\tilde{y}}_1,\dots ,{\tilde{y}}_m$) of the following linearized problem

$\begin{array}{llll}\hfill & {\partial }_t{\tilde{y}}_j\left(t,x\right)-D_j\mathrm{\Delta }{\tilde{y}}_j\left(t,x\right)={\partial }_y{R}_j^{\mathrm{\Omega }}\left(y\left(t,x\right),\mathrm{\lambda }\right){\tilde{y}}_j& \hfill & \\ \hfill & -{\partial }_{\mathrm{\lambda }}{R}_j^{\mathrm{\Omega }}\left(y\left(t,x\right),\mathrm{\lambda }\right)e_i\text{for}\left(t,x\right)\in \left(0,T\right)\times \mathrm{\Omega },\end{array}$ (7a)

together with the boundary conditions

$\begin{array}{llll}\hfill & -D_j\nabla {\tilde{y}}_j\left(t,x\right)\cdot \nu \left(x\right)=-{\partial }_y{R}_j^b\left(y\left(t,x\right),\mathrm{\lambda }\right){\tilde{y}}_j& \hfill & \\ \hfill & +{\partial }_{\mathrm{\lambda }}{R}_j^b\left(y\left(t,x\right),\mathrm{\lambda }\right)e_i\text{for}\left(t,x\right)\in \left(0,T\right)\times {\mathrm{\Gamma }}_b,\end{array}$ (7b)

$\begin{array}{llll}\hfill & -D_j\nabla {\tilde{y}}_j\left(t,x\right)\cdot \nu \left(x\right)=& \hfill & \\ \hfill & -{\theta }^e\left(R_j^e\right)\prime \left(y\left(t,x\right)\right){\tilde{y}}_j\text{for}\left(t,x\right)\in \left(0,T\right)\times {\mathrm{\Gamma }}_e,\end{array}$ (7c)

$\begin{array}{llll}\hfill & -D_j\nabla {\tilde{y}}_j\left(t,x\right)\cdot \nu \left(x\right)=& \hfill & \\ \hfill & -{\theta }^i\left(R_j^i\right)\prime \left(y\left(t,x\right)\right){\tilde{y}}_j\text{for}\left(t,x\right)\in \left(0,T\right)\times {\mathrm{\Gamma }}_i,\end{array}$ (7d)

$\begin{array}{ll}\hfill & -D_j\nabla {\tilde{y}}_j\left(t,x\right)\cdot \nu \left(x\right)=0\text{for}\left(t,x\right)\in \left(0,T\right)\times {\mathrm{\Gamma }}_0,\end{array}$ (7e)

and the initial condition

${\tilde{y}}_j\left(0,x\right)=0\text{for}x\in \mathrm{\Omega },$ (lf)

where ∂_y and ∂_λ denote partial derivatives and ()′ [e.g., ($R_j^i$)′] denotes total derivatives of the corresponding kinetic functions, respectively, and where y is the solution of equations (1a)–(1f). The solution of this linearized problem can be computed simultaneously with the direct simulation using the same discretization method as described in Section 2.2. Compared to the evaluation of the objective function (5), the additional computational effort to obtain the derivative is solving n_p linear equation systems in each time step. Instead of using the linearized problem (7), we can also obtain a derivative from an adjoint problem, see again, e.g., Hinze et al. (2009). The latter approach is more efficient when n_p is large.

With the derivative at hand, we can then solve the reduced optimization problem subject to further parameter constraints for instance with primal-dual interior-point algorithms using a quasi-Newton approximation of the Hessian. These methods are known to have both good theoretical and practical behavior, see, e.g., Forsgren et al. (2002). Given some initial parameter λ⁽⁰⁾, they compute a sequence of parameters λ^(k), k = 1, 2, 3, …, converging to an optimal λ until for some k* first order conditions or stationarity in the objective function are satisfied within a tolerance tol_X or tol_fun, respectively. We set λ* = λ^(k*), J* = J(y(λ*)), and J⁰ = J(y(λ⁽⁰⁾)).

In order to determine an optimal stoichiometry λ for our sucrolytic chamber, we have the three parameters λ_inv, λ_hk, and λ_pgi, but we may eliminate for example λ_pgi by the condition (2). This yields n_p = 2, together with a linear inequality constraint λ_inv + λ_hk ≤ 1 and box constraints λ_inv, λ_hk ∈ [0,1]. For our numerical results in Section 3.2 for the SE-scenario, we have used the interior-point algorithm implemented in the Optimization Toolbox in MATLAB (2014) with a BFGS quasi-Newton approximation of the hessian, where equations (5) and (6) are computed using the methods described in Section 2.2. Again we use the software package M++ for our implementation. As initial parameter λ⁽⁰⁾, we either use the stoichiometry from equal activity λ^E given by equation (3) or the uniform stoichiometry ${\mathrm{\lambda }}^{\left(0\right)}={\mathrm{\lambda }}^U=\left(\frac{1}{3},\frac{1}{3},\frac{1}{3}\right)$ depending on which has a larger value J⁽⁰⁾. Furthermore, we use tol_X = 1.0e–06 and tol_fun = 1.0e–04 for all our optimization runs.

The derivative computation and hence the above method can directly be adapted to other real parameters entering the kinetic functions such as, e.g., the permeability θ. Combinatorial problems from (additional) discrete parameters such as the number of beads can be solved using enumeration or, approximately, by derivative-based optimization techniques based on relaxation methods, even for time dependent parameters as in Hante and Sager (2013). Continuously varying geometry parameters such as bead diameters or shape and size of the chamber Ω_c can be handled similarly by shape derivatives as in Leugering et al. (2011). For the SE-scenario, we enumerate optimal λ* for various combinations of number of beads n_b, different choices of hexokinases in the sucrolytic pathway and different time horizons T.

3. Results

We apply our simulation and optimization methods to the situation of the microreactor in the SE-scenario for realistic parameters. Hereby, we consider two-dimensional quadratic reaction chambers Ω_c = (0, L)² with L = 500 and L = 3000 μm, and different numbers of uniformly distributed beads, each of diameter of 10 μm. The outflow boundary is Γ_e = {L} × (0, L). The rest of the chambers’ boundary is impermeable. The parameters of the mathematical model are given in Table 3. The final time T varies between 0.5 and 4 h. This is appropriate, due to the fact that the outflow-membrane Γ_e is in general unstable and bursts after several hours. We stress that, due to the need of small time step sizes and the ratio of bead and chamber size, larger time horizons and chamber sizes lead to a higher computational complexity.

For our numerical investigations, we consider two types of hexokinase, namely HsHK2 and ScHK2. We also have two stoichiometry values which play a role in our calculations: the uniform stoichiometry ${\mathrm{\lambda }}^U=\left(\frac{1}{3},\frac{1}{3},\frac{1}{3}\right)$, and the equal activity stoichiometry λ^E described in equation (3). More precisely, λ^E is used for the numerical simulations, whereas either of the two is used as initial parameter in the optimization, depending on which leads to a larger value of the total glucose-6-phosphate outflow.

3.1. Simulation Results

We consider a quadratic chamber Ω_c with a length of 3000 μm, a fixed run time of T = 3 h, a fixed enzyme stoichiometry λ = λ^E and a number of beads n_b ∈ {0², 1², 2², … ,6², 12², 18², … ,60²}, where n_b = 0 refers to a scenario with enzymes distributed in the bulk solution.

The meshes used for spatial discretization consist of between 15,190 triangles and 22,941 edges for 1 bead, and 93,184 triangles and 151,592 edges for 2916 beads. Detailed information about the meshes is available in the Supplementary Material. We used a time step size of 0.25 s. All calculations have been repeated for T = 1.5 h on uniform refined meshes and half of the time step size. No significant differences between the results on the coarse and on the fine meshes were observed.

First results are concerned with the numerical computation of the spatial and temporal distribution of the metabolites. Figure 4 shows the distribution of G6P for a scenario with 1296 beads and the hexokinase HsHK2 after 3 h for completely and partially loaded beads.

FIGURE 4

Figure 4. G6P concentration in $\frac{\mathit{\text{mol}}}{m^{\text{2}}}$ in a scenario with n_b = 1296 beads and the hexokinase HsHK2, for a completely loaded beads (left) and partially loaded beads, where the total number of enzymes distributed to the beads corresponds to that on one completely loaded bead (right), after 3 h.

Next results are concerned with the computation of quantities, which are particularly relevant for our microreactor, namely the export of G6P at time t ∈ [0,T], i.e., the outflow rate at time t given by

${\int }_{{\mathrm{\Gamma }}_e}-D_{\text{G}6\text{P}}\nabla y_{\text{G}6\text{P}}\left(t,x\right)\cdot \nu \left(x\right)d\mathrm{\sigma }\left(x\right),$ (8)

and the production rate for the metabolites at the beads at time t given by

${\int }_{{\mathrm{\Gamma }}_b}{D}_j\nabla y_j\left(t,x\right)\cdot \nu \left(x\right)d\mathrm{\sigma }\left(x\right),$ (9)

for j = G, …, G6P. In computing these quantities, we consider the following arrangements:

(S.1) Increase the number of beads n_b, where beads are completely loaded.

(S.2) The total number of enzymes in the domain is fixed (and given by the number of binding sites on one completely loaded bead) and enzymes are distributed uniformly on n_b beads or in bulk. The latter case, we characterize by n_b = 0.

Both arrangements are simulated for the two hexokinases HsHK2 and ScHK2, and the beads are distributed periodically in the bulk, like, e.g., in Figure 4. For the arrangement (S.2), we also consider the case without beads where all enzymes are distributed in the bulk. Then the reaction terms $R_j^{\mathrm{\Omega }}$ in equation (1a) are equal to $R_j^b$ from Table 2. In (S.1), the maximal velocity $v_{\max }^i=k_{\text{cat}}^i{\left[E\right]}_0$ corresponding to enzyme i, i ∈ {inv, hk, pgi} is calculated with [E]₀ = [E]_full where [E]_full is the concentration of binding sites on one bead. These values are denoted by $v_{\max }^{i,\text{full}}$, and are given in Table 3. In (S.2), we fix the total concentration of enzyme to be given by [E]_full. However, this concentration is now distributed on n_b beads. Thus, the values for $v_{\max }^i,i\in \left\{\text{inv},\text{hk},\mathit{\text{pgi}}\right\}$ are now obtained by dividing the values used in (S.1) by the number of beads n_b.

The production and outflow rate for G6P are plotted in Figure 5 for (S.1) and Figure 6 for (S.2). For fully loaded beads the production rate increases with increasing number of beads, and reaches saturation after approximately 1000 s in the case of HsHK2 and after approximately 10,800 s in the case of ScHK2 (Figure 5, upper pictures). Concerning the outflow rate, in case of HsHK2 (Figure 5 right bottom) an upper bound is asymptotically reached, the value of which seams to be the same for all n_b ≥ 576 (simulations with longer runtime suggest that the value seams to be the same for n_b ≥ 16).

FIGURE 5

Figure 5. G6P production rate (top) and outflow rate (bottom) in $\frac{\mathit{\text{mol}}}{s}$ for arrangement (S.1) with n_b ∈ {1, 144, 576, 1296, 2916} and the enzymes ScHK2 (left) and HsHK2 (right).

FIGURE 6

Figure 6. G6P production rate (top) and outflow rate (bottom) in $\frac{\mathit{\text{mol}}}{s}$ for the arrangement (S.2) n_b ∈ {0, 1, 4, 36, 900} and the enzymes ScHK2 (left) and HsHK2 (right).

In case of partially loaded beads corresponding to arrangement (S.2), the production rate of G6P at the beads up to a time of 3 h first decreases with increasing number of beads and then stays nearly constant (Table 4; Figure 6, upper figures). Please note, that for the arrangement (S.2) the outflow rates are not reaching saturation, which implies that the transporters activity are not limiting the outflow for the chosen enzyme concentration [E]_full. We also remark that for the chosen stoichiometry λ^E, in case of arrangement (S.2), the production rate of G6P for HsHK2 is by a factor 10³ to 10⁴ larger than for ScHK2, i.e., in (S.2) the selected hexokinase has a greater influence on the conversion process as in (S.1).

TABLE 4

Table 4. Production and outflow rates at T = 3h for scenario (S.2) with the hexokinase HsHK2 and n_b ∈ {1², 2², 6², 18², 30², 42², 54², 0}, and times until half of these rates are reached.

In order to quantify how fast production and outflow rate for various number of beads increase in time in the scenario (S.2) with the hexokinase HsHK2, in Table 4, we give the values of these rates at T = 3 h, and the times the reactor needs to reach half of these values.

3.2. Optimization Results

We aim at finding the optimal stoichiometry λ* for the SE-scenario. Since our focus is on computing a variety of different scenarios with a reasonable computational effort, we have chosen the smaller chamber Ω_c = (0.500 μm)² for these experiments. The meshes used for the spatial discretization of this chamber consist of between 820 triangles and 1270 edges for 1 bead, and 2394 triangles and 3883 edges for 64 beads. Further information about the meshes is available in the Supplementary Material.

For our study, we consider the following arrangements:

(O.1) For a fixed runtime of T = 1 h, we increase the number of completely loaded beads in the chamber with n_b ∈ {1, 4, 9, 16, 25, 36, 49, 64}.

(O.2) For a sucrolytic chamber with one completely loaded bead, we optimize the G6P outflow for different runtimes T ∈ {0.5, 1, 1.5, 2, … 4} hours.

We study these arrangements for the hexokinases ScHK2 and HsHK2.

For the enzyme ScHK2, we find in case of arrangement (O.1) that the optimal stoichiometry for a runtime of 1 h is between 3 and 12% for invertase, between 82 and 94% for hexokinase and between 3 and 6% for phosphoglucose isomerase. With an increasing number of beads, the proportion of ScHK2 increases and the proportions of invertase and phosphoglucose isomerase decrease. The exact results are listed in Table 5, where for each number of beads n_b, the optimal stoichiometry λ* as well as the corresponding total outflow J* are displayed. From the arrangement (O.2), we find that the optimal stoichiometry for one bead is between 8 and 13% for invertase, between 80 and 87% for hexokinase and between 5 and 6% for phosphoglucose isomerase, again with an increase in the proportion of binding sites loaded with ScHK2 in expense of proportion for ScINV, but with an almost constant proportion of binding sites loaded with ScPGI when the runtime is successively increased. The exact results are listed in Table 6, where for each run time T, the optimal stoichiometry λ* as well as the corresponding total outflow J* are displayed.

TABLE 5

Table 5. Optimal parameters λ* for the SE-scenario with ScHK2 (top) and HsHK2 (bottom), for a fixed runtime T = 3600, and different numbers n_b of beads.

TABLE 6

Table 6. Optimal parameters λ* for the SE-scenario with ScHK2 (top) and HsHK2 (bottom), for one bead and different runtimes T on a mesh with n_h = 820.

In the presence of the hexokinase HsHK2, we find from the first arrangement that the optimal stoichiometry for a runtime of 1 h is between 5 and 9% for invertase, between 91 and 95% for hexokinase and approximately 0% for the phosphoglucose isomerase. With an increasing number of beads, the proportion of hexokinase increases and the proportion for invertase decreases. The exact results are listed in Table 5. From the arrangement (O.2), we find that the optimal stoichiometry for one bead is between 7 and 10% for invertase and between 90 and 93% for the hexokinase again with an increase in the proportion of binding sites loaded with HsHK2 in expense of the proportions for ScINV. As in arrangement (O.1), the optimal proportion of phosphoglucose isomerase is approximately 0%. The exact results are again listed in Table 6. We stress that the optimization results in Table 5 correspond to local maxima.

4. Discussion

In this paper, we developed mathematical techniques for the simulation and optimization of metabolic processes in biological microreactors with membrane-bounded subcompartments. The results from our models will be incorporated into the design of a microreactor and will be validated by experiments with this microreactor. The microreactor is based on a microfluidic system, and consists of chambers separated by membranes which carry specific transporters for input of substrates and export of products. Inside the chambers, magnetic nano-beads carrying multienzyme-complexes are distributed, and immobilized by an outer magnetic field.

We have set up a mathematical model describing the spatio-temporal dynamics of the metabolite concentrations involved in the assembled metabolic pathway, and develop efficient and accurate approaches for the numerical simulation of this model. Furthermore, we show that model-based optimization for such systems is feasible by the methods presented. These approaches are applied to the proof-of-concept microreactor carrying out the synthetic pathway for the conversion of sucrose to glucose-6-phosphate. The obtained results shed light on the following important questions linked to the design and functionality of microreactors.

4.1. Are Spatially Resolved Models Necessary in the Description of Biochemical Microreactors?

The simulation of the distribution of metabolites, e.g., that of G6P in the presence of the enzyme HsHK2 (Figure 4), shows spatial patterns, which differ for different experimental arrangements. In case of completely loaded beads (Figure 4, left) the highest concentration of G6P can be found around the beads. If the beads are just partially loaded (Figure 4, right), a concentration gradient across the chamber prevails. In this case, the total number of enzymes distributed to the beads corresponds to the number of enzymes on one completely loaded bead. These spatial patterns in the metabolite distributions confirm the demand of spatially resolved models in the description of metabolic processes carried out in synthetic microreactors and, ultimately in living cells, which are more complex in architecture. Furthermore, the limiting effect of diffusion also shows that spatial effects have to be taken into account. This is visible in arrangement (S.2), where in the initial phase of the simulation the production rate of G6P at the beads decreases with increasing number of beads (Figure 6, upper figures); however, the outflow rate of G6P is higher for a large number of beads. This results from the fact that for high number of beads there are more beads close to the outflow boundary Γ_e, and consequently G6P diffusion to the boundary Γ_e needs less time. We mention, however, that for longer runtimes, the higher production rate of G6P for lower numbers of beads is predominating, and the outflow rates of G6P are higher for less number of beads (Figure 6, lower pictures).

4.2. Has Spatial Organization and Microcompartmentation of Enzymes the Potential to Influence the Efficiency of Metabolic Pathways?

Comparing the outflow rate of the product G6P for a fixed total enzyme quantity which is distributed on different numbers of beads, it turns out that the production rate can be increased with decreasing number of beads (Table 4; Figure 6, lower figures), especially, assembling all enzymes on one bead leads to the maximal outflow of G6P. This remains true also after comparing the upper scenarios with the scenario, in which enzymes are distributed in the bulk fluid within the reactor chamber. This last scenario turns out to lead to the lowest G6P outflow. For the scenario with one bead and the hexokinase HsHK2, the outflow of G6P is approximately 60% higher than the total outflow for the scenario with enzymes distributed within the bulk fluid with the same hexokinase (Table 4). These results suggest that microcompartmentation of enzymes increases the total output of G6P, not only in the model but also in vitro.

4.3. What Are the Limitations for the Productivity of the Modeled (Synthetic) Membrane-Bounded Microreactors?

The simulation results for the arrangement (S.1) (with completely loaded beads) show that the outflow rates for G6P reach saturation, in spite of the fact that the activity inside the chamber increases with increasing number of beads. This can directly be seen for the hexokinase HsHK2 (Figure 5, right bottom), whereas for ScHK2 this effect does not occur up to the run time of 3 h (Figure 5, left bottom). However, our computations show that at later time, this limiting effect also occurs for this hexokinase. This suggests that the limiting factor of the G6P outflow is the activity of the transporters in the membrane Γ_e. Finally, we mention that for the arrangement (S.2) the outflow rates are not reaching saturation, which implies that the activity of the transporters is not limiting the outflow for the chosen enzyme concentration [E]_full.

4.4. Can Optimal Design Help to Improve the Yield of a Biosynthetic Pathway?

Our optimization results for the microreactor in the SE-scenario with the hexokinase ScHK2, show that after a run time T = 3600 s, the total outflow of G6P, denoted by J*, for the optimal stoichiometry λ* is on average twice as much as the total outflow J⁰ for the uniform stoichiometry λ^U, and even much more than for the stoichiometry λ^E, e.g., 2746% in case of one bead. For one bead, this improvement is also visualized in Figure 7 (left). When run time is successively increased (for a reactor with one bead) starting from T = 1800 s, again the optimal stoichiometries enhance the total outflow of the desired product G6P on average twofold compared to the uniform stoichiometry and even much more compared to the stoichiometry λ^E. The comparison of the G6P outflow over time for stoichiometries computed for different runtimes shows that a larger proportion of binding sites loaded with ScINV leads to a higher G6P outflow at the initial phase while a larger proportion of binding sites loaded with ScHK2 leads to a higher G6P outflow on the long run (Figure 8).

FIGURE 7

Figure 7. G6P outflow for the SE-scenario with n_b = 1 compared for the uniform (λ^U), equal activity (λ^E) and optimal (λ*) stoichiometry for T = 7200 and ScHK2 (left) and HsHK2 (right) computed with Δt = 0.25 and n_h = 820.

FIGURE 8

Figure 8. Predicted difference of in the G6P outflow for the SE-scenario with n_b = 1 up to T = 7200 s comparing the optimal stoichiometry λ* obtained for T = 7200 and that for T = 1800 both for ScHK2 (dash-dotted line) and HsHK2 (solid line) from Table 6.

Optimization of the microreactor setup with the hexokinase HsHK2 suggests that the phosphoglucose isomerase can be eliminated for the run times considered in the optimization scenarios. Furthermore, in contrast to the results for ScHK2, the degree of improvement of the optimal stoichiometries is between 50% (for one bead) and 1% (for 64 beads) compared to the stoichiometry computed from equal activity and slightly higher compared to the uniform stoichiometry. For one bead, the improvement is also visualized in Figure 7 (right). When time is successively increased (for a reactor with one bead) from T = 1800–14,400 s, the optimal stoichiometries now enhance the total outflow of G6P. This enhancement is 20% up to 78% compared to the stoichiometry computed from equal activity and slightly more compared to the uniform stoichiometry, e.g., 147% for T = 1800. The observation obtained for ScHK2 concerning the role of the invertase and hexokinase for the initial phase and on the long run is now even more pronounced (Figure 8).

To understand the optimization result obtained for the proportion of phosphoglucose isomerase (the enzyme catalyzing the reaction G6P ⇌ F6P), in case of the microreactor with HsHK2, namely that this enzyme should be better eliminated, we simulated the production rates for all metabolites, for one completely loaded bead (Figure 9). We see that, with respect to our considered time interval, the production rates for G and F increase very fast to a maximum-level around 0.3 $\frac{\mathit{mol}}{s}$. The reaction rates in Table 2 imply that the conversion of S into F is constant, since the concentration of S is constant. Furthermore, Figure 9 shows that the production rate of F6P increases in time, but the production rate of F remains nearly constant. Hence, we conclude that there is a metabolic flux from G6P to F6P, what contradicts the reactor to obtain a maximal outflow of G6P, and should be avoided by considering a small activity of phosphoglucose isomerase.

FIGURE 9

Figure 9. Production rates in $\frac{\mathit{\text{mol}}}{s}$ for one fully loaded bead in the SE-scenario with the hexokinase HsHK2.

We stress that none of the reported maxima can be guaranteed to be global optima. Nevertheless, concerning the quality of the maxima, we note that all search paths from λ⁽⁰⁾ to λ* transverse a large part of the parameter space and are associated with a significant decrease in the costs. Moreover, we have tested various initial points λ⁽⁰⁾ for the scenario with one bead, HsHK2 and T = 3600. In all cases, we have observed convergence to the same λ*. We conclude that the reported local maxima are reasonable candidates for an optimized stoichometry for the respective microreactor.

Finally, we emphasize that the simulation and optimization approaches developed in this work can be repeated with minor adaptions for more complex biochemical microreactors and for other optimization parameters, like, e.g., other real parameters entering the reaction kinetics.

Author Contributions

MN-R and LV gave substantial contributions to the Abstract and Introduction. MG and MN-R gave substantial contributions to the conception of the Derivation of the Mathematical Model, see Section 2.1. TE and PK gave substantial contributions to the conception of the Numerical Analysis and Simulations, see Section 2.2 and 3.1. FH and GL gave substantial contributions to the conception of the Optimization Methods and Results, see Section 2.3 and 3.2. LV gave substantial contributions to the acquisition, analysis, and interpretation of data for the work. All authors contributed essentially in the discussion of the results, see Section 4. Also, all authors have revisited the work critically for important intellectual content, approved the version to be published, and agreed to be accountable for all aspects of the work.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The work of the first and the second author was supported by the Emerging Fields Initiative (EFI) for Synthetic Biology at the Friedrich-Alexander-Universität Erlangen-Nürnberg.

Supplementary Material

The Supplementary Material for this article can be found online at http://journal.frontiersin.org/article/10.3389/fbioe.2016.00013

Footnote

1. ^Gahn, M., Neuss-Radu, M., and Knabner, P. (2015). Homogenization of reactiondiffusion processes in a two-component porous medium with nonlinear flux conditions at the interface. SIAM J. Appl. Math. (Under Review).

References

Asanomi, Y., Yamaguchi, H., Miyazaki, M., and Maeda, H. (2011). Enzyme-immobilized microfluidic process reactors. Molecules 16, 6041–6059. doi: 10.3390/molecules16076041

PubMed Abstract | CrossRef Full Text | Google Scholar

Bashor, C. J., Horwitz, A. A., Peisajovich, S. G., and Lim, W. A. (2010). Rewiring cells: synthetic biology as a tool to interrogate the organizational principles of living systems. Annu. Rev. Biophys. 39, 515. doi:10.1146/annurev.biophys.050708.133652

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyle, P. M., and Silver, P. A. (2009). Harnessing nature’s toolbox: regulatory elements for synthetic biology. J. R. Soc. Interface 6(Suppl. 4), S535–S546. doi:10.1098/rsif.2008.0521.focus

CrossRef Full Text | Google Scholar

Chen, A. H., and Silver, P. A. (2012). Designing biological compartmentalization. Trends Cell Biol. 22, 662–670. doi:10.1016/j.tcb.2012.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Forsgren, A., Gill, P. E., and Wright, M. H. (2002). Interior methods for nonlinear optimization. SIAM Rev. 44, 525–597. doi:10.1137/S0036144502414942

CrossRef Full Text | Google Scholar

Gao, H., and Leary, J. A. (2003). Multiplex inhibitor screening and kinetic constant determinations for yeast hexokinase using mass spectrometry based assays. J. Am. Soc. Mass Spectrom. 14, 173–181. doi:10.1016/S1044-0305(02)00867-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Gloyn, A., Odili, S., Zelent, D., Buettger, C., Castleden, H., Steele, A., et al. (2005). Insights into the structure and regulation of glucokinase from a novel mutation (v62m), which causes maturity-onset diabetes of the young. J. Biol. Chem. 280, 14105–14113. doi:10.1074/jbc.M413146200

PubMed Abstract | CrossRef Full Text | Google Scholar

Hante, F. M., and Sager, S. (2013). Relaxation methods for mixed-integer optimal control of partial differential equations. Comput. Optimization Appl. 55, 197–225. doi:10.1007/s10589-012-9518-3

CrossRef Full Text | Google Scholar

Hinze, M., Pinnau, R., Ulbrich, M., and Ulbrich, S. (2009). Optimization with PDE Constraints, Volume 23 of Mathematical Modelling: Theory and Applications. New York: Springer.

Google Scholar

Kräutle, S., and Knabner, P. (2005). A new numerical reduction scheme for fully coupled multicomponent transport-reaction problems in porous media. Water Resour. Res. 41, 09414. doi:10.1029/2004WR003624

CrossRef Full Text | Google Scholar

Lafraya, A., Sanz-Aparicio, J., Polaina, J., and Marin-Navarro, J. (2011). Fructo-oligosaccharide synthesis by mutant versions of Saccharomyces cerevisiae invertase. Appl. Environ. Microbiol. 77, 6148–6157. doi:10.1128/AEM.05032-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Leugering, G., Novotny, A. A., Perla Menzala, G., and Sokołowski, J. (2011). On shape optimization for an evolution coupled system. Appl. Math. Optimization 64, 441–466. doi:10.1007/s00245-011-9148-7

CrossRef Full Text | Google Scholar

Lin, H.-Y., Kao, Y.-H., Chen, S.-T., and Meng, M. (2009). Effects of inherited mutations on catalytic activity and structural stability of human glucose-6-phosphate isomerase expressed in Escherichia coli. Biochim. Biophys. Acta 1794, 315–323. doi:10.1016/j.bbapap.2008.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Marguet, M., Bonduelle, C., and Lecommandoux, S. (2013). Multicompartmentalized polymeric systems: towards biomimetic cellular structure and function. Chem. Soc. Rev. 42, 512–529. doi:10.1039/c2cs35312a

PubMed Abstract | CrossRef Full Text | Google Scholar

MATLAB and Statistic Toolbox Release. (2014). The Mathworks, Inc., Natick, MA, USA.

Google Scholar

Nomura, A., Shin, S., Mehdi, O. O., and Kauffmann, J.-M. (2004). Preparation, characterization, and application of an enzyme-immobilized magnetic microreactor for flow injection analysis. Anal. Chem. 76, 5498–5502. doi:10.1021/ac049489v

PubMed Abstract | CrossRef Full Text | Google Scholar

Peters, R. J., Louzao, I., and van Hest, J. C. (2012). From polymeric nanoreactors to artificial organelles. Chem. Sci. 3, 335–342. doi:10.1039/C2SC00803C

CrossRef Full Text | Google Scholar

Purnick, P. E., and Weiss, R. (2009). The second wave of synthetic biology: from modules to systems. Nat. Rev. Mol. Cell Biol. 10, 410–422. doi:10.1038/nrm2698

PubMed Abstract | CrossRef Full Text | Google Scholar

Raviart, P. A., and Thomas, J. M. (1977). A mixed finite element method for second-order elliptic problems. Math. Aspects Finite Elem. Method 606, 292–315. doi:10.1007/BFb0064470

CrossRef Full Text | Google Scholar

Roodbeen, R., and van Hest, J. (2009). Synthetic cells and organelles: compartmentalization strategies. Bioessays 31, 1299–1308. doi:10.1002/bies.200900106

PubMed Abstract | CrossRef Full Text | Google Scholar

Segel, I. H. (1975). Enzyme Kinetics. New York: John Wiley & Sons, Inc.

Google Scholar

Somarowthu, S., Brodkin, H., D’Aquino, J., Ringe, D., Ondrechen, M., and Beuning, P. (2011). A tale of two isomerases: compact versus extended active sites in ketosteroid isomerase and phosphoglucose isomerase. Biochemistry 50, 9283–9295. doi:10.1021/bi201089v

PubMed Abstract | CrossRef Full Text | Google Scholar

Wieners, C. (2005). “Distributed point objects. a new concept for parallel finite elements,” in Domain Decomposition Methods in Science and Engineering, Volume 40 of Lecture Notes in Computational Science and Engineering, eds T. Barth, M. Griebel, D. Keyes, R. Nieminen, D. Roose, T. Schlick, et al. (Berlin, Heidelberg: Springer), 175–182.

Google Scholar