{'fraction_non_alphanumeric': 0.05810043668122271, 'fraction_numerical': 0.018799126637554585, 'mean_word_length': 3.830819533804451, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Chemotaxis on networks: Analysis and numerical approximation

We consider the Keller-Segel model of chemotaxis on one-dimensional networks. Using a variational characterization of solutions, positivity preservation, conservation of mass, and energy estimates, we establish global existence of weak solutions and uniform bounds. This extends related results of Osaki and Yagi to the network context. We then analyze the discretization of the system by finite elements and an implicit time-stepping scheme. Mass lumping and upwinding are used to guarantee the positivity of the solutions on the discrete level. This allows us to deduce uniform bounds for the numerical approximations and to establish order optimal convergence of the discrete approximations to the continuous solution without artificial smoothness requirements. In addition, we prove convergence rates under reasonable assumptions. Some numerical tests are presented to illustrate the theoretical results.


Introduction
Back in the 1970s, Keller and Segel [23] introduced their celebrated model of chemotaxis describing the collective movement of cellular organisms in response to the distribution of a chemical substance. According to [8,19], the minimal system given by ∂ t u − div(α∇u − χu∇c) = 0, ∂ t c − div(β∇c) + γc = δu, serves as a prototype for studying many interesting mathematical aspects of chemotaxis. In the context of biological applications, u denotes the density of the population of interest and c is the concentration of the chemoattractant. The differential equations are usually augmented by homogeneous Neumann boundary conditions ∂ n u = ∂ n c = 0 which leads to global conservation of the population and to preservation of positivity in both variables.
The analysis of this parabolic-parabolic model of chemotaxis is well understood by now. Local existence of a unique solution was proven by Yagi in [36] and global existence was established by Nagai, Senba, and Yoshida [25] under a smallness condition on the initial data. For large initial mass, blow-up in finite time was demonstrated in [15,22] for dimension d ≥ 2. Finite-time blow-up cannot occur in dimension d = 1 and thus solutions for the nonlinear parabolic system exist globally in time [18,28]. We refer to [19,21] for an overview about various models of chemotaxis and further theoretical results.
Apart from the analysis, also the numerical approximation of chemotaxis models has attracted significant interest in the literature. Nakaguchi and Yagi [26] proposed fully discrete schemes obtained by Galerkin approximation in space and Runge-Kutta methods in time, and they established convergence rates under appropriate smoothness assumptions on the true solution. Filbet [13] studied finite volume approximations and proved convergence also for non-smooth solutions. Higher order schemes have been studied in [7,11]. In order to retain the conservation of mass and the positivity of the numerical solution, Saito [31] proposed and analyzed an upwind finite element scheme for multidimensional problems; a slightly different approach was considered [32]. Let us also mention recent work [3,14,27] concerning the approximation of hyperbolic models of chemotaxis.
In this paper, we study chemotaxis on one-dimensional networks modeled by a system of partial differential-algebraic equations on finite metric graphs [24]. A one-dimensional version of the minimal system is imposed on every edge of the graph and complemented by algebraic coupling conditions at the vertices to ensure continuity of the solution and conservation of mass across junctions. A corresponding system has been proposed and investigated by Borsche et al. [2], who considered a positivity preserving finite volume discretization and established the well-posedness of their scheme. Let us also mention recent work by Camilli and Corrias [6], who considered problems with constant coefficients and proved well-posedness by perturbation arguments and using the mapping properties of the heat semi-group on networks. Numerical methods for hyperbolic models of chemotaxis on networks were investigated by Bretti et al. [5].
Before we proceed, let us briefly discuss the main contributions of our manuscript: By extending the functional analytic framework of Osaki and Yagi [28], we consider the chemotaxis problem on networks as a semilinear parabolic system which allows us to establish existence of a local solution by Galerkin approximation, energy estimates, and perturbation arguments. Following the ideas of [20,28], we further show that the solution remains positive, provided that the initial values are positive, and we prove that the total mass of the population is conserved for all time. This yields uniform a-priori estimates for the L 1 -norm of the density u and allows us to derive sharper energy estimates by which we can show that the solution can at most grow polynomially in time and hence exists globally. These results can be seen as a natural generalization of those in [20,28] to one-dimensional networks. However, we use somewhat different energy estimates in our proofs which allows us to apply our analysis also to problems with discontinuous model parameters and to networks of rather general topology. Our method of proof also differs from that in [6] and our results are more general, in particular, our analysis covers the case of non-constant and discontinuous coefficients. As preparation for the second part of the manuscript, we also establish higher regularity of solutions.
After having proved the global existence and uniqueness of solutions, we turn to their systematic numerical approximation. For the discretization, we here consider a Galerkin approximation in space by finite elements combined with an implicit timestepping scheme. In order to ensure positivity of the discrete solutions, we employ a mass-lumping strategy and an upwind discretization for the convective term. The resulting scheme has a similar structure as that considered by Saito [31] for chemotaxis problems in multiple dimensions, but the formulation of our scheme is closer to that of the continuous problem, which facilitates the analysis substantially. Some alternative but related approaches can be found in [13] and [32]. Using similar methods of proof as on the analytical level, we derive uniform bounds for the discrete approximations and we establish convergence of the numerical solution and order optimal convergence rates under reasonable smoothness assumptions. Our analysis is somewhat sharper and more general than that presented in [31]. In particular, we do not require a strong restriction on the time step to guarantee the stability of our fully discrete scheme and we obtain convergence in the general case without artificial smoothness assumptions.
The remainder of the manuscript is organized as follows: In Section 2, we introduce our notation and the problem under investigation, and we give a variational characterization of solutions which will be the basis for the rest of the manuscript. In Section 3, we establish the existence and uniqueness of solutions and derive uniform bounds that grow at most polynomially in time. In addition, we prove higher regularity of solutions under natural smoothness and compatibility conditions on the initial data. The numerical approximation is introduced in Section 4 and we establish uniform global bounds for the discrete solutions. In Sections 5 and 6, we prove the convergence of discrete solutions to the true solution and we establish order optimal convergence rates under reasonable smoothness assumptions. For illustration of our theoretical findings, we present some numerical tests in Section 7 and we close the presentation with a short summary and a discussion of possible directions for future research.

Preliminaries
We start by introducing our notation and then formally state the chemotaxis problem on the network to be considered for the rest of the paper.
2.1. Network topology. Let (V, E) be a finite directed and connected graph [1] with vertices v ∈ V and edges e ∈ E. To any edge e = (v 1 , v 2 ) pointing from vertex v 1 to v 2 , we set n e (v 1 ) = −1 and n e (v 2 ) = 1. The matrix with entries N ij = n e j (v i ) then is the incidence matrix of the graph. For any v ∈ V , we denote by E(v) = {e ∈ E : e = (v, ·) or e = (·, v)} the set of edges starting or ending at v, and for e ∈ E we define V(e) = {v ∈ V : e = (v, ·) or e = (·, v)}. We further denote by V b = {v ∈ V : |E(v)| = 1} the set of boundary vertices and call V 0 = V \ V b the set of interior vertices. A small example illustrating our notation is presented in Figure 2.1.

Function spaces.
To any edge e ∈ E we associate a positive length e > 0 and with some abuse of notation, we always identify the topological edge e with the geometric interval [0, e ] in the sequel. Let be the vector with entries e . Following the notation of [24], we call the triple G = (V, E, ) a metric graph. We further denote by The non-zero entries of the incidence matrix are The corresponding norm is given by v L 2 (E) = v, v 1/2 E and the spaces L p (E), 1 ≤ p ≤ ∞ are defined accordingly. In addition, we will also make use of the function space consisting of continuous functions with square integrable weak derivatives. This space is complete when equipped with the norm defined by v 2 Note that by continuity and density, the scalar product ·, · E can be extended to the duality product on H 1 (E) × H 1 (E), for which we use the same symbol.
For T > 0 and some Banach space X, we denote by L p (0, T ; X) the space of measurable functions with values v(t) ∈ X and with finite norm v p L p (0,T ;X) = T 0 v(t) p X . Spaces of differentiable functions in time are denoted by W k,p (0, T ; X) and equipped with their natural norms. As usual, we write H k (0, T ; X) = W k,2 (0, T ; X) for convenience. Let us recall that the embedding of the energy space into C([0, T ]; L 2 (E)) is continuous, which can be proven with similar arguments as on single intervals. As a consequence, the evaluation v(t) is well-defined for functions v ∈ W (0, T ) and one has a uniform bound with a constant C independent of v; we refer to [12] for details and further references.
2.3. Problem statement. Let G = (V, E, ) be a finite directed metric graph as introduced above. On every edge e ∈ E, the chemotactic movement shall be described by e ∈ E, t > 0, (2.2) with model parameters α, β, γ, δ, χ to be specified below. Recall that f e = f | e denotes the restriction of a function f onto the edge e. In addition to the above equations, we assume that the solution is continuous across vertices, i.e., and we require that the population and concentration are conserved at all vertices, i.e., These conditions imply that no mass is gained or lost at interior vertices v ∈ V 0 or across the boundary v ∈ V b of the network. To complete the definition of our model problem, we finally assume to have knowledge of the initial values u e (0) = u e,0 , c e (0) = c e,0 , e ∈ E. (2.6) Any pair of sufficiently regular functions (u, c), for instance, continuously differentiable in time and twice continuously differentiable in space on every edge, will be called a regular solution of (2.1)-(2.6) on [0, T ], if it satisfies the equations in a pointwise sense.

2.4.
Variational characterization of solutions. Throughout our analysis, we will make use of the following weak characterization of regular solutions.
for all test functions v, q ∈ H 1 (E) and all points t ∈ [0, T ] in time.
Proof. Let us start with the second identity. Multiplication of (2.2) by a test function q e on every edge e, integration over e, and summation over all edges e ∈ E leads to Via integration-by-parts on every edge e, the last term can be transformed to .
. Exchanging the order of summation and using the continuity condition (2.3), which implies that q e (v) = q(v) for some q(v) and all e ∈ E(v), the last term can be further evaluated as For the last equality, we made use of the coupling condition (2.5). A combination of the above formulas already yields the second identity of the lemma, and the first one can be derived with very similar arguments.
Remark 2.2. The equations (2.7)-(2.8) also make sense for less regular functions, e.g., The particular choice of these spaces will become clear from our analysis. Such a pair of functions (u, c) which satisfies (2.7)-(2.8) for a.a. t ∈ [0, T ] will be called a weak solution of problem (2.1)-(2.5). Note that the first term in (2.7) has to be interpreted as a duality product here. By standard embedding results [12], one can see that u, c ∈ C([0, T ]; L 2 (E)) which allows to satisfy the initial values in a reasonable way.

Well-posedness
In order to guarantee the existence and uniqueness of solutions of problem (2.1)-(2.6), we make the following assumptions on the parameters and the initial values.
and the norm of the solution can be bounded by the norm of the initial data.
A detailed proof is given in the appendix. Let us mention already here that the particular functional analytic setting allows to consider (2.1)-(2.6) as a semilinear parabolic system and the result can thus be proven by a fixed-point argument. The positivity of the initial values in assumption (A2) is not required for the local existence.
3.2. Global solutions. As a next step, we now show that the norm of the solution does not blow up in finite time and, therefore, the solution exists globally.
Here P (t) is a polynomial in t with coefficients that depend only on the bounds for the parameters in (A1), on the geometry of the graph, and on u 0 L 2 (E) and c 0 H 1 (E) .
Similarly as in [20,28], our proof is based on conservation and positivity preservation of solutions, which we state explicitly as a preparatory result.
Proof. The first assertion follows by testing (2.7) with v ≡ 1. Now let u − = min(u, 0) denote the negative part of u. Testing equation (2.7) with v = u − , we conclude that 1 2 From Theorem 3.1, we deduce that χ ∂ x c(t) L 2 (E) ≤ C for a.a. 0 ≤ t ≤ T . Using Lemma A.1, we thus obtain , where C only depends on and C G . Choosing = α/(2C) allows to absorb the first term in the left hand side of the energy estimate, and by Lemma A.2, we deduce that in the last identity we used that u(0) = u 0 ≥ 0. This shows that u(t) ≥ 0 on its domain of existence. The non-negativity of c can then be derived with similar arguments.
Proof of Theorem 3.2. The positivity of the solution is already guaranteed by Lemma 3.3. The remaining proof of the a-priori estimate then proceeds in several steps.
Step 1: As a direct consequence of Lemma 3.3, we obtain that for all 0 ≤ t ≤ T . Using this identity, we further deduce from (2.8), by testing with q = 1 and integration over time, that Note that P 1 (t) is a polynomial in t whose coefficients depend continuously on the data.
Step 2: Testing (2.8) with q = c(t) yields 1 2 Here we used (3.1) and Lemma A.1 for the second estimate, and employed (3.2) and Young's inequality for the third estimate. Note that P 2 (t) is again a polynomial of t with coefficients depending continuously on the problem data. By some elementary computations and integration with respect to time, we further obtain with polynomial P 3 (t) depending only on P 2 (t) and β.
Step 3: Testing equation (2.8) with q = ∂ t c(t) yields the estimate with C 2 = C 2 (C G , M, δ) and C 3 = C 3 (C G , M, δ). Here we used Lemma A.1 and (3.1) in the second estimate, and Young's inequality for the third. By integration in time, we get with constant C 4 = C 4 (C G , M, β, β, δ) and polynomial P 4 (t) whose coefficients again depend continuously on the data. By squaring the previous estimate, we further get with P 5 (t) = 2P 2 4 (t) and C 5 = 2C 2 4 . In the derivation of this estimate, we used Hölder's inequality to bound the term Step 4: Testing equation (2.7) with v = u(t) leads to Together with the estimate (3.3) and using Young's inequality, one can then see that The coefficients of the polynomial P 6 (t) again depend continuously on the problem data.
Step 5: Inserting the last expression in the first estimate of Step 4, slightly rearranging the terms, and integrating with respect to time now yields with polynomial P 7 (t) whose coefficients depend continuously on the data. This is the required estimate for u. A combination with (3.3) and (3.4) yields the bounds for c.

3.3.
Higher regularity. We now state some bounds for the solution in stronger norms which are obtained under the assumption that the initial values have higher regularity and satisfy the usual compatibility conditions. We thus assume in the following that : w| e ∈ H k (e) ∀e ∈ E} the space of piecewise smooth functions with appropriate regularity on the individual edges e ∈ E. Under these assumptions, one can show that the solution (u, c) in fact enjoys higher regularity.

