Synchronising and non-synchronising dynamics

# Synchronising and Non-synchronising dynamics for a two-species aggregation model

Casimir Emako-Kazianou Sorbonne Universités, UPMC Univ Paris 06, UMR 7598, Laboratoire Jacques-Louis Lions, F-75005, Paris, France
CNRS, UMR 7598, Laboratoire Jacques-Louis Lions, F-75005, Paris, France
INRIA-Paris-Rocquencourt, EPC MAMBA, Domaine de Voluceau, BP105, 78153 Le Chesnay Cedex
Jie Liao Department of Mathematics, East China University of Science and Technology, Shanghai, 200237, P. R. China  and  Nicolas Vauchelet Sorbonne Universités, UPMC Univ Paris 06, UMR 7598, Laboratoire Jacques-Louis Lions, F-75005, Paris, France
CNRS, UMR 7598, Laboratoire Jacques-Louis Lions, F-75005, Paris, France
INRIA-Paris-Rocquencourt, EPC MAMBA, Domaine de Voluceau, BP105, 78153 Le Chesnay Cedex
July 3, 2019
###### Abstract.

This paper deals with analysis and numerical simulations of a one-dimensional two-species hyperbolic aggregation model. This model is formed by a system of transport equations with nonlocal velocities, which describes the aggregate dynamics of a two-species population in interaction appearing for instance in bacterial chemotaxis. Blow-up of classical solutions occurs in finite time. This raises the question to define measure-valued 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 measure-valued 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 measure-valued 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 blow-up 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, two-species chemotaxis, aggregate dynamics.
###### 2010 Mathematics Subject Classification:
35D30, 35Q92, 45K05, 65M08, 92D25

## 1. Introduction

Aggregation phenomena for a population of individuals interacting through an interaction potential are usually modelled by the so-called 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 one-dimensional case, the system of equations writes:

 (1.1) ∂tρα+χα∂x(a(ρ)ρα)=0,for α=1,2,

with

 a(ρ):=∫R∂xK(x−y)ρ(t,dy),ρ:=θ1ρ1+θ2ρ2,

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:

1. .

2. .

3. .

4. is -concave with i.e.,

 ∀x,y∈R∗,(∂xK(x)−∂xK(y))(x−y)≤λ(x−y)2.

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  where the velocity of pedestrians is influenced by the distribution of neighbours. This equation can also be applied to model opinion formation (see ) 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 global-in-time existence for weak measure solutions has been investigated. In , existence of weak solutions for single species model has been obtained as a gradient flow. This technique has been extended to the two-species model at hand in . 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 one-dimensional space in . 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 one-species 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 non-synchronising 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 one-species case, which corresponds to the particular case of (1.1) when setting , have been investigated by several authors. In  the authors propose a finite volume method consistent with the gradient flow structure of the equation, but no convergence result has been obtained. In , a Lagrangian method is proposed (see also the recent work ). For the dynamics after blow up, a finite volume scheme which converges to the theoretical solution is proposed in [19, 6]. In the two-species 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 non-synchronising 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 .

System (1.1) can be derived from a hyperbolic limit of a kinetic chemotaxis model. In the case of two-velocities and in one space dimension, the kinetic chemotaxis model is given by

 (1.2) ⎧⎪⎨⎪⎩∂tfεα+v∂xfεα=1ε∫V(Tα[S](v′,v)fεα(v′)−Tα[S](v,v′)fεα(v))dv′,α=1,2,v∈V={±1},−∂xxSε+Sε=θ1(fε1(1)+fε1(−1))+θ2(fε2(1)+fε2(−1)),

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 

 (1.3) Tα[S](v,v′)=ψα(1+χαv∂xS),

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  to model the run-and-tumble process. Existence of solutions to this two species kinetic system has been studied in .

Summing and substracting equations (1.2) with respect to for yields

 (1.4) ∂tρεα+∂xJεα=0,
 (1.5) ∂tJεα+∂xρεα=2ψαε(χα∂xSερεα−Jεα),α=1,2,

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) ∂tρ0α+χα∂x((∂xS0)ρ0α)=0,

