Nonlinear flux-limited models for chemotaxis on networks

In this paper we consider macroscopic nonlinear moment models for the approximation of kinetic chemotaxis equations on a network. Coupling conditions at the nodes of the network for these models are derived from the coupling conditions of kinetic equations. The results of the different models are compared and relations to a Keller-Segel model on networks are discussed. For a numerical approximation of the governing equations an asymptotic well-balanced schemes is extended to directed graphs. Kinetic and macroscopic equations are investigated numerically and their solutions are compared for tripod and more general networks.


1.
Introduction. The Keller-Segel equations have been widely used to model chemotaxis. These equations and in particular, the properties of their solutions have been intensively investigated, see for example [9,26,32,36]. For adapted and expanded versions of the Keller-Segel model we refer to [7,13,14,29,32,33,12]. Moreover, in recent literature improved flux-limited Keller-Segel models, taking into account the finite speed of propagation, have been developed in [1]. For surveys and extended reference lists, see, for example [2,3].
Our starting point is the classical kinetic chemotaxis equation [11]. Scaling it with the so called diffusive scaling leads to the Keller-Segel equation, see [11]. In general, the derivation of Keller-Segel type models, including flux-limited diffusion models and Fokker-Planck type models, from underlying kinetic or microscopic models is discussed for example in [2,12,11]. In particular, using moment closure approaches one may obtain macroscopic equations intermediate between kinetic and Keller-Segel equations, see the above mentioned references or [19,28,6]. Linear moment closure models called the P1-model (in the radiative transfer literature) or the Cattaneo equations are in many situations resonable approximations of the kinetic equation [6]. However, a major drawback of these approximations is that they do not guarantee the positivity of the cell density. Thus, in this paper we consider a non-linear closure method to deal with this problem. Additionally, we also introduce half-moment non-linear closures for the kinetic equations and obtain associated hydrodynamic macroscopic equations similar to [17].
To model chemotaxis on a network, the crucial point is to define suitable coupling conditions. In previous work [5,6] coupling conditions for the scalar Keller-Segel equation and kinetic equation have been discussed. An analytical investigation can be found in [10]. These investigations deal with the case of a network problem for a parabolic system of equations. Conservation of mass and continuity of density and chemo-attractant at the nodes are used as coupling conditions. Moreover, in [8,25] the hyperbolic-parabolic Cattaneo equations on a graph are studied. For general work on kinetic equations on networks, for example for traffic flow equations, we refer to [18] or [27].
In the present work we will derive the coupling conditions for nonlinear halfmoment models from the coupling condition for the kinetic model. These coupling conditions guarantee on the one hand the conservation of mass through nodes, and on the other hand, they satisfy a positivity condition for the density. In the present case, as in [8,25], the hyperbolic equations require another approach to the coupling conditions as in the Keller-Segel context. However, in the diffusive limit, when the scaling parameter goes to 0, all coupling conditions converge to those of the Keller-Segel model, i.e. the conservation of mass through nodes and a continuity condition, compare for example [5].
Numerical methods for these nonlinear models have to deal with issues like the numerical moment realizability and (as in the linear case) with a large range of the scaling parameter describing different regimes. Thus, it is desirable to use an asymptotic preserving well-balanced (APWB) scheme. These numerical methods work uniformly with respect to the scaling parameter and one obtains in particular a scheme suitable for computations near the diffusive limit. We refer to [22] for the basic idea of such schemes. Combining this with a numerical scheme to solve the coupling problem we obtain results for the full network problem.
The paper is organized as follows. In section 2 we discuss the kinetic model and its nonlinear hydrodynamic full-and half-moment approximations which are derived using non-linear closure function. Section 3 considers the coupling conditions for the kinetic equations and derived conditions for the macroscopic models. Section 4 and 5 contain details about the numerical schemes on an interval and the numerical treatment of the coupling conditions. Finally, the numerical results for tripod and more general networks are shown in Section 6.
2. Kinetic and macroscopic models for chemotaxis.
is the density of the chemoattractant. Let α ∈ R + denote the rate of turning into the favourable direction given by the chemoattractant, let λ ∈ R + be the turning rate of the cells and D, γ ρ , γ m ∈ R + be the diffusivity, production and decay rate of the chemoattractant, respectively. As in [6], we consider the kinetic equation Here, we have defined We assume λ ≥ α in order to guarantee positivity of the turning kernel for all v, compare [11]. We introduce a parameter ∈ [0, ∞) and scale the kinetic equation (1) in the following way   This equation will be used as the starting point for the macroscopic models considered in the next section. The kinetic equation (3) can be transformed using the even-and odd-parities [30] r for positive velocities (v ≥ 0). This leads to the system When → 0, the limit of the kinetic equation (3) is the Keller-Segel equation with flux-limited chemoattractant. The procedure is shortly revisited. Start with the scaled version of the kinetic equation (6). Integrating the equation for r over v ≥ 0 we obtain Moreover we observe r = ρ 2 + O( 2 ) and j = α∂xm . Plugging the expansions into (8) gives the diffusion limit of the kinetic equation, the modified or flux-limited Keller-Segel equation In addition, from the scaled version for kinetic chemotaxis equation (3), using moment closure approaches one may obtain macroscopic, hyperbolic equations intermediate between kinetic and Keller-Segel equations, see [19,29]. There are two kinds of moment models for chemotaxis, full-moment models and half-moment models, that are corresponding to full-and half-moment closure methods respectively.

