On Computational Modeling in Tumor Growth

Abstract

The paper addresses modeling of avascular and vascular tumor growth within the framework of continuum mechanics and the adopted numerical solution strategies. The models involve tumor cells, both viable and necrotic, healthy cells, extracellular matrix (ECM), interstitial fluid, neovasculature and co-opted blood vessels, nutrients, waste products, and their interaction and evolution. Attention is focused on the more recent models which are much richer than earlier ones, i.e. they address more aspects of this complicated problem. An important element is how the governing equations are obtained and how the many interfaces between the above listed components are dealt with. These considerations suggest the definition of different classes of models comprised of diffusion, single phase flow and multiphase flow models with or without a solid phase. A multiphase flow model in a deforming porous medium (ECM) is chosen as reference model since it appears to invoke the least number of simplifying assumptions and has the largest potential for further development. The strategies adopted in the choice of the many model dependent constitutive relationships are discussed in detail. Two applications referring to two different model classes conclude the paper.

Appendices

Appendix A: Linear Momentum Balance Equation for a Fluid Phase

The general conservation of momentum Eq. (15) will be denoted for the fluid phase using the letter f as a qualifier. In the paper the f will be specified as either l, h or t.

$$ \begin{aligned}[b] &\frac{\partial ( \varepsilon^{f}\rho^{f}\mathbf{v}^{\overline{f}} )}{\partial t} + \nabla \cdot \bigl( \varepsilon^{f}\rho^{f}\mathbf{v}^{\overline{f}} \mathbf{v}^{\overline{f}} \bigr) - \nabla \cdot \bigl( \varepsilon^{f} \mathbf{t}^{\overline{\overline{f}}} \bigr) \\ &\quad {}- \varepsilon^{f}\rho^{f}\mathbf{g}^{\overline{f}} - \sum _{\kappa \in \Im_{cf}} \biggl( \sum_{i \in \Im_{s}} \mathop{M_{v}}\limits ^{i\kappa \to if}\mathbf{v}^{\overline{f}} + \mathop{\mathbf{T}}\limits^{\kappa \to f} \biggr) = 0 \end{aligned} $$

(A.1)

where $\mathbf{g}^{\overline{f}}$ is the body force, $\mathop{M_{v}}\limits^{i\kappa \to if}\mathbf{v}^{\overline{f}}$represents the momentum exchange from the κ to the f phase due to mass exchange of species i, $\mathop{\mathbf{T}}\limits^{\kappa \to f}$ is the interaction force between phase f and the adjacent interfaces, and $\mathbf{t}^{\overline{\overline{f}}}$ is the stress tensor. If the inertial terms are considered to be negligible, as is the case for slow flow in a porous medium, the first two terms in Eq. (A.1) can be neglected. Additionally, the momentum exchange due to mass transfer, $\mathop{M_{v}}\limits^{i\kappa \to if}\mathbf{v}^{\overline{f}}$may also be considered small since this term is of the same order of magnitude as the inertial terms. Thus the momentum equation simplifies to

$$ - \nabla \cdot \bigl( \varepsilon^{f} \mathbf{t}^{\overline{\overline{f}}} \bigr) - \varepsilon^{f}\rho^{f} \mathbf{g}^{\overline{f}} - \sum_{\kappa \in \Im_{cf}} \mathop{\mathbf{T}}\limits^{\kappa \to f} = 0 $$

(A.2)

The TCAT method of closure involves arranging terms in the entropy inequality into force-flux pairs. At equilibrium each member of the force-flux pair will be zero. This equilibrium constraint guides closure of the conservation system for near equilibrium situations. In the case here where the flows are slow, the near-equilibrium state assumption is appropriate. Based on the TCAT procedure, the elements of the entropy inequality relating to flow velocity that arise in the entropy inequality are

$$ \begin{aligned}[b] &{-} \frac{1}{\theta^{\overline{\overline{f}}}} \biggl[ \varepsilon^{f}\rho^{f}\mathbf{g}^{\overline{f}} + \varepsilon^{f}\rho^{f}\nabla \bigl( \varsigma^{\overline{f}} + \psi^{\overline{f}} \bigr) - \nabla \bigl( \varepsilon^{f}p^{f} \bigr) \\ &\quad {}+ \sum_{\kappa \in \Im_{cf}} \mathop{\mathbf{T}}\limits^{\kappa \to f} \biggr] \cdot \bigl( \mathbf{v}^{\overline{f}} - \mathbf{v}^{\overline{s}} \bigr) + \frac{\varepsilon^{f}}{\theta^{\overline{\overline{f}}}} \bigl( \mathbf{t}^{\overline{\overline{f}}} + p^{f} \mathbf{1} \bigr):\mathbf{d}^{\overline{\overline{f}}} \ge 0 \end{aligned} $$

