Splitting Schemes&Segregation In Reaction-(Cross-)Diffusion Systems

One of the most fascinating phenomena observed in reaction-diffusion systems is the emergence of segregated solutions, i.e. population densities with disjoint supports. We analyse such a reaction cross-diffusion system. In order to prove existence of weak solutions for a wide class of initial data without restriction about their supports or their positivity, we propose a variational splitting scheme combining ODEs with methods from optimal transport. In addition, this approach allows us to prove conservation of segregation for initially segregated data even in the presence of vacuum.


Introduction.
In this work we consider the following reaction-diffusion system for the evolution on an interval x ∈ Ω of two species with population densities ρ, η ≥ 0 : where χ : [0, ∞) → [0, ∞) is a C 1 super-linear function modelling nonlinear diffusion and F i and G i , i = 1, 2 model the reaction phenomena. Systems of this type appear naturally in mathematical biology. A fundamental biological phenomenon in interactions among different biological species is the inhibition or activation of growth whenever two populations occupy the same habitat. One species may promote or suppress the proliferation of the other species. In models involving cells or bacteria, the limited growth of different cell types can be attributed to volume or size constraints of the individual cells forming the different populations. The diffusive part in (1) was originally introduced in the seminal papers [30] and [8,7] and exhibits an intriguing phenomenon: segregated densities remain separated at all times. In fact, nonlinear diffusions are natural ways to include volume filling effects into mathematical biology models, see [42,20] in the case of the classical Keller-Segel system. They help to avoid blow-up in these aggregation models in a biologically meaningful way and lead generically to asymptotic stabilisation. In the absence of reactions, the system leads to the nonlinear diffusion equation (2) ∂ t σ = ∂ x (σ∂ x χ ′ (σ)) = ∂ 2 x β(σ).
with σ = ρ+η, in which χ ′ (σ) models the resistance to compression of the whole group of individuals σ. The natural assumption on (2) in order to be a diffusion equation is β ′ (s) > 0 for s > 0, possibly degenerating at s = 0, or equivalently χ ′′ (s) > 0 for s > 0. The particular relevant case of χ(s) = s 2 /2 can be understood as the mean-field limit of interacting particles with very localised repulsion, see [40,16,12,35]. Related reaction-diffusion models to (1) appear in tissue growth models where cell adhesion and volume effects are important factors determining cell sorting in heterogeneous cell populations, see [38,23,17,22], and zebrafish lateral line patterning [51]. They are also basic building bricks for a variety of cancer invasion models in the literature [26,43,47,32,6,3,4,5,29,28,45] in which the coupling with other biologically meaningful modelling factors such as extracellular matrix, enzymatic activators and other substances are taken into account. These works usually involve drift terms due to long range attraction and/or repulsion between individuals leading to related mathematical difficulties with respect to (1), see for instance [34].
The nonlinear diffusion equation (2) is well-studied, see [48], and it can be understood as a gradient flow, see [41,2,49,44], of the energy functional in the metric space of probability measures endowed with a suitable topology induced by the euclidean transport distance denoted by d 2 . The wellposedness of solutions to the nonlinear diffusion equation (2) was obtained in [41] by means of the so called JKOscheme, cf. [33], which is a particular case of the minimising movement scheme by De Giorgi, see [1] and the references therein. The idea of such a scheme is to recursively construct a sequence by solving a minimisation problem in a certain metric space (X, d). Given some initial conditionσ for Eq. (2) and a fixed time step 0 < τ < 1 we set σ 0 τ =σ, and then recursively define for n ∈ N. The seminal work of R. J. McCann [37] shows that E(σ) is displacement convex or geodesically convex on the metric space of probability measures on the line endowed with d 2 as soon as β is nondecreasing on (0, +∞), called the McCann's condition in short. We also refer to [2,Chap. 9], [24, p. 26], or [50,Chap. 17] for this classical notion, and [13] for related issues. Upon choosing a proper time interpolation σ τ , it can be proven that the sequence {σ τ } τ converges to a weak solution of Eq. (2), see [49,2,50,44] and the references therein.
In this work, we propose a variational splitting scheme in order to construct weak solutions to the reaction-diffusion system (1). More precisely, we solve in an inner time stepping the diffusive part of the system by the JKO scheme related to the nonlinear diffusion equation (2), and then we transport both densities ρ, η through the flow generated by the equation for the total population σ. Note that in this step the total and individual masses of the populations are unchanged in time. In the second inner time stepping of the splitting scheme, we solve the system of ODEs parameterised by the spatial variable x ∈ Ω leading to the final approximation of our new population densities after a time step. This variational splitting scheme will be written in details in Section 2. The splitting between reaction and diffusion steps is natural from the numerical analysis view point as it has already been used for variations of Keller-Segel models where the diffusion step is solved by the JKO scheme [11,25] in the case of a single population density coupled with a system of reaction-diffusion equations.
Our main result shows the convergence of the splitting variational scheme towards weak solutions of the system (1). The main mathematical difficulty here arises from the cross-diffusion term allowing for segregation fronts to form in the solutions. This phenomenon was proven in [8] in the case of initial data with separated supports for the populations. More precisely, while [30] constructs a source solution of the system without reactions similar to the well-known Barenblatt-Pattle profiles [48] for nonlinear diffusions, [8] constructs a solution to the system without reactions by formulating it as a free boundary problem for a single effective equation, and by characterising the segregation front through this free boundary. This approach can only work in case the support of both populations are at a positive distance to each other initially. Later [6,9] combined both the nonlinear diffusion and the reaction to obtain a system similar to (1) showing similar segregation phenomena by regularisation techniques. However, their approach heavily relies on the absence of vacuum as they assume that σ 0 is bounded below by a positive constant.
These remarkable results have severe consequences -initially smooth solutions lose their regularity when both densities meet each other. In fact, they become discontinuous at the contact interface immediately. This phenomenon legitimises our functional space choice as bounded functions of bounded variation, see the precise notion of weak solution and assumptions on the initial data in the next section.
In contrast to [8,6,9], we show the convergence of our variational splitting scheme for general initial data even in the presence of vacuum and for general nonlinearities. Moreover, we recover their result in [8,6] about segregation fronts by showing that initial data which are initially segregated remain segregated for all times. In fact, we are even able to drop their restrictive assumption of strict boundedness away from zero, i.e. absence of vacuum. An important technical point in our proof relies on displacement convexity of an auxiliary functional that allows to obtain further regularity on the approximate solutions in order to pass to the limit in the nonlinear diffusion terms. This auxiliary functional imposes a slightly more restricted set of nonlinear diffusions satisfying some integrability condition at the origin, see the precise conditions in the next section.
The rest of the paper is organised as follows: In Section 2 we introduce the variational splitting scheme, present the main result, and explain the strategy of the proof. Section 3 is dedicated to deriving all estimates necessary for proving the existence theorem as well as the segregation theorem, and finally, in Section 4 we conclude by illustrating the result with some numerical examples.
2. Preliminaries and main result. As already mentioned, our main aim is to study the existence of weak solutions for the following one-dimensional two species cross-diffusion and reaction system: in Ω, where Ω ⊂ R is an open bounded interval, T > 0. Moreover, χ denotes an internal energy density and F i and G i , i = 1, 2 model the reaction phenomena. As mentioned before the space of bounded functions with bounded variation is a natural functional setting.
Definition 2.1 (Space of functions of bounded variations). Let f :Ω → R. We define its variation with respect to a partition P := {x 1 < x 2 < · · · < x |P | } ⊂ Ω by We call f a function of bounded variation if its total variation sup P V P (f ) < ∞ is finite. Here the supremum is taken over all partitions of Ω. We denote by BV (Ω) the set of functions whose variation is bounded. Equipped with the norm the set BV (Ω) is a vector space.
We shall see in the remainder of this section that the vector space of bounded function with bounded variation is a good choice to construct solutions. In our analysis we will exploit the following property.
is a real algebra.
The proof of the previous result is standard. Notice that in one dimension, BVregularity implies boundedness. We prefer to write BV (Ω) ∩ L ∞ (Ω) for the sake of clarity and possible future generalisations.

