Astrocytic tracer dynamics estimated from [1-11C]-acetate PET measurements

Mathematical Medicine and Biology Advance Access published November 24, 2014 Mathematical Medicine and Biology (2014) Page 1 of 16 doi:10.1093/imammb/dqu021

Astrocytic tracer dynamics estimated from [1-11 C]-acetate PET measurements Andrea Arnold Center for Quantitative Sciences in Biomedicine, North Carolina State University, Campus Box 8213, 2700 Stinson Drive, 308 Cox Hall, Raleigh, NC 27695-8213, USA and Department of Mathematics, North Carolina State University, Campus Box 8205, 2108 SAS Hall, 2311 Stinson Drive, Raleigh, NC 27695-8205, USA

Albert Gjedde Department of Neuroscience and Pharmacology, University of Copenhagen, Blegdamsvej 3B, DK-2200, Copenhagen, Denmark and Department of Nuclear Medicine and PET Centre, Aarhus University Hospitals, Aarhus, Nørrebrogade 44 DK-8000, Denmark Peter Iversen Department of Nuclear Medicine and PET Centre, Aarhus University Hospitals, Aarhus, Nørrebrogade 44 DK-8000, Denmark and Erkki Somersalo Department of Mathematics, Applied Mathematics and Statistics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, USA [Received on 13 March 2014; revised on 24 August 2014; accepted on 8 October 2014] We address the problem of estimating the unknown parameters of a model of tracer kinetics from sequences of positron emission tomography (PET) scan data using a statistical sequential algorithm for the inference of magnitudes of dynamic parameters. The method, based on Bayesian statistical inference, is a modification of a recently proposed particle filtering and sequential Monte Carlo algorithm, where instead of preassigning the accuracy in the propagation of each particle, we fix the time step and account for the numerical errors in the innovation term. We apply the algorithm to PET images of [1-11 C]-acetate-derived tracer accumulation, estimating the transport rates in a three-compartment model of astrocytic uptake and metabolism of the tracer for a cohort of 18 volunteers from 3 groups, corresponding to healthy control individuals, cirrhotic liver and hepatic encephalopathy patients. The distribution of the parameters for the individuals and for the groups presented within the Bayesian framework support the hypothesis that the parameters for the hepatic encephalopathy group follow a significantly different distribution than the other two groups. The biological implications of the findings are also discussed. Keywords: parameter estimation; tracer kinetics; PET imaging; particle filters; sequential Monte Carlo (SMC).

c The authors 2014. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. 

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

Daniela Calvetti∗ Department of Mathematics, Applied Mathematics and Statistics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, USA ∗ Corresponding author: Email: [email protected], [email protected]

2 of 16

A. ARNOLD ET AL.

1. Introduction

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

Understanding the neurochemistry and metabolic processes in brain relies on indirect measurements such as dynamic positron emission tomography (PET) records. Although PET imaging, i.e. the passage from photon counting data to activity map in the brain, is a well-established practice, the quantitative imaging of kinetic model parameters from imaging data is less well developed and understood. In vivo neurodynamics is usually described in terms of parametric compartment models (Gjedde et al., 2011), and the distribution of parameter values carries important information about potential disease of the brain; see, e.g. Gjedde et al. (2010) and Iversen et al. (2009). The estimation of kinetic parameters can be approached either from an optimization point of view, by fitting the model to the observations, or by using statistical methods. The advantage of the latter approach is that it naturally provides a framework for including uncertainties in the modelling paradigm, including model inputs such as the arterial tracer concentration, and measurement errors. The Bayesian statistical framework (Kaipio and Somersalo, 2004; Calvetti and Somersalo, 2009) allows qualitative complementary information to be encoded in the prior distribution, and it produces an output naturally equipped with a reliability estimate. In this article, we analyse the quantitative PET imaging problem using sequential Bayesian estimation and, in particular, particle filtering sequential Monte Carlo (PF-SMC) methods. In this work, we consider a kinetic model of [1-11 C]-acetate imaging. Acetate is a metabolite preferentially taken up into and metabolized by astrocytes, not neurons, which makes acetate tracers particularly suited for studying astrocytic metabolism in brain (Wanievski and Martin, 1998; Hertz et al., 2007; Wyss et al., 2011). The particular question that is addressed in this paper is related to the brain function during hepatic encephalopathy. It is known that serious liver malfunctioning affects, among other things, the cognitive functions, but the detailed mechanisms are still to be discovered. In particular, a question of interest is whether the elevated level of ammonium in plasma during hepatic encephalopathy affects primarily the neurons or astrocytes in brain. Acetate passes from the blood into the brain tissue through the blood–brain barrier and is further taken up by the astrocytes. Inside the cell, acetate enters the TCA cycle in the mitochondria, where the carbons are released to form carbon dioxide that is consequently cleared from the tissue. The labelling of the first carbon by the 11 C isotope, as opposed to labelling the second carbon of the acetate carbon skeleton, is motivated by the fact that in this way, the isotope is cleared already during the second turn of the TCA cycle, thus minimizing the potential contaminate of the neuron through the glutamineglutamate cycle (Wyss et al., 2011; Lanz et al., 2012). Experimental data collection and the inverse problem can be summarized as follows: a radioactive tracer is administered into the circulation of the subject, and the tracer concentration in arterial blood is estimated from direct arterial measurements of the radioactivity. The blood flow carries the tracer to the brain where it enters the tissue, and the time dependent positron emission activity, proportional to the tracer mass, is measured with a PET device and mapped to a series of tomographic images, yielding an estimate of the time dependent activity distribution in the tissue. More specifically, dynamic PET images track the concentration of radioactivity as a function of time in each image voxel, and voxels in the region of interest are averaged, yielding the mean tissue activity curves which are part of the data for the present study. The interpretation of the data requires a model describing the tracer’s kinetics. The tracer’s concentration in plasma, on the other hand, which constitutes the input of the model, is measured directly from blood samples. A standard way of describing the marker dynamics is to use a compartment model in which three compartments are identified: a blood compartment, a precursor compartment corresponding to acetate taken up by the tissue and a product compartment corresponding to acetate in the astrocyte as a metabolite. A schematic picture of the compartment model is shown in Fig. 1. The mathematical

