G Bretti

R Natalini
email: roberto.natalini@cnr.it&g.bretti@cnr.it

M Ribot
email: ribot@unice.fr

A HYPERBOLIC MODEL OF CHEMOTAXIS ON A NETWORK: A NUMERICAL STUDY

Keywords: 1991 Mathematics Subject Classification. 65M06, 35L50, 92B05, 92C17, 92C42 hyperbolic system on network, initial-boundary value problem, transmission conditions, asymptotic behavior, finite difference schemes, chemotaxis

In this paper we deal with a semilinear hyperbolic chemotaxis model in one space dimension evolving on a network, with suitable transmission conditions at nodes. This framework is motivated by tissue-engineering scaffolds used for improving wound healing. We introduce a numerical scheme, which guarantees global mass densities conservation. Moreover our scheme is able to yield a correct approximation of the effects of the source term at equilibrium. Several numerical tests are presented to show the behavior of solutions and to discuss the stability and the accuracy of our approximation.

Introduction

The movement of bacteria, cells or other microorganisms under the effect of a chemical stimulus, represented by a chemoattractant, has been widely studied in mathematics in the last two decades, see [START_REF] Horstmann | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences[END_REF][START_REF] Murray | Mathematical biology. I. An introduction[END_REF][START_REF] Perthame | Transport equations in biology[END_REF], and numerous models involving partial differential equations have been proposed. The basic unknowns in these chemotactic models are the density of individuals and the concentrations of some chemical attractants. One of the most considered models is the Patlak-Keller-Segel system [START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF], where the evolution of the density of cells is described by a parabolic equation, and the concentration of a chemoattractant is generally given by a parabolic or elliptic equation, depending on the different regimes to be described and on authors' choices. The behavior of this system is quite well known now: in the one-dimensional case, the solution is always global in time, while in two and more dimensions the solutions exist globally in time or blow up according to the size of the initial data. However, a drawback of this model is that the diffusion leads to a fast dissipation or an explosive behavior, and prevents us to observe intermediate organized structures, like aggregation patterns.

By contrast, models based on hyperbolic/kinetic equations for the evolution of the density of individuals, are characterized by a finite speed of propagation and have registered a growing consideration in the last few years [START_REF] Gamba | Percolation, morphogenesis, and Burgers dynamics in blood vessels formation[END_REF][START_REF] Perthame | Transport equations in biology[END_REF][START_REF] Filbet | Derivation of hyperbolic models for chemosensitive movement[END_REF][START_REF] Hillen | Hyperbolic models for chemosensitive movement[END_REF][START_REF] Dolak | Cattaneo models for chemosensitive movement. Numerical solution and pattern formation[END_REF]. In such models, the population is divided in compartments depending on the velocity of propagation of individuals, giving raise to kinetic type equations, either with continuous or discrete velocities.

Here we consider an hyperbolic-parabolic system which arises as a simple model for chemotaxis:

(1.1)

   u t + v x = 0, v t + λ 2 u x = φ x u -v, φ t -D φ xx = au -bφ.
Such kind of models were originally considered in [START_REF] Segel | A theoretical study of receptor mechanisms in bacterial chemotaxis[END_REF], and later reconsidered in [START_REF] Greenberg | Stability results for a diffusion equation with functional drift approximating a chemotaxis model[END_REF]. They are based on an adaptation to the chemotactic case of the so-called hyperbolic heat or Cattaneo or telegraph equation, adding a source term accounting for the chemotactic motion in the equation for the flux. The function u is the density of cells in the considered medium, v is their averaged flux and φ denotes the density of chemoattractant. The individuals move at a constant speed λ, changing their direction along the axis during the time. The positive constant D is the diffusion coefficient of the chemoattractant; the positive coefficients a and b, are respectively its production and degradation rates.

These equations are expected to behave asymptotically as the corresponding parabolic equations, but displaying a different and richer transitory regime, and this is what is known to happen at least without the chemotactic term. Analytically, these models have been studied in [START_REF] Hillen | Hyperbolic models for chemotaxis in 1-D. Nonlinear Anal[END_REF][START_REF] Hillen | Existence of weak solutions for a hyperbolic model of chemosensitive movement[END_REF] and more recently in [START_REF] Guarguaglini | Stability of constant states and qualitative behavior of solutions to a one dimensional hyper-bolic model of chemotaxis[END_REF], where the analytical features were almost completely worked out, at least around constant equilibrium states, where it is proved that, at least for the Cauchy problem, the solutions of the hyperbolic and parabolic models are close for large times.

The novelty of this paper is to consider this one dimensional model on a network. More precisely, we consider system in the form (1.1) on each arc of the network, and so we have to consider one set of solutions (u, v, φ) for each arc. Functions on different arcs are coupled using suitable transmission conditions on each node of the network. Conservation laws or wave equations on networks have already been studied, for example in [START_REF] Garavello | Traffic flow on networks[END_REF] for traffic flows or in [START_REF] Dáger | Wave propagation, observation and control in 1-d flexible multi-structures[END_REF][START_REF] Valein | Stabilization of the wave equation on 1-D networks[END_REF] for flexible strings distributed along a planar graph. However, here we consider different types of transmission conditions, which impose the continuity of the fluxes rather than the continuity of the densities. Therefore, in this article, a particular care will be given to the proper setting and the numerical approximation of the transmission conditions at nodes, both for the hyperbolic and the parabolic parts of (1.1). In particular, some conditions have to be imposed on the approximation of the boundary conditions, in order to ensure the conservation of the total mass of the system. Let us also mention that a first analytical study of system (1.1) on a network, coupled through transmission conditions of this type, is carried out in [START_REF] Guaraldo | Some analytical results for hyperbolic chemotaxis model on networks[END_REF].

The study of this system is motivated by the tissue-engineering research concerning the movement of fibroblasts on artificial scaffolds [START_REF] Harley | Microarchitecture of Three-Dimensional Scaffolds Inuences Cell Migration Behavior via Junction Interactions[END_REF][START_REF] Mandal | Cell proliferation and migration in silk broin 3D scaffolds[END_REF][START_REF] Spadaccio | A G-CSF functionalized PLLA scaffold for wound repair: an in vitro preliminary study[END_REF], during the process of dermal wound healing. The natural process of healing of a damaged tissue occurs through a first phase in which fibroblasts, the stem cells to be in charge of to the reparation of dermal tissue, create a new extracellular matrix, essentially made by collagen, and, driven by chemotaxis, migrate to fill the wound. In recent years, tissue-engineering research has developed some new techniques, which aim at accelerating the wound healing. Actually, cellular migration on an injured body zone is stimulated and improved by using artificial scaffolds constituted by a network of crossed polymeric threads inserted within the wound, which mimic the extracellular matrix. The fibroblasts's reparation action is accelerated, since they already have a support and also they are constrained to move along the network, with less degrees of freedom, and it is believed that this approach could be effective in minimizing scarring [START_REF] Spadaccio | A G-CSF functionalized PLLA scaffold for wound repair: an in vitro preliminary study[END_REF]. Therefore, our simple model of chemotaxis on a network, which can be obtained by reducing the kinetic model of cell movement on a 3D extracellular matrix proposed in [START_REF] Chauviere | Mathematical framework to model migration of cell population in extracellular matrix, Chapter in "Cell mechanics. From single scale-based models to multiscale modeling[END_REF] to the case of a network, see [START_REF] Guaraldo | Some analytical results for hyperbolic chemotaxis model on networks[END_REF], is a good candidate for reproducing this configuration: the arcs of the network stand for the fibers of the scaffold and the transport equations give the evolution of the density of fibroblasts on each fiber. However, in this paper, we only address the numerical aspects of this problem. More direct applications of this framework to the real biomedical problem will be explored in future research. As reported in [START_REF] Harley | Microarchitecture of Three-Dimensional Scaffolds Inuences Cell Migration Behavior via Junction Interactions[END_REF], the present understanding of the critical biochemical and biophysical parameters that affect cell motility in three-dimensional environments is quite limited. Nevertheless, it has been observed that junction interactions affect local directional persistence as well as cell speed at and away from the junctions, so providing a new mechanism to control cell motility by using the extracellular microstructure. Therefore, mathematical modeling and simulations could play a crucial role in providing a better understanding of these phenomena, and an optimization tool for designing improved scaffolds.

