Multi-scale modeling of Snail-mediated response to hypoxia in tumor progression
Giulia Chiari Martina Conte Marcello Delitala





  
  
   



 

 



bbold
AMSb












 























theoremTheorem[section]
corollaryCorollary
*mainMain Theorem
lemma[theorem]Lemma
propProposition
conjectureConjecture
*problemProblem
thmTheorem
oss*Observation
*rem*Remark
definition
defin[theorem]Definition
remRemark
*notationNotation





























|mid
tr
T





sistem
{[ .




,1,2,3]Giulia Chiari[Corresponding authors: martina.conte@polito.it, giulia.chiari@polito.it. These authors contributed equally to the work.]
*,1,4]Martina Conte
1]Marcello Delitala
[1]Department of Mathematical Sciences "G. L. Lagrange", Politecnico di Torino Corso Duca degli Abruzzi 24, 10129 Torino, Italy
[2]Department of Mathematics “G. Peano”, Università di Torino Via Carlo Alberto 10, 10124 Torino, Italy
[3]Department of Mathematics, Swinburne University of Technology John St, Hawthorn VIC 3122, Australia
[4]Division of Mathematical Oncology and Computational Systems Biology Department of Computational and Quantitative Medicine, Beckman Research Institute City of Hope National Medical Center, Duarte, CA, USA
                     


