On the numerical discretization of a tumor progression model driven by competing migration mechanisms

1.   Introduction

1.1.   Background

The phenomenon of lymph node metastasis has been recognized for a long time. However, the underlying mechanism by which malignant tumor cells leave the primary tumor site, invade the lymphatics and metastasize to lymph nodes are unclear ^[11]. In ^[12,24] it is suggested that interstitial flow (IF) caused by lymphatic drainage directs tumor cell migration by utilizing interstitial flow to create and amplify gradients of a chemical agent (chemokine) and thus move by chemotaxis toward the lymphatic in a process termed autologous chemotaxis. Polacheck et al. ^[16] extended the study by Shields et al. ^[24], demonstrating that the strength of flow as well as the cell seeding density affected the migration direction. Their work provides further evidence that autologous chemotaxis is the mechanism that leads to cancer cell migration with the flow, but there is another mechanism that causes strain-induced migration against the flow. It was shown that a low cell seeding density culture migrated with the flow in accordance with the behavior reported in ^[24]. However, for a high cell seeding density the migration was dominated by upstream migration. In addition, this migration against the flow direction was sensitive to the fluid velocity. Higher fluid velocity increased the upstream migration. Moreover, when chemotaxis was blocked, both high and low density cultures migrated upstream. Polacheck et al. further investigated the effects of interstitial fluid flow stresses imparted on cells ^[17]. A first attempt to describe the phenomenon of autologous chemotaxis theoretically in a continuum setting where both cancer cells and fluid are represented through separate mass and momentum equations ^[2], can be found in ^[30]. The model was based on a biphasic cell-fluid-matrix formulation where matrix is assumed to be stagnant occupying a certain volume $1-\phi$ of the bulk volume whereas cells and fluid compete for space in the remaining void volume described by the porosity $\phi$. Various model parameters were determined to comply with the experimental results in ^[12,24] regarding a realistic fluid velocity and cell migration velocity. In particular, efforts were made to explain the role of a governing component of the cell velocity reflecting chemotaxis-driven motion toward a chemical agent (chemokine) that is skewed in the flow direction. A step further was taken in ^[31] where the biphasic model was extended to account for the mechanically induced upwind migration component reported in ^[16]. (Further modeling of this phenomenon can be found in ^[23]). The resulting model was shown to capture many of the trends reported in ^[16] concerning the internal competition between the downstream chemotaxis migration and the fluid-stress related upstream migration. However, an open question is to what extent these competing mechanisms can play a role as a possible driver for metastatic behavior? That is, to what extent are tumor cells at the periphery of the primary tumor able to detach from the primary tumor, a behavior that might pave the way for lymph node metastasis where cancer cells escape to nearby lymphatic vessels? For further discussion and some preliminary investigations of that issue we refer the reader to ^[8].

1.2.   Main objectives

The purpose of the current work is to take a closer look at different aspects of the biphasic cell-fluid model and try to put the model into a larger context. We start reviewing a slightly simplified version of the model studied in ^[8,31] where the chemical agents that are involved have been reduced from three to one. This is done while maintaing the same competing migration mechanisms. The first question we are interested in is: How can we make use of a general and robust discretization method to obtain numerical reliable approximations of the true solution? The idea we employ is to borrow techniques from a recent study of a cell-fluid model, which represents a biphasic description of the model introduced in the seminal work ^[29] of cell swimming, that was explored in ^[19]. Secondly, we use the novel discretization and explore in a one-dimensional tumor setting to what extent the model possesses a detachment mechanism, as explained above. The third question we pay attention to is: Can we derive a version of the cell-fluid-matrix model that resembles the Tuval model presented in ^[29]? This is partly motivated by the fact that this model appears to be ideal for systematical studies of the interplay between fluid dynamics, as described by the incompressible (Navier-)Stokes formulation, and chemotactic-driven cell motion. The model has released a rather intense study of finer details, as revealed through mathematical analysis, involved in the competing forces that ultimately dictate the cell-fluid dynamics ^{[1,3,6,13,27,28,33]}. The fourth question is a natural consequence of this last model: How can we exploit the fact that the fluid dynamics apparently depend weakly on the explicit presence of cancer cells, whereas cancer cells are more sensitive to the fluid behavior?

2.   A cell-fluid-matrix model

2.1.   General multiphase formulation

In the multiphase modeling framework, the tumor-host environment is considered as a mixture of three interacting continua ^[2,5,9,22]: the stagnant matrix, whose main component is extracellular matrix (ECM), that is described by a volume fraction $1-\phi$; the cellular phase comprises tumor cells represented by a volume fraction $\alpha_c$ and density $\rho_c$ moving within the pore space represented by porosity $\phi$ with an interstitial velocity $\mathbf{u}_c^p$; and the extracellular fluid phase represented by the volume fraction $\alpha_w$ and density $\rho_w$ moving with a velocity $\mathbf{u}_w^p$.

