Synchronising and Nonsynchronising dynamics for a twospecies aggregation model
Abstract.
This paper deals with analysis and numerical simulations of a onedimensional twospecies hyperbolic aggregation model. This model is formed by a system of transport equations with nonlocal velocities, which describes the aggregate dynamics of a twospecies population in interaction appearing for instance in bacterial chemotaxis. Blowup of classical solutions occurs in finite time. This raises the question to define measurevalued solutions for this system. To this aim, we use the duality method developed for transport equations with discontinuous velocity to prove the existence and uniqueness of measurevalued solutions. The proof relies on a stability result. In addition, this approach allows to study the hyperbolic limit of a kinetic chemotaxis model. Moreover, we propose a finite volume numerical scheme whose convergence towards measurevalued solutions is proved. It allows for numerical simulations capturing the behaviour after blow up. Finally, numerical simulations illustrate the complex dynamics of aggregates until the formation of a single aggregate: after blowup of classical solutions, aggregates of different species are synchronising or nonsynchronising when collide, that is move together or separately, depending on the parameters of the model and masses of species involved.
Key words and phrases:
hydrodynamic limit, duality solution, twospecies chemotaxis, aggregate dynamics.2010 Mathematics Subject Classification:
35D30, 35Q92, 45K05, 65M08, 92D251. Introduction
Aggregation phenomena for a population of individuals interacting through an interaction potential are usually modelled by the socalled aggregation equation which is a nonlocal nonlinear conservation equation. This equation governs the dynamics of the density of individuals subject to an interaction potential . In this work, we are interested in the case where the population consists of two species which respond to the interaction potential in different ways. In the onedimensional case, the system of equations writes:
(1.1) 
with
where for are positive constants.
In this work, we are interested in the case where the interaction potential in (1.1) is pointy i.e. satisfies the following assumptions:

.

.

.

is concave with i.e.,
The aggregation equation arises in several applications in biology and physics. In fact, it is encountered in the modelling of cells which move in response to chemical cues. The velocity of cells depending on the distribution of nearby cells represents the gradient of the chemical substance which triggers the motion. Cells gather and form accumulations near regions more exposed to oxygen as observed in [20, 24]. We can also describe the movement of pedestrians using the aggregation equation as in [16] where the velocity of pedestrians is influenced by the distribution of neighbours. This equation can also be applied to model opinion formation (see [25]) where interactions between different opinions can be expressed by a convolution with the kernel .
From the mathematical point of view, it is known that solutions to the aggregation equation with a pointy potential blow up in finite time (see e.g [11, 5, 2]). Then globalintime existence for weak measure solutions has been investigated. In [5], existence of weak solutions for single species model has been obtained as a gradient flow. This technique has been extended to the twospecies model at hand in [11]. Another approach of defining weak solution for such kind of model has been proposed in [18, 17] for the single species case. In this approach, the aggregation equation is seen as a transport equation with a discontinuous velocity . Then solutions in the sense of duality have been defined for the aggregation equation.
Duality solutions has been introduced for linear transport equations with discontinuous velocity in the onedimensional space in [3]. Then it has been adapted to the study of nonlinear transport equations in [4, 18, 17]. In [18, 17], the authors use this notion of duality solutions for the onespecies aggregation equation. Such solutions are constructed by approximating the problem with particles, i.e. looking for a solution given by a finite sum of Dirac delta functions. Particles attract themselves through the interacting potential , when two particles collide, they stick to form a bigger particle.
In this work, we extend this approach to the two species case. To do so, we need to modify the strategy to the problem at hand. Indeed, collisions between particles of different species are more complex: particles can move together or separately after collision. This synchronising or nonsynchronising dynamics implies several difficulties for the treatment of the dynamics of particles. In fact, particles of different species can not stick when they collide. Then an approximate problem is constructed by considering the transport equation with the a regularized velocity. Then measure valued solutions are constructed by using a stability result.
An important advantage of this approach is that it allows to prove convergence of finite volume schemes. Numerical simulations of the aggregation equation for the onespecies case, which corresponds to the particular case of (1.1) when setting , have been investigated by several authors. In [8] the authors propose a finite volume method consistent with the gradient flow structure of the equation, but no convergence result has been obtained. In [9], a Lagrangian method is proposed (see also the recent work [7]). For the dynamics after blow up, a finite volume scheme which converges to the theoretical solution is proposed in [19, 6]. In the twospecies case, the behaviour is more complex since the interaction between the two species can occur and they may synchronise or not i.e. move together or separately depending on the parameters of the models and the masses of species. A numerical scheme illustrating this interesting synchronising or nonsynchronising dynamics is provided in Section 6. In addition, a theoretical result on the convergence of the numerical approximation obtained with our numerical scheme towards the duality solution is given. Such complex interactions phenomena have been observed experimentally in [13].
System (1.1) can be derived from a hyperbolic limit of a kinetic chemotaxis model. In the case of twovelocities and in one space dimension, the kinetic chemotaxis model is given by
(1.2) 
where stands for the distribution function of th species at time , position and velocity , is the concentration of the chemical substance, is the tumbling kernel from direction to direction and is a small parameter. This tumbling kernel being affected by the gradient of the chemoattractant, is chosen as in [12]
(1.3) 
where is a positive constant called the natural tumbling kernel and is the chemosensitivity to the chemical . This kinetic model for chemotaxis has been introduced in [21] to model the runandtumble process. Existence of solutions to this two species kinetic system has been studied in [14].
Summing and substracting equations (1.2) with respect to for yields
(1.4) 
(1.5) 
where and . Taking formally the limit in (1.5), we deduce that in the sense of distributions. Injecting in (1.4), we deduce formally that satisfies the limiting equation:
(1.6) 
where satisfies the elliptic equation:
This latter equation can be solved explicitly on and is given by
(1.7) 
Then we recover system (1.1). This formal computation can be made rigorous.
The rigorous derivation of (1.6) from system (1.2) will be proved in this work.
The paper is organized as follows. We first recall some basic notations and notions about the duality solutions and state our main results. Section 3 is devoted to the derivation of the macroscopic velocity used to define properly the product and duality solutions. Existence and uniqueness of duality solutions are proved in Section 4, as well as its equivalence to gradient flow solutions. The convergence of the kinetic model (1.2) as towards the aggregation model (1.6)(1.7) is shown in Section 5. Finally, a numerical scheme that captures the synchronising and nonsynchronising behaviour of the aggregate equation is studied in Section 6, as well as several numerical examples showing the synchronising and nonsynchronising dynamics.
2. Notations and main results
2.1. Notations
We will make use of the following notations. Let , we denote

