Federica Bubba

Benoît Perthame

Daniele Cerroni
email: daniele.cerroni@polimi.it

Pasquale Ciarletta
email: pasquale.ciarletta@polimi.it

Paolo Zunino
email: paolo.zunino@polimi.it

Benoit Perthame
email: benoit.perthame@sorbonne-universite.fr

A coupled 3D-1D multiscale Keller-Segel model of chemotaxis and its application to cancer invasion

Keywords: 2020 Mathematics Subject Classification. Primary: 92B05, 62P10; Secondary: 92D25 Chemotaxis, Keller-Segel, multiscale model, cancer

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1. Introduction. In many complex biological systems the macroscopic, visible dynamics is the result of exchanges and communications between features that act at different scales of time and space. Working like microscopic lens, mathematical modeling allows to capture the variety of phenomena that characterize multiscale systems. For this reason, in recent years the development, analysis and numerical resolution of multiscale models has become a fertile field of research in applied mathematics. Multiscale models have been successful in the description of a large variety of phenomena, in particular in the context of biomedical sciences. For instance, cellular growth and migration are regulated by the much smaller proteins inside the cells [START_REF] Preziosi | Mathematical models of the interaction of cells and cell aggregates with the extracellular matrix[END_REF]; metastatic invasion is influenced by the reciprocal communications between the cells and the surrounding tissue [START_REF] Trucu | Multiscale computational modelling and analysis of cancer invasion[END_REF]; in microcirculation and cardiac tissue perfusion, oxygen and nutrients can be exchanged between the tissue and the network of small capillary vessels [START_REF] Angelo | On the coupling of 1d and 3d di↵usionreaction equations. Application to tissue perfusion problems[END_REF][START_REF] Possenti | A computational model for microcirculation including Fahraeus-Lindqvist e↵ect, plasma skimming and fluid exchange with the tissue interstitium[END_REF]. In the latter field of application, the multiscale feature arises from exchanges between embedded geometries: small-scale networks of capillaries are embedded in a surrounding bulk structure (the cardiac tissue). Despite the growing performance of computers, the numerical resolution of realistic multiscale problems could remain very demanding from a computational viewpoint, as it would require very fine grids to solve the smallest scales. As a consequence, significant e↵orts have been put into the development and analysis of novel models and numerical methods that allow to reduce the computational complexity of multiscale problems. Here, we focus on the topological model reduction method [START_REF] Angelo | Multiscale modelling of metabolism and transport phenomena in living tissues[END_REF][START_REF] Laurino | Derivation and analysis of coupled PDEs on manifolds with high dimensionality gap arising from topological model reduction[END_REF][START_REF] Quarteroni | Geometric multiscale modeling of the cardiovascular system, between theory and practice[END_REF] that, based on spatial averaging techniques, enables to reduce a fully threedimensional problem involving embedded multiscale geometries into one that features coupled partial di↵erential equations with di↵erent dimensions, i.e., 3D-1D. The formulation and the analysis of coupled 3D-1D problems was performed by D'Angelo [START_REF] Angelo | Finite element approximation of elliptic problems with Dirac measure terms in weighted spaces: Applications to one-and three-dimensional coupled problems[END_REF][START_REF] Angelo | Multiscale modelling of metabolism and transport phenomena in living tissues[END_REF][START_REF] Angelo | On the coupling of 1d and 3d di↵usionreaction equations. Application to tissue perfusion problems[END_REF] and applied in a large variety of contexts, for instance microcirculation [START_REF] Possenti | A computational model for microcirculation including Fahraeus-Lindqvist e↵ect, plasma skimming and fluid exchange with the tissue interstitium[END_REF] and drug delivery for tumor treatments [START_REF] Cattaneo | A computational model of drug delivery through microcirculation to compare di↵erent tumor treatments[END_REF].

In this paper we present, derive and simulate a novel version of the Keller-Segel model of chemotaxis on embedded multiscale geometries, whereby a one-dimensional network is embedded in a three-dimensional bulk domain. In biology, chemotaxis refers to the ability of cells, bacteria or other organisms to orient their movement in response to gradients of chemical signals present in the environment. Chemotaxis is known to play a fundamental role in many biological processes, such as bacterial self-organization, wound healing, cancer growth and metastasis formation [START_REF] Kay | Changing directions in the study of chemotaxis[END_REF][START_REF] Roussos | Chemotaxis in cancer[END_REF][START_REF] Wadhams | Making sense of it all: bacterial chemotaxis[END_REF]. Mathematically, chemosensitive movement of biological populations can be described by the well-known Keller-Segel model [START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF][START_REF] Keller | Model for chemotaxis[END_REF], a system of coupled partial di↵erential equations describing the chemotaxis-driven self-organization of unicellular organisms orchestrated by the means of chemical signals emitted by the organisms themselves. The Keller-Segel system has attracted great attention thanks to its ability to faithfully reproduce chemotactic movement and numerous versions have been proposed. In fact, other than the original continuum population-based version, di↵erent scales of details can be considered, i.e., microscopic individual-based and mesoscopic descriptions [START_REF] Buttenschoen | A space-jump derivation for non-local models of cellcell adhesion and non-local chemotaxis[END_REF][START_REF] Hillen | The di↵usion limit of transport equations derived from velocity-jump processes[END_REF][START_REF] Othmer | The di↵usion limit of transport equations II: Chemotaxis equations[END_REF][START_REF] Painter | Volume-filling and quorum-sensing in models for chemosensitive movement[END_REF]. Recently, the Keller-Segel system has been investigated numerically and analytically on networkshaped domains [START_REF] Borsche | The scalar Keller-Segel model on networks[END_REF][START_REF] Bretti | On modeling maze solving ability of slime mold via a hyperbolic model of chemotaxis[END_REF][START_REF] Bretti | A hyperbolic model of chemotaxis on a network: a numerical study[END_REF] which intrinsically incorporate multidimensional e↵ects as they allow to enrich plain one-dimensional dynamics with multidimensional features. Furthermore, in [START_REF] Macdonald | A computational method for the coupled solution of reaction-di↵usion equations on evolving domains and manifolds: Application to a model of cell migration and chemotaxis[END_REF][START_REF] Mackenzie | Local modulation of chemoattractant concentrations by single cells: dissection using a bulk-surface computational model[END_REF] the authors consider coupled multidimensional bulk-surface (i.e., 3D-2D) versions of the system. We refer the interested reader to the surveys [START_REF] Bellomo | Towards a mathematical theory of Keller-Segel models of pattern formation in biological tissues[END_REF][START_REF] Othmer | The mathematical analysis of biological aggregation and dispersal: progress, problems and perspectives[END_REF] for a complete overview of the multiscale versions of the Keller-Segel model and its numerous variants. In this work, we apply the topological model reduction technique to derive a coupled 3D-1D formulation of the Keller-Segel model on embedded multiscale geometries, whereby a chemosensitive cellular population is constrained to move along a one-dimensional network embedded in a three-dimensional domain in which are expressed the chemical signals that influence cellular migration. Compared to previous results, the novelty of the present paper is to incorporate for the first time embedded multiscale features by considering one-dimensional network-shaped domains embedded in a three-dimensional bulk domain, therefore with a dimensional gap greater than 1.

The ultimate goal of our study stems from the willing to employ an embedded multiscale Keller-Segel-type system to model the invasion of potentially metastatic cancer cells in the extracellular matrix (ECM), the highly organized structure of fibers, proteins and other macromolecules that provides the sca↵old for cellular migration. In fact, biological evidences highlight that invasive cancer cells are able to remodel their microenvironment, locally degrade the ECM and exploit the network of fibers to spread in the surrounding tissue. Therefore, the intrinsic multiscale character of cancer invasion requires suitable mathematical models that take into account the dynamics acting on di↵erent scales. The problem of describing cellular migration and cancer invasion in the ECM has been tackled by several authors and the need to employ a multiscale approach is evident. Whilst some works have investigated individual-based models [START_REF] Drasdo | A single-cell-based model of tumor growth in vitro: monolayers and spheroids[END_REF][START_REF] Schlüter | Computational modeling of single-cell migration: The leading role of extracellular matrix fibers[END_REF] or continuum ones [START_REF] Bitsouni | Mathematical modelling of cancer invasion: The multiple roles of TGF-pathway on tumour proliferation and cell adhesion[END_REF][START_REF] Chaplain | Mathematical modelling of cancer cell invasion of tissue: dynamic heterogeneity[END_REF][START_REF] Chaplain | Mathematical modelling of cancer cell invasion of tissue: the role of the urokinase plasminogen activation system[END_REF][START_REF] Chaplain | Mathematical modelling of cancer invasion: the importance of cell-cell adhesion and cell-matrix adhesion[END_REF], multiscale and hybrid discrete-continuum approaches seem predominant; see, for instance, the works of Anderson et al. [START_REF] Anderson | A hybrid mathematical model of solid tumour invasion: The importance of cell adhesion[END_REF][START_REF] Anderson | Mathematical modelling of tumour invasion and metastasis[END_REF], Chaplain and collaborators [START_REF] Ramis-Conde | Mathematical modelling of cancer cell invasion of tissue[END_REF][START_REF] Sfakianakis | A hybrid multiscale model for cancer invasion of the extracellular matrix[END_REF][START_REF] Trucu | Multiscale computational modelling and analysis of cancer invasion[END_REF], Chauvière and Preziosi [START_REF] Chauvière | Mathematical framework to model migration of cell population in extracellular matrix[END_REF][START_REF] Chauvière | Modeling the motion of a cell population in the extracellular matrix[END_REF], Loy and Preziosi [START_REF] Loy | Kinetic models with non-local sensing determining cell polarization and speed according to independent cues[END_REF], Painter [START_REF] Painter | Modelling cell migration strategies in the extracellular matrix[END_REF], Rocha et al. [START_REF] Rocha | A hybrid three-scale model of tumor growth[END_REF] and references therein. Some of the mentioned works explore the importance of the orientation of matrix fibers on cellular infiltration and incorporate this factor in the model as a parameter or a variable. However, the multiscale 3D-1D approach we propose here paves the way towards more realistic models of cancer invasion, as it allows to investigate how the infiltration of cancer cells in the tissue is influenced by both the topography of the network of ECM fibers on which they migrate and the interplay between local and three-dimensional e↵ects.