2.2.
Full-moment non-linear models. Consider the following averaged quantities where ρ represents the density and q the flow of the cells. The fact that f (x, t, v) ≥ 0, it leads to the realizability conditions for full-moment models, see [16,37]: This leads to the equations In order to close the model (11), a moment closure is specified by using a suitable approximation of f in terms of a function of ρ and q to obtain expressions for vq(x, t) leads to the P1 model for chemotaxis as in [8,25,6]. However this approach does not guarantee the realizability condition (10) and the cell density f can be negative. Therefore, one may assume that the cell density f (x, t, v) is approximated by an exponential function f (x, t, v) = a exp(vb) (12) with a ∈ R + and b ∈ R. This closure function is chosen using an entropy minimization principle, see for example [23] . This leads to guaranteeing the realizability condition. The function h(u) is called the Eddington factor. Without refering to an underlying closure function the function h(u) can also be chosen directly such that the realizability condition is fulfilled. In [35,23] models with different forms for this function have been proposed, for example, the Kershaw closure h(u) = 1 3 + 2 3 u 2 or the Levermore-Lorentz Eddington factor The resulting full-moment macroscopic model for chemotaxis derived by the nonlinear closure (12) reads

M1 MODEL ON NETWORKS 385
Scaling and denoting q := q it can be rewritten as Under the realizability condition, the model (14) is hyperbolic and converges to the Keller-Segel equations when → 0.

