Technical Report:
Time-Activity-Curve Integration in Lu-177 Therapies in Nuclear Medicine

from

Dr. rer. nat. Theresa Ida Götz

from Weiden (Germany)

SPECT/CT

     Single-Photon Emission Computed Tomography / CT
PET/CT

     Positron-Emission Tomography / CT
SPECT

     Single-Photon Emission Computed Tomography
PET

     Positron-Emission Tomography
CT

     Computed Tomograpy
MIRD

     Medical Internal Radiation Dose Committee
DVK

     Dose Voxel Kernel
DPK

     Dose Point Kernel
MC

     Monte Carlo Simulation
NET

     Neuroendocrine tumor
ADT

     Androgen Deprivation Therapy
PRRT

     Peptide Radionuclide Receptor Therapy
mCSPC

     metastativc Castration-Sensitive Prostate Cancer
SSTR

     Somatostatin Receptor
TTD

     Total Tumor Dose
PC

     Prostate Cancer
PSMA

     Prostate Specific Membrane Antigen
TAC

     Time Activity Curve
EM

     Expectation Maximization
TIA

     Time Integrated Activity
VOI

     Volume Of Interest
OLINDA

     Organ Level Internal Dose Assessment with Exponential Modeling
SE

     Standard Error
SD

     Standard Deviation
SS

     Sample Size
PF

     Particle Filter
SNR

     Signal to Noise Ratio
DNN

     Deep Neural Networks
EMD

     Empirical Mode Decomposition
EEMD

     Ensemble Empirical Mode Decomposition
MEEMD

    Multi-dimensional Ensemble Empirical Mode Decomposition
BEMD

     Bi-Dimensional EMD
CNN

     Convolutional Neural Network
IMF

     Intrinsic Mode Function
BIMF

     Bi-dimensional Intrinsic Mode Function
BEEMD

     Bi-dimensional Ensemble Empirical Mode Decomposition (EEMD)
GiT-BEMD

     ????
BIM

     Bi-dimensional Intrinsic Mode
MLE

     Maximum Likelihood Estimation
ReLU

     Leaky Rectified Linear Units
MLP

     Multi-Layer Perceptrons
LCN

     Local Contrast Normalization
SGD

     Stochastic Gradient Descent
KL

     Kullback-Leibler divergence
CgA

     Chromogranin A level
PSA

     Prostate Specific Antigen
AUC

     Area Under the Curve
SEM

     State Evolution Model
OM

     Observation Model
PFF

     Particle Fit
SF

     Simple Fit
HU

     Hounsfield Units
STS

     Soft Tissue Scaling
CS

     Center Scaling
DS

     Density Scaling
PS

     Percentage Scaling
DVK-NN

     DVK via Neural Networks
NN

     Neural Network
DE-NN-EMD

     Dose Estimation via Neural Network and EMD
DE-NN

     Dose Estimation via Neural Network
LOOCV

     Leave-One-Out Cross-Validation
RLT

     Radioligand Therapy
FBP

    Filtered Back-Projection
SIR

    Sampling Importance Resampling
PCA

    Principal Component Analysis
CMSC

    Comparable Minimal Scale Combination Principle
GT

    Ground truth

##  Chapter 1 Introduction

Therapy with Lu-177-labelled somatostatin analogues and PSMA inhibitors are
established tools in the management of patients with inoperable and
metastasised neuroendocrine tumours (NET) and prostate cancer (PCa) [27, 8].
Dosimetry could provide information for assessing the patient-individual dose-
burden due to radiotherapy with Lu-177-DOTATOC and Lu-177-PSMA as
radiopharmaceutical [32]. Most commonly, the pharmaceutical is delivered over
several (3-6) therapy cycles, where at the first cycle a dosimetry is
performed. For this, several SPECT images can be collected at multiple time
points. The resulting organ-wise time-activity-curve (TAC) is then integrated
to estimate the cumulated activity, which, together with the so called S-value
determines the energy dose, according to the Medical Internal Radiation Dose
(MIRD) scheme. From this, the dose limiting organ is identified (often being
kidneys), and the total amount of injected activity which can be applied over
multiple subsequent radiotherapy cycles is determined.
Commonly, the dose assessment is limited to mean-organ doses, since per-voxel
3-D dose-distributions are often regarded as too noisy for estimating absorbed
energy doses. Nevertheless, some studies from literature regarding this topic
are available.

###  1.1 Literature on Lu-177 Dosimetry

In the following, a brief overview of the literature regarding dosimetry based
on organ- or voxel-wise time-activity curves is provided:

Patient-specific dosimetry for I-131 tyroid tumor therapy is reported in [25].
Such personalized dosimetry better accounts for radionuclide distribution and
patient anatomy. However, its accuracy depends on the quality of the cumulated
activity which enters the dosimetric calculations and which is derived from
SPECT or PET imaging studies. Three to four image acquisitions were performed
in $15$ patients with metastatic tyroid carcinomas. A voxel-wise integration
over time of the PET images provided a 3D spatial distribution of accumulated
activity, and absorbed dose maps were obtained employing the software package
3D Internal Dosimetry. The study indicated a high inter-tumor and inter-
patient variability of absorbed dose distribution underlining the need for a
patient-specific dosimetry.

A similar study is reported in [30], where the dosimetry in I-131 internal
emitter therapy is discussed for follicular lymphoma patients being treated
with I-131-tositumomab. The study employs voxel-wise integration of time-
activity-curves derived from SPECT/CT images at multiple time points, which
were subsequently registered at multiple time points on a dual-modality
scanner. Maps of integrated TACs have been used together with CT images to
determine absorbed energy dose spatial distributions by Monte Carlo (MC)
simulations. Results are compared to state-of-the-art dosimetric calculations.

A PRRT with Lu-177-labeled somatostatin analogues has been reported already by
[29]. In their study, the authors compared the dosimetric parameter uptake,
half-life (kinetics), and mean absorbed organ and tumor doses of Lu-177-DOTA-
NOC and Lu-177-DOTA-TATE in $69$ patients suffering from a NET with high
somatostatin receptor expression. Absorbed doses were assessed using the MIRD
scheme. DOTA-TATE showed a higher uptake in tumor lesions as well as a higher
mean absorbed dose. However, the high inter-patient variability deemed an
individualized dosimetry obligatory.

Individualized dosimetry in patients undergoing therapy with Lu-177-DOTA-D-
Phe-Ryr-octreotate was studied in [21]. SPECT data were analyzed considering
whole organ volumes of interest (VOIs) to obtain the amount of radioactivity
in the organ. In addition, small VOIs were employed to estimate the
radioactivity concentration in the respective tissues. With a cohort of $24$
patients SPECT images were acquired at $1$h, $24$h, $96$h, and $168$h after
infusion of the radio-pharmaceutical and absorbed doses were calculated in
non-tumor affected organs like kidney, spleen and liver. Authors concluded
that SPECT dosimetry based on small VOIs has been more reliable and delivered
more accurate results than conventional planar dosimetry or dosimetry based on
whole organ VOIs.

In a later study of Lu-177-labeled peptides for NET therapy, [12] observed
that side effects are few and remain mild, if kidney-protective agents were
used, thus increasing the progression-free survival to more than $40$ months.
Apart from kidneys, bone marrow is another organ at risk in radionuclide
therapy.

In [22] an individualized dosimetry protocol for the bone marrow was
developed. Furthermore, patients, undergoing fractionated therapy with $7.4$
GBq/cycle, were identified who first reached an accumulated dose of either $2$
Gy to the bone marrow or $23$ Gy to the kidneys. Planar whole-body and
SPECT/CT images over the abdomen were acquired at $24$h, $96$h, and $168$h
after administering Lu-177-octreotate. In addition, blood samples were
collected six times during the first $24$h, with additional blood samples
taken at $96$h and $168$h from several patients. Also urine samples were taken
from several patients during the first $24h$. The absorbed radiation dose was
computed from blood and urine activity curves as well as an organ-based
analysis of whole body images. The absorbed dose to the kidney was calculated
from the pharmacokinetic data obtained from SPECT/CT. The study revealed that
the kidney was the dose limiting organ in $197$ out of $200$ patients. The
number of cycles until the limiting dose was reached showed a high variability
($2-10$ cycles), thus suggesting the need for an individualized radiation
treatment.

Another recent study by [10] proposed an automatized voxel-wise dosimetry tool
for radiotherapy based on multiple SPECT/CT images. The latter were acquired
at $4,24,72$ hours after Lu-177-DOTATATE PRRT for $17$ patients suffering from
advanced NETs. Deformable image registration and in-house developed algorithms
for pharmacokinetic uptake interpolation at voxel level were employed.
Absorbed energy dose was computed from cumulated activity based on MC-
determined voxel S-values and compared to results obtained using a
conventional MIRD protocol (OLINDA).

Dosimetry of organs and tumors is often avoided because of the additional
burden for patients. Therefore, in [7] an approximation has been advocated
which allows the calculation of absorbed energy doses from a single
measurement of the abdominal radioactivity distribution. In the study, the
activity kinetics were assessed in $29$ patients suffering from NETs or
meningioma after administering Lu-177-DOTATATE or -DOTATOC. Mono-or bi-
exponential decay curves were fit to SPECT/CT imaging data from kidneys,
livers, spleens and NET lesions. The time-integrated activity was compared to
an approximation given by the product of three factors: the time $t$ of a
single measurement, the expected reading at time $t$ and the factor $2/\ln 2$.
Tissue-specific deviations of the approximation from the time-integrated
activity were computed for several time points ranging from $t=24$h to
$t=144$h. Results were inconclusive in that best accordance was achieved for
$t=96$h, while a reduced accuracy resulted for shorter and longer time spans.

##  Chapter 2 Dataset

###  2.1 Patient collectiv

Table 2.1: Patient data concerning the DOTATOC radiotherapy. Here $m$ denotes the weight of each person, $A_{0}$ the initially injected radioactivity, $g$ the grading of the NET, $CgA_{1}$ denotes the CgA before therapy and $CgA_{2}$ after therapy, $\Delta t$ the time betweent the $CgA$ values, $n$ the number of therapy cycles and the sex is coded in $1$ for female and $0$ for male. Age | sex | $g$ | $m$ | $A_{0}$ | $CgA_{1}$ | $CgA_{2}$ | $\Delta t$ | $n$
---|---|---|---|---|---|---|---|---
$[y]$ |  |  | $(kg)$ | $(MBq)$ |  ($\mu g/l$) |  ($\mu g/l$) | $(m)$ |
45 | 1 | 3 | 64 | 6077 | $494$ | $1026$ | $5$ | $4$
66 | 0 | 3 | 79 | 6375 | $1121$ | $706$ | $4$ | $4$
54 | 0 | 1 | 77 | 7054 | $661$ | $1956$ | $4$ | $2$
67 | 1 | 2 | 70 | 6821 | $1352$ | $1473$ | $6$ | $3$
78 | 0 | 2 | 88 | 5773 | $34$ | $31$ | $4$ | $3$
54 | 0 | 2 | 56 | 7023 | $54$ | $77$ | $4$ | $3$
67 | 0 | 2 | 75 | 5809 | $597$ | $1748$ | $2$ | $2$
54 | 0 | 2 | 71 | 6644 | $1319$ | $682$ | $4$ | $4$
77 | 1 | 2 | 49 | 6378 | $224$ | $246$ | $1$ | $1$
52 | 0 | 3 | 78 | 6233 | $46$ | $83$ | $5$ | $2$
52 | 0 | 3 | 103 | 7265 | $76$ | $47$ | $4$ | $3$
71 | 0 | 2 | 118 | 6757 | $51$ | $51$ | $3$ | $2$
49 | 0 | 3 | 80 | 6328 | $139$ | $214$ | $2$ | $2$
$66.5$ |  | $2.3$ | $77.5$ | 6503 | $589.8$ | $768.1$ | $3.9$ | $2.8$