The resulting model is summed up below, and we refer to ^[7,30,31] for more details:

with additional algebraic equations and correlations given by

The first two equations of (2.1) describe the mass balance for cancer cells and fluid, respectively, as they compete for the same pore space. Equation (2.2)$_1$ expresses that cells and fluid completely occupy the pore space. The blood circulation associated with the leaky and aberrant tumor vasculature is expressed through $Q = Q_v-Q_l$ involved in (2.1)$_2$. Generally, we may consider the phases as compressible fluids by specifying appropriate pressure-density closure relations $\rho_i = \rho_i(P_i)$, $i = c, w$. The two subsequent equations (2.1)$_{3, 4}$ describe the corresponding momentum balance. In particular, the pressure gradient $\nabla P_c$ (cell phase) is balanced against the cell-matrix resistance force $-\hat{\zeta}_c \mathbf{u}_c^p$ whereas the interaction term $\hat{\zeta}(\mathbf{u}_w^p- \mathbf{u}_c^p)$ accounts for a force generated by the cancer cells in response to the stress from the flowing fluid. The last term describes internal viscous effects. Similarly, (2.1)$_{4}$ describes how the fluid pressure gradient $\nabla P_w$ is balanced by a fluid-matrix resistance force $-\hat{\zeta}_w \mathbf{u}_w^p$ followed by the cell-fluid interaction term and internal viscosity. In Section 4 we provide further information how the parameters involved in $\hat{\zeta}_w$, $\hat{\zeta}_c$, and $\hat{\zeta}$ in (2.1)$_{3, 4}$ are set. The last equation (2.1)$_5$ accounts for the chemical agent $C$ which is secreted from the cancer cells and move by diffusion and fluid velocity $\mathbf{u}_w^p$.

The cell phase pressure $P_c$ feels additional stress due to cell-cell interaction through $\Delta P(\alpha_c)$ and migration-related stress due to chemotaxis through $\Lambda(C)$, as expressed by the relation $(2.2)_2$. Similar to ^[7,30,31], we use

and

where $C_M$ represents a maximal concentration. The reaction term $R_C$ in (2.2)$_{3}$ accounts for production of chemokine by the cancer cells themselves (the first bracket) and consumption and absorption through the lymphatic drainage (second bracket). The last line (2.2)$_{4}$ describes that the fluid produced by the blood vasculature in the tumor through $Q_v$ is controlled by the internal vascular pressure $\widetilde{P}_v^*$ and the conductivity associated with the blood vessel wall through $T_v$. Similarly for the lymphatic absorption of fluid through $Q_l$. We refer to Section 4 for more information how parameters are set, including (2.3) and (2.4).

Assuming a constant porosity $\phi$ where we rescale the time variable $t\rightarrow t/\phi$ and introduce a phase velocity $\mathbf{u}_c = \phi \mathbf{u}_c^p$ and $\mathbf{u}_w = \phi \mathbf{u}_w^p$, we may write the model in the form

Here we have also implicitly absorbed $\phi$ in the coefficients $\hat{\zeta}_c$, $\hat{\zeta}_w$, $\hat{\zeta}$, and $\varepsilon_c, \varepsilon_w$. We have also used the approximation $C \alpha_w \approx C$ for the transport equation describing $C$.

2.2.   Incompressible and inviscid version

Following ^[30] we assume fluids are incompressible. This allows us to derive an explicit expression for the interstitial cell velocity $\mathbf{u}_c$ which reflects the competition between different migration mechanisms as well the role played by the different interaction terms $\hat{\zeta}_w$, $\hat{\zeta}_c$, and $\hat{\zeta}$. From (2.1) we obtain the following dimensionless simplified version:

The model is combined with the boundary conditions

where $\nu$ is the outward normal on $\partial \Omega$ and $P^*$ a known pressure at the boundary. Corresponding initial data are

From the two momentum equations (2.6)$_{3, 4}$ we can compute explicit expressions for the cell and fluid velocity, respectively, $\mathbf{u}_c$ and $\mathbf{u}_w$, by ignoring the viscous terms, i.e., we set $\varepsilon_c = \varepsilon_w = 0$. We refer to ^[30] for details. The following expressions are found:

where $\mathbf{U}_T = \alpha_w \mathbf{u}_w+ \alpha_c \mathbf{u}_c$ and with fractional flow functions $\hat{f}_c(\alpha_c)$ and $\hat{f}_w(\alpha_c)$ given by

where the coefficients $\hat{\lambda}_c$, $\hat{\lambda}_w$, and $\hat{\lambda}_T$ (so-called mobility functions ^[37]) are given by

and total mobility $\hat{\lambda}_T$ given by

The model (2.6) then takes the more compact form:

where $\mathbf{u}_c$ is given by (2.9)$_1$. Note that we also need $\mathbf{u}_w$ given by (2.9)$_2$ to solve (2.13)$_{2}$. Moreover, in order to compute $\mathbf{U}_T$ we first solve the elliptic problem for $P_w$ (see ^[30] for details)

Knowing $P_w$, we can compute $\mathbf{U}_T$ needed in (2.9). We have the following expression for $\mathbf{U}_T$ ^[30]:

Previous discretization ^[30,31] has been based on the formulation (2.13)–(2.15) combined with the expressions (2.9)–(2.12). In the following, we describe another approach which is based on the original formulation (2.6) and which is motivated by more recent investigations of incompressible and compressible multiphase flow models for creeping fluids in a porous media ^{[18,20,21,25]}. These models do not involve any chemotaxis tranport mechanisms. The approach we explore in the current work builds on the discretization of a cell-fluid model presented in ^[19]. However, that model is nevertheless different from the tumor model we study here since it only involves a downstream chemotactic migration and not the mechanical upstream migration mechanism. The main advantage of this approach is that it can naturally deal with the more general case where both viscous terms (internal viscosity) are included as well as compressible fluids. In the following, for the sake of simplicity, we assume incompressibility.

2.3.   A general discretization method based on the original full model (2.6)

The dimensionless domain $\Omega = [0, 1]$ is considered for the one-dimensional model. We introduce the approximate cell volume fraction and chemokine concentration $\{\alpha_{cj}(t)\}_{j = 1}^{N}$ and $\{C_{j}(t)\}_{j = 1}^{N}$ and fluid pressure $\{P_{wj}(t)\}_{j = 1}^{N}$ associated with the nodes $\{x_j\}_{j = 1}^{N}$ whereas the approximate velocities $\{u_{w, j+1/2}\}_{j = 0}^{N}$ and $\{u_{c, j+1/2}\}_{j = 0}^{N}$ are associated with the grid block interfaces $\{x_{j+1/2}\}_{j = 0}^{N}$. That is, grid block center positions are given by

and corresponding grid block interfaces $x_{j+1/2}$ are given by

where $\Delta x = 1/N$.

We assume that we have an approximate solution $(\alpha_{c, j}^k, C_j^k, P_{w, j}^k, u_{c, j+1/2}^k, u_{w, j+1/2}^k)$ at time $t^k$. Next, we describe how to obtain a new approximate solution $(\alpha_{c, j}^{k+1}, C_j^{k+1}, P_{w, j}^{k+1}, u_{c, j+1/2}^{k+1}, u_{w, j+1/2}^{k+1})$ at time $t^{k+1}$ through a 2-step algorithm.

Step 1: Mass transport

We solve for $(\alpha_{c, j}^{k+1}, C_j^{k+1})$ by using the following discrete version of (2.6)$_1$ and (2.6)$_5$, respectively:

where

where

and

Then, we compute $\alpha_{w, j}^{k+1}$ from (2.2)$_1$. Next, we compute pressure and velocities simultaneously at time level $k+1$.

Step 2: Computation of velocities and pressure

We solve for $P_{w, j}^{k+1}$, $u_{c, j+1/2}^{k+1}$ and $u_{w, j+1/2}^{k+1}$ by considering the following algebraic system, which is a discrete version of the sum of (2.6)$_1$ and (2.6)$_2$, combined with (2.6)$_{3, 4}$:

which is combined with the momentum balance equations

For the cell-fluid interaction term $\hat{\zeta}(\alpha_c)$ we discretize as follows:

For cell volume fraction $\alpha_c$ at interface $j+1/2$ we use

For water volume fraction $\alpha_w$ interface $j+1/2$ we use

These terms appear in product with the pressure gradients appearing on the left-hand-side of (2.22).

The discretization of $[\hat{\zeta}_c^{k+1}]_{j+1/2}$ and $[\hat{\zeta}_w^{k+1}]_{j+1/2}$ appearing in (2.22) on the right-hand-side is based on upwind similar to the treatment of the terms $[\alpha_c]_{j+1/2}^{k+1}$ (2.24) and $[\alpha_w]_{j+1/2}^{k+1}$ (2.25), respectively.

3.   A reduced cell-fluid-matrix model

Now, we want to derive a reduced version of the model (2.6). The key element is that we seek a model where we have decoupled to some extent the fluid dynamics from the cancer cell dynamics. More precisely, we seek to take advantage of the fact that fluid velocity field largely is dictated by the fluid-matrix resistance force $-\hat{\zeta}_w \mathbf{u}_w$ in (2.6)$_4$ and the leaky vasculature and collecting lymphatic vessels through $Q = Q_v-Q_l$ in (2.6)$_2$. Cancer cells, on the other hand, apparently are more sensitive to the fluid. The resulting model may be considered as a counterpart of the Tuval's model that have been used for cell-fluid dynamics ^[29] but now in a porous media setting where interaction with the stagnant matrix also is involved.