2.3.
Half-moment non-linear models. We will construct a macroscopic, halfmoment model for (3) by applying a half-moment nonlinear closure, which was introduced in [17], see [6] for a linear version. Consider the following quantities Using These conditions again guarantee the existence of a non-negative density for these equations [16,37].
Integrating the kinetic equation (3) leads to the half-moment model The nonlinear closure for the half-moment models is derived similarly as for the full-moment case. One uses the ansatz Plugging (18) in (15) leads to Condition (19) guarantees the realizability condition (16). We note that an explicit maximum-entropy Eddington factor approximating the above closure is obtained by using the Kershaw closure [34] h Finally, the half-moment model for chemotaxis reads the half-moment model can be rewritten as with P (ρ, q ,ρ,q) : The half-moment model (23) is hyperbolic and has again the Keller-Segel equations as macroscopic diffusive limit when goes to 0.
3. Asymptotic preserving well-balanced (APWB) schemes for the nonlinear models on an interval. We concentrate in the following on the discretization of the equations for the movements of the cell density. The equation for the chemoattractant m is always solved with a forward central difference scheme as in [5]. We construct an APWB numerical scheme for the nonlinear chemotaxis model following closely the strategies developed in [21] for the linear Cattaneo model and [4,15] for radiation hydrodynamics. Moreover, we refer to [24,20,31] for further details on AP and well-balanced schemes.
3.1. Framework for APWB scheme. First, let us recall the framework of the APWB scheme in [24,20,21] for the Cattaneo model By introducing the diagonal variables (25) is rewritten as Using the well-balanced Godunov scheme for (26) is given as with l := √ a ∆t ∆x and ∂ x m n The APWB numerical solution (28) can be rewritten explicitly as It is easy to check that for small enough values of ≤ √ a ∆t ∆x ≤ 1, the APWB scheme (29) is positivity preserving for U and V , if the following parabolic CFL restriction holds ∆t ≤ λ∆x 2 4 (a + ∆x) .
3.2. APWB scheme for the nonlinear full-moment model. Based on the above considerations and the schemes developed in the above cited references we consider the following relaxation system for (14) For β → 0 the system converges to (14). Note that for any value of β the system (30) converges to the Keller-Segel equation with a diffusion coefficient equal to a. System (30) is split into two into separate sub systems. The right hand side treated with a simple implicit numerical scheme and we obtain for β = 0 which is used as the initial condition for the system given by the left hand side. By introducing the new variables the left hand side of system (30) can be rewritten as Now, we can apply the APWB schemes presented in 3.1 for each subsystem in (33).
Remark that for a := max h q ρ the quantities in (32) are nonnegative and the scheme converges formally to the Keller-Segel equations as → 0.
Remark 1. We note that under similar conditions as in the linear case, it can be shown, that the scheme preserves positivity in U and V . However, this does not mean that we can prove that the density ρ remains positive. In the linear case, this is not an issue to be investigated, since the density ρ in the Cattaneo model is not necessarily positive. In contrast, in the nonlinear case, it is an important open issue to be investigated.
3.3. APWB scheme for half-moment M1 model. Following the same procedure as for the full-moment model in the previous section, we consider the following relaxation system for the half-moment model (23) with h ± := h ± q ±q ρ± ρ and P = P (ρ, q ,ρ,q) ,P = P (ρ, q ,ρ,q) defined in section 2.3. To obtain the correct asymptotic behaviour using a splitting scheme, we introduce a parameter σ < λ, and rewrite the system (34) as Using again an implicit scheme for the right hand side of (35) and setting β = 0, one obtains for the first step of the splitting scheme This relaxation scheme ensures that system (35) has the correct asymptotic limit as → 0. Thus (35) converges to the Keller-Segel equation (9). Moreover, the solution in the relaxation part (36) preserves the realizability condition i.e., | q (1) | ≤ ρ (1) .
The transport part in (35) includes four independent subsystems with (36) as the initial data. Introducing the characteristic variables ( we use the APWB scheme for the left hand side of (35). With a := max w ρ we have again the positivity for the characteristic variables (37). Moreover, the half-moment model (34) converges formally to the Keller-Segel equation (9) when β → 0 and → 0.

4.
Coupling conditions for chemotaxis models on networks. In this section we state a hierarchy of coupling conditions for the kinetic and the nonlinear halfmoment model proposed in the previous section. In [6] coupling conditions have been derived for the kinetic model and linear macroscopic approximations. These conditions can be used as well for the nonlinear half moment model discussed in the previous section. Finding conditions for the nonlinear full moment model is much harder since the signs of the eigen values can change with the states of the system. First the coupling conditions for the kinetic model are recalled to derive thereof the nonlinear half moment coupling conditions.

4.1.
Coupling condition for kinetic equations [6]. Consider the kinetic equation As stated in [6], the coupling conditions for (38) should assign a value to all f (v) for v ∈ [0, 1]. A possible choice is for v ∈ [0, 1] and i = 1, . . . , N . In order to conserve the total mass in the system the matrix A ∈ R N ×N has to fulfill N i=1 a i,j = 1 ∀j = 1, . . . , N .
Moreover, positivity requires a i,j ≥ 0, i, j = 1, . . . , N . If all edges are treated equally, the coupling conditions for three edges are  This can be generalized to the cases N > 3 by choosing a i,j = 1 N −1 i = j and a i,i = 0.
Remark 2. Proving existence, uniqueness or stability of the kinetic equation with the above coupling conditions on the graph are interesting further issues, which would deserve a deeper analysis.

4.2.
Coupling condition for the nonlinear half-moment model. Since the kinetic coupling conditions (39) are linear and independent of v we obtain directly by integrating the coupling conditions for the half moment model (21)  Recall that all properties of (39) are inherited by the above coupling conditions. For example the total mass in the system is conserved since In particular, the coupling conditions for the half-moment model yield the positivity of the densities, i.e. ρ + i ≥ 0 ∀i = 1, . . . , N if ρ − j ≥ 0 ∀j = 1, . . . , N .

4.3.
Numerical coupling of the nonlinear half-moment model. Here only the case of a junction coupling three edges will be treated, but all procedures easily extend to arbitrary junctions. Consider the left hand side of (35) which includes four independent subsystems of eight variables For the additional numerical quantities we need further coupling conditions. Note that for i = 1, ..., N from the relaxation part of (35) we have: z i = q i ,ẑ i =q i and w i = P i ,ŵ i =P i : Thus, conditions for the states z andẑ are found from those prescribed for q ∓ . To obtain conditions for w andŵ we have to find conditions for P ± i . These are easily derived from the kinetic coupling conditions (39). Multiply (39) by v 2 and integrate over positive velocities to obtain  Inserting (22) into (40), we express the coupling conditions in the variables of (41)  and with (42) we obtain  This leads to  Since on each edge there are four waves with positive speed, we need exactly twelve equations (43),(45) for three edges, as provided.
Remark 3. We have not investigated the positivity preserving property of the scheme for the density ρ, see Remark 1 for the case of a single line. It would be even more interesting to investigate this property on the graph taking into account the discretization of the coupling conditions. 5. Numerical tests. In this section we investigate the proposed models in several numerical test cases for different values of . For simplicity we use for the nonlinear half-moment model the explicit Kershaw closure, i.e. h ± given by (20). In all the considered examples we will use the following values for the parameters λ = α = 1, D = 1, γ ρ = 1 and γ m = 0.1. The spatial resolution is ∆x = 0.02 and the time step is chosen according to the CFL condition. At end points where no coupling conditions are imposed zero Neumann boundary conditions are applied. For the kinetic model we discretize the velocity space V = [−1, 1] with N v = 50 cells.

Numerical solutions on an interval.
Test 1. We consider the interval [0, 2] and ∆x = 0.005. As initial conditions for the kinetic equation we use The remaining initial values for the hydrodynamic equations can be derived from the kinetic initial condition. At t = 0 no chemoattractant m(x, 0) = 0 is present.
In the Figures 2 and 3 the densities at t = 0.2 for the values = 1, 0.5, 0.1 and = 10 −6 are shown. For the kinetic model we observe that diffusion depends strongly on the value of as expected. The closest approximation to the kinetic equations are the half moment models, in particular, the nonlinear half-moment model. For the half moment P 1 model for = 1 and = 0.5 one can clearly observe the four waves generated by the advective part. In the full moment P1 model only two waves are used to approximate the kinetic solution. The flux-limited Keller-Segel equation is evolving too fast for large . In general, the half-moment models are clearly superior to the full moment models and the nonlinear models provide better approximations than the linear models. For small all models converge to the solution of the Keller-Segel equation.

Test 2.
In the second test we investigate the positivity preserving of the nonlinear moment models which is not guaranteed in the case of the linear models. We use non-equilibrium initial conditions for the kinetic equation (   with The initial conditions for the full-moment models are The solutions are plotted in Figure 4, for ∆x = 0.02 at t = 0.5. Even though the initial conditions satisfy the realizability conditions ρ ± ≥ |q ± | (for half-moment models) and ρ ≥ |q| (for full-moment models), only the nonlinear full-and halfmoment models preserve the positivity of the cell density.  We study a junction connecting three outgoing edges, as depicted in Figure 1. The results of the nonlinear half-moment model with the coupling conditions derived in section 4 are compared with the results of the kinetic, linear half and full moment and Keller-Segel model with coupling conditions described in [6]. On each of the edges we consider the interval [0, 1]. As initial conditions we choose which is consistent with the following values for the kinetic equation All other quantities are initially zero.
In is again the best approximation to the kinetic result. As the value of decreases all solutions approach the solution of the Keller-Segel equation.

Numerical solutions on a larger network.
In this last test case we consider a larger network of 31 edges and 23 nodes as shown in Figure 9. The length of the short edges is 0.5, for the longer ones we have 1 and √ 2 respectively. Note that this network does not only contain three way junctions, but also nodes connecting up to five edges. At the open ends of the network Dirichlet boundary conditions are imposed. For the kinetic model we prescribe f (0, v, t) = 1 2 for v > 0. The boundary values for the other models can be derived thereof. As initial conditions all values are set to zero except the density in the edges at the outer boundaries, which is set to 1. For the numerical computation we used 15, 30 and 42 cells for the spatial discretization of the edges, respectively . The velocity space in the kinetic model is resolved with 30 points. The scaling parameter is chosen as = 1.
In Figure 9 the density at time t = 5 for all four models is shown. There is an inflow from both sides of the network. As observed in the previous tests, the states of the Keller-Segel model propagate faster, such that the network is filled at     almost coincide. Concerning the computation times, the least expensive model is the full moment P 1 model with approximately 15% of the computation time of the kinetic model. The Keller-Segel model needs 18%, the half moment P 1 model needs 22% and the half moment M 1 model needs 23% of the kinetic computation time. The overall CPU time for the kinetic computation with 30 spatial and 30 velocity cells per edge has been approximately 210 min on an Intel Core i5-3230M with 2.6 GHz .
6. Conclusions. In this paper we developed nonlinear half-moment models for chemotaxis models which approximate very efficiently the solution of a network problem governed by kinetic equations. An APWB scheme is investigated for these nonlinear moment models. We derive from coupling conditions for the kinetic models corresponding conditions for the macroscopic models. Note that this procedure  also might be applied to other kinetic models. In the numerical tests we investigated the dynamics for different coupling conditions and for varying values of . A simulation on a larger network showed the applicability to more complicated structures and differences in behavior of macroscopic and kinetic models. Finally, we remark that there are several open issues. Most important would be analytical investigations of the kinetic network problem and a thorough investigation of the positivity preserving property of the APWB scheme for the density ρ in the nonlinear half moment model on the graph.