Bayesian calibration, validation, and uncertainty quantification of diffuse interface models of tumor growth

Abstract

The idea that one can possibly develop computational models that predict the emergence, growth, or decline of tumors in living tissue is enormously intriguing as such predictions could revolutionize medicine and bring a new paradigm into the treatment and prevention of a class of the deadliest maladies affecting humankind. But at the heart of this subject is the notion of predictability itself, the ambiguity involved in selecting and implementing effective models, and the acquisition of relevant data, all factors that contribute to the difficulty of predicting such complex events as tumor growth with quantifiable uncertainty. In this work, we attempt to lay out a framework, based on Bayesian probability, for systematically addressing the questions of Validation, the process of investigating the accuracy with which a mathematical model is able to reproduce particular physical events, and Uncertainty quantification, developing measures of the degree of confidence with which a computer model predicts particular quantities of interest. For illustrative purposes, we exercise the process using virtual data for models of tumor growth based on diffuse-interface theories of mixtures utilizing virtual data.

Appendix A: Discretization of the models

In this section, we consider the discretization and numerical solution of the three models described in Sect. 4.2. We consider a mixed finite element approximation in space combined with a semi-implicit time-stepping algorithm. We focus on model $\mathcal M _3$ given in (26). The discretization of $\mathcal M _1$ can be derived from $\mathcal M _3$ as a special case. Of course, $\mathcal M _2$ is equivalent to $\mathcal M _3$ up to time-independent parameters.

1.1 A.1 Weak formulation

Anticipating possible low regularity in solutions and in preparation for general mixed finite element approximations, we first consider a weak-mixed formulation of (26). Let $\mathsf{V},\,\mathsf{V}_g$ and $\mathsf{W}_{u_0}$ denote the spaces

$$\begin{aligned} \mathsf{V}&:= L^2\left(0,T; H^1\left(\varOmega \right)\right)\!, \quad \mathsf{V}_g := \Big \{ v\in \mathsf{V} : v = g \text{ on} \, \partial \varOmega ,\, \text{ a.e.} t \Big \},\\ \mathsf{W}_{u_0}&:= \Big \{ v\in \mathsf{V} \cap L^{\infty }(0,T;L^\infty (\varOmega )),\, v_t \in \mathsf{V}^{\prime }=L^2\big (0,T;(H^1(\varOmega ))^{\prime }\big ),\, v(0) = u_0 \Big \}. \end{aligned}$$

In particular, note that $\mathsf{V}_0 = L^2\left(0,T; H^1_0(\varOmega )\right)$ and $\mathsf{V}_1 = \{1\} + \mathsf{V}_0$. We shall now define a weak solution of (25) to be a triple $\mathbf{u}:=(u,\mu ,c) \in \mathsf{W}_{u_0} \times \mathsf{V} \times \mathsf{V}_1$ such that

$$\begin{aligned} \begin{aligned} \int \limits _0^T \Big (\left\langle u_t, \varphi \right\rangle + B(\mathbf{u};\mathbf{w}) \Big ) \, \mathrm{d} t = 0 \quad \quad \forall \mathbf{w} := (\varphi ,\eta ,\xi )\in \mathsf{V} \times \mathsf{V} \times \mathsf{V}_0, \end{aligned} \end{aligned}$$

(40)

where $B\left(\cdot ; \cdot \right)$ is the semi-linear form,

$$\begin{aligned} B\left(\mathbf{u};\mathbf{w}\right)&= \left( M u^2 \nabla \mu ,\nabla \varphi \right) - \left(P c u - Au,\varphi \right)\nonumber \\&\quad +\left(\mu - f^{\prime }\left(u \right), \eta \right) + \epsilon \left( \chi \,c, \eta \right) - \epsilon ^2 \left(\nabla u,\nabla \eta \right)\nonumber \\&\quad + \left(D\, \nabla c,\nabla \xi \right) + \left( c u,\xi \right), \end{aligned}$$

(41)