This study comprises $26$ patients suffering from NETs ($13$) or Prostate
Cancer (PC) ($13$) and which underwent a ${}^{177}Lu$-DOTATOC or
${}^{177}Lu$\- Prostate Specific Membrane Antigen (PSMA) therapy. The patient
cohort consists of $23$ male and $3$ female patients with an average age of
$63.8\pm 10$ years at time of therapy. The on average injected activity was
$A_{0}=6643\pm 421$ MBq.

The information about the patients suffering on NET is collected in table 2.1
and for the patients with prostate cancer in 2.2.

Table 2.2: Patient data concerning the PSMA radiotherapy. Here $m$ denotes the weight of each person, $A_{0}$ the initially injected radioactivity, $Gleason$ the Gleason Score of the PC, $PSA_{1}$ denotes the PSA level before therapy and $PSA_{2}$ after therapy, $\Delta t$ the time betweent the $PSA$ values and $n$ the number of therapy cycles. Age | Gleason | $m$ | $A_{0}$ | $PSA_{1}$ | $PSA_{2}$ | $\Delta t$ | $n$
---|---|---|---|---|---|---|---
$[y]$ |  | $(kg)$ | $(MBq)$ |  ($ng/ml$) |  ($ng/ml$) | $(m)$ |
80 | $9$ | 88 | 6699 | $1277$ | $594$ | $2$ | $1$
67 | $8$ | 104 | 6997 | $96$ | $120$ | $6$ | $2$
75 | $9$ | 74 | 6314 | $22$ | $5$ | $2$ | $1$
55 | $6$ | 113 | 5931 | $122$ | $523$ | $2$ | $1$
61 | $9$ | 85 | 6899 | $32$ | $74$ | $5$ | $2$
66 | $9$ | 97 | 7139 | $92$ | $199$ | $2$ | $1$
63 | $9$ | 94 | 7103 | $484$ | $289$ | $2$ | $1$
65 | $9$ | 80 | 7030 | $58$ | $0.1$ | $3$ | $2$
66 | $9$ | 86 | 6545 | $469$ | $445$ | $1$ | $1$
78 | $9$ | 106 | 6768 | $93$ | $86$ | $3$ | $2$
79 | $9$ | 75 | 6915 | $138$ | $25$ | $2$ | $1$
62 | $8$ | 68 | 6824 | $2251$ | $1946$ | $2$ | $1$
55 | $9$ | 94 | 7006 | $148$ | $53$ | $2$ | $2$
$65.7$ | $8.4$ | $90$ | $6782$ | $406.3$ | $336.0$ | $2.6$ | $1.4$

###  2.2 Image Acquisition

####  2.2.1  Single-Photon Emission Computed Tomography (SPECT)

Data was acquired on a Single-Photon Emission Computed Tomography / CT
(SPECT/CT) system (Siemens Symbia T2) at time points $4$h, $24$h, $48$h, and
$72$h after administering the radiopharmaceutical. The acquisition of the $24$
h p.i. image was carried out as part of a hybrid SPECT/CT acquisition, After
manual co-registration, the Computed Tomograpy (CT) part was used for
attenuation correction of the subsequent SPECT. Acquisition was done based on
an inhouse standard quantitative ${}^{177}Lu$ \- protocol, which is described
in detail in [20]. Therefore, it is herein only briefly outlined:

  * •

SPECT using medium energy collimators, $3^{\circ}$ angular sampling, $15$ min
total dwell time

  * •

Iterative ordered-subset conjugate-gradient reconstruction of the $208$ keV
photopeak data with $24$ iterations, $1$ subset, matrix $256$x$256$

  * •

Point-spread-function modelling in reconstruction

  * •

Triple energy window based scatter correction

  * •

CT-based attenuation correction

  * •

No post-reconstruction smoothing

The reconstruction, using the ordered-subset conjugate-gradient algorithm, was
carried out on a Siemens research workstation. The algorithm outputs fully
quantitative reconstructed SPECT images providing activities per voxel
measured in $[Bq/ml]$. Furthermore the SPECT data is alined with the
corresponding CT data.

####  2.2.2 CT

Additionally, a low-dose CT was carried out as part of the multi-modal
SPECT/CT acquisition to enable necessary attenuation corrections of the SPECT
images. The CT covered the same field-of-view as the SPECT and was acquired
and reconstructed using the following parameters:

  * •

Slice collimation of 2 x 5 mm, pitch of 1.8, time per rotation of 0.8 s, tube
voltage of 130 kVp, tube current of 30 mAs effective

  * •

Filter-Back-Projection reconstruction with B08s and B41s Kernels, $512$x$512$
matrix, $2.5$ mm slice thickness

  * •

B08s image was used for attenuation correction of the SPECT data

  * •

B41s image was used for defining organ and tumor volumes-of-interest

Reconstructed SPECT images were partitioned into voxels of size $V_{vox}=4.79\
mm^{3}$ and were cropped to $82^{3}$ voxels. Thereby large portions of the
reconstructed volume were removed which predominantly contained air.

###  2.3 Common Preprocessing

Figure 2.1: CT-image, SPECT-image and hyprid SPECT/CT-image as shown in ITK-
snap. In the right bottom corner the segmentation of the right and left
kidney, the spleen and a liver tumor can be seen.

In figure 2.1 a CT-image, a SPECT-image and a SPECT/CT-image are presented.
For all patients, the kidneys, spleen and tumors were segmented by a medical
expert. The tumors were sorted to bone lesion, lymph nodes, visceral and
primary tumor. The segmentation was done using the software ITK-snap [31],
which was carried out on the original CT datasets. The resulting segmentation
was subsequently down-sampled to match the voxel-size of $(4.7mm)^{3}$ of the
SPECT image. For each patient, the follwoing information is available:

  * •

SPECT-Image acquired 4h after injection

  * •

SPECT-Image acquired 24h after injection

  * •

SPECT-Image acquired 48h after injection

  * •

SPECT-Image acquired 72h after injection

  * •

CT-Image acquired 24h after injection

  * •

Organ and tumor segmentation

Exemplarly, for one patient and one slice the eight images are illustrated in
figure 2.2.

Figure 2.2: The four SPECT-images, the CT-image and the corresponding
segmentation map illustrated for one slice of one patient.

##  Chapter 3 Sequential Monte Carlo Techniques

Non-linear dynamical systems evolve either continuously or via transitions
between discrete states. The latter approximation is only reasonable if the
mean residence time in any state is longer than the transition time. Particle
filters represent discrete state space models which estimate the states
$\mathbf{X}(\mathbf{r},t)$ of a dynamical system from noisy and disturbed
observations $\mathbf{Z}(\mathbf{r},t)$. This is achieved by combining
available measured data $\mathbf{Z}(\mathbf{r},t)$ with prior knowledge about
the underlying physical phenomena driving the dyanmics and the measurement
process. The underlying physics allow to construct a state evolution model
while the known measurement process provides information to construct an
observation model. With these ingredients and observations of the system
dynamics, estimates of the unobservable dynamic state variables can be
produced sequentially.

###  3.1 Background

State evolution models are commonly based on a Markov model, while observation
models need to reflect the specifics of the physical measurement process.
Both, evolution and observation models are most often described in terms of
conditional probability density functions. Hence we have (throughout we use
the short hand notation
$\mathbf{z}(t_{k})\equiv\mathbf{z}_{k},\\{\mathbf{z}(t_{0}),\ldots,\mathbf{z}(t_{k})\\}\equiv\mathbf{z}_{0:k}$
etc.)

$\displaystyle\mathbf{x}_{k+1}$ $\displaystyle\sim$ $\displaystyle
p\left(\mathbf{x}_{k+1}|\mathbf{x}_{0:k}\right)\approx
p\left(\mathbf{x}_{k+1}|\mathbf{x}_{k}\right)$ (3.1)
$\displaystyle\mathbf{z}_{k}$ $\displaystyle\sim$ $\displaystyle
p\left(\mathbf{z}_{k}|\mathbf{x}_{0:k}\right)\approx
p\left(\mathbf{z}_{k}|\mathbf{x}_{k}\right)$ (3.2)

where $\sim$ means “represented through”, $\mathbf{x}_{k}$ denotes an
unobservable state vector of the dynamical system and $\mathbf{z}_{k}$ a
corresponding observation vector of the system through a measurement at
discrete time $t_{k}$ [15]. Note that the conditional probability densities
only depend on the last state of the system. This is in accord with the Markov
assumption which expresses common intuition that any prediction of a future
state should be based primarily on the recent history of the evolving system.

Now given the task that we want to predict unobservable states of an evolving
system from sequential noisy observations, Bayesian inference [17] allows to
integrate available information, i. e. observations $\mathbf{z}_{0:k}$ and an
associated physical model of the dynamical system, parametrized through
dynamic states $\mathbf{x}_{0:k}$. As new information $\mathbf{z}_{k}$ is
available, it can be combined with previous information
$p(\mathbf{x}_{0:k-1},\mathbf{z}_{0:k-1})$ via the Bayes Theorem to deduce the
underlying system state $\mathbf{x}_{k}$. Employing the Markov assumption we
have

$p(\mathbf{x}_{k}|\mathbf{z}_{0:k})=p(\mathbf{x}_{k}|\mathbf{z}_{k},\mathbf{z}_{0:k-1})\approx\frac{p(\mathbf{z}_{k}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{z}_{k-1})}{\sum_{k}\left(p(\mathbf{z}_{k}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{z}_{k-1})\right)}$
(3.3)

The latter forms the basis for statistical estimates of the posterior density
$p(\mathbf{x}_{k}|\mathbf{z}_{0:k})$ of the unobservable state vector
$\mathbf{x}_{k}$ given the likelihood of the observation
$p(\mathbf{z}_{k}|\mathbf{x}_{k})$, some prior knowledge
$p(\mathbf{x}_{k}|\mathbf{z}_{k-1})$ about the observation process and the
data evidence $p(\mathbf{z}_{k})$.

Note that the term posterior density denotes different densities depending on
the application:

  * •

Filtering: $p(\mathbf{x}_{k}|\mathbf{z}_{0:k})$

  * •

Prediction: $p(\mathbf{x}_{k+\tau}|\mathbf{z}_{0:k})$

  * •

Smoothing: $p(\mathbf{x}_{k-\tau}|\mathbf{z}_{0:k})$

where $\tau>0$ denotes the prediction or smoothing lag.