Here, we consider the dynamics of a cellular population of density u(t, x) that emits a chemical signal of concentration c(t, x), called chemoattractant, described by the minimal Keller-Segel model:

8 > < > : @u @t r • (D 1 ru u v) = 0, t>0, x 2 D, @c @t D 2 c = ↵ u I D c, t > 0, x 2 ⌦ ⇢ R 3 , (1)
where ⌦ ⇢ R 3 is a regular bounded three-dimensional domain, D represents the domain on which the cellular dynamics takes place and v models the velocity field of the cellular population, proportional to chemical gradients. Starting from a fully three-dimensional Keller-Segel system, whereby D ⌘ ⌦ ⇢ R 3 and v = rc in System (1), in Section 2 we assume that the cellular population is restricted to move in a generalized three-dimensional cylinder ⌃ ⇢ R 3 entirely embedded in the bulk domain ⌦, therefore D ⌘ ⌃ ⇢ R 3 . We refer to this model as the restricted 3D Keller-Segel model. In Section 3, we let the cylinder ⌃ shrink to its one-dimensional centerline ⇤ and we apply the topological model reduction technique [START_REF] Angelo | Multiscale modelling of metabolism and transport phenomena in living tissues[END_REF][START_REF] Laurino | Derivation and analysis of coupled PDEs on manifolds with high dimensionality gap arising from topological model reduction[END_REF] to derive the coupled 3D-1D version of System (1), whereby D ⌘ ⇤ ⇢ R and v is the spatial average of the gradient of the chemoattractant concentration c. As a consequence, the reduced 3D-1D system features two partial di↵erential equations with di↵erent dimensions: a one-dimensional equation for u on the network ⇤ and a three-dimensional equation for c in the bulk domain ⌦.

We numerically solve the coupled 3D-1D Keller-Segel model in Section 4, where we also present the finite element scheme employed for the discretization. We consider three di↵erent one-dimensional manifolds: a segment in Subsection 4.2, a bifurcated network in Subsection 4.3 and more complex networks with several branches in Subsection 4.4.

In Sections 5 and 6 we apply the coupled 3D-1D framework to a model of cancer invasion in the tissue. More in details, in Section 5 we briefly summarize the background notions on cellular invasion in biological tissues and, inspired by the Keller-Segel-type model presented in [START_REF] Anderson | Mathematical modelling of tumour invasion and metastasis[END_REF], we derive a multiscale multidimensional model of cancer invasion. In this model, cells migrate along the branches of a one-dimensional network and their movement is influenced by the matrix density and the chemical signal produced by the cells, both three-dimensional quantities.

In Section 6, we perform numerical experiments of the model on several network topologies and we discuss the numerical results in light of the biological findings. Finally, in Section 7, we conclude the paper, present possible numerical advancements and provide insights on further researches.

2.

The restricted 3D Keller-Segel model. In a regular, bounded domain ⌦ ⇢ R 3 , let u(t, x) be the density of a cellular population at time t > 0 and at position x 2 ⌦ and let c(t, x) be the concentration of the chemical signal, called chemoattractant. We consider the fully three-dimensional Keller-Segel system 8 > > > < > > > :

@u @t r • (D 1 ru u rc) = 0, @c @t D 2 c = ↵ u c, t > 0, x 2 ⌦ ⇢ R 3 , (2)
with assigned nonnegative initial conditions u(0, x) = u 0 (x) 0, and c(0, x) = c 0 (x) 0, in ⌦. Moreover, we prescribe no-flux boundary conditions on the boundary @⌦ of ⌦ for the cellular density

( D 1 ru + u rc) • n ⌦ = 0, on @⌦,
and Dirichlet and mixed-Robin boundary conditions for the chemical concentration c on @⌦ DIR and @⌦ MIX respectively:

c = c ext , on @⌦ DIR , D 2 rc • n ⌦ = c (c c ? ) , on @⌦ MIX ,
where n ⌦ denotes the outward unit vector normal to @⌦ and where @⌦ DIR [@⌦ MIX = @⌦.

The first equation in System (2) describes the evolution in time of cells as the superposition of two e↵ects: random, unbiased cellular movement modeled according to the Fickian law of di↵usion with di↵usivity D 1 > 0, and chemotactic movement in the direction of greater concentration of the chemoattractant c, modeled as an advection term with velocity field rc and chemotactic sensitivity > 0. The second equation describes random di↵usive movement for c, with di↵usion rate D 2 > 0, depletion of the chemoattractant at a rate > 0 and production by the cells at a rate ↵ > 0.

In this section, we present the restricted version of the fully three-dimensional Keller-Segel model (2), whereby we suppose that cells are constrained to move in a generalized three-dimensional cylinder embedded in the bulk domain ⌦.

We shall first set up the geometrical configuration of the three and one-dimensional problems. This will allow us to apply the topological model reduction technique to derive the reduced 3D-1D Keller-Segel model in the subsequent section.

2.1. Geometry. Let (s) = [ 1 (s), 2 (s), 3 (s)], be a C 2 regular curve in the three-dimensional space parametrized by the arc-length s 2 (0, S), S > 0. We assume that 0 (s) = 1 such that the arc-length and the coordinate coincide. We call [⇠(s), ⌫(s), ⇣(s)] the Frenet's frame related to , namely the tangent, normal and binormal unit vectors, respectively. We define a generalized cylinder ⌃ as the volume obtained by sweeping a twodimensional surface (i.e., the cross-section) along . For each s 2 (0, S) we define the cross-section through the parametrization D(s) = [x(r, ✓), y(r, ✓)] : (0, R(s)) ⇥ (0, ⇥(s)) ! R 2 and we assume that D(s) is convex for any s. Let @D(s) := [x(R(s), ✓), y(R(s), ✓)] : (0, ⇥(s)) ! R 2 be the boundary of the cross-section D and assume it is a piecewise C 2 regular curve. Then, the generalized cylinder ⌃ can be defined as the volume Moreover, we denote by 0 = { (0) + D(0)} and S = { (S) + D(S)} respectively the top and bottom boundaries of ⌃. We assume that ⌃ is fully embedded in ⌦, i.e., the distance between @⌃ and @⌦ is strictly positive. Moreover, we let the cross-sections D(s) change size but not shape with s.

⌃ = { (s) + x(r, ✓) ⌫(s) + y(r, ✓) ⇣(s), r 2 (0, R(s)), ✓ 2 (0, ⇥(s)), s 2 (0, S)}, (3)
We can exploit the structure of the generalized cylinder to decompose the integrals on ⌃ as follows: for a su ciently regular function w defined on ⌃, we write

Z ⌃ w d! = Z ⇤ Z D(s) w d ds = Z ⇤ |D(s)| w ds, (4)
where w(s) : ⇤ ! R denotes the mean value of w on the cross-section D(s) and is defined as

w(s) := 1 |D(s)| Z D(s) w d .

2.2.

The Keller-Segel model with restricted cell motility. We now assume that the cells move preferentially along the generalized cylinder ⌃ defined in (3). This assumption is motivated by the observation that cells in tissues move by contact guidance, which means that they have the tendency to follow fibers orientation [START_REF] Painter | Modelling cell migration strategies in the extracellular matrix[END_REF]. To enforce this assumption in the original three-dimensional Keller-Segel model (2), we let ⇢ ⌃ be a regular weight function, entirely supported in ⌃, that represents the probability of a cell to travel away from the centerline ⇤ of the cylinder. Then, we restrict the motility of cells to the cylinder ⌃ by considering the following equation for the cellular density:

@u @t r • ( ⇢ ⌃ ) = 0, := D 1 ru u rc, in ⌦. (5)
We now want to prove that the assumption that the cellular flux is confined in the cylinder ⌃ leads to a restricted 3D Keller-Segel model, in which the cellular density u is nonvanishing only in ⌃. In other words, we prove that if the initial density u 0 is supported in ⌃ then

Z ⌦ u d! ⌘ Z ⌃ u d!, for all t > 0. (6)
In fact, since ⌦ is bounded, we can apply the theorem on the partitions of the unity to find two functions ⇧ ⌃ and ⇧ ? ⌃ that verify

• 0  ⇧ ⌃ , ⇧ ? ⌃  1, • ⇧ ⌃ + ⇧ ? ⌃ = 1, in ⌦. In particular, ⇧ ⌃ = 1 and ⇧ ? ⌃ = 0 in ⌃. Then Z ⌦ @u @t d! = Z ⌦ @u @t ⇧ ⌃ d! + Z ⌦ @u @t ⇧ ? ⌃ d!.
Let us consider first the second term. We use the definition of ⇧ ? ⌃ and Equation (5) to discover

Z ⌦ @u @t ⇧ ? ⌃ d! = Z ⌦/⌃ @u @t ⇧ ? ⌃ d! = Z ⌦/⌃ r • ( ⇢ ⌃ ) ⇧ ? ⌃ d! ⌘ 0, since ⇢ ⌃ is entirely supported in ⌃.
As for the first term, we write

Z ⌦ @u @t ⇧ ⌃ d! = Z ⌃ @u @t ⇧ ⌃ d! + Z ⌦/⌃ @u @t ⇧ ⌃ d!.
Using again Equation ( 5) and arguing as before, we infer R

⌦/⌃ @u @t ⇧ ⌃ d! ⌘ 0. Fi- nally, recalling that ⇧ ⌃ = 1 in ⌃, we find that Z ⌦ @u @t ⇧ ⌃ d! = Z ⌃ @u @t ⇧ ⌃ d! ⌘ Z ⌃ @u @t d!,
which leads to [START_REF] Borsche | The scalar Keller-Segel model on networks[END_REF] after integration in time and taking into account the assumption that

R ⌦ u 0 d! ⌘ R ⌃ u 0 d!.
The above arguments lead to consider the following restricted version of the threedimensional Keller-Segel system (2):

8 > < > : @u @t r • ( ⇢ ⌃ ) = 0, t>0, x 2 ⌃ ⇢ R 3 , @c @t D 2 c = ↵ u I ⌃ c, t > 0, x 2 ⌦ ⇢ R 3 , (7)
where I ⌃ denotes the indicator function of ⌃. On the boundary of ⌃, we prescribe no-flux boundary conditions for u: This might deserve a short motivation: assuming cells can move only when they are enough to the fiber?

( D 1 ru + u rc) • n ⌃ = 0, on @⌃,
where n ⌃ is the normal unit vector pointing outward the boundary of the cylinder ⌃.

In the remaining of the paper, we assume that the transversal diameter of ⌃ is much smaller then the diameter of ⌦, i.e., max s2(0,S)

diam(D(s)) ⌧ diam(⌦).
When the domain ⌃ consists of many cylinders, such as in the context of the biological applications mentioned in the introduction, the numerical resolution of the three-dimensional problem (7) could become very demanding. For this reason, we let ⌃ shrink to its one-dimensional centerline ⇤ and derive a multiscale multidimensional version of problem [START_REF] Bretti | On modeling maze solving ability of slime mold via a hyperbolic model of chemotaxis[END_REF].

3. The coupled 3D-1D Keller-Segel model. We are now in position to apply the topological model reduction to the restricted fully three-dimensional Keller-Segel system (5). We discuss separately the derivation of the equation for the cellular density u and for the chemoattractant c when the three-dimensional cylinder ⌃ shrinks to its centerline ⇤. We then consider the generalization of the model to a one-dimensional network with several branches.

3.1.

Topological model reduction of the problem on ⌃. Let P be an arbitrary portion of the generalized cylinder ⌃, delimited by the lateral surface P and the two cross-sections D(s 1 ) and D(s 2 ), with 0 < s 1 < s 2 < S.

Integrating the equation for u on P and applying the divergence theorem we get

Z P @u @t d! Z P r • ( ⇢ ⌃ ) d! = Z P @u @t d! Z @P ( ⇢ ⌃ ) • n ⌃ d .
Thanks to the integral decomposition formula in (4), we simplify the first term as

Z P @u @t d! = Z s2 s1 Z D(s) @u @t d ds = d dt Z s2 s1 |D(s)| u ds.
As for the second term, we decompose the boundary of P to get

Z @P ( ⇢ ⌃ ) • n ⌃ d = Z P ( ⇢ ⌃ ) • n ⌃ d + Z D(s1) ( ⇢ ⌃ ) • n ⌃ d + Z D(s2) ( ⇢ ⌃ ) • n ⌃ d .
The boundary condition for the cellular flux on @⌃, together with the fundamental theorem of the integral calculus gives

Z @P ( ⇢ ⌃ ) • n ⌃ d = Z D(s1) ( ⇢ ⌃ ) • n ⌃ d + Z D(s2) ( ⇢ ⌃ ) • n ⌃ d = Z D(s1) ✓ D 1 @u @s u rc • ◆ ⇢ ⌃ d + Z D(s2) ✓ D 1 @u @s u rc • ◆ ⇢ ⌃ d = Z s2 s1 d ds Z D(s) ✓ D 1 @u @s u rc • ◆ ⇢ ⌃ d ds.
Similarly to (4), we now define the weighted mean value -with weight ⇢ ⌃ -for a general regular function w as

w ⌃ (s) := 1 |D(s)| Z D(s) w ⇢ ⌃ d .
As a consequence, we get

Z @P ( ⇢ ⌃ ) • n ⌃ d = Z s2 s1 d ds " |D(s)| ✓ D 1 @u @s u rc • ◆ ⌃ # ds.
We now suppose that u is uniform on each cross-section, i.e., u(t, s, r, ✓) = U (t, s).

As a consequence, U = U and, recalling that all parameters are constants, we obtain ✓

D 1 @u @s u rc • ◆ ⌃ = D 1 ⇢ @U @s U v c • ,
where

⇢ := 1 |D(s)| R D(s) ⇢ ⌃ d and v c := (rc) ⌃ .
Thanks to the generality of P, the problem on ⌃ can be reduced to the problem on its centerline ⇤. Hence, we get

|D(s)| @U @t @ @s ✓ |D(s)| D 1 ⇢ @U @s |D(s)| U v c • ◆ = 0, on ⇤, (8)
complemented with the following conditions on the boundary of ⇤:

|D(s)| D 1 ⇢ @U @s + |D(s)| U v c • = 0, for s = 0, S.
Finally, we consider the variational formulation of the reduced one-dimensional problem [START_REF] Bretti | A hyperbolic model of chemotaxis on a network: a numerical study[END_REF]. After multiplying Equation ( 8) by a test function 2 H 1 (⇤) and integration by parts, we get

✓ @U @t , ◆ ⇤,|D(s)| + ✓ D 1 ⇢ @U @s U v c • , @ @s ◆ ⇤,|D(s)| = 0, (9)
where we used the short-hand notation

(f, g) ⇤,|D(s)| := Z ⇤ |D(s)| f (s) g(s) ds.
3.2. Topological model reduction of the problem on ⌦. Let ' 2 H 1 @⌦ DIR (⌦). Multiplying the second equation in ( 7) by ' and integrating by parts yields

Z ⌦ @c @t ' d! + Z ⌦ D 2 rcr' d!+ Z @⌦ MIX c c ' d + Z ⌦ c ' d! = Z @⌦ MIX c c ? ' d + Z ⌦ ↵ u ' I ⌃ d!.
The last integral can be simplified applying the topological model reduction. In order to achieve this, we write the solution u and the test function ' as the sum of their average on D(s) and some perturbation:

u = u + ũ, ' = ' + '
, where we suppose that ũ and ' have zero-average on each cross-section D(s). Thanks to this assumption, we get

Z ⌦ ↵ u ' I ⌃ d! = Z ⇤ ↵ Z D(s) (u + ũ) (' + ') d ds = Z ⇤ ↵ Z D(s) u v d ds + Z ⇤ ↵ Z D(s) ũ ' d ds.
Assuming that the product of the perturbations is small, namely that Z

D(s) ũ ' d ⇡ 0,
and recalling that u = U = U we get

Z ⌦ ↵ u ' I ⌃ d! = Z ⇤ ↵|D(s)| U ' ds.
As a consequence, the variational formulation for the problem on ⌦ becomes ✓ @c @t , '

◆ ⌦ + (D 2 rc, r') ⌦ + ( c c, ') @⌦ MIX + ( c, ') ⌦ = ( c c ? , ') @⌦ MIX + ↵ U, ' ⇤,|D(s)| , (10)
for all ' 2 H 1 @⌦ DIR (⌦), where as in the previous subsection we have introduced the short-hand notations

(f, g) := Z f g d!, with 2 ⌦, @⌦ MIX . ( 11
)
3.3. Extension of the model to a network. We now discuss the extension of the coupled 3D-1D Keller-Segel model ( 9)-( 10) to the case where ⇤ is a network. We define the network as a set of branches

i (s) = [ i 1 (s), i 2 (s), i 3 (s)], s 2 (0, S i ), i = 1, . . . , N, N 2
N, that are C 2 regular three-dimensional curves connected by the vertexes j 2 J . We call K j the set of the indexes i of the branches i (s) that are connected to the vertex j 2 J . We decompose K j into two subsets: K j , the set of indexes of the branches originating in j, and K + j , the set of indexes of the branches that end in the vertex j. Finally, we define B as the set of indexes of dead-end branches, that can similarly be split up into the subsets B and B + . Generalizing the geometrical setting presented in Section 2, we can define ⌃ as the union of the generalized cylinders obtained by sweeping the two-dimensional cross-sections D(s), s 2 (0, S i ), along i . The centerlines of ⌃ are now the one-dimensional sets

⇤ i = { i (s), s 2 (0, S i )}, i = 1, . . . , N,
of lengths S i . The topological model reduction technique presented for the single-branch case can also be applied to the network setting branch by branch. For i = 1, . . . , N, we set

U i := U ⇤i ,
and we formulate the coupled 3D-1D Keller-Segel model on the network ⇤ = [ N i=1 ⇤ i as follows:

8 > > < > > : |D(s)| @U i @t @ @s ✓ |D(s)| D 1 ⇢ @U i @s |D(s)| U i v c • i ◆ = 0, on ⇤ i , @c @t D 2 c = ↵ U ⇤ c, in ⌦. (12a) (12b)
In addition to the boundary conditions on @⇤ and @⌦

|D(s)| ✓ D 1 ⇢ @U i @s + U i v c • i ◆ s=0 = 0, 8i 2 B , |D(s)| ✓ D 1 ⇢ @U i @s + U i v c • i ◆ s=Si = 0, 8i 2 B + , c = c ext , on @⌦ DIR , D 2 rc • n ⌦ = c (c c ? ) , on @⌦ MIX , (13)
we also prescribe Neumann-Kirchho↵ conditions to ensure flux conservation at each vertex of the network that is not dead-end: for j 2 J we impose

X i2K j |D(s)| ✓ D 1 ⇢ @U i @s + U i v c • i ◆ s=0 = X i2K + j |D(s)| ✓ D 1 ⇢ @U i @s + U i v c • i ◆ s=Si . (14)
3.4. Variational formulation. To write the variational formulation of the coupled 3D-1D Keller-Segel model, we define the Sobolev space H 1 (⇤) as the space of continuous functions V : ⇤ ! R, such that their restrictions V i to ⇤ i is in [START_REF] Laurino | Derivation and analysis of coupled PDEs on manifolds with high dimensionality gap arising from topological model reduction[END_REF]. Thus, for U, 2 H 1 (⇤) we multiply the density equation in ( 12) by and integrate by parts to get

H 1 (⇤ i ),
N X i=1 " ✓ @U i @t , i ◆ ⇤i,|D(s)| + ✓ D 1 ⇢ @U i @s U i v c • i , @ i @s ◆ ⇤i,|D(s)| # = N X i=1 |D(s)| ✓ D 1 ⇢ @U i @s + U i v c • i ◆ i s=0 N X i=1 |D(s)| ✓ D 1 ⇢ @U i @s + U i v c • i ◆ i s=Si .
Using the fact that is continuous on ⇤, we can factorize the test functions. Moreover, the no-flux boundary conditions on the dead-end branches B (13) and the flux-continuity Neumann-Kirchho↵ conditions [START_REF] Chauvière | Mathematical framework to model migration of cell population in extracellular matrix[END_REF] at the junctions make the righthand side vanish. Thus, the variational formulation of the problem ( 12) is

✓ @U @t , ◆ ⇤,|D(s)| + ✓ D 1 ⇢ @U @s U v c • , @ @s ◆ ⇤,|D(s)| = 0, (15a) coupled with ✓ @c @t , ' ◆ ⌦ + (D 2 rc, r') ⌦ + ( c c, ') @⌦ MIX + ( c, ') ⌦ = ( c c ? , ') @⌦ MIX + ↵ U, ' ⇤,|D(s)| , (15b)
for all 2 H 1 (⇤) and ' 2 H 1 @⌦ DIR (⌦).

4. Numerical simulations.

4.1. Numerical scheme. Let us consider a quasi-uniform partition T h ⌦ of ⌦ and an admissible partition T h ⇤ of ⇤ and let V ⌦ h and V ⇤ h be two finite element spaces of piecewise linear functions defined on T h ⌦ and T h ⇤ respectively.

We consider the following semi-implicit spatio-temporal discretization of (15a)-(15b): for k = 0, 1, . . . , N T , find

U k+1 h 2 V ⇤ h and c k+1 h 2 V ⌦ h such that 8 > > > > > > > > > > < > > > > > > > > > > : 1 t U k+1 h , h ⇤,|D(s)| + D 1 ⇢ dU h ds k+1 U k+1 h v k c,h • , d h ds ! ⇤,|D(s)| = 1 t U k h , h ⇤,|D(s)| , 1 t c k+1 h , ' h ⌦ + D 2 rc k+1 h , r' h ⌦ + c k+1 h , ' h ⌦ + c c k+1 h , ' h @⌦ MIX = 1 t c k h , ' h ⌦ + ( c c ? h , ' h ) @⌦ MIX + ↵ U k h , ' h ⇤,|D(s)| , ( 16
) for all h 2 V ⇤ h and ' h 2 V ⌦ h .
In System ( 16), we have introduced the discretized velocity v k c,h :=

