Simulation of Multispecies Desmoplastic Cancer Growth via a Fully Adaptive Non-linear Full Multigrid Algorithm

A fully adaptive non-linear full multigrid (FMG) algorithm is implemented to computationally simulate a model of multispecies desmoplastic tumor growth in three spatial dimensions. The algorithm solves a thermodynamic mixture model employing a diffuse interface approach with Cahn-Hilliard-type fourth-order equations that are coupled, non-linear, and numerically stiff. The tumor model includes extracellular matrix (ECM) as a major component with elastic energy contribution in its chemical potential term. Blood and lymphatic vasculatures are simulated via continuum representations. The model employs advection-reaction-diffusion partial differential equations (PDEs) for the cell, ECM, and vascular components, and reaction-diffusion PDEs for the elements diffusing from the vessels. This study provides the details of the numerical solution obtained by applying the fully adaptive non-linear FMG algorithm with finite difference method to solve this complex system of PDEs. The results indicate that this type of computational model can simulate the extracellular matrix-rich desmoplastic tumor microenvironment typical of fibrotic tumors, such as pancreatic adenocarcinoma.

Introduction

The process of cancer progression is driven by the communication between tumor cells and their surroundings. A dynamic tumor microenvironment typically consists of highly proliferating neoplastic cells of different phenotypes, necrotic tumor cells, infiltrating innate and adaptive immune inflammatory cells, cancer-associated fibroblasts, cancer stem cells, extracellular matrix (ECM), blood and lymphatic vessels, pericytes, healthy host cells, and a variety of soluble molecules (Hanahan and Weinberg, 2000; de Visser and Coussens, 2006; Tlsty and Coussens, 2006; Whiteside, 2008; Perez-Moreno, 2009). These cellular and molecular elements dictate the tumor progress from unregulated neoplastic growth to potential metastasis. In its heterogeneous milieu with complex tumor-induced interactions and mechanical stress, an evolving tumor mass also undergoes transient morphological changes arising from cell motility and cell-cell/cell-ECM interactions. Mathematical modeling of cancer progression including its associated microenvironment may be a useful tool for predicting tumor dynamics and cancer response to therapy.

To study the desmoplastic tumor microenvironment, we recently presented a tumor model (Ng and Frieboes, 2017) as a continuum scale multicomponent-multispecies system consisting of heterogeneous cell types and ECM. This thermodynamic mixture model, inspired by the one derived in Wise et al. (2008) and Frieboes et al. (2010), includes metabolic reactions, tumorigenic factors, desmoplastic response, as well as tumor-induced angiogenesis and lymphangiogenesis. The diffuse interface approach is implemented as derived in Wise et al. (2008), where thermodynamically consistent Darcy velocities and Fickian diffusive terms are determined from the energy variation. In the Helmholtz free energy equation, the square gradient model (Cahn and Hilliard, 1958; Cahn, 1959; Yang et al., 1976; Rowlinson, 1979) is used to describe interfaces arising from the adhesive properties of cells and ECM components, and an elastic energy term by Leo et al. (1998) is added to represent the elastic properties of the ECM component. Continuous blood and lymphatic vessel densities are modified from cell fluxes employed in Anderson and Chaplain (1998), Chaplain (1996), and Mantzaris et al. (2004). Sprout initiation conditions of vessels, as well as interactions between angiogenic factors, proteolytic enzymes, and the ECM component are inspired by Levine et al. (2000), Levine et al. (2001a,b). Lymphangiogenesis is assumed to behave similarly to angiogenesis. Nutrients and waste products from cell metabolism are governed by fluxes and reaction rates modified from Casciari et al. (1992).

Continuum models represent cell populations and molecular species that influence the cell cycle events as continuous variables (see recent reviews Roose et al., 2007; Preziosi and Tosin, 2009a; Tracqui, 2009; Byrne, 2010; Cristini and Lowengrub, 2010; Edelman et al., 2010; Kreeger and Lauffenburger, 2010; Lowengrub et al., 2010; Osborne et al., 2010; Rejniak and McCawley, 2010; Vineis et al., 2010; Andasari et al., 2011; Chaplain, 2011; Deisboeck et al., 2011; Frieboes et al., 2011; Michor et al., 2011; Rejniak and Anderson, 2011; Bachmann et al., 2012; Oden et al., 2015 and references therein), These models typically implement ODE or PDE approaches to describe an advection-diffusion-reaction system. Continuum multiphase/mixture mechanochemical models include chemical and mechanical interactions between phases (cell types or species) (see Araujo and McElwain, 2004; Hatzikirou et al., 2005; Quaranta et al., 2005; Byrne et al., 2006; Graziano and Preziosi, 2007; Roose et al., 2007; Astanin and Preziosi, 2008; Preziosi and Tosin, 2009a; Tracqui, 2009; Lowengrub et al., 2010) and associated references). These models typically introduce a stress tensor, velocity, and pressure for each phase by enforcing mass, momentum, and energy balances (Ambrosi and Preziosi, 2002; Breward et al., 2002, 2003; Byrne and Preziosi, 2003; Byrne et al., 2003; Araujo and McElwain, 2005a,b; Graziano and Preziosi, 2007; Astanin and Preziosi, 2008; Galle et al., 2009; Preziosi and Tosin, 2009b; Bresch et al., 2010; Preziosi et al., 2010; Preziosi and Vitale, 2011; Sciumé et al., 2013; Klika, 2014). Related to these mixture models is the diffuse interface approach (Oden et al., 2010; Hawkins-Daarud et al., 2012; Chen et al., 2014), for which the square gradient theory can be used to describe smooth transitions within thin interfacial regions. The gradient contributes to the Helmholtz free energy, from which the component velocities, pressures, and diffusive terms can be derived (Wise et al., 2008; Chen and Lowengrub, 2014). Continuum single- or multi-phase models that consider the effects of cell-cell and/or cell-ECM adhesion have included (Frieboes et al., 2007, 2013; Ambrosi and Preziosi, 2009; Bearer et al., 2009; Kuusela and Alt, 2009; Chatelain Clément et al., 2011; Escher and Matioc, 2013). In Gerisch and Chaplain (2008), Preziosi and Tosin (2009b), Psiuk-Maksymowicz (2013), Wu et al. (2013), Sciumé et al. (2014a,b), and Arduino and Preziosi (2015), the ECM is represented as a key component of the tumoral tissue.

This paper presents the fully adaptive non-linear FMG algorithm and numerical solution of the model presented in Ng and Frieboes (2017). This diffuse interface model is characterized by non-linear fourth-order Cahn-Hilliard type PDEs and narrow transition layers. A fully adaptive block-structure Cartesian mesh is used on adaptively refined grid levels to address the need of transient locally-refined mesh regions for the narrow transition layers. To avoid severe time-step restrictions of explicit methods and fine resolution of transition layers, a fully adaptive non-linear finite difference multigrid method inspired by the pioneering work of Wise et al. (2007) and Wise et al. (2011) is implemented to solve the system of equations. The equations are discretized in time by the semi-implicit Crank-Nicolson method. The error smoothing steps employ the non-linear Full Approximation Scheme, and iterations are carried out in V-cycles. Modified from the numerical solution and method given in Wise et al. (2007, 2011) we implement a fully adaptive non-linear full multigrid (FMG) algorithm with finite difference method to solve this complex set of PDEs.

Materials and Methods

Desmoplastic Tumor Model

Consider a tumor growing in a tissue domain Ω ⊂ ℝ³ where tumor and healthy cells, as well as the ECM, are tracked with continuous volume fractions. In the diffuse interface model, adhesive forces hold the tumor cells together, creating a boundary layer of finite thickness between the tumoral and healthy regions. An appropriate distribution of ECM across the domain can be achieved by a carefully chosen free energy term. Dimensionless governing equations for multispecies tumor growth systems are presented here.

We assume that the liquid (extracellular fluid) volume fraction stays constant, as well as the total solid (ECM, tumor and healthy cells) volume fraction, and set the densities of all components to unity. Dimensionless variables ${\tilde{\varphi }}_V$, ${\tilde{\varphi }}_D$, ${\tilde{\varphi }}_E$, and ${\tilde{\varphi }}_H$ are the normalized volume fraction (normalized by the total volume fraction) of viable tumor cells, dead tumor cells, ECM, and healthy host cells, respectively. Their transient diffusion-convection governing equations are given as Ng and Frieboes (2017):

Cells and ECM Components

$\begin{array}{ll}\begin{array}{l}\frac{\partial {\tilde{\varphi }}_V}{\partial \tilde{t}}+\nabla ·({\tilde{\varphi }}_V{\tilde{\text{u}}}_{\alpha })=\nabla ·({\tilde{M}}_V\nabla {\tilde{\mu }}_T)+{\tilde{S}}_V\end{array}& (2.1)\end{array}$

$\begin{array}{ll}\begin{array}{l}\frac{\partial {\tilde{\varphi }}_D}{\partial \tilde{t}}+\nabla ·({\tilde{\varphi }}_D{\tilde{\text{u}}}_{\alpha })=\nabla ·({\tilde{M}}_D\nabla {\tilde{\mu }}_T)+{\tilde{S}}_D\end{array}& (2.2)\end{array}$

$\begin{array}{ll}\begin{array}{l}\frac{\partial {\tilde{\varphi }}_E}{\partial \tilde{t}}+\nabla ·({\tilde{\varphi }}_E{\tilde{\text{u}}}_{\alpha })=\nabla ·({\tilde{M}}_E\nabla {\tilde{\mu }}_E)+{\tilde{S}}_E\end{array}& (2.3)\end{array}$

$\begin{array}{ll}\begin{array}{l}{\tilde{\varphi }}_H=1-{\tilde{\varphi }}_V-{\tilde{\varphi }}_D-{\tilde{\varphi }}_E\end{array}& (2.4)\end{array}$

The total tumor volume fraction is ${\tilde{\varphi }}_T={\tilde{\varphi }}_V+{\tilde{\varphi }}_D$ and ${\tilde{M}}_i=\tilde{M}{\tilde{\varphi }}_i$ is the positive non-constant mobility of component i. Chemical potential terms are given below:

Chemical Potentials

$\begin{array}{ll}\begin{array}{l}{\tilde{\mu }}_T=\frac{\partial {\tilde{F}}_b}{\partial {\tilde{\varphi }}_T}-{{\widetilde{ϵ}}_T}^2{\nabla }^2{\tilde{\varphi }}_T-{{\widetilde{ϵ}}_{TE}}^2{\nabla }^2{\tilde{\varphi }}_E\end{array}& (2.5)\end{array}$

$\begin{array}{ll}\begin{array}{l}{\tilde{\mu }}_E=\frac{\partial {\tilde{F}}_b}{\partial {\tilde{\varphi }}_E}+\frac{\partial \tilde{{W}}}{\partial {\tilde{\varphi }}_E}-{{\widetilde{ϵ}}_E}^2{\nabla }^2{\tilde{\varphi }}_E-{{\widetilde{ϵ}}_{TE}}^2{\nabla }^2{\tilde{\varphi }}_T\end{array}& (2.6)\end{array}$

$\begin{array}{ll}\begin{array}{l}\frac{\partial {\tilde{F}}_b}{\partial {\tilde{\varphi }}_T}=2A_1{\tilde{\varphi }}_T(1-{\tilde{\varphi }}_T-{\tilde{\varphi }}_E)(1-2{\tilde{\varphi }}_T-{\tilde{\varphi }}_E)\\ +(A_5-A_3)(2{\tilde{\varphi }}_E-A_5-A_3)\end{array}& (2.7)\end{array}$

$\begin{array}{ll}\begin{array}{l}\frac{\partial {\tilde{F}}_b}{\partial {\tilde{\varphi }}_E}=2({\tilde{\varphi }}_T+A_2)({\tilde{\varphi }}_E-A_3)-2A_1{({\tilde{\varphi }}_T)}^2(1-{\tilde{\varphi }}_T-{\tilde{\varphi }}_E)\\ +({\tilde{\varphi }}_E-A_5)\left[2(1-{\tilde{\varphi }}_T+A_4)-3{\tilde{\varphi }}_E+A_5\right]\end{array}& (2.8)\end{array}$

$\begin{array}{ll}\begin{array}{l}\frac{\partial \tilde{{W}}}{\partial {\tilde{\varphi }}_E}={\widetilde{ϵ}}_e\left[6{\tilde{\varphi }}_E(1-{\tilde{\varphi }}_E)\right]{\displaystyle \sum _{i,j=1}^3}\left[\frac{1}{2}{({\tilde{{E}}}_T)}_{ij}{\widetilde{𝕋}}_{ij}^{{}^{*}}-{({\tilde{{E}}}_T^{{}^{*}})}_{ij}{\widetilde{𝕋}}_{ij}\right]\end{array}& (2.9)\end{array}$

The bulk free energy term ${\tilde{F}}_b$ is adapted from a tertiary semi-immiscible system described by Kim and Lowengrub (2005), where A₁ – A₅ are constants. The dimensionless elastic strain energy term $\tilde{{W}}$ used in Equation (2.9), which follows the form given by Leo et al. (1998) and Garcke (2005), is computed from the elements given below.

Elastic Energy

$\begin{array}{ll}\begin{array}{l}{\widetilde{𝕋}}_{mn}=2{\tilde{L}}_2{({\tilde{{E}}}_T)}_{mn}+{\tilde{L}}_1{\delta }_{mn}{\displaystyle \sum _{s=1}^3}{({\tilde{{E}}}_T)}_{ss}\end{array}& (2.10)\end{array}$

$\begin{array}{ll}\begin{array}{l}{({\tilde{{E}}}_T)}_{mn}={\tilde{{E}}}_{mn}-{\tilde{{E}}}_{mn}^{{}^{*}}\end{array}& (2.11)\end{array}$

$\begin{array}{ll}\begin{array}{l}{({\tilde{{E}}}_T^{{}^{*}})}_{mn}={({\tilde{{E}}}_E^{{}^{*}})}_{mn}-{({\tilde{{E}}}_C^{{}^{*}})}_{mn}\end{array}& (2.12)\end{array}$

$\begin{array}{ll}\begin{array}{l}{\tilde{{E}}}_{mn}^{{}^{*}}=Q_3({\tilde{\varphi }}_E){({\tilde{{E}}}_T^{{}^{*}})}_{mn}+{({\tilde{{E}}}_C^{{}^{*}})}_{mn}\end{array}& (2.13)\end{array}$

$\begin{array}{ll}\begin{array}{l}{\widetilde{𝕋}}_{mn}^{{}^{*}}=2(1-{\tilde{L}}_2^C){({\tilde{{E}}}_T)}_{mn}\\ +({\tilde{L}}_1^E-{\tilde{L}}_1^C){\delta }_{mn}{\displaystyle \sum _{k=1}^3}{({\tilde{{E}}}_T)}_{kk}\end{array}& (2.14)\end{array}$