The most widely known Bayesian filter method is the Kalman filter [11], which,
however, relies on linear models with additive Gaussian noise, i. e. it
assumes the distribution $p(\mathbf{x}_{k}|\mathbf{z}_{0:k})$ to be a multi-
variate Gaussian. Extensions of the Kalman filter concern sequential Monte
Carlo methods [4] which represent the posterior density in terms of random
samples (particles) and associated weights. Such sequential Monte Carlo
methods are variously known as Particle Filter, bootstrap filter [5], survival
of the fittest [13], condensation algorithm [9] and can be applied to non-
linear models with non-Gaussian errors [2]. Thus, given any dynamical state
representation $f(\mathbf{x}_{k})$, one can estimate the related expected
value via

$\displaystyle\mathbb{E}(f(\mathbf{x}_{k}))$ $\displaystyle=$
$\displaystyle\int
f(\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{z}_{0:k})d\mathbf{x}_{k}$ (3.4)
$\displaystyle=$ $\displaystyle\frac{\int
f(\mathbf{x}_{k})p(\mathbf{z}_{k}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{z}_{k-1})d\mathbf{x}_{k}}{\int
p(\mathbf{z}_{k}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{z}_{k-1})d\mathbf{x}_{k}}$
$\displaystyle\approx$
$\displaystyle\sum_{l=1}^{L}w_{k}^{(l)}f(\mathbf{x}_{k}^{(l)})$

Here $\mathbf{x}_{k}^{(l)}$ denotes a set of samples, called particles, drawn
from the posterior distribution $p(\mathbf{x}_{k}|\mathbf{z}_{k-1})$. The
sampling weights are given by

$w_{k}^{(l)}=\frac{p\left(\mathbf{z}_{k}|\mathbf{x}_{k}^{(l)}\right)}{\sum_{l^{\prime}=1}^{L}p\left(\mathbf{z}_{k}|\mathbf{x}_{k}^{(l^{\prime})}\right)}$
(3.5)

where the same samples are used in the numerator and the denominator. Hence a
particle filter approximates the posterior distribution
$p(\mathbf{x}_{k}|\mathbf{z}_{0:k})$ by a set of particles
$\mathbf{x}_{k}^{(l)}$, sampled from $p(\mathbf{x}_{k}|\mathbf{z}_{k-1})$, and
their weights $\\{w_{k}^{(l)}\\}$. The latter are normalized to $0\leq
w_{k}^{(l)}\leq 1$ and obey the closure relation $\sum_{l}w_{k}^{(l)}=1$.

In practice one is concerned with a sequential sampling scheme, which works as
follows:

  * •

Suppose, a set of samples $\\{\mathbf{x}_{k}^{(l)}\\}$ and related weights
$w_{k}^{(l)}$ have been obtained at time $t_{k}$. These samples and weights
represent the posterior distribution $p(\mathbf{x}_{k}|\mathbf{z}_{k})$ at
time step $t_{k}$.

  * •

Suppose, an observation $\mathbf{z}_{k+1}$ has been made subsequently

How can one estimate corresponding samples $\mathbf{x}_{k+1}^{(l)}$ and
weights $w_{k+1}^{(l)}$?

  * •

First, new samples are drawn from the distribution

$\displaystyle p(\mathbf{x}_{k+1}|\mathbf{z}_{0:k})$ $\displaystyle=$
$\displaystyle\int
p(\mathbf{x}_{k+1}|\mathbf{x}_{k},\mathbf{z}_{0:k})p(\mathbf{x}_{k}|\mathbf{z}_{0:k})d\mathbf{x}_{k}$
(3.6) $\displaystyle=$ $\displaystyle\int
p(\mathbf{x}_{k+1}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{z}_{0:k})d\mathbf{x}_{k}$
$\displaystyle=$ $\displaystyle\int
p(\mathbf{x}_{k+1}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{z}_{k},\mathbf{z}_{0:k-1})d\mathbf{x}_{k}$
$\displaystyle=$ $\displaystyle\frac{\int
p(\mathbf{x}_{k+1}|\mathbf{x}_{k})p(\mathbf{z}_{k}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{z}_{k-1})d\mathbf{x}_{k}}{\int
p(\mathbf{z}_{k}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{z}_{k-1})d\mathbf{x}_{k}}$
$\displaystyle=$
$\displaystyle\sum_{l=1}^{L}w_{k}^{(l)}p\left(\mathbf{x}_{k+1}|\mathbf{x}_{k}^{(l)}\right)$

where the conditional independence property has been used

$\displaystyle p(\mathbf{x}_{k+1}|\mathbf{x}_{k},\mathbf{z}_{0:k})$
$\displaystyle=$ $\displaystyle p(\mathbf{x}_{k+1}|\mathbf{x}_{k})$ (3.7)
$\displaystyle p(\mathbf{z}_{k}|\mathbf{x}_{k},\mathbf{z}_{0:k-1})$
$\displaystyle=$ $\displaystyle p(\mathbf{z}_{k}|\mathbf{x}_{k})$ (3.8)

The above distribution is a mixture distribution. Hence, samples
$\mathbf{x}_{k}^{(l)}$ can be drawn by choosing a mixture component $l$ with
probability $w_{k}^{(l)}$ and drawing a sample from the corresponding mixture
component of the distribution.

  * •

Second, for each sample the new observation $\mathbf{z}_{k+1}$ is used to
estimate the related weights

$w_{k}^{(l)}\propto p(\mathbf{z}_{k+1}|\mathbf{x}_{k+1}^{(l)})$ (3.9)

Thus using a probabilistic description and the Markov assumption, one is
concerned with the following densities:

  * •

The likelihood of the state $p(\mathbf{x}_{k}|\mathbf{z}_{k-1})$

  * •

The Markov transition model $p(\mathbf{x}_{k+1}|\mathbf{x}_{k})$ $\rightarrow$
state evolution model

  * •

The sensor model $p(\mathbf{z}_{k}|\mathbf{x}_{k})$ $\rightarrow$ observation
model

###  3.2 State Estimation

While theoretical derivations are based on the general Markov model,
algorithms are based on evolution and observation models as detailed next.

####  3.2.1 The model equations

Let the dynamic state of a system be described by a state vector
$\mathbf{x}(t)\in\mathcal{R}^{N}$ whose components represent the dynamical
variables of the system. The dynamic evolution of the system, called the state
evolution model $\mathbf{f}(...)$, is described by the following equation for
the state vector

$\mathbf{x}(t_{k})\equiv\mathbf{x}_{k}=\mathbf{f}\left(\mathbf{x}_{k-1},\boldsymbol{\epsilon}^{(x)}\right)$
(3.10)

where $k=1,2,.......$ denotes discrete time instances $t_{k}$ and
$\boldsymbol{\epsilon}^{(x)}\in\mathcal{R}^{L}$ denotes the state noise.

If measurements are available at $t_{k}$, they are collected in the
observation vector $\mathbf{z}$, and its relation to the state vector is given
by the observation model $\mathbf{h}(...)$

$\mathbf{z}(t_{k})\equiv\mathbf{z}_{k}=\mathbf{h}\left(\mathbf{x}_{k},\boldsymbol{\epsilon}^{(m)}\right)$
(3.11)

where $\boldsymbol{\epsilon}^{(m)}\in\mathcal{R}^{M}$ denotes the measurement
noise. The various noise distributions $p(\boldsymbol{\epsilon}^{(x)})$ and
$p(\boldsymbol{\epsilon}^{(m)})$ are assumed to be known. In some
applications, also additional known control variables $\mathbf{u}_{k}$ enter
the problem, but are omitted here.

The state estimation problem aims at obtaining information about the state
vector $\mathbf{x}_{k}$ based on the state evolution model $\mathbf{f}(...)$
and based on the measurements $\mathbf{z}_{1:k}$, given by the observation
model $\mathbf{h}(...)$.

####  3.2.2 The simplifying assumptions

The state evolution and observation models rely on the following simplifying
assumptions:

  * •

The sequence of state vectors $\mathbf{x}_{k}$ forms a Markov process [3] and
does not depend on the sequence of observations

$p(\mathbf{x}_{k}|\mathbf{x}_{0:k-1},\mathbf{z}_{1:k-1})=p(\mathbf{x}_{k}|\mathbf{x}_{k-1})$
(3.12)

  * •

The observations $\mathbf{z}_{k}$ form a Markov process with respect to the
history of the state vector $\mathbf{x}_{k}$, i. e.

$p(\mathbf{z}_{k}|\mathbf{x}_{0},\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{k-1})=p(\mathbf{z}_{k}|\mathbf{x}_{k})$
(3.13)

  * •

The probability of an observation only depends on the current state

$p(\mathbf{z}_{k}|\mathbf{x}_{0:k},\mathbf{z}_{1:k-1})=p(\mathbf{z}_{k}|\mathbf{x}_{k})$
(3.14)

  * •

For $i\neq j$ the noise vectors $\boldsymbol{\epsilon}^{(x)}_{i}$ and
$\boldsymbol{\epsilon}^{(x)}_{j}$, as well as
$\boldsymbol{\epsilon}^{(m)}_{i}$ and $\boldsymbol{\epsilon}^{(m)}_{j}$, are
mutually independent

  * •

The noise vectors $\boldsymbol{\epsilon}^{(x)}_{i}$ and
$\boldsymbol{\epsilon}^{(m)}_{j}$ are mutually independent for all
$i,j=1,2,....$

  * •

All noise vectors are also mutually independent of the initial state
$\mathbf{x}_{0}$

####  3.2.3 The problems to be solved

Several problems can be tackled with this approach

  * •

The prediction problem determines $p(\mathbf{x}_{k}|\mathbf{z}_{1:k-1})$, i.
e. it predicts the state of the system at time $t_{k}$ based on observations
made at times $t_{1},\ldots,t_{k-1}$

  * •

The filtering problem determines $p(\mathbf{x}_{k}|\mathbf{z}_{1:k})$, i. e.
it filters out from the observations made up to time point $t_{k}$ the state
of the system at time $t_{k}$.

  * •

The fixed lag smoothing problem determines
$p(\mathbf{x}_{k}|\mathbf{z}_{1:k+\tau})$, i. e. it determines the state of
the system at time $t_{k}$ given observations over a time span
$t_{1},\ldots,t_{k+\tau}$, where $\tau\geq 1$ denotes the fixed time lag.

  * •

The whole domain smoothing problem determines
$p(\mathbf{x}_{k}|\mathbf{z}_{1:K})$, where
$\mathbf{z}_{1:K}=\mathbf{z}_{i},i=1,\ldots,K$ denotes the complete
observation sequence.

In this study, only the filtering problem is of interest here. By assuming
that $p(\mathbf{x}_{0}|\mathbf{z}_{0})=p(\mathbf{x}_{0})$ is available, the
posterior probability density $p(\mathbf{x}_{k}|\mathbf{z}_{1:k})$ is then
obtained with Bayesian filters in two steps: prediction and update.

###  3.3 Bayesian Filters

The Bayesian solution to compute the posterior distribution,
$p(\mathbf{x}_{k}|\mathbf{z}_{1:k})$, of the state vector $\mathbf{x}(t_{k})$,
given past observations $\mathbf{z}_{1:k}$, is given by the general Bayesian
update recursion