⇣ rc k h ⌘ ⌃
and, for k = 0, we suppose U 0 h and c 0 h given nonnegative initial conditions, projections of u 0 and c 0 on the respective discrete spaces.

Setting N ⇤ h := dimV ⇤ h and N ⌦ h := dimV ⌦ h , we let { h l } N ⇤ h l=1 and {' h l } N ⌦ h l=1
be two independent sets of finite elements basis respectively of V ⇤ h and V ⌦ h . We consider the finite element approximations of the solutions

U k h = N ⇤ h X l=1 U k l h l (s), s 2 ⇤, c k h := N ⌦ h X l=1 c k l ' h l (x), x 2 ⌦,
and collect their degrees of freedom in the vectors of unknowns

U k := [U k 1 , • • • , U k N ⇤ h ] T
, and

c k := [c k 1 , • • • , c k N ⌦ h ] T .
Then, we define the finite element matrices for the problem on ⇤

[M ⇤ ] n,l := h l , h n ⇤,|D(s)| , n , l = 1, . . . , N ⇤ h , ⇥ A k ⇤ ⇤ n,l := ✓ D 1 ⇢ d h l ds h l v k c,h • , d h n ds ◆ ⇤,|D(s)| , n ,l= 1, . . . , N ⇤ h , for the problem on ⌦ [M ⌦ ] n,l := ' h l , ' h n ⌦ , n , l = 1, . . . , N ⌦ h , [A ⌦ ] n,l := D 2 r' h l , r' h n ⌦ + ' h l , ' h n ⌦ + c ' h l , ' h n @⌦ MIX , n,l= 1, . . . , N ⌦ h , [F ⌦ ] n,l := c c ? h , ' h n @⌦ MIX , n = 1, . . . , N ⌦ h ,
and the matrix deriving from the coupling of the 3D with the 1D problem

[B] n,l := ⇣ ↵ h l , ' h n ⌘ ⇤,|D(s)| , l = 1, . . . , N ⇤ h , n = 1, . . . , N ⌦ h .
At this point, we can rewrite the variational formulation of the coupled 3D-1D Keller-Segel model in the following matrix form: for k = 0, 1, 2, . . . , we search for

U k+1 2 R N ⇤ h and c k+1 2 R N ⌦ h such that 2 6 4 1 t M ⇤ + A k ⇤ • • 1 t M ⌦ + A ⌦ 3 7 5  U k+1 c k+1 = 2 6 4 1 t U k 1 t c k + F ⌦ + BU k 3 7 5 .
Remark 1. For the sake of simplicity, in the following simulations we let the crosssection be constant along the network, in such a way that |D(s)| ⌘ ⇡ R 2 , where R is the radius of the three-dimensional cylinder of the non reduced problem, that we assume to be the same for all the branches of the network. Moreover, we choose

⇢ ⌃ ⌘ 1 in ⌃.
The simulations are obtained using the software GetFEM [START_REF] Renard | GetFEM: Automated FE modeling of multiphysics problems based on a generic weak form language[END_REF]. The code has been developed starting from the project MANWorks, available on Github (https: //github.com/stefano-brambilla-853558/MANworks).

4.2.

A network with a single branch. As a first numerical test, we solve the coupled 3D-1D Keller-Segel system ( 12)-( 13)-( 14) on a single-branch network ⇤ completely embedded in a unit cubic domain ⌦ = (0, 1) 3 (cf. Figure 1). In this simplified context, the segment -of endpoints [x 0 , y 0 , z 0 ] = [0.1, 0.5, 0.5] and [x 1 , y 1 , z 1 ] = [0.9, 0.5, 0.5] -represents the one-dimensional manifold along which the cellular population can move. We assume that there is an initial nodule of cells on the branch concentrated around the leftmost endpoint x = 0.1, of density

u 0 = e |x 0.1| 2 0.01 , x 2 [0.1, 0.9], (17)
and that the initial concentration of chemoattractant is increasing with x, namely

c 0 (x, y, z) = 1 e |x| 2 0.1 . ( 18
)
We simulate the evolution of System [START_REF] Chaplain | Mathematical modelling of cancer cell invasion of tissue: the role of the urokinase plasminogen activation system[END_REF] with the following parameters Parameter values in Figure 2 Eq. (12a) D 1 = 0.01 = 5 Eq. (12b) D 2 = 1 ↵ = 20 = 0 and with no-flux boundary conditions on the tissue boundary, i.e., with c ⌘ 0 on @⌦. We let the di↵usion rate of the chemoattractant be much bigger that the one of the cellular population, coherently to what has been observed experimentally (see, e.g., Murray [START_REF] Murray | An introduction[END_REF]). Moreover, we let the chemotactic sensitivity be much stronger than cellular di↵usivity D 1 in such a way that the chemotactic aggregation prevails over cellular di↵usion in the long time. For this simulation, we employ a uniform mesh of 21 points in each direction for the discrete three-dimensional domain, while for the one-dimensional segment we employ a uniform mesh of 200 points. The time step is fixed to t = 0.01 and the radius of the original three-dimensional segment is R = 0.05. Figure 2 shows the evolution in time of the solutions to [START_REF] Chaplain | Mathematical modelling of cancer cell invasion of tissue: the role of the urokinase plasminogen activation system[END_REF] at four subsequent times. The nodule of cells initially elongates due to di↵usion and chemoattraction. As expected, cells travel along the segment in the direction of higher concentrations of the chemoattractant. This process is relatively fast. In fact, all the cells aggregate near the rightmost endpoint x = 0.9 of the segment in about ten iterations (t ⇡ 0.1). Once they reach this position in space, the cellular aggregate concentrates into a single peak (cf. Figure2(c)). At this point, the chemical concentration reaches a (nearly) uniform state due to the e↵ects of di↵usion. In Figure 3 we display the cellular density and the chemical concentration along the segment (thus for x 2 [0.1, 0.9] and y = z = 0.5) at four di↵erent times. In particular, in Figures 3(c)-(d) one can remark the e↵ects of the production of chemoattractant by the cellular population in the equation for c: the quasi-constant chemical concentration is perturbed around the cellular aggregate. These e↵ects are small, a consequence of the choice of the radius R and of the chemoattractant production rate ↵. We also remark that, once it has reached the endpoint x = 0.9, the cellular aggregate slowly travels back, in the direction of the endpoint x = 0.1, see Figures 3(c)-(d). This e↵ect is a consequence of the asymmetry of the initial problem. The radius of the non reduced three-dimensional network is constant everywhere in the three branches of the bifurcation, R = 0.05. As in the previous section, we center the initial density of cells at the leftmost endpoint of the bifurcation, setting for x 2 [0.1, 0.9] u 0 = e (x 0.1) 2 0.01 , and create the chemotactic gradient in the direction of increasing values of x by letting c 0 (x, y, z) = 1 e |x| 2 0.1 . We discretize ⌦ with a uniform mesh of 21 points in each direction while we choose a uniform discretization of 100 points for each branch of the network. The time step is set to t = 0.01.

We display the time evolution of the solutions in Figure 4, at t = 0, 0.04, 0.16 and 1. The initial nodule of cells concentrated in x = 0.1 follows the chemical gradient in the direction of increasing values of x. At the junction, located at x = 0.6, the nodule splits into two equal aggregates. Once the cells have entered the two branches as two equal aggregates, since the chemical gradient is weaker, the di↵usive e↵ects prevail over chemotactic forces, and the cells mainly spread in the two rightmost branches. We remark that the state showed in Figure 4 at the final time t = 1 is not a steady state: the cells are slowly di↵using and we expect a uniform final state.

4.4.

Networks with multiple branches. Here, we investigate the dynamics of the solutions of the coupled 3D-1D Keller-Segel system on more complex networks composed of multiple branches. We consider first a network of 8 equal branches of radius R = 0.05 that lays on the plane (x, y) and is symmetric with respect to the lines y = 0.5 and x = 0.5. We simulate the dynamics of System [START_REF] Chaplain | Mathematical modelling of cancer cell invasion of tissue: the role of the urokinase plasminogen activation system[END_REF] with the parameter values listed in the following table :  Parameter values in Figure 5 Eq. (12a) D 1 = 0.01 = 5 Eq. (12b) D 2 = 1 ↵ = 1 = 0.01 and prescribing no-flux boundary conditions for the chemical concentration on @⌦. We employ a uniform discretization of 31 points in each direction for ⌦ and a uniform mesh with 21 points in each branch of the network. The time step is t = 0.01.

Figure 5 shows the evolution in time of the corresponding solutions of the multiscale Keller-Segel model, at t = 0, 0.04, and 0.3. The small initial nodules of cells are initially given by ( 17), concentrated near x = 0.1 on each of the three left branches of the network. The chemical gradient, of initial density given by ( 18), creates an external bias that transports symmetrically the cells on the opposite branches of the network. The chemoattractant di↵uses and uniformly degrades. Let us notice that the invasion of the cells in the tissue is not significantly perturbed by the presence of the two central branches, transversal with respect to the direction of the chemical gradient (cf. Figure 5(b)). This is a remarkable di↵erence with the behavior of the solutions of the system we will simulate in Section 6.

In the second case, inspired by [START_REF] Bretti | On modeling maze solving ability of slime mold via a hyperbolic model of chemotaxis[END_REF], we consider a slightly more complex network, displayed in Figure 6. Here, we study how a small aggregate of cells, placed at the leftmost branch of the network, of initial density u 0 = e (x 0.1) 2 +(y 0.9) 2 0.025 , evolves through a series of multifurcations, guided by a chemical gradient of initial concentration c 0 (x, y, z) = x + (1 y). Figure 6 displays the time evolution at six subsequent times of the cellular density and the corresponding chemoattractant concentration, solutions to the 3D-1D Keller-Segel model whereby the parameters have been set as follows:

Parameter values in Figure 6 Eq. (12a) D 1 = 0.2 = 17 Eq. (12b) D 2 = 0.001 ↵ = 1 = 0.0 Since we aim at observing the e↵ects of the chemical attraction on the spatial distribution of cells over time, we choose a much smaller di↵usivity for the chemoattractant, compared to the one of the cellular population, and a relatively strong chemotactic sensitivity. As displayed in Figure 6, the cells spread through the network by filling most importantly the shortest path, on which the chemical gradients are stronger. In the long time, all cells aggregate in the rightmost branch of the network, in a single-peak stationary state. 5. A multiscale model of metastatic invasion. Building on the results presented so far, we propose a multiscale model of migration and invasion of aggressive cancer cells in the extracellular matrix. In Subsection 5.1 we provide the biological background on the structure of the extracellular matrix as well as on its relevance in cellular migration. We highlight that the ability of aggressive cancerous cells to invade the matrix by exploiting its network of fibers is one of the hallmarks of metastatic progression. In Subsection 5.2, we consider the model of cancer invasion in the extracellular matrix proposed by Anderson et al. in [START_REF] Anderson | Mathematical modelling of tumour invasion and metastasis[END_REF] and generalize its formulation to the present context of embedded multiscale geometries. We then derive a coupled 3D-1D model whereby the one-dimensional dynamics of cells is influenced by two three-dimensional quantities: the density of the extracellular matrix and the concentration of matrix degradative enzymes produced by the cells themselves.

5.1. Biological background. The extracellular matrix (ECM) is a complex meshwork of proteins and other macromolecules that provides the physical architecture for cellular adhesion and migration. The ECM is composed of two main components: fibrous proteins (such as collagens, elastins, fibronectins and laminins) and proteoglycans that, like an hydrated a gel, fills the extracellular interstitial space [START_REF] Frantz | The extracellular matrix at a glance[END_REF].

Collagens are the most present proteins in the ECM: they determine the physical three-dimensional structure of the matrix and support cell adhesion, migration and chemotaxis, thanks to their ability to assemble into supramolecular structures, such as fibrils and networks [START_REF] Eble | The extracellular matrix in tumor progression and metastasis[END_REF][START_REF] Frantz | The extracellular matrix at a glance[END_REF]. The topography and sti↵ness of the extracellular matrix plays a fundamental role in cell migration and metastasis invasion. In fact, in a healthy tissue the ECM is a relaxed network of fibers and gel-like molecules that dynamically resists to cellular stresses. However, cancerous cells are able to secrete ECM-modifying enzymes that alter the matrix structure. In this altered situations, the collagens fibers and networks undergo a reorientation and linearization that increases the sti↵ness of the extracellular matrix [START_REF] Eble | The extracellular matrix in tumor progression and metastasis[END_REF][START_REF] Egeblad | Dynamic interplay between the collagen sca↵old and tumor evolution[END_REF]. This abnormal restructuring of the collagen fibers facilitates cancerous invasion into the interstitial tissue and promotes metastases formation. In fact, the fibers work as "highways" for the invading cells by facilitating cancer dissemination in the tissue [START_REF] Egeblad | Dynamic interplay between the collagen sca↵old and tumor evolution[END_REF]. On the other hand, the ECM meshwork can constitute a barrier to cellular migration. A hallmark of malignancy in tumors is the ability of cells to breach the dense extramatricial barriers. In fact, in order to further invading or invade? the surrounding tissue, malignant tumoral cells secret matrix-degrading enzymes (MDEs) and create tunnels in the surrounding interstitial tissue to make their way to the blood vessels [START_REF] Egeblad | Dynamic interplay between the collagen sca↵old and tumor evolution[END_REF][START_REF] Painter | Modelling cell migration strategies in the extracellular matrix[END_REF].

The architecture and topology of the extracellular matrix plays a relevant role on the cellular migration and the first phases of metastasis formation. In the following of this section, we present, study and simulate a coupled 3D-1D version of a Keller-Segel-type model on networks of the metastatic invasion and migration in the ECM. 5.2. The model. In order to model malignant cancer cell invasion in the ECM and to study the interactions between the cellular population and the surrounding tissue, inspired by the model presented in [START_REF] Anderson | Mathematical modelling of tumour invasion and metastasis[END_REF], we let u(t, x) be the density of cancer cells, m(t, x) be the density of the extracellular matrix and c(t, x) be the concentration of matrix-degrading enzymes (MDEs). In a regular bounded three-dimensional domain ⌦ representing the tissue, the model proposed in [START_REF] Anderson | Mathematical modelling of tumour invasion and metastasis[END_REF]

is 8 > > > > > > < > > > > > > : @u @t r • (D 1 c ru u rm) = 0, @c @t D 2 c = ↵ u c, @m @t =  c m, ( 19a
) (19b) (19c)
for x 2 ⌦ and t > 0, with nonnegative initial conditions u(0, x) = u 0 (x) 0, c(0, x) = c 0 (x) 0, and m(0, x) = m 0 (x) 0, in ⌦ and equipped with no-flux boundary conditions for u and c on @⌦

( D 1 c ru + u rm) • n ⌦ = 0, D 2 rc • n ⌦ = 0.
In System [START_REF] Drasdo | A single-cell-based model of tumor growth in vitro: monolayers and spheroids[END_REF], cancer cells move by random, di↵usive movement with di↵usion rate D 1 c, with D 1 > 0, thus proportional to the concentration of MDEs: cellular motility increases in regions of the domain where the concentration of chemical factors is higher, modeling chemokinesis. Moreover, the random chemokinetic movement is biased in the direction of higher density of the ECM: cells move up matrix gradients in order to invade the tissue. Since the macromolecules and proteins that constitute the ECM are non-di↵usible, this cellular response is called haptotaxis. The MDEs concentration di↵uses with constant rate D 2 > 0, undergoes a natural decay with rate > 0 and is produced by cells at a rate ↵ > 0. Finally, we let the ECM be degraded by the MDEs at a rate  c, with  > 0 constant. It is straightforward to adapt the topological model order reduction technique detailed in Section 3 to the three-dimensional model [START_REF] Drasdo | A single-cell-based model of tumor growth in vitro: monolayers and spheroids[END_REF]. In this context, we represent the extracellular matrix composed of a three-dimensional bulk region ⌦ containing proteoglycans and unstructured, degraded collagens in which a one-dimensional network ⇤ of linearized collagen fibrils is embedded. We then let m(t, x) represent the density of the three-dimensional structure of the ECM and suppose that the cellular motility is restricted to the network ⇤. We make the assumption that the MDEs can di↵use in the entire domain ⌦ and only degrade the interstitial region of the extracellular matrix. As in Section 3, we consider a one-dimensional network composed of N branches ⇤ i and call U i := U ⇤i the restrictions of the density u to the i-th branch. The multiscale 3D-1D version of ( 19) can be obtained following the methodology described in Section 3:

8 > > > > > > > < > > > > > > > : |D(s)| @U i @t @ @s ✓ |D(s)| D 1 c ⌃ @U i @s |D(s)| U i v m • i ◆ = 0, on ⇤ i , @c @t D 2 c = ↵ U ⇤i c, in ⌦, @m @t =  c m, in ⌦, (20a) (20b) (20c)
for t > 0, with v m := (rm) ⌃ . We assign no-flux and mixed boundary conditions on @⇤ and @⌦, respectively

|D(s)| ✓ D 1 c ⌃ @U i @s + U i v m • i ◆ s=0 = 0, 8i 2 B , |D(s)| ✓ D 1 c ⌃ @U i @s + U i v m • i ◆ s=Si = 0, 8i 2 B + , c = c ext , on @⌦ DIR , D 2 rc • n ⌦ = c (c c ? ) , on @⌦ MIX , (21)
and Neumann-Kirchho↵ conditions at each vertex of the network that is not deadend: for j 2 J X

i2K j |D(s)| ✓ D 1 c ⌃ @U i @s + U i v m • i ◆ s=0 = X i2K + j |D(s)| ✓ D 1 c ⌃ @U i @s + U i v m • i ◆ s=Si . (22)
We highlight that the ability of cells to degrade locally the matrix through the action of the MDEs allows them to create small tunnels in the ECM density that facilitate their migration along the fibers. This process, known in biology as contact guidance, is particularly relevant in the coupled 3D-1D formulation of the model, as we will see in Section 6.

On the parameter values. As mentioned in [START_REF] Anderson | Mathematical modelling of tumour invasion and metastasis[END_REF] and references therein, from experimental evidences the cell di↵usivity is estimated to be of order D 1 ⇠ 10 10 cm 2 s 1 and the haptotactic strength ⇠ 10 3 cm 2 s 1 , while the di↵usion rate of the chemicals is in the range D 2 ⇠ 10 8 10 10 cm 2 s 1 .

In the simulations we will perform, the cross-section of the network is always a circular section of constant radius R and we will always consider the case where R ⌧ 1. Therefore, even when

D 1 ⇠ D 2 , the ratio ⇡ R 2 D 1 D 2
⌧ 1 will remain coherent with biological measurements. However, regarding the advection strength -for which biological evidences indicate a value of at least 13 orders of magnitude greater that D 1 -computational limits constrain to lower values of which, in any case, will always satisfy c ⌃ D 1 .

6. Numerical simulations.

6.1. Numerical scheme. Let us consider the same partitions T h ⌦ of ⌦ and T h ⇤ of ⇤ and the same finite element spaces V ⌦ h and V ⇤ h presented in Section 4. For the multiscale model [START_REF] Drasdo | A single-cell-based model of tumor growth in vitro: monolayers and spheroids[END_REF] we introduce the following semi-implicit spatio-temporal discretization: for k = 0, 1, . . . , N T , find

U k+1 h 2 V ⇤ h and c k+1 h , m k+1 k 2 V ⌦ h such that 8 > > > > > > > > > < > > > > > > > > > : 1 t U k+1 h , h ⇤,|D(s)| + D 1 ⇢ c k h⌃ dU h ds k+1 U k+1 h v k m,h • , d h ds ! ⇤,|D(s)| = 1 t U k h , h ⇤,|D(s)| , 1 t c k+1 h , ' h ⌦ + D 2 rc k+1 h , r' h ⌦ + c k+1 h , ' h ⌦ + c c k+1 h , ' h @⌦ MIX = 1 t c k h , ' h ⌦ + ( c c ? h , ' h ) @⌦ MIX + ↵ U k h , ' h ⇤,|D(s)| , 1 t m k+1 h , ' h ⌦ +  c k h m k+1 h , ' h ⌦ = 1 t c k h , ' h ⌦ , (23) for
all h 2 V ⇤ h and ' h 2 V ⌦ h , with v k m,h := ⇣ rm k h ⌘ ⌃ . We let { h l } N ⇤ h l=1 and {' h l } N ⌦
h l=1 be two independent finite elements basis of V ⇤ h and V ⌦ h , respectively, and write