Metric structure.
Consider Ω ⊆ R an open bounded interval, and denote by M + (Ω) the set of positive and finite measures. Throughout this paper we will make use of the following notation P m (Ω) := µ ∈ M + (Ω) µ(Ω) = m , that is, the set of positive measures with mass m > 0. Consider a measure µ ∈ P m (Ω) and a Borel map T : R → R. We denote by ν = T # µ ∈ P m (Ω) the push-forward measure of µ through T , defined by for all Borel functions f on Ω. We call T a transport map pushing µ to ν. We endow the space P m (Ω) with the p−Wasserstein distance, p ≥ 1, Here, Π m (µ 1 , µ 2 ) is the set of all transport plans between µ 1 and µ 2 , that is the set of positive measures of fixed mass, γ ∈ P m (Ω × Ω), defined on the product space such that π i # γ = µ i , for i = 1, 2, where π i denotes the projection operator on the i-th component of the product space. If µ 1 is absolutely continuous with respect to the Lebesgue measure the optimal transport map, γ, is unique and can be written as γ = (id, T ) # µ. In addition there exists a Kantorovich potential, ϕ, that is linked to the transport map, T , in the following way We refer to [49,2,50,44] and the references therein for a good account of the properties of transport distances and the state of the art in gradient flows/steepest descents of functionals in metric spaces of probability measures. While transport distances are an incredibly powerful tool for dealing with transport PDEs exhibiting a gradient flow structure, it is not applicable in the presence of source terms. This is owing to the fact that it is only defined for two measures of the same mass.
To resolve this shortcoming we will make use of the Bounded-Lipschitz distance d BL , classically used for the derivation of the Vlasov equation, see [39,46,19,18,31] and the references therein. The Bounded-Lipschitz distance d BL , also frequently called flat metric, is defined as follows Since our problem is posed on the product space we extend the metric setting to the space M + × M + , the product space of nonnegative measures in the canonical way.
(Ω) for i = 1, 2 as well as µ 2 (Ω) = µ 3 (Ω) and ν 2 (Ω) = ν 3 (Ω). Then Throughout two quantities are crucial for our analysis -the sum of the two densities, σ, and the ratio, r, between one density, say ρ, and the sum where we assume no vacuum, i.e. σ > 0. A straightforward computation shows these functions satisfy the following system of PDEs where to denote the reaction terms in the transformed variables. In order to simplify the analysis in Section 3, let us introduce the more concise notation Note that these functions are Lipschitz and bounded as they are linear combinations of BV ∩ L ∞ functions. Thus the transformed system (7) can be rewritten in the more compact form Let us note here, that taking F i = G i = 0 in (1), yields the following nonlinear cross-diffusion system where the sum σ satisfies the nonlinear diffusion equation (2).
In order to be useful for our purposes, we need some properties on the internal energy density. The role of these properties will be clearer after the statement of our main result.
As mentioned in the introduction, the space of probability measures endowed with d 2 has proven to be an exceptional choice of a metric space. In this case the minimiser in Eq. (3) satisfies the so-called optimality condition cf. [44,Proposition 7.20], for instance. Here ϕ denotes the associated Kantorovich potential [44,Theorem 1.17]. Notice that this optimality is only obtained in those references for non-degenerate problems, where β ′ (0 + ) > 0. However, similar approximation techniques as those developed in [41, p. 156] and [24, p. 27] allow to overcome the difficulties associated to degenerate diffusions, where β ′ (0 + ) = 0. In the rest of the paper, we will proceed as if we were dealing with nonlinearities leading to nondegenerate diffusions since by this standard approximation procedure, the same result can be obtained for the degenerate ones. We are now ready to introduce our notion of weak solutions.
> 0, and is therefore displacement convex. This fact will be crucial in Section 3.3 in order to prove Lemma 3.10. Let us just state here that it is necessary to obtain additional regularity from the dissipation of this functional on the nonlinear diffusion term, thus allowing us to pass to the limit in the approximating sequence.
2. Let us note that we can also allow for nonlocal reaction terms, i.e.
and similarly for G i , for i = 1, 2. The only assumption we need to impose on the kernels is that they are smooth and integrable. These models are found in modelling pattern formation, as for instance the kernel-based Turing pattern system [36] or the nonlinear aggregation-diffusion system [51]. 3. Similarly, we can -at least formally -interpret the system (9) as a gradient flow of the functional in the product Wasserstein space.