$\begin{array}{ll}\begin{array}{l}{\tilde{{E}}}_{mn}=\frac{1}{2}\left(\frac{\partial {ũ}_m}{\partial {\tilde{x}}_n}+\frac{\partial {ũ}_n}{\partial {\tilde{x}}_m}\right)\end{array}& (2.15)\end{array}$

$\begin{array}{ll}\begin{array}{l}{\tilde{L}}_i=Q_3({\tilde{\varphi }}_E)({\tilde{L}}_i^E-{\tilde{L}}_i^C)+{\tilde{L}}_i^{C}wherei=1,2\end{array}& (2.16)\end{array}$

where $\tilde{{E}}$ is the infinitesimal strain, ${\tilde{\text{u}}}_{\text{d}}$ is the displacement vector, ${\tilde{{E}}}_{\text{E}}^{*}$ and ${\tilde{{E}}}_{\text{C}}^{*}$ are the Eigenstrain tensors for ECM and cells, respectively. In Equations (2.10) and (2.14), δ_mn = 1 for m = n and δ_mn = 0 for m≠n. ${\tilde{L}}_1^E$, ${\tilde{L}}_2^E$, ${\tilde{L}}_1^C$, and ${\tilde{L}}_2^C$ are Lamé constants for the ECM and cell components. The displacement vector ${\tilde{\text{u}}}_{\text{d}}$ is solved by setting $\nabla ·{\widetilde{𝕋}}_i=0$, where ${\widetilde{𝕋}}_i=\begin{array}{c}\left[{\widetilde{𝕋}}_{1i}{\widetilde{𝕋}}_{2i}{\widetilde{𝕋}}_{3i}\right]\\ \end{array}$ and i = 1, 2, 3.

The dimensionless cell-ECM phase pressure $\tilde{p}$ and interstitial fluid phase pressure $\tilde{q}$, with their corresponding velocities ${\tilde{\text{u}}}_{\alpha }$ and ${\tilde{\text{u}}}_{\beta }$ are computed using the following equations:

Pressures and Velocities

$\begin{array}{ll}\begin{array}{l}\nabla ·\left[{\tilde{k}}_{\alpha }\left(\nabla \tilde{p}-\frac{{\tilde{\gamma }}_T}{{\widetilde{ϵ}}_T}{\tilde{\mu }}_T\nabla {\tilde{\varphi }}_T-\frac{{\tilde{\gamma }}_E}{{\widetilde{ϵ}}_E}{\tilde{\mu }}_E\nabla {\tilde{\varphi }}_E\right)\right]\\ =-({\tilde{S}}_V+{\tilde{S}}_D+{\tilde{S}}_E)\end{array}& (2.17)\end{array}$

$\begin{array}{ll}\begin{array}{l}{\nabla }^2\tilde{q}=\frac{R_{\alpha ,\beta }}{{\tilde{k}}_{\beta }}({\tilde{S}}_V+{\tilde{S}}_D+{\tilde{S}}_E)\end{array}& (2.18)\end{array}$

$\begin{array}{ll}\begin{array}{l}{\tilde{\text{u}}}_{\alpha }=-{\tilde{k}}_{\alpha }\left[\nabla \tilde{p}-\frac{{\tilde{\gamma }}_T}{{\widetilde{ϵ}}_T}{\tilde{\mu }}_T\nabla {\tilde{\varphi }}_T-\frac{{\tilde{\gamma }}_E}{{\widetilde{ϵ}}_E}{\tilde{\mu }}_E\nabla {\tilde{\varphi }}_E\right]\end{array}& (2.19)\end{array}$

$\begin{array}{ll}\begin{array}{l}{\tilde{\text{u}}}_{\beta }=-{\tilde{k}}_{\beta }\nabla \tilde{q}\end{array}& (2.20)\end{array}$

$\begin{array}{ll}\begin{array}{l}{\tilde{\text{u}}}_E={\tilde{\text{u}}}_{\alpha }-\tilde{M}\nabla ({\tilde{\mu }}_E)\end{array}& (2.21)\end{array}$

where ${\tilde{k}}_{\alpha }=f({\tilde{\varphi }}_T,{\tilde{\varphi }}_E)$ and ${\tilde{k}}_{\beta }$ are motilities of the solid and liquid phase respectively, and ${\tilde{\text{u}}}_E$ is the ECM component velocity.

Dimensionless nutrients and waste products concentrations ñ, $\tilde{g}$, $\tilde{w}$, $\tilde{\updownarrow }$, $\tilde{b}$, ã, $\tilde{s}$, and $\tilde{r}$ represent O₂, glucose, CO₂, lactate, bicarbonate, H⁺, Na⁺, and Cl⁻ respectively. We use the following quasi-steady state governing equations for tracking nutrients and waste products within the tissue domain:

Nutrients and Waste Products

$\begin{array}{ll}\begin{array}{l}\nabla ·({\tilde{D}}_n\nabla ñ)+{\tilde{k}}_{n1}{ñ}_C-({\tilde{k}}_{n1}+{\tilde{k}}_{n2})ñ=0\end{array}& (2.22)\end{array}$

$\begin{array}{ll}\begin{array}{l}\nabla ·({\tilde{D}}_g\nabla \tilde{g})+{\tilde{k}}_{g1}{\tilde{g}}_C-({\tilde{k}}_{g1}+{\tilde{k}}_{g2})\tilde{g}=0\end{array}& (2.23)\end{array}$

$\begin{array}{ll}\begin{array}{l}\nabla ·({\tilde{D}}_w\nabla \tilde{w})+{\tilde{k}}_{n2}ñ+{\tilde{k}}_r\tilde{b}ã+{\tilde{k}}_w{\tilde{w}}_C-({\tilde{k}}_f+{\tilde{k}}_w)\tilde{w}=0\end{array}& (2.24)\end{array}$

$\begin{array}{ll}\begin{array}{l}\nabla ·\{{\tilde{D}}_{\ell }[\nabla \tilde{\ell }\\ -{\tilde{z}}_{\ell }\tilde{\ell }\left(\frac{{\tilde{z}}_{\ell }{\tilde{D}}_{\ell }\nabla \tilde{\ell }+{\tilde{z}}_b{\tilde{D}}_b\nabla \tilde{b}+{\tilde{D}}_a\nabla \tilde{a}+{\tilde{z}}_s{\tilde{D}}_s\nabla \tilde{s}+{\tilde{z}}_r{\tilde{D}}_r\nabla \tilde{r}}{{\tilde{z}}_{\ell }^2{\tilde{D}}_{\ell }\tilde{\ell }+{\tilde{z}}_b^2{\tilde{D}}_b\tilde{b}+{\tilde{D}}_a\tilde{a}+{\tilde{z}}_s^2{\tilde{D}}_s\tilde{s}+{\tilde{z}}_r^2{\tilde{D}}_r\tilde{r}}\right)]\}\\ +2R_{g,n}\left({\tilde{k}}_{g2}\tilde{g}\right)-\frac{1}{3}\left({\tilde{k}}_{n2}\tilde{n}\right)+{\tilde{k}}_{\ell }{\tilde{\ell }}_C-{\tilde{k}}_{\ell }\tilde{\ell }=0\end{array}& (2.25)\end{array}$

$\begin{array}{ll}\begin{array}{l}\nabla ·\{{\tilde{D}}_b[\nabla \tilde{b}\\ -{\tilde{z}}_b\tilde{b}\left(\frac{{\tilde{z}}_{\ell }{\tilde{D}}_{\ell }\nabla \tilde{\ell }+{\tilde{z}}_b{\tilde{D}}_b\nabla \tilde{b}+{\tilde{D}}_a\nabla \tilde{a}+{\tilde{z}}_s{\tilde{D}}_s\nabla \tilde{s}+{\tilde{z}}_r{\tilde{D}}_r\nabla \tilde{r}}{{\tilde{z}}_{\ell }^2{\tilde{D}}_{\ell }\tilde{\ell }+{\tilde{z}}_b^2{\tilde{D}}_b\tilde{b}+{\tilde{D}}_a\tilde{a}+{\tilde{z}}_s^2{\tilde{D}}_s\tilde{s}+{\tilde{z}}_r^2{\tilde{D}}_r\tilde{r}}\right)]\}\\ +{\tilde{k}}_f\tilde{w}-{\tilde{k}}_r\tilde{b}\tilde{a}=0\end{array}& (2.26)\end{array}$

$\begin{array}{ll}\begin{array}{l}\nabla ·\{{\tilde{D}}_a[\nabla \tilde{a}\\ -{\tilde{z}}_a\tilde{a}\left(\frac{{\tilde{z}}_{\ell }{\tilde{D}}_{\ell }\nabla \tilde{\ell }+{\tilde{z}}_b{\tilde{D}}_b\nabla \tilde{b}+{\tilde{D}}_a\nabla \tilde{a}+{\tilde{z}}_s{\tilde{D}}_s\nabla \tilde{s}+{\tilde{z}}_r{\tilde{D}}_r\nabla \tilde{r}}{{\tilde{z}}_{\ell }^2{\tilde{D}}_{\ell }\tilde{\ell }+{\tilde{z}}_b^2{\tilde{D}}_b\tilde{b}+{\tilde{D}}_a\tilde{a}+{\tilde{z}}_s^2{\tilde{D}}_s\tilde{s}+{\tilde{z}}_r^2{\tilde{D}}_r\tilde{r}}\right)]\}\\ +2R_{g,n}\left({\tilde{k}}_{g2}\tilde{g}\right)\text{}-\text{}\frac{1}{3}\left({\tilde{k}}_{n2}\tilde{n}\right)\text{}+\text{}{\tilde{k}}_f\tilde{w}\text{}-\text{}{\tilde{k}}_r\tilde{b}\tilde{a}\text{}+\text{}{\tilde{k}}_{\ell }{\tilde{\ell }}_C\text{}-\text{}{\tilde{k}}_{\ell }\tilde{\ell }\text{}=\text{}0\end{array}& (2.27)\end{array}$

$\begin{array}{ll}\begin{array}{l}\nabla ·\{{\tilde{D}}_s[\nabla \tilde{s}\\ -{\tilde{z}}_s\tilde{s}\left(\frac{{\tilde{z}}_{\ell }{\tilde{D}}_{\ell }\nabla \tilde{\ell }+{\tilde{z}}_b{\tilde{D}}_b\nabla \tilde{b}+{\tilde{D}}_a\nabla \tilde{a}+{\tilde{z}}_s{\tilde{D}}_s\nabla \tilde{s}+{\tilde{z}}_r{\tilde{D}}_r\nabla \tilde{r}}{{\tilde{z}}_{\ell }^2{\tilde{D}}_{\ell }\tilde{\ell }+{\tilde{z}}_b^2{\tilde{D}}_b\tilde{b}+{\tilde{D}}_a\tilde{a}+{\tilde{z}}_s^2{\tilde{D}}_s\tilde{s}+{\tilde{z}}_r^2{\tilde{D}}_r\tilde{r}}\right)]\}\\ =0\end{array}& (2.28)\end{array}$

$\begin{array}{ll}\begin{array}{l}\tilde{r}=-\frac{1}{{\tilde{z}}_r}\left({\tilde{z}}_{\ell }\tilde{\ell }+{\tilde{z}}_b\tilde{b}+{\tilde{z}}_aã+{\tilde{z}}_s\tilde{s}\right)\end{array}& (2.29)\end{array}$

The flux terms of charged species follow those given by Casciari et al. (1992). The terms ${\tilde{k}}_{n1}$, ${\tilde{k}}_{n2}$, ${\tilde{k}}_{g1}$, ${\tilde{k}}_{g2}$, ${\tilde{k}}_w$, and ${\tilde{k}}_{\updownarrow }$ are combined rate constants used in source terms of nutrients and waste products, whereas ${\tilde{k}}_f$ and ${\tilde{k}}_r$ are the forward and reverse rate constants of the reaction $C{O}_2+H_{2}O\leftrightarrow HC{O}_3^{-}+H^{+}$.

Dimensionless concentrations of tumor growth factors, tumor angiogenic factors, matrix degrading enzymes, and myofibroblastic cells are represented by $t\tilde{g}f$, $t\tilde{a}f$, $\tilde{m}$, and ${\tilde{F}}_E$ respectively. Quasi-steady state equations for $t\tilde{g}f$ and $t\tilde{a}f$, as well as transient governing equations for $\tilde{m}$ and ${\tilde{F}}_E$ are given below:

Tumorigenic Species

$\begin{array}{ll}\begin{array}{l}\nabla ·({\tilde{D}}_{tgf}\nabla t\tilde{g}f)+{\tilde{\lambda }}_{tgf}\\ -({\tilde{\lambda }}_{tgf}+{\tilde{\lambda }}_{de,tgf}+{\tilde{\lambda }}_{U,tgf})t\tilde{g}f=0\end{array}& (2.30)\end{array}$

$\begin{array}{ll}\begin{array}{l}\nabla ·({\tilde{D}}_{taf}\nabla t\tilde{a}f)+{\tilde{\lambda }}_{taf}\\ -({\tilde{\lambda }}_{taf}+{\tilde{\lambda }}_{de,taf}+{\tilde{\lambda }}_{U,taf})t\tilde{a}f=0\end{array}& (2.31)\end{array}$

$\begin{array}{ll}\begin{array}{l}\frac{\partial \tilde{m}}{\partial \tilde{t}}=\nabla ·({\tilde{D}}_m\nabla \tilde{m})+{\tilde{S}}_m\end{array}& (2.32)\end{array}$

$\begin{array}{ll}\begin{array}{l}\frac{\partial {\tilde{F}}_E}{\partial t}+\nabla ·({\tilde{F}}_E{\tilde{\text{u}}}_E)=-\nabla ·({\tilde{D}}_F{\tilde{F}}_E\nabla t\tilde{g}f)+{\tilde{S}}_{FE}\end{array}& (2.33)\end{array}$

Here, ${\tilde{F}}_E$ reflects the concentration within the ECM phase, assuming the volume of myofibroblastic cells is negligible and the ECM phase is continuous throughout the domain. Similarly, dimensionless ECM based concentrations of blood and lymphatic vessels, ${\tilde{B}}_n^E$ and ${\tilde{L}}_n^E$ respectively, with their corresponding diffusive flux terms are given by.

Blood and Lymphatic Vessels

$\begin{array}{ll}\begin{array}{l}\frac{\partial {\tilde{B}}_n^E}{\partial t}+\nabla ·({\tilde{B}}_n^E{\tilde{\text{u}}}_E)=-\nabla ·{\tilde{\text{J}}}_{BnE}+{\tilde{S}}_{BnE}\end{array}& (2.34)\end{array}$