ESTIMATION OF PET TRACER DYNAMICS

3 of 16

compartment model describing the mass balance in the brain tissue is governed by a system of differential equations (Gjedde et al., 2011) dm1 (t) = K1 c(t) − (k2 + k3 )m1 (t) dt dm2 (t) = k3 m1 (t) − k5 m2 (t) dt

(1.1) (1.2)

where c(t) is the arterial tracer concentration, or arterial input function (AIF), estimated from measured blood samples drawn during the measurement process, and m1 (t) and m2 (t) represent the mass of tracer content in the precursor and product compartments, respectively. The coefficients k2 , k3 and k5 are transport rates between the compartments, and K1 is the clearance of the tracer from blood to tissue, weighted with the capillary compartment’s virtual volume V0 . In this model, k4 , which would correspond to back-transport to the precursor compartment, is assumed to be vanishingly small because of the unidirectionality of the biochemical reaction represented by k3 . The initial values for the differential equations are m1 (0) = m2 (0) = 0, and the total tracer content in the tissue volume of interest is m(t) = V0 c(t) + m1 (t) + m2 (t),

(1.3)

where we assume that the PET image provides an estimate of its time course. This estimate, which we refer to as the tissue’s time–activity curve (TAC), is the data for the present problem. The above model is a version of the standard compartment models used to study tracer kinetics (Gjedde et al., 2011). The model does not address questions such as tracer recirculation and reuptake, thus the effect of the labelled CO2 on the PET data is incompletely addressed (Iida et al., 1986, 1993;

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

Fig. 1. A schematic representation of the compartment model of tracer kinetics. The dashed module represents the reuptake of CO2 , which is not included in the model.

4 of 16

A. ARNOLD ET AL.

Keiding and Pavese, 2013). Since the contribution of 11 CO2 to the input function was not corrected, the values for 11 C-acetate will be overestimated due to contributions of oxidation of acetate by other body tissues. This may be a confound for physiological interpretation of data in disease states involving the liver, but not necessarily for the mathematical analysis of a test data set. It is not the aim of this work to improve the existing model, but rather to develop a methodology to estimate the model parameters V0 , K1 and kj , where j = 2, 3, 5, from the PET image data using SMC methods. The approach relies on a Bayesian paradigm of inverse problems, in which the unknowns are modelled as random variables, and the data are used to infer the posterior probability distribution. We refer to Kaipio and Somersalo (2004) and Calvetti and Somersalo (2009) as general references for the Bayesian formalism, and to Liu and West (2001), Doucet et al. (2001, Chapter 10) and Arnold et al. (2013) as references for the SMC methodology.

We start by introducing some notations. The unknown parameters describing the tracer dynamics in tissue are collected in a single vector that we denote by θ , θ = (V0 , K1 , k2 , k3 , k5 ) ∈ R5 . The input function c(t) obtained from plasma radioactivity measurements is assumed to be known. Integrating the system (1.1)–(1.2) with the vanishing initial values gives the pair (m1 (t, θ ), m2 (t, θ )) and the idealized model for the errorless observation, m(t, θ ) = θ (1)c(t) + m1 (t, θ ) + m2 (t, θ ),

(2.1)

where θ (1) refers to the first component of the vector θ , while the observed data consist of noisy samples at discrete times tj , dj = m(tj , θ ) + ej , 1  j  T, (2.2) where, for simplicity, the noise is modelled as additive error. It is common to assume that ej represents only the measurement uncertainty, while in reality it should also account for the modelling error. This question will be addressed in the following section. In the Bayesian framework, randomness represents lack of information. Therefore, the unknown θ is modelled as a random variable, and the information about it is encoded in the respective probability distribution. We assume here that probability distributions can be expressed in terms of probability densities in the parameter space. The prior probability density, encoding all information about the unknown, is denoted by πpr (θ ). The probability density of the data dj , if θ were known, gives the likelihood density π(dj | θ ), and encodes the information about the noise ej . Finally, Bayes’ formula gives the posterior density, which is an update of the prior density based on the additional information that the data provide, π(θ | d) ∝ πpr (θ )π(d | θ ), where ‘∝’ stands for ‘equal up to a multiplicative normalizing constant’. 2.1

Error model and hierarchical update

Consider the TAC data model (2.2), where the vector ej represents the aggregation of measurement errors, imprecisions in the model, as well as model discrepancy. When referring to measurement error in

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

2. Bayesian framework

ESTIMATION OF PET TRACER DYNAMICS

5 of 16

