BACKGROUND
ordinary differential equations  are commonly used for modeling biological and biochemical systems.
ode models are often subject to considerable uncertainty and/or variability in both initial conditions and parameters  <cit> .
particularly in the case of nonlinear odes, it is essential to have efficient and accurate techniques for analyzing the effects of uncertainty and variability on the dynamical behavior.
the effect of variations in the input on model behavior , the model sensitivity, can be analyzed in various ways.
most numerical approaches address the problem either by computing local sensitivity indices   <cit> , by solving the ode for a statistically large ensemble of random or quasi-random input values  <cit> , or by approximating the functional relationship of the input and output  <cit> .
when uncertainty can be narrowed down to 'small' perturbations, it is often sufficient to study its effects locally.
it is, however, difficult to determine a priori if the uncertainty is small, and in many biological applications the assumption of small perturbations either is questionable  or has been shown to be wrong  <cit> .
to study those effects globally the sensitivity analysis problem can be formulated in terms of an ode with random initial conditions.
the task then is to determine the probability density function , or some features of the pdf, at a time t >  <dig>  approaches to global sensitivity analysis typically rely on monte carlo-type random or quasi-random approximations to the output distribution  <cit> .
in such settings, information on the probability distribution is not directly accessible but encoded in the position and denseness of the sampling points.
the quality of the approximation at a point will thus depend on the coverage of a surrounding region of that point, and involves some density estimation steps that might require problem specific knowledge  <cit> .
investigating regions of low probability with high accuracy generally requires prohibitively large numbers of sample points .
in  <cit> , a global sensitivity analysis of odes was proposed based on recasting the problem in terms of a first-order partial differential equation  that describes the evolution of the associated pdf .
knowledge of the entire pdf, in contrast to some specific observables , fully characterizes the impact of uncertainty and variability.
the pde view also facilitates model assessment and parameter estimation studies, since the model likelihood naturally offers a way of assessing deviations of the model output from experimental data.
the availability of the likelihood then gives access to a wide range of well-established numerical tools for model assessment and parameter estimation  <cit> .
to numerically solve the particular class of pdes that arise in odes with random initial data, finite difference schemes are commonly applied  <cit> .
these methods typically become computationally prohibitive beyond three dimensions.
moreover, the numerical treatment of pdes is generally less accessible to practitioners.
the unscented kalman filter  for time-continuous systems  <cit> , originally designed for the approximate solution of the closely related fokker-planck equation, has also been applied to solve this pde and to obtain sensitivity estimates of ode models  <cit> .
the ukf yields normal approximations of the output pdf, where the estimated mean and variance are second-order approximations of the true output mean and variance.
as we illustrate in our examples, non-linear ode models can easily give rise to strongly non-normal output distributions, even if the initial distribution is normal, such that in this case an analysis using the ukf will give misleading results.
it is however well-known that first-order pdes can be solved using the method of characteristics  <cit> .
with the method of characteristics, the pdf can be computed along trajectories of the system by solving the original ode with an additional dimension for the density.
this method thus provides a bridge from the intricate pde description back to the ode setting, where numerical solutions are readily accessible to practitioners.
propagating points along the trajectories of the system is a common feature of the method of characteristics, monte carlo methods and the unscented kalman filter.
in contrast to the other methods, however, the method of characteristics also propagates information  along the trajectories, and these values are exact up to the accuracy of the ode solver.
the additional computational costs of solving the extended ode are negligible as the extra dimension does not necessitate additional sample/discretization points.
in fact, since the pdf at a given state value is directly computable using the method of characteristics and its accuracy at a point does not depend on the denseness of sampling points in the surrounding region, considerably fewer discretization points are sufficient to obtain accurate estimates as compared to monte carlo-based methods  <cit> .
although the method of characteristics is known and widely used in other fields such as meteorology, e.g.
<cit> , its applicability and potential has not been adequately recognized among the biological modeling community.
in this article we review the formulation of odes with random inputs in terms of the pde description and its solution via the method of characteristics.
we demonstrate the use of the method of characteristics for sensitivity analysis as well as for likelihood-based model assessment and parameter estimation and illustrate its benefits and limitations by means of numerical case studies.
RESULTS
methodology
odes with uncertain or variable input
consider an ode of the form
x˙=f,x=x <dig>
with x ∈ ℝd.
typically, x belongs to some extended state space comprising the state variable z ∈ ℝn and the parameters p ∈ ℝm of the system with d = n + m. that is
x:=andf:=0),
where f describes the dynamics of z given the parameter values p, which are assumed to remain constant in time.
for example, z might denote the concentration of a molecular species of a metabolic network or signalling pathway, and p the associated reaction rate constants.
under the assumption that f : ℝd → ℝd is continuously differentiable with respect to x, the initial value problem  is known to have a unique solution x for t ≥  <dig>  <cit> .
uncertainty and variability in the model input can be modeled by assuming that x <dig> = x <dig> is a random variable with pdf u <dig> : ℝd → ℝ. consequently, the solution {xt}t≥ <dig> of the initial value problem  is a stochastic  process.
for any t ≥  <dig>  let us denote by u = u, u:ℝ0+×ℝd→ℝ, the pdf of the random variable xt, i.e., ℙ=∫−∞xudy.
in this setting, the sensitivity of the ode  at time t >  <dig> with respect to the initial density u <dig> amounts to computing the density u at t = t. for continuously differentiable f and u, the density u satisfies the first-order linear pde  <cit>
∂∂tu=−div=−∑i=1d∂∂xi, u=u <dig>
where fi denotes the i-th component of f , and divf is the sum of the partial derivatives of fi .
note that  is the fokker-planck equation corresponding to a stochastic differential equation with zero diffusion  <cit> .
computing the pdf along solutions of the ode
it is well-known that first-order pdes of the form  can be solved using the method of characteristics  <cit> , which in our case is identical to the solution of  along the solution of the initial value problem .
define ρ = u), where x denotes the solution of the initial value problem .
applying the chain rule, ρ obeys the ode:
ρ˙=∂∂tu)+∑i=1d∂∂xiu)⋅x˙i=−div))+grad u)⋅x˙.
noting that x˙=f and div) = divf  · u + grad u · f , where grad u= is the gradient of u, and using  we finally obtain from
ρ˙=−divf)⋅ρ.
the pde  can thus be solved pointwise for each x <dig> by solving the original ode  together with an extra dimension for the density ρ, i.e.,
x˙=fρ˙=−divf⋅ρ,
with initial conditions x = x <dig> and ρ = u <dig>  since x ∈ ℝd, the new system  has d +  <dig> dimensions.
the computational effort for solving the extra dimension is negligible compared to the information gained.
the distribution u at time t >  <dig> can be approximated by the following two steps:
discretize the region of interest in the state space ℝd, resulting in discretization points ξi, i =  <dig> ...,n.
for the initial values , u0)) ∈ ℝd+ <dig>  i =  <dig> ...n, solve the extended ode  to compute the density u)) at t = t.
this procedure directly yields the density values u) along the trajectories ξi, t ∈  up to the accuracy of the ode solver.
in the subsequent examples we will further show how to modify step  in order to investigate the distribution on particular regions of the state and parameter space.
in comparison, monte carlo-based methods require density estimation subsequent to solving the ode for the sample points:
sample the initial distribution u <dig>  resulting in sampling points ξi ∈ ℝd, i =  <dig> ...,n.
for each sampling point ξi solve the original ode  to obtain ξi  at t = t.
estimate the density from the propagated points ξi, i =  <dig> ...,n; for example, by considering a neighborhood b of a point x and approximating the pdf by the relative frequency, i.e., u ≈ # {ξi ∈ b}/n.
in contrast to the method of characteristics, the quality of the approximation u at a single point x depends on the total number of sample points  <cit> .
to obtain good estimates of the pdf in regions with low probability, monte carlo-based methods typically require very large sample sizes.
the quality of the approximation may moreover depend on the structure of the pdf to be estimated, and might therefore necessitate problem specific knowledge.
an example is given below, where the pdf is concentrated on a non-linear manifold due to fast, contracting directions.
we illustrate the advantages of the method of characteristics for sensitivity analysis of ode models by two examples in gene expression.
we will see that descriptors such as mean and variance may provide only poor information about the pdf.
our first example demonstrates the benefits of the method of characteristics in terms of an efficient computation of the pdf in regions with low probability.
in the second example, the pdf contracts onto a lower-dimensional manifold of the state space, and the method of characteristics provides the density values directly along that manifold.
in a third example we further show how the method can be used to compute the likelihood of an ode model and thus facilitates comparison to experimental data.
for the sake of simplicity, we choose normal initial distributions to describe state and parameter uncertainty and variability.
the method of characteristics, however, provides a general strategy to compute model uncertainty/variability and likelihood for arbitrary distributions of initial values and parameters with the only restriction that the associated pdf u must be continuously differentiable.
the choice of normal initial distributions also illustrates that the assumption of normal output distributions--underlying the ukf and least squares parameter estimation--is easily violated for non-linear odes, even if the initial distribution is normal.
examples
sensitivity analysis and the impact of variability
example  <dig>  consider a protein x activating its own expression by cooperatively binding to the promoter that positively regulates its own expression .
assuming that x is diluted due to cell-growth with rate constant kd >  <dig>  the concentration x of the protein is modeled by the ode  <cit>
x˙=f=vmax⋅xβkβ+xβ−kd⋅x.
the first summand describes the activation of gene expression in terms of a hill function  <cit> , where vmax >  <dig> denotes the maximal expression rate, k >  <dig> the protein concentration corresponding to the half-maximal expression rate, and β >  <dig> describes the cooperativity of promoter activation.
a saturable synthesis was chosen to reflect a finite gene copy number encoding for x. the second summand describes the dilution of x. an initial density u <dig> may, for example, represent the abundance of the protein x in individuals of a population of cells.
most of the distribution is translated linearly, since for large x the hill-type activation term in  is approximately constant.
for smaller values of x, however, the right hand side of  is strongly nonlinear.
for values of x close to zero the dynamics are very slow, which causes the formation of a heavy tail.
this is a characteristic feature of bimodality associated with positive feedback loops.
due to only few discretization points on the interval x ∈  <cit>  at t =  <dig> the structure of the heavy tail is only poorly resolved.
in contrast to conventional monte carlo methods, the method of characteristics allows for a refined computation of the density on any sub-interval of interest by the following two steps:
discretize the sub-interval of interest, and solve the original ode  for each discretization point ξi, i =  <dig> ...n, backward in time from t = t to  <dig> to obtain ξi .
for the initial values , u0)), i =  <dig> ...n, solve  forward in time to compute the density along the trajectories.
the two-step procedure is illustrated in figure  <dig>  this way the distribution can be studied on arbitrary regions of the state and parameter space, and no subsequent normalization is required.
we used the two-step procedure to obtain an improved resolution of the heavy tail ) and observe the formation of a second mode close to the origin.
as an interpretation, this may imply that for this part of the population, x does not reach a certain threshold concentration within the given time interval .
in other applications, such as toxicological risk assessment studies  <cit> , the information may analogously be used to determine the percentage of a population that exceeds or remains below a certain toxicological threshold.
to illustrate that the method of characteristics can be applied to any continuously differentiable pdf, not only normal pdfs, we repeated the above computations for an initial exponential distribution with mean μ =  <dig>  &2).
with standard monte carlo-based sensitivity approaches such localized information is difficult to obtain, when the region of interest has only low probability.
the heavy tail in figure  <dig> has a total probability of approximately  <dig> , which means that in expectation only  <dig> % of the monte carlo sampling points will lie in the interval  <cit> .
with the two-step procedure we used  <dig> discretization points to approximate the heavy tail.
compared to our approach, it would require  <dig>  <dig> monte carlo sample points to expect the same coverage on the interval  <cit> , and the subsequent step of density estimation required by monte carlo methods further impacts the approximation quality  <cit> .
apart from variability in the initial concentration x, we can additionally account for variability in the parameter values.
we computed the density for a state space extended according to  by the cooperativity β, by the maximal expression rate vmax and by both β and vmax, i.e., with extended state space variables ', ' and ', respectively.
the initial distribution was assumed to be a joint normal distribution, where x had mean and variance as before, and the means of vmax and β were set to  <dig> and  <dig>  respectively, each with variance  <dig> .
using the above two-step procedure, we computed the densities at t =  <dig>  figure  <dig> depicts the marginal distributions of protein concentration x under the different scenarios of variable/uncertain parameters obtained by subsequent integration over the parameter dimensions that are not considered.
for example, the marginal density in x is obtained from the joint density in  by
u=∫udβ.
numerically, we discretized the above integral using the midpoint rule.
comparison of the distributions indicates that variability in the cooperativity β has minor impact on the final variability of the protein concentration x .
the corresponding joint distributions of ' and ' are shown in figure  <dig> as well as the two-dimensional marginal distributions for the three-dimensional case obtained by integrating only β or vmax, respectively.
same as in the one-dimensional marginal distribution, it can be seen that the variability in protein concentration is mainly dominated by the variability in the maximal expression rate vmax.
the numerical integration ), necessary to visualize multivariate densities for d ≥  <dig> or to compute observables such as mean and variance, is currently the computationally limiting step in the application of the method of characteristics, since it requires a regular  discretization of the state and parameter space, which becomes prohibitive in high dimensions.
for low- and moderate-dimensional systems of odes, however, the method of characteristics provides a more efficient and equally simple alternative to conventional approaches for studying uncertainty and variability of ode systems: the density information obtained via the method of characteristics is expectedly more accurate than estimates obtained with monte carlo methods  <cit>  and considerably richer than a simple characterization of the variability/uncertainty by means of certain indicators .
it moreover facilitates further statistical analysis such as model assessment and parameter estimation, e.g., by means of information theoretical approaches  <cit> , as illustrated later in example  <dig> for the case of parameter estimation.
example  <dig>  consider the genetic toggle switch  <cit> , where two proteins x <dig> and x <dig> mutually repress the other protein's expression .
in  <cit> , the concentrations x <dig> and x <dig> of x <dig> and x <dig>  respectively, are modeled by the two-dimensional system
x˙1=α11+β1−kd⋅x1x˙2=α21+β2−kd⋅x <dig>
the parameter α <dig> >  <dig> represents the effective expression rate of protein x <dig>  and β <dig> >  <dig> describes the repression cooperativity of the promoter that regulates the expression of x <dig> by x <dig>  analogously, α <dig> and β <dig> describe the effective expression rate of protein x <dig> and its promoter's cooperativity of repression by x <dig>  respectively.
the parameter q <dig>  corresponds to the concentration of x <dig>  that represses the promoter activity of x <dig>  by 50%.
the parameter kd again denotes the dilution rate constant.
assume that the concentrations x <dig> and x <dig> at t =  <dig> have a joint normal distribution with mean μ = '  and covariance matrix Σ = diag{ <dig> ,  <dig> } .
the initial density u <dig> is shown in figure  <dig> along with the vector field defined by  with symmetric parameters α <dig> = α <dig> = α =  <dig> , β <dig> = β  <dig> = β =  <dig> , q <dig> = q <dig> =  <dig>  and kd =  <dig> .
we computed the density u solving  with f as in .
at t =  <dig>  the majority of the probability is concentrated on the slow manifold of the vector field with large variance along the manifold ).
such steep distributions on lower-dimensional manifolds pose problems to many other methods: methods based on the estimation of mean and co-variance fail to describe the true shape of the distribution.
most pde solvers have numerical problems with such steep gradients.
monte carlo-based density estimation may yield a too coarse-grained approximation of the true density if knowledge of the manifold is not provided a priori ; that is, to arrange the bins of a histogram or for kernel density estimation the centers of the kernel functions along the manifold.
using the method of characteristics we avoid these problems: the steep gradients do not pose any problems to the independently computed trajectories, which directly yield the density values on the attractor manifold.
the method thus does not require problem specific ingenuity.
deterministic models in a likelihood setting for comparison with experimental data
so far we have discussed the use of the method of characteristics to study the sensitivity of ode models.
it however also offers benefits when comparing model output with experimental data for model assessment, such as validation/falsification and selection between different models, or parameter estimation.
an exact match of the deterministic model with the data is unlikely, and a quantification of the mismatch remains a critical issue.
some numerical approaches are based on verifying that the experimental data lies within regions of the state space that are reachable with certain parameter sets of the  ode model  <cit> .
most other approaches assess the mismatch based on root mean square deviations and select a model or parameters based on a minimization of these errors.
such least squares approaches are based on the assumption that deviations are due to additive gaussian noise, usually assumed to be identical and independent at all points in time  <cit> .
as we have seen in the previous examples, a normal distribution can typically not be expected for general, nonlinear ode models, even if the initial distribution is assumed to be normal.
in addition, classical least squares approaches do not allow taking into account prior information on the initial condition or on the parameter to be estimated.
stochastic approaches offer a natural way of assessing deviations of the model output from data based on the likelihood function.
given a set of data points d = {ξ <dig> ...,ξn }, the likelihood of a model is defined as the probability that the model predicts this data.
for continuous state space models m, the likelihood ℒ given d is defined via the pdf by the product of its values at the data points
ℒ=∏i=1nu,
where u denotes the density function of the output distribution of the model m. based on the likelihood, there are many methods available for model assessment and parameter estimation  <cit> .
we can directly use the method of characteristics to efficiently compute the likelihood of an ode model for given data points.
we first consider the case that parameters are known and only initial conditions are affected by uncertainty.
given prior knowledge of this uncertainty defined by the initial density u <dig>  and further given data points d = {ξ <dig> ...,ξn} at a time t >  <dig>  the likelihood for each single data point can be computed analogously to the two-step procedure in example 1:
solve the ode  for each of the n data points ξi ∈ ℝd, i =  <dig> ...,n, backward in time from t = t to  <dig> to obtain ξi.
for the initial values , u0)) ∈ ℝd+ <dig>  i =  <dig> ...n, solve the extended ode  forward in time to obtain the likelihood values for each data point u).
in accordance with , the likelihood ℒ of the model given the data is the product of the likelihood for each single data point, i.e., ℒ=∏i=1nu).
as the number n of data points is in most cases comparatively small, the above two-step procedure is an efficient way to obtain the likelihood of an ode model.
example  <dig>  reconsider the ode model  of autoregulated gene expression from example  <dig>  assume that we want to estimate the maximal expression rate vmax of the protein x based on an observation ξ =  <dig>  at time t =  <dig>  for sake of clarity we only consider one data point.
for several data points the same procedure as described below applies, and the final likelihood is given by the product of the single-data likelihood values.
as we are interested in the likelihood of different values of vmax, we consider the autoregulation model  extended according to  by vmax.
we apply the above two-step procedure to a representative ensemble {υ <dig> ...,υn} of values of vmax.
for each pair , υi) the backward-solution yields a different value , υi).
given prior knowledge in terms of a joint pdf u <dig> for x <dig> and vmax, the forward-solution of  with initial conditions , υi), u <dig>  υi)) then yields the likelihood values u, υi)) associated with each vi, i =  <dig> ...,n.
we computed the likelihood of a set of equidistant values of vmax ∈  <cit>  using the same parameter values as in example  <dig>  for two different scenarios of prior information:  a joint normal distribution of x <dig> and vmax with parameters μ = ' and Σ = diag{ <dig> ,  <dig> } , and  a joint distribution of x <dig> and vmax, where x <dig> is normally distributed with n and vmax is independently uniformly distributed on the interval  <cit>  .
the first scenario accounts for prior knowledge of x <dig> and vmax, where a more or less precise knowledge of vmax is given  is small).
accordingly, the maximum-likelihood estimate is vmax*≈ <dig> close to the prior mean of vmax.
in the second scenario no prior information on vmax was imposed .
the maximum-likelihood estimate of vmax is therefore solely determined by the value of vmax that yields the initial value closest to μ =  <dig>  since the data point ξ =  <dig> is relatively unlikely for larger vmax ), scenario  yields a smaller maximum-likelihood estimate of vmax*≈ <dig> .
CONCLUSIONS
studying the effects of uncertainty and variability in initial values or parameters of ode models can be computationally intensive, since it generally involves solving the system a large number of times.
the method of characteristics offers a simple yet accurate alternative to conventional approaches for small- and moderate-dimensional systems.
the approach does not assume a particular shape of neither input nor output distribution, it only requires the pdf to be sufficiently smooth  and yields density values that are exact up the accuracy of the ode solver used.
our first two examples illustrate how a precise characterization of the model uncertainty/variability can be obtained with only few trajectories.
in this context we also demonstrated that the analysis can be efficiently restricted to certain sub-regions of the state/parameter space.
one limitation of the two-step procedure used for the latter analysis is that for chaotic models the backward-forward solution of the ode is ill-conditioned  <cit> .
in such cases one may resort to approximate techniques as the ukf  <cit> , but care must be taken that the assumption of the kalman filter, that the distribution remains approximately normal, is satiesfied.
another limitation, currently the main one, is the need for uniform grids when dimensions are to be integrated from the final density.
but we anticipate that the method of characteristics will prove useful in the context of error control for approximate solution methods of eq.
or  such as monte carlo or the unscented kalman filter in higher dimensions by providing exact values of the pdf at particular points in state space.
as another application we considered the comparison of model results with experimental data.
for deterministic models numerical approaches typically rely on root mean squared errors to quantify deviations.
their minimization can be interpreted as the maximum-likelihood estimate based on the assumption that deviations are normally  distributed at all times.
while being a simple assumption, for general nonlinear ode models it is rarely expected to hold.
in the third example we described a simple framework, where the method of characteristics was applied to maximum-likelihood parameter estimation based on a distribution that accounts for prior knowledge of parameters and initial values and for the system dynamics.
we provide matlab files illustrating the method of characteristics in additional file  <dig>
authors' contributions
ayw, wh and rhm planned and performed the research, ayw performed the numerical simulations, all authors contributed to the design and the writing of the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
matlab files illustrating the usage of the method of characteristics.
click here for file
acknowledgements
this work was supported in part by the science foundation ireland, grant sfi07/rpr/i <dig>  we thank ken r. duffy  for his critical reading of the manuscript.

 ordinary differential equations are commonly used for modeling biological and biochemical systems. ode models are often subject to considerable uncertainty and/or variability in both initial conditions and parameters. particularly in the case of nonlinear odes, it is essential to have efficient and accurate techniques for analyzing the effects of uncertainty and variability on the dynamical behavior. the effect of variations in the input on model behavior, the model sensitivity, can be analyzed in various ways. most numerical approaches address the problem either by computing local sensitivity indices, by solving the ode for a statistically large ensemble of random or quasi-random input values, or by approximating the functional relationship of the input and output. when uncertainty can be narrowed down to'small' perturbations, it is often sufficient to study its effects locally. it is, however, difficult to determine a priori if the uncertainty is small, and in many biological applications the assumption of small perturbations either is questionable or has been shown to be wrong. to study those effects globally the sensitivity analysis problem can be formulated in terms of an ode with random initial conditions. the task then is to determine the probability density function, or some features of the pdf, at a time t > approaches to global sensitivity analysis typically rely on monte carlo-type random or quasi-random approximations to the output distribution. in such settings, information on the probability distribution is not directly accessible but encoded in the position and denseness of the sampling points. the quality of the approximation at a point will thus depend on the coverage of a surrounding region of that point, and involves some density estimation steps that might require problem specific knowledge. investigating regions of low probability with high accuracy generally requires prohibitively large numbers of sample points. in, a global sensitivity analysis of odes was proposed based on recasting the problem in terms of a first-order partial differential equation that describes the evolution of the associated pdf. knowledge of the entire pdf, in contrast to some specific observables, fully characterizes the impact of uncertainty and variability. the pde view also facilitates model assessment and parameter estimation studies, since the model likelihood naturally offers a way of assessing deviations of the model output from experimental data. the availability of the likelihood then gives access to a wide range of well-established numerical tools for model assessment and parameter estimation. to numerically solve the particular class of pdes that arise in odes with random initial data, finite difference schemes are commonly applied. these methods typically become computationally prohibitive beyond three dimensions. moreover, the numerical treatment of pdes is generally less accessible to practitioners. the unscented kalman filter for time-continuous systems, originally designed for the approximate solution of the closely related fokker-planck equation, has also been applied to solve this pde and to obtain sensitivity estimates of ode models. the ukf yields normal approximations of the output pdf, where the estimated mean and variance are second-order approximations of the true output mean and variance. as we illustrate in our examples, non-linear ode models can easily give rise to strongly non-normal output distributions, even if the initial distribution is normal, such that in this case an analysis using the ukf will give misleading results. it is however well-known that first-order pdes can be solved using the method of characteristics. with the method of characteristics, the pdf can be computed along trajectories of the system by solving the original ode with an additional dimension for the density. this method thus provides a bridge from the intricate pde description back to the ode setting, where numerical solutions are readily accessible to practitioners. propagating points along the trajectories of the system is a common feature of the method of characteristics, monte carlo methods and the unscented kalman filter. in contrast to the other methods, however, the method of characteristics also propagates information along the trajectories, and these values are exact up to the accuracy of the ode solver. the additional computational costs of solving the extended ode are negligible as the extra dimension does not necessitate additional sample/discretization points. in fact, since the pdf at a given state value is directly computable using the method of characteristics and its accuracy at a point does not depend on the denseness of sampling points in the surrounding region, considerably fewer discretization points are sufficient to obtain accurate estimates as compared to monte carlo-based methods. although the method of characteristics is known and widely used in other fields such as meteorology, e.g., its applicability and potential has not been adequately recognized among the biological modeling community. in this article we review the formulation of odes with random inputs in terms of the pde description and its solution via the method of characteristics. we demonstrate the use of the method of characteristics for sensitivity analysis as well as for likelihood-based model assessment and parameter estimation and illustrate its benefits and limitations by means of numerical case studies.  methodology odes with uncertain or variable input consider an ode of the form x˙=f,x=x with x ∈ ℝd. typically, x belongs to some extended state space comprising the state variable z ∈ ℝn and the parameters p ∈ ℝm of the system with d = n + m. that is x:=andf:=0), where f describes the dynamics of z given the parameter values p, which are assumed to remain constant in time. for example, z might denote the concentration of a molecular species of a metabolic network or signalling pathway, and p the associated reaction rate constants. under the assumption that f: ℝd → ℝd is continuously differentiable with respect to x, the initial value problem is known to have a unique solution x for t ≥. uncertainty and variability in the model input can be modeled by assuming that x = x is a random variable with pdf u: ℝd → ℝ. consequently, the solution {xt}t≥ of the initial value problem is a stochastic process. for any t ≥ let us denote by u = u, u:ℝ0+×ℝd→ℝ, the pdf of the random variable xt, i.e., ℙ=∫−∞xudy. in this setting, the sensitivity of the ode at time t > with respect to the initial density u amounts to computing the density u at t = t. for continuously differentiable f and u, the density u satisfies the first-order linear pde ∂∂tu=−div=−∑i=1d∂∂xi, u=u where fi denotes the i-th component of f, and divf is the sum of the partial derivatives of fi. note that is the fokker-planck equation corresponding to a stochastic differential equation with zero diffusion. computing the pdf along solutions of the ode it is well-known that first-order pdes of the form can be solved using the method of characteristics, which in our case is identical to the solution of along the solution of the initial value problem. define ρ = u), where x denotes the solution of the initial value problem. applying the chain rule, ρ obeys the ode: ρ˙=∂∂tu)+∑i=1d∂∂xiu)⋅x˙i=−div))+grad u)⋅x˙. noting that x˙=f and div) = divf · u + grad u · f, where grad u= is the gradient of u, and using we finally obtain from ρ˙=−divf)⋅ρ. the pde can thus be solved pointwise for each x by solving the original ode together with an extra dimension for the density ρ, i.e., x˙=fρ˙=−divf⋅ρ, with initial conditions x = x and ρ = u since x ∈ ℝd, the new system has d + dimensions. the computational effort for solving the extra dimension is negligible compared to the information gained. the distribution u at time t > can be approximated by the following two steps: discretize the region of interest in the state space ℝd, resulting in discretization points ξi, i =...,n. for the initial values, u0)) ∈ ℝd+ i =...n, solve the extended ode to compute the density u)) at t = t. this procedure directly yields the density values u) along the trajectories ξi, t ∈ up to the accuracy of the ode solver. in the subsequent examples we will further show how to modify step in order to investigate the distribution on particular regions of the state and parameter space. in comparison, monte carlo-based methods require density estimation subsequent to solving the ode for the sample points: sample the initial distribution u resulting in sampling points ξi ∈ ℝd, i =...,n. for each sampling point ξi solve the original ode to obtain ξi at t = t. estimate the density from the propagated points ξi, i =...,n; for example, by considering a neighborhood b of a point x and approximating the pdf by the relative frequency, i.e., u ≈ # {ξi ∈ b}/n. in contrast to the method of characteristics, the quality of the approximation u at a single point x depends on the total number of sample points. to obtain good estimates of the pdf in regions with low probability, monte carlo-based methods typically require very large sample sizes. the quality of the approximation may moreover depend on the structure of the pdf to be estimated, and might therefore necessitate problem specific knowledge. an example is given below, where the pdf is concentrated on a non-linear manifold due to fast, contracting directions. we illustrate the advantages of the method of characteristics for sensitivity analysis of ode models by two examples in gene expression. we will see that descriptors such as mean and variance may provide only poor information about the pdf. our first example demonstrates the benefits of the method of characteristics in terms of an efficient computation of the pdf in regions with low probability. in the second example, the pdf contracts onto a lower-dimensional manifold of the state space, and the method of characteristics provides the density values directly along that manifold. in a third example we further show how the method can be used to compute the likelihood of an ode model and thus facilitates comparison to experimental data. for the sake of simplicity, we choose normal initial distributions to describe state and parameter uncertainty and variability. the method of characteristics, however, provides a general strategy to compute model uncertainty/variability and likelihood for arbitrary distributions of initial values and parameters with the only restriction that the associated pdf u must be continuously differentiable. the choice of normal initial distributions also illustrates that the assumption of normal output distributions--underlying the ukf and least squares parameter estimation--is easily violated for non-linear odes, even if the initial distribution is normal. examples sensitivity analysis and the impact of variability example consider a protein x activating its own expression by cooperatively binding to the promoter that positively regulates its own expression. assuming that x is diluted due to cell-growth with rate constant kd > the concentration x of the protein is modeled by the ode x˙=f=vmax⋅xβkβ+xβ−kd⋅x. the first summand describes the activation of gene expression in terms of a hill function, where vmax > denotes the maximal expression rate, k > the protein concentration corresponding to the half-maximal expression rate, and β > describes the cooperativity of promoter activation. a saturable synthesis was chosen to reflect a finite gene copy number encoding for x. the second summand describes the dilution of x. an initial density u may, for example, represent the abundance of the protein x in individuals of a population of cells. most of the distribution is translated linearly, since for large x the hill-type activation term in is approximately constant. for smaller values of x, however, the right hand side of is strongly nonlinear. for values of x close to zero the dynamics are very slow, which causes the formation of a heavy tail. this is a characteristic feature of bimodality associated with positive feedback loops. due to only few discretization points on the interval x ∈ at t = the structure of the heavy tail is only poorly resolved. in contrast to conventional monte carlo methods, the method of characteristics allows for a refined computation of the density on any sub-interval of interest by the following two steps: discretize the sub-interval of interest, and solve the original ode for each discretization point ξi, i =...n, backward in time from t = t to to obtain ξi. for the initial values, u0)), i =...n, solve forward in time to compute the density along the trajectories. the two-step procedure is illustrated in figure this way the distribution can be studied on arbitrary regions of the state and parameter space, and no subsequent normalization is required. we used the two-step procedure to obtain an improved resolution of the heavy tail ) and observe the formation of a second mode close to the origin. as an interpretation, this may imply that for this part of the population, x does not reach a certain threshold concentration within the given time interval. in other applications, such as toxicological risk assessment studies, the information may analogously be used to determine the percentage of a population that exceeds or remains below a certain toxicological threshold. to illustrate that the method of characteristics can be applied to any continuously differentiable pdf, not only normal pdfs, we repeated the above computations for an initial exponential distribution with mean μ = &2). with standard monte carlo-based sensitivity approaches such localized information is difficult to obtain, when the region of interest has only low probability. the heavy tail in figure has a total probability of approximately, which means that in expectation only % of the monte carlo sampling points will lie in the interval. with the two-step procedure we used discretization points to approximate the heavy tail. compared to our approach, it would require monte carlo sample points to expect the same coverage on the interval, and the subsequent step of density estimation required by monte carlo methods further impacts the approximation quality. apart from variability in the initial concentration x, we can additionally account for variability in the parameter values. we computed the density for a state space extended according to by the cooperativity β, by the maximal expression rate vmax and by both β and vmax, i.e., with extended state space variables',' and', respectively. the initial distribution was assumed to be a joint normal distribution, where x had mean and variance as before, and the means of vmax and β were set to and respectively, each with variance. using the above two-step procedure, we computed the densities at t = figure depicts the marginal distributions of protein concentration x under the different scenarios of variable/uncertain parameters obtained by subsequent integration over the parameter dimensions that are not considered. for example, the marginal density in x is obtained from the joint density in by u=∫udβ. numerically, we discretized the above integral using the midpoint rule. comparison of the distributions indicates that variability in the cooperativity β has minor impact on the final variability of the protein concentration x. the corresponding joint distributions of' and' are shown in figure as well as the two-dimensional marginal distributions for the three-dimensional case obtained by integrating only β or vmax, respectively. same as in the one-dimensional marginal distribution, it can be seen that the variability in protein concentration is mainly dominated by the variability in the maximal expression rate vmax. the numerical integration ), necessary to visualize multivariate densities for d ≥ or to compute observables such as mean and variance, is currently the computationally limiting step in the application of the method of characteristics, since it requires a regular discretization of the state and parameter space, which becomes prohibitive in high dimensions. for low- and moderate-dimensional systems of odes, however, the method of characteristics provides a more efficient and equally simple alternative to conventional approaches for studying uncertainty and variability of ode systems: the density information obtained via the method of characteristics is expectedly more accurate than estimates obtained with monte carlo methods and considerably richer than a simple characterization of the variability/uncertainty by means of certain indicators. it moreover facilitates further statistical analysis such as model assessment and parameter estimation, e.g., by means of information theoretical approaches, as illustrated later in example for the case of parameter estimation. example consider the genetic toggle switch, where two proteins x and x mutually repress the other protein's expression. in, the concentrations x and x of x and x respectively, are modeled by the two-dimensional system x˙1=α11+β1−kd⋅x1x˙2=α21+β2−kd⋅x the parameter α > represents the effective expression rate of protein x and β > describes the repression cooperativity of the promoter that regulates the expression of x by x analogously, α and β describe the effective expression rate of protein x and its promoter's cooperativity of repression by x respectively. the parameter q corresponds to the concentration of x that represses the promoter activity of x by 50%. the parameter kd again denotes the dilution rate constant. assume that the concentrations x and x at t = have a joint normal distribution with mean μ =' and covariance matrix Σ = diag{, }. the initial density u is shown in figure along with the vector field defined by with symmetric parameters α = α = α =, β = β = β =, q = q = and kd =. we computed the density u solving with f as in. at t = the majority of the probability is concentrated on the slow manifold of the vector field with large variance along the manifold ). such steep distributions on lower-dimensional manifolds pose problems to many other methods: methods based on the estimation of mean and co-variance fail to describe the true shape of the distribution. most pde solvers have numerical problems with such steep gradients. monte carlo-based density estimation may yield a too coarse-grained approximation of the true density if knowledge of the manifold is not provided a priori; that is, to arrange the bins of a histogram or for kernel density estimation the centers of the kernel functions along the manifold. using the method of characteristics we avoid these problems: the steep gradients do not pose any problems to the independently computed trajectories, which directly yield the density values on the attractor manifold. the method thus does not require problem specific ingenuity. deterministic models in a likelihood setting for comparison with experimental data so far we have discussed the use of the method of characteristics to study the sensitivity of ode models. it however also offers benefits when comparing model output with experimental data for model assessment, such as validation/falsification and selection between different models, or parameter estimation. an exact match of the deterministic model with the data is unlikely, and a quantification of the mismatch remains a critical issue. some numerical approaches are based on verifying that the experimental data lies within regions of the state space that are reachable with certain parameter sets of the ode model. most other approaches assess the mismatch based on root mean square deviations and select a model or parameters based on a minimization of these errors. such least squares approaches are based on the assumption that deviations are due to additive gaussian noise, usually assumed to be identical and independent at all points in time. as we have seen in the previous examples, a normal distribution can typically not be expected for general, nonlinear ode models, even if the initial distribution is assumed to be normal. in addition, classical least squares approaches do not allow taking into account prior information on the initial condition or on the parameter to be estimated. stochastic approaches offer a natural way of assessing deviations of the model output from data based on the likelihood function. given a set of data points d = {ξ...,ξn }, the likelihood of a model is defined as the probability that the model predicts this data. for continuous state space models m, the likelihood ℒ given d is defined via the pdf by the product of its values at the data points ℒ=∏i=1nu, where u denotes the density function of the output distribution of the model m. based on the likelihood, there are many methods available for model assessment and parameter estimation. we can directly use the method of characteristics to efficiently compute the likelihood of an ode model for given data points. we first consider the case that parameters are known and only initial conditions are affected by uncertainty. given prior knowledge of this uncertainty defined by the initial density u and further given data points d = {ξ...,ξn} at a time t > the likelihood for each single data point can be computed analogously to the two-step procedure in example 1: solve the ode for each of the n data points ξi ∈ ℝd, i =...,n, backward in time from t = t to to obtain ξi. for the initial values, u0)) ∈ ℝd+ i =...n, solve the extended ode forward in time to obtain the likelihood values for each data point u). in accordance with, the likelihood ℒ of the model given the data is the product of the likelihood for each single data point, i.e., ℒ=∏i=1nu). as the number n of data points is in most cases comparatively small, the above two-step procedure is an efficient way to obtain the likelihood of an ode model. example reconsider the ode model of autoregulated gene expression from example assume that we want to estimate the maximal expression rate vmax of the protein x based on an observation ξ = at time t = for sake of clarity we only consider one data point. for several data points the same procedure as described below applies, and the final likelihood is given by the product of the single-data likelihood values. as we are interested in the likelihood of different values of vmax, we consider the autoregulation model extended according to by vmax. we apply the above two-step procedure to a representative ensemble {υ...,υn} of values of vmax. for each pair, υi) the backward-solution yields a different value, υi). given prior knowledge in terms of a joint pdf u for x and vmax, the forward-solution of with initial conditions, υi), u υi)) then yields the likelihood values u, υi)) associated with each vi, i =...,n. we computed the likelihood of a set of equidistant values of vmax ∈ using the same parameter values as in example for two different scenarios of prior information: a joint normal distribution of x and vmax with parameters μ =' and Σ = diag{, }, and a joint distribution of x and vmax, where x is normally distributed with n and vmax is independently uniformly distributed on the interval. the first scenario accounts for prior knowledge of x and vmax, where a more or less precise knowledge of vmax is given is small). accordingly, the maximum-likelihood estimate is vmax*≈ close to the prior mean of vmax. in the second scenario no prior information on vmax was imposed. the maximum-likelihood estimate of vmax is therefore solely determined by the value of vmax that yields the initial value closest to μ = since the data point ξ = is relatively unlikely for larger vmax ), scenario yields a smaller maximum-likelihood estimate of vmax*≈.  studying the effects of uncertainty and variability in initial values or parameters of ode models can be computationally intensive, since it generally involves solving the system a large number of times. the method of characteristics offers a simple yet accurate alternative to conventional approaches for small- and moderate-dimensional systems. the approach does not assume a particular shape of neither input nor output distribution, it only requires the pdf to be sufficiently smooth and yields density values that are exact up the accuracy of the ode solver used. our first two examples illustrate how a precise characterization of the model uncertainty/variability can be obtained with only few trajectories. in this context we also demonstrated that the analysis can be efficiently restricted to certain sub-regions of the state/parameter space. one limitation of the two-step procedure used for the latter analysis is that for chaotic models the backward-forward solution of the ode is ill-conditioned. in such cases one may resort to approximate techniques as the ukf, but care must be taken that the assumption of the kalman filter, that the distribution remains approximately normal, is satiesfied. another limitation, currently the main one, is the need for uniform grids when dimensions are to be integrated from the final density. but we anticipate that the method of characteristics will prove useful in the context of error control for approximate solution methods of eq. or such as monte carlo or the unscented kalman filter in higher dimensions by providing exact values of the pdf at particular points in state space. as another application we considered the comparison of model results with experimental data. for deterministic models numerical approaches typically rely on root mean squared errors to quantify deviations. their minimization can be interpreted as the maximum-likelihood estimate based on the assumption that deviations are normally distributed at all times. while being a simple assumption, for general nonlinear ode models it is rarely expected to hold. in the third example we described a simple framework, where the method of characteristics was applied to maximum-likelihood parameter estimation based on a distribution that accounts for prior knowledge of parameters and initial values and for the system dynamics. we provide matlab files illustrating the method of characteristics in additional file authors' contributions ayw, wh and rhm planned and performed the research, ayw performed the numerical simulations, all authors contributed to the design and the writing of the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 matlab files illustrating the usage of the method of characteristics. click here for file acknowledgements this work was supported in part by the science foundation ireland, grant sfi07/rpr/i we thank ken r. duffy for his critical reading of the manuscript.
both, uncertainty and variability can be quantified in terms of a probability density function on the state and parameter space. in many applications, ordinary differential equation models are subject to uncertainty or variability in initial conditions and parameters. the partial differential equation that describes the evolution of this probability density function has a form that is particularly amenable to application of the well-known method of characteristics. this leads to simple methods for studying uncertainty, variability and likelihood, with significant advantages over more traditional monte carlo and related approaches especially when studying regions with low probability. the value of the density at some point in time is directly accessible by the solution of the original ode extended by a single extra dimension.
 in many applications, ordinary differential equation models are subject to uncertainty or variability in initial conditions and parameters. both, uncertainty and variability can be quantified in terms of a probability density function on the state and parameter space.  the partial differential equation that describes the evolution of this probability density function has a form that is particularly amenable to application of the well-known method of characteristics. the value of the density at some point in time is directly accessible by the solution of the original ode extended by a single extra dimension. this leads to simple methods for studying uncertainty, variability and likelihood, with significant advantages over more traditional monte carlo and related approaches especially when studying regions with low probability.  while such approaches based on the method of characteristics are common practice in other disciplines, their advantages for the study of biological systems have so far remained unrecognized. several examples illustrate performance and accuracy of the approach and its limitations. 