Modeling Interactions among Migration, Growth and Pressure in Tumor Dynamics

: What are the biomechanical implications in the dynamics and evolution of a growing solid tumor? Although the analysis of some of the biochemical aspects related to the signaling pathways involved in the spread of tumors has advanced notably in recent times, their feedback with the mechanical aspects is a crucial challenge for a global understanding of the problem. The aim of this paper is to try to illustrate the role and the interaction between some evolutionary processes (growth, pressure, homeostasis, elasticity, or dispersion by ﬂux-saturated and porous media) that lead to collective cell dynamics and deﬁnes a propagation front that is in agreement with the experimental data. The treatment of these topics is approached mainly from the point of view of the modeling and the numerical approach of the resulting system of partial differential equations, which can be placed in the context of the Hele-Shaw-type models. This study proves that local growth terms related to homeostatic pressure give rise to retrograde diffusion phenomena, which compete against migration through ﬂux-saturated dispersion terms.


Introduction
This paper seeks to address the analysis of the interaction among pressure, growth and cell migration in avascular tumors. We propose a thermodynamically consistent approach that links the theory of finite growth with a novel migration term based on flow-saturated mechanisms that allows controlling the dispersion front of the tumor, both from a qualitative point of view defining the characteristics of the front, and quantitatively since the speed of the tumor can be regulated from experimental data.
Cancer is the leading cause of death worldwide. According to data published in [1], 19.3 million new cases were detected and 10.0 million deaths occurred in 2020. Due to the relevance of the issue together with the technical difficulties of obtaining data from in vivo experiments, mathematical modeling appears to be a useful tool to predict and anticipate tumor development. Some papers [2][3][4][5][6] highlight the challenges surrounding the dynamics of cancer and review mathematical models that attempt to respond to these concerns. From one side, research has been focused, on a large degree, on models that control growth and treatments through biochemical interactions that identify morphogens and target genes involved in deregulation and growth associated with tumor processes [4,7]. In addition, the processes of migration and growth entail a reorganization of the collective distribution of cells, even a substantial change in their form when, for example, the epithelial-mesenchymal transition occurs [8]. These processes might produce a substantial change in the elastic properties and aggregation forces and, therefore, in the intracellular and extracellular pressure that affects the cell structure and the environment, both from the biomechanical and biochemical points of view.
In the last few decades, the biology-mechanics interaction has been confirmed as a key player in tumor growth. One of the first mechanical pieces of evidence was supported in 1997 when Helmlinger et al. reported measurements of adenocarcinoma spheroids embedded in agarose gels matrices with different concentrations (different levels of stress). They found that spheroids proliferation was inhibited according to increasing gel concentrations [9]. Further experiments proved stress distribution also affected the shaping and patterning of tumor spheroids [10,11]. Nonetheless, the inhibitory effect that pressure produces on spheroids has been demonstrated to be reversible, disappearing if stress differences vanish [12]. This fact may be explained because pressure seems to cause a quiescent state of cells in which there is a balance between duplication and a certain degree of pressure until stress ceases [13].
Thus, there is clear evidence that cells sense and respond to mechanical forces which regulate biochemical cascade and tumor fate through mechanotransduction [14]. Mathematical models need to consider both biochemical and biomechanics agents that interact in the processes of cellular rearrangement, growth and migration in order to be realistic and to have an accurate predictive character.
The relationship between mathematical models and biology is relatively recent, although prominent authors such as Fick, Fisher, Keller, Segel, and Murray have contributed to drawing attention to the mutual scientific interest for both sciences. This field opens new frontiers to research, particularly in a problem of social relevance such as cancer disease. The mathematical biomechanical approach emerged strongly in the mid-1990s by the name of finite strain theory to consider the growth-stress interplay [15][16][17]. The latest studies have investigated how external pressure of the surrounding tissue could limit tumor evolution through an empirical law, a mechanical function that affects mobility and proliferation [18][19][20]. Moreover, poroelasticity arose to explain in part the cooperation of two main regions of tumor tissue: the solid tumor and the extracellular matrix. This connection has been mainly modeled by Darcy's law, as proposed in [21][22][23], giving rise to Hele Shaw-type models.
Recently, the idea of Hele-Shaw models has been extended to cell competition [24], where it is proposed a poroelastic model based on a porous media equation where pressure regulates both migration and growth. The continuity equation for cell density evolution ρ(t, x), where x ∈ R n , n = 2 or 3 and t ∈ R + , is described by: where ν is the viscosity coefficient, and P = P(ρ) represents the pressure (given, as usual, by the equation of state in terms of the density), which is often used as Darcy's law [25][26][27], and G(ρ, P I ) is the growth term correlated with proliferation and internal pressure, P I . Of course, in (1) there is an underlying dimensioning process to make each of the terms involved consistent. In this sense, depending on the nonlinear flux that we use to define the migration dynamics (in (1) it corresponds to the term F = ρ∇P), the role of the constant ν has different nuances. Our interest will be to adapt the non-linear operator associated with the migration through a flux generically written in form F (ρ, ∇P). Therefore, depending on this non-linear relationship, the constant ν (which generally represents the viscosity or kinematic viscosity) will require an adaptation in its units so that Equation (1) or its modified one remains dimensionless. The growth function has been proposed to be limited by P I , whereas the pressure is homeostatic P h , it is equal to P I . Then, there is a balance between apoptosis and proliferation and G(ρ, P I )| P h = 0. This equilibrium can be local in space, as shown in Section 3. In [25][26][27], a functional framework is proposed for function G to provide a well-posed system of equations. However, modeling G is still challenging. Under these premises, Shraiman proposed a theoretically mechanical feedback mechanism that regulates the homeostatic pressure. In this approach, even if cells were in the exact mechanic environment (same mechanical properties), cells could mutate and grow faster or slower than their surrounding tissue causing cell competition and mechanical stresses. In that case, each cell could compare its proliferation and growth rate with that of the microenvironment and regulate tissue growth [28]. The mechanism proposed to growth function includes proliferation (how cells duplicate) and the velocity of cell growth (the rate at which cells acquire mass): where S(ρ) is the net proliferation of cells and v g the rate of mass acquisition, in which mechanics is considered [29]. However, he also claimed that if cells can move within the tissue, pressure induced by non-uniform growth could be relieved and the mechanism would not work in the whole tissue. The equation proposed by Shraiman for the evolution of the internal pressure takes the form where γ (t, x) is the average growth rate over the tissue, and γ(t, x) is a function depending on the evolution of the density and the dynamics volume or area, which represents the local tumor growth and satisfies γ = αρ [18,19,28]. The area increase rate is α(t, x).
Mechanical constants related to elasticity are κ and µ (bulk and shear modulus), and χ is associated with boundary conditions (see Section 2 for a deduction and an extension of the model). As we will comment below, the above equation leads to a retrograde diffusion term in density (see Section 2.2), a process in which the driving forces are directed towards their source in a movement more typical of aggregation forces than of dispersion. This type of phenomenon causes extreme density concentrations that can produce extrusion or even rupture of the bonds that maintain the aggregated cells (such as those associated with E-Cadherine, among others). Stress and growth phenomena cause a collective reorganization that leads to part of the traction forces generated by each cell being transmitted to the growing cells themselves. This process influences the growth itself and stress is transferred to the border of the growing tissue. However, what happens when movement is added to a migration process? Traction forces that drive the migration of collective cells are generated by some cell units behind the displacement front [30] and have a local influence on their environment, both through biochemical and biomechanical signaling pathways that interact with each other by modifying their respective gradients. Our proposal for migration is to modify the flux that determines the non-linear diffusion free-boundary model (1) by a saturated-flux of the type: where c ρ is the maximum speed of propagation of the solution support and ν ρ the respective viscosity. As we have pointed out previously, the units of the constants ν and ν ρ are not the same and depend on the dimensionless process of the equation, although both are identified as viscosities (or kinematic viscosity in the second case). Migration is a controlled speed process in which a propagation front appears, where pressure plays an essential role since mobile cells communicate the stress state to adjacent cells [31]. Moreover, internal growth pressure continues to generate with a dual effect of aggregation and homeostasis. Then, proliferation, retrograde diffusion, and flux-saturated dispersion are mechanisms that interact and compete, regulating growth and motility.
Flux-saturated models are essentially equations in divergence form so that their flux saturates at a constant value as long as the size of the gradients is large enough. These models frequently appear in some areas of mathematical physics (radioactive transport theory and astrophysics, for example) and are gaining relevance in mathematical biology (morphogenesis and tumor dynamics). This type of model allows us to combine both the diffusion of porous media and flux-saturated mechanisms to obtain a deeper understanding of each of them and their mutual interaction that opens the possibility of new emerging behaviors. The opportunities for this type of flow increase if we also add the growth terms and the influence of pressure (see Section 2.2 below). In the resulting model, new mathematical problems arise that are of great interest and constitute challenges to analyze in the near future. The essential advantage of including the effects of saturated flow compared to considering only the evolution in a porous medium is that saturated flow terms introduce a new biological parameter: the speed of the tumor propagation front, which enables us to control and regulate the profile of tumor progression. This speed can be measured experimentally and allows the propagation front to be defined with more precision, which in the case of porous media is only characterized by internal pressure and not by the characteristics of each agent involved (in our case tumor cells) [32]. This limited speed also affects the evolution of the rest of the terms involved in the process: growth and pressure. Then, the strength of this study lies in joining some of the most relevant items of the current growth models: controlling the tumor propagation front [25,27] and including the role of mechanics [15][16][17]28,29]. In particular, we propose to combine these approaches, incorporating the mechanical feedback of (2) as growth function suggested by the model (1) with a flux-saturated mechanism to explain the tumor front of propagation, but also the internal tumor dynamics and the mechanical feedback phenomena competition through the following system: where A(t, x) denotes N-dimensional volume per unit of the tumor mass. However, the term area is used generically in the paper to describe the two-dimensional form. The carrying capacity is K and β(t, x) is the proliferation rate. The pressure of the flux verifies Darcy's Law P(ρ) = cρ m , with m ≥ 1. An alternative to Darcy's law would be to add second-order terms in the derivative using Brinkman's law, see [33]. Note that the dependency of A with respect to the space variable is parametric and inherited from the dependency of A with respect to ρ(t, x), which is considered as a continuous distribution.
To summarize, in this paper, we analyze the system (3) from the modeling perspective and study the qualitative phenomena resulting from the competition of the mechanisms previously introduced. Section 2 establishes the basis for the deduction of the model from elastic energy minimization principles and analyzes the mathematical properties. Finally, in Section 3, we study various cases of interest with the help of high-resolution numerical methods typical of systems of conservation laws.

