Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Comput Biol Chem. Author manuscript; available in PMC 2010 October 1.
Published in final edited form as:
PMCID: PMC2753754

Algebraic Methods for Inferring Biochemical Networks: a Maximum Likelihood Approach


We present a novel method for identifying a biochemical reaction network based on multiple sets of estimated reaction rates in the corresponding reaction rate equations arriving from various (possibly different) experiments. The current method, unlike some of the graphical approaches proposed in the literature, uses the values of the experimental measurements only relative to the geometry of the biochemical reactions under the assumption that the underlying reaction network is the same for all the experiments. The proposed approach utilizes algebraic statistical methods in order to parametrize the set of possible reactions so as to identify the most likely network structure, and is easily scalable to very complicated biochemical systems involving a large number of species and reactions. The method is illustrated with a numerical example of a hypothetical network arising form a “mass transfer”-type model.

Keywords: Biochemical reaction network, law of mass action, algebraic statistical model, polyhedral geometry

1. Introduction

In modern biological research, it is very common to collect detailed information on time-dependent chemical concentration data for large networks of biochemical reactions (see survey papers Crampin et al., 2004; Maria, 2004). Often, the main purpose of collecting such data is to identify the exact structure of a network of chemical reactions for which the identity of the chemical species present in the network is known but a priori no information is available on the species interactions (see e.g. Margolini and Califano, 2007). The problem is of interest both in the setting of classical theoretical chemistry, as well as, more recently, in the context of molecular and systems biology problems and as such has received a lot of attention in the literature over last several decades as evidenced by multiple papers devoted to the topic (Bansal et al., 2007; Fay and Balogh, 1968; Himmeblau et al., 1967; Hosten, 1979; Karnaukhov et al., 2007; Rudakov, 1960, 1970; Schuster et al., 2002; Vajda et al., 1986).

In general, two very different reaction networks might generate identical mass-action dynamical system models, making it impossible to discriminate between them, even if one is given experimental data of perfect accuracy and unlimited temporal resolution. Sometimes this lack of uniqueness is referred to as the “fundamental dogma of chemical kinetics”, although it is actually not a well known fact in the biochemistry or chemical engineering communities (Crampin et al., 2004; Erdi and Toth, 1989; Epstein and Pojman, 2002). Necessary and sufficient conditions for two reaction networks to give rise to the same deterministic dynamical system model (i.e., the same reaction rate equations) are described in Craciun and Pantea (2008), where the problem of identifiability of reaction networks given high accuracy data was analyzed in detail. The key observation is that, if we think of reactions as vectors, it is possible for different sets of such vectors to span the same positive cones, or at least to span positive cones that have nonempty intersection (see Figure 1 for an example).

Figure 1
Identifiability of reaction networks given experimental data. Note that, for a deterministic mass-action model, even if we can estimate the vector K of parameter values with great accuracy, we cannot determine if the “correct” reaction ...

On the other hand, it is often the case that experimental measurements for the study of a specific reaction network or pathway are being collected under many different experimental conditions, which affect the values of reaction rate parameters. Almost always, the reactions of interest are not “elementary reactions”, for which the reaction rates parameters must be constant, but they are so called “overall reactions” that summarize several elementary reaction steps. In that case the reaction rates parameters may reflect the concentrations of biochemical species which have not been included explicitly in the model. In such circumstances the reaction rate parameters are not constant, but rather depend on specific experimental conditions, such as concentrations of enzymes and other intermediate species. Therefore, the estimated vector of reaction rate parameters will not be the same for all experimental conditions, but each specific experimental setting will give rise to one such vector of parameters. However, the set of all these vectors should span a specific cone, whose extreme rays should identify exactly the set of reactions that gave rise to the data.