the context of dynamic PET imaging, we mean the photon emission of the volume element in question, calculated from the PET counting data. By imprecisions in the model, we refer to uncertainties that have an effect on reliability of the estimates, most notably uncertainties related to the AIF; see, e.g. Iida et al. (1986, 1993). In this work, we address one possible error source, related to an incompletely known offset of the AIF and TAC curves: The blood sample data are collected independently of the PET scan, and the delay between the time when the tracer arrives at the brain and at the blood sample site varies. This delay may affect, in particular, the estimates based on data at the very beginning of the estimation process. Therefore, assume that the convection time from the brain to the site where the blood sample is drawn contains an unknown estimation error τ . To include this error in the model, we replace (2.2) by 1  j  T.

(2.3)

The offset τ could be added to the model as an unknown that needs to be estimated from the data. Here, we treat it as a nuisance parameter and, to simplify the discussion, we approximate m(t, θ ) by its first-order Taylor polynomial and write dj ≈ m(tj , θ ) + mt (tj , θ )τ + ej = m(tj , θ ) + nj + ej ,

(2.4)

where mt denotes the derivative of m with respect to t. As is common for the modelling errors, the offset error nj depends on the unknown θ itself. To overcome this difficulty, we propose here a hierarchical procedure. Assume that θˆ represents a typical value of the parameter θ . We model the prior density of θ in terms of the typical value, i.e. we define a conditional prior for the components θ (k) of the parameter vector, πpr (θ (k) | θˆ (k)) ∼ Uniform(α θˆ (k), Aθˆ (k)) for some scalars 0 < α < 1 < A. Furthermore, we approximate the modelling error term by nj = mt (tj , θ )τ ≈ mt (tj , θˆ )τ , such that the only source of randomness is the random offset τ , which we will model as a zero mean Gaussian random variable. To find a typical value θˆ , we consider the model (2.2), assuming that the offset vanishes. We assume for simplicity Gaussian measurement error with constant variance σ 2 , which leads to a likelihood model   1 π(dj | θ ) ∝ exp − 2 |dj − m(tj , θ )|2 . 2σ Further, by assuming conditionally independent measurements, the likelihood of the full data vector d is ⎞ ⎛ T T   1 π(dj | θ ) ∝ exp ⎝− 2 |dj − m(tj , θ )|2 ⎠ . π(d | θ ) = 2σ j=1 j=1

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

dj = m(tj + τ , θ ) + ej ,

6 of 16

A. ARNOLD ET AL.

Assuming a weakly informative prior imposing the positivity of the parameters θ (k) and some large upper bound, 0  θk  Θ(k), πpr (θ ) ∝

K 

1[0,Θ(k)] (θ (k)),

k=1

where 1[0,Θ(k)] is the indicator function of the interval [0, Θ(k)], the posterior density becomes ⎞ T  1 π(θ | d) ∝ 1[0,Θ(k)] (θ (k)) exp ⎝− 2 |dj − m(tj , θ )|2 ⎠ . 2σ j=1 k=1 ⎛

K 

(2.5)

j=1

The constrained minimization problem (2.6) is straightforward to implement. To estimate the sensitivity of estimated θ to modelling errors, a more comprehensive analysis is necessary. An extensive Markov chain Monte Carlo run for sampling the posterior is one option; a more approximate but computationally lighter alternative is to use parametric bootstrapping, which gives an idea of the posterior density assuming that the prior is non-informative. In this article we use a PF-SMC algorithm combined with the error remodelling approach explained as follows: 1. Estimate θ using the simple likelihood model, and compute the typical value θˆ from (2.6). 2. Write a conditional prior and conditional likelihood for the components θ (k), ˆ ∼ Uniform(α θˆ (k), Aθˆ (k)) πpr (θ (k) | θ(k)) for some scalars 0 < α < 1 < A, and the likelihood π(dj | θ , θˆ ) ∼ N (m(t, θ ), σ 2 + mt (tj , θˆ )2 var(τ )). 3. Run a PF-SMC algorithm to explore the conditional posterior density, π(θ | d, θˆ ) ∝ πpr (θ | θˆ )π(d | θ , θˆ ). The error variances σ 2 and var(τ ) are estimated heuristically from the data: a systematic and automatic way of assigning their values would require information about the measurement procedure that was not available. The state-of-the-art methodology in analysing dynamic PET data is to use non-linear curve fitting, see, e.g. Keiding and Pavese (2013). This amounts to solving the least squares optimization problem (2.6), which is naturally imbedded in the proposed empirical Bayes formulation. Hence, the proposed algorithm extends the existing methodology towards a more complete analysis of the variability of the estimates due to model uncertainties and imprecisions in data. In the next section, we briefly review the PF-SMC algorithm proposed in Arnold et al. (2013).

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

Maximizing the posterior distribution is tantamount to solving for the constrained least squares solution, ⎧ ⎫ T ⎨ ⎬ θˆ = argmin |dj − m(tj , θ )|2 subject to 0  θ (k)  Θ(k). (2.6) ⎩ ⎭

ESTIMATION OF PET TRACER DYNAMICS

2.2

7 of 16

SMC and survival of the fittest

In summarizing the PF-SMC algorithm, assume that the time interval corresponding to the data acquisition is discretized by time instances t0 < t1 < · · · < tT , some of which, but not necessarily all, coincide with the data acquisition. However, to simplify the notation, and to avoid introduction of new indices, we assume that data arrive at each tj , denoting dj = ∅ if the measurement is missing. We denote by Dj = {d1 , d2 , . . . , dj } the data acquired up to time t = tj . The posterior density of θ conditional on data Dj , π(θ | Dj , θˆ ), is summarized by the sample

Let θ¯j and Cj denote the sample-based mean and covariance of θ , θ¯j =

N 1  n θ , N n=1 j

1  n (θ − θ¯j )(θjn − θ¯j )T . N − 1 n=1 j N