$\displaystyle\mbox{measurement
update}:p\left(\mathbf{x}_{k}|\mathbf{z}_{1:k}\right)$ $\displaystyle=$
$\displaystyle\frac{p(\mathbf{z}_{k}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{z}_{1:k-1})}{p(\mathbf{z}_{k}|\mathbf{z}_{1:k-1})}$
(3.15)
$\displaystyle\mbox{normalization}:p(\mathbf{z}_{k}|\mathbf{z}_{1:k-1})$
$\displaystyle=$
$\displaystyle\int_{\mathds{R}^{N}}p(\mathbf{z}_{k}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{z}_{1:k-1})d\mathbf{x}_{k}$
(3.16) $\displaystyle\mbox{time update}:p(\mathbf{x}_{k+1}|\mathbf{z}_{1:k})$
$\displaystyle=$
$\displaystyle\int_{\mathds{R}^{N}}p(\mathbf{x}_{k+1}|\mathbf{x}_{k})p(\mathbf{x}_{k}|\mathbf{z}_{1:k})d\mathbf{x}_{k}$
(3.17)

The posterior distribution is the primary output from a nonlinear filter, from
which standard measures as the minimum mean square estimate
$\hat{\mathbf{x}}_{k|k}^{mms}$ and its covariance $\mathbf{C}_{k|k}^{(x),mms}$
can be extracted:

$\displaystyle\hat{\mathbf{x}}_{k|k}^{mms}$ $\displaystyle=$
$\displaystyle\int\mathbf{x}_{k}p(\mathbf{x}_{k}|\mathbf{z}_{1:k})d\mathbf{x}_{k}$
(3.18) $\displaystyle\mathbf{C}_{k|k}^{(x),mms}$ $\displaystyle=$
$\displaystyle\int\left(\mathbf{x}_{k}-\hat{\mathbf{x}}_{k}^{mms}\right)\left(\mathbf{x}_{k}-\hat{\mathbf{x}}_{k}^{mms}\right)^{T}p(\mathbf{x}_{k}|\mathbf{z}_{1:k})d\mathbf{x}_{k}$
(3.19)

For a linear Gaussian model, the Kalman filter recursion relations provide a
solution to this Bayesian filtering problem. But for non-linear or non-
Gaussian models no finite-dimensional representations of the posterior density
exist, and numerical approximations are needed.

###  3.4 Particle Filter

If linear Gaussian evolution - observation models are inadequate, the use of
the Kalman filter does not result in optimal solutions. In such cases the
posterior density of underlying states $\mathbf{x}_{k}$, given the sequence of
observations $\mathbf{z}_{0:k}$ is not analytic, implying that the
distribution $p(\mathbf{x}_{k}|\mathbf{z}_{0:k})$ cannot be expressed in a
simple form. The application of Monte Carlo techniques then appears as the
most general and robust approach to analyze non-linear system dynamics and/or
non-Gaussian state distributions. This is the case despite the availability of
the so-called extended or unscented Kalman filter, which generally involves a
linearization of the problem.

The particle filter (PF) is also known under various names as there are:

  * •

the bootstrap filter,

  * •

the condensation algorithm,

  * •

interacting particle approximations,

  * •

survival of the fittest.

The particle filter approximates the posterior distribution of unobservable
states, $p(\mathbf{x}_{k}|\mathbf{z}_{0:k})$, with a discrete density which
will be evaluated at a dynamic stochastic grid in this work.

####  3.4.1 Sampling on a dynamic stochastic grid

The Particle Filter (PF) is based on a direct application of the Bayesian
recursion relations and approximates the posterior distribution,
$p(\mathbf{x}_{k}|\mathbf{z}_{0:k})$, with a discrete density sampled on a
dynamic stochastic grid. This kind of sampling has some noteworthy properties:

  * •

A dynamic stochastic grid $\mathbf{x}^{i}_{k}$ changes over time and
represents a very efficient representation of the state space.

  * •

The PF generates and evaluates a set $\\{\mathbf{x}^{i}_{1:k}\\}_{i=1}^{N}$ of
$N$ different trajectories rather than current states $\mathbf{x}^{i}_{k}$
only. This affects the time update of the Bayesian recursion relations as
follows:

$\displaystyle p\left(\mathbf{x}^{i}_{1:k+1}|\mathbf{z}_{1:k}\right)$
$\displaystyle=$ $\displaystyle
p\left(\mathbf{x}^{i}_{1:k}|\mathbf{z}_{1:k}\right)p\left(\mathbf{x}^{i}_{k+1}|\mathbf{x}^{i}_{1:k},\mathbf{z}_{1:k}\right)$
(3.20) $\displaystyle=$ $\displaystyle
w^{i}_{k|k}p\left(\mathbf{x}^{i}_{k+1}|\mathbf{x}^{i}_{k}\right)$

  * •

The new grid is obtained by sampling from

$p\left(\mathbf{x}_{k+1}|\mathbf{z}_{1:k}\right)=\sum_{i=1}^{N}p\left(\mathbf{x}^{i}_{1:k+1}|\mathbf{z}_{1:k}\right)=\sum_{i=1}^{N}w^{i}_{k|k}p\left(\mathbf{x}^{i}_{k+1}|\mathbf{x}^{i}_{k}\right)$
(3.21)

Employing importance sampling amounts to introducing an auxiliary proposal
density $q\left(\mathbf{x}_{k+1}|\mathbf{x}_{k},\mathbf{z}_{k+1}\right)$, from
which it is easy to sample. The advantage of using an auxiliary proposal
distribution can be seen, if one re-writes the Bayesian recursion as follows:

$\displaystyle p\left(\mathbf{x}_{k+1}|\mathbf{z}_{1:k}\right)$
$\displaystyle=$
$\displaystyle\int_{\mathcal{R}^{N}}p\left(\mathbf{x}_{k+1}|\mathbf{x}_{k}\right)p\left(\mathbf{x}_{k}|\mathbf{z}_{1:k}\right)d\mathbf{x}_{k}$
$\displaystyle=$
$\displaystyle\int_{\mathcal{R}^{N}}q\left(\mathbf{x}_{k+1}|\mathbf{x}_{k},\mathbf{z}_{k+1}\right)\frac{p\left(\mathbf{x}_{k+1}|\mathbf{x}_{k}\right)}{q\left(\mathbf{x}_{k+1}|\mathbf{x}_{k},\mathbf{z}_{k+1}\right)}p\left(\mathbf{x}_{k}|\mathbf{z}_{1:k}\right)d\mathbf{x}_{k}$

Now for each particle a random sample from this auxiliary proposal
distribution $\mathbf{x}^{i}_{k+1}\sim
q\left(\mathbf{x}_{k+1}|\mathbf{x}^{i}_{k},\mathbf{z}_{k+1}\right)$ is drawn,
and the posterior probability is adjusted for each particle with the
importance weight

$\displaystyle p\left(\mathbf{x}_{1:k+1}|\mathbf{z}_{1:k}\right)$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\frac{p\left(\mathbf{x}^{i}_{k+1}|\mathbf{x}^{i}_{k}\right)}{q\left(\mathbf{x}^{i}_{k+1}|\mathbf{x}^{i}_{k},\mathbf{z}_{k+1}\right)}w^{i}_{k|k}\delta\left(\mathbf{x}_{1:k+1}-\mathbf{x}^{i}_{1:k+1}\right)$
(3.23) $\displaystyle=$ $\displaystyle
w^{i}_{k+1|k}\delta\left(\mathbf{x}_{1:k+1}-\mathbf{x}^{i}_{1:k+1}\right)$

Note that the proposal distribution
$q\left(\mathbf{x}^{i}_{k+1}|\mathbf{x}^{i}_{k},\mathbf{z}_{k+1}\right)$
depends on the last state in the particle trajectory $\mathbf{x}^{i}_{k}$ as
well as the next measurement $\mathbf{z}_{k+1}$. The simplest choice of
proposal distribution is to use the dynamic model itself, i. e.
$q\left(\mathbf{x}^{i}_{k+1}|\mathbf{x}^{i}_{k},\mathbf{z}_{k+1}\right)=p\left(\mathbf{x}^{i}_{k+1}|\mathbf{x}^{i}_{k}\right)$,
leading to $w^{i}_{k+1|k}=w^{i}_{k|k}$.

Thus, the Particle Filter is a sequential Monte Carlo technique for the
solution of the state estimation problem. It is a sampling method which
approximates the posterior distribution by making use of its temporal
structure. Again, the key idea is to represent the required posterior density
function by a set of random samples (particles) with associated weights
yielding a particle representation of the posterior density according to

$p(\mathbf{x}_{k}|\mathbf{z}_{0:k})\approx\sum_{i=1}^{N}w^{i}_{k-1}\delta\left(\mathbf{x}_{k}-\mathbf{x}^{i}_{k-1}\right)$
(3.24)

where $w^{i}_{k}$ represents the weight of particle $\mathbf{x}^{i}_{k}$.
Estimates are then computed based on these samples and weights. As the number
of samples becomes very large, this Monte Carlo characterization becomes an
equivalent representation of the posterior probability function, and the
solution approaches the optimal Bayesian estimate.

Note that for the posterior density a recursive representation exists
according to

$p\left(\mathbf{x}_{k}|\mathbf{z}_{0:k}\right)=\alpha\cdot
p\left(\mathbf{z}_{k}|\mathbf{x}_{k}\right)\int
d\mathbf{x}_{k-1}p\left(\mathbf{x}_{k-1}|\mathbf{z}_{0:k-1}\right)p\left(\mathbf{x}_{k}|\mathbf{x}_{k-1}\right)$
(3.25)

With a particle representation, this recursive relation simplifies to

$p\left(\mathbf{x}_{k}|\mathbf{z}_{0:k}\right)\approx\alpha\cdot
p\left(\mathbf{z}_{k}|\mathbf{x}_{k}\right)\sum_{i=1}^{N}w^{i}_{k-1}p\left(\mathbf{x}_{k}|\mathbf{x}^{i}_{k-1}\right)$
(3.26)

How do we create a proper set of particles for representing the posterior
distribution $p\left(\mathbf{x}_{k}|\mathbf{z}_{0:k}\right)$? The answer is
importance sampling which will be explained next.

####  3.4.2 Sequential importance sampling

The technique of importance sampling is a method for generating fair samples
of a distribution $p(\mathbf{x}(t))$. Suppose $p(\mathbf{x}(t))$ is a density
from which it is difficult to draw samples, but it is easy to evaluate
$p(\mathbf{x}^{i}(t))$ for some particular instances $\mathbf{x}^{i}$, i. e.
on a grid. Then, an approximation to $p(\mathbf{x})$ can be given by:

$\displaystyle p(\mathbf{x})$ $\displaystyle\approx$
$\displaystyle\sum_{i=1}^{N}w^{i}\delta\left(\mathbf{x}-\mathbf{x}^{i}\right)$
(3.27) $\displaystyle w^{i}$ $\displaystyle=$
$\displaystyle\frac{p(\mathbf{x})}{q\left(\mathbf{x}^{i}\right)}$ (3.28)