The purpose of the current paper is to propose a statistical method based on the above geometric considerations, which allows one to take advantage of the inherent stochasticity in the data, in order to determine the unique reaction network that can best account for the results of all the available experiments pooled together. The idea is related to the notion of an algebraic statistical model (as described in Pachter and Sturmfels, 2005, Chapter 1), and relies on mapping the estimated reaction parameters into an appropriate convex region of the span of reaction vectors of a network, using the underlying geometry to identify the reactions which are most likely to span that region. As shown below, this approach reduces the network identification problem to a statistical inference problem for the parameters of a multinomial distribution, which may then be solved for instance using the classical likelihood methods.

2. Maximum Likelihood Inference for a Biochemical Reaction Network

In this section we develop a formal way of inferring a most likely subnetwork of a given conic network (i.e., network represented by a cone like the one in Figure 1) of the minimal spanning dimension. For the inference purpose, in the network of m reactions and d species we assume that the empirical data D = {Ki, i = 1 … , k} [subset or is implied by] Rd is available in the form of (multiple) estimates of the parameters of the system of differential equations corresponding to a hypothesized biochemical network. As illustrated in Figure 1, such networks are in general “unidentifiable” in the sense that different chemical reaction networks may give rise to the same system of differential equations. However, in the stochastic or “statistical” sense it is possible to identify the “most likely” (i.e., maximizing the appropriate likelihood function) network as indicated by the data D.

Throughout this paper we will use a common abuse of notation, where we think of the species A1, … , Ad as the formal generators of a d-dimensional space, and also identify them with the canonical basis of Rd. If y and y′ are non-negative linear combinations of the species, then the ordered pair yy′ is a reaction and its reaction vector is the formal difference y′y. Also, y is called the source complex and y′ is called the target complex of this reaction. Note that the identication of the species with the canonical basis of Rd allows us to think of the reaction vectors as elements of Rd.

2.1. Multinomial model

Consider d species, and m possible reactions with reaction vectors R = {R1, … , Rm} [subset or is implied by] Rd among the species. For more details about how each reaction generates a reaction vector see Craciun and Pantea 2008. Let Rd denote the collection of all (md) positive cones spanned by subsets of d reactions in R.

Denote by cone(R) the positive cone generated by the reaction vectors in R. Let S be the partition of cone(R) obtained by all possible intersections of non-degenerate cones in Rd. Suppose S contains n full-dimensional regions S1, … , Sn; throughout we shall refer to these regions as building blocks, and to n as the number of building blocks.

Let Δm–1 be a probability simplex in Rm and let θ [sm epsilon] Δm–1 be a vector of probabilities associated with the reactions that give rise to R. We assume that these m reactions have the same source complex (i.e., form a conic network), since, as explained in Craciun and Pantea (2008), the identifiability of a network can be addressed one source complex at a time. Define the polynomial map




for i = 1, … , n.

We take1 vol(CSi)vol(C)=0ifvol(C)=0. Define s(θ)=Σσθσ(1)θσ(d) and


In this setting p [sm epsilon] Rn is our statistical model for the data, after we substitute θm=1Σj=1m1θj. Note that we may interpret the monomials θσ(1) (...) θσ(d) in (1) as the probabilities of a given data point being generated by the d-tuple of reactions σ(1), … σ(d). With this interpretation the coordinate pi of the map p in (2) is simply the conditional probability that the data point is observed in Si given that it was generated by a d-tuple of reactions. Note that the map p is rational but, as we shall see below, the model may be re-parametrized into an equivalent one involving only the multilinear map (1).

Let ui denote the number of data points in Si. The log-likelihood function corresponding to a given data allocation is


. Our inference problem is to find

θ^=argmaxθl(θ)subject toΣi=1mθi=1andθi0.

Example 2.1. Consider the two reaction networks described in Figure 1. The model has d = 3 species and a total of m = 5 possible reactions R = {R1, … , R5} = {A1 → 2A1, A1A1 + A2, A1 → 2A3, A1 → 2A2, A1A1 + A3}. In this case there are n = 5 building blocks S1, … , S5 defined by the intersections of all non-trivial reaction cones generated by reaction triples, illustrated in Figure 2.