Internal Pressure and Growth
Growth is defined here as the process of adding mass and area that occurs through cell division. Although soft tissues have been commonly modeled as hyperelastic materials [34][35][36], the aim of this study does not focus on the constitutive model of the tumor rather the interaction between its mechanical parameters and speed of propagation. Specifically, we have simplified finite growth theory by assuming infinitesimal strains, infinitesimal displacements, elastic, isotropic and linear behavior of the growing tissue. Even though the main equations are also formulated in tensorial notation, the index notation is the one adopted.
Consider a continuous medium at initial configuration B o , which after growth and move lead to the configuration B t at time t ∈ R + . We refer to the total strain tensor that tissue experiments as ε ij . Similar to thermoelasticity, the strain tensor is an additive decomposition (4) of growth deformation ε g ij and mechanical deformation ε e ij , which guarantees mechanical equilibrium accounting for grown tissue [18,19,29,37].
Mechanical deformation is described by the Green-Lagrange tensor (also called Cauchy's or small strain tensor) ε e ij = 1 x) has been previously defined in Section 1.
The following process might be useful as a framework for other linear constitutive laws in small strains. For non-linearities, we refer to finite growth theory studies [15][16][17]29,38]. We have adapted the method reported by [28] to include logistic proliferation function and to account for flux-saturated dispersion. In this work, we have assumed the classical linear elastic isotropic energy: The first term on the right represents pure shear, and the second term is the hydrostatic compression. Analogously to the thermoelastic problem [37,39], the volumetric term is affected by local tumor growth γ(t, x). Mechanical parameters used are the shear µ and the bulk κ modulus [37]. Strain energy is widely accepted in mechanics in its tensorial form that, including growth, can be written as: Constitutive equations are obtained by deriving strain elastic energy function from strains. Then, the Euler-Lagrange associated with the stress-strain relations leads to or in its tensorial form where S refers to the Second Piola-Kirchoff stress tensor. In small strains, S σ e σ e σ e , being σ e ij = (σ σ σ e ) ij and σ e ij the Cauchy stress tensor.
According to hydrostatic compression, pressure is affected only by the volumetric elastic deformation, P I = κε e kk [37]. Taking into account the elastic strain tensor (4), then the internal local pressure is: Furthermore, the relative area changes relate to the local growth rate through the strains (see (7) below), which are derived from the minimization of the strain energy. The process is accurate as long as the time scale of growth is slower than the elastic response of cells [28], and it is particularly useful to study the dependence between mechanical and dispersive biological parameters.
In the mimimization of the functional (5), boundary conditions must be included. These conditions lead to the emergence of an additive scalar term χ 1 (x) as follows which satisfies ∂ 2 χ 1 ∂x 2 = 0. If the tissue is unbounded, then it is deduced that χ 1 = 0. Combining Equations (6) and (7), the local pressure could be expressed by For uniform growth rate and free boundary conditions, the elastic response of tissue readjusts tissue deformation growth by uniform dilatation without inducing additional stresses. With non-homogenous growth, a cluster of cells could grow slower or faster than its surroundings cells causing the cluster strain or compression and the surrounding compression or strain. Depending on stress level, cells may lose their adhesion and become disaggregated from tissue ( Figure 1a) even they could cause buckling in adjacent cells as well (Figure 1b). Then, considering the growth of the surroundings cells, pressure is described as follows: in which we have modified the corresponding expression of the boundary term by using now χ, which includes the new boundary conditions related to the environment. There is also the possibility of restricting the influence of the internal pressure using a relaxation term: The relaxation term limits the growth of temporary "memory" [28], characteristic of the rearrangement process, where the coefficient τ may depend on particular cell adhesion proteins.
The question now is how to modulate this local average growth rate γ (t, x). Assuming that cytonemes (membrane nanotubules) not only mediate cellular communication [40][41][42][43][44][45] but additionally mechanosense the closest cells, it is reasonably accepted that cells feel the pressure and growth of adjacent cells differently depending on how they are connected. Cytonemes usually have an extension of approximately 3-7 cell units [44] (about 60 µm) and cell-to-cell communication is mediated by cytomeme-cytoneme or cytoneme-membrane interaction as presented in Figure 2. We have applied this theory of cell communication to model the local interaction of growth and its influence on internal homeostatic pressure. Then, the average growth is a regular function for each cell, which has support on the range within cytonemes operate: where φ is a regularizing function with compact support on (− , ), being the cytoneme zone of influence and sensitivity of each cell to its environment.
Taken into account the non-uniform growth proposed in (8) and the definition of P I in terms of strains (6), the relative deformation is: When we work with problems defined in all the space, we must require that the trace of the tensor strain tends to zero at infinity, which leads to the boundary term being zero, χ = 0, without more than integrating the previous equation. Given that area change is proportional to the divergence of cell growth velocity v g [29] and the trace of the strain rate tensor, the above expression can be reformulated as: Assuming the previous definition of local tumor growth γ(t, x), the evolution of area reads: Regarding the mechanisms of proliferation, we have considered that the number of cells N divide at rate β(t, x) by means of a logistic function until the carrying capacity K of the medium is reached: where, again, the parametric dependency of N on the space variable is inherited from the dependency on β. Although the proliferation rate has been widely described as constant in many models, it is well known that not all parts of a tumor proliferate equally. In many cases, the rims related to the propagation fronts have a higher proliferation ratio [31].
Therefore, we have considered β variable both spatially and temporally. The number of cells relates to cell density by the relation ρ(t, x) = N A −1 . Then, deriving this relationship of (10) and (11), the equation for density giving rise to The above deduction is independent of the spatial dimension considered.