where $q(\mathbf{x}^{i})$ denotes any auxiliary proposal distribution, also
called importance density. In particular, a uniform sampling of the state
space could be used but would lead to sample depletion, or sample degeneracy,
or sample impoverishment. A more direct auxiliary proposal distribution would
be an approximation to the posterior density
$p\left(\mathbf{x}_{k}|\mathbf{z}_{1:k}\right)$. This is achieved by re-
sampling, which introduces the required information feedback from the
observations, so trajectories that perform well, will survive the re-sampling.

####  3.4.3 A generic Particle Filter algorithm

  * •

Choose

    * –

a proposal distribution
$q\left(\mathbf{x}_{k+1}|\mathbf{x}_{1:k},\mathbf{z}_{k+1}\right)$,

    * –

a re-sampling strategy and

    * –

the number of particles $N$.

  * •

Initialization: Generate $\mathbf{x}^{i}_{1}\sim p_{x_{0}},\;i=1,\ldots,N$ and
set $w^{i}_{1|0}=\frac{1}{N}$.

  * •

Iteration: for $k=1,2,......$ do

    * –

Measurement update: For $i=1,2,\ldots,N$ compute

$w^{i}_{k|k}=\frac{w^{i}_{k|k-1}p\left(\mathbf{z}_{k}|\mathbf{x}^{i}_{k}\right)}{\sum_{i=1}^{N}w^{i}_{k|k-1}p\left(\mathbf{z}_{k}|\mathbf{x}^{i}_{k}\right)}$
(3.29)

    * –

Estimation: The filtering posterior density is approximated by

$\hat{p}\left(\mathbf{x}_{1:k}|\mathbf{z}_{1:k}\right)=\sum_{i=1}^{N}w^{i}_{k|k}\delta\left(\mathbf{x}_{1:k}-\mathbf{x}^{i}_{1:k}\right)$
(3.30)

and the mean is approximated by

$\langle\mathbf{x}\rangle\approx\sum_{i=1}^{N}w^{i}_{k|k}\mathbf{x}^{i}_{1:k}$
(3.31)

    * –

Re-sampling: Optionally at each time, take $N$ samples with replacement from
the set $\\{\mathbf{x}^{i}_{1:k}\\}_{i=1}^{N}$ where the probability to take
sample $i$ is $w^{i}_{k|k}=1/N$

    * –

Time update: Generate predictions according to the chosen proposal density

$\mathbf{x}^{i}_{k+1}\sim
q\left(\mathbf{x}_{k+1}|\mathbf{x}^{i}_{k},\mathbf{z}_{k+1}\right)$ (3.32)

and compensate for the importance weight

$w^{i}_{k+1|k}=w^{i}_{k|k}\frac{p\left(\mathbf{x}^{i}_{k+1}|\mathbf{x}^{i}_{k}\right)}{q\left(\mathbf{x}^{i}_{k+1}|\mathbf{x}^{i}_{k},\mathbf{z}_{k+1}\right)}$
(3.33)

Note that the algorithm exhibits the fundamental structure of the Bayesian
recursion relations as detailed above. Most common forms of PF algorithms
combine the weight updates into one equation according to

$w^{i}_{k|k}\propto
w^{i}_{k-1|k-1}\frac{p\left(\mathbf{z}_{k}|\mathbf{x}^{i}_{k}\right)p\left(\mathbf{x}^{i}_{k}|\mathbf{x}^{i}_{k-1}\right)}{q\left(\mathbf{x}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{z}_{k}\right)}$
(3.34)

The PF algorithm outputs an approximation of the trajectory posterior density
$p\left(\mathbf{x}_{1:k}|\mathbf{z}_{1:k}\right)$. For a filtering problem,
the simplest solution would be just to extract the last state
$\mathbf{x}^{i}_{k}$ from the trajectory $\mathbf{x}^{i}_{k}$ and use the
particle approximation

$\hat{p}\left(\mathbf{x}_{k}|\mathbf{z}_{1:k}\right)=\sum_{i=1}^{N}w^{i}_{k|k}\delta\left(\mathbf{x}_{k}-\mathbf{x}^{i}_{k}\right)$
(3.35)

However, this is incorrect as in general all paths $\mathbf{x}^{j}_{1:k-1}$
can lead to the state $\mathbf{x}^{i}_{k}$. The correct solution, taking all
paths leading to $\mathbf{x}^{i}_{k}$ into account, leads to an importance
weight

$w^{i}_{k+1|k}=\sum_{j=1}^{N}w^{j}_{k|k}\frac{p\left(\mathbf{x}^{i}_{k+1}|\mathbf{x}^{j}_{k}\right)}{q\left(\mathbf{x}^{i}_{k+1}|\mathbf{x}^{i}_{k},\mathbf{z}_{k+1}\right)}$
(3.36)

which replaces the time update of the normalized weights given in the generic
PF algorithm above. This solution is called the marginal Particle Filter (mPF)
which has many applications in system identification and robotics.
Unfortunately, the complexity is now of order $O(N^{2})$.

Prediction to get $p\left(\mathbf{x}_{1:k+m}|\mathbf{z}_{1:k}\right)$ can be
implemented by repeating the time update in the generic PF algorithm $m$
times.

####  3.4.4 Importance re-sampling

The basic particle filter suffers from sample depletion where all but a few
particles will have negligible weights. Re-sampling solves this problem, but
inevitably destroys information and thus increases uncertainty by the random
sampling. In this work the bootstrap PF, also called Sampling Importance
Resampling (SIR), is employed which applies re-sampling each time. To avoid
sample depletion, an auxiliary importance density needs to be selected
appropriately as the prior density
$p\left(\mathbf{x}_{k}|\mathbf{x}^{i}_{k-1}\right)$.

####  3.4.5 The SIR algorithm

The SIR algorithm uses re-sampling at every iteration. It can be summarized in
the following steps, as applied to the system evolution from $t_{k-1}$ to
$t_{k}$:

  * Step 1

    * –

For $i=1,\ldots,N$ draw new particles $\mathbf{x}^{i}_{k}$ from the importance
density by employing the transition model

$q(\mathbf{x}_{k})=\sum_{i=1}^{N}w^{i}_{k-1}p\left(\mathbf{x}_{k}|\mathbf{x}^{i}_{k-1}\right)$
(3.37)

To do so, choose a random number $r$ uniformly from $[0,1]$ and choose
particle $i=r$, then sample from the prior density
$p\left(\mathbf{x}_{k}|\mathbf{x}^{i}_{k-1}\right)$.

    * –

Use the corresponding likelihood to calculate corresponding weights

$w^{i}_{k}=p\left(\mathbf{z}_{k}|\mathbf{x}^{i}_{k}\right)$ (3.38)

The samples $\mathbf{x}^{i}_{k}$, employed above, are fair samples from
$p\left(\mathbf{x}_{k}|\mathbf{z}_{0:k-1}\right)$ and re-weighting them
accounts for the evidence of the observations $\mathbf{z}_{k}$.

  * Step 2

Calculate the total weight $W_{k}=\sum_{i=1}^{N}w^{i}_{k}$.
Normalize the particle weights $w^{i}_{k}=W_{k}^{-1}w^{i}_{k}\;\forall
i=1,....,N$

  * Step 3

Re-sample the particles by doing

    * –

Compute the cumulative sum of weights $W^{i}_{k}=W^{i}_{k-1}+w^{i}_{k};\forall
i=1,\ldots,N,\;W^{0}=0$

    * –

Let $i=1$ and draw a starting point $u_{1}$ from a uniform distribution
$U[0,N^{-1}]$

    * –

For $j=1,\ldots,N$ do the following

      * *

move along the cumulative sum of weights by setting $u_{j}=u_{1}+N^{-1}(j-1)$

      * *

while $u_{j}>W^{i}$ set $i=i+1$

      * *

assign samples $\mathbf{x}^{j}_{k}=\mathbf{x}^{i}_{k}$

      * *

assign weights $w^{j}_{k}=N^{-1}$

The resampling procedure just described avoids to have many degenerate
particles with vanishing weights but it also leads to a loss of diversity in
the sense that the resulting samples may contain many redundant particles.
This phenomenon is called sample impoverishment and is often observed in case
of small process noise. In this situation, all particles collapse to a single
particle within few instants $t_{k}$.

####  3.4.6 Effective number of samples

An indicator of the degree of depletion is the effective number of samples,
defined in terms of the coefficient of variation $\sigma_{\epsilon_{x}}$ as

$N_{eff}=\frac{N}{1+\sigma_{\epsilon_{x}}^{2}(w^{i}_{k|k})}$ (3.39)

The effective number of samples is thus at its maximum $N_{eff}=N$ when all
weights are equal $w^{i}_{k|k}=\frac{1}{N}$, and the lowest value it can
attain is $N_{eff}=1$, which occurs when $w^{i}_{k|k}=1$ with probability
$1/N$ and $w^{i}_{k|k}=0$ with probability $(N-1)/N$.

In practical applications this number could be approximated by

$\hat{N}_{eff}=\frac{1}{\sum_{i=1}^{N}\left(w^{i}_{k|k}\right)^{2}}$ (3.40)

when we have $1\leq\hat{N}_{eff}\leq N$. Again the upper bound $N_{eff}=N$ is
attained when all particles have the same weight, and the lower bound
$N_{eff}=1$ when all the probability mass is devoted to a single particle. The
resampling condition in the PF can now be defined as $N_{eff}<N_{th}$, and the
threshold can for instance be chosen as $N_{th}=\frac{2N}{3}$.

####  3.4.7 Choice of proposal/importance distribution

The choice of proposal distribution clearly influences the depletion problem.
The most general proposal distribution has the form
$q\left(\mathbf{x}_{1:k}|\mathbf{z}_{1:k}\right)$. This means that the whole
trajectory needs be sampled at each iteration, which in real-time applications
is not realistic. But the general proposal can be factorized as

$q\left(\mathbf{x}_{1:k}|\mathbf{z}_{1:k}\right)=q\left(\mathbf{x}_{k}|\mathbf{x}_{1:k-1},\mathbf{z}_{1:k}\right)q\left(\mathbf{x}_{1:k-1}|\mathbf{z}_{1:k}\right)$
(3.41)

The most common approximation in applications is to reduce the path
$\mathbf{x}_{1:k-1}$ and only sample the new state $\mathbf{x}_{k}$, so the
proposal distribution $q\left(\mathbf{x}_{1:k}|\mathbf{z}_{1:k}\right)$ is
replaced by
$q\left(\mathbf{x}_{k}|\mathbf{x}_{1:k-1},\mathbf{z}_{1:k}\right)$, which, due
to the Markov assumption, can be written as

$q\left(\mathbf{x}_{k}|\mathbf{x}_{1:k-1},\mathbf{z}_{1:k}\right)\approx
q\left(\mathbf{x}_{k}|\mathbf{x}_{k-1},\mathbf{z}_{k}\right).$ (3.42)

The approximate proposal density predicts good values of the current state
$\mathbf{x}_{k}$ only, not of the whole trajectory $\mathbf{x}_{1:k}$. For
further insight, one needs to discuss the dependence of this proposal
distribution on the Signal to Noise Ratio (SNR). Here, the SNR is defined as
the ratio of the maximal value of the likelihood and prior, respectively,