The main focus of this paper is on the construction of an effective numerical scheme for computing the solutions to this problem, which is not an easy task, even for the case of a single arc. In that case, non constant highly concentrated stationary solutions are expected and schemes which are able to capture these large gradients in an accurate way are needed. The main problem is to balance correctly the source term with the differential part, in order to avoid an incorrect approximation of the density flux at equilibrium, as first observed in [START_REF] Guarguaglini | Stability of constant states and qualitative behavior of solutions to a one dimensional hyper-bolic model of chemotaxis[END_REF]. Asymptotic High Order schemes (AHO) were introduced in [START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF], inspired by [START_REF] Aregba-Driollet | Asymptotic high-order schemes for 2 × 2 dissipative hyperbolic systems[END_REF], to deal with this kind of inaccuracies. These schemes are based on standard finite differences methods, modified by a suitable treatment of the source terms, and they take into account for the behavior of the solutions near non constant stationary states. An alternative approach, inspired by the well-balanced methods, has been proposed in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF][START_REF] Gosse | Well-balanced numerical approximations display asymptotic decay toward Maxwellian distributions for a model of chemotaxis in a bounded interval[END_REF], with similar results. However the methods in [START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF] seem easier to be generalized to the present framework.

Regarding the problem considered in this paper, the main difficulty is in the discretization of the transmission conditions at node, also enforcing global mass conservation at the discrete level. Therefore, in Section 2 we explain some analytical properties of problem (1.1), with a particular emphasis on boundary and transmission conditions. Section 3 is devoted to the numerical approximation of the problem based on a AHO scheme with a suitable discretization of the transmission and boundary conditions ensuring the mass conservation. In the present paper, we have chosen to consider only the second order version of the scheme, which is enough for our purposes, but it is easy to adapt also the third order schemes proposed in [START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF]. Remark that here, unlike the single interval case, we are forced, for any given time step, to fix the space step on each arc using relation (3.14) introduced in Section 3, to obtain consistency on the boundary. Numerical tests (not shown) confirm the necessity of this supplementary constraint.

Finally, in Section 4, we report some numerical experiments, to show the behavior and the stability of our scheme. A special attention is given to the stability of the scheme near nodes and the correct behavior of the approximation for large times and near asymptotic states. It has to be mentioned that during this research we observed, in contrast with what happens for the diffusive models, the appearance of blow-up phenomena even for data of relative moderated size. Even if, up to now, there are no rigorous results, which can help to decide if these singular events are really occurring, or they are just a numerical artifact, our close investigation in Subsection 4.3 gives a strong indication towards the first alternative.

Analytical background

Let us define a network or a connected graph G = (N , A), as composed of two finite sets, a set of P nodes (or vertices) N and a set of N arcs (or edges) A, such that an arc connects a pair of nodes. Since arcs are bidirectional the graph is non-oriented, but we need to fix an artificial orientation in order to fix a sign to the velocities. The network is therefore composed of "oriented" arcs and there are two different types of intervals at a node p ∈ N : incoming ones -the set of these intervals is denoted by I p -and outgoing ones -whose set is denoted by O p . For example, on the network depicted in Figure 1, 1, 2 ∈ I and 3, 4 ∈ O. We will also denote in the following by I out and O out the set of the arcs incoming or outgoing from the outer boundaries. The N arcs of the network are parametrized as intervals a i = [0, L i ], i = 1, . . . , N , and for an incoming arc, L i is the abscissa of the node, whereas it is 0 for an outgoing arc. 2.1. Evolution equations for the problem. We consider system (1.1) on each arc and rewrite it in diagonal variables for its hyperbolic part by setting (2.1)

u ± = 1 2 u ± v λ .
Here u + and u -are the Riemann invariants of the system and u + (resp. u -) denotes the density of cells following the orientation of the arc (resp. the density of cells going in the opposite direction). This transformation is inverted by u = u + + u - and v = λ(u +u -), and yields:

(2.2)

           u + t + λu + x = 1 2λ (φ x -λ)u + + (φ x + λ)u -, u - t -λu - x = - 1 2λ (φ x -λ)u + + (φ x + λ)u -, φ t -Dφ xx = a(u + + u -) -bφ.
We can also denote by T ± = 1 2λ (φ x ∓ λ) the turning rates (namely the probabilities of cells to change direction) and a(u + + u -)bφ represents the production and degradation of the chemoattractant. We assume that all the cells are moving along an arc with the same velocity λ (in modulus), which may depend however on the characteristics of the arc. For the moment, we omitted the indexes related to the arc number since no confusion was possible. From now on, however, we need to distinguish the quantities on different arcs and we denote by u ± i , u i , v i and φ i the values of the corresponding variables on the i-th arc. On the outer boundaries, we could consider general boundary conditions:

(2.3) u + i (0, t) = α i (t)u - i (0, t) + β i (t), if i ∈ I out , u - i (L i , t) = α i (t)u + i (L i , t) + β i (t), if i ∈ O out .
For α i (t) = 1 and β i (t) = 0, we just recover the standard no-flux boundary condition (2.4) u + i (., t) = u - i (., t) (which is equivalent to v(., t) = 0). On the outer boundaries, we also consider no-flux (Neumann) boundary conditions for φ, which read

(2.5) ∂ x φ i (., t) = 0.
The no-flux boundary conditions mean that, on the boundary, the fluxes of cells and chemoattractants are null. This condition could be generalized, for example in the case when we assume that there is a production of fibroblasts on the boundary.

2.2.

Transmission conditions at a node. Now, let us describe how to define the conditions at a node; this is an important point, since the behavior of the solution will be very different according to the conditions we choose. Moreover, let us recall that the coupling between the densities on the arcs are obtained through these conditions. At node p ∈ N , we have to give values to the components such that the corresponding characteristics are going out of the node. Therefore, we consider the following transmission conditions at node:

(2.6)          u - i (L i , t) = j∈Ip ξ i,j u + j (L j , t) + j∈Op ξ i,j u - j (0, t), if i ∈ I p , u + i (0, t) = j∈Ip ξ i,j u + j (L j , t) + j∈Op ξ i,j u - j (0, t), if i ∈ O p ,
where the constant ξ i,j ∈ [0, 1] are the transmission coefficients: they represent the probability that a cell at a node decides to move from the i-th to the j-th arc of the network, also including the turnabout on the same arc. Let us notice that the condition differs when the arc is an incoming or an outgoing arc. Indeed, for an incoming (resp. outgoing) arc, the value of the function u + i (resp. u - i ) at the node is obtained through the system and we need only to define u - i (resp. u + i ) at the boundary.

These transmission conditions do not guarantee the continuity of the densities at node; however, we are interested in having the continuity of the fluxes at the node, meaning that we cannot loose nor gain any cells during the passage through a node. This is obtained using a condition mixing the transmission coefficients ξ i,j and the velocities of the arcs connected at node p. Fixing a node and denoting the velocities of the arcs by λ i , i ∈ I p ∪ O p , in order to have the flux conservation at node p, which is given by: (2.7)

i∈Ip λ i (u + i (L i , t) -u - i (L i , t)) = i∈Op λ i (u + i (0, t) -u - i (0, t)),
it is enough to impose the following conditions:

(2.8)

i∈Ip∪Op λ i ξ i,j = λ j , j ∈ I p ∪ O p .
Notice that, condition (2.7), can be rewritten in the uv variables as (2.9)

i∈Ip v i (L i , t) = i∈Op v i (0, t).
This condition ensures that the global mass µ(t) of the system is conserved along the time, namely:

(2.10)

µ(t) = N i=1 Li 0 u i (x, t)dx = µ 0 := N i=1
Li 0 u i (x, 0)dx, for all t > 0.

2.3. Dissipative transmission coefficients for the hyperbolic problem. It is sometimes useful to restrict our attention to the case of positive transmission coefficients of dissipative type, in the sense that they ensure energy decay of the solutions to the linear version of system (1.1), namely:

