Generating patient-specific virtual tumor populations with reaction-diffusion models and molecular imaging data

1.   Introduction

Broadly speaking, mathematical oncology seeks to develop a collection of in silico models that accurately predict the growth and progression of malignant tumors and their response to treatment ^[1]. While the basic scientific goal is to develop a rigorous quantitative science of cancer ^[2], a more readily attainable goal is to use existing models and data to inform clinical practice in some substantive way ^[3]. One approach to clinical mathematical oncology, which we call mathematical model-based precision medicine, is to use mathematical models to make patient-specific prognostic and therapeutic predictions which are subsequently employed to direct treatment choices and predict patient outcomes. In our view, such predictions must be conditioned on patient data in order for the strategy to qualify as precision medicine ^[4,5,6].

Two essential features of oncology are heterogeneity and uncertainty. Heterogeneity can manifest as interpatient (between-patient), intrapatient (within-patient) and intratumor (within-tumor); see ^[7] for a more complete discussion of the clinical impact of heterogeneity in cancer and ^[8] for a review of the mathematical strategies to address heterogeneity. From a modeling perspective, heterogeneity implies a statistical approach where model parameters are treated as spatially distributed random variables with both population and patient-specific probability distributions ^[5]. Heterogeneity also highlights the issue of uncertainty: For a given patient, most if not all model parameters are unknown prior to the acquisition of patient-specific data, and even with such data, uncertainties may remain, since model parameters may only be partially identifiable and the data acquisition system may introduce additional irreducible noise. Image data provides a clear example of this issue, as imaging introduces information loss due to detector noise, low resolution and the possibility of null functions ^[9]. This uncertainty must be accounted for if model-based predictions are to be used reliably in the clinical context.

Once a relevant mathematical model has been selected, prediction takes place in silico, whereby a simulated patient is produced and Quantities-of-Interest (QoIs) are calculated. In this context, we say that the parameter values and resulting simulations correspond to a virtual patient. Owing to the heterogeneity and uncertainty issues discussed above, in both population and patient-specific applications of mathematical oncology it is desirable to generate suitably randomized Virtual Patient Populations (VPPs). If the VPP represents an untreated group, it is called a Virtual Control Population (VCP) ^[4,10], while if the VPP represents a treatment group, it is called a Virtual Treatment Population (VTP). With appropriate VPPs, an in silico Virtual Clinical Trial (VCT) can be performed to assess treatment choices, again in either the population or patient-specific case. Population VPPs and VCTs can assist in making general prognosis and treatment recommendations or simulate a regulatory trial ^[11], while a patient-specific VPP and/or VCT can be used to assess uncertainty in a particular patient's disease progression or treatment response ^[6,12]. While population VCTs have found usage in device design and regulatory contexts ^[13], patient-specific VPPs and VCTs have seen limited usage.

In previous work ^[5,6], we have highlighted the application of spatiotemporal random field modeling to address heterogeneity and uncertainty. In this work, we discuss specific VPP-based strategies towards the practical implementation of the framework presented in ^[5,6]. Specifically, we describe a method for generating VCPs by pairing a spatiotemporal reaction-diffusion equation (RDE) model for avascular, pre-metastatic tumor growth with spatially inhomogenous random field coefficients. We discuss the RDE model in sections 2.1 and 2.2 and the random field models in section 2.3. Together, these models provide a solution of the 'forward problem' for VPP generation. Two important inverse problems then arise. In the first, we consider in sections 4.1 and 4.2 the problem of estimating the RDE coefficient fields for a specific patient from noisy molecular imaging data; proof-of-concept simulation results are given in sections 5.1 and 5.2. For the second inverse problem, we must estimate the statistical parameters that define the random field coefficient models from a collection of many patients' data, hence solving the statistical calibration problem for virtual populations; we discuss this problem in section 4.3. Both maximum likelihood and Bayesian techniques are presented, and both require an accurate likelihood model connecting model parameters to available noisy data. To this end, in section 3 we discuss statistical models for molecular Emission Computed Tomography (ECT) data. We make this emphasis both because well-validated likelihood models for ECT are available ^[9], and because it is possible to model a direct connection between physiological processes such as tumor cell density and spatial growth rate, and certain emission imaging techniques and tracers.

Note that the results given in sections 5.1 and 5.2 are generated in a pure simulation setting, where both ground-truth tumors and image data are simulated. This allows for methodological comparisons with known ground truth and tight control of all nuisance parameters that influence the data. A summary of future work that extends these results is provided in section 6.

2.   Tumor growth modeling with heterogeneity and uncertainty

In this section, we introduce the tumor growth model that we employ for the remainder of the paper. We emphasize that the models we present are simplifications of real tumor biology, which is complex and not yet fully understood. These models are intended to be more 'clinically relevant' ^[3] than biologically accurate, due to their relative simplicity and amenability to patient-specific inference from spatiotemporal imaging data. In section 2.2 we discuss the parameter-to-solution map, in section 2.3 we discuss random field models for the tumor growth model coefficients, and in section 2.4 discuss the notions of virtual populations and virtual clinical trials.

2.1.   Reaction-diffusion tumor growth models

A wide variety of Partial Differential Equation (PDE) models have been employed to model spatiotemporal tumor growth, and the derivation and mathematical analysis of such models is adequately described elsewhere ^{[14,15,16,17,18,19,20,21]}. Here, we are interested in computational techniques for generating virtual patient populations, and thus consider a relatively simple PDE model (2.1) that has gained traction in the context of GlioBlastoma Multiforme (GBM). The randomization and statistical calibration techniques we present are independent of the tumor growth model, however, so the choice of Eq (2.1) is largely for illustration purposes.

We denote a spatiotemporal density of tumor cells by $\mathit{\boldsymbol{n}} = n(\mathit{\boldsymbol{x}}, t)$, with $\mathit{\boldsymbol{x}}\in V$ being a $d-$dimensional position vector in the spatial domain $V$ and $0\leq t\leq T$, where $T$ is the maximum number of days simulated. The units of $\mathit{\boldsymbol{n}}$ are cells per unit volume, and we take $(\mathit{\boldsymbol{x}}, t)$ to have units of cm and days, respectively. For in vivo tumors $d = 3$, though for simulation purposes and modeling of in vitro and certain in vivo experimental setups we also consider spatial domains in $d = 2$ as a reasonable approximation. The general reaction-diffusion equation for $\mathit{\boldsymbol{n}}$ is

The spatially varying scalar diffusion coefficient $\mathit{\boldsymbol{D}} = D(\mathit{\boldsymbol{x}})$ describes the rate of apparent cell diffusion, while the growth function $\mathit{\boldsymbol{\mathcal{G}}}(\mathit{\boldsymbol{n}})$ describes the growth and competition between the tumor and its environment. The boundary condition is taken to be the no-flux (Neumann) condition $\mathit{\boldsymbol{\hat{\mathit{\boldsymbol{n}}u}}}\cdot\nabla\mathit{\boldsymbol{n}} = 0$, where $\mathit{\boldsymbol{\hat{\mathit{\boldsymbol{n}}u}}}$ is the outward boundary normal of the domain $V$. Note that a treatment function $\mathit{\boldsymbol{\mathcal{T}}}(\mathit{\boldsymbol{n}})$ can be added to Eq (2.1) to model an intervention such as chemotherapy, radiotherapy or surgical resection, but we postpone discussion of such models to a future work.

While in the scalar tumor growth context the Gompertz function $\mathit{\boldsymbol{\mathcal{G}}}(\mathit{\boldsymbol{n}}) = -\mathit{\boldsymbol{\rho}}\mathit{\boldsymbol{n}}\ln(\mathit{\boldsymbol{n}}/\mathit{\boldsymbol{\kappa}})$ has gained plausible mechanistic support via the notion of quiescence ^[22], to our knowledge the spatial RDE growth function is largely a phenomenological choice. The Fisher-KPP growth function is frequently selected in the GBM context, and is the one we will use in this paper:

The growth $\rho(\mathit{\boldsymbol{x}})$ and carrying capacity $\kappa(\mathit{\boldsymbol{x}})$ functions, which have respective units of cells per day and cells, will be discussed further below.

2.2.   Patient-specific parameters and the parameter-to-solution map

Note that in Eqs (2.1) and (2.2), several additional parameter functions are required to fully specify the solution to the PDE. These include the diffusion coefficient $D(\mathit{\boldsymbol{x}})$, the initial cell density $n_0(\mathit{\boldsymbol{x}})$, the growth $\rho(\mathit{\boldsymbol{x}})$ and carrying capacity $\kappa(\mathit{\boldsymbol{x}})$. These parameter functions are certainly patient-specific, since they correspond to local physiology and nutrient concentrations. For notational simplicity, we will use the symbol $\mathit{\boldsymbol{\beta}}$ to denote the collection of all parameters necessary to specify a solution to Eq (2.1). In the context of this paper, a virtual patient then corresponds to a particular choice of $\mathit{\boldsymbol{\beta}}$, which we will indicate by saying that patient $j$ has parameter $\mathit{\boldsymbol{\beta}}_j$. We will assume for the remainder of this article the initial condition is given by a Gaussian

