Lyapunov stability of competitive cells dynamics in tumor mechanobiology

Poromechanics plays a key role in modelling hard and soft tissue behaviours, by providing a thermodynamic framework in which chemo-mechanical mutual interactions among fluid and solid constituents can be consistently rooted, at different scale levels. In this context, how different biological species (including cells, extra-cellular components and chemical metabolites) interplay within complex environments is studied for characterizing the mechanobiology of tumor growth, governed by intratumoral residual stresses that initiate mechanotransductive processes deregulating normal tissue homeostasis and leading to tissue remodelling. Despite the coupling between tumor poroelasticity and interspecific competitive dynamics has recently highlighted how microscopic cells and environment interactions influence growth-associated stresses and tumor pathophysiology, the nonlinear interlacing among biochemical factors and mechanics somehow hindered the possibility of gaining qualitative insights into cells dynamics. Motivated by this, in the present work we recover the linear poroelasticity in order to benefit of a reduced complexity, so first deriving the well-known Lyapunov stability criterion from the thermodynamic dissipation principle and then analysing the stability of the mechanical competition among cells fighting for common space and resources during cancer growth and invasion. At the end, the linear poroelastic model enriched by interspecific dynamics is also exploited to show how growth anisotropy can alter the stress field in spherical tumor masses, by thus indirectly affecting cell mechano-sensing.


