A method for zooming of nonlinear models of biochemical systems

Basic Definitions and Assumptions

is then commonly used to reduce the model. We refer to a state as fast in the time interval ≤ T₀ ≤ t < T₁ if Eq. (4) is valid in this time interval, and holds for the class of all considered inputs to the system. Note that T₁ = ∞ in the case that the systems remains in QSS, which may for example not be the case for models with switches (where, e.g., the values of a subset of the state variables may change when a certain condition is fulfilled) [34, 35].

Note that Eq. (3) (steady state) necessitates that Eq. (4) (quasi-steady state) is fulfilled, but not vice versa. Although QSS implies that (some of) the terms of the right-hand side of the ODE are large and leaves the left-hand side (derivative term) negligible, the derivative term may still be large enough for the state variable in QSS to change considerably during the time-span of a simulation; the key is that these changes mainly occur on a slow manifold.

Zooming of Linear Models

where x _r ∈ ℝ^{(n - w + 1)}, A _r ∈ ℝ^{(n - w + 1) × (n - w + 1)}, B _r ∈ ℝ^{(n - w + 1) × m}, and C _r ∈ ℝ^{l × (n - w + 1)}.

where $k_{jL}^{\prime }$ is the rate parameter in the reaction from state variable $x_{lL}^{\prime }$ to state variable x _j in the reduced model, and k _ji is the rate parameter in the reaction from state variable x _i to state variable x _j in the original model.

Finally, note that Eqs. (5) and (6) provide a link between ${ℳ}_o$ and ${ℳ}_r$, which constitute two different levels of granularity. It is this link between the models that make us consider them as two different degrees of zooming, and the primary goal of this paper is to establish such a link also for nonlinear models.

Extension to Nonlinear Models - Initial Observations

We will now present some key observations that are used in the derivation of the method for zooming of nonlinear models.

where $M\in {\mathbb{N}}^{n-n_r\times n}$.

where Eqs. (9) and (10) were used.

since the term S _f r _f (x _f , p) dominates the term S _s r _s (x, p). Note that Eq. (4) and consequently Eq. (13) only hold in T₀ ≤ t < T₁, since the state is only known to be fast in this time span.

Eq. (12) defines the ODEs of the reduced model (the lumped state variables), and Eqs. (11) and (13) can in principle be used to calculate back-translation formulae, as is demonstrated with the small example model in the next section. However, as we shall see, this approach requires the explicit algebraic solution to a system of nonlinear equations, which is typically an infeasible task. Furthermore, there is not a clear one-to-one mapping between the state variables of the original and reduced models as in the case of proper lumping [27].

A Small Example Model

We will now present a small example model, with three fast state variables, which is reduced with the approach discussed above. An alternative approach is then demonstrated with the advantage that it scales better to larger models.

where the lumped state variables L₁ and L₂ are introduced. Note that Eq. (12) defines the dynamics of the lumped state variables.

where $K_1=\frac{k_{-1}}{k_1}$.

We can now employ Eq. (12) to solve the ODEs for the lumped state variables L₁ and L₂, and use Eqs. (19)-(21) as back-translation formulae to compute the trajectories of the original state variables A, B, and C. However, note that for even slightly larger clusters of fast species than the one discussed here it would not be possible to calculate algebraic expressions of the original state variables with this approach, since it builds on the explicit solution of a system of nonlinear equations, which quickly becomes infeasible with growing problem size.

The reduced model consists of the two state variables L₁ and B (note that L₂ does not appear in the reduced model), and the dynamics is described by Eqs. (12) and (25), respectively. Note that the state variables A and C can be back-translated from the reduced model with Eqs. (23) and (24). This approach is a bit more intricate than the first, but comes with the advantage that we do not need to solve a system of nonlinear equations.

A Method for Zooming of Nonlinear Models

We will now step-by-step present a method that can be used to construct zoomable nonlinear biochemical models. This involves two sub-goals: i) to identify a reduced model that shares important characteristics with the original model, ii) to derive back-translation formulae that can be used to compute the original state variables and parameters from the reduced model.