with $\mathit{\boldsymbol{x}}_0 = [0.5, 0.5]^T$, $A = 5/2\pi\sigma^2$ and $\sigma^2 = 1$$e$$-4$; this represents an initial collection of 5 tumor cells, well-localized around the center of our computational domain. The collection of all remaining unknown parameter fields is thus $\mathit{\boldsymbol{\beta}} = (\mathit{\boldsymbol{D}}, \mathit{\boldsymbol{\rho}}, \mathit{\boldsymbol{\kappa}})$. To avoid pathological situations, we also assume the technical condition that $\mathit{\boldsymbol{\beta}}\in\mathbb{B}$, where $\mathbb{B}$ is a set of parameter functions selected to guarantee that a positive, bounded solution $n(\mathit{\boldsymbol{x}}, t)$ of Eq (2.1) exists for all $\mathit{\boldsymbol{\beta}}\in\mathbb{B}$ and for $0\leq t\leq T$. Selecting such a $\mathbb{B}$ requires careful technical analysis of Eq (2.1); see ^[23].

While we assume that each patient corresponds to a unique set of parameters $\mathit{\boldsymbol{\beta}} = \mathit{\boldsymbol{\beta}}_j$, it is likely that $\mathit{\boldsymbol{\beta}}_j$ is effectively unknown or uncertain for a given patient. To address this, we treat $\mathit{\boldsymbol{\beta}}$ as a function space-valued random vector in both the population and patient-specific cases, and provide corresponding random field models in section 2.3 below. In a simulation setting, we must choose a computational representation for $\mathit{\boldsymbol{\beta}}_j$, which we indicate with a superscript $v$ for 'virtual'. So, $\mathit{\boldsymbol{\beta}}_j^{(v)}$ is a virtual (i.e., computational) sample corresponding to a sample for $\mathit{\boldsymbol{\beta}}$.

Assuming Eq (2.1) is classically well-posed for every (or almost every) $\mathit{\boldsymbol{\beta}}\in\mathbb{B}$, we can construct the deterministic parameter-to-solution map

where $\mathit{\boldsymbol{n}}_j^{(model)} = n_j^{(model)}(\mathit{\boldsymbol{x}}, t)$ solves Eq (2.1) with the parameter vector $\mathit{\boldsymbol{\beta}} = \mathit{\boldsymbol{\beta}}_j$, corresponding to patient $j$. In a simulation setting, the forward map (2.4) is implemented via a numerical differential equation solver. We have employed a finite difference method-of-lines approach with a five-point spatial stencil and Runge-Kutta time stepper, as discussed in ^[24]. The result is an approximate solution to Eq (2.1), where the virtual coefficient sample $\mathit{\boldsymbol{\beta}}_j^{(v)}$ is processed through computational algorithm $\mathit{\boldsymbol{\mathcal{F}}}^{(v)}$ to produce the virtual tumor cell density $\mathit{\boldsymbol{n}}_j^{(v)}$:

The choice of discretization specifies the format for the virtual cell density $\mathit{\boldsymbol{n}}_j^{(v)}$; with our finite difference scheme, we can assume that $\mathit{\boldsymbol{n}}_j^{(v)}$ takes the form of an $N_x\times N_y\times N_t$ voxel array (for $d = 2$). The specification of $\mathit{\boldsymbol{\beta}}_j^{(v)}$ is discussed below; we will assume a voxel-free form (Eq (2.8) below) that can be evaluated on an arbitrary grid when needed. Note that the virtual cell density will differ both from a patient's true cell density and from a measured cell density by discretzation, model and measurement errors; see ^[8] for a discussion of these errors. A demonstration of our scheme implementing (2.5) in $d = 2$ is shown in Figure 1.

In the next section, we will introduce random field models that allow Eqs (2.4) and (2.5) to be used to generate virtual tumor populations. Note briefly that we have made assumptions that allow us to analyze and simulate Eq (2.1) as a Random PDE (RPDE) model, that is as a classically well-posed PDE whose coefficients are randomized. Note that this is in contrast to a Stochastic PDE (SPDE) model, which requires more advanced analysis and simulation techniques. See e.g., ^[5] for discussion of a technique based on characteristic functionals (introduced briefly below), and ^[25,26,27] for further discussion of the distinction between SPDE and RPDE.

2.3.   Random field models for RDE coefficients

As discussed in the introduction, any mathematical model that is employed in a clinical precision medicine context must account for the possibility of heterogeneity and uncertainty by assuming that model parameters are sampled from a probability distribution. Since the model parameters that specify Eqs (2.1) and (2.2) are spatial functions, we employ the theory of random fields. We will avoid technical discussions, sufficing to say that one can treat $\mathit{\boldsymbol{\beta}}$ as a function space-valued random vector, written formally as $\mathit{\boldsymbol{\beta}}:\Omega\rightarrow \mathbb{B}$, where $(\Omega, \mathscr{F}, \mathbb{P})$ is a probability space and $\mathbb{B}$ is a set of functions mentioned previously, selected to guarantee the existence of the forward map (2.4). We assume that $\mathbb{B}\subset\mathbb{X}$, where $\mathbb{X}$ is a Hilbert space with inner product $(\mathit{\boldsymbol{\beta}}, \mathit{\boldsymbol{\varphi }})_\mathbb{X}$, taken to be the $L^2(V)$ inner product for vector-valued functions. In the most general setting, the complete statistics of a second-order random process are fully specified by its marginal distributions (through the Kolmogorov construction), its abstract probability distribution $\mathbb{P}_\mathit{\boldsymbol{\beta}}$, or its characteristic functional: ^[5,9,28,29]

Frequently, a process is completely specified by its first few moments e.g., its mean function $\mathit{\boldsymbol{\bar{\beta}}}(\mathit{\boldsymbol{x}}) = \mathbb{E}[\mathit{\boldsymbol{\beta}}]$ and its covariance function $\mathit{\boldsymbol{k}}(\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{x}}^\prime) = \mathbb{E}[(\mathit{\boldsymbol{\beta}}(\mathit{\boldsymbol{x}}) -\mathit{\boldsymbol{\bar{\beta}}}(\mathit{\boldsymbol{x}}))(\mathit{\boldsymbol{\beta}}(\mathit{\boldsymbol{x}}^\prime) -\mathit{\boldsymbol{\bar{\beta}}}(\mathit{\boldsymbol{x}}^\prime))^{\text{T}}]$. While its use requires functional calculus, the characteristic functional (2.6) is typically a compact and efficient way to describe the complete statistics of a random field, and its use in the biological context is a topic of current investigation ^[5]. Refer to the recent books ^[27] and ^[30] for a complete discussion of the theory of Hilbert space-valued random vectors.

In this work, we will simulate random fields by first assuming a convenient functional form for the realizations of the process, then randomizing the corresponding (finitely many) parameters to achieve desired statistics. This method allows us to guarantee that realizations of $\mathit{\boldsymbol{\beta}}$ satisfy conditions for Eq (2.1) to be well-posed, and allows for straightforward simulations using standard code packages. In particular, the random field models we use will take the form of a randomized synthesis map, which for a scalar-valued function $\varphi _j(\mathit{\boldsymbol{x}})$ takes the form:

where ${\bf{\Phi }}(\cdot):\mathbb{R}^N\rightarrow \mathbb{X}$ is a synthesis map that maps a (finite-dimensional) patient-specific parameter $\mathit{\boldsymbol{\theta}}_j$ to the function $\varphi _j^{(v)}(\mathit{\boldsymbol{x}})$. By selecting a probability distribution $\mathbb{P}(\mathit{\boldsymbol{\theta}}_j\vert\mathit{\boldsymbol{\theta}}_p)$ that depends on the population parameter $\mathit{\boldsymbol{\theta}}_p\in{\bf{\Theta }}\subset\mathbb{R}^p$, the synthesis map (2.7) produces a random field. In sections 4.1 and 4.2, we discuss ways to estimate a particular $\mathit{\boldsymbol{\theta}}_j$ from image data, while in section 4.3 we discuss ways to estimate the population parameter $\mathit{\boldsymbol{\theta}}_p$ from a collection of images.

Because the coefficients in the RDE (2.1) must be non-negative and bounded to be physiologically realistic, it is advantageous to begin with a random field model that always produces non-negative, bounded realizations. For our simulations, we employ 'lumpy-type' random field models ^[31,32], a realization of which has the functional form

If the lump function $\ell(\mathit{\boldsymbol{x}}; \mathit{\boldsymbol{\gamma}})$ and the lump amplitudes $b_l$ are non-negative and bounded, and the number of lumps is finite, then Eq (2.8) is guaranteed to produce non-negative and bounded realizations. For simplicity, we use isotropic unnormalized Gaussian lumps:

Random fields with this choice of lump function were proposed to model soft tissues with smoothly varying characteristics ^[31], while more sophisticated choices of $\ell(\mathit{\boldsymbol{x}})$, such as the clustered lump technique ^[32], can model more complex tissues such as breast or brain tissue. Ultimately, the choice of $\ell(\mathit{\boldsymbol{x}})$ and its shape parameter $\mathit{\boldsymbol{\gamma}}$ requires a statistical calibration step such as that discussed in section 4.3.

We randomize Eq (2.8) as follows. We first select a lump shape parameter $\mathit{\boldsymbol{\gamma}} = \sigma^2 = \sigma_0^2$, a constant lump amplitude $b_0$, and a mean number of lumps $\bar{L}\ll L_{max}$. We then draw $L\sim {\text{Poi}}(\bar{L})$, and set $b_\ell = b_0$ for $1\leq \ell\leq L$, and $b_\ell = 0$ otherwise. Lastly, we draw $\mathit{\boldsymbol{x}}_1\ldots, \mathit{\boldsymbol{x}}_L$, I.I.D. uniform in the domain $V$ and form the sum (2.8). Note that the resulting random field realization can be written in terms of a nonlinear synthesis map (2.7) with $N = (d+1)L_{max}$. We use the notation $LB(\bar{L}, b_0, \sigma_0^2)$ to indicate a lumpy-type field with mean number of lumps $\bar{L}$, amplitude $b_0$ and isotropic lump variance $\sigma_0^2$. A demonstration of the behavior of lumpy-type random fields of the form given in Eq (2.8) is shown in Figure 1, top row.

We reiterate that a single realization of the random field specified in Eq (2.8) is fully specified by the parameters $\sigma^2_0, b_0$ and the lump center list $[\mathit{\boldsymbol{x}}_l]$, while the complete statistics of the random field are fully specified by $\bar{L}, \sigma^2_0$ and $b_0$. This distinction is important when we discuss parameter estimation problems in sections 4.1–4.3.

2.4.   Virtual populations and quantitative biomarkers

Suppose that $\mathit{\boldsymbol{q}}\in\mathbb{R}^n$ is a low-dimensional vector quantifying patient status or prognosis, for instance (informally) $\mathit{\boldsymbol{q}} = [\text{metric of tumor burden}, \text{metric of normal tissue function}]$. In the model-based precision medicine context, such a vector is generally called a quantitative biomarker ^[6,33,34]. The purpose of a virtual patient population, and more generally a virtual clinical trial, is to predict the distribution of $\mathit{\boldsymbol{q}}$ for either a control or treatment group. We will assume for convenience that $\mathit{\boldsymbol{q}} = q\in \mathbb{R}$ is a scalar functional of $\mathit{\boldsymbol{\beta}}$ and $\mathit{\boldsymbol{n}}$:

Given a virtual population of simulated parameters and tumors, $\left(\mathit{\boldsymbol{\beta}}_1^{(v)}, \mathit{\boldsymbol{n}}_1^{(v)}\right), \ldots, \left(\mathit{\boldsymbol{\beta}}_J^{(v)}, \mathit{\boldsymbol{n}}_J^{(v)}\right)$, we estimate the biomarker $q$ for each virtual patient using a method $\mathcal{M}^{(v)}$, resulting in a vector

If $\mathit{\boldsymbol{Q}}^{(v)}$ is generated from a population distribution, that is, the random vector $\mathit{\boldsymbol{\beta}}$ models inter-patient heterogeneity, then we say that $\mathit{\boldsymbol{Q}}^{(v)}$ is a population virtual biomarker sample. If $\mathit{\boldsymbol{\beta}}$ models patient-specific uncertainties, for instance if $\mathit{\boldsymbol{\beta}}$ is sampled from the conditional random vector $\mathit{\boldsymbol{\beta}}_j\vert \mathit{\boldsymbol{g}}_j$ where $\mathit{\boldsymbol{g}}_j$ is some patient-specific data, then we say that $\mathit{\boldsymbol{Q}}^{(v)}$ is a patient-specific virtual biomarker sample. Statistical properties of $q$ such as its mean, standard deviation and distribution can be estimated from $\mathit{\boldsymbol{Q}}^{(v)}$. By generating virtual control and treatment groups, Eq (2.11) can be used to simulate a clinical trial. If the VPP employed in Eq (2.11) represents a treatment group, decisions regarding the intervention can be made based on these statistics. For instance, a rigorous decision-theoretic approach to patient-specific treatment selection would define an additional function $U(\mathit{\boldsymbol{q}})$ which measures the utility of achieving the biomarker $\mathit{\boldsymbol{q}}(\mathit{\boldsymbol{\tau}})$ via the intervention $\mathit{\boldsymbol{\tau}}$. Since $\mathit{\boldsymbol{q}}$ is a random vector, the typical approach is to maximize the expected utility, that is, we seek the intervention $\mathit{\boldsymbol{\tau}}^*$ which solves

While numerous patient-specific therapy optimization schemes have been suggested ^[35], to our knowledge a scheme based on Eq (2.12) has not seen practical implementation in this context. The VPP methodologies presented here may prove applicable to a future application of Eq (2.12), since the expected value in Eq (2.12) can be approximated via Eq (2.11).

An example of a quantity-of-interest that we will use to illustrate our methods is the tumor burden, which is the total number of tumor cells as calculated from $n(\mathit{\boldsymbol{x}}, t)$:

If the goal is to predict $N(t)$ for some $t$ in the future, a VPP allows one to assess the uncertainty in this value as a result of uncertainty in the patient-specific parameter $\mathit{\boldsymbol{\beta}}_j$ being 'propagated' through the models $\mathit{\boldsymbol{\mathcal{F}}}$ and $\mathit{\boldsymbol{\mathcal{M}}}$ ^[25,30]. An illustration of uncertainty propagation is demonstrated in Figure 2, where several choices of the population parameter $\mathit{\boldsymbol{\theta}}_p$ are chosen and the resulting (estimated) probability distribution of $N(100)$ is shown. Further examples are given in section 5 below. Another possible QoI, the integrated log-kill, is similar to Eq (2.13) but with the integrand $n(\mathit{\boldsymbol{x}}, t)$ replaced by $\ln(n(\mathit{\boldsymbol{x}}, t)/n_0(\mathit{\boldsymbol{x}}))$. This metric is derived from a mass action drug effect and is applicable in chemotherapy evaluation ^[5]. Note that we do not discuss the important topic of exactly which scalar QoIs have real clinical impact in terms of patient outcomes; refer to e.g., ^[33,34].

3.   Molecular emission imaging data

In a pure simulation setting, a virtual biomarker sample of the form shown in Eq (2.11) can be computed without reference to any particular patient-specific data, so long as the corresponding probability distributions are well-calibrated. However, the model-based precision medicine context requires that patient-specific data be used to create individualized VPPs, hence making uncertainty-informed decisions about patient $j$ possible. We will focus on the case of molecular Emission Computed Tomography (ECT) data, for reasons discussed below.

3.1.   Preclinical and clinical imaging modalities

In the clinic, image data typically takes the form of reconstructed tomographic imaging studies such as MRI, CT, or $^{18}$F-FDG Positron Emission Tomography (PET). Note that while it can be difficult to precisely relate these standard tomographic data to either $\mathit{\boldsymbol{n}}_j$ or $\mathit{\boldsymbol{\beta}}_j$, it is nonetheless possible to postulate such relationships and thereby extract useful information from these images. An example of such an approach is the usage of Apparent Diffusion Coefficient (ADC) maps, generated using Diffusion-Weighted MRI (DWMRI), to estimate $\mathit{\boldsymbol{n}}_j$ through the application of a model which relates tissue cellularity to ADC ^[3,36,37]. The basic difficulty in such efforts is that malignancy is not directly measured by any of these standard techniques: For example, ADC maps estimate the apparent diffusion coefficient of water molecules and CT measures X-Ray attenuation, neither of which is uniquely altered by the presence of a malignant tumor. While increased $^{18}$F-FDG uptake implies increased cellular metabolism along the glycolytic pathway, a byproduct of malignant growth, this effect is not unique to malignant cells, since benign tumors can also display an increased metabolic rate ^[38].

In this article we consider a wide class of imaging techniques known as Emission Computed Tomography (ECT), broadly defined as any imaging modality (such as PET) which detects a molecular tracer distribution via optical or nuclear detection techniques. We choose to focus on ECT data because the wide variety of available tracers allows for more direct modeling relating these data to tumor cell density $\mathit{\boldsymbol{n}}_j$ and the RDE coefficient fields $\mathit{\boldsymbol{\beta}}_j$. As an example, genetically engineered cancer cells can be made to express Green Fluorescent Protein (GFP), which acts as a photon activity distribution when excited ^[39], and thus ECT, typically in the form of intravital microscopy (e.g., window chambers ^[40]) can be used to image $\mathit{\boldsymbol{n}}_j$ directly in vivo, at least in the small animal setting. Meanwhile, the growth factor $\mathit{\boldsymbol{\rho}}_j$ can potentially be measured using the PET tracer $^{18}$F-FLT, which is a radioactive version of the molecule thymadine used in DNA synthesis ^[41,42]. Ex vivo, the Ki-67 immunohistochemistry stain would typically be employed to measure proliferation, but this requires either a biopsy or sacrifice of the animal. Current progress in small animal imaging is working towards producing in vivo imaging agents for common immunohistochemistries, but to our knowledge an in vivo Ki-67 analog is not available. The other advantages of ECT are the possibility to perform polyscopic imaging (that is, measurement of several distinct physiological processes simultaneously) ^[43], and the wealth of literature on the statistics of raw and reconstructed ECT data ^[9].