$\begin{array}{ll}\begin{array}{l}\frac{\partial {\tilde{L}}_n^E}{\partial t}+\nabla ·({\tilde{L}}_n^E{\tilde{\text{u}}}_E)=-\nabla ·{\tilde{\text{J}}}_{LnE}+{\tilde{S}}_{LnE}\end{array}& (2.35)\end{array}$

$\begin{array}{ll}\begin{array}{l}{\tilde{J}}_{BnE}={\tilde{\chi }}_{che,BnE}{{A}}_{che,BnE}{\tilde{B}}_n^E\nabla t\tilde{a}f\\ +{\tilde{\chi }}_{hap,BnE}{{A}}_{hap,BnE}{\tilde{B}}_n^E\nabla {\tilde{\varphi }}_E-{\tilde{D}}_{BnE}\nabla {\tilde{B}}_n^E\end{array}& (2.36)\end{array}$

$\begin{array}{ll}\begin{array}{l}{\tilde{\text{J}}}_{LnE}={\tilde{\chi }}_{che,LnE}{{A}}_{che,LnE}{\tilde{L}}_n^E\nabla t\tilde{a}f\\ +{\tilde{\chi }}_{hap,LnE}{{A}}_{hap,LnE}{\tilde{L}}_n^E\nabla {\tilde{\varphi }}_E-{\tilde{D}}_{LnE}\nabla {\tilde{L}}_n^E\end{array}& (2.37)\end{array}$

Tissue effective diffusivities ${\tilde{D}}_n$, ${\tilde{D}}_g$, ${\tilde{D}}_w$, ${\tilde{D}}_{\updownarrow }$, ${\tilde{D}}_b$, ${\tilde{D}}_a$, ${\tilde{D}}_s$, ${\tilde{D}}_{tgf}$, ${\tilde{D}}_{taf}$, ${\tilde{D}}_m$, ${\tilde{D}}_{FE}$, ${\tilde{D}}_{BnE}$, ${\tilde{D}}_{LnE}$, and tissue effective mass transfer coefficients ${\tilde{\lambda }}_{B,n}$, ${\tilde{\lambda }}_{B,g}$, ${\tilde{\lambda }}_{B,w}$, ${\tilde{\lambda }}_{B,{l}}$, are represented by $\tilde{\psi }$ and determined as follow:

$\begin{array}{ll}\begin{array}{l}\tilde{\psi }={\tilde{\psi }}_E{Q}_3({\tilde{\varphi }}_E)\\ +\left[1-Q_3({\tilde{\varphi }}_E)\right]\left\{{\tilde{\psi }}_T{Q}_3(\frac{{\tilde{\varphi }}_T}{{\tilde{\varphi }}_C})+{\tilde{\psi }}_H\left[1-Q_3(\frac{{\tilde{\varphi }}_T}{{\tilde{\varphi }}_C})\right]\right\}\end{array}& (2.38)\end{array}$

where ${\tilde{\psi }}_E$, ${\tilde{\psi }}_T$, and ${\tilde{\psi }}_H$ are dimensionless diffusivity or transfer coefficient in the ECM, tumor, and healthy-host cell domain respectively. The total cell volume fraction is ${\tilde{\varphi }}_C={\tilde{\varphi }}_T+{\tilde{\varphi }}_H$. The full list of parameters, constants, factors used in non-dimensionalization, source terms, rate terms, and adjustment factors can be found in Supplementary Tables 1–6.

We define the following Neumann boundary conditions for the cell and ECM volume fractions, and Dirichlet boundary conditions for solid cell pressure and chemical potentials at all external boundaries:

$\begin{array}{ll}\begin{array}{l}n·\nabla {\tilde{\varphi }}_{\text{V}}=n·\nabla {\tilde{\varphi }}_{\text{D}}=n·\nabla {\tilde{\varphi }}_{\text{E}}=0,\\ \tilde{p}={\tilde{\mu }}_T={\tilde{\mu }}_E={\tilde{\mu }}_H=0\end{array}& (2.39)\end{array}$

where n is the outward normal of a boundary. For nutrients and waste products, as well as tumorigenic species, Dirichlet boundary conditions are imposed, with the exception of myofibroblastic cell species, where Neumann boundary conditions are applied:

$\begin{array}{ll}\begin{array}{l}\tilde{n}=\tilde{g}=1,\\ \tilde{w}=\tilde{\ell }=\tilde{b}=\tilde{a}=\tilde{s}=t\tilde{g}f=t\tilde{a}f=\tilde{m}=0,\\ n·\nabla {\tilde{F}}_E=0.\end{array}& (2.40)\end{array}$

Blood and lymphatic vessels are assumed to be at their corresponding far-field values at external boundaries:

$\begin{array}{ll}\begin{array}{l}{\tilde{B}}_n^E={\tilde{B}}_{\infty },{\tilde{L}}_n^E={\tilde{L}}_{\infty },\end{array}& (2.41)\end{array}$

where ${\tilde{B}}_{\infty }={\tilde{L}}_{\infty }=0.2$ is used here.

Numerical Methods

The mathematical model is first discretized in time and space. After reorganizing the discretized equations, a non-linear relaxation procedure is employed.

Computational Domain

We consider a rectangular 3D domain Ω = (L_x, R_x) × (L_y, R_y) × (L_z, R_z). Let the domain be discretized into N_x × N_y × N_z cells. Bounded by boundary cells, the domain is covered by the following sets of cell centers, cell edge points, and cell grid points:

$\begin{array}{ll}\begin{array}{l}C=\left\{(x_i,y_j,z_k)|0\le i\le N_x+1,0\le j\le N_y+1,0\le k\le N_z+1\right\}\end{array}& (\mathrm{2.2.1})\end{array}$

$\begin{array}{ll}\begin{array}{l}{E}_{ew}=\left\{(x_{i+\frac{1}{2}},y_j,z_k)|0\le i\le N_x,1\le j\le N_y,1\le k\le N_z\right\},\end{array}& (\mathrm{2.2.2})\end{array}$

$\begin{array}{ll}\begin{array}{l}{E}_{ns}=\left\{(x_i,y_{j+\frac{1}{2}},z_k)|1\le i\le N_x,0\le j\le N_y,1\le k\le N_z,\right\},\end{array}& (\mathrm{2.2.3})\end{array}$

$\begin{array}{ll}\begin{array}{l}{E}_{tb}=\left\{(x_i,y_j,z_{k+\frac{1}{2}})|1\le i\le N_x,1\le j\le N_y,0\le k\le N_z,\right\},\end{array}& (\mathrm{2.2.4})\end{array}$

$\begin{array}{ll}\begin{array}{l}G=\left\{(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}},z_{k+\frac{1}{2}})|0\le i\le N_x,0\le j\le N_y,0\le k\le N_z,\right\}\end{array}& (\mathrm{2.2.5})\end{array}$

where E_ew, E_ns, and E_tb represent east-west, north-south, and top-bottom sets of cell edges respectively, and G represents cell corners. Let grid spaces be

$\begin{array}{ll}\begin{array}{l}\Delta x_i=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}},\\ \Delta y_j=\Delta y_{j+\frac{1}{2}}-y_{j-\frac{1}{2}},\\ \Delta z_k=z_{k+\frac{1}{2}}-z_{k-\frac{1}{2}},\end{array}& (\mathrm{2.2.6})\end{array}$

and if partitions in the three directions are uniform and equal:

$\begin{array}{ll}\begin{array}{l}\Delta x_i=\Delta x_j=\Delta x_k=\eta ,\end{array}& (\mathrm{2.2.7})\end{array}$

$\begin{array}{ll}\begin{array}{l}\eta =\frac{R_x-L_x}{N_x}=\frac{R_y-L_y}{N_y}=\frac{R_z-L_z}{N_z},\end{array}& (\mathrm{2.2.8})\end{array}$

for 1 ≤ i ≤ N_x, 1 ≤ j ≤ N_y, and 1 ≤ k ≤ N_z, we get the following coordinates for cell centers and cell edges or corners defined in Equations (2.2.1)–(2.2.5):

$\begin{array}{ll}\begin{array}{l}{x}_i=L_x+\left(i-\frac{1}{2}\right)\eta ,y_j=L_y+\left(j-\frac{1}{2}\right)\eta ,\\ z_k=L_z+\left(k-\frac{1}{2}\right)\eta ,\end{array}& (\mathrm{2.2.9})\end{array}$

$\begin{array}{ll}\begin{array}{l}{x}_{i+\frac{1}{2}}=L_x+i\eta ,y_{j+\frac{1}{2}}=L_y+j\eta ,\\ z_{k+\frac{1}{2}}=L_z+k\eta ,\end{array}& (\mathrm{2.2.10})\end{array}$

following their corresponding ranges indicated in Equations (2.2.1)–(2.2.5).

The differential, Laplacian and flux terms associated with the model are in Supplementary Materials.

Model Discretization

The model consists of a set of stiff differential equations that are fourth-order in space. At time step a with time step size θ, they are discretized in time using the Crank-Nicolson Method as in Wise et al. (2011). Details of the discretization as well as the multigrid V-cycle iterations and non-linear Gauss–Seidel relaxations are given in Supplementary Materials.

The tumor model solved using the adaptive full multigrid V-cycle was coded in C and simulations were performed on a node equipped with 768 GB of RAM and 32 Intel Xeon 3.3 GHz cores running CentOS 6.7 x86_64. The algorithms were partially parallelized using OpenMP to achieve higher performance.

Results

Self-Adaptive Full Multigrid Full Approximation Scheme V-Cycle Algorithm

The self-adaptive full multigrid process combines FMG with the self-adaptive approach. The process involves constructing a multilevel block-structured mesh and seeking the solution on the new mesh structure. Mesh construction and refinement, as well as FAS multigrid in an adaptive FMG setting were developed as follows.

Multilevel Mesh Refinement

We start with a rectangular computational domain as described in Methods. Each level κ covers a domain of Ω_κ with mesh size η_κ. Grid levels are numbered as κ = κ_min, …, 0, …, κ_max, where κ_min represents the coarsest mesh level and κ_max the finest. Global grid levels, κ = 0 − κ_min, cover the entire computational domain Ω_κ = Ω. The finest global grid level, κ = 0, is referred as the root level. Grid levels κ = 1 to κ = κ_max are levels with refined mesh covering a domain of Ω_κ+1 ⊆ Ω_κ ⊆ Ω, each consists of n_κ rectangular blocks, B_{κ, 1}, …, B_κ,nκ, of uniform grids, G_{κ, 1}, …, G_κ,nκ. Note that n_κ = 1 for all global grid levels. Cell-centered discretization is used with grid spacing ${\eta }_{\kappa +1}=\frac{{\eta }_{\kappa }}{2}$. The hierarchy of levels and meshes are shown in Figure 1.

FIGURE 1

Figure 1. For simplification and ease of visualization, a two dimensional example is given here. The multigrid algorithm consists of a hierarchy of grid levels κ = −2, …, 2, shown here on the top rowwith their corresponding subdomains/blocks. Level κ = 0 is the root level and levels κ = 1 and 2 are refinement levels. Global levels κ = 0, −1, and −2 cover the entire computational domain, whereas refinement levels κ = 1 and 2 cover subdomains where refinement is needed. Shown in the example here, the locally refined block-structured κ = 1 consists of one block, B_{1, 1}, and κ = 2 consists of two blocks, B_{2, 1} and B_{2, 2}. The middle row shows the corresponding meshes for each level. For refinement levels in the example here, mesh G_{1, 1} for block B_{1, 1}, G_{2, 1} and G_{2, 2} for blocks B_{2, 1} and B_{2, 2}. The composite grids of all levels is shown in the last row.

After solution is obtained for the current finest level, we use the undivided gradient test (Wise et al., 2007) to tag cells for refinement. The undivided gradient test is used to capture the diffuse interface region (tumor-host cell domain). Let the critical value of some indicative variable ψ_κ(x, y, z) be C_κ, cell-centered coordinates on level κ be x_{i, j, k}, the following test is performed on cells on level κ:

$\begin{array}{ll}\begin{array}{l}{F}_{\kappa }=\begin{array}{c}\left\{x_{i,j,k}\in {\Omega }_{\kappa }|\sqrt{\begin{array}{c}{\left[{\psi }_{i+1,j,k}-{\psi }_{i-1,j,k}\right]}^2\\ +{\left[{\psi }_{i,j+1,k}-{\psi }_{i,j-1,k}\right]}^2\\ +{\left[{\psi }_{i,j,k+1}-{\psi }_{i,j,k-1}\right]}^2\end{array}}>C_{\kappa }\right\}\end{array}\end{array}& (\mathrm{3.1.1})\end{array}$

where F_κ are coordinates of flagged cells. Another criteria that can be used to tag cells for refinement is the relative truncation error test (Trottenberg et al., 2001). In the relative truncation error test, the relative truncation error with respect to Ω_κ and Ω_κ−1,

$\begin{array}{ll}\begin{array}{l}{\tau }_h^{2h}=L_{\kappa -1}(R_{\kappa }^{\kappa -1}{\psi }_{\kappa }^{a,r,{\nu }_1})-R_{\kappa }^{\kappa -1}{L}_{\kappa }({\psi }_{\kappa }^{a,r,{\nu }_1})\end{array}& (\mathrm{3.1.2})\end{array}$

as used in Equation (3.1.20) in the FAS cycle, is used to flag cells for refinement. A subroutine is created and called to perform this task:

$\begin{array}{ll}\begin{array}{l}{F}_{\kappa }^{a,n}=\mathrm{\text{FLAG}}({\psi }_{\kappa }^{a,r})\end{array}& (\mathrm{3.1.3})\end{array}$

where the current finest-grid data is used and $F_{\kappa }^{a,n}$ is a list of flagged coordinates.

The flagged cells are then rearranged into patches of rectangular refined mesh in the next finer level using a clustering algorithm from Bell et al. (1994) and Berger and Rigoutsos (1991) with minor modification. Create a list of blocks and start with just one block containing all flagged cells given by $F_{\kappa }^{a,n}$. The following procedure is used to divide flagged cells into blocks:

(1) Compute efficiency of the block, which is the ratio of flagged cells to total cells. If the efficiency is below the threshold efficiency, and, if the size of the block is more than twice the threshold size, then continue. Else, accept and add the current block to the list, return to Step (1) for the next block on the list.

(2) Compute signatures, which is the number of flagged cells in each slice along each direction. Checking for gaps along all dimensions starting from the longest edge, find a gap closest to the center of that direction. If gaps exist and an optimum gap location is found, slice the block into two along the gap, and go to Step (5). If no gaps are found, continue.