Splitting scheme.
We are now ready to introduce our splitting scheme for equation (4). Let some initial data ρ 0 , η 0 ∈ BV (Ω) ∩ L ∞ (Ω) be given such that there exists a function r 0 ∈ BV (Ω) such that and 0 ≤ r 0 ≤ 1. Furthermore we assume F i and G i , i = 1, 2 are bounded and Lipschitz with respect to ρ and η and we impose G 1 (0, ·) ≥ 0 and G 2 (·, 0) ≥ 0 to ensure positivity of solutions.

Combination of both steps -construction of a solution.
Throughout this paper we refer to Eq. (11a) as reaction step and to Eq. (11c) as diffusion step, respectively.
We will say that two densities ρ, η ∈ BV (Ω) ∩ L ∞ (Ω) are segregated if the intersection of the interior of their supports is empty. We are now ready to state our main result.
3. Proof of the main result. This section is dedicated to proving the main result of the paper -the convergence of the approximation obtained by the splitting scheme to a solution of the system. It is organised as follows: in Section 3.1 and 3.2 we establish the crucial BV-estimates and L ∞ -bounds. In Section 3.3 we combine the estimates from the previous sections in order to get uniform estimates for a whole iteration. Finally, in Section 3.4 we show how to extract a convergent subsequence and identify its limit as a weak solution to system (4).