3.1.   An incompressible version

We seek to take advantage of the fact that while the cancer cells are quite sensitive to the fluid flow, both through the downstream chemotaxis migration as well as the upstream fluid-stress driven migration, the interstitial fluid depends weakly on the presence of cancer cells. We may assume incompressible phases and consider an approximation and modification of (2.6) motivated by the rewritten form (2.13). The rationale behind the approximation is the following three points:

(i) Fluid velocity typically is a 100-fold larger than cell migration velocity ^[16,24]. This suggests that $\mathbf{U}_T\approx \alpha_w \mathbf{u}_w$.

(ii) We want to ensure that the basic relation (2.15), which determines the fluid velocity field through a Darcy-like equation, is taken into account.

(iii) Experimental work has revealed that interstitial flow may cause shear stress ranging from pulsatile and turbulent (near vasculature network) to primarily laminar convection, with fluid velocity influenced by interstitial porosity and pressure as well as vasculature density, permeability (conductivity), and viscoelastic properties ^[4,11]. This suggests that we use a full momentum equation for the fluid with inclusion of acceleration terms.

Taking these aspects into account, we suggest a cell-fluid-matrix model of the form

Here $\mathbf{u} = \mathbf{u}_w$ and $\rho$ refers to the fluid density (constant). In (3.1) we see the appearance of functions $\widetilde{f_c}$, $\widetilde{f_w}$, and $\widetilde{h}$ instead of $\hat{f_c}$, $\hat{f_w}$, and $\hat{h}$ given by (2.10) signalling that there is a slight difference. Similarly for $\widetilde{\lambda}_c$ and $\widetilde{\lambda}_T$. This is explained in the next subsection. A more precise justification of using (3.1) to represent the biphasic model (2.6) is provided through the numerical experiments presented in Section 4. We may rewrite (3.1) in the form (essentially using $\nabla\cdot (\alpha_w\mathbf{u}) = Q$)

This model bears similarity to the chemotaxis-Navier-Stokes model ^[29] when $\alpha_w\approx 1$ and $Q = 0$. A typical variant of this model is given by

Here $n$ is cell density, $D_n$ cell diffusion coefficient, $\mathbf{\Phi}$ gravitational potential, and $S(n, C)$ represents production/consumption of the chemical agent, possibly, under the influence of the cells $n$. Comparing (3.3)$_2$ and (3.1)$_2$, we see that the fluid field is affected by the cells through the buoyancy term $n\nabla \mathbf{\Phi}$ in the first model, whereas the latter is under influence of the cell population through the resistance force term $-\frac{\alpha_w}{\widetilde{\lambda}_T} \mathbf{u}$ on the right-hand-side and cell mobility-related term $\frac{\widetilde{\lambda}_c}{\widetilde{\lambda}_T}\nabla(\Delta P +\Lambda)$ of the momentum balance, as well as the equation $\nabla\cdot (\alpha_w\mathbf{u}) = Q$.

Remark 1. Numerical computations are carried out in Section 4 to explore features possessed by the reduced model (3.1). Some natural questions are: What is the difference between (3.1) and (2.6)? Can we ignore the term $\frac{\widetilde{\lambda}_c}{\widetilde{\lambda}_T}\nabla(\Delta P(\alpha_c)+\Lambda(C))$ in (3.1)$_2$? Can we replace $\alpha_w\mathbf{u}$ by $\mathbf{u}$ appearing in (3.1)$_1$ and (3.1)$_2$ and through that provide a further decoupling of cell and fluid dynamics? Can we take advantage of the looser coupling between fluid phase and cell phase and improve the efficiency of computing numerical solutions by doing fewer calculations of the fluid velocity $\mathbf{u}_w$ and pressure $P_w$?

3.2.   Derivation of the reduced model (3.1) from the biphasic model (2.6)

We want to build a numerical scheme for the reduced model (3.1) which makes use of the principles of the discrete scheme presented in Section 2.3 for the model (2.6). For that purpose we need to revisit calculations involved in Section 2.2 that give us expressions for the cell and fluid velocity $\mathbf{u}_c$ and $\mathbf{u}_w$ as well as the expressions for $\hat{f}_c$, $\hat{f}_w$, and $\hat{h}$. In particular, since the dicretization of $\alpha_c$ and $\alpha_w$ given by (2.24) and (2.25) not necessarily obey (2.2)$_1$ at the interface $j+1/2$, we should avoid using that.

From (2.6)$_{3, 4}$ we have the momentum balance when restricting ourselves to the one-dimensional setting (setting $\varepsilon_c = \varepsilon_w = 0$):

We solve for $u_c$ and $u_w$ from the two equations in (3.4) and find expressions of the form

Having these two phase velocities, the total velocity $U_T = \alpha_c u_c+ \alpha_w u_w$ can be obtained:

Therefore, from (3.6)

where

