Tumor boundary instability induced by nutrient consumption and supply

Yu Feng Yu Feng, Beijing International Center for Mathematical Research, Peking University, No. 5 Yiheyuan Road Haidian District, Beijing, P.R.China 100871 [email protected] , Min Tang Min Tang: School of Mathematical Sciences, Institute of Natural Sciences, MOE-LSC, Shanghai Jiao Tong University, Shanghai, 200240, P. R. China [email protected] , Xiaoqian Xu Xiaoqian Xu: Zu Chongzhi Center for Mathematics and Computational Sciences, Duke Kunshan University, China [email protected] and Zhennan Zhou Zhennan Zhou, Beijing International Center for Mathematical Research, Peking University, No. 5 Yiheyuan Road Haidian District, Beijing, P.R.China 100871 [email protected]

Abstract.

We investigate the tumor boundary instability induced by nutrient consumption and supply based on a Hele-Shaw model derived from taking the incompressible limit of a cell density model. We analyze the boundary stability/instability in two scenarios: 1) the front of the traveling wave; 2) the radially symmetric boundary. In each scenario, we investigate the boundary behaviors under two different nutrient supply regimes, in vitro and in vivo. Our main conclusion is that for either scenario, the in vitro regime always stabilizes the tumor’s boundary regardless of the nutrient consumption rate. However, boundary instability may occur when the tumor cells aggressively consume nutrients, and the nutrient supply is governed by the in vivo regime.

1. Introduction

Tumor, one of the major diseases threatening human life and health, has been widely concerned. The mathematical study of tumors has a long history and constantly active. We refer the reader to the textbook [cristini2010multiscale, cristini2017introduction] and review articles [araujo2004history, lowengrub2009nonlinear, byrne2006modelling, roose2007mathematical]. Previous studies and experiments indicate that the shape of tumors is one of the critical criteria to distinguish malignant from benign. Specifically, malignant tumors are more likely to form dendritic structures than benign ones. Therefore, it is significant to detect and predict the formation of tumor boundary instability through mathematical models. Before discussing the mathematical studies of tumor morphology, we review relevant mathematical models as follows.

The first class of model was initiated by Greenspan in 1976 [greenspan1976growth], which further inspired a mass of mathematical studies on tumor growth (e.g., [byrne1996growth, chaplain1996avascular, zheng2005nonlinear, friedman2007bifurcationA]). The tumor is regarded as an incompressible fluid satisfying mass conversation. More precisely, these free boundary type models have two main ingredients. One is the nutrient concentration $\sigma$ governed by a reaction-diffusion equation, which considers the consumption by the cells and the supplement by vessels. The other main component is the internal pressure $p$ , which further induces the cell velocity $v$ via different physical laws (e.g., Darcy’s law [greenspan1976growth, byrne1996growth, friedman1999analysis, cristini2003nonlinear], Stokes law [friedman2007bifurcationA, friedman2002quasistatic, friedman2002quasi], and Darcy&Stokes law [franks2003interactions, franks2009interactions, zheng2005nonlinear, king2006mathematical, pham2011predictions]). Finally, the two ingredients are coupled via the mass conservation of incompressible tumor cells, which yields the relation $\nabla\cdot v=\lambda(\sigma)$ , with the cell proliferation rate $\lambda$ depending on $\sigma$ . To close the model, the Laplace-Young condition ( $p|_{\partial\Omega}=\gamma\kappa$ , where $\kappa$ is the mean curvature, and $\gamma$ stands for the surface tension coefficient) is imposed on the tumor-host interface. For some variant models, people replace the Laplace-Young condition with other curvature-dependent boundary conditions (see, e.g., [turian2019morphological, lu2019nonlinear, pham2018nonlinear]). More sophisticated models were also investigated recently. In particular, we mention the studies based on the two-phase models [pham2018nonlinear, turian2019morphological, lu2019nonlinear], and the works involve chemotaxis [he2022incompressible, kim2022density].

Most studies on the stability/instability of tumor boundary are based on the above class of models and have been investigated from different points of view. Among them, for different models (e.g., Darcy [fontelos2003symmetry, friedman2001symmetry, friedman2006bifurcation, friedman2006asymptotic, friedman2008stability]; and Stokes [friedman2006free, friedman2007bifurcationA, friedman2007bifurcationB]), Friedman et al. proved the existence of non-radially symmetric steady states analytically and classified the stability/instability of the boundaries from the Hopf bifurcation point of view. Specifically, in their studies, the bifurcation parameter is characterized by the cell proliferation rate or ratio to cell-cell adhesiveness. Then the authors showed that the boundary stability/instability changes when the parameter crosses a specific bifurcation point. On the other hand, Cristini et al. in [cristini2003nonlinear], as the pioneers, employ asymptotic analysis to study and predict the tumor evolution. Their work is of great significance to the dynamic simulation of tumors and nurtured more related works in this direction [macklin2007nonlinear, pham2018nonlinear, turian2019morphological, lu2020complex, lu2022nonlinear]. All the research demonstrated that many factors could induce the tumor’s boundary instability, including but not limited to vascularization [cristini2003nonlinear, pham2018nonlinear, lu2019nonlinear, lu2022nonlinear], proliferation [friedman2001symmetry, friedman2008stability, friedman2006free, friedman2007bifurcationB, cristini2003nonlinear, lu2022nonlinear], apoptosis [friedman2001symmetry, friedman2008stability, friedman2006free, friedman2007bifurcationB, cristini2003nonlinear, turian2019morphological, lu2022nonlinear, pham2018nonlinear, lu2022nonlinear], cell-cell adhesion [friedman2001symmetry, friedman2008stability, friedman2006free, friedman2007bifurcationB, cristini2003nonlinear, pham2018nonlinear], bending rigidity [turian2019morphological, lu2019nonlinear], microenvironment [turian2019morphological, pham2018nonlinear, macklin2007nonlinear, lu2020complex], chemotaxis [lu2022nonlinear, lu2020complex].

In recent decades, tumor modeling from different perspectives has emerged and developed. In particular, one could consider the density model proposed by Byrne and Drasdo in [byrne2009individual], in which the tumor cell density $\rho$ is governed by a porous medium type equation, and the internal pressure $p$ is induced by the power rule $p=\rho^{m}$ with the parameter $m>1$ . The power rule enables $p$ naturally vanish on the tumor boundary. Moreover, the boundary velocity $v$ is governed by Darcy’s law $v=-\nabla p|_{\partial\Omega}$ . Previous research indicates that the porous media type equations have an asymptote concerning the parameter m tending to infinity [aronson1998limit, gil2001convergence, igbida2002mesa, kim2003uniqueness, kim2009homogenization]. Motivated by this, Perthame et al. derived the second kind of free boundary model in [perthame2014hele] by taking the incompressible limit (sending $m$ to infinity), or equivalently mesa-limit of the density models. An asymptotic preserving numerical scheme was designed by J.Liu et al. in [liu2018accurate], the scheme naturally connects the numerical solutions to the density models to that of the free boundary models.

In the mesa-limit free boundary models proposed in [perthame2014hele], the limit density $\rho_{\infty}$ can only take value in $[0,1]$ , and the corresponding limit pressure $p_{\infty}$ is characterized by a monotone Hele-Shaw graph. More specifically, $p_{\infty}$ vanishes on the unsaturated region where $\rho_{\infty}<1$ (see equation (2.7)). The Hele-Shaw graph representation of pressure brings the following advantages. Firstly, in the Hele-Shaw type model, the formation of a necrotic core can be described by an obstacle problem [guillen2022hele], which leads $\rho_{\infty}$ to decay exponentially in the necrotic core. Due to the Hele-Shaw graph, the pressure $p_{\infty}$ naturally vanishes there instead of taking negative values. Secondly, a transparent regime called "patch solutions" exists, in which $\rho_{\infty}$ remains in the form of $\chi_{D(t)}$ , i.e., the indicator function of the tumor region. Again, to satisfy the corresponding Hele-Shaw graph, $p_{\infty}$ has to vanish on the tumor’s interface (where $\rho_{\infty}$ drops from $1$ to $0$ ), which is significantly different from the first kind of free boundary models developed from [greenspan1976growth], in which the internal pressure relies on the boundary curvature $\kappa$ as mentioned previously. Moreover, in the mesa-limit free boundary models, the boundary velocity is still induced by Darcy’s law $v_{\infty}=-\nabla p_{\infty}|_{\partial\Omega}$ . For completeness, the derivation of the mesa-limit model is summarized in Section 2.1. Albeit various successful explorations based on such mesa-limit free boundary models [mellet2017hele, david2021free, guillen2022hele, perthame2015incompressible, david2021convergence, kim2003uniqueness, kim2016free, kim2018porous, kim2022density, kim2019singular, david2021incompressible, perthame2014traveling, jacobs2022tumor, liu2021toward, dou2022tumor], the study on its boundary stability/instability is yet thoroughly open.

The primary purpose of this paper is to investigate whether boundary instability will arise in the mesa-limit free boundary models, which should shed light on the boundary stability of the cell density models when $m$ is sufficiently large. To simplify the discussion, we consider tumors in the avascular stage with saturated cell density so that the density function $\rho_{\infty}$ is a patch solution, and the tumor has a sharp interface. As the first attempt in this regard, we explore the instability caused by nutrient consumption and supply. A similar mechanism can induce boundary instability in other biological systems, see [ben1994generic] for nutrient induce boundary instability in bacterial colony growth models. The role of nutrition in tumor models has been widely studied, and we refer the reader to the latest article in this direction [jacobs2022tumor]. Inspired by [perthame2014traveling], we divide the nutrient models into two kinds, in vitro and in vivo, according to the nutrient supply regime. In either regime, the nutrient is consumed linearly in the tumor region with a rate $\lambda>0$ . However, in the in vitro model, we assume that a liquid surrounds the tumor with nutrient concentration $c_{B}$ . Mathematically, the nutrient concentration remains $c_{B}$ at the tumor-host interface. For the in vivo model, the nutrient is transported by vessels outside the tumor and reaches $c_{B}$ at the far field. Correspondingly, we assume the exchange rate outside the tumor is determined by the concentration difference from the background, i.e., $c_{B}-c$ . The two nutrient models will be specified more clearly in Section 2.1.2.