In an initialization step of the method for a model ${ℳ}_o$ we first formulate mathematical equations for all conservation relations Eq. (8), state variables in steady-state Eq. (3), and quasi-steady state assumptions Eq. (4). If additional properties of the system are known, we also formulate the corresponding equations.

Step 1

The first step of the method is to identify the apparent conservation relations in the model.

Definition 1: Let S _f be the stoichiometric matrix for the reactions r _f (x _f , p) as defined in Eq. (9). Each subset of state variables for which the corresponding rows of S _f are linearly dependent constitutes an apparent conservation relation. Hence the apparent conservation relations lie in the left null space of S _f and the dimension of this space is n - rank(S _f ).

Note that the apparent conservation relations are defined in Eq. (11). It is trivial to identify the set of all linearly dependent rows of S _f with a mathematical computing software (e.g., SBtoolbox for Matlab [36]).

Step 2

The second step of the method is to define the state variables of the reduced model, which we refer to as modified lumped state variables.

Definition 2: Let x be a lumped state variable corresponding to a subset of the state variables in an apparent conservation relation. Then x is a modified lumped state variable if the lumping scheme with respect to the state variables of the original model is proper.

where M _m is a × n _f matrix with elements equal to 0 or 1 and column sums equal to 1, and l _m denotes the modified lumped state variables. We typically have a large freedom in the choice of M _m . The number of state variables is maximally reduced if all exact conservation relations in the model are retained as modified lumped state variables (and replaced by constants).

Step 3

where we used Eq. (27) in the last step.

A modified lumped state variable for which an insufficient number of linear and linearly independent equations are available may still be used in the reduced model. However, the back-translation of the modified lumped state variable to the original state variables is then not possible, and step 3 of the method is ignored.

Step 4

where the matrix M _m is defined in Eq. (26), and $g_i\left(l_m,p\right)\equiv x_{f_i}$ is introduced to simplify the notation.

Proof:

where Eq. (12) was used in the last step. □

The matrix I + J(l _m ; p) is symbolically invertible, but may in general contain singularities for particular combinations of parameters values and state variable values. However, the matrix is always invertible for the models discussed in this paper, since the corresponding determinants are strictly positive.

Step 5

The final step of the method is to back-translate the modified lumped state variables to the original state variables with the fraction parameters derived in step 3. This allows a comparison between the predictions by the reduced model to those of the original model.

The implementation of the method is straight-forward, and we have used Matlab (R2008b) together with the SBtoolbox [36] as computing software for the models in this paper.

Enzyme Kinetics Model

Note that the complexes C _s and C _p are formed by S bound to E, and P bound to E, respectively. This model is frequently occurring as part of larger models of biological systems, although the reaction from C _S to C _P is sometimes neglected, or reversible. The ODEs for the model are listed in Appendix A.1, where the three reactions are defined as: r₁ = k₁SE - k_-1C _S , r₂ = k₂C _S , and r₃ = k₃C _P - k_-3PE.

which is constant since the total amount of the enzyme E is conserved in the system.

Reduction of the Enzyme Kinetics Model

The first step of the method is to identify the apparent conservation relations from the matrix S _f . Since r₂(x, p) is dominated by r₁(x, p) and r₃(x, p) the model ODEs can be written on the form of Eq. (9)

$\begin{array}{ll}\hfill \left(\begin{array}{c}\hfill Ṡ\hfill \\ \hfill Ė\hfill \\ \hfill Ṗ\hfill \\ \hfill \stackrel{°}{C_S}\hfill \\ \hfill \stackrel{°}{C_P}\hfill \end{array}\right)& =S_s{r}_s\left(x,p\right)+S_f{r}_f\left(x_f,p\right)=\\ =\left(\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \\ \hfill -1\hfill \\ \hfill 1\hfill \end{array}\right)r_2+\left(\begin{array}{cc}\hfill -1\hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill \end{array}\right)\left(\begin{array}{c}\hfill r_1\hfill \\ \hfill r_3\hfill \end{array}\right),\end{array}$