Cj =

The probability density is approximated by a Gaussian mixture model, π(θ | Dj , θˆ ) ≈

N 

wnj N (θ | θ¯jn , s2 Cj ),

n=1

where the means θ¯jn , called auxiliary particles, are defined via a shrinkage towards the sample mean, θ¯jn = aθjn + (1 − a)θ¯j ,

0 < a < 1,

and s is a scaling factor to counterbalance the shrinkage, s2 = 1 − a2 . The weights wnj represent probabilities that are updated at each step recursively as described below. The initial sample, when j = 0, is ˆ and the weights are set equal, wn0 = 1/N. drawn from the conditional prior π(θ | θ), Let X (t) = (m1 (t), m2 (t)) denote the state variable at time t, satisfying the system of differential equations (1.1)–(1.2). We need to propagate the solution X from tj to tj+1 using a numerical scheme, denoted by Ψ . If Xj is the approximate value of the current state, Xj ≈ X (tj ), we have X (tj+1 ) ≈ Xj+1 = Ψ (Xj , θ ), where we indicate the dependency of the propagated value on the parameter vector θ . In the particle filter methods, we assume that at each t = tj we have a sample of state solutions Mj = {Xj1 , . . . , XjN }, which corresponds to the parameter sample Sj . The initial state sample is defined as X01 = · · · = X0N = 0, corresponding to the known initial value.

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

Sj = {θj1 , . . . , θjN }.

8 of 16

A. ARNOLD ET AL.

The PF-SMC algorithm thus proceeds through the following steps: 1. Survival of the fittest: (a) Propagate Mj using the auxiliary particles as parameters, n X¯ j+1 = Ψ (Xjn , θ¯jn ),

1  n  N.

n Denote by m ¯ nj the predicted observation (2.1) computed with X¯ j+1 and θ¯jn .

Scale the fitness parameters gjn by their sum over the particle index n. (c) Draw indices 1 , . . . , N from the index set {1, 2, . . . , N} with replacement, where gjn is the probability for n to be drawn. n ¯ n n (d) Reshuffle the samples, replacing (X¯ j+1 , θ¯jn ) with (X¯ j+1 , θj ), and Xjn with Xj n . In this new n sample, the most fit auxiliary particles, for which gj is large, appear more often than the less fit particles.

2. Proliferation: n (a) Draw the parameter realizations θj+1 , 1  n  N, from the Gaussian N (θ | θ¯jn , s2 Cj ).

(b) Repropagate each resampled particle with the new parameters and add innovation, n n Xj+1 = Ψ (Xjn , θj+1 ) + vnj+1 ,

1  n  N,

where vnj+1 is drawn from a zero mean Gaussian distribution with variance set to correspond to the estimated numerical approximation error of the selected numerical scheme Ψ . (c) Update the weights to equal the likelihood ratio, wnj+1 =

n n π(dj+1 | Xj+1 , θj+1 , θˆ ) . n n π(dj+1 | X¯ j+1 , θ¯j+1 , θˆ )

(d) Recompute the parameter ensemble mean and covariance θ¯j+1 and Cj+1 . The PF-SMC algorithm recently proposed in Arnold et al. (2013) uses linear multistep method (LMM) numerical schemes with a fixed time step to solve the systems of differential equations. We refer to the original article for the details on how higher order LMM solvers can be expressed in the proper Markovian framework and the assignment of the innovation variance in Step 2(b) using error control estimates.

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

(b) Compute the fitness of the auxiliary particle θ¯jn by evaluating the weighted likelihood functions,   1 n n n 2 ¯ j | , σj2 = σ 2 + mt (tj θˆ )2 var(τ ). gj = wj exp − 2 |dj − m σj

ESTIMATION OF PET TRACER DYNAMICS

9 of 16

3. Results and discussion The algorithm described in the previous section was applied to data collected from a cohort of 18 volunteers from 3 different groups: a healthy control group of 5 individuals, a group of 7 individuals suffering from cirrhotic liver, and a group of 6 individuals with the condition of hepatic encephalopathy. As pointed out in the introduction, the motivation for the study was to find whether the liver condition of the subjects in the last two groups, through a perturbed ammonium level in their circulation, significantly affected the astrocytic metabolism at the level of acetate uptake and clearance. The details about the cohort of subjects, the scanning procedures and the reconstruction of the data can be found in Iversen et al. (2009). The [1-11 C]-acetate PET protocol lasted 45 min and the data corresponds to 50 time frames, with the interval between them increasing in time: every 5 s for the first 18 frames, every 10 s for the next 12 frame, then every 30 s for the subsequent 7 frames, followed by one frame after 120 s and finishing with 12 frame every 180 s. During the [1-11 C]-acetate recording, arterial blood samples of 0.5 ml each were collected manually at the following intervals: every 5 s for the first 18 times, every 10 s for the next 9 times, then every 60 s for the subsequent 12 frames, and finally 300 s for the last 6 times. Total blood 11 C radioactivity concentrations were measured using a well counter (Packard Instrument Company, Meriden, CT, USA), cross-calibrated with the tomograph, and corrected for radioactive decay to the start of the scan. The volume of interest in our computed examples is that of the whole-brain, and the corresponding TACs were extracted from the dynamic PET scan. The study was performed in accordance with the Helsinki II Declaration and approved by the ethics committee of Aarhus County. Informed consent was obtained from each subject. No complications were observed during the procedures. The 18 subjects are labelled with an index p, 1  p  18, and we identify the index sets with the corresponding group, writing HC = {1  p  5},