Equation (3.7) corresponds to the momentum equation in (3.1)$_2$ (ignoring the acceleration term and viscous term).

Next, let us see how we can identify an expression for the cell velocity $u_c$ that explain the form of the cell mass balance (3.1)$_1$. We can eliminate the phase pressure gradient $P_{wx}$ in (3.4) by multiplying by $\alpha_w$ and $\alpha_c$ in (3.4)$_1$ and (3.4)$_2$, respectively, and subtracting. Using that $u_w = \frac{U_T- \alpha_c u_c}{ \alpha_w}$ this calculation yields

where cell migration $u_c$ has three components: advective $u_{c1}$, diffusive $u_{c2}$ and autologous chemotaxis $u_{c3}$; and $\widetilde{f}_c$ and $\widetilde{h}$ take the form

Remark 2. It should be noted that (3.10) is equivalent to (2.10) in light of (2.2)$_1$. Similarly for (3.8) versus (2.11) and (2.12). In the derivation of (3.8) and (3.10) we avoid using the relation (2.2)$_1$. This ensures that the discretization of the reduced model (3.1), which is presented in the next section, is consistent with the discretization of the more general biphasic model (2.6).

3.3.   Numerical scheme implemented for the reduced model (3.1)

We assume that we have an approximate solution $(\alpha_{c, j}^k, C_j^k, P_{w, j}^k, u_{c, j+1/2}^k, u_{w, j+1/2}^k)$ at time $t^k$. Next, we describe how to obtain a new approximate solution $(\alpha_{c, j}^{k+1}, C_j^{k+1}, P_{w, j}^{k+1}, u_{c, j+1/2}^{k+1}, u_{w, j+1/2}^{k+1})$ at time $t^{k+1}$ through a 3-step algorithm.

Step 1: Mass transport

We express the cell mass balance equation (3.1)$_1$ in terms of the cell velocity (3.9) and (3.10). Then we solve for $(\alpha_{c, j}^{k+1}, C_j^{k+1})$ by using the following discrete version of (3.1)$_1$ and (3.1)$_3$, respectively:

where

where

and

Then, we compute $\alpha_{w, j}^{k+1}$ from (2.2)$_1$. Next, we compute pressure and fluid velocity simultaneously at time level $k+1$.

Step 2: Computation of fluid velocity and pressure

We solve for $P_{w, j}^{k+1}$ and $u_{w, j+1/2}^{k+1}$ by considering the following algebraic system which is a discrete version of (3.1)$_2$:

which is combined with the momentum balance equation

where

for cell

for water

Step 3: Computation of cell velocity

After obtaining the water velocity $u_{w, j+1/2}^{k+1}$, the cell velocity is calculated based on (3.9) and (3.10).

Herein, $\alpha_{c, j+1/2}^{k+1}$ and $\alpha_{w, j+1/2}^{k+1}$ appearing in (3.21) correspond to (3.19) and (3.20).

4.   Results

A 1D spatial domain is considered with a realistic length 0.01m. We first present the numerical solution of the full two-phase cell-fluid-matrix model (2.6). A main purpose is to illustrate the ability of the model to create isolated groups of cancer cells as a result of proper tuning of the downstream and upstream mechanism. We also explore the role of parameters that characterize the cancer cell phenotype (i.e., strength of downstream versus upstream migration through $\hat{\zeta}$), as well as parameters that affect the fluid drainage through $T_v(x)$, $T_l(x)$, and internal vasculature pressure $\widetilde{P}_l^*$. Next, we compare the solution of the reduced model (3.1) with the same input data as the original model (2.6) where we also ignore the acceleration terms $\rho(\mathbf{u}_t+ (\mathbf{u}\cdot\nabla) \mathbf{u})$ in the fluid momentum equation. For the cell-fluid regime we consider here, this part of the model plays a minor role. We also assess the possibility for further simplification of $(3.1)$ by ignoring more terms in the fluid momentum equation (3.1)$_2$ as mentioned in Remark 1. Finally, we check how to accelerate the computations of the reduced model by reducing the number of updates of $u_w$ and $P_w$ in $(3.1)_2$.

For input data required to solve the model we refer to Table 1. The choice of these parameters are largely set as in ^[31] where the model was aligned to the experimental behavior as reported in ^[16]. We compute solutions after a time T = 600 (dimensionless) with reference time $T^* = 10^4s$, which amounts to approximately 70 days.

Initial cell volume fraction and chemokine concentration

We consider an initial cell volume fraction (primary tumor) corresponding to

and initial chemokine is set to zero

Boundary conditions

The model is combined with the boundary conditions

where $\nu$ is the outward normal on $\partial \Omega$ and $P^* = 0$ bar relatively the reference pressure 1 bar.

Choice of input parameters

The following values are used for parameters related to the vascular flow, $Q_v$, given by $(2.2)_4$ and involved in the continuity equation for the fluid $(2.6)_2$:

