FINITE ELEMENT APPROXIMATION OF A POPULATION SPATIAL ADAPTATION MODEL

In [18], Sighesada, Kawasaki and Teramoto presented a system of partial differential equations for modeling spatial segregation of interacting species. Apart from competitive Lotka-Volterra (reaction) and population pressure (cross-diffusion) terms, a convective term modeling the populations attraction to more favorable environmental regions was included. In this article, we study numerically a modification of their convective term to take account for the notion of spatial adaptation of populations. After describing the model, in which a time non-local drift term is considered, we propose a numerical discretization in terms of a mass-preserving time semi-implicit finite element method. Finally, we provied the results of some biologically inspired numerical experiments showing qualitative differences between the original model of [18] and the model proposed in this article.

1. Introduction.In [18], Shigesada, Kawasaki and Teramoto introduced a model for analyzing the spatial segregation patterns arising in the evolution of populations of two species which are ruled by • competition for similar resources, • population pressure, and • environmental quality.These biological interactions are realized mathematically in the form of a time evolution drift-cross diffusion system of partial differential equations, for i = 1, 2, where u i denotes population density, U is the environmental potential, modeling areas where the environmental conditions are more or less favorable [18,17], the non-negative diffusion coefficients c i and a ij model the random diffusion of individuals and the inter-and intra-specific population pressure, respectively, and d i are usually assumed to be real constants determining the attraction strength of the potential maxima.Function f i is a competition Lotka-Volterra type function, where α i ≥ 0 is the intrinsic growth rate of the i−species, and β ij ≥ 0 are the coefficients of inter-and intra-specific competition.
In this article, for simplicity, we assume Equations (1) to be satisfied in the bounded domain Ω × (0, T ), with Ω ⊂ R an open interval, and T > 0, although the multi-dimensional case Ω ⊂ R N with N ≤ 3 could be also treated.The problem is completed by prescribing non-flux boundary conditions and non-negative initial data: for i = 1, 2, where ν denotes the exterior unit normal to ∂Ω.Problem ( 1)-( 3) has received much attention since its introduction due to the interesting spatial pattern formation of solutions, referred to as segregation.These patterns do not arise in the linear diffusion model, i.e. for a ij = 0, i, j = 1, 2, where, if in addition d i = 0, then the steady state solutions are constants determined by the zeros of the Lotka-Volterra terms.These constant solutions correspond, in general terms, to two kind of competitions: weak, which implies coexistence, and strong which implies extinction of at least one population.Lou and Ni [15,16] analyzed the steady state problem corresponding to (1) (with d i = 0) and proved the existence of non-constant solutions for some parameter combinations including weak and strong competition.Their results seem to indicate that while the intensity of diffusion (c i ) and self-diffusion (a ii ) tend to suppress pattern formation, those of cross-diffusion (a 12 , a 21 ) seems to help create segregation patterns.
The first results on the existence of solutions of problem ( 1)-( 3) were proven under certain restrictions on the self-and cross-diffusion coefficients.For instance, for sufficiently small cross-diffusion terms (or small initial data) and vanishing selfdiffusion coefficients a 11 = a 22 = 0, Deuring proved the global existence of solutions in [4].For the case c 1 = c 2 a global existence result in one space dimension was obtained by Kim [14].Furthermore, under the condition 8a 11 > a 12 , 8a 22 > a 21 , ( Yagi [19] showed the global existence of solutions in two space dimensions assuming a 12 = a 21 .A global existence result for weak solutions in any space dimension under assumption (4) can be found in [8].Condition (4) can be easily understood by observing that in this case, the diffusion matrix is positive definite, hence yielding an elliptic operator.If the condition (4) does not hold, there are choices of c i , a ij , u i ≥ 0 for which the diffusion matrix is not positive definite.In [9] (see also [10] for some extensions of the result) the existence of global weak solutions for any a ij > 0 was proven by using a suitable entropy functional.However, the proof uses the embedding H 1 (Ω) ⊂ L ∞ (Ω) in a crucial way, restricting the result to one space dimension.The one-dimensional result was later generalized by Chen and Jüngel [2] to up to three space dimensions without any additional restrictions than those given in [9].We refer to [15,16] for the corresponding stationary problem and notice that related models appear, among other fields, in ecology, chemotaxis, granular material and semiconductor theories [5,12,7,11,3].On the numerical side, a first approach based on a time Euler semi-discrete scheme was proven to be convergent in [9], in the one dimensional spatial case.More recently, other numerical approaches have been introduced in the context of Euler-Galerkin approximations (N ≤ 3) by Chen and Jüngel [2], of finite element methods (N ≤ 3) by Barret and Blowey [1] and of particle methods (N = 1) by Gambino et al. [13].However, for the 1D case, the results of all these methods seems to be similar.
In all these works the focus is set either on finding conditions on the diffusion coefficients c i and a ij , and the Lotka-Volterra coefficients α i and β ij which ensure the existence of solutions of problem (1)-(3), allowing to define a convergent numerical method to approximate them, either on finding conditions which imply qualitative properties such as the co-existence or the extinction of populations in the steady state problem.However, the drift term responsible of directing the populations towards the maxima of the environmental potential has been always assumed to be linear and depending on constant coefficients (d i ) which express the strength of attraction of the different populations to these maxima.In [6] we introduced a new drift coefficient dependence which allow us to modelize spatial adaptation by means of a memory mechanism which strengthen the attraction of population to a point if the population density in such point has been high in the past.More explicitly, we assume d i to have the form In addition, another modification of the original problem was introduced in [6], assuming the intrinsic growth rate coefficients α i to be non-constant.Indeed, once that we consider an heterogeneous space domain in which populations are driven to the environmental potential maxima, it seems reasonable to assume a dependence of α i on U such that larger growth of populations takes place in better environmental regions.We set with α i ≥ 0. From the analytical point of view, both modifications may be treated introducing minor changes in the proof of Theorem 1 of [9] (see also [2,1], for N ≤ 3) allowing to prove the existence of weak solutions of problem ( 1)-(3), with α i defined by (6) and d i by (5), see [6] for details.
The main contribution in this paper relates to the numerical method employed to simulate particular solutions of Problem (1)-(3) under assumptions ( 5) and ( 6).In [6], a finite differences scheme for the spatial discretization was considered.Although, as mentioned, in the one-dimensional case we did not find any decisive practical advantage on using any particular method of approximation, the finite differences discretization introduces some complications involving the boundary data, specially for variable coefficients α i , implying that a discrete version of (4) must be assumed to hold on the boundary to ensure the uniquenes of the discrete solution.However, as it is well known, non-flow boundary conditions are implemented in a simple way in the finite element method, not requiring any additional assumption on the coefficients of the problem.Moreover, the advantages of the finite element method against the finite differences method are well known in the space multidimensional case.

Numerical discretization and examples.
In this section we present numerical simulations illustrating differences between the behaviors of solutions corresponding to constant or variable convective coefficients d i .In the first example, of qualitative nature, we simulate the situation in which a catastrophic natural event changes abruptly the spatial location of the maxima of the environmental potential, i.e., the more favored environmental region.We see that in the case of spatial adaptation, represented by d i given by ( 5), the extinction of the population more intensely adapted to the initial potential maximum is possible as a result of the very low population density left in areas far from this maximum, which implies a bad competitive positioning near the new potential maximum after the catastrophic event.However, for constant d i the population is able to recover and dominate again in the new favored region.
The second example, in which the potential maximum is kept time independent, shows two interesting biological properties.First, that the segregation of populations is more intense in the case of spatial adaptation than in the case of constant convective coefficients.Second, that intense spatial adaptation may lead to coexistence in cases where the constant convection coefficients lead to extinction.

Finite element approximation.
As mentioned in the Introduction, we take the spatial dimension to be n = 1, e.g.Ω ⊂ R is chosen to be an open interval.For a more detailed description of the following finite element approximation, including the case of higher spatial dimensions, see [1].On the interval Ω, we consider a family of quasi-uniform partitionings T h , h > 0, consisting of disjoint and open subintervals I with h I = |I| and h = max I∈T h h I , so that ΩD = ∪ I∈T h Ī. Associated with T h is the finite element space Let J be the set of nodes of T h and {p j } j∈J the coordinates of these nodes.Let {ϕ j } j∈J be the standard basis functions for S h , that is ϕ j ∈ S h , ϕ ≥ 0 in Ω, and ϕ j (p i ) = δ ij for all i, j ∈ J.The following functions were considered in [1] to obtain a discrete analogue of the entropy inequality which allows to control the possible non-positivity of discrete approximate solutions.We define Λ ε : with F ε : R → [0, ∞) given by and For the time discretization, we consider a partitioning 0 = t 0 < t 1 < . . .< t N −1 < T N = T of [0, T ] into possibly variable time steps τ n = t n − t n−1 , n = 1, . . ., N .We set τ = max n τ n .For any given ε ∈ (0, 1), we then consider the following finite element approximation of problem ( 1)-( 3 where π h : C( ΩD ) → S h is the usual interpolation operator, with (π h η)(p j ) = η(p j ) for all j ∈ J, U n stands for U (•, t n ), and u 0 ε,i ∈ S h is an approximation of u 0 i , for instance its L 2 projection on S h .Trivial modifications of the proof of Theorem 2.1 of [1], allows to ensure the existence of solutions of problem (7).More concretely, if ) and τ n is small enough, then there exists a ) ∈ (S h ) 2 to the n-th step of problem (7).In addition, if τ → 0 with either τ 1 ≤ Ch 2 or u 0 i ∈ H 1 (Ω), and if εh −1/2 → 0 then a subsequence (not relabeled) of may be extracted such that (u ε,1 , u ε,2 ) → (u 1 , u 2 ) in a suitable sense, being (u 1 , u 2 ) a weak solution of problem ( 1)-(3).

