A key problem in systems biology is estimating dynamical models of gene regulatory networks. The mathematical modeling of expression dynamics, combined with model parameter estimation, has been crucial to unraveling complex regulatory programs [1], to recognizing the robustness of the regulatory architecture of the segment polarity genes to variations in initial conditions and parametric variation [2–4], to studying mechanisms of robustness and evolution of the control of the cell cycle in yeast [5, 6], to identifying surprising shifts in the expression domains of the gap genes and the regulatory interactions responsible [7, 8], and to numerous other studies (e.g., [9–16]).

Methods for estimating dynamical models depend on the form of the model and of the data available. We focus on the problem of estimating differential equation models of gene network dynamics based on time series data. Assuming one notion of expression is associated to each gene–for example, mRNA or protein expression level, but not both–then a generic ordinary differential equation (ODE) model for

*N* genes can be formulated as

where **x** is the vector of expression levels of the *N* genes, and **f** produces a vector of time derivatives of expression depending on the current expression levels and on some adjustable parameters *θ.* These parameters typically encode such features as kinetic rates for the production and decay of gene products, and regulatory influences between the genes. The regulatory architecture of the system–that is, which genes' expression derivatives depend on which other genes' expression–may be made explicit in the function **f** (e.g., [2]), or it may be implicit in the parameters (e.g., [7, 8]), in which case optimizing the parameters implicitly determines network architecture. Extensions of our work to modeling both mRNA and protein levels of expression, for example, are straightforward, as would be extensions to functions **f** that depend on time or to delay differential equations, where the derivatives depend on the state of the system in the past. We will also assume that the expression data is collected from the wild type network, though initial conditions may vary. Knock-out or over-expression data has also proven useful in genetic network inference, both in theory [17] and in practice [18]. However, wild-type data is far more common and easier to generate than genetic perturbation data.

To introduce the dynamics estimation approach we investigate, suppose for simplicity that we have access to a single time series

**y** (

*t*
_{0}),

**y** (

*t*
_{1}),…

*,*
**y** (

*t*
_{
T
})

*,* where each vector contains possibly-noisy observed expression values for all

*N* of the genes. Suppose further that we have chosen the form of our model,

**f** (

**x**,

*θ*)

*.* Most often, parameters of an ODE model are estimated by minimizing the squared error between modeled and observed expression

where **x** (*t*
_{
i
}) denotes the solution to the ODE (Eq. 1) with parameters *θ.* The initial conditions are often taken from the observed data, **x** (*t*
_{0}) = **y** (*t*
_{0}) *,* or they may be part of the parameters *θ.* Even when the ODE model is linear, so that
for some matrix *A,* this optimization is not trivial. The solution to such an ODE is given by the matrix exponential **x** (*t*) = *e*
^{
At
}
**x** (*t*
_{0}), so that the dependence of the error on the parameters (*A*) is not straightforward. Still, linear differential equation models have been fit efficiently to expression data by various means, most prominently by recasting the problem into other more convenient forms [9, 19, 20]. When **f** is nonlinear, as is typically the case when trying to make more detailed models of network dynamics, then solving the minimization (Eq. 2) is all the more difficult.

There is another major approach to fitting ODE models, however, via functional data analysis (FDA) [21]. The fundamental idea of FDA is to transform a data series (e.g., the time series **y** (*t*
_{
i
})) into a continuous function (**ŷ** (*t*)). This transformation often involves "denoising" the data, using smoothing splines or some other basis function approximation. Various estimation problems can then be solved in terms of these functions. FDA approaches have made some inroads in the literature on gene expression analysis. It has particularly appeared in papers on dimensionality reduction, clustering or classification for microarray expression time-series data (e.g., [22–24]). More relevantly to the present paper, several works have proposed estimating linear dynamical models from (microarray) expression data [25, 26]. As will be explained in greater detail below, this approach to dynamical modeling allows the estimation problem to be reduced to one of regression, which carries both statistical and computational advantages.

There are several approaches to using FDA ideas in estimating differential equations [

21,

27]. For the general problem of estimating differential equation parameters with nonlinear (or linear) dynamics function

**f**, the most direct approach is to create the smooth of the expression series

**ŷ** (

*t*) and then to differentiate that to produce

. The model parameters can then be fit so that they recreate the estimated derivates as accurately as possible, rather than recreating the observed trajectory as accurately as possible.

This error criterion is different from Equation

2. We will call that one trajectory-based error, and Equation

3 the derivated-based error. The FDA approach thus changes the problem being solved, rather than being an alternative method for solving the traditional formulation of ODE fitting. The derivative-based error has several major computational advantages that allow it to be optimized much more efficiently. First, evaluating the derivatived-based error for any particular parameter set

*θ* is more efficient than for trajectory-based error. It does not require computation of a solution to the ODE (Eq.

1), but only evaluating the dynamics function

**f** along the estimated trajectory

**ŷ** (

*t*). Depending on how

**ŷ** (

*t*) is represented, its derivative,

, may be efficiently calculable as well. A second major advantage is that for typical models, the parameters

*θ* can be partitioned into subsets

*θ*
_{
g
} that are specific to each gene

*g* 's dynamics. In this case, the derivative-based error decomposes into a sum of terms for each gene.

Here,
denotes the element of the vector
corresponding to gene *g,* and similarly *f*
_{
g
} denotes the element of **f** pertaining to gene *g*. Thus, if there are *N* genes and, say, *K* free parameters per gene, the fitting problem is reduced from a single nonlinear *N* × *K* -dimensional problem to a set of *N* independent *K* -dimensional optimization problems. Such a reduction in dimensionality is typically very favorable when solving nonlinear optimization problems. Finally, although the derivative-based error criterion is still a nonlinear optimization problem for arbitrary dynamics functions **f**, informally, the optimization tends to be "less nonlinear" than for the trajectory-based error. In part, this is because the error involves only the evaluation of the dynamics function rather than solutions to the dynamics equation. Typically, **f** is not taken to be anything more complicated than a generalized linear model [28], so that minimizing derivated-based error is a generalized-linear least squares regression problem–a type of problem routinely solved in statistical analyses.

Despite the potential advantages of the FDA approach, we believe it has not been seriously evaluated on the problem of estimating nonlinear models of gene expression dynamics. In particular, neither its efficiency nor its ability to correctly estimate regulatory network architecture have been evaluated. Here, we formulate and test FDA approaches on data from two different real networks, the gap gene system of *Drosophila* melanogaster [7, 8] and the synthetic IRMA network [29], and on simulated data generated by the GeneNetWeaver program [30]. We show that the FDA approach is extremely efficient at fitting nonlinear dynamical models of these data sets. In fact, it is so efficient that we can explicitly enumerate and test all possible regulatory architectures, which, to our knowledge, has never been achievable before for this type of modeling. These enumerations clarify the key regulatory factors, as well as interactions between factors, that explain the observed expression dynamics. We also assess the accuracy with which regulatory relationships are correctly extracted from the data, and compare it to other state-of-the-art fitting approaches. In general, the approach seems as successful as any other at determining which genes regulate which, and is very successful at discriminating the types of regulatory interactions–activation or repression.