where satisfies the elliptic equation:

 −∂xxS0+S0=θ1ρ01+θ2ρ02.

This latter equation can be solved explicitly on and is given by

 (1.7) S0=K∗(θ1ρ01+θ2ρ02),K=12e−|x|.

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 non-synchronising behaviour of the aggregate equation is studied in Section 6, as well as several numerical examples showing the synchronising and non-synchronising 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:

 Mb(R)={μ∈Mloc(R),|μ|(R)<+∞}.
• is the space of time-continuous bounded Borel measures endowed with the weak topology.

• is the Wasserstein space of order 2:

 P2(R)={μnonnegative borel measures in% Rs.t |μ|(R)=1,∫R|x|2μ(dx)<∞}.
• For , we define :

 ˆH={H(x),for x≠0,0,else.

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) ∀x,y∈R,(ˆ∂xK(x)−ˆ∂xK(y))(x−y)≤λ(x−y)2.

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  for one dimensional linear scalar conservation law with discontinuous coefficients. Let us then consider the linear conservation equation:

 (2.2) ∂tρ+∂x(b(t,x)ρ)=0,in ]0,T[×R,

with . We assume weak regularity of the velocity field and satisfies the so-called one-sided Lipschitz (OSL) condition:

 (2.3) ∂xb≤γ(t),γ∈L1(]0,T[),in the sense of % distributions.

In order to define duality solutions, we introduce the related backward problem

 (2.4) {∂tp+b∂xp=0,p(T,⋅)=pT∈Liploc(R).

We define the set of exceptional solutions as follows

 E:={p∈Liploc(]0,T[×R)solution % to (???) with pT=0}.
###### Definition 2.1 (Reversible solutions to (2.4)).

We say that is a reversible solution to (2.4) if and only if satisfies (2.4) and is locally constant on , where is defined by

 Ve:={(t,x)∈]0,T[×R;∃pe∈Epe(t,x)≠0}.
###### Definition 2.2 (Duality solutions to (2.2), see ).

We say is a duality solution to (2.2) in if for any , and any reversible solution to (2.4) compactly supported in , the function is constant on .

The following result shows existence and weak stability for duality solutions provided that the velocity field satisfied the one-sided-Lipschitz condition.

###### Theorem 2.3 (Theorem 2.1 in ).
1. Given . Under the assumption (2.3), there exists a unique , duality solution to (2.2), such that .

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

 ∂tρ+∂x(^bρ)=0,in the distributional % sense.
3. Let be a bounded sequence in , with in . Assume that , where is bounded in . Consider a sequence of duality solutions to

 ∂tρn+∂x(bnρn)=0,in ]0,T[×R,

such that is bounded in , and . Then in , where is the duality solution to

 ∂tρ+∂x(bρ)=0,in ]0,T[×R,ρ(0,⋅)=ρini.

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 .

###### Definition 2.4.

(Duality solutions for system (1.1)) We say that is a duality solution to (1.1) if there exists and satisfying in the sense of distributions, such that for all ,

 (2.5) ∂tρα+χα∂x(^a(ρ)ρα)=0,for α=1,2,ρ=θ1ρ1+θ2ρ2,

in the sense of duality on and a.e. We emphasize that the final datum for (2.5) should be instead of .

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 14, there exists a unique duality solution to (1.1) in the sense of Definition 2.4 with such that

 (2.6) ^a(ρ):=∫Rˆ∂xK(x−y)ρ(t,dy),ρ=θ1ρ1+θ2ρ2.

This duality solution is equivalent to the gradient flow solution defined in .

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).

Assume that for . Let and be a solution to the kinetic-elliptic equation (1.2) such that and and .
Then, as , converges in the following sense:

 ρεα:=fεα(1)+fεα(−1)⇀ραweakly in SM,for α=1,2, Sε⇀SinC([0,T],W1,∞(R))−weak,

where is the unique duality solution of (1.6) and given in Theorem 2.5.

The condition in the previous theorem is needed to guarantee that the tumbling kernel defined in (1.3) is positive.

To conclude this Section on the main results, we emphasize that, a finite volume scheme to simulate (2.5) is proposed in Section 6 and its convergence towards duality solutions is stated in Theorem 6.3.

## 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

 ∂xKn∈C0(R),∀x∈R,∂xKn(−x)=−∂xKn(x),

and

 ∥∂xKn∥∞≤∥∂xK∥∞,∂xxKn≤λin the distributional sense.
###### 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

 n∂xK(1n)≤λ.

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) ∀ρ∈SM,an(ρ):=∫R∂xKn(x−y)ρ(t,dy).