U k h = N ⇤ h X l=1 U k l h l (s), s 2 ⇤, c k h := N ⌦ h X l=1 c k l ' h l (x), m k h := N ⌦ h X l=1 m k l ' h l (x), x 2 ⌦.
Defining the vectors containing the degrees of freedom of the unknowns

U k := [U k 1 , • • • , U k N ⇤ h ] T , c k := [c k 1 , • • • , c k N ⌦ h ] T , and
m k := [m k 1 , • • • , m k N ⌦ h
] T , we can rewrite the variational formulation of the coupled 3D-1D Keller-Segel model in the following matrix form: for k = 0, 1, 2, . . . , we search for

U k+1 2 R N ⇤ h and c k+1 , m k+1 2 R N ⌦ h such that 2 4 1 t M ⇤ + A k ⇤ • • • 1 t M ⌦ + A ⌦ • • • 1 t M ⌦ + D k ⌦ 3 5 2 4 U k+1 c k+1 m k+1 3 5 = 2 4 1 t U k 1 t c k + F ⌦ + BU k 1 t m k 3 5 .
where

⇥ A k ⇤ ⇤ n,l := ✓ D 1 ⇢ c k h⌃ d h l ds h l v k m,h • , d h n ds ◆ ⇤,|D(s)| , n,l= 1, . . . , N ⇤ h , ⇥ D k ⌦ ⇤ n,l :=  c k h ' h l , ' h n ⌦ , n , l = 1, . . . , N ⌦ h ,
and the remaining matrices have been defined in Section 4.

