In the past decade, there has been a rapid development in systems biology approaches driven by high-throughput data characterizing regulations of genetic networks, interactions among proteins, and reactions in metabolic pathways. These data usually provide a specific scenario of a biological system which may be compared with an alternative system, for instance, expression levels of biomarkers associated with a disease pattern versus healthy controls. Extending the snapshot-type data to condensed data in a time sequence, which is more suitable for profiling the temporal dynamics, provides insights into the functions and underlying regulating mechanisms of the biological system. To gain these insights, mathematical representations of biological systems established with temporal data are highly desirable.

Establishment of proper mathematical representations requires identification of a suitable model with an adequate framework and structure to determine the parameters in the framework. For the structure identification of a model, extensive research has been carried out and mathematical models have been developed to represent the instantaneous rate of a process as an explicit function of all the state variables

*x* ∈

*R*
^{
n
} that have a direct influence on the process. In these representations, the rate of change in each variable

*x*
_{
i
}, is determined by the difference between influx and efflux of the variable (equation

1) and each flux

*v*
_{
i
} is approximated by a product of power-law functions as shown in equation

2:

where

*γ*
_{
i
} represents the rate constant and

*ρ*
_{
ij
} represents the kinetic order of

*v*
_{
i
} with respect to

*x*
_{
j
} [

1–

3]. These approximations have been widely applied to modeling and analysis of biochemical systems for allowing computational simulation of dynamics and parameter estimation for unknown

*γ*
_{
i
} and

*ρ*
_{
ij
}. However, the representations have a low range of accuracy when saturation and cooperativity are represented. To address the reaction rate of enzyme-catalyzed reactions with cooperativity and saturation, a preferred mathematical function is the Hill equation which is described as:

where *V* denotes the maximal rate, *u*
_{
i
} denotes the reaction factor, *K* represents the saturation constant, and *H* denotes the Hill coefficient. The Hill coefficient *H* corresponds to the number of binding sites in the molecule that catalyzes the process [4]. Though there is some disagreement regarding how accurate it is to determine the Hill coefficient *H* by the number of binding sites [5], *H* is generally assumed to be a known constant that can be estimated from experimental data. However, it remains problematic to determine the parameters of *V* and *K*.

For linearly parameterized systems, the least squares method generally gives the optimal estimate of parameters. In addition to the least square approach, an adaptive estimation algorithm serves as a powerful tool to estimate the unknown parameters in ordinary differential equations (ODE) [6–9][10]. For nonlinearly parameterized dynamics, Cao and colleagues have studied the conditions for parameter convergence if the nonlinearly parameterized function satisfies the Lipschitz condition [11]. Qu and colleagues have proposed an adaptive control algorithm for a nonlinearly parameterized system with specific input/control in lieu of requiring the Lipschitz condition with respect to parameters [12]. However, a Hill equation in an ODE does not satisfy Lipschitz condition with respect to the parameter *K* and, generally, it is difficult to apply a continuous control determined by estimated parameters and states of biological systems due to the lack of real time measurements. While estimating the parameters of the Hill equation in ODEs provides accurate prediction of the reactions, it is very difficult to incorporate the continuous evaluation of the states that is needed to better understand the regulation of the biological system. Additionally, it is even more challenging when there is sparse experimental data in discrete time sequences, as is often the case.

Bayesian approaches have been widely used for machine learning, adaptive filters, information theory, and pattern recognition [13–16][17][18]. Specifically, Markov chain Monte Carlo (MCMC) methods have demonstrated to be a powerful inference tool for complex systems raised in behavioral science and computational biology [19, 20][21]. MCMC gains its popularity due to its easy implementation, ability to generate statistically samples from a target high dimensional distribution, good inference performance, and convenience for statistical analysis of results. Therefore, it is very promising to apply MCMC methods to estimate parameters in nonlinearly parameterized dynamics.

The aim of this study is to estimate the unknown parameters

**θ** using a Bayesian approach in nonlinear ODEs representing a biological system as equation (

4):

In this representation, *x* ∈ *R*
^{
n
} denotes the system's state variables, for instance, the concentrations of biochemical factors, and *x*
_{0} is the initial state, *f*(·) is a set of nonlinear functions describing the dynamical property of the biological systems, *u*(*t*) ∈ *R*
^{
l
} is the systems input denoting concentrations of stimuli, and θ ∈ *R*
^{
p
} are parameters that characterize dynamic reactions, *y* ∈ *R*
^{
n
} represents the observed data subject to a Gaussian white noise *ε*(*t*) ~ *N*(0,*σ*
^{2}), *g*(·) represents a measurement function and atypical format will be an identical matrix. We assume we have discrete time series of *y*(*t*), *and*
*u*(*t*) and all parameters in **θ** are positive.

We applied our Bayesian algorithm to estimate unknown parameters in the biological pathways involved in the left ventricle (LV) response to myocardial infarction, which involves inflammatory and fibrotic components typical of a wound healing response. Macrophages begin to infiltrate the LV at day 3 post-MI and are stimulated by interleukin-10 (IL-10) to release transforming growth factor β (TGF-β). In turn, TGF-β stimulates fibroblasts to secrete extracellular matrix components that are necessary for an adequate scar to form. Estimates of the parameters were close to their true value with considerably small estimatiom errors, particularly with regards to small noise variances.