Multiphase modeling and qualitative analysis of the growth of tumor cords

In this paper a macroscopic model of tumor cord growth is developed, relying on the mathematical theory of deformable porous media. Tumor is modeled as a saturated mixture of proliferating cells, extracellular fluid and extracellular matrix, that occupies a spatial region close to a blood vessel whence cells get the nutrient needed for their vital functions. Growth of tumor cells takes place within a healthy host tissue, which is in turn modeled as a saturated mixture of non-proliferating cells. Interactions between these two regions are accounted for as an essential mechanism for the growth of the tumor mass. By weakening the role of the extracellular matrix, which is regarded as a rigid non-remodeling scaffold, a system of two partial differential equations is derived, describing the evolution of the cell volume ratio coupled to the dynamics of the nutrient, whose higher and lower concentration levels determine proliferation or death of tumor cells, respectively. Numerical simulations of a reference two-dimensional problem are shown and commented, and a qualitative mathematical analysis of some of its key issues is proposed.


Introduction
Mathematical modeling of tumor growth dates back to many years, as the considerable literature on the subject demonstrates. Recent surveys by Araujo and McElwain [3] and by Preziosi and Graziano [19] highlight the increasing interest by biologists and mathematicians toward the qualitative and quantitative comprehension of this complex phenomenon, with the common aim of devising suitable descriptive tools to predict the behavior of the system and to virtually test the effect of different, possibly new therapies.
As a matter of fact, it is all but an easy task to conceive a mathematical model capable to take into account the wide variety of biomechanical and biochemical factors, in most part intimately correlated to each other, that contribute to the overall dynamics of the system. This is particularly evident if one considers that tumor growth is essentially a multiscale phenomenon, in which subcellular, cellular, and tissue levels are strongly interconnected in terms of cause and effect. In many cases, the reason for some macroscopic outcome relies on the activation of specific internal mechanisms to the cells. These usually trigger particular collective behaviors of groups of cells, resulting finally in an observable effect at the macroscopic scale. In addition, when dealing with living matter, as cells and tissues of human body are, classical laws of Newtonian mechanics are not in principle the only rules that the evolution of the system obeys to, due to some sort of intelligence that makes cells able to self-organize according to possible changes in their state or in the environment they live in.
Without claiming to be thorough, and referring instead to the above-cited papers for a comprehensive review on the state of the art, we briefly outline here some basic facts about mathematical modeling of tumor growth at the macroscopic scale, in order to sketch the context which the present work fits in.
A first thread of models of tumor growth to be mentioned (see e.g., Byrne [12] and the main references listed therein) is based upon a particular geometry of the tumor mass, the so-called spheroid, which essentially consists in a spherically-shaped three-dimensional aggregate of cells, primarily produced in vitro for experimental purposes. In this specific configuration, growth is addressed by focusing on the time evolution of the external radius R(t) of the spheroid. An integrodifferential equation is obtained by equating the time variation of the volume of the spheroid to the overall cell proliferation, under the assumption that the latter is somehow related to the concentration of a chemical, usually oxygen, which provides the cells with the nutrient needed for their vital functions. This way it is unnecessary to explicitly track the evolution of the cell density, and the model is closed by simply adding a stand-alone reaction-diffusion equation for the oxygen concentration, however posed in a domain which changes in time according to the evolution of R(t). The resulting mathematical problem is an integro-differential free boundary problem, which has been proved to admit solutions with suitable regularity properties and to predict, in some proper ranges of the parameters, the evolution of the system toward a steady state characterized by a finite steady growth size of the spheroid (Bueno et al. [10], Friedman and Reitich [16]).
A slight variant of this framework entails the introduction of two regions in the spheroid, namely an outer viable zone in which cells proliferate and a central necrotic core in which cells starve and die due to an insufficient delivery of oxygen from the periphery of the tumor. In this case, specific reaction-diffusion equations for the nutrient concentration are set up in either region, accounting for the fact that living cells take up oxygen from the environment, while necrotic cells only receive a small amount of oxygen by means of diffusion. In addition, an integro-differential equation is written for the evolution of the outer proliferating shell, in the same spirit as the corresponding equation of the previous model. The resulting problem describes once again the evolution of the nutrient concentration across the spheroid, but it is characterized by two free boundaries delimiting the necrotic core and the external periphery of the spheroid, respectively. From the mathematical point of view, the problem has been extensively studied by Cui and Friedman [14], who addressed the existence of a stationary solution and the convergence to it of the time dependent model under suitable assumptions on the parameters.
An alternative class of mathematical models of tumor growth relies instead on the theory of deformable porous media. In such a context, the tumor is regarded as a mixture of various mutually interacting components, which obey standard coupled mass and momentum balances. One of the main novelties with respect to the above-described spheroid models is that now no particular geometry is in principle required to write the equations, so that the global conformation of the tumor is in turn an unknown of the problem which can be duly studied by the model itself. However, the assumption of specific geometrical properties, like e.g., spherical or axial symmetry, or even the reference to one-dimensional settings are customary in the literature also in this case, in order to address simplified explorative models.
Mixture-based models usually include one or more state variables specifically devoted to track the evolution in time and space of the cell population. In the early models of spheroids this was not a major issue since, as recalled above, the growth of the tumor was described at an essentially geometrical level by invoking a volume balance. Ambrosi and Preziosi [2] pointed out that this amounts to assuming a constant cell volume ratio within the tumor, so that the growth of the spheroid would be essentially dictated by the necessity to preserve such a constraint in spite of cell proliferation. Therefore, those models can be viewed as particular approximations of a more general multiphase setting. In many cases, the tumor is modeled as a biphasic mixture consisting of cells and extracellular material. Frequently, the former are further divided into several subpopulations of proliferating, quiescent, and necrotic cells, with possible interchanges among the populations determined by the availability of some nutrients that cells need in order to carry out their vital functions. The extracellular material is commonly identified with extracellular fluid or extracellular matrix (stroma). It fills the interstices among the cells, so that no voids are left within the mixture and the latter can be regarded as saturated. In some models of vascular tumors, the extracellular material also includes a distribution of blood vessels mimicking the coopted vasculature under angiogenic effects (Breward et al. [9]).
In the theory of deformable porous media, the velocities of the constituents of the mixture are determined from suitable momentum equations accounting for the internal stress of each phase, as well as for the mutual external interactions among different phases. Therefore, as usual in continuum mechanics, suitable constitutive relations for the components of the tumor have to be specified, in order to render at the tissue level the proper mechanical behavior. It is plain that this poses the major modeling problems, because living tissues can be assimilated to classical materials only on first approximation. For a detailed study of the mechanics involved in tumor growth we refer the interested reader to Ambrosi and Mollica [1]. Here we simply remark that in some particular geometries mass balance equations may be sufficient to determine, at a purely kinematic level, the velocities of the components, at least under suitably reinforced assumptions on the saturation of the mixture (Bertuzzi et al. [5]).
A rigorous derivation of model equations for avascular tumor growth according to the principles of the theory of mixtures can be found in Byrne and Preziosi [11] and in Breward et al. [8]. Following analogous theoretical guidelines, in this paper we are concerned instead with vascular tumor structures, that develop mainly in the axial direction along the wall of a blood vessel. Due to the particular shape they take, they are named tumor cords. For such tumors no particular form of angiogenesis is needed, as they directly catch the nutrient by diffusion from the vessel which they grow around. Mathematical modeling of tumor cord growth has already been addressed by a number of papers by Bertuzzi and coworkers [5,6,7]. In particular, they use a mixture theory approach supplemented by proper kinematic constraints that allow to express the velocity of the components of the mixture without invoking any momentum balance. In addition, they assume cylindrical symmetry of the cords around the blood vessels, further neglecting axial variations of the state variables or performing suitable averages in the axial direction when needed, so as to reduce the domain of the problem to an annulus representing the two-dimensional cross-section of the cord. They also add a sophisticated kinematic constraint to describe the evolution of an inner viable core in interaction with an outer necrotic ring, which are modeled as two distinct zones by means of specific equations for each of them. Finally, they propose a detailed mathematical analysis of both the stationary and the evolution problems generated by their model. Conversely, we are interested in modeling tumor cord growth starting from sufficiently general principles of the theory of mixtures, so as to provide a model able in principle to face many different physical and geometrical configurations. A basic (i.e., minimal) version of the model is developed and studied here. For further extensions and improvements, as well as for more applications to different systems, we refer to Preziosi and Tosin [20].
The paper is organized into six more sections that follow this Introduction. In Sect. 2 the equations of the model and the proper boundary conditions are derived, then in Sect. 3 the mathematical problem is stated in non-dimensional form. In Sect. 4 numerical simulations are shown and commented, with the aim of testing the ability of the model to reproduce some basic relevant features of the system. Sections 5, 6 address the qualitative mathematical analysis of some key issues of the problem, and Sect. 7 finally draws conclusions and briefly sketches research perspectives.