To employ the statistical estimation strategies discussed in section 4, we require statistical models for ECT data. We emphasize the case of raw, un-reconstructed ECT data because the statistics are better understood: Most reconstructed images produced in the clinic are generated by a proprietary 'black box' method provided by the imaging system manufacturer. Our framework can still be applied to such data by carefully modifying the $\mathit{\boldsymbol{\mathcal{H}}}$ operator defined in Eq (3.3) to account for the reconstruction step, however it can be difficult to calibrate an appropriate noise model for likelihood-based estimation. If the MLEM algorithm is known to have been employed for reconstruction, see ^[44] for a discussion of the resulting noise characteristics.

3.2.   Statistics of ECT data

In an ECT system, a chemical tracer is engineered to emit some measurable particle, such as photons, charged particles or nuclear decay products ^[45]. A tracer corresponds to a particle activity distribution $f_j(\mathit{\boldsymbol{x}}, t)$, having units of particles emitted per unit volume per unit time. For the case of GFP imaging discussed above, a reasonable model is $f_j^{(GFP)}(\mathit{\boldsymbol{x}}, t) \propto n_j(\mathit{\boldsymbol{x}}, t)$, that is, the photon activity is proportional to the actual tumor cell density, with constant of proportionality related to the photon yield per cell (i.e., the number of photons detected, on average, per cell). Similarly, we can assume $\mathit{\boldsymbol{f}}_j^{(^{18}\text{F-FLT})}\propto \mathit{\boldsymbol{\rho}}_j$, the growth factor. The energy produced by $f_j(\mathit{\boldsymbol{x}}, t)$ propagates through the tissue medium and is detected by specialized ECT imaging hardware, which can range from optical CCD and CMOS detectors to scintillation-type detectors for high-energy particles. The statistics of raw ECT data are discussed at length elsewhere ^[9,46], so we summarize by noting that two data formats are commonly available, known as binned-mode data and particle-processed data. Because the presentation of particle-processing requires a more technical discussion of imaging detectors, we only consider binned-mode data in this work; see e.g., ^[45,46,47] for a discussion of particle processing. In a future work we will address the usage of particle processing data in mathematical oncology.

Binned-mode ECT data arises when detected particles are sorted into spatial bins corresponding to a discretization of the detector. These spatial bins need not correspond to physical detector pixels, which many ECT imaging detectors do not have. The physics-based statistics of particle counting nearly always lead to a Poisson data model for particle count data, which can be written in terms of a noisy linear functional of the activity $\mathit{\boldsymbol{f}}_j$ as follows. Suppose that the imaging system has $M$ detector bins, with data for patient $j$ labeled $g_{j1}, \ldots, g_{jM}$. The statistical model for binned-mode ECT data is that $\mathit{\boldsymbol{g}}_{j}$ is a Poisson random vector with mean

where $h_m(\mathit{\boldsymbol{x}}, t)$ is called the $m$-th sensitivity function of the imaging system (assumed to be fixed for all patients), $f_j(\mathit{\boldsymbol{x}}, t)$ is the particle activity distribution corresponding to the tracer administered to patient $j$, and $T$ is the total imaging time. For tomographic systems, $h_m(\mathit{\boldsymbol{x}}, t)$ must account for both the propagation within the patient, typically following a Radiative Transport Equation (RTE), as well as the geometry and blur characteristics of the imaging detector ^[6,9]. For this work, we will consider a simplified imaging geometry, intended to model intravital microscopy techniques such as mouse window chambers, where optical and/or gamma-ray images are collected with a planar detector above a pseudo-2D activity distribution. To a reasonable approximation, a Gaussian blur model

is then appropriate. This blur accounts for both imaging elements such as lenses and collimators as well as position estimation blur ^[9,46]. More sophisticated imaging system models, for instance 3D tomographic systems, will be discussed in a future work (see section 6).

Expressing Eq (3.1) as a linear continuous-to-discrete operator $\mathit{\boldsymbol{\mathcal{H}}}:\mathbb{X}\rightarrow \mathbb{R}^M$ allows us to write a binned-mode system model succinctly as

where $\mathit{\boldsymbol{\bar{g}}}_j = [\bar{g}_{j1}, \ldots, \bar{g}_{jM}]^\intercal = \mathit{\boldsymbol{\mathcal{H}f}}_j\in\mathbb{R}^M$ is defined via Eq (3.1). Considering $\mathit{\boldsymbol{f}}_j$ as a parameter, the probability distribution for the count vector $\mathit{\boldsymbol{g}}_j$ is thus

If $\mathit{\boldsymbol{f}}_j$ has been parameterized using a finite-dimensional parameter $\mathit{\boldsymbol{\theta}}_j$, for instance via a synthesis map of the form shown in Eq (2.8), we can define a nonlinear function $\mathit{\boldsymbol{H}}:\mathbb{R}^N\rightarrow\mathbb{R}^M$, where for image reconstruction typically $N < M$, as follows:

This function computes the mean image data $\mathit{\boldsymbol{\bar{g}}}_j$ for the synthesized object $\mathit{\boldsymbol{f}}_j = {\bf{\Phi }}(\mathit{\boldsymbol{\theta}}_j)$, and is implemented computationally via numerical quadrature methods applied to Eq (3.1). The component functions of $\mathit{\boldsymbol{H}}$ are denoted $H_m(\mathit{\boldsymbol{\theta}}_j):\mathbb{R}^N\rightarrow \mathbb{R}$. Note that if the synthesis method is nonlinear, as in Eq (2.8), then $\mathit{\boldsymbol{H}}$ is a nonlinear function of the vector $\mathit{\boldsymbol{\theta}}_j$, while if ${\bf{\Phi }}$ is linear, then $\mathit{\boldsymbol{H}}$ is linear. In either case, we can write a probability distribution for $\mathit{\boldsymbol{g}}_j\vert\mathit{\boldsymbol{\theta}}_j$ as follows:

We will make use of Eqs (3.4)–(3.6) for likelihood-based estimation methods in section 4 below.

4.   Statistical estimation and inversion techniques

Assume that for patient $j$, ECT data $\mathit{\boldsymbol{g}}_j\in \mathbb{R}^M$ has been collected, where $\mathit{\boldsymbol{g}}_j$ can either be monoscopic (a single imaging study) or polyscopic (consisting of multiple imaging studies, corresponding to multiple tracers). In the previous section, we provided statistical models for $\mathit{\boldsymbol{g}}_j$ that relate it to the activity distribution $\mathit{\boldsymbol{f}}_j$ via Eq (3.3). As discussed there, ECT tracers allow $\mathit{\boldsymbol{f}}_j$ to be directly related to $\mathit{\boldsymbol{n}}_j$ and $\mathit{\boldsymbol{\rho}}_j$. In this section, we discuss statistical techniques for estimating $\mathit{\boldsymbol{f}}_j$ (hence $\mathit{\boldsymbol{n}}_j$ and $\mathit{\boldsymbol{\rho}}_j$) from $\mathit{\boldsymbol{g}}_j$. In section 4.1, we discuss point estimation techniques for an individual patient, where a single estimate of $\mathit{\boldsymbol{f}}_j$ is produced, then in section 4.2 we discuss patient-specific Bayesian strategies where the solution consists of samples from a posterior distribution. In section 4.3, we discuss methods for estimating a population parameter $\mathit{\boldsymbol{\theta}}_p$ from a database $\mathit{\boldsymbol{G}} = [\mathit{\boldsymbol{g}}_1, \ldots, \mathit{\boldsymbol{g}}_J]$ of images collected for $J$ patients.

4.1.   Patient-specific maximum likelihood parameter estimation