Experiments.
As in [1], we use the following fixed point algorithm for solving the system of nonlinear algebraic equations for (u n ε,1 , u n ε,2 ) arising at each time level from the approximations (7).For t = t 0 = 0, set u 0 i = u 0 ε,i .For t = t n , let u n−1 ε,i be given and set u n,0 ε,i = u n−1 ε,i .Then, for k ≥ 1 find u n,k ε,i such that for i, j = 1, 2, with j = i, and for all ϕ ∈ S h ) ∇u n,k ε,i , ∇ϕ We adopted the stopping criteria max i=1,2 with tol = 10 −7 in the experiments, and set In all experiments we integrated in time until a numerical stationary solution, u S i , was achieved.This was determined by max Unless otherwise stated, in all the experiments we use the data given in Table 1.

Experiment 1. Intensive adaptation may lead to extinction after a catastrophic environmental event.
In this example we explore the effects that sudden environmental changes may have on the extinction of populations which have adapted intensively to some region.For this experiment we use the following Lotka-Volterra functions with α 1 (U ) = 320U , α 2 (U ) = 300(0.99U+ 0.01), and β ij = 150 for i, j = 1, 2. For the case of variable d i , we define them by (5).When running for constant d i , we take which is constant due to the election of the initial datum, see Table 1.The environmental potential, given also in Table 1, is first set with the maximum at x 0 = 0.2.We run the simulation until t = 0.3 is reached and a sudden change of the potential maximum, to x 0 = 0.8, is produced.Then we continue till the steady state is nearly reached, which we assume to be when For the case of constant d i this happens at t = 4.29 while for d i variable it takes till t = 12.86.Notice that the only difference between equations and data for populations 1 and 2 is the definition of the growth rate coefficients α i .
In Fig. 1 we show time slices of the evolution (left to right and up to down) of both populations (u 1 continuous line, u 2 dotted line) for the case of variable d i (first row) and constant d i (second row).In the first column, we plotted the population distribution just before the sudden change of maximum of the potential that takes place in t = 0.3.A notorious concentration and growing of population 1 in the neighborhood of x 0 = 0.2 is accompanied by an almost extinction of population 2, due to α 1 > α 2 , for both choices of d i .However, although hardly visible from the plots, we checked that u 2 > u 1 in regions far from the potential maximum.After the catastrophic event, at t = 0.3, the behavior of populations differs absolutely, leading to steady states in which the dominance is switched (d i variable) or mantained (d i constant).
In Fig. 2 we plot the spatial adaptation terms produced in the case of variable d i , which rule the strength of the convection term.It also interesting to note that the biological notion of spatial adaptation represented by these terms is time varying and that intense adaptation to some region, x 0 = 0.2 in this example, may be weakened and practically disappear if the region is left uninhabitated for a long time.
In Figs.1-2, the approximation to the steady state solution satisfying ( 8) is reached around time t = 3. Experiment 2. Adaptation in a stable environment may enhance segregation and promote coexistence.In this example we compare the segregation magnitudes for the cases of variable and constant d i .We run two experiments, one with zero Lotka-Volterra functions, implying mass conservation for the continuous model, and another with similar Lotka-Volterra functions than those of Experiment 1.The environmental potential, given in Table 1, is time independent and with the maximum at x 0 = 0.5, to check the symmetry preserving property of the discretization scheme.In order to have some distinction between populations, we set, for the case of variable d i , with ε 1 = 2 and ε 2 = 1.The convection coefficients for the corresponding problem with constant d i are given by d i = ε i u i0 .We run the simulation until the steady state is nearly established, using the same criterium than in Experiment 1.Notice that the only difference between equations and data for populations 1 and 2 is in the parameters ε i .
In Fig. 3 we show the steady state for constant d i (left) and variable d i (right), which is reached for t * ≈ 0.714 and t * ≈ 8.186, respectively.Continuous line corresponds to population 1 and dotted line to population 2. The mass conservation property is well captured by the discrete model, being the relative difference lower than 8 × 10 −4 for d i constant and an order greater for d i variable.The symmetry of the solution is also conserved, being the relative difference max i=1,2 of the order of 10 −15 , for both d i cases, where ũi (x, t) = u i (1 − x, t).The concentration of mass in the neighborhood (0.45, 0.55) of x 0 = 0.5 is, for the variable case 0.55 0.45 where the percentage points are in terms of the total mass, while, for the constant case is The qualitative differences between the solutions of constant or variable d i are clearly seen in Fig. 3.In the case of variable d i , the population with more intensive adaptation capacity (population 1 due to ε 1 > ε 2 ) concentrates in the surroundings of the potential maximum x 0 while population 2 reaches, in fact, its minimum density at this point.However, in the case of constant d i , both populations reach their maximum value at the potential maximum x 0 = 0.5.
For the second example of this experiment we use the Lotka-Volterra functions given in (9) with α i (U ) = 30U , for i = 1, 2 and β ij = 15 for i, j = 1, 2. In Fig. 4 we show the steady state for constant d i (left) and variable d i (right).Continuous line corresponds to population 1 and dotted line to population 2. Notably, the effect of rapid concentration of population 1 around the maximum, for the case of variable d i , leads the system to a coexistence steady state, although with a high degree of segregation in the surroundings of x 0 , where population 2 attains a minimum.However, for the case of constant d i , population 2 is extincted.Experiment 3. Preservation of segregation.In this example we consider a environmental potential Ũ with two maxima located at x 0 = 0.2 and x 1 = 0.8, given by Ũ (x) = 0.3U (x − x 0 ) + 0.7U (x − x 1 ) for x ∈ (0, 1), with U given in Table 1.We use the Lotka-Volterra functions given in (9) with α i (U ) = 2U , for i = 1, 2 and β ij = 1.5 for i, j = 1, 2. In addition, the initial distributions of populations is assumed to be segregated one from each other and occupying a region containing one of the maxima.More explicitly, we take Observe that the differential equations satisfied by u 1 and u 2 are formally the same, so that the only difference among them must be due to the asymetrical initial data.
In Figure 5 we see the results of considering variable d i (left) and constant d i = 0.5 (right) defined as in the previous experiments.In both cases, the formal steady state takes longer than in the previous examples to be achieved The most interesting phenomenum arising in this example is the preservation of segregation in the case of variable d i , i.e., in the case in which spatial adaptation takes place.This is explained by the fact that nonlinear non-local transport effects drive intensely each population to the local maximum they initially occupy, weakening in this way the transport to the other maximum.However, in the case of linear transport the final distributions are almost identical due to the constant attractive strength populations feel towards both maxima.