Introduction
Cancer disease occurs when irreversible pathological changes cause alterations of cells natural programs. Mutations of cell cycles induce abnormal proliferation generating internal forces that are continuously counteracted by the elastic resistance of the host tissue. As a result, solid tumor microenvironment becomes a dynamical scenario in which mechanical stresses, fluid pressure and elastic properties rapidly coevolve with the development of the mass, by promoting competitive interactions between the resident cell populations through a process known as mechanical cell competition [1]. Together with the ascertained stiffening of solid masses with respect to the host surroundings, the role Executive Editor: Xi-Qiao Feng B Massimiliano Fraldi fraldi@unina.it 1 Department of Structures for Engineering and Architecture, University of Napoli Federico II, Napoli, Italy that solid stress plays in many phenomena accompanying the initiation and the progression of cancer has been noticeably recognized [2][3][4]. The study of inhomogeneous residual stresses has in fact become a crucial aspect in oncophysics in order to understand some key mechanisms of tumor pathophysiology that could help both diagnostics and therapeutic interventions [5]. In fact, growth-induced stresses alter the homeostatic properties of the host tissue and cause different mechanosensitive cellular responses that may lead to dramatic degenerations, including central necrosis due to the abnormal intratumoral pressures compressing lymphatic vessels and the malignant progression of leading cancer cells driven by anomalous stress gradients at the host-tumor interface and peritumoral vascularization [6,7]. Also, solid stress acts as a mechanical barrier that contrasts drug inflow by compromising the effectiveness of cancer therapies [8][9][10]. For these reasons, the study of the effects of solid stress on soft tissue homeostasis has been investigated both theoretically and experimentally, by observing the in vitro growth of confined multicellular spheroids [11][12][13][14] as well as by means of magnetic or mechanical actuators to apply forces on in vivo systems [2,15]. The underlying mechanisms transducing a mechanical cue into a biochemical signal are explained on the basis of some conformational changes and molecular pathways that modify cancer single-cell properties, which exhibit a different cytoskeletal stiffness as well as altered adhesion and motility capabilities [11,[16][17][18][19][20]. Therefore, the understanding of the mechano-responses of cancer systems, from single-cells to entire masses, could actually lead to innovative targeting and therapeutic strategies that, by exploiting the mechanical differences between tumor and host as a selectivity principle, provide the possibility to mechanically attack cancer cells by preserving the healthy surroundings [21][22][23]. At the tissue level, the synergistic interplay between bio-chemical and physical signals significantly affect the biomechanical response of tumor masses, their study thus requiring enriched models that incorporate the most significant physiological events characterizing the mechanobiology of solid tumor growth [24][25][26][27]. As a consequence, the prediction of the influence that mechanical stresses can have on multiscale heterogeneous systems and the characterization of the dynamics to which tumor growth obeys have been investigated with increasing attention in oncophysics. In this framework, a great variety of modelling strategies have been proposed, which have from a side contributed to give some answers to open issues, on the other hand also highlighting the difficulties to include all the relevant aspects characterizing the complexity of the mechanobiology at the basis of the tumor development. In particular, how tumor invasion governed by cell-cell and cellenvironment dynamics and mechanical stresses do interact is not yet completely understood [28]. From a mechanical point of view, the macroscopic behaviour of growing tumors has been from time described by means of linearly elastic [29] and poroelastic approaches [7,25,30], as well as by means of hyperelastic models providing finite deformations [9,31,32]. In particular, poroelasticity is crucial to follow the development of solid stress and interstitial fluid pressure (IFP) within the tissue, this also allowing to trace the role of interstitial flows within the solid extracellular matrix (ECM) in cellular growth and migration [33][34][35]. However, most of these mechanical approaches include the effect of volumetric growth in the momentum equation by considering a priori prescribed growth functions [36] or by uploading phenomenological macroscopic tumor laws [9], in this way inevitably losing the key coupling with the underlying microscopic dynamics. In fact, while the growthassociated strains affect the stress, often the feedback of both the environmental stress and the related pressure gradients on the cells proliferation rate is not taken into account. More recently, coupled modelling approaches to characterize the elastic growth of prostate benign hyperplasia provided the influence of a mechanical feedback in the diffusion and the proliferation potential of the solid phase, by additionally including the interaction with chemicals [37,38]. The interactions among tissue constituents has been traced, at the macroscale, by means of multiphase models based on reaction-diffusion equations by accounting for the presence of microstress and macrostress fields [39][40][41][42][43]. At the microscale level, several works suggest the adoption of models analysing the mutual effects in terms of exchanged forces to simulate cell-cell interaction [43][44][45], while tensegritybased single-cell models have been recently used to also furnish a way for biomechanically discriminating among tumor and healthy cells on the basis of the redistribution of internal pre-stretch [17,46]. At the level of cell aggregates, statistical mechanics arguments have been used to demonstrate that cells interactions represent a functional constraint to the macroscopic tumor growth law [47] and the use of evolutionary models, in the framework of population ecology, has represented a particularly successful strategy to study the dynamics of cells' collective interactions through competitive/cooperative mechanisms, in which however the feedback from the environmental on cell behaviours seems to be still absent [48][49][50][51][52][53][54]. In the light of these considerations, the full coupling of population dynamics with the tissue biomechanical response and stress-driven factors in modelling tumor growth has been the crucial point developed in very recent works by some of the present authors [27,55]. Therein, the theoretical prediction of cell interspecific behaviour and the cells-extracellular matrix interaction have been described by means of a Volterra-Lotka (VL)-like model following a predator-prey logic, by assuming that cancer and healthy cells compete for the shared space and the available resources. As a matter of fact, this idea has also recently found confirmation, since the mechanical cell competition for the common space has been demonstrated to be a process that actually triggers the preferential elimination of one cell population with respect to others through apoptosis, with the generating competitive interactions regulated by short-range biochemical inductions, mechanical constraints and compressive forces [1,[56][57][58]. Despite the adoption of a geometrically and constitutively nonlinear framework allows to follow the evolution of the tumor mass by accounting for finite motion and large deformations [55,59], some crucial information concerning the quality of the dynamical system are somehow blurred by the high degree of mathematical coupling between mechanical and species variables in the system of nonlinear partial differential equations governing the problem. Therefore, to better elucidate how solid species behave and to gain further insights about the interplay between mechanical and biological factors, we here recover the linear elasticity in order to benefit of a reduced mathematical complexity for analysing the stability of the mechanobiological systems. This is accomplished by showing that the well-known Lyapunov stability criterion can be derived from the thermodynamic dissipation principle that defines the constitutive framework of the poroelastic model enriched with competitive cells dynamics. By following this way, we then study the asymptotic behaviour of the species also with the help of an uncoupling strategy based on the difference between the elastic and biological characteristic times involved in the continuity equations. Finally, under the assumption of linear poroelaticity, some results of numerical simulations that involve possible anisotropic growth, which depends on tumor cell density, intratumoral stress distribution and on how neighbouring cells interact each other during proliferation by straining the matrix, are shown to point out the consequences of non-stereologically growing masses in terms of change in material properties (remodelling) and tumor stiffening.

Thermodynamic framework and stability of motion in growing systems
The theoretical model used to predict tumor growth combines the classical poroelastic field equations with inelastic growth through the direct (full) coupling of the poromechanical balances with a suitable dynamical system obeying a Volterra-Lotka logic. In this way, the inelastic deformation due to volumetric growth can be seen as the net result of the interaction among the biological solid species that inhabit the tissue. In particular, we start from some well-established thermodynamic considerations in order to derive via dissipation principle both the constitutive assumptions of the linear poroelasticity and the configurational thermodynamic forces conjugated with growth deformation. In addition, in the light of the coupling between the mechanical problem and species dynamics, the dissipation principle also allows to find a condition that ensures both the stability associated with the motion of the solid constituents and the thermodynamic consistency of the model. To this purpose, the growing continuum can be seen as a multi-component system, in which the solid constituents are denoted by a subscript i ∈ S , while the fluid phase is indicated by the subscript i = F. However, the high water content in the involved biological species lets to assume a constant true density (say ) for each constituent, so that the total density of the system remains unchanged: the sum of all the volumetric fractions giving the unity (in what follows, the subscript i will exclusively refer to the solid constituents and the hypothesis of constant density will be hold true). Growth is an inelastic process that inevitably implies dissipation. For this reason, the free energy density ψ is here expressed as the sum of a reversible energy aliquot, which is linked to the purely poroelastic contribution and classically depends on the elastic strain, say E e , and the fluid content φ F [60], and a growth-associated contribution that is instead a function of the solid species φ i . In explicit: Under the hypothesis of an isothermal process, the balance of energy can be written by introducing a specific metabolic energy contribution, say g , which is directly responsible of the volumetric growth g. Hence, the first principle of thermodynamics reads as: where U and V are the internal energy per unit mass and the volume measure, respectively, while D =Ė = sym(u ⊗ ∇) is the symmetrical part of the velocity gradient. In fact, in the linear regime, the total deformation E and the displacement field u are related each other by the compatibility equation E = sym(u ⊗ ∇), which can be decomposed in the sum of the elastic strain and the growth strain, i.e. E = E e + γ g (see Refs. [61,62]), where γ is a tensor of anisotropy that distributes volume growth in the different directions through specific weight coefficients. By following well-known approaches of the Thermodynamics of multicomponent bodies, the entropy imbalance must account for the energy contribution carried out by the species transport and internal supply [63]. In the light of these considerations, the second principle is written down by separately including a thermodynamic force f g conjugated to the rateġ (representing a volumetric rate of entropy supply due to growth of the solid species) and a thermodynamic tensor force f γ that is instead associated to the internal remodelling variables that describe the structural remodelling caused by deformation and growth [64][65][66]. Furthermore, by adding up a conductive rate of entropy supply due to the fluid transport [63], the entropy inequality reads as: where q F is the fluid flux vector and μ F = ρ −1 F ( p − p 0 ) is the associated fluid chemical potential that allows to describe energy exchanges due to fluid transport in the control volume. By subtracting Eq. (4) from Eq. (3), and identifying the free energy per unit volume defined in Eq. (2) as ψ = ρ(U − S), the following dissipation inequality can be obtained: To express the former equation in a local form, the Gauss theorem is applied to the second member of the left side of Eq. (5), while the chain rule on the right side allows to express the rate of ψ according to Eq. (2). In this way, it results: By both invoking the additive decomposition of the strain and taking into account the fluid continuity equation: the terms in Eq. (6) can be collected to have: from which one readily gets the standard constitutive assumptions for the poroelastic variables Equation (8) can be further reduced by taking a Darcy type flux vector q F = −K∇μ F (with K being a definite positive matrix) and by properly introducing a fluid source term driven by the static pressure gradient Γ F ∝ p 0 − p. In order to make explicit constitutive assumptions and obtain a suitable expression for the linear elastic potential as a function of the poroelastic fields (the strain tensor and the interstitial pressure, according to most of the classical models [27,60]), a Legendre transform is performed on ψ e in a way to have a dual energy density: By then considering its differential form and accounting Eq. (9b) one has the identity: (11) from which Equation (12b) together with Eq. (10) let to obtain an explicit constitutive relation for the fluid content, which reads as linearizing the incremental fluid fraction dφ F , thus writing in which A = −∂ 2 ψ * e /(∂E e ∂ p) and M −1 = −∂ 2 ψ * e /∂ p 2 , φ F0 being the initial fluid fraction at zero elastic strain and at the reference pressure p 0 . Hence, the substitution of Eq. (13) into Eq. (12b) and a direct integration with respect to p gives: where the energy density termψ resulting from the integration can be interpreted as the effective energy aliquot associated to the solid deformation. Within a linear isotropic framework, this coincides with the quadratic elastic St. Venant-Kirchhoff strain energy density, i.e.: with C = ∂ 2 ψ * e /∂E 2 e = 2μI + λI ⊗ I being the isotropic fourth order stiffness tensor and thus, together with the relation (9a)a, the Terzaghi stress principle can be immediately derived: (16) In the light of the specific constitutive assumptions here introduced, the dissipation inequality (8) reduces to: These terms are evidently not independent, since the coupling suggests a relation between the volumetric growth strain g and the vector of the biological species. However, they carry different information that allow to discuss them separately. The first two terms are in fact associated to kinematical effects of growth and let to assume an expression for the thermodynamical forces conjugated to growth and remodelling. In particular, the growth volumetric strain is driven by both mechanical and metabolic factors and the configurational force conjugated to the growth can be expressed as f g = σ : γ + g . This relation (that also resembles the Eshelbian stress-energy tensor classically derived in nonlinear formulations [27,67]) highlights the role of a chemical sources in triggering volume growth, mediated by the weighted sum of the stresses that can invite or inhibit some growth kinetics in the different directions. In case of isotropic growth γ = 1/3 I, so that the growth results to be driven by the hydrostatic pressure σ hyd . Analogously, also stress-driven anisotropy coefficients could be in general assumed by considering the Kirchhoff-type stress f γ = gσ as conjugated pair [65,68,69]. These relations qualitatively suggest that growth takes place in presence of available metabolic resources, in a way to be sustained also in a stressed (e.g. compressed) environment, where mechanical impairments could occur. In fact, the respect of Eq. (17) would imply that the lack of metabolites and unfavourable compression would either arrest the growth or cause resorption processes, by orienting the growth towards less stressed compartments. Also, these driving forces allow the introduction of suitable feedback mechanisms thanks to which it is possible to set the coupling between the mechanical variables and the fitness functions of biological species. In particular, one can assume a direct influence of the mechanical stress on the growth, by introducing non-constant VL coefficients able to describe how cells intrinsic rates feel the environmental stress through specific mechano-sensing parameters. The reduced dissipation (17) finally gives the condition ∇ φ ψ g ·φ ≤ 0. This dissipation inequality can be connected to the stability of the species dynamics. In fact, the emerged inequality represents the well-known Lyapunov stability criterion [70], which relates the stability of the biological dynamics to the existence of a positive-definite Lyapunov function in the positive octant of the phase field that vanishes at the equilibrium points and such that its time derivative is locally negative semi-definite in a neighbourhood of the candidate attractor points. Inequality (17) then actually states a substantial equivalence between the thermodynamical and stability of the multi-species dynamics, in the light of the fact that the entropy principle can be attributed to the study of the stability of the Lyapunov function of the generalized VL system at hand, a key parallel for the present application that has been also deeply investigated (without considering mechanical coupling) in the field of ecological models in Ref. [71]. By virtue of the Lyapunov stability theorem assessments, one can suppose that, in the neighbourhood of an equilibrium point φ * , the dissipation potential can assume the following expression: in which the general form of the evolution equations of the solid constituentsφ i = φ i Γ i has been involved, the function Γ i representing the fitness function of the i-th species, whose expression will be specified in what follows for both the interacting cells and the ECM with reference to the particular problem. Also, Eq. (18) recalls the Boltzmann entropy formula. In order to define a specific dissipation function for the species, the study of the stationary points is then required.

Tumor growth resulting from poroelasticity enriched by predator-prey cell dynamics
Cancer disease involves a cascade of pathological genomic alterations that compromise, often irreversibly, ordinary cellular programs and functionalities. These mutations are accompanied by the progressive modification of the extracellular environmental conditions, defined by immune response, matrix metabolism as well as stiffness, mechanical and biochemical gradients. The growth of solid tumors can be treated physically as a mechanical process in which a heterogeneous mass expands within (or against) a surrounding host tissue. Tumor expansion is in fact controlled by internal driving stresses, which are counterbalanced by the mechanical (elastic) resistance provided by the surrounding environment. Internal stresses are mostly originated by the abnormal proliferation of cancer cells, whose viability and distribution depend on the intratumoral availability of nutrients. This implies that the physical forces pushing the tumor ahead do not involve the sole surface tension and the pressure of the surrounding medium, but also explicit active cellular forces in the momentum balance that, in turn, activate mechanosensing cellular responses. With the aim to model the complex machine of the host-tumor interaction in growing solid tumors, we here adopt an enriched poroelasticity approach in which the mechanically activated stress fields, fluid pressure and nutrient walkway are all coupled with spatially inhomogeneous and time-varying bulk growth. This growth is induced by competitive dynamics occurring at the microscopic scale level among healthy cells, cancer cells and ECM (projected at the macroscale through their corresponding volume fractions) and modelled by introducing ad hoc non-linear Lotka/Volterra-like equations. The basic idea is that cancer and healthy cell species do not compete directly, as it would happen in a pure predator-prey logic, but instead contend the common resources and the available space. Cells competition leads to the preferential elimination of a cell population by the other, mediated by specific environmental conditions and interspecific communication pathways [1,47]. The introduction of these effects enriches previously proposed elastic and poroelastic models [28,30] by explicitly introducing the actual competition among cell species and by reproducing the experimentally observed coupled dynamics in which the presence of one species tends to somehow limit the development of the other one. This aspect is here described in explicit as a behavioural phenomenon occurring in the cells community without prescribing macroscopic tumor growth laws. In fact, in the proposed model, mutual interactions in turn modify the intrinsic growth rates of the cell populations and leads to spatially inhomogeneous elastic and residual stresses as well as non-uniform IFP distributions. In the present treatment, however, in order to elucidate in the simplest way the key aspects of the dynamics at hand, further elements that would imply a direct competition between cancer and healthy cells depending on additional factors, such as the anti-oncogenic potentials of some immune cells or the aggressiveness of pre-malignant cells which become malignant as a result of mutation processes, will be voluntarily neglected.

Governing equations of linear poroelasticity incorporating growth
By recalling the constitutive equations of linear poroelasticity [60], the strain tensor E, defined on a closed subset Ω ⊂ R 3 in presence of volumetric growth, can be written as where γ g is the growth strain tensor previously introduced, with g ∈ R being the volumetric growth strain function and the matrix γ = diag {γ k } containing the anisotropic multipliers for each principal direction so that tr (γ ) = k γ k = 1. The matrix S represents the drained compliance fourth rank elasticity tensor (the superscript (d) will be avoided) while σ e f f is the effective stress tensor resulting from Eq. (16). Its definition includes the Cauchy stress tensor σ and the interstitial fluid pressure p (IFP), while A is commonly referred as the Biot effective stress coefficient symmetric tensor, equal to A = I − C · S (m) : I, in which C = S −1 is the drained stiffness tensor and S (m) is the matrix-associated compliance fourth-order tensor, while I and I respectively indicate the fourth-rank identity tensor and the second-rank identity vector (in Voigt notation). By deriving the stress tensor σ from Eq. (19) and in absence of body forces and neglecting inertia terms, the stress equations of motions in three dimensions read as in which E e = (E − γ g) is the elastic part of the deformation. By combining Eq. (20b) with compatibility Eq. (19b), together with the hypothesis of considering an elastic isotropic material (also implying that A = AI), the quasistatic balance of linear momentum can be written as: with are the Lamé moduli. The other basic field variable of poroelasticity is the dimensionless (under the hypothesis of constant density) variation in fluid content ζ that, according to the constitutive Eq. (13), can be linearly related to the elastic strain and to the pore pressure [60]: This relation allows to rewrite the fluid mass conservation equation as: where q F = −(1/υ F )K∇ p, with υ F denoting the fluid viscosity and K the intrinsic permeability symmetric tensor, whereas Γ F is a source/sink term (fluid mass supply per unit volume) introduced to account for the fluid interchange from the leaky capillaries to the absorbing lymphatic vessels within the interstitial space at the microcirculation level, modelled according to the Starling's theory. In particular, fluid movement within the microvascular beds is governed by the hydrostatic gradient between the IFP and the capillary pressure ( p v − p) [10], contrasted by the difference between capillary and interstitial osmotic pressures (π v − π ι ) weighted by a reflection coefficient . In addition, the effect of lymphatic drainage in the opposite direction is included, driven by the difference between the IFP and lymphatic vessels pressure p l [72][73][74]. Pressure gradients in Γ F are mediated by two conductivity coefficients, named κ v and κ l respectively. Other characteristic poroelastic constants appearing in Eq. (22a) and used to define the inverse of the Biot modulus M −1 are the effective compressibility factor C e f f , in which the porosity ϕ and the fluid bulk modulus K (F) are also included (the fluid is supposed to be incompressible, so that 1/K (F) → 0), and the Skempton compliance difference tensor B, linked by the following relationships: In Eq. (24a), the Skempton tensor reduces to B = (B/3) I because of the isotropy assumption, with B = 1/K − 1/K (m) / C e f f , tending to the unity for fluidsaturated materials, as well as for K (m) approaching to K (F) , whether or not the fluid and the matrix are assumed incompressible, in this case being the porosity ϕ actually uneffective. By further introducing the undrained elastic constants as and by also invoking isotropy, the scalar Biot effective stress coefficient A obtained in Eq. (21) can be expressed in terms of undrained and drained constants and the Skempton coefficient B Then, by focusing on a spherically symmetric case to have the deformations Eq. (19) be written as E = diag {∂u/∂r , u/r , u/r } and the multiplier γ = diag {γ r , (1 − γ r ) /2, (1 − γ r ) /2}, the stress-strain-pore pressure constitutive Eq. (20a) take the form: in agreement with literature formulations regarding the characterization of porous materials and the anisotropic expansion of linear elasticity by means of prescribed growth laws [28,75]. By starting from the definition of the isotropic Skempton coefficient B in Eq. (24b), analogous passages let to find C e f f = A/(K B) = 9 ν (u) − ν / B 2 E 1 + ν (u) , so that the constitutive relation for fluid content ζ in Eq. (22a) combined with the expression of the Biot modulus given by Eq. (24a) becomes with = trE, the latter relation showing how the growth g counteracts the dilation of pore spaces that could be occupied by the growing solid phase, which is also directly responsible of the nutrient consumption carried by the fluid. The adopted framework then allows to express the governing equations of the poromechanical problem (26) and (28) by involving a minimal set of measurable constitutive parameters, i.e. the drained moduli E and ν, the undrained Poisson ratio ν (u) and the Skempton coefficient B.

Mechanical and interspecific coupling and specialization of the competitive interactions
By recalling that the growth strain field appearing in the poroelastic constitutive Eqs. (27) and (28) is an explicit function of the solid species, i.e. g = g (φ), where the vector φ collects the volumetric fractions of the constituents involved in the dynamics, in this case represented by the healthy and tumor cells and the ECM, say φ H , φ T and φ M , the poroelastic field Eqs. (20b) and (23) are finally coupled with the dynamical system leading to write the mass balances in the form Herein, Γ g = {Γ i }, i = (T , H , M) represents the vector of the net rates of the solid species responsible of the growth term g, whose expression explains the particular way in which the species interact each other within the environment. To this aim, by focusing on the competition between cancer and healthy cells interacting with the ECM in a shared environment by challenging for the common resources, VL equations have been here employed to describe the dynamics among the species. Evolutionary game theory and population ecology are widely used in the literature to model some cell-cell and cell-immune system interactions [51,53,54,76]. Here, the proposed strategy enriches these dynamics with the possibility of providing the coevolution of the environmental changes through non-constant parameters and poromechanical variables. In agreement with the previous nonlinear formulations of the model [27,55], whose outcomes are also confirmed by recent works demonstrating that tumor heterogeneity and cell-cell interactions within a shared space constrain cellular proliferation by affecting the macroscopic tumor progression [1,47], mutual inhibiting inter-species parameters are introduced to trace back the competition that each species exerts on the other one by affecting the growth rate. In particular, with reference to the particular case of spherically symmetrical solid tumors, system (29) can be specialized by providing functional couplings in the system of equations: Cell-specific Eqs. (30c) and (30d) recall the classical VL-like ones. However, differently from the classical VL equations in which coefficients are all constant, in order to realistically describe the competitive dynamics among cell species and ECM in tumor spheroids, two main aspects here suggest to adopt varying coefficients in the VL equations. In fact, since cells activity takes place in presence of nutrients, supposed to be homogeneously dissolved within the fluid, the proliferation rates are here proportional to the variable fluid content φ F = φ F0 + ζ . Additionally, the role that mechanical stress has in controlling cell proliferation is also taken into account. Several literature studies in fact investigate with increasing interest the importance of anomalous mechanical compression in controlling tumor progression via apoptosis, due to the increased cell packing and reduced proliferation, as well as by regulating cancer cells migration towards less compressed districts through the activation of specific biochemical pathways [1][2][3][4]11,13,14]. In order to include mechanotransduction processes in the system, specific η T and η H are modelled as inhomogeneous parameters that result to be a function of the hydrostatic stress perceived by the cells σ hyd = trσ /3. Hence, given that hydrostatic compression inhibits the intrinsic cell rates η T and η H in cells experience quiescence, and by assuming that both cancer and healthy cells are virtually indistinguishable in such a state, the cell proliferation rates are supposed to reduce to a common lower growth rate η q . In agreement with recent literature works proposing a role for a certain homeostatic pressure in controlling competitive cell interactions [1,58], the mechano-sensing behaviour is modelled here by adopting the transition function displayed in Fig. 1, whose expression reads as follows: where K = {T , H }, χ σ is a constant and the coefficients η T 0 and η H 0 represent the nominal replication rates, depending upon the doubling times T T and T H of tumor and healthy cell species, respectively, lying in the range 17-40 hours for the most of cell types, as also experimentally determined [55,77]. The quiescent metabolic rate η q is a reduced non zero metabolic coefficient accounting for the fact that quiescent cells preserve a basal metabolic activity until favourable environmental conditions are restored [78]. To account for competition, the growth rates η T and η H are then penalized through the VL-coefficients α I J , with {I , J } = {T , H , M}, the case I = J denoting self-inhibiting terms while the coefficients with I = J providing the constrains due to species mutual interactions. Accordingly, the mutual competition terms α H T and α T H weight the way in which a cell species is influenced by the presence of the other one, while the coefficients α H M and α T M measure the occupation degree of the ECM in the space contended by cells. The coefficients α T T and α H H are self-competition terms accounting for the carrying capacity of each cellular species. These coefficients can be considered constant and not affected by the mechanical stress, this because quiescence does not likely influence the aggressive mutagenic phenotype of cancer cells against normal cells. All these considerations let to describe a sufficiently complex population dynamics with respect to the standard VL models, being the game of interaction affected by the environmental parameters as well as by the possibility of changing the game strategy through the modulation of the intrinsic rates in order to find the fittest way to survive. The third interspecific Eq. (30e) represents the balance of ECM in which biochemical differences between healthy and tumor ECM components are neglected and therefore the resulting overall ECM fraction, φ M , dynamically depends on synthesis and degradation processes promoted by cells, through the production characteristic rates β T and β H , and the loss rate δ M , weighted by the coefficients α MT and α M H . In addition, to also take into account that, within a tumor, abnormal IFP rise can occur because ordinary lymphatic drainage mechanisms are altered due to enzymes-mediated degradation and stress-induced collapse [74], the conductivity parameter κ l is assumed to decrease with the tumor fraction according to the law: with κ ln being an assigned constant denoting the nominal value of drainage conductivity in normal tissues. Finally, the net growth strain term (assuming a common density) can be defined in terms of the solid constituents as: the additional subscript 0 indicating the volumetric fractions of tumor and healthy cells and matrix, respectively, at the conventional initial time t = 0. The values of the parameters used in the simulations are all summarized in Table 1. Before presenting case studies, it is worth to discuss some important considerations regarding the stability of the dynamics describing cancer cell invasion, by studying the related equilibrium points and by applying the stability principle to the coupled system (30).