Some Mathematical Properties of the Growth Model with Internal Pressure
In this subsection, we study and try to give meaning to the problem of evolution (10) and (12).
Assume first that β, κ, µ and K are positive constants, and α is a positive and bounded measurable function. The integral operator is which is well defined for any measurable and bounded function f , where φ (x) = −N φ x , being > 0 the size of the local distance of influence of a cell on its environment and φ is a measurable, positive and bounded function with integral one, supported in the unit ball in R N . Let L ∞ m (R N ) be the vector space of the measurable and bounded functions endowed with the norm which is a Banach space. Note that unlike L ∞ , we have not taken the essential supreme and, therefore, we are not going to identify the functions that are the same almost everywhere. Furthermore, the subscript m highlights the hypothesis that the functions of this space are also measurable. With this definition, it becomes evident that the operator · : L ∞ m (R N ) → L ∞ m (R N ) is continuous, as can be deduced from the classical property of convolution For each ρ 0 , A 0 ∈ L ∞ m (R N ), the initial value problem (10) and (12) associated with has a unique maximal solution, for t ∈ I w , I w being the maximal interval of existence. This result follows directly from the application of the Picard-Lindelöf theorem in Banach spaces for which is locally Lipschitz in L ∞ m (R N ) 2 . Furthermore, because of the continuity in t of the system, function (t, x) → (ρ(t, x), A(t, x)) is a differentiable function in the variable t, measurable in x, and the system of Equations (10) and (12), is verified at all points. This is the basic idea of the proof of the following result in which anx is fixed and the existence in the variable t of the resulting ordinary differential equation can be proved.
• A(t,x) = 0, for somet ∈ I w . Then, A(t,x) = 0, for every t ∈ I w , where I w is the maximal interval of existence.
This result has two important consequences.
Then, the solutions of the system (10) and (12) have a maximal interval of existence I w that contains the entire positive real line [0, ∞) and verify for any t ≥ 0.
Proof. From the previous lemma, we have that ρ(t, x) ≥ 0 and A(t, x) ≥ 0, for t ≥ 0, t ∈ I w . We use the fact that for any ρ ≥ 0. This implies that ∂ ∂t ρ ≤ βρ, from where the first estimate is complete. For the second, observe that As a consequence of the Lipschitzianity of the functions of a variable included in T , the operator is Lipschitz in bounded sets of L ∞ m (R N ) 2 . Then, as in the theory of R N the following result can be proved Lemma 2. Let B be a Banach space and f : B → B a Lipschitzian functional in bounded sets. Let I w = (w − , w + ) and z : I w → B a maximal solution of z = f (z). Then, if w + < ∞, At w − the result is symmetric.
As a consequence of the Lemma 2 and the estimates obtained previously, it follows that I w = (w − , ∞).
We are now going to introduce a result that specifies the domain of definition of the solutions.