(A.3)

In this equation, $\theta^{\overline{\overline{f}}}$ is the macroscale temperature of the f phase, $\psi^{\overline{f}}$ is the gravitational potential, $\varsigma^{\overline{f}}$ is the chemical potential, p ^f is the fluid pressure, $\mathbf{v}^{\overline{s}}$ is the velocity of the solid phase and $\mathbf{d}^{\overline{\overline{f}}}$ the rate of strain tensor of the phase $f (\mathbf{d}^{\overline{\overline{f}}} = \frac{1}{2} [ \nabla \mathbf{v}^{\overline{f}} + ( \nabla \mathbf{v}^{\overline{f}} )^{T} ])$. All of these quantities are macroscale averages.

Consider the variability in volume fraction of the f phase to be small. For this situation, $\nabla \psi^{\overline{f}} + \mathbf{g} = 0$. Additionally, consider an isothermal case such that the Gibbs-Duhem equation provides $\rho^{f}\nabla \varsigma^{\overline{f}} - \nabla p^{f} = 0$. Application of these two conditions to Eq. (A.3) reduces it to

$$ \begin{aligned}[b] &\frac{\varepsilon^{f}}{\theta^{\overline{\overline{f}}}} \bigl( \mathbf{t}^{\overline{\overline{f}}} + p^{f}\mathbf{1} \bigr):\mathbf{d}^{\overline{\overline{f}}} \\ &\quad {}- \frac{1}{\theta^{\overline{\overline{f}}}} \biggl[ - p^{f}\nabla \varepsilon^{f} + \sum _{\kappa \in \Im_{cf}} \mathop{\mathbf{T}}\limits^{\kappa \to f} \biggr] \cdot \bigl( \mathbf{v}^{\overline{f}} - \mathbf{v}^{\overline{s}} \bigr) \ge 0 \end{aligned} $$

(A.4)

This equation contains two independent force-flux products. The stipulation that both elements of each product pair must be zero at equilibrium and the requirement that the grouping of terms must be non-negative suggests the linear relations

$$ \sum_{\kappa \in \Im_{cf}} \mathop{\mathbf{T}}\limits^{\kappa \to f} - p^{f}\nabla \varepsilon^{f} = - \mathbf{R}^{f} \cdot \bigl( \mathbf{v}^{\overline{f}} - \mathbf{v}^{\overline{s}} \bigr) $$

(A.5)

and

$$ \mathbf{t}^{\overline{\overline{f}}} + p^{f}\mathbf{1} = \mathbf{A}^{f}:\mathbf{d}^{\overline{\overline{f}}} $$

(A.6)

In the first relation, R ^f is a symmetric, positive, semi-definite tensor accounting for the resistance to flow. In the second relation, A ^f is fourth order tensor that accounts for the dependence of the stress tensor on the rate of strain. At the macroscale for slow flow, this tensor is taken to be zero such that

$$ \mathbf{t}^{\overline{\overline{f}}} = - p^{f}\mathbf{1} $$

(A.7)

is the resulting form of the stress tensor. We note that this does not imply that the fluid is inviscid. The effects of viscosity are accounted for at the macroscale by the momentum exchange term $\mathop{\mathbf{T}}\limits^{\kappa \to f}$.

Substitution of the closure relations Eqs. (A.5) and (A.7) into Eq. (A.2) provides the momentum equation in the form

$$ \varepsilon^{f}\nabla p^{f} - \varepsilon^{f}\rho^{f}\mathbf{g}^{\overline{f}} + \mathbf{R}^{f} \cdot \bigl( \mathbf{v}^{\overline{f}} - \mathbf{v}^{\overline{s}} \bigr) = 0 $$

(A.8)

Typically this relation is expressed as

$$ - \mathbf{K}^{f} \cdot \bigl( \nabla p^{f} - \rho^{f}\mathbf{g}^{\overline{f}} \bigr) = \varepsilon^{f} \bigl( \mathbf{v}^{\overline{f}} - \mathbf{v}^{\overline{s}} \bigr) $$