$SNR\propto\frac{max_{x_{k}}p(\mathbf{z}_{k}|\mathbf{x}_{k})}{p(\mathbf{x}_{k}|\mathbf{x}_{k-1})}$
(3.43)

Note that for a linear Gaussian model this yields

$SNR\propto\sqrt{\frac{\det(\mathbf{C}^{\epsilon_{x}})}{\det(\mathbf{C}^{\epsilon_{m}})}}$
(3.44)

Thus, the SNR is high if the measurement noise is small compared to the signal
noise. Given this measure, one could consider different sampling strategies,
but here only the one used in this work is explained.

Likelihood sampling: For medium or high SNR, samples are drawn from the
likelihood $p\left(\mathbf{z}_{k}|\mathbf{x}_{k}\right)$.

The proposal distribution
$q\left(\mathbf{x}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{z}_{k}\right)$ can be
factorized as follows

$\displaystyle
q\left(\mathbf{x}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{z}_{k}\right)\approx
p\left(\mathbf{x}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{z}_{k}\right)$
$\displaystyle=$ $\displaystyle
p\left(\mathbf{x}_{k}|\mathbf{x}^{i}_{k-1}\right)\frac{p\left(\mathbf{z}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{x}_{k}\right)}{p\left(\mathbf{z}_{k}|\mathbf{x}^{i}_{k-1}\right)}$
(3.45) $\displaystyle=$ $\displaystyle
p\left(\mathbf{x}_{k}|\mathbf{x}^{i}_{k-1}\right)\frac{p\left(\mathbf{z}_{k}|\mathbf{x}_{k}\right)}{p\left(\mathbf{z}_{k}|\mathbf{x}^{i}_{k-1}\right)}$

Now consider the case that the likelihood
$p\left(\mathbf{z}_{k}|\mathbf{x}_{k}\right)$ is much more peaky than the
prior $p(\mathbf{x}_{k}|\mathbf{x}_{k-1})$ and if it is integrable in state
space, $\mathbf{x}_{k}$, then one can set

$d^{2}(\mathbf{x}_{n},\mathbf{x}_{m})p\left(\mathbf{x}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{z}_{k}\right)\propto
p\left(\mathbf{z}_{k}|\mathbf{x}_{k}\right)$ (3.46)

Thus, a suitable proposal distribution for the high SNR case is based on a
scaled likelihood function

$q\left(\mathbf{x}_{k}|\mathbf{x}^{i}_{k-1},\mathbf{z}_{k}\right)\propto
p\left(\mathbf{z}_{k}|\mathbf{x}_{k}\right)$ (3.47)

This choice then yields the following weight update

$w^{i}_{k|k}=w^{i}_{k-1|k-1}p\left(\mathbf{x}^{i}_{k}|\mathbf{x}^{i}_{k-1}\right)$
(3.48)

Sampling from the likelihood requires that the likelihood function
$p\left(\mathbf{z}_{k}|\mathbf{x}_{k}\right)$ is integrable with respect to
$\mathbf{x}_{k}$. This is not the case when $N>L$, i. e. the number of
unobservable states is larger than the number of observations (measurements).
The interpretation in this case is that for each value of $\mathbf{z}_{k}$,
there is an infinite-dimensional manifold of possible $\mathbf{x}_{k}$ to
sample from, each one equally likely.

###  3.5 Resampling Schemes

The method of sampling from the likelihood has the drawback that it becomes
unstable as $k$ increases. There is a discrepancy between the weights,
therefore the algorithm can be stabilized by performing a re-sampling
sufficiently often. The weighted approximate density is modified by each re-
sampling step. This is done by eliminating particles with low weights and
multiplying particles with an important weight. The filtering posterior
density then becomes:

$\hat{p}\left(\mathbf{x}_{0:k},\mathbf{z}_{1:k}\right)=\sum_{i=1}^{N}\frac{n_{i}}{N}\delta(\mathbf{x}_{0:k}-\mathbf{x}_{0:k}^{i})$
(3.49)

where $n_{i}$ is the number of copies of particle $\mathbf{x}_{0:k}^{i}$. In
the following, the three main re-sampling schemes will be described.

####  3.5.1 Multinomial Resampling

The multinomial resampling scheme generates $N$ ordered uniform random numbers
on the interval $(0,1]$. The new particle set is selected according to the
multinomial distribution [1].

$\displaystyle u_{k}$ $\displaystyle=u_{k+1}\tilde{u}_{k}^{\frac{1}{k}},\
with\ \tilde{u}_{k}\sim U[0,1)$ $\displaystyle u_{N}$
$\displaystyle=\tilde{u}_{N}^{\frac{1}{N}}$ $\displaystyle\mathbf{x}_{k}$
$\displaystyle=x(F^{-1}(u_{k}))$ $\displaystyle=x_{i}\ with\ i\ so\ that\
u_{k}\in\big{[}\sum_{s=1}^{i-1}\omega_{s},\sum_{s=1}^{i}\omega_{s}\big{)},$

where the generalized inverse of the cumulative probability distribution is
$F^{-1}$. The weights are normalized.

####  3.5.2 Systematic resampling

With systematic re-sampling, the particle set is selected according to the
multinomial distribution of $N$ generated ordered numbers $u_{k}$.

$u_{k}=\frac{(k-1)+\tilde{u}}{N},\ with\ \tilde{u}\sim U[0,1)$ (3.50)

####  3.5.3 Residual resampling

With residual re-sampling, the variance can be decreased. We define:

$\displaystyle N^{i}=\lfloor n\omega^{i}\rfloor+\bar{N}^{i},\ \mbox{for}\
i=1,\ldots,m,$

where the $\bar{N}^{1},...,\bar{N}^{n}$ are multinomially distributed, and
$\lfloor\rfloor$ means the integer part. The multinomial distribution
$M(n-R;\bar{\omega}^{1},...,\bar{\omega}^{n})$ has the following arguments:

$\displaystyle R$ $\displaystyle=\sum_{i=1}^{m}\lfloor n\omega^{i}\rfloor$
$\displaystyle\bar{\omega}^{i}$ $\displaystyle=\frac{n\omega^{i}-\lfloor
n\omega^{i}\rfloor}{n-R},\ i=1,...,m$

The multinomial counts $\bar{N}^{1},...,\bar{N}^{n}$ are generated in the same
way as in the multinomial re-sampling described above.

##  Chapter 4 Time-Integrated Activity Map (TIA)

In this chapter the four SPECT-images as shown in figure 2.2 will be
integrated voxelwise to obtain the number of decays per voxel. The activity
values will be denoised by the particle filter. After the denoising, a mono-
exponential function is fit to the data and integrated over time.

###  4.1 Particle Filter Method

For dosimetry it is more important to have a patient specific model than an
on-line tracking method. Therefore all SPECT images were used to estimate a
state evolution model. First of all, Volume Of Interest (VOI)s were segmented
in a CT image by a medical expert. To this contour a safety margin of 2 voxels
was added to account for partial volume effects as well as deviations from the
registration or inter-organ movements. By summing up the activity values of
all voxels enclosed by the contour, an organ-specific total activity was
obtained for each time point. An additional time point was added at
$t_{k}=600\ h$, at which the specific activity was set to $0\ Bq/ml$
deliberately. The time-dependent activity (see figure 4.2) was expressed by a
parametrized analytic function which was fit to the observed time-dependent
activities according to

$\ln\left(\frac{a(t)}{a(0)}\right)=-\left(\frac{ln(2)}{T_{1/2}}\cdot
t\right)=-\frac{t}{\tau_{eff}}$ (4.1)

Note that the effective lifetime

$\tau_{eff}=\left(\frac{1}{\tau_{phys}}+\frac{1}{\tau_{biol}}\right)^{-1}=\frac{T_{1/2}}{\ln
2}$ (4.2)

Figure 4.1: Time activity curve for one kidney: $f(0)=1\%IA,T_{1/2}=77.9h$.

is composed of the physical lifetime of the radionuclide and its biological
lifetime due to the metabolic turnover of the radio-pharmaceutical. While the
biological half-life of the isotope should be identical across an organ, the
amplitude $a(0)$ of the Time Activity Curve (TAC) is different for every
voxel. Note that $a(0)$ is different from the initially injected activity
$A_{0}$ and generally unkown.

####  4.1.1 State evolution model (SEM)

For each voxel of the kidney an own State Evolution Model (SEM) is computed.
At the beginning, the four activity values of the same voxel plus an
additional value of $a(t_{k})=0Bq/ml$ at $t_{k}=600\ h$ are approximated by a
mono-exponential decay. Hereby, the half-life is assumed to be equal to the
previously computed half-life of the whole kidney. From this fit, the
amplitude $\langle a(0)\rangle$ of the state evolution model can be obtained.
Moreover, the inevitable noise on the data is modeled as Gaussian white noise
with an amplitude corresponding to $0.1\langle a(0)\rangle$.

####  4.1.2 Observation model (OM)

The Observation Model (OM) $\mathbf{h}$ represents the distribution of
acquisition times of the SPECT image and is modeled by a Normal distribution
around the estimated particle.

####  4.1.3 Similarity measure

The aim of the voxelwise integration of the activity is to obtain the number
of nuclear disintegrations $N^{Lu}_{n}$ for each voxel $n\in S$ according to:

$N^{Lu}_{n}=\int_{0}^{t_{K}}a_{n}(t^{\prime})dt^{\prime}\approx\sum_{k=1}^{K}a_{n}(t_{k})\cdot\Delta
t$ (4.3)

where $K=5$ in this study. The activity distribution is assumed to be
identical for each voxel, hence it can be estimated for each organ at every
time point. This fact can be explored to find a good parameter configuration.
As similarity measure the Kullback-Leibler divergence (KL) was used. An
averaged histogram of the $k=4$ time points was calculated, for which the
histograms were centered. A histogram of the voxel activities $h(\langle
a_{n}\rangle)$ was calculated, where $\langle
a\rangle_{n}=(1/K)\sum_{k}a_{n}(t_{k})$ denotes the activity of voxel $n$
averaged over the four time points. The similarity between different
histograms has been computed employing the KL.

Figure 4.2: Shown are two histogramms of voxelactivities: Left side with a
Kullback-Leibler divergence of $0.67$ and on the right side with $0.12$.

###  4.2 Simple Fit Method

Instead of fitting the activity values as smoothed by the particle filter
(’particle fit’), one can use another approach to calculate a voxel-wise decay
map. One can fit a mono-exponential function to the four activity values per
voxel. In the follwoing this method is called ’simple fit’.

For each method the Time Integrated Activity (TIA) in a specific source region
$S$ was calculated by integrating, for every voxel, the fit function over time
and sum over all voxels.

$\displaystyle TIA_{S}$ $\displaystyle=$ $\displaystyle A_{0}\sum_{n\in
S}\int_{0}^{\infty}\cdot e^{-\frac{t}{\tau_{eff,n}}}dt$

