Email updates

Keep up to date with the latest news and content from BMC Systems Biology and BioMed Central.

Open Access Highly Accessed Methodology article

A continuous-time adaptive particle filter for estimations under measurement time uncertainties with an application to a plasma-leucine mixed effects model

Annette Krengel1*, Jan Hauth1, Marja-Riitta Taskinen2, Martin Adiels34 and Mats Jirstrand5

Author affiliations

1 Fraunhofer-Institut für Techno- und Wirtschaftsmathematik (ITWM; Fraunhofer Institute for Industrial Mathematics), Kaiserslautern, Germany

2 Department of Medicine, University of Helsinki, Helsinki, Finland

3 Department of Mathematical Sciences, University of Gothenburg, Gothenburg, Sweden

4 Wallenberg laboratory, Sahlgrenska Center for Cardiovascular Research, Department of Medicine, University of Gothenburg, Gothenburg, Sweden

5 Fraunhofer Chalmers Centre (FCC), Gothenburg, Sweden

For all author emails, please log on.

Citation and License

BMC Systems Biology 2013, 7:8  doi:10.1186/1752-0509-7-8


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1752-0509/7/8


Received:31 October 2011
Accepted:7 January 2013
Published:19 January 2013

© 2013 Krengel et al.; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

When mathematical modelling is applied to many different application areas, a common task is the estimation of states and parameters based on measurements. With this kind of inference making, uncertainties in the time when the measurements have been taken are often neglected, but especially in applications taken from the life sciences, this kind of errors can considerably influence the estimation results. As an example in the context of personalized medicine, the model-based assessment of the effectiveness of drugs is becoming to play an important role. Systems biology may help here by providing good pharmacokinetic and pharmacodynamic (PK/PD) models. Inference on these systems based on data gained from clinical studies with several patient groups becomes a major challenge. Particle filters are a promising approach to tackle these difficulties but are by itself not ready to handle uncertainties in measurement times.

Results

In this article, we describe a variant of the standard particle filter (PF) algorithm which allows state and parameter estimation with the inclusion of measurement time uncertainties (MTU). The modified particle filter, which we call MTU-PF, also allows the application of an adaptive stepsize choice in the time-continuous case to avoid degeneracy problems. The modification is based on the model assumption of uncertain measurement times. While the assumption of randomness in the measurements themselves is common, the corresponding measurement times are generally taken as deterministic and exactly known. Especially in cases where the data are gained from measurements on blood or tissue samples, a relatively high uncertainty in the true measurement time seems to be a natural assumption. Our method is appropriate in cases where relatively few data are used from a relatively large number of groups or individuals, which introduce mixed effects in the model. This is a typical setting of clinical studies. We demonstrate the method on a small artificial example and apply it to a mixed effects model of plasma-leucine kinetics with data from a clinical study which included 34 patients.

Conclusions

Comparisons of our MTU-PF with the standard PF and with an alternative Maximum Likelihood estimation method on the small artificial example clearly show that the MTU-PF obtains better estimations. Considering the application to the data from the clinical study, the MTU-PF shows a similar performance with respect to the quality of estimated parameters compared with the standard particle filter, but besides that, the MTU algorithm shows to be less prone to degeneration than the standard particle filter.

Keywords:
Particle filter; Sequential Monte Carlo methods; Nonlinear filtering; Parameter estimation; Measurement time uncertainties; PK/PD; Mixed effects; Leucine kinetics

Background

Measurement time uncertainties

Uncertainty in the time at which a measurement is taken is an often neglected source of random error. While in many application areas, this kind of error is generally small and indeed neglectable (due to automated measurements and precise timings), in others it may be of real influence, especially in the life sciences. As a prominent example, one may consider pharmacokinetic and pharmacodynamic (PK/PD) models which are used to describe the metabolic interactions and the effects of a chemical agent (like a drug or a labelled substance) over time inside an organism, respectively.

A typical population experiment in the PK/PD context consists in the analysis of the contents of the blood plasma of several individuals with respect to concentrations of certain molecules of interest. For this purpose, blood probes have to be taken from each individual at certain (fixed) time points after a certain event has occurred (e.g. a drug or a labelled substance has been applied). It is clear from the setting of the experiments that there is some variation in the real point in time when the blood probe has been taken: the true time when the measurement value has been obtained might be shortly before or after the intended time, and this true measurement time is not known to us. Since the inclusion of those time uncertainties in the model usually makes the analysis more difficult, it is standard to lump the time uncertainties with the measurement error. But especially at early times when concentrations change quickly, this may easily lead to wrong estimations, even if one assumes very high variances of the measurement error (we will demonstrate this later on a simple example). On the other hand, the inclusion of measurement time uncertainties (MTU) in algorithms aiming at inference making in complex models is not straightforward. In this article, we will present a modification of the Particle Filter (PF) algorithm (which we call MTU-PF) which is able to fully include a statistical model of the time uncertainties.

Inference in complex systems

The assessment of the effectiveness of a drug in a clinical study has been done in the past by the direct computation of relatively simple statistical values. The enormous increase in complexity of the underlying models, due to present developments in medicine and biology, for instance in the areas of personalized medicine or systems biology, increases also the need for more sophisticated model-based inference methods.

The estimation of unobservable internal variables or model parameters from data which have been obtained from blood or tissue samples at several time points can reveal information on the concentrations and effectiveness of the substance under question. If these data come from individuals which belong to two different (or even more) groups, e.g. test and control group, mixed effects are introduced in the underlying models. The inherent non-linearity and high variability of biological processes adds considerably to the difficulties one faces during the inference step. Inference in connection with dynamic models plays a major role in many other application areas. State and parameter estimation as well as model discrimination and validation are most common, but also optimal control problems should be mentioned.

It is often not enough to consider (independent) measurement noise [1]. Correlations between residuals are not uncommon, and the violation of this statistical assumption may lead to wrong estimates. A natural way to include correlated noise is to model two different types of noise: the dynamic (process or system) noise which is present in the dynamics of the system states and originates either from true random fluctuations in the system or from unmodelled dynamics in the system, and the measurement noise which is introduced by the measurement procedure or equipment and modelled by independent residuals. One possible approach is to use state space models which consist of a time-continuous model for the system states, e.g. based on Stochastic Differential Equations (SDEs), and a separate model for the time-discrete measurements.

Parameter estimation with Maximum Likelihood approach

Parameter estimation in state space systems is a difficult problem. In a context where the system dynamics are modelled by Ordinary Differential Equations (ODEs) without correlated noise, the problem is most often considered as a (deterministic) optimization problem based on a Maximum Likelihood (ML) formulation. An overview of these approaches can be found in [2] and [3]; see also [4], which consider also other aspects like identifiability. A generalization of the ML approach including more flexible cost functions is given by the prediction error estimation method ([5]). The introduction of system noise in the state variables leads to optimization problems with SDE constraints. In this case, internal system states which cannot be directly observed need to be estimated jointly with the parameters, given the data. For this purpose, the parameter estimation methods must be augmented by appropriate state filtering methods. An overview of ML parameter estimation in these types of models is given in [6]. If the SDEs are non-linear, linearizations to the Kalman Filter, like the Extended Kalman Filter (EKF) or the Unscented Kalman Filter (UKF), are used to establish approximations to the means and covariances of the filter distributions over time. All those approximations suffer from the fact that they approximate the filtering distributions of the states by a Gaussian distribution at all time points and cannot adequately approximate skewed or multimodal distributions. Better approximations are provided by simulation based methods like Sequential Monte Carlo (SMC) algorithms where good convergence results have been established ([7]). Nevertheless, they suffer from several drawbacks when applied to the joint estimation of dynamic states and fixed parameters ([8-10], see also [11]).

Parameter estimation in a Bayesian context

In a Bayesian context, in contrast to the “classical” ML approach, a prior distribution is assigned to the parameter vector, hence the parameters can be treated as random variables. In this sense, parameter estimation is done by evaluating the so-called posterior distribution which can be computed (at least theoretically) by Bayes’ theorem given the observations (measurements) and the prior distribution. In the context of high-dimensional spaces, this requires the computation of high-dimensional integrals which is not possible to do analytically. For this purpose, Markov Chain Monte Carlo (MCMC) methods provide powerful tools for the computation of simulation-based approximations to the posterior distribution. Again, in the context of the joint estimation of dynamic states and fixed parameters, the design of good proposal densities is a very difficult problem which renders the use of standard MCMC methods like the Metropolis-Hastings sampler impractical for the purposes of parameter estimation in state space systems.

It has long been a wish to combine both (dynamic) SMC and (static) MCMC methods to provide a general tool for the joint estimation of dynamic states and static parameters. Only recently, Andrieu et al. [11] proposed a very promising combination of both types of Monte Carlo approaches called Particle Markov Chain Monte Carlo (PMCMC) which is generally applicable and where also convergence has been proved.

In the present article, even though the PMCMC approach might be the preferred method for parameter estimation in state space systems, we will concentrate solely on the SMC methods, since our modification affects only this part. However, to be able to do parameter estimation in a pure SMC context, we rely on an approach that is very often used to avoid problems with the estimation of constant parameters. This approach consists in the introduction of artificial dynamics in the parameters, that means the parameters are allowed to slightly change their values over time. In this way, and in a Bayesian context, the parameters can be treated exactly in the same way as the system states. After building an augmented system state by concatenating the parameter vector and the state vector, the joint estimation of states and parameters reduces to filtering of the augmented state vector which makes SMC methods directly applicable to the problem.

Particle filters for state and parameter estimation

Particle filters ([12-14]) belong to the class of SMC methods for state filtering in state space models. Using the state augmentation approach, the method is also capable of estimating system parameters. The standard particle filter is designed for discrete, non-linear, and non-Gaussian models and can routinely be adapted to the continuous case with measurements at discrete times. The idea of the particle filter is that, at each time point, there is a sample based representation (the weighted particles) of the current estimate of the inner states and parameters which is based on the measurements that have been obtained up to the current time point. The particle cloud is then propagated through time, and the particles and weights are updated accordingly at each time point where measurements are available.

Non-Linear Mixed Effects models

Estimation in a Non-linear Mixed Effects model (NLME) involves the estimation of both global and individual parameters. With classical maximum likelihood estimation, the individual parameters are random variables equipped with a distribution while the global parameters remain constants with a “true” but unknown value. If the underlying model equations are non-linear, this leads to likelihood functions which are not analytically accessible and one has to rely on approximations. In the context where the system dynamics are modelled by ODEs, the most popular algorithm for NLME parameter estimation in the PK/PD context is the tool NONMEM ([15]). In [1] an estimation algorithm for NLME models based on Stochastic Differential Equations (SDEs) was proposed that uses the First-Order Conditional Estimation (FOCE) method to approximate the likelihood in combination with the EKF estimation in the SDEs. This has been added to NONMEM ([16]). In [17], a comparison between ODE and SDE based parameter estimation has been performed which showed that the interindividual variabilities were in general estimated to be smaller for the SDE model. Donnet and Samson ([18]) proposed a stochastic version of the Expectation-Maximization (SAEM) algorithm (for the estimation of the global parameters) in combination with MCMC methods (for the estimation of states and individual parameters). However, since MCMC exhibits slow mixing properties in the context of the estimation of states and parameters in state space models, in [19] MCMC has been replaced by the more promising PMCMC approach of Andrieu et al. ([11]).