The mathematical model
In this section we propose a mathematical model for the description of the growth of a tumor mass along a source of nutrient, working in the framework of the theory of mixtures. In doing this, one of our main goals will be the minimality of the model, meaning that we will constantly aim at taking into account only those mechano-chemical mechanisms strictly relevant for an essential representation of the macroscopic phenomenon. Of course, we are aware that in this way many interesting aspects, whose inclusion would make the model closer to physical reality and perhaps less academic, are left out of the study. However, we believe that a minimal model along with a related qualitative mathematical analysis are necessary starting points before tackling more complicated situations, and that such a minimal model may be eventually used as the core of future more accurate models.
For the sake of definiteness, we will deduce the model having in mind the two-dimensional geometry depicted in Fig. 1. Actually, since the model is intrinsically multidimensional, any further extension to different geometries and possibly to higher dimensions is mostly technical, requiring in principle the same ideas up to a suitable reinterpretation of the symbolic form of the differential equations that we will write.  Figure 1. The two-dimensional domain of the problem. A spatial region Q is divided into two subregions Ω and Ω c , representing the tumor cord and the host tissue respectively, separated by a free boundary S. The blood vessel coincides with the horizontal z axis.
Let us then consider a spatial region Q ⊂ R 2 occupied by a mixture of cells, extracellular fluid, and extracellular matrix (ECM), that we suppose saturated. Denoting by φ c , φ , φ m the respective volume ratios of the constituents, this condition is expressed by The region Q is internally divided into two subregions Ω, Ω c separated by a one-dimensional manifold S, so that Q = Ω ∪ Ω c and Ω ∩ Ω c = S. We assume that Ω represents the growing tumor cord and Ω c the surrounding host tissue, while S plays in this context the role of a free boundary delimiting the cord. As we will see, the model guarantees that tumor cells located in Ω and normal cells located in Ω c never mix, hence the sole volume ratio φ c is sufficient to track the evolution of the whole cell population. The blood vessel around which the tumor cord develops is schematized by the horizontal z axis, and its presence will be incorporated into the equations of the model in terms of boundary conditions. This is possible, and is indeed useful, because we are not interested in the interactions of the vessel with the tumor, so that its geometry and mechanics can be duly ignored.

Momentum equations.
It is assumed that normal cells and tumor cells differ only in that the latter undergo an unregulated proliferation, while mechanical properties are the same for both. Hence, one can write a unique momentum equation of the form where: • p is introduced as a Lagrange multiplier to satisfy the saturation constraint (1), and is then interpreted as the interstitial pressure of the extracellular fluid; • T c is the excess stress tensor of the cellular matter, accounting for cell-cell internal stress; • m c is the resultant of the actions on the cellular matter due to the interactions with the extracellular fluid and the extracellular matrix. Notice that in Eq. (2) inertia is neglected, in view of the relatively slow dynamics of the cells during their growth process.
Assuming that cells behave like an elastic fluid, we introduce a scalar, possibly nonlinear function Σ = Σ(φ c ) such that I being the identity matrix. Positive values of Σ for high φ c denote compression of the cells, while negative values for low φ c imply a packing tendency. In addition, we postulate the existence of a value φ 0 > 0 such that Σ(φ 0 ) = 0, identifying a configuration in which cells get in touch without Cell stress function Σ(φ c ) and its first order approximation (tangent line) at the cell stress-free volume ratio φ c = φ 0 .
exerting any stress on each other. In the sequel, we will occasionally refer to φ 0 as the stress-free cell volume ratio (see Fig. 2). Concerning the external actions, it is convenient to split m c into two contributions m c , m cm specifically related to the interactions of the cells with the fluid and the ECM. More in detail, we can imagine a simple viscous friction among the components of the mixture, that is where v c , v , v m denote the velocities of the constituents. In particular, in the classical Darcy-like theory of porous media it is customary to express the interaction with the fluid phase as with K > 0 representing the permeability of the mixture. Regarding the interaction of the cells with the extracellular matrix, the theory is instead by far less classical and consolidated, and certainly viscous friction can be possibly accepted only as a first approximation. Indeed, cells are known to attach to the ECM via suitable adhesion sites, and to detach only in presence of sufficiently strong forces (Baumgartner et al. [4], Canetta et al. [13], Sun et al. [21]), therefore such an interaction is definitely not a pure sliding. However, to keep things simple we set for a certain constant Λ cm > 0, referring the interested reader to Preziosi and Tosin [20] for a more detailed mathematical modeling of this term of the equations. Introducing Eqs. (3), (4), (5) into Eq. (2) we find, after some standard manipulations, For the extracellular fluid, an equation similar to Eq. (2) can be written: with m = m c +m m , that is the sum of the interactions with the cells and the ECM, respectively. However, due to the interpretation of p provided above, the excess stress tensor T is commonly neglected and the pressure is regarded as the main internal stress of the fluid, so that Eq. (7) actually simplifies as Relying again on a Darcy-like framework, we think of the interaction terms as proportional to the relative velocities of the interacting constituents, with specifically due to an action-reaction principle, and (10) m m ≈ 0 to emphasize the fact that the interaction of the extracellular fluid with the extracellular matrix is actually negligible with respect to the corresponding one of the cells (cf. Eq. (5)). It might be inferred that a more accurate modeling of the latter may allow for a better expression of the former. However, such topics are beyond the scope of a minimal model of tumor growth and will therefore not be addressed in the present work. Putting Eqs. (9), (10) into Eq. (8) yields the relation which represents the classical form of the well-known Darcy's law.
Finally, an equation similar to Eqs. (2), (7) holds in principle also for the ECM: with m m = m mc + m m = −m cm − m m in view of the action-reaction principle. However, it requires to specify the excess stress tensor T m , which appears to be a very complicated object due to the fibrous nature of the extracellular matrix and to the continuous remodeling it undergoes when cells move within it. Therefore, wishing to focus mainly on the growth of the tumor cord in interaction with a source of nutrient, we choose to consider the ECM as a rigid non-remodeling scaffold in which cells live and grow and extracellular fluid flows. As a consequence, we replace Eq. (12) by (13) v m = 0 in Q and we simply observe that T m plays now the role of a (tensor) Lagrangian multiplier to satisfy the constraint (13). In addition, we take the volume ratio φ m as a constant, say so that the saturation constraint (1) specializes as It is interesting to distinguish in the momentum equation (6) for the cells the contributions due to the interaction with the extracellular fluid, represented by the terms φ c ∇p and m c , from the remaining ones. In most cases, one can assume that the former are negligible, in terms of order of magnitude, with respect to the latter, i.e., (16) so that the main balance is between the intercellular stress and the interaction force with the extracellular matrix. Owing to this, Eq. (6) along with Eq. (13) provides a particularly simple expression for the velocity v c of the cells: By consequence, Eqs. (11), (15) give in turn