8)
A complete proof is presented in the appendix. Let us note that this is exactly the additional regularity that can be expected, see [12,Sec 7.1.3], and which allows us to establish order optimal convergence rates for the numerical approximation in Section 6.

Discretization
We now turn to the systematic discretization of the chemotaxis problem on networks by a finite element method in space and an implicit time stepping scheme.
and P k (T h (e)) = {v ∈ L 2 (e) : v| T ∈ P k (T ), T ∈ T h (e)} the spaces of piecewise polynomials on T h (E) and T h (e), respectively, and by P k (T ) the space of polynomials of degree ≤ k on the element T . Note that P k (T h (e)) ⊂ L 2 (e), but in general P k (T h (e)) ⊂ H 1 (e).
For the approximation of the population density u and the concentration c in space, we consider the finite element space of continuous and piecewise linear functions over the mesh T h (E). We further denote by Recall that π h is uniformly bounded on L 2 (E) and on H 1 (E). Moreover, These estimates are readily proven by the standard approximation error estimates for the L 2 -projection, inverse inequalities, and interpolation arguments; see [4] for details.
Let us now turn to the time discretization. We choose a time step size τ > 0 and set t n = nτ for n ≥ 0. For any given sequence {a n } n≥0 , we denote by d τ a n = 1 τ (a n − a n−1 ) (4.3) the backward difference quotient with respect to the given time discretization. Let us note that for any scalar product ·, · with associated norm · , there holds 1 2 d τ a n 2 = 1 2τ a n 2 − 1 2τ a n−1 2 = d τ a n , a n − τ 2 d τ a n 2 , which can be verified by some basic calculations. As a direct consequence, one has 1 2τ a n 2 ≤ 1 2τ a n−1 2 + d τ a n , a n , (4.4) which will be frequently used to derive discrete energy estimates in our analysis below.