is the space of nonnegative functions of .

is the space of continuous functions that vanish at infinity.

is the set of local Borel measures, those whose total variation is finite:

is the space of timecontinuous bounded Borel measures endowed with the weak topology.

is the Wasserstein space of order 2:

For , we define :
We notice that if satisfies 2 and 4, we have by taking in 4 and using 2 that, , . We deduce that 4 holds for i.e.:
(2.1) 
We recall a compactness result on . If there exists a sequence of bounded measures such that their total variations are uniformly bounded, then there exists a subsequence of that converges weakly to in .
2.2. Duality solutions
For the sake of completeness, we recall the notion of duality solutions which has been introduced in [3] for one dimensional linear scalar conservation law with discontinuous coefficients. Let us then consider the linear conservation equation:
(2.2) 
with . We assume weak regularity of the velocity field and satisfies the socalled onesided Lipschitz (OSL) condition:
(2.3) 
In order to define duality solutions, we introduce the related backward problem
(2.4) 
We define the set of exceptional solutions as follows
Definition 2.1 (Reversible solutions to (2.4)).
The following result shows existence and weak stability for duality solutions provided that the velocity field satisfied the onesidedLipschitz condition.
Theorem 2.3 (Theorem 2.1 in [3]).

There exists a bounded Borel function , called universal representative of such that a. e., and for any duality solution ,