Proposition 2.
Let ρ 0 ,ρ 0 , A 0 ,Ã 0 ∈ L ∞ m (R N ) such that ρ 0 −ρ 0 and A 0 −Ã 0 are zero almost everywhere. Let ρ, A andρ,Ã be the corresponding solutions of the initial value problem in I w and I w , respectively. Then, there exists a full measure set R such that for any (t, x) ∈ I w ∩Ĩ w × R N .
Returning to the general context of Lemma 2, it can be verified that in every compact interval J, where the solution of z = f (z), z(0) = z 0 is well defined, each sequence of approximate solutions defined by a polygonal line, called the Euler polygon, with partition norm tending to 0 converges in norm and uniformly in J to the corresponding solution of the problem.
Proof. Let J be a compact interval containing t = 0 that is common in the definition of ρ,ρ, A,Ã. It is immediately verified that given a partition ∆ of J the corresponding polygonal verify ρ Therefore, the result of the convergence of the polygonal (Euler) lines, and taking into account that the uniform limit of measurable functions is measurable, completes the proof, we refer to [46] for more details.
Let ρ 0 ≥ 0 and A 0 ≥ 0. Returning to the relationship between the area and the density, As is well known for the logistic equation, N(t, x) converges point-wise in x towards the carrying capacity K, when t → ∞. We define the support of the initial condition We have the following result then N(t, x) tends uniformly and exponentially towards the stationary value K.
Note that if ρ 0 or A 0 has a jump discontinuity atx, then ρ(t,x) or A(t,x) maintain a discontinuity for small t, respectively. In the same way, if the initial data ρ 0 and A 0 are continuous but α has a jump discontinuity, then ρ(t, ·) and A(t, ·) keep the discontinuity for t small.
Let us now see in this growth and internal pressure model ( (10) and (12)) how the tissue area evolves. To do this, we are going to make the change of variable B(t, x) = 1 Assume that inf{α(x); x ∈ J} > 0.
Then, under hypotheses (14) and (15), αN is uniformly also positive. Therefore, we As a consequence, the following property A(t, x) −→ +∞ holds uniformly in J, provided that (14) and (15) are verified. Thus, under these conditions, the density decreases as since N(t, x) is uniformly bounded in J (see Proposition 3). Note that C could be very small and, therefore, there may be a transition time until the denominator of the previous estimate begins to be greater than one.
We are going to see intuitively how the equation for ρ is equivalent to a retrograde diffusion equation, which involves a process of aggregation and concentration, as well as the propagation of singularities. To do that, we consider α is constant and make a Taylor expansion of ρ(x − y) in the term corresponding to the mean value · . For simplicity, we consider the φ function, which determines · in (13), defined by where B N is the unit ball centered at zero in R N , |B N | is its measure, and ℵ B N (x) is the usual characteristic function on B N . Then, for every regular function f , we have From here, we have deduced 4 ) is the remainder in the Taylor expansion, which depends on the fourth derivative of f . Therefore, we find Consequently, the term of the internal pressure produces a retrograde diffusion effect on the density in (12), while in the equation of evolution of the area (10), the sign is the opposite.