A basis of the left null space of S _f is given by the row vectors of M _f , which is defined by

$\begin{array}{lll}\hfill l& =\left(\begin{array}{c}\hfill L_S\hfill \\ \hfill L_P\hfill \\ \hfill L_E\hfill \end{array}\right)=& \hfill \text{(1)}\\ =\left(\begin{array}{ccccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill \end{array}\right)\left(\begin{array}{c}\hfill S\hfill \\ \hfill E\hfill \\ \hfill P\hfill \\ \hfill C_S\hfill \\ \hfill C_P\hfill \end{array}\right)=& \hfill \text{(2)}\\ =M_f{x}_f,& \hfill \text{(3)}\end{array}$

(37)

where L _S and L _P are apparent conservation relations and LE is an exact conservation relation.

The second step is to define the modified lumped state variables on the form of Eq. (26). The number of state variables is maximally reduced if L _E is retained as a state variable in the reduced model (i.e., since L _E , unlike L _S and L _P , can be replaced by a constant). The vector of modified lumped state variables then is defined

$\begin{array}{ll}\hfill l_m& =\left(\begin{array}{c}\hfill S\hfill \\ \hfill P\hfill \\ \hfill L_E\hfill \end{array}\right)=\left(\begin{array}{ccccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill \end{array}\right)\left(\begin{array}{c}\hfill S\hfill \\ \hfill E\hfill \\ \hfill P\hfill \\ \hfill C_S\hfill \\ \hfill C_P\hfill \end{array}\right)=\\ =M_f{x}_f,\end{array}$

In the third step of the method we calculate fraction parameters for the modified lumped state variable L _E . There are five equations (Eqs. (34)-(35) and (37)) that are linear w.r.t. the state variables E, C _S , and C _P , which are lumped into the state variable L _E . Note that only n _m = 3 equations are required to derive fraction parameters, and we use Eqs. (34)-(36) to formulate an equation system as in Eq. (28) with the solution

$E={\eta }_E{L}_E\approx \frac{1}{1+M_{1}S+M_{3}P}{L}_E\equiv g_2\left(l_m,p\right),$

(38)

$C_S={\eta }_{C_S}{L}_E\approx \frac{M_{1}S}{1+M_{1}S+M_{3}P}{L}_E\equiv g_4\left(l_m,p\right),$

(39)

$C_P={\eta }_{C_P}{L}_E\approx \frac{M_{3}P}{1+M_{1}S+M_{3}P}{L}_E\equiv g_5\left(l_m,p\right),$

(40)

where $M_1=\frac{k_1}{k-1}$ and $M_3=\frac{k_3}{k-3}$. The two remaining modified lumped state variables correspond to S and P in the original model, so we define that $l_{m_1}=S\equiv g_1\left(l_m,p\right)$ and $l_{m_2}=P\equiv g_3\left(l_m,p\right)$.

In the fourth step we derive the rate of change of the modified lumped state variables. Eq. (32) gives that ${\stackrel{°}{L}}_E=0$, which is replaced by a constant, and

$Ṡ\approx \frac{-k_2{M}_{1}S{L}_E\left({\left(1+M_{1}S+M_{3}P\right)}^2+M_3{L}_E\right)}{\varphi \left(S,P,L_E,M_1,M_3\right)},$

(41)

$Ṗ\approx \frac{k_2{M}_{1}S{L}_E\left({\left(1+M_{1}S+M_{3}P\right)}^2+M_1{L}_E\right)}{\varphi \left(S,P,L_E,M_1,M_3\right)}.$

(42)

where

$\varphi \left(S,P,L_E,M_1,M_3\right)=\left({\left(1+M_{1}S+M_{3}P\right)}^3+\left(M_1+M_3+M_1{M}_3\left(P+S\right)\right)\left(1+M_{1}S+M_{3}P\right)L_E+M_1{M}_3{L}_E^2\right).$

The two ODEs in Eqs (41)-(42) define the dynamics of the state variables in the reduced model. We finally note that the exact conservation relation for the substrate, L _T = S + P + C _S + C _P , together with Eqs. (39)-(40) can be used to reduce the model further to a single state.

In the fifth step of the method we use the fraction parameters, defined in Eqs. (38)-(40), to back-translate the modified lumped state variables to the state variables of the original model. A comparison between predictions of the original state variables, from simulations of the original model and the reduced model, is presented in Figure 1. Implementations of the original model (Additional file 1), the reduced model (Additional file 2), and a script for simulation with SBtoolbox2 for MATLAB [36] (Additional file 3), are available in Additional files.

The only assumption that was used in the derivation of the reduced model is that the reaction terms r₁(x, p) and r₃(x, p) dominate the reaction term r₂(x, p), which results in that all state variables are in QSS. To assess the impact of these assumptions on the reduced model we compute the relative difference between the state variables in the original and in the reduced model

${\epsilon }_i\left(t\right)=\frac{|x_i^o\left(t\right)-x_i^r\left(t\right)|}{x_i^o\left(t\right)},i=1,\dots ,n,$

where $x_i^o$ is state variable i in the original model, $x_i^r$ is the corresponding back-translated state variable in the reduced model, and |x| denotes the elementwise absolute values of x. The maximal mean and infinity norm of ε _i (t) in Eq. (43) over time is presented in Table 1 for parameter values over five orders of magnitude. In general, the reduced model appears to be robust to changes in the parameter values, although slightly more sensitive to some parameters (e.g., small values of k_-1, large values of k₁, or large values of k₂, which violate the assumptions used in the reduction). However, note that the validity of the QSS assumption may also depend on the state variables, for example the total concentration of the enzyme. Interestingly, we observed that the reduced model can well approximate the original model over several order of magnitudes around the nominal enzyme concentration (L _E = 1). It is well-known that the QSS approximation is only valid for sufficiently small enzyme concentrations, and as expected the performance of the reduced model starts to decrease for immense enzyme concentrations.

Param./Factor	10-2	10-1	100	101	102
k 1	0.00053/0.0082	0.0018/0.0071	0.0010/0.0059	0.019/0.18	0.19/1.8
k -1	0.19/1.9	0.0061/0.079	0.0010/0.0059	0.0017/0.0025	0.00024/0.0028
K 2	0.00019/0.0060	0.00082/0.0055	0.0010/0.0059	0.035/0.23	0.21/0.39
k 3	0.0068/0.014	0.0053/0.010	0.0010/0.0059	0.0067/0.044	0.035/0.29
K -3	0.020/0.30	0.0048/0.032	0.0010/0.0059	0.0051/0.0093	0.0069/0.015
All	0.012/1.2	0.0014/0.062	0.0010/0.0059	0.025/0.41	0.010/0.045

The nominal parameter values (k₁ = 1000, k_-1 = 2000, k₂ = 1, k₃ = 1000, and k_-3 = 3000) are modified by a multiplicative factor, and the maximal (for any state variable) time average/infinity norm of the relative difference between the original and the reduced model is presented above. Note that only concentrations larger than 10^-6 are considered in the analysis above, due to potential numerical inaccuracies.

Note that all the state variables of the original model have a direct biological interpretation also in the reduced model, and that Eqs. (38)-(40) can be used to back-translate the state variables. The reduced model may be depicted

$S+{\eta }_E{L}_E\underset{k_1}{\stackrel{k_{-1}}{\leftrightarrow }}{\eta }_{C_S}{L}_E\underset{k_2}{\overset{}{\to }}{\eta }_{C_P}{L}_E\underset{k_3}{\stackrel{k_{-3}}{\leftrightarrow }}P+{\eta }_E{L}_E$

where the fraction parameters specify the distribution of the enzyme among the corresponding original state variables.

Glucose Transport in Budding Yeast

A model for the transport of glucose into a cell of baker's yeast (S. cerevisiae), which constitutes the first step of glycolysis, is presented in [37]. The inflow of glucose is modeled as a facilitated diffusion process, in which a carrier enzyme is responsible for the transport between the inner and outer regions of the cellular membrane. It is assumed that glucose 6-phosphate (G6P) has an inhibitory role in the glucose transport process by binding to the transporter. A graphical representation of the model is shown in Figure 2, and the ODEs for the state variables are listed in Appendix A.2.

In [27] we described how the calculation of fraction parameters, based on a set of assumptions, leads to the same reaction rates in the reduced model as were reported in [37]. The assumptions are that state variables participating in reactions for uptake and release of glucose and G6P across the cell membrane are in QSS, that the transporter is conserved, and that the concentrations of the transporter in the inner and outer regions of the cellular membrane are constant.

It is not clear how equations for back-translation of the state variables in the reduced model in [27, 37] can be derived. The reduced model has three state variables; external- and internal glucose and G6P, but only two differential equations for the in- and outflow of glucose, since the ODE for G6P is replaced by a representative function that is inferred from the G6P data. Our method does not rely on that such information is available, although it would in principle be possible to utilize data fitted functions for the state variables. Also note that the equations in the reduced model describe the total influx and efflux of glucose across the membrane [27], which cannot be interpreted w.r.t. the state variables of the original model. Other assumptions in [37] that complicates a comparison with our method are that the efflux of glucose is negligible, that the concentration of glucose in the cytosol is negligible, and that the concentrations of the transporter are constant in the inner and outer regions of the cell membrane (Eq. (50)).

It is not possible to generate a reduced model by direct substitution of the fraction parameters that were derived in [27] into the ODEs of the original model, since this would lead to the prediction that the state variables are constant (as discussed in [27]). We will now instead illustrate how our method can be used to derive a reduced and zoomable version of the glucose transport model.

Reduction of the Glucose Transport Model

Before applying our method to the glucose transport model we tried an alternative approach. Eqs. (43)-(46) were solved w.r.t. the state variables in QSS, and the resulting expressions were then substituted into the remaining ODEs. The details of the derivation of the reduced model are presented in Appendix A.3. The reduced model does not produce satisfactory predictions for any other state variable than $x_{G6P}^i$, which remains approximately constant during the simulation. Implementations of the original model (Additional file 4), the reduced model (Additional file 5), and a script for simulation with SBtoolbox2 for MATLAB [36] (Additional file 6), are available in Additional files.

Since the first approach turned out to be insufficient for reduction of the glucose transport model we applied our method to the same model. Following [37], we initially assumed constant transporter concentrations in the inner and outer regions of the cellular membrane, as defined by Eq. (50). For the details on the derivation of the reduced model we refer to Appendix A.4. The reduced model clearly performs better than the model resulting from the first approach, but it is still not satisfactory. However, the assumption of constant regional concentrations of the transporter may not be valid since the transport of glucose across the cell membrane is a rate limiting step in the model, and appears to be important for the state variable dynamics. We therefore decided to neglect Eq. (50) in the reduction process.

Implementations of the original model (Additional file 4), the reduced model (Additional file 7), and a script for simulation with SBtoolbox2 for MATLAB [36] (Additional file 8), are available in Additional files. In the first step of the method we identify the following apparent conservation relations

$\begin{array}{lll}\hfill l=& \left(\begin{array}{c}\hfill L_{Gl{c}_1}\hfill \\ \hfill L_{E_1}\hfill \\ \hfill L_{Gl{c}_2}\hfill \\ \hfill L_{G6P}\hfill \\ \hfill L_{E_2}\hfill \end{array}\right)=Mx=& \hfill \text{(1)}\\ \hfill =& \left(\begin{array}{ccccccccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right)x.\end{array}$

(51)

We note that there are two disjoint clusters of fast reactions in the model, corresponding to the outer- and inner parts of the cell membrane.

In the second step we define the modified lumped state variables. We decide to keep the lumped state variables $L_{E_1}$ and L_{G 6P}as modified lumped state variables. The choice to keep L_{G 6P}leads to the largest possible reduction in the number of state variables, since the conservation of L_{G 6P}is exact (which is not true for any other state variable in l). The modified lumped state variables are defined true for any other state variable in l). The modified lumped state variables are defined

$\begin{array}{lll}\hfill l_m=& \left(\begin{array}{c}\hfill x_{Glc}^e\hfill \\ \hfill L_{E_1}\hfill \\ \hfill x_{Glc}^i\hfill \\ \hfill L_{G6P}\hfill \\ \hfill L_{E_3}\hfill \end{array}\right)=M_{m}x=& \hfill \text{(1)}\\ \hfill =& \left(\begin{array}{ccccccccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1\hfill \end{array}\right)x,\end{array}$

(52)

where

$x=\left(\begin{array}{c}\hfill x_{Glc}^e\hfill \\ \hfill x_{Glc}^i\hfill \\ \hfill x_{E-G6P}^i\hfill \\ \hfill x_{E-Glc-G6P}^i\hfill \\ \hfill x_{G6P}^i\hfill \\ \hfill x_{E-Glc}^e\hfill \\ \hfill x_{E-Glc}^i\hfill \\ \hfill x_E^e\hfill \\ \hfill x_E^i\hfill \end{array}\right).$

Note that two of the state variables in the original model, $x_{Glc}^e$ and $x_{Glc}^i$, are also modified lumped state variables.

In the third step of the method we calculate fraction parameters for the modified lumped state variables that correspond to more than one of the original state variables (i.e., $L_{E_1}$, L_{G 6P}, and $L_{E_3}$). All of the modified lumped state variables satisfy the requirement that at least n _m of Eqs. (43)-(49) and Eq. (52) are linear, and linearly independent, with respect to the corresponding original state variables. Eqs. (43) and (52) form a nonlinear equation system with the solution

$\begin{array}{lll}\hfill \left(\begin{array}{c}\hfill x_E^e\hfill \\ \hfill x_{E-Glc}^e\hfill \end{array}\right)=& \left(\begin{array}{c}\hfill {\eta }_E^e\hfill \\ \hfill {\eta }_{E-Glc}^e\hfill \end{array}\right)L_{E_1}\approx & \hfill \text{(1)}\\ \hfill \approx & \frac{1}{\left(K_1+x_{Glc}^e\right)}\left(\begin{array}{c}\hfill K_1\hfill \\ \hfill x_{Glc}^e\hfill \end{array}\right)L_{E_1}.\end{array}$

(53)

Let us define that $g_6\equiv x_{E-Glc}^e$ and $g_8\equiv x_E^e$. Similarly, the fractions of the two carrier state variables in the inner regions of the cell to the lumped state variable $L_{E_3}$ can be computed from Eqs. (44) and (52) (Eq. (28))

$\begin{array}{lll}\hfill \left(\begin{array}{c}\hfill x_{E-Glc}^i\hfill \\ \hfill x_E^i\hfill \end{array}\right)=& \left(\begin{array}{c}\hfill {\eta }_{E-Glc}^i\hfill \\ \hfill {\eta }_E^i\hfill \end{array}\right)L_{E_3}\approx & \hfill \text{(1)}\\ \approx \frac{1}{\left(x_{Glc}^i+K_2\right)}\left(\begin{array}{c}\hfill x_{Glc}^i\hfill \\ \hfill K_2\hfill \end{array}\right)L_{E_3}.\end{array}$

(54)

We define that $g_7\equiv x_{E-Glc}^i$ and $g_9\equiv x_E^i$. The fraction parameters for the G6P-state variables can be computed from Eqs. (45)-(46) and Eq. (28) with the solution

$\begin{array}{ll}\left(\begin{array}{c}\hfill x_{E-G6P}^i\hfill \\ \hfill x_{E-Glc-G6P}^i\hfill \\ \hfill x_{G6P}^i\hfill \end{array}\right)=\left(\begin{array}{c}\hfill {\eta }_{E-G6P}^i\hfill \\ \hfill {\eta }_{E-Glc-G6P}^i\hfill \\ \hfill {\eta }_{G6P}^i\hfill \end{array}\right)L_{G6P}\approx & \hfill \text{(1)}\\ \approx \frac{1}{\xi }\left(\begin{array}{c}\hfill K_3{x}_E^i\hfill \\ \hfill K_4{x}_{E-Glc}^i\hfill \\ \hfill K_3{K}_4\hfill \end{array}\right)L_{G6P}.\end{array}$

(55)

where $\xi =K_3{K}_4+K_4{x}_{E-Glc}^i+K_3{x}_E^i$. We define that $g_3\equiv x_{E-G6P}^i$, $g_4\equiv x_{E-Glc-G6P}^i$, and $g_5\equiv x_{G6P}^i$. For the two original state variables that are kept as modified lumped state variables in the reduced model we define that $g_1\equiv x_{Glc}^e$ and $g_2\equiv x_{Glc}^i$.

The fourth step of the method is to derive rate equations for the modified lumped state variables. Since the apparent conservations are separated into two disjoint clusters of fast state variables, we can treat the model for the inner and outer regions of the membrane separately. Let the modified lumped state variables corresponding to the outer region of the cell membrane be denoted by $l_{m_1}={\left(x_{Glc}^e{L}_{E_1}\right)}^T$, and the variables in the inner region of the cell membrane by $l_{m_2}={\left(x_{Glc}^i{L}_{E_3}\right)}^T$. Note that the state variable L_{G 6P}can be replaced by a constant in the model, since ${\stackrel{°}{L}}_{G6P}=0$. The ODEs of the modified lumped state variables are derived with Eq. (31). In the inner region the ODEs are

$\begin{array}{lll}\hfill {\stackrel{°}{l}}_{m_1}=& \left(\begin{array}{c}\hfill {ẋ}_{Glc}^e\hfill \\ \hfill {\stackrel{°}{L}}_{E_1}\hfill \end{array}\right)={\left(I+J\right)}^{-1}{\stackrel{°}{l}}_1=& \hfill \text{(1)}\\ ={\left(\begin{array}{cc}\hfill 1+J_{11}\hfill & \hfill J_{12}\hfill \\ \hfill 0\hfill & \hfill 1\hfill \end{array}\right)}^{-1}{\stackrel{°}{l}}_1,\end{array}$

(56)

where

$\begin{array}{lll}\hfill {\stackrel{°}{l}}_1=& \left(\begin{array}{c}\hfill {\stackrel{°}{L}}_{Gl{c}_1}\hfill \\ \hfill {\stackrel{°}{L}}_{E_1}\hfill \end{array}\right)\approx & \hfill \text{(1)}\\ \approx \left(\begin{array}{c}\hfill -\alpha \left(x_{EGlc}^e-x_{EGlc}^i\right)\hfill \\ \hfill -\alpha \left(x_{EGlc}^e-x_{EGlc}^i\right)-\beta \left(x_E^e-x_E^i\right)\hfill \end{array}\right).\end{array}$

(57)

In the larger outer region of the cell membrane the ODEs take the form

$\begin{array}{lll}\hfill {\stackrel{°}{l}}_{m_2}\approx & \left(\begin{array}{c}\hfill {ẋ}_{Glc}^i\hfill \\ \hfill {\stackrel{°}{L}}_{E_3}\hfill \end{array}\right)={\left(I+J\right)}^{-1}{\stackrel{°}{l}}_2=& \hfill \text{(1)}\\ \hfill =& {\left(\begin{array}{cc}\hfill 1+J_{11}\hfill & \hfill J_{12}\hfill \\ \hfill J_{21}\hfill & \hfill 1+J_{22}\hfill \end{array}\right)}^{-1}{\stackrel{°}{l}}_2,\end{array}$

(58)

where

$\begin{array}{lll}\hfill {\stackrel{°}{l}}_2=& \left(\begin{array}{c}\hfill {\stackrel{°}{L}}_{Gl{c}_2}\hfill \\ \hfill {\stackrel{°}{L}}_{E_2}\hfill \end{array}\right)\approx & \hfill \text{(1)}\\ \hfill \approx & \left(\begin{array}{c}\hfill -\alpha \left(x_{EGlc}^i-x_{EGlc}^e\right)\hfill \\ \hfill -\alpha \left(x_{EGlc}^i-x_{EGlc}^e\right)-\beta \left(x_E^i-x_E^e\right)\hfill \end{array}\right).\end{array}$

(59)

In the fifth step of the method the reduced model, which is defined by the four ODEs in Eqs. (56)-(59), is simulated. The trajectories of the state variables of the reduced model can be back-translated to the original state variables with the fraction parameters defined in Eqs. (53), (54), and (55). The simulation results are shown in Figure 3 and Figure 4. All the state variables can be back-translated properly, which shows that the model properties that are important for recovery of the state variables are retained in the reduction. Implementations of the original model (Additional file 4), the reduced model (Additional file 9), and a script for simulation with SBtoolbox2 for MATLAB [36] (Additional file 10), are available in Additional files.

If we use $L_{E_2}$ instead of L_{G 6P}as a modified lumped state variable, the reduced model will have the same state variables as the reduced model in [27, 37] (i.e., $x_{Glc}^e,x_{Glc}^i$, and $x_{G6P}^i)$ two additional state variables for the transporter. This gives a reduced model with five state variables, but equally many parameters as in the previous case. A comparison between the original model and the reduced model, w.r.t. the original state variables, is shown in Figure 5 and Figure 6. As can be seen the comparison is very good, in fact it is even slightly better than for the reduced model with four state variables. The details of the derivation of the reduced model are presented in Appendix A.5. Implementations of the original model (Additional file 4), the reduced model (Additional file 11), and a script for simulation with SBtoolbox2 for MATLAB [36] (Additional file 12), are available in Additional files.

Note that the only assumption used to derive the reduced model is that states that are involved in reactions at the membrane are in QSS. To investigate the parameter space region in which the QSS assumptions are valid we use the measure defined in Eq. (43). The maximal mean and infinity norm of the relative difference between the original and the reduced model in Eq. (43) over time is presented in Table 2. The reduced model appears to be relatively robust to changes in the parameters, although sensitive to small values of k_-4 and to large values of k₄. This is mainly due to that a large proportion of the transporter E is absorbed in $x_{E-G6P}^i$, which leads to that some of the QSS assumptions are invalid. We also observed that relative difference between the models is insensitive to the total concentration of the transporter for several orders of magnitude around the nominal value (L _E = 0.01). However, note that this observation is specific to the studied model and may not be generalizable to other similar biochemical models.

Param./Factor	10-2	10-1	100	101	102
k 1	0.079/0.091	0.079/0.10	0.081/0.091	0.083/0.10	0.084/0.11
k -1	0.085/0.23	0.084/0.093	0.081/0.091	0.078/0.10	0.078/0.089
k 2	0.16/0.93	0.080/0.55	0.081/0.091	0.063/0.083	0.054/0.062
k -2	0.054/0.10	0.062/0.074	0.081/0.091	0.097/0.10	0.10/0.10
k 3	0.081/0.091	0.082/0.091	0.081/0.091	0.079/0.089	0.075/0.082
k -3	0.070/0.080	0.079/0.089	0.081/0.091	0.082/0.091	0.081/0.091
k 4	0.23/0.32	0.21/0.30	0.081/0.091	6.36/6.93	310/336
k -4	310/336	6.36/6.93	0.081/0.091	0.21/0.30	0.23/0.32
α	0.10/0.10	0.095/0.17	0.081/0.091	0.080/0.089	0.083/0.088
β	0.076/0.088	0.091/0.096	0.081/0.091	0.075/0.15	0.074/0.19
All	0.16/0.97	0.081/0.69	0.081/0.091	0.060/0.091	0.054/0.081

The nominal parameter values (k₁ = 1000, k_-1 = 1100, k₂ = 1000, k_-2 = 1200, k₃ = 1000, k_-3 = 7000, k₄ = 1000, k_-4 = 1100, α = 4.2, and β = 1) are modified by a multiplicative factor, and the maximal (for any state variable) time average/infinity norm of the relative difference between the original and the reduced model is presented. Note that due to potential numerical inaccuracies only concentrations larger than 10^-6 are considered in the analysis.