Deformations of acid-mediated invasive tumors in a model with Allee effect

Abstract

We consider a Gatenby-Gawlinski-type model of invasive tumors in the presence of an Allee effect. We describe the construction of bistable one-dimensional traveling fronts using singular perturbation techniques in different parameter regimes corresponding to tumor interfaces with, or without, an acellular gap. By extending the front as a planar interface, we perform a stability analysis to long wavelength perturbations transverse to the direction of front propagation and derive a simple stability criterion for the front in two spatial dimensions. In particular we find that in general the presence of the acellular gap indicates transversal instability of the associated planar front, which can lead to complex interfacial dynamics such as the development of finger-like protrusions and/or different invasion speeds.

Keywords: 34D15; 35B35; 35B36; 35C07; 92C50.

Fig. 1

Results of a direct numerical simulation of (1.1) exhibiting the transversal long wavelength instability for parameter values taken from Gatenby and Gawlinski (1996) $(a,\kappa ,{\delta }_1,{\delta }_2,{\delta }_3,\rho ,\epsilon )=(0.35,0.1,12.5,0.1,70.0,1.0,0.0063)$. In the left panel, the underlying (longitudinally stable) one-dimensional traveling wave is shown (in which the normal cell density U is plotted in blue, the tumor cell density V in red, and the acid concentration W in yellow. The six right panels depict the results of two-dimensional simulations. The initial conditions of the two-dimensional run were constructed by trivially extending this 1D-stable profile in the y-direction and adding a small amount of noise. The simulations were performed in a co-moving frame $\xi =x+ct$; in the laboratory frame the tumor would travel with speed $c=0.0401$ to the left, corresponding to the wave speed of the initial front profile (so that in the absence of any instability, the fronts would appear stationary). The corresponding tumor U-profiles (top three panels) and V-profiles (bottom three panels) are plotted (for $t=20000,30000,40000$ from left to right); for all profiles yellow indicates high density of cells, and blue depicts low density of cells

Fig. 2

Bistable traveling front solutions in the benign (left) and malignant (right) cases. Note that in the malignant case, an acellular gap may appear in between the interfaces formed by the U and V profiles and that this case corresponds to the situation of Fig. 1

Fig. 3

Schematic of the singular heteroclinic orbit in the benign case where $u\ne 0$ after the invasion of the traveling wave (left) and the malignant case where $u=0$ after the invasion (right)

Fig. 4

(Left) The fast orbit $(v_1,q_1)(\xi ;w)$ in the subspace $u=1-{\delta }_{1}w$ with $0\le w<min\{\rho {(1-a)}^2/(4{\delta }_2),1/{\delta }_1\}$. (Right) The fast orbit $(v_0,q_0)(\xi ;w)$ in the subspace $u=0$ with $w>0$

Fig. 5

The slow orbits on ${\mathcal{M}}_{0,1}^0$ and ${\mathcal{M}}_{0,1}^{+}$, corresponding to ${\mathcal{W}}^{\text{u}}(0,0)$ and ${\mathcal{W}}^{\text{u}}(V^{+},0)$, respectively, in the benign (left), malignant no-gap (center), and malignant gap (right) cases. The manifolds ${\mathcal{W}}^{\text{s,u}}(0,0)$ on ${\mathcal{M}}_{0,1}^0$ are depicted in solid red, while the manifolds ${\mathcal{W}}^{\text{s,u}}(V^{+},0)$ on ${\mathcal{M}}_{0,1}^{+}$ are depicted in dashed red, and the vertical dashed line indicates the subspace $w={\delta }_1^{-1}$ at which the reduced dynamics transition from ${\mathcal{M}}_0^{0,+}$ to ${\mathcal{M}}_1^{0,+}$; see also Fig. 6

Fig. 6

Structure of the singular heteroclinic orbit in the benign (left), malignant no-gap (middle), and malignant gap (right) cases. The subspace $\{u=0\}$ is depicted in purple, while the subspace $\{u=1-{\delta }_{1}w\}$ is depicted in blue

Fig. 7

Plot of gap width versus ${\delta }_1\in (0.1,15)$ obtained by numerical continuation of the traveling wave equation (2.2) for the parameter values $(a,\kappa ,{\delta }_2,{\delta }_3,\rho )=(0.1,0.1,0.1,70,1.0)$ and $\epsilon =0.0063$ (red), $\epsilon ={10}^{-4}$ (yellow). The gap width was computed by measuring the spatial width where both the u and v profiles of the corresponding front solution were below a threshold value of $10\epsilon$. Also plotted (blue) is the singular limit gap width obtained by solving (2.7) for $w_{\ast }^{}$ using Mathematica and integrating (2.5) to obtain the time spent along ${\mathcal{M}}_0^0$ between $w={\delta }_1^{-1}$ and $w=w_{\ast }^{}$

Fig. 8

1D traveling wave profiles obtained for the parameter values $(a,\kappa ,{\delta }_1,{\delta }_2,{\delta }_3,\rho ,\epsilon )=(0.1,0.1,12.5,0.1,70,1,0.0063)$ (first row), $(a,\kappa ,{\delta }_1,{\delta }_2,{\delta }_3,\rho ,\epsilon )=(0.25,0.1,12.5,0.1,70,1,0.0063)$ (second row), and $(a,\kappa ,{\delta }_1,{\delta }_2,{\delta }_3,\rho ,\epsilon )=(0.35,0.1,12.5,0.1,70,1,0.0063)$ (third row), $(a,\kappa ,{\delta }_1,{\delta }_2,{\delta }_3,\rho ,\epsilon )=(0.25,0.05,11.5,3,1,15,0.05)$ (fourth row). The u, v, w profiles are plotted in blue, red, yellow, respectively. Profiles were obtained by solving the traveling wave equation (2.1) in MATLAB. Also shown are the 1D spectra, providing numerical evidence that all four solutions are 1D-stable, as well as a continuation of the critical eigenvalue ${\lambda }_{\text{c}}(\ell )$ for small, positive values of the wavenumber $\ell$

Fig. 9

Results of numerical continuation in AUTO07p (Doedel et al. 2007) for the parameter values $(a,\kappa ,{\delta }_2,{\delta }_3,\rho ,\epsilon )=(0.1,0.1,0.1,70,1.0,0.0063)$ for values of ${\delta }_1\in (0.05,15)$: wave speed c versus ${\delta }_1$ (left), ${\lambda }_{\text{c},2}$ versus ${\delta }_1$ (right)

Fig. 10

Results of numerical continuation for the parameter values $(\kappa ,{\delta }_2,{\delta }_3,\rho ,\epsilon )=(0.1,0.1,70,1.0,$0.0063) for a range of a-values for ${\delta }_1=0.6$ (blue) and ${\delta }_1=12.5$ (red): wave speed c versus a (left), ${\lambda }_{\text{c},2}$ versus a (right)

Fig. 11

Results of numerical continuation of the curve ${\lambda }_{\text{c},2}=0$ for the parameter values $(a,\kappa ,{\delta }_3,\rho )=(0.1,0.1,70,1.0)$ in the $({\delta }_1,{\delta }_2)$-plane for $\epsilon =0.0063$ (yellow) and $\epsilon ={10}^{-5}$ (blue). Below each curve we have ${\lambda }_{\text{c},2}<0$, while ${\lambda }_{\text{c},2}>0$ above. Plotted in red is the malignant/benign boundary ${\delta }_1{V}_{+}=1$

Fig. 12

Results of direct numerical simulations exhibiting different manifestations of the long wavelength instability. The simulations were performed in a co-moving frame corresponding to the wave speed of the initial front profile (so that in the absence of any instability, the fronts would appear stationary), with Neumann boundary conditions in $\xi$, and periodic boundary conditions in y. Finite differences were used for spatial discretization, and MATLAB’s ode15s routine was used for time stepping. The corresponding v-profile is depicted, where yellow indicates high density of tumor cells, and blue depicts low density of tumor cells. (First row) Simulation for the parameter values $(a,\kappa ,{\delta }_1,{\delta }_2,{\delta }_3,\rho ,\epsilon )=(0.25,0.1,12.5,0.1,70.0,1.0,0.0063)$ at the times $t=5000,10000,15000,20000$ from left to right, with wave speed $c=0.2211$. (Second row) Simulation for the parameter values $(a,\kappa ,{\delta }_1,{\delta }_2,{\delta }_3,\rho ,\epsilon )=(0.35,0.1,12.5,0.1,70.0,1.0,0.0063)$ at the times $t=10000,20000,30000,40000$ from left to right, with wave speed $c=0.0401$. (Third row) Simulation for the parameter values $(a,\kappa ,{\delta }_1,{\delta }_2,{\delta }_3,\rho ,\epsilon )=(0.25,0.05,11.5,3,1,15,0.05)$ at the times $t=250,350,450,550$ from left to right, with wave speed $c=0.3296$. In each of the three simulations, initial conditions were constructed by trivially extending the 1D-stable profiles from Fig. 8 (rows 2 through 4) in the y-direction and adding a small amount of positive noise