Stability of cancer invasion
The linear formulation of the model allows for noticeable simplifications with respect to the nonlinear dynamics proposed in previous works by some of the present authors [27,55], because the degree of coupling between mechanical and species variables is here reduced. This is particularly advantageous for discussing the equilibria and stability of the system (30). In order to study the stability of the dynamics in a slender way, it is here proposed to separate the first two poroelastic nonhomogeneous PDEs from the multispecies equations on the basis of their characteristic times. In fact, the first two poroelastic Eqs. (30a), (30b) are intrinsically gov-erned by the speed of the elastic medium (that is proportional to √ K / ) and by the rise/decay times of the interstitial pressure τ + p ≈ (κ v M) −1 and τ − p ≈ (κ l M) −1 , all of the order of the seconds or fraction of seconds. On the other hand, cells dynamics processes will adapt nearby their attractors presumably in a slower manner, because of both the greater biological intrinsic proliferation times (of the order of several hours or days) and the VL mutual interactions. On these bases, by condensing the poroelastic field variables within the vector S = {u, p}, the dynamic-coupled differential equations (30) can be rewritten in the following way that highlights the interplay between fast variables S and slow variables φ, since, regardless of how time can be scaled, the relative rates will verify the condition Ṡ ∞ φ ∞ , while L(·) and Π(·) collect the linear operators of equations (30a), (30b). This implies that, for rescaled times t comparable with the characteristic times of S, one can obtain a fast subsystem (FS) corresponding to where the species variables responsible for the macroscopic growth do not vary, while, for rescaled times t compatible with the (slowest) characteristic time of the cells, a slow subsystem (SS) can be instead written down as The dynamics of the original system can be explained in terms of the interlaced evolution of the respective fast and slow subsystems. This simplifying assumption also supports the previously made hypothesis according to which uncoupled poroelastic and species potentials were considered in Eq. (2), the functions ψ e and ψ g depending on the poroelastic variables and the species vector, respectively, and their respective rates being de facto separable since they can be performed by mutually fixing the growth level or the stress conditions. By then considering a constant φ in the FS, the problem (35) reduces to a classic poroelastic problem with a prescribed growth-associated inelastic strain, whose stability can be ascribed to the well-established ellipticity of the governing equations of standard poroelasticity. At the cell characteristic time scales, Eq. (36) let to study the stability of the species for prescribed microenvironmental conditions. In fact, under this hypothesis, when species oscillate in the neighbourhood of their stable attractors, stress and fluid do not significantly vary in a way that the stability condition (17) of the growth trajectories can be satisfied by holding the state variables E e and φ F at a certain constant value, say {E e , φ F } (or, equivalently, {E e , p}). In the light of the theorem proposed by Tuljapurkar (see Ref. [79]), this condition is equivalent to study the local stability of the species motion around a given equilibrium point φ * in the phase space delimited by the tetrahedron , p)}. In addition, it is noticed that the localization of the entropy inequality requires that the stability of the species trajectories can be expressed (at each point P of the open domain Ω of the problem) by means of the ordinary condition: Therefore, by focusing on the radially symmetric problem set by Eq. (30), let us introduce a dimensionless variable defined as z = r /b + η T 0 t to so obtain the following ordinary differential equation (ODE) subsystem of three-species LV equations: in which tilde denotes the initial LV coefficients divided by η T 0 , while the function q > 0 takes into account the influence of the poromechanical coupling term, in which both fluid pressure and elastic strain are fixed at the slow time scale of species oscillation. The positivity of the coupling function q derives by observing that both the mechanosensing transition function (31) and the fluid phase assume non negative values, the complete saturation of the pore space being avoidedin the light of the linear constitutive assumption providing small porosity variations. The three-species system (39) let to determine in a more direct way both the physically and biologically consistent equilibrium points and the trajectories lying in the domain I. This is also guaranteed by considering positive initial conditions for the species, which ensures that the positive octant is invariant and that all the trajectories of motion keep bounded in this (0, R + )×(0, R + )×(0, R + ), as shown by Itik and Banks [53] through the help of two lemmas that are here reported for completeness: Lemma 1 With all positive initial conditions, the solutions of the system (39) lie in (0, R + ) × (0, R + ) × (0, R + ) . This follows by observing that each coordinate hyper-plane is invariant.