4.2.
Auxiliary results. In order to guarantee the positivity of solutions also on the discrete level, some modifications of the standard Galerkin approach will be required. In the following, we introduce the main ingredients and present some basic results.
Quasi-interpolation. For the approximation of the initial values, we will use a quasiinterpolation operator π h : where [4,9] for details. Moreover, the operator is positivity-preserving, i.e., which follows directly from the construction (4.5) via local averaging.
Mass-lumping. It is well-known that some sort of mass lumping is required to ensure a discrete maximum principle for the finite element approximation of parabolic problems, see e.g. [33]. To this end, we define for u, v ∈ H 1 pw (E) the lumped scalar product where I T (w) ∈ P 1 (T ) denotes the linear interpolation of the function w ∈ H 1 (T ). This corresponds to using numerical integration on every element T ∈ T h by the trapezoidal rule. One can verify by simple calculations that the induced norm Upwinding. As a final ingredient for our discretization scheme, we now describe the upwind technique for the convective term. Let χ∂ x c n−1 Note that for any n ≥ 1, the operator π n−1 h : Moreover, the piecewise constant function π n−1 h u interpolates u at the vertices x j and, therefore, standard Taylor estimates yield for all u ∈ H 1 (E) with a uniform constant C independent of the mesh.

4.3.
Definition of the discretization scheme. We are now in the position to formulate our numerical approximation scheme for the weak formulation of problem (2.1)-(2.6).
We will show below that the discrete solution is well-defined, that the scheme is positivity preserving, and that the total population is conserved for all time. For our analysis, we will need some auxiliary results which we state next. 4.4. Algebraic properties. Let us start with summarizing some algebraic properties of the discretization scheme (4.12)-(4.13). We denote by {φ i : i = 1, . . . , N } the nodal basis of the finite element space V h defined by (4.14) Now Moreover, the discretized system (4.12)-(4.13) can be written in algebraic form as with system matrices defined by From the definition of the matrices, one can immediately deduce the following properties.
Note that the entries of all system matrices are defined by integrals over the network E which can be split into integrals over individual elements T ∈ T h (E) respectively.
The properties in (a) then follow directly from (4.14). (b) In a similar manner, we can decompose Let u n , c n denote the coordinate vectors for the functions u n h and c n h . Using (A2), (4.7), (4.15), and the definition of the initial values, one can see that u 0 ≥ 0 and c 0 ≥ 0. Moreover, the algebraic system defining the solution at iteration n can be rewritten as deduce that M (1) and M (δ) have positive entries and that the system matrices S u and S c are M-matrices; see [29] for details. By the inverse positivity of M-matrices one can thus infer that u n ≥ 0 and c n ≥ 0, if u n−1 ≥ 0 and c n−1 ≥ 0. Positivity of the discrete solution (u n h , c n h ) for all n ≥ 0 then follows by induction and noting that u n ≥ 0, c n ≥ 0; see above. The conservation property finally follows by testing (4.12) with v h = 1.
Remark 4.4. A related algebraic argument was used by Saito in [31] to establish positivity of the discrete solution. In that work, however, a strong smallness condition on the time step size τ was required to ensure that the system matrices are row-wise diagonally dominant. Here we use the fact that the matrices S u and S c are column-wise diagonally dominant for any τ > 0 which, together with the sign conditions on the entries, is sufficient to prove the M-matrix property. This allows us to avoid the strong smallness condition on τ . The same argument should allow to improve the analysis for the method of [31]. 4.6. Uniform a-priori bounds. With similar reasoning as on the continuous level, we are now able to derive uniform bounds also for the discrete solutions.
with a polynomial P (t) whose coefficients only depend on the problem data.
Proof. The proof of Theorem 3.2 applies almost verbatim. In particular, the constants and polynomial functions appearing in the proof are of the same form as those appearing in Theorem 3.2. For convenience of the reader, we repeat the main steps. All estimates will hold uniformly for all time steps t n ≤ T .
Step 1: As a direct consequence of Lemma 4.3, we obtain Without loss of generality we may set M h = M by altering either of the two constants in the respective proofs. Testing (4.13) with q h = 1 and using the positivity of c n h , we get . From the uniform bound on u n h L 1 (E) and induction, we thus deduce that c n h L 1 (E) ≤ c 0 h L 1 (E) + t n δM h =: P 1 (t n ). Note that the polynomial P 1 (t) has the same form as that in the proof of Theorem 3.2.
Step 2: Testing (4.13) with q h = c n h and using (4.4) leads to , with the same polynomial P 2 (t) as in the proof of Theorem 3.2. Summation over n yields Step 3: By testing equation (4.13) with q h = d τ c n h and proceeding with the same arguments as in the proof of Theorem 3.2, one arrives at Squaring this estimate and some elementary estimates further lead to The constant C 5 and the polynomial P 5 are again of the same form as those in Theorem 3.2.
Step 4: We now turn to the estimates for u n h . By testing (4.12) with v h = u n h , using (4.4), and proceeding in the same manner as in the proof of Theorem 3.2, we get From the previous estimates for c n h , we may then further deduce that and this bound holds with the same polynomial P 6 (t) as in the proof of Theorem 3.2.
Step 5: A combination of the estimates in Step 4 now leads to which yields the desired bound for u n h . The corresponding estimates for c n h then follow by inserting this bound into the estimates of Step 2 and Step 3.
Let us emphasize that all arguments used for the analysis of the problem on the continuous level carry over to the discrete setting almost verbatim. The reason for this is that by its variational character, the proposed method with mass lumping, upwinding, and implicit time integration inherits all important structures of the continuous problem.