3.
Conclusions.The mathematical model ( 1)-( 3) introduced by Shigesada, Kawasaki and Teramoto (1979) [18] to reproduce the behaviour of interacting species which are affected not only by competition or random displacement but by population pressure and attraction to favorable environmental regions has been a source of interesting mathematical and biological discussion.One of the more remarkable properties of the model is the formation of segregation patterns, which are observed in the field, and which lead to non-trivial steady state configurations.In this article we proposed the consideration of a new term in the equations, the time non-local term (5), which may be interpreted as a spatial adaptation intensity or capacity of the populations.From the mathematical point of view, the introduction of this new term do not pose additional difficulties for the achievement of results on existence and regularity of solutions.Moreover, the numerical discretization seems to behave well in a similar range of parameters and data than the original model since the way in which the new term induces the concentration of one of the populations in a narrow region is always bounded.However, the introduction of this new term produces important quantitative and qualitative differences with respect to the original model.We showed numerical experiments in which the behavior of solutions is qualitatively different in terms of coexistence and extinction.We also showed that the segregation-concentration effect already present in the original model is enhanced quantitatively with the introduction of the spatial adaptation term.

Figure 1 .
Figure 1.Experiment 1. Population evolution for variable di (above)and constant di (below).Time slices just before the sudden change of environmental potential maximum at t = 0.3, which relocates from x = 0.2 to x = 0.8, and for the computed steady state.The horizontal lines correspond to the initial populations densities.

Figure 2 .
Figure 2. Experiment 1. Evolution of the adaptation intensity terms, di.Above: time slices before the sudden change of environmental potential maximum at t = 0.3, which relocates from x = 0.2 to x = 0.8.Below: time slices after change of environmental potential.

Figure 3 .
Figure 3. Experiment 2. Lotka-Volterra terms set to zero and environmental potential maximum at x0 = 0.5.Case of variable di (left) and constant di (right).The horizontal line corresponds to the initial populations densities.Notice the change of scale between figures.

Figure 4 .Figure 5 .
Figure 4. Experiment 2. Competitive Lotka-Volterra terms and environmental potential maximum at x0 = 0.5.Case of variable di (left) and constant di (right).The horizontal line corresponds to the initial populations densities.Notice the change of scale between figures.

Table 1 .
Parameter values common for all the experiments