Lemma 2 The solutions of Eq. (39) with initial values in
Then, starting from Eq. (39), the Jacobi matrix J Γ = ∇ φ Γ of the linearized system results to be (tildes are hereinafter omitted for the sake of simplicity): Steady-state points are found by setting Γ = 0. The physically admissible points which can potentially occur describe the following exclusive situations of interest.

Cells fractions extinction. This first situation is depicted
by three equilibrium points, which can be written as The value φ * M1 can be determined in the two first cases respectively as 1/α H M and 1/α T M , while it is undetermined in the last case. Then, the first points can be included in the last one so that a sole equilibrium φ * 1 ⊂ I can be analyzed. Basically, in the particular case in which cells approach the extinction, this state describes a situation in which the solid matrix remains unchanged in absence of cell activity. This is also due to the fact that, for example, no self-degradation had been provided in the model through an intrinsic decay rate. 2. Healthy fraction dominance. This situation is associated to the steady state φ * 2 with expression (41) and corresponds to a system in its healthy stage, without invading tumor species.
3. Tumor fraction dominance. The equilibrium point (say φ * 3 ) is given by (42) and it is associated to a completely tumor stage, in which space is fully invaded by cancer cell. 4. Solid species coexistence. These (two) points have a longer expression in terms of model coefficients, and are here simply denoted as φ * } for the sake of simplicity.
With reference to parameters adopted (see table Table 1), a situation of strong dominance occurs for the tumor species. The introduction of the aggressive tumor species with a fitness function that, in relation to the modified LV coefficients, prevails on the healthy counterpart, turns the ECM-healthy cell coexistence (i.e. in absence of tumor species) into an unstable equilibrium. In particular, the insurgence of a predator species within the ecosystem possessing the given characteristics heavily affects a condition of normal homeostasis. Provided that points in 1) and 4) have no physical interest for the consideration at hand (in particular, equilibria 4) lead to not acceptable solutions for the adopted parameters), this corresponds to observe that the tumor presence generates an unstable direction in correspondence of point 2), for which the perturbation of the condition φ T = 0 brings the system in a new equilibrium state characterized by cancer dominance. The eigenvalues associated to points 2) and 3) (say y 2 and y 3 ) are: The projections in the phase space of Fig. 2a,b respectively show equilibria φ * 2 and φ * 3 in their planes, i.e. φ T = 0 and φ H = 0. In particular, focusing on Fig. 2a, the projections of phases trajectories in the plane φ T = 0 exhibit a convergent trend towards the point φ * 2 . This dynamics thus connotes normal homeostasis in the healthy districts, the plane φ T = 0 constituting a stable subspace for state 2), described by the two dimensional stable subspace E 2 = span{v 21 , v 22 }, where v 21 and v 22 are the eigenvectors corresponding to the eigenvalues y 2 with negative real part (illustrated by the blue arrows of Fig. 2a). Out-of-plane perturbations, i.e. the appearance of tumor cells, will intercept new trajectories that lead the system to a new malignant equilibrium. In particular, the third eigenvalue y 23 has positive real part for the adopted LV coefficients and also given that q is strictly positive, and thus generates an unstable out-of-plane direction v 23 according to the first Lyapunov stability theorem (represented in Fig. 2a  by means of a red arrow). Therefore, the activation of the tumor fitness function, as also clearly highlighted in Fig. 2c in which the cell phase space for φ M = φ * M2 is represented, strongly diverts the trajectories of the biological species from the healthy subspace (the φ H -axis, in this case) by pointing towards the tumor stage. This means that the healthy condition cannot be maintained within the body points at incipient cancer invasion. On the other hand, in the tumor regions of the domain, the equilibrium φ * 3 results a stable attractor with {y 3 } < 0. As shown in Fig. 2b in the plane φ H = 0, trajectories converge to a stable tumor node, out-of-plane perturbations not producing escaping trajectories. This can be also observed in the cells phase space of Fig. 2d , where trajectories towards the healthy axis are absent and the tumor attractor keeps stable. The influence of intratumoral compression on species intrinsic dynamics is reported in Fig. 3. In particular, by considering a fixed amount of fluid, that, in a linear regime, does not exhibit significant variations, the evolution of species trajectories for the healthy and cancer dominances has been analysed for three selected values of the mechanosensing curve S K , representative of high, intermediate and low stress conditions, respectively. From results in Fig. 3, it clearly emerges that the velocity of cell species is strongly reduced by higher hydrostatic pressures, the flow lines appearing practically orthogonal to the cells axis in this situation, while trajectories progressively point to their respective nodes with increasing speed when environmental stress is relieved.
Naturally, all these considerations are related to the chosen set of coefficients. The variation of the coefficients α I J , δ I and β I may generate different scenarios, in which altered dynamics may occur exhibiting different stable attractors (like the ones related to the cells coexistence that additionally include oscillating patterns that are not obtained here, or the one related to the tumor extinction in which the healthy species is able to survive to a given ad hoc treatment) also in presence of a tumor-type fitness function in the species rates. All the coefficients in fact reflect some cell-cell and cell-ECM interactions occurring at the microscopic scale level and may be actually interpreted as the net result of multiple reactions-based processes. There could be therefore the possibility of directly modelling these interactions by means of a more complex multiscale approach in which chemical balances are included, these ones working as intermediate regulators of the cross-talks among populations at the tissue scale. In this way, the eventual competition between tumor and immune system could be also taken into account, by so envisaging a perturbation of these mechanisms, in particular by explicitly adding antagonist species like drugs in order to analyze the variations of the tumor strategy coefficients and capture possible changes into the system fate through the analysis of its asymptotic behaviour.

