Tumor growth towards lower extracellular matrix conductivity regions under Darcy’s Law and steady morphology

Abstract

We study a classic Darcy’s law model for tumor cell motion with inhomogeneous and isotropic conductivity. The tumor cells are assumed to be a constant density fluid flowing through porous extracellular matrix (ECM). The ECM is assumed to be rigid and motionless with constant porosity. One and two dimensional simulations show that the tumor mass grows from high to low conductivity regions when the tumor morphology is steady. In the one-dimensional case, we proved that when the tumor size is steady, the tumor grows towards lower conductivity regions. We conclude that this phenomenon is produced by the coupling of a special inward flow pattern in the steady tumor and Darcy’s law which gives faster flow speed in higher conductivity regions.

Appendix

1.1 Basics of Darcy’s law and mass transport in porous media

This section extracts some information useful for our discussion from a classic book by Bear (Bear 1972). Above all, notice that in the continuum approach of fluid dynamics in porous media, all the kinematic and dynamic variables are defined as averaged values in a proper volume (representative elementary volume), such as porosity, velocity, pressure. Darcy’s law for homogeneous fluid flowing through inhomogeneous isotropic porous media can be written as

$$\begin{aligned} \mathbf {q}=-K \nabla \varphi , \end{aligned}$$

(34)

where $\mathbf {q}$ is the specific flux with dimension length/time, measuring the volume of fluid discharged through porous media per unit area per unit time. Let $\phi$ be the porosity (the fraction of pores or void space in a unit volume) and $\mathbf {u}$ be the flow velocity of the fluid in the pores, then relation between specific flux and velocity is expressed with the following Dupuit-Forchheimer equation,

$$\begin{aligned} \mathbf {q} = \phi \mathbf {u}. \end{aligned}$$

(35)

The value of $\phi$ is between 0 and 1 and often expressed by a percentage. The porosity is defined in such a way that around any point in the porous media, the void space fraction is $\phi$. The coefficient $K = K(\mathbf {x})$, $\mathbf {x}\in {\mathbb {R}}^3$, is the hydraulic conductivity and a scalar variable for inhomogeneous isotropic media. Its expression is

$$\begin{aligned} K=\frac{\rho g k}{\mu }, \end{aligned}$$

(36)

where $\rho$ is the fluid density, g is the gravity acceleration, k is the intrinsic permeability (dependent solely on the property of the solid matrix) with dimension $length^2$, $\mu$ is the dynamic viscosity of the fluid. Thus, the dimension of K is length/time. The quantity $\varphi$ is the piezometric or hydraulic head and given by

$$\begin{aligned} \varphi =\frac{p}{\rho g} + z. \end{aligned}$$

(37)

Here, p is hydrostatic pressure and z is the height. All these quantities with their dimensions are listed in Table 3. Using relation (35), the Darcy’s law (34) can be rewritten as

$$\begin{aligned} \mathbf {u} = - \frac{K}{\phi } \nabla \varphi . \end{aligned}$$

(38)

Table 3 Quantities used in Darcy’s law for porous media. Here the units are the ones usually used in groundwater hydrology

The mass balance equation for the tumor cell component according to (Bear 1972) (equation 4.3.1 therein) is

$$\begin{aligned} \frac{ \partial (\rho \phi )}{\partial t} + \nabla \cdot (\rho \phi \mathbf {u}) = I_t, \end{aligned}$$

(39)

where $I_t$ is the tumor cell source term. Under the assumption of constant density $\rho$ and fixed porosity in time, $\frac{\partial (\rho \phi )}{\partial t}=0$. The source of tumor cells given in (Cristini et al. 2003) is $\rho G(n-A)$, where G is the division rate and A is the apoptotic threshold. Because the tumor cell source term must be only present in the void space, therefore, $I_t=\phi \rho G(n-A)$. Then the mass balance equation is reduced to

$$\begin{aligned} \nabla \cdot (\phi \mathbf {u}) = \phi G(n-A). \end{aligned}$$

(40)

In (Cristini et al. 2003), the mass equation is simply

$$\begin{aligned} \nabla \cdot \mathbf {u} = G(n-A). \end{aligned}$$

(41)

Compared with (40), this must be done by assuming $\phi$ is a constant in space. For the sake of simplicity and comparison, we accept the assumption here.

Define ${\tilde{p}}=\frac{\varphi }{\phi }$ and use the definition of K in (36), then (38) becomes

$$\begin{aligned} \mathbf {u} = -K \nabla {\tilde{p}}. \end{aligned}$$

(42)

Ignoring z in (37) and dropping the tilde for ${\tilde{p}}$, we get (1).