On the other hand, in a Bayesian context, also the global parameters are equipped with a (prior) probability distribution, and the conceptual difference between global and individual parameters vanishes. The mixed effects model can then be considered as a hierarchical model with dependent parameters ([20,21], see also [22] for a more recent population-based Bayesian approach to PK/PD modelling). Simulation-based (Monte Carlo) methods can easily be adapted to this case. Nevertheless, the above mentioned challenges to both SMC and MCMC methods are even higher due to the increased number of states and parameters in NLME models (the number of states and individual parameters has to be multiplied by the number of individuals).

Aim of the article

Our goal is two-fold: Firstly, we want to show that the particle filter algorithm is applicable (with our modifications) also to more complex models when time uncertainties are formulated explicitly. Secondly, we want to show that the modification may even provide the possibility for further enhancement of the performance of the algorithm by presenting an adaptive time-stepping scheme which is only possible in the context of the new algorithm.

We do not claim that our MTU algorithm generally performs better or worse than the standard filter, nor that it should be the preferred method for estimation in non-linear mixed effects models. Rather, we provide a method which is usable for models where time uncertainties may play a major role. In these cases, it may indeed lead to better estimations. On the other hand, our method transfers the time-discrete particle filter approach, where updates based on the measurements very strictly depend on the measurement times, to a truly time-continuous approach, where updates to the filtering distributions can be performed at every point on the time-scale. Since we want to focus on the time uncertainties, we neglect discussing further issues like identifiability, model evaluation and model discrimination. In our application to the model of plasma-leucine kinetics, we try to avoid these issues by providing ad-hoc values to some of the parameters (especially to the variances of the system states).

Motivating example

Let us have a look at an example for illustrating the benefits of a separate modelling of measurement time uncertainties. Let us consider a state space system given by the ODE

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M1">View MathML</a>

with parameters <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M2">View MathML</a>. Here <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M3">View MathML</a> denotes the state of the system at time t. We call the state trajectories obtained by this deterministic system the nominal evolutions of the states. We add noise to the system in a standard way by introducing an additional term <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M4">View MathML</a>, with a standard Wiener process <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M5">View MathML</a> and a diffusion parameter <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M6">View MathML</a>. This leads to the following SDE describing the evolution of the state q over time:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M7">View MathML</a>

Furthermore, let the initial state q(0) be given by a log-normal distribution with parameters <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M8">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M9">View MathML</a> (mean and variance of the logarithm of q(0), respectively). The parameters chosen in our implementation of this example are shown in Table 1.

Table 1. Parameters for the motivating example

We assume that M measurements of the state q(t) will be taken at times Tj and that each measurement j, j=1,…,M is disturbed by normal noise with mean q(tj) and with a fixed variance <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M14">View MathML</a>, i.e. measurement j is distributed according to <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M15">View MathML</a>. Usually, the times Tj are assumed to be known. In contrast, we will assume that in addition to the measurement value error, there is some uncertainty about the exact times where the measurements have been taken. If we attempt to take the jth measurement at the intended (or nominal) time <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M16">View MathML</a>, the measurement in fact takes place at time Tj which may be shortly before or after the intended time <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M17">View MathML</a>. A natural way to model these uncertainties is to assume that the measurement time Tj is given as a realization of a random variable Tj. In our example, we assume that Tj follows a truncated normal distribution given by the density

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M18">View MathML</a>

with normalizing constant

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M19">View MathML</a>

and given intended measurement times <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M20">View MathML</a>. Figure 1 shows the different distributions for one measurement j for all possible intended times <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M21">View MathML</a>. In each of the subfigures (a)-(d), the shaded green area gives an impression of the “density” of the distribution of the measurement value yj in dependence of the intended measurement time <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M22">View MathML</a> on the x-axis, while the dark-green dashed line depicts the nominal evolution of the state q over time. Subfigure (a) shows the distribution of the measurement values with time uncertainties, while (b)-(d) depict the distribution of the measurement values with known measurement times (in this case <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M23">View MathML</a>), for several standard deviations σy: in (b), the original standard deviation is used, while in (c) and (d), higher standard deviations are used which correspond to the cases with lumped value and time variations.

thumbnailFigure 1. Assumed measurement distributions for the motivating example. Measurement distribution resulting from (a) separate modelling of measurement time uncertainties and measurement value uncertainties, with σy=0.005, and (b-d) lumped time-and-value uncertainties with several different assumed lumped measurement variances σy. The dashed dark-green line depicts the nominal evolution of the state q over time. The green shaded area depicts the region where the measurements are expected.

Comparing Figures 1(a) and 1(b-d), we observe that the distributions of the measurements exhibit clearly different shapes. For the “true” model depicted in Figure 1(a), if we consider a single point in time that lies in a time segment where the state values change quickly, the distribution of the measurement at this certain point in time is quite broad. The variance in the measured value is very high, whereas it is small in time segments where the state values change slowly. In contrast, for the standard particle filter, the measurement variance is constant and hence the assumed measurement distributions differ remarkably from the “true” distributions, howsoever the value of <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M24">View MathML</a> is chosen. It must be expected that this leads to difficulties when inference on the states and parameters needs to be done based on these models. We will resume our example after having presented the MTU particle filter and will show that this is indeed the case.

Methods

We divide this section into three subsections. In the first subsection, we fix the state and observation model we want to consider. In the second subsection entitled “Standard case” we outline the standard particle filter algorithm in the context of time-continuous states with time-discrete measurements, and the various probability distributions involved. Although nothing is new in this subsection, it serves several purposes. Firstly, the time-continuous case is relatively rarely considered in the literature; secondly, the derivation of our modification needs a slightly more general formulation than it is standard for the discrete-time filter; and lastly, the comparison of our modified version with the standard case might more clearly reveal the differences between the two approaches. In the third subsection entitled “MTU particle filter”, we present our new modification of the particle filter. In the following section “Results and Discussion”, we compare the new MTU particle filter to the standard particle filter and to an alternative Maximum Likelihood estimation method on a simple artificial example. We also present an application of our MTU-PF method to a PK/PD study in a non-linear mixed-effects setting in direct comparison with the standard particle filter.

Note: a list of all used symbols with a short explanation can be found at the end of this paper.

The model

State process

Let <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M25">View MathML</a> be a probability space and for each t∈ [t0,) with <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M26">View MathML</a> let <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M27">View MathML</a> be an arbitrary measurable space. For each t∈ [t0,) let further <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M28">View MathML</a> be an <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M29">View MathML</a> measurable random variable such that <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M30">View MathML</a> is a continuous-time Markov process with general state space

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M31">View MathML</a>

For each t∈ [t0,), denote by <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M32">View MathML</a> the pushforward measure of P under Xt, i.e. <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M33">View MathML</a> for all <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M34">View MathML</a>. Further, denote by <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M35">View MathML</a> the pushforward measure of P under <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M36">View MathML</a> (with the corresponding product algebra). Analogously, denote by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M37">View MathML</a>

the state space restricted to the interval [t0,t], and denote by <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M38">View MathML</a> the corresponding pushforward measure. For each s and t with t>st0, let Ks,t(xs, dxt) be the Markov kernel of the process <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M39">View MathML</a> from time s to time t.

An important special case for <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M40">View MathML</a> is given by a multidimensional Itô process on <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M41">View MathML</a> (equipped with the corresponding Borel σ-algebra) defined through a stochastic differential equation (SDE)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M42">View MathML</a>

with drift a(x,t), diffusion matrix B(x,t), multidimensional standard Wiener process <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M43">View MathML</a>, and initial variable <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M44">View MathML</a>. In this case, it is possible to sample directly (at least approximately) from the kernels Ks,t when a suitable discretization method is applied, for instance the Euler-Maruyama method.

Observations / measurements

Let the process <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M45">View MathML</a> be observed via M random variables Y1:M with values in measurable spaces <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M46">View MathML</a>. Each single observation (measurement) yj depends on the state variable <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M47">View MathML</a> at some time Tj and on the observation time (measurement time) Tj itself. We assume that, given the observation time Tj and the state <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M48">View MathML</a>, the variable yj is independent of all other variables, and the conditional measure can be expressed via some conditional probability density <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M49">View MathML</a> with respect to a reference measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M50">View MathML</a> on <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M51">View MathML</a>. We do not require any further restrictions on g such as linear dependence on the states or Gaussianity.

Observation / measurement times

The observation times (measurement times) Tj for j=1,…,M are usually assumed to be deterministically given and known. Our variant of the particle filter will be based on the assumption that the observation times Tj are themselves realizations of random variables Tj. These variables model the uncertainty about exact observation times. In contrast to the observation variables yj, the observation times Tj are never observed (measured). We assume that all information available to us is their probability distribution on the half axis [t0,), while in the case of the observations yj, we know both the densities <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M52">View MathML</a>and the observed values yj.

In this article, we will only consider the simplest case where each variable Tj is independent of all others. Dependencies between the Tj’s, especially concerning the order of the observation times, may be considered natural but would lead to more complicated algorithms. However, order dependencies can easily be introduced via restrictions on the support of the variables. In general, the probability distribution of every single variable Tj shall be given by a density γj(tj) with respect to the Lebesgue measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M53">View MathML</a> on the interval [t0,).

In the following, we will consider the two cases mentioned, where either all Tj are deterministic and known or all Tj are random and unknown. Note that the first case formally coincides with the case that Tj is random but observed. We will therefore stick to the notation <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M54">View MathML</a> for the observation densities in both cases.

Standard case: measurement times deterministic and known

We will first consider the standard case, where the observation times Tj are known. For simplicity, we assume here that the observation times t1:M are strictly ordered increasingly, i.e. t0<t1⋯<tM.

The standard case of the particle filter is usually formulated for discrete-time Markov processes <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M55">View MathML</a> with general state space where the state variables are only defined at the initial time t0 and at the times t1,…,tM when measurements occur. Nevertheless, this case is included in our more general framework where Xt is defined for all tt0. One just focuses on the state variables for those times only. In view of the later generalization to random observation times, we will consider the fixed values Tj as realizations of random variables Tj and condition all occurring densities on them. As mentioned above this assumption leads to the same results as if we assumed the values Tj to be given deterministically.

Full model and filter model

The full model is given by the joint density of the variables <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M56">View MathML</a> and Y1:M (conditioned on the observation times T1:M=t1:M) with respect to the product measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M57">View MathML</a>:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M58">View MathML</a>

(1)

The filter at a given time tk is based on a reduced model. This model is given by the joint density of the variables <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M59">View MathML</a> and Y1:k (conditioned on T1:M=t1:M) with respect to the product measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M60">View MathML</a>:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M61">View MathML</a>

(2)

This density is based on the state sequence <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M62">View MathML</a>. In contrast, we can focus on the single state <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M63">View MathML</a> by considering the joint density of the variables <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M64">View MathML</a> and Y1:k (given T1:M=t1:M) with respect to <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M65">View MathML</a>. It can be computed by marginalization as follows:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M66">View MathML</a>

(3)

and the filter density at time tk with respect to <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M67">View MathML</a> can then be computed with Bayes’ theorem:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M68">View MathML</a>

(4)

with

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M69">View MathML</a>

(5)

For general (non-linear) models, the practical computation of the filter density is very difficult. Nevertheless, the particle filter computes a Monte Carlo approximation using the fact that the filter densities <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M70">View MathML</a> can be computed recursively. This is done in two steps. We consider the filter distribution at time tk-1 given by the probabilities

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M71">View MathML</a>