(A.9)

where K ^f=(ε ^f)²(R ^f)⁻¹is called the hydraulic conductivity.

The hydraulic conductivity depends on the properties of both the flowing fluid and the solid porous material. For an isotropic medium, K ^f=K ^f 1. The morphology and topology of the solid media are important in determining the hydraulic conductivity of the cellular solid phases. The conductivity is influenced by the cell size distribution, shape of the cells, tortuosity of passages, specific surface area, and porosity (the sum of the fluid volume fractions). It also depends on the density and viscosity of the fluid. Neglecting gravity in Eq. (A.8) yields Eq. (17).

Appendix B: Coefficients of the Matrices Appearing in Eq. (42)

In the following equations K ^s is the Bulk modulus of the solid skeleton and $\frac{\partial \mathbf{e}_{sw}^{s}}{\partial t} = \frac{\mathbf{1}}{3K^{s}}\frac{\partial p^{s}}{\partial t}$

$$ \mathbf{C}_{nn} = \int_{\varOmega} \mathbf{N}_{n}^{T} \bigl( \varepsilon S^{l} \mathbf{N}_{n} \bigr) d\varOmega $$

(B.1)

$$ \begin{aligned}[b] \mathbf{C}_{tt} &= \int_{\varOmega} \mathbf{N}_{t}^{T} \biggl[ \varepsilon \mathbf{N}_{t} + \frac{S^{t}}{K^{s}} \biggl( S^{t}\frac{\partial p^{t}}{\partial S^{t}} + S^{h} \frac{\partial p^{h}}{\partial S^{t}} \\ &\quad {}+ p^{t} - p^{l} \biggr) \mathbf{N}_{t} \biggr] d\varOmega \end{aligned} $$

(B.2)

$$ \begin{aligned}[b] &\mathbf{C}_{th} = \int_{\varOmega} \mathbf{N}_{t}^{T} \biggl[ \frac{S^{t}}{K^{s}} \biggl( S^{t}\frac{\partial p^{t}}{\partial S^{h}} + S^{h}\frac{\partial p^{h}}{\partial S^{h}} + p^{h} - p^{l} \biggr)\mathbf{N}_{h} \biggr] d \varOmega \end{aligned} $$

(B.3)

$$ \mathbf{C}_{tl} = \int_{\varOmega} \mathbf{N}_{t}^{T} \biggl[ \frac{S^{t}}{K^{s}} \biggl( S^{t}\frac{\partial p^{t}}{\partial p^{l}} + S^{h}\frac{\partial p^{h}}{\partial p^{l}} + S^{l} \biggr)\mathbf{N}_{l} \biggr] d\varOmega $$

(B.4)

$$ \begin{aligned}[b] &\mathbf{C}_{ht} = \int_{\varOmega} \mathbf{N}_{h}^{T} \biggl[ \frac{S^{h}}{K^{s}} \biggl( S^{t}\frac{\partial p^{t}}{\partial S^{t}} + S^{h}\frac{\partial p^{h}}{\partial S^{t}} + p^{t} - p^{l} \biggr)\mathbf{N}_{t} \biggr] d \varOmega \end{aligned} $$

(B.5)

$$ \begin{aligned}[b] \mathbf{C}_{hh} &= \int_{\varOmega} \mathbf{N}_{h}^{T} \biggl[ \varepsilon \mathbf{N}_{h} + \frac{S^{h}}{K^{s}} \biggl( S^{t}\frac{\partial p^{t}}{\partial S^{h}} + S^{h} \frac{\partial p^{h}}{\partial S^{h}} \\ &\quad {}+ p^{h} - p^{l} \biggr) \mathbf{N}_{h} \biggr] d\varOmega \end{aligned} $$

(B.6)

$$ \mathbf{C}_{hl} = \int_{\varOmega} \mathbf{N}_{h}^{T} \biggl[ \frac{S^{h}}{K^{s}} \biggl( S^{t}\frac{\partial p^{t}}{\partial p^{l}} + S^{h}\frac{\partial p^{h}}{\partial p^{l}} + S^{l} \biggr)\mathbf{N}_{l} \biggr] d\varOmega $$

(B.7)