Given data $\mathit{\boldsymbol{g}}_j \in\mathbb{R}^M$ for patient $j$, we consider the estimation of the underlying activity distribution $\mathit{\boldsymbol{f}}_j$ using the imaging system model $\mathit{\boldsymbol{g}}_j = \mathit{\boldsymbol{\mathcal{H}f}}_j + \mathit{\boldsymbol{\eta}}_j$, given in Eq (3.3). To approach this infinite-dimensional inverse problem, we select a finite-dimensional parameterization $\mathit{\boldsymbol{f}}_j = {\bf{\Phi }}(\mathit{\boldsymbol{\theta}}_j)$, where $\mathit{\boldsymbol{\theta}}_j\in \mathbb{R}^N$, and consider the problem of estimating $\mathit{\boldsymbol{\theta}}_j$ from $\mathit{\boldsymbol{g}}_j$, that is, solving the inverse problem $\mathit{\boldsymbol{g}}_j = \mathit{\boldsymbol{H}}(\mathit{\boldsymbol{\theta}}_j) + \mathit{\boldsymbol{\eta}}_j$. Note that moving from $\mathit{\boldsymbol{f}}_j$ to $\mathit{\boldsymbol{\theta}}_j$ might introduce an approximation error, owing to the fact that $\mathit{\boldsymbol{f}}_j$ need not be expressible in same form as the synthesis used for the reconstruction; in ECT imaging, $\mathit{\boldsymbol{f}}_j$ is typically quite smooth, so this approximation error is typically insignificant when compared to the influence of detector resolution and noise. While typical image reconstruction algorithms assume a voxel-based parameterization of $\mathit{\boldsymbol{f}}_j$, we elect to use a lumpy-type synthesis of the form given in Eq (2.8). As discussed, the activity distribution can correspond to the cell density $\mathit{\boldsymbol{n}}_j$, or the growth factor $\rho_j(\mathit{\boldsymbol{x}})$. We denote a generic estimate of $\mathit{\boldsymbol{\theta}}_j$ produced from $\mathit{\boldsymbol{g}}_j$ by a hat, i.e., $\mathit{\boldsymbol{\hat{\theta}}}_j$ is an estimate of $\mathit{\boldsymbol{\theta}}_j$ produced from the data $\mathit{\boldsymbol{g}}_j$. While a wide array of strategies exist to solve this inverse problem, we consider only likelihood-based methodologies.

Recall from section 3.2 that the probability distribution for binned-mode ECT data is Poisson. For an activity distribution $\mathit{\boldsymbol{f}}_j$ parameterized using the synthesis map ${\bf{\Phi }}$, i.e. $\mathit{\boldsymbol{f}}_j = \mathit{\boldsymbol{f}}_j(\mathit{\boldsymbol{\theta}}_j)$, we thus have a likelihood function

where $P(\mathit{\boldsymbol{g}}_j\vert\mathit{\boldsymbol{f}}_j)$ is the Poisson distribution given in Eq (3.4), and $H_m(\mathit{\boldsymbol{\theta}}_j)$ are the components of the function defined in Eq (3.5). Note the abuse of notation in Eq (4.1): The two functions $P(\mathit{\boldsymbol{g}}_j\vert\mathit{\boldsymbol{\theta}}_j)$ and $P(\mathit{\boldsymbol{g}}_j\vert\mathit{\boldsymbol{f}}_j)$ have different parameter types, but we take Eq (4.1) to define the former in terms of the latter. Taking its logarithm (and ignoring the $\mathit{\boldsymbol{\theta}}_j$-independent constant) results in the log-likelihood function

We now define the Maximum Likelihood Estimate (MLE) as the maximum of either Eq (4.1) or Eq (4.2) over the set ${\bf{\Theta }}$ of allowed parameters:

A variety of optimization algorithms exist to compute $\mathit{\boldsymbol{\hat{\theta}}}_j^{(ML)}$, including the Expectation Maximization (EM) algorithm and Newton-type methods ^[48,49]. If the imaging operator $\mathit{\boldsymbol{H}}(\mathit{\boldsymbol{\theta}}_j)$ is linear, the EM algorithm applied to Poisson data leads to an iterative method known in the imaging community as the MLEM algorithm. Starting from $\mathit{\boldsymbol{\hat{\theta}}}^{(0)}_j$, the MLEM algorithm performs the iterative update

We will employ the method (4.4) to compute some MLEs for linear imaging models in section 5, while for nonlinear models we use the Matlab algorithm $\tt{fmincon}$, which is a based on an interior point method with BFGS quasi-Newton steps ^[50]. A demonstration of algorithm (4.4) applied to a reconstruction of $n_j(\mathit{\boldsymbol{x}}, t)$ from image data of the form (3.1) is shown in Figure 3. Note that the reconstruction displays visual artifacts due to the function $\mathit{\boldsymbol{H}}(\mathit{\boldsymbol{\theta}})$ having a nontrivial null space. We have found that these artifacts have minimal impact on the task at hand, which is to predict the progression of $n(\mathit{\boldsymbol{x}}, t)$ and $N(t)$. Future work will analyze in more detail the question of task-based image quality, as defined in ^[9,51], in the mathematical oncology context.

We note briefly that under certain technical conditions, an MLE is asymptotically consistent and efficient, meaning for well-specified likelihood models, the MLE converges to the true parameter and attains the Cramér-Rao lower bound. Asymptotically, the distribution of the estimate $\mathit{\boldsymbol{\hat{\theta}}}_j^{(ML)}$ (as a random variable depending on the data $\mathit{\boldsymbol{g}}_j$) converges to a normal with covariance given by the inverse of the Fisher information ^[9]:

The Fisher information matrix $\mathit{\boldsymbol{F}} = \mathit{\boldsymbol{F}}(\mathit{\boldsymbol{\theta}}_j)$ can be used to design more efficient measurement systems by 'tuning' the design parameters of the system to modify the likelihood in a desirable manner ^[52]. Note also that the Fisher information matrix can be evaluated 'off-line', by computing the expected values in Eq (4.5), which allows one to analyze the uncertainty in the MLE in a frequentist setting. While $\mathit{\boldsymbol{F}}$ depends on the true $\mathit{\boldsymbol{\theta}}_j$, for many ECT systems this dependence is minimal (i.e., the Fisher information matrix is nearly constant with respect to $\mathit{\boldsymbol{\theta}}_j$), which allows one to compute a single (system-dependent) $\mathit{\boldsymbol{F}}$ to be used for all MLEs derived from the system. One can also estimate confidence intervals for scalar quantities of interest derived from estimates of $\mathit{\boldsymbol{\theta}}_j$ using transformation rules for the Fisher information matrix ^[53].

To use the MLE procedure to generate a patient-specific virtual population, we assume that $\mathit{\boldsymbol{g}}_j = \mathit{\boldsymbol{\mathcal{H}f}}_j$, where either $\mathit{\boldsymbol{f}}_j = n_j(\mathit{\boldsymbol{x}}, t_0)$, $\mathit{\boldsymbol{f}}_j = \mathit{\boldsymbol{\rho}}_j$, or $\mathit{\boldsymbol{f}}_j = [n_j(\mathit{\boldsymbol{x}}, t_0), \rho_j(\mathit{\boldsymbol{x}})]^T$, that is, we image either the cell density at a single time point, or the growth factor, or both the cell density and the growth factor. From this data, we take the following steps:

(i) Use MLE to obtain $\mathit{\boldsymbol{\hat{\theta}}}_j^{(\mathit{\boldsymbol{n}})}$ (and $\mathit{\boldsymbol{\hat{\theta}}}_j^{(\mathit{\boldsymbol{\rho}})}$, if measured) from which we synthesize $\mathit{\boldsymbol{\hat{n}}}_j$ (and $\mathit{\boldsymbol{\hat{\rho}}}_j$, if measured).

(ii) Randomize the remaining parameters $\mathit{\boldsymbol{D}}$, $\mathit{\boldsymbol{\kappa}}$ (and $\mathit{\boldsymbol{\rho}}$, if not measured) using a well-calibrated random field model such as Eq (2.8), resulting in $\mathit{\boldsymbol{D}}_j^{(1)}, \ldots, \mathit{\boldsymbol{D}}_j^{(J)}$, $\mathit{\boldsymbol{\kappa}}_j^{(1)}, \ldots, \mathit{\boldsymbol{\kappa}}_j^{(J)}$ (and $\mathit{\boldsymbol{\rho}}_j^{(1)}, \ldots, \mathit{\boldsymbol{\rho}}_j^{(J)}$, if required).

(iii) Generate a patient-specific virtual population by using $\mathit{\boldsymbol{\hat{n}}}_j$ as the initial condition, either $\mathit{\boldsymbol{\hat{\rho}}}_j$ or $\mathit{\boldsymbol{\rho}}_j^{(1)}, \ldots, \mathit{\boldsymbol{\rho}}_j^{(J)}$ for the growth factor, and $\mathit{\boldsymbol{D}}_j^{(1)}, \ldots, \mathit{\boldsymbol{D}}_j^{(J)}$, $\mathit{\boldsymbol{\kappa}}_j^{(1)}, \ldots, \mathit{\boldsymbol{\kappa}}_j^{(J)}$ for the diffusion coefficient and carrying capacity, and solving the forward problem (2.5) for each sample.

(iv) With the resulting VPP, compute any desired quantitative biomarkers as in Eq (2.11).