(6)

for each set <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M72">View MathML</a>. We then get first the prediction distribution, i.e. the distribution of <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M73">View MathML</a> given the data Y1:k-1 (and t1:M), by use of the kernel <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M74">View MathML</a>:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M75">View MathML</a>

(7)

for each set <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M76">View MathML</a>. Then we use Bayes’ theorem to get the filter distribution at time tk:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M77">View MathML</a>

(8)

for each set <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M78">View MathML</a>, with normalizing constant

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M79">View MathML</a>

(9)

Importance sampling

Another ingredient for the particle filter is sequential importance sampling. We assume that a second Markov chain <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M80">View MathML</a> on the same state space is given with pushforward measures <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M81','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M81">View MathML</a> and kernels <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M82','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M82">View MathML</a> for j=1,…,M. We assume that for each <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M83','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M83">View MathML</a>, the measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M84','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M84">View MathML</a> is absolutely continuous with respect to the measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M85','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M85">View MathML</a>. It follows that the Radon-Nikodym derivative (written as conditional density)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M86">View MathML</a>

exists. We further assume that the pushforward measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M87">View MathML</a> under <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M88">View MathML</a> is absolutely continuous with respect to the corresponding pushforward measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M89','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M89">View MathML</a> under <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M90','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M90">View MathML</a> with Radon-Nikodym derivative

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M91">View MathML</a>

For sequential importance sampling, we need to be able to sample from the initial measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M92','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M92">View MathML</a> and from the kernels

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M93">View MathML</a>

for each <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M94','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M94">View MathML</a>, and to compute <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M95','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M95">View MathML</a> as well as <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M96">View MathML</a> pointwise.

Using

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M97','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M97">View MathML</a>

we can then write the recursive formula (8) for the filter distribution at time tk as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M98','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M98">View MathML</a>

(10)

for each <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M99','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M99">View MathML</a>. The direct computation of the normalizing constants <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M100','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M100">View MathML</a> (while Y1:M is assumed to be fixed) is not necessary. Sequential importance sampling is performed as follows. Draw a number N of realizations <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M101','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M101">View MathML</a> from <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M102','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M102">View MathML</a> and compute the corresponding unnormalized weights

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M103','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M103">View MathML</a>

Then, for all k=1,…,M, sample realizations <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M104','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M104">View MathML</a> from the kernel <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M105','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M105">View MathML</a> for each i=1,…,N and compute the unnormalized weights

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M106','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M106">View MathML</a>

For suitable integrable functions h (e.g. fulfilling some mild restrictions on how fast h may increase with x, see [23] for details), the expectation of h with respect to the filter density conditioned on the observations Y1:k=y1:k, given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M107','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M107">View MathML</a>

(11)

can then be approximated by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M108','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M108">View MathML</a>

(12)

where N is the number of particles. In fact, it can be shown that as N approaches infinity, these empirical expectations converge to the filter expectations:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M109','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M109">View MathML</a>

(13)

Note that if we can sample from the Markov kernels of <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M110','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M110">View MathML</a>, we can choose <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M111','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M111">View MathML</a> (at least in law), whence <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M112','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M112">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M113','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M113">View MathML</a>. This is a standard choice, but in terms of efficiency of the particle filter algorithm not always the best one. On the other hand, finding a suitable Markov chain <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M114','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M114">View MathML</a> different from <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M115','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M115">View MathML</a> is not an easy task.

Resampling

If the number N of samples through time is fixed, the samples obtained by sequential importance sampling quickly degenerate since most of the normalized weights decrease rapidly towards 0. The degree of degeneracy is often measured by an estimate of the so-called effective sample size (ESS). This estimate at time t is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M116','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M116">View MathML</a>

(14)

where

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M117','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M117">View MathML</a>

(15)

are the normalized weights. It obtains its maximal value N if all weights are equal, and it approaches 1 if the variance of the weights and thus the degree of degeneracy increases. To avoid this degeneration of the samples, a resampling step needs to be done when the ESS drops below a threshold NThreshold (which is usually chosen to be N/2).

Resampling at some time s is based on given non-negative (unnormalized) selection weights <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M118','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M118">View MathML</a> for each particle index i: One repeatedly selects particles with probabilities <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M119','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M119">View MathML</a> given by the normalized selection weights

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M120','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M120">View MathML</a>

(16)

This is multinomial resampling. There exist procedures where each single particle is still selected with probability <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M121','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M121">View MathML</a>, but with reduced overall variance, for instance stratified resampling or systematic resampling which should be preferred [24,25]. In any case, resampling defines a selection function ι:II on the index set I: ={1,…,N}. Resampling is then done by replacing the state samples <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M122','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M122">View MathML</a> by the selected state samples <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M123','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M123">View MathML</a>. Since before selection the probability that the particle i will be chosen is <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M124','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M124">View MathML</a> for each draw, the expected number of times that particle i has been chosen after N draws is <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M125','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M125">View MathML</a>. To correct for the introduced bias, the normalized weight <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M126','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M126">View MathML</a> for each selected particle i needs then to be corrected by replacing it by the weight

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M127','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M127">View MathML</a>

(17)

(using (16)). The necessary correction is therefore achieved if the unnormalized weights <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M128','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M128">View MathML</a> are replaced by the corrected unnormalized weights <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M129','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M129">View MathML</a>.

Note that in the original particle filter, the selection weights <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M130','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M130">View MathML</a> at time s are chosen to be the particle weights (before the replacement), i.e.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M131','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M131">View MathML</a>

such that after the resampling step the unnormalized weights are all equal to 1. Nevertheless, in general their choice is free and may be based on the observations (which is used in the so-called auxiliary particle filter [26]).

Particle filter algorithm

The particle filter computes the state realizations and weights recursively through time. In its standard form, the particle filter can be stated as in algorithm 1.

Note that if one chooses <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M132','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M132">View MathML</a> (in law), then <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M133','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M133">View MathML</a> and the update of the weights simplifies to

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M134','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M134">View MathML</a>

Data Likelihood

Model validation or discrimination is generally based on the data likelihood

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M135','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M135">View MathML</a>

(18)

Algorithm 1 Standard particle filter

for given observations Y1:k. Without resampling, the data likelihood could be approximated by the empirical mean of the unnormalized weights, i.e. by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M137','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M137">View MathML</a>

(19)

because this is the empirical estimate for the above expectation. After a resampling step, this is not valid any longer. Nevertheless, in any case, the data likelihood can be computed recursively by the following estimate of the ratio <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M138','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M138">View MathML</a>:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M139','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M139">View MathML</a>

(20)

with initial estimate <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M140','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M140">View MathML</a> (see e.g. [27]).

MTU particle filter: Uncertain measurement times

We now assume that each observation time Tj is a realization of a random variable Tj. Its distribution is expressed via densities γj with respect to the Lebesgue measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M141','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M141">View MathML</a>. The observation times Tj themselves are not observed.

Full model

The full model in this case will include complete continuous state paths, since the observation times are now distributed over the complete time axis [t0,), and the observations may potentially depend on every state <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M142','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M142">View MathML</a> for tj∈[t0,). Consider therefore the joint density of the variables <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M143','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M143">View MathML</a>, Y1:M and t1:M, with respect to the product measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M144','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M144">View MathML</a>:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M145','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M145">View MathML</a>

(21)

Filter model

The filter at a given time tt0 is again based on a reduced model. This model is given by the joint density of the following variables: <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M146','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M146">View MathML</a>, denoting the state paths until time t; further only those variables yj for which Tjt; and finally t1:M. This density is given with respect to the product measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M147','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M147">View MathML</a> by:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M148','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M148">View MathML</a>

(22)

Note that we cannot use the simple notation of the standard case where for filtering only the first k observations are taken into consideration at time tk, since neither the observations are ordered in time nor the times Tj are fixed in advance. For this reason we have to include all measurements Y1:M also into the filter model. Note that even though we use the complete data Y1:M=y1:M in the notation, only those yj have to be known at time t for which Tjt holds. To avoid confusion, we mark all densities connected to the filter model at time t by a hat superscript (and by the index t).

We will now derive formulas for the filter density. Since we assume that the observation times t1:M are not observed, we use marginalization to get the joint density for <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M149','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M149">View MathML</a> and Y1:M only, which is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M150','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M150">View MathML</a>

(23)

with respect to the product measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M151','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M151">View MathML</a>. We will further simplify this density. If we define

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M152','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M152">View MathML</a>

then

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M153','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M153">View MathML</a>

and further

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M154','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M154">View MathML</a>

(24)

where the last step is possible because the factor indexed by j does not depend on <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M155','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M155">View MathML</a> for jj. For each j, we can split the integration by Tj at the time point t into two parts and get

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M156','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M156">View MathML</a>

where the last step follows from the fact that γj is a probability density and therefore

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M157','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M157">View MathML</a>

holds. Inserting this into (24), we get

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M158','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M158">View MathML</a>

With a further marginalization, we get the joint density of Xt and Y1:M for the filter model,

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M159','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M159">View MathML</a>

(25)

which is with respect to the product measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M160','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M160">View MathML</a>. From this density, we finally can compute the filter density with respect to <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M161','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M161">View MathML</a> by applying Bayes’ theorem:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M162','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M162">View MathML</a>

(26)

where

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M163','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M163">View MathML</a>

(27)

is the data likelihood with respect to the measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M164','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M164">View MathML</a>.

Effective computation of the filter distributions

In the following paragraph, we will show how the densities of the filter distributions given by (26) can be effectively computed. This is the basis for the formulation of our MTU particle filter method.

Let the observations <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M165','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M165">View MathML</a> with <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M166','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M166">View MathML</a> be given. For each time t∈ [t0,) and for each j∈{1,…,M}, we define random variables

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M167','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M167">View MathML</a>

by the following system of ODEs

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M168','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M168">View MathML</a>

(28)

for each ωΩ with initial values

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M169','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M169">View MathML</a>

(29)

and by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M170','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M170">View MathML</a>

(30)

We will show that for each set A in the σ-algebra generated by the variable Xt, it holds that

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M171','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M171">View MathML</a>

(31)

where <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M172','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M172">View MathML</a> is the filter density. That means we can use the processes Wj,t and Wt to compute the filter distributions through time. From this, it follows immediately that we can also compute filter expectations. Indeed, for any real-valued measurable function h on <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M173','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M173">View MathML</a> such that E[|h(Xt)|]<, it holds that the expectation of h(Xt) given Y1:M=y1:M with respect to the filtered state Xt defined by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M174','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M174">View MathML</a>

(32)

is given by the following equation:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M175','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M175">View MathML</a>

(33)

To show our assertion, we consider the processes Wj,t for j=1,…,M. According to (28) and (29), each Wj,t is defined as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M176','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M176">View MathML</a>

(34)

so, with (30),

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M177','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M177">View MathML</a>

(35)

holds. Thus, to prove (31), we have to show that for each set A from the σ-algebra generated by Xt,

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M178','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M178">View MathML</a>

holds (see (25) and (26)). It is enough to show the equality for the numerator, i.e.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M179','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M179">View MathML</a>

since the equality of the denominator follows then immediately from the special case A=Ω and from the fact that

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M180','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M180">View MathML</a>