Effect of growth anisotropy in spherical tumor masses
By focusing again on a spherically symmetric case, the model can be studied by analyzing the evolution of a growing heterogeneous sphere with radius b in which internal distributions of the solid species are assigned at the initial time, by finding numerical solutions to the system (30) in the five unknowns {u, p, φ T , φ H , φ M }, being the sole radial non zero displacement component and all the functions depending upon r and time t. An internal inclusion with radius r i → 0 is considered in order to avoid singularities, by there prescribing a vanishing displacement. At the external radius r = b, traction-free conditions are considered, i.e. σ r (b, t) = 0. In addition, no Darcy fluxes occur at r = b and r = r i . With reference to the initial conditions, solid species distributions, pressure and displacement have been prescribed at the starting time t i = 0 + as follows: Herein D i (r ) = {1 + exp [χ r (r − a) /b]} −1 is the initial distribution function of the tumor species, where χ r is a smoothing coefficient and a is the initial tumor front, while the constants φ T 0 and φ H 0 are the amplitudes of the initial tumor and healthy phases. These conditions imply that no residual stress is provided in the initial configuration. Numerical simulations, performed by means of the NDSolve package in Wolfram Mathematica ® [80], were carried out by varying the value of the growth anisotropy coefficient γ r . This anisotropy multiplier, which also affects the growth in the circumferential directions because of symmetry, is a coupling parameter taking into account the fact that growth along a given direction occurs by exerting mechanical strain and somehow translates the cellular patterns assumed by the interacting cells during the growth process. By here assuming γ r to be constant, its value is taken as variable between 1/3, which characterizes a case of isotropic growth in which the strain induced by the mass change equally deforms the solid matrix in the principal directions, and the limit case of γ r = 0. The latter value would instead correspond to a situation in which cells grow by mainly straining in the circumferential direction, the neighbouring cells mechanically interacting each other within each spherical layer and the tumor accretion occurring through the progressive formation of so called cell cords that radially invade the space, as exhibited by specific cellular patterns in some multicellular tumor spheroids (see Ref. [81]). These two limit situations clearly produce different stress distributions within the tumor microenvironment. In particular, by looking at Fig. 4, the anisotropization of the growth induces higher intratumoral compressions, whose magnitudes result to be up to two-fold Radial (upper row) and circumferential (lower row) stresses at different times and for varying anisotropy coefficient γ r of those ones obtained from the isotropic case classically assumed. This difference can be explained observing that layered growth is characterized by a stronger circumferential mechanical interaction between cells belonging to the same stratum, which in turn generates a tightening effect that elastically wraps the inner districts of the spheroid by so intensifying the level of intratumoral compression. It is worth noting that, although cell mechanosensing is here assumed to exclusively depend on the hydrostatic part of the stress, possible growth anisotropies indirectly affect cells proliferation through the level of induced residual stresses, thus representing a further indirect mechanical feedback in the light of the established fully coupled framework. Also, an increased hydrostatic compression would also affect the elastic properties of the tumor. In fact, the residual stress imprisoned in the tumor spheroid is known to promote the material remodelling of the host tissue stiffness [59]. In particular, the hydrostatic pressure influences the local bulk modulus, with a consequent effect on the evaluation of the effective bulk response of the mass as well. In the present case, the tightening associated to the anisotropy patterns would deepen the stiffening of the tumor, which is one of the most apparent features characterizing solid tumor growth in soft tissues. Such an effect would be crucial to be investigated for the implications that the accurate prediction of intratumoral solid stresses and mechanical stiffness have on both anti-cancer interventions and diagnostic procedures. Actually, the variation in mechanical properties is in fact a determinant marker in hardness-based cancer detection and, although durometry is still a qualitative and not completely engineered method, the clinical palpation remains one the most adopted practices to individuate the presence and to assess the grade of non metastatic masses and tumor nodules [82,83]. The magnitude of the stresses also affect the species development through the mechanical feedback previously described. As shown in Fig. 5, the tumor species fast grow within in the tumor region and slightly starts to invade the surrounding host tissue. Also, by comparing the two limit situations, tumor volumetric fraction results to be more inhibited by the pressure in the anisotropic case, while no relevant changes occur in the outer healthy surroundings. In the presented simulation, the growth of the tumor amount results greater than its migration potential, meaning that the tumor cells (for the assigned initial state) tend to "bulky" proliferate in the internal zone by penalizing the healthy constituents and by saturating the available space before diffusing. Experimental observations report that tumor migration is strongly enhanced by the level of mechanical compression, which induces border cancer cells to augment their motility [84]. In the present case, an initially unstressed environment was assigned, this implying that the stress within the spheroid can be considered still sustainable by the internal proliferating cells.
The abnormal proliferation of tumor cells also affects the level of interstitial fluid pressure. As displayed in Fig. 6, the IFP rises within the tumor interior somehow proportionally to the tumor density with respect to the isotropic and anisotropic limit coefficients. IFP gradients generate an outward flux at the tumor-host interface that causes the centrifugal diversion of the fluids carrying nutrients and thus contributes to slow down the proliferative potential of inner cells. However, by looking at Fig. 6, the variation in fluid content does not follow the IFP and exhibits an inverse behaviour, the intratumoral fluid loss for γ r = 0 greater than the one observed when γ r = 1/3. Such an impairment is due to the fact that, recalling the poroelastic constitutive relation connecting pore strain, pressure and fluid content, in the former case the pores are much more compressed by the higher stresses and increase the fluid walkway throughout the solid matrix by accelerating cells inactivation in the tumor core. In the isotropic case, stress-driven flow diversion is lower, this also ensuring a more prolonged supply of nutrients to the internal cell layers. Fig. 7 a Calibration of the mechano-sensing feedback on proliferative and apoptotic cell rates on the basis of available literature measurements (data re-adapted from Ref. [13]). b Comparison between proliferative and apoptotic tumor fractions at two weeks, in case of unconfined and confined growth, respectively (with η K = 1.42 −1 and η a = 1.95 −1 ). c Comparison of the theoretically predicted and experimentally obtained variations of apoptotic index as a function of the applied pressure level, calculated for colon carcinoma multicellular spheroids derived from mouse CT26 cell lines (experimental data adapted from Ref. [85])