where $S$ is the source region, $n$ indicates one voxel, $A_{0}$ denotes the
injected activity and $\tau_{eff,n}$ the half-life, which can differ between
the voxels for both, the ’simple fit’ and ’particle fit’ method. Note that the
approximation $a_{n}(t_{0})\approx A_{0}\forall\;n\in S$ accounts for the fact
that the complex metabolic processes which distribute the radioactive
pharmaceutical in the body is generally unknown. Finally, also the
distribution of half-lives across all voxels in the VOIs was calculated.

##  Chapter 5 Results

###  5.1 Parameter evaluation

Finding a good parameter configuration, which yields similar histograms of the
number of decays within an organ, and which also yields a good match to the
measured activity distribution of all activity values within a VOI, the
standard deviations for the state evolution model and the observation model
were considered adjustable parameters. The configurations were computed for
all $26$ patients, and for each configuration, the histogram of the number of
decays per voxel was determined and compared with the averaged activity
histogram. The KL was computed for all configurations. The $\sigma_{obs}$ was
varied between $0.35\leq\sigma_{obs}\leq 0.99$. The considered standard
deviations are rather large because the estimated precision of the activity
measurement through SPECT imaging is in this range according to [19].

In the following table, the KLs, averaged over all patients, are collected:

Table 5.1: Kullback - Leibler divergences, averaged over the entire patient cohort model$\diagdown$ observation | 0.35 | 0.5 | 0.75 | 0.99
---|---|---|---|---
0.01 | 0.12 | 0.17 | 0.17 | 0.17
0.10 | 0.16 | 0.17 | 0.17 | 0.17
0.25 | 0.16 | 0.16 | 0.16 | 0.16
0.35 | 0.58 | 0.54 | 0.48 | 0.46

The values in each row are very similar to each other, indicating that the KL
is rather robust against uncertanties in the observation model. The smallest
value for the KL is marked in blue.

Table 5.2: Relative deviations of the sum of voxelwise integrated activities from the corresponding whole organ integrated activity in percent model $\diagdown$ observation | 0.35 | 0.5 | 0.75 | 0.99
---|---|---|---|---
0.01 | $-34.2\%$ | $-32.6\%$ | $-32.6\%$ | ${\color[rgb]{0,0,1}-32.6\%}$
0.10 | $-35.8\%$ | $-33.7\%$ | $-33.4\%$ | $-33.4\%$
0.25 | $-44.4\%$ | $-40.3\%$ | $-39.0\%$ | $-38.5\%$
0.35 | $-64.4\%$ | $-59.3\%$ | $-56.3\%$ | $-55.9\%$

As a second decision criterion for an optimal parameter configuration, the
relative deviation of the total number of decays estimated either by
integrating the time-dependent activity at every voxel and then summing over
all voxels within the VOI or by integrating the whole organ (here the kidney)
activity over time was calculated.

In the following table the resulting relative deviations are collected. The
deviations for the total number of decays are large because of a high
statistical uncertainty in the voxel activity values. To correct for these
large deviations, the voxelwise integrated activity map can be normalized with
the total number of decays per organ. For the following results, we choose the
parameters as $\sigma_{observation}=0.35$ and $\sigma_{model}=0.10$.

To obtain stable results for the single activity values, the number of
particles for each time step was varied. In figure 5.1 the determined
difference between the total number of decays per organ for different numbers
of particles and the total number of decays for $2000$ particles is
illustrated. With only $150$ particles, the difference to the reference value
almost disappears.

Figure 5.1: Variation of the number of particles.

Moreover, while applying the particle filter, different resampling schemes
were evaluated, but there was no difference in the total number of decays per
organ. Therefore we decided to perform no resampling. The possibility to loose
too much diversity is very low with only four time steps.

###  5.2 Voxelwise TAC

For each voxel, a state evolution model was defined through four measured time
points. According to the state evolution model and the observation model,
particles were generated and an estimation of the true activity value was
obtained through its likelihood. An analytic function was then adapted to
either the four estimates of the true activity values or the observed activity
values. In figure 5.2, the actual measurements and the estimates of the state
evolution model as represented by properly adapted analytic functions, whereby
the adaptation to the observed values is drawn in black and the adaptation to
the estimates of the true activity values is colored in grey. Through particle
filtering, the outlier at the second time point could be eliminated very well.

Figure 5.2: Time activity curves for two different voxels at four time steps.

###  5.3 Whole patient

In figure 2.2, one slice of the SPECT images at four time points for one
patient have been illustrated. Now the voxelwise results for the integration
of the TACs are shown in figure 5.3. The red voxels in the right image,
showing the results from the simple fit method, are those with a diverging
integral, the calculated number of decays is equal to infinity.

Figure 5.3: One slice of the results for the particle filtering method (left)
and the simple fit method (right).

The voxelwise integration of the particle filtered activity values was
performed for the whole patient. For the segmented organs (kidney, spleen and
tumor), a corresponding half-life was calculated and used for the state
evolution model. For the voxels outside the VOIs, an over all voxels from the
same tissue class averaged half-life was calculated and used for the model.
For each VOI, a histogram of the number of decays per voxel was obtained and
is shown in figure 5.4. A second histogram is also shown, were the number of
decays per voxel is presented, which is calculated by a simple voxel-by-voxel
integration. Both histograms are almost indistinguishable.

Figure 5.4: Histograms of number of decays per voxel for the particle filter
method (blue) and the voxelwise simple fit of the four activity values (red).

Considering the histograms of all voxel-specific half-lives, which were used
in the analytic function for the determination of the TIA, the particle filter
method yielded histograms centered at the half-life from the state evolution
model, while the simple fit method very often yielded half-lives larger than
the physical half-life of Lutetium (see figure 5.5). Also the widths of the
two distributions differ considerably. The much narrower width of the
distribution for the particle filter method depends on the $\sigma$ parameters
of the SEM and the OM.

Figure 5.5: Histograms of halflives per voxel for the particle filter method
(blue) and the voxelwise simple fit of the four activity values (red). The
halflife for the whole organ is $81.6h$.

##  Chapter 6 Discussion

In this work, we proposed a method for processing voxel-wise TACs prior to
fitting them to mono-exponentials and subsequently integrating them
temporally. The TACs were obtained by sequential, quantitative SPECT/CT
imaging for dosimetry of $26$ patients under Lu-177 targeted radiotherapy. The
proposed pre-processing method is a Particle Filter, which belongs to the
class of sequential Monte-Carlo methods. The aim of the method is a de-noising
of the observed TACs. The goals were the following:

  * •

The whole-organ TIA map, as obtained on basis of the filtered data, should be
as similar as possible to the conventional evaluation. The latter obtains TIA
maps by fitting and integrating the organ-averaged time-activity curve.

  * •

The histograms of the TIA maps, resulting from the PF, should resemble the
TIA-histograms, as initially estimated from SPECT images.

We have shown that Particle Filters are suitable for a reduction in the
variance of the effective half-lives of voxels in kidneys, spleen, and tumors.
Most notably, the number of voxels which exhibit half-lives longer than the
physical half-life of the radioisotope is significantly reduced.

For two voxel-specific methods to estimate whole organ TIAs, namely particle
fit (PFF) and simple fit (SF), we quantified deviations to whole-organ TIAs
estimated by time-integration of the whole organ activity.

When comparing the resulting half-lives, we found that application of particle
filtering results in a strong reduction of variance, compared to SF. As
intended, the voxel-specific half-lives for PFF were distributed around the
half-life used in the state evolution model. Importantly, the number of voxels
which exhibit a longer-than-physical half-life was greatly reduced with PFF,
compared to SF. In summary, the PFF method does not fix the half-life to a
single value, rather it exhibits a narrow distribution and reduces the number
of implausible half-lives to almost zero in our study.

Unfortunately, literature for comparing different methods for fitting voxel-
specific TACs is scarce so far. Thus, no generally accepted standard-method
exists. Guidelines recommend either not modelling at all, rather performing
trapezoidal integration of TACs, or modelling by exponential functions [28,
26], although these recommendations are usually in the context of whole-organ
dosimetry.

In the literature of voxel-based dosimetry in targeted radiotherapy, most
groups applied no modelling at all for the time period covered by imaging, and
mono-exponential extrapolation after the last imaging [25, 18, 30, 6, 23].
Jackson et al. [10] used tri-exponential modelling of the TAC, but applied
regularization in order to limit the number of voxels with implausible TIA. In
the study of Kost et al. [14], the user could choose between mono- or bi-
exponential modelling of pharmacokinetics. However, the focus of their study
was not on the comparison of the different models and consequently they did
not report on this.

Marcatili et al. [16] implemented an algorithm which automatically selects the
best model for each voxel, either pure mono-exponential decay or linear uptake
followed by mono-exponential decay, and no modelling and trapezoidal
integration as fall-back if modelling fails. They did not provide a systematic
comparison between the different models but stated that the fall-back was used
in less than $5\%$ of voxels for their application (I-124 Positron-Emission
Tomography (PET)). However, it could be hypothesized that for SPECT imaging of
typical Lu-177 activities, the image noise would be comparably higher than for
PET, which potentially could affect fitting of per-voxel TACs negatively.

To the best of our knowledge, the study of Sarrut et al. [24] is the only
study that systematically compared different models for voxel-wise TAC
fitting. They focused on an algorithm called VoMM, that automatically selected
the best model for fitting the TAC of individual voxels. The models considered
were mono- or bi-exponential functions with variable degrees-of-freedom of the
fit ranging from 2 to 4 parameters. They compared the results obtained by VoMM
and by using fixed models for all voxels in a numerical experiment with
varying noise levels and also in one dosimetry patient. In any case, 6 imaging
time points were used. For the numerical experiment, they evaluated the Root
Mean Squared Error between ground-truth and modelled TIAs and also the
percentage of voxels which could not be fit successfully ($\%$FF). They found
that the Root Mean Squared Error and $\%$FF was lower for VoMM as compared to
the fixed models. Notably, they also showed that a low count number and thus
high noise level in a voxel’s TAC would result in an increased Root Mean
Squared Error and $\%$FF. For their patient data, they found that $\%$FF was
lower for the VoMM than for most fixed models. However, one fixed model did
have an even lower $\%$FF.

Altogether, literature on this topic shows that voxel-wise fitting of TACs,
subsequent calculation of TIAs, and thus voxel-specific dosimetry is feasible.
However, it also turns out that especially the TAC fitting is most likely
negatively affected by noise. Most groups apply some technique to either
circumvent modelling of especially the early part of the TACs by using
numerical integration or implement some kind of regularization in order to
limit the bias introduced by voxel noise. Most groups did only assess if
voxel-wise fitting could be done successfully from an algorithmic point of
view, but did not evaluate if a successful fit resulted in feasible TIAs.
Based on our experience and from our results obtained by the simple fit
method, we suspect that many voxels’ TIAs, although obtained by successful
fitting, are still incorrect. This problem could be alleviated by application
of Particle Filters, which could be understood as de-noising based on a
plausible assumption/model. Altogether, we think that particle filters offer
an interesting alternative to methods reported from literature and should be
further evaluated in forthcoming studies.

###  6.1 Conclusion