Convergence
Based on the uniform bounds of the discrete solution provided by Theorem 4.5, we now establish the convergence of the discretization scheme towards the unique solution of the continuous problem without additional regularity assumptions. For a convenient presentation of our results, we will interpret the discrete functions (w n h ) n≥0 as (piecewise) continuous functions of time. Such extensions will be denoted with double subscripts w h,τ . Proof. The proof is based on standard arguments for the analysis of nonlinear parabolic problems, see e.g. [30]. We therefore only sketch the main arguments.
Step 1a: Let (u h,τ , c h,τ ) be piecewise constant in time with values given by u h,τ (t) = u n h and c h,τ (t) = c n h for t n−1 < t ≤ t n . Furthermore, let c h,τ and π h,τ be piecewise constant in time with values c h,τ (t) = c n−1 h and π h,τ (t) = π n−1 h for t n−1 < t ≤ t n . Then from the variational characterization (4.12)-(4.13), one can deduce that for all test functions v h ∈ V h and q h ∈ V h , and for a.a. 0 ≤ t ≤ T .
Step 2a: Recall that π h : , the approximation properties of the L 2 -projection π h , and the quasi-uniformity of the mesh T h (E), one can see that Then by the uniform bounds for the discrete solution (u n h , c n h ) n≥0 and the quasi-uniformity of the mesh T h (E), one can see that Recall that H k pw = {v ∈ L 2 (E) : v| e ∈ H k (e)} is the space of piecewise smooth functions. Thus, the family {L h,τ } h,τ >0 of linear operators is uniformly bounded on H 1 (E) and L h,τ v → 0 for v ∈ H 2 pw (E) ∩ H 1 (E) which is dense in H 1 (E). Hence, L h,τ v → 0 for all v ∈ H 1 (E) by the Banach-Steinhaus theorem [16]. In summary, we thus have shown that for all 0 ≤ t < t ≤ T and all v ∈ H 1 (E) as h, τ → 0. It then follows from Lebesgue's differentiation theorem [12] that Step 2b: In a similar manner, one can show that for all v ∈ H 1 (E) and for a.a. 0 ≤ t ≤ T . This shows that (u * , c * ) satisfies equation (2.8).
We next turn to the verification of identity (2.7).
Step 3a: With the very same arguments as above, one can show that for all v ∈ H 1 (E) and for a.a. 0 ≤ t ≤ T as h, τ → 0. It thus remains to establish the convergence for the convective term, which we do next.
Step 3b: Let us define linear operators B h,τ : Then from the uniform bounds for the discrete solution (u n h , c n h ) n≥0 , the H 1 -stability of the L 2 -projection π h on the quasi-uniform mesh T h (E), and (4.10), one can deduce that This shows that the family of operators {B h,τ } h,τ >0 is uniformly bounded for all h, τ > 0. As a next step, we decompose By the Aubin-Lions lemma [30], we know that u h,τ → u * in L 2 (0, T ; L 2 (E)). Moreover, for all 0 ≤ t ≤ t ≤ T and for any test function v ∈ H 1 (E). From the estimate (4.11), we therefore deduce that π h,τ u h,τ − u h,τ L 2 (t ,t ;L ∞ (E)) ≤ Ch 1/2 ∂ x u h,τ L 2 (t ,t ;L 2 (E)) . Using the H 1 -stability of the L 2 -projection π h and the uniform a-priori bounds for the discrete solution (u n h , c n h ) n≥0 , one can then see that In summary, we thus have shown that With the assertions of Step 3a, this shows that (u * , c * ) solves (2.7).
Step 4: From the estimates for the quasi-interpolation operator π h stated in Section 2, one can deduce that π h v → v in L 2 (E) for all v ∈ L 2 (E). Together with the continuity of the trace mapping for the space W (0, T ), this yields u * (0) = u 0 and c * (0) = c 0 .
Step 5: In summary, we have shown that (u * , c * ) is a weak solution of (2.1)-(2.6), and from the uniqueness stated in Theorem 3.2, we infer that u * = u and c * = c.

