A note on the replicator equation with explicit space and global regulation

A replicator equation with explicit space and global regulation is considered. This model provides a natural framework to follow frequencies of species that are distributed in the space. For this model, analogues to classical notions of the Nash equilibrium and evolutionary stable strategy are provided. A sufficient condition for a uniform stationary state to be a spatially distributed evolutionary stable state is presented and illustrated with examples.

Arguably, one of the simplest replicator equations takes the forṁ Here v = v(t) = (v 1 (t), . . . , v n (t)) ∈ R n is a vector-function of n variables, A is a constant n × n matrix with elements a ij ∈ R, (Av) i is the i-th element of the vector Av, (Av) i = n i=1 a ij v j (t), f loc (t) is a function, which is determined later, and dot, as usual, is used to denote differentiation with respect to the time variable t. It is customarily supposed that v i (t) describes a relative abundance of the i-th species such that the total concentration n i=1 v i (t) is kept constant and often equal, without loss of generality, to 1. Therefore, the state space of (1) is the simplex S n = {v : n i=1 v i (t) = 1, v i (t) ≥ 0, i = 1, . . . , n}, which is invariant under (1) if we set f loc (t) = Av, v = n i=1 (Av) i v i , so that · , · denotes the usual scalar product in R n . Since (Av) i gives the net rate of growth of the i-th species in our system, the quantity (Av) i is termed as the Malthusian fitness, and hence f loc (t) represents the average fitness of the population at time t.
System (1) is a very well studied object, see, e.g., [8,18,17]. Very briefly, the rest points of (1) are given by the solutions to the following system (Av) 1 = (Av) 2 = . . . = (Av) n = Av, v = β, v ∈ S n . ( In general, (2) can have no solutions in int S n = {v : . . , n}, a unique solution, or infinitely many solutions. The necessary condition to have a solution to (2) is the linear dependence of the rows of A, which we denote A i , and the vector 1 n = (1, 1, . . . , 1): A natural way to derive the replicator equation (1) from the first principles is to start with a selection system of the formẏ where y = (y 1 (t), . . . , y n (t)) ∈ R n + is a vector of absolute sizes, and F i (v) denotes the Malthusian fitness (the per capita birth rate) of species y i , which can depend on the structure of the total population at the time t. Assuming that y does not tend to zero, the change of the variables v i = y i /( n i=1 y i ) leads to the replicator equation (1), if F i (v) = (Av) i . Therefore, it is equivalent to study the selection system (3) or the replicator equation (1) (see [20,21] for more details).
As it was mentioned, the replicator equation (1) naturally arises in the evolutionary game theory (for the origin, see [22,23]). There exists a parallel between the concepts of the game theory with a payoff matrix A and the behavior of the solutions to the replicator equation (1) [17]. In particular, one of the central notions of the game theory, the Nash equilibriumv, is defined as suchv ∈ S n for which v, Av ≤ v, Av for any v ∈ S n ; and an evolutionary stable state for any v =v in a neighborhood ofv ∈ S n . It can be shown (see [17] for the details and proofs) that ifv is a Nash equilibrium of the game with payoff matrix A, thenv is a rest point of (1). Moreover, ifv is a rest point of (1) and Lyapunov stable, then it is a Nash equilibrium; and ifv is ESS, then it is an asymptotically stable rest point of (1). We remark that model (1) is a system of ordinary differential equations (ODE), i.e., it is a mean-field model. A significant attention was drawn to the replicator equation in the case when heterogeneous spatial structure can be included into the model formulation [11]. One of the suggested solution was spatially explicit models (see [1,2,10] for the models of molecular evolution). In general, there are several different approaches to include spatial structure into the replicator equation. The solution to the problem when all the diffusion rates are equal is straightforward: in this case, following the ecological approach, we can just add the Laplace operator to the right hand sides of (1) (this was used, e.g., in [14,16]). However, the assumption of the equal diffusion rates would be too stringent in the general situation. To overcome this problem, Vickers et al. [8,19,28] introduced a special form of the population regulation to allow for different diffusion rates. In these works a nonlinear term is used that provides a local regulation of the populations under question to keep the total population size constant, although no particular biological mechanism is known that lets individuals adapt their per capita birth and death rates to local circumstances [13]. A more straightforward approach, in our view, would be to start with a spatially explicit selection system of the form . . , n, are diffusion coefficients. This approach is similar to the lines how the replicator equation (1) can be obtained from (3). In this way, assuming impenetrable boundary of the area Ω, we automatically obtain the condition of the global population regulation which was considered in some earlier works on the mathematical models of the prebiotic molecular evolution [3,4,5,6,29]. Here we extend this approach to the general replicator equation.
The rest of the paper is organized as follows. In Section 2 the model formulation is presented, together with some additional definitions and notation. Section 3 is devoted to stability analysis of spatially homogeneous equilibria of the distributed replicator equation. In Section 4 we formulate possible extensions of the notions of the Nash equilibrium and ESS for our spatially explicit model and present some consequences of the new definitions. In Section 5 we derive sufficient conditions for a distributed ESS along with some illustrative examples.

Replicator equation with explicit space
Let Ω be a bounded domain, Ω ∈ R m , m = 1, 2, or 3, with a piecewise-smooth boundary Γ. In the following we assume, without loss of generality, that the volume of Ω is equal to 1, i.e., Ω dx=1. A spatially explicit analogue to (1) is given by the following reaction-diffusion system Here u i = u i (x, t), x ∈ Ω, ∂ t = ∂ ∂t , ∆ is the Laplace operator, in Cartesian coordinates ∆ = m k=1 A is a given constant matrix, and d i > 0, i = 1, . . . , n are the diffusion coefficients.
The initial conditions are u i (x, 0) = ϕ i (x), i = 1, . . . , n, and the form of f sp (t) will be determined later.
It is natural to assume that we consider closed systems (see also [29]), i.e., we have the boundary conditions where n is the normal vector to the boundary Γ.
As it was discussed in Section 1, the global regulation of the total species concentrations occurs in the system, such that for any time moment t. This condition is analogous to the condition for the constant concentration in the finite-dimensional case [17]. From the boundary condition (7) and the integral invariant (8) the expression for the function f sp (t) follows: Suppose that for any fixed x ∈ Ω each function u i (x, t) is differentiable with respect to variable t, and belongs to the Sobolev space W 1 2 (Ω) if m = 1 or W 2 2 (Ω) if m = 2, 3 as the function of x for any fixed t > 0. Here W k 2 , k = 1, 2 is the space of functions, which have square integrable derivatives with respect to x ∈ Ω up to the order k. Note, that from the embedding theorem (e.g., [7,24]) we have that any function from the space W k 2 (Ω) is continuous, except possibly on the set of measure zero, and this continuity is used in some proofs below.
Denote Ω t = Ω × [0, ∞) and consider the space of functions B(Ω t ) with the norm Hereinafter we denote S n (Ω t ) the set of non-negative vector-functions u(x, t), u i (x, t) ∈ B(Ω t ), i = 1, . . . , n, which satisfy (8), and we use the notation int S n (Ω t ) for the set of functions We consider weak solutions to (6), i.e., the solutions should satisfy the integral identity for any function η = η(x, t) on compact support, which is differentiable on [0, ∞) with respect to t and belongs to the Sobolev space W 1 2 (Ω) for any fixed t > 0. Remark that generally system (6) is not a "system of differential equations" because its right-hand side contains functional (9).
The steady state solutions to (6) can be found as the solutions to the following elliptic problem with the boundary conditions The set of all non-negative vector-functions (11) is denoted S n (Ω). Using (11) and (10) we obtain that i.e., f sp is a constant. The rest points of (1) are spatially homogeneous solutions to (10); the converse also holds: any spatially homogeneous solution to (10) is a rest point of (1). In [5,6] we proved that for sufficiently small values of the diffusion coefficients d i there exist nonhomogeneous solutions to (10) in the case (Au) i = k i u i and (Au) i = k i u i−1 for arbitrary positive constants k i (these are autocatalytic and hypercyclic systems, respectively, for more details see [5,6]). Here we undertake a task to investigate the stability of spatially homogeneous solutions to (6), which can differ from the stability of the rest points of the local system (1) due to the explicit space in the model. We also consider the conditions necessary for spatially non-homogeneous solutions to appear. And yet another purpose is to transfer the definitions of the Nash equilibria and ESS for the distributed system (6)- (8) and apply these concepts to study the asymptotic behavior of (6).

Stability of spatially homogeneous equilibria
We start with the standard definition that is used extensively throughout the text.
ofŵ(x) that for any initial data of (6) from U δ it follows that for any t ≥ 0. If in (13) the left hand side tends to zero, thenŵ(x) is asymptotically stable.
Consider the following eigenvalue problem The eigenfunction system of (14) is given by and forms a complete system in the Sobolev space W 2 2 (Ω) (e.g., [24]), additionally where δ ij is the Kronecker symbol. The corresponding eigenvalues satisfy the condition The following theorem gives a necessary condition for a spatially homogeneous solution to (6) be asymptotically stable, if the corresponding rest point of (1) is asymptotically stable. Theorem 1. Letv ∈ int S n be an asymptotically stable rest point of (1). Then for this point to be an asymptotically stable homogeneous stationary solution to (6) it is necessary that where λ 1 is the first non-zero eigenvalue of (14).
Proof. Remark that ifv is Lyapunov stable thenv is a Nash equilibrium of the game with the payoff matrix A. Consider the solutions to (6)- (8) assuming that the Cauchy data are perturbed: which is always possible since the eigenfunctions ψ i (x) of (14) form a complete system in W k 2 (Ω). c i k (t), k = 1, 2, 3, . . . , i = 1, . . . , n are smooth functions of t. Note that w 0 i (x) = w i (x, 0), i = 1, . . . , n. The spatially homogeneous equilibriumv is stable if c i k (t) → 0, t → ∞ for all i and k.
Due to the fact that we have from (14) that Substituting (17) into (6) and retaining in the usual way only linear terms with respect to w i we obtain the following equations: for i = 1, . . . , n. Here c 0 = (c 1 0 (t), . . . , c n 0 (t)), and τ denotes transposition. From (2) it follows that the last term in (19) is zero. From (4) and (18) we also have that Therefore, we obtain that System (20) is linear, with the matrix Q = q ij i,j=1,...,n with the elements which coincides with the Jacobi matrix of (1) at the rest pointv. Due to the assumptions the rest point of (1) is asymptotically stable, therefore the trivial stationary point of (20) is also asymptotically stable. Multiplying equations (6) consequently by ψ k (x), k = 1, 2, . . . and retaining only linear terms, after substituting (17) into (6), we obtain where c k = (c 1 k (t), . . . , c n k (t)), λ k are the eigenvalues of (14). Linear system (22) has the matrix R k with the elements r k ij = a ijvi − λ k d i δ ij , where δ ij is the Kronecker symbol. From (21) and the last expressions we have that System (20) is stable, hence tr Q = n i=1 q ii < 0. On the other hand therefore tr R k is negative if (16) holds, which completes the proof. Denoting System (22) has the matrix For the equilibria to be stable we need first that tr R k < 0, which is equivalent to Condition (23) is a generalization of (16) since c 1 c 2 /(c 1 + c 2 ) < 0. The second condition is det R k > 0, or (24) is satisfied. The last inequality is equivalent to In the case ∆ < 0 the condition for the replicator equation (6) to have asymptotically stable equilibrium is that the largest root of the quadratic equation Therefore the necessary and sufficient conditions for the homogeneous rest pointû ∈ int S 2 (Ω) be asymptotically stable are the conditions (23) and (25) or the conditions (23) and (26). Consider for example the simplest hypercycle equation with A = 0 a 12 a 21 0 .
Here we have that condition (23) is always satisfied. The condition det R k > 0 can be written as d 1 d 2 λ 2 k > a 12 a 21v1v2 . If this condition does not hold then, as was shown in [6] for the general n-dimensional case, the inner rest point becomes unstable.
If Ω is one-dimensional, then for the cases (Au) i = k i u i and (Au) i = k i u i−1 there exist nonuniform steady state solutions to (10)-(12) (see [4,6]). More precisely, a necessary condition for such a solution to exist is We note that π 2 is the first non-zero eigenvalue of (14) when Ω = [0, 1]. Remark that this condition is a particular case of the general condition (16), and we conjecture that the same situation occurs in the general case. The natural question is whether non-uniform stationary solutions appear in the general case (6). Here we show that, at least for some matrices A, system (6) possesses non-uniform steady state solutions.
We rewrite (10) in the form for any i = 1, . . . , n. The initial data are We assume that matrix A is such that Aw > 0 for any w > 0. System (27) is conservative, its rest points (ŵ,p) can be found from p i = 0, (Aŵ) i = β, β = Aŵ,ŵ , i = 1, . . . , n, thereforeŵ is a rest point of (1). Consider the Jacobi matrix of (27) evaluated at the rest point (0,ŵ): where 0 is the n × n zero matrix, I is the n-dimensional identity matrix, and J l d is the Jacobi matrix of (1) at the rest pointv =ŵ when i-th row divided by d i for each row. The eigenvalues of J are the roots of the equation where det J l is the determinant of the Jacobi matrix of (1) atv. if det J l < 0 then for any n ≥ 2 the last equation has pure imaginary roots. If det J l > 0 then the same is possible when n = 2k + 1, k = 1, 2, . . . Let us introduce the following sets in the phase space: From (28) it follows that at the initial "time" the orbit of the system (27) belongs to Σ. Suppose that w(0) ∈ U − . Then the second equation in (27) yields that functions p i (x) increase as x increases, and therefore p i (x) > 0, x > 0, i = 1, . . . , n. From the condition on A it follows that (Aw(0)) i < (Aw(x)) i , x > 0 and therefore the values (Aw(x)) i − f sp decrease as x increases, hence w(x) has to cross the hyperplane Π where (Aw(x)) i = f sp , i = 1, . . . , n. After that, the solution gets into the set U + , which implies that p i (x), i = 1, . . . , n are decreasing and w(x) are increasing. Therefore there exists such value x * i such that p i (x * i ) = w ′ i (x * i ) = 0, i = 1, . . . , n (see Fig. 1).
The diffusion coefficients d i characterize (n + i)−th component of the movement speed along the phase trajectories. An increase (decrease) in values d i corresponds to the decrease (increase, respectively) in the movement speed of (n + i)−th component of the speed vector. Therefore, it is possible to find such values of d i , i = 1, . . . , n such that all x * i = 1. We remark that if we reduce values d i twice this would mean that the phase orbit again would reach Π (one cycle). Reducing d i four times, we obtain the orbit that makes two cycles in the state space, and so on. Therefore, from the discussion above, it follows that system (10) can have non-uniform stationary solutions which may possess an arbitrary number of oscillations (see a similar discussion in [6], where some examples are given).

Dynamics of the distributed replicator equation
Definition 2. We shall say that the vector functionŵ( for any vector-function u(x, t) ∈ S n (Ω t ), u(x, t) =ŵ(x).

Remark 2.
Letŵ ∈ S n satisfy (29). Thenŵ is a Nash equilibrium in the game with payoff A. Indeed, (8) implies thatū(t) ∈ S n for any t. Therefore we have ū, Aŵ ≤ Aŵ,ŵ , which means thatŵ is a Nash equilibrium and therefore a rest point of (1).
Proof. Due to the fact thatŵ(x) is Lyapunov stable, then, for any initial data from a neighborhood U δ ofŵ(x) in the space W k 2 (Ω), the corresponding solution satisfies (13). Suppose thatŵ(x) is not a distributed Nash equilibrium. Then, using continuity of the scalar product, there exists an index i and constant ξ > 0 such that for all u(x, t) ∈ int S n (Ω t ) in a neighborhood ofŵ(x) ∈ int S n (Ω).
From Definitions 1 and 4 it follows that stability in the mean integral sense is weaker than Lyapunov stability. For instance, consider a simple example: let g(x, t) ∈ W 1 2 , x ∈ [0, 1] and can be represented as Let us suppose that c 0 (t) → 0 when t → ∞.
does not necessarily tend to zero.
Theorem 3. Letŵ ∈ int S n be a spatially homogeneous solution to (10) (i.e.,ŵ ∈ int S n is a rest point of (1)). Ifŵ is DESS thenŵ is an asymptotically stable solution to (6) in the sense of the mean integral value.
Proof. Consider the set of function u(x, t) ∈ S n (Ω t ) belonging to a neighborhood U δ ofŵ in the space B(Ω t ). Define the functional We can always choose U δ such that u i (x, t) > 0 for all i becauseŵ i > 0. Using (32) we obtain where we suppress the dependence on t and x for simplicity. From (34) we have that The functional V is bounded. Indeed, applying Jensen's inequality for sums Hence, SinceV (u) > 0 for all u ∈ U δ then V is a strict Laypunov functional for (6)   Using (34) we obtain (36).

Sufficient conditions for DESS
Application of the results obtained so far for particular distributed replicator systems reduces to the problem of checking conditions for DESS for spatially homogeneous stationary solution w ∈ int S n Let us introduce the following function If the stationary pointŵ ∈ int S n meets the conditions for DESS then M (t) < 0 in a neighborhood ofŵ in the space W k 2 (Ω) for functions u(x, t) ∈ S n (Ω t ) from neighborhood ofŵ(x) in the space S n (Ω t ), such that u(x, t) =ŵ(x), x ∈ Ω, t ≥ 0.
Example 1. Consider system (6) with the matrix We have Ac 0 (t), c 0 (t) = −(c 0 1 ) 2 (b − a + c − d), and hence condition (42) is satisfied if In Fig. 2 (16) is not satisfied and there are spatially heterogeneous stationary solutions to the system (6). The dynamics of the solutions is shown in Fig. 2(a), (b). Spatially non-uniform solutions are stable in this case. In Fig. 3 all possible stable spatially non-uniform solutions are presented, each of which has its own basin of attraction.
Example 3 (Hypercycle equation). Consider a hypercycle system with three members and the matrix which yields (42). x w i (x) (d) Figure 3: Asymptotically stable stationary non-uniform solutions to the problem as in Example 1. See text for details. The numerical scheme is presented in [5] Matrix A is circulant, and its eigenvalues are Assume a < (b + c)/2, b > a > c > 0, bc > a 2 , then Re λ 2,3 < 0. λ 1 has the eigenvector (1, 1, 1). This vector is orthogonal to any c 0 that satisfies (39), hence A is negatively determined on this set.
Let a 1 + a 3 > a 2 . The eigenvalues of A are given by In the same vein as in Example 4, λ 1 corresponds to the eigenvector u = (1, 1, 1, 1), which is orthogonal to any c 0 satisfying (39). The other eigenvalues have negative real parts, therefore A is negatively determined. Consider the following matrix  Fig. 4). In Fig. 5 possible asymptotically stable non-uniform solutions are shown. The details of the numerical scheme are given in [5].
Here the conditions for the spatially homogeneous solutions to be DESS are fulfilled, therefore, this solution is stable in the sense of the mean integral value, which can be seen from the pictures.