Figure 2
Configuration of the five building blocks for the possible reactions R = {R1, … , R5} = {A1 → 2A1, A1A1 + A2, A1 → 2A3, A1 → 2A2, A1A1 + A3} of example 2.1.

Thus denoting Cjkl = cone(Rj, Rk, Rl) for any triple {j, k, l} [set membership] {1, … , 5} we have


. Note that the cones C124 and C135 are degenerate and are not involved in the definitions of the Si's. Denoting further vjkl(i) = vol(CjklSi)/vol(Cjkl) for any triple {j, k, l} [set membership] {1, … , 5}, we see that the the map (1) becomes


where the coefficients satisfy Σivjkl(i) = 1 for any triple {j, k, l} appearing on the right-hand-side in the formulas above. The rational map (2) is therefore given by


where the sum in the denominator extends over all distinct triples {j, k, l} excluding {1, 2, 4} and {1, 3, 5}, i.e., the ones corresponding to degenerate cones.

2.2. Multilinear representation

The model representation via a rational map (2) may be equivalently described in terms of a simpler polynomial map (1) as follows. Let us substitute [theta w/ tilde]i = θis−1/d for i = 1, … m and define


Note that gi : >0mRn and l([theta w/ tilde]) = l(θ). Thus we may consider a following more convenient version of (4). Find

θ^=argmaxθ~l(θ~)subject toΣσθ~σ(1)θ~σ(d)=Σig~i(θ~)=Σi=1npi(θ)=1,iθ~i0.


Consider a fixed d-tuple of reactions (say, σ1) and in the formulas for gi (i = 1 … , n) substitute [theta w/ tilde]σ1(1)[theta w/ tilde]σ1(d) = 1 – Σσσ1 [theta w/ tilde]σ(1)[theta w/ tilde]σ(d). Note that the resulting algebraic statistical map is multilinear i.e, linear in one parameter [theta w/ tilde]k when all others are fixed. For instance, as a function of [theta w/ tilde]1 we have


where Σi ai = 0 and Σi bi = 1 and ai, bi are given in terms of [theta w/ tilde]l for l > 2.

By Varchenko's theorem (see Pachter and Sturmfels, 2005, chapter 1) the conditional, one dimensional version of problem (*) may be now solved iteratively for each pi ([theta w/ tilde]1|·), i = 1, … , n by finding a unique root of the score equations in the regions bounded by the ratios −bi/ai.

2.3. Maximization algorithm

Due to the conditional convexity of the one dimensional problems the above considerations suggest that the following algorithm for (local) maximization of l([theta w/ tilde]) should be valid (see also Pachter and Sturmfels, 2005, Example 1.7, page 11):

Algorithm 2.1.

  1. Pick initial vectors [theta w/ tilde] and [theta w/ tilde]old [set membership] Rm.
  2. While |l([theta w/ tilde]) − l([theta w/ tilde]old)| > [sm epsilon]
    • [theta w/ tilde]old[theta w/ tilde]
    • for k=1 to m do
    • compute ai, bi (as functions of [theta w/ tilde]j, j ≠ k)
    • – identify the bounded interval as determined by Varchenko's fromula which is statistically meaningful (there is only one).
    • – use a simple hill-climbing algorithm to find an optimal θ~kopt in that interval
    • – update [theta w/ tilde]kθ~kopt
  3. Recover θ from [theta w/ tilde] by taking θk = [theta w/ tilde]ki[theta w/ tilde]i.

The advantage of the algorithm above is that it reduces a potentially very complicated multivariate optimization problem in which d and m are large to iteratively solving of a simple, univariate one. The disadvantage is that due to its dimension-iterative character the algorithm is seen to be slow and for smaller networks perhaps less efficient than some off-the-shelf optimization algorithms available in commercial software (e.g., some modified hill-climbing methods with random restarts). For that reason in our numerical example below we used the standard Matlab optimization package rather than Alg. 2.1.