Let be a bounded sequence in , with in . Assume that , where is bounded in . Consider a sequence of duality solutions to
such that is bounded in , and . Then in , where is the duality solution to
Moreover, weakly in .
2.3. Main results
Up to a rescaling, we can assume without loss of generality that the total mass of each species is normalized to 1. Then we will work in the space of probabilities for densities.
The first theorem states the existence and uniqueness of duality solutions for system (1.1) and its equivalence with the gradient flow solution considered in [11].
Definition 2.4.
Then, we have the following existence and uniqueness result:
Theorem 2.5 (Existence, uniqueness of duality solution and equivalence to gradient flow solution).
Let and . Under assumptions 1–4, there exists a unique duality solution to (1.1) in the sense of Definition 2.4 with such that
(2.6) 
This duality solution is equivalent to the gradient flow solution defined in [11].
In our second main result, we prove the convergence of the kinetic model (1.2) towards the aggregation model.
Theorem 2.6 (Hydrodynamical limit of the kinetic model).
The condition in the previous theorem is needed to guarantee that the tumbling kernel defined in (1.3) is positive.
3. Macroscopic velocity
In this section, we find the representative of for which existence and uniqueness of duality solutions hold. To this end, we consider the similar system of transport equations to (1.1) associated to the velocity which converges to . Next, the limit of the product for is computed.
3.1. Regularisation
We build a sequence which converges to by considering the sequence of regularised kernels approaching .
Lemma 3.1.
Let be the sequence of regular kernels defined by
Then
and
Proof.
From 1, and since is continuous at , we conclude that . From 2, we deduce that is an odd function. Using the definition of and 3, we get that . From the construction of , we have that outside the interval and from 4 one sees in in the sense of distributions. In addition, if we take and in 4, we have that
Since in , we conclude that in in the sense of distributions. Finally, we obtain that in the sense of distributions. ∎
In the rest of the paper, the notation will refer to the regularised kernels of Lemma 3.1. Given , the velocity is defined similarly to (2.6) as
(3.1) 
In the following lemma, we show that if and admit weak limits and respectively in , then the limit of the product is . Contrary to [22] where the twodimensional case is considered, this limiting measure does not charge the diagonal.
Lemma 3.2.
Proof.
Before starting the proof of the lemma, we first introduce some notations which simplify the computations
(3.2)  
For , we denote
Step 1: Convergence almost everywhere in time of .
Since , we have
where and are defined by
From the definition of in Lemma 3.1, it follows that
The estimate on in Lemma 3.1 and 3 imply that
with and defined in (3.2).
Let . Since the set converges to the empty set, there exists such that ,
(3.3) 
For all , we observe that , we have
(3.4) 
From the weak convergence of , , we note that the sequence converges weakly to . Since the total variation of is constant in , the tight convergence is achieved. Then, there exists such that
From (3.4) and (3.3), we conclude that ,
Hence, for all , we get
(3.5) 
We deduce that .
Next, we show that tends to zero.
where is an integer which will be fixed later. From the construction of in Lemma 3.1, we get
Therefore, one has
Let . Using (3.4), by the same token as previously, there exists such that for all ,
Setting , we conclude that for all ,
For , we notice that is a continuous function that vanishes on the diagonal and we have
The tight convergence of to implies that there exists such that for all
Therefore for all , one has
(3.6) 
This implies that converges to 0.
Combining (3.5) and (3.6), we deduce that for almost every , converges to 0.
Step 2: Lebesgue’s dominated convergence theorem
For all , we have that
Since converges almost everywhere to 0, converges to zero from Lebesgue’s dominated convergence theorem. ∎
3.2. OSL condition on the macrosocopic velocity
Proposition 3.3.
4. Existence and uniqueness of duality solutions
4.1. Proof of the existence of duality solutions in Theorem 2.5
The proof is divided into several steps. First, we construct an approximate problem for which the existence of duality solutions holds.
Then, we pass to the limit in the approximate problem to get the existence of duality solutions thanks to the weak stability of duality solutions stated in Theorem 2.3
and recover Equation (2.5) from Lemma 3.2.
Finally, we recover the bound on the second order moment.
Step 1: Existence of duality solutions for the approximate problem
The macroscopic velocity is replaced by an approximation defined in (3.1) and the following system is considered:
(4.1) 
Since is not Lipschitz continuous, we first consider an approximation of obtained by a convolution with a molifier. The solution to the following equation is investigated.
(4.2) 
where is given by
Applying Theorem 1.1 in [10] gives the existence of solutions in and . Since the velocity field is Lipschitz, is a duality solution. In addition, for we have for the following estimate:
Then,
(4.3) 
Using (4.3) and the density of in , we deduce that .
Moreover, the equicontinuity of in follows from (4.3) and the density of in .
Since , Ascoli Theorem gives the existence of a subsequence in of which converges to a limit named in .
We pass to the limit when tends to infinity in Equation (4.2) and obtain that satisfies (4.1).
Step 2 : Extraction of a convergent subsequence of and existence of duality solutions.
As above, there exists a subsequence of in such that
Let us find the equation satisfied by in the distributional sense. Let be in . Since satisfies (4.1) in the distributional sense, we have
Using Lemma 3.2, we can pass to the limit in the latter equation and obtain,
Thus satisfies (2.5) in the sense of distributions. From Proposition 3.3, the macroscopic velocity satisfies an uniform OSL condition. Then, by weak stability of duality solutions in (see Theorem 2.3 (3)), we deduce that
Step 3 : Finite second order moment.
From Equation (2.5), we deduce that the first and second moments satisfy in the sense of distributions