$$ \begin{aligned}[b] &\mathbf{C}_{lt} = \int_{\varOmega} \mathbf{N}_{l}^{T} \biggl[ \frac{1}{K^{s}} \biggl( S^{t}\frac{\partial p^{t}}{\partial S^{t}} + S^{h}\frac{\partial p^{h}}{\partial S^{t}} + p^{t} - p^{l} \biggr)\mathbf{N}_{t} \biggr] d \varOmega \end{aligned} $$

(B.8)

$$ \begin{aligned}[b] &\mathbf{C}_{lh} = \int_{\varOmega} \mathbf{N}_{l}^{T} \biggl[ \frac{1}{K^{s}} \biggl( S^{t}\frac{\partial p^{t}}{\partial S^{h}} + S^{h}\frac{\partial p^{h}}{\partial S^{h}} + p^{h} - p^{l} \biggr)\mathbf{N}_{h} \biggr] d \varOmega \end{aligned} $$

(B.9)

$$ \mathbf{C}_{ll} = \int_{\varOmega} \mathbf{N}_{l}^{T} \biggl[ \frac{1}{K^{s}} \biggl( S^{t}\frac{\partial p^{t}}{\partial p^{l}} + S^{h}\frac{\partial p^{h}}{\partial p^{l}} + S^{l} \biggr)\mathbf{N}_{l} \biggr] d\varOmega $$

(B.10)

$$ ( \mathbf{C}_{uu} )_{ij} = - \int _{\varOmega} \mathbf{B}^{\operatorname{T}} \mathbf{D}_{s} \mathbf{B} d\varOmega $$

(B.11)

$$ \mathbf{K}_{nn} = \int_{\varOmega} ( \nabla \mathbf{N}_{n} )^{T} \bigl( \varepsilon S^{l}D_{eff}^{\overline{nl}}\nabla \mathbf{N}_{n} \bigr)d\varOmega $$

(B.12)

$$ \mathbf{K}_{tt} = \int_{\varOmega} ( \nabla \mathbf{N}_{t} )^{T} \biggl( \frac{k_{rel}^{t}\mathbf{k}^{ts}}{\mu^{t}} \frac{\partial p^{t}}{\partial S^{t}}\nabla \mathbf{N}_{t} \biggr)d\varOmega $$

(B.13)

$$ \mathbf{K}_{th} = \int_{\varOmega} ( \nabla \mathbf{N}_{t} )^{T} \biggl( \frac{k_{rel}^{t}\mathbf{k}^{ts}}{\mu^{t}} \frac{\partial p^{t}}{\partial S^{h}}\nabla \mathbf{N}_{h} \biggr)d\varOmega $$

(B.14)

$$ \mathbf{K}_{tl} = \int_{\varOmega} ( \nabla \mathbf{N}_{t} )^{T} \biggl( \frac{k_{rel}^{t}\mathbf{k}^{ts}}{\mu^{t}} \frac{\partial p^{t}}{\partial p^{l}}\nabla \mathbf{N}_{l} \biggr)d\varOmega $$

(B.15)

$$ \mathbf{K}_{ht} = \int_{\varOmega} ( \nabla \mathbf{N}_{h} )^{T} \biggl( \frac{k_{rel}^{h}\mathbf{k}^{hs}}{\mu^{h}} \frac{\partial p^{h}}{\partial S^{t}}\nabla \mathbf{N}_{t} \biggr)d\varOmega $$

(B.16)

$$ \mathbf{K}_{hh} = \int_{\varOmega} ( \nabla \mathbf{N}_{h} )^{T} \biggl( \frac{k_{rel}^{h}\mathbf{k}^{hs}}{\mu^{h}} \frac{\partial p^{h}}{\partial S^{h}}\nabla \mathbf{N}_{h} \biggr)d\varOmega $$

(B.17)

$$ \mathbf{K}_{hl} = \int_{\varOmega} ( \nabla \mathbf{N}_{h} )^{T} \biggl( \frac{k_{rel}^{h}\mathbf{k}^{hs}}{\mu^{h}} \frac{\partial p^{h}}{\partial p^{l}}\nabla \mathbf{N}_{l} \biggr)d\varOmega $$

(B.18)

$$ \begin{aligned}[b] &\mathbf{K}_{lt} = \int_{\varOmega} ( \nabla \mathbf{N}_{l} )^{T} \biggl( \frac{k_{rel}^{t}\mathbf{k}^{ts}}{\mu^{t}} \frac{\partial p^{t}}{\partial S^{t}}\nabla \mathbf{N}_{t} \,{+}\, \frac{k_{rel}^{h}\mathbf{k}^{hs}}{\mu^{h}} \frac{\partial p^{h}}{\partial S^{t}}\nabla \mathbf{N}_{t} \biggr)d\varOmega \end{aligned} $$