Error estimates
Under suitable smoothness assumptions on the true solution, we can now also derive quantitative convergence rates. The aim of this section is to prove the following result.
with constant C that only depends on the bounds for the coefficients, the time horizon T , the geometry of the network, and the norm of the solution (u, c) in (3.7)-(3.8).
In the usual way [34,35], we decompose the errors in the two solution components via into approximation and discrete error components. For our analysis, we choose the operator R h : In the following two sections, we derive bounds for the two error contributions of the individual solution components, and we then complete the proof of the above theorem.
6.1. Approximation error. We start with summarizing some elementary properties of the H 1 -projection R h : H 1 (E) → V h defined above. These are well-known for a single edge e, see e.g. [4], and can be generalized easily to the network setting.
Here ∂ xx u L 2 (E) = ( T ∂ xx u 2 L 2 (T ) ) 1/2 denotes the norm of the broken derivative ∂ xx u. As a direct consequence of these estimates, we obtain the following bounds for the approximation error contributions in the above error splitting.
For the proof of the corresponding estimates for the discrete error components we require a number of auxiliary results which are stated and proved in the next section.
6.2. Auxiliary results. As a preliminary step, we now state estimates for some terms that will arise in the analysis of the discrete error components e n h , d n h below. Lemma 6.4. Let u h , v h ∈ V h , then for any > 0, Proof. Let us note that the numerical integration is exact, if the product u h v h is piecewise linear. By the Bramble-Hilbert lemma and scaling arguments, one can then see that Via an inverse inequality, the second term can be bounded by ∂ The result then follows by combination of the last two inequalities, summation over all elements, and application of the Cauchy-Schwarz inequality.
Using the properties of the H 1 -projection R h , we can derive the following bounds.
Lemma 6.5. Let w ∈ H 1 (E). Then for any > 0, Proof. By application of Lemma 6.2, the Cauchy-Schwarz inequality, and the norm equivalence estimates (4.8), we obtain The assertion of the lemma now follows via Young's inequality.
As a next step, we derive an estimate for the errors introduced through mass lumping and the approximation of the time derivatives by finite differences. Lemma 6.6. Let w ∈ L ∞ (0, T ; L 2 (E)) ∩ L 2 (0, T ; H 1 (E)). Then for any > 0, Proof. We start with splitting the error by By Lemmas 6.4 and 6.5, we readily obtain In the second step, we used the fundamental theorem of calculus the Cauchy-Schwarz, and Young's inequality. The second term can be further estimated by Here we used Taylor expansion and the Cauchy-Schwarz inequality for the first term, the norm equivalence (4.8) for the second, and applied Young's inequality to split the product. A combination of the two estimates for (i) and (ii) yields the assertion.
As a last step in our preliminary considerations, we now derive a bound for the error introduced through the upwinding strategy. Lemma 6.7. Let c n−1 h ∈ V h be given, and let π n−1 h denote the corresponding upwind operator as defined in (4.9). Then for all > 0, Recall that e n h = R h u(t n ) − u n h and d n h = R h c(t n ) − c n h are the discrete error contributions.
Proof. By application of the Cauchy-Schwarz and Young's inequality, we obtain Using triangle inequalities, the first term can be further estimated by Before turning to the estimation of the individual terms, let us make a preliminary observation. From the embedding of H 1 (E) in L ∞ (E), the bounds for u(t) H 1 (E) in (3.7), and the estimates of Theorem 4.5, we obtain and these bounds hold uniformly for all 0 ≤ t ≤ T, n ≥ 0. With the aid of Hölder's inequality, Taylor estimates, and the Cauchy-Schwarz inequality, we can then estimate the first term in the above error expansion by For the second term, we introduce temporal averages c n = 1 τ t n t n−1 c(t)dt for the function c. Then by the Hölder and triangle inequalities, and using the bound (6.4), we get L 2 (t n−1 ,t n ;H 1 (E)) + h 2 /τ ∂ xx c 2 L 2 (t n−1 ,t n ;L 2 (E)) ) 1/2 . For the third term in the above estimate, we get Using Hölder's inequality, the uniform bounds for ∂ x c n h L 2 (E) provided by Theorem 4.5, and similar reasoning as in the estimate of the term (ii), we obtain L 2 (t n−1 ,t n ;H 1 (E)) + h 3 /τ ∂ xx u 2 L 2 (t n−1 ,t n ;L 2 (E)) ) 1/2 . In the last step, we employed Taylor estimates, the Cauchy-Schwarz inequality, and the properties of the H 1 -projection. By the estimate (4.11) for the upwind projection and similar arguments as before, the fifth term can be estimated by L 2 (t n−1 ,t n ;H 1 (E)) + h 2 /τ ∂ x u 2 L 2 (t n−1 ,t n ;H 1 pw (E)) ) 1/2 .
Using the interpolation inequality (A.2) and the uniform bounds (6.4) for the discrete solution, the last term can finally be bounded by h ) 1/2 . The assertion of the lemma now follows by squaring the estimates of the terms (i)-(vi), and combination with the first inequality of the proof. 6.3. Estimates for the discrete error. We are now in the position to derive the following bounds for the discrete error contributions.
The constants only depend on the parameters of the problem, the geometry of the graph, and the norms of the solution components appearing in (3.7)-(3.8).
Proof. From the variational characterization of u and u n h , one can see that As in the proof of the preceding lemma, we have for all > 0, C τ ∂ t u 2 L 2 (t n−1 ,t n ;H 1 (E)) + h 2 /τ ∂ xx u 2 L 2 (t n−1 ,t n ;L 2 (E)) .
The remaining terms on the right hand side can be estimated by the previous lemmas. Testing with v h = e n h and using some elementary computations, one thus obtains Choosing = α/2 allows to absorb the first term on the right hand side by the left hand side of the inequality. Multiplication by 2τ , summation over n, and using the bounds From the variational equations characterizing c and c n h , we obtain that As before, we denote by c n = 1 τ t n t n−1 c(t)dt the temporal averages of the function c on the subinterval [t n−1 , t n ). We then have for all > 0, Here, we have used the equivalence of the norms (4.8) for the first term, Lemma 6.4 for the second term, as well as Taylor estimates and the properties of the H 1 -projector. In the same way, we obtain for all > 0, The remaining terms on the right hand side can now be estimated by the results of the previous section. Choosing q h = d n h and some elementary computations, we thus arrive at L 2 (t n−1 ,t n ;L 2 (E)) + h 2 τ ∂ xx c 2 L 2 (t n−1 ,t n ;L 2 (E)) + ∂ x c 2 L 2 (t n−1 ,t n ;L 2 (E)) + ∂ x u 2 L 2 (t n−1 ,t n ;L 2 (E)) . With = β/2, the first term on the right hand side can be absorbed in the left hand side. Multiplying by 2τ , summing over n, and using the bounds (3.7)-(3.8), we obtain Inserting this estimate into (6.5) and applying the discrete Gronwall lemma (A.3) once more, we can conclude that This is the required estimate for the first component of the discrete error. Using this estimate in (6.6) yields the bound for the second component of the discrete error.
6.4. Proof of Theorem 6.1. The error splitting (6.1)-(6.2), we directly obtain The assertion of the theorem now follows by simply estimating the individual terms with the help of the bounds provided by Lemma 6.3 and Lemma 6.8