Multi-scale modeling of Snail-mediated response to hypoxia in tumor progression
    [
    April 29, 2024
===============================================================================





Tumor cell migration within the microenvironment is a crucial aspect for cancer progression and, in this context, hypoxia has a significant role. An inadequate oxygen supply acts as an environmental stressor inducing migratory bias and phenotypic changes. In this paper, we propose a novel multi-scale mathematical model to analyze the pivotal role of Snail protein expression in the cellular responses to hypoxia. Starting from the description of single-cell dynamics driven by the Snail protein, we construct the corresponding kinetic transport equation that describes the evolution of the cell distribution. Subsequently, we employ proper scaling arguments to formally derive the equations for the statistical moments of the cell distribution, which govern the macroscopic tumor dynamics. Numerical simulations of the model are performed in various scenarios with biological relevance to provide insights into the role of the multiple tactic terms, the impact of Snail expression on cell proliferation, and the emergence of hypoxia-induced migration patterns. Moreover, quantitative comparison with experimental data shows the model's reliability in measuring the impact of Snail transcription on cell migratory potential. Through our findings, we shed light on the potential of our mathematical framework in advancing the understanding of the biological mechanisms driving tumor progression.

 
  

 
 


§ INTRODUCTION

Migration of tumor cells into the normal tissue under the influence of biochemical and biophysical components of the micro-environment is one of the hallmarks of cancer <cit.>. However, because of the highly complex biology at the cellular and molecular level and in the interactions with the surrounding environment, the exact dynamics driving cell migration are still not completely well understood. Tissue oxygenation is one of the prominent traits in this context. It has been suggested that oxygen concentration highly influences the switch between migrating and proliferating cell behavior, the invasiveness, and aggressiveness of the tumor cells. Moreover, the deprivation of oxygen acts as an environmental stressor, promoting a long series of mutations that strongly impact the tumor dynamics <cit.>.
 
Hypoxia, defined as an insufficient supply of oxygen, has long been recognized as a contributing factor to the tumor microenvironment (TME) <cit.>. Not only it can induce a pronounced migratory bias of the cells towards favorable areas, but it can even determine their phenotype and interaction strategies. It has been clinically observed that, in solid tumors, the oxygen distribution is heterogeneous with oxygen levels ranging from normal to hypoxic and severe hypoxic <cit.>. Under hypoxic conditions, tumor cells undergo morphological and molecular changes to adjust their behaviors and acquire the abilities to adapt to hypoxia and escape apoptosis <cit.>. There are several adaptive responses of tumor cells to hypoxia <cit.>, which may involve the secretion of specific transcription factors, like hypoxia-inducible factor 1 (HIF-1), the upregulation of hypoxia-inducible angiogenic factors, sustaining new vessels formation <cit.>, or the glycolysis activation <cit.>. Among others, HIF-1 activation controls the expression of Snail transcription. 

The Snail superfamily of transcription factors includes Snail1, Slug, and Scratch proteins <cit.>. It is well documented that Snail protein directly represses E-cadherin and thus it is a key inducer of epithelial-mesenchymal transition, a biological process defining progression from a polarised epithelial cell phenotype to a mesenchymal phenotype<cit.>. In addition to regulating epithelial-to-mesenchymal transition (EMT) and cell migration, overexpression of Snail induces resistance to apoptosis and tumor recurrence <cit.>. The Snail-mediated survival of epithelial cells may thus enhance the ability of tumor cells to invade and metastasize. Overexpression of Snail has been reported to be a sufficient inducer of EMT as well as a predictor for an aggressive tumor phenotype. It has recently been demonstrated that Snail expression is induced by hypoxic conditions and is regulated by HIF-1α expression at the transcriptional level <cit.>. Up-regulation of Snail-1 correlates with metastasis and poor prognosis, whereas silencing of Snail-1 is critical for reducing tumor growth and invasiveness <cit.>. Since such complex processes and their mutual conditioning scenarios are difficult to assess experimentally, mathematical models can help understanding the underlying biological mechanisms, test hypotheses, and even make predictions.

Cancer progression is a complex process, involving several factors, and taking place as both an individual and a collective process. Microscopic intracellular dynamics, occurring at the individual level, influence the mesoscopic behavior of the cells, which determines the macroscopic evolution for the cell population density. Previous models for tumor invasion have been proposed in discrete or continuous frameworks. The former are based on the description of individual cell dynamics moving on a lattice <cit.> and can involve continuous equations for the evolution of external factors (e.g., chemoattractant concentration, density of ECM fibers, low pH levels) - the so-called hybrid models <cit.>. Concerning the latter, different classes of fully continuous models for tumor cell migration have been developed. Many of them are versions or extensions of a classical reaction-diffusion model proposed by Murray <cit.>, while more recent works take into account the advection bias of tumor cells describing motility adjustment to extracellular signals. Some of these are directly set on the macroscopic scale and rely on balance equations for mass, flux, or momentum <cit.>, or on integro-differential equations accounting for the development of specific intra-tumor structure <cit.> (see also the review in <cit.> for settings with multiple taxis in the larger context of cell migration). More recently, the use of kinetic transport equations (KTEs) in the kinetic theory of active particles (KTAP) framework has been largely applied to the study of cell migration, in general, <cit.>, and in the specific context of tumor evolution <cit.>. Kinetic models are intrinsically multi-scale models that characterize the dynamics of distribution functions of tumor cells which may depend, besides time and position, on several kinetic variables, such as microscopic velocity or activity variables. These models use Boltzmann-type equations for the cell population density and scaling arguments to derive the macroscopic setting. Among those models, <cit.> accounts for effects of hypoxia or hypoxia-driven acidity on the migration and invasion process of tumor cells. Concerning the taxis terms obtained in this class of models, in <cit.> they are derived from the description of subcellular dynamics for receptor binding, which leads in the mesoscopic KTE to transport terms w.r.t. the activity variables. In <cit.>, instead, the use of turning rates depending on the pathwise gradient of some chemotactic signal leads to various types of taxis at the macroscopic level. Moreover, in <cit.>, forces and stress, acting on the cells and depending on the chemical and physical composition of the environment, translates into transport terms w.r.t. the velocity variable in the corresponding KTE. In particular, in <cit.>, the authors consider a flux-limited description of the transport terms. Flux limitations have been introduced in the modeling of cell motility to reduce the infinite speed of propagation triggered by linear diffusion and the excessive influence of the latter on the spread of cells. Their derivation from KTEs has been provided formally in <cit.> and rigorously in <cit.>. In both cases, the derivation is based on an appropriate choice of the signal response function involved in the turning operator and depends on the directional derivative of the signal. In <cit.>, instead, an alternative approach based on characterizing velocity dynamics at the single-cell level is proposed. Lastly, concerning the mathematical modeling of Snail dynamics, several works have been proposed for theoretically studying Snail's role in the epithelial-to-mesenchymal transition process, especially looking at the interactions among microRNAs and transcription factors - miR-34, miR-200, Zeb, and Snail - at the single-cell level <cit.>, while its connection to cellular motility has not been largely investigated. 
 
In the present work, we propose a multi-scale mathematical modeling approach for describing tumor invasion in response to tissue hypoxia, investigating the interplay between molecular signaling pathways and cell dynamics. The model connects single-cell behavior driven by Snail expression with macroscopic scale dynamics describing tumor migration in the tissue. Starting from the approach proposed in <cit.>, we introduce a novel description of the internal variable dynamics for Snail expression and we account for flux-limited operators in the single-cell velocity dynamics. At the mesoscopic level, cell evolution is described in terms of a classical kinetic transport equation with different formulations for the proliferative operator. From this description, using proper up-scaling arguments, we derive the macroscopic setting. Moreover, we numerically investigate the model's capability of capturing different biologically relevant scenarios concerning hypoxia and Snail effects on cell migration and proliferation. We show how the model is reliable in replicating different experimental results and offers new perspectives for interpreting experimental findings.

The paper is organized as follows. Section <ref> provides the setup of microscopic and mesoscopic equations for the dynamics of tumor cells. Section <ref> contains the derivation of the macroscopic equation for tumor cell density evolution, which features flux-limited chemotaxis towards increasing oxygen concentrations and self-diffusion, as well as a proliferation term modeling the inverse correlation between moving and proliferating cell capability in two possible manner. In Section <ref>, four numerical experiments are proposed to show both the qualitative behavior of the proposed model in different scenarios and its capability to qualitatively replicate experimental data concerning Snail impact of tumor cell migration in different oxygen conditions. Finally, Section <ref> provides a discussion of the main outcome of our model, along with some perspectives.





§ MODELING

 In this note, we propose a multi-scale model for describing tumor progression in response to tissue hypoxia, whose influence on the cancer cells is mediated by Snail dynamics. Following the well-established literature regarding multi-scale models for tumor invasion <cit.>, the model setting proposed here is built using the classical tools and methods of kinetic theory. Our main aim is to obtain a detailed description of the tumor cell dynamics, taking into account the effect of microscopic signaling pathways in the mechanisms of tumor response to hypoxia. 
 
 Starting from the microscopic level of interaction between cells and oxygen, we consider the dynamics of the Snail signaling pathway which is involved in the cell response to hypoxic microenvironmental conditions. Moreover, following <cit.>, we provide a microscopic description of velocity dynamics that depend on oxygen and macroscopic cell density tactic gradient, with Snail expression influencing cell motility. Then, we set up the corresponding kinetic transport equation describing the evolution of the cell distribution
 in relation to the prescribed microscopic dynamics. Performing a proper model upscaling, we obtain the equations for the statistical moments of the cell distribution. These describe the dynamics of tumor cells, which are driven by limited diffusion and oxygen-mediated drift, and the evolution of the average Snail expression of the tumor population.




 §.§ Microscopic scale

At this level, we model the dynamics of the microscopic variable y∈, describing the expression of the Snail protein, and the microscopic cell velocity ∈, both influenced by the oxygen levels. We assume oxygen to be time-independent, i.e., O_2=O_2() has a fixed distribution that does not evolve in time. Possible extensions of this approach are discussed in Section <ref>. 



  §.§.§ Dynamics of Snail protein expression

Concerning Snail protein expression, we model the process of protein synthesization from gene transcription and its regulation depending on the oxygen dynamics. Snail proteins are transcription factors involved in the regulation of hypoxia-driven cell migration and invasion <cit.>. It has been shown in several different types of human cancer that overexpression of Snail induces invasion and metastasis <cit.>. The expression of Snail protein is controlled by the oxygen levels, decreasing when the tumor mass is properly oxygenated. Relying on the description of the temporal evolution of the total level of Snail proposed in <cit.>, here we model its dynamics with the following equation:

    dydt=g_sH(y)H(O_2)-γ_s y ,

where g_s and γ_s represent the basal transcription and degradation rates, respectively, while the functions H(·) model the transcription activation/inhibition mechanisms. We recall here that O_2=O_2(). Generally, the functions H(G) can be described as

    H(G)=1+λ_GG/G_01+G/G_0,

where G is a generic agent influencing the transcription and G_0 its reference value. In particular, λ_G>1 models an activation mechanism supported by the agent G, while λ_G<1 refers to an inhibition mechanism driven by G. In the case of Snail, it has been shown that it has a self-regulatory (inhibition) mechanism <cit.>, which we describe as

    H(y)=11+y/y_0 ,
 
with y_0 maximum Snail expression.

The assumption λ_G=0 is made to avoid overloading the model with further parameters. However, other choices of this function can be included in the model. 

Concerning O_2, it exerts an inhibitory mechanism on the transcription of Snail. However, for the definition of the microscopic dynamics and the derivation of the macroscopic model, we keep a general expression for the function H(O_2). We then specify the value of the parameter λ_O_2 used for the numerical simulations in Table <ref>. Rescaling y/y_0 ⇝ y and g_s/y_0 ⇝ g_s, we simplify the notation as

    dydt=g_s11+yH(O_2)-γ_s y:=G(y,O_2)

with y∈=(0,1). Looking at its quasi-steady state solutions y^*, we observe that

    dydt=0   ⟺ g_s11+yH(O_2)-γ_s y=0
       ⟺[g_sH(O_2)-γ_s y(1+y)]11+y=0
       ⟺ y(1+y)=g_sγ_sH(O_2)
       ⟺ y^2 +y-g_sγ_sH(O_2)=0
       ⟺ y=-12(1±√(1+4g_sγ_sH(O_2))) .

As y represents a biological quantity accounting for the expression of Snail, no negative values are admitted. Thus, the only acceptable steady-state solution for our system is given by

    y^*=12(√(1+4g_sγ_sH(O_2))-1) .


To ensure that the equilibrium distribution y^* belongs to , we have to require that 
    g_sγ_smax_∈ℝ^n{H(O_2())}< 2 .
 If we consider an inhibitory function H(O_2) such that λ_O_2=0, than the condition reads g_s< 2γ_s.




  §.§.§ Dynamics of cell velocity

Concerning the microscopic velocity , we model the mechanism by which cells tend to migrate by aligning to two different gradients. Precisely, increasing gradients of oxygen attract tumor cells towards better oxygenated areas, while tumor cells tend to avoid crowded regions with high cell densities. In both cases, the smaller the amount of Snail expression, the lower the cell's tendency to move along the directions of these gradients. Under these assumptions, the preferred direction of a cell can be modeled by a weighted sum of the two gradients. Thus, velocity dynamics are modeled with the following equation

    ddt=g(y,O_2,M)-a_2 ,

where the function g(y,O_2,M) describes cell acceleration, while the second term models cell deceleration, with a_2 a positive constant scaling cell deceleration. Cells, in fact, tend to slow down or randomly move in the absence of external signals. Concerning cell acceleration, we set

    g(y,O_2,M)=a_1(y,O_2,M) .

Here, (y,O_2,M) is the vector gradient modeling cell alignment along the directions given by the gradient of oxygen O_2 and the gradient of macroscopic tumor cell density M, while a_1 is a positive constant scaling cell acceleration. As introduced above, we assume that the cell's tendency to follow oxygen gradient is enhanced by high Snail expression, which is one of the mechanisms of cell response to hypoxia. At the same time, since low levels of Snail promote high levels of E-cadherin expression <cit.>, which is responsible for cell-cell adhesion, we assume that cell tendency to avoid high cell density region is also positively regulated by y. In fact, high levels of Snail would promote less adhesion between cells and, thus, a more enhanced tendency to escape from the tumor core. Thus, we choose

    (y,O_2,M)=y (β∇_ O_2√((O_2,0X)^2+|∇_ O_2|^2)-(1-β)∇_ M√((K_MX)^2+|∇_ M|^2)) .
 
Here, K_M and O_2,0 are reference values for tumor cells and oxygen, X>0 is a constant to be selected in correspondence to appropriate time and length scales, while the parameter β∈(0,1) weights the contributions of the two tactic terms, depending on the main microenvironmental cue. Thus, (<ref>) can be written as

    ddt=a_1y (β∇_ O_2√((O_2,0X)^2+|∇_ O_2|^2)-(1-β)∇_ M√((K_MX)^2+|∇_ M|^2))-a_2=:S(,y,O_2,M) .

We observe that g(y,O_2,M) is bounded

    |g(y,O_2,M)|=|a_1 |<a_1

and, thus, the speed s=||<a_1a_2:=s_ub, with s_ub an upper bound for cell speed. Finally, we complete the microscopic level system by modeling the changes in the cell position ∈ℝ^n as

    ddt= .

Thus, collecting equations (<ref>), (<ref>), and (<ref>), the complete system for the microscopic level dynamics reads

    ddt= ,
    ddt=a_1 y (β∇_ O_2√((O_2,0X)^2+|∇_ O_2|^2)-(1-β)∇_ M√((K_MX)^2+|∇_ M|^2))-a_2 ,
    dydt=g_s11+yH(O_2)-γ_s y .




 §.§ Mesoscopic scale

At this level, we consider the cell density distribution p(t,,,y):[0,T]×ℝ^n××→ℝ, depending on time t, position ∈ℝ^n, microscopic velocity ∈, and internal variable for Snail protein expression y∈. In particular, the microscopic velocity vector  can be written as =sϑ with cell speed s∈(0,s_ub) and cell direction ϑ∈𝕊^n-1. For describing the mesoscopic dynamics of tumor cells, we consider the following kinetic transport equation:

    ∂ p∂ t +∇_· ( p) +∂∂ y(G(y,O_2)p)+∇_· (S(,y,O_2,M) p)=𝒫[p] .

 Here, the functions G(y,O_2) and S(,y,O_2,M) are given by (<ref>) and (<ref>), respectively, while the operator 𝒫[p] describes the proliferation process. We generally describe it as
 
    𝒫[p]=μ_1(M,O_2,s)∫_μ_2(y')χ(t,,y,y')p(t,,,y')dy' .

 Here, the coefficient function μ_1(M,O_2,s) accounts for the possible effect of cell speed and oxygen level changes, as well as overcrowding of the environment, on cell proliferation. Instead, the integral operator, involving the coefficient function μ_2(y) and the kernel χ(t,,y,y') describes the role of Snail expression in the proliferation process. In particular, χ(t,,y,y') is the probability kernel representing the likelihood of cells to receive a Snail expression regime y after the division of a cell with Snail expression y'. We propose two possible choices for the proliferation term expression, both based on the assumption that cell capabilities of moving and proliferating are inversely correlated (go-or-grow hypothesis <cit.> and tradeoffs <cit.>). 
 
     
  ∙ In the first case, we assume that

    μ_1(M,O_2,s):=μs_ub-s^*s_ub(1-MK_M)O_2O_2,0+O_2    and    μ_2(y):=1     ∀ y∈.

Moreover, we assume that the level of Snail expression in a daughter cell is equal to the one of its mother, i.e., the kernel χ(t,,y,y')=δ(y-y'). With these assumptions, the proliferative operator reads 

    𝒫_1[p]=μs_ub-s^*s_ub(1-MK_M)O_2O_2,0+O_2p(t,,,y)

With this choice, the impact of proliferation-migration tradeoffs is taken into account by relating μ_1(M,O_2,s) to cell speed in a decreasing manner. This allows to account for the impact of y only in an indirect manner, namely through the variation of s. Moreover, as the adaptation of speed to the surrounding environment is faster than the proliferation time, the
velocity is approximated by its steady-state ^* and the corresponding speed is denoted by s^*=|^*|.

  ∙ For the second case, we assume that 

    μ_1(M,O_2,s):=μ(1-MK_M)O_2O_2,0+O_2    and    μ_2(y):=1-y     ∀ y∈.

With this choice, the mentioned tradeoffs are taken into account in the integral term, which models a reduced proliferation for high levels of Snail expression. In particular, we assume that the kernel does not depend on the level of Snail expression of the mother cell, i.e., χ=χ(t,,y), and distribution of the level of Snail expression in a daughter cell is  symmetrical around the quasi-steady-state y^*, i.e.,

    ∫_(y-y^*)χ(t,,y)dy=0 .

With these assumptions, the proliferative operator reads 

    𝒫_2[p]:=μ(1-MK_M)O_2O_2,0+O_2∫_(1-y')χ(t,,y)p(t,,,y')dy'

 
In both cases, the macroscopic cell density M is given by 

    M:=∫_∫_ p(t,,,y)dy d .

The introduced descriptions of the proliferation process are such that proliferation is reduced for highly motile cells, with a direct or indirect effect of Snail expression. In both case, the parameter μ represents the constant tumor proliferation rate and we also include a direct effect of oxygen on tumor cell proliferation, which is enhanced at higher oxygen levels. We will compare, then, the corresponding macroscopic models in the numerical experiment 2 in Section <ref>. 



§ DERIVATION OF MACROSCOPIC SYSTEM

Due to the high dimension of (<ref>), solving directly this kinetic equation has to face several challenges, especially related to its complexity and to a high computational cost. Therefore, in this Section we aim at deducing a macroscopic counterpart of (<ref>). 



 §.§ Assumptions
 
To obtain a closed system of macroscopic equations from the integration of (<ref>) w.r.t y and , we need to make the following assumptions on the moments of the distribution function:


    ∇_·∫_(-^*)(v_i-v_i^*)pd≈ 0 ,      ∫_∫_(v_i-v_i^*)(y-y^*)pdyd≈ 0 ,      ∫_∫_(y-y^*)^2pdyd≈ 0 .
 
Precisely, with v_i we indicate the i-th component of the velocity vector , while y^* and ^* are the steady-state solutions of microscopic equations (<ref>) and (<ref>), respectively. With these assumptions, we state that some of the second-order moments for the tumor cell distribution w.r.t. deviations of  and y from their steady-states are negligible, as well as the second-order moment w.r.t y. These are reasonable choices since the microscopic dynamics of protein expression and velocity changes happen faster in comparison to the kinetic behavior of tumor cells.

Considering the rescaling described in the above Sections, the domains  and  are given by =(0,1) and =B_s_ub^n(0)=(0,s_ub)×𝕊^n-1. Following the approach proposed in <cit.>, we assume 
the distribution p to be compactly supported in the × space. Precisely, for equation (<ref>), boundary conditions w.r.t. these variables need to be prescribed at the inflow boundary of  and . Considering the dynamics in (<ref>), a protein expression state y ∈∂ is part of the inflow boundary if G(y,O_2)· n≤ 0 for n outward normal on the boundary. Given ∂={0,1}, it holds

    G(0,O_2)· n=g_sH(O_2)· (-1)<0
       G(1,O_2)· n=(2 g_sH(O_2)-γ_s)· (1)≤ 0   if condition (<ref>) holds.

Thus, the inflow boundary of  coincides with ∂ and boundary conditions can be prescribed on the whole ∂. Instead, considering the dynamics in (<ref>) and a velocity vector ∈∂, we have that ||=s_ub and the outward normal n=/s_ub. The velocity vector  is part of the inflow boundary if S(,y,O_2,M)· n≤ 0, i.e.,

    S(,y,O_2,M)· n   =a_1s_ub(  b(y,O_2,M)·) -a_2s_ub( ·)
       ≤a_1s_ub| b(y,O_2,M)| ||-a_2s_ub||^2
    
         ≤ a_1-a_2 s_ub=0 .

Thus, the inflow boundary of  coincides with ∂ and boundary conditions can be prescribed on the whole ∂. Precisely, (<ref>) and (<ref>) allow to conclude that the characteristics of the transport part of equation (<ref>) that start in ℝ^n×× do not leave this set.
 


 §.§ Equation for the moments

Here, we upscale (<ref>) to obtain the equation for the macroscopic cell density M(t,). Firstly, we rescale the quantities introduced above as p/K_M ⇝ p, M/K_M ⇝ M, O_2/O_2,0⇝ O_2, and s^*/s_ub⇝ s^*. Then, we introduce a small parameter ε≪ 1 to rescale time and space as


    t̂=ε^κ t ,  =ε ,
       ĝ_s=ε^-νg_s ,  γ̂_s=ε^-νγ_s ,

with κ, ν >0. The rescaling of the reaction rates g_s and γ_s means the rescaling of dy/dt, while the negative power of ε is chosen to reflect the fact that these dynamics are the fastest among all included processes. Moreover, assuming that 1/X is of order ε, we observe that ^*=^*. For simplicity of writing, we drop the hat symbol from all variables and, thus, equation (<ref>) reads

    ε^κ∂ p∂ t +ε∇_· ( p) +ε^-ν∂∂ y(G(y,O_2)p)+∇_· (S(,y,O_2,M)p)=ε^κ𝒫_k[p]

where 

    G(y,O_2)=g_s11+yH(O_2)-γ_sy ,
       S(,y,O_2,M)=a_1y (β∇_ O_2√(1+|∇_ O_2|^2)-(1-β)∇_ M√(1+|∇_ M|^2))-a_2 ,

and 𝒫_k[p] would be given by either

    𝒫_1[p]=μ(1-s^*)(1-M)O_21+O_2p(t,,,y)

or

    𝒫_2[p]:=μ(1-M)O_21+O_2∫_(1-y')χ(t,,y)p(t,,,y')dy' ,

which corresponds to the rescaled versions of (<ref>) or (<ref>), respectively. Together with the introduced macroscopic tumor density (<ref>), let consider the following notations for the moment of the distribution function p:


    m(t,,y):=∫_ p(t,,,y)d ,             m^(t,,y):=∫_ p(t,,,y)d , 
       m_i^(t,,y):=∫_ v_i p(t,,,y)d ,
                   M^_i(t,)=∫_∫_ v_i p(t,,,y)d dy ,
       M^(t,)=∫_∫_ p(t,,,y)d dy ,        M^y(t,)=∫_∫_ yp(t,,,y)d dy .


Considering (<ref>) and integrating w.r.t  we derive the equation for m(t,,y), i.e.,

    ε^κ∂ m∂ t +ε∇_· m^ +ε^-ν∂∂ y(G(y,O_2)m)+∫_∇_· S(,y,O_2,M) p d=ε^κ∫_𝒫_k[p]d ,

Here

    ∫_∇_·  S(,y,O_2,M) p  d=0

for the boundary conditions imposed on , while for the proliferative operator we have either

    ∫_𝒫_1[p]d=∫_μ(1-s^*)(1-M)O_21+O_2p(t,,,y)d= μ(1-s^*)(1-M)O_21+O_2 m(t,,y)=𝒫_1[m]

or

    ∫_𝒫_2[p]d   =∫_μ(1-M)O_21+O_2∫_(1-y')χ(t,,y)p(t,,,y')dy'd
    
       =μ(1-M)O_21+O_2∫_(1-y')χ(t,,y)m(t,,y')dy'=𝒫_2[m] .

Thus, the equation for m(t,,y) reads

    ε^κ∂ m∂ t +ε∇_· m^ +ε^-ν∂∂ y(G(y,O_2)m)=ε^κ𝒫_k[m] .

Considering again (<ref>), multiplying it by v_i and integrating w.r.t  we derive the equation for m_i^(t,,y), which is the i-th component of the vector m^, i.e.,


    m^:=
    [ m_1^; m_2^;  ...; m_n^ ]
    .

The equation for m_i^ reads


    ε^κ∂ m_i^∂ t +ε∫_ v_i∇_· ( p)d +ε^-ν∫_ v_i ∂∂ y(G(y,O_2)p)d+∫_ v_i ∇_· (S(,y,O_2,M) p)d=ε^κ∫_ v_i 𝒫_k[p]d .

Here,

    ∫_ v_i∇_· ( p)d   =∇_·∫_ v_i ( p)d=∇_·∫_ (v_i-v_i^*)(-^*)pd+∇_·∫_ (v_i^*+v_i^*-v_i^*^*)pd
       =∇_· (^*m_i^+v_i^*m^-v_i^*^*m)

for the assumption in Section <ref>, while


    ∫_ v_i ∇_· (S(,y,O_2,M) p)d   =∫_ v_i ∂_v_i(S_ip)d+∑_j=1
     j i^N ∫_ v_i∂_v_j(S_jp)d
       =∫_⋃_j=1
    j i^N_j∫__i v_i ∂_v_i(S_ip) dv_i  d+∑_j=1
     j i^N  ∫_⋃_k=1
    k j^N_k v_i∫__j∂_v_j(S_jp) dv_j  d

where d is (n-1)-th components vector such that 


    =(v_i,)∈_i×⋃_j=1
    j i_j= .
 
Under the boundary conditions in Section <ref>, the second term on the right-hand-side is equal to 0, while using the chain rule the first term reduces to


    ∫_⋃_j=1
    j i^N_j∫__i v_i ∂_v_i(S_ip) dv_i  d   =∫_⋃_j=1
    j i^N_j[v_i (S_i p)_|∂_i-∫__iS_ipdv_i]d=-∫_ S_ipd
    
       =- a_1 y (β(∇_ O_2)_i√(1+|∇_ O_2|^2)-(1-β)(∇_ M)_i√(1+|∇_ M|^2))m+a_2m_i^
       =-g_i(y,O_2,M)m+a_2m_i^ ,

where g_i(y,O_2,M) represents the i-th component of the vector function g(y,M,O_2) defined in (<ref>). Concerning the proliferative operator, we have either

    ∫_ v_i𝒫_1[p]d=∫_ v_i μ(1-s^*)(1-M)O_21+O_2p(t,,,y)d= μ(1-s^*)(1-M)O_21+O_2 m_i^(t,,y)=𝒫_1[m_i^]

or

    ∫_ v_i𝒫_2[p]d   =∫_ v_iμ(1-M)O_21+O_2∫_(1-y')χ(t,,y)p(t,,,y')dy'd
    
       =μ(1-M)O_21+O_2∫_(1-y')χ(t,,y)m_i^(t,,y')dy'=𝒫_2[m_i^] .

Thus, the equation for m_i^ reads

    ε^κ∂ m_i^∂ t +ε∇_· (^*m_i^+v_i^*m^-ε v_i^*^*m) +ε^-ν∂∂ y(G(y,O_2)m_i^)-(g_i(y,O_2,M)m-a_2m_i^)=ε^κ𝒫_k[m_i^] .

Therefore, the system for the n+1 variables (m,m_1^,m_2^,..,m_n^) is given by

    ε^κ+ν∂ m∂ t(t,,y) +ε^1+ν∇_· m^(t,,y) +∂∂ y(G(y,O_2)m(t,,y))=ε^κ+ν𝒫_k[m](t,,y) ,
    ε^κ+ν∂ m_i^∂ t(t,,y)    +ε^1+ν∇_·[^*m_i^+v_i^*m^-v_i^*^*m](t,,y) +∂∂ y[G(y,O_2)m_i^(t,,y)]
       =ε^ν[g_i(y,O_2,M)m(t,,y)-a_2m_i^(t,,y)]+ε^κ+ν𝒫_k[m_i^](t,,y)      ∀ i=1...n .

We remark that ∇_· m^=∑_j=1^n∂_x_jm_j^. Then, integrating (<ref>) w.r.t y, for the assumptions in Section <ref>, we immediately get

    ε^κ∂ M∂ t(t,)+ε∇_· M^(t,)=ε^κ∫_𝒫_k[m]

where for the proliferative operator we have either

    ∫_𝒫_1[m]dy=∫_μ(1-s^*)(1-M)O_21+O_2m(t,,y)dy= μ(1-s^*)(1-M(t,))O_21+O_2M(t,)

or

    ∫_𝒫_2[m]dy   =∫_μ(1-M(t,))O_21+O_2∫_(1-y')χ(t,,y)m(t,,y')dy'dy
    
       =μ(1-M)O_21+O_2∬_(1-y')χ(t,,y)m(t,,y')dy'dy
       =μ(1-M(t,))O_21+O_2(M(t,)-M^y(t,)) .

Instead, integrating (<ref>) w.r.t. y, we obtain 


    ε^κ∂ M_i^∂ t   +ε∇_·(^*M_i^+v_i^*M^-v_i^*^*M) +ε^-ν∫_∂∂ y(G(y,O_2)m_i^)dy-∫_ (g_i(y,O_2,M)m-a_2m_i^)dy
       =ε^κ∫_𝒫_k[m_i^]dy ,

where the first integral vanishes for the boundary conditions in Section <ref>. Concerning the second integral on the left-hand-side, we have

    ∫_ (g_i(y,O_2,M)m-a_2m_i^)dy=g̃_i(O_2,M)M^y-a_2M_i^ ,

where

    g̃_i(O_2,M)=a_1 (β(∇_ O_2)_i√(1+|∇_ O_2|^2)-(1-β)(∇_ M)_i√(1+|∇_ M|^2)) .

Instead, the proliferative operator reads either

    ∫_𝒫_1[m_i^]dy=∫_μ(1-s^*)(1-M)O_21+O_2m_i^(t,,y)dy= μ(1-s^*)(1-M)O_21+O_2 M_i^(t,)

or

    ∫_𝒫_2[m_i^]dy   =∫_μ(1-M)O_21+O_2∫_(1-y')χ(t,,y)m_i^(t,,y')dy'dy
    
       =μ(1-M)O_21+O_2(M_i^-∬_ y'χ(t,,y)m_i^(t,,y')dy'dy)
       =μ(1-M)O_21+O_2(M_i^-∫_∫_ (v_i-v_i^*)y p(t,,,y)dyd-v_i^*∫_∫_ y p(t,,,y)dyd)
       =μ(1-M)O_21+O_2(M_i^-∫_∫_ (v_i-v_i^*)(y-y^*) p(t,,,y)dyd
           -y^*∫_∫_(v_i-v_i^*)p(t,,,y)dyd-v_i^*∫_∫_ y p(t,,,y)dyd)
       =μ(1-M(t,))O_21+O_2(M_i^(t,)-y^*M_i^(t,)-v_i^*M^y(t,)+y^*v_i^*M(t,)) .

Therefore, we obtain

    ε^κ∂ M_i^∂ t(t,)    +ε∇_·[^*M_i^+v_i^*M^- v_i^*^*M](t,) -[g̃_i(O_2,M)M^y(t,)-a_2M_i^(t,)]=ε^κ∫_𝒫_k[m_i^]dy .

Finally, multiplying (<ref>) by y and integrating w.r.t. y we obtain 


    ε^κ∂ M^y∂ t+ε∇_·∫_ y M^ dy +ε^-ν∫_ y∂∂ y(G(y,O_2)m) dy=ε^κ∫_ y 𝒫_k[m] (t,,y)dy .

Here,


    ∇_·∫_ y M^ dy=∇_·∫_∫_ y ( p)d dy   =∇_·(∫_∫_ (-^*)(y-y^*)pd dy +∫_∫_ (^*y+ y^*-^*y^*)pd dy) 
       =∇_·(^*M^y+y^*M^-^*y^*M)

thanks to the assumption in Section <ref>, while


    ∫_ y∂∂ y(G(y,O_2)m) dy   =(y(Gm)_|∂-∫_G(y,O_2)mdy)=-∫_G(y,O_2)mdy=-∫_(g_s11+yH(O_2)-γ_sy)mdy
       =γ_s M^y-g_sH(O_2) ∫_11+ym dy .

Considering the Taylor expansion of 11+y around y^* we get

    11+y=11+y^*+y-y^*≈ (1+y^*)^-1-(1+y^*)^-2(y-y^*)+(1+y^*)^-3(y-y^*)^2+O((y-y^*)^2) .

Thus, ignoring the higher-order terms of this expansion, we obtain


    ∫_11+ym dy   =∫_((1+y^*)^-1-(1+y^*)^-2(y-y^*)+(1+y^*)^-3(y-y^*)^2)m dy
       =(1+y^*)^-1 M- (1+y^*)^-2 M^y+y^*(1+y^*)^-2M ,

i.e.,


    ∫_ y∂∂ y(G(y,O_2)p) dy   =γ_sM^y-g_sH(O_2)1+y^*(M-11+y^*(M^y-y^*M))
       =-g_sH(O_2)(1+y^*)^2(1+2y^*)M+(γ_s+g_sH(O_2)(1+y^*)^2)M^y .

The proliferative operator, instead, reads either

    ∫_ y𝒫_1[m]dy=∫_μ(1-s^*)(1-M)O_21+O_2ym(t,,y)dy= μ(1-s^*)(1-M)O_21+O_2 M^y(t,)

or


    ∫_ y 𝒫_2[m]dy   =∫_μ(1-M)O_21+O_2∫_ y(1-y')χ(t,,y)m(t,,y')dy'dy
       =μ(1-M)O_21+O_2(∫_ yχ(t,,y)dy∫_ (1-y')m(t,,y')dy')
       =μ(1-M)O_21+O_2(M-M^y)(∫_ (y-y^*)χ(t,,y)dy+y^*∫_χ(t,,y)dy)
       =μ(1-M(t,))O_21+O_2y^*[M(t,)-M^y(t,)] .

Thus, the equation for M^y reads

    ε^κ∂ M^y∂ t(t,)    +ε∇_·[^*M^y+y^*M^-^*y^*M](t,) -ε^-νg_sH(O_2)(1+y^*)^2(1+2y^*)M(t,)
       +ε^-ν(γ_s+g_sH(O_2)(1+y^*)^2)M^y(t,)=ε^κ∫_ y𝒫_k[m](t,,y)dy .

Thus, the system for M(t,), M^(t,), and M^y(t,) reads

    ε^κ-1∂ M∂ t(t,)+∇_· M^(t,)=ε^κ-1∫_𝒫_k[m](t,,y)dy ,
       ε^κ∂ M_i^∂ t(t,)+ε∇_·[^*M_i^+v_i^*M^- v_i^*^*M](t,) -(g̃_i(O_2,M)M^y(t,)-a_2M_i^(t,))
       =ε^κ∫_𝒫_k[m_i^](t,,y)dy ,                ∀ i=1...n ,
       ε^κ+ν∂ M^y∂ t(t,)  +ε^1+ν∇_·[^*M^y+y^*M^-^*y^*M](t,) -g_sH(O_2)(1+y^*)^2(1+2y^*)M(t,)
       +( γ_s+g_sH(O_2)(1+y^*)^2)M^y(t,)=ε^κ+ν∫_ y𝒫_k[m](t,,y)dy .

Now, we consider the expansion of the previously introduced moments in the form


    M=M_0+ε M_1+O(ε^2) ,
       M^=M_0^+ε M_1^+O(ε^2) ,
       M^y=M^y_0+ε M^y_1+O(ε^2) ,

such that M_0=lim_ε→ 0M. As the focus of this study is on the tumor response to hypoxia and, thus, the influence of the environmental chemotactic cue given by the oxygen levels, cell movement has a very clear directional component. Therefore, following the well-established literature <cit.>, we consider here a hyperbolic limit of the moment system (<ref>), i.e., we choose κ=1. Passing formally to the limit ε→ 0 in (<ref>), from the first equation we get

    ∂ M_0∂ t(t,)+∇_· M_0^(t,)=∫_𝒫_k[m_0](t,,y)dy ,

where m_0 is the zero-order term in the expansion of the moment m(t,,y). Then, from the equation for M_i^ in (<ref>), we get

    g̃_i(O_2,M_0)M_0^y-a_2M_i,0^=0   ⟺    M_i,0^=g̃_i(O_2,M_0)a_2M_0^y

that, considering the vector M_0^ with i-th component M_i,0^, means 

    M_0^=g̃(O_2,M_0)a_2M_0^y

with 

    g̃(O_2,M_0)=a_1 (β∇_ O_2√(1+|∇_ O_2|^2)-(1-β)∇_ M_0√(1+|∇_ M_0|^2)) .

Finally, from the last equation in (<ref>), we obtain

    -g_sH(O_2)(1+y^*)^2(1+2y^*)M_0+(γ_s+g_sH(O_2)(1+y^*)^2)M_0^y=0 ,  ⟺   M_0^y=g_sH(O_2)(1+2y^*)γ_s(1+y^*)^2+g_sH(O_2)M_0 .

Thus, plugging (<ref>) into (<ref>), and the latter in (<ref>), we obtain the following equation for the evolution of the macroscopic cell density M_0

    ∂ M_0∂ t(t,)+∇_·[ F(y^*,O_2())M_0(t,)(β∇_ O_2()√(1+|∇_ O_2()|^2)-(1-β)∇_ M_0(t,)√(1+|∇_ M_0(t,)|^2)) ]=ℙ_k[M_0](t,)

where we set the velocity field  (or chemotactic sensitivity)

    F(y^*,O_2()):=a_1a_2g_sH(O_2())(1+2y^*)γ_s(1+y^*)^2+g_sH(O_2()) ,

while the proliferative operator is given by either

    ℙ_1[M_0](t,)=∫_𝒫_1[m_0](t,,y)dy=μ(1-s^*)(1-M_0(t,))O_2()1+O_2()M_0(t,)

or

    ℙ_2[M_0](t,)=∫_𝒫_2[m_0](t,,y)dy=μ(1-M_0(t,))O_2()1+O_2()(1-a_2a_1F(y^*,O_2()))M_0(t,) .

Since we are interested in the impact of the microscopic Snail protein expression of a cell on the overall macroscopic evolution of the population, for the simulation described in Section <ref> we also analyze the evolution of M_0^y that is linked to the average expression of the Snail protein in the cell population. Thus, discarding the subscripts _0, the resulting macroscopic system reads

    ∂ M∂ t(t,)+∇_·[ F(y^*,O_2())M(t,)(β∇_ O_2()√(1+|∇_ O_2()|^2)-(1-β)∇_ M(t,)√(1+|∇_ M(t,)|^2)) ]=ℙ_k[M](t,) ,
    
    M^y(t,)=a_2a_1F(y^*,O_2())M(t,) .





§ NUMERICAL EXPERIMENTS

We perform 2D numerical simulations of the resulting macroscopic system (<ref>) to analyze in-silico scenarios of tumor progression under varied oxygen conditions and Snail expressions. Numerical simulations are conducted using a self-developed code in Python 3.10.12, whose details are provided in Section <ref>. The parameter values used for the simulations are reported in Table <ref>. Sections <ref>, <ref>, <ref>, and <ref> present the results of the four main numerical experiments.

    
  
     Experiment 1 - In Section <ref>, we analyze the differences in cell migration and Snail distribution over time between a chemotactic dominated scenario and an anti-crowding dominated scenario by varying the value of the weighting parameter β and, thus, its impact on cell motion.
    
    
  
    Experiment 2 - In Section <ref>, we investigate the impact of Snail expression on cell proliferation, 
    and we compare the model's evolution using the two different choices of the proliferative operator given in (<ref>) and (<ref>).
    
    
   
   Experiment 3 -  In Section <ref>, we show our model capabilities to quantitatively replicate experimental results from two different studies on human cancer cells. Firstly, we consider the findings of <cit.> regarding the effect of Snail over-expression or knockdown on the migration capability of human breast cancer cells. Secondly, we refer to the comparison presented in <cit.>, where the effects of hypoxia and a combination of hypoxia and Snail knockdown on the motility of human hepatocarcinoma cells are studied.
    
    
   
    Experiment 4 - In Section <ref>, we analyze the hypoxia-driven spatial distribution of Snail expression within a tumor and we show its consistency with the results shown in <cit.>.







 §.§ Numerical method

To perform numerical simulations of the model, we adapt the numerical method presented in <cit.> to our problem structure. In detail, we rewrite the first equation of (<ref>) as

    ∂ M∂ t(t,)= - T(y^*,O_2,M)(t,) +ℙ_k[M](t,)

where

    (y^*,O_2,M)(t,) = ∇_·[ F(y^*,O_2())M(t,)(β∇_ O_2()√(1+|∇_ O_2()|^2)-(1-β)∇_ M(t,)√(1+|∇_ M(t,)|^2)) ]

rules the movement of cells, while ℙ_k[M](t,) is the proliferation term defined in (<ref>) (k=1) and (<ref>) (k=2). Setting =(x_1,x_2), we consider the geometric domain Ω = [x_1,min, x_1,max]×[x_2,min,x_2,max]⊆ℝ^2, where we introduce a uniform Cartesian mesh consisting of the cells
C_j,k := [x_1,j-1/2, x_1,j+1/2] × [x_2,k-1/2, x_2,k+1/2], for j=0,…, N_x_1 and for k=0,…, N_x_2, of size Δ x_1   ×Δ x_2. We adopt a splitting method, accounting first for the conservative part (y^*,O_2,M) and, then, for the reaction term ℙ_k[M]. Precisely, defining 


    M_j,k(t)= 1/Δ x_1Δ x_2∫_C_j,k M(t,x_1,x_2) d𝐱 ,
 for the conservative part we adopt the general semi-discrete finite-volume scheme given by


    T_j+1/2,k = u^+_j+1/2,kM^E_j,k + u^-_j+1/2,kM^W_j+1,k


    T_j,k+1/2 = v^+_j,k+1/2M^N_j,k + v^-_j,k+1/2M^S_j,k+1 .

Here:

    
  * (·)^+ and (·)^- indicate the positive and negative part of their arguments, respectively, i.e., (·)^+= max{0,·} and (·)^-= min{0,·};
    
    
  * the apices E, W, N, S indicate East, West, North, and South and correspond to the evaluation of the piecewise linear reconstruction using the following first-order truncation of Taylor expansion
    
    M(x_1,x_2) = M_j,k+ (∂_x_1 M)_j,k(x_1-x_1,j) + (∂_x_2 M)_j,k(x_2-x_2,k),    (x_1,x_2)∈ C_j,k

    at the cell interfaces (x_1,j+1/2,x_2,k), (x_1,j-1/2,x_2,k), (x_1,j,x_2,k+1/2),
    (x_1,j,x_2,k-1/2), respectively;
    
    
  * defining
    
    U_j,k = F(y^*(x_1,j,x_2,k),O_2(x_1,j,x_2,k))(β∇_ O_2(x_1,j,x_2,k)√(1+|∇_ O_2(x_1,j,x_2,k)|^2)-(1-β)∇_ M(t,x_1,j,x_2,k)√(1+|∇_ M(t,x_1,j,x_2,k)|^2)) ,

    
    then u := U_x_1 and v := U_x_2 are the components of U along the horizontal and vertical axis respectively.

Note that the derivatives in the middle points are evaluated as


    (∂_x_1 M)_j+1/2,k = M_j+1,k-M_j,k/Δ x_1 ,     (∂_x_2 M)_j,k+1/2 = M_j,k+1-M_j,k/Δ x_2
 ,
while the derivatives in the nodes are initially evaluated as


    (∂_x_1 M)_j,k = M_j+1,k-M_j-1,k/2Δ x_1,     
    (∂_x_2 M)_j,k = M_j,k+1-M_j,k-1/2Δ x_2

and then corrected using a generalized minmod limiter in order to preserve the positivity of the linear reconstruction M (for further details see <cit.>). 
For the time discretization, we use the forward Euler method. We denote with apex l the discretized time step, i.e., 

    t^l = t_0 +∑_i=1^l-1Δ t_i .

To optimize the performances, we use adaptive time steps obtained by imposing the positivity-preserving CFL 

    Δ t_l ≤Δ𝒯_l:=min{Δ x_14a,Δ x_24b}

where a=j,kmax(|u^l_j+1/2,k|)
and b=j,kmax(|v^l_j,k+1/2|). Therefore, starting from the discretized initial condition M_j,k^0 provided for each j=0,…, N_x_1 and for k=0,…,N_x_2, given that M_j,k^l is the numerical approximation of M_j,k(t^l), the numerical scheme reads

    M^l+1/2_j,k =
    M^l_j,k-Δ t_lΔ x_1(T^l_j+1/2,k - T^l_j-1/2,k) - Δ t_lΔ x_2(T^l_j,k+1/2 - T^l_j,k-1/2) 
    
    M^l+1_j,k  = M^l+1/2_j,k + Δ t_l (ℙ_k)^l+1/2_j,k

for l=1,…,N_l. In the proposed experiments, we set the spatial domain Ω =  [0,50] × [0,50]  mm^2 and we consider the time t ∈ [0,T], with T>0. Dealing with a limited domain, we set no entry flux boundary conditions.



 §.§ Experiment 1: chemotactic or anti-crowding dominated motion

In this first experiment, we simulate our model with the aim of comparing cell behaviors in two different scenarios. We consider a first scenario in which cell movement is dominated by a chemotactic attraction toward increasing oxygen concentrations, strongly reducing the impact of the natural cell anti-crowding mechanism. Instead, in a second scenario, we emphasize the role of anti-crowding dynamics, which helps cells to avoid highly dense regions.

Considering the macroscopic system (<ref>), we note that the parameter β∈(0,1) impacts the motility behavior of cells, weighting the influence of oxygen and cell density gradients on the direction of tumor cell migration. Higher values of β imply a stronger impact of chemotaxis compared to anti-crowding dynamics. Thus, to capture the two described scenarios, we choose β=0.98 to emphasize the role of chemotactic movement toward increasing oxygen concentrations, while β=0.8 to account for a stronger effect of the anti-crowding mechanism. In our experiment, we consider a fixed oxygen source located in the top-right corner of the domain Ω and whose expression is given by
 


    O_2(𝐱) = A_O_2 e^-(𝐱-𝐱_O_2)^2θ^2_O_2


with A_O_2=0.8, 𝐱_O_2 = [45, 45] mm, and θ_O_2=34 mm. 
For the tumor cells, we assume that the initial tumor mass M(t=0,):=M_0() is located in the opposite (bottom-left) corner of the domain Ω and it is defined as


    M_0(𝐱) = A_M e^-(𝐱-𝐱_M)^2θ^2_M

with A_M=0.9, 𝐱_M = [9, 9] mm, and θ_M=3 mm. The initial configuration of this setting is shown in Figure <ref>.


To better characterize the dynamics we divide the domain into four different zones depending on the relative oxygen conditions, which trigger cell motility and proliferation. Figure <ref> and the related Table <ref> summarize the combination of the environmental conditions in each of these areas along the 1D section of the domain.

   
 
The results of this first experiment are illustrated in Figure <ref>. Precisely, the left column refers to the simulations obtained for β=0.98 (chemotactic dominated scenario), while the right column to β=0.8 (anti-crowding dominated scenario). 

The simulations track the evolution of the tumor mass over a time range of approximately T=21 days.
In the first row of Figure <ref>, we illustrate the progression of the tumor mass in the domain Ω at four different temporal steps. The contour plot illustrates the density map, while the contour lines represent the tumor's edge (defined by a density threshold corresponding to 10% of the carrying capacity). To better visualize the differences in the tumor evolution in the two scenarios we consider a 1D section of the domain Ω illustrating the tumor profile along the bisecting line. Moreover, in order to account for the epigenetic trait information, for every value y ∈ and every time t, we evaluate the fraction of tumor mass having an equilibrium Snail expression y^* lower than y, i.e.,

    C_y^*(t, y) = ∫_D_yM̅(t,) d

where D_y={∈Ω:y^*()∈(-∞,y)} and M̅(t,) is the normalized tumor density distribution, and we introduce the quantity

    c_y^*(t, y)=∂_y C_y^*(t,y) ,

which provides an indication of the fraction of mass that, at a given time t, has a certain Snail expression y ∈. For this quantity, we evaluate mean η(t), mode ψ(t), and standard deviation σ(t) as


    η(t) = ∫_ y c_y^*(t,y) dy  ,
       ψ(t) = _y∈(c_y^*(t,y)) , 
       σ(t) = √(∫_ y^2 c_y^*(t,y) dy - η(t)^2) .

 Integrating the information about the environmental conditions (Figure <ref>) with the data concerning the expression of Snail (third row of Figure <ref>), we can better grasp the features characterizing the two depicted scenarios. Depending on β, cell behaviors towards more oxygenated regions show evident differences. Initially, low oxygen concentration corresponds to high Snail expression and they collectively contribute both to a high tactic sensitivity F(y^*,O), driving cell drift, and a low proliferation rate. Thus, as the mass is initially situated in an area characterized by low oxygenation, this oxygen deprivation triggers cell motility. As the mass progresses towards the upper right corner (temporal step t_1), the low oxygen levels still induce a more motile than proliferative cell phenotype, enforced by a high mean η(t) of c_y^*(y,t) (depicted in the third row of Figure <ref>). By the time t_1, the mass has spread enough to develop a smoother profile, resulting in a reduced density gradient in favor of stronger chemotactic motion due to the oxygen gradient. This is particularly evident for β=0.98. Cells are, thus, accelerated and move more compactly towards the oxygen source, and this is evidenced in both columns of Figure <ref> by the reduction in mass width orthogonal to the chemotactic gradient between times t_1 and t_2, confirming cell convergence towards the oxygen Gaussian distribution. As the cells reach areas closer to the source (temporal step t_2), a still strong oxygen gradient is countered by increased oxygenation levels, and, thus, inhibition of motility in favor of an enhanced proliferation rate. This results also in a lower level of η(t) (third row of Figure <ref>). This combination results in masses developing a distinct tail to the left, with some cells remaining outside the region of orderly motility due to lower oxygen levels. Many motile cells originating from less oxygenated areas "push" against a slowing front where cells are less motile, but contribute to the increase of the density due to their proliferative capability. At the final time t_3, the mass has reached a region characterized by high oxygenation levels with an almost negligible oxygen gradient. At this stage, ordered motile dynamics become nearly absent and the mean η(t) has reached a lower stable value. This is particularly evident in the scenario with β=0.98, where the mass tends to regain a more radial symmetry and slightly expands under the influence of density pressure driven by proliferative dynamics and anti-crowding.

From a qualitative viewpoint, in both scenarios the dynamics are initially characterized by a diffusing mass that expands orthogonally to the direction of chemotactic motion (more evident for β = 0.8) and then by a direct movement towards the oxygen source. Moreover, in both cases, the presence of flux-limited operators in the drift term determines steep and well-defined invasion fronts, reducing the typical artificial tails characterizing linear diffusion and, thus, its excessive influence on cell spread. However, comparing the experiments conducted for the two scenarios makes it evident how conditions favoring anti-crowding dynamics can lead to significant changes in tumor shape even for small variations in β. In fact, when anti-crowding is the dominant mechanism, cells tend to move toward the location of the oxygen source with a large spread in the domain, with respect to the chemotactic dominated scenario, which shows cells compactly migrating towards more oxygenated regions. Looking at the tumor profiles at time t_1 and t_2 (second row in Figure <ref>), the differences in height and the size of the mass support are evident. This is reflected also in the evolution of mean, mode, and standard deviation of c_y^*(y,t). In fact, for higher values of β, mean η and mode ψ show a similar trend over time. However, the lower β, the greater the differences in their evolution. This is because the larger spread of the tumor mass and the slower cell movement in the domain keep the value of the mean higher for a longer time. Moreover, this determines the wider variety of values covered by the distribution and, thus, a larger standard deviation. 





 §.§ Experiment 2: impact of Snail expression on cell proliferation

In the second experiment, we focus on the proliferative dynamics characterizing tumor cells. Here, we concentrate on the scenario where chemotaxis drives cell motility and compare two proliferative models: ℙ_1, introduced in (<ref>), and ℙ_2, introduced in (<ref>). It is worth recalling that the shared elements in these two modeling choices are: (i) a proliferative rate inversely proportional to cell density; (ii) an increase in proliferative activity correlating with higher levels of available oxygen; (iii) the assumption that cells capabilities of moving and proliferating
are inversely correlated. What distinguishes these approaches is the epigenetic or phenotypic characterization of the duality between cells' motility and proliferative dynamics. In the case of ℙ_1, the factor (1-s^*) ensures a direct correlation between higher Snail expression and lower proliferative activity. Instead, the environmental factor has an indirect impact on the trade-off characterization, as oxygen density influences proliferative activity only indirectly by determining Snail expression. In contrast, for ℙ_2, the factor 1-a_2/a_1 F(y^, O_2()) ensures that the proliferative activity decreases as F(y^*, O_2()) increases, which is directly proportional to Snail expression and inversely proportional to oxygen density.

To compare these two proliferative choices, we conduct two simulations under identical environmental and initial conditions for tumor and oxygen, as well as for the same values of the model parameters. The only difference between the simulations lies in the formulation of the proliferative term. To quantify the results, we consider the 1D section of the domain Ω along the bisecting line (similar to the previous experiment) and we show the difference between the tumor densities resulting from (<ref>) with ℙ_1 (M_ℙ_1(t,)) or ℙ_2 (M_ℙ_2(t,)) at four equally spaced time points: t_1=5.2d, t_2=10.4d, t_3=15.6d, and t_4=20.8d. Moreover, defining 

    Q(t) = ∫_ΩM(t,) d

as the total amount of tumor cells in the spatial domain, we introduce the percentage mass increment of tumor mass from the initial configuration as 


    ℐ_%(t)= 100·(Q(t)-Q(0)Q(0))

with Q(0):=Q(t=0), we compute the center of mass of the tumor 

    Υ(t) = ∫_ΩM̅(t,) d

with M̅(t,) = M(t,) / Q(t), and the velocity of the center of mass

    𝐕(t) = ∂_t Υ(t) .

Results of this experiment are shown in Figure <ref>. Precisely, the top row illustrates the evolution of the difference between M_ℙ_1(t,) and M_ℙ_2(t,) along the 1D section. In the bottom row, we depict the temporal evolution of the mass increase ℐ_%, defined in (<ref>), (left plot) and the velocity of the center of mass, defined in (<ref>), along the bisecting line (right plot). 

 

Analysis of the plot in the top row reveals that the difference is consistently positive, supporting the intuitive notion that a stronger trade-off is determined in the case of ℙ_2, where both the epigenetic trait and environmental factor directly contribute to module proliferation. This observation is further corroborated by the bottom-row plots of Figure <ref>. In fact, we notice that the final percentage increment is six times higher for the epigenetic-driven duality (ℙ_1) with respect to the case in which there is a direct contribution of both epigenetic trait and environmental factor (ℙ_2). Considering, instead, the evolution for the center of mass, from the bottom-right plot of Figure <ref> we observe a perfect overlap of its velocity dynamics throughout the experiment. This confirms the fact that any observed differences can be solely attributed to the proliferative dynamics resulting from the distinct trade-offs under analysis, while the spatial dynamics are not affected.



 §.§ Experiment 3: impact of Snail expression and hypoxia on cancer cell migration

In this third experiment, we aim to assess the impact of Snail expression and exposure to hypoxia on cancer cell migration. Our goal is to validate our model by replicating experimental results that investigate the role of hypoxia in cell migration and determine whether motility can be triggered by various factors, including inhibition or up-regulation of Snail expression. 
To achieve this, we specify the parameter values and the environmental conditions such that they replicate different experimental scenarios, and we compare the resulting outcomes with empirical observations.

We consider the parameter g_s, accounting for Snail transcription and, starting from the reference value of g_s =1.5  (molecules·min^-1), used for the experiments in Section <ref> and <ref>, we define up-regulation and down-regulation of Snail expression by setting as g_s =2.1  (molecules·min^-1) and g_s =0.9 (molecules·min^-1), respectively. This corresponds to variations of 0.6 above and below the reference value. To ensure that the condition (<ref>) holds true in all the scenarios, we consider a slightly higher value for Snail degradation rate with respect to the previous experiments, i.e., we set  γ_s = 0.03  (min^-1). For the environmental conditions, by referring to <cit.> we consider levels of oxygenation compatible with normoxia (7%) and pathological hypoxia (1%) and we set the scaling factor O_2,0 such that, in our model, these conditions are represented by O_2=0.7 and O_2=0.1.

We aim at qualitatively replicating the experimental findings proposed in <cit.> and <cit.>. Specifically, in <cit.> the authors investigate human breast cancer cells (cell lineage MCF-7). In their experiment, they analyze fold change in tumor cell migration by migration assays using Transwell migration chambers. Precisely, cells are suspended in upper Transwell chambers in serum-free media and allow to migrate towards a serum gradient (10%) in the lower chamber for 6 hours. The experiment is repeated in normoxic conditions by transiently overexpressing and silencing Snail protein expression. Instead, in <cit.>, the authors employ a similar methodology with the human hepatocellular carcinoma (cell lines HepG2). They assess cell motility with the same migration assays comparing experiments conducted in normoxic and hypoxic conditions.

We replicate the chamber setup by considering our square domain Ω intersected by a vertical membrane aligned parallel to the x_2 axis and positioned at x_1=25 mm. A 1D section of the chamber setup is illustrated in Figure <ref>. 
 
Here, the left part of the domain (for x_1<25 mm) represents the upper chamber, where all cells are initially distributed following


    M_0()= 
        
          A_M                  x_1 < x_s
    
          A_Me^-(x_1-x_s)^2θ^2_M     x_1 ≥ x_s

where A_M = 0.9, x_s = 10 mm and θ_M = 3 mm. We consider the following linear expression for the oxygen distribution

    O() = O_min + (O_max - O_min)/50 x_1

with O_min=0.7 and O_max=1.0 in normoxic conditions, while O_min=0.1 and O_max=0.4 in hypoxic conditions. This choice establishes a fixed oxygen gradient along the chamber, which is consistent with the biological setting. We conduct five experiments, which are summarized in Table <ref>.


Under the aforementioned conditions, we allow cells to move in response to the environmental stimuli for a duration of T=6 hours. Subsequently, we measure the quantity of tumor mass that has passed through the membrane as

    Q(t) = ∫_Ω_BM(t,) d ,

where Ω_B = [25,50] × [0,50] mm^2 represents the bottom chamber.

We designate the results obtained in scenario (A) as the control case and we use them to normalize the outcomes of the other experiments. Figure <ref> collects all the results of the five tests. In the top row, we show the results related to the scenarios (A), (B), and (C) and we compare them with the data taken from <cit.>, while in the bottom row, we refer to scenarios (A), (D), and (E) and we compare the results with the data taken from <cit.>. Each row comprises two columns. The left column provides a map of the values taken by F(y^*, O_2) against various levels of oxygenation and Snail transcription, while the right column provides histograms comparing the results of in-vitro (black) and in-silico (red) experiments. 



As clearly shown in Figure <ref>, for both the breast cancer and the hepatocarcinoma cases, the model is able to effectively replicate the experimental data. Specifically, in the case of breast cancer, the in-silico results closely resemble those obtained in the in-vitro experiments for the scenario (B). In scenario (C), any discrepancy appears to be merely apparent, as the in-silico results fall within the error band of the in-vitro experiments, which notably exhibits a wider range of data. For hepatocarcinoma, there is a remarkably high correspondence between the in-vitro and in-silico data, and experiments (D) and (E) show a notable match. Furthermore, it is worth noticing how the level curves of F(g_s,O_2) (left column of Figure <ref>) provide insight into the experimental observations. Specifically, previous experimental works have noticed that the knockdown of Snail nullifies the motility advantage gained under hypoxic conditions (D) compared to normoxia (A), bringing the motility to a level comparable with the control case (as observable by comparing scenarios (E) and (A)). In our case, this can be observed a priori by looking at the location of the corresponding bullets on the level plot of F(g_s,O_2): (A) and (E) are located, in fact, almost on the same level curve. This implies that, given equal cell density (ensured at least initially with identical initial conditions) and a consistent oxygen gradient (maintained at the same value across the domain and for all oxygenation conditions in this experiment), the term governing cell speed assumes comparable values in both experiments (slightly lower for (E)). Therefore, comparable levels of fold changes in cell migration are expected as well as the slight discrepancy in the number of cells passing through the membrane between the two cases.




 §.§ Experiment 4: hypoxia-driven ring structure in tumor and Snail distributions

As last experiment, we refer to additional results shown in <cit.>, where the authors 
investigate the expression of Snail in non-invasive ductal carcinoma in situ (DCIS), an early breast cancer, considering a model system of hypoxia in-vivo. Considering a central necrotic area, their analysis of several DCIS samples revealed a typical pattern of HIF-1α expression, with increasing staining intensity approaching the areas of necrosis, and similar spatial distribution for the nuclear expression of Snail, gradually increased approaching the necrosis (see Figure 6 in <cit.>). In particular, the authors show that hypoxia induces Snail expression independently of other EMT markers. 

To qualitatively reproduce these observations, we simulate an initial tumor mass located in the center of the domain. We assume an oxygen distribution that decreases towards the center of the domain, leading to highly hypoxic (or necrotic) areas due to higher consumption in regions with higher cell density. The initial condition for cancer cells is given by (<ref>) setting A_M=1, 𝐱_M = [25,25] mm, and θ_M=5 mm, while for the fixed oxygen distribution we consider 1-O(𝐱), with O() defined by (<ref>) and A_O_2=0.8, 𝐱_O_2 = [25,25] mm, and θ_O_2=13 mm. The initial conditions for tumor cells and oxygen are illustrated in Figure <ref>.



Figure <ref> collects the results of this experiment at four time steps: t_0=0h, t_1=5h, t_2=20h, and t_3=35h. The first row depicts a 2D representation of the tumor mass, including a density map and contour lines highlighting the tumor's edge (defined by a density threshold corresponding to 10% of the carrying capacity). From system (<ref>), defined the average expression of the Snail protein in the cell population as

    M̅^y(t,) = M^y(t,)M(t,)= a_2a_1 F(y^*,O_2)  for ∈supp(M) ,

we illustrate its evolution in the second row of Figure <ref>. Finally, in the third row, the 1D section of the tumor mass density along the bisecting line at the four specified time steps is shown.


We observe that, initially, both anti-crowding and chemotactic stimuli point in the same direction, which tends toquickly move the cells away from the central hypoxic region, where cell density is high and oxygen concentration is low. Consequently, in this initial phase, cell migration is rapid, leading to a fast transformation of the peaked initial Gaussian for cell distribution into a smoother bubble profile (as shown at time t_1). As time progresses, the prevalence of chemotactic motion results in a depletion of cells from the central mass, gradually giving rise to a ring-like structure (times t_2 and t_3). During this phase, movement starts to slow down due to two main factors. Firstly, comparing the position of the ring with respect to the oxygen profile (shown in Figure <ref>), we observe a decrease in motility caused by both high levels of oxygenation (reducing the tactic sensitivity F(y^*,O_2)) and low oxygen gradients (reducing the chemotactic stimulus). Secondly, the slowdown is due to the anti-crowding mechanism, which, once the void forms at the center of the mass, would induce cells from the inner part of the ring to move towards the center, conflicting with the chemoattractant-driven movement. These observations are also consistent with the plots in the second row of Figure <ref>, where M̅^y(t,) is shown. They illustrate how, initially, the average expression of Snail is high, inducing rapid cell migration, and it increases approaching the central hypoxic core. Then, while the tumor mass moves outward, this expression decreases as cells reach more oxygenated areas, still maintaining higher values around the inner border of the ring. It is interesting how the model is able to qualitatively capture the two main dynamics shown in the data from <cit.>. In fact, the model reproduces both the experimentally observed ring shape of the tumor mass and the spatial distribution of the average Snail expression, mirroring the findings of the experimental study.




§ CONCLUSION


The migration of tumor cells in response to oxygen concentration gradients remains a critical area of investigation in cancer biology. While the role of hypoxia in promoting tumor aggressiveness and metastasis is well recognized, the exact mechanisms driving cell migration in response to oxygen levels are still an area of investigation and understanding these mechanisms may be crucial for developing effective therapeutic strategies.

In this study, we developed a novel mathematical model to investigate the interplay between hypoxia, molecular signaling pathways, and tumor cell migration. Specifically, we proposed a multi-scale model that naturally integrates single-cell behavior driven by Snail expression with macroscopic scale dynamics describing tumor migration in the tissue. Our approach employs tools and methods from the kinetic theory of active particles to construct a kinetic transport equation that describes the evolution of the tumor cell distribution based on detailed microscopic dynamics. By employing proper scaling arguments, we formally derived the equations for the statistical moments of the cell distribution. These capture cell density dynamics, influenced by limited non-linear diffusion and oxygen-mediated drift, and the evolution of the average Snail expression within the tumor population, which directly relates to tumor migratory capability. Overall, our model offers a detailed description of macroscopic tumor cell dynamics, considering the effect of microscopic Snail signaling pathways in the mechanisms of tumor response to hypoxia. This modeling approach represents a promising way to integrate molecular signaling pathways with cell migration dynamics.

We validated the model in different scenarios with biological relevance, focusing on the role of chemotactic dominated motion and anti-crowding effects, and analyzing the effect of Snail expression on cell migration and proliferation. We also showed the reliability of our approach by testing its ability to quantitatively replicate experimental results from two different studies published in the literature. We investigated the effect of hypoxia and Snail knockdown on the motility of cancer cells, comparing our results with those presented by <cit.> on human hepatocarcinomas. Moreover, we considered the findings of <cit.> and we replicate in-silico the results regarding the effect of Snail over-expression and Snail knockdown on the migration capability of human breast cancer cells in normoxic conditions. We also analyzed the spatial distribution of Snail expression within the tumor mass in response to hypoxia, showing how the model is able to reproduce the patterns observed experimentally in <cit.>. These results support the idea that our mathematical framework can offer new perspectives for interpreting experimental data and understanding the underlying biological mechanisms driving tumor migration.

Moving forward, it will be important to explore the implications of our findings in the context of clinical outcomes and therapeutic interventions. Particularly, our results highlight the importance of considering the dynamic regulation of Snail expression in response to hypoxia. This finding underscores the potential significance of developing strategies to target Snail as a therapeutic approach to control tumor cell migration and metastasis. Furthermore, incorporating heterogeneous and dynamic environmental factors, such as a non-stationary oxygen distribution with its consumption by tumor cells, could improve the predictive power of our model and enhance the quantitative fitting of the experimental data, ultimately leading to a better understanding of tumor invasion.

In summary, the proposed mathematical modeling approach is a novel and valuable tool to integrate detailed descriptions of microscopic cell dynamics with cell evolution at a macroscopic (tissue) level. In particular, the multi-scale modeling approach allows to properly derive the macroscopic terms driving cell evolution from a detailed description of the single-cell dynamics, instead of phenomenological stating them directly at the macroscopic level. Our findings offer interesting interpretations of the complex dynamics underlying tumor progression and motility, providing new perspectives for interpreting experimental data and understanding the biological mechanisms driving tumor development. This also paves the way for personalized medicine approaches tailored to individual tumor characteristics.




 §.§ Acknowledgement

This paper has been partially supported by the Italian Ministry of Education, Universities and Research, through the MIUR grant Dipartimento di Eccellenza 2018-2022, project E11G18000350001 (MC, GC, MD), by the National Group of Mathematical Physics (GNFM-INdAM) through the INdAM–GNFM Project (CUP E53C22001930001) "From kinetic to macroscopic models for tumor-immune system competition" (MC), and by the Modeling Nature Research Unit, Grant QUAL21-011 (MC), funded by Consejería de Universidad, Investigación e Innovación (Junta de Andalucía). MC also acknowledges support from City of Hope’s Global Scholar Program.
  



 §.§ Author contribution statement

Conceptualization: G.C., M.C., M.D. Methodology: G.C., M.C., M.D. Formal analysis: G.C., M.C., M.D. Resources: G.C., M.C., M.D. Software: G.C. Visualization: G.C., M.C. Writing-original draft preparation: G.C., M.C., M.D.  Writing—review and editing: G.C., M.C., M.D. Supervision: M.D. Project administration: M.D. Funding acquisition: M.C., M.D. All authors contributed to the article and approved the submitted version.




 §.§ Competing interests

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.



* 
abbrv ]