Biochemical Reaction Systems
Most biological processes of interest to systems biology are modeled by means of open biochemical reaction systems; i.e., systems that exchange mass with their surroundings [14]. Living cells, for example, are open systems maintaining themselves by exchanging materials with their environment. Mass exchange is modeled either explicitly or implicitly. Examples of explicit modeling include: clamped species, reactions with null species as reactants or products, and irreversible reactions [15]. Clamped species are chemicals whose concentrations are held fixed. They are often used to model molecular species whose concentrations are affected by unknown reactions. It is apparent that these chemicals must be supplied or removed from the system at appropriate rates to ensure that their concentrations do not deviate from their fixed values. On the other hand, reactions with null reactants or products model mass transfer in and out of the system, respectively. Finally, the use of an irreversible reaction is based on the assumption that the concentration of at least one of its products is clamped to zero (otherwise, the reaction can be reversed at sufficiently high product concentrations), which implies mass transfer as well. An example of implicitly modeling mass exchange is a reaction with reactants or products not modeled by the system. A common case would be the phosphorylation of a protein without explicitly modeling conversion of ATP to ADP [5].
An open biochemical reaction system is comprised of N molecular species X_{1}, X_{2},..., X_{
N
} that interact through M coupled reactions, given by
where , ℳ: = {1, 2,..., M}, and ν_{
nm
}, are the stoichiometries of the reactants and products. We can characterize an open biochemical reaction system at time t ≥ 0 by the concentrations of all molecular species at t. Clamped molecular species have concentrations that do not vary with time, whereas, the concentrations of the remaining 'dynamic' species evolve as a function of time. We will assume that the system characterizes reactions in an ideal and wellstirred (homogeneous) mixture at constant temperature and volume and that the concentrations vary continuously in time. In this case, we can describe the dynamic evolution of the molecular concentrations in the system by the following chemical kinetic equations:
initialized by setting x_{
n
}(0) = q_{
n
}, , for some initial concentrations q_{
n
}, , where ϕ m(t, k) is the net flux of the m^{th} reaction, is the net stoichiometry coefficient of the n^{th} molecular species associated with the m^{th} reaction, is an observation time window of interest, and k is a vector of kinetic parameters , where . The parameters k characterize the biochemical reaction system at hand and are independent of the molecular concentrations . Moreover, we assume that these parameters do not vary with time.
By appropriately pruning and modifying an open biochemical reaction system, we can derive a closed subsystem (i.e., a system that does not exchange mass with its surroundings) that lends itself to thermodynamic analysis since external or unknown thermodynamic forces no longer exist. To do so, we first remove all null and irreversible reactions, as well as partially modeled reactions (i.e., reactions with incomplete stoichiometries). Next, we remove clamping of molecular species involved in reversible reactions. Subsequently, we keep only reactions that are thermodynamically independent. Thermodynamic independence among reactions means that a reaction is only driven by its own thermodynamic force, which implies that the affinity of the reaction will be zero if and only if there is no change in its degree of advancement. This condition is usually fulfilled if we keep only elementary reactions (i.e., reactions that take place in one single step). As a consequence, we obtain a closed reaction set ℳ_{0} ⊆ ℳ. Finally, we remove all species that are no longer involved with any of the reactions in ℳ_{0}, leaving only the molecular species associated with the closed subsystem.
The main rationale behind the second step is that the kinetic parameters k considered in this paper are assumed to be independent of the molecular concentrations. As a consequence, the values of these parameters will not change if the concentrations of the clamped species are allowed to vary. Therefore, we can construct a (possibly artificial) situation in which the concentrations of the clamped species vary as if they were dynamic species. Because our goal in creating the closed subsystem is to discover and enforce thermodynamic constraints on the kinetic parameters, we must include the clamped species in our model. This is contrary to simply removing all clamped species and the associated reactions, since this approach will not allow us to determine thermodynamically feasible values for the kinetic parameters of the removed reactions.
The third step is due to a simplification imposed on us by the current state of nonequilibrium thermodynamics. Thermodynamically dependent reactions influence each other, since the thermodynamic force of one reaction may drive the other reaction and vice versa. Unfortunately, it is not clear at this point how to deal with thermodynamically coupled reactions. Future research may be necessary to address this issue.
The resulting closed biochemical reaction subsystem is comprised of N_{0} molecular species that interact through M_{0} coupled reversible reactions. The dynamic evolution of the molecular concentrations in this system is governed by:
initialized by x_{
n
} (0) = q_{
n
} , for . We will characterize all reactions in ℳ_{0} by the generalized massaction rate law [16]. In this case, the net fluxes are given by
for m ∈ ℳ_{0}, where r_{2m1}, r_{2m}are the (generalized) rate constants of the m^{th} forward and reverse reactions, respectively. The quantity f_{
m
} [x(t), π ] can be any positive and finite function of the concentrations x(t) and may depend on a set of kinetic parameters π. For usual mass action kinetics, f_{
m
} [x(t), π ] = 1. However, for more complex schemes, this function usually takes a rational or polynomial form. It is known that all reversible reaction rate laws in ideal mixtures (including reversible MichaelisMenten kinetics) can be described by (3) [5]. Note that k includes the kinetic parameters π as well as the rate constants {r_{2m1}, r_{2m}, m ∈ ℳ_{0}}. It also includes the kinetic parameters of all reactions in ℳ \ ℳ_{0}.
It is a direct consequence of thermodynamic analysis that a closed biochemical reaction system will asymptotically reach a unique nonzero state of chemical equilibrium at which all concentrations become stationary (assuming, of course, we have nonzero initial conditions), satisfying the following detailed balance equations:
As a consequence,
These constraints must be satisfied by the rate constants in order for the closed biochemical reaction system to be thermodynamically feasible.
The constraints implied by (5) correspond to the reaction 'cycles' in the system. A reaction cycle is comprised of those reactions corresponding to the nonzero elements of a vector in the null space null () of the stoichiometry matrix of the closed system. Clearly, (2) will be at a fixed point whenever the net fluxes of the underlying reactions are set equal to the corresponding elements of a vector in null (). If we denote by s the N_{0}× 1 vector whose n^{th} element is the log steadystate concentration ln and by z the M_{0}× 1 vector whose m^{th} element is the log equilibrium constant z_{
m
} = ln(r_{2m1}/r_{2m}), then we can write (5) in a matrixvector form as . Now, if b is a vector in the null space of , then and , or
which are the wellknown Wegscheider conditions [6]. These conditions express necessary and sufficient constraints on the reaction rate constants of a closed biochemical reaction system to be thermodynamically feasible. Thus, if we denote the set of all thermodynamically feasible parameters k by , then any satisfies (6) and, likewise, any k that satisfies (6) is a member of .
We want to emphasize that, in open biochemical reaction systems, the rate constants of the reversible reactions must also be constrained by the Wegscheider conditions, even if the system is far from equilibrium. To identify these constraints, we need to prune an open biochemical reaction system into a closed subsystem, by employing the technique discussed previously, and use the resulting stoichiometry matrix to calculate the constraints given by (6).
An equally important observation is that the rate constants of the reactions pruned from an open system are not constrained by the Wegscheider conditions, since (4) must only be satisfied by the reactions in the closed subsystem. Furthermore, if a reaction m in a closed system is not part of a cycle, then b_{
m
} = 0, for every b∈ null (), and its forward and reverse rate constants will not be thermodynamically constrained, since these constants trivially satisfy (6).
It can be shown (see Additional file 1) that the entropy production rate of an open biochemical reaction system at chemical equilibrium in which the net fluxes of the reactions in ℳ_{0} equal to the elements of a vector in the null space of the stoichiometry matrix , is given by
where A = 6.02214084(18) × 10^{23}mol^{1} is the Avogadro number, V is the system volume, and k_{
B
} = 1.3806504 × 10^{23}JK^{1} is the Boltzmann constant. According to the second law of thermodynamics, the entropy production rate must always be greater than or equal to zero, with equality if and only if the system is at thermodynamic equilibrium. It is therefore clear from (7) that the Wegscheider conditions imply that the entropy production rate must be zero in this case (i.e., the system must be at thermodynamic equilibrium). As a consequence, the chemical motive force (which is the amount of energy added to the system per unit time due to mass exchange through its boundary) and the heat dissipation rate must also be zero. This makes intuitive sense, since a reaction cycle leaves all molecular concentrations unchanged and, therefore, there is no change in internal energy or mass flow through the system boundary. Clearly, we can think of the Wegscheider conditions as being a direct consequence of the thermodynamic requirement that σ (b) = 0, for every b∈ null ().
Linear constraints
Unfortunately, (6) imposes a possibly infinite number of nonlinear constraints on the rate constants of a closed biochemical reaction system. However, it is sufficient to satisfy (6) for M_{2} = M_{0}  M_{1} basis vectors {b^{(i)}, i = 1, 2,..., M_{2}} of the null space of , where M_{0} is the number of reactions and M_{1} = rank () (see Additional file 1). By using this observation and by taking logarithms on both sides of (6), we obtain the following linear constraints on the lograte constants :
where is the m^{th} component of the i^{th} basis vector b^{(i)}. In Additional file 1 we derive an analytical formula for the basis vectors of the null space of [see Equation (S.1)]. As a consequence, the possibly infinite nonlinear Wegscheider conditions given by (6) are equivalent to much more manageable finite linear constraints on the logvalues of the parameters of the closed subsystem, given by
where κ:= ln(k) and is an M_{2} × J matrix that can be easily constructed from knowledge of using (8). Hence, if and only if ln(k) = 0.
We should note that there might be additional linear constraints that we may wish to impose on the logarithms of the kinetic parameters of a biochemical reaction system. Here are some examples:

By using an appropriate experimental procedure, we may be able to accurately measure the equilibrium constant R_{
m
}of the m^{th} reversible reaction. For a (generalized) massaction reaction, we have R_{
m
}:= r_{2m1}/r_{2m}and thus we obtain a linear constraint on the log rate constants of the m^{th} reaction, where is the log value of the measured equilibrium constant.

For a reversible MichaelisMenten reaction, the Haldane relationship implies a linear constraint between the logarithms of the kinetic parameters of a reversible MichaelisMenten reaction and [17].

By using experimental techniques, such as plasmon resonance or atomic force microscopy, it may be possible to obtain a highly accurate measurement of an individual kinetic parameter k_{
j
}. In this case, we must impose the (trivial) constraint , where e_{
j
}is the j^{th} column of the J × J identity matrix.

To reduce the dimensionality of parameter estimation, we may employ a sensitivity analysis approach, such as the one proposed in [18, 19], to identify parameters that do not appreciably influence the cost of estimation. Determining accurate values for these parameters is inconsequential to the behavior of the biochemical reaction system at hand. Therefore, we can fix these parameters to some nominal values throughout model calibration, resulting again in linear constraints of the form .

We may want to expand an existing (validated) thermodynamically feasible model to include additional reactions and molecular species. We can do this by fixing the parameters of the existing model using linear constraints , where is the j^{th} parameter value of the existing model. Then estimation takes place only on the parameters associated with the new reactions.
For a given biochemical reaction system, we can combine all possible linear constraints on the logarithms κ of the kinetic parameters k into a single matrix equation of the form:
where is an appropriately constructed L × J matrix and c is an L × 1 vector of known values determined by the constraints, with L being the number of constraints. Note that if there are no constraints other than the Wegscheider conditions, we would simply have and c= 0.
Model calibration
We will now assume that we have obtained noisy measurements of the concentration dynamics of selected molecular species in an open biochemical reaction system of known stoichiometry, obtained by a set of distinct experiments at discrete time points, within the observation time window . The problem of model calibration we consider in this paper is to determine thermodynamically consistent values for the kinetic parameters , so that the concentration dynamics produced by the estimated system match y in some welldefined sense. Note that, for a given , the dynamics are computed via (1) using initial conditions that correspond to the p^{th} experimental condition. Data y may be obtained by appropriately designed in vivo or in vitro experiments, or by simulating an established and experimentally validated biochemical reaction model whose kinetic parameters are thermodynamically infeasible. Instead of focusing on the quality of the estimated values of the kinetic parameters, it has been argued in [20] that matching predicted and experimental observations of molecular concentrations is the right thing to do due to the 'sloppiness' of biochemical reaction systems (different combinations of parameter values may essentially lead to the same concentration dynamics).
There are many estimation techniques we can use to address the previous problem, such as maximum likelihood or Bayesian inference. The final product of these techniques is a cost function C(κy) used to quantify the overall error between the predicted concentration measurements x, obtained by simulating the biochemical reaction system with kinetic parameter values k = exp{κ}, and the noisy observations y. In an effort to reduce the typically large dynamic range of kinetic parameter values, it is customary to estimate the log values κ instead of k. As a consequence, the problem of interest here is to compute an estimate of the log kinetic parameters κ, so that
where is the set of all κ's satisfying the linear constraints given by (10); i.e.,. For simplicity, we consider in this paper the leastsquares error cost criterion, given by
This error criterion is a consequence of a maximum likelihood approach to parameter estimation under the assumption of normally distributed observation errors. Note that the cost C depends on the log kinetic parameters κ through the molecular concentrations .
In this paper, we refer to the constrained optimization problem given by (11) as Thermodynamically Consistent Model Calibration (TCMC). The importance of TCMC lies on the formulation of the model calibration problem as one of constrained optimization via (11), with constraints that ensure at least the thermodynamic feasibility of the resulting model. A useful observation is that TCMC is agnostic to the choice of the cost function used and the algorithm employed for optimization. Moreover, we can easily transform the constrained optimization problem given by (11) to a standard unconstrained problem. Indeed, a well known result of linear algebra implies that , where , κ_{0} is a J × 1 vector that satisfies is a J × d matrix whose columns form a basis for the null space of matrix , and ℜ ^{d} is the space of all d × 1 real valued vectors. Thus, if we can find a particular solution κ_{0} to the constraints , we can reformulate the constrained optimization problem given by (11) as the following lower dimensional unconstrained problem:
where
Note that we assume here that there is more than one solution to . If only one solution exists, optimization is not necessary, since this solution will be our parameter estimate. On the other hand, if has no solution, then we cannot find a κ that will simultaneously satisfy all necessary constraints, indicating that we must reformulate the problem.
The objective function C_{0} is nonconvex with possibly many local minima. As a consequence, a gradientbased optimization algorithm for solving (13) may prematurely terminate at a local minimum with much larger cost than the globally minimum cost. To ameliorate this problem, we have decided in this paper to use a stochastic optimization algorithm, namely simulated annealing (SA) [21]. Stochastic optimization algorithms can move away from premature local minima, thus resulting in better solutions to optimization problems than when using deterministic techniques. Although many choices exist for optimization, such as the simultaneous perturbation stochastic approximation (SPSA) method employed in our previous work [11] and genetic algorithms, SA is a stochastic search optimization algorithm that enjoys several advantages over other algorithms. In particular, the most important features of SA are ease of implementation and the ability to avoid premature convergence by jumping away from local minima en route to finding a global minimum. In SA, a new value of v is proposed nearby the current value. The proposed value becomes the new value with a certain probability based on cost improvement. If the proposed value is not accepted, then the current value is used. The proposed value is usually drawn from an appropriately chosen probability distribution around the current value (e.g., a Gaussian distribution centered at the current value). See Additional file 1 for a detailed description of the SA algorithm used.
A natural question that arises here is whether different choices for κ_{0} and affect the final result of optimization. If we had an algorithm that could always find the global solution to a nonconvex optimization problem, then the choice of κ_{0} and would have no effect on the solution. Since however global minima are difficult to find, we expect that different choices for κ_{0} and will have some impact on the final solution. Note that it would be advantageous to choose κ_{0} as close as possible to the globally optimal solution. We attempt to do so in our subsequent example by taking κ_{0} to be a solution to that is closest (in the leastsquares sense) to published values. On the other hand, we expect that the choice of will have only a minor effect on optimization, since different matrices amount to scaling or rotating the axes of the parameter space being searched. Good optimization algorithms, such as SA, are expected to be robust to such alterations.
Example
We now demonstrate the proposed TCMC method by reestimating the kinetic parameters of a classical model of the EGF/ERK signaling cascade [13]. This model consists of three compartments (extracellular space, cytoplasm, and endosomal volume), N = 100 biochemical species, and M = 125 reactions. Moreover, it is characterized by 90 different kinetic parameter values, a number that is smaller than the total number of individual reactions, due to the fact that some reactions share the same kinetic parameters, whereas, other reactions are not associated with any parameters.
Although the Schoeberl model has provided valuable insights into the biological mechanisms underlying EGF/ERK signaling, the values of the kinetic parameters published in the literature are thermodynamically infeasible. As a consequence, the concentration dynamics produced by the published model are physically impossible and could not occur in nature. By using TCMC to recompute thermodynamically feasible values for the kinetic parameters, we can construct a physically realizable model whose dynamics are expected to reflect the true behavior of EGF/ERK signaling more accurately than the dynamics produced by the published model.
We use the version of the Schoeberl model published in the BioModels database [22]. Moreover, we employ the same experimental time series data of ERKPP activity used for creating the original model [13]. To simplify implementation of TCMC, we assume here that the Schoeberl model is characterized by J = 2M = 250 kinetic parameters (two parameters per reaction). For an irreversible reaction, we constrain the rate constant of the reverse reaction to be equal to zero. For those reactions not associated with any kinetic parameters, we assign two artificial parameters (one for the forward and one for the reverse reaction) and constrain both their values to be zero.
We implement the following TCMC procedure to reestimate the values of the kinetic parameters in a thermodynamically consistent manner. First, we find the closed subsystem of the Schoeberl model (see Additional file 1). This subsystem consists of N_{0} = 93 molecular species and M_{0} = 83 reversible elementary (monomolecular and bimolecular) reactions. Then, we determine the thermodynamic constraints by using the 93 × 83 stoichiometry matrix of the closed subsystem. It turns out that the rank of the stoichiometry matrix is 65. As a consequence, the closed subsystem contains M_{2} = 83  65 = 18 independent reaction cycles, determined by the columns of matrix , given by Equation (S.1) of Additional file 1 (see the file for details on the reaction cycles). Therefore, the rate constants of the closed subsystem must satisfy 18 independent Wegscheider conditions. It turns out that the published Schoeberl model satisfies only 10 of these conditions. As a matter of fact, the entropy reaction rates, given by (7), associated with 5 independent reaction cycles are negative, in direct violation of the second law of thermodynamics, whereas, the entropy reaction rates of 3 reaction cycles are positive, with the remaining being equal to zero (see Table S.2 in Additional file 1).
Next, we construct matrix and vector c by combining the 18 thermodynamic constraints with 167 linear equality constraints originally built into the model that relate various parameters across reactions (see Additional file 1 for details on these constraints), thus producing the linear constraints . In this case, is a 185 × 250 matrix, whereas, c is a 250 × 1 vector with elements 0 or ∞ (this corresponds to kinetic parameters whose values are fixed to zero). It turns out that rank () = 171, which implies that the dimension of the null space of is d = 79. Subsequently, we find a basis for the null space of and form the 250 × 79 matrix . Moreover, we find a particular solution κ_{0} of that is closest, in the leastsquares sense, among all other solutions to the published thermodynamically infeasible parameter values. We accomplish this by using a wellknown approach for solving constrained leastsquares problems [23]. Finally, by using (12) on the aforementioned experimental time series data and SA, we calculate the 79 × 1 vector by minimizing the cost function C_{0}(vy), given by (14), and set .
We take the thermodynamically feasible log kinetic parameter values κ_{0} to be close to the published kinetic parameter values, since the published values already produce a good match between the experimentally available and predicted molecular dynamics. In this case, TCMC provides a thermodynamically consistent correction of κ_{0}, by means of , that reduces the cost of estimation, thus further improving the match between the experimentally obtained and predicted dynamics.
In Figure 1, we depict the concentration dynamics of ERKPP obtained by the published model (dotted curves), for two different input EGF concentrations, namely 50 ng/mL and 5 ng/mL, whereas, in Figure S1 of Additional file 1 we depict the concentration dynamics of ERKPP obtained by the published model for three additional concentrations of input EGF, namely 0.5 ng/mL, 0.125 ng/mL, and 0.0625 ng/mL. Clearly, the resulting dynamics match the available densitometric data (indicated by the circles) rather poorly.
We should note here that the ERKPP dynamics originally published in [13] seem to match the available data pretty well  see Figure 2F in [13]. However, the model detailed in the original publication cannot be used to reproduce the results. On the other hand, the model available in the BioModels database does not reproduce the results published in the Schoeberl paper but produces dynamics that are very similar to the ones reported in that paper.
The solid curves in Figure 1 and Figure S1 of Additional file 1 depict the ERKPP concentration dynamics estimated by TCMC. TCMC results in a thermodynamically consistent model of the EGF/ERK signalling cascade that produces ERKPP concentration dynamics which match the available experimental data noticeably better than the dynamics produced by the published model. As a matter of fact, model fit between predicted and experimentally measured ERKPP concentration dynamics, quantified by the leastsquares error given by (12), is reduced by 69%. TCMC simultaneously adjusts the values of the kinetic parameters in order to minimize the cost of fitting the ERKPP response to the available densitometric data. This 'collective fitting' strategy has been recognized in the literature [20] as being more desirable than constructing a biochemical reaction system model from individual parameter estimates in a piecewise fashion, which is the case with the published Schoeberl model.
To separate the effect of collective fitting versus imposing the underlying thermodynamic constraints on the kinetic parameters, we use the same simulated annealing algorithm employed by TCMC to estimate the kinetic parameters without including the thermodynamic constraints. The resulting (thermodynamically infeasible) estimated dynamics are depicted by the dashed curves in Figure 1 and Figure S1 of Additional file 1. As expected, these dynamics fit the densitometric data better than the dynamics obtained by the published model. As a matter of fact, model fit between predicted and experimentally measured ERKPP concentration dynamics is reduced by 70% in this case, which is slightly better than the one obtained by TCMC. It will shortly become clear however that the solution obtained without imposing the thermodynamic constraints leads to unrealistic system behavior. In the following, we discuss a number of advantages gained by TCMC over estimating the kinetic parameters using a collective fitting approach that does not consider the underlying thermodynamic constraints.