(B.19)

$$ \begin{aligned}[b] \mathbf{K}_{lh} &= \int_{\varOmega} ( \nabla \mathbf{N}_{l} )^{T} \biggl( \frac{k_{rel}^{t}\mathbf{k}^{ts}}{\mu^{t}} \frac{\partial p^{t}}{\partial S^{h}}\nabla \mathbf{N}_{h} \\ &\quad {}+ \frac{k_{rel}^{h}\mathbf{k}^{hs}}{\mu^{h}} \frac{\partial p^{h}}{\partial S^{h}}\nabla \mathbf{N}_{h} \biggr)d\varOmega \end{aligned} $$

(B.20)

$$ \begin{aligned}[b] \mathbf{K}_{ll} &= \int_{\varOmega} ( \nabla \mathbf{N}_{l} )^{T} \biggl( \frac{k_{rel}^{t}\mathbf{k}^{ts}}{\mu^{t}} \frac{\partial p^{t}}{\partial p^{l}}\nabla \mathbf{N}_{l} + \frac{k_{rel}^{h}\mathbf{k}^{hs}}{\mu^{h}} \frac{\partial p^{h}}{\partial p^{l}}\nabla \mathbf{N}_{l} \\ &\quad {}+ \frac{k_{rel}^{l}\mathbf{k}^{ls}}{\mu^{l}}\nabla \mathbf{N}_{l} \biggr)d\varOmega \end{aligned} $$

(B.21)

$$ \begin{aligned}[b] &\mathbf{f}_{n} = \int_{\varOmega} \mathbf{N}_{n}^{T} \biggl( \frac{1}{\rho} \Bigl( \omega^{\overline{nl}}\mathop{M}\limits ^{l \to t} - \mathop{M}\limits ^{nl \to t} \Bigr) - \varepsilon S^{l}\mathbf{v}^{\overline{l}} \cdot \nabla \omega^{\overline{nl}} \biggr) d\varOmega \end{aligned} $$

(B.22)

$$ \begin{aligned}[b] \mathbf{f}_{t} &= \int_{\varOmega} \mathbf{N}_{t}^{T} \biggl[ \frac{1}{\rho} \mathop{M} \limits_{growth}^{l \to t} - S^{t}\operatorname{tr} \biggl( \frac{\partial \mathbf{e}^{s}}{\partial t} - \frac{\partial \mathbf{e}_{sw}^{s}}{\partial t} \biggr) \\ &\quad {}- \nabla \bigl( S^{t} \bigr) \cdot \biggl( \varepsilon \frac{\partial \mathbf{u}^{s}}{\partial t} \biggr) \biggr] d\varOmega \end{aligned} $$

(B.23)

$$ \begin{aligned}[b] \mathbf{f}_{h} &= \int_{\varOmega} \mathbf{N}_{h}^{T} \biggl[ - S^{h} \operatorname{tr} \biggl( \frac{\partial \mathbf{e}^{s}}{\partial t} - \frac{\partial \mathbf{e}_{sw}^{s}}{\partial t} \biggr) \\ &\quad {}- \nabla \bigl( S^{h} \bigr) \cdot \biggl( \varepsilon \frac{\partial \mathbf{u}^{s}}{\partial t} \biggr) \biggr] d\varOmega \end{aligned} $$

(B.24)

$$ \mathbf{f}_{l} = \int_{\varOmega} \mathbf{N}_{l}^{T} \biggl[ - \operatorname{tr} \biggl( \frac{\partial \mathbf{e}^{s}}{\partial t} - \frac{\partial \mathbf{e}_{sw}^{s}}{\partial t} \biggr) \biggr] d\varOmega $$

(B.25)

$$ \mathbf{f}_{u} = \int_{\varOmega} \mathbf{B}^{T} \biggl( \mathbf{D}_{s}\frac{\partial \mathbf{e}_{vp}^{s}}{\partial t} \biggr) d\varOmega + \int_{\varOmega} \mathbf{B}^{T} \biggl( \mathbf{D}_{s}\frac{\partial \mathbf{e}_{sw}^{s}}{\partial t} \biggr) d\varOmega $$

(B.26)