Our study of boundary stability/instability consists of two scenarios. We begin with a relatively simple case, the front of traveling waves, in which quantitative properties can be studied more explicitly. In this case, the unperturbed tumor region corresponds to a half plane with the boundary being a vertical line propagating with a constant normal velocity. Then we test the boundary stability/instability by adding a perturbation with frequency $l\in\mathbb{R}^{+}$ and amplitude $\delta$ . Our analysis shows that in the in vitro regime, $\delta$ always decreases to zero as time propagates. In other words, the boundary is stable for any frequency perturbation. In contrast, in vivo regime, there exists a threshold value $L$ such that the perturbation with a frequency smaller than $L$ becomes unstable when the nutrient consumption rate, $\lambda$ , is larger than one.

The above case corresponds to the boundary stability/instability while the tumor is infinitely large. In order to further explore the influence of the finite size effect on the boundary stability/instability, we consider the perturbation of radially symmetric boundary with different wave numbers $l\in\mathbb{N}$ and radius $R$ . Our analysis shows that the in vitro regime still suppresses the increase of perturbation amplitude and stabilizes the boundary regardless of the consumption rate, perturbation wave number, and tumor size. For the in vivo regime, when the consumption rate $\lambda$ is less than or equal to one, the boundary behaves identically the same as the in vitro case. However, when $\lambda$ is greater than one, the continuous growth of tumor radius will enable perturbation wave number to become unstable in turn (from low to high). Further more, as $R$ is approaching infinity, the results in the radial case connect to the counterparts in the traveling wave case.

The main contribution of this work is to show that tumor boundary instability can be induced by nutrient consumption and supply. As a by-product, our results indicate that the cell apoptosis and curvature-dependent boundary conditions present abundantly in previous studies (e.g., [cristini2003nonlinear, friedman2007bifurcationA]) are unnecessary for tumor boundary instability formation.

The paper is organized as follows. In Section 2, we first derive our free boundary models by taking the incompressible limit of density models characterized by porous medium type equations in Section 2.1. Besides that, we also introduce the in vitro and in vivo nutrient regimes in this subsection. Furthermore, the corresponding analytic solutions are derived in Section 2.2. Section 3 is devoted to introducing the linear perturbation technique in a general framework. Then, by using the technique in Section 3, we study the boundary stability of the traveling wave and the radially symmetric boundary under the two nutrient regimes, respectively, in Section 4 and Section 5 (with main results in Section 4.1 and Section 5.1). Finally, we summarize our results and discuss future research plans in Section LABEL:Conclusion.

2. Preliminary

2.1. model introduction

2.1.1. The cell density model and its Hele-Shaw limit

To study the tumor growth under nutrient supply, let $\rho(x,t)$ denote the cell population density and $c(x,t)$ be the nutrient concentration. We assume the production rate of tumor cells is given by the growth function $G(c)$ , which only depends on the nutrient concentration. On the other hand, we introduce

(2.1)

D(t)=\left\{\rho(x,t)>0\right\}

to denote the support of $\rho$ . Physically, it presents the tumoral region at time $t$ . We assume the tumoral region expands with a finite speed governed by the Darcy law $v=-\nabla p$ via the pressure $p(\rho)=\rho^{m}$ . Thus, the cell density $\rho$ satisfies the equation:

(2.2)

\frac{\partial}{\partial t}\rho-\nabla\cdot\left(\rho\nabla p(\rho)\right)=\rho G(c),\quad x\in\mathbb{R}^{2},\quad t\geqslant 0.

For the growth function $G(c)$ , we assume

(2.3)

G(c)=G_{0}c,\quad\text{with }\quad G_{0}>0,

note that in contrast to the nutrient models in [tang2014composite, perthame2014hele], we eliminate the possibility of the formation of a necrotic core by assuming that $G(\cdot)$ is always positive and linear (for simplicity), since this project aims to study the boundary instability induced by the nutrient distribution itself.

Many researches, e.g. [perthame2014hele, david2021free, david2021incompressible, kim2016free, kim2018porous, guillen2022hele], indicate that there is a limit as $m\rightarrow\infty$ which turns out to be a solution to a free boundary problem of Hele-Shaw type. To see what happens, we multiply equation (2.2) by $m\rho^{m-1}$ on both sides to get

(2.4)

\frac{\partial}{\partial t}p(\rho)=\lvert\nabla p(\rho)\rvert^{2}+mp(\rho)\Delta p(\rho)+mG_{0}p(\rho)c.

Hence, if we send $m\rightarrow\infty$ , we formally obtain the so called complementarity condition (see [perthame2014hele, david2021free] for a slight different model):

(2.5)

p_{\infty}(\Delta p_{\infty}+G_{0}c)=0.

On the other hand, the cell density $\rho(x,t)$ converges to the weak solution (see [perthame2014hele]) of

(2.6)

\frac{\partial}{\partial t}\rho_{\infty}-\nabla\cdot\left(\rho_{\infty}\nabla p_{\infty}\right)=\rho_{\infty}G(c),

and $p_{\infty}$ compels the limit density $\rho_{\infty}$ only take value in the range of $\left[0,1\right]$ for any initial date $\rho_{0}\in\left[0,1\right]$ (see Theorem 4.1 in [perthame2014hele] for a slightly different model). Moreover, the limit pressure $p_{\infty}$ belongs to the Hele-Shaw monotone graph:

(2.7)