Cell Density Evolution
The aim of this section is to add movement to the cell density since the system (10) and (12), introduced in the previous paragraph, determined solely by growth and internal pressure does not allow cell dispersal. The idea is to analyze the competition between these two mechanisms in opposition with a dispersion term.
Then, the continuity equation of cell density that takes into account motility, proliferation and growth of the cells takes the form ∂ρ ∂t = ν ρ ∇ · J (ρ, ∇P) + ρG(ρ, P I ) (17) where the first term on the right refers to the flux of cells and the second refers to the function that incorporates growth and internal pressure G(ρ, P I ) by means of proliferation S(ρ) and velocity of cell growth v g , which has been obtained in (9). The flux of cells J has been commonly modeled as linear diffusion by choosing the chemical potential J = ∇ ln(ρ). Nonetheless, the core issue of linear diffusion is the questionable control of the front of propagation which produces an immediately invasion of the cells throughout the tissue, causing the loss of a well-defined tumor front propagation [31], which is in contrast with experiments. Several studies have suggested an overcome to this problem through a relaxation term [47,48]. To the best of the author's knowledge, this term is not currently used in oncological mathematical models, in particular, because Rubin demonstrated that the Cattaneo model could violate the second principle of thermodynamics [49].
As we have commented in the introduction, a possibility that provides a finite speed of propagation in the cell spreading is to choose Darcy's law in J to model the pressure. This law leads to diffusion in porous media and the study of the Hele-Shaw model (extensively analyzed in [24][25][26][27] and some references therein).
The speed of cell propagation is a measurable experimental biological parameter that may be incorporated into the model. We aim to clear control the speed of propagation of the front with the so-called flux-saturated equations, see [32,[50][51][52][53] for an introduction to this type of dispersive systems. The purpose is to combine two non-linear diffusion mechanisms: porous media and the flux-saturated terms, which results in a saturated flow as long as the gradient size is large enough. The proposed flux is defined by the equation: where we have considered a Darcy-law-type model P(ρ) = cρ m , with m ≥ 1, being m an empirical parameter related to the propagation and the medium porosity. The finite speed of propagation that limits the speed of the front and which could depend also on the internal pressure is represented by c ρ . In the previous references, it is proved that the solutions of the flux-saturated equations preserve the migration fronts and provide a qualitative and quantitative fit with the experimental data [31,45], in the sense of maintaining compactness of the support and the finite speed of propagation. It also preserves the possible discontinuities, allowing the appearance of emerging invasion profiles. This class of concepts has often been used in similar contexts, such as cell communication [45,54] and morphogenesis processes associated with the equations of Keller-Segel [55,56]. Finally, to complete our model (17) and (18), we define the growth function expressed in terms of proliferation S(ρ) and area acquisition: which, taking (9) and (10) into account, leads to the system of coupled differential Equations (3)- (18). Thus, we propose a model in which migration based on flux-saturated and porous medium dispersion phenomena compete with (local) proliferation and internal pressure, which regulate growth by retrograde diffusion. This combination of mechanisms is discussed in the next section.

Numerical Treatment of the Model
All the equations have been solved in the one-dimensional case with the corresponding adaptation of the radial case in the superior dimension. For the numerical treatment, we have implemented a self-developed code in Matlab (MathWorks Inc., Natick, MA, USA). Time discretization is based on third order explicit Runge-Kutta method. We have reached 3.8 · 10 6 points at t = 67.5 h. For the spatial discretization of saturated flux we have adapted the fifth order WENO method [57,58] with 1000 points. The need to preserve the tumor propagation front and not introduce numerical viscosity implies that the methods based on porous media operators or those with saturated flux require high approximation orders [32,57]. This high order scheme is also necessary due to the strong discontinuities caused by internal (homeostatic) pressure changes. These have been influenced by its retrograde diffusion process together with the flux-saturated mechanisms. The method provides accurate high order stability resolution keeping non-oscillating, stable and sharp discontinuity transitions. Simulations have been carried out for different ranges of parameters (see Table 1). This work is focused on the interaction between migration, growth, and pressure. For this reason, we consider a range of parameter values to analize computationally the different type of patterns caused by the proliferation-pressure-propagation interaction. In order to asses the proposed model, a set of parameters have been included from the literature: (i) growth rate of human colon adenocarcinoma obtained experimentally [59]; (ii) dynamics in front of glioblastomas [31]; (iii) data of LN229 glioma for mechanical parameters reported by [60], where Elastic Young modulus is considered E ≈ 1KPa. Assuming the Poisson coefficient v = 0.45, the bulk and shear modulus are calculated as follows:

Growth and Internal Pressure Dependence without Migration
Let us start by considering growth without migration. Assuming that α and β are constant parameters, cell density evolution is slowed down by the α parameter effect compared with the simple logistic growth, in which α = 0. We have also found that a higher α implies more pressure effect, which decreases the cell density ( Figure 3a). Interestingly, cells of the front (corner regions of the support) feel the difference of pressure caused by the non existence of adjacent cells, which results in a corner peak in cell density that increases with time ( Figure 3b). If pressure is extremely high, which corresponds with values values α 1 · 10 −1 mm 2 cell −1 h −1 , the system is not able to readjust the stresses and the density is concentrated in the pressure jump zones. This phenomenon is a consequence of the retrograde diffusion process that governs the density dynamics causing an aggregation process when there are pressure differences, increasing this effect with growth processes. All concentrations increase with time, while the density is uniformly bounded for each choice of α, which agrees with the results of the previous section with Equation (16). This fact is manifested in the pressure jump zones or in the corners of the domain (Figure 3), and its biological basis is the homeostatic readjustment of the stresses driven by pressure and growth [10].
To take account of more realistic situation and more logical approach, we have considered heterogenous growth where a cluster of cells grow slower (Figure 3c) or faster than its surroundings (Figure 3d), which also implies a difference of pressure (higher or lower α in the core region). Then, instabilities in density are also generated and the dynamics is not able to assume them by attenuating its effect. Moreover, these instabilities increase over time up to a maximum value (bounded by the carrying capacity) unless the system is endowed with a migration process that can disperse these concentrations or instabilities in density, which will be studied in the next section. These instabilities have been tested by comparing them with a solution in meshes of different sizes, in particular finer than the instabilities due to the pressure difference, verifying that they are not consequences of the numerical schemes used. In the results, we have not scaled with respect to carrying capacity. Note that in this case, the carrying capacity is associated with the product of the density ρ times the area A and, therefore, its growth and decrease can be mutually compensated in relation to the carrying capacity. Figure 3. Growth and internal pressure dependence without migration. (a) Cluster of cells growing slower than its surroundings (α = 1 · 10 −1 and α = 1 · 10 −2 mm 2 cell −1 h −1 in the core and the inner region, respectively). (b) Cluster of cells growing faster than its surroundings (α = 1 · 10 −3 and α = 1 · 10 −2 mm 2 cell −1 h −1 in the core and the inner region, respectively). For the four cases, proliferation rate is assumed constant, β = 3.45 · 10 −2 h −1 . (c) Cluster of cells growing slower than its surroundings (α = 1 · 10 −1 and α = 1 · 10 −2 mm 2 cell −1 h −1 in the core and the inner region, respectively). (d) Cluster of cells growing faster than its surroundings (α = 1 · 10 −3 and α = 1 · 10 −2 mm 2 cell −1 h −1 in the core and the inner region, respectively). For the four cases, proliferation rate is assumed constant, β = 3.45 · 10 −2 h −1 and the figures have the same axis label. The results correspond to a tumor with a radius of 1 cm.

Migration, Growth and Pressure Interactions
In this subsection, we are going to test the complete model (3) in which cells are allowed to move through a porous medium by means of a flux-saturated dispersion term. This dispersion process regulates the effects in the density of the internal pressure shape. The velocity of the tumor front is independent of the parameter or function α, even for high α values; see Figure 4a, in which the speed of the front is compared with different heats of α. However, there are significant qualitative modifications in the evolution of the front as a function of internal pressure. These modifications are reflected in the tendency to lose the convex shape of the tumor density (characteristic of flux-saturated [50][51][52][53]), which, as shown in Figure 4b-d, may not occur when the internal pressure difference is high.
Migration seems to regulate stress differences caused by a cluster of cells growing slower than its surroundings for m = 1, which did not happen in the absence of flux as it is reported in Figure 3. Here, it is observed how mechanical feedback leads to an equilibrium situation of the cluster that can be identified with a quiescent state regulated by internal pressure and local growth (Figure 4b). Remarkably, if proliferation rate β is diminished in time (at t = 30 h and 50 h), then the mechanism is delayed, which means that local homeostatic pressure regulation is achieved asymptotically in time (Figure 4c).
As occurs in the case in which the internal pressure is not considered, the excess growth itself can cause a singularity of jump in the evolution of the density since the proliferation of the tumor support is not capable of assimilating its internal growth rate, see [50], in the case of study of the traveling waves associated with the mechanism of flux-saturation. Therefore, extreme values of growth or internal pressure could lead to eliminate the damping of the instabilities in density, as occurs in Figure 4, or even to the emergence of new jump discontinuities.
Furthermore, the poroelastic parameter m contributes to the regulation of the density redistribution process redirected by internal pressure. Indeed, despite the front propagation speed being limited by c p and is exactly c p for m = 1, the propagation speed decreases with the increase of poroelastic coefficient, which could create cell density instabilities (jump or concentrations) not to be dampened. As we have proved in our analysis (with the same values of the parameters of the case m = 1 of Figure 5), considering m = 2, the flow mechanism is not capable of dissipating the internal pressure differences and their consequent oscillations in density within the tumor. For large values of m, in particular for m > 1 the internal pressure differences could dissipate increasing the speed of tumor propagation c p , not the viscosity ν ρ , although the increase in speed may not be compatible with experimental values [31].

Discussion and Conclusions
In this paper, we have proposed how internal pressure (homeostatic), proliferation, and dispersal migration play a relevant role as competitors regulating the reorganization of tumor cell density in its evolution.
Numerical analysis of the model provides further evidence that increasing internal pressure, which may be caused by growth factors or other forces, triggers a decrease in proliferation. We have also examined the local effect in space and time of these phenomena on tumor dynamics. In the absence of migration, the growth mechanisms combined with those of internal pressure give rise to retrograde diffusion Equation (12) that cause large oscillations in cell density. However, as we have detailed analytically, these oscillations are uniformly bounded and depend on the rates of proliferation, internal pressure, and the carrying capacity of the logistic model.
Our model proposes a migration process governed by a flux-saturated mechanism. This dispersive movement is affected by the pressure in different ways: (i) the rate of propagation of the tumor invasion front is independent of internal pressure and growth, but not of the poroelasticity of the medium given by Darcy's law and (ii) the inner part of the tumor tries to rearrange the stresses until a homeostatic state emerges. With this in mind, the model can reproduce patterns of cell quiescence and shows how this process can be reversed, modifying in time the parameters of the growth and internal pressure functions.
Taken together, our findings suggest that retrograde diffusion mechanisms compete against those of non-linear diffusion through flux-saturation in porous media that attempts to counteract these retrograde diffusion oscillations. Then, the whole system incorporates a rich variety of patterns that reproduce experimentally proven phenomena, but that lacked an adequate mathematical context, which provides valuable possibilities for research, both from an analytical and modeling perspective. For instance, from a strictly mathematical modeling approach, including internal pressure terms alters the tendency towards convexification that flux-saturated processes have [32] and adds great scope of discontinuity in density.
In conclusion, the proposed system can reproduce a large part of the evolutionary patterns associated with tissue growth (tumors in particular) evidenced by the experiments, which adds value to this model. The α parameter that defines the internal evolution of the pressure can be measured experimentally, which also opens the possibility of altering or manipulating it externally (for example using mechanical waves) to control tumor proliferation.
Finally, some broad opportunities start with the analysis of this model, and they will be the subject of forthcoming work. Specifically, it is possible to incorporate the aforementioned processes governed by Brinkman's law and the analytical studies on the qualitative aspects of the solutions (oscillations or non-convexification of shape in density, among others). We are also interested in incorporating in a short time the biochemical interactions with the biomechanics studied here, which will modulate both growth and cell adhesion properties. Within the next few years, the reengineering of tumors (via directed drugs or oncotripsy) is likely to become a challenging component in the analysis of tumor dynamics [61,62]. Therefore, future mathematical models on the current topic are required to throw up many questions about the action and consequences of external agents, especially through the feedback between mechanics and biochemical processes involved, on tumor progression.