CL = {6  p  12},

HE = {13  p  18},

where HC refers to healthy control, CL to cirrhotic liver and HE to hepatic encephalopathy. For each individual p, we record the optimized least squares estimate and the parameter sample with corresponding weights as n θˆ(p) , (θ(p) , wn(p) ), 1  n  N.

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

The basic PF-SMC algorithm was modified for the current application to meet the challenges posed by the initial, very rapid ramping of the tracer’s mass and the 100-fold increase in the intervals between data measurements from the beginning to the end of the procedure. More specifically, instead of keeping the time step fixed for each particle for the entire duration of the experiment, we adjusted the time step according to the frequency of the data collection, keeping it the same for all particles. The details about how the time step was selected and how the LMM had to be modified can be found in Arnold (2014). During the up-ramping phase of the experiment, the particle retention rate (PRR), defined as the percentage of fit particles, may become prohibitively low: to avoid such an occurrence, if the PRR drops below a certain threshold, which was set to 60% in the current application, we resample until the PRR is above the threshold. Because of the complexity of the problem and the need for a large number of particles, the PFSMC algorithm is computationally demanding. It was recently shown in Arnold et al. (2014) that it is possible to reformulate the problem in so as to take advantage of vectorized and parallelized computing environment, significantly reducing the computing time. Our computed examples were implemented using a vectorized version of the PF-SMC algorithm.

10 of 16

A. ARNOLD ET AL.

V0

K1

0.200

0.015

0.150 0.010 0.100 0.005 0.050 0.000

HC

CL

0.000

HE

HC

0.009

0.020

0.006

0.010

0.003

HC

CL

HE

k3

0.030

0.000

CL

0.000

HE

HC

CL

HE

k5 0.007 0.005 0.003 0.001 0.000

HC

CL

HE

Fig. 2. Box plots indicating the 50 and 90% credibility intervals for each subject with respect to the five parameters. Box plots for the healthy control group (labelled as HC) are shown in black, for the cirrhotic liver group (CL) in blue and for the hepatic encephalopathy group (HE) in red, respectively. The red dot on each box plot marks the optimized least squares estimate of the corresponding parameter and subject.

To obtain the parameter samples, for each individual p we run the hierarchical PF-SMC algorithm, filtering a sample size of N = 100, 000 particles, using the third-order Adams–Moulton (AM) scheme for the time propagation and the fourth-order AM scheme for the error estimation. The resulting parameter distributions are summarized in Fig. 2, where the box plots indicate the 50 and 90% credibility intervals of the conditional distributions π(θ(p) | θˆ(p) , d(p) ), where d(p) is the data vector

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

k2

11 of 16

ESTIMATION OF PET TRACER DYNAMICS

Table 1 LMM PF-SMC parameter results for each subject p in the healthy control group, labelled by subject number 1  p  5. The optimized parameter value (computed via least squares) as well as the posterior mean of the filtered sample distribution and the upper (75%) and lower (25%) quartiles are given V0

K1

k2

k3

k5

8.8914e−03 8.9489e−03 1.1578e−02 6.3371e−03

9.1126e−03 8.8374e−03 9.7043e−03 7.9528e−03

2.4863e−02 2.3401e−02 2.5913e−02 2.0803e−02

5.7240e−03 5.3588e−03 6.7686e−03 3.8626e−03

5.0372e−03 4.7153e−03 5.4572e−03 3.8935e−03

Subject 2 Optimized value Posterior mean Upper quartile Lower quartile

2.7843e−02 3.2851e−02 3.8617e−02 2.8162e−02

8.8373e−03 8.4330e−03 9.2887e−03 7.5562e−03

2.3209e−02 2.2079e−02 2.4475e−02 1.9564e−02

3.4826e−03 3.2498e−03 4.0786e−03 2.3685e−03

3.0598e−03 2.7650e−03 3.2039e−03 2.2597e−03

Subject 3 Optimized value Posterior mean Upper quartile Lower quartile

9.6508e−02 1.1467e−01 1.2852e−01 1.0360e−01

9.7631e−03 9.7486e−03 1.0940e−02 8.5796e−03

1.8257e−02 1.8837e−02 2.1130e−02 1.6624e−02

2.2531e−03 2.4727e−03 3.0072e−03 2.0034e−03

2.9420e−03 3.0761e−03 3.5728e−03 2.5951e−03

Subject 4 Optimized value Posterior mean Upper quartile Lower quartile

8.3579e−02 9.7014e−02 1.0644e−01 8.8837e−02

5.6257e−03 5.3257e−03 5.9777e−03 4.6361e−03

1.3614e−02 1.3399e−02 1.5076e−02 1.1686e−02

1.7553e−03 1.6429e−03 2.0786e−03 1.1827e−03

1.6701e−03 1.4641e−03 1.6888e−03 1.2124e−03

Subject 5 Optimized value Posterior mean Upper quartile Lower quartile

7.4945e−02 8.8266e−02 9.5864e−02 8.1943e−02

3.4092e−03 2.7094e−03 3.0568e−03 2.2705e−03

1.2487e−02 9.8134e−03 1.1424e−02 7.6249e−03

4.3138e−03 3.6135e−03 4.5033e−03 2.5508e−03

1.3619e−03 1.0863e−03 1.2022e−03 9.4770e−04