(2.11)

u t + v x = 0, v t + λ 2 u x = -v,
on a general network, with no-flux conditions (2.4) on the external nodes, and transmission conditions (2.6) at the internal nodes, always assuming the flux conservation condition (2.8) at nodes. To obtain the decay in time of the energy, which is defined by

E(t) = N i=1 Li 0 u 2 i (x, t) + v 2 i (x, t) λ 2 i dx 1/2
, it is sufficient to impose some equalities on the coefficients, as proved in [START_REF] Guaraldo | Some analytical results for hyperbolic chemotaxis model on networks[END_REF].

Proposition 1 ([12]

). The energy associated with the solutions to system (2.11), with no-flux conditions (2.4) on the external nodes, and transmission conditions (2.6) at the internal nodes, assuming condition (2.8), is decreasing if the transmission coefficients ξ i,j belong to [0, 1], and at every node p ∈ N , we have:

(2.12)

j∈Ip∪Op ξ i,j = 1 for all i ∈ I p ∪ O p .
Actually, in [START_REF] Guaraldo | Some analytical results for hyperbolic chemotaxis model on networks[END_REF], it is proved that under the assumptions of Proposition 1, it is possible to define a monotone generator of semigroup, and then a contraction semigroup, in the Sobolev space H 1 , for the linear transmission problem (2.11) on a network. Let us remark also that in the simplest case of a network composed by two arcs (one incoming and one outgoing, see next Figure 2), these conditions are also necessary in order to have the dissipation property. In such a case we have that dissipativity is given iff:

(2.13) max 0, λ 1 -λ 2 λ 1 ≤ ξ 1,1 ≤ 1, λ 2 (1 -ξ 2,2 ) = λ 1 (1 -ξ 1,1 ).
Using the previous relations and conditions on the coefficients ξ i,j given by (2.8), we obtain the values for the two missing coefficients:

(2.14)

ξ 1,2 = 1 -ξ 1,1 , ξ 2,1 = λ 1 λ 2 (1 -ξ 1,1 )
, so, we have only one degree of freedom.

2.4.

Transmission conditions for φ. Now let us consider the transmission conditions for φ in system (1.1). We complement conditions (2.3), (2.5), and (2.6) with a transmission condition for φ. As previously, we do not impose the continuity of the density of chemoattractant φ, but only the continuity of the flux at node p ∈ N . Therefore, we use the Kedem-Katchalsky permeability condition [START_REF] Kedem | Thermodynamic analysis of the permeability of biological membrane to non-electrolytes[END_REF], which has been first proposed in the case of flux through a membrane. For some positive coefficients κ i,j , we impose at node

(2.15) D i ∂ n φ i = j∈Ip∪Op κ i,j (φ j -φ i ), i ∈ I p ∪ O p .
The condition (2.16) κ i,j = κ j,i , i, j = 1, . . . , N yields the conservation of the fluxes at node p, that is to say

i∈Ip∪Op D i ∂ n φ i = 0.
Let us also notice that we can assume that κ i,i = 0, i = 1, . . . , N , which does not change condition (2.15). Finally, notice that the positivity of the transmission coefficients κ i,j , guarantees the energy dissipation for the equation for φ in (1.1), when the term in u is absent.

2.5. Stationary solutions. First we consider stationary solutions, which are known to drive the asymptotic behavior of the system. Let us consider the case of stationary solutions of system (1.1), complemented with boundary conditions (2.4), (2.5), (2.6), and (2.15). In the general case, we find on each arc the following solution :

(2.17)

       v i = constant, u i = exp(φ i /λ 2 i ) C i - v i λ 2 i x 0 exp(-φ i (y)/λ 2 i )dy , -D i φ i,xx = a i u i -b i φ i ,
which leads to solve, on each arc, the scalar non-local equation:

(2.18) -D i φ i,xx = a i exp(φ i /λ 2 i ) C i - v i λ 2 i x 0 exp(-φ i (y)/λ 2 i )dy -b i φ i ,
which has to be coupled at each node by the boundary conditions (2.4), (2.5), (2.6), and (2.15). We can prove easily that in the case of dissipative coefficients ξ i,j satisfying (2.8), (2.12) and the condition ξ i,j > 0, if all the fluxes v i are null, then the density u is continuous at a node, namely at a node p, the functions u i , i ∈ I p ∪ O p have all the same values. However, this is not the general case. For the simplest network composed of one incoming I = {1} and one outgoing O = {2} arc, represented in Fig. 2, we find on each interval that v 1 = v 2 = 0 from condition (2.4), and so we obtain the following local system for φ 1 and φ 2 :

(2.19) -D 1 φ 1,xx = a 1 C 1 exp(φ 1 /λ 2 1 ) -b 1 φ 1 , -D 2 φ 2,xx = a 2 C 2 exp(φ 2 /λ 2
2 )b 2 φ 2 , with boundary conditions (2.5) and (2.15) for φ 1 and φ 2 , which reads

∂ x φ 1 (L 1 ) = ∂ x φ 2 (0) = κ 1,2 (φ 2 (0) -φ 1 (L 1 )),
and

∂ x φ 1 (0) = ∂ x φ 2 (L 2 ) = 0.
We have also to take into account the following condition given by transmission condition (2.6) :

λ 2 ξ 2,1 C 1 exp(φ 1 (L 1 )/λ 2 1 ) = λ 1 ξ 1,2 C 2 exp(φ 2 (0)/λ 2 2
). Solving the corresponding system for φ 1 and φ 2 is a difficult task, even numerically, since an infinite number of solutions exist both for φ 1 and φ 2 , as in the case of a single interval [START_REF] Guarguaglini | Stability of constant states and qualitative behavior of solutions to a one dimensional hyper-bolic model of chemotaxis[END_REF], and it should be necessary to make them verify the above conditions at node. In order to simplify our study, we limit ourselves to state a result in the case of constant (in space) stationary solutions to system (1.1). Proposition 2. Let us consider a general network G = (N , A) and system (1.1) set on each arc of the network, complemented with boundary and transmission conditions (2.4), (2.5), (2.6), and (2.15).

(i) For general values of transmission coefficients ξ i,j satisfying (2.8), there is no non trivial constant stationary solution.

(ii) For the special case of transmission coefficients ξ i,j satisfying the dissipation relations (2.8) and (2.12) and of the ratios a i /b i being equal to the same constant on each arc, there exists a one-parameter stationary solution, which is constant by arc.

Proof. Take a constant (in space) stationary solution to system (1.1). This means that on each arc of the network, we have three constant values (u i , v i , φ i ), which satisfy v i = 0, since v i = u i φ ix = 0, a i u i = b i φ i , and boundary conditions (2.6), (2.15), which become in that case (2.20)

u i = j∈Ip∪Op ξ i,j u j , and
(2.21) 0 = j =i κ i,j (φ j -φ i ).
We remark that conditions (2.4) and (2.5) are automatically satisfied. (i) Denoting by N the number of arcs of the network, we have to fix therefore N unknowns to determine the stationary solution. Conditions (2.20) -(2.21) impose 4 equations by arc, unless the arc is connected to an outer node. In that case, there are only 2 conditions. To sum up, if we denote by N out the number of outer nodes, we need to satisfy 4N -2N out conditions. Taking into account relations (2.8), we obtain that equations (2.6) are linked and the system can be reduced to a system of 4N -2N out -N in conditions, where N in is the number of inner nodes, which is, generally speaking, greater than the number of unknowns. Therefore, unless some particular sets of coefficients κ i,j and ξ i,j , the only solution for previous system is the null one on each arc.

(ii) Now, let us consider transmission coefficients ξ i,j satisfying relations (2.8) and (2.12). We also assume that there exists a constant α such that, for all i, we have a i = αb i . In that case, we can find a stationary solution defined on each arc by (U, 0, αU ). Such kind of solution satisfies clearly the transmission condition (2.21), but satisfies also condition (2.20) with relations (2.12).