(3) Compute second derivatives of the signatures along all directions. Checking for inflection points starting from the longest edge, find an inflection point closest to the center of that direction. If inflection points are found and an optimum inflection location is found, slice the block into two along the inflection point, and go to Step (5). If no inflection points are found, continue.

(4) If the efficiency of the block is above a minimum allowed, accept and add the current block to the list, returning to Step (1) for the next block on the list. Else, divide the block into two along the mid-point of the longest dimension and continue.

(5) Delete all empty slices along the edges of each block, and add the two trimmed blocks to the list. Repeat from Step (1) until all blocks created are checked.

The above task is assigned to the following subroutine of block generation:

$\begin{array}{ll}\begin{array}{l}{B}_{\kappa +1}=\mathrm{\text{BLOCKGEN}}(threshold\text{\_}eff,threshold\text{\_}size,min\text{\_}eff,F_{\kappa }^{a,n})\end{array}& (\mathrm{3.1.4})\end{array}$

where B_κ+1 is an array of blocks on the refined level κ + 1, with rows of coordinates corresponding to corners of each block.

Following the generation of new refined grids on level κ + 1, the new grids are populated with data from level κ and the old level κ + 1. All cell-centered data for variables on the newly refined grids must be generated by higher order interpolation, such as cubic interpolation given by Equations (3.1.28) and (3.1.29), from the coarse grids below. If a new level κ + 1 grid cell overlaps any old level κ + 1 grid cell, the previous time-step data for the new grid cell is copied from the old. However, for any new level κ + 1 grid cell that does not overlap with an old κ + 1 grid cell, its previous time-step data has to be obtained from the coarse grid below. For all old level κ + 1 grid cells that do not overlap with any new κ + 1 grid cells, data are averaged and stored in the coarse level κ grid. Refer to Figure 2 for illustration.

FIGURE 2

Figure 2. Populating a newly refined grid. For simplification and ease of visualization, a two dimensional example is illustrated here. An existing block, with its boundary marked by - -, is overlapped by a newly refined block with boundary marked by . All cell-centered variable data on the newly refined grid are obtained via cubic interpolation from the coarse grid data below. The previous time-step data for the new grid points (•), which overlap old grid points, are copied directly from the old grid points they overlap; for new grid points (▴) that do not overlap with any old grid, their previous time-step data are obtained from the coarse grid below. The remaining data on the old grid (°) that do not overlap any of the new grid points are averaged and copied to the coarse grid data points (×).

After populating the newly refined grid, ghost cells surrounding each grid patch must be constructed. We first compute data for the ghost cells using the Π quadratic interpolation given by Colella et al. (2009). Illustrated in Figure 3, quadratic interpolation is performed twice. First, to get data from the coarse-grid for the intermediate points a and b:

$\begin{array}{ll}\begin{array}{ll}{I}_{\kappa -1}^{\kappa }{\psi }_{\kappa -1}\left(x,y\right)=& \frac{1}{32}[30{\psi }_{\kappa -1}\left(x,y-\frac{{\eta }_{\kappa }}{2}\right)\\ & +5{\psi }_{\kappa -1}\left(x,y+\frac{3{\eta }_{\kappa }}{2}\right)]\\ & -3{\psi }_{\kappa -1}\left(x,y-\frac{5{\eta }_{\kappa }}{2}\right)],\end{array}& (\mathrm{3.1.5})\end{array}$

$\begin{array}{ll}\begin{array}{ll}{I}_{\kappa -1}^{\kappa }{\psi }_{\kappa -1}\left(x,y\right)=& \frac{1}{32}[30{\psi }_{\kappa -1}\left(x,y+\frac{{\eta }_{\kappa }}{2}\right)\\ & +5{\psi }_{\kappa -1}\left(x,y-\frac{3{\eta }_{\kappa }}{2}\right)]\\ & -3{\psi }_{\kappa -1}\left(x,y+\frac{5{\eta }_{\kappa }}{2}\right)],\end{array}& (\mathrm{3.1.6})\end{array}$

respectively. Then, to compute fine-grid boundary data on ghost cell c, perform quadratic interpolation using points a, e, and g. Similarly, interpolate using points b, f, and h to obtain data at ghost cell d. Following Wise et al. (2007), all ghost cells of each fine-grid patch are first filled using the quadratic interpolation above. If a ghost cell falls on a cell on any neighboring fine-grid patches of the same level, then the ghost cell data is replaced by the more accurate cell-centered data from the cell of the neighboring patch.

FIGURE 3

Figure 3. Interpolation at fine-grid ghost layers. For simplification and ease of visualization, a two dimensional illustration is presented. The coarse grid mesh is shown using thick solid lines, and the finer grid mesh is shown using thin solid lines. The internal cell-centered grid points of the finer grid are depicted by (▴) and its external/ghost grid points are represented by (°). Only relevant cell-centered grid points (♦) on the coarse grids are shown. First, data on intermediate points (×) are obtained using quadratic interpolation from their nearest three coarse grid points (♦). Corner ghost cells (°) are obtained by quadratic interpolation using the nearest coarse-grid data (♦) and two fine-grid data (▴) on the diagonal, shown by green dashed lines. Remaining ghost cells (°) are obtained from the intermediate points (×) and two fine-grid interior points (▴), shown by red dashed lines.

Adaptive FAS V-Cycle

Consider a case with κ grid levels with κ = 0 being the root level. The levels κ = κ_max and κ = κ_min represent the level with the finest and coarsest grid, respectively. Each level of κ > 0 contains one or more blocks of refined rectangular grids covering a subdomain Ω_κ, whereas levels κ ≤ 0 cover the entire computational domain Ω_κ = Ω.

Averaging is used for the restriction operator $R_{\kappa }^{\kappa -1}$:

$\begin{array}{ll}\begin{array}{ll}{R}_{\kappa }^{\kappa -1}{\psi }_{\kappa }\left(x,y,z\right)=& \frac{1}{8}[{\psi }_{\kappa }\left(x-\frac{{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2},z-\frac{{\eta }_{\kappa }}{2}\right)\\ & +{\psi }_{\kappa }\left(x-\frac{{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2},z-\frac{{\eta }_{\kappa }}{2}\right)\\ & +{\psi }_{\kappa }\left(x+\frac{{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2},z-\frac{{\eta }_{\kappa }}{2}\right)\\ & +{\psi }_{\kappa }\left(x+\frac{{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2},z-\frac{{\eta }_{\kappa }}{2}\right)\\ & +{\psi }_{\kappa }\left(x-\frac{{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2},z+\frac{{\eta }_{\kappa }}{2}\right)\\ & +{\psi }_{\kappa }\left(x-\frac{{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2},z+\frac{{\eta }_{\kappa }}{2}\right)\\ & +{\psi }_{\kappa }\left(x+\frac{{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2},z+\frac{{\eta }_{\kappa }}{2}\right)\\ & +{\psi }_{\kappa }\left(x\text{}+\text{}\frac{{\eta }_{\kappa }}{2},y\text{}+\text{}\frac{{\eta }_{\kappa }}{2},z\text{}+\text{}\frac{{\eta }_{\kappa }}{2}\right)]\end{array}& (\mathrm{3.1.7})\end{array}$

which reduces to a four-point average in a 2D case. Linear interpolation is used for the prolongation operators $P_{\kappa -1}^{\kappa }$ in error correction steps within a V-cycle. In the 2D cell-centered discretization cases as depicted in Figure 4, the bilinear interpolation operators are given by Trottenberg et al. (2001):

$\begin{array}{ll}\begin{array}{ll}{P}_{\kappa -1}^{\kappa }{\psi }_{\kappa -1}\left(x,y\right)=& \frac{1}{16}[9{\psi }_{\kappa -1}\left(x-\frac{{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2}\right)\\ & +3{\psi }_{\kappa -1}\left(x-\frac{{\eta }_{\kappa }}{2},y-\frac{3{\eta }_{\kappa }}{2}\right)\\ & +3{\psi }_{\kappa -1}\left(x+\frac{3{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2}\right)\\ & +{\psi }_{\kappa -1}\left(x+\frac{3{\eta }_{\kappa }}{2},y-\frac{3{\eta }_{\kappa }}{2}\right)]\end{array}& (\mathrm{3.1.8})\end{array}$

$\begin{array}{ll}\begin{array}{ll}{P}_{\kappa -1}^{\kappa }{\psi }_{\kappa -1}\left(x,y\right)=& \frac{1}{16}[9{\psi }_{\kappa -1}\left(x+\frac{{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2}\right)\\ & +3{\psi }_{\kappa -1}\left(x+\frac{{\eta }_{\kappa }}{2},y-\frac{3{\eta }_{\kappa }}{2}\right)\\ & +3{\psi }_{\kappa -1}\left(x-\frac{3{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2}\right)\\ & +{\psi }_{\kappa -1}\left(x-\frac{3{\eta }_{\kappa }}{2},y-\frac{3{\eta }_{\kappa }}{2}\right)]\end{array}& (\mathrm{3.1.9})\end{array}$

$\begin{array}{ll}\begin{array}{ll}{P}_{\kappa -1}^{\kappa }{\psi }_{\kappa -1}\left(x,y\right)=& \frac{1}{16}[9{\psi }_{\kappa -1}\left(x-\frac{{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2}\right)\\ & +3{\psi }_{\kappa -1}\left(x-\frac{{\eta }_{\kappa }}{2},y+\frac{3{\eta }_{\kappa }}{2}\right)\\ & +3{\psi }_{\kappa -1}\left(x+\frac{3{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2}\right)\\ & +{\psi }_{\kappa -1}\left(x+\frac{3{\eta }_{\kappa }}{2},y+\frac{3{\eta }_{\kappa }}{2}\right)]\end{array}& (\mathrm{3.1.10})\end{array}$

$\begin{array}{ll}\begin{array}{ll}{P}_{\kappa -1}^{\kappa }{\psi }_{\kappa -1}\left(x,y\right)=& \frac{1}{16}[9{\psi }_{\kappa -1}\left(x+\frac{{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2}\right)\\ & +3{\psi }_{\kappa -1}\left(x+\frac{{\eta }_{\kappa }}{2},y+\frac{3{\eta }_{\kappa }}{2}\right)\\ & +3{\psi }_{\kappa -1}\left(x-\frac{3{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2}\right)\\ & +{\psi }_{\kappa -1}\left(x-\frac{3{\eta }_{\kappa }}{2},y+\frac{3{\eta }_{\kappa }}{2}\right)]\end{array}& (\mathrm{3.1.11})\end{array}$

for points marked a, b, c, and d, respectively. For 3D cases as shown in Figure 5, trilinear interpolation produces the following

$\begin{array}{ll}\begin{array}{ll}{P}_{\kappa -1}^{\kappa }{\psi }_{\kappa -1}\left(x,y,z\right)=& \frac{1}{64}[27{\psi }_{\kappa -1}\left(x+\frac{{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2}\right)\\ & +9{\psi }_{\kappa -1}\left(x+\frac{{\eta }_{\kappa }}{2},y+\frac{3{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2}\right)\\ & +9{\psi }_{\kappa -1}\left(x+\frac{{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2}\right)\\ & +9{\psi }_{\kappa -1}\left(x-\frac{3{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2}\right)\\ & +3{\psi }_{\kappa -1}\left(x+\frac{{\eta }_{\kappa }}{2},y+\frac{3{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2}\right)\\ & +3{\psi }_{\kappa -1}\left(x-\frac{3{\eta }_{\kappa }}{2},y+\frac{3{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2}\right)\\ & +3{\psi }_{\kappa -1}\left(x-\frac{3{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2}\right)\\ & +{\psi }_{\kappa -1}\left(x-\frac{3{\eta }_{\kappa }}{2},y+\frac{3{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2}\right)]\end{array}& (\mathrm{3.1.12})\end{array}$

$\begin{array}{ll}\begin{array}{ll}{P}_{\kappa -1}^{\kappa }{\psi }_{\kappa -1}(x,y,z)=\hfill & \frac{1}{64}[27{\psi }_{\kappa -1}(x+\frac{{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2})\hfill \\ \hfill & +9{\psi }_{\kappa -1}(x+\frac{{\eta }_{\kappa }}{2},y-\frac{3{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2})\hfill \\ \hfill & +9{\psi }_{\kappa -1}(x+\frac{{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2})\hfill \\ \hfill & +9{\psi }_{\kappa -1}(x-\frac{3{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2})\hfill \\ \hfill & +3{\psi }_{\kappa -1}(x+\frac{{\eta }_{\kappa }}{2},y-\frac{3{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2})\hfill \\ \hfill & +3{\psi }_{\kappa -1}(x-\frac{3{\eta }_{\kappa }}{2},y-\frac{3{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2})\hfill \\ \hfill & +3{\psi }_{\kappa -1}(x-\frac{3{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2})\hfill \\ \hfill & +{\psi }_{\kappa -1}(x-\frac{3{\eta }_{\kappa }}{2},y-\frac{3{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2})]\hfill \end{array}& (\mathrm{3.1.13})\end{array}$

$\begin{array}{ll}\begin{array}{ll}{P}_{\kappa -1}^{\kappa }{\psi }_{\kappa -1}(x,y,z)=\hfill & \frac{1}{64}[27{\psi }_{\kappa -1}(x-\frac{{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2})\hfill \\ \hfill & +9{\psi }_{\kappa -1}(x-\frac{{\eta }_{\kappa }}{2},y+\frac{3{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2})\hfill \\ \hfill & +9{\psi }_{\kappa -1}(x-\frac{{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2})\hfill \\ \hfill & +9{\psi }_{\kappa -1}(x+\frac{3{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2})\hfill \\ \hfill & +3{\psi }_{\kappa -1}(x-\frac{{\eta }_{\kappa }}{2},y+\frac{3{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2})\hfill \\ \hfill & +3{\psi }_{\kappa -1}(x+\frac{3{\eta }_{\kappa }}{2},y+\frac{3{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2})\hfill \\ \hfill & +3{\psi }_{\kappa -1}(x+\frac{3{\eta }_{\kappa }}{2},y-\frac{{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2})\hfill \\ \hfill & +{\psi }_{\kappa -1}(x+\frac{3{\eta }_{\kappa }}{2},y+\frac{3{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2})]\hfill \end{array}& (\mathrm{3.1.14})\end{array}$

$\begin{array}{ll}\begin{array}{ll}{P}_{\kappa -1}^{\kappa }{\psi }_{\kappa -1}(x,y,z)=\hfill & \frac{1}{64}[27{\psi }_{\kappa -1}(x-\frac{{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2})\hfill \\ \hfill & +9{\psi }_{\kappa -1}(x-\frac{{\eta }_{\kappa }}{2},y-\frac{3{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2})\hfill \\ \hfill & +9{\psi }_{\kappa -1}(x-\frac{{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2})\hfill \\ \hfill & +9{\psi }_{\kappa -1}(x+\frac{3{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2})\hfill \\ \hfill & +3{\psi }_{\kappa -1}(x-\frac{{\eta }_{\kappa }}{2},y-\frac{3{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2})\hfill \\ \hfill & +3{\psi }_{\kappa -1}(x+\frac{3{\eta }_{\kappa }}{2},y-\frac{3{\eta }_{\kappa }}{2},z\pm \frac{{\eta }_{\kappa }}{2})\hfill \\ \hfill & +3{\psi }_{\kappa -1}(x+\frac{3{\eta }_{\kappa }}{2},y+\frac{{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2})\hfill \\ \hfill & +{\psi }_{\kappa -1}(x+\frac{3{\eta }_{\kappa }}{2},y-\frac{3{\eta }_{\kappa }}{2},z\mp \frac{3{\eta }_{\kappa }}{2})]\hfill \end{array}& (\mathrm{3.1.15})\end{array}$

for points a and c, b and d, e and g, f and h, respectively. The coarse grid operator L_κ−1 is generated in a way that is analogous to the computation of L_κ on the fine grid. Alternatively, the Galerkin coarse grid operator can be used:

$\begin{array}{ll}\begin{array}{l}{L}_{\kappa -1}=R_{\kappa }^{\kappa -1}{L}_{\kappa }{P}_{\kappa -1}^{\kappa }\end{array}& (\mathrm{3.1.16})\end{array}$

where $R_{\kappa }^{\kappa -1}$ and $P_{\kappa -1}^{\kappa }$ are appropriately chosen transfer operators (Trottenberg et al., 2001).

FIGURE 4

Figure 4. Arrangement of unknowns in 2D cell-centered discretization.

FIGURE 5

Figure 5. Arrangement of unknowns in 3D cell-centered discretization.

A FAS V-cycle consists of two iterating components. The outer time iteration travels through one V–loop, starting from the finest mesh level, looking for the fixed point solution of $\text{L}({\psi }_{i,j,k}^{a,r},{\psi }_{i,j,k}^{a,r-1})=\text{R}({\psi }_{i,j,k}^{a-1},{\psi }_{i,j,k}^{a,r-1})$ as derived in Supplement (Multigrid V–Cycle Iterations). Starting from ψ^{a, 0}, each successive V–cycle produces an approximation that is converging toward the fixed point solution, ψ^a, at the current time step a. After r* outer time iterations (after r* V–loops), if the approximation ψ^{a, r}* results in an error within a tolerable bound, an approximated solution for the current time step is established by ψ^a = ψ^{a, r}*. The inner iteration is the non-linear Gauss-Seidel relaxation as described in Supplement (Non-linear Gauss–Seidel Relaxations) and given in Supplementary Equation (1.4.39). This procedure occurs twice on each level. Relaxation first takes place during the error smoothing step while traveling down from the finest to the coarsest grid level, and again after coarse grid correction while going up from the coarsest to the finest grid level. Hence, the starting approximation for the error smoothing step is obtained from the previous iteration ${\psi }_{\kappa }^{a,r,0}={\psi }_{\kappa }^{a,r-1}$ or restricted from the relaxed fine-grid approximation ${\psi }_{\kappa }^{a,r,0}=R_{\kappa }^{\kappa -1}{\psi }_{\kappa +1}^{a,r,{\nu }_1}$, whereas the coarse-grid-corrected approximation ${\psi }_{\kappa }^{a,r,0}={\psi }_{\kappa }^{a,r,CGC}$ is used during the correction step.

The computation of a new time iterate ${\psi }_{\kappa }^{a,r}$ using an adaptive FAS multigrid-cycle, starting with any of the refined or global grid levels described in the Multilevel Mesh Refinement subsection, can be done recursively using a two–level model as summarized in the operator below (Brandt, 1977; Trottenberg et al., 2001; Wise et al., 2007).

Recursive adaptive cycle operator:

$\begin{array}{ll}\begin{array}{l}{\psi }_{\kappa }^{a,r}=\mathrm{\text{ADAPFAS}}(\kappa ,\gamma ,{\nu }_0,{\nu }_1,{\nu }_2,{\psi }_{\kappa }^{a,r-1},{\psi }_{\kappa -1}^{a,r-1},L_{\kappa },R_{\kappa })\end{array}& (\mathrm{3.1.17})\end{array}$

I. Pre-smoothing – Compute a smoothed approximation ${\psi }_{\kappa }^{a,r,{\nu }_1}$ by applying ν₁ smoothing steps to ${\psi }_{\kappa }^{a,r-1}$ on Ω_κ:

$\begin{array}{ll}\begin{array}{l}{\psi }_{\kappa }^{a,r,{\nu }_1}=\mathrm{\text{SMOOTH}}({\nu }_1,{\psi }_{\kappa }^{a,r-1},L_{\kappa },R_{\kappa })\end{array}& (\mathrm{3.1.18})\end{array}$

II. Coarse-grid correction: – Initialize the coarse-grid iterate:

$\begin{array}{ll}{\psi }_{\kappa -1}^{a,r,0}=\{\begin{array}{cc}{R}_{\kappa }^{\kappa -1}{\psi }_{\kappa }^{a,r,{\nu }_1}& on{\Omega }_{\kappa -1}\cap {\Omega }_{\kappa }\\ {\psi }_{\kappa -1}^{a,r-1}& on{\Omega }_{\kappa -1}-{\Omega }_{\kappa }\end{array}& (\mathrm{3.1.19})\end{array}$

– Update the ghost cells on κ − 1 level using interpolation and exchange for ${\psi }_{\kappa -1}^{a,r,0}$ on neighboring patches.

– Compute the coarse-grid RHS:

$\begin{array}{ll}{R}_{\kappa -1}=\{\begin{array}{ll}{R}_{\kappa }^{\kappa -1}\left[R_{\kappa }-L_{\kappa }\left({\psi }_{\kappa }^{a,r,{\nu }_1}\right)\right]\hfill & \hfill \\ +L_{\kappa -1}\left(R_{\kappa }^{\kappa -1}{\psi }_{\kappa }^{a,r,{\nu }_1}\right)\hfill & on{\Omega }_{\kappa -1}\cap {\Omega }_{\kappa }\hfill \\ R_{\kappa -1}\hfill & on{\Omega }_{\kappa -1}-{\Omega }_{\kappa }\hfill \end{array}& (\mathrm{3.1.20})\end{array}$

– Compute an approximate solution ${\overline{\psi }}_{\kappa -1}^{a,r}$ of the following coarse-grid equation on Ω_κ−1:

$\begin{array}{ll}\begin{array}{l}{L}_{\kappa -1}({\overline{\psi }}_{\kappa -1}^{a,r})=R_{\kappa -1}\end{array}& (\mathrm{3.1.21})\end{array}$

If κ = κ_min+1, employ a direct solver or perform ν₀ smoothing steps:

$\begin{array}{ll}\begin{array}{l}{\overline{\psi }}_{{\kappa }_{min}}^{a,r}=\mathrm{\text{SMOOTH}}({\nu }_0,{\psi }_{{\kappa }_{min}}^{a,r,0},L_{{\kappa }_{min}},R_{{\kappa }_{min}})\end{array}& (\mathrm{3.1.22})\end{array}$

If κ > κ_min + 1, solve Equation (3.1.21) by employing γ adaptive FAS cycle using ${\psi }_{\kappa -1}^{a,r,0}$ as initial approximation:

$\begin{array}{ll}\begin{array}{l}{\overline{\mathbf{\psi }}}_{\kappa -1}^{a,r}=\text{ADAPFAS}(\kappa -1,\gamma ,{\nu }_0,{\nu }_1,{\nu }_2,{\mathbf{\psi }}_{\kappa -1}^{a,r,0}\\ {\mathbf{\psi }}_{\kappa -2}^{a,r-1}{\mathit{L}}_{\kappa -1},{\mathit{R}}_{\kappa -1})\end{array}& (\mathrm{3.1.23})\end{array}$

– Arrive at approximation for time iteration r on Ω_κ−1−Ω_κ:

$\begin{array}{ll}\begin{array}{l}{\psi }_{\kappa -1}^{a,r}={\overline{\psi }}_{\kappa -1}^{a,r}\end{array}& (\mathrm{3.1.24})\end{array}$

– Compute coarse-grid correction on Ω_κ−1∩Ω_κ:

$\begin{array}{ll}\begin{array}{l}{\text{e}}_{\kappa -1}^{a,r}={\overline{\psi }}_{\kappa -1}^{a,r}-{\psi }_{\kappa -1}^{a,r,0}\end{array}& (\mathrm{3.1.25})\end{array}$

– Interpolate the correction and compute the coarse-grid corrected approximation on Ω_κ:

$\begin{array}{ll}\begin{array}{l}{\psi }_{\kappa }^{a,r,CGC}={\psi }_{\kappa }^{a,r,{\nu }_1}+P_{\kappa -1}^{\kappa }{\text{e}}_{\kappa -1}^{a,r}\end{array}& (\mathrm{3.1.26})\end{array}$

– Update the ghost cells on κ level using interpolation and exchange for ${\psi }_{\kappa }^{a,r,CGC}$ on neighboring patches.

III. Post-smoothing:

Compute ${\psi }_{\kappa }^{a,r}$ by applying ν₂ smoothing steps to ${\psi }_{\kappa }^{a,r,CGC}$ on Ω_κ:

$\begin{array}{ll}\begin{array}{l}{\psi }_{\kappa }^{a,r}=\mathrm{\text{SMOOTH}}({\nu }_2,{\psi }_{\kappa }^{a,r,CGC},L_{\kappa },R_{\kappa })\end{array}& (\mathrm{3.1.27})\end{array}$

– Update the ghost cells on κ level using interpolation and exchange for ${\psi }_{\kappa }^{a,r}$ on neighboring patches.

V– or W–cycle is defined via the cycle index γ = 1 and γ = 2, respectively. The initial value used in V-cycles may come from the converged solution at previous time ψ^{a, 0} = ψ^a−1 or from coarse-grid approximations as in FMG method described in the next section. The SMOOTH routine in Equations (3.1.18), (3.1.22), and (3.1.27) uses Red-Black (or odd-even) ordering for the relaxation process, with ν_b extra smoothing steps for d_κ cells within the boundaries.

Adaptive FMG Structure

In FMG, a better initial approximation for finer-grid iteration is obtained from a coarser-grid level (Kronsjö and Dahlquist, 1972; Brandt, 1977). The procedure begins with a set of solutions on the coarsest grid, followed by interpolation of the solution set to a fine-grid level, providing an initial guess for the fine-grid multigrid approximation. After a few multigrid cycles on the fine-grid level, the solution is again interpolated to the next fine-grid level. This proceeds until the finest grid level has been reached, if refinement is needed for each level κ > 0. Discretization accuracy can be reached on each grid level within a few multigrid cycles in this manner. The structure of full multigrid is illustrated in Figure 6.

FIGURE 6

Figure 6. FMG on self-adaptive grids. Beginning on a series of global coarse grids, finer grid is introduced a layer per multigrid cycle. ⊗ represents solution on the coarsest grid; ⊗ represents smoothing on global grid levels; ∙ represents smoothing on refined grid levels; indicates FMG interpolation; ↘ indicates restriction; and ↗ indicates prolongation.

A higher order scheme is normally employed for an FMG interpolation. Cubic Lagrangian interpolation is one such method used in the FMG interpolation steps to compute an approximation at the fine grid. In 1D as depicted in Figure 7, approximation ψ_κ(x) at the fine grid point can be computed as

$\begin{array}{ll}\begin{array}{l}{\prod }_{\kappa -1}^{\kappa }{\mathbf{\psi }}_{\kappa -1}(x)=\frac{1}{128}[-7{\mathbf{\psi }}_{\kappa -1}(x-\frac{5{\eta }_{\kappa }}{2})\hfill \\ +105{\mathbf{\psi }}_{\kappa -1}(x-\frac{{\eta }_{\kappa }}{2})\hfill \\ +35{\mathbf{\psi }}_{\kappa -1}(x+\frac{3{\eta }_{\kappa }}{2})\hfill \\ -5{\mathbf{\psi }}_{\kappa -1}(x+\frac{7{\eta }_{\kappa }}{2})]\hfill \end{array}& (\mathrm{3.1.28})\end{array}$

$\begin{array}{ll}\begin{array}{l}{\prod }_{\kappa -1}^{\kappa }{\mathbf{\psi }}_{\kappa -1}(x)=\frac{1}{128}[-7{\mathbf{\psi }}_{\kappa -1}(x-\frac{5{\eta }_{\kappa }}{2})\hfill \\ +105{\mathbf{\psi }}_{\kappa -1}(x-\frac{{\eta }_{\kappa }}{2})\hfill \\ +35{\mathbf{\psi }}_{\kappa -1}(x-\frac{3{\eta }_{\kappa }}{2})\hfill \\ -5{\mathbf{\psi }}_{\kappa -1}(x+\frac{7{\eta }_{\kappa }}{2})]\hfill \end{array}& (\mathrm{3.1.29})\end{array}$

for point a and b respectively. Cubic interpolation in the y- and z-directions for bicubic and tricubic interpolation, in 2D and 3D cases, respectively, are analogous to the 1D calculation shown.

FIGURE 7

Figure 7. Arrangement of unknowns in 1D cell-centered discretization.

FMG can be performed on locally refined grids which are defined a priori. Here, we use self-adaptive FMG, where refined grids may or may not be generated, one layer at a time, based on some predefined criteria. Self-adaptive FMG advances as follows:

(1) Perform FMG on a series of global grids until a satisfactory accuracy has been obtained on grid level κ = 0.

(2) Determine if refining the finest grid would be beneficial, based on the criteria outlined in the Multilevel Mesh Refinement subsection. If no, iteration stops. If yes, refine grid in local regions.

(3) Use cubic interpolation on the current fine-grid (κ) approximation to compute an initial guess for the refined grid level (κ + 1), proceed by performing a number of multigrid cycles on the refined level until the required accuracy has been satisfied.

(4) If the finest level κ_max has been reached, the time iteration stops. Otherwise, return to Step (2).

Let tol_κ, γ_κ, and r_κ be the residual tolerance, cycle index, and number of multigrid cycles performed on level κ. The self-adaptive FMG-FAS iteration is given by the following:

Initialize ${\psi }_{\kappa }^{a=0}$, κ = κ_min, …, κ_max

Set γ_κ, r_κ, and tol_κ, κ = κ_min, …, κ_max

   For a = 1, a_max

      Set ${\psi }_{\kappa }^{a,r=0}={\psi }_{\kappa }^{a-1},\kappa ={\kappa }_{min},\dots ,{\kappa }_{max}$

$\begin{array}{ll}{\mathit{\psi }}_{\kappa }^a=\text{FMG}\left(\mathit{\kappa },\mathit{\gamma },{\mathit{\nu }}_{\mathit{0}},{\mathit{\nu }}_{\mathit{1}},{\mathit{\nu }}_{\mathit{2}},{\text{tol}}_{\text{k}}\right),\mathit{\kappa }={\mathit{\kappa }}_{\mathbf{min}},\cdots ,{\mathit{\kappa }}_{\mathbf{max}}& (\mathrm{3.1.30})\end{array}$

   End For

The recursive FMG routine from Rude (1993) is modified for Equation (3.1.30) above and shown as the following:

 For κ = κ_min, κ_max

If κ = κ_min, employ a direct solver or perform ν₀ smoothing steps at the coarsest level and return the solution:

    ${\overline{\psi }}_{{\kappa }_{min}}^{a,0,{\nu }_0}=\mathrm{\text{SMOOTH}}({\nu }_0,{\psi }_{{\kappa }_{min}}^{a,r=0},L_{{\kappa }_{min}},R_{{\kappa }_{min}})$

Else

If κ < 0, perform r_κ multigrid V-cycles at level κ:

r = r + 1

${\psi }_{\kappa }^{a,r}=\mathrm{\text{ADAPFAS}}(\kappa ,{\gamma }_{\kappa },to{l}_{\kappa },{\nu }_0,{\nu }_1,{\nu }_2,{\psi }_{\kappa }^{a,r-1},{\psi }_{\kappa -1}^{a,r-1},$ L_κ, R_κ)

Else perform multigrid V-cycles at level κ until the approximation yields an error within bound:

Do

r = r + 1

${\psi }_{\kappa }^{a,r}=\mathrm{\text{ADAPFAS}}(\kappa ,{\gamma }_{\kappa },to{l}_{\kappa },{\nu }_0,{\nu }_1,{\nu }_2,{\psi }_{\kappa }^{a,r-1},{\psi }_{\kappa -1}^{a,r-1},$ L_κ, R_κ)

While $\left|\right|R_{\kappa }-L_{\kappa }({\psi }_{\kappa }^{a,r})\left|\right|>to{l}_{\kappa }$

End If

If κ < 0, interpolate the coarse-grid solution to obtain fine-grid initial guess:

${\psi }_{\kappa +1}^{a,r-1}={\Pi }_{\kappa }^{\kappa +1}{\psi }_{\kappa }^{a,r-1}$

Else

If 0 ≤ κ < κ_max, flag cells on level κ to determine if refinement is needed:

$F_{\kappa }^{a,r-1}=\mathrm{\text{FLAG}}({\psi }_{\kappa }^{a,r-1})$

If any cells are flagged to be refined, generate blocks on level κ + 1 accordingly:

$\begin{array}{l}{B}_{\kappa +1}=\text{BLOCKGEN }(threshold\text{\_}eff\text{ , }threshold\text{\_}size\text{, }min\text{\_}eff\text{ ,}\\ F_{\kappa }^{a,r-1})\end{array}$

Then interpolate the coarse-grid solution to obtain fine-grid initial guess:

${\psi }_{\kappa +1}^{a,r-1}={\Pi }_{\kappa }^{\kappa +1}{\psi }_{\kappa }^{a,r-1}$

End If

   End For

A parameter σ is used to adapt the accuracy to the different level of meshes. Tolerance is reduced by σ for each finer mesh level, $to{l}_{\kappa }=\frac{to{l}_{\kappa -1}}{\sigma }$. The parameter σ and the initial value of tolerance have to be selected with care in order to ensure the overall accuracy of the final grid structure.

System Solution

Next, we applied the fully adaptive non-linear FMG algorithm to solve the desmoplastic tumor model. The set of tumor model related parameters, constants, source terms, rate expressions, and adjustment factors are listed in Supplementary Tables 1–6. We assume that the tissue consists of a fixed fraction of liquid interstitial-fluid phase (ϕ_β = 0.2) and solid cell-ECM phase (ϕ_α = 0.8). We also assume that the solid cell-ECM phase is only composed of viable tumor cells (${\tilde{\varphi }}_V$), dead tumor cells (${\tilde{\varphi }}_D$), extracellular matrix ECM (${\tilde{\varphi }}_E$), and healthy host cells (${\tilde{\varphi }}_H$). Nutrients, waste products, and tumorigenic species are assumed to be carried in the liquid phase, while the myofibroblastic cell species (${\tilde{F}}_E$), blood (${\tilde{B}}_n^E$) and lymphatic (${\tilde{L}}_n^E$) microvessels are assumed to take up negligible volume and exist within the continuous ECM component.

We assume that the source and sink of viable tumor cells comes from cell proliferation and losses to necrosis respectively. The dead tumor cell species hence comes from the necrosed viable tumor cells and undergoes lysis. The ECM is assumed to be secreted by the myofibroblastic cell species and degraded by the matrix degrading enzyme species ($\tilde{m}$). The healthy host cells species is assumed to be homeostatic. Nutrients such as oxygen (ñ) and glucose ($\tilde{g}$) enter the tissue via blood vessels and are consumed mainly by the viable tumor species. Tumor growth factors ($t\tilde{g}f$) and tumor angiogenic factors ($t\tilde{a}f$) are assumed to be secreted by viable tumor cells and can undergo degradation. Tumor angiogenic factors also have an uptake term by proliferating vessels. Myofibroblastic cell species is assumed to go through proliferation, apoptosis, and necrosis, while the blood and lymphatic vessels undergo restructuring and potential loss due to crushing by the surrounding tissue pressure. To slightly perturb the symmetry of tumor progression, we let the fraction of blood vessels that are sprouting be different in different regions of the domain, ranging from 0.1 to 0.8.

As described in Methods, all PDEs were discretized using the Crank-Nicolson method. To solve the set of cell-centered discretized equations, we applied full multigrid algorithm and adaptive full approximation scheme V-cycle with Gauss-Seidel Red-Black smoothing to the multiple grid level system. Here, we use a total of five grid levels as depicted in Figure 1, with three increasingly refined global levels (κ = −2, −1, 0) covering the entire domain Ω and two adaptively refined levels (κ = 1, 2). The domain is Ω = (0, 40)³ and the mesh sizes for the finest global level to the finest adaptively refined level are η₀ = 40/32, η₁ = 40/64, and η₂ = 40/128. The time step-size used is θ = 1 × 10⁻² and minimum tolerance is set at 5 × 10⁻⁴. Key solver related parameters are listed in Table 1.

TABLE 1

Table 1. Solver related parameters.

The center of the computational domain is simply seeded with viable tumor cells (${\tilde{\varphi }}_V=0.65$) at the beginning of the simulation. We assume that there are no dead tumor cells (${\tilde{\varphi }}_D=0$) initially and the ECM is distributed evenly (${\tilde{\varphi }}_E=0.35$) across the domain. The healthy host cells thus take up the remaining volume in the solid phase of the tissue. The initial tumor shape is therefore rectangular with a sharp interface. The simplified source term for the viable tumor cell species used in Equation (2.1), where only mitosis and necrosis are considered, is given by:

$\begin{array}{ll}\begin{array}{l}{\tilde{S}}_v=\{{\tilde{\lambda }}_{M,V}\tilde{n}\left(1+t\tilde{g}f\right){\rm H}\left(\tilde{n}-{\tilde{n}}_h\right)\\ -{\tilde{\lambda }}_{N,V}\left[1-{\rm H}\left(\tilde{n}-{\tilde{n}}_{v,V}\right){\rm H}\left(\tilde{g}-{\tilde{g}}_{v,V}\right)\right]\}{\tilde{\varphi }}_V,\end{array}& (5.1)\end{array}$

where ${H}$ is the Heaviside function, ñ_h = 0.3 is the hypoxic level of oxygen, ${\tilde{\lambda }}_{M,V}$ = 1 and ${\tilde{\lambda }}_{N,V}=3$ are the rate constants of mitosis and necrosis, respectively, for viable tumor cells. ñ_{v, V} = 0.21 and ${\tilde{g}}_{v,V}=0.1$ are the oxygen and glucose viability limits, respectively, for the viable tumor cell species. The mitosis rate is assumed to be upregulated by the level of tumor growth factors. The center of the tumor mass experiences a significant drop in nutrient levels. When the oxygen level drops below the hypoxic level, viable tumor cells cease to reproduce. Necrosis takes place if the oxygen level falls below the viability limit, where viable tumor cells necrosed to dead tumor cells. The source term for the dead tumor cell species in Equation (2.2) is hence given by

$\begin{array}{ll}\begin{array}{l}{\tilde{S}}_D={\tilde{\lambda }}_{N,V}\left[1-{H}(ñ-{ñ}_{v,V}){H}(\tilde{g}-{\tilde{g}}_{v,V})\right]{\tilde{\varphi }}_V-{\tilde{\lambda }}_{L,D}{\tilde{\varphi }}_D,\end{array}& (5.2)\end{array}$

where ${\tilde{\lambda }}_{L,D}=1$ is the lysis rate constant for the dead tumor cell species. The lysed dead tumor cells are assumed loss to the interstitial fluid. The other major solid component of the tissue is the ECM, and its source term used in Equation (2.3) is

$\begin{array}{ll}\begin{array}{l}{\tilde{S}}_E={\tilde{\lambda }}_{F,E}{{A}}_{F,E}\tilde{F}-{\tilde{\lambda }}_{de,E}\tilde{m}{\tilde{\varphi }}_E,\end{array}& (5.3)\end{array}$

where ${\tilde{\lambda }}_{F,E}=5$ is the secretion rate constant of ECM by the myofibroblastic cell species $\tilde{F}$ and ${\tilde{\lambda }}_{de,E}=5$ is the degradation rate constant of ECM by the matrix degrading enzyme species $\tilde{m}$. Factors affecting the rate constant of ECM secretion by the myofibroblastic cell species are included in the adjustment factor ${{A}}_{F,E}$ given by

$\begin{array}{ll}\begin{array}{l}{{A}}_{F,E}=\begin{array}{c}(1-{\tilde{\varphi }}_T-{\tilde{\varphi }}_E)(1+t\tilde{g}f)\left[1+{{F}}_{n,E}^F\frac{{ñ}_h-ñ}{{ñ}_h-{ñ}_{v,F}}{H}({ñ}_h-ñ)\right]\\ {H}(ñ-{ñ}_{v,F}){H}(t\tilde{g}f-{t\tilde{g}f}_{F,E}){H}(1-{\tilde{\varphi }}_T-{\tilde{\varphi }}_E)\end{array}\end{array}& (5.4)\end{array}$

where ${{F}}_{n,E}^F=2$ is the effective factor of hypoxia on upregulating the production of ECM by myofibroblastic cells and ñ_{v, F} = 0.21 is the oxygen viability limit of the myofibroblastic cell species. The constant ${t\tilde{g}f}_{F,E}=0.2$ is a lower threshold of tumor growth factors, below which the production of ECM by the myofibroblastic cell species is assumed negligible. A complete list of source terms, rate expressions, and adjustment factors is shown in Supplementary Table 6.

The transient tumor progression is shown in Figure 8, via tracking the evolution of the tumor isosurface ${\tilde{\varphi }}_T=0.15$. With the specific set of constants and parameters used in this case, the tumor undergoes regression in the time frame shown. Nutrients levels drop below the viability limit in the center of the tumor mass, causing dead tumor cells to accumulate from necrosed viable tumor cells, which are then lysed from the tumor mass. The increase in tumor growth factors in the tumor region upregulates the secretion of ECM by the myofibroblastic cell species, resulting in a high ECM environment within the tumor mass (as in Figure 9).

FIGURE 8

Figure 8. The transient ϕ_T = 0.15 isosurface at time = 10, 30, and 50. The total tumor volume fraction ϕ_T = ϕ_V + ϕ_D. The center of the domain is initially seeded with viable tumor cells ϕ_V = 0.65. The starting extracellular matrix volume fraction is set to be homogenously distributed across the domain at ϕ_E = 0.35. The remaining volume fraction is thus made up by healthy host cells, denoted by ϕ_H.

FIGURE 9

Figure 9. 1-D Profiles of species with trilinear (left side) vs. cubic (right side) interpolation used in FMG-interpolation at time = 10 (sliced at j = k = 58). First row: tumor viable species ${\tilde{\varphi }}_V$, dead species ${\tilde{\varphi }}_D$, and ECM species ${\tilde{\varphi }}_E$. The overall tumor pressure is labeled by p_α. Second row: diffusible substances driving the tumor evolution, including oxygen (O₂), glucose (Glu), carbon dioxide (CO₂), bicarbonate (Bic), lactate (Lac), and hydrogen ions (H⁺). Third row: Concentration of myfibroblasts (myF), tumor growth factors (TGF), and matrix degrading enzymes (MDE). Fourth row: corresponding density of blood vasculature (Bn), lymphatic vasculature (Ln), and tumor angiogenic factors (TAF). The myofibrobrastic cell species density in tissue (shown in the third row of plots as myF) is computed by $\tilde{F}={\tilde{F}}_E{\tilde{\varphi }}_E$. Similarly, the blood and lymphatic vessel densities in tissue (shown in the last row of plots as Bn and Ln) are given by ${\tilde{B}}_n={\tilde{B}}_n^E{\tilde{\varphi }}_E$ and ${\tilde{L}}_n={\tilde{L}}_n^E{\tilde{\varphi }}_E$ respectively.

As mentioned in Results (Adaptive FMG Structure), a higher order scheme is normally used in the FMG interpolation step when interpolating the solution to a finer grid level. Trilinear interpolation as shown in Equations (3.1.12–3.1.15) and cubic interpolation given in Equations (3.1.28) and (3.1.29) are used in the FMG interpolation step and the simulation results for time = 10 are compared in Figure 9. The ECM and viable tumor levels are slightly higher in the case where cubic interpolation is used, resulting in higher solid cell–ECM phase pressure. Since the myofibroblastic cell species and vessels are assumed to reside within the ECM component, the higher ECM volume fraction in the high order scheme case also results in higher myofibroblastic cell species (presented in green as myF) and vessel densities (presented in red as Bn for blood and blue as Ln for lymphatic) in tissue. Concentration profiles of interstitial fluid based species show no visible difference in the two cases.

The adaptive block-structured mesh system for the simulation in Figures 8, 9 at time = 50 is shown in Figure 10. For global level κ = 0, there is one block (shown in black) covering the entire domain with mesh size η₀ = 40/32. In the first level of refinement at level κ = 1, there are four blocks (shown in green) with mesh size η₁ = 40/64. In the next level of refinement at level κ = 2 and within the coarser domain Ω₁, there are six blocks (shown in red) with mesh size η₂ = 40/128. As shown in Table 1, the critical value of C_κ = 0.05 in Equation (3.1.1) is used in the undivided gradient test to flag cells for refinement. We also buffer each flag cell by flagging four cells surrounding it in each direction, creating a cube of 9 × 9 × 9 flagged cells, or less if it is near any external or ghost boundaries. The thresholdsize sets the minimum block size to 10³, resulting in overall bigger blocks generated by the BLOCKGEN routine in Equation (3.1.4).

FIGURE 10

Figure 10. ϕ_T = 0.15 isosurface at time = 50 shown with the blocks of the block-structured mesh in the global and refined levels. There is one block (black) covering the entire global domain on the global level = 0 with the mesh spacing η₀ = 40/32. There are two levels of refinement. Blocks covering the refined level = 1 (green) correspond to mesh spacing η₁ = 40/64, whereas blocks on the refined level = 2 (red) correspond to mesh spacing η₂ = 40/128.

In Figure 11, the transient degree of freedom (DOF) for the simulation in Figures 8, 9 is plotted. The degree of freedom is represented by the total cell-centered grid points in the 3D domain. The mesh with just the three global levels (κ = −2, −1, 0) has only 37,376 DOF. In the first iteration, both refinement levels (κ = 1, 2) begin to contain flagged cells and the mesh now has 224,640 DOF. At time = 50, there are 554,568 DOF, increased by a factor of 2 from the first iteration.

FIGURE 11

Figure 11. Degree of freedom, represented by total cell-centered grid points in the 3D domain, from time = 0 to 50. Here, external boundary and internal (ghost) boundary grid points are not accounted.

A 3D convergence test as described in Wise et al. (2007) was done by varying grid and time-step sizes. Three simulations were performed where grid spacings used for the root level κ = 0 are 16³, 32³, and 64³, respectively. We focused on the convergence of the tumor volume fraction, ${\tilde{\varphi }}_T$, and in all three cases, data were interpolated to the cell centers of their global uniform grids corresponding to their finest mesh η_κmax. With κ_max = 2, the grid sizes η₂ and the corresponding time step sizes θ, which follow the linear refinement path θ = 0.032η, are outlined in Table 2. The errors were calculated by comparing cell center data to the averaged cell center data of the finer grid set, and the l₂ norms of the errors were used to obtain the convergence rates as shown in Wise et al. (2007). The errors and the rate of convergence are also listed in Table 2. The results show that the algorithm presented herein is first-order accurate. Even though the discretization is second order in both time and space, the first-order rate of convergence attained, as also reported by Wise et al. (2011), is expected since the surface adhesion terms in Supplementary Equation (1.2.6) were treated explicitly. Non-smooth functions used in adjustment factors and first order interpolation function used in prolongation are other potential contributing factors to the first-order convergence.

TABLE 2

Table 2. Errors and convergence rate.

Discussion

This paper illustrates the application of a fully adaptive, non-linear full multigrid, finite-difference algorithm to solve a diffuse interface desmoplastic tumor system (Ng and Frieboes, 2017). The set of PDEs in the model is discretized in time using the Crank-Nicolson method. A Non-linear Full Approximation Scheme is used in the full multigrid V-cycle iterations, and Red-Black ordering is used in the Gauss-Seidel relaxation. A block-structure multilevel Cartesian mesh isused consisting of three global levels with mesh sizes η₋₂ = 40/8, η₋₁ = 40/16, η₀ = 40/32, and two adaptively refined levels with mesh sizes η₁ = 40/64, η₂ = 40/128. A numerical simulation of desmoplastic tumor progression in 3D which shows that the algorithm is capable of simulating an ECM-rich tumor environment and the handling of morphological changes during tumor progression (Ng and Frieboes, 2017).

The diffuse interface model is chosen to bypass the need to track the transient position of sharp interfaces, as well as enforcing complicated boundary conditions across these interfaces, such as those required in a sharp interface model. The adaptive multigrid algorithm used can efficiently handle the narrow transition layers as well as larger morphological evolutions. As opposed to lexicographic ordering, the Red-Black sweep used here in relaxations/error smoothings allows for easy parallelization to minimize the computational expense.

In the future, 3D meshes with more refinement levels and finer grids will be evaluated. Different criteria, such as the relative truncation error test or a simple volume fraction test, will be used to flag cells for refinement, and the minimum block size allowed during block generation will also be reduced. The mathematical model may be augmented to include additional tumor cell species as well as immune cell species, together with their corresponding functions and effects. Tumorigenic species such as hormonal growth factors and chemoattractants could be added. An anticancer drug species may be added to study various tumor responses to different therapeutics. A more expansive role and interaction of the lymphatic system with the microenvironment (e.g., Scianna et al., 2013; Swartz, 2014) could be implemented. The current numerical model is coded in C and is partially parallelized. Full parallelization of the program will be explored to speed up and extend the capabilities of the code to handle models of larger scale.

Author Contributions

HF: study conception and scientific oversight; CN and HF: model definition and design; CN: initial implementation; CN and HF: code debugging; CN and HF: manuscript writing and revision. Both authors reviewed and approved the final manuscript.

Funding

This work was supported in part by the National Institutes of Health/National Cancer Institute (grant numbers U54CA143907, R01CA180149, and R15CA203605), and by the Department of Bioengineering and the James Graham Brown Cancer Center at the University of Louisville.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer LP and handling Editor declared their shared affiliation.

Acknowledgments

This work was conducted in part using the resources in the research computing group and the Cardinal Research Cluster (CRC) at the University of Louisville. The authors are grateful to Harrison Simrall for computing assistance and to Dylan Goodin for proofreading the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys.2018.00821/full#supplementary-material

References

Ambrosi, D., and Preziosi, L. (2002). On the closure of mass balance models for tumor growth. Math. Mod. Methods Appl. Sci. 12, 737–754. doi: 10.1142/S0218202502001878

CrossRef Full Text | Google Scholar

Ambrosi, D., and Preziosi, L. (2009). Cell adhesion mechanisms and stress relaxation in the mechanics of tumours. Biomech. Model. Mechanobiol. 8, 397–413. doi: 10.1007/s10237-008-0145-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Andasari, V., Gerisch, A., Lolas, G., South, A. P., and Chaplain, M. A. J. (2011). Mathematical modeling of cancer cell invasion of tissue: biological insight from mathematical analysis and computational simulation. J. Math. Biol. 63, 141–171. doi: 10.1007/s00285-010-0369-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, A. R. A., and Chaplain, M. A. J. (1998). Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull. Math. Biol. 60, 857–899. doi: 10.1006/bulm.1998.0042

PubMed Abstract | CrossRef Full Text | Google Scholar

Araujo, R. P., and McElwain, D. L. S. (2004). A history of the study of solid tumour growth: the contribution of mathematical modelling. Bull. Math. Biol. 66, 1039–1091. doi: 10.1016/j.bulm.2003.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Araujo, R. P., and McElwain, D. L. S. (2005a). A mixture theory for the genesis of residual stresses in growing tissues I: a general formulation. SIAM J. Appl. Math. 65, 1261–1284. doi: 10.1137/040607113

CrossRef Full Text | Google Scholar

Araujo, R. P., and McElwain, D. L. S. (2005b). A mixture theory for the genesis of residual streses in growing tissues II: Solutions to the biphasic equations for a multicell spheroid. SIAM J. Appl. Math. 66, 447–467. doi: 10.1137/040607125

CrossRef Full Text | Google Scholar

Arduino, A., and Preziosi, L. (2015). A multiphase model of tumour segregation in situ by a heterogeneous extracellular matrix. Int. J. Non Linear Mech. 75, 22–30. doi: 10.1016/j.ijnonlinmec.2015.04.007

CrossRef Full Text | Google Scholar

Astanin, S., and Preziosi, L. (2007). “Multiphase models of tumour growth,” in Selected Topics on Cancer Modelling: Genesis - Evolution - Immune Competition - Therapy, eds N. Bellomo, M. Chaplain, and E. DeAngelis (Boston, MA: Birkhäuser), 1–31.

Google Scholar

Bachmann, J., Raue, A., Schilling, M., Becker, V., Timmer, J., and Klingmüller, U. (2012). Predictive mathematical models of cancer signalling pathways. J. Intern. Med. 271, 155–165. doi: 10.1111/j.1365-2796.2011.02492.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bearer, E. L., Lowengrub, J. S., Frieboes, H. B., Chuang, Y.-L., Jin, F., Wise, S. M., et al. (2009). Multiparameter computational modeling of tumor invasion. Cancer Res. 69, 4493–4501. doi: 10.1158/0008-5472.CAN-08-3834

PubMed Abstract | CrossRef Full Text | Google Scholar

Bell, J., Berger, M., Saltzman, J., and Welcome, M. (1994). 3-Dimensional adaptive mesh refinement for hyperbolic conservation-laws. SIAM J. Sci. Comp. 15, 127–138. doi: 10.1137/0915008

CrossRef Full Text | Google Scholar

Berger, M., and Rigoutsos, I. (1991). An algorithm for point clustering and grid generation. Ieee Trans. Syst. Man Cybernet. 21, 1278–1286. doi: 10.1109/21.120081

CrossRef Full Text | Google Scholar

Brandt, A. (1977). Multi-Level Adaptive Solutions to Boundary-Value Problems. Math. Comp. 31, 333–390. doi: 10.1090/S0025-5718-1977-0431719-X

CrossRef Full Text | Google Scholar

Bresch, D., Colin, T., Grenier, E., Ribba, B., and Saut, O. (2010). Computational modeling of solid tumor growth: the avascular stage. SIAM J. Sci. Compt. 32, 2321–2344. doi: 10.1137/070708895

CrossRef Full Text | Google Scholar

Breward, C. J. W., Byrne, H. M., and Lewis, C. E. (2002). The role of cell-cell interactions in a two-phase model for avascular tumour growth. J. Math. Biol. 45, 125–152. doi: 10.1007/s002850200149

PubMed Abstract | CrossRef Full Text | Google Scholar

Breward, C. J. W., Byrne, H. M., and Lewis, C. E. (2003). A multiphase model describing vascular tumour growth. Bull. Math. Biol. 65, 609–640. doi: 10.1016/S0092-8240(03)00027-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Byrne, H. M. (2010). Dissecting cancer through mathematics: from the cell to the animal model. Nat. Rev. Cancer 10, 221–230. doi: 10.1038/nrc2808

PubMed Abstract | CrossRef Full Text | Google Scholar

Byrne, H. M., Alarcon, T., Owen, M. R., Webb, S. D., and Maini, P. K. (2006). Modelling aspects of cancer dynamics: a review. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 364, 1563–1578. doi: 10.1098/rsta.2006.1786

PubMed Abstract | CrossRef Full Text | Google Scholar

Byrne, H. M., King, J. R., McElwain, D. L. S., and Preziosi, L. (2003). A two-phase model of solid tumour growth. Appl. Math. Lett. 16, 567–573. doi: 10.1016/S0893-9659(03)00038-7

CrossRef Full Text | Google Scholar

Byrne, H., and Preziosi, L. (2003). Modelling solid tumour growth using the theory of mixtures. Math. Med. Biol. 20, 341–366. doi: 10.1093/imammb/20.4.341

PubMed Abstract | CrossRef Full Text | Google Scholar

Cahn, J. W. (1959). Free energy of a nonuniform system. II. Thermodynamic basis. J. Chem. Phys. 30, 1121–1124. doi: 10.1063/1.1730145

CrossRef Full Text | Google Scholar

Cahn, J. W., and Hilliard, J. E. (1958). Free energy of a nonuniform system. I. Interfacial free energy. J. Chem. Phys. 28, 258–267. doi: 10.1063/1.1744102

CrossRef Full Text | Google Scholar

Casciari, J. J., Sotirchos, S. V., and Sutherland, R. M. (1992). Mathematical-modeling of microenvironment and growth in emt6/ro multicellular tumor spheroids. Cell Prolif. 25, 1–22. doi: 10.1111/j.1365-2184.1992.tb01433.x

CrossRef Full Text | Google Scholar

Chaplain, M. A. J. (1996). Avascular growth, angiogenesis and vascular growth in solid tumours: The mathematical modelling of the stages of tumour development. Math. Comp. Model. 23, 47–87. doi: 10.1016/0895-7177(96)00019-2

CrossRef Full Text | Google Scholar

Chaplain, M. A. J. (2011). Multiscale mathematical modelling in biology and medicine. IMA J. Appl. Math. 76, 371–388. doi: 10.1093/imamat/hxr025

CrossRef Full Text | Google Scholar

Chatelain Clément, C., Ciarletta, P., and Ben Amar, M. (2011). Morphological changes in early melanoma development: influence of nutrients, growth inhibitors and cell-adhesion mechanisms. J. Theor. Biol. 290, 46–59. doi: 10.1016/j.jtbi.2011.08.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., and Lowengrub, J. S. (2014). Tumor growth in complex, evolving microenvironmental geometries: a diffuse domain approach. J. Theor. Biol. 361, 14–30. doi: 10.1016/j.jtbi.2014.06.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Wise, S. M., Shenoy, V. B., and Lowengrub, J. S. (2014). A stable scheme for a nonlinear, multiphase tumor growth model with an elastic membrane. Int. J. Numer. Method. Biomed. Eng. 30, 726–754. doi: 10.1002/cnm.2624

PubMed Abstract | CrossRef Full Text | Google Scholar

Colella, P., Graves, D. T., Keen, N. D., Ligocki, T. J., Martin, D. F., McCorquodale, P. W., et al. (2009). “Chombo software package for AMR Applications Desgin Document,” in Applied Numerical Algorithms Group CRD, Berkeley, CA.

Google Scholar

Cristini, V., and Lowengrub, J. (2010). Multiscale Modeling of Cancer: An Integrated Experimental and Mathematical Modeling Approach. Cambridge: Cambridge University Press.

Google Scholar

Deisboeck, T. S., Wang, Z. H., Macklin, P., and Cristini, V. (2011). “Multiscale cancer modeling,” in Annual Review of Biomedical Engineering, Vol 13. Annual Review of Biomedical Engineering, eds M. L. Yarmush, J. S. Duncan, and M. L. Gray (Palo Alto, CA: Annual Reviews), 127–55.

Google Scholar

de Visser, K. E., and Coussens, L. M. (2006). The inflammatory tumor microenvironment and its impact on cancer development. Contrib. Microbiol. 13, 118–137. doi: 10.1159/000092969

PubMed Abstract | CrossRef Full Text | Google Scholar

Edelman, L. B., Eddy, J. A., and Price, N. D. (2010). In silico models of cancer. Wiley Interdiscip Rev Syst. Biol. 2, 438–459. doi: 10.1002/wsbm.75

PubMed Abstract | CrossRef Full Text | Google Scholar

Escher, J., and Matioc, A.-V. (2013). Analysis of a two-phase model describing the growth of solid tumors. Eur. J. Appl. Math. 24, 25–48. doi: 10.1017/S0956792512000290

CrossRef Full Text | Google Scholar

Frieboes, H. B., Chaplain, M. A. J., Thompson, A. M., Bearer, E. L., and Lowengrub, J. S. (2011). Physical oncology: a bench-to-bedside quantitative and predictive approach. Cancer Res. 71, 298–302. doi: 10.1158/0008-5472.CAN-10-2676

PubMed Abstract | CrossRef Full Text | Google Scholar

Frieboes, H. B., Jin, F., Chuang, Y. L., Wise, S. M., Lowengrub, J. S., and Cristini, V. (2010). Three-dimensional multispecies nonlinear tumor growth-II: tumor invasion and angiogenesis. J. Theor. Biol. 264, 1254–1278. doi: 10.1016/j.jtbi.2010.02.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Frieboes, H. B., Lowengrub, J. S., Wise, S., Zheng, X., Macklin, P., Bearer, E. L., et al. (2007). Computer simulation of glioma growth and morphology. NeuroImage 37(Suppl. 1), S59–S70. doi: 10.1016/j.neuroimage.2007.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Frieboes, H. B., Smith, B. R., Chuang, Y.-L., Ito, K., Roettgers, A. M., Gambhir, S. S., et al. (2013). An integrated computational/experimental model of lymphoma growth. PLoS Comp. Biol. 9:e1003008. doi: 10.1371/journal.pcbi.1003008

PubMed Abstract | CrossRef Full Text | Google Scholar

Galle, J., Preziosi, L., and Tosin, A. (2009). Contact inhibition of growth described using a multiphase model and an individual cell based model. Appl. Math. Lett. 22, 1483–1490. doi: 10.1016/j.aml.2008.06.051

CrossRef Full Text | Google Scholar

Garcke, H. (2005). On a Cahn-Hilliard model for phase separation with elastic misfit. Ann. Institut. Henri Poincare Anal. Non Lineaire 22, 165–185. doi: 10.1016/j.anihpc.2004.07.001

CrossRef Full Text | Google Scholar

Gerisch, A., and Chaplain, M. A. J. (2008). Mathematical modelling of cancer cell invasion of tissue: Local and non-local models and the effect of adhesion. J. Theor. Biol. 250, 684–704. doi: 10.1016/j.jtbi.2007.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Graziano, L., and Preziosi, L. (2007). “Mechanics in tumor growth,” in Modeling of Biological Materials, eds F. Mollica, L. Preziosi, and K. R. Rajagopal (Basel: Birkhäuser), 263–321.

Google Scholar

Hanahan, D., and Weinberg, R. A. (2000). The hallmarks of cancer. Cell 100, 57–70. doi: 10.1016/S0092-8674(00)81683-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Hatzikirou, H., Deutsch, A., Schaller, C., Simon, M., and Swanson, K. (2005). Mathematical modelling of glioblastoma tumour development: a review. Math. Models Meth. Appl. Sci. 15, 1779–1794. doi: 10.1142/S0218202505000960

CrossRef Full Text | Google Scholar

Hawkins-Daarud, A., van der Zee, K. G., and Oden, J. T. (2012). Numerical simulation of a thermodynamically consistent four-species tumor growth model. Int. J. Numer. Method. Biomed. Eng. 28, 3–24. doi: 10.1002/cnm.1467

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J., and Lowengrub, J. (2005). Phase field modeling and simulation of three-phase flows. Interface Free Bound. 7, 435–466. doi: 10.4171/IFB/132

CrossRef Full Text | Google Scholar

Klika, V. (2014). A guide through available mixture theories for applications in critical reviews in solid state and material science. Critic. Rev. Solid State Materials Sci. 39, 154–174. doi: 10.1080/10408436.2012.719132

CrossRef Full Text | Google Scholar

Kreeger, P. K., and Lauffenburger, D. A. (2010). Cancer systems biology: a network modeling perspective. Carcinogenesis 31, 2–8. doi: 10.1093/carcin/bgp261

PubMed Abstract | CrossRef Full Text | Google Scholar

Kronsjö, L., and Dahlquist, G. (1972). On the design of nested iterations for elliptic difference equations. BIT 12, 63–71. doi: 10.1007/BF01932674

CrossRef Full Text | Google Scholar

Kuusela, E., and Alt, W. (2009). Continuum model of cell adhesion and migration. J. Math. Biol. 58, 135–161. doi: 10.1007/s00285-008-0179-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Leo, P. H., Lowengrub, J. S., and Jou, H. J. (1998). A diffuse interface model for microstructural evolution in elastically stressed solids. Acta Mater. 46, 2113–2130. doi: 10.1016/S1359-6454(97)00377-7

CrossRef Full Text | Google Scholar

Levine, H. A., Pamuk, S., Sleeman, B. D., and Nilsen-Hamilton, M. (2001a). Mathematical modeling of capillary formation and development in tumor angiogenesis: penetration into the stroma. Bull. Math. Biol. 63, 801–863. doi: 10.1006/bulm.2001.0240

PubMed Abstract | CrossRef Full Text | Google Scholar

Levine, H. A., Sleeman, B. D., and Nilsen-Hamilton, M. (2000). A mathematical model for the roles of pericytes and macrophages in the initiation of angiogenesis. I. The role of protease inhibitors in preventing angiogenesis. Math. Biosci. 168, 77–115. doi: 10.1016/S0025-5564(00)00034-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Levine, H. A., Sleeman, B. D., and Nilsen-Hamilton, M. (2001b). Mathematical modeling of the onset of capillary formation initiating angiogenesis. J. Math. Biol. 42, 195–238. doi: 10.1007/s002850000037

PubMed Abstract | CrossRef Full Text | Google Scholar

Lowengrub, J. S., Frieboes, H. B., Jin, F., Chuang, Y. L., Li, X., Macklin, P., et al. (2010). Nonlinear modelling of cancer: bridging the gap between cells and tumours. Nonlinearity 23, R1–R91. doi: 10.1088/0951-7715/23/1/R01

PubMed Abstract | CrossRef Full Text | Google Scholar

Mantzaris, N. V., Webb, S., and Othmer, H. G. (2004). Mathematical modeling of tumor-induced angiogenesis. J. Math. Biol. 49, 111–187. doi: 10.1007/s00285-003-0262-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Michor, F., Liphardt, J., Ferrari, M., and Widom, J. (2011). What does physics have to do with cancer? Nat. Rev. Cancer 11, 657–670. doi: 10.1038/nrc3092

PubMed Abstract | CrossRef Full Text | Google Scholar

Ng, C. F., and Frieboes, H. B. (2017). Model of vascular desmoplastic multispecies tumor growth. J. Theor. Biol. 430, 245–282. doi: 10.1016/j.jtbi.2017.05.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Oden, J. T., Hawkins, A., and Prudhomme, S. (2010). General diffuse-interface theories and an approach to predictive tumor growth modeling. Math. Models Meth. Appl. Sci. 20, 477–517. doi: 10.1142/S0218202510004313

CrossRef Full Text | Google Scholar

Oden, J. T., Lima, E. A. B. F., Almeida, R. C., Feng, Y., Rylander, M. N., Fuentes, D., et al. (2015). Toward predictive multiscale modeling of vascular tumor growth: computational and experimental oncology for tumor prediction. Arch. Comp. Methods Engineer. 23, 735–779. doi: 10.1007/s11831-015-9156-x

CrossRef Full Text

Osborne, J. M., Walter, A., Kershaw, S. K., Mirams, G. R., Fletcher, A. G., Pathmanathan, P., et al. (2010). A hybrid approach to multi-scale modelling of cancer. Philos. Trans. 368, 5013–5028. doi: 10.1098/rsta.2010.0173

PubMed Abstract | CrossRef Full Text | Google Scholar

Perez-Moreno, M. (2009). When neighbourhood matters: tumour microenvironment. Clin. Transl. Oncol. 11, 70–74. doi: 10.1007/s12094-009-0316-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Preziosi, L., Ambrosi, D., and Verdier, C. (2010). An elasto-visco-plastic model of cell aggregates. J. Theor. Biol. 262, 35–47. doi: 10.1016/j.jtbi.2009.08.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Preziosi, L., and Tosin, A. (2009a). Multiphase and multiscale trends in cancer modelling. Math. Model. Nat. Phenom. 4, 1–11. doi: 10.1051/mmnp/20094301

CrossRef Full Text | Google Scholar

Preziosi, L., and Tosin, A. (2009b). Multiphase modelling of tumour growth and extracellular matrix interaction: mathematical tools and applications. J. Math. Biol. 58, 4–5. doi: 10.1007/s00285-008-0218-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Preziosi, L., and Vitale, G. A. (2011). Multiphase model of tumor and tissue growth including cell adhesion, and plastic reorganization. Math. Models Meth. Appl. Sci. 21, 1901–1932. doi: 10.1142/S0218202511005593

CrossRef Full Text | Google Scholar

Psiuk-Maksymowicz, K. (2013). Multiphase modelling of desmoplastic tumour growth. J. Theor. Biol. 329, 52–63. doi: 10.1016/j.jtbi.2013.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Quaranta, V., Weaver, A. M., Cummings, P. T., and Anderson, A. R. A. (2005). Mathematical modeling of cancer: The future of prognosis and treatment. Clin. Chim. Acta 357, 173–179. doi: 10.1016/j.cccn.2005.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Rejniak, K. A., and Anderson, A. R. A. (2011). Hybrid models of tumor growth. Wiley Interdiscip Rev. Syst. Biol. 3, 115–125. doi: 10.1002/wsbm.102

PubMed Abstract | CrossRef Full Text | Google Scholar

Rejniak, K. A., and McCawley, L. J. (2010). Current trends in mathematical modeling of tumor-microenvironment interactions: a survey of tools and applications. Exp. Biol. Med. 235, 411–423. doi: 10.1258/ebm.2009.009230

PubMed Abstract | CrossRef Full Text | Google Scholar

Roose, T., Chapman, S. J., and Maini, P. K. (2007). Mathematical models of avascular tumor growth. SIAM Rev. 49, 179–208. doi: 10.1137/S0036144504446291

CrossRef Full Text | Google Scholar

Rowlinson, J. S. (1979). Translation of vanderwaals,jd the thermodynamic theory of capillarity under the hypothesis of a continuous variation of density. J. Stat. Phys. 20, 197–244. doi: 10.1007/BF01011513

CrossRef Full Text | Google Scholar

Rude, U. (1993). Fully adaptive multigrid methods. SIAM J. Numer. Anal. 30, 230–248. doi: 10.1137/0730011

CrossRef Full Text | Google Scholar

Scianna, M., Bell, C. G., and Preziosi, L. (2013). A review of mathematical models for the formation of vascular networks. J. Theor. Biol. 333, 174–209. doi: 10.1016/j.jtbi.2013.04.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Sciumé, G., Gray, W. G., Hussain, F., Ferrari, M., Decuzzi, P., and Schrefler, B. A. (2014a). Three phase flow dynamics in tumor growth. Comp. Mech. 53, 465–484. doi: 10.1007/s00466-013-0956-2

CrossRef Full Text | Google Scholar

Sciumé, G., Santagiuliana, R., Ferrari, M., Decuzzi, P., and Schrefler, B. A. (2014b). A tumor growth model with deformable ECM. Phys. Biol. 11:065004. doi: 10.1088/1478-3975/11/6/065004

PubMed Abstract | CrossRef Full Text | Google Scholar

Sciumé, G., Shelton, S., Gray, W. G., Miller, C. T., Hussain, F., Ferrari, M., et al. (2013). A multiphase model for three-dimensional tumor growth. New J. Phys. 15:015005. doi: 10.1088/1367-2630/15/1/015005

PubMed Abstract | CrossRef Full Text | Google Scholar

Swartz, M. A. (2014). Immunomodulatory roles of lymphatic vessels in cancer progression. Cancer Immunol. Res. 2, 701–707. doi: 10.1158/2326-6066.CIR-14-0115

PubMed Abstract | CrossRef Full Text | Google Scholar

Tlsty, T. D., and Coussens, L. M. (2006). Tumor stroma and regulation of cancer development. Ann. Rev. Pathol. 1, 119–150. doi: 10.1146/annurev.pathol.1.110304.100224

PubMed Abstract | CrossRef Full Text | Google Scholar

Tracqui, P. (2009). Biophysical Models of Tumour Growth. Reports on Progress in Physics.

Trottenberg, U., Oosterlee, C. W., and Schüller, A. (2001). Multigrid. San Diego, CA: Academic Press.

Vineis, P., Schatzkin, A., and Potter, J. D. (2010). Models of carcinogenesis: an overview. Carcinogenesis 31, 1703–1709. doi: 10.1093/carcin/bgq087

PubMed Abstract | CrossRef Full Text | Google Scholar

Whiteside, T. L. (2008). The tumor microenvironment and its role in promoting tumor growth. Oncogene 27, 5904–5912. doi: 10.1038/onc.2008.271

PubMed Abstract | CrossRef Full Text | Google Scholar

Wise, S., Kim, J., and Lowengrub, J. (2007). Solving the regularized, strongly anisotropic Cahn-Hilliard equation by an adaptive nonlinear multigrid method. J. Comp. Phys. 226, 414–446. doi: 10.1016/j.jcp.2007.04.020

CrossRef Full Text | Google Scholar

Wise, S. M., Lowengrub, J. S., and Cristini, V. (2011). An adaptive multigrid algorithm for simulating solid tumor growth using mixture models. Math. Comp. Model. 53, 1–20. doi: 10.1016/j.mcm.2010.07.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Wise, S. M., Lowengrub, J. S., Frieboes, H. B., and Cristini, V. (2008). Three-dimensional multispecies nonlinear tumor growth - I - Model and numerical method. J. Theor. Biol. 253, 524–543. doi: 10.1016/j.jtbi.2008.03.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, M., Frieboes, H. B., McDougall, S. R., Chaplain, M. A. J., Cristini, V., and Lowengrub, J. (2013). The effect of interstitial pressure on tumor growth: coupling with the blood and lymphatic vascular systems. J. Theor. Biol. 320, 131–151. doi: 10.1016/j.jtbi.2012.11.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, A. J. M., Fleming, P. D., and Gibbs, J. H. (1976). Molecular theory of surface-tension. J. Chem. Phys. 64, 3732–3747. doi: 10.1063/1.432687

CrossRef Full Text | Google Scholar