Note that in the above procedure, we make no attempt to estimate $\mathit{\boldsymbol{D}}$ and $\mathit{\boldsymbol{\kappa}}$ by performing the nonlinear parameter estimation task implied by Eq (2.4), nor do we attempt to estimate $n(\mathit{\boldsymbol{x}}, t)$ for any $t < t_0$. From single time point data such as this, an identifiability (i.e., uniqueness) issue arises preventing such estimation from succeeding. This procedure is used in sections 5.1 and 5.2 to produce virtual populations.

Note that one issue that is not addressed in the above technique is uncertainty inherent in the noisy data $\mathit{\boldsymbol{g}}_j$, namely that $\mathit{\boldsymbol{g}}_j$ being random means the MLE is also random. While Fisher information analysis can be employed to analyze this uncertainty in an off-line frequentist setting, this analysis does not provide a clear method for producing alternates to $\mathit{\boldsymbol{\hat{\theta}}}_j$, since one only has access to a single realization of $\mathit{\boldsymbol{g}}_j$. The Bayesian methodology discussed below provides an avenue towards producing many plausible estimates from a single realization of $\mathit{\boldsymbol{g}}_j$.

4.2.   Patient-specific Bayesian solution

In a patient-specific Bayesian solution to the inverse problem $\mathit{\boldsymbol{g}}_j = \mathit{\boldsymbol{\mathcal{H}f}}_j + \mathit{\boldsymbol{\eta}}_j$, we assume that $\mathit{\boldsymbol{f}}_j$ were sampled from a prior distribution $\mathbb{P}_0(\mathit{\boldsymbol{f}})$. Note that $\mathbb{P}_0(\mathit{\boldsymbol{f}})$ need not correspond to a well-calibrated model of the true population in order to proceed with the inference, but the method is typically much more reliable if it does. Then, since $\mathit{\boldsymbol{g}}_j$ is statistically coupled to $\mathit{\boldsymbol{f}}_j$ through Eq (3.3), we can consider the conditional random vector $\mathit{\boldsymbol{g}}_j\vert\mathit{\boldsymbol{f}}_j$, whose statistics are described by the likelihood function $\mathcal{L}(\mathit{\boldsymbol{f}}_j\vert \mathit{\boldsymbol{g}}_j) = \mathbb{P}(\mathit{\boldsymbol{g}}_j\vert\mathit{\boldsymbol{f}}_j)$ as discussed previously. A Bayesian inverse solution to the inference problem is to consider the posterior random vector $\mathit{\boldsymbol{f}}_j\vert \mathit{\boldsymbol{g}}_j$ instead of a single point estimate $\mathit{\boldsymbol{\hat{f}}}_j$. Samples from the posterior can be used to generate a virtual population as described in the procedure above, where now $\mathit{\boldsymbol{\hat{n}}}_j$ and $\mathit{\boldsymbol{\hat{\rho}}}_j$ are replaced by samples from their corresponding posteriors to account for uncertainty in the data $\mathit{\boldsymbol{g}}_j$.

Bayesian personalization of tumor growth has been considered in ^[21,54], for example. The former considers an RDE similar to Eq (2.1) but requires two segmented MRI images to perform the personalization. The latter generates ensemble members (i.e., virtual patients) by using fixed RDE coefficients ($\mathit{\boldsymbol{\beta}}$ in our notation) and imposing an ad hoc additive Gaussian noise model on the system evolution, in order to leverage the ensemble Kalman filter technique.

The distribution of the posterior $\mathit{\boldsymbol{f}}_j\vert\mathit{\boldsymbol{g}}_j$ is given via Bayes rule, which states informally that $\mathbb{P}_{post}(\mathit{\boldsymbol{f}}_j\vert\mathit{\boldsymbol{g}}_j) \propto \mathbb{P}_0(\mathit{\boldsymbol{f}})\mathcal{L}(\mathit{\boldsymbol{f}}_j\vert \mathit{\boldsymbol{g}}_j)$. More rigorously ^[55], we would say that

Practically speaking, in this paper we consider (via Eqs (2.7) and (2.8)) finite-dimensional parameterizations of $\mathit{\boldsymbol{f}}$, so that when coupled with binned-mode ECT data, the result is a finite-dimensional Bayesian inverse problem ^[56]. Taking as before a lumpy-type parameterization $\mathit{\boldsymbol{f}}_j = {\bf{\Phi }}(\mathit{\boldsymbol{\theta}}_j)$, and with the data described by the nonlinear function $\mathit{\boldsymbol{\bar{g}}}_j = \mathit{\boldsymbol{H}}(\mathit{\boldsymbol{\theta}}_j) = \mathit{\boldsymbol{\mathcal{H}}}{\bf{\Phi }}(\mathit{\boldsymbol{\theta}}_j)$ defined in section 3.2, we now assign a prior on the vector $\mathit{\boldsymbol{\theta}}_j$, which we assume has PDF $p_0(\mathit{\boldsymbol{\theta}}_j)$. The Poisson likelihood $\mathcal{L}(\mathit{\boldsymbol{\theta}}_j\vert\mathit{\boldsymbol{g}}_j)$ is given in Eq (3.6), leading to a description of the posterior as

If either the prior is infinite-dimensional or the ECT data is assumed to be of particle-processing form, a fully infinite-dimensional Bayesian approach must be considered ^[55].

To draw samples from the posterior (4.6), we apply a Markov Chain Monte Carlo (MCMC) algorithm based on the one presented in ^[57]. First, we assume for each tracer $\mathit{\boldsymbol{f}}_j$ an expansion of the form Eq (2.8), so that $\mathit{\boldsymbol{f}}_j = {\bf{\Phi }}(\mathit{\boldsymbol{\theta}}_j)$ with $\mathit{\boldsymbol{\theta}}_j = [b_1, \ldots, b_{L_{max}}, \mathit{\boldsymbol{x}}_1, \ldots, \mathit{\boldsymbol{x}}_{L_{max}}, \mathit{\boldsymbol{\gamma}}_1, \ldots, \mathit{\boldsymbol{\gamma}}_{L_{max}}]$. Note that in some cases, we reduce the computational burden by fixing the $\mathit{\boldsymbol{x}}_l$ and/or $\mathit{\boldsymbol{\gamma}}_l$ and only sample the remaining parameters. The MCMC procedure is then as follows:

(i) Given the current sample $\mathit{\boldsymbol{\theta}}_j^{(n)}$, we propose a new sample $\mathit{\boldsymbol{\theta}}_j^\prime \sim \pi\left(\mathit{\boldsymbol{\theta}}_j^\prime\vert \mathit{\boldsymbol{\theta}}_j^{(n)}\right)$ by randomly selecting $L_{mcmc}$ lumps and perturbing their amplitudes $b_l$, location $\mathit{\boldsymbol{x}}_l$ and shape $\mathit{\boldsymbol{\gamma}}_l$ by a Gaussian with diagonal covariance $\Sigma$. This proposal is symmetric.

(ii) From $\mathit{\boldsymbol{\theta}}_j^\prime$, the activity $\mathit{\boldsymbol{f}}_j^\prime$ is synthesized using Eq (2.8) and the proposed mean image $\mathit{\boldsymbol{\bar{g}}}^\prime = \mathit{\boldsymbol{\mathcal{H}f}}_j^\prime$ is simulated using Eq (3.5).

(iii) The acceptance probability is computed using the Metropolis-Hastings rule ^[25]:

after which we accept $\mathit{\boldsymbol{\theta}}_j^{(n+1)} = \mathit{\boldsymbol{\theta}}_j^\prime$ with probability $Q(\mathit{\boldsymbol{\theta}}_j^\prime, \mathit{\boldsymbol{\theta}}_j^{(n)})$, and take $\mathit{\boldsymbol{\theta}}_j^{(n+1)} = \mathit{\boldsymbol{\theta}}_j^{(n)}$ with probability $1-Q(\mathit{\boldsymbol{\theta}}_j^\prime, \mathit{\boldsymbol{\theta}}_j^{(n)})$.

By construction of the Metropolis-Hastings algorithm, the Markov chain samples $\mathit{\boldsymbol{\theta}}_j^{(1)}, \mathit{\boldsymbol{\theta}}_j^{(2)}, \ldots$ have stationary distribution given by the posterior $p(\mathit{\boldsymbol{\theta}}_j\vert \mathit{\boldsymbol{g}}_j)$, and can thus be used to synthesize virtual patients for the calculation of any quantitative biomarkers as discussed in section 2.4. An example of this approach is discussed in section 5.1. Note that the standard Metropolis-Hastings algorithm performs best for $\mathit{\boldsymbol{\theta}}$ with relatively small dimension; various modifications for high-dimensional problems have also been proposed ^[25,30,55].

4.3.   Calibration of population distributions

As illustrated in Figure 2, the population parameter $\mathit{\boldsymbol{\theta}}_p$ which specifies the random field models used to generate the virtual population has a substantial influence on the statistics of any biomarkers calculated using Eq (2.11). The question of how to estimate $\mathit{\boldsymbol{\theta}}_p$ from data thus arises. Such calibration problems in the context of VCTs have been discussed in ^[58,59], but their methodology accounts only for ODE models arising in systems pharmacology and does not directly apply to our context of PDE tumor growth models, random fields and imaging data. In this section, we discuss briefly a technique to estimate the population parameters which specify a lumpy-type model (2.8) from a database of ECT images $\mathit{\boldsymbol{G}} = [\mathit{\boldsymbol{g}}_1, \ldots, \mathit{\boldsymbol{g}}_J]$.