Using the variable transformation <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M181','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M181">View MathML</a>, we get

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M182','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M182">View MathML</a>

This is what we wanted to show.

Weights

Since for each t and for each j=1,…,M the random variables Wj,t and Wt depend only on the process <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M183','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M183">View MathML</a> until time t, we can define functions <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M184','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M184">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M185','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M185">View MathML</a> by setting for each <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M186','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M186">View MathML</a>:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M187','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M187">View MathML</a>

It follows from (34) that

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M188','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M188">View MathML</a>

(36)

for each j and from (35) that

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M189','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M189">View MathML</a>

(37)

The values of Wt will serve as weights in the MTU particle filter. We will call Wj,t the partial weights. Since in each discretization scheme which is applied to solve the integral in the formula (36) for Wj,t the integrand has to be evaluated, we may run into practical problems if we use it as it is written in the formula. If the density <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M190','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M190">View MathML</a> evaluates to a very small value and if we work with fixed-precision numbers, the subtraction of 1 will result in a value which may be practically equal to -1. If this error accumulates over time, we may end up with wrong values for <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M191','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M191">View MathML</a>. To reduce the computational error, the integral could be split up in the following way:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M192','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M192">View MathML</a>

(38)

where the cumulative distribution function

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M193','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M193">View MathML</a>

(39)

is independent of <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M194','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M194">View MathML</a> and where the part

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M195','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M195">View MathML</a>

(40)

depends on the path <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M196','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M196">View MathML</a>. It is even more convenient to compute <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M197','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M197">View MathML</a> by evaluating the antiderivative of γj, if it is computationally available.

Note that the definition of the filter distribution is dependent on the reference measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M198','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M198">View MathML</a>. A suitable change of this measure may help to further increase the efficiency of the algorithm. This issue still has to be explored.

Resampling

Special attention is needed for the computation of the weights after resampling steps have been applied. As mentioned earlier, resampling at time s is done by randomly generating a selection function ι:II (with index set I={1,…,N}) based on given non-negative (unnormalized) selection weights <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M199','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M199">View MathML</a> for each particle index i. The state samples <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M200','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M200">View MathML</a> have then to be replaced by the selected state samples <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M201','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M201">View MathML</a>, and the unnormalized weights <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M202','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M202">View MathML</a> by the corrected weights <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M203','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M203">View MathML</a>. We assume that resampling steps at times s1,…,s with t0s1⋯<st have occurred, states and weights have been replaced at these times, and within this paragraph, we denote them by <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M204','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M204">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M205','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M205">View MathML</a>, respectively, for each particle i. By definition, we have at t=s

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M206','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M206">View MathML</a>

and

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M207','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M207">View MathML</a>

For each time ts and for each particle i, the corrected weights are then recursively given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M208','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M208">View MathML</a>

(41)

with

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M209','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M209">View MathML</a>

Since the process Wt computes the uncorrected weights <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M210','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M210">View MathML</a>, we have to correct the weights explicitly in the algorithm by dividing them by the cumulative product

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M211','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M211">View MathML</a>

Note that we have to select these products during a resampling step similarly to the selection of the states, i.e. at t=s,

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M212','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M212">View MathML</a>

The MTU particle filter algorithm

A practical MTU particle filter can be obtained by any discretization scheme based on an (arbitrary) time discretization t0 = τ0τ1τD. Similar to the standard particle filter, sampling need not necessarily be done from the state process <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M213','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M213">View MathML</a> directly. One can instead sample from another process <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M214','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M214">View MathML</a>, provided that the Radon-Nikodym derivatives

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M215','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M215">View MathML</a>