$\left\langle \cdot ,\cdot \right\rangle$ denotes the duality pairing between $\left(H^1 \left(\varOmega \right)\right)^{\prime }$ and $H^1\left(\varOmega \right),\,\left(\cdot ,\cdot \right)$ denotes the $L^2$ inner product over $\varOmega$, and $\left(\cdot ,\cdot \right)_{\partial \varOmega }$ is the $L^2$ inner product over $\partial \varOmega$.

The well-posedness of the above weak formulation (40) is an open question not considered here. The major difficulties are the nonlinear mobility $M u^2$ and the nonlinear coupling between $c$ and $u$, i.e. the reaction terms $cu$. On the other hand, the simpler model $\mathcal M _1$ is a compact perturbation of the Cahn–Hilliard equation, the well-posedness of which can be proven; see for instance Elliott (1989) and Elliott and Garcke (1996). The semilinear form (41) is bounded above for $u$ bounded [in $\mathsf{W}_{u_0} \subset L^\infty (0,T;L^\infty (\varOmega ))$], but this may not be true for certain cases (see Caffarelli and Muler (1995) for a related result for the Cahn–Hilliard equation).

1.2 A.2 Semi-discrete finite element approximations

We now address the construction of semi-discrete conforming finite element approximations of the system (40). Such a formulation involves a discretization of the spatial variations of $\left(u,\mu ,c\right)$ keeping, for the moment, the time variations continuous in time.

Let $\mathcal T ^h$ denote a member of a family of partitions of domain $\varOmega$ into meshes of non-overlapping convex finite elements $\varOmega _k$ such that $\overline{\varOmega } = \bigcup _{k=1}^{N\left(h\right)} \overline{\varOmega }_k$ and $\varOmega _k\cap \varOmega _j = \emptyset ,\, k\ne j$. Each $\varOmega _k$ is the image of a master element $\tilde{\varOmega }$ under an invertible (generally affine) map $F_k$. We construct finite-dimensional subspaces of $H^1\left(\varOmega \right)$ and $H^1_0(\varOmega )$ defined by

$$\begin{aligned} \mathcal S ^h := \left\{ v^h \in H^1\left(\varOmega \right) : v^h|_{\varOmega _k} = \hat{v}\circ F_k^{-1}, 1\le k\le N\left(h\right), \hat{v} \in \mathbb P ^K\big (\overline{\tilde{\varOmega }}\big )\right\} , \end{aligned}$$

(42)

and $\mathcal S ^h_0 := \mathcal S ^h \cap H_0^1(\varOmega )$, respectively. We further set $\mathcal S ^h_1 := \{1\} + \mathcal S ^h_0$. Here $\mathbb P ^K(\tilde{\varOmega })$ is either the space of polynomials of degree $\le K$ defined on the closure of $\tilde{\varOmega }$ or the space of tensor products of polynomials of degree $K$ on the closure of $\tilde{\varOmega }$. With these notations and conventions in place, the semi-discrete approximation $\mathbf{u}^h := (u^h,\mu ^h,c^h):[0,T]\rightarrow \mathcal S ^h\times \mathcal S ^h \times \mathcal S ^h_1$ is defined as follows:

$$\begin{aligned} \big ( u^h_t(t), \varphi ^h \big ) \!+\! B(\mathbf{u}^h(t);\mathbf{w}^h) \!=\! 0 \quad \quad \forall \mathbf{w}^h \!:=\! (\varphi ^h,\eta ^h,\xi ^h)\in \mathcal S ^h\times \mathcal S ^h \times \mathcal S ^h_0,\qquad \end{aligned}$$

(43)

for a.e. $t\in (0,T)$. The initial condition for $u^h$ is set by $u^h(0) = \mathcal I ^h \, u_0$, where $\mathcal I ^h: H^1(\varOmega )\rightarrow \mathcal S ^h$ is a suitable interpolation or projection.