In the case (i) of the previous proposition, since the total initial mass is strictly positive and is preserved in time, we cannot expect the system to converge asymptotically to a stationary state which is constant on each arc and so nonconstant asymptotic solutions are expected. In the case (ii), the constant state can be reached, and U is determined by the total mass of the initial data.

Numerical schemes

Here we introduce our numerical schemes. We first give some details about schemes for system (1.1) on a single interval and the discretization of boundary conditions presented in [START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF]. Therefore, our main goal will be to generalize these schemes to the case of a network. In the two first subsections, we will concentrate on the discretization of the hyperbolic part, whereas the discretization of the parabolic part will be treated in subsection 3.3. [START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF] about AHO schemes for system (3.1) on a single interval. Let us consider a fixed single interval [0, L]. We define a numerical grid using the following notations: h is the space grid size, k is the time grid size and (x j , t n ) = (jh, nk) for j = 0, . . . , M + 1, n ∈ N are the grid points. In this subsection, we denote by w n,j the discretization of function w on the grid at time t n and at point x j for j = 0, . . . , M + 1 and n ≥ 0. We also use the notation f n,j for f (x j , t n ), where f is an explicitly known function depending on (x, t). Here we describe the discretization of system (1.1) with no-flux boundary conditions v(0, t) = v(L, t) = 0, denoting by f = φ x u and omitting the parabolic equation for φ. Since we also work with Neumann boundary conditions for the φ function, the function f will satisfy the following conditions on the boundary : f (0, t) = f (L, t) = 0. We therefore consider the following system (3.1)

Short review of the results from

u t + v x = 0, v t + λ 2 u x = f -v
and rewrite it in a diagonal form, using the usual change of variables (2.1),

     u - t -λu - x = 1 2 (u + -u -) - 1 2λ f, u + t + λu + x = 1 2 (u --u + ) + 1 2λ f. (3.2)
Set ω = u - u + , so that we can rewrite the system in vector form

(3.3) ω t + Λω x = Bω + F, with Λ = -λ 0 0 λ , B = 1 2 -1 1 1 -1 and F = 1 2λ
-f f . As shown in [START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF], to have a reliable scheme, with a correct resolution of fluxes at equilibrium, we have to deal with Asymptotically High Order schemes in the following form :

(3.4) ω n+1,i -ω n,i k + Λ 2h ω n,i+1 -ω n,i-1 - λ 2h (ω n,i+1 -2ω n,i +ω n,i-1 ) = B ℓ ω n,i+ℓ + ℓ=-1,0,1 D ℓ F n,i+ℓ .
With the following choice of the matrices (3.5)

B 0 = 1 4 -1 1 1 -1 , B 1 = 1 4 -1 1 0 0 , B -1 = 1 4 0 0 1 -1 , D 0 = 1 2 1 0 0 1 , D -1 = 1 2 0 0 0 1 , D 1 = 1 2
1 0 0 0 , we have a second-order AHO scheme on every stationary solutions, which is enough to balance the flux of the system at equilibrium. This means that the scheme is second order when evaluated on stationary solutions. Monotonicity conditions are satisfied if h ≤ 4λ and k ≤ 4h h + 4λ

, see [START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF] for more details. Let us mention that it should be easy to consider third-order AHO schemes, but for simplicity (these schemes require a fourth-order AHO scheme for the parabolic equation with a fivepoints discretization for φ x ), we prefer to limit our presentation to the second-order case.

Boundary conditions for scheme (3.4) have to be treated carefully, to enforce mass-conservation. In [START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF], the following boundary conditions were used :

(3.6) v n+1,0 = v n+1,M+1 = 0, u n+1,0 = 1 -λ k h u n,0 + λ k h u n,1 -k 1 h - 1 2λ v n,1 - k 2λ f n,1 , u n+1,M+1 = 1 -λ k h u n,M+1 + λ k h u n,M + k 1 h - 1 2λ v n,M + k 2λ f n,M .
These boundary conditions have been obtained by calculating the difference of the discrete mass at two successive computational times and defining u n+1,0 and u n+1,M+1 as a function of the discrete quantities computed at time t n in order to cancel exactly this difference. Consequently, the discrete mass will be preserved in time as the continuous mass L 0 u(x, t)dx is conserved for system (3.1) with boundary conditions v(0, t) = v(L, t) = 0, at the continuous level. This technique will be generalized in this paper to the case of a network.

3.2. The AHO scheme for system (3.1) in the case of a network. Let us consider a network as previously defined in Section 2. Each arc a i ∈ A, 1 ≤ i ≤ N , is parametrized as an interval a i = [0, L i ] and is discretized with a space step h i and discretization points x j i for j = 0, . . . , M i + 1. We still denote by k the time step, which is the same for all the arcs of the network. In this subsection, we denote by w n,j i the discretization on the grid at time t n and at point x j i of a function w i , i = 1, . . . , N on the i-th arc for j = 0, . . . , M i + 1 and n ≥ 0. Now, we consider the AHO scheme (3.4) on each interval, and we rewrite it in the u-v variables thanks to the change of variables (2.1), in order to define the discrete boundary and transmission conditions. We keep the possibility to use different AHO schemes on different intervals and therefore the coefficients of the scheme will be indexed by the number of the arc. Let R = 1 1 -λ λ be the matrix associated to the change of variables (2.1), namely such that u v = R u - u + . We rewrite (3.4) in the variables u and v as : (3.7)

u n+1,j i = u n,j i - k 2h i v n,j+1 i -v n,j-1 i + λ i k 2h i (u n,j+1 i -2u n,j i + u n,j-1 i ) + k 2 ℓ=-1,0,1 β ℓ u,u,i u n,j+ℓ i + 1 λ i ℓ=-1,0,1 β ℓ u,v,i v n,j+ℓ i + 1 λ i ℓ=-1,0,1 γ ℓ u,i f n,j+ℓ i , v n+1,j i = v n,j i - λ 2 i k 2h i u n,j+1 i -u n,j-1 i + λ i k 2h i (v n,j+1 i -2v n,j i + v n,j-1 i ) + k 2 λ i ℓ=-1,0,1 β ℓ v,u,i u n,j+ℓ i + ℓ=-1,0,1 β ℓ v,v,i v n,j+ℓ i + ℓ=-1,0,1 γ ℓ v,i f n,j+ℓ i , with coefficients β ℓ u,u , β ℓ u,v , β ℓ v,u , β ℓ v,v and γ ℓ u , γ ℓ v defined by (3.8) RB ℓ R -1 = 1 2 β ℓ u,u β ℓ u,v /λ λβ ℓ v,u β ℓ v,v , RD ℓ R -1 = 1 2 * γ ℓ u /λ * γ ℓ v .
Now, we define the numerical boundary conditions associated to this scheme. As before for equation (3.6), we need four boundary or transmission conditions to implement this scheme on each interval. Considering an arc and its initial and end nodes, there are two possibilities: either they are external nodes, namely nodes from the outer boundaries linked to only one arc, or they are internal nodes connecting several arcs together. The boundary and transmission conditions will therefore depend on this feature. Below, we will impose two boundary conditions (3.9)-(3.12) at outer nodes, and two transmission conditions (3.10)-(3.13) at inner nodes.

The first type of boundary conditions will come from condition (2.4) at outer nodes :

(3.9) v n+1,0 i = 0, if i ∈ I out , v n+1,Mi+1 i = 0, if i ∈ O out ,
where I out (resp. O out ) means that the arc is incoming from (resp. outgoing to) the outer boundary. The second one will come from a discretization of the transmission condition (2.6) at node p, that is to say

(3.10)          u n,Mi+1 -,i = j∈Ip ξ i,j u n,Mj+1 +,j + j∈Op ξ i,j u n,0 -,j , if i ∈ I p , u n,0 +,i = j∈Ip ξ i,j u n,Mj+1 +,j + j∈Op ξ i,j u n,0 -,j , if i ∈ O p .
However, these relations link all the unknowns together and they cannot be used alone. An effective way to compute all these quantities will be presented after equation (3.13) below. We still have two missing conditions per arc, which can be recovered by imposing the exact mass conservation between two successive computational steps. The discrete total mass is given by

I n tot = N i=1
I n i , where the mass corresponding to the arc i is defined as:

(3.11)

I n i = h i   u n,0 i 2 + Mi j=1 u n,j i + u n,Mi+1 i 2   .
Computing I n+1 tot -I n tot , we find:

I n+1 tot -I n tot = N i=1 h i k 2 1 k (u n+1,0 i -u n,0 i ) + 1 h i (v n,1 i + v n,0 i ) + λ i h i (u n,0 i -u n,1 i ) + β -1 u,u,i u n,0 i -β 1 u,u,i u n,1 i - 1 λ i β 1 u,v,i v n,1 i + 1 λ i β -1 u,v,i v n,0 i - 1 λ i γ 1 u,i f n,1 i -γ -1 u,i f n,0 i + h i k 2 1 k (u n+1,Mi+1 i -u n,Mi+1 i ) - 1 h i (v n,Mi i + v n,Mi+1 i ) + λ i h i (u n,Mi+1 i -u n,Mi i ) + β 1 u,u,i u n,Mi+1 i -β -1 u,u,i u n,M i - 1 λ i β -1 u,v,i v n,Mi i + 1 λ i β 1 u,v,i v n,Mi+1 i + 1 λ i γ 1 u,i f n,Mi+1 i -γ -1 u,i f n,Mi i .
We are going to impose boundary conditions such that the right-hand side in the previous difference is exactly canceled. On the outer boundaries we obtain the following type of boundary conditions, following equation (3.6) :

(3.12)                              u n+1,0 i = 1 -λ i k h i -kβ -1 u,u,i u n,0 i + k λ i h i + β 1 u,u,i u n,1 i -k 1 h i - β 1 u,v,i λ i v n,1 i + k λ i γ 1 u,i f n,1 i -γ -1 u,i f n,0 i , if i ∈ I out , u n+1,Mi+1 i = 1 -λ i k h i -kβ -1 u,u,i u n,Mi+1 i + k λ i h i + β 1 u,u,i u n,Mi i + k 1 h i - β -1 u,v,i λ i v n,Mi i - k λ i γ 1 u,i f n,Mi+1 i -γ -1 u,i f n,Mi i , if i ∈ O out ,
where I out and O out have the same meaning as previously. These expressions correspond to boundary conditions (3.6) in the case of a more general AHO scheme [START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF]. Then, using the conditions (3.12) to simplify the computation of I n+1 tot -I n tot and summing with respect to the nodes instead of the arcs, we can rewrite the remaining difference of mass in u ± variables as:

I n+1 tot -I n tot = p∈N i∈Op h i k 2 1 k u n+1,0 +,i + 1 k u n+1,0 -,i + u n,0 +,i - 1 k + 2 λ i h i + β -1 u,u,i + β -1 u,v,i + u n,0 -,i - 1 k + β -1 u,u,i -β -1 u,v,i -u n,1 +,i β 1 u,u,i + β 1 u,v,i + u n,1 -,i -2 λ i h i -β 1 u,u,i + β 1 u,v,i - 1 λ i (γ 1 u,i f n,1 i -γ -1 u,i f n,0 i ) + p∈N i∈Ip h i k 2 1 k u n+1,Mi+1 +,i + 1 k u n+1,Mi+1 -,i + u n,Mi+1 +,i - 1 k + β 1 u,u,i + β 1 u,v,i + u n,Mi+1 -,i - 1 k + 2 λ i h i + β 1 u,u,i -β 1 u,v,i -u n,Mi +,i ( 2
λ i h i + β -1 u,u,i + β -1 u,v,i ) + u n,Mi -,i -β -1 u,u,i + β -1 u,v,i + 1 λ i (γ 1 u,i f n,Mi+1 i -γ -1 u,i f n,Mi i ) .
Therefore, using the transmission conditions (3.10) for u n+1,Mi+1

-,i if i ∈ I p and for u n+1,0 +,i if i ∈ O p , we can split the equation interval by interval and obtain the following numerical boundary conditions: (3.13)

u n+1,Mi+1 +,i = h i   h i + j∈Ip∪Op h j ξ j,i   -1 × u n,Mi+1 +,i (1 -kβ 1 u,u,i -kβ 1 u,v,i ) + u n,Mi+1 -,i 1 -2k λ i h i -kβ 1 u,u,i + kβ 1 u,v,i + ku n,Mi +,i 2 λ i h i + β -1 u,u,i + β -1 u,v,i + ku n,Mi -,i (β -1 u,u,i -β -1 u,v,i ) - k λ i γ 1 u,i f n,Mi+1 i -γ -1 u,i f n,Mi i , if i ∈ u n+1,0 -,i = h i   h i + j∈Ip∪Op h j ξ j,i   -1 × u n,0 +,i (1 -2k λ i h i -kβ -1 u,u,i -kβ -1 u,v,i ) + u n,0 -,i 1 -kβ -1 u,u,i + kβ -1 u,v,i + ku n,1 +,i (β 1 u,u,i + β 1 u,v,i ) + ku n,1 -,i 2
λ i h i + β 1 u,u,i -β 1 u,v,i + k λ i (γ 1 u,i f n,1 i -γ -1 u,i f n,0 i ) , if i ∈ O p .
Once these quantities are computed, we can use equations (3.10) at time t n+1 , to obtain u n+1,Mi+1

-,i if i ∈ I p and u n+1,0

+,i if i ∈ O p .
In conclusion, we have imposed four boundary conditions (3.9), (3.10), (3.12), and (3.13) on each interval. Conditions (3.9) and (3.12) deal with the outer boundary and are written in the u-v variables, whereas conditions (3.10) and (3.13) deal with the node and are written in the u ± variables. Under these conditions, the total numerical mass is conserved at each step. Now, we have to discuss the consistency of all these conditions. First, conditions (3.9), (3.10) are imposed exactly. Besides, it has been proved in [START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF] that conditions (3.12), set on the outer boundary, are generally of order one and of order two on stationary solutions. Finally, we need to consider the consistency of the conditions (3.13) at node. We present here only the case i ∈ O p . Expanding in Taylor series up to order one, we get:

u n+1,0 -,i -1 + j∈Ip∪Op h j h i ξ j,i -1 × u n,0 +,i (1 -2k λ i h i -kβ -1 u,u,i -kβ -1 u,v,i ) + u n,0 -,i 1 -kβ -1 u,u,i + kβ -1 u,v,i + ku n,1 +,i (β 1 u,u,i + β 1 u,v,i ) + ku n,1 -,i 2
λ i h i + β 1 u,u,i -β 1 u,v,i + k λ i (γ 1 u,i f n,1 i -γ -1 u,i f n,0 i ) = u n,0 -,i 1 -1 + j∈Ip∪Op h j h i ξ j,i -1 1 + 2k λ i h i -u n,0 +,i 1 + j∈Ip∪Op h j h i ξ j,i -1 1 -2k λ i h i + O(k + i∈Ip∪Op h i ).
Now, to have consistency, namely to cancel the last two terms on the R.H.S., we need to impose the following condition linking the space and the time step on each arc :

(3.14)

h i = 2kλ i ,
which implies, thanks to (2.8):

  1 + j∈Ip∪Op h j h i ξ j,i   -1 = 1 2 .
Under this condition and using equations (3.2), expanding in Taylor series up to order three we find:

u n+1,0 -,i -1 + j∈Ip∪Op h j h i ξ j,i -1 × u n,0 +,i (1 -2k λ i h i -kβ -1 u,u,i -kβ -1 u,v,i ) + u n,0 -,i 1 -kβ -1 u,u,i + kβ -1 u,v,i + ku n,1 +,i (β 1 u,u,i + β 1 u,v,i ) + ku n,1 -,i 2
λ i h i + β 1 u,u,i -β 1 u,v,i + k λ i (γ 1 u,i f n,1 i -γ -1 u,i f n,0 i ) = k 2 u n,0 -,i (β 1 u,v,i + β -1 u,u,i -β 1 u,u,i -β -1 u,v,i ) + k∂ t u n,0 -,i + k 2 2 ∂ tt u n,0 -,i -kλ i 1 + k(β 1 u,u,i -β 1 u,v,i ) ∂ x u n,0 -,i -k 2 λ 2 i ∂ xx u n,0 -,i + k 2 u n,0 +,i (β -1 u,u,i + β -1 u,v,i -β 1 u,u,i -β 1 u,v,i ) -k 2 λ i β 1 u,u,i + β 1 u,v,i ∂ x u n,0 +,i - k 2λ i (γ 1 u,i -γ -1 u,i )f n,0 i -k 2 γ 1 u,i ∂ x f n,0 i + O(k 3 ) = k 2 u n,0 -,i (-1 + β 1 u,v,i + β -1 u,u,i -β 1 u,u,i -β -1 u,v,i ) + u n,0 +,i (1 + β -1 u,u,i + β -1 u,v,i -β 1 u,u,i -β 1 u,v,i ) - 1 λ i (1 + γ 1 u,i -γ -1 u,i )f n,0 i + k 2 1 2 ∂ tt u n,0 -,i -λ i ∂ tx u n,0 -,i -λ i ( 1 2 + β 1 u,u,i -β 1 u,v,i )∂ x u n,0 -,i -λ i (- 1 2 + β 1 u,u,i + β 1 u,v,i )∂ x u n,0 +,i -(γ 1 u,i + 1 2 )∂ x f n,0 i + O(k 3 ).
Thanks to this development we can state our general result of consistency.

Proposition 3. Given a general scheme in the form (3.4), the conditions (3.13) at node are consistent only if on each arc the condition (3.14) is verified. To have the second order accuracy at node the following conditions on the coefficients of the scheme have to be verified:

(3.15) β 1 u,u,i = β -1 u,u,i , β 1 u,v,i -β -1 u,v,i = 1, γ -1 u,i -γ 1 u,i = 1.
Moreover, to have a third order accuracy for stationary solutions, we need :

(3.16) β 1 u,u,i = β -1 u,u,i = 0, β 1 u,v,i = -β -1 u,v,i = 1 2 , γ 1 u,i = -γ -1 u,i = - 1 2 .
Notice that, all these conditions are satisfied for the Roe scheme defined by (3.5).

3.3.

Discretization of the parabolic equation for φ in system (2.2). Now, let us explain how to compute the approximations f n+1,j i of the function f on the arc i at discretization point x j i and time t n+1 needed for computing (3.7), (3.12), and (3.13). Referring to system (2.2), we have f = φ x u, where φ satisfies the parabolic equation φ t -D φ xx = aubφ on each arc. Boundary conditions for φ are given by equations (2.5) on the outer boundary and (2.15) at a node.

We solve the parabolic equation, using a finite differences scheme in space and a Crank-Nicolson method in time, namely an explicit-implicit method in time.

Therefore, we will have the following equation for φ n,j i , 1 ≤ j ≤ M i , (3.17)

φ n+1,j i = φ n,j i - D i k 2h 2 i -φ n,j+1 i + 2φ n,j i -φ n,j-1 i - D i k 2h 2 i -φ n+1,j+1 i + 2φ n+1,j i -φ n+1,j-1 i + a i k 2 (u n+1,j i + u n,j i ) - b i k 2 (φ n+1,j i + φ n,j i
). Now, let us find the two boundary conditions needed on each interval. As in subsection 3.2, the boundary conditions will be given in the case of an outer node and in the case of an inner node. On the outer boundary, condition (2.5) for φ is discretized using a second order approximation, which is

(3.18)      φ n,0 i = 4 3 φ n,1 i - 1 3 φ n,2 i , if i ∈ I out , φ n,Mi+1 i = 4 3 φ n,Mi i - 1 3 φ n,Mi-1 i , if i ∈ O out .
Let us now describe our numerical approximation for the transmission condition (2.15) which, as the transmission condition for the hyperbolic part (2.6), couples the φ functions of arcs having a node in common. Condition (2.15) is discretized using the same second-order discretization formula as before, namely we have at node p,

φ n,Mi+1 i = 4 3 φ n,Mi i - 1 3 φ n,Mi-1 i + 2 3
h i D i j∈Ip κ i,j (φ n,Mj +1 j -φ n,Mi+1 i ) + 2 3 h i D i j∈Op κ i,j (φ n,0 j -φ n,Mi+1 i ), if i ∈ I φ n,0 i = 4 3 φ n,1 i - 1 3 φ n,2 i + 2 3
h i D i j∈Ip κ i,j (φ n,Mj +1 j -φ n,0 i ) + 2 3 h i D i j∈Op κ i,j (φ n,0 j -φ n,0 i ), if i ∈ O p .
These relations can be rewritten as :

(3.19)   1 + 2 3 h i D i j∈Ip∪Op κ i,j   =η p i φ n,Mi+1 i = 4 3 φ n,Mi i - 1 3 φ n,Mi-1 i + 2 3
h i D i j∈Ip κ i,j φ n,Mj +1 j + 2 3
h i D i j∈Op κ i,j φ n,0 j , if i ∈   1 + 2 3 h i D i j∈Ip∪Op κ i,j   =η p i φ n,0 i = 4 3 φ n,1 i - 1 3 φ n,2 i + 2 3
h i D i j∈Ip κ i,j φ n,Mj +1 j + 2 3
h i D i j∈Op κ i,j φ n,0 j , if i ∈ O p .
Let us remark that the previous discretizations are compatible with relations (3.18) considering that for outer boundaries the coefficients κ i,j are null. Therefore, in this case, the value of η out i is just equal to 1. Since equations (3.19) are coupling the unknowns of all arcs altogether, we have to solve a large system which contains all the equations of type (3.17) and also the discretizations of transmission conditions (3.19).

Once the values of φ n+1,j i are known, we can compute a second-order discretization of the derivatives of φ which gives the values of the f function, namely :

φ n+1,j x,i =                1 2 h i φ n+1,j+1 i -φ n+1,j-1 i , 1 ≤ j ≤ M i , 1 2 h i -φ n+1,2 i + 4φ n+1,1 i -3φ n+1,0 i , j = 0, 1 2 h i φ n+1,Mi-1 i -4φ n+1,Mi i + 3φ n+1,Mi+1 i , j = M i + 1.
The discretization of f needed at equations (3.7),(3.12), and (3.13) is therefore given by f n+1,j i = φ n+1,j

x,i u n+1,j i .

Numerical tests

Here we present some numerical experiments for system (1.1) on networks, with the use of the methods introduced in Section 3, namely the second-order AHO scheme for the hyperbolic part, complemented with the Crank-Nicolson scheme for the parabolic part. We start with a simple test for the AHO scheme on the hyperbolic part of Section 3 in the case of a simplified system, where φ x is equal to a constant α on each arc, for which we know the exact stationary states. 4.1. Case φ x constant. For this example, we omit the equation for φ so that the system becomes (4.1)

     u + t + λu + x = 1 2λ (α -λ)u + + (α + λ)u -, u - t -λu - x = - 1 2λ (α -λ)u + + (α + λ)u -.
This system is suitable to test the accuracy of the numerical approximation, since it is easy to compute its asymptotic stationary solutions. We also rewrite the previous system (4.1) using the usual change of variables (2.1) which gives (4.2)