and for lymphatic absorption $Q_l$ in (2.2)$_4$ we use

Here the values of $\widetilde{P}_v^*$ and $\widetilde{P}_l^*$ are given relatively the reference pressure 1 bar (752 mmHg). Regarding $\Delta P(\alpha_c)$ given by (2.3), we set $\gamma = 1000$ Pa. For the various interaction forces $\hat{\zeta}_c$, $\hat{\zeta}_w$, and $\hat{\zeta}$ involved in (2.6)$_{3, 4}$ we use the following correlations ^[15,26,31]:

Parameters for the fluid-matrix interaction is set to

and corresponding values for the cell-matrix interaction are set as

whereas the following set of parameters are used for the cell-fluid interaction term $\hat{\zeta}$:

The negative sign of $I$ reflects a cell-generated mechanical force against the fluid flow direction. In Figure 1 we have plotted the fractional flow function $\hat{f}_c(\alpha_c)$ given by (3.10)$_1$ and coefficient function $\hat{h}(\alpha_c)$ given by (3.10)$_2$. Due to the negative coefficient $I$ in the cell-fluid term $\hat{\zeta}$, the fractional flow function $\hat{f}_c(\alpha_c)$ contains a negative dip that gives rise to an upstream advective transport through $u_{c1}$, see (3.9) and (3.10). The upstream migration becomes larger with larger values of $u_w$ and increases with higher volume fraction of cells.

4.1.   The cell-fluid-matrix model (2.6)

We first illustrate a cell progression pattern after a time $T = 600$ by using the model (2.6) and the corresponding numerical scheme described in Section 2.3. Figure 2(A) shows the relevant position of the leaky intratumoral vasculature (through $T_v$) and collecting peritumoral lymphatic vessels (through $T_l$). As seen from Figure 2(A), lymphatic vessels are essentially non-functional in the initial cell region and the fluid is forced to flow towards the peritumoral region, which can be seen from the velocity $u_w$ profile in Figure 2(D). Note that the positive value of velocity $u_w$ represents flow from left to right in the domain, vice versa. This outward flow of fluid gives rise to the fluid pressure $P_w$ (blue) shown in Figure 2(F). Chemokine $C$ is secreted by the cancer cells and goes with the fluid flow towards the lymphatic region where it tends to accumulate and create higher concentrations, see panel (C). This triggers a chemotactic-driven outgoing cell migration. At the same time, the strong fluid velocity found at the tumor periphery generates a relatively strong mechanical-driven upstream migration. What is interesting is that the cells' upstream migration due to fluid flow is almost offset by the velocity components due to diffusion and chemotaxis in the middle region of domain where there seems to be close to zero total velocity, see Figure 2(E). However, the velocity caused by chemotaxis dominates the total cell velocity in the region outside the primary tumor where the upstream migration effect is relatively low due to the low cell volume fraction, thus leading to detachment of cancer cells on both sides of the primary tumor, see Figure 2(B). The net effect is, as seen from panel (B), that the primary tumor is largely kept intact while islands of cancer cells are formed in the peritumoral region. This demonstrates how the experimentally observed migration mechanisms reported in ^[16,17,24] may serve as an effective means for lymph node metastases.

What is the role of the upstream migration for generating cancer cell detachment

In this example, we modify the cell-fluid interaction term in (4.9) and reduce the upstream effect by a factor 2 by setting $\hat{k} = 2\hat{k}_w$. Figure 3 illustrates the corresponding result. In panel (B) we see that most of the cells move towards the lymphatic region since the upstream velocity component is reduced. Most strikingly, there is no detachment mechanism that can create cell islands, as shown in Figure 2(B). The chemokine concentration profile in panel (C) is quite smooth due to a higher consumption rate with higher cell volume fraction. As a consequence, as seen in panel (E), both the upstream component and the chemotatic component of cell velocity are rather small compared to the case in Figure 2(E). This shows how the upstream migration of cells caused by fluid flow plays a vital role in forming detached islands of cells.

What is the effect of the vasculature and lymphatic conductivities $T_v$ and $T_l$

In this example, we use a heterogeneous distribution of the vasculature and lymphatic conductivities $T_v$ and $T_l$, see panel (A) in Figure 4, rather than a constant value as in the example in Figure 2(A). A first observation from Figure 4(B) is that a new isolated island is formed on each side of the primary tumor, in addition to the old one shown in Figure 2(B). This is consistent with the non-uniform distribution of the lymphatic vessels through $T_l$ shown in Figure 4(A). Looking at the cell velocity components in Figure 4(E), it is clear that the upstream mechanical-driven cell velocity (pink) and downstream chemical-driven (black) have opposite sign. The detachment mechanism observed for the base case (in Figure 2) is at work again and is able to generate a second island, by the help of the non-uniform absorption of fluid described through $T_l$. The heterogeneous lymphatic vessels influence the source term in the chemokine equation (2.20) and fluid velocity (Figure 4(D)) through the mass transport equation (2.21).