In the next section we revert to the notation of Section 1 and the original problem (4). Based on (*) in this section we may thus extend map g to >0m, take s(θ) = 1 in (3) and re-cast the original likelihood maximization problem (4) as

θ^=argmaxθΣiuiloggi(θ)subject toΣσθσ(1)θσ(d)=1andθi0

where the gi's are given by (1).

3. Simulated Numerical Example

In this section we illustrate the ideas discussed above by analyzing a specific numerical example in detail.

If we have d chemical species and data of the form D = {Ki, i = 1 … , k}, then we would hope that the statistical algorithm described above should recover the most likely d reactions out of a given list of md possible reactions, by finding the maximizing vector [theta w/ hat] of the corresponding log-likelihood function. In what follows the setup of the problem is that of (4′). To this end, consider the following four-dimensional example:

equation image

where Ai, i = 0, 1, 2, 3, denote four chemical species. We shall use the above reaction network to simulate “experimentally measured” data and to test the performance of our method outlined in Section 2. To this end we shall augment the above network by including one or more “incorrect” reactions, and shall check whether our likelihood-based algorithm (4′) is able to identify the original “correct” set of four reactions.

3.1. Data generation

Note that the (deterministic) dynamics of the chemical reaction network (5) is governed by linear differential equations of the form


where the parameters γi, i = 0 … 3 are linear combinations of the rate constants, given by the stoichiometry of the true network. In our example each data point Ki [set membership] D (i = 1, … , k) was generated by estimating the set of parameters (γi) of the true reaction network (5). For each one of the k data points the set of parameters (γi) was obtained from rate constants drawn independently from a gamma distribution G(α, λ) with parameters α = 1.5 and λ = 1. In order to identify the coordinates of the points in D, the estimated parameters [gamma with circumflex]i, i = 0, 1, 2, 3, were calculated each time by fitting the trajectories (6) to the time series data points generated from the stochastic process tracing (6) (see Ethier and Kurtz, 1986). The Gillespie algorithm (see Rempala et al., 2006) was used to generate the 20 equally-spaced values of the trajectory of random process on the interval (0, 1) with the fixed initial condition. An example of three random trajectories with independently generated reaction constants values is given in Figure 3. These three trajectories would give rise to three independently estimated sets of values (γi) and consequently to three data points Ki [set membership] D. The fitting was based on the least-squares criterion which is statistically justified for estimation purpose of (γi) in this particular case by an appropriate central limit theorem (see e.g. Ethier and Kurtz, 1986, chapter 11).

Figure 3
An example of generation of the data points Ki [set membership] D for i = 1, 2, 3 via a two-step process of simulation and estimation. Three stochastic trajectories of the reaction network (5) were simulated via Gillespie algorithm with propensity (reaction) ...

We view the single resulting ([gamma with circumflex]0, [gamma with circumflex]1, [gamma with circumflex]2, [gamma with circumflex]3) as coordinates of a point K in the species coordinate system [A0, A1, A2, A3]. This representation of data points is not related to a choice of the reactions; note, however, that if the estimation error is sufficiently small2, most data points lie inside the open convex cone generated by the true reactions (5). As shown in Craciun and Pantea (2008) the coordinates of Ki in the basis given by the reaction vectors in (5) are precisely the estimates of the true rate constants.

The data set D = {Ki, i = 1 … , k} used in the simulation described above was based on multiple batches of k = 10 data points. The first three data points from the first batch are summarized in Table 1.

Table 1
Three sets of parameters (γ0, γ1, γ23) of system (6) corresponding to the trajectories depicted in Figure 3 along with their estimated values (obtained via least-squares fitting) and standard errors of the estimates.

In order to test our method, let us first add one incorrect reaction, A0A2, and from this point on suppose we have no a priori knowledge of the true chemistry; therefore, the five possible reactions are as shown in (7).