Suppose we assume a lumpy-type model for each of the RDE coefficient fields, $\mathit{\boldsymbol{\beta}} = [\mathit{\boldsymbol{D}}, \mathit{\boldsymbol{\rho}}, \mathit{\boldsymbol{\kappa}}]$, say each is of the form given in Eq (2.8) with $\ell(\mathit{\boldsymbol{x}})$ given by a Gaussian (2.9). While the lump function $\ell(\mathit{\boldsymbol{x}})$ can itself be the target of calibration (see ^[60], for instance), we assume Eq (2.9) for simplicity. The statistics of the random field $LB(\bar{L}, b_0, \sigma^2)$ are thus fully specified by the three scalars $\bar{L}, b_0$ and $\sigma^2$, which we write as $\mathit{\boldsymbol{\theta}}_p = [\bar{L}, b, \sigma^2]\in {\bf{\Theta }} = (0, \infty)^3$. The problem of calibrating a population distribution such as this is to take a database of noisy image data $\mathit{\boldsymbol{G}} = [\mathit{\boldsymbol{g}}_1, \ldots, \mathit{\boldsymbol{g}}_J]$ for a population of $J$ patients and produce an estimate $\mathit{\boldsymbol{\hat{\theta}}}_p = \left[\widehat{\bar{L}}, \hat{b}, \widehat{\sigma^2}\right]$. We can again apply either MLE or the Bayes approach as discussed in sections 4.1 and 4.2.

To develop a likelihood model for $\mathit{\boldsymbol{\theta}}_p$, we consider first the generic case where $\mathit{\boldsymbol{g}}_j$ consists of direct images of some activity $\mathit{\boldsymbol{f}}_j$, for instance $\mathit{\boldsymbol{f}}_j \propto \mathit{\boldsymbol{\rho}}_j$ as discussed in section 5.1. This case corresponds to that considered in ^[61]. Assuming inter-subject independence, we have

Using the law of total probability and the fact that the distribution of the image data is independent of the object statistics, we can write

In (binned-mode) ECT imaging, the probability distribution $P(\mathit{\boldsymbol{g}}_j\vert\mathit{\boldsymbol{f}}_j)$ is given by the Poisson (3.4), which depends on the imaging operator $\mathit{\boldsymbol{\mathcal{H}}}$. The expression (4.8) states that the probability of observing a given image $\mathit{\boldsymbol{g}}_j$ if the population parameter is $\mathit{\boldsymbol{\theta}}_p$ is the probability of observing $\mathit{\boldsymbol{g}}_j$ given the object $\mathit{\boldsymbol{f}}_j$, averaged over of all possible random field realizations $\mathit{\boldsymbol{f}}_j$ with statistics specified by the parameter $\mathit{\boldsymbol{\theta}}_p$. The expectation in Eq (4.8) is ostensibly infinite dimensional (since the $\mathit{\boldsymbol{f}}_j$ are functions), though it can be approximated using Monte Carlo:

Equation (4.9) specifies a computational method for approximating the likelihood (4.7): For fixed $\mathit{\boldsymbol{\theta}}_p$, compute Eq (4.9) for each image in the database, then compute the product (or its logarithm) in Eq (4.7). The MLE of $\mathit{\boldsymbol{\theta}}_p$ is then defined as before in Eq (4.3), with $\mathit{\boldsymbol{\theta}}_p\in{\bf{\Theta }}_p$, where ${\bf{\Theta }}_p$ is the set of admissible population parameters. Note that the estimate (4.9) can be computed in parallel using graphics processors, which drastically reduces the time required for each likelihood calculation, which would have previously been a bottleneck in the application of a traditional MLE. Despite this, the MLE optimization problem (4.3) applied to the likelihood (4.7) poses a challenge for standard numerical optimization routines due to scaling and non-convexity issues. See section 6 for a brief discussion of future work aimed at addressing these issues.

5.   Results

In this section, we discuss two in silico experiments that illustrate the methodologies discussed above. For both, a common 'ground truth' cell density path $n(\mathit{\boldsymbol{x}}, t)$ is generated by sampling a random field for each coefficient and simulating until $t = T = 365$ days. In the first experiment, we use (simulated) ECT data for $f_j(\mathit{\boldsymbol{x}}) = n(\mathit{\boldsymbol{x}}, 100)$ to compute both an MLE and a collection of posterior samples for the tumor cell density at $t = t_0 = 100$ days. Then, a virtual population is generated by using the resulting $\mathit{\boldsymbol{\hat{n}}}_j$ (or posterior samples) as the initial condition in Eq (2.1), randomizing the remaining coefficients and simulating until $t = T$. The resulting populations are then compared to the 'true' $n(\mathit{\boldsymbol{x}}, t)$, using the tumor burden $N(t)$ defined in Eq (2.13) as the quantity-of-interest. In the second experiment, we use ECT data and MLE for both the growth factor $\rho(\mathit{\boldsymbol{x}})$ and the cell density $n(\mathit{\boldsymbol{x}}, t_0)$, again generating a VPP by using the estimated cell density as the initial condition $n_0(\mathit{\boldsymbol{x}})$ and $\mathit{\boldsymbol{\rho}}_j = \mathit{\boldsymbol{\hat{\rho}}}_j$, then randomizing both $\mathit{\boldsymbol{D}}$ and $\mathit{\boldsymbol{\kappa}}$.

For both experiments, we use the following random fields, both to generate the initial 'true' tumor path $n(\mathit{\boldsymbol{x}}, t)$ and to randomize unmeasured coefficients. We assume $\mathit{\boldsymbol{D}}\sim LB(20, 1\text{e$-7$}, 0.04)$, $\mathit{\boldsymbol{\rho}}\sim LB(200, 0.25, 0.002)$, and $\mathit{\boldsymbol{\kappa}}\sim LB(100, 5\text{e7}, 0.1)$. These parameters were selected to achieve an average tumor size of approximately $10^8$ cells after the first year; obviously a calibration procedure such as the one outlined in section 4.3 would be required to estimate population parameters for a real population of tumors. The exact sample used is shown in Figure 4.

5.1.   Experiment 1: Maximum likelihood estimation and posterior sampling of cell density $\mathit{\boldsymbol{n}}$

In this experiment, we suppose that $\mathit{\boldsymbol{g}}_j^{(n)} = \mathit{\boldsymbol{\mathcal{H}}}^{(n)}\mathit{\boldsymbol{n}}_j + \mathit{\boldsymbol{\eta}}_j^{(n)}$ consists of ECT measurements of the cell density $n_j(\mathit{\boldsymbol{x}}, t_0)$, with $t_0 = 100$ days, and that no other measurements are available. We assume that the experimental setup is confined to a thin slab so that a 2D simulation described in section 3.1 is appropriate, and that the imaging system corresponds to Eq (3.1) with blur kernel given by Eq (3.2). The blur is taken to be $\sigma^2_{blur} = 0.0212$, which corresponds to an imaging resolution (Full-Width-Half-Max, FWHM) of 500 $\mu$m, and we assume that the photon yield is $A = 1$e$-3$ (detected photons per cell). Once the mean image data $\mathit{\boldsymbol{\bar{g}}}_j = \mathit{\boldsymbol{\mathcal{H}f}}_j$ is computed using Eq (3.1), a randomized image is produced using $\tt{poissrnd}$.

To perform the image reconstruction for $n_j(\mathit{\boldsymbol{x}}, t_0)$, we apply the linear MLEM algorithm (4.4) with a reconstruction of the type (2.8), with fixed lump centers and lump variance, so that $\mathit{\boldsymbol{\theta}} = [b_1, \ldots, b_{L_{max}}]$. Specifically, we fix the lump centers on a grid of size $25\times 25$, uniformly spaced on $[0.35, 0.65]^2$, and set the lump variance to $\sigma^2 = 1$e$-2$$/256$. The matrix $\mathit{\boldsymbol{H}}$ with columns given by $\mathit{\boldsymbol{\bar{g}}}_l = \mathit{\boldsymbol{\mathcal{H}}}(\ell(\mathit{\boldsymbol{x}}-\mathit{\boldsymbol{x}}_{l}))$ is then computed using Eq (3.1), and the MLEM algorithm (4.4) applied for $K_{max} = 5000$ iterations, with initial guess $\mathit{\boldsymbol{\hat{\theta}}}_j^{(0)} = \mathit{\boldsymbol{H}}^T \mathit{\boldsymbol{g}}_j$. From the final $\mathit{\boldsymbol{\hat{\theta}}}_j^{(K_{max})}$, we synthesize the estimate $\hat{n}_j(\mathit{\boldsymbol{x}}, t_0)$ using Eq (2.8). The results of this reconstruction are shown in Figure 5.