In the following lemma, we show that if and admit weak limits and respectively in , then the limit of the product is . Contrary to  where the two-dimensional case is considered, this limiting measure does not charge the diagonal.

###### Lemma 3.2.

For , let be a sequence such that . Suppose that there exists in such that

 ρnα⇀ραweakly in SM,

Then, we have

 an(ρnα)ρnβ⇀^a(ρα)ρβweakly in Mb([0,T]×R),α,β=1,2,

where and are defined in (2.6),(3.1) respectively. That is for ,

 ∫T0∫R2∂xKn(x−y)ρnα(t,dx)ρnβ(t,dy)ϕ(t,x)dt→∫T0∫R2ˆ∂xK(x−y)ρα(t,dx)ρβ(t,dy)ϕ(t,x)dt.
###### Proof.

Before starting the proof of the lemma, we first introduce some notations which simplify the computations

 (3.2) μn(t,dx,dy) :=ρnα(t,dx)⊗ρnβ(t,dy),En:={(x,y)∈R2,x≠y,|x−y|≤1n}, μ(t,dx,dy) :=ρα(t,dx)⊗ρβ(t,dy).

For , we denote

 An(t):=∫R2∂xKn(x−y)μn(t,dx,dy)ϕ(t,x)−∫R2ˆ∂xK(x−y)μ(t,dx,dy)ϕ(t,x),

Step 1: Convergence almost everywhere in time of .
Since , we have

 An(t) =∫R2ˆ∂xKn(x−y)μn(t,dx,dy)ϕ(t,x)−∫R2ˆ∂xK(x−y)μ(t,dx,dy)ϕ(t,x), =In(t)+IIn(t),

where and are defined by

 In(t) :=∫R2(ˆ∂xKn(x−y)−ˆ∂xK(x−y))μn(t,dx,dy)ϕ(t,x), IIn(t) :=∫R2ˆ∂xK(x−y)(μn(t,dx,dy)−μ(t,dx,dy))ϕ(t,x).

From the definition of in Lemma 3.1, it follows that

 In(t)=∫En(∂xKn(x−y)−∂xK(x−y))μn(t,dx,dy)ϕ(t,x).

The estimate on in Lemma 3.1 and 3 imply that

 ∣∣In(t)∣∣≤2∥ϕ∥L∞∥∂xK∥L∞μn(t,En),

with and defined in (3.2).
Let . Since the set converges to the empty set, there exists such that ,

 (3.3) μ(t,En)≤ε.

For all , we observe that , we have

 (3.4) μn(t,En)≤μn(t,EN)≤(μn−μ)(t,EN)+μ(t,EN).

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

 |μn−μ|(t,EN)≤ε.

From (3.4) and (3.3), we conclude that ,

 μn(t,En)≤2ε.

Hence, for all , we get

 (3.5) ∣∣In(t)∣∣≤2∥ϕ∥L∞∥∂xK∥L∞μn(t,En)≤4∥ϕ∥L∞∥∂xK∥L∞ ε.

We deduce that .

Next, we show that tends to zero.

 IIn(t) =∫R2(ˆ∂xK(x−y)−ˆ∂xKR(x−y))ϕ(t,x)(μn(t,dx,dy)−μn(t,dx,dy)) +∫R2ˆ∂xKR(x−y)ϕ(t,x)(μn(t,dx,dy)−μ(t,dx,dy)), :=II1n(t)+II2n(t),

where is an integer which will be fixed later. From the construction of in Lemma 3.1, we get

 II1n=∫ER(∂xK(x−y)−∂xKR(x−y))ϕ(t,x)(μn(t,dx,dy)−μ(t,dx,dy)).