2.2.
Mass balance equations. For the system constituted by cells and extracellular fluid, one can also write a balance of mass under the assumption that they form a closed mixture, that is a mixture which does not exchange matter with the external environment. In more detail, we have to take into account that tumor cells, unlike normal cells of the host tissue, are characterized by an intense proliferation, which is mainly responsible for the macroscopic growth of the tumor cord. In order to focus specifically on this phenomenon, we can reasonably disregard the standard reproduction activity of the host cells, so that our mass balance equations read (19) ∂φ c ∂t in Ω c .
In Eqs. (19), G(φ c , c) represents a source/sink of cell mass due to reproduction and death of tumor cells in connection with the availability of some nutrient c. No contribution of ECM to the cell growth is explicitly accounted for, though it is known that the extracellular matrix contains growth factors. The reason is that, in the present context, ECM is not supposed to play a primary role in the dynamics of the system, hence its inclusion in the function G would result at most in constant coefficients related to the volume ratio φ m (cf. Eq. (14)). It is assumed that dead tumor cells are degraded into extracellular fluid, and conversely that the latter is consumed whenever a tumor cell duplicates to originate new cells. We incidentally notice that Eqs. (19), (20) contain the implicit hypothesis that the density of the cells and of the extracellular fluid is the same, and moreover that it is constant in the whole mixture.
χ Ω (t, x) = 1 if x ∈ Ω at time t 0 otherwise, we can substitute Eq. (17) into the first of Eqs. (19) to get a single equation for the cell volume ratio φ c : The fluid volume ratio φ is then determined algebraically a posteriori from Eq. (15) as φ = φ * − φ c . In addition, summing term by term Eqs. (19), (20) and recalling again Eq. (15) we further discover which, together with Eqs. (17), (18), implies This means that, in view of assumption (16), we allow for a nonconstant pressure field p that can be in turn recovered a posteriori as a by-product of the integration of Eq. (21). Finally, once also p is known, Eq. (18) yields the velocity v of the extracellular fluid.
2.3. Cell growth and nutrient diffusion. The function G = G(φ c , c) appearing in Eqs. (19), (21) describes gain (if G > 0) or loss (if G < 0) of cell mass caused by reproduction or death of tumor cells, respectively. These are triggered by the presence of some nutrient within the mixture, whose concentration per unit volume is denoted by c. As a matter of fact, such a nutrient can be mainly identified with the oxygen carried by the blood, which diffuses through the mixture from the wall of the vessel around which tumor cells grow. It is assumed that oxygen molecules occupy no space within the mixture.
It can be inferred that the basic mechanism by which the available amount of oxygen affects the vital functions of the cells consists in the existence of a critical threshold c 0 > 0 of the concentration c. When the latter is above this value, cells are sufficiently fed and can duplicate. Conversely, when c falls below c 0 , cells starve and die. A more sophisticated process may include also a quiescency stage, taking place when oxygen concentration is below c 0 but above a second threshold c 1 < c 0 . In this state, tumor cells neither reproduce nor die, they simply enter a survival state waiting for being reactivated, if c grows above c 0 , or for definitely dying, in case c falls further below c 1 .
¿From the modeling point of view, it is convenient to express the function G in the form , so that the new function g carries the information about the specific growth dynamics of the cells, while the function Γ accounts for growth regulation by oxygen concentration.
We want to point out that a major requirement on the cell volume ratio φ c is that it should be bounded between 0 and φ * , in order to guarantee physical consistency of the solutions returned by the model (cf. Eq. (15)). Cell growth mechanisms play undoubtedly a fundamental role in this respect, and the choice of g turns out to be crucial in order to generate physically significant mathematical models. Therefore, it should rely not only upon physical intuition but also on a qualitative knowledge of the mathematical properties of the model itself. Specific features of g leading to mathematically well-posed problems will be discussed in detail in Sect. 5. Here we simply anticipate that a logistic term (see Fig. 3 is a possible, physically reasonable choice that indeed works.
Concerning Γ, it will be clear at the end of Sect. 5 that it mainly influences the distribution of tumor cells far from the blood vessel at the steady state. Specifically, it determines the maximum size reached by the cord in the transverse x direction. For the moment, we just observe that the two scenarios previously discussed may correspond to respectively (see Fig. 3, right).
As already mentioned, oxygen reaches cells by diffusion from the vessel wall. Actually, to some extent it is also advected through the mixture by the extracellular fluid, but it is known that advection is globally quite inefficient with respect to diffusion, due to the relatively slow dynamics of the constituents of the mixture compared to the high diffusivity of oxygen molecules (Netti and Jain [17]). Therefore we can write for c a standard reaction-diffusion equation of the form where D > 0 is the oxygen diffusion coefficient within the mixture, and the term at the right-hand side models the uptake of nutrient by the sole tumor cells (this is the reason for the indicator function χ Ω ) as proportional to the product φ c c via a phenomenological constant α > 0.

Remark 1.
It is worth noting that in the cell growth function G, as well as in the right-hand side of the first of Eqs. (20), we have not included a term related to the natural death of cells (e.g., a term of the form −δφ c , δ > 0). Similarly, in Eq. (23) neither oxygen absorption by normal host cells nor possible chemical degradation of the nutrient (e.g., by a term like τ −1 c, with τ > 0 proportional to the characteristic half-life of oxygen molecules) has been taken into account. Of course, these effects could be straightforwardly incorporated in the model, but in the present context they seem negligible with respect to the main dynamics we are concerned with.

2.4.
Boundary and interface conditions. The system of equations (21), (23) needs to be supplemented by suitable conditions for the two main state variables of the problem, i.e., φ c and c, at the boundary of the domain as well as on the interface S between the regions Ω and Ω c . With reference to Fig. 1, we see that the domain Q has four boundaries, whose just one is physical, namely the edge l 1 representing the blood vessel, while the remaining three serve uniquely to confine geometrically the domain of the problem. In particular, l 2 can be considered as a symmetry boundary, in the sense that one can imagine a symmetric evolution of the system for z < 0, while the upper and right boundaries l 3 , l 4 should be thought of as 'sufficiently far' (ideally, x, z → +∞) in the host tissue to be unaffected by the dynamics of the growing tumor cord.
In the sequel, the symbol n will denote a generic outward normal unit vector, which has to be referred from time to time to the appropriate boundary or interface under consideration.
At the vessel wall we assume no detachment of both tumor and normal cells, which essentially amounts to setting their normal velocity to zero. Moreover, we prescribe the typical oxygen concentration c b found in the blood. Taking Eq. (17) into account, we have therefore the following conditions: On the symmetry boundary we assign canonical symmetry conditions on φ c and c: which in our specific case take the form ∂ z φ c = ∂ z c = 0 for z = 0. Finally, at the 'far' upper and right boundaries we assume that the host tissue is unperturbed, i.e., unstressed, and we set zero flux of oxygen: Note in particular that the second of these conditions is interpreted in the present context as ∂ x c = 0 on l 3 , ∂ z c = 0 on l 4 . According to the classical theory of deformable porous media (see e.g., Preziosi [18]), at the interface S between the tumor cord and the host tissue suitable conditions on the continuity of the internal cell stress and of the flux of oxygen have to be imposed: where · denotes the jump across the interface S. In particular, given the constitutive relation (3), the condition on the cell stress corresponds to φ c Σ(φ c ) = 0, and, if one assumes that Σ is a continuous function of φ c , this further reduces to φ c = 0, that is the continuity of the cell volume ratio across S. In addition, the motion of S as a material interface for the cellular matter is regulated by the following ordinary differential equation: where v c is given by Eq. (17), which simply accounts for the fact that, in a Lagrangian framework, the points of S move with the normal velocity of the cells themselves.

Nondimensional statement of the problem
Let L, τ , Φ, C denote characteristic values of length, time, cell volume ratio, and nutrient concentration respectively, that we choose such that Let us moreover introduce the non-dimensional variablesx,t,φ c ,c defined by the relations Dropping from now on the tildes from the dimensionless variables and the subscript 'c' from φ c to simplify the notation, the problem originated by Eqs.
plus the free boundary condition The choice of φ * , namely the per cent amount of free space left by the extracellular matrix, as reference value for the cell volume ratio implies that at a non-dimensional level we should expect 0 ≤ φ ≤ 1 in Q in view of Eq. (15), as well as 0 < φ 0 < 1 coherently with the fact that in the dimensional formulation of the model φ 0 is implicitly assumed to be positive and lower than φ * for physical consistency reasons. Similarly, since we expect c b to be the maximum value of nutrient concentration found in the system, we consider 0 ≤ c ≤ 1 in Q.

Numerical results
Equations (26) have been numerically integrated over the reference domain Q = [0, 2.5]×[0, 10] via a standard Finite Element discretization, coupled with a level set technique to track the evolution of the free boundary S. Specifically, the following functions have been used: along with this set of non-dimensional parameters: We refer in particular the reader to Subsect. 5.3 for further motivations on the choice of Σ.  At early times (t = 100) the cord begins to grow keeping a round homogeneous shape, due to the high availability of nutrient across the domain. However, the progressive consumption of nutrient by the cells makes oxygen concentration fall below the critical threshold c 0 at the top of the cord (t = 325), that is in the farthest zone from the blood vessel. As a consequence, the transverse width of the cord starts decreasing (t = 650) and the whole structure takes a more elongated shape. At later times (t = 900), a front and a rear regions can be clearly distinguished along the cord. The first one is a sort of growing rounded head, where most living tumor cells are concentrated because they are globally fed by a sufficient amount of nutrient. The second one is a kind of straight tail of nearly constant transverse width, in which the distribution of cells and oxygen depends uniquely on the distance from the blood vessel and has substantially reached a steady state.
The dynamics of the system depicted by the model highlights two main issues to be investigated, namely the axial and the transverse growth of the cord. In more detail, the former is concerned with the determination of the shape and the growth speed of the head in its motion along the blood vessel, possibly in connection with traveling waves describing the evolution of the front part of the interface S. The latter is instead aimed at characterizing the typical transverse width of the tail, taking advantage of the above consideration that the rear part of the cord is basically stationary with an evident special dependence of φ and c on the space variables x and z. In this paper, specifically in Sect. 5 below, we address this second problem, leaving the first one for a possible forthcoming work.
We begin from a purely descriptive point of view, observing that in the tail of the cord the cell volume ratio attains its maximum value at the vessel wall (x = 0), as it can be expected for in that zone cell proliferation is constantly sustained by a direct nutrient supply. Conversely, at the interface S and in the corresponding upper region inside the host tissue it decreases toward the stress-free value φ 0 . Hence the expansion of the tail toward a transverse stationary configuration causes tumor and host cells to achieve an equilibrium spatial distribution characterized by a null mutual stress as well as a null internal stress of the surrounding tissue. We incidentally notice that an analogous unstressed state cannot be reached at the head of the cord, which is never stationary for external cells are continuously compressed and pushed away by growing tumor cells that penetrate into the host tissue.
With respect to oxygen distribution, we note that the tail can be ideally divided into two stripes: An inner one, extending from the vessel to a certain depth inside the cord, where c ≥ c 0 , and an outer one, reaching the periphery of the cord, in which c < c 0 . It is straightforward to assimilate them to the so-called viable and necrotic regions, respectively, whose existence is experimentally confirmed and explicitly described by means of suitable equations in the mathematical model of tumor cords by Bertuzzi et al. [5]. However, we want to point out that in the present context the formation of these two distinct regions is not postulated a priori as a modeling assumption, but is recovered a posteriori as a result of the specific dynamics entailed by the model itself. Oxygen concentration, which is maximum at its source, namely the blood vessel, decreases while diffusing through the cord due to the absorption by tumor cells. The farther from the vessel cells are the lower the quantity of nutrient they are reached by is, so that some of them starve and can no longer proliferate. The process continues until in a certain portion of the cord the net proliferation rate is zero: Then that portion of cord is in equilibrium, it stops growing transversally and finally reaches a steady state.

The stationary problem
This part of the paper is devoted to the analysis of the mathematical model presented in the previous Sects. 2 and 3. In particular, referring to the nondimensional form (26) of the equations, we consider a steady state in which the cord is supposed to be infinitely long in the longitudinal z direction, and to have grown up to a width w > 0 in the transverse x direction (see Fig. 6). In such a configuration, the state variables φ and c depend only on x in Ω.
In this geometrical setting, the free boundary S coincides with the line x = w, while the tumor cord and the host tissue are represented by the stripe 0 < x < w and the half-plane x > w, respectively. According to Eqs. (26), an admissible steady solution (φ, c) in Ω c satisfies φ = φ 0 , ∇c = 0, hence we can confine ourselves to the region Ω assuming the continuity conditions φ = φ 0 , The objectives of the present study are twofold: (a) As a first step (cf. Subsect. 5.1), to establish existence of physically significant solutions with appropriate regularity. In this stage, the width w of the cord is regarded to all purposes as a parameter of the model. Specifically, we consider the following class of boundary value problems 1 : in the unknowns φ, c :Ī → R + , where: (i) I = (0, 1) is the rescaled domain after the substitution x → wx, which reduces the free boundary problem to a fixed boundary problem. Notice that w appears now among the coefficients of the equations. (ii) F = F (φ) is a 'generalized' stress function linked to the cell stress function Σ = Σ(φ) by the relation where the superscript stands for derivation with respect to φ. This way it results in accordance with the first of Eqs. (26). (iii) g = g(φ) is the specific cell growth function, which determines how cells proliferate on the basis of their current distribution. (iv) Γ = Γ(c) is the growth regulation function, which expresses promotion or inhibition of cellular proliferation as a consequence of the availability of nutrient c. (v) α > 0, φ 0 ∈ (0, 1) are phenomenological parameters related to the consumption rate of nutrient by the cells and to the stress-free cell density, respectively.
1 The condition ∂xφ = 0 for x = 0 of Problem (32) is equivalent to the original boundary condition ∂x(φΣ(φ)) = 0 of Problem (26) in a sense that will be made precise later in Remark 4.

Remark 2.
Given the dimensionless cell stress function Σ, condition F (0) = 0 is readily obtained by properly choosing the integration constant in Eq. (33). Conversely, we observe that the fulfillment of condition F (1) = 1 may simply require a suitable rescaling of Σ.
Assumption 2. We assume that g is Lipschitz continuous on [0, 1], with Lipschitz constant L g > 0: Moreover, we require Assumption 3. We assume that Γ is Lipschitz continuous on [0, 1], with Lipschitz constant L Γ > 0: From Assumption 1 it follows that F is monotonically increasing and invertible on [0, 1], and that its inverse function, henceforth denoted by f for brevity, is in turn increasing and differentiable on every compact subset of (0, 1]. We agree to denote by L f,ε the Lipschitz constant of f on every interval of the form [F (ε), 1], 0 < ε < 1. Notice that the assumptions on F amount mainly to a request of ellipticity of the nonlinear equation for φ, which is however potentially degenerate at φ = 0.
On the other hand, Assumptions 2, 3 entail the boundedness of g, Γ on [0, 1], hence there exist For modeling reasons discussed in Sect. 3, we are interested in nonnegative continuous solutions φ, c bounded from above by 1. In more detail, in view of the possible degeneracy of the equation for φ in zero, we require that φ be strictly positive in the domainĪ, therefore in the sequel it will be instrumental to refer to the following sets: where ε > 0 is a fixed parameter. Notice that we cannot expect solutions φ > φ 0 onĪ, for they would not match the boundary condition of Problem (32) at x = 1, hence we assume ε < φ 0 .
Remark 3. We anticipate that, as far as the solution to the free boundary problem is concerned, the specific value of the parameter ε in the range (0, φ 0 ) is irrelevant, because the cell density φ can be shown to satisfy automatically φ ≥ φ 0 on the whole intervalĪ (cf. Theorem 7). The introduction of ε is due to the necessity to obtain existence in V ε × U of solutions to Problem (32), among which to look later for possible solutions to the free boundary problem. However, the theory we are going to develop will provide as by-product an optimal criterion on whose basis to select the value of ε.
In the following, we will consider the main results of our analysis without going into the details of their proofs, so as to give an overall picture of the mathematical problem and of its solution. The interested reader can find all technical proofs postponed in Sect. 6. We simply point out here that the constant C P , occasionally appearing in several forthcoming formulas, is the Poincaré constant of the domain I.

5.1.
Existence and regularity of the solutions. We begin by addressing the question of existence and regularity of solutions to Problem (32) for fixed positive w. Under suitable assumptions on the parameters of the equations, we will identify a family of solutions (φ, c) ∈ V ε × U attached to each w ranging in an appropriate set of values.
Fix then w > 0. The strategy we adopt to solve Problem (32) consists of the following steps: (1) We choose any ϕ ∈ V ε and study the boundary value problem , establishing existence and uniqueness of a solution c ∈ U .
(2) We then choose any σ ∈ U and study the boundary value problem , finding conditions on the parameters that guarantee existence and uniqueness of a solution φ ∈ V ε .
(3) As a consequence of the previous steps, we can define an operator A : V ε → V ε that to every function ϕ ∈ V ε associates the solution φ ∈ V ε of Problem (35) through the solution c ∈ U of Problem (34). By showing that A satisfies the hypotheses of Schauder Fixed Point Theorem, we finally find a function φ ∈ V ε , and consequently also a function c ∈ U , solving Problem (32) as desired. (4) At last, we prove that under the same conditions on the parameters established in the previous steps it is possible to control the distance between any solution φ ∈ V ε to Problem (32) and φ 0 , so that a perturbative expansion of φ about φ 0 be justifiable in addressing next the free boundary problem.

Proposition 1.
To every ϕ ∈ V ε there exists a unique solution c ∈ U of Problem (34), which is further twice continuously differentiable and monotonically decreasing onĪ, with a maximum value equal to 1 attained for x = 0.
Step 3: Solution to Problem (32). We will henceforth always assume βw < 1, even when not explicitly stated. Let us introduce the operator A 1 : V ε → U that to every ϕ ∈ V ε associates the unique solution c = A 1 (ϕ) ∈ U of Problem (34). Analogously, let A 2 : U → V ε be the operator that to such c ∈ U associates the unique solution φ = A 2 (c) ∈ V ε of Problem (35) with σ = c. By composition, we define the operator A = A 2 • A 1 : V ε → V ε that to every ϕ ∈ V ε associates the unique solution φ ∈ V ε of the problem then the pair (φ, c = A 1 (φ)) ∈ V ε × U is a solution to Problem (32). Our task is therefore to show that A admits such a fixed point.
For this, it is first useful to know that Moreover, since in applying fixed point techniques a key feature is usually the compactness of the involved operator, the following result is particularly welcome.

Proposition 4. The operator
Thanks to Propositions 3, 4, we can apply Schauder Fixed Point Theorem to solve Problem (32). For the sake of completeness, we explicitly recall in the statement of the theorem all hypotheses that bring to the result.
Theorem 5 gives existence but not uniqueness of solutions (φ, c) ∈ V ε × U to Problem (32). Observe however that, owing to Proposition 1, to any fixed point φ ∈ V ε of the operator A there corresponds a unique c = A 1 (φ) ∈ U such that the pair (φ, c) ∈ V ε × U solves Problem (32).
Step 4: Controlling the distance between φ and φ 0 . For any solution (φ, c) ∈ V ε × U to Problem (32), we are interested now in controlling the distance between φ and the constant φ 0 , in view of the application of a perturbative technique to solve the free boundary problem. We find out that the quantity β 2 w plays a fundamental role, in fact we have: Theorem 6. Let φ ∈ V ε be any solution to Problem (32). Then ¿From this result we infer that, since β 2 w ≤ βw, for even moderately small values of βw (with respect to 1) the solution φ becomes close to φ 0 provided g M (maximum cell proliferation) and L g (maximum proliferation rate) satisfy g M /L g = O(C P ), C P 0.64 (cf. Sect. 6).  To this end, we formulate additional assumptions on Γ: Assumption 4. In addition to Assumption 3, we assume that Γ is nondecreasing on [0, 1]: Notice that Assumptions 3 and 4 imply the existence of a point c 0 ∈ (0, 1) such that Γ(c 0 ) = 0, with moreover Γ(ξ) ≤ 0 for ξ ∈ [0, c 0 ), Γ(ξ) > 0 for ξ ∈ (c 0 , 1].
We state first of all a crucial result that clarifies the role of the parameter ε.
In view of Theorem 7, the specific value of ε ∈ (0, φ 0 ) is not relevant to solve the free boundary problem, indeed any possible solution, provided it exists, is 'naturally' bounded from below by φ 0 . Considering that Theorem 5 guarantees existence of solutions if βw < 1 and that, on the other hand, it results β → +∞ when both ε → 0 + and ε → φ − 0 (cf. Proposition 2), we deduce that the optimal criterion for the choice of ε aims at minimizing β(ε), so as to get the widest possible range of admissible values of w.
Let us look for a first order expansion of φ in a neighborhood of the constant stress-free value φ 0 : , where ν > 0 is a 'small' parameter to be identified from Problem (32). Correspondingly, also c is expanded with respect to ν in a similar way: . Inserting Eqs. (40), (41) into the second of Eqs. (32), and retaining only the terms up to the zeroth order in ν, yields the following problem for c (0) : where we have emphasized the dependence of c (0) on w as a parameter. It is an easy task to check that c (0) w ∈ U , indeed Problem (42) falls into Proposition 1 with the special choice ϕ ≡ φ 0 ∈ V ε . In order to get an equation satisfied by φ (1) , and to detect at the same time the parameter ν, we proceed in a similar fashion, introducing expressions (40), (41) into the first of Eqs. (32). Assuming formal differentiability of g, Γ we find, after some standard algebra, . As Γ is not identically zero, another term of zeroth order in ν should be found in this equation, besides the one at the right-hand side. Estimating the order of magnitude of Γ(c As anticipated above, the parameter ν must be 'small' in order for the perturbative expansion to be meaningful. Equation (44) shows that this requirement corresponds in essence to the same condition prescribed by Theorem 5 for the existence of solutions to Problem (32). Therefore, the existence theory and the perturbative approach are consistent with one another, as they impose the same conditions on the parameters β, w of the model. In other words, for a certain value of the quantity βw, they either apply or fail both.
Retaining only the terms up to the zeroth order in ν, the problem for the perturbation φ (1) reads: is an appropriate positive constant. In particular, boundary conditions for φ (1) are deduced from those prescribed on φ taking Eq. (40) into account.
Unlike the previous case for c (0) w , Problem (45) cannot be solved in a closed form for a generic function Γ. However, classical theory for linear elliptic equations guarantees existence and uniqueness of a solution φ (1) ∈ H 1 0,1 (I), which is actually continuous due to Sobolev embedding theorems for I ⊆ R. Furthermore, since c (0) w ∈ U and Γ is continuous on [0, 1], we even get ∂ 2 x φ (1) ∈ C(Ī), hence φ (1) is a classical solution to Problem (45). This allows us to evaluate its derivative for x ∈Ī by integrating once the differential equation in (45): whence, owing to the boundary conditions for x = 0, x = 1, Solving the (approximate) free boundary problem corresponds therefore to finding the roots w > 0 of Eq. (46). After some preliminary considerations about the dependence of c (iii) For all c 0 ∈ (0, 1) there exists w * > 0 such that if w ≥ w * the equation c (0) w (x) = c 0 in the unknown x has exactly one solutionx w ∈Ī. (iv) For any fixed c 0 ∈ (0, 1) one hasx w → 0 when w → +∞.
we have all we need to solve our free boundary problem. The following result should be compared to the analogous one proved by Bueno et al. [10]. Notice that Theorem 9 does not guarantee that the solution w 0 actually satisfies βw 0 < 1. In practice, this condition has to be checked a posteriori, after solving Eq. (46) explicitly. Concerning this, we observe that, once the function Γ has been specified, Eq. (46) often specializes in an algebraic equation in the unknown w, that can be solved with arbitrary precision, possibly with the aid of suitable numerical techniques. After getting the solution w 0 claimed by Theorem 9, two situations may arise: Condition βw 0 < 1 is either satisfied, hence, in view of the reasonings above, one can take w 0 as a good approximation of the solution to the free boundary problem, or it is violated, in which case one should reject w 0 and conclude that, for the specific model parameters at hand, no information on the existence of a solution satisfying the free boundary condition can be obtained from the present theory.
We show an example of application of this procedure in the next Sect. 5.3.

A noticeable example.
In order to illustrate how the previous theory applies to a specific model, we consider for µ ≥ 1 the following class of intercellular stress functions Σ (cf. Eq. (28)): whose behavior for several µ is illustrated in Fig. 7. Notice that for µ > 2 the function Σ grows unboundedly with φ, and is thus compatible with the expected behavior of an intercellular pressure-like stress, while for 1 ≤ µ ≤ 2 it is bounded from above, and if µ = 2 it even tends to zero for φ → +∞. From the modeling point of view, the meaningful cases correspond to µ > 2. However, the theory developed in the previous sections also covers the range of values 1 ≤ µ ≤ 2, which may be of some theoretical interest due to the particular form of the equations they originate (for instance, for µ = 1 the differential part of the equation for φ is linear). According to Eq. (33), the corresponding generalized stress function F is which gives rise to the stationary porous medium equation for the cell volume ratio φ. It is straightforward to check that such an F complies with Assumption 1 for all µ ≥ 1. Moreover, we choose for g and Γ the expressions given by Eqs. (29), (30): where 0 < c 0 < 1 and γ > 0 are constant, whence also Assumptions 2, 3, 4 are satisfied.
The inverse function f of F on [0, 1] is Notice that f is actually Lipschitz on every compact subset of (0, 1] for all µ ≥ 1, but for µ > 1 it is not Lipschitz on the whole interval [0, 1] due to the vertical tangent at φ = 0. However, we know that the interesting parameter is its Lipschitz constant L f,ε on [F (ε), 1], which exists for any µ ≥ 1 and equals Moreover, For all w > 0 such that βw < 1, that is (48) w < min 1 Theorem 5 guarantees the existence of a solution (φ, c) ∈ V ε × U to the problem Let us now look for a solution to the free boundary problem (39) via the small parameter approximation discussed in Subsect. 5.2. The specific form of the function Γ we are using allows for an easy estimate of the solutions to Eq. (46). Inserting the expression of c For the sake of definiteness, let us fix the same parameters already used in Eq. (31), then let us determine the optimal ε by solving the min-max problem This gives ε 0.5 and correspondingly β 1 β 2 0.55, whence it can be easily computed that the estimate (51) agrees with the previous bound (48). Solving Eq. (50) via an iterative numerical procedure yields w 0 1.45. Since βw 0 0.80, the small parameter assumption and the related perturbative approach are consistent with the problem at hand, as it is definitely confirmed by the value of w returned by the Finite Element solution (see Fig. 4, 8) of the exact problem, i.e., w 1.44. The behaviors of the relative errors on φ and c (see Figure 9): where φ (1) is analytically computed for this particular case from Eq. (45) using C = Γ −1 M , and where we have set ν = (β 2 w) 2 , further prove the good quality of the small parameter approximation with respect to the exact solution.

Technical proofs
Here we collect, without specific comments, the proofs of all results stated in the previous Sect. 5. Before going into technical details, however, it is convenient to fix some basic facts and notations that we will extensively use in the sequel.
We denote by H 1 0,0 (I) (H 1 0,1 (I), respectively) the closed linear subspace of the Sobolev space H 1 (I) = W 1,2 (I) consisting of all functions u ∈ H 1 (I) whose trace vanishes at x = 0 (x = 1, respectively). We recall that Poincaré inequality holds true for functions u ∈ H 1 0,0 (I), or u ∈ H 1 0,1 (I): where the Poincaré constant C P equals 2 π . As a consequence, the quantity u H(I) := ∂ x u L 2 (I) defines a norm in both H 1 0,0 (I) and H 1 0,1 (I) equivalent to the usual H 1 -norm. Notice that in referring to such a norm we use generically the short subscript H(I) for either of the spaces, the meaning being recoverable each time from the context.
When dealing with continuous functions u ∈ C(Ī), we indicate by u ∞ their ∞-norm overĪ: Finally, we use the symbols u + , u − for the positive and negative part of a function u: In view a theorem due to G. Stampacchia, it is known that u ∈ H 1 (I) implies u + , u − ∈ H 1 (I) with moreover Let us start by a mainly technical Lemma, which will be sometimes referenced in the forthcoming proofs.
Lemma 0. Let h, k ∈ C(Ī) with moreover h(x) ≥ 0 for all x ∈Ī. If u ∈ C(Ī) solves the equation in I with either of the following sets of boundary conditions: Proof. Using standard techniques for linear elliptic problems, it can be shown that the proposed equation admits a unique solution u ∈ H 1 0,0 (I) in case (i), and a unique solution u ∈ H 1 0,1 (I) in case (ii). In both cases, u ∈ C(Ī) due to Sobolev embedding theorems for I ⊆ R. Moreover, Lax-Milgram Theorem entails the following a priori estimate: which, considering that and, from Morrey inequality, that u ∞ ≤ u H(I) , yields u ∞ ≤ k ∞ and thus the thesis. Proposition 1. To every ϕ ∈ V ε there exists a unique solution c ∈ U of Problem (34), which is further twice continuously differentiable and monotonically decreasing onĪ, with a maximum value equal to 1 attained for x = 0.

Proof.
(1) By means of the substitutionc = c−1, Problem (34) is converted into the following linear elliptic problem with homogeneous boundary conditions: (2) We check now that c ≥ 0 onĪ. For this, observe first that c − ∈ H 1 0,0 (I), then multiply Eq. (34) by c − and integrate by parts over I: We deduce c − H(I) = 0, which yields c − = 0 almost everywhere in I and, by continuity, c(x) ≥ 0 for all x ∈Ī.
(3) Finally, we show that c ≤ 1 onĪ. Since c ∈ C(Ī), the second order derivative ∂ 2 x c exists as a distribution. But from Eq. (34) we discover that it actually coincides with a continuous function: x c = αw 2 ϕc ∈ C(Ī), so that c turns out to be a classical solution of Problem (34). This allows us to compute, for every x ∈Ī, whence, using the boundary condition at x = 1, Therefore we have that c is nonincreasing onĪ, and consequently we deduce max x∈Ī c(x) = c(0) = 1, which completes the proof.
(1) First we rewrite Problem (35) making the substitution which in turn implies φ = f (u + u 0 ) for u 0 := F (φ 0 ) ∈ (0, 1): Notice that for φ ∈ V ε the boundary conditions at x = 0 of Problems (35) and (54) agree, (2) Now we prove that for all w > 0 Problem (35) The functiong(u) being continuous, if P (u) denotes any of its antiderivatives on R then Problem (55) can be viewed as the Euler-Lagrange equation associated to the energy functional E : H 1 0,1 (I) → R, Due to the assumptions on the sign of g and the monotonicity of f , it resultsg(u) = 0 for u ≤ −u 0 or u ≥ 1 − u 0 , hence P (u) is bounded: There exists a constant C > 0 such that |P (u)| ≤ C for all u ∈ R. This entails first that E is finite on H 1 0,1 (I): E(u) ≤ 1 2 u 2 H(I) + w 2 CΓ M < +∞, and moreover that it is coercive: where C > 0 is a new suitable constant. Thus any minimizing sequence is bounded in H 1 0,1 (I), and has therefore weakly convergent subsequences. If we can show that E is sequentially weakly lower semicontinuous on H 1 0,1 (I), we will deduce the existence of a minimizer u ∈ H 1 0,1 (I), that is a solution to the auxiliary problem (55). For this, notice that E belongs to the class of functionals of the form where L : R 2 × I → R is the function given by Clearly, L is convex in p for each z ∈ R, x ∈ I, and moreover it is bounded from below as it results L(p, z, x) ≥ −w 2 CΓ M for all (p, z, x) ∈ R 2 × I. These two features are sufficient in order for E to be sequentially weakly lower semicontinuous on H 1 (I) (see Evans [15], Chapter 8, p. 446 for further details), hence also on H 1 0,1 (I) as desired. We conclude that the auxiliary problem (55) has solutions u ∈ H 1 0,1 (I) which, as a consequence of Sobolev embedding theorems for I ⊆ R, are actually continuous functions u ∈ C(Ī).
(b) We claim that any solution u ∈ H 1 0,1 (I) to Problem (55) satisfies the inequalities −u 0 ≤ u(x) ≤ 1 − u 0 for all x ∈Ī. To see this, let us multiply Eq. (55) by (u + u 0 ) − ∈ H 1 0,1 (I) and integrate by parts over I: Sinceg ( We observe now that for u ∈ [−u 0 , 1 − u 0 ] one hasg(u) = g(f (u + u 0 )), thus any solution to the auxiliary problem is also a solution to Problem (54). Therefore, the latter admits solutions u ∈ C(Ī) bounded between −u 0 and 1 − u 0 , which, via Eq. (53), originate solutions φ ∈ C(Ī) to Problem (35) such that 0 ≤ φ ≤ 1 onĪ. (3) Finally, we show that it is possible to get solutions φ ∈ V ε to Problem (35) provided w is sufficiently small. For this, let φ ∈ C(Ī), 0 ≤ φ ≤ 1, be a solution to Problem (35). Since F (φ 0 ) is constant, we can rewrite the equation for φ in the equivalent form Thus, multiplying both sides of the above equation by F (φ) − F (φ 0 ) and integrating by parts over I, we discover where we have used the fact that, for φ ∈ [0, 1], the function g is bounded by g M . In view of Morrey inequality we further deduce .
(4) Regarding uniqueness, suppose, under the hypothesis β 1 w ≤ 1, that φ 1 , φ 2 ∈ V ε are two solutions to Problem (35). Subtracting the respective equations we find Multiplying both sides by F (φ 2 ) − F (φ 1 ) and integrating by parts over I we get where we have used the Lipschitz continuity of g on [0, 1]. Setting for all x ∈Ī, and the uniqueness follows from the bijectivity of F . (5) In view of the previous results and given the definition of β, we conclude that if βw < 1 then there exists a unique solution φ ∈ V ε to Problem (35) as desired.
Proposition 3. The operators A 1 : V ε → U , A 2 : U → V ε are Lipschitz continuous, and so is Proof.
(1) Continuity of A 1 is implied by Proposition 3, hence in order to have compactness we need to prove that A 1 (V ε ) is relatively compact in U . We do it by showing that for any {ϕ n } ⊂ V ε the sequence {c n } = {A 1 (ϕ n )} ⊂ U contains a convergent subsequence.
Due to ϕ n ∞ , c n ∞ ≤ 1, Eq. (34) yields ∂ 2 x c n ∞ ≤ αw 2 for all n, hence, using the boundary condition at x = 1, In view of these estimates, we conclude that the sequence {c n } is uniformly bounded and equicontinuous. Owing to Ascoli-Arzelà compactness criterion, we can therefore extract a subsequence {c n k } converging uniformly onĪ to some c ∈ A 1 (V ε ): whence the compactness of A 1 follows. (2) To show the compactness of A we rely on that A = A 2 • A 1 . Continuity of A is implied by Proposition 3. Moreover, if {ϕ n } is a sequence in V ε and we let c n = A 1 (ϕ n ), φ n = A 2 (c n ) = A(ϕ n ), the compactness of A 1 previously proved allows us to assume, passing to a subsequence if necessary, that {c n } converges in U . But then the continuity of A 2 implies that also {φ n } converges in V ε , and we have the thesis.
Theorem 5. Let Assumptions 1, 2, 3 hold with βw < 1, where β = β(ε) is the constant defined by Eq. (37). Then there exists a solution (φ, c) ∈ V ε × U to Problem (32) satisfying the a priori estimate Proof. We know from Proposition 4 that A : V ε → V ε is compact. We claim now that V ε is convex and closed in C(Ī).
The a priori estimate on φ, c readily follows from Lax-Milgram Theorem applied to Problem 1 with ϕ = φ. Theorem 6. Let φ ∈ V ε be any solution to Problem (32). Then Proof. ¿From Proposition 2 we know Rearranging the coefficients according to the definition of β 2 given by Eq. (36) we get whence the thesis follows using Morrey inequality at the left-hand side and then considering that Theorem 7. Let Assumption 4 hold and assume that a solution (φ, c) ∈ V ε × U exists to Problem (32), such that φ fulfills the free boundary condition (39). Then φ(x) ≥ φ 0 for all x ∈Ī.
(iii) For all c 0 ∈ (0, 1) there exists w * > 0 such that if w ≥ w * the equation in the unknown x has exactly one solutionx w ∈Ī.
Proof. Some of these properties are more easily proved by referring to Problem (42) rather than to its explicit solution (43).
(1) To show continuity of the mapping w → c so that, owing to Lemma 0, we have w1 (x) for all x ∈Ī as w 2 → w 1 , and continuity follows.
We observe that, under the hypothesis w ≥ w * , it results c w (x) < c 0 for x ∈ (x w , 1]. (5) Finally, let us consider statement (iv). First we check that the limit ofx w for w → +∞ exists by observing that for w 2 ≥ w 1 ≥ w * it results , where we have used the monotonicity of w → c (0) w and the definition ofx w1 ,x w2 respectively. Since x → c (0) w1 (x) is nonincreasing, this saysx w2 ≤x w1 , hence w →x w is monotone nonincreasing and admits therefore a limit at +∞. Also notice that Assume now by contradiction that lim w→+∞xw = ξ > 0. Fix then w ≥ w * and consider the point x = ξ ∈ I. Owing to (57) we have ξ ≤x w , whence c which is in contrast with (ii). Thusx w → 0 and we obtain (iv). (1) We claim that C is continuous on R + . Indeed, let w 1 , w 2 ≥ 0. Using Assumption 3 and Proposition 8 we get Due to the positivity of C(0), existence of a positive root of Eq. (46) is simply obtained by showing that C(w) < 0 for w large enough: Continuity will do the rest.
For c 0 ∈ (0, 1) fixed by Assumption 4, let w * > 0 be the value of w mentioned by Proposition 8(iii): We prove that there exists w * * > w * such that C(w) < 0 for all w ≥ w * * . To this end, we will henceforth assume w ≥ w * and rewrite C as (2) Notice that C 1 (w) > 0 for all w ≥ w * , indeed Γ(c (0) w (x)) > 0 for x ∈ [0,x w ). The Mean Value Theorem implies the existence of x w ∈ (0,x w ) such that w (x w )). By consequence, using Proposition 8(iv) we get C 1 (w) ≤ Γ Mxw → 0 for w → +∞.

Conclusions and perspectives
In this paper we have developed a macroscopic model which describes the evolution of a tumor cord using the theory of mixtures. The tumor mass has been modeled as a three-phasic porous medium composed by cells, extracellular fluid, and extracellular matrix, growing along a blood vessel whence cells get the necessary nutrient for their vital functions. As a matter of fact, since the main interest was in investigating the growth process in connection with the availability of nutrient, the extracellular matrix has been regarded as a rigid non-remodeling scaffold in which cells live and proliferate and extracellular fluid flows, therefore it has not played a primary role in the dynamics of the system. By means of suitable principles of mixture theory, in particular the reference to a Darcy-like framework to model interactions among different components, the problem has been reduced to a system of two partial differential equations. The first, for the cell volume ratio, is derived from the cell mass balance in which the velocity is explicitly expressed in terms of the internal stress state of the cells. It also includes a source/sink term describing proliferation or death of the cells according to the current size of the cell population itself, as well as to the available amount of oxygen present in the mixture. The second, accounting for the dynamics of the nutrient, is a linear parabolic equation modeling the diffusion of oxygen from the vessel wall toward the mixture and its contemporaneous uptake by the cells. This simple setting has to be understood in the light of the experimental evidence that the contribution of a possible transport by the extracellular fluid to the motility of oxygen molecules is definitely negligible with respect to their high diffusivity, hence the corresponding advection terms can be rightly dropped from the equations. Interaction of the growing tumor cord with a surrounding host tissue composed by normal non-proliferating cells has also been considered.
The mathematical formalization of these ideas results in a free boundary problem constituted by a system of two partial differential equations supplemented by suitable boundary and interface conditions. The free boundary is represented in this context by the interface separating the cord from the host tissue, which changes in time according to the growth of the tumor in interaction with the external normal cells and has to be determined as a further unknown of the problem.
The model is intrinsically multidimensional. In addition, it is not linked to a particular geometry because it does not assume any special configuration of the domain. In view of this, it is in principle applicable to a wide range of different scenarios. In this work, a two-dimensional problem has been specifically considered, focusing on the axial growth of the cord along the blood vessel and on its simultaneous expansion outward the vessel. Numerical simulations have shown ability of the model to reproduce specific features of vascular tumors, like the formation of an inner vital zone, constituted by sufficiently fed cells that duplicate, and an outer necrotic zone at the periphery of the cord, where conversely cells are not reached by a proper quantity of nutrient and cannot hence proliferate nor survive. It is worth pointing out that some different models of tumor growth assume a priori the existence of such zones and describe them invoking specific equations. In the present model, their formation is instead recovered a posteriori as a consequence of the overall dynamics predicted by a unique set of equations, by appealing to relatively general guidelines. Another interesting feature made evident by the numerical results concerns the evolution of the system toward a steady state. The cord clearly exhibits two different behaviors in the front part, the head, which is vital since always globally fed by a sufficient amount of oxygen, and in the rear part, the tail, which elongates as the head moves forward along the blood vessel, keeping a rectilinear shape with a nearly constant width in the transverse direction.
These considerations suggest that, as far as axial growth is concerned, two main questions have to be addressed, namely the dynamics of the head in terms of shape and velocity of propagation and the equilibrium of the tail in terms of maximum penetration inside the host tissue. In the present work, we have mathematically investigated this second issue, studying existence and regularity of physically significant solutions that describe the distribution of cells and nutrient across the tail at the steady state. By addressing furthermore the free boundary problem via a perturbative expansion of the solution, we have also been able to characterize qualitatively and quantitatively the steady growth width of the cord.
The model developed in this paper should be regarded as minimal, in the sense that it is deduced from few outline principles of a general well-coded theory, and takes into account only those really fundamental aspects to obtain a qualitative mathematical and biological description of the macroscopic phenomenon under consideration. Of course, it is liable to many improvements, which still at a basic level may involve, for instance, a more accurate description of both host cells and extracellular matrix, that here are essentially passive and lacking in their own physiology, and a more sophisticated coupling between the dynamics of the tumor cells and the nutrient, accounting for different uptake rates on the basis of the specific functions (survival, proliferation) carried out by cells. However, we believe that a minimal model, for which both a theoretical comprehension and validation are possible in view of its relatively simple structure, constitutes a solid foundation, hence a safe starting point, to tackle more complicated situations.