What is the role of the internal lymphatic vasculature pressure $\widetilde{P}_l^*$

The role of lymph nodes for cancer cell migration to and through lymphatic vessels, survival in draining lymph nodes and further spread to other distant organs, is currently intensely debated ^[11,32,38]. In Figure 5 we run the base again but the value of the internal lymphatic vasculature pressure $\widetilde{P}_l^*$ in $(4.5)$ take two different values: for the left region $\widetilde{P}_l^* = -7.5\;$mmHg and for the right region $\widetilde{P}_l^* = -3\;$mmHg. As we can expect, the fluid pressure $P_w$ in Figure 5(F) takes a lower value in the left region (i.e., a stronger absorption of fluid) compared to the right side. Moreover, the maximal water velocity $u_w$ in the right domain is lower than the one in the left side. This affects the distribution of chemokine $C$ within the fluid phase seen in panel (C). As a result, the cell velocity component due to chemotaxis seen in Figure 5(E) is smaller in the peritumoral region on the right side. In turn, the developed isolated island of cells on the right side of the primary tumor becomes smaller, as compared to the left hand side island.

4.2.   The reduced cell-fluid-matrix model (3.1)

Now our question is: can we use the reduced cell-fluid model (3.1) to reproduce the same behavior as seen for the full model in Figure 2? The potential benefits of this investigation is that it can not only help us to speed up our simulations, especially when there is a need for huge amounts of output data, but it can also provide further insight into the cell-fluid migration mechanisms. We use the same input data as the base case in Figure 2 and make use of the reduced model (3.1) with its corresponding numerical scheme developed in Sec.3.3. Figure 6 shows the result for the reduced model and the plots are compared with those from the full two-fluid model (2.6). Clearly, a nice match is seen between the reduced model and the full model. This verifies that the key elements of the migration mechanisms are kept in the reduced model.

Features of the reduced model (3.1)

In view of the above comparison of the reduced model (3.1) and the original (2.6) two natural questions are: (i) Can we take a step further in decoupling the fluid flow dynamics and cell migration dynamics by ignoring the extra pressure terms $\frac{\widetilde{\lambda}_c}{\widetilde{\lambda}_T}\nabla(\Delta P(\alpha_c)+\Lambda(C))$ in the momentum equation in $(3.1)_2$, as commented in Remark 1? (ii) Is it possible to speed up the computation efficiency of the reduced model? To answer the first question, we run the reduced model where we have ignored this term, i.e., setting $\Delta P$ and $\Lambda$ as 0. A comparion is shown in Figure 7(A1, B1). Clearly, the results are quite similar. The reason can be explained by the fact that the coefficient ($\frac{\hat{\lambda}_c}{\hat{\lambda}_T}$, cell fraction flow) takes a rather small value due to the low mobility of cancer cells. Consequently, for the current fluid-cell-matrix "flow regime" it makes perfect sense to remove the impact from cancer cells on the fluid through these pressure-related terms. It might be tempting to suggest a further decoupling between fluid and cells by replacing $\nabla\cdot (\alpha_w\mathbf{u}) = Q$ by $\nabla\cdot \mathbf{u} = Q$ as indicated in Remark 1. Numerical investigations indicate that this completely changes the tumor progression behavior. This is quite reasonable since the volume fraction of cancer cells $\alpha_c$ typically can be relatively high at the tumor periphery where the fluid velocity $u_w$ is strongest, and the presence of cancer cell phase cannot be ignored.

Finally, what about point (ii)? We want to check the possibility of accelerating the computations of the reduced model by decreasing the number of updates of $u_w$ and $P_w$ in Step 2 (Section 3.3) of the numerical algorithm. Based on our tests, we find that at the early stage of the simulation, the cell-fluid dynamics is quite fast. So, at that stage it is necessary with a tight coupling, as described by the algorithm in Section 3.3. However, after a time $t = 10$, when the fluid and pressure field have been established, we stop updating $u_w$ and $P_w$. The result of this is shown in Figure 7(A2, B2). Clearly, it works perfectly well to strongly reduce the updates of $u_w$ and $P_w$ (after an initial transient time period till $t = 10$). We can see hardly no difference between this simplified calculation of $u_w$ and $P_w$ and the scheme where $u_w$ and $P_w$ are updated for every $\Delta t$.

5.   Concluding remarks

The main findings from our investigations in a one-dimensional tumor setting are:

(i) We have verified that a novel discretization of the general biphasic cell-fluid-matrix model (2.6) reliably can deal with the competing migration mechanisms, one that triggers downstream chemotactic migration, another which gives rise to a mechanical driven upstream. In particular, we find that these two opposite migration mechanisms, when properly tuned, do not cancel out each other but instead can be a means to form isolated groups of cancer cells that are situated in the region outside the primary tumor and close to lymphatic vessels. This detachment mechanism reflects an aggressive behavior that might be linked to lymp node metastases where cancer cells are able to escape through lymphatic vessels ^[4,11,38].