Therefore, one has

 ∣∣II1n(t)∣∣≤2∥∂xK∥L∞∥ϕ∥L∞(μn(t,ER)+μ(t,ER)).

Let . Using (3.4), by the same token as previously, there exists such that for all ,

 μn(t,EN)≤2ε,μ(t,En)≤ε,

Setting , we conclude that for all ,

 ∣∣II1n(t)∣∣≤6ε∥∂xK∥L∞∥ϕ∥L∞.

For , we notice that is a continuous function that vanishes on the diagonal and we have

 ∫R2ˆ∂xKN(x−y)ϕ(t,x)(μn−μ)(t,dx,dy)=∫R2∂xKN(x−y)ϕ(t,x)(μn−μ)(t,dx,dy).

The tight convergence of to implies that there exists such that for all

 ∣∣∣∫R2ˆ∂xKM(x−y)ϕ(t,x)(μn−μ)(t,dx,dy)∣∣∣≤ε.

Therefore for all , one has

 (3.6) ∣∣IIn(t)∣∣≤ε(1+6∥∂xK∥L∞∥ϕ∥L∞).

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

 |An(t)|≤2∥ϕ∥L∞∥∂xK∥L∞MαMβ.

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.

Let be positive constants and be a positive measure such that . Let be such that assumption 4 hold. Let and be defined in (2.6) and (3.1) respectively. Then, there exists such that

 ∂x^a(t,x)≤κ(t),∂xan(t,x)≤κ(t), in the sense of distributions.
###### Proof.

For , we compute:

 (^a(ρ)(t,x)−^a(ρ)(t,y))(x−y)=∫R(ˆ∂xK(x−z)−ˆ∂xK(y−z))(x−y)ρ(t,dz).

Using the -concavity of , we deduce from (2.1)

 (^a(ρ)(t,x)−^a(ρ)(t,y))(x−y)≤λ(x−y)2∫Rρ(t,dz)≤λM(x−y)2.

Since is also concave from the proof of Lemma 3.1, we get the one-sided Lipschitz estimate on by the same token as for . ∎

## 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) ∂tρnα+χα∂x(an(θ1ρn1+θ2ρn2)ρnα)=0,for % α=1,2.

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) ∂tρn,mα+χα∂x(an,m(θ1ρn,m1+θ2ρn,m2)ρn,mα)=0,for α=1,2,

where is given by

 an,m(ρ):=∫R∂xKn,m(x−y)ρ(t,dy).

Applying Theorem 1.1 in  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:

 ddt(∫Rρn,mα(t,dx)ϕ(x))=∫R2∂xKn,m(x−y)(θ1ρn,m1(t,dy)+θ2ρn,m2(t,dy))ρn,mα(t,dx)∂xϕ(x).

Then,

 (4.3) ∣∣∣ddt(∫Rρn,mα(t,dx)ϕ(x))∣∣∣≤∥∂xϕ∥L∞∥∂xK∥L∞(θ1+θ2).

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

 ρnα⇀ραweakly in SM,for α=1,2.

Let us find the equation satisfied by in the distributional sense. Let be in . Since satisfies (4.1) in the distributional sense, we have

 ∫T0∫R∂tϕ(t,x)ρnα(t,dx)dt+χα∫T0∫Ran(θ1ρn1+θ2ρn2)ρnα(t,dx)ϕ(t,x)dt=0.

Using Lemma 3.2, we can pass to the limit in the latter equation and obtain,

 ∫T0∫R∂tϕ(t,x)ρα(t,dx)dt+χα∫T0∫R^a(θ1ρ1+θ2ρ2)ρα(t,dx)ϕ(t,x)dt=0.

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

 ∂tρα+χα∂x(^a(ρ)ρα)=0,for α=1,2, in the sense of duality in ]0,T[×R.

Step 3 : Finite second order moment.
From Equation (2.5), we deduce that the first and second moments satisfy in the sense of distributions

 ddt(∫|x|ρα(t,dx)) =−∫sgn(x)^a(ρ)ρα(t,dx), ddt(∫|x|2ρα(t,