p_{\infty}(\rho_{\infty})=\left\{\begin{array}[]{rcr}0,\qquad 0\leqslant\rho_{\infty}<1,\\ \left[0,\infty\right),\qquad\rho_{\infty}=1.\\ \end{array}\right.

The incompressible limit and the complementarity condition of a fluid mechanical related model have been rigorously justified in [perthame2014hele, david2021free]. And the incompressible limit of (2.2) (coupled with nutrient models that will be introduced in the next section) was verified numerically in [liu2018analysis].

We define the support of $p_{\infty}$ to be

(2.8)

D_{\infty}(t)=\left\{p_{\infty}(x,t)>0\right\},

then (2.5) and (2.7) together yield


(2.9a)	$\displaystyle-\Delta p_{\infty}=G_{0}c\qquad$	$\displaystyle\text{for}\quad x\in D_{\infty}(t),$
(2.9b)	$\displaystyle p_{\infty}=0,\qquad$	$\displaystyle\text{for}\quad x\in\mathbb{R}^{2}\setminus D_{\infty}(t),$
and $\rho_{\infty}=1$ in $D_{\infty}$ .

Therefore, once the nutrient concentration $c(x,t)$ is known one can recover $p_{\infty}$ from the elliptic equation above.

Now we justify the relationship between $D(t)$ and $D_{\infty}(t)$ . Observe that when $m$ is finite, $\rho$ and $p(\rho)$ have the same support $D(t)$ , whereas as $m$ tends to infinity, $\rho_{\infty}$ may have larger support than $p_{\infty}$ . However, a large class of initial data, see e.g. [perthame2016some], enable the free boundary problem (correspond to (2.6) and (2.9)) possess patch solutions, i.e., $\rho_{\infty}=\chi_{D_{\infty}}$ , where $\chi_{A}$ stands for the indicator function of the set $A$ . In this case, the support of $p_{\infty}$ coincides with that of $\rho_{\infty}$ . Moreover, the boundary velocity $v$ is governed by Darcy law $v=-\nabla p_{\infty}$ . Further, the boundary moving speed along the normal direction at the boundary point $x$ , denote by $\sigma(x)$ , is given by:

(2.10)

\sigma(x)=-\nabla p_{\infty}\cdot\hat{n}(x),

where $\hat{n}(x)$ is the outer unit normal vector at $x\in\partial D_{\infty}(t)$ . The boundary speed for more general initial data was studied in [kim2018porous].

As the end of this subsection, we emphasize that in our free boundary model, as the limit of the density models, the pressure $p_{\infty}$ always vanishes on $\partial D_{\infty}$ . However, as mentioned previously, in the first kind free boundary models, the internal pressure $\tilde{p}$ is assumed to satisfy the so-called Laplace-Young condition (or some other curvature dependent boundary condition). Mathematically, the boundary condition (2.9b) is replaced by

(2.11)

\tilde{p}(x)=\gamma\kappa(x),

where $\gamma>0$ is a constant coefficient, and $\kappa(x)$ denotes the curvature at the boundary point $x$ . In the related studies, the curvature condition (2.11) plays an essential role (e.g., [cristini2003nonlinear, friedman2007bifurcationA]).

2.1.2. Two nutrient models

Regarding the nutrient, it diffuses freely over the two dimensional plane. However, inside the tumoral region, the cells consume the nutrient. While outside the tumor, the nutrient exchanges with the far field concentration $c_{B}$ provided by the surrounding environment or vasculature. It follows that the following reaction-diffusion equation can govern the consumption, exchange, and diffusion of the nutrient in general:

(2.12)

\tau\partial_{t}c-\Delta c+\Psi(\rho,c)=0,

where $\tau$ is the characteristic time scale of the nutrient change, and $\Psi(\rho,c)$ describes the overall effects of the nutrient supply regime outside the tumor and the nutrient consumption by cells inside the tumor. To simplify the mathematical analysis, we drop the time derivative in (2.12) and consider following elliptic formulation instead

(2.13)

-\Delta c+\Psi(\rho,c)=0.

This is reasonable because $\tau\ll 1$ (see, e.g., [greenspan1972models, adam2012survey, byrne1996growth]). As in [perthame2014traveling], two specific developed and widely studied models are the in vitro and the in vivo model.

For the in vitro model, we assume that the tumor is surrounded by a liquid in which the exchange rate with the background is so fast that the nutrient concentration can be assumed to be the same constant $c_{B}$ as that of the surrounding liquid, while inside the tumoral region, the consumption function is bi-linear in both the concentration $c$ and the cell density $\rho$ with consumption rate $\lambda>0$ . The boundary instability was observed in models where tissues aggressively consume nutrients [maury2014congestion]. Therefore, in our models, we expect boundaries are more likely to be unstable when $\lambda$ is large. When the in vitro is coupled with the cell density model (2.2), equation (2.13) writes


(2.14a)	$\displaystyle-\Delta c+\lambda\rho c=0,\qquad$	$\displaystyle\text{for}\quad x\in D(t),$
(2.14b)	$\displaystyle c=c_{B},\qquad$	$\displaystyle\text{for}\quad x\in\mathbb{R}^{2}\setminus D(t).$

By considering the incompressible limit of the density model (sending $m\rightarrow\infty$ ), and concern patch solutions $\rho_{\infty}=\chi_{D_{\infty}}$ . Equation (2.14) tends to:


(2.15a)	$\displaystyle-\Delta c+\lambda c=0,\qquad$	$\displaystyle\text{for}\quad x\in D_{\infty}(t),$
(2.15b)	$\displaystyle c=c_{B},\qquad$	$\displaystyle\text{for}\quad x\in\mathbb{R}^{2}\setminus D_{\infty}(t).$

For the in vivo model, the consumption of nutrients within the tumor region (where $\rho>0$ ) remains the same as in the in vitro model. However, in the in vivo model, the nutrients are provided by vessels of the healthy tissue surrounding the tumor, while the healthy tissue consumes nutrients as well. This leads to the nutrient supply outside the tumor being determined by the concentration difference from the background, $c_{B}-c$ , with a positive coefficient $\tilde{\lambda}$ . Mathematically, the overall function $\Psi(\rho,c)$ is written as $\Psi(\rho,c)=\lambda\rho c\cdot\chi_{D}-\tilde{\lambda}(c_{B}-c)\cdot\chi_{D^{c}}$ . For simplicity, we set $\tilde{\lambda}=1$ and $\lambda>0$ . Note that this expression guarantees the nutrient concentration reaches $c_{B}$ at the far field. A more detailed discussion of this issue can be found in [chatelain2011emergence, jagiella2012parameterization].

With the same reason as the previous case, by taking $m\rightarrow\infty$ in the density model and concerning patch solutions, we get the in vivo nutrient equations for the limit free boundary model,


(2.16a)	$\displaystyle-\Delta c+\lambda c=0,\qquad$	$\displaystyle\text{for}\quad x\in D_{\infty}(t),$
(2.16b)	$\displaystyle-\Delta c=c_{B}-c,\qquad$	$\displaystyle\text{for}\quad x\in\mathbb{R}^{2}\setminus D_{\infty}(t).$

Moreover, we need to emphasize that the in vivo we refer to is different from the previous articles (see, e.g., [cristini2003nonlinear]) in which in vivo corresponds to the vascularization inside the tumor.

The uneven growth phenomena in the tumor models are conjectured due to the non-uniform distribution of nutrients [maury2014congestion]. More precisely, in contrast to the fingertips region, the nutrient is inadequate around the valley since more cells consume nutrients there. Consequently, the tissue around the tips grows faster than the valleys, and therefore instability occurs. In the in vitro model, the concentration of the nutrient will match the background concentration $c_{B}$ at the boundary regardless of the regions. However, for the in vivo model, the nutrient is directly available only from healthy tissue; this regime will enlarge the concentration difference at the tips and valleys. Therefore, we expect tumor borders are more prone to grow unevenly in the in vivo models, in particular when the consumption rate $\lambda$ is relatively large.

2.2. Analytic solutions

Starting from this section, we focus on the mesa limit free boundary models. Therefore, for simplicity of the notations, we drop the free boundary models’ subscripts and use $D(t)$ , $\rho$ , and $p$ to denote the tumoral region, cell density, and pressure in the limit model. On the other hand, through this paper, we use $I_{j}$ and $K_{j}$ ( $j\in\mathbb{N}$ ) to denote the second kind of modified Bessel functions, their definitions and basic properties are reviewed in Appendix LABEL:sec:Properties_of_Bessel_functions.

The models introduced in Section 2.1 have been studied in [liu2018analysis] when $\lambda=1$ . In particular, the authors derived $2D$ radially symmetric solutions for the free boundary models, which are coupled with either the in vitro or the in vivo model. Moreover, their computation yields that as the radius of the tumor tends to infinity, the boundary velocity tends to be a finite constant. In other words, the radially symmetric solutions converge to traveling wave solutions.

For self-consistency, we recall the derivation of the radially symmetric solutions in [liu2018analysis] in this section. Besides that we also derive the traveling wave solutions for the two nutrient models and verify that they are indeed the limit of the radially symmetric solutions as radius goes to infinity. The analytical solutions in this section will serve as the cornerstone of subsequent perturbation analysis. Now, we begin with the traveling wave scenario.

2.2.1. traveling plane solution for the in vitro model

For solving two-dimensional traveling wave solutions, we fix the traveling front at $\xi=x-\sigma t=0$ , where $\sigma$ stands for the traveling speed and will be determined later. Without loss of generality, let the tumoral region be the left half plane, that is $D(t)=\left\{(\xi,y)|\xi\leqslant 0\right\}$ . One can easily see that in the unperturbed two-dimensional problem, to find its solution reduces to solve a one-dimensional problem. Moreover, we disclose that the variable $y$ will serve as the perturbation parameter in the perturbation problems, which will be investigated later. The one dimensional problem writes:


(2.17a)	$\displaystyle-\partial^{2}_{\xi}c+\lambda c=0,\qquad$	$\displaystyle\text{for}\quad\xi\leqslant 0,$
(2.17b)	$\displaystyle c=c_{B},\qquad$	$\displaystyle\text{at}\quad\xi=0,$
in addition, we also assume the concentration of nutrient vanish at the center of tumor, that is
(2.17c)	$c(\xi)=0,\qquad\text{for}\quad\xi=-\infty.$

And the equations for pressure $p(\xi,y)$ , i.e., (2.9) and (2.10) reads


(2.18a)	$\displaystyle-\partial^{2}_{\xi}p=G_{0}c,\qquad$	$\displaystyle\text{for}\quad\xi\leqslant 0,$
(2.18b)	$\displaystyle p=0,\qquad$	$\displaystyle\text{for}\quad\xi\geqslant 0,$
and traveling speed is given by
(2.18c)	$\sigma=-\partial_{\xi}p(0).$
Since the gradient of the pressure is always equal to zero at the center of the tumor, we also have
(2.18d)	$\partial_{\xi}p(\xi)=0,\qquad\text{for}\quad\xi=-\infty.$

By solving (2.17) we get

(2.19)

c=c_{B}e^{\sqrt{\lambda}\xi},\qquad\text{for}\quad\xi\leqslant 0,

plug the above expression into (2.18) to solve for $p$ and get:

(2.20)

p(\xi)=-\frac{G_{0}c_{B}}{\lambda}e^{\sqrt{\lambda}\xi}+\frac{G_{0}c_{B}}{\lambda}\qquad\text{for}\quad\xi\leqslant 0.

Then, we can further find the traveling speed

(2.21)

\sigma=-\partial_{\xi}p(0)=\frac{G_{0}c_{B}}{\sqrt{\lambda}}.

2.2.2. traveling plane solution for the in vivo model

For the in vivo model, the only difference from the in vitro model is the equations for $c(\xi,y)$ are replaced by


(2.22a)	$\displaystyle-\partial^{2}_{\xi}c+\lambda c=0,\qquad$	$\displaystyle\text{for}\quad\xi\leqslant 0,$
(2.22b)	$\displaystyle-\partial^{2}_{\xi}c=c_{B}-c,\qquad$	$\displaystyle\text{for}\quad\xi\geqslant 0,$
(2.22c)	$\displaystyle c(\xi)=0,\qquad$	$\displaystyle\text{for}\quad\xi=-\infty.$
in addition, $c$ and $\partial_{\xi}c$ are both continuous at the boundary of the tumor, that is
(2.22d)	$c(0^{-})=c(0^{+})\quad\text{and}\quad\partial_{\xi}c(0^{-})=\partial_{\xi}c(0^{+}).$

And the pressure $p$ still satisfies (2.18).

By solving (2.22) we get

(2.23)

c(\xi)=\left\{\begin{array}[]{rcr}\frac{c_{B}}{\sqrt{\lambda}+1}e^{\sqrt{\lambda}\xi}\stackrel{{\scriptstyle\scriptscriptstyle\textup{def}}}{{=}}c^{\text{(i)}}(\xi)&\text{for}\quad\xi\leqslant 0,\\ -\frac{\sqrt{\lambda}c_{B}}{\sqrt{\lambda}+1}e^{-\xi}+c_{B}\stackrel{{\scriptstyle\scriptscriptstyle\textup{def}}}{{=}}c^{\text{(o)}}(\xi)&\text{for}\quad\xi\geqslant 0.\\ \end{array}\right.

and plug the above expression into (2.18) to derive $p$ and get

(2.24)

p(\xi)=-\frac{G_{0}c_{B}}{\lambda(\sqrt{\lambda}+1)}e^{\sqrt{\lambda}\xi}+\frac{G_{0}c_{B}}{\lambda(\sqrt{\lambda}+1)}\qquad\text{for}\quad\xi\leqslant 0.

And the boundary speed is given by

(2.25)

\sigma=-\partial_{\xi}p(0)=\frac{G_{0}c_{B}}{\lambda+\sqrt{\lambda}}.

By now, we have finished the derivation for the traveling wave solutions. In the next two subsections, we recall the derivation for the radially symmetric scenario in [liu2018analysis] and verify that the boundary speeds converge to the traveling waves’ for the corresponding nutrient regime.

2.2.3. 2D radially symmetric model for the in vitro model

For the radially symmetric free boundary model, the tumoral region becomes $D(t)=\mathbb{B}_{R(t)}(0)$ (a disk centered at origin with radius $R$ ). In this case, we employ polar coordinates $(r,\theta)$ , and we can conclude that the solutions are independent of $\theta$ by symmetry. However the variable $\theta$ will play an important role in the perturbed problem, which will be seen in the later sections. Thus, for the free boundary model with nutrients governed by the in vitro model, equation (2.15) can be further written as


(2.26a)	$\displaystyle-\frac{1}{r}\partial_{r}(r\partial_{r}c)+\lambda c=0,\qquad$	$\displaystyle\text{for}\quad r\leqslant R(t),$
(2.26b)	$\displaystyle c=c_{B},\qquad$	$\displaystyle\text{for}\quad r\geqslant R(t).$

And the equations for pressure $p$ (2.9) and (2.10) reads


(2.27a)	$\displaystyle-\frac{1}{r}\partial_{r}(r\partial_{r}p)=G_{0}c\qquad$	$\displaystyle\text{for}\quad r\leqslant R(t),$
(2.27b)	$\displaystyle p=0,\qquad$	$\displaystyle\text{for}\quad r\geqslant R(t),$
(2.27c)	$\displaystyle\sigma(R(t))=-\partial_{r}p(R(t)),\qquad$	$\displaystyle\text{on}\quad\partial\mathbb{B}_{R}(0).$
And by symmetry, we also require
(2.27d)	$\partial_{r}p(0)=0.$

By solving (2.26) we get

(2.28)

c(r,t)=c_{B}\frac{I_{0}(\sqrt{\lambda}r)}{I_{0}(\sqrt{\lambda}R)}\qquad\text{for}\quad r\leqslant R(t).

Plug the above expression into (2.27) to solve for $p$ , and we get:

(2.29)

p(r,t)=-\frac{G_{0}c_{B}}{\lambda I_{0}(\sqrt{\lambda}R(t))}I_{0}(\sqrt{\lambda}r)+\frac{G_{0}}{\lambda}c_{B}\qquad\text{for}\quad r\leqslant R(t).

And the boundary velocity is given by

(2.30)

\dot{R}=\sigma(R(t))=-\frac{\partial p}{\partial r}(R(t))=\frac{G_{0}c_{B}I_{1}(\sqrt{\lambda}R)}{\sqrt{\lambda}I_{0}(\sqrt{\lambda}R)}.

Note that as $R(t)\rightarrow\infty$ the speed limit is $\frac{G_{0}c_{B}}{\sqrt{\lambda}}$ , which recovers the speed for the traveling wave solution (2.21).

2.2.4. 2D radially symmetric model for the in vivo model

The computation is similar to the previous case, except that the equations for nutrient are replaced by


(2.31a)	$\displaystyle-\frac{1}{r}\partial_{r}(r\partial_{r}c)+\lambda c=0,\qquad$	$\displaystyle\text{for}\quad r\leqslant R(t),$
(2.31b)	$\displaystyle-\frac{1}{r}\partial_{r}(r\partial_{r}c)=c_{B}-c,\qquad$	$\displaystyle\text{for}\quad r\geqslant R(t).$

By solving above two equations and using the continuity of both $c$ and $\partial_{r}c$ at $R(t)$ , we get

(2.32)

c(r,t)=\left\{\begin{array}[]{rcr}c_{B}a_{0}(R)I_{0}(\sqrt{\lambda}r)\stackrel{{\scriptstyle\scriptscriptstyle\textup{def}}}{{=}}c^{\text{(i)}}(r,t),&\text{for}\quad r\leqslant R(t),\\ c_{B}(1+b_{0}(R)K_{0}(r))\stackrel{{\scriptstyle\scriptscriptstyle\textup{def}}}{{=}}c^{\text{(o)}}(r,t),&\text{for}\quad r\geqslant R(t),\\ \end{array}\right.

where $a_{0}$ and $b_{0}$ are given by


(2.33a)	$\displaystyle a_{0}(R)$	$\displaystyle=\frac{K_{1}(R)}{\sqrt{\lambda}K_{0}(R)I_{1}(\sqrt{\lambda}R)+K_{1}(R)I_{0}(\sqrt{\lambda}R)}\stackrel{{\scriptstyle\scriptscriptstyle\textup{def}}}{{=}}\frac{K_{1}(R)}{C(R)},$
(2.33b)	$\displaystyle b_{0}(R)$	$\displaystyle=-\frac{\sqrt{\lambda}I_{1}(R)}{\sqrt{\lambda}K_{0}(R)I_{1}(\sqrt{\lambda}R)+K_{1}(R)I_{0}(\sqrt{\lambda}R)}\stackrel{{\scriptstyle\scriptscriptstyle\textup{def}}}{{=}}-\frac{\sqrt{\lambda}I_{1}(R)}{C(R)}.$

Then from the pressure equations (2.27), we can solve and get

(2.34)

p(r,t)=-\frac{G_{0}c_{B}}{\lambda}a_{0}(R)I_{0}(\sqrt{\lambda}r)+\frac{G_{0}c_{B}}{\lambda}a_{0}(R)I_{0}(\sqrt{\lambda}R),\qquad\text{for}\quad r\leqslant R(t).

And the velocity of the boundary is given by

(2.35)

\dot{R}=\sigma(R(t))=\frac{G_{0}c_{B}K_{1}(R)I_{1}(\sqrt{\lambda}R)}{\lambda K_{0}(R)I_{1}(\sqrt{\lambda}R)+\sqrt{\lambda}K_{1}(R)I_{0}(\sqrt{\lambda}R)}\leqslant\frac{G_{0}c_{B}I_{1}(\sqrt{\lambda}R)}{\sqrt{\lambda}I_{0}(\sqrt{\lambda}R)},

which implies that the speed in the in vivo model is slower than that in the in vitro model. Again, by sending $R\rightarrow\infty$ , we get the limiting speed for the in vivo model is $\frac{c_{B}G_{0}}{\lambda+\sqrt{\lambda}}$ , which recovers the speed for the traveling wave in (2.25).

3. Framework of the perturbation analysis

We devote this section to establishing the general framework of our asymptotic analysis. Such analysis involves classical techniques which was originally developed by Mullins et al. in [mullins1963morphological] and widely used in [cristini2003nonlinear, macklin2007nonlinear, pham2018nonlinear, turian2019morphological, lu2020complex, lu2022nonlinear], whereas we present it as generic methodology which in theory can be applied to other problems as well.

We divide our analysis into three parts as follows.

3.1. Perturbation of the boundary

We study the perturbation of two kinds of boundaries, the radial boundary and the front of traveling waves, and the relationship between them. In either case, we have a proper coordinate system denoted as $(\zeta,\vartheta)$ . For simplicity, we assume the boundary profile is a curve $\mathcal{O}_{t}\subseteq\mathbb{R}^{2}$ , which can be parameterized by the variable $\vartheta$ in the following form:

(3.1)

\mathcal{O}_{t}(\vartheta)=\left\{(\zeta,\vartheta)|\zeta=\mathcal{Z}(t,\vartheta),\vartheta\in\mathcal{R}\right\}

with some contour index function $\mathcal{Z}(t,\vartheta)$ and range $\mathcal{R}$ .

For the radial case, the unperturbed tumor region at time $t$ is given by a disk with radius $R(t)$ , that is $D(t)=\mathbb{B}_{R(t)}$ . In this case, equations and functions are naturally presented in terms of the polar coordinate. Therefore, $(\zeta,\vartheta)=(r,\theta)$ and $\mathcal{R}=\left[-\pi,\pi\right)$ . Further more, the tumor boundary at time $t$ can be written as:

(3.2)

\mathcal{B}_{t}(\theta)=\left\{(r,\theta)|r=R(t),\theta\in\left[-\pi,\pi\right)\right\}.

For the traveling wave case, we employ the Euler coordinate $(\xi,y)$ (where $\xi=x-\sigma t$ ). In this case, the tumor region is a half plane with a moving front. We fix the front (propagate to the right) at $\xi=0$ with traveling speed $\sigma$ , and the tumor region, therefore, become $D(t)=\left\{(\xi,y)|\xi\leqslant 0\right\}$ . Then, we write the traveling front more clearly in the parameter curve form:

(3.3)

\mathcal{B}(y)=\left\{(\xi,y)|\xi=0,y\in\mathbb{R}\right\}.

For the purpose of introducing perturbation method in a general framework, we combine the two scenarios above in the following unified notations. Let $D(t)$ still presents the tumor region at time $t$ ; and the boundary curve writes

(3.4)

\mathcal{B}_{t}(\vartheta)=\left\{(\zeta,\vartheta)|\zeta=Z(t),\vartheta\in\mathcal{R}\right\}.

Moreover, any point $B\in\mathcal{B}_{t}$ can be presented as $B(Z,\vartheta_{*})$ for some $\vartheta_{*}\in\mathcal{R}$ . Note that in either case above, the index function $Z(t)$ is independent on the parameter variable $\vartheta$ . More precisely, for the radial case $Z(t)=R(t)$ , and (3.4) stands for (3.2); for the traveling wave case, (3.4) stands for (3.3) with $Z(t)$ takes constant value $0$ .

Next, we add a small perturbation to the two kinds of boundaries. From the parameterization representation point of view, the perturbation replaces the boundary curve (3.4) by:

(3.5)

\tilde{\mathcal{B}}_{t}(\vartheta)=\left\{(\zeta,\vartheta)|\zeta=Z(t)+\delta(t)\mathcal{P}(\vartheta),\vartheta\in\mathcal{R}\right\},

where $\delta(t)\ll 1$ stands for the amplitude of the perturbation, and $\mathcal{P}(\vartheta)$ characterizes the perturbation profile. Thus, the perturbed boundary at time $t$ is still parameterized by the variable $\vartheta$ . Intuitively, (3.5) means that the perturbation will push the point $(Z,\vartheta_{*})\in\mathcal{B}_{t}$ to $(Z+\delta\mathcal{P}(\vartheta_{*}),\vartheta_{*})\in\tilde{\mathcal{B}}_{t}$ for any $\vartheta_{*}\in\mathcal{R}$ . Note that the perturbation form $\eqref{eqn:boundary perturbation}$ enables the evolution of the perturbation term to reduce to the evolution of the amplitude function $\delta(t)$ while its spatial profile persists. Such an ansatz with temporal and spatial degrees of freedom separated makes sense only when the profile function represents a typical model of a general classical of contours. In the next, we explain how to choose the perturbation profiles in the two cases.

In the radial symmetry case, the profile $\mathcal{P}(\theta)$ is parameterized by $\theta\in[-\pi,\pi)$ and it can be expressed as a Fourier expansions in general. In particular, for the single wave perturbation with wave number $l$ , $\mathcal{P}(\theta)$ takes the form of:

(3.6)

\mathcal{P}(\theta)=C_{1}\cos{l\theta}+C_{2}\sin{l\theta},\quad\text{with}\quad l\in\mathbb{N}^{+},

where $C_{1},C_{2}$ are constant coefficients. Note that by rotating the coordinate system and rescaling on $\delta(t)$ , without loss of generality we can simply take

(3.7)

\mathcal{P}(\theta)=\cos{l\theta}\stackrel{{\scriptstyle\scriptscriptstyle\textup{def}}}{{=}}\mathcal{P}_{l}(\theta).

For the traveling wave case, the profile is parameterized by $y\in\mathbb{R}$ . By a similar reason to the radial case, we can simply consider

(3.8)

\mathcal{P}_{l}(y)=\cos{ly},\quad\text{with}\quad l\in\mathbb{R}^{+},

otherwise we can just shift the profile along $y$ -axis.

It is important to note that for the perturbation of the traveling wave, we are actually allowed to take $\mathcal{P}(y)=\cos{ly}$ with $l\in\mathbb{R}^{+}$ . However, only integer frequencies perturbation are reasonable for the radial case, since $\mathcal{P}(\theta)$ has to be a $2\pi$ -periodic function.

3.2. Solutions after perturbation

Let $\tilde{D}(t)$ , enclosed by $\tilde{\mathcal{B}}_{t}$ , denote the tumoral region after the perturbation. Then the perturbed functions $(\tilde{c},\tilde{p},\tilde{\rho})$ satisfy the equations (boundary conditions will be specified in the next subsection):


(3.9a)	$\displaystyle-\Delta\tilde{c}+\Psi(\tilde{\rho},\tilde{c})=0,\qquad$	$\displaystyle\text{on}\quad\mathbb{R}^{2},$
(3.9b)	$\displaystyle-\Delta\tilde{p}=G_{0}\tilde{c},\qquad$	$\displaystyle\text{in}\quad\tilde{D}(t),$
recall that $\Psi(\tilde{\rho},\tilde{c})$ reads (2.15) in the in vitro model and (2.16) in the in vivo model.

When the boundary perturbation vanishes, (3.9) reduce to the the unperturbed problem, where the solutions are given in a closed-form. In the presence of the boundary perturbation, we still have $\tilde{\rho}=\chi_{\tilde{D}}$ since it remains as a patch, but the solution to $\tilde{c}$ and $\tilde{p}$ are no longer available. However, we can alternatively seek asymptotic solutions of $\tilde{c}$ and $\tilde{p}$ with respect to $\delta$ , while the condition $\tilde{\rho}=\chi_{\tilde{D}}$ help to linearize the calculation. We elaborate the asymptotic analysis procedures as follows.

Firstly, corresponding to the small perturbation (3.5), we have the following asymtotic expansion with respect to the small value $\delta$ :


(3.10a)	$\displaystyle\tilde{c}(\zeta,\vartheta,t)$	$\displaystyle=c_{0}(\zeta,t)+\delta c_{1}(\zeta,\vartheta,t)+O(\delta^{2}),$
(3.10b)	$\displaystyle\tilde{p}(\zeta,\vartheta,t)$	$\displaystyle=p_{0}(\zeta,t)+\delta p_{1}(\zeta,\vartheta,t)+O(\delta^{2}).$

Since the perturbation scale is assumed to be very small, i.e., $\delta\ll 1$ , the behavior of the perturbed solutions are dominated by the unperturbed ones. Thus, the leading order terms $c_{0}(\zeta,t)$ and $p_{0}(\zeta,t)$ take the same expression as the solutions without perturbation, which have been solved in Section 2.2. On the other hand, the main response corresponding to the perturbation are captured by the first-order terms $c_{1}(\zeta,\vartheta,t)$ and $p_{1}(\zeta,\vartheta,t)$ . Note that besides variable $\zeta$ they depend on $\vartheta$ as well.

We continue to investigate the structures of $c_{1}$ and $p_{1}$ when the perturbation profile (3.5) is given by $\mathcal{P}(\vartheta)=\mathcal{P}_{l}(\vartheta)$ , here $\mathcal{P}_{l}(\vartheta)$ presents (3.7) or (3.8) in the respective case. In either case, the perturbed tumoral region $\tilde{D}$ still possess a symmetry, or periodicity, respect to the parameter $\vartheta$ ( $\theta$ for the radial case and $y$ for the traveling wave case). Then we have following conclusion for the perturbed solutions $(\tilde{c},\tilde{p})$ .

Lemma 3.1.

If the perturbed solutions are unique, then they must process the same periodicity as the boundary geometry.

Proof.

For either scenario, the front of traveling wave or radially symmetric boundary, we assume the boundary has periodicity $\vartheta^{*}$ . Then, with respect to (3.5) we have:

(3.11)

\zeta(\vartheta)=\zeta(\vartheta+\vartheta^{*}),

where $\zeta(\vartheta)\stackrel{{\scriptstyle\scriptscriptstyle\textup{def}}}{{=}}Z(t)+\delta\mathcal{P}(\vartheta)$ . For $\mathcal{P}(\vartheta)=\cos{l\vartheta}$ , $\vartheta^{*}$ is given by $\vartheta^{*}=\frac{2\pi}{l}$ . We define the translation operator $\tau_{\vartheta^{*}}(\zeta,\vartheta):(\zeta,\vartheta)\mapsto(\zeta,\vartheta+\vartheta^{*})$ . One can easily observe that the nutrient equations (for either in vitro or in vivo) and the pressure equation are both invariant under $\tau_{\vartheta^{*}}$ since the operator simply corresponds to a translation or a rotation, and diffusion operator is isotropic. Moreover, the boundary geometry and boundary conditions remain the same under the operator $\tau_{\vartheta^{*}}$ as well. Thus, the uniqueness of the solution yield that the unique solutions $c^{*}$ and $p^{*}$ must possess the same periodicity as the boundary geometry. That is,


(3.12a)	$\displaystyle c^{*}(\zeta,\vartheta)$	$\displaystyle=c^{}(\tau_{\vartheta^{}}(\zeta,\vartheta)),$
(3.12b)	$\displaystyle p^{*}(\zeta,\vartheta)$	$\displaystyle=p^{}(\tau_{\vartheta^{}}(\zeta,\vartheta)).$

∎

According to the above lemma, to be consistent with the boundary’s periodicity, we expand $c_{1}(\zeta,\vartheta,t)$ and $p_{1}(\zeta,\vartheta,t)$ as Fourier series, and (3.10) can be further written as:


(3.13a)	$\displaystyle\tilde{c}(\zeta,\vartheta,t)$	$\displaystyle=c_{0}(\zeta,t)+\delta(t)\Sigma_{k=1}^{\infty}c_{1}^{k}(\zeta,t)\mathcal{P}_{l}^{k}(\vartheta)+O(\delta^{2}),$
(3.13b)	$\displaystyle\tilde{p}(\zeta,\vartheta,t)$	$\displaystyle=p_{0}(\zeta,t)+\delta(t)\Sigma_{k=1}^{\infty}p_{1}^{k}(\zeta,t)\mathcal{P}_{l}^{k}(\vartheta)+O(\delta^{2}),$
where $\mathcal{P}_{l}^{k}(\vartheta)=\cos{kl\vartheta}$ .

In the above expansions, $c_{1}^{k}(\zeta,t)$ and $p_{1}^{k}(\zeta,t)$ (with $k\in\mathbb{N}^{+}$ ) serve as the Fourier coefficients with $O(1)$ . From the calculation in the later sections (Section 4.2 and Section LABEL:sec:Radial_The_detailed_calculations_for_the_two_nutrient_regimes), we will see that only $c_{1}^{1}$ and $p_{1}^{1}$ , the coefficients of the wave number that is the same as the perturbation, do not vanish. Therefore, it suffices to keep the first term in the series (3.13) and drop the superscript in $c_{1}^{1}$ , $p_{1}^{1}$ and, $\mathcal{P}_{l}^{1}$ . Thus, (3.13) writes


(3.14a)	$\displaystyle\tilde{c}(\zeta,\vartheta,t)$	$\displaystyle=c_{0}(\zeta,t)+\delta(t)c_{1}(\zeta,t)\mathcal{P}_{l}(\vartheta)+O(\delta^{2}),$
(3.14b)	$\displaystyle\tilde{p}(\zeta,\vartheta,t)$	$\displaystyle=p_{0}(\zeta,t)+\delta(t)p_{1}(\zeta,t)\mathcal{P}_{l}(\vartheta)+O(\delta^{2}).$

In the traveling wave case, the dependency of $t$ can be removed for the terms $c_{j}$ and $p_{j}$ ( $j=\left\{0,1\right\}$ ), since the unperturbed tumor boundary do not evolve in time. Finally, by plugging the expansion (3.14) into (3.9) and collect the first order terms we get


(3.15a)	$\displaystyle-\Delta(c_{1}(\zeta,t)\mathcal{P}_{l}(\vartheta))+\lambda c_{1}(\zeta,t)\mathcal{P}_{l}(\vartheta)=0,\qquad$	$\displaystyle\text{in}\quad\tilde{D}(t),$
(3.15b)	$\displaystyle-\Delta(p_{1}(\zeta,t)\mathcal{P}_{l}(\vartheta))=G_{0}(c_{1}(\zeta,t)\mathcal{P}_{l}(\vartheta)),\qquad$	$\displaystyle\text{in}\quad\tilde{D}(t),$
for either nutrient regime. In addition, for the in vivo model $c_{1}$ also satisfies
(3.15c)	$-\Delta(c_{1}(\zeta,t)\mathcal{P}_{l}(\vartheta))+c_{1}(\zeta,t)\mathcal{P}_{l}(\vartheta)=0,\qquad\text{in}\quad\mathbb{R}^{2}\setminus\tilde{D}(t).$

where we used the fact that the zero order terms satisfy (3.9). By solving (3.15), one can get the solutions of $c_{1}$ and $p_{1}$ for the respective models. Note that (3.15) implies the expression of $c_{1}$ and $p_{1}$ depend on the wave number $l$ . The detailed computation will be carried out for the specific cases in the later sections.

3.3. Match the boundary condition

In the last part of this section, we explain how to determine the particular solutions of $c_{1}$ and $p_{1}$ by matching the boundary conditions. We also show that by using the expression of $p_{1}$ , one can determine the evolution of the perturbation magnitude.

In this section, we always assume the perturbation profile $\mathcal{P}(\vartheta)$ is given by $\mathcal{P}_{l}(\vartheta)$ . And note that given $t$ for any $\vartheta\in\mathcal{R}$ , $(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta)$ presents a point on the perturbed boundary $\tilde{\mathcal{B}}_{t}$ , which is originally at the position $(Z,\vartheta)\in\mathcal{B}_{t}$ . Recall that $c_{0}$ and $c_{1}$ (similarly for $p_{0}$ and $p_{1}$ ) only depend on the variable $\zeta$ in space, and the unperturbed boundary $\mathcal{B}_{t}$ is characterized as the contour of $\zeta$ with level set index $Z(t)$ (see (3.4)).

Since the analytical solutions are not available for the perturbed problem, it is not practical to enforce the boundary conditions in the precise way. Instead, since we seek the first order perturbation solutions due to the boundary variation, we can approximately match the the perturbed solutions at the perturbed boundary up to $O(\delta^{2})$ error with the their evaluations at the unperturbed boundary.

For the in vitro model, the perturbed solution $\tilde{c}$ satisfies the boundary condition:

(3.16)

\tilde{c}=c_{B},\quad\text{at}\quad\tilde{\mathcal{B}}_{t}.

Thus by using expansion (3.14a), we can evaluate $\tilde{c}$ at the perturbed boundary point $(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta)$ to get

(3.17)		$\displaystyle\tilde{c}(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta,t)$	$\displaystyle=c_{0}(Z+\delta\mathcal{P}_{l}(\vartheta),t)+\delta c_{1}(Z+\delta\mathcal{P}_{l}(\vartheta),t)\mathcal{P}_{l}(\vartheta)+O(\delta^{2})$
		$\displaystyle=c_{0}(Z,t)+\delta\partial_{\zeta}c_{0}(Z,t)\mathcal{P}_{l}(\vartheta)+\delta c_{1}(Z,t)\mathcal{P}_{l}(\vartheta)+O(\delta^{2}),$

where we used the Taylor expansions for $c_{0}(Z+\delta\mathcal{P}_{l},t)$ and $c_{1}(Z+\delta\mathcal{P}_{l},t)$ . By the boundary conditions of the perturbed and unperturbed problems, we have $\tilde{c}(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta,t)=c_{0}(Z,t)=c_{B}$ for arbitrary $\vartheta\in\mathcal{R}$ . Thus the zero order terms in (3.17) are canceled out, and by balancing the first order terms in (3.17) we get

(3.18)

\partial_{\zeta}c_{0}(Z,t)+c_{1}(Z,t)=0.

While for the in vivo model, $\tilde{c}$ and its normal derivative are both continuous at $\tilde{\mathcal{B}}_{t}$ . And in either kind of boundary, the normal derivative of $\tilde{c}(\zeta,\vartheta,t)$ on $\tilde{\mathcal{B}}_{t}$ is given by $\partial_{\zeta}\tilde{c}(\zeta,\vartheta,t)$ for any $\tilde{B}(\zeta,\vartheta)\in\tilde{\mathcal{B}}_{t}$ . And if we decompose the solution $\tilde{c}$ according to the regions as $\tilde{c}=\tilde{c}^{\text{(i)}}\chi_{\tilde{D}}+\tilde{c}^{\text{(o)}}\chi_{\mathbb{R}^{2}\setminus\tilde{D}}$ , i.e., let $\tilde{c}^{\text{(i)}}$ denotes the nutrient solution inside the tumor, and $\tilde{c}^{\text{(o)}}$ presents the outside solution. Then, the continuity at the boundary yields


(3.19a)	$\displaystyle\tilde{c}^{\text{(i)}}(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta,t)=\tilde{c}^{\text{(o)}}(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta,t),\qquad$	$\displaystyle\forall\vartheta\in\mathcal{R}$
(3.19b)	$\displaystyle\partial_{\zeta}\tilde{c}^{\text{(i)}}(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta,t)=\partial_{\zeta}\tilde{c}^{\text{(o)}}(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta,t),\qquad$	$\displaystyle\forall\vartheta\in\mathcal{R},$

With the same spirit of (3.17), for $\partial_{\zeta}\tilde{c}(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta,t)$ we have

(3.20)

\partial_{\zeta}\tilde{c}(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta,t)=\partial_{\zeta}c_{0}(Z,t)+\partial^{2}_{\zeta}c_{0}(Z,t)\delta\mathcal{P}_{l}(\vartheta)+\partial_{\zeta}c_{1}(Z,t)\delta\mathcal{P}_{l}(\vartheta)+O(\delta^{2}).

Since $c_{0}$ is the solution to the unperturbed problem (given by (2.32) or (2.23) for the respect case), $c_{0}$ and $c_{0}^{\prime}$ are both continuous at the unperturbed boundary $\mathcal{B}_{t}=\left\{\zeta=Z(t)\right\}$ . More precisely,

(3.21)		$\displaystyle c_{0}^{\text{(i)}}(Z,t)$	$\displaystyle=c_{0}^{\text{(o)}}(Z,t),$
(3.22)		$\displaystyle\partial_{\zeta}c_{0}^{\text{(i)}}(Z,t)$	$\displaystyle=\partial_{\zeta}c_{0}^{\text{(o)}}(Z,t).$

Thus, by using the expansions (3.17) and (3.20), the boundary condition (3.19) yields


(3.23a)	$\displaystyle c_{1}^{\text{(i)}}(Z,t)$	$\displaystyle=c_{1}^{\text{(o)}}(Z,t),$
(3.23b)	$\displaystyle\partial^{2}_{\zeta}c_{0}^{\text{(i)}}(Z,t)+\partial_{\zeta}c_{1}^{\text{(i)}}(Z,t)$	$\displaystyle=\partial^{2}_{\zeta}c_{0}^{\text{(o)}}(Z,t)+\partial_{\zeta}c_{1}^{\text{(o)}}(Z,t).$

By using (3.18) or (3.23) as the boundary condition for $c_{1}$ , we can work out the particular solution of $c_{1}$ in the respective cases. We mention that when the boundary is the traveling front, to carry out the full expression of $c_{1}$ , we also need to use the boundary condition $c_{1}(-\infty,y)=c_{1}(+\infty,y)=0$ for any $y\in\mathbb{R}$ , which is derived from $\tilde{c}(-\infty,y)=0$ and $\tilde{c}(+\infty,y)=c_{B}$ for any $y\in\mathbb{R}$ . The detail calculations will be carried out for each specific case later (Section 4.2 and Section LABEL:sec:Radial_The_detailed_calculations_for_the_two_nutrient_regimes).

In either nutrient model, the perturbed pressure solution $\tilde{p}$ satisfies the boundary condition

(3.24)

\tilde{p}=0,\quad\text{at}\quad\tilde{\mathcal{B}}_{t}.

Similar to the previous calculations. By using the expansion (3.14b) to evaluate $\tilde{p}$ at $(Z+\delta\mathcal{P}_{l}(\vartheta))\in\tilde{\mathcal{B}}_{t}$ , we get

(3.25)

\tilde{p}(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta,t)=p_{0}(Z,t)+\partial_{\zeta}p_{0}(Z,t)\delta\mathcal{P}_{l}(\vartheta)+p_{1}(Z,t)\delta\mathcal{P}_{l}(\vartheta)+O(\delta^{2}).

The perturbed and unperturbed boundary condition yield that $\tilde{p}(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta,t)$ and $p_{0}(Z,t)$ both equal to zero. In particular, $p_{0}$ as the unperturbed solution has already been solved in the Section 2.2. Thus we get

(3.26)

\partial_{\zeta}p_{0}(Z,t)+p_{1}(Z,t)=0.

Then by using the expression of $c_{1}$ (see (3.15b)) and (3.26), we can further determine the particular solution of $p_{1}$ . For the traveling wave case, we also use the condition $\partial_{\zeta}p_{1}(-\infty,y)=0$ for $\forall y\in\mathbb{R}$ . Finally, the normal boundary speed (2.10) yields:

(3.27)

\frac{d(Z(t)+\delta(t)\mathcal{P}_{l}(\vartheta))}{dt}=-\partial_{\zeta}\tilde{p}(Z+\delta\mathcal{P}_{l}(\vartheta),\vartheta,t).

By plugging the expression of $\tilde{p}$ into (3.27) and taking Taylor expansion for the $\zeta$ variable, we get

(3.28)

\frac{dZ}{dt}+\frac{d\delta}{dt}\mathcal{P}_{l}(\vartheta)=-\left(\partial_{\zeta}p_{0}(Z,t)+\partial^{2}_{\zeta}p_{0}(Z,t)\delta\mathcal{P}_{l}(\vartheta)+\partial_{\zeta}p_{1}(Z,t)\delta\mathcal{P}_{l}(\vartheta)+O(\delta^{2})\right).

Since the unperturbed problem yields $\frac{dZ}{dt}=-\partial_{\zeta}p_{0}(Z,t)$ , the above identity can be further simplified into

(3.29)

\delta^{-1}\frac{d\delta}{dt}=-\left(\partial^{2}_{\zeta}p_{0}(Z,t)+\partial_{\zeta}p_{1}(Z,t)+O(\delta)\right).

In the end, we determine the evolution of the perturbation magnitude by the sign of $\delta^{-1}\frac{d\delta}{dt}$ . If it is positive, then it implies the growing of the magnitude. For the radial case, $Z(t)$ is given by the unperturbed tumor radius $R(t)$ , therefore $\delta^{-1}\frac{d\delta}{dt}\sim-(\partial^{2}_{r}p_{0}(R(t),t)+\partial_{r}p_{1}(R(t),t))$ . Note that the leading order of $\delta^{-1}\frac{d\delta}{dt}$ is independent on $\theta$ , which parameterize the boundaries. While, for the traveling wave case, $Z(t)=0$ , thus $\delta^{-1}\frac{d\delta}{dt}\sim-(\partial^{2}_{\xi}p_{0}(0)+\partial_{\xi}p_{1}(0))$ , which is independent of $y$ . Furthermore, under the same nutrient regime, we expect that the boundary instability of the radius case will coincide with that of the traveling wave when $R$ increase to infinity.

4. Stability of traveling waves in the two nutrient models

In this section, we study the boundary stability of the traveling wave front under two nutrient regimes. In Section 4.1, we establish the set up and main conclusions. In Section 4.2, we work out the expression of $\delta^{-1}\frac{d\delta}{dt}$ for the two nutrient models, which serves as the proof of Theorem 4.1. And in Section 4.3, we prove the mathematical properties of $\delta^{-1}\frac{d\delta}{dt}$ summarized in Corollary 4.2, and these properties further yield the boundary behaviors summarized in Remark 4.3.

4.1. Setup and main results

As presented in Section 3.1, in the traveling wave case we employ the Euler coordinate system $(\xi,y)$ . In the absence of perturbation, the tumor boundary is defined by (3.3) with the level set index $Z(t)=0$ , and the tumor region is the left half space $D(t)=\left\{(\xi,y)|\xi\leqslant 0\right\}$ . Then following the framework of Section 3, we consider the perturbation by a single wave mode:

(4.1)

\mathcal{P}_{l}(y)=\cos{ly}\qquad\text{with}\quad l\in\mathbb{R}^{+},

thus the perturbed boundary (3.5) writes:

(4.2)

\tilde{\mathcal{B}}_{t}(y)=\left\{(\xi,y)|\xi=\delta(t)\cos{ly},y\in\mathbb{R}\right\},

and the perturbed tumor region becomes

(4.3)

\tilde{D}(t)=\left\{(\xi,y)|\xi\leqslant\delta(t)\cos{ly},y\in\mathbb{R}\right\}.

Then correspond to the above perturbation, the perturbed solutions $c$ and $p$ solves (3.9). Note that we dropped the tilde of the perturbed solutions for simplicity. Further more, the perturbed solutions possess the asymptotic expansions:


(4.4a)	$\displaystyle c(\xi,y,t)$	$\displaystyle=c_{0}(\xi)+\delta(t)c_{1}(\xi,y)+O(\delta^{2}),$
(4.4b)	$\displaystyle p(\xi,y,t)$	$\displaystyle=p_{0}(\xi)+\delta(t)p_{1}(\xi,y)+O(\delta^{2}),$
where the leading order terms $c_{0}$ and $p_{0}$ correspond to the solution of the unperturbed problems, which have been solved in Section 2.2.1 and Section 2.2.2 for the respective nutrient model. And the first order terms $c_{1}(\xi,y,t)$ and $p_{1}(\xi,y,t)$ can be further expanded as

(4.4c)	$\displaystyle c_{1}(\xi,y)=\Sigma_{k=1}^{\infty}c_{1}^{k}(\xi)\cos{kly},$
(4.4d)	$\displaystyle p_{1}(\xi,y)=\Sigma_{k=1}^{\infty}p_{1}^{k}(\xi)\cos{kly}.$
with $l\in\mathbb{R}^{+}$ , so that $c_{1}$ has the same periodicity as the boundary geometry.

However, from the calculation later we will see that $c_{1}^{k}(\xi)=p_{1}^{k}(\xi)=0$ for any $k\neq 1$ .

For the in vivo model, we use the superscript (i) or (o) to denote the solution inside or outside the tumor region $\tilde{D}(t)$ . Then according to (3.15a) and (3.15c) we have


(4.5a)		$\displaystyle-\Delta c_{1}^{\text{(i)}}(\xi,y)+\lambda c_{1}^{\text{(i)}}(\xi,y)=0,$
(4.5b)		$\displaystyle-\Delta c_{1}^{\text{(o)}}(\xi,y)+c_{1}^{\text{(o)}}(\xi,y)=0.$

And by using the expansion in (3.17) and (3.20), the series form of $c_{1}(\xi,y)$ in $\eqref{eqn:TW expansion of c1}$ , the boundary condition (3.19) yields


(4.6a)	$\displaystyle c_{1}^{\text{(i)},k}(0)$	$\displaystyle=c_{1}^{\text{(o)},k}(0),\quad\forall k\in\mathbb{N}^{+},$
(4.6b)	$\displaystyle\partial_{\xi}c_{1}^{\text{(i)},k}(0)$	$\displaystyle=\partial_{\xi}c_{1}^{\text{(o)},k}(0),\quad\forall k\geqslant 2,$
(4.6c)	$\displaystyle\partial^{2}_{\xi}c_{0}^{\text{(i)}}(0)+\partial_{\xi}c_{1}^{\text{(i)},1}(0)$	$\displaystyle=\partial^{2}_{\xi}c_{0}^{\text{(o)}}(0)+\partial_{\xi}c_{1}^{\text{(o)},1}(0).$
Further more, the assumptions $c(-\infty,y)=0$ and $c(+\infty,y)=c_{B}$ for any $y\in\mathbb{R}$ yields

(4.6d)	$\displaystyle c_{1}^{\text{(i)},k}(-\infty)=0,$
(4.6e)	$\displaystyle c_{1}^{\text{(o)},k}(+\infty)=0,$
for any $k\in\mathbb{N}^{+}$ .

For the in vitro model, $c$ presents the nutrient solution inside the tumor and equation (3.15a) writes

(4.7)

-\Delta c_{1}(\xi,y)+\lambda c_{1}(\xi,y)=0,\quad\text{in}\quad\tilde{D}(t).

By using (3.17) and the series expansion of $c_{1}$ in (4.4c), the boundary condition (3.16) yields:


(4.8a)		$\displaystyle\partial_{\zeta}c_{0}(0)+c_{1}^{1}(0)=0.$
(4.8b)		$\displaystyle c_{1}^{k}(0)=0,\quad\forall k\geqslant 2.$
and similar to the in vivo model, the assumption $c(-\infty,y)=0$ for any $y\in\mathbb{R}$ gives
(4.8c)		$c_{1}^{k}(-\infty)=0,\qquad\forall k\in\mathbb{N}^{+}.$

Once $c_{1}(\xi,y)$ is determined, we can move on to the study of the first order term of pressure, i.e., $p_{1}(\xi,y)$ . Under either nutrient regime, $p_{1}(\xi,y)$ satisfies the equation:

(4.9)

-\Delta p_{1}(\xi,y)=G_{0}c_{1}(\xi,y),\quad\text{in}\quad\tilde{D}(t).

By using the expansion (3.25), the series form of $p_{1}$ in (4.4d), the boundary condition (3.24) yields


(4.10a)		$\displaystyle\partial_{\zeta}p_{0}(0)+p_{1}^{1}(0)=0,$
(4.10b)		$\displaystyle p_{1}^{k}(0)=0,\quad\forall k\geqslant 2.$
On the other hand, for the traveling wave case we require $\partial_{\zeta}p(-\infty,y)=0$ , which further yields $\partial_{\xi}p_{1}(-\infty,y)=0$ . Therefore,
(4.10c)		$\partial_{\xi}p_{1}^{k}(-\infty)=0,\qquad\text{for}\quad\forall k\in\mathbb{N}^{+}.$

Once the expression of $p_{1}(\xi,y)$ is determined, we can further work out the expression of $\delta^{-1}\frac{d\delta}{dt}$ for the two nutrient regimes, which determines the evolution of the perturbation amplitude. Now we establish the main conclusions, and the details of the calculation will be left to the next subsection.

Theorem 4.1.

Given growing rate $G_{0}>0$ , background concentration $c_{B}>0$ , nutrient consumption rate $\lambda>0$ , and perturbation frequency $l\in\mathbb{R}^{+}$ . The perturbation evolution function, $\delta^{-1}\frac{d\delta}{dt}$ , of the in vitro model is given by:

(4.11)

\delta^{-1}\frac{d\delta}{dt}=\frac{G_{0}c_{B}}{\sqrt{\lambda}}\cdot(\sqrt{\lambda}-\sqrt{\lambda+l^{2}})+O(\delta)\stackrel{{\scriptstyle\scriptscriptstyle\textup{def}}}{{=}}F_{1}(\lambda,l)+O(\delta).

For the in vivo model, $\delta^{-1}\frac{d\delta}{dt}$ is given by:

(4.12)		$\displaystyle\delta^{-1}\frac{d\delta}{dt}$	$\displaystyle=\frac{G_{0}c_{B}}{\sqrt{\lambda}}\left(\frac{\sqrt{\lambda}-l}{\sqrt{\lambda}+1}+\frac{l-\sqrt{\lambda+l^{2}}}{\sqrt{\lambda+l^{2}}+\sqrt{1+l^{2}}}\right)+O(\delta)$
		$\displaystyle\stackrel{{\scriptstyle\scriptscriptstyle\textup{def}}}{{=}}F_{2}(\lambda,l)+O(\delta).$

Refer to caption — Figure 1. $F_{2}(\lambda,l)$ for $G_{0}=1$ , $c_{B}=100$ , left one: $\lambda=0.8$ , right one: $\lambda=2$ .

Note that in either nutrient regime, the value of $G_{0},c_{B}>0$ serve as scaling parameters, therefore do not influence the quantitative behavior of $\delta^{-1}\frac{d\delta}{dt}$ . For the in vitro model, one can easily check that $F_{1}(\lambda,l)$ is always negative. For the in vivo model, $F_{2}(\lambda,l)$ is plotted in Figure 1 for different choice of $\lambda$ . Base on the expression of $F_{1}$ , $F_{2}$ , and the observations from Figure 1, we prove following mathematical properties for the evolution equations.

Corollary 4.2.

Fix $G_{0}>0$ and $c_{B}>0$ . For any $\lambda>0$ and $l>0$ , $F_{1}(\lambda,l)$ in (4.11) is always negative, therefore the perturbation amplitude always decreases in the in vitro model. For the in vivo model, $\delta^{-1}\frac{d\delta}{dt}$ is given by $F_{2}(\lambda,l)$ as in (4.12). When $l$ approaches zero, $F_{2}$ has the asymptote:


(4.13a)		$\displaystyle F_{2}(\lambda,l)\sim\frac{(\lambda-1)\cdot l^{2}+O(l^{3})}{2\lambda(\sqrt{\lambda+1})(\sqrt{\lambda+l^{2}}+\sqrt{1+l^{2}})},$
and the limit at infinity:
(4.13b)		$\displaystyle\lim_{l\rightarrow+\infty}F_{2}(\lambda,l)\rightarrow-\infty.$

Further more, for $\lambda>1$ there exists $L>0$ such that $F_{2}(\lambda,l)>0$ for $l\in(0,L)$ . And for $0<\lambda\leqslant 1$ , we have $F_{2}(\lambda,l)<0$ for any $l>0$ . And in particular, when $\lambda=1$ , $F_{2}$ can be further simplified into:

(4.14)

F_{2}(1,l)=\frac{l(1-\sqrt{1+l^{2}})}{2\sqrt{1+l^{2}}}<0,\qquad\text{for}\quad l>0.

Remark 4.3.

The mathematical properties of $F_{2}(\lambda,l)$ in Corollary 4.2 imply the following boundary behaviors:

(1)

When the consumption rate is relatively large, the amplitude of low-frequency perturbations can grow in time, and the boundary propagation, therefore, becomes unstable. However, the amplitude of high-frequency perturbation decays. Correspondingly, the perturbed, wave like, boundary degenerates to the vertical line.
(2)

When the nutrient consumption rate is relatively small, the perturbation amplitude decreases for perturbation of any frequency, i.e., the wave like boundary always evolve to a vertical line in this regime.

Remark 4.4.

For either nutrient model, $\delta^{-1}\frac{d\delta}{dt}\rightarrow 0$ as $l\rightarrow 0$ , we claim that this relate to the single wave perturbation of the radially symmetric solution as its radius $R$ goes to infinity. This relationship will be further discussed in Section LABEL:sec:Relationship_between_the_Radial_boundary_and_the_traveling_wave_boundary.

4.2. The detailed calculations for the two nutrient regimes

In this subsection we work out the expression of $\delta^{-1}\frac{d\delta}{dt}$ in Theorem 4.1 for the two nutrient models.

For the in vivo case. Plugging the expansion (4.4c) of $c_{1}(\xi,y)$ into (4.5), together with the conditions (4.6d) and (4.6e), for any $k\in\mathbb{N}^{+}$ we have:


(4.15a)	$\displaystyle c^{\text{(i)},k}_{1}(\xi)$	$\displaystyle=a_{1}^{k}e^{\sqrt{\lambda+k^{2}}\xi}\qquad$	$\displaystyle\text{for}\quad\xi\leqslant 0,$
(4.15b)	$\displaystyle c^{\text{(o)},k}_{1}(\xi)$	$\displaystyle=b_{1}^{k}e^{-\sqrt{1+k^{2}}\xi}\qquad$	$\displaystyle\text{for}\quad\xi\geqslant 0.$

Recall that the leading order terms $c_{0}^{\text{(i)}}(\xi)$ and $c_{0}^{\text{(o)}}(\xi)$ are given by (2.23). Then (4.6a)-(4.6c) yield $a_{1}^{k}=b_{1}^{k}=0$ for any $k\neq 1$ , for $k=1$ we get nontrivial solution:

(4.16)

c^{\text{(i)},1}_{1}(\xi)=c^{\text{(o)},1}_{1}(\xi)=-\frac{\sqrt{\lambda}c_{B}}{\sqrt{\lambda+l^{2}}+\sqrt{1+l^{2}}}e^{\sqrt{\lambda+l^{2}}\xi}.

By now, $c_{1}^{\text{(i)},k}(\xi)$ and $c_{1}^{\text{(o)},k}(\xi)$ are determined for any $k$ . Therefore, $c_{1}(\xi,y)$ is determined. Then by solving equation (4.9) together with boundary conditions (4.10) (with $p_{0}$ given by (2.24)), we get $p_{1}^{k}(\xi)=0$ for any $k\neq 1$ , and:

(4.17a)		$p_{1}^{1}(\xi)=Ae^{l\xi}-\frac{G_{0}}{\lambda}c^{\text{(i)},1}_{1}(\xi)=Ae^{l\xi}+\frac{G_{0}c_{B}}{\sqrt{\lambda}(\sqrt{\lambda+l^{2}}+\sqrt{1+l^{2}})}e^{\sqrt{\lambda+l^{2}}\xi},$
with $A$ given by:
(4.17b)		$A=\frac{G_{0}c_{B}}{\sqrt{\lambda}}\left(\frac{1}{\sqrt{\lambda}+1}-\frac{1}{\sqrt{\lambda+l^{2}}+\sqrt{1+l^{2}}}\right).$

By using the expression of $p_{0}(\xi)$ and $p_{1}^{1}(\xi)$ , (3.29) yields that up to an error of $O(\delta)$ :

(4.18)		$\displaystyle\delta^{-1}\frac{d\delta}{dt}$	$\displaystyle=-\left(\partial^{2}_{\xi}p_{0}(0)+\partial_{\xi}p_{1}^{1}(0)\right)$
		$\displaystyle=\frac{G_{0}c_{B}}{\sqrt{\lambda}}\left(\frac{\sqrt{\lambda}-l}{\sqrt{\lambda}+1}+\frac{l-\sqrt{\lambda+l^{2}}}{\sqrt{\lambda+l^{2}}+\sqrt{1+l^{2}}}\right)=F_{2}(\lambda,l).$

For the in vitro model. Plugging the expansion of $c_{1}(\xi,y)$ (4.4c) into (4.7), together with the conditions (4.8c), for any $k\in\mathbb{N}^{+}$ we have:

(4.19)

c^{k}_{1}(\xi)=a_{1}^{k}e^{\sqrt{\lambda+k^{2}}\xi}\qquad\text{for}\quad\xi\leqslant 0.

And the leading order term $c_{0}(\xi)$ for this case is given by (2.19). Then by using boundary condition (4.8), we get $c^{k}_{1}(\xi)=0$ for any $k\neq 1$ , and for $k=1$ :

(4.20)

c^{1}_{1}(\xi)=-c_{B}\sqrt{\lambda}e^{\sqrt{\lambda+l^{2}}\xi}.

Then similar to the previous case, by solving equation (4.9) together with boundary conditions (4.10) (with $p_{0}$ given by (2.20)), we get $p_{1}^{k}(\xi)=0$ for any $k\neq 1$ . And for $k=1$ :

(4.21)

p^{1}_{1}(\xi)=\frac{G_{0}c_{B}}{\sqrt{\lambda}}e^{\sqrt{\lambda+l^{2}}\xi}.

Finally, by using the expression of $p_{0}$ and $p_{1}^{1}$ , (3.29) yields up to an error of $O(\delta)$ :

(4.22)		$\displaystyle\delta^{-1}\frac{d\delta}{dt}$	$\displaystyle=-\left(\partial^{2}_{\xi}p_{0}(0)+\partial_{\xi}p_{1}^{1}(0)\right)$
		$\displaystyle=\frac{G_{0}c_{B}}{\sqrt{\lambda}}(\sqrt{\lambda}-\sqrt{\lambda+l^{2}})=F_{1}(\lambda,l).$

By now we complete the proof of Theorem 4.1.

4.3. Boundary stability analysis for the two nutrient models

In this subsection, we prove the mathematical properties of $F_{1}$ and $F_{2}$ in Corollary 4.2, which further yields the boundary behaviors in Remark 4.3.

For the in vitro model, one can easily check that $F_{1}(\lambda,l)\leqslant 0$ for any frequency $l\geqslant 0$ . Thus, the amplitude of the perturbation decays as time evolves for any single frequency perturbation. Therefore, the proof of the argument for the in vitro model in Corollary 4.2 is completed.

For the in vivo model, $\delta^{-1}\frac{d\delta}{dt}$ is given by (4.12) and plotted in Figure 1. The limit (4.13b) is obviously for checking. For the asymptote of $l$ approaches $0$ , note that

(4.23)		$\displaystyle F_{2}(\lambda,l)$	$\displaystyle=\frac{(\sqrt{\lambda}-l)(\sqrt{\lambda+l^{2}}+\sqrt{1+l^{2}})+(l-\sqrt{\lambda+l^{2}})(\sqrt{\lambda}+1)}{\sqrt{\lambda}(\sqrt{\lambda}+1)(\sqrt{\lambda+l^{2}}+\sqrt{1+l^{2}})}$
		$\displaystyle=\frac{N(\lambda,l)}{2\lambda(\sqrt{\lambda}+1)(\sqrt{\lambda+l^{2}}+\sqrt{1+l^{2}})}.$

By using the Taylor expansion $\sqrt{1+x}=1+\frac{x}{2}+O(x^{2})$ , when $l$ approaches to $0$ we have

(4.24)

N(\lambda,l)=\frac{\lambda-1}{2\sqrt{\lambda}}l^{2}+O(l^{3}).

Therefore, if $\lambda>1$ then $F_{2}(\lambda,l)>0$ for $l$ close to zero, and $F_{2}(\lambda,l)<0$ for $l$ sufficiently large. Thus, by the intermediate value theorem and the continuity of $F_{2}(\lambda,l)$ in $l$ , there exists $L>0$ such that $F_{2}(\lambda,l)>0$ for $l\in(0,L)$ .

Now, we prove $F_{2}(\lambda,l)<0$ for any $0<\lambda\leqslant 1$ and $l>0$ . From the expression in (4.12), $F_{2}(\lambda,l)<0$ hold obviously for $l\geqslant\sqrt{\lambda}$ . On the other hand, the denominator in (4.23) is always positive. Therefore, it is sufficient for us to show the numerator $N(\lambda,l)$ is negative for $0<l<\sqrt{\lambda}\leqslant 1$ . Indeed, by taking the derivative of $N(\lambda,l)$ with respect to $l$ , we get

(4.25)	$\displaystyle\frac{\partial N(\lambda,l)}{\partial l}$	$\displaystyle=-(\sqrt{\lambda+l^{2}}+\sqrt{1+l^{2}})+(\sqrt{\lambda}+1)+\frac{-l(l+1)}{\sqrt{\lambda+l^{2}}}+\frac{l(\sqrt{\lambda}-l)}{\sqrt{1+l^{2}}}$
	$\displaystyle\leqslant-\frac{l(l+1)}{\sqrt{\lambda+l^{2}}}+\frac{l(\sqrt{\lambda}-l)}{\sqrt{1+l^{2}}}$
	$\displaystyle\leqslant\frac{l(\sqrt{\lambda}-1-2l)}{\sqrt{\lambda+l^{2}}}<0.$

for $0<l<\sqrt{\lambda}\leqslant 1$ . Finally, combine with the fact $N(\lambda,0)=0$ , we can conclude $N(\lambda,l)<0$ for $0<l<\sqrt{\lambda}\leqslant 1$ . By now, we complete the proof of Corollary 4.2.

5. Stability of radially symmetric boundary in the two nutrient models

In this section, we study the boundary stability of the radially symmetric solution under the two nutrient regimes. The structure of this section is arranged as follow. We establish the setups and main conclusions in Section 5.1. After that we carry out the calculation of $\delta^{-1}\frac{d\delta}{dt}$ for the two nutrient models in Section LABEL:sec:Radial_The_detailed_calculations_for_the_two_nutrient_regimes, which serves as the proof of Theorem LABEL:thm:main_thm_for_Radial. Then, we proof the mathematical properties of the perturbation evolution functions summarized in Corollary LABEL:cor:Radial in Section LABEL:sec:Radial_Boundary_stability_analysis_for_the_two_nutrient_models. Finally, we discuss the relationship between the perturbation of the radial boundary and the traveling wave boundary in Section LABEL:sec:Relationship_between_the_Radial_boundary_and_the_traveling_wave_boundary.

5.1. Setup and main results

For the radial case, we employ the polar coordinate system $(r,\theta)$ . Before the perturbation, the tumor boundary is defined by (3.3) with the level set index $Z(t)=R(t)$ , and the tumor region corresponds to the disk with radius $R(t)$ , $D(t)=\left\{(r,\theta)|r\leqslant R(t)\right\}$ . Following the framework of Section 3, we consider the perturbation by a single wave mode, i.e. $\mathcal{P}_{l}(\theta)=\cos{l\theta}$ with $\theta\in[-\pi,\pi)$ . Then, the perturbed boundary (3.5) writes:

(5.1)

\tilde{\mathcal{B}}_{t}(\theta)=\left\{(r,\theta)|r=R(t)+\delta(t)\cos{l\theta},\theta\in[-\pi,\pi)\right\}.

Then the perturbed solutions (drop the tilde) $c$ and $p$ solves (3.9), with the perturbed tumor region $\tilde{D}(t)$ enclosed by $\tilde{\mathcal{B}}_{t}(\theta)$ . Further more, $c$ and $p$ have the asymptotic expansions:


(5.2a)
(5.2b)