u t + v x = 0, v t + λ 2 u x = α u -v,
with α a constant. To satisfy the subcharacteristic condition in [START_REF] Natalini | Convergence to equilibrium for the relaxation approximation of conservation laws[END_REF], we also assume that (

Let us explain how to find the stationary states in the case of the two-arcs network of Figure 2. The method can be easily generalized to more complex networks. In that case, the stationary solutions satisfy the following equations on the intervals I 1 and I 2 :

v i,x = 0, λ 2 i u i,x = α i u i -v i , that is to say (4.4) v i = constant, u i = C i exp(α i x/λ 2 i ) + v i /α i .
Since both intervals are connected to the outer boundary, due to boundary condition (2.4), we have v 1 = v 2 = 0. Therefore we obtain non constant solutions on each arc, given by u

± i = u i 2 = C i 2 exp(α i x/λ 2 i )
and the constants C i are computed thanks to condition (2.6). Remark that, in that case, we do not expect to have asymptotic states given by constant stationary solutions, since the only possible constant solution is the null one, which will be unsuitable, due to the constraint of the conservation of mass. Set (4.5)

C 1 = C 1 λ 1 exp(α 1 L 1 /λ 2 1 ), C 2 = C 2 λ 2 .
These constants solve the following system :

(4.6) M C = λ 1 (ξ 1,1 -1) λ 2 ξ 1,2 λ 1 ξ 2,1 λ 2 (ξ 2,2 -1)
C 1 C 2 = 0.
According to (2.8), Ker M = {0}, and so we have at most one equation and two unknowns. Therefore, there exists at least one family of non trivial stationary solutions to system (4.2) and exactly one family when dim Ker M = 1. Remark that in the general case of a single node with an arbitrary number of incoming and outcoming arcs, assuming that all coefficients ξ i,j are strictly positive -or more generally, that the matrix formed by these coefficients is irreducible, which is somewhat meaningful in the biological context -, we can prove that we have exactly dim Ker M = 1, thanks to the classical Perron-Frobenius theorem.

In the case we are looking for an asymptotic state as a stationary state of the system, we can also take into account the conservation for mass. In that case, the stationary state we compute should have the same mass as the initial datum. More precisely, according to equation

µ 0 = 2 i=1 Li 0 C i exp αx λ i 2 dx = 2 i=1 C i λ 2 i α exp αL i λ i 2 -1 ,
we have that the free parameter is fixed by the mass conservation.

In particular we set L 1 = 4, L 2 = 1, α i = α = 0.5, λ 1 = 2, λ 2 = 1 and take the dissipative transmission coefficients ξ 1,1 = 0.8, ξ 2,1 = 0.4, ξ 1,2 = 0.2, ξ 2,2 = 0.6. If µ 0 = 250, the system is solved by C 1 ∼ 28.13 and C 2 ∼ 56.25, so that the stationary solutions are u 1 = C 1 exp(x/8) and u 2 = C 2 exp(x/2), with C 1 ∼ 34.12 and C 2 ∼ 56.25. The numerical simulations provide the asymptotic densities plotted in Fig. 3 and we notice a nice agreement with the stationary solutions computed analytically. Remark that densities are continuous at the node as explained in Section 2.5 for dissipative coefficients and vanishing fluxes.

In Fig. 4 we present the log-log plot of the error in the L 1 -norm computed using the formula (4.9) of Section 4, between the approximated and the asymptotic solutions to system (4.2). The results in Fig. 4 show that the AHO approximation scheme provides the stationary solutions of the simplified hyperbolic model (4.2) with an accuracy of first order, and the error for the flux function v tends clearly to zero, faster than for the function u.

More examples and results showing the asymptotic behavior of solutions to the simple problem (4.2) on larger networks can be found in [START_REF] Bretti | A numerical scheme for a hyperbolic relaxation model on networks[END_REF], while some analytical results are given in [START_REF] Guaraldo | Some analytical results for hyperbolic chemotaxis model on networks[END_REF]. 4.2. Asymptotic solutions to the full system (2.2). Next, we deal with the full system (2.2), which now include the chemotaxis equation. First, we consider again a network with only two arcs. We take the following data: the total mass µ 0 = 160 distributed as a small perturbation of the value 20 on two arcs of length L 1 = 6 and L 2 = 2, see Fig. 5, a i = b i = 1, u i (x, 0) = φ i (x, 0) and v i (x, 0) = 0, i = 1, 2 and λ 1 = 5, λ 2 = 4. In the next figures we represent the asymptotic stable solutions to system (2.2) on the two-arcs network, produced by our scheme. All the solutions are plotted at a time where the stationary state is already reached. In particular, in Fig. 6 we plot a constant solution obtained using the dissipative transmission coefficients of Section 2.3. In that case we can observe what was explained in Section 2.5, namely that in the case of two arcs and one node, there exist particular dissipative transmission coefficients, such that the asymptotic stationary solutions are constants on all the arcs. In Fig. 7 we plot the more common case of non-constant solutions, obtained using different parameters and non dissipative coefficients. In both cases the limit flux function v is equal to error for u error for v order 1

Figure 4. Log-log plot of the error in L 1 norm between the approximated and the the asymptotic solutions, as a function of the space step, for the solutions to system (4.2). The error is displayed in blue for u, and in green for v. Initial data are distributed on the network as a symmetric perturbation of the value 50. We used different space steps satisfying condition (3.14), with

λ 1 = 2, λ 2 = 1, L 1 = 4, L 2 = 1, µ 0 = 250, T = 100.
zero everywhere, since for the stationary solution the flux is constant, the flux on the external nodes is zero, and all the arcs are connected to external nodes.

Let us now consider a larger network composed of twelve nodes and four arcs, see Fig. 8. We choose some non dissipative transmission coefficients, given in Table 1, in order to satisfy condition (2.8). Let us consider as initial condition on the incoming arc 5, the function plotted in Fig. 9, where we put a small symmetric perturbation of the constant state u = 110.

In this case it is hard to compute analytically the stationary solutions. We only know that non-constant solutions are generally expected, according to the discussion in Section 2.5. In Fig. 10 we plot the asymptotic densities on the network node by node, starting from North-East and proceeding in a clockwise direction. Notice that most of the arcs are repeated in the different figures. In Fig. 11 the asymptotic fluxes are represented, and again our scheme is able to stabilize them correctly. We notice that the fluxes of arcs connected to outer boundaries vanish, whereas the fluxes of inner arcs, even if they are constant, are different from zero.

4.3.

Instabilities: the appearance of numerical blow-up. Let us consider some cases that present a strong asymptotical instability. Indeed, for some values of the parameters of the problem, namely of the arc's length L and the cell velocity λ, in connection with the total mass distributed on the arcs of the network, we can observe increasing oscillations, which eventually may cause the blow-up of solutions. It is important to notice that the blow-up can be already observed for this model even for a single arc, see Example 1 below, when the total mass µ 0 is large with respect to the characteristic parameters L and λ. However, here the presence of more arcs, and so, a greater total length and total mass, makes this kind of phenomenon much more frequent.

Example 1. Here we assume that we have only one interval with L = 1 and λ = 10 and we take, as initial condition for the density and the chemoattractant, a

λ i = λ = 10, L i = 1, a i = b i = D i = 1.
The transmission coefficients can be found in Table 1. The total initial mass µ 0 = 1320 is distributed as a perturbation of the constant state 110 on arc 5 as in Fig. 9 and as the constant density 110 on the other arcs, with h i = h = 0.01, k = 0.0005. symmetric perturbation of a constant state C 0 = 9000. The total mass is µ 0 = 9000, as shown in Figure 12. The solution presents a clear blow-up at time T = 0.1, see Fig. 13. This blow-up seems associated to non physical negative values of the density function u, and it is observed in the same way even for refined meshes (see Table 2 for the case of two arcs). This is not surprising, since the quasimonotonicity of the system, see again [START_REF] Natalini | Convergence to equilibrium for the relaxation approximation of conservation laws[END_REF], is violated when the gradient φ x is larger than λ.

= 0.2, ξ 2,2 = 0.2, ξ 1,1 = 0.1, ξ 2,1 = 0.3, ξ 8,1 = 0.3, ξ 7,1 = 0.3, Node N-E ξ 1,2 = 0.2, ξ 2,2 = 0.2, ξ 8,2 = 0.3, ξ 7,2 = 0.3, ξ 1,8 = 0.2, ξ 2
Example 2. Here we take two arcs of length L 1 = 6 and L 2 = 2 and the initial density as in Fig. 5, with a i = b i = 1, u i (x, 0) = φ i (x, 0), and v i (x, 0) = 0, i = 1, 2. Then we change the values of velocities λ 1 and λ 2 in order to see how they influence the behavior of solutions to system (1.1). At the junction we assume transmission and dissipative coefficients, taking ξ 1,1 = 0.96 and then satisfying equations (2.13)- (2.14). What we observe is that solutions blow up in finite time or not according to the relative values of λ 1 and λ 2 , as it is shown in Figure 14. More precisely, we can observe three different regimes. If λ 2 is large with respect to 1 λ1-2 , solutions stay bounded and converge to stationary solutions (green "x" in Figure 14). If λ 1 is small with λ 2 large enough, then solutions blow up in finite time (red "+" in Figure 14). Finally, there is a small region in between, λ 1 around the value 3 and λ 2 small enough, such that solutions present a large spike at the boundaries (marked by blue asterisks "*"). . Schematization of the regions describing the behavior of solution for µ 0 =160 and the velocities λ 1 and λ 2 varying: blowup (marked by red crosses "+"), solutions with a spike at the boundaries (marked by blue asterisks "*") and stable stationary solutions (marked by green "x").

Let us now focus on the blow-up behavior. Referring to Fig. 14, we can choose a pair of velocities belonging to the blow-up region marked by red crosses "+", to say λ 1 = 1 and λ 2 = 2. The time step just before the numerical blow-up time of corresponding solutions, starting from initial data as in Fig. 5, is plotted in Fig. 15. Even if apparently we are close to the transmission point, there are many grid points separating it from the blow-up point. To show that the blow-up is not just a numerical artifact, we perform the same simulation with the same data, but on refined grids. In Table 2 we report the blow-up time of solutions to system (1.1) for a fixed global mass µ 0 when either the CFL condition ν = k h λ or h go to zero. Out of the case of ν = 1, which appears to be more unstable, the blow-up time is independent of the meshes and has to be considered to occur in the analytical solutions. Table 2. Blow-up times of the solutions to system (1.1) when either the CFL condition ν = k h λ or h go to zero, with transmission coefficients of dissipative type, L 1 = 6, L 2 = 2, λ 1 = 1, λ 2 = 2, µ 0 = 160. where w T j (h) denotes the numerical solution obtained with the space step discretization equal to h, computed in x j at the final time T . The total L 1 -error is (4.10)

Blow-up time

T OT err = N i=1 e i (h i ).

Table 3 shows the L 1 -error (4.9) on the asymptotic solutions u, φ and v and order of convergence (4.7) of the approximation scheme applied to the considered network.

The results in Table 3 show the effectiveness of AHO approximation scheme in the solution of the transmission problem represented by the hyperbolic model (1.1). We notice indeed that even in this more general case the scheme still keeps a formal accuracy of first order, although the interactions at the boundaries could deteriorate its accuracy.

Figure 1 .

 1 Figure 1. First example of network.

Figure 2 .

 2 Figure 2. One incoming and one outgoing arc connected at a node.

Figure 3 .

 3 Figure 3. Comparison between the densities of the exact and the numerical stationary solutions on arcs 1 and 2 obtained for λ 1 = 2, λ 2 = 1, α i = α = 0.5, initial mass µ 0 = 250 distributed on the network as a symmetric perturbation of the value 50, L 1 = 4, L 2 = 1, dissipative coefficients ξ 1,1 = 0.8, ξ 2,1 = 0.4, ξ 1,2 = 0.2, ξ 2,2 = 0.6 and time T = 28.

19

 19

Figure 5 .

 5 Figure 5. Initial data corresponding to the total mass µ 0 = 160.

Figure 6 .

 6 Figure 6. Asymptotic solution for λ 1 = 5, λ 2 = 4, dissipative coefficients ξ 1,1 = 0.8, ξ 2,1 = 0.25, ξ 1,2 = 0.2, ξ 2,2 = 0.75, T = 7.7.

Figure 7 .Figure 8 .

 78 Figure 7. Asymptotic solution at time T = 30 for λ 1 = 5, λ 2 = 4, in case of non-dissipative coefficients ξ 1,1 = 0.8, ξ 2,1 = 0.25, ξ 1,2 = 0.24, ξ 2,2 = 0.7.

Figure 9 .Figure 10 .

 910 Figure 9. Initial condition for u and φ on arc 5 of the network presented in Fig. 8.

ξ 12 ,

 12 12 = 0.1, ξ 11,12 = 0.3, ξ 3,12 = 0.3, ξ 4,12 = 0.3, Node S-W ξ 12,11 = 0.2, ξ 11,11 = 0.2, ξ 3,11 = 0.3, ξ 4,11 = 0.3, ξ 12,3 = 0.2, ξ 11,3 = 0.2, ξ 3,3 = 0.4, ξ 4,3 = 0.2, ξ 12,4 = 0.5, ξ 11,4 = 0.1, ξ 3,4 = 0.2, ξ 4,4 = 0.2, ξ 3,3 = 0.1, ξ 10,3 = 0.3, ξ 9,3 = 0.3, ξ 2,3 = 0.3, Node S-E ξ 3,10 = 0.2, ξ 10,10 = 0.2, ξ 9,10 = 0.3, ξ 2,10 = 0.3, ξ 3,9 = 0.2, ξ 10,9 = 0.2, ξ 9,9 = 0.4, ξ 2,9 = 0.2, ξ 3,2 = 0.5, ξ 10,2 = 0.1, ξ 9,2

Figure 11 .

 11 Figure 11. The asymptotic fluxes of the arcs of the network composed of 12 arcs and 4 nodes at time T = 30, with the same data as Fig. 10.

Figure 12 .Figure 13 .

 1213 Figure 12. The initial condition u 0 (x) is a symmetric perturbation of a constant state C 0 = 9000, the total mass is µ 0 = 9000.

2 Figure 14

 214 Figure 14. Schematization of the regions describing the behavior of solution for µ 0 =160 and the velocities λ 1 and λ 2 varying: blowup (marked by red crosses "+"), solutions with a spike at the boundaries (marked by blue asterisks "*") and stable stationary solutions (marked by green "x").

Figure 15 . 4 . 4 . 2 e i (h i ) e i hi 2 , i = 1 , 4

 15442214 Figure 15. Blow-up at time T = 4, for initial data as in Fig.5, with L 1 = 6, L 2 = 2, λ 1 = 1 and λ 2 = 2, dissipative coefficients with ξ 1,1 = 0.96, the total mass is equal to µ 0 = 160: on the left the density u and the concentration φ, on the right the flux v. The space steps are equal to h 1 = 0.001, h 2 = 0.002.

Table 3 .

 3 Acknowledgements : The research leading to these results has received funding from the European Union Seventh Framework Programme [FP7/2007-2013] under 8.87206e-05 0.982631 8.87207e-05 -0.058657 1.44060e-07 0.00625 0.980243 4.41941e-05 0.990856 4.41954e-05 0.666605 1.49949e-07 0.003125 0.986317 2.20550e-05 0.992983 2.20651e-05 0.863690 9.43741e-08 0.0015625 0.937936 1.10172e-05 0.937109 1.10280e-05 0.955806 5.17981e-08 Orders and errors of the approximation scheme for the solutions to system (1.1), L i = 1, λ i = 4, i = 1, 2, µ 0 = 120.056, T = 25. grant agreement n257462 HYCON2 Network of excellence. This work has also been partially supported by the PRIN project 2008-2009 "Equazioni iperboliche non lineari e fluidodinamica" and by the ANR project MONUMENTALG, ANR-10-JCJC 0103.

Table 1 .

 1 Transmission coefficients used for the numerical simulations of Figures 10 and 11 given node by node.

		,8 = 0.2,	ξ 8,8 = 0.4, ξ 7,8 = 0.2,
	ξ 1,7 = 0.5,	ξ 2,7 = 0.1,	ξ 8,7 = 0.2, ξ 7,7 = 0.2,
	ξ 5,5 = 0.1,	ξ 4,5 = 0.3,	ξ 1,5 = 0.3, ξ 6,5 = 0.3,
	Node N-W ξ 5,4 = 0.2,	ξ 4,4 = 0.2,	ξ 1,4 = 0.3, ξ 6,4 = 0.3,
	ξ 5,1 = 0.2,	ξ 4,1 = 0.2,	ξ 1,1 = 0.4, ξ 6,1 = 0.2,
	ξ 5,6 = 0.5,	ξ 4,6 = 0.1,	ξ 1,6 = 0.2, ξ 6,6 = 0.2.

ℓ=-1,0,1