3.1.
Estimates for reaction step. Since the right-hand sides, Σ(σ, r), R(σ, r) are Lipschitz continuous in both components we note that the solution of the reaction system is unique.
Next we address the BV-estimates during the reaction step.
Proposition 3.2 (Bounded variation of r n+1/2 and σ n+1/2 ). Let us consider (r n , σ n ) as initial data for our splitting scheme. Then, the reaction step is BV-stable in the following sense.
for some positive K, depending only on the Lipschitz constants of F i , G i and the L ∞ -bounds on F i , G i , for i = 1, 2.

Upon integrating in time we get
Now, let P ⊂ Ω be an arbitrary partition. We compute the variation of σ and r with respect to P and obtain where c only depends on the L ∞ -bounds and the Lipschitz-continuity of A i , for ∈ {1, 2, 3} and the L ∞ -bounds of σ and r. Applying Gronwall's lemma we finally obtain Passing to the supremum on the right-hand side and then on the left-hand side yields the result.

Estimates for diffusive step.
This section is devoted to establishing BVestimates and L ∞ -bounds in the diffusive step. To this end we will make use of the following remark.
Proof. Choose x 0 ∈ argmax(σ n+1 ). Then, by the optimality condition, Eq. (10), we have x 0 ∈ argmin(ϕ) where we used the fact that χ ′′ ≥ 0, cf. Definition 2.5, (NL − iii). Hence ϕ ′′ (x 0 ) ≥ 0 and consequently, by passing to the derivative in Eq. (5) we get where T is the transport map from σ n+1 to σ n+1/2 . After a change of variables we get for any x ∈ Ω. For the non-negativity we observe that T ′ ≥ 0 by Brenier's theorem [14,15,49,50]. Thus Finally the bounds for r n+1 follow from its definition, cf. Eq. (11d), as the composition with a monotone function does not change the infimum and the supremum of a function. This concludes the proof.
Proof. The result for the BV-norm of the minimiser, σ n+1 , is shown analogously to the proof of Theorem 1.1., cf. [27]. Now we need to show that the BV-norm of the ratio r does not increase. Recall the definition of r n+1 , cf. Eq. (11d), as where T is the transport map such that ρ n+1/2 = T # ρ n+1 and σ n+1/2 = T # σ n+1 . Note that it is indeed the same function and there holds T ′ ≥ 0, by Brenier's theorem [14,15,49,50] . Now, let P ⊂ Ω be any partition of Ω. There holds where P ′ is another partition induced by the monotone map T . Taking the supremum over all partitions P , we get Finally note that, indeed, on supp(σ n+1 ) as we shall see now. According to Remark 3.3, the same transport map T pushes ρ n+1 onto ρ n+1/2 and σ n+1 onto σ n+1/2 . As a consequence the densities satisfy 3.3. Combined estimates for an entire splitting step. We have now garnered all information necessary to pass to the limit. Let us combine the estimates from the previous section in the following lemma. Proof. The uniform L ∞ -bound is a consequence of combining Propositions 3.1 and 3.4. We use these uniform L ∞ -bounds in the estimates for the BV -norm, cf. Propositions 3.2 and 3.5. Combining both for the reaction and diffusion step we also obtain the uniform BV -bounds.
As a result of Lemma 3.6 and Lemma 2.2 we obtain the following corollary. Finally, we need to prove an estimate on the cross-diffusion term to be able to pass to the limit later. This estimate is achieved in Lemma 3.10 which is preceded by two technical lemmas. We exploit the existence of an auxiliary functional guaranteed by Definition 2.5. Note that in the absence of the reaction part this would indeed be an entropy in the classical sense, i.e. that it is decayed along solutions. Since we are interested in a uniform estimate we shall begin by proving a control of this functional during each reaction phase. Proof. A straight forward computation yields using the uniform L ∞ -bound on σ and the fact that |Ω| < ∞. Hence, where c is independent of n.
Next, we address the diffusion step. As mentioned earlier the auxiliary functional, K, is an entropy for the diffusive part and from its dissipation we obtain the necessary regularity, as asserted in the following lemma. Lemma 3.9 (H 1 -bound for σ n+1 ). The minimiser of the JKO step satisfies the following estimate Let σ s = (T s ) # σ n+1 be the geodesic interpolation between σ s | s=0 = σ n+1 and σ s | s=1 = σ n+1/2 , given by for the associated Kantorovich potential, ϕ, cf. Eq. (5). As a consequence the velocity field is given by v s = (T − id) • T −1 s , satisfying the following continuity equation, We differentiate the entropy along the geodesic and obtain Thus, at s = 0, the evolution of the entropy Eq. (8) becomes Using the optimality condition, Eq. (10), we obtain where the last inequality is a consequence of the geodesic convexity of the entropy, cf. Remark 2.7. This concludes the proof.
We combine the previous lemmas to obtain the desired estimate for a full iteration and finally for the piecewise constant interpolation, σ τ .
for some positive constant depending only on T .
Proof. This statement is a consequence of combining Lemma 3.8 and Lemma 3.9 to get Summing over n = 0 . . . N − 1 gives which yields the statement of the lemma.
Lemma 3.11 (Total-square 2−Wasserstein distance estimates). For every n ∈ {1, . . . , N } consider two consecutive steps for (11c), U n+1/2 = (ρ n+1/2 , η n+1/2 ) and U n+1 = (ρ n+1 , η n+1 ), then there exists a constant C such that Proof. Using the minimising property of U n+1 we have Adding and subtracting E(U n ) on the right-hand side, and considering that due to a Taylor expansion of χ, where ξ ∈ [σ n , σ n+1/2 ]. Using the L ∞ -bounds on σ, r and F i , G i , we obtain for a full time step. Finally, summing over n gives the result: 3.4. Convergence. We now prove the convergence result. Proof. The proof is based on the application of a generalised version of the Ascoli-Arzelà theorem, cf. Ref. [2], Section 3.
We begin by establishing 'almost continuity' of the approximation. To this end let 0 ≤ s < t ≤ T be two time instances. Then there exist two uniquely determined integers, m, k, satisfying where U τ (t) = (ρ τ (t), η τ (t)) as defined in Definition 2.8. It becomes apparent that we need to address the bounded Lipschitz term next. To this end we use the triangulation established in Corollary 2.4 to estimate it by the L 1 -distance in the reaction step and the W 1 -distance in the diffusion step. Hence For the reaction step an argument similar to Eq. (14) yields U n − U n+1/2 2 L 1 (Ω) ≤ Cτ 2 , and, using that the p-Wasserstein distances are ordered, we even have where we used the total-square estimate, Lemma 3.11. Therefore Eq. (15) becomes Thus we get the 'almost 1 2 -Hölder continuity' for the curve U τ (t) and we obtain the uniform narrow compactness on compact time intervals by using the refined version of Ascoli-Arzelà's theorem, cf. [2], Section 3.
Proof. Note that it suffices to show the result for (ρ τ ) τ >0 as the same argument applies for (η τ ) τ >0 . By Proposition 3.12 we can extract a subsequence, still denoted the same, such that  Proof. The proof consists of two parts -the diffusion part and the reaction part. We write them in their respective weak formulation and combine them to obtain the complete approximation of the weak formulation. Here we only show the argument for ρ as the corresponding result for η is shown analogously. We begin with the diffusion part.
Multiplying Eq. (18) by τ and summing from m to k − 1, we obtain We are now ready to pass to the limit. The strong convergence obtained in Corollary 3.13 allows us to pass to the limit in the nonlinear reaction terms. Moreover the crossdiffusion term converges due to the weak-strong L 2 (0, T ; L 2 (Ω))-duality, by Lemma 3.13 and Corollary 3.13. Passing to the limit τ → 0 we obtain the weak formulation for ρ. The same argument for η yields the statement.
We end the paper with a stunning result which can be seen as a generalisation of the result of Bertsch et al., cf. [8,7]. In their papers they prove that initially segregated species stay segregated at all times. We can drop their assumption that i.e. that both species are ordered and prove the following, more general theorem. Then, there exists a solution such that both species stay segregated at all times.
Proof. It suffices to show this property at the level of the discrete scheme since, once it is established, we use the strong L 2 -convergence of ρ τ and η τ to show segregation is also kept in the limit: Thus let us now show the property at the level of the approximation. It is clear that in the reaction step segregation is kept as it does not change the support of both ρ and η. Hence we only need to prove that segregation is kept in the diffusion step. This is done by contradiction. Let us assume there exists an 1 ≤ n ≤ N such that while Ω ρ n+1 (x)η n+1 (x) dx > 0.
4. Numerical simulations. In this section we illustrate our main result by showing some numerical simulations using a different splitting scheme. This splitting scheme is based on the finite volume scheme introduced in [10,21,23]. We solve the cross-diffusion inner time step with this finite volume scheme instead of the optimal mass transport approach. Then we solve the reaction step along the ODEs. We refer to [25] for the numerical difficulties related to the implementation of the mass transport approach in a splitting fashion.   4.1. Lotka-Volterra type reaction. Let us consider the system ∂ t ρ = ∂ x (ρ∂ x (σ)) + ρ(1 − σ), modelling two species avoiding overcrowding due to the nonlinear cross-diffusion term as well. In addition the growth term is of Lotka-Volterra type, cf. Figures 1 & 2. Let us consider now a little modification of system (20), as in [6], Using the same parameters and domain as in [6] Figures 3 & 4 show complete segregation in agreement with their result and moreover numerically validate our main result in Theorem 3.15.