Let us denote by $\left\{ \varphi _i\left(\mathbf{x}\right)\right\} _{i=1}^N$ a set of basis functions of $\mathcal S ^h$ generated on a partition $\mathcal T ^h$ by the finite element approximations: $\mathcal S ^h = \text{ span}\left\{ \varphi _i\right\}$. Then each of the component approximations is of the form

$$\begin{aligned} u^h\left(t,\mathbf{x}\right)&= \sum _{i=1}^N u_i\left(t\right)\varphi _i\left(\mathbf{x}\right), \quad \mu ^h\left(t,\mathbf{x}\right) = \sum _{i=1}^N \mu _i\left(t\right) \varphi _i\left(\mathbf{x}\right),\nonumber \\&\quad c^h\left(t,\mathbf{x}\right) = \sum _{i=1}^N c_i\left(t\right)\varphi _i\left(\mathbf{x}\right)\!. \end{aligned}$$

(44)

Here $u_i\left(t\right),\,\mu _i\left(t\right)$ and $c_i\left(t\right)$ are continuous functions in time. Upon introducing these into (43), we obtain a system of nonlinear ordinary differential equations (ODE’s) in the unknown discrete variables, $u_i\left(t\right),\,\mu _i\left(t\right)$ and $c_i\left(t\right)$. We proceed by developing a time discretization of these ODE’s.

1.3 A.3 Time discretization

An attractive class of robust schemes for time-integration of Cahn–Hilliard-type problems are semi-implicit schemes Eyre (1998); Wise et al. (2009); Wise (2010). Semi-implicit schemes are preferred for their superior stability properties allowing large time-steps to be taken. These schemes split the free energy $f$ in a contractive and expansive part, $f(u) = f_c(u) - f_e(u)$, and accordingly treat $f_c$ implicitly and $f_e$ explicitly. The splitting is such that $f_c$ and $f_e$ are both convex functions, at least in the domain of interest (for $u\in [0,1]$). We shall employ the splitting $f_c(u) = \frac{3}{2}\, u^2$ and $f_e(u):=f_c(u)-f(u)$.

Let $\varDelta t$ denote the time step and $t_n = n \varDelta t$, for $n=0,\ldots , N:=T/\varDelta t$ corresponding discrete times. We denote the approximation to $\mathbf{u}^h(t_n)$ by $\mathbf{u}_n$. Then, the fully discrete approximation $\mathbf{u}_n:=(u_n,\mu _n,c_n)\in \mathcal S ^h\times \mathcal S ^h \times \mathcal S ^h_1$ is defined recursively by

$$\begin{aligned} \left(\frac{u_{n} - u_{n-1}}{\varDelta t}, \varphi ^h \right) + \left( M u_{n-1}^2 \nabla \mu _n,\nabla \varphi ^h\right) - \left(P c_{n-1} u_{n} -A u_{n},\varphi \right)&= 0 \quad&\forall \varphi ^h\in S^h,\end{aligned}$$

(45a)

$$\begin{aligned} \left(\mu _n - f_c^{\prime }\left(u_n \right) + f_e^{\prime }(u_{n-1}), \eta ^h \right) + \epsilon \left( \chi (t_n)\, c_n, \eta ^h \right) - \epsilon ^2 \left(\nabla u_n,\nabla \eta ^h\right)&= 0 \quad&\forall \eta ^h\in S^h,\end{aligned}$$

(45b)

$$\begin{aligned} \left(D(t_n)\,\nabla c_n,\nabla \xi ^h\right) + \left( c_n u_{n-1} ,\xi ^h\right)&= 0 \quad&\forall \xi ^h\in S^h_0, \end{aligned}$$

(45c)

for $n=1,\ldots , N$. The scheme is initialized using the initial condition $u_0 = u^h(0) = \mathcal I ^h\,u_0$.

Note that analogous to $f$, we have treated the other nonlinear terms $(M u^2 \nabla \mu ,\,cu)$ in a semi-implicit fashion. Since our splitting of $f(u)$ resulted in a linear contractive part $f^{\prime }_c(u)$, the above scheme requires the solution of a linear system of equations at each time step.