We first describe a perspective of our method, and then the two models are analyzed to illustrate its performance. The two models were chosen from representative kinetic models for biological phenomena at the molecular level: one model (Model 1) is composed of two variables, analogous to molecular binding and dissociation, such as affinity binding in an antibody cross-link, and the other model (Model 2) is composed of four variables, analogous to a molecular reaction cascade, such as phosphorylation in signal transduction. Notably, we assumed that only one variable is measured among the variables in the two models.
Overview of present method
The key point of this study is the introduction of new constraints obtained by differential elimination into the objective function, to improve the parameter accuracy. Following an explanation of differential elimination, the method of introducing the constraints is briefly described.
Differential algebra aims at studying differential equations from a purely algebraic point of view [6, 7]. Differential elimination theory is a sub theory of differential algebra [3], based on Rosenfeld-Gröbner [4]. The differential elimination rewrites the inputted system of differential equations to another equivalent system according to ranking (order of terms). Here, we provide an example of differential elimination, as shown below, according to Boulier [3, 5].
Assume a model of two variables, x1 and x2, in Fig. 1, which is described by the following system of parametric ordinary differential equations,
,(1)
where k12, k21, ke and Ve are some constants. Here, two molecules are assumed to bind according to Michaelis-Menten kinetics. The differential elimination then produces the following two equations equivalent to the above system.
(2)
When we define the left sides of the above system as C1,t and C2,t, C2,t is composed of x1, its derivatives, and the parameters obtained by eliminating x2, and C1,t is composed of x1, its derivatives, the parameters and x2. Note that x2 in C1,t can be expressed by x1, its derivatives and the parameters in C2,t. Then, the values of C1,t and C2,t can be calculated, if we have time-series data of x1, and they would be zero, if all parameters were exactly estimated. Thus, C1,t and C2,t can be regarded as a kind of error function that expresses the difference between the measured and estimated data.
In general, the typical objective function for evaluating the reproducibility of an experimentally measured time-series for a parameter set is the total relative error, E. The parameter set is then estimated when the total relative error falls below a given threshold. However, the immense searching space of parameter values frequently hinders correct parameter estimation. To overcome this problem, we introduce the constraint between the estimate obtained by differential elimination (DE constraints), C, into the objective function, i.e.,
,(3)
where α is a weighting factor, which is approximately estimated by Pareto optimal solutions for E and C, and then is manually modified (see details in Methods).
Model 1
We analyzed a network model for the binding and dissociation of two molecules (Figure 2A). According to the kinetics of the model (see also Methods), the reference curve of one variable, xAB, was generated (Figure 2B), and two optimization techniques, genetic algorithm (GA) and particle swarm optimization (PSO), were applied to it to evaluate the effect of the introduction of differential elimination constraints (DE constraints) (see details in Additional File 1) into the objective function.
Overall, the introduction of DE constraints into the objective function was highly effective for correctly estimating the parameter values in both GA and PSO (Figure 3). By using GA (Figure 3A), k
p
and km were correctly estimated with the introduction, while the estimation of k
p
failed without the introduction. Indeed, the most frequent values estimated with the introduction (right side of Figure 3A) were found in the bins corresponding to the range between 0.045 and 0.055 for k
p
and between 0.45 and 0.55 for k
m
. In contrast, the most frequent values estimated without the introduction (left side of Figure 3A) were found in the range between 0.065 and 0.075 for k
p
, while those for k
m
were correctly estimated. By using PSO (Figure 3B), k
m
was correctly estimated with the introduction, but k
p
failed, while the estimations of both parameters failed without the introduction. Furthermore, another difference between the estimations with and without the introduction is the distribution form of the estimated values, although the numbers of trial successes in the optimization were different with and without the introduction (see details in Methods). As seen in Figure 3A and 3A, the values with the introduction were sharply distributed, while those without the introduction were widely distributed. The introduction of DE constraints contracted the parameter space to facilitate the estimation of the correct values. As a result, the parameter accuracy was improved by the new objective function with the introduction of DE constraints in Model 1.
Figure 4 clarifies the contraction of parameter space with the introduction of DE constraints into the objective function. Indeed, the estimated values with the introduction by using GA and PSO were concentrated around the correct values (right side of Figure 4). In contrast, the estimated values without the introduction by using the two optimization techniques were broadly distributed (left side). Although the numbers of estimated parameter sets were different with and without the introduction (see details in Methods), the distributions by using the two techniques without the introduction show weak positive correlations. This indicates that the ratio of estimated parameter sets was approximately kept, but the correct estimations failed, without the introduction.
Model 2
We analyzed a network model for the molecular cascade reaction of four molecules (Figure 5A). According to the kinetics of the model (see also Methods), the reference curve of one variable, x1, was generated (Figure 5B), and GA and PSO were applied to it to evaluate the effect of the introduction of DE constraints (see details in Additional File 2) into the objective function.
In Model 2, the introduction of DE constraints into the objective function was also highly effective for correctly estimating the parameter values in both GA and PSO (Figure 6). By using GA (Figure 6A), all three parameters were correctly estimated with the introduction (right side), while the estimations of k31 and k41 failed without the introduction (left side). By using PSO (Figure 6B), all three parameters were also correctly estimated with the introduction (right side), while the estimations of k41 failed without the introduction (left side). Furthermore, the features of the distribution forms of the estimated values were similar to those in Model 1 (Figure 3). As seen in Figure 6A and 6B, the distribution of the estimated values with the introduction was sharp (right side), while that without the introduction was wide (left side). As a result, the parameter accuracy was also improved by the new objective function to contract the parameter space with the introduction of DE constraints in Model 2.
The contraction of parameter space with the introduction of DE constraints into the objective function is shown more clearly in Figure 7. The features of the parameter space in Figure 7 are similar to those in Figure 4. Indeed, the estimated values with the introduction by using GA and PSO were concentrated around the correct values (right side of Figure 7), while the estimated values without the introduction were broadly distributed (left side). In addition, the distributions by using the two techniques without the introduction also show weak positive correlations, similar to the case in Figure 4. Without the introduction, the ratio of estimated parameter sets was approximately maintained, but the correct estimations failed.