Numerical illustration
In this section, we illustrate the theoretical results of the paper by some numerical tests.  In order to verify the theoretical convergence rates obtained in section 6, we proceed as follows: Since no analytical solution is available for this test case, we use the difference between the numerical solutions on two different meshes with meshsize H and h = H/2 as an approximation for the actual error. Since the method has been proven to converge with a given rate, this yields accurate approximations for the true error. The results obtained in our numerical tests are presented in Table 7  As predicted by our theoretical results in Section 6, we observe first order of convergence in both solution components and norms used for our analysis. 7.2. Block network. As a second test case, we consider an example proposed in [2], namely the block network depicted in Figure 7 The initial values for the two solution components are simply chosen as In order to get an interesting behavior, the boundary conditions at the two ports of the network are chose as follows: .
This means that the bacteria with density u enter the network at node v 0 and they are expected to move towards v 26 where the chemoattractant c is added. Due to the different properties of the individual pipes, we expect the bacteria to move along the red path  As can be seen from the images, both solution components are smooth and the system behavior seems to be computed correctly as expected.

Summary and discussion
In this paper, we considered a chemotaxis model on networks described by a system of partial differential-algebraic equations, which can be seen as a natural generalization of the minimal model to the network context. By extending the results of Osaki and Yagi, we were able to establish global existence and uniqueness of solutions based on perturbation arguments for semilinear evolution problems, conservation of mass, and positivity of the solutions. In addition, we derived regularity estimates for the solutions.
The arguments of the proofs were developed in a way that allowed us to prove the global existence and uniqueness of solutions also for numerical approximations obtained by a finite element discretization with mass lumping and upwinding and an implicit Euler method for time integration. The discrete approximations could be shown to converge to the unique global solution of the chemotaxis problem without artificial smoothness requirements on the solution. In addition, we could establish order optimal convergence rates under minimal smoothness assumptions.
The arguments used for the analysis of the finite element method presented in this paper may be used to improve also the convergence results of the method of Saito [31], in particular, to get rid of the strong stepsize restrictions needed there. Also a generalization to chemotaxis models with nonlinear coefficients seems possible to some extent. A class of problem that would certainly deserve further considerations, also from a numerical point of view, are models of haptotaxis, where no diffusion is present in the equation governing the chemoattractant.
The following embedding inequalities are used several times in our proofs.
Lemma A.1. Let (V, E, ) be a finite metric graph. Then for any > 0, for all f ∈ H 1 (E) with a constant C G only depending on the geometry of the graph.
Proof. The assertions follow by applying the classical Gagliardo-Nirenberg inequality on every edge, summation over all edges, using the equivalence of the l p -norms on finite sequences, and application of Young's inequality.
We also require the following continuous and discrete version of Gronwall's inequality.
Further suppose that τ γ n < 1 for n = 1, . . . , N and set σ n := (1 − τ γ n ) −1 . Then Proof. Existence of a unique weak solution follows by Galerkin approximation and standard arguments; see [10,12] for details. To keep track of the constants, we give a short proof of the a-priori estimates. Testing the variational equation (B.2) with q = c(t) yields 1 2 Here we used the bounds for the coefficients for the first step and Young's inequality with > 0 for the second. An application of Lemma A.2 with = 1/t further yields ) . The first estimate then follows by noting that z 2 L 2 (0,t;L 2 (E)) ≤ t z 2 L ∞ (0,t;L 2 (E)) , which follows by the Cauchy-Schwarz inequality, and some further elementary computations. For the second bound of the lemma, we test (B.2) with q = ∂ t c(t) which yields . By rearranging the terms, integration in time, and using assumption (A1), we obtain . This yields the second estimate and completes the proof of the lemma.
As a next step, we now derive a-priori estimates for the second solution component.
Proof. Testing (B.1) with v = u(t), we obtain via assumption (A1) that 1 2 , where we used Lemma A.1 to estimate u(t) L ∞ (E) . Choosing = α/(4χC) and applying Young's inequality yields with a constant C = C(χ, α, C, C G ). Inserting this expression in the above estimate, rearranging some of the terms, and applying Lemma A.2 further yields The result then follows by taking the maximum over t ∈ [0, T ] on both sides.
By combination of the two previous results, we directly obtain the following assertion.
Based on the previous results, we can define a mapping Lemma B.4. For any T > 0, u 0 ∈ H 1 (E), and c 0 ∈ L 2 (E), the mapping Φ T is well defined on X T = L ∞ (0, T ; L 2 (E)). Moreover, for any R > u 0 L 2 (E) there exists T (R) > 0 such that for all 0 < T ≤ T (R), Φ T maps B T,R = {z ∈ X T : z X T ≤ R} into itself.
Proof. The assertion follows directly from the previous lemmas.
With similar arguments as employed for the proof of the assertions stated in the previous section, we can show that Φ T is Lipschitz continuous on X T = L ∞ (0, T ; L 2 (E)).
Lemma B.5. Let (A1) hold and let u 0 , c 0 and T (R) be given as in Lemma B.4. Then for any 0 < T ≤ T (R), we have with Lipschitz constant L(T ) = C 5 T e C 6 T and constants C 5 , C 6 independent of T .
Proof. Let (u, c) and ( u, c) denote the solutions of (B.1)-(B.2) with the same initial values u 0 ∈ L 2 (E) and c 0 ∈ H 1 (E) but with different data z and z ∈ B T,R . Then ∂ t c(t) − ∂ t c(t), q E + β∂ x c(t) − β∂ x c(t), ∂ x q E + γc(t) − γ c(t), q E = δz(t) − δ z(t), q E for all suitable test functions q and 0 ≤ t ≤ T . In addition, c(0) − c(0) = 0. With similar arguments as in the proof of Lemma B.1, one can see that Next observe that u − u satisfies u(0) − u(0) = 0 and, in addition, there holds ∂ t u(t) − u(t), v E + α∂ x u(t) − α u(t), ∂ x v E = χ(u(t) − u(t))∂ x c(t), ∂ x v E + χ u(t)(∂ x c(t) − ∂ x c(t)), ∂ x v E for any suitable test function v and a.a. 0 ≤ t ≤ T . By choosing v = u(t) − u(t) and applying some elementary manipulations, one can deduce that 1 2 With similar arguments as were used to bound the term ( * ) in the proof of Lemma B.2, the first term can be further estimated by with C depending on the graph, on the bounds in assumption (A1), on the norms c 0 H 1 (E) and u 0 L 2 (E) of the initial values, as well as on R and T (R) as required. Using Lemma B.2, Lemma A.1, and (B.4), the second term can be estimated by with C depending on the problem data like C as required. A combination of these estimates and an application of Lemma A.2 finally yields the result.
Problem (2.1)-(2.6) is equivalent to the fixed-point problem u = Φ T (u). As a direct consequence of the previous Lemmas, one can see that for all 0 < T ≤ T (R) with T (R) chosen sufficiently small, Φ T maps B R,T = {z ∈ L ∞ (0, T ; L 2 (E)) : z L ∞ (0,T ;L 2 (E)) ≤ R} into itself and is a contraction. Hence, Banach's fixed-point theorem guarantees the existence of a unique fixed-point u ∈ B R,T with u = Φ T (u). An application of Lemma B.3 then yields the assertion of the theorem.
Appendix C. Proof of Theorem 3.4 By formally differentiating the system (2.1)-(2.2) we obtain where w = ∂ t u and d = ∂ t c are the derivatives of the solution (u, c). Similarly as in the previous section, we will prove existence of local solutions via Banach's fixed-point theorem. To this end, we consider the following linearized system for all v, q ∈ H 1 (E) and 0 ≤ t ≤ T with z ∈ L ∞ (0, T ; L 2 (E)) given.
Proof. It suffices to prove the a-priori estimate. Existence of a unique solution can then be obtained by Galerkin approximation. According to Lemma B.1, we have ∂ x d L ∞ (0,t;L 2 (E)) ≤ C( d 0 H 1 (E) + t 1/2 z L ∞ (0,t;L 2 (E)) ) for a.a. t ∈ [0, T ]. Testing (C.3) with v = w(t) we obtain the differential inequality Integrating between 0 and t, using that c ∈ L ∞ (0, T ; H 1 (E)), u ∈ L 2 (0, T ; H 1 (E)), applying the interpolation inequality (A.2) and Gronwall's lemma, we arrive at The bound of the lemma is then obtained by taking the supremum over t on both sides.
Similarly to the previous section, we can now define the fixed-point map where X T = L ∞ (0, T ; L 2 (E)) is chosen as before and where w denotes the first component of the solution of system (C.3)-(C.4) with given initial values w 0 ∈ L 2 (E) and d 0 ∈ H 1 (E). By the previous lemma, Ψ T is well-defined and maps X T into itself. As a next step, we verify that Ψ T is a contraction, if T is chosen sufficiently small.
Lemma C.2. Let the assumptions of Lemma C.1 be valid. Then with Lipschitz constant L(T ) = C P (T )T 1/2 e CQ(T ) , where C , C are constants and P (T ), Q(T ) are polynomials of T that only depend on the problem data and the geometry of the graph.
Proof. Let (w, d) and ( w, d) be two solutions of (C.3)-(C.4) with the same initial values w 0 ∈ L 2 (E) and d 0 ∈ H 1 (E), but with different data z, z ∈ X T . With the same arguments as in the proof of Lemma B.5, we obtain Moreover, the difference w − w satisfies w(0) − w(0) = 0 and By choosing v = w(t) − w(t), one can then deduce that