(ii) The cancer cells are responsible for formation of a leaky and aberrant intratumoral vasculature modeled through the source term $Q_v$ as well as remodelling of the matrix structure and density, which affect the tissue conductivity as represented through the fluid-matrix resistance force coefficient $\hat{\zeta}_w$. Hence, implicitly the cancer cells affect the fluid velocity field $\mathbf{u}_w$ and fluid pressure $P_w$, as expressed by the pressure equation (2.14). Here we note that total mobility $\hat{\lambda}_T$ (2.12) is dominated by fluid mobility $\hat{\lambda}_w$ in the presence of both phases since from (2.11) it follows that $\frac{\hat{\lambda}_c}{\hat{\lambda}_w}\sim \frac{ \alpha_c^2 \hat{\zeta}_w}{ \alpha_w^2\hat{\zeta_c}}\sim \frac{I_w}{I_c} \sim 0.005$ in view of (4.8). In the current version of the model the vasculature and tissue conductivity through $Q_v$ and $\hat{\zeta}_w$ are assumed to be stationary. By taking advantage of the fact that explicit presence of cancer cells in the model (2.6) has a weaker impact on the fluid, we have identified a chemotaxis-Navier-Stokes version (3.1) which bears similarity to the much studied cell-fluid model for cell swimming (3.3) ^[29]. It has been verified that the behavior of the full model (2.6) largely carries over to the reduced model (3.1).

(iii) As a byproduct we have found that the numerical discretization we employ can be greatly accelerated by taking advantage of the fact that the fluid pressure and velocity possibly change slowly as compared to the cell migration.

What is the potential importance of the findings presented in this work? The usefulness of the proposed numerical scheme for (2.6) and (3.1) lies in the the fact that it can naturally be extended to include compressible phases and viscous terms as well as more cell types that also are involved. For example, other cells like stromal cells and immune cells play important roles, depending on the special application of the model. In order to use the model for a more real-life situation where one has access to data through MRI images, etc, it will necessary to study the model in 2D and 3D setting. Methods used to relate the model to available data require many simulations. Hence, the possibility to accelerate computations is essential. The reduced model and its discrete counterpart opens for such applications. Combination with novel machine learning based frameworks ^[14] can potentially lead to even more efficient calculations.

Lymph node metastasis is a main reason why cancer becomes a deadly disease ^[11]. That is, instead of just having a growing coherent tumor, the tumor evolves into a heterogeneous mix of aggressive cancer cell phenotypes that can make use of different migration mechanisms for detachment. The model we explore contains a specific mechanism for creating isolated groupes of cancer cells that are located in the peritumoral region. A natural question to raise in this context is: What are critical parameters that determine the relative strength of the various migration mechanisms involved in the model (2.6) and its simplified version (3.1)? As explored in Section 4, whether the model gives rise to aggressive behavior or not depends on a fine-tuned balance between the magnitude of different forces that are involved, as expressed through the momentum balance equations. Mathematical analysis can detect finer details pertaining to such issues, insight that may not easily be detected from numerical simulations. Long time behavior and asymptotic behavior represents particular areas of interest. Singular behavior (in a broad sense), e.g., formation of isolated islands that are able to persist over long time, are of interest. A simplified version of a model similar to (2.5) has been investigated in a one-dimensional setting in ^[9]. Similarly, the competing downstream and upstream migration were explored by mathematical methods in ^[10] and precise statements of the long-time behavior were detected under appropriate constraints. An interesting question is whether such results can be taken over to more general models as discussed in this manuscript.

The chemotaxis-Navier-Stokes model (3.3) has attracted extensive interest the last two decades. A main question is to what extent various results known to hold for a pure chemotaxis model still hold under the influence of a fluid velocity field described by the (Navier-)Stokes model ^[34] influenced by the cell population through buoyancy. This include results on global solvability and blow-up behavior (critical parameters), e.g., ^{[1,3,6,13,27,28,35,36]} and references therein, as well as long-time behavior ^[33]. An interesting question is whether the model includes the possibility to represent spatial-dependent patterns over long time ^[29]. The long-time result in ^[33] ensures in a precise sense that the cell population will converge to a constant state without spatial variations. In the context of the tumor model, this amounts to a situation where isolated groups of cancer cells cannot persist over time. In that light an interesting question is whether it is possible to build on the experience with the analysis of the Tuval model to obtain qualitative insight into the solutions generated by the tumor model (3.1) or an appropriate modification/simplification of it.

Acknowledgments

The first author thanks the University of Stavanger for funding of the postdoc position.

Conflict of interest

The authors declare no conflict of interest.