of the pth subject. The results are also given in tabular form, organized by subject group: healthy control subjects in Table 1, cirrhotic liver subjects in Tables 2 and 3, and hepatic encephalopathy subjects in Table 4. These results indicate that while there is no significant difference in the estimated parameters between the healthy control and cirrhotic liver groups, in hepatic encephalopathy subjects the parameters related to acetate uptake, in particular, i.e. V0 and K1 , show a clear offset. Likewise, the k5 values in the hepatic encephalopathy group seem to be consistently lower than in the other groups. To quantify the differences between the groups beyond a mere visual inspection, we performed the following statistical Bayes ratio test. Let X denote a chosen reference group (healthy control, cirrhotic liver or hepatic encephalopathy) and Y a selected test group. We start by fitting a Gaussian distribun | p ∈ X , 1  n  N}. Given the parameter θ(p) for an tion N (θ | μX , CX ) to the combined sample {θ(p)

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

Parameter Subject 1 Optimized value Posterior mean Upper quartile Lower quartile

12 of 16

A. ARNOLD ET AL.

Table 2 LMM PF-SMC parameter results for subjects 6  p  10 in the cirrhotic liver group, labelled by subject number. The optimized parameter value (computed via least squares) as well as the posterior mean of the filtered sample distribution and the upper (75%) and lower (25%) quartiles are given V0

K1

k2

k3

k5

6.2424e−02 6.4943e−02 7.9495e−02 5.0891e−02

8.9456e−03 8.7107e−03 9.7667e−03 7.6447e−03

2.2924e−02 2.2542e−02 2.5332e−02 1.9674e−02

2.5699e−03 2.5062e−03 3.1902e−03 1.8031e−03

3.2486e−03 2.9698e−03 3.4833e−03 2.3962e−03

Subject 7 Optimized value Posterior mean Upper quartile Lower quartile

8.6474e−02 8.1059e−02 9.1732e−02 6.9716e−02

6.3450e−03 6.3240e−03 7.0845e−03 5.5673e−03

1.3867e−02 1.3928e−02 1.5591e−02 1.2283e−02

2.0623e−03 2.0990e−03 2.5934e−03 1.6111e−03

2.0401e−03 1.9654e−03 2.2537e−03 1.6681e−03

Subject 8 Optimized value Posterior mean Upper quartile Lower quartile

1.2194e−01 1.2634e−01 1.4114e−01 1.1194e−01

6.8327e−03 6.9997e−03 7.9292e−03 6.1090e−03

1.5541e−02 1.6248e−02 1.8245e−02 1.4336e−02

1.3862e−03 1.3980e−03 1.7431e−03 1.0566e−03

2.0483e−03 2.0432e−03 2.3962e−03 1.6836e−03

Subject 9 Optimized value Posterior mean Upper quartile Lower quartile

7.9554e−02 8.6554e−02 9.6469e−02 7.7539e−02

4.6056e−03 4.4903e−03 5.0371e−03 3.9271e−03

1.2285e−02 1.2433e−02 1.4032e−02 1.0840e−02

2.4207e−03 2.4304e−03 3.0483e−03 1.8147e−03

1.3916e−03 1.3011e−03 1.4936e−03 1.0996e−03

Subject 10 Optimized value Posterior mean Upper quartile Lower quartile

4.5701e−02 5.5152e−02 6.3375e−02 4.8836e−02

6.6120e−03 6.1530e−03 6.9295e−03 5.3202e−03

2.1954e−02 2.0826e−02 2.3644e−02 1.7818e−02

3.5441e−03 3.4431e−03 4.3590e−03 2.4984e−03

2.5974e−03 2.3465e−03 2.7102e−03 1.9483e−03

individual from the test group, p ∈ Y , we can evaluate the mean likelihood, averaged over the conditional posterior, L(p | X , θˆ(p) ) = ≈



N (θ(p) | μX , CX )π(θ(p) | θˆ(p) , d(p) )dθ(p)

N 

n wn(p) N (θ(p) | μX , CX ).

n=1

This quantity expresses the average likelihood of getting the observation θ(p) if it would come from the distribution calculated from the reference group data. To average over the test group, we fit the sample {θˆ(p) | p ∈ Y } with a uniform distribution, which we assume represents the distribution of the hyperparameter θˆ over the test group Y , denoted π(θˆ | Y ).

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

Parameter Subject 6 Optimized value Posterior mean Upper quartile Lower quartile

13 of 16

ESTIMATION OF PET TRACER DYNAMICS

Table 3 LMM PF-SMC parameter results for subjects p = 11 and p = 12 in the cirrhotic liver group, labelled by subject number. The optimized parameter value (computed via least squares) as well as the posterior mean of the filtered sample distribution and the upper (75%) and lower (25%) quartiles are given V0

K1

k2

k3

k5

1.6253e−01 1.5490e−01 1.6371e−01 1.4573e−01

3.7216e−03 4.0900e−03 4.6239e−03 3.6094e−03

1.1151e−02 1.2413e−02 1.4211e−02 1.0790e−02

1.5924e−03 1.7288e−03 2.1220e−03 1.3669e−03

1.3767e−03 1.4476e−03 1.6551e−03 1.2471e−03

Subject 12 Optimized value Posterior mean Upper quartile Lower quartile

6.8065e−02 8.2760e−02 9.2506e−02 7.5241e−02

6.7645e−03 6.4842e−03 7.2070e−03 5.7296e−03

1.6687e−02 1.6449e−02 1.8144e−02 1.4722e−02

1.9993e−03 1.9803e−03 2.5240e−03 1.4388e−03

2.2676e−03 2.2527e−03 2.6489e−03 1.8481e−03

The uniform distribution is justified by the small sample size and the rejection of a multivariate normal distribution via Mardia’s test for skewness and kurtosis; see Mardia (1970, 1974) for test details. The likelihood averaged over the combined posterior density  L(Y , X ) = L(p | X , θˆp )π(θˆp | Y )d θˆp ≈

nY 1  L(p | X , θˆ(p) ), nY p=1

can be thought of as an average measure for the probability of the sample of group Y as a whole, assuming that it came from the distribution of the reference group X. To obtain meaningful relative measure of how well the groups fit, we compute a version of a Bayes ratio, defining R(Y , X ) =

L(Y , X ) , L(X , X )

where the reference value R(Y , X ) = 1 would indicate a perfect match. Table 5 lists the Bayes ratios R for all possible pairs X and Y, while Table 6 gives the Bayes ratios comparing healthy control with the combined diseased group, cirrhotic liver + hepatic encephalopathy. While these ratios do not suggest any stark differences between healthy control and cirrhotic liver or the combined diseased group, Table 5 shows a clear distinction between hepatic encephalopathy and the other groups, both when hepatic encephalopathy cohort is treated as the control group and as the test group. From the physiological point of view, the lowered values of V0 and K1 in the hepatic encephalopathy group point towards a lowered cerebral blood flow (CBF), which is in line with previous findings (Keiding and Pavese, 2013). On the other hand, it has been previously found (Gjedde et al., 2010) that oxygen delivery is not limited in hepatic encephalopathy, while a decrease in the oxygen cerebral metabolic rate (CMRO2 ) was observed in this cohort. Therefore, we may conclude that the low CMRO2 must be due to

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

Parameter Subject 11 Optimized value Posterior mean Upper quartile Lower quartile

14 of 16

A. ARNOLD ET AL.

Table 4 LMM PF-SMC parameter results for each subject p in the hepatic encephalopathy group, labeled by subject number 13  p  18. The optimized parameter value (computed via least squares) as well as the posterior mean of the filtered sample distribution and the upper (75%) and lower (25%) quartiles are given V0

K1

k2

k3

k5

4.4654e−02 5.3271e−02 6.0117e−02 4.7824e−02

4.4032e−03 4.2031e−03 4.6359e−03 3.7525e−03

1.1005e−02 1.0523e−02 1.1689e−02 9.3077e−03

1.3250e−03 1.2412e−03 1.5647e−03 9.0020e−04

1.2227e−03 1.1120e−03 1.2983e−03 9.0731e−04

Subject 14 Optimized value Posterior mean Upper quartile Lower quartile

1.6110e−02 2.0956e−02 2.3107e−02 1.9617e−02

4.2687e−03 4.1881e−03 4.4592e−03 3.9188e−03

1.0831e−02 1.0793e−02 1.1625e−02 9.9625e−03

1.6558e−03 1.5888e−03 2.0094e−03 1.1517e−03

1.9719e−03 1.8064e−03 2.0693e−03 1.5186e−03

Subject 15 Optimized value Posterior mean Upper quartile Lower quartile

3.1157e−02 4.0486e−02 4.3950e−02 3.8186e−02

4.6739e−03 4.5001e−03 4.7844e−03 4.2048e−03

1.1737e−02 1.1340e−02 1.2243e−02 1.0404e−02

2.1741e−03 1.9457e−03 2.4494e−03 1.3962e−03

1.8991e−03 1.6556e−03 1.8711e−03 1.4190e−03

Subject 16 Optimized value Posterior mean Upper quartile Lower quartile

2.7991e−02 3.5249e−02 3.9306e−02 3.2354e−02

2.7334e−03 2.3085e−03 2.5910e−03 1.9808e−03

1.6681e−02 1.3433e−02 1.5423e−02 1.0935e−02

5.5866e−03 4.8445e−03 5.9660e−03 3.5671e−03

9.8447e−04 8.2596e−04 9.1559e−04 7.2420e−04

Subject 17 Optimized value Posterior mean Upper quartile Lower quartile

3.9232e−02 4.4772e−02 4.8628e−02 4.1349e−02

3.3517e−03 2.9977e−03 3.2354e−03 2.7418e−03

1.0546e−02 9.1167e−03 9.9861e−03 8.1459e−03

3.5964e−03 3.0473e−03 3.6583e−03 2.3450e−03

1.6859e−03 1.4430e−03 1.5816e−03 1.2868e−03

Subject 18 Optimized value Posterior mean Upper quartile Lower quartile

4.6142e−02 5.8406e−02 6.3940e−02 5.4438e−02

4.4364e−03 4.2594e−03 4.6723e−03 3.8351e−03

1.0484e−02 1.0239e−02 1.1291e−02 9.1618e−03

1.5789e−03 1.5104e−03 1.9177e−03 1.0871e−03

1.8681e−03 1.7735e−03 2.0431e−03 1.4834e−03

a not specified metabolic inhibition, and the low CBF is a consequence, rather than a cause, of the low CMRO2 . The natural question that arises is whether the metabolic inhibition occurs in astrocytes, neurons or both. The relatively insignificant difference in the parameter k3 in the three cohorts suggests that the oxidative metabolism of the astrocyte in hepatic encephalopathy is not compromised, thus shifting the attention to the oxidative metabolism in neuron. Methodologically, the current approach for the parameter estimation accounts for the intrinsic variability of the parameter values within individuals, not only within a cohort. It also provides a natural way to assess the sensitivity of the data to the parameters, removing the need for a separate sensitivity

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

Parameter Subject 13 Optimized value Posterior mean Upper quartile Lower quartile

ESTIMATION OF PET TRACER DYNAMICS

15 of 16

Table 5 Bayes ratios for all possible pairs of control groups X and test groups Y X Y HC CL HE

HC 1.0000 0.9591 0.2205

CL 0.5821 1.0000 0.0444

HE 0.0042 0.0034 1.0000

X Y HC CL + HE

HC 1.0000 0.6182

CL + HE 0.5285 1.0000

analysis based on linearization or bootstrapping. This approach is in line with the current interest in uncertainty quantification within the mathematical sciences. Moreover, the estimation of the distribution of the parameters in sample form makes it possible to carry out correlational studies between pairs of parameters to uncover hidden dependencies. Acknowledgement This work was partly supported by NSF DMS project number 1312424 (Erkki Somersalo) and by grant number 246665 from the Simons Foundation (Daniela Calvetti). In addition to co-authors Albert Gjedde and Peter Iversen, we acknowledge the following members of the ‘liver-brain group’ at Aarhus University and the University of Copenhagen who made the study possible: Kim Mouridsen PhD, Mikkel B. Hansen PhD, Svend B. Jensen PhD, Michael Sørensen MD PhD, Lasse K. Bak PhD, Helle Waagepetersen PhD, Arne Schousboe DSc, Peter Ott MD DSc, Hendrik Vilstrup MD DSc, and Susanne Keiding MD DSc. References Arnold, A. (2014) Sequential Monte Carlo parameter estimation for differential equations. Ph.D. Thesis, Case Western Reserve University, Cleveland, OH. Arnold, A., Calvetti, D. & Somersalo, E. (2013) Linear multistep methods, particle filtering and sequential Monte Carlo. Inverse Probl., 29, 085007. Arnold, A., Calvetti, D. & Somersalo, E. (2014) Vectorized and parallel SMC parameter estimation for stiff ODEs (submitted for publication). Preprint: arXiv:1411.1702 [stat.CO]. Calvetti, D. & Somersalo, E. (2009) Introduction to Bayesian Scientific Computing. New York: Springer. Doucet, A., de Freitas, N. & Gordon, N. (eds) (2001) Sequential Monte Carlo Methods in Practice. New York: Springer.

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

Table 6 Bayes ratios comparing the healthy control group with the combined diseased group, cirrhotic liver + hepatic encephalopathy

16 of 16

A. ARNOLD ET AL.

Downloaded from http://imammb.oxfordjournals.org/ by guest on January 31, 2015

Gjedde, A., Bauer, W.R. & Wong, D.F. (2011) Neurokinetics: The Dynamics of Neurobiology In Vivo. New York: Springer. Gjedde, A., Keiding, S., Vilstrup, H. & Iversen, P. (2010) No oxygen delivery limitation in hepatic encephalopathy. Metab. Brain Dis., 25, 57–63. Hertz, L., Peng, L. & Dienel, G.A. (2007) Energy metabolism in astrocytes: high rate of oxidative metabolism and spatiotemporal dependence on glycolysis/glycogenolysis. J. Cereb. Blood Flow Metab., 27, 219–249. Iida, H., Jones, T. & Miura, S. (1993) Modeling approach to eliminate the need of separate arterial plasma in oxygen-15 inhalation positron emission tomography. J. Nucl. Med., 34, 1333–1340. Iida, H., Kanno, I., Miura, S., Murakami, M., Takahashi, K. & Uemura, K. (1986) Error analysis of quantitative cerebral blood flow measurement using H2 1 5O autoradiography and positron emission tomography, with respect to the dispersion of the input function. J. Cereb. Blood Flow Metab., 6, 536–545. Iversen, P., Sørensen, M., Bak, L.K., Waagepetersen, H.S., Vafaee, M.S., Borghammer, P., Mouridsen, K., Jensen, S.B., Vilstrup, H., Schousboe, A., Ott, P., Gjedde, A. & Keiding, S. (2009) Low cerebral oxygen consumption and blood flow in patients with cirrhosis and an acute episode of hepatic encephalopathy. Gastroenterology, 136, 863–871. Kaipio, J.P. & Somersalo, E. (2004) Statistical and Computational Inverse Problems. New York: Springer. Keiding, S. & Pavese, N. (2013) Brain metabolism in patients with hepatic encephalopathy studied by PET and MR. Arch. Biochem. Biophys., 536, 131–142. Lanz, B., Uffmann, K., Wyss, M.T., Weber, B., Buck, A. & Gruetter, R. (2012) A two-compartment mathematical model of neuroglial metabolism using [1-11 C] acetate. J. Cereb. Blood Flow Metab., 32, 548–559. Liu, J. & West, M. (2001) Combined parameter and state estimation in simulation-based filtering. Sequential Monte Carlo Methods in Practice (A. Doucet, J.F.G. de Freitas & N.J. Gordon eds). New York: Springer, pp. 197–223. Mardia, K.V. (1970) Measures of multivariate skewness and kurtosis with applications. Biometrika, 57, 519–530. Mardia, K.V. (1974) Applications of some measures of multivariate skewness and kurtosis in testing normality and robustness studies. Sankhy¯a: Indian J. Stat., Ser. B, 36, 115–128. Wanievski, R.A. & Martin, D.L. (1998) Preferential utilization of acetate by astrocytes is attributable to transport. J. Neurosci., 18, 5225–5233. Wyss, M.T., Magistretti, P.J., Buck, A. & Weber, B. (2011) Mini review: labeled acetate as a marker of astrocytic metabolism. J. Cereb. Blood Flow Metab., 31, 1668–1674.