Some experimental evidences on stress-dependent apoptosis in the competitive interactions
The presented model can be enriched by introducing stressdependent apoptosis terms into the interspecific interplays, which allows to directly compare our theoretical predictions with available literature findings in which it is observed in vitro how compressive stress induces a reduction of cell proliferative potential and, in turn, drives cell towards apoptotic stages [12,13,85]. Therefore, by assuming that stress influences apoptotic rates through a mechanism complementary to the way in which it affects division rates, we re-parametrized the mechano-sensing function S K in Eq. (31) on the basis of the relationship between pressure and cell apoptotic rates provided by Montel et al. [13] for colon carcinoma cell spheroids derived from mouse CT26 cell lineages, in which collected data varied with an analogous trend, as also shown in Fig. 7a. In the competitive system, cell rates can be modified by considering in explicit the mutation of tumor and healthy species into an apoptotic state as intratumoral compression increases. By then indicating with the superscript a the updated functions, each cell specific rate in Eq. (30) becomes Γ K denote the apoptotic counterparts of each phenotype obeying the additional bal- K . By using proper cell-specific intrinsic rates, numerical simulations were carried out by considering a case of isotropic growth and by imposing different pressures at the external boundaries, i.e. σ r (b, t) = −p ext , so as to reproduce analogous stress conditions with respect to those described in the case of multicellular spheroids grown within gel environments under controlled pressures. Theoretical predictions in Fig. 7b show that, inside the tumor, isotropic compression induces a prevalence of the apoptotic fraction, while the active tumor cells move towards the periphery by forming a proliferative ring at the tumor-host interface. These outcomes reproduce the morphological pattern exhibited by a vast variety of solid tumor spheroids that grow in vitro under non-zero stress states by forming apoptotic cores, in which caspase-3 activity is often over-expressed with respect to the surrounding rim, which is instead mainly occupied by dividing cells [86]. The relative apoptotic percentage obtained from computations well replicates the experimental data reported in Fig. 7c and referred to the case of colon carcinoma cell spheroids subjected to different levels of external pressure [13,85]. These results demonstrate that the presented strategy, by employing a minimal set of constitutive parameters and characteristic cell rates, can help to shed light on the way of consistently modelling the interlaced roles of mechanics and interspecific cell competition within tumor microenvironments to predict the development of solid tumors and related residual stress fields.