In this study, the challenge was to calculate the number of decays per voxel
from the SPECT-images. Therefore, a statistical method called particle filter
was used. With this method the voxelwise noise of the SPECT images could be
removed by matching the measuremenst to a physical model of the radioactive
decay. This method yields a homogeneous distribution of the decays per voxel
by taking into account different organs. The particle filter method was
compared to a simple fit method. The latter generally resulted in
inhomogeneous decay maps where certain voxels contained diverging numbers of
decays. For $26$ patients suffering from a prostate cancer or a neuroendocrine
tumor, the corresponding decay maps could be determined reliably.

## Bibliography

  * [1] R. Douc and O. Cappe.  Comparison of resampling schemes for particle filtering.  pages 64 – 69, 2005.
  * [2] A. Doucet, N. de Freitas, and N. (Eds.) Gordon.  Sequential Monte Carlo Methods in Practice.  Springer, 2001.
  * [3] A. Eberle.  Markoc Processes, 2017.
  * [4] G.S. Fishman.  Monte Carlo - Concepts, Algorithms and Applications.  Springer, 1996.
  * [5] N.J. Gordon, D.J. Salmond, and A.F.M. Smith.  Novel approach to nonlinear/non-gaussianbayesian state estimation.  Radar and Signal Processing, IEEE Proceedings, 140(2):107–113, 1993\.
  * [6] E. Grassi, F. Fioroni, V. Ferri, E. Mezzenga, M.A. Sarti, T. Paulus, N. Lanconelli, A. Filice, A. Versari, and M. Iori.  Quantitative comparison between the commercial software stratos® by philips and a homemade software for voxel-dosimetry in radiopeptide therapy.  Physica Medica: European Journal of Medical Physics, 31(1):72–79, 2015.
  * [7] Heribert Hänscheid, Constantin Lapa, Andreas K Buck, Michael Lassmann, and Rudolf A Werner.  Dose mapping after endoradiotherapy with 177lu-dotatate/dotatoc by a single measurement after 4 days.  Journal of nuclear medicine: official publication, Society of Nuclear Medicine, 59(1):75–81, 2018.
  * [8] Michael S Hofman, John Violet, Rodney J Hicks, Justin Ferdinandus, Sue Ping Thang, Tim Akhurst, Amir Iravani, Grace Kong, Aravind Ravi Kumar, Declan G Murphy, Peter Eu, Price Jackson, Mark Scalzo, Scott G Williams, and Shahneen Sandhu.  [177lu]-psma-617 radionuclide treatment in patients with metastatic castration-resistant prostate cancer (lupsma trial): a single-centre, single-arm, phase 2 study.  The Lancet Oncology, 19(6):825 – 833, 2018.
  * [9] M. Isard and A. Blake.  CONDENSATION– conditional density propagation of visual tracking.  International Journal of Computer Vision, 29(1):5–28, 2001.
  * [10] P. Jackson, J.M. Beauregard, M.S. Hofman, T. Kron, A. Hogg, and R.J. Hicks.  An automated voxelized dosimetry tool for radionuclide therapy based on serial quantitative spect/ct imaging.  Medical physics, 40(11), 2013.
  * [11] R.E. Kalman.  A new approach to linear filtering and prediction problems.  J. Basic Engineering, pages 35–45, 1960.
  * [12] BLR Kam, JJM Teunissen, Eric P Krenning, Wouter W de Herder, Saima Khan, EI Van Vliet, and Dirk Jan Kwekkeboom.  Lutetium-labelled peptides for therapy of neuroendocrine tumours.  European journal of nuclear medicine and molecular imaging, 39(1):103–112, 2012.
  * [13] K. Kanazawa, D. Koller, and S. Russell.  Stochastic simulation algorithms for dynamic probabilistic networks.  In Proceedings of the Eleventh Conference on Uncertainty in Artificial Intelligence, UAI’95, pages 346–351, San Francisco, CA, USA, 1995\. Morgan Kaufmann Publishers Inc.
  * [14] S. D. Kost, Y.K. Dewaraja, R. G. Abramson, and M. G. Stabin.  Vida: A voxel-based dosimetry method for targeted radionuclide therapy using geant4.  Cancer Biotherapy and Radiopharmaceuticals, 30(1):16–26, 2015.  PMID: 25594357.
  * [15] D.J.C. MacKay.  Information Theory, Inference and Learning Algorithms.  Cambridge University Press, 2003.
  * [16] S. Marcatili, C. Pettinato, S. Daniels, G. Lewis, P. Edwards, S. Fanti, and E. Spezi.  Development and validation of raydose: a geant4-based application for molecular radiotherapy.  Physics in Medicine and Biology, 58(8):2491, 2013.
  * [17] A. O’Hagan and J. Forster.  Kendall’s Advanced Theory od Statistics, Vol. 2B: Bayesian Inference.  Arnold Publishers, 1999.
  * [18] A. Prideaux, H. Song, R. Hobbs, B. He, E. Frey, P. Ladenson, R. Wahl, and Sgouros G.  3D radiobiologic dosimetry: application of radiobiologic modelling to patient-specific 3D imaging-based internal dosimetry.  J. Nucl. Med,, 48:1008–1016, 2007.
  * [19] P. Ritt, J. Hornegger, and T. Kuwert.  Technik und physikalische aspekte der spect/ct.  Der Nuklearmediziner, 34(01):9–20, 2011.
  * [20] J.C. Sanders, T. Kuwert, J. Hornegger, and Ph. Ritt.  Quantitative SPECT/CT imaging of 177-lu with In Vivo Validation in Patients Undergoing Peptide Receptor Radionuclide Therapy.  Mol. Imaging Biol, 2014.
  * [21] Mattias Sandström, Ulrike Garske, Dan Granberg, Anders Sundin, and Hans Lundqvist.  Individualized dosimetry in patients undergoing therapy with 177 lu-dota-d-phe 1-tyr 3-octreotate.  European journal of nuclear medicine and molecular imaging, 37(2):212–225, 2010.
  * [22] Mattias Sandström, Ulrike Garske-Román, Dan Granberg, Silvia Johansson, Charles Widström, Barbro Eriksson, Anders Sundin, Hans Lundqvist, and Mark Lubberink.  Individualized dosimetry of kidney and bone marrow in patients undergoing 177lu-dota-octreotate treatment.  Journal of Nuclear Medicine, 54(1):33–41, 2013.
  * [23] D. Sarrut, J.N. Badel, A. Halty, G. Garin, D. Perol, P. Cassier, J.Y. Blay, D. Kryza, and A.L. Giraudet.  3d absorbed dose distribution estimated by monte carlo simulation in radionuclide therapy with a monoclonal antibody targeting synovial sarcoma.  EJNMMI Physics, 4(1):6, Jan 2017.
  * [24] D. Sarrut, A. Halty, J.N. Badel, L. Ferrer, and M. Bardis.  Voxel-based multimodel fitting method for modeling time activity curves in spect images.  Medical Physics, 44(12):6280–6288, 2017.
  * [25] G. Sgouros, K. S. Kolbert, A. Sheikh, K. S. Pentlow, et al.  Patient-specific dosimetry for 131-i thyroid cancer therapy using 124-i pet and 3-dimensional-internal dosimetry (3d-id) software.  The Journal of Nuclear Medicine, 45(8):1366, 2004.
  * [26] J.A. Siegel, S.R. Thomas, J.B. Stubbs, M.G. Stabin, M.T. Hays, K.F. Koral, J.S. Robertson, R.W. Howell, D. R. Wessels, B. W .and Fisher, et al.  Mird pamphlet no. 16: techniques for quantitative radiopharmaceutical biodistribution data acquisition and analysis for use in human radiation dose estimates.  Journal of Nuclear Medicine, 40(2):37S, 1999.
  * [27] Jonathan Strosberg, Ghassan El-Haddad, Edward Wolin, Andrew Hendifar, James Yao, Beth Chasen, Erik Mittra, Pamela L. Kunz, Matthew H. Kulke, Heather Jacene, David Bushnell, Thomas M. O Dorisio, Richard P. Baum, Harshad R. Kulkarni, Martyn Caplin, Rachida Lebtahi, Timothy Hobday, Ebrahim Delpassand, Eric Van Cutsem, Al Benson, Rajaventhan Srirajaskanthan, Marianne Pavel, Jaime Mora, Jordan Berlin, Enrique Grande, Nicholas Reed, Ettore Seregni, Kjell Ã–berg, Maribel Lopera Sierra, Paola Santoro, Thomas Thevenet, Jack L. Erion, Philippe Ruszniewski, Dik Kwekkeboom, and Eric Krenning.  Phase 3 trial of 177lu-dotatate for midgut neuroendocrine tumors.  New England Journal of Medicine, 376(2):125–135, 2017.  PMID: 28076709.
  * [28] P. Vicini, A. B. Brill, and A. Stabin, M.and Rescigno.  Kinetic modeling in support of radionuclide dose assessment.  Seminars in nuclear medicine, 38(5):335–346, 2008.
  * [29] Christiane Wehrmann, Stefan Senftleben, Carolin Zachert, Dirk Müller, and Richard P Baum.  Results of individual patient dosimetry in peptide receptor radionuclide therapy with 177lu dota-tate and 177lu dota-noc.  Cancer biotherapy & radiopharmaceuticals, 22(3):406–416, 2007\.
  * [30] S.J. Wilderman, A.M. Avram, J. Kritzman, R. Ackerman, and Y.K. Dewaraja.  Dosimetry in 131i internal emitter therapy using voxel dependent integrated time-activities derived from multiple, registered spect and ct images.  In Nuclear Science Symposium Conference Record, 2006. IEEE, volume 6, pages 3492–3496. IEEE, 2006.
  * [31] P.A. Yushkevich, J. Piven, H. Cody Hazlett, R. Gimpel Smith, S. Ho, J.C. Gee, and G. Gerig.  User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability.  Neuroimage, 31(3):1116–1128, 2006.
  * [32] John J. Zaknun, L. Bodei, J. Mueller-Brand, M. E. Pavel, R. P. Baum, D. Hörsch, M. S. O’Dorisio, T. M. O’Dorisiol, J. R. Howe, M. Cremonesi, and D. J. Kwekkeboom.  The joint iaea, eanm, and snmmi practical guidance on peptide receptor radionuclide therapy (prrnt) in neuroendocrine tumours.  European Journal of Nuclear Medicine and Molecular Imaging, 40(5):800–816, May 2013.

  *[EEMD]: Ensemble Empirical Mode Decomposition
  *[NET]: Neuroendocrine tumor
  *[CgA]: Chromogranin A level
  *[PC]: Prostate Cancer
  *[PSMA]: Prostate Specific Membrane Antigen
  *[PSA]: Prostate Specific Antigen
  *[SPECT]: Single-Photon Emission Computed Tomography
  *[SPECT/CT]: Single-Photon Emission Computed Tomography / CT
  *[CT]: Computed Tomograpy
  *[PF]: Particle Filter
  *[SIR]: Sampling Importance Resampling
  *[SNR]: Signal to Noise Ratio
  *[VOI]: Volume Of Interest
  *[TAC]: Time Activity Curve
  *[SEM]: State Evolution Model
  *[OM]: Observation Model
  *[KL]: Kullback-Leibler divergence
  *[TIA]: Time Integrated Activity
  *[PET]: Positron-Emission Tomography