6.2.

A network with a single branch. We solve numerically the multiscale model ( 20) on a single-branch network ⇤ embedded in a unit cubic domain ⌦, as showed in Figure 1 in Section 4. The original three-dimensional cylinder has a constant radius R = 0.05. We simulate System [START_REF] Eble | The extracellular matrix in tumor progression and metastasis[END_REF] with the following values of the parameters: Parameter values in Figure 9 Eq. (20a) D 1 = 0.01 = 0.5 Eq. (20b) D 2 = 0.001 ↵ = 10 = 0 Eq. (20c)  = 0.5 and we assign no-flux boundary conditions on @⌦ for the MDEs concentration c. Following [START_REF] Anderson | Mathematical modelling of tumour invasion and metastasis[END_REF], we set the initial densities and concentrations to

u 0 (x) = e (x 0.1) 2 0.01 , on ⇤, c 0 (x, y, z) = 0.5, in ⌦, m 0 (x, y, z) = 1 1 2 e |x| 2 0.01 , in ⌦.
We employ a uniform grid of 21 points in each direction for ⌦ and a uniform discretization of 200 points for the single-branch network ⇤. The time-step is fixed, t = 0.01. Figures 8 and9 display the time evolution of the cellular density U , of the ECM density m and of the MDEs concentration c. In Figure 8, U and m are plotted on the one-dimensional segment (thus for x 2 [0.1, 0.9]) while in Figure 9 U is plotted on the one-dimensional branch and m and c in the corresponding three-dimensional domain at di↵erent times of the evolution.

