NUMERICAL STUDY OF TWO-SPECIES CHEMOTAXIS MODELS

We first conduct a comparative numerical study of two recently proposed two-species chemotaxis models. We show that different scenarios are possible: depending on the initial masses, either one or both cell densities may blow up, or a global solution may exist. In particular, our numerical results indicate answers on some open questions of possible blow up stated in [4, 7]. We then introduce two regularizations of the studied models and demonstrate that their solutions are capable of developing spiky structure without blowing up.


1.
Introduction. Chemotaxis is one of the most important mechanism in movement of biological species. It describes a collective movement of cells or biological species that is oriented towards the chemoattractant gradient. A classical PDEbased model of chemotaxis, the Patlak-Keller-Segel (PKS) system, was first proposed in [27] and [18,19]. The classical PKS system as well as its more recent modifications (see, e.g., [2,14,16,17,28] and references therein) are capable of describing the aggregation of biological species and consequently a mechanism for self-organization of biological systems. In some of the above models, the cell aggregation may lead to a finite time blow-up of the solution provided the total initial cell mass is above a certain threshold (see [3,11,24]). However, the formation of singularities in the PKS model may be viewed as a purely mathematical artifact. Many regularized PKS-type systems admit global classical (yet spiky) solutions, which may better describe the biological aggregation phenomenon (as, e.g., in the regularized models studied in [2,23,26,29,33], see also a recent review [14]).
The aim of this paper is twofold. We first numerically investigate two models proposed and analytically studied in [4,5,6,7,8,9,20,34]. The first system reads    (ρ 1 ) t + χ 1 ∇ · (ρ 1 ∇c) = µ 1 ∆ρ 1 , (ρ 2 ) t + χ 2 ∇ · (ρ 2 ∇c) = µ 2 ∆ρ 2 , c t = D∆c + α 1 ρ 1 + α 2 ρ 2 − βc, where ρ 1 and ρ 2 denote the cell densities of the first and second species, respectively, and c is the concentration of the chemoattractant. The positive constants µ i , χ i , α i (i = 1, 2) and β are parameters of cell diffusion, chemotactic sensitivities, production and consumption rates, respectively. Finally, D is the chemoattractant diffusion coefficient. Throughout the paper we will assume that the second species has larger chemotactic sensitivity than the first one, that is, Since the molecular diffusion is typically much faster than the cell diffusion, that is, both µ 1 D and µ 2 D, the system (1) can be simplified by assuming that µ i /D ≈ 0, i = 1, 2, which lead to the following system:    (ρ 1 ) t + χ 1 ∇ · (ρ 1 ∇c) = µ 1 ∆ρ 1 , (ρ 2 ) t + χ 2 ∇ · (ρ 2 ∇c) = µ 2 ∆ρ 2 , ∆c + γ 1 ρ 1 + γ 2 ρ 2 − ζc = 0, The models (1) and (2) can be viewed as direct extensions of the PKS system to the case of the chemotaxis motion of two noncompeting species that produce the same chemoattractant. As in the case of the classical PKS model, the main question we numerically investigate is: what are the conditions on initial masses and chemotactic parameters that determine whether the solution remain smooth and bounded or it blows up in a finite time? In the latter case, we would also like to determine whether the cell densities of both species blow up simultaneously or the cell density of the species with a larger chemosensitivity constant blows up first. These questions were first raised in [34] and further studied in [4,5,7]. More precisely, it has been proven analytically that depending on a particular set of parameter values and initial cell densities the solution may be globally regular or it may blow up within a finite time. However, in some cases, the theory fails to predict the behavior of the solution. Our numerical simulations indicate answers to some of these open questions.
The second main goal of our paper is to present and study two regularizations of (1). The first regularized system reads where the functions Q 1 and Q 2 are smooth saturated chemotactic fluxes Q i (u 1 , . . . , i (u) that satisfy the following properties: where C j i are constants. This regularization is similar to the one proposed in [2] for the PKS system. It is based on a fundamental biological property of the chemotactic flux-its boundedness (this feature is almost always lost in weakly nonlinear, small gradients expansions, underlying the derivation of most continuum models). The synthesized form of the saturated fluxes a function, which has the following universal features: It is linear at small gradients of c and is bounded at large gradients of c. There is a certain arbitrariness in the choice of the chemotactic flux functions Q i . A typical example of saturated chemotactic fluxes, which is used in all of our numerical experiments, is for i = 1, 2 where s * i are switching parameters, which define small gradient values, for which the system (3) reduces to the original system (1) (or (2)) so that the effect of saturated chemotactic flux is felt at large gradient regimes only. This is expected to result in solutions which are spiky but yet bounded for all times. The second regularized system, is similar to the density-dependent regularization from [31,32]. Here, κ is a (small) regularization parameter and κ → 0 leads to the original system (1). As in the single species case [13,31], one expects a global classical solution of (6) to exist. At the same time, we will numerically demonstrate that the solutions of (6) typically have spiky structure and thus the system (6) can be used to model the aggregation phenomenon.
The present paper will be organized as follows. In Section 2, we prove a-priori estimates for the regularized system (3). Section 3 is devoted to the description of the numerical method, which is a modification of the second-order positivity preserving upwind scheme from [1,2]. Finally, Section 4 is devoted to the presentation and discussion of our numerical simulations.
2. L ∞ bounds via Moser-Alikakos iteration. Let us consider the original initial-boundary value problem (IBVP) for the system (3): where ∂Ω is a Lipschitz continuous boundary with the outer normal n.
In this section, we will prove a-priori estimates for positive solutions of (7). The following result is a generalization of analogous result for the one-species chemotaxis system from [2]. Theorem 2.1. Let (ρ 1 (x, t), ρ 2 (x, t), c(x, t)) be a positive classical solution of the IBVP (7) with bounded nonnegative initial data. Then, for all x ∈Ω and t ≥ 0, c(x, t) ≤ c 0 where C = C(d, Ω) is a constant, which depends on d and Ω only, and C i = max j C (j) i , i = 1, 2. Proof. We begin by multiplying the density equations in (7) by ρ s−1 1 and ρ s−1 2 (s ≥ 2), respectively. Then, integrating by parts, applying the chain rule, using the boundedness of |Q i | and the inequality (see [2,21]) with a suitable ε yields the following estimates: The last term in (12) is estimated using the inequality (11) with u = ρ s 2 i and ε such that Substituting (13) into (12), we obtain

TWO-SPECIES CHEMOTAXIS 135
We then fix T ∈ (0, ∞), multiply both sides of the last inequality by the integrating factor e κt , where κ : Let us now define the function which satisfies (from (14)): where the constant K depends on d and Ω only. Taking s = 2 k , k = 1, 2, . . . and applying the above inequality recursively we obtain Finally, we note the total mass of the cells remains constant in time (this can be verified by integrating the density equations in (7) over Ω), and therefore and the estimate (8) for the cell densities ρ i (x, t) follows from (16), (17). Let us denote by M (1) := max i=1,2 (M 1 (1), M 2 (1)) . In order to obtain a bound on chemoattractant concentration c(x, t), we compare it with the solution of the following initial value problem (IVP): The comparison principle then yields which completes the proof of Theorem 2.1.

Remark 1.
An alternative L ∞ bounds on ρ 1 and ρ 2 can be obtained using the results from [15]. However, the estimates (8) and (9) explicitly show the dependence of the bounds on the number of dimensions d and the viscosity coefficients µ 1 and µ 2 .
Remark 2. Note that using the L ∞ -bounds established in Theorem 2.1, parabolic boundary L p -estimates and Schauder estimates (see, e.g., [21]), one can obtain that (ρ 1 ) t , (ρ 2 ) t , c t and all spatial partial derivatives of ρ 1 , ρ 2 and c up to order two are bounded onΩ × [0, ∞). This will lead to a global existence result similar to the one established in [2] for the one-species chemotaxis system. As it has been illustrated by our numerical experiments (shown in Section 4), the regularized solution while being bounded for all times, may develop spiky (even multi-spiky) structures that model aggregation phenomena.
3. Numerical scheme. The numerical results presented in this paper are obtained using a second-order positivity preserving upwind scheme, which is a rather straightforward extension of the hybrid finite-volume-finite-difference scheme developed in [1,2] for the single species chemotaxis models. In this section, we briefly describe the scheme for the two-species chemotaxis system where ε = 0 or ε = 1 and the functions g, Q 1 and Q 2 may be either linear or nonlinear so that the system (19) reduces to either (1), (2), (3) or (6). We introduce a Cartesian mesh consisting of the uniform cells of the size ∆x∆y centered at (x j , y k ). The computed quantities are the cell averages of cell densities ρ i , and the point values of the chemoattractant concentration c, c j,k (t) ≈ c(x j , y k , t), which are evolved in time according to the semi-discrete scheme: Here, (H i ) x j+ 1 2 ,k and (H i ) y j,k+ 1 2 are the following upwind numerical fluxes: where Q (1) i and Q (2) i are the components of the vector function Q i (see (4)) and The point values (ρ i ) are obtained using the piecewise linear reconstruction with the slopes ((ρ i ) x ) j,k and ((ρ i ) y ) j,k calculated using the minmod2 limiter (see, e.g., [22,25,30]): Thus, we have for i = 1, 2 . Remark 3. Notice that in the above formulae, the quantities (ρ i ) j,k , c j,k , (H i ) x , ((ρ i ) x ) j,k , ((ρ i ) y ) j,k and the functions ( ρ i )(x, y), i = 1, 2 depend on time, but we suppress this dependence for brevity.
If ε = 1, then the semi-discrete scheme (20) is a system of time-dependent ODEs, which has to be integrated numerically using a stable and accurate ODE solver. In this paper, we have used the third-order strong stability preserving (SSP) Runge-Kutta method from [10]. The efficiency of the fully discrete method can be improved by applying an SSP implicit-explicit Runge-Kutta method (see, e.g., [12] and references therein), as discussed in [1].
If ε = 0, then the last equation in (19) becomes an elliptic equation for c, and consequently the last equation in (20) becomes a system of linear algebraic equations. This system has to be solved using a proper linear algebra solver. One time step of the resulting algorithm will then consist of an explicit time advance of ρ 1 and ρ 2 followed by solving the last equation in (20) for c using the values of ρ 1 and ρ 2 from the new time level. 4. Numerical experiments. In this section, we present the results of our numerical experiments that clarify the behavior of the solutions of the studied two-species chemotaxis systems in two space dimensions. We restrict our consideration to the two-dimensional (2-D) case since the theoretical results/open questions in [4,5,7] were obtained/formulated for the 2-D version of the system (2).

4.1.
Parabolic-elliptic systems. We first consider the parabolic-elliptic system (2) and its regularization with a bounded chemotaxis fluxes (with Q 1 and Q 2 satisfying (4)): Without loss of generality we set µ 2 = 1. We denote by θ 1 and θ 2 the initial masses Following [4,7], we split the (θ 1 , θ 2 )-plane into the following four regions, outlined in Figure 1: • Region B: • Region C: • Region D: In [4,7], the following results were proved for the 2-D IVP for the parabolicelliptic system (2) with γ 1 = γ 2 = ζ = µ 2 = 1: • There is a global classical solution in Region A (the proof is based on the energy functions that provide a-priori estimates for the entropy of (2)); • In Region C, ρ 2 blows up faster than ρ 1 ; • In Region D, ρ 1 and ρ 2 blow up at the same rate.
The question on the solution behavior in Region B remains open and we investigate it numerically. We study the systems (2) and (22) (2). Example 1. Global existence in Region A.. We first consider the system (2) with χ 1 = 1, χ 2 = 10, µ 1 = 1, and subject to the following initial data: As one can see in Figure 2, the magnitude of both ρ 1 and ρ 2 decays and the solution remains smooth and bounded. . We now consider the system (2) with χ 1 = 1, χ 2 = 20, µ 1 = 1, and subject to the same, radially symmetric initial data (23). Figure 3 suggests that ρ 2 blows up while ρ 1 remains bounded. Moreover, the magnitude of ρ 1 seems to decay in time. However, this would contradict the analytical results on simultaneous blow-up in Region B proved in [5,7] for radially symmetric initial data. We therefore perform a mesh refinement study to carefully monitor the behavior of max Ω ρ 1 (x, y, 0.05) and max Ω ρ 2 (x, y, 0.05) as a function of N , where the computational grid is N × N . The dependence of max Ω ρ 1 (x, y, 0.05) on N together with the algebraic function ξ 1 (N ) = 1.266(2N + 10) 1/4 are shown in Figure 4 (left). This results indicate that ρ 1 still blows up, but does not develop a δ-type singularity and therefore its blow-up is extremely hard to verify numerically.
In contrast, ρ 2 collapses to a δ-function as indicated in Figure 4 (right), where we plot max Ω ρ 2 (x, y, 0.05) as a function of N together with the quadratic function ξ 2 (N ) = 0.0266(N − 8) 2 . Note that this quadratic increase indeed reflects a δtype singularity since using either a finite-volume, finite-difference or finite-element method 2-D δ-functions can only be resolved so that max j,k (ρ 2 ) j,k ∼ 1 ∆x∆y .  Example 3. Blow-up of ρ 1 and ρ 2 with non-radial initial data in Region B.. The theoretical blow-up results for Region B reported in [5,7] only apply to radially symmetric initial data. However, behavior of solutions with non-radially symmetric initial data is still an open problem. In order to numerically investigate this case we take the following initial data: ρ 1 (x, y, 0) = 12.5 e −100(x 2 /16+y 2 ) , ρ 2 (x, y, 0) = 12.5 e −100(x 2 +y 2 /16) , and numerically solve the IVP (2)  To better understand the difference in the behavior of ρ 1 and ρ 2 , we perform a mesh refinement study similar to the one conducted in Example 2. The results shown in Figure 6 support the conjecture that non-radially symmetric initial data from Region B lead to the same different types of blow-up as in the radially symmetric case. Example 4. ρ 2 blows up faster than ρ 1 in Region C.. Next, we consider the system (2) with χ 1 = 6, χ 2 = 100, µ 1 = 1, and subject to the following initial data: ρ 1 (x, y, 0) = 10 e −100(x 2 +y 2 ) , ρ 2 (x, y, 0) = 90 e −100(x 2 +y 2 ) . Figure 7 shows that both ρ 2 and ρ 1 blow up, while c stays bounded. One can also observe that ρ 2 seems to blow up faster than ρ 1 . To verify this, we perform the mesh refinement study and plot the results obtained on the 200 × 200 and 400 × 400 uniform grids at the same time t = 0.007. As one can see, the magnitude of ρ 2 increases by a factor of about 4 (from 8.9258 · 10 3 to 3.2590 · 10 4 ), which clearly indicates that by this time ρ 2 has already blown up. At the same time, ρ 1 increases only by a factor of about 2 (from 55.9119 to 105.1929), which means that ρ 1 is going to blow up a little later. Notice that this numerical experiment confirms the analytical result from [4,5,7].

Remark 4.
We have conducted more numerical experiments (not reported here for the sake of brevity) that confirm the analytical results from [4,7]: global existence of the solution in Region A as well as simultaneous blow-up in Region D. 4.1.2. The regularized system (22). Example 5. Region B, regularized solution. We now consider the system (22), (5) with s * 1 = s * 2 = 20 and the same data as in Example 2. A mesh refinement study presented in Figure 8 clearly demonstrates that saturated chemotaxis flux prevents blow-up of ρ 2 though a spiky structure is developed. Example 6. Region C, regularized solution. Here, we consider the regularized system (22), (5) with s * 1 = s * 2 = 20 and the same data as in Example 4. The obtained solution is shown in Figure 9, where one can see a spiky structure in both density components. Note that while the magnitude of ρ 2 has increased, the magnitude of ρ 1 has slightly decreased. A mesh refinement study (not presented here for the sake of brevity) indicates that unlike the solution of the original system (2) (shown in Example 4), the solution of the regularized system (22) does not blow up. Moreover, by the time t = 0.05 the regularized solution has already reached its (numerical) steady state.

Parabolic systems.
In this section, we numerically study behavior of solutions of the parabolic system (1) and its two regularizations (3) and (6). In the parabolic case, no analytical results that could have split the (θ 1 , θ 2 )-plane into particular regions (as it has been done in the parabolic-elliptic case in Section 4.1) are available. However, one can expect the solutions of the parabolic system to behave rather similarly to the solutions of the parabolic-elliptic systems, especially when µ 1 /D and µ 2 /D are small (in all of our numerical experiments, we take µ 1 = µ 2 = 1 and D = 10).
A snapshot of the computed solution at time t = 0.01 is plotted in Figure 10. We then double the chemotactic sensitivity of the second species and take χ 2 = 10. This leads to the nonmonotone behavior of ρ 2 : Its magnitude first increases, but then starts decreasing and by time t = 0.05 (see Figure 11) the solution looks similar to the one obtained with χ 2 = 5 (the only difference is that the ratio between the maximum values of ρ 2 and ρ 1 is now about 3 times larger.  Example 8. ρ 2 blows up faster than ρ 1 . Next, we consider the system (1) subject to the same initial data (25), but with much larger chemotactic sensitivity constants χ 1 = 5 and χ 2 = 60. Figure 12 shows that both ρ 2 and ρ 1 blow up, but ρ 2 blows up  faster than ρ 1 as it is confirmed by the performed mesh refinement study. As one can see, when the mesh size is doubled, the magnitude of ρ 2 increases by a factor of about 4 (from 4.5158 · 10 4 to 1.8027 · 10 5 ), which clearly indicates that by this time ρ 2 has already blown up. At the same time, ρ 1 increases only by a factor of less than 2 (from 1.5263 · 10 3 to 2.8660 · 10 3 ), which means that ρ 1 is going to blow up a little later. Example 9. Blow-up of large initial data. In this example, we take much larger initial data, ρ 1 (x, y, 0) ≡ ρ 2 (x, y, 0) = 5000 e −100(x 2 +y 2 ) , c(x, y, 0) ≡ 1, (26) and the same chemotactic sensitivity constants χ 1 = 5 and χ 2 = 60. As it can be seen in Figure 13, both ρ 1 and ρ 2 blow up now much faster than in Example 8.  (3) and (6). Example 10. Spiky solutions-no blow-up. We now study the behavior of the solutions of the regularized systems (3) and (6) subject to the same initial data (26) as in Example 9. We first compute the solution of the system (3), (5) with s * 1 = s * 2 = 20. As one can see in Figure 14, both ρ 1 and ρ 2 increase and the spikes are formed (notice that they have about the same magnitude even though χ 2 is much larger than χ 1 : this is an effect of the regularization). However, the mesh refinement study clearly demonstrates that the solution has not blown up. Our further numerical studies indicate that the obtained spiky solution is a steady state: it does not change as the time increases.
The solution of the second regularized system (6) with κ = 0.01 is shown in Figure 15. As one can clearly see, this regularized solution does not blow up as well. However, it behaves differently: by the time t = 0.01, max Ω ρ 2 has increased, while max Ω ρ 1 has decreased. As in the previous case, the obtained solution is a numerical steady state.

5.
Conclusions. In this paper, we have presented a comprehensive numerical study of several two-species chemotaxis Patlak-Keller-Segel type models. We have considered both the parabolic-elliptic as well as the fully parabolic systems. The simplest (yet very challenging for rigorous mathematical analysis) parabolic-elliptic case has been analytically studied in [4,7]. It has been proven there that under certain conditions on the initial cell densities and chemotactic sensitivity coefficients the system admits global regular solutions, while under a different set of conditions the densities of both species will simultaneously blow up within a finite time. If none of those conditions is satisfied, the question of global existence vs. finite time blow-up remains open. The aim of the present paper has been to present an extensive numerical study of possible configurations and indicate answers to some open questions posted in [4,7]. More precisely, we have demonstrated that for the parabolic-elliptic system the following 3 scenarios are possible: a global solution may exist, the density of one species may blows up faster than the density of the second species, both densities may blow up simultaneously.
In the fully parabolic case, the situation is more complicated and a complete identification of the corresponding conditions for the initial data and parameters is not yet available. Nevertheless, we have demonstrated that the same 3 scenarios are still possible.
Since blow up of the solution can be viewed as a purely mathematical artifact of the classical Patlak-Keller-Segel type models, we have also studied two different regularizations that yield spiky but bounded solutions. We have derived a-priori estimates that confirm that the solutions of the regularized system indeed remain bounded. We have also conducted a number of numerical experiments to study behavior of the obtained spiky solutions.