To demonstrate the usage of this estimate to generate a virtual tumor population as discussed in section 4.1, we use the resulting $\hat{n}_j(\mathit{\boldsymbol{x}}, t_0)$ as the initial condition $n_0(\mathit{\boldsymbol{x}})$ in Eq (2.1), and randomize the remaining coefficient fields using the lumpy-type random field models described in the introduction to this section to produce a virtual population of size $J = 64$. In Figure 6, the quantity-of-interest $N(t)$ is shown for the resulting virtual population, while in Figure 7, a collection of 10 virtual tumors is shown for $t = 150$ (i.e., 50 days post-imaging).

To demonstrate the usage of the Bayesian methods discussed in section 4.2, we generated samples from the posterior $\mathit{\boldsymbol{n}}_j\vert\mathit{\boldsymbol{g}}_j$ to use as part of the virtual population generation procedure. For the MCMC procedure, we use the algorithm outlined in section 4.2, moving $L_{mcmc} = 3$ lumps for each proposal with a proposal distribution of $\mathcal{N}(0, \sigma^2_{mcmc}\mathit{\boldsymbol{I}})$, taking $\sigma_{mcmc} = 1$e4. The prior for $\mathit{\boldsymbol{\theta}}_j = [b_1, \ldots, B_{L_{max}}]$ is chosen to be I.I.D. uniform on the interval $[0, 1\text{e}10]$. Starting from $\mathit{\boldsymbol{\theta}}_j^{(0)} = \mathit{\boldsymbol{\hat{\theta}}}_j^{(ML)}$, i.e., starting at the maximum likelihood estimate, we generate $J = 512$ samples, then synthesize $\mathit{\boldsymbol{n}}_j^{(1)}, \ldots, \mathit{\boldsymbol{n}}_j^{(J)}$. Note that the subscript remains $j$ as these samples correspond to patient $j$, while the superscript indicates that these are posterior samples, generated using MCMC. Starting at the maximum likelihood estimate (which for our prior corresponds to the MAP estimate) reduces the need for an expensive burn-in period in which the first $J_0$ samples are thrown out to ensure that the chain is sampling from the posterior. The acceptance rate for our chain (number of proposals accepted divided by total number of samples) is 0.495. Using each posterior sample as an initial condition, we again generate a virtual population by solving Eq (2.1) with the remaining coefficients $\mathit{\boldsymbol{D}}, \mathit{\boldsymbol{\rho}}, \mathit{\boldsymbol{\kappa}}$ randomized according to the lumpy-type field models specified previously. The results of this virtual population are shown in Figure 8.

For both the maximum likelihood and Bayesian solutions, the resulting virtual populations under-estimate, on average, the true tumor burden both over the short- and long-term, which indicates that additional data may be needed to estimate other model coefficients to gain a more accurate patient-specific prediction; we discuss on possibility in section 5.2 below. While the Bayesian technique accounts for uncertainty inherent in the image data $\mathit{\boldsymbol{g}}_j^{(n)}$, this uncertainty appears relatively small as it relates to the prediction of $N(t)$; the variance of $N(t)$ for the virtual population generated using the Bayesian technique is only slightly larger than the variance for the MLE population.

5.2.   Experiment 2: Maximum likelihood estimation of both $\mathit{\boldsymbol{n}}$ and $\mathit{\boldsymbol{\rho}}$

In Experiment 1, we saw that the virtual population under-estimated the true tumor path in both the MLE and Bayesian examples, suggesting that additional data may be required. In this experiment, we assume that the image data $\mathit{\boldsymbol{g}}_j$ consists of (simulated) polyscopic data of the following form. We assume that we have access first to cell density measurements of the form $\mathit{\boldsymbol{g}}_j^{(n)} = \mathit{\boldsymbol{\mathcal{H}}}^{(n)}\mathit{\boldsymbol{n}}_j + \mathit{\boldsymbol{\eta}}_j^{(n)}$, where $\mathit{\boldsymbol{\mathcal{H}}}^{(n)}$ is the same as in Experiment 1. The second dataset is of the form $\mathit{\boldsymbol{g}}_j^{(\rho)} = \mathit{\boldsymbol{\mathcal{H}}}^{(\rho)}\mathit{\boldsymbol{\rho}}_j + \mathit{\boldsymbol{\eta}}_j^{(\rho)}$, where we now have $\sigma^2_{blur} = 0.0425$ (1 mm FWHM). As in the first experiment, we assume that imaging is performed at $t = t_0 = 100$. The image reconstruction for $\mathit{\boldsymbol{n}}_j$ is performed using MLEM as in experiment 1 (Figure 5). For $\rho_j(\mathit{\boldsymbol{x}})$, we use MLE with a nonlinear synthesis map as follows. The lumpy-type expansion (2.8) with Gaussian lump (2.9) is assumed where now $\mathit{\boldsymbol{\theta}}_j = [\mathit{\boldsymbol{x}}_1, \ldots, \mathit{\boldsymbol{x}}_{L_{max}}, b, \sigma^2]$, that is, we fit the lump positions and a common lump amplitude and variance. The result is a reconstruction of the form

where we have chosen $L_{max} = 60$ for this experiment. We assume further the constraints that $\mathit{\boldsymbol{x}}_l\in [0, 1]^2$, $b\in (0, \infty)$, and $\sigma^2\in (0, 0.1]$, and use the Matlab algorithm $\tt{fmincon}$ to perform the MLE optimization (4.3). The results of this reconstruction are shown in Figure 9.

In Figure 10, we compare the tumor burden predicted using the virtual population to the 'true' tumor path. Over the first 30 day post-imaging period, the ensemble mean $\bar{N}(t)$ for the virtual population tracks the true tumor path closely, indicating that if both $\mathit{\boldsymbol{n}}$ and $\mathit{\boldsymbol{\rho}}$ are estimated, the result appears unbiased over this short time period. This is also demonstrated by looking at the individual virtual population samples in Figure 11: Visual inspection reveals a closer overall appearance to the 'true' tumor than was apparent in Figure 7. Over a longer time period, the ensemble still trends towards a lower overall growth rate than the 'true' tumor path, as seen previously in Figure 6.

6.   Discussions

In this article, we have presented a rigorous, end-to-end statistical framework for addressing both inter- and intra-patient heterogeneity for the commonly used reaction-diffusion tumor growth model (2.1), emphasizing the usage of raw molecular emission imaging data to estimate unknown model parameters. We demonstrated these methods in two in silico experiments, using a simulated ground truth tumor path and simulated imaging data to illustrate patient-specific Virtual Patient Population (VPP) generation.

In section 2.3, we presented a family of random field models that are both easy to simulate and demonstrate a wide variety of spatial heterogeneity; in section 4.3, we discussed how the parameters defining the statistics of these random field models can be estimated from a collection of molecular images, and discussed some of the computational challenges inherent in this task. Future work will discuss algorithms for solving this MLE problem, as well as additional estimation strategies to calibrate more general random field models from noisy data.

In section 3, we discussed molecular Emission Computed Tomography (ECT) imaging systems and their associated statistical models. For simplicity, we emphasized the case of performing direct, possibly polyscopic, 2D imaging of a thin sample. While the 2D approximation is very limiting in general, it is applicable to the experimental setting of mouse window chambers, a topic of future work in our group. In another future work, we will address fully 3D tomographic ECT systems with particle-processing data and the associated computational and image quality questions that arise. Understanding image data quality in a task performance setting can help design better imaging systems ^[51], and to our knowledge, there is currently minimal effort formally evaluating the quality of image data in the mathematical oncology context.

In section 4, we presented general likelihood-based statistical procedures for estimating both patient-specific and population parameters from noisy molecular imaging data, and in section 5, we presented two experiments to illustrate these procedures. In the first experiment, measurements of $n_j(\mathit{\boldsymbol{x}}, 100)$ are used to construct a maximum likelihood estimate and posterior samples, which are subsequently used to generate a patient-specific VPP. In the second experiment, we add an additional imaging study, using only maximum likelihood to estimate both $n_j(\mathit{\boldsymbol{x}}, 100)$ and $\rho_j(\mathit{\boldsymbol{x}})$. The resulting virtual population demonstrated a more accurate prediction of the true $N_j(t)$, over a 30 day post-imaging prediction period. In a future work, we will investigate the question of optimizing imaging acquisition schemes for the task of predicting $N(t)$. We will also address the addition of treatment models and evaluation techniques discussed in ^[5] to our virtual population methodology, working towards providing a general strategy to solving the optimization problem (2.12).

Acknowledgments

The author was supported by the NIBIB grant R01EB000803 and travel funding from the Statistical and Applied Mathematical Sciences Institute. The author also acknowledges very useful comments from reviewers, discussions with the SAMSI working group on tumor heterogeneity and with Harrison H. Barrett and Luca Caucci from the University of Arizona.

Conflict of interest

The authors declare no conflict of interest for this paper.