The initial nodule of cells invades and degrades the matrix in the form of an asymmetric pulse that, as expected, travels in the direction of higher density of the ECM, thus for increasing values of x. The matrix slowly degrades under the local e↵ects of the degrading enzymes. At about the time t ⇡ 0.24 (cf. Figure 8(c)), a small aggregate of cells detaches from the main peak, guided by the non-monotonicity in the ECM density. Compared to the corresponding simulation performed in Section 4, this is a novel e↵ect, consequence of the local matrix degradation. In the long time, most of the cells reach the rightmost endpoint of the segment, at x = 0.9, while the extracellular matrix degrades up to complete depletion. Notice that the way in which the ECM degrades depends not only on the rate of production of MDEs, the parameter ↵, and on the degradation rate , but also on the di↵usion rate of the MDEs. In fact, for bigger values of D 2 , the matrix degrades more uniformly in the three-dimensional domain ⌦. For this simulation, D 2 is chosen small compared to D 1 in order to enhance the locality of matrix degradation. This e↵ect can be seen clearly in Figure 8(f ). Notice that D 2 ⇡ R 2 D 1 remains true, as discussed previously. We also remark that in a single-branch network, the ability of invasion of the cellular population is not significantly a↵ected by the choice of D 2 . Influences of network topology and multidimensionality on invasion. Next, we investigate how the topology of the network or the three-dimensional spatial effects modify the invasive abilities of an initial nodule of cells along the single-branch network. Considering the same initial conditions, values of the parameters of the model and numerical setting, we first slightly modify the topology of the network by adding a small branch of length 0.1 to the main segment. This alteration creates a junction at the node x = 0.5. The corresponding solutions are displayed in Figure 10(a)-(b) at two di↵erent times. The initial dynamics remains the same until the aggregate of cells reaches the junction: there, some cells leave the main aggregate and fill the small transversal branch (cf. Figure 10(a)). As a result, the aggregate reaches the rightmost endpoint of the main segment later: the invasion is slowed down. In fact, at time t = 1.2, the peak has not yet reached the endpoint x = 0.9, see Figure 10(b). This behavior is easily explained: since some of the cells have escaped from the main segment to fill the smallest one, the term that accounts for production of MDEs in the equation for c becomes weaker. As a consequence, the degradation of the matrix is less e↵ective.

As a second case, we investigate how the multidimensionality of the model affects the dynamics of the solutions. To this end, we add a small branch of length 0.4, parallel to the main segment and at a constant distance of 0.2 from it. The corresponding dynamics is displayed in Figures 10(c)-(d). This alteration speeds up the invasion of the initial nodule of cells: the rightmost endpoint is reached at the time t ⇡ 1.12, about 8 iterations before compared to the simulation displayed in Figure 9. In fact, in this case, the cells on the smallest branch contribute to the degradation of the matrix, aiding the invasion of the tissue. In this simulation, the discretization of the three-dimensional domain is uniform, with 21 points in each direction, while we employ 50 points for each branch of the network. The time step is t = 0.01.

The two initial aggregates of cells travel in the direction of increasing values of x, guided by larger densities of the matrix. Initially, the invasion is quick and the transversal branches (i.e., the branches that are orthogonal to the direction of increasing gradients of m) do not seem to interfere with the invasion. However, when the cells enter the part of the domain where the gradient of m is nearly zero for the choice of m 0 (cf. Figure 11) and for the e↵ects of di↵usion and boundary conditions, at x = 0.6 (see Figure 12(d)), the cellular invasion slows down and stops The invasion is also influenced by the choice of the radius. In fact, for values of R bigger than 0.025 on the one hand the degradation of the matrix would be more e↵ective, leading to more invasive dynamics but, on the other hand, the loss of mass in the transversal branches would be higher, leading to smaller aggregates for which the invasion would be more di cult. Influence of fibers orientation on cellular invasion. Invasive cancer cells are able to remodel the meshwork of collagen fibers that constitute the extracellular matrix. In particular, from biological evidences it comes to light that sti↵er matrices (whereby the collagen fibers have been linearized by the enzymes secreted by the cancerous cells and are mostly oriented in the direction of blood vessels) facilitate cellular invasion. We now employ the multidimensional Keller-Segel-type system (20) to qualitatively reproduce di↵erences in the invasive abilities of cancerous cells depending on the topology of the network, and in particular on the orientation of its branches. Let us consider the following set of parameters in System [START_REF] Eble | The extracellular matrix in tumor progression and metastasis[END_REF] Parameter values in Figures 14 Eq. (20a) D 1 = 0.01 = 0.5 Eq. (20b) D 2 = 0.01 ↵ = 10 = 0 Eq. (20c)  = 0.5 with R = 0.025 and no-flux boundary conditions on the boundary @⌦ of the unit cubic domain ⌦. We study the evolution in time of three small nodules of cells on three equal parallel branches ⇤ 1 , ⇤ 2 , ⇤ 3 that lie on the plane z = 0.5, for x 2 [0.1, 0.9] and of three di↵erent ordinates y 1 = 0.48, y 2 = 0.5, y 3 = 0.52, respectively. The initial cellular density on these branches is

u 0 = e (x 0.1) 2 0.01
, on ⇤ 1 , ⇤ 2 , ⇤ 3 , and we take the following initial distributions for the matrix m 0 and the chemical c 0 :

m 0 (x, y, z) = 1 1 2 e x 2
0.01 , c 0 (x, y, z) ⌘ 0.1, in ⌦.