Conclusions
With the aim to predict the evolution of growth-induced stresses, stress gradients and pressures inside growing solid tumors, mechanical models need to be coupled with dynamical systems able to trace back the most relevant interactions among the biological constituents determining tissue development by accounting for the presence of competition terms and chemo-mechanical feedbacks. In this regard, the Authors proposed a coupled nonlinear theory in which poroelasticity and interspecific Volterra-Lotka dynamics allowed to characterize the growth of spheroidal solid tumors, by providing experimental confirmations of the predicted stresses and remodelling [27,55]. By starting from the proposed approach, the present work proposed a mechanically linear framework in order to reduce the degree of coupling and study the stability of the dynamics related to tumor growth in the threespecies model. The model was also exploited to highlight that the eventual presence of growth anisotropy, depending on how neighbouring cells interact during proliferation by deforming the matrix, can lead to higher intratumoral stresses with consequences on the remodelling of the mass. This result could be further investigated also by resorting to nonlinear simulations and a more complex coupling that provides the use of stress-governed evolution laws of the anisotropic multipliers in order to account for their spatial and temporal inhomogeneities. Moreover, the macroscopic description of tumor progression provided in the present paper could be enriched by including inhomogeneous interspecific coefficients that actually depend on the main biochemical species involved in the signalling path-ways of the mechanotransductive processes that are initiated by physical events such as stress, flow pressure and substrate rigidity [15,16,87]. In particular, by following a multiscale approach, the interplay between molecular and mechanical signals could be modelled by considering ad hoc introduced reaction-based equations tracing the activation of selected biochemical mediators that govern, at the microscale level, internal cell reconfigurations as well as the cell-cell and the cell-environment interactions. Additionally, the same proposed strategy might be applied to test in silico the effectiveness of drug agents in competing dynamics with cells and tissues as well as to design drugs at the scales of interest. All these aspects would lead to a more accurate characterization of tumor mechanical microenvironment that is one of the most challenging aspects of mechanobiology. In fact, the prediction (and understanding) of the response of tumor systems, from the single-cell scale to the mass continuum level, still represents a crucial breakthrough for implementing engineered strategies to selectively target cancer cells as well as for conceiving and optimizing mechanically-based protocols of anti-cancer treatments in the field of the precise medicine for oncology. To this purpose, the modelling of mechanosensing via mechanical cell competition, by including stress-dependent interplays, could contribute to elucidate still unclear dynamics or even to envisage new ways for altering the tumor fate through the use of ad hoc mechanical stimuli, by so reorienting the species evolution towards controlled states, as also shown through comparison of theoretical predictions and experimental evidences. For all these reasons, it is felt that the coupling between Continuum Mechanics, cell competitions and evolutionary laws could help to unveil still concealed mechanisms regulating tumor growth as well as other different multiscale/multiphysics biological processes in which the evolution of living systems within their environment is governed by the interplay between mechanical and biochemical events [61,88,89].