for each s,t∈[t0,t] with s<t exist and can be evaluated pointwise. (In fact, it suffices that ϱt | s(xt | xs) exists for all states <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M216','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M216">View MathML</a> which are reachable from some initial state <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M217','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M217">View MathML</a> with <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M218','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M218">View MathML</a>). In the special case that we sample from the Markov kernels of <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M219','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M219">View MathML</a> directly, ϱt | s≡1 for all s<t. Our MTU particle filter is described in algorithm 2. Here we suppress the indices s1,…,s in the notations of states and weights (see last paragraph).

Algorithm 2 MTU particle filter

The algorithm as it is written here has to be enriched with concrete discretization methods for the sampling and update steps. For instance, if the process <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M221','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M221">View MathML</a> is a multidimensional Itô process defined through a stochastic differential equation (SDE)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M222','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M222">View MathML</a>

with drift a(x,t), diffusion matrix B(x,t), multidimensional standard Wiener process <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M223','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M223">View MathML</a>, and initial variable <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M224','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M224">View MathML</a>, then the Euler-Maruyama method can be used for discretization. We set Δτd: =τd-τd-1. The sampling is then done by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M225','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M225">View MathML</a>

(42)

where ηi is a sample from a standard normal distribution (with mean 0 and variance 1). Further, the update step can be done in the simplest case using Euler discretization:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M226','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M226">View MathML</a>

(43)

and

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M227','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M227">View MathML</a>

(44)

Of course, better discretization schemes are possible. As we have already mentioned, if the antiderivative Gj with Gj(t0)=0 (which is in fact the distribution function) of γj is available, then one should rather use

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M228','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M228">View MathML</a>

for the computation of the values <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M229','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M229">View MathML</a>.

Adaptive stepsize

To be able to fully exploit our MTU particle filter method, the discretization stepsize must be chosen appropriately. One simple possibility is to use a very small stepsize throughout the complete procedure. A quite high computation time will result from that. This can be reduced if an adaptive stepsize is chosen. We propose to determine the stepsize Δτd online depending on the ESS estimate. The stepsize should decrease when the ESS drops rapidly, and it should increase again if the ESS estimate changes only marginally. In detail, the following procedure can be applied.

In each step of the algorithm, we obtain an initial guess of the stepsize by a linear interpolation between a maximal stepsize if the ESS had not changed since the last step, and a minimal stepsize if the ESS had dropped by the number N of samples (actually, the maximal difference that can be obtained is N-1). From this initial guess, we compute the increments of the partial weights, and from them we predict the ESS in the next step based on the current stepsize guess. If the difference between this predicted ESS and the current ESS drops by more than a certain relative amount (we use 10%), then a new guess of the stepsize is computed by dividing the current guess by 2. With this new guess, a new predicted ESS is computed, and the test can be applied again. This procedure will be applied iteratively until either the difference between predicted and current ESS drops by less than the prescribed amount, or the stepsize guess falls below a prescribed minimal stepsize. The current stepsize or the minimal stepsize, respectively, is then accepted, and the algorithm proceeds with this stepsize. See algorithm 3 for a formal description of the stepsize determination.

Algorithm 3 Determination of the adaptive stepsize

If we sample from the Markov kernels of <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M231','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M231">View MathML</a> directly, then <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M232','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M232">View MathML</a> and the update of the weights does not depend on the new states <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M233','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M233">View MathML</a>. Hence we only need to compute the weight update and the corresponding ESS estimate until we find an adequate stepsize. In this case, it is not necessary to sample the new states in each iteration which renders the algorithm computationally more effective.

Note that this procedure cannot be performed when the measurement times are fixed (i.e. in the standard particle filter). In this case the ESS does not depend on the stepsize, and a reduction of it will not improve the ESS. The application of the MTU particle filter with distributed measurement times is essential to be able to use this adaptive stepsize procedure.

Data Likelihood

As mentioned earlier for the standard case of the particle filter algorithm, without resampling, the data likelihood could be approximated by the empirical mean of the unnormalized weights (see (19)). In our case, this would be

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M234','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M234">View MathML</a>

After each resampling step, we have to correct this formula. Using the same notation as in the paragraph on resampling, for each time ts and for each particle i, the corrected estimate is given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M235','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M235">View MathML</a>

(45)

with

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M236','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M236">View MathML</a>

(see (41)). This can be seen by considering recursively the correction at the time of the resampling step . Since

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M237','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M237">View MathML</a>

is the expected number of times the particle i has been selected after N draws, we have to correct the weights <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M238','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M238">View MathML</a> in selection step by dividing them by this number. The following computation then proves the correctness of the formula via induction:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M239','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M239">View MathML</a>

We get algorithm 4 for the computation of the empirical estimate of the data likelihood, which needs to be done in parallel to the MTU particle filter algorithm.

Algorithm 4 Estimation of the data likelihood

Offline and online estimation

Two main cases of estimation procedures may be distinguished: the offline estimation procedure used for example for parameter estimation, and the online estimation procedure for control purposes. While in the offline case the measurements are completely available before estimation begins, we know only some of the measurement values at some certain time t during the online procedure. In this latter case, an online computation where all measurement times have been modelled with probability densities with infinite support is impossible, because in this situation all measurement values must already be known at time t0. Online estimation is nevertheless possible if the support of the measurement time densities is finite. The online estimation must then be delayed by the diameter of the respective supports.

Implementation

We have implemented the proposed algorithm in Mathematica as part of a Parameter Estimation Toolbox for Systems Biology developed by the Systems Biology group at the Fraunhofer Chalmers Centre (FCC) in Göteborg (Sweden). Furthermore, we have implemented it also in the statistical computing language R [28]. All figures in this article have been created using this R implementation.

Results and Discussion

Motivating example - results

In this section, we resume our motivating example and use it in a parameter estimation setting to compare our MTU filter to both the standard particle filter and to a state-of-the-art Maximum Likelihood (ML) method which is not based on Monte Carlo techniques. To this aim, we first create virtual “measurements” at four intended measurement times <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M241','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M241">View MathML</a>, j=1,…,4 by simulation runs with our true model (see Table 1). This is done as follows: We simulate a single state path (q(t))t∈[0,10] based on the “true” parameter values α=1 and β=3. Then, for each intended measurement time <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M242','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M242">View MathML</a>, we sample an actual measurement time Tj according to the density γj. Now, for each Tj, we sample a measurement value yj from the distribution <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M243','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M243">View MathML</a>. Figure 2 shows one set of measurement values obtained in this way. The intended measurement times are 0.5, 1, 2, 4 and the measurement values we got from one simulation run are 1.083346, 2.550290, 2.700863, 2.949450. We will use these values in the following estimation runs. Our aim is to estimate the parameter vector θ=(α,β) with the following three differrent estimation procedures and compare the results:

1. ML estimation on “lumped” measurements,

2. Standard PF on “lumped” measurements, and

3. MTU-PF.

thumbnailFigure 2. Measurement distribution with a set of possible measurements. The dashed dark-green line depicts the nominal evolution of the state q over time. The green shaded area depicts the region where the measurements are expected. The yellow circles depict a set of possible measurement values.

Note that in both the standard PF and the ML estimation method, it is not possible to include the uncertainties in the measurement times directly. Rather, we have to lump these uncertainties with the uncertainties in the measurement values by increasing their variances <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M244','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M244">View MathML</a> and to treat the <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M245','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M245">View MathML</a> like true measurement times. In both cases, we tested several values of “lumped” variances <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M246','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M246">View MathML</a>.

We have implemented the standard particle filter and the MTU particle filter in the statistical computing language R [28]. The ML estimations have been done using the Parameter Estimation Toolbox for Systems Biology developed by the Systems Biology group at the Fraunhofer Chalmers Centre (FCC) in Göteborg (Sweden).

ML estimation

The ML estimation is a standard method based on a Maximum Likelihood approach, that is, on the maximization of the likelihood function L(θ;y1:M) with respect to the parameter vector θ. The computation of the likelihood requires the use of estimates of the state mean and covariance matrix which are routinely obtained by the use of a non-linear derivation of the Kalman filter.

Given measurements Y1:M at a total of M time points Tj which are assumed to be fixed and known, the likelihood function can be written as

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M247','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M247">View MathML</a>

(46)

where <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M248','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M248">View MathML</a> is the conditional probability density function of the j-th observation yj given all previous observations y1:j-1 as given in (9), here explicitly based on the given parameter vector θ. As already mentioned, analytical solutions are not generally available. An exception to this is the case that all measurements are conditionally Gaussian, that means that all conditional densities <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M249','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M249">View MathML</a> are Gaussian. In this case, the Kalman filter yields the correct solution. The assumption that the measurements are conditionally Gaussian is commonly used as an approximation even in those cases where this is not true. If one accepts that this assumption gives an approximation which is close enough to the true model, then it is only necessary to propagate the means and covariances of the measurements, since the Gaussian probability distribution is completely characterized by the first two moments. We introduce the notation

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M250','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M250">View MathML</a>

and

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M251','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M251">View MathML</a>

for the prediction of the mean and covariance of the observation variables, respectively. The computation of these values can be achieved by some derivates of the Kalman filter, commonly used are the Extended Kalman Filter (EKF) or the Unscented Kalman Filter (UKF), more exactly those versions of these filters which are time-continuous in the states and time-discrete in the measurements (continuous/discrete EKF and UKF, respectively). The residuals ϵj are then defined by the differences between the measurements yj and their estimations <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M252','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M252">View MathML</a>:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M253','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M253">View MathML</a>

The assumption of normally distributed observations leads to an approximation of the likelihood function of the following form

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M254','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M254">View MathML</a>

where l is the number (dimension) of the observation variables. The negative logarithm of this approximated likelihood function is

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M255','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M255">View MathML</a>

The problem to finding maximum likelihood estimates of the model parameters takes the form of a nonlinear optimization problem:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M256','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M256">View MathML</a>

Roughly, the estimation with this approach is done as follows:

1. Choose an initial guess θ for the parameter vector θ.

2. Based on the current parameter vector θ and for each (intended) measurement time Tj, compute residuals ϵj and estimates of the covariance matrices Rj|j-1 using (an approximation to) the Kalman filter.

3. Compute the likelihood L based on the state estimates ϵj and Rj|j-1.

4. Use one step of some local optimization technique to find an improved parameter vector <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M257','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M257">View MathML</a> that decreases the negative log-likelihood function; replace θ by this improved parameter vector <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M258','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M258">View MathML</a>.

5. Repeat steps 2 to 4 until the parameter vector shows convergence.

The results of several runs of parameter estimations for different values of the lumped measurement variances <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M259','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M259">View MathML</a> are shown in Figure 3. For each <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M260','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M260">View MathML</a>, we performed 100 runs with different initial state and parameter values. Note that the choice of an initial parameter value is mandatory for all approaches based on a local optimization method, and that different choices may lead to different results. On the other hand, in our example, we assume a log-normal prior distribution for the initial states (see Table 1), while the Kalman filter based approximations always assume normally distributed initial state values. Conditioned on the intial state (and current parameter), our example system is Gaussian linear which means that the Kalman filter yields exact estimates. We therefore decided the use the following procedure: Before each run, we sample the starting values for the parameters from the same priors as are available to the particle filter (see Table 2), and we also sample initial values q0 for the states from the correct distributions (see Table 1). Then we start the estimation procedure based on these initial parameter and state samples, and we keep the resulting estimated parameters as one “sample”. We repeat this procedure for all of the 100 runs. In this way, the standard procedure uses exact the same information as the particle filters, with the only exception of the measurement time distributions γj. Figure 3 shows the distributions of the estimated parameter values as box plots.

thumbnailFigure 3. Parametersα andβ estimated by ML estimator. Box plots of estimated parameters of 100 runs each. The medians are depicted with circles. The bottom and top of the box are the 0.25-quantile and the 0.75-quantile, i.e. 50% of the values lie within the box. The whiskers mark the 0.025-quantile and the 0.975-quantile, i.e. 95% of the values lie between the whiskers. The horizontal red lines denote the true parameter values.

Table 2. Choices for the PFs for the motivating example

The estimated parameters of all estimation runs are far away from the true values; even quite large variances <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M263','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M263">View MathML</a> do not improve the performance of the estimations. (We also tested the ML approach on a different data set of 4 measurements which was created without randomness in the measurement times: the ML approach performed much better in this case since the estimated parameters varied around the true values; data not shown).

Estimation with particle filters

Tests with the particle filters (both standard and MTU) have been performed in the following way: First we perform an estimation run where we use these measurements to estimate the parameters α and β of the system equation <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M264','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M264">View MathML</a>. We perform this estimation with the MTU particle filter and with the standard particle filter under the same conditions and with the same seed for the random number generator. Thus the initial distribution of the particles is the same in both cases. For the standard filter, we use different levels of variance <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M265','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M265">View MathML</a> for the measurement noise. Afterwards, the empirical medians of the final parameter distributions are used to perform a simulation run in each case, where we estimate the state q of the system based on the measurements and these fixed parameters. For both runs, we compute an estimate of the data likelihood over time. The choices of stepsizes and artificial parameter dynamics can be found in Table 2 (the values for the artificial parameter dynamics are computed by interpolation in the same way as we will describe later in the application of the MTU-PF to the plasma-leucine model).

Table 3 shows a comparison of the estimated parameters and data log-likelihood values. The estimated parameters (empirical medians) of the MTU particle filter (α=1.012 and β=3.010) are very close to the true values (α=1 and β=3) whereas the estimated values of the standard particle filter are significantly worse. The data log-likelihood in both estimation and simulation run is considerably higher for the MTU filter, compared to any of the standard filter runs.

Table 3. Estimated parameters and data log-likelihood values for the MTU particle filter and the standard particle filter with different standard deviationsσy for the measurement noise

In Figure 4, we show the distributions of the filtered states obtained by simulation runs with the estimated parameters. The violet shaded area in Figure 4(a) depicts the estimated state distributions of the state q for the MTU particle filter with the light-blue line marking its median. The dashed dark-green line depicts the nominal evolution of the state q over time which should be approximated by the filter. Figures 4(b-d) show the estimated state distributions of the state q for the standard particle filter with different assumed lumped measurement variances. As can be seen from these figures, the approximation by the MTU particle filter is very close to the evolution of the state with the true parameters. On the contrary, however we choose the measurement variance in the standard case, the algorithm is not able to adequately approximate the correct state evolution.

thumbnailFigure 4. Filtered state distributions based on simulation runs with estimated parameters for the motivating example. MTU particle filter (a) and standard particle filter with several different assumed lumped measurement variances σy (b-d). Violet shaded area: Filtered state distributions based on empirical quantiles. Solid blue line: Median of filtered states. Dashed dark-green line: nominal evolution of the state q over time.

Figure 5 displays the simulated measurement distributions corresponding to the filtered state distributions based on the simulation runs with estimated parameters. Here we notice again that the MTU particle filter can adapt well to the situation, whereas the measurement distributions which can be realized by the standard filter do not fit well with the data points.

thumbnailFigure 5. Simulated measurements corresponding to the filtered state distributions for the motivating example. MTU particle filter (a) and standard particle filter with several different assumed lumped measurement variances <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M266','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M266">View MathML</a> (b-d). Violet shaded area: Measurement distributions based on empirical quantiles. Solid blue line: Median of simulated measurements. Yellow circles: Measurements.

In addition to the significant improvement in the parameter and state estimation, the MTU particle filter has another benefit. Figure 6 shows the development of the effective sample size estimate during the estimation runs. In the standard particle filter runs where the lumped measurement variance is assumed to be small, the predicted values for states and parameters do not fit well with the measurements. At the measurement times, when the predicted states are compared with the measurements and weighted accordingly, most weights decrease rapidly. This leads to a high variance of the weights and the effective sample size estimate drops severely, indicating that only very few particles effectively contribute to the estimation. This degeneration of the particle cloud is prevented by the MTU particle filter, as can be seen in Figure 6(a). For the standard filter runs where <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M267','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M267">View MathML</a> is assumed to be larger, the ESS estimate does not drop as severely as with a small <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M268','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M268">View MathML</a>. This is due to the fact that, with a measurement variance which is assumed to be very high, more or less all possible simulated measurements are considered to fit well with the true measurements. However, these filters are not able to establish reasonably good estimates of the states and parameters and are clearly outperformed by the MTU particle filter.

thumbnailFigure 6. Estimated effectice sample size during estimation run. MTU particle filter (a) and standard particle filter with several different assumed lumped measurement variances σy (b-d). The dotted blue line denotes the resampling threshold.

Application to the Plasma-Leucine Model with Population Data

In this section, we apply our MTU particle filter algorithm to a leucine kinetics model (Demant et al. [29] based on Cobelli et al. [30]) with data taken from a clinical study on diabetes patients ([31,32], see Additional file 1). We perform Bayesian population-based (i.e. NLME) parameter estimation with this model. The model and the data have been previously used in a maximum likelihood estimation context in [17]. It should be noted that the original ODE model needs to be extended by some kind of stochastic process variability in order to turn it into an SDE model. The approach taken in [17] is different from our approach in the way that in [17] the stochastic fluctuations are assumed to be in the tracee (plasma leucine) while we here assume the variability to be in the tracer (labelled leucine). Both assumptions are plausible; a final decision on the best way to model the process dynamics has not yet been made.

Additional file 1. Data used for estimation and simulation. The data consist of one record for each patient, each record consisting of 4 lines in the following format:

patient id [COMMA] group (control or diabetes) [COMMA] leucine initial value [NEW LINE]measurement 1 time [COMMA] measurement 2 time [COMMA] …[COMMA] measurement np time [NEW LINE]measurement 1 value[COMMA] measurement 2 value [COMMA] …[COMMA] measurement np value [NEW LINE][NEW LINE]

where np is the number of measurements for patient p, time is given in hours, and each measurement value is the measured ratio of isotope labeled leucine.

Format: CSV Size: 9KB Download fileOpen Data

The Leucine model

In [31] (see also the thesis [33]), a new combined multicompartmental model for apolipoprotein B-100 (apoB) and triglyceride metabolism in very low density lipoprotein (VLDL) subfractions has been developed, see Figure 7. VLDL are transporters of triglycerides and cholesterol from the liver to the periphery, and elevated levels are associated with increased risk for cardiovascular events. Each VLDL particle has exactly one apoB molecule attached which makes apoB a suitable marker for triglyceride transport.

thumbnailFigure 7. Multicompartmental model for apolipoprotein B-100 (apoB) and triglyceride (TG) metabolism in very low density lipoprotein (VLDL) subfractions. This multicompartmental model was developed in [31]. Circles depict compartments and arrows depict fluxes between compartments. The model includes separate modules for leucine (yellow) and glycerol (red). The free leucine plasma kinetics is modeled by two pools (3 and 4) and a plasma compartment (1), which interchange materials with an intrahepatic compartment (2). Compartment 2 feeds the apoB synthetic machinery. For glycerol, the plasma compartment (13) is connected to a pooling compartment (12) and feeds TG synthesis, which consists of a fast pathway (14) and a slow pathway (21). The assembly of lipoprotein is modeled by separate delays for apoB (11) and TG (22). The plasma kinetics of apoB and TG is modeled by a four-compartment hydrolysis chain, consisting of compartments 5, 6, 8, and 10 for apoB and compartments 15, 16, 18, and 20 for TG. Compartments 5/15 and 6/16 are associated with VLDL1, together with a slowly decaying compartment 7/17. Compartments 8/18 and 10/20 together with the slowly decaying compartment 9/19 form the VLDL2 module. Lipolysis of TG is modeled to take place in the transfer between two compartments. For details on the full model, see [31,34], or [33]. For our purposes, we use a restricted model consisting of the leucine pool, e.g. compartments 1 to 4, see also Figure 8.

thumbnailFigure 8. Schematic depiction of the restricted model (leucine pool) [[33]]. This scheme is a subscheme of Figure 7. Circles depict compartments. Arrows depict fluxes between compartments and are labelled with the corresponding fractional transfer coefficients. Compartment 1 is the plasma-leucine compartment where the leucine is injected. Compartment 2 is an intrahepatic compartment which is the source for apoB synthesis. Compartments 3 and 4 are body protein pools. The output is at compartment 1. Compartment 11 is a delay compartment used only as output from compartment 2.

The secreted particles become denser and denser as triglycerides are delivered to target organs such as muscles and adipose tissue and the relative protein content is increased. As the density increases, the VLDL becomes an intermediate density lipoprotein (IDL) and finally a low density lipoprotein (LDL).

For our purposes we use only the part of the model concerning the leucine pool (compartments 1-4), see Figure 8. This subsystem was first used for apoB kinetic studies by Demant et al. [29] as an implementation of the work by Cobelli et al. [30]. The output is at compartment 1 and at compartment 2. The influx into compartment 1 will be denoted U1.

The data are obtained from tracer/tracee experiments. Here the tracee (i.e. the concentration we are actually interested in) is the leucine amino acids in the apoB molecules. Additional labelled leucine (tracer) is injected as a bolus infusion. Knowledge about the kinetics (fluxes between compartments) of the tracee can be gained by studying the kinetics of the tracer. In the restricted model four compartments are considered: plasma leucine (1), intra-hepatic leucine (2), and two plasma protein pools (3 and 4). In the full model, additional compartments represent VLDL subfractions (compartments 5-11 in Figure 7). The four-compartment system is linked to compartments 5-11 through compartment 2.

For each compartment i where i=1,…,4, let Qi and Qi denote the mass of the tracee and the tracer, respectively. Similarly, let Ui and Ui denote the input for the tracee and the tracer, respectively. For tracer/tracee experiments, Qi is assumed to be in steady state. If the concentration level of the labelled injection is small compared with the overall concentration levels, and if the model is linear, then approximately

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M269','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M269">View MathML</a>

where <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M270','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M270">View MathML</a>, u(t)=(u1(t),0,0,0)T and

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M271','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M271">View MathML</a>

where kj,i for ij is the transfer coefficient of the tracers from compartment i to compartment j (compartments 0 and 11 are considered to be output compartments), and for each i=1,…,4

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M272','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M272">View MathML</a>

Throughout this paper, the time unit used is hours, all fractional transfer coefficients are given in the unit h-1, and the amount of material in compartments is given in mg. In our model, only k0,1, k1,2, k1,3, k2,1, k3,1, k3,4, k4,3 and k11,2 are assumed to be non-zero, while additionally the following dependencies of the transfer coefficients are assumed to be valid:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M273','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M273">View MathML</a>

The fractional transfer coefficient k11,2 has to be fixed for the system to be identifiable. We set k11,2=0.01h-1, as an estimated average from current results. We build stochastic differential equations (SDEs) from the resulting ordinary differential equations by adding noise terms which are given by standard Wiener processes <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M274','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M274">View MathML</a> multiplied by diffusion parameters σ1,…,σ4, respectively. The leucine pool subsystem which we consider here (compartments 1-4) interacts with the surroundings only via an initial input into compartment 1 at time t=0, an output from compartment 1, and another output from compartment 2 (towards compartment 11). All other flows are processes acting inside the subsystem and hence should follow the principle of mass conservation. Therefore we add the stochastic terms to the ODE system in the following way:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M275','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M275">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M276','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M276">View MathML</a>

We fix the diffusion parameters to be σ1=σ2=σ3=σ4=3. Initial conditions are given by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M277','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M277">View MathML</a>

The patients get a bolus injection and therefore the input u1(t) will be modelled as a delta distribution at time t=0h,

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M278','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M278">View MathML</a>

Practically, this means that only the initial condition q1 is affected by it, and we can replace the initial condition for q1 by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M279','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M279">View MathML</a>

and set u1(t)=0 in the differential equations.

The same differential equations, without noise terms, are assumed to be valid for the states Qi and the input U1 of the tracee:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M280','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M280">View MathML</a>

where <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M281','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M281">View MathML</a>

, u(t)=(u1(t),0,0,0)T. It is assumed that the tracee input U1(t)=U1 is constant but unknown. We will therefore estimate it together with the transfer parameters. Since for the tracee steady state is assumed (i.e. dQi(t)/ dt=0), it is possible to solve those equations for Q1(t), and we get:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M282','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M282">View MathML</a>

The output is a measurement proportional to the ratio between the tracer and the tracee, disturbed by log-normal noise:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M283','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M283">View MathML</a>

where we assume the value of the variance parameter to be <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M284','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M284">View MathML</a> (this is the variance of logξt). The parameter p1 denotes the unknown proportion of plasma leucine that actually is in the plasma. The parameters p1 and U1 are not jointly identifiable, therefore we fix p1=0.65. More details concerning the deterministic model (without noise terms) can be found in [31] and [33]. Note that the stochastic disturbances are not part of the original model, they are rather augmentations of the model used in this article.

The mixed effects model

The model as presented in the last paragraph contains only flux parameters kj,i which are assumed to be the same for every individual. Neither does it account for individual differences between several persons, nor does it account for possible changes of the flux parameters when the persons under consideration are affected by a disease or a treatment. To be able to treat these differences in an appropriate way, we introduce group and patient specific parameters in the model; namely the transfer coefficient k0,1 will be split into a group dependent and a patient dependent part. In this way, we introduce so-called mixed effects into the model. Mixed effects generally increase the difficulty in inference making. In the following test runs, we will use measurement data previously reported as individual data for a total of 34 subjects [31,32]. Among them, 15 patients belonged to the diabetes group, and 19 to the control group. From experiments, it can be observed that the degradation rate k0,1 of plasma leucine is significantly different for people with and without diabetes. We therefore assume that the expected value of k0,1 in each group differs and may be either <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M285','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M285">View MathML</a> or <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M286','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M286">View MathML</a> corresponding to the diabetes group and the control group, respectively. Additionally, we assume that we have patient-dependent random factors ζp modelling the parametric uncertainties among individuals, such that finally

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M287','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M287">View MathML</a>

where all ζp’s are static and independently log-normally distributed:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M288','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M288">View MathML</a>

for p = 1,…,34. As a consequence, each of the states q1,…,q4 has to be considered separately for each patient p. We indicate this by writing <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M289','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M289">View MathML</a>, p=1,…,34.

The aim of the estimation runs is, apart from estimating the remaining parameters, to show that the group dependent parameters <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M290','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M290">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M291','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M291">View MathML</a> are indeed different. We want to apply Bayesian estimation to the parameters. For this reason, we principally have to treat them like state variables in the particle filter. Since the estimation of static variables is problematic with particle filter methods, it is standard to introduce small artificial stochastic dynamics to the parameters consisting of normal increments with decaying variances [35]. The same is true for the static noise parameters ηp which also are to be estimated. Our process Xt is then given as an augmented state vector

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M292','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M292">View MathML</a>

The overall model is thus a non-linear mixed effects model with three levels of effects (parameters), namely global parameters, group dependent parameters (<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M293','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M293">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M294','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M294">View MathML</a>), and personal parameters (ζp). Nevertheless, since the core of the model is linear, i.e. the states q1,…,q4 conditioned on all parameters, a Rao-Blackwellization concerning the linear parts of the model is possible, and the Kalman filter applied to this linear partial model can be used in combination with particle filtering for the non-linear parts [36]. Since the model is used for demonstration purposes only, we did not use this technique although in principle it is possible.

Estimation runs

Estimation and simulation runs have been performed with data from all 34 patients (19 from control group and 15 from diabetes group). The computer experiments have been done as follows. We first estimate parameters with the MTU particle filter. Separately, we estimate parameters with the standard particle filter under the same conditions and with the same seed for the random number generator. The initial distribution of the particles is then the same in both cases. For both runs, we compute an estimate of the effective sample size (ESS) and the data likelihood over time. Both estimates allow a performance comparison of the MTU versus the standard particle filter. Secondly, the empirical medians of the final parameter distributions are afterwards used to perform simulation runs in both cases. Both versions of the particle filter, this time with parameters fixed to the estimated values, are used to perform these simulations. In these simulation runs, the data are used for the computation of the data likelihood conditioned on the estimated parameters. The resulting distributions of the simulated measurements can be compared to the true measurements, both visually and by observing the data likelihood.

We used 10,000 particles and a resampling threshold of 7,500. Stepsizes in the MTU filter are between 10-7h and 10-3h, adaptively computed based on the ESS estimate. In the standard filter, we use a fixed stepsize of 10-3h. The data contain measurements until time t=8h, but we perform estimations and simulations only until time t=1h, mainly to reduce computation time; anyhow, after t=1h, the tracer concentrations are relatively low and do not change considerably, and therefore the measurements are not expected to improve the estimations significantly. In our implementation, we directly sample from the states Xt (i.e. <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M295','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M295">View MathML</a> in law) using the Euler-Maruyama scheme for discretization [37].

As mentioned earlier, Bayesian parameter estimation with particle filters requires the introduction of artificial dynamic noise for the parameters. It is standard to use normal increments with decaying variances [35]. Since in our case all parameters with exception of the ηp’s are assumed to be positive, the application of a “log-normal” noise (based on a geometric Brownian motion) in place of the standard normal noise seems to be more appropriate for these parameters. In detail, it has been done as follows. The priors and the distributions of the artificial noise are chosen to be log-normal for all parameters with exception of the individual parameters ηp which have normal priors and noise. The prior for the parameters <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M296','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M296">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M297','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M297">View MathML</a>, k1,2, k1,3, k3,1, k4,3 is <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M298','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M298">View MathML</a>, the prior for U1 is <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M299','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M299">View MathML</a>. The prior for each ηp is <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M300','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M300">View MathML</a>. The parameter update (“artificial noise”) is generally performed according to

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M301','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M301">View MathML</a>

for <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M302','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M302">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M303','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M303">View MathML</a>, k1,2, k1,3, k3,1, k4,3, U1, and according to

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M304','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M304">View MathML</a>

for θ=ηp. It is standard to decrease the variance of the artificial noise over time. In our case we have chosen the diffusion parameter σθ of each artificial noise variable to be dependent on time via a quadratic function

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M305','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M305">View MathML</a>

with parameter dependent coefficients aθ and bθ. Practically, aθ and bθ have been determined by fixing two interpolation points (t0,σθ(t0)) and (t1,σθ(t1)). It holds:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M306','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M306">View MathML</a>

We found best performance with the following choices: For all parameters, we have chosen t0=0h, t1=2h and σθ(t1)=σθ(t0)/10 which means that the diffusion parameter has dropped to 10% of its initial value at time t=2h. The initial values are σθ(t0)=0.5h for all parameters with exception of the ηp’s which have higher initial diffusion σθ(t0)=1h.

As mentioned earlier, we have performed two different estimation and simulation runs, one with the MTU particle filter with distributed measurement times, and for comparison one run with the standard particle filter. In the MTU particle filter, the distributions of the measurement times are truncated normal distributions with mean equal to the nominal value of the measurement time and with variance 0.0012. The normal distribution is truncated at the time point 0.01h left and right from the mean value. In Figure 9, the development of the estimated effective sample size and the estimated data likelihood is shown, both with respect to time t. The development of predictions of the measurements during a simulation run with the final estimated parameters is shown in Figures 10 and 11 for the run with the standard particle filter, and in Figures 12 and 13 for the run with the MTU particle filter. Finally, in Figures 14 and 15, box plots of the final estimated global and individual parameters are shown, respectively.

thumbnailFigure 9. Development of estimated effective sample size and estimated data likelihood during estimation runs. Standard particle filter (top) and MTU particle filter (bottom).

thumbnailFigure 10. Predicted measurement distribution over time (standard particle filter) during simulation run. Patients from diabetes group. Development of predicted measurement distributions over time during simulation run with parameters estimated by the standard particle filter. Circles: Measurements. Solid line: Median of simulated measurements. Violet shaded area: Measurement distributions based on empirical quantiles.

thumbnailFigure 11. Predicted measurement distribution over time (standard particle filter) during simulation run. Patients from control group. Development of predicted measurement distributions over time during simulation run with parameters estimated by the standard particle filter. Circles: Measurements. Solid line: Median of simulated measurements. Violet shaded area: Measurement distributions based on empirical quantiles.

thumbnailFigure 12. Predicted measurement distribution over time (MTU particle filter) during simulation run. Patients from diabetes group. Development of predicted measurement distributions over time during simulation run with parameters estimated by the MTU particle filter. Circles: Measurements. Solid line: Median of simulated measurements. Violet shaded area: Measurement distributions based on empirical quantiles.

thumbnailFigure 13. Predicted measurement distribution over time (MTU particle filter) during simulation run. Patients from control group. Development of predicted measurement distributions over time during simulation run with parameters estimated by the MTU particle filter. Circles: Measurements. Solid line: Median of simulated measurements. Violet shaded area: Measurement distributions based on empirical quantiles.

thumbnailFigure 14. Estimated global and group parameters. Box plots of estimated distributions for the global parameters and the group dependent parameters. The medians are depicted with triangles (standard particle filter) or circles (MTU particle filter). The bottom and top of the box are the 0.25-quantile and the 0.75-quantile, i.e. 50% of the values lie within the box. The whiskers mark the 0.025-quantile and the 0.975-quantile, i.e. 95% of the values lie between the whiskers. The values for U1 have been scaled by a factor of 0.01 in order to fit into the plot.

thumbnailFigure 15. Estimated individual parameters. Box plots of estimated distributions for the individual parameters. The medians are depicted with triangles (standard particle filter) or circles (MTU particle filter). The bottom and top of the box are the 0.25-quantile and the 0.75-quantile, i.e. 50% of the values lie within the box. The whiskers mark the 0.025-quantile and the 0.975-quantile, i.e. 95% of the values lie between the whiskers.

A comparison of the results of MTU-PF and standard particle filter shows that both algorithms exhibit very similar performance with respect to the quality of the estimated parameters, since in both cases the development of the data likelihood is very similar both for estimation and simulation. The estimated log-likelihood of the data at the final time is 137.239 for the MTU particle filter and 136.207 in the standard case, which is practically equal. The computation time of the MTU-PF is only slightly higher than the one of the standard filter. Visual inspection of the simulation runs shows that the simulated measurements of the model with parameters estimated by both filters fit the data in a similar way. This impression is supported by the values of the estimated data likelihood. The MTU particle filter has a final log-likelihood value of 157.622, similar to the one of the standard filter with a value of 155.952. The difference is insignificant.

In contrast to the insignificant differences concerning the resulting likelihoods between the MTU-PF and the standard PF, the development of the ESS estimate in the estimation runs differs remarkably. With the MTU particle filter, the ESS estimate shows a high value at all times and does not drop below a value of 7032.661 during estimation. This is only slightly lower than the resampling threshold of 7500 (see Figure 9, upper part). The standard particle filter shows a much worse performance. As can be observed from the lower part of Figure 9, the ESS drops several times to very low values, reaching a minimum of 101.102. The ESS measures in some sense the variance of the normalized particle weights (low ESS means high variance) and highly varying weights indicate that the particle cloud is in a bad condition, since in this case a few particles with very high weights dominate the majority of the particles with very low weights. The extreme case is ESS=1 where only one particle has a significant weight. After a resampling step, only the information carried by this particle is present in the current particle cloud. Results obtained by estimation runs where such low ESS values have occurred cannot be trusted very much. If in contrast a resampling step is performed while the ESS is only slightly below the resampling threshold, which is the case in the MTU-PF algorithm, then many different particles will be chosen during resampling and a high percentage of the information contained in the particle cloud will be carried over to subsequent steps. In this sense, the MTU particle filter avoids the degeneration of the particle cloud by controlling the value of the ESS and holding it at a high value at every time. As a consequence, the estimation results should be considered to be more reliable than with the standard algorithm.

A look at the estimated values of the group parameters <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M307','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M307">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M308','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M308">View MathML</a> (see Figure 14) shows that in both estimation runs (MTU-PF and standard PF), the rate <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M309','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M309">View MathML</a> for diabetes patients is only about 60 percent of the rate <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M310','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M310">View MathML</a> of the control group (standard particle filter: 0.337h-1 vs 0.557h-1, MTU particle filter: 0.346h-1 vs 0.577h-1). The good performance especially of the MTU particle filter strengthens the confidence in the obtained result and leads to the conclusion that the secretion rate k0,1 is indeed lower for the diabetes patients than for the people from the control group.

Conclusions and future work

We proposed a new modification of the particle filter algorithm which works in continuous-time settings. It allows the direct inclusion of measurement time uncertainties in the underlying model. The modifications additionally allow the use of time-stepping strategies to improve the performance of the algorithm. The assumption of a random distribution of measurement times is natural in many applications.

The MTU-PF method is generally applicable. Even when measurement times may be assumed to be concentrated on single time points, our method can be used as a kind of regularization of the standard particle filter method if artificial distributions with highly concentrated masses around the measurement points are introduced.

We compared the performance of the MTU-PF to the standard PF and to an alternative Maximum Likelihood estimation method on a small artificial example. The results clearly show the advantage of the application of the MTU-PF in cases of uncertain measurement times.

We believe that our MTU particle filter is especially suitable for biological/medical applications where — compared to technical applications — the variance of the measurement values is relatively high due to biological variation and because relatively few consecutive measurements are possible. We provided an illustrative application from a PK/PD study. A comparison of our MTU particle filter with the standard filter showed that in this case our method is able to avoid weight degeneracy measured in terms of the Effective Sample Size (ESS) estimate. Even though the estimations in this application with both standard and MTU particle filters are reasonably good, the fit to the measurements is still not perfect. Whether that is due to our choices of model and parameters, or due to the known weaknesses of the general SMC method (which also our modification necessarily suffers from), remains to be evaluated in greater detail. In subsequent studies, we plan to apply the new algorithm to the complete liver-plasma model using additional measurements to be able to draw conclusions of greater medical relevance.

Another topic may also prove interesting for future work. Our experiments showed that the results are highly dependent on the choice of the development of the diffusion coefficients of the artificial noise necessary for Bayesian parameter estimation. While there is a general agreement that these coefficients should decrease over time, there is currently a lack of automated methods for making appropriate choices of the diffusion coefficients, both for initial values and dynamic development. However, to provide a really practical method for parameter estimation in non-linear mixed effects models (or even in models which adhibit only global parameters), our approach could also be combined with methods better suited to the estimation of fixed parameters, a good candidate being the PMCMC methods proposed in [11]. This is future work.

In summary, we believe that the method presented in this article opens the door to even more efficient and reliable state sampling and parameter estimation methods based on the particle filter algorithm operating on continuous-time stochastic state space systems.

Notation

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M311','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M311">View MathML</a> probability space

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M312','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M312">View MathML</a> arbitrary measurable space (for each t∈[t0,) with <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M313','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M313">View MathML</a>)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M314','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M314">View MathML</a><a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M315','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M315">View MathML</a> measurable random variable

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M316','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M316">View MathML</a> continuous-time Markov process with general state space <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M317','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M317">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M318','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M318">View MathML</a> state space of <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M319','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M319">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M320','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M320">View MathML</a> the pushforward measure of P under Xt, i.e. <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M321','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M321">View MathML</a><a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M322','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M322">View MathML</a> for all <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M323','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M323">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M324','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M324">View MathML</a> the pushforward measure of P under <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M325','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M325">View MathML</a> (with the corresponding product algebra)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M326','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M326">View MathML</a> the state space restricted to the interval [t0,t]

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M327','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M327">View MathML</a> the corresponding pushforward measure

N number of particles

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M328','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M328">View MathML</a> state samples at time t

Ks,t(xs, dxt) the Markov kernel of the process <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M329','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M329">View MathML</a> from time s to time t

a(x,t) drift vector

B(x,t) diffusion matrix

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M330','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M330">View MathML</a> multidimensional standard Wiener process

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M331','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M331">View MathML</a> measurable space

Y1:M observation random variable with values in measurable spaces <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M332','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M332">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M333','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M333">View MathML</a> conditional probability density with respect to a reference measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M334','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M334">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M335','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M335">View MathML</a> reference measure <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M336','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M336">View MathML</a> on<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M337','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M337">View MathML</a>

tj (j=1,…,M) observation times

Tj random variables modelling the uncertainty about exact observation times

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M338','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M338">View MathML</a> Lebesgue measure on the interval [t0,)

γj(tj) probability density of Tj with respect to <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M339','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M339">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M340','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M340">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M341','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M341">View MathML</a>, and Markov chain, pushforward

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M342','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M342">View MathML</a> measure, and kernel for importance sampling

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M343','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M343">View MathML</a> Radon-Nikodym derivative<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M344','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M344">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M345','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M345">View MathML</a> Radon-Nikodym derivative at start time t0

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M346','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M346">View MathML</a> unnormalized weight of particle i at time t

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M347','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M347">View MathML</a> normalized weight of particle i at time t

s-th resampling time

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M348','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M348">View MathML</a> (unnormalized) selection weight

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M349','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M349">View MathML</a> normalized selection weight (probability that particle i will be selected during resampling)

ι:II selection function on the index set I: ={1,…,N} for the -th resampling step

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M350','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M350">View MathML</a> data likelihood

ft full density at time t

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M351','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M351">View MathML</a> filter density at time t, i.e. only those observations yj are included for which tjt

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M352','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M352">View MathML</a> conditional probability density <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M353','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M353">View MathML</a> if Tjt, 1 otherwise

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M354','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M354">View MathML</a> stochastic process given by dWj,t(ω)=(gj(yj | Xt(ω),t) -1)γj(t) dt

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M355','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M355">View MathML</a> product process of the Wj,t, where j=1,…,M

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M356','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M356">View MathML</a> expectation of h(Xt) given Y1:M=y1:M with respect to the filtered state Xt

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M357','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M357">View MathML</a> partial weight at time t (w.r.t. the j-th observation)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M358','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M358">View MathML</a> weight at time t

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M359','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M359">View MathML</a> cumulative distribution function <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M360','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M360">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M361','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M361">View MathML</a><a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M362','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M362">View MathML</a>

t0=τ0<τ1…<τD time discretization

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M363','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M363">View MathML</a> cumulative product of the selection weights for particle i

Δτd=τd-τd-1 stepsize

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M364','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M364">View MathML</a> correction factor for the data likelihood at time τd

Qi mass of the tracee in compartment i

qi mass of the tracer in compartment i

Ui input for the tracee in compartment i (e.g. U1 denotes the influx into compartment 1)

ui input for the tracer in compartment i

kj,i transfer coefficient of the tracers from compartment i to compartment j

σi diffusion parameter

ξt log-normal noise (<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M365','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M365">View MathML</a>)

<a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M366','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M366">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M367','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M367">View MathML</a> degradation rate of plasma leucine for people in the diabetes group and the control group, respectively

ζp, ηp patient-dependent random factors modelling the parametric uncertainties among individuals (ζp= exp(ηp) with <a onClick="popup('http://www.biomedcentral.com/1752-0509/7/8/mathml/M368','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/7/8/mathml/M368">View MathML</a>)

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

JH and AK contributed to the underlying ideas of the method and its elaboration and implemented the algorithm. Both authors contributed equally to the initial manuscript. MJ and MA commented and refined the manuscript. MA provided the compartmental model and the data. MT acquired the data. MJ contributed to the discussion and supervised the project. All authors read and approved the final manuscript.

Acknowledgements

Thanks go to Michaela Höhne from the University of Kaiserslautern for her help in preparing the figures. Data was generated by Martin Adiels, Professor Marja-Riitta Taskinen, Helsinki University, and Professor Jan Borén, Göteborg University, and has been published elsewhere ([31,32]). We would like to thank the anonymous reviewers for their helpful suggestions and feedback. This work was partially supported by the Fraunhofer Internal Programs under Grant No. MAVO 819 620. The work was also supported by grants from the European Commission 7th Framework Programme (UNICELLSYS, grant No 201142), the Swedish Foundation for Strategic Research through the Gothenburg Mathematical Modelling Centre (GMMC), the Sahlgrenska Center for Cardiovascular and Metabolic Research (CMR), the Sigrid Juselius Foundation, and the Helsinki University Central Hospital Research Foundation, Helsinki, Finland.

References

  1. Overgaard RV, Jonsson N, Tornøe CW, Madsen H: Non-linear mixed-effects models with stochastic differential equations: implementation of an estimation algorithm.

    Journal of Pharmacokinetics and Pharmacodynamics 2005, 32:85-107. PubMed Abstract | Publisher Full Text OpenURL

  2. Moles CG, Mendes P, Banga JR: Parameter estimation in biochemical pathways: a comparison of global optimization methods.

    Genome research 2003, 13(11):2467-2474. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  3. Chou ICC, Voit EO: Recent developments in parameter estimation and structure identification of biochemical and genomic systems.

    Mathematical biosciences 2009, 219(2):57-83. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  4. Ashyraliyev M, Fomekong-Nanfack Y, Kaandorp JA, Blom JG: Systems biology: parameter estimation for biochemical models.

    FEBS Journal 2009, 276(4):886-902. PubMed Abstract | Publisher Full Text OpenURL

  5. Bohlin T: Practical Grey-Box Process Identification: Theory and Applications. Advances in Industrial Control, London: Springer; 2010. OpenURL

  6. Nielsen J, Madsen H, Young P: Parameter Estimation in Stochastic Differential Equations: An Overview.

    Annual Reviews in Control 2000, 24:83-94. OpenURL

  7. Crisan D, Doucet A: A survey of convergence results on particle filtering methods for practitioners.

    IEEE Transactions on Signal Processing 2002, 50(3):736-746. Publisher Full Text OpenURL

  8. Andrieu C, De Freitas N, Doucet A: Sequential MCMC for Bayesian model selection. In Proceedings of the IEEE Signal Processing Workshop on Higher-Order Statistics, 1999. Caesarea, Israel: IEEE; 1999:130-134. OpenURL

  9. Fearnhead P: MCMC, sufficient statistics and particle filters.

    Journal Of Computational And Graphical Statistics 2002, 11(4):848-862. Publisher Full Text OpenURL

  10. Storvik G: Particle filters for state-space models with the presence of unknown static parameters.

    IEEE Transactions on Signal Processing 2002, 50(2):281-289. Publisher Full Text OpenURL

  11. Andrieu C, Doucet A, Holenstein R: Particle Markov chain Monte Carlo methods.

    Journal of the Royal Statistical Society - Series B: Statistical Methodology 2010, 72(3):269-342. Publisher Full Text OpenURL

  12. Gordon N, Salmond D, Smith A: Novel approach to nonlinear/non-Gaussian Bayesian state estimation.

    IEE-Proceedings-F 1993, 140(2):107-113. OpenURL

  13. Cappé O, Godsill S, Moulines E: An overview of existing methods and recent advances in Sequential Monte Carlo.

    Proceedings of the IEEE 2007, 95(5):899-924. OpenURL

  14. Doucet A, de Freitas N, Gordon N (Eds): Sequential Monte Carlo Methods in Practice. Statistics for Engineering and Information Science, New York: Springer Verlag; 2001. OpenURL

  15. Beal S, Sheiner L: NONMEM User’s Guides. NONMEM Project Group, University of California, San Francisco: ; 1994. OpenURL

  16. Tornøe CW, Overgaard RV, Agersø H, Nielsen HA, Madsen H, Jonsson EN: Stochastic differential equations in NONMEM: implementation, application, and comparison with ordinary differential equations.

    Pharmaceutical Research 2005, 22(8):1247-1258. PubMed Abstract | Publisher Full Text OpenURL

  17. Berglund M, Sunnåker M, Adiels M, Jirstrand M, Wennberg B: Investigations of a compartmental model for leucine kinetics using non-linear mixed effects models with ordinary and stochastic differential equations.

    Mathematical Medicine and Biology: A Journal of the IMA 2011, 29(4):361-384. PubMed Abstract | Publisher Full Text OpenURL

  18. Donnet S, Samson A: EM algorithm coupled with particle filter for maximum likelihood parameter estimation of stochastic differential mixed-effects models. [http://hal.archives-ouvertes.fr/hal-00519576 webcite]. [Preprint version 2 – 21 jul 2011]

  19. Donnet S, Samson A: Parametric inference for mixed models defined by stochastic differential equations.

    ESAIM Probability and Statistics 2008, 12:196-218. OpenURL

  20. Racine-Poon A, Wakefield J: Statistical methods for population pharmacokinetic modelling.

    Statistical Methods in Medical Research 1998, 7:63-84. PubMed Abstract | Publisher Full Text OpenURL

  21. Sheiner L, Wakefield J: Population modelling in drug development.

    Statistical Methods in Medical Research 1999, 8(3):183. PubMed Abstract | Publisher Full Text OpenURL

  22. Andersen KE, Hojbjerre M: A population-based Bayesian approach to the minimal model of glucose and insulin homeostasis.

    Statistics in Medicine 2005, 24(15):2381-2400. PubMed Abstract | Publisher Full Text OpenURL

  23. Hu XL, Schön T, Ljung L: A basic convergence result for particle filtering.

    IEEE Transactions on Signal Processing 2008, 56(4):1337-1348. OpenURL

  24. Douc R, Cappé O, Moulines E: Comparison of Resampling Schemes for Particle Filtering. In Proceedings of the 4th International Symposium on Image and Signal Processing and Analysis, 2005. ISPA, Zagreb, Croatia: IEEE; 2005:64-69. OpenURL

  25. Hol JD, Schön TB, Gustafsson F: On resampling algorithms for particle filters. In Proceedings of Nonlinear Statistical Signal Processing Workshop (NSSPW). Cambridge, UK: IEEE; 2006:79-82. OpenURL

  26. Pitt MK, Shephard N: Filtering via simulation: auxiliary particle filter.

    Journal of the American Statistical Association 1999, 94:590-599. Publisher Full Text OpenURL

  27. Del Moral P, Doucet A, Jasra A: Sequential Monte Carlo samplers.

    Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2006, 68(3):411-436. Publisher Full Text OpenURL

  28. R Development Core Team: R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria: ; 2011.

    [http://www.R-project.org webcite]. [ISBN 3-900051-07-0]

    OpenURL

  29. Demant T, Packard CJ, Demmelmair H, Stewart P, Bedynek A, Bedford D, Seidel D, Shepherd J: Sensitive methods to study human apolipoprotein B metabolism using stable isotope-labeled amino acids.

    Am J Physiol 1996, 270(6 Pt 1):E1022-36. PubMed Abstract | Publisher Full Text OpenURL

  30. Cobelli C, Saccomani MP, Tessari P, Biolo G, Luzi L, Matthews DE: Compartmental model of leucine kinetics in humans.

    Am J Physiol 1991, 261(4 Pt 1):E539-50. PubMed Abstract | Publisher Full Text OpenURL

  31. Adiels M, Packard C, Caslake MJ, Stewart P, Soro A, Westerbacka J, Wennberg B, Olofsson SO, Taskinen MR, Borén J: A new combined multicompartmental model for apolipoprotein B-100 and triglyceride metabolism in VLDL subfractions.

    Journal of Lipid Research 2005, 46:58-67. PubMed Abstract | Publisher Full Text OpenURL

  32. Adiels M, Borén J, Caslake MJ, Stewart P, Soro A, Westerbacka J, Wennberg B, Olofsson SO, Packard C, Taskinen MR: Overproduction of VLDL1 Driven by Hyperglycemia Is a Dominant Feature of Diabetic Dyslipidemia.

    Arteriosclerosis, Thrombosis, and Vascular Biology 2005, 25(8):1697-1703. Publisher Full Text OpenURL

  33. Adiels M: A compartmental model for kinetics of apolipoprotein B-100 and triglycerides in VLDL1 and VLDL2 in normolipidemic subjects.

    Licentiate thesis,

    Chalmers University of Technology, Göteborg 2002

    OpenURL

  34. Adiels M, Olofsson SO, Taskinen MR, Borén J: Overproduction of very low-density lipoproteins is the hallmark of the dyslipidemia in the metabolic syndrome.

    Arteriosclerosis, thrombosis, and vascular biology 2008, 28(7):1225-1236. Publisher Full Text OpenURL

  35. Hürzeler M, Künsch HR: Approximating and Maximizing the Likelihood for a General State-Space Model. In Sequential Monte Carlo Methods in Practice. New York: Springer; 2001. OpenURL

  36. Doucet A, de Freitas N, Murphy KP, Russell SJ: Rao-Blackwellised Particle Filtering for Dynamic Bayesian Networks. In UAI ’00: Proceedings of the 16th Conference on Uncertainty in Artificial Intelligence. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc.; 2000:176-183. OpenURL

  37. Kloeden PE, Platen E: Numerical solution of stochastic differential equations. Berlin: Springer Verlag; 1999. OpenURL