We then consider two slightly di↵erent configurations of the network, let us call them N 1 and N 2 . In both cases, to the three parallel branches ⇤ 1 , ⇤ 2 and ⇤ 3 , we add 12 branches of di↵erent lengths and orientations on which the initial cellular density is constant and equal to u 0 = 0.2. The di↵erence between N 1 and N 2 lies in the fact that in N 2 most of the branches are oriented in the direction of increasing vales of x, thus of increasing values of ECM density. The total length of the branches of N 1 and N 2 remains unaltered. Therefore, the two configurations have the same cellular mass. We compare the evolution in time of the three initial aggregates on ⇤ 1 , ⇤ 2 and ⇤ 3 in Figures 14 and15 in the three-dimensional domain and on the line x 2 [0.1, 0.9], respectively, for three subsequent times. The di↵erence between the two networks N 1 and N 2 is small (only the orientation of 4 branches was modified), and so is the di↵erence between the corresponding solutions. Nevertheless, it is evident from Figure 15 how for the network N 2 , in which more branches are oriented in the direction of increasing values of x, the invasion of cells is deeper. This can be noticed most importantly in the final times t = 0.92, 1.46. This di↵erence is greater for the cells on the branch ⇤ 3 , since a bigger aggregate has advanced on the branch at the final time of the simulation (cf. Figures 15(e) and (f ), blue line). A di↵erence can be noticed also for the cells on the branch ⇤ 1 (cf. Figures 15(e 7. Discussion and future work. Multiscale mathematical models and the related numerical methods have been intensively studied in the context of biomedical applications. In most cases, the di↵erent scales are either treated via separate models, e.g., by deriving PDE-based models from individual-based ones [START_REF] Drasdo | A single-cell-based model of tumor growth in vitro: monolayers and spheroids[END_REF][START_REF] Painter | Volume-filling and quorum-sensing in models for chemosensitive movement[END_REF][START_REF] Hillen | The di↵usion limit of transport equations derived from velocity-jump processes[END_REF][START_REF] Othmer | The di↵usion limit of transport equations II: Chemotaxis equations[END_REF], or by integrating di↵erent approaches, i.e., individual-based and PDE-based, as in [START_REF] Ramis-Conde | Mathematical modelling of cancer cell invasion of tissue[END_REF][START_REF] Rocha | A hybrid three-scale model of tumor growth[END_REF][START_REF] Sfakianakis | A hybrid multiscale model for cancer invasion of the extracellular matrix[END_REF][START_REF] Trucu | Multiscale computational modelling and analysis of cancer invasion[END_REF].

In this work we have presented a novel advancement in the study of chemotaxis models by representing the multiscale features as embedded geometries of di↵erent dimensions, namely one-dimensional network-shaped manifolds embedded in threedimensional domains. More in details, we have generalized the topological model reduction technique [START_REF] Angelo | Multiscale modelling of metabolism and transport phenomena in living tissues[END_REF][START_REF] Laurino | Derivation and analysis of coupled PDEs on manifolds with high dimensionality gap arising from topological model reduction[END_REF] to chemotaxis models, deriving a coupled 3D-1D version of the Keller-Segel system whereby the chemosensitive population moves along the one-dimensional network following three-dimensional signals. Moreover, we have investigated the dynamics of the solutions via several numerical experiments involving networks of di↵erent structures: a single-segment, a bifurcation and networks with several branches and multifurcations. The multiscale Keller-Segel model presented here, to the best of our knowledge, improves previous results. In fact, on the one hand, previous works study chemotaxis models where both the cellular population and the chemical concentration evolve on one-dimensional networks [START_REF] Bretti | A hyperbolic model of chemotaxis on a network: a numerical study[END_REF][START_REF] Bretti | On modeling maze solving ability of slime mold via a hyperbolic model of chemotaxis[END_REF][START_REF] Borsche | The scalar Keller-Segel model on networks[END_REF] and, hence, multidimensionality is not integrated in the model. On the other hand, Mackenzie and collaborators [START_REF] Mackenzie | Local modulation of chemoattractant concentrations by single cells: dissection using a bulk-surface computational model[END_REF][START_REF] Macdonald | A computational method for the coupled solution of reaction-di↵usion equations on evolving domains and manifolds: Application to a model of cell migration and chemotaxis[END_REF] have considered coupled multidimensional bulk-surface models of chemotaxis. However, in their work the multidimensionality gap is only 1 (i.e., 3D for the bulk and 2D for the surface) and the intriguing dynamics arising from network-shaped domains are not included. Additionally, we have proposed the coupled 3D-1D Keller-Segel system as a model of invasion of potentially metastatic cancer cells in the extracellular matrix, represented as a three-dimensional domain containing a structured network of fibers. In this context, cells move along the one-dimensional network towards regions where the ECM is more dense (haptotaxis) and emit a three-dimensional di↵usible enzyme that degrades the matrix. We investigate how both the topography of the fibers and the multidimensional coupling influence the invasive potentiality of cells. In particular, we find that a network in which the fibers are oriented in the same direction facilitates cellular invasion, confirming biological evidences. We believe that our results open the way towards more realistic models of cancer invasion in tissues. However, numerous improvements can be considered. From a modeling point of view, we aim at implementing much more realistic networks of ECM fibers, such as those presented in [START_REF] Schlüter | Computational modeling of single-cell migration: The leading role of extracellular matrix fibers[END_REF][START_REF] Chauvière | Mathematical framework to model migration of cell population in extracellular matrix[END_REF], to improve results on the crucial role of fibers orientation in the invasion of cells. From a numerical point of view, more robust numerical schemes should be considered. In fact, the finite element method employed in this paper does not preserve exactly steady states and the mass of solutions. Moreover, it does not allow to consider parameter values in realistic ranges. Limiting ourselves to finite elements methods, schemes based on discontinuous Galerkin discretization [START_REF] Yekaterina | New interior penalty discontinuous Galerkin methods for the Keller-Segel Chemotaxis model[END_REF][START_REF] Epshteyn | Discontinuous Galerkin methods for the chemotaxis and haptotaxis models[END_REF][START_REF] Li | Local discontinuous Galerkin method for the Keller-Segel chemotaxis model[END_REF] or on scalar auxiliary variable techniques [START_REF] Poulain | Scalar auxiliary variable finite element scheme for the parabolicparabolic Keller-Segel model: positivity preserving and energy stability[END_REF] could be more reliable. It is also worth noticing that the modeling framework presented in this work can be adapted to take into account several other applications whereby multiscale features naturally arise, for instance wound healing, tissue-engineering, and neurodegenerative diseases, such as Alzheimer.

The coupled 3D-1D systems proposed here might also deserve deeper theoretical understanding. For instance energy dissipation is a remarkable of the standard Keller-Segel system, which lead to elaborate dissipative numerical schemes, see [START_REF] Almeida | Energy and implicit discretization of the Fokker-Planck and Keller-Segel type equations[END_REF] and the references therein. This leads to understand how our reduction procedure behaves with respect to energy properties. Finite time blow-up is another theoretical feature of the Keller-Segel system which has been widely studied in 2D and 3D, see the recent paper [START_REF] Souplet | Blow-up profiles for the parabolicelliptic Keller-Segel system in dimensions n 3[END_REF], and rises questions in the coupled formalism.

  of centerline ⇤ = { (s), s 2 (0, S)} ,and lateral boundary = { (s) + x(R(s), ✓) ⌫(s) + y(R(s), ✓) ⇣(s), ✓ 2 (0, ⇥(s)), s 2 (0, S)}.

Figure 1 .

 1 Figure 1. Outline of the unit cubic domain ⌦ (in black) and the embedded single-branch network ⇤ (red line).

2 Figure 2 .

 22 Figure 2. Evolution of the cellular density U on a single-branch network and of the corresponding three-dimensional chemoattractant concentration c at four di↵erent times, t = 0.0, 0.08, 0.26 and 1.2.

2 Figure 3 .

 23 Figure 3. Initial conditions and subsequent states of the evolution in time of the cellular density U and of the concentration of the chemoattractant c along the single-branch network, thus for s 2 [0, 0.8] (or, equivalently, for x 2 [0.1, 0.9]).

4. 3 .

 3 A bifurcated network. We consider now the case where the cellular population moves along a bifurcated network symmetric with respect to the line y = 0.5 and embedded in a unit cubic domain ⌦ as in the previous section. Considering no-flux boundary conditions for c ( c ⌘ 0) on @⌦, we set in[START_REF] Chaplain | Mathematical modelling of cancer cell invasion of tissue: the role of the urokinase plasminogen activation system[END_REF] the parameter values summarized in the following table:Parameter values in Figure 4 Eq. (12a) D 1 = 0.01 = 10 Eq. (12b) D 2 = 1 ↵ = 1 = 0.01 (a) t = 0 (b) t = 0.04 (c) t = 0.16 (d) t = 1

Figure 4 .

 4 Figure 4. Evolution in time of the solution of the 3D-1D Keller-Segel model, on a one-dimensional bifurcated network embedded in a unit three-dimensional domain, at four di↵erent times of the simulation: t = 0, 0.04, 0.16, 1.

3 Figure 5 .

 35 Figure 5. (a)-(c) Time evolution of three initial aggregates of cells on a symmetric network and the corresponding chemical concentration at three di↵erent times, t = 0, 0.04, and 0.3. (d) Internal visualization of the chemoattractant concentration (on a di↵erent scale of values).

Figure 6 .

 6 Figure 6. Evolution of the cellular density and of the chemical concentration, solutions to the multiscale Keller-Segel model, at six di↵erent times.

Figure 7 .

 7 Figure 7. Schematic representation of the metastatic cascade, adapted from [20]. Malignant tumor cells escape from the primary tumor site by breaching the extracellular barriers (such as the basement membrane) (1), take advantage of the collagen fibrils and networks to invade the tissue (2) and enter the blood circulation (3).

2 Figure 8 .

 28 Figure 8. Time evolution of the cellular density U (blue line) andthe corresponding ECM density m (red line) on the one-dimensional single-branch network, for x 2 [0.1, 0.9] and at four di↵erent times t = 0, 0.12, 0.24, 1.2.

2 Figure 9 .

 29 Figure 9. Solutions of the multiscale Keller-Segel-type model (20) in the single-branch test case at three di↵erent times. The evolution of the cellular density U on the one-dimensional segment is displayed together with the density of the extracellular matrix m (left panels) and with the MDEs concentration c (right panels).

12 Figure 10 .

 1210 Figure 10. The time evolution of the cellular density U and the density of the extracellular matrix m for two di↵erent networks. The invasion is slowed down ((a)-(b)) or sped up ((c)-(d)) compared to the simulation in Figure 9.

6. 3 . 1 ,

 31 Networks with multiple branches. Let us consider the unit cubic threedimensional domain for the tissue ⌦ and study the one-dimensional dynamics on more complex networks.We first consider the evolution of two small nodules of cells, of initial density u 0 (x) = e (x 0.1) 2 0.01, for x < 0.2, an extracellular matrix of initial density (cf. Figure11)m 0 (x, y, z)and a small concentration of MDEs constant in all the domain, namely c 0 ⌘ 0.1 in ⌦.

Figure 11 .

 11 Figure 11. Initial density m0 of the extracellular matrix for x 2 [0, 1]. System (20) is solved by setting Parameter values in Figure 12 Eq. (20a) D 1 = 0.01 = 0.5 Eq. (20b) D 2 = 0.25 ↵ = 10 = 0 Eq. (20c)  = 0.5 with no-flux boundary conditions on @⌦ and a constant radius R = 0.025. The corresponding evolution in time of the solutions is displayed in Figure 12.In this simulation, the discretization of the three-dimensional domain is uniform, with 21 points in each direction, while we employ 50 points for each branch of the network. The time step is t = 0.01.The two initial aggregates of cells travel in the direction of increasing values of x, guided by larger densities of the matrix. Initially, the invasion is quick and the transversal branches (i.e., the branches that are orthogonal to the direction of increasing gradients of m) do not seem to interfere with the invasion. However, when the cells enter the part of the domain where the gradient of m is nearly zero for the choice of m 0 (cf. Figure11) and for the e↵ects of di↵usion and boundary conditions, at x = 0.6 (see Figure12(d)), the cellular invasion slows down and stops

6 Figure 12 .

 612 Figure 12. Evolution in time of the cellular density U and the corresponding density of the ECM m, solutions of the multiscale Keller-Segeltype model (20), at six subsequent times of the simulations.

6 Figure 13 .

 613 Figure 13. Evolution in time of the cellular density U and the corresponding density of the ECM m, solutions of the multiscale Keller-Segeltype model[START_REF] Eble | The extracellular matrix in tumor progression and metastasis[END_REF]. The invasion is more e↵ective without the presence of transversal branches, compared to the solutions displayed in Figure12.

  ) and (f ), red line): most of the cells have reached the rightmost endpoint of the branch at time t = 1.46.

46 Figure 14 .

 4614 Figure 14. Comparison of the evolution in time of the cellular density U for the network N1 (left panels) and N2 (right panels) at three subsequent times of the simulations.

46 Figure 15 .

 4615 Figure 15. Comparison of the evolution in time of the cellular density U on the three parallel branches ⇤1 (red line), ⇤2 (green line) and ⇤3 (blue line) for the networks N1 (left panels) and N2 (right panels).

Acknowledgement. We are grateful to the family of Federica Bubba for their support in this work. Her memory as brilliant scholar and generous human being will remain indelible among all of us.

This work is partly supported by MIUR grant PRIN grant 2017KL4EF3. B.P. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 740623).