equation image

Later in this section we also consider the case where we add not just one, but several incorrect reactions.

3.2. Calculation of the log-likelihood function

In order to obtain an estimate [theta w/ hat] via (4′) one needs to be able to evaluate the map (1), i.e., in addition to the data counts vector u [set membership] Rn in (3) one also needs to know the values of the coefficients of the polynomial map. Whereas the calculation of the exact values is difficult for d > 2, one may typically resort to Monte-Carlo approximations (see e.g. Huggins and Yoshida, to appear). In our current example, for a non-degenerate (i.e., 4-dimensional) cone C, we have computed the approximate relative volumes vol(C ∩ Si)/vol(C) using the following Monte Carlo method. For each cone C we generated N = 2000 points inside C with the corresponding conical coordinates randomly drawn from the four independent gamma G(1.5, 1) random variable and then counted the proportion of the total points falling into CSi i.e., used the approximation

vol(CSi)vol(C)(#points inCSi)Ni=1,n.

With the coefficient values determined as above, the coordinate polynomial maps gi in (1) were easily calculated now by identifying the cones that contained the appropriate building block regions Si.

3.3. Visualization of the chemical network

The geometry in (7) can be visualized in the 3-dimensional subspace W [subset or is implied by] Rn generated by {A1,A2A3}. This follows as all the reaction targets are in this sub-space, and we can understand the configuration of relevant four-dimensional cones by looking at their intersections with W. Each four-dimensional cone with vertex X0 intersects W along a tetrahedron. The intersections of all these tetrahedra cut out the building blocks corresponding to our example (7), as illustrated in Figures Figures44 and and5.5. There are five vertices labeled by numbers corresponding to the five target reactions in (7); they form a six-faced convex polyhedron P. Let C be the intersection of line passing through points 1 and 2, denoted (12), with the plane (345). Then all building blocks are tetrahedra with a vertex at C and the opposite face being one of the six faces of the polyhedron P. For example, the building block (C245) is depicted in Figure 4.

Figure 4
Geometry of building blocks for reaction network (7).
Figure 5
Faces of polyhedron corresponding to reaction network (7).

Not surprisingly, the 10 data points generated in our example were distributed among the building blocks that compose the tetrahedron (1234) corresponding to the true reactions; 6 data points fell inside the building block (C234) and 4 inside (C134). The log-likelihood function was found in this case as


3.4. Maximization of log-likelihood

In order to maximize l(θ) or, equivalently, to minimize −l(θ), rather than using the conditional Algorithm 2.1 we used the Matlab function fmincon for constrained optimization. As in (4′), the constraint is given by the condition s(θ) = Σσ θσ(1)θσ(4) = 1 and comes from the fact that there are 4 reactions in the true network. This constraint also assures that g maps into the probability simplex Δn–1, i.e., defines an algebraic variety (polynomial map) which corresponds to a valid statistical model.

The optimization is repeated 2m–l times (i.e., 16 times for the example for network (7)) with random initial conditions satisfying the constrain. A list of (local) minima was created and entries were merged if they were sufficiently close. The point θ that realized the smallest local minimum was reported together with the percentage of time the algorithm ended up at that particular point (success rate). The output for example (7) given by the customized Matlab function was

  • Minimum of negative log-likelihood: 6.94
  • Theta:
  • 1 1 1 1 0.
  • Hits: 16 out of 16, 100%.

As we may see from these results, in the notation of Figures Figures44 and and55 the algorithm identified the true reactions (targets) 1, 2, 3 and 4 and discarded the incorrect reaction 5.

3.5. More numerical comparisons

We also ran example (5) with, respectively, two, three, four and five incorrect reactions added to the set of four correct ones. The true network was always identified and the success rate (percentage of correct hits for various random initial guess) was high. The results of these experiments are summarized in Table 2. Aditionally, in order to investigate the robustness of our procedure against the lack of precision in the estimates of the γ's, we have added the zero-mean Gaussian term with varying variance to the multiple batches (3000) of k = 10 sets of estimators in the network (7) with one erroneous reaction. The estimated probability of assigning the lowest probability to the extraneous reaction vector as a function of increased noise added to the mean-squared estimated values is presented in Figure 6. As seen from the figure with small to moderate random noise added to the values of the reaction constants the likelihood method is still able to recover the correct set of reactions at remarkably high rate. Not surprisingly, for very noisy data (above one SD) the recovery rate decreases significantly.

Figure 6
The rate of recovery of the correct reactions out of the network (7) as a function of the size of the Gaussian noise added to the least-squares estimates of the (γi) parameters. The solid graph is a smoothed representation of the empirical values ...
Table 2
Summary of average numerical results obtained from multiple experiments, each with k = 10 data points, N = 2000 rays in the Monte Carlo relative volume computation, and 2m−1 optimizations with random initial guess. The analysis was performed on ...

4. Summary and Discussion

We have proposed herein a statistical method for inferring a biochemical reaction network given several sets of data that originate from “noisy” versions of the reaction rate equations associated with the network. As illustrated in the earlier work of some of the authors (Craciun and Pantea, 2008), in the usual deterministic sense such networks are in general unidentifiable, i.e., different chemical reaction networks may give rise to exactly the same reaction rate equations. In practice, the situation is further complicated since the coefficients of the reaction rate equations are estimated from available experimental data, and hence are subject to measurement error and, moreover, their actual values may differ at different experimental conditions, i.e. at different data points. The statistical approach described here is largely unaffected by these problems, as it only relies on the geometry of the network relative to the data distribution, in order to identify the sets of most likely reactions. Hence, the method takes advantage of the algebraic and geometric representation of the network rather than merely the observed experimental values of the network species, as is commonly the case in network inference models based on graphical methods, like e.g. Bayesian or probabilistic boolean networks (see T. Richardson, 2002). Still, in order to use the proposed multinomial parametrization of a biochemical network, the method does require a valid way of mapping the experimentally estimated rate coefficients into the networks' appropriate convex regions, and with very large measurement errors is likely to perform poorly as illustrated in Figure 6. On the other hand, precisely because of the need for the experimental data mapping, the method has a very attractive feature of being able to potentially combine a variety of different data sets obtained by various methods into one set of experimental points placed in the convex hull of the network building blocks. This potential universality of the method's applicability requires further study and possibly a development of additional statistical methodology beyond the scope of our present work. In the current paper our main goal was to present a proof-of-concept example based on simulated data, with a purposefully straightforward but non-trivial model discrimination problem. For the example provided in this paper the method was seen to perform very well, with almost perfect discrimination against incorrect models even as the complexity of the model selection problem increased.

Nonetheless, further study and development is needed to assess how well the method may perform on more challenging and realistic data sets. In particular, one of the aspects of the methodology which was not pursued here, and which could improve its computational scalability, is the utilization of techniques from computational algebra in order to increase the efficiency and further automate the proposed maximization algorithm.


The authors would like to thank Peter Huggins and Ruriko Yoshida for very helpful discussions, and additionally thank Peter Huggins for making available his Matlab script for volume calculations. The research was partially sponsored by FRG grant NSF–DMS 0840695 (Rempala) and NSF–DMS 0553687 (Craciun) as well as by the NIH grant 1R01DE019243-01 (Rempala) and the NIH grant 1R01GM086881 (Craciun).


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

1In general, it may be beneficial to consider various measures vol(·) which are absolutely continuous w.r.t. the usual Lebesque measure. For instance in Section 3 we describe an example where this measure is defined via gamma densities.

2For our current example a brief inspection of the Table 1 indicates a reasonably good agreement between the estimates and the true values of the parameters (γ0, γ1, γ2, γ3) both in terms of actual values as well as the corresponding SE's.

Contributor Information

Gheorghe Craciun, Department of Mathematics and Department of Biomolecular Chemistry, University of Wisconsin-Madison. ude.csiw.htam@nuicarc, Phone: (608)265-3391, Fax: (608)263-8891.

Casian Pantea, Department of Mathematics, University of Wisconsin-Madison. ude.csiw.htam@aetnap.

Grzegorz A. Rempala, Department of Biostatistics, Medical College of Georgia, Augusta, GA 30912. ude.gcm@alapmerg.


  • Bansal M, Belcastro V, Ambesi-Impiombato A, di Bernardo D. How to infer gene networks from expression profiles. Mol. Syst. Biol. 2007;78(3) doi:10.1038/msb4100120. [PMC free article] [PubMed]
  • Craciun G, Pantea C. Identifiability of chemical reaction networks. J. Math. Chem. 2008;44(1):244–259.
  • Crampin E, Schnell S, McSharry PE. Mathematical and computational techniques to deduce complex biochemical reaction mechanisms. Prog. Biophys. Mol. Biol. 2004;86:77–112. [PubMed]
  • Epstein I, Pojman J. An Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos. Oxford University Press; Oxford, UK: 2002.
  • Erdi P, Toth J. Mathematical Models of Chemical Reactions: Theory and Applications of Deterministic and Stochastic Models. Princeton University Press; Princeton, NJ: 1989.
  • Ethier S, Kurtz T. Wiley Series in Probability and Mathematical Statistics: Probability and Mathematical Statistics. John Wiley & Sons Inc.; New York, NY: 1986. Markov Processes.
  • Fay L, Balogh A. Determination of reaction order and rate constants on the basis of the parameter estimation of differential equations. Acta Chim. Acad. Sci. Hun. 1968;57(4):391–401.
  • Himmeblau D, Jones C, Bischoff K. Determination of rate constants for complex kinetic models. Ind. Eng. Chem. Fundam. 1967;6(4):539–543.
  • Hosten L. A comparative study of short cut procedures for parameter estimation in differential equations. Comput. Chem. Eng. 1979;3:117–126.
  • Huggins P, Yoshida R. First steps toward the geometry of cophylogeny arXiv:0809.1908
  • Karnaukhov A, Karnaukhova E, Williamson J. Numerical matrices method for nonlinear system identification and description of dynamics of biochemical reaction networks. Biophys. J. 2007;92:3459–3473. [PubMed]
  • Margolini A, Califano A. Theory and limitations of genetic networks inference from microarray data. Ann. N.Y. Acad Sci. 2007;1115:51–72. [PubMed]
  • Maria G. A review of algorithms and trends in kinetic model identification for chemical and biochemical systems. Chem. Biochem. Eng. Q. 2004;18(3):195–222.
  • Pachter L, Sturmfels B. Algebraic Statistics for Computational Biology. Cambridge University Press; USA: 2005.
  • Rempala G, Ramos K, Kalbfleisch T. A stochastic model of gene transcription: An application to l1 retrotransposition events. J. Theor. Biol. 2006;242(1):101–116. [PubMed]
  • Rudakov E. Differential methods of determination of rate constants of non-complicated chemical reactions. Kinet. Catal. 1960;1:177–187.
  • Rudakov E. Determination of rate constants. method of support function. Kinet. Catal. 1970;11:228–234.
  • Schuster S, Hilgetag C, Woods J, Fell D. Reaction routes in biochemical reaction systems: algebraic properties, validated calculation procedure and example from nucleotide metabolism. J. Math. Biol. 2002;45:153–181. [PubMed]
  • Richardson T, S. P. Ancestral graph markov models. Ann. Statist. 2002;30:962–1030.
  • Vajda S, Valko P, Yermakova A. A direct-indirect procedure for estimating kinetic parameters. Computers Chem. Eng. 1986;10(1):49–58.