Gradient flow formulations of discrete and continuous evolutionary models: a unifying perspective

TL;DR

This paper unifies classical biological evolution models by reformulating the Moran process, Kimura Equation, and Replicator Equation as gradient flows, revealing their geometric relationships and convergence properties.

Contribution

It introduces a gradient flow framework for these models, connecting discrete and continuous evolutionary dynamics through geometric and variational principles.

Findings

Reformulation of Moran process and Kimura Equation as gradient flows.

Conditions for convergence of gradient structures from Moran to Kimura, and Kimura to Replicator.

Geometric characterization of evolutionary processes as energy minimization problems.

Abstract

We consider three classical models of biological evolution: (i) the Moran process, an example of a reducible Markov Chain; (ii) the Kimura Equation, a particular case of a degenerated Fokker-Planck Diffusion; (iii) the Replicator Equation, a paradigm in Evolutionary Game Theory. While these approaches are not completely equivalent, they are intimately connected, since (ii) is the diffusion approximation of (i), and (iii) is obtained from (ii) in an appropriate limit. It is well known that the Replicator Dynamics for two strategies is a gradient flow with respect to the celebrated Shahshahani distance. We reformulate the Moran process and the Kimura Equation as gradient flows and in the sequel we discuss conditions such that the associated gradient structures converge: (i) to (ii) and (ii) to (iii). This provides a geometric characterisation of these evolutionary processes and provides a reformulation of the above examples as time minimization of free energy functionals.

Equations

$$M_{ij}=\left\{\begin{array}[]{ll}\frac{N-j}{N}\mathsf{s}_{j},&\ i=j+1\ ,\\ \frac{j}{N}\mathsf{s}_{j}+\frac{(N-j)}{N}(1-\mathsf{s}_{j}),&\ i=j,\\ \frac{j}{N}(1-\mathsf{s}_{j}),&\ i=j-1,\\ 0,&\text{otherwise.}\end{array}\right.$$

$$\mathbf{p}(t+\Delta t)=\mathbf{Mp}(t),$$

$$\frac{\mathrm{d}\mathbf{p}}{\mathrm{d}t}=\left(\mathbf{M}-\mathbf{I}\right)\mathbf{p},$$

$${\mathsf{s}}_i=\frac{i}{N}\left[1-\frac{2}{\kappa N}\left(1-\frac{i}{N}\right)V^{\prime }\left(\frac{i}{N}\right)\right]$$

$${\partial }_{t}p=\frac{\kappa }{2}{\partial }_{xx}^2\left(x(1-x)p\right)+{\partial }_x\left(x(1-x)V^{\prime }(x)p\right)$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^{1}p(t,x)\mathrm{d}x=0\text{mbox}and\frac{\mathrm{d}}{\mathrm{d}t}{\int }_0^1\phi (x)p(t,x)\mathrm{d}x=0,$$

$$\dot{X}=-X(1-X)V^{\prime }(X).$$

$${\mathbf{R}}^{†}\mathbf{MR}=\left(\begin{array}{cc}\tilde{\mathbf{M}}& 0\\ \mathbf{A}& \mathbf{I}\end{array}\right)$$

$$\mathbf{M}=\left(\begin{array}{cc}\tilde{\mathbf{M}}& 0\\ \mathbf{A}& \mathbf{I}\end{array}\right)$$

$$\tilde{\mathbf{p}}=(p_i{)}_{i=0,\dots ,N-k},\frac{\mathrm{d}\tilde{\mathbf{p}}}{\mathrm{d}t}=(\tilde{\mathbf{M}}-\mathbf{I})\tilde{\mathbf{p}}.$$

$$\mathbf{M}=\left(\begin{array}{ccc}1& {\mathbf{\alpha }}^{†}& 0\\ 0& \tilde{\mathbf{M}}& 0\\ 0& {\mathbf{\beta }}^{†}& 1\end{array}\right)$$

$$⟨\mathbf{p}(t),{\mathbf{F}}_l⟩=⟨\mathbf{p}(0),{\mathbf{F}}_l⟩,\forall t>0\text{ and }l=1,\dots ,k.$$

$$\exp \left(t(\mathbf{M}-\mathbf{I})\right)=\left(\begin{array}{cc}\exp \left(t(\tilde{\mathbf{M}}-\mathbf{I})\right)& 0\\ \mathbf{A}(\tilde{\mathbf{M}}-\mathbf{I}{)}^{-1}\left(\exp \left(t(\tilde{\mathbf{M}}-\mathbf{I})\right)-\mathbf{I}\right)& \mathbf{I}\end{array}\right),$$

$$\mu :=\rho (\tilde{\mathbf{M}})\in (0,1).$$

$${\tilde{w}}_i{\tilde{M}}_{ij}{\tilde{z}}_j={\tilde{w}}_j{\tilde{M}}_{ji}{\tilde{z}}_i,\forall i,j=0,\dots ,N-k.$$

$$\mathrm{Diag}[\tilde{\mathbf{w}}]\tilde{\mathbf{M}}\mathrm{Diag}[\tilde{\mathbf{z}}]=\mathrm{Diag}[\tilde{\mathbf{w}}\circ {\mathbf{d}}^{-1}]\mathbf{T}\mathrm{Diag}[\mathbf{d}\circ \tilde{\mathbf{z}}],$$

$$\mathbf{K}:=\frac{1}{\mu }\mathrm{Diag}[\tilde{\mathbf{w}}{]}^{-1}{\tilde{\mathbf{M}}}^{†}\mathrm{Diag}[\tilde{\mathbf{w}}].$$

$$\mathbf{\pi }:=\tilde{\mathbf{w}}\circ \tilde{\mathbf{z}}.$$

$${\mathbf{\pi }}^{†}\mathbf{K}$$

$$={\mu }^{-1}{\tilde{\mathbf{z}}}^{†}{\tilde{\mathbf{M}}}^{†}\mathrm{Diag}[\tilde{\mathbf{w}}]={\tilde{\mathbf{z}}}^{†}\mathrm{Diag}[\tilde{\mathbf{w}}]={\mathbf{\pi }}^{†}.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}⟨\tilde{\mathbf{w}},\tilde{\mathbf{p}}(t)⟩=⟨\tilde{\mathbf{w}},(\tilde{\mathbf{M}}-\tilde{\mathbf{I}})\tilde{\mathbf{p}}(t)⟩=(\mu -1)⟨\tilde{\mathbf{w}},\tilde{\mathbf{p}}(t)⟩$$

$$⟨\tilde{\mathbf{w}},\tilde{\mathbf{p}}(t)⟩=e^{(\mu -1)t}⟨\tilde{\mathbf{w}},\tilde{\mathbf{p}}(0)⟩.$$

$$\mathbf{q}(t):=\frac{\tilde{\mathbf{w}}\circ \tilde{\mathbf{p}}(\text{nicefrac}t\mu )}{⟨\tilde{\mathbf{w}},\tilde{\mathbf{p}}(\text{nicefrac}t\mu )⟩}\text{and}\mathbf{Q}(t):=\frac{\mathrm{d}\mathbf{q}}{\mathrm{d}\mathbf{\pi }}(t)={\left(\frac{q_i(t)}{{\pi }_i}\right)}_i.$$

$$H(\mathbf{q}|\mathbf{\pi }):=\sum _{i=0}^{N-k}\frac{q_i}{{\pi }_i}\log \left(\frac{q_i}{{\pi }_i}\right){\pi }_i=\sum _{i=0}^{N-k}{q}_i\log \left(\frac{q_i}{{\pi }_i}\right).$$

$$\beta(x,y):=\left\{\begin{array}[]{ll}\frac{x-y}{\log x-\log y}&\mbox{if }x\neq y\\ x&\mbox{otherwise}\end{array}\right..$$

$${\mathcal{W}}_N^2({\mathbf{Q}}^0,{\mathbf{Q}}^1):=\underset{\mathbf{Q},\mathbf{\psi }}{\inf }\left\{\frac{1}{2}{\int }_0^1\sum _{i,j}|{\psi }_i^{\tau }-{\psi }_j^{\tau }{|}^2{K}_{ij}\beta (Q_i^{\tau },Q_j^{\tau }){\pi }_i\mathrm{d}\tau \right\},$$

$$\left\{\begin{array}[]{l}\frac{\mathrm{d}}{\mathrm{d}\tau}Q^{\tau}_{i}+\sum\limits_{j}(\psi^{\tau}_{j}-\psi^{\tau}_{i})K_{ij}\,\beta(Q^{\tau}_{i},Q^{\tau}_{j})=0\quad\,\text{a.e. }\tau\in[0,1],\,\forall i\ ,\\ \mathbf{Q}^{\tau}|_{\tau=0}=\mathbf{Q}^{0},\qquad\mathbf{Q}^{\tau}|_{\tau=1}=\mathbf{Q}^{1}.\end{array}\right.$$

$${\left[{\mathrm{grad}}_{{\mathcal{W}}_N}F(\mathbf{Q})\right]}_i=-\sum _j{K}_{ij}\beta (Q_i,Q_j)\left(f^{\prime }(Q_j)-f^{\prime }(Q_i)\right).$$

$$\frac{\mathrm{d}\mathbf{q}}{\mathrm{d}t}=-{\mathrm{grad}}_{{\mathcal{W}}_N}H(\mathbf{q}|\mathbf{\pi })$$

$$\frac{\mathrm{d}}{\mathrm{d}\tau }F({\mathbf{Q}}^{\tau })={⟨{\mathrm{grad}}_{{\mathcal{W}}_N}F({\mathbf{Q}}^{\tau }),\frac{\mathrm{d}{\mathbf{Q}}^{\tau }}{\mathrm{d}\tau }⟩}_{{\mathbf{Q}}^{\tau }}.$$

$${\mathrm{grad}}_{{\mathcal{W}}_N}F(\mathbf{Q})=-{\mathrm{div}}_{\mathbf{Q}}(\nabla f^{\prime }(\mathbf{Q})).$$

Full text

Keywords: Gradient Flow Structure; Optimal Transport; Replicator Dynamics; Shahshahani Distance; Reducible Markov Chains; Kimura Equation.

1 Introduction

1.1 Background

From a contemporary perspective, evolution can be conveniently described as being the product of changes in allele frequencies within a population — cf. [43]. Albeit with an apparently simple definition, evolution is actually a complex phenomenon and, as such, it comprises many different mechanisms: natural selection, mutation, genetic-drift, to name only a few.

The need to understand these different mechanisms and, more recently, their interplay, led to the development of a plethora of models in evolutionary dynamics: discrete time Markov chains were used in the early 20th century to study genetic drift (e.g., the Wright-Fisher [110, 48] and the Moran [89] processes); continuous time stochastic processes in the mid 20 century geared towards molecular evolution, as is the case of Kimura Partial Differential Equation [70]; and, finally, systems of Ordinary Differential Equations (ODE) that are used to model natural selection, cf. [105]. These three classes of basic models can be considered as a classical triad [19, 21], and they will be the focus of this work. It should be also noticed that, more recently, new modelling paradigms have been used — most notably Individual Based Models [59, 34] and kinetic models [10, 106].

The study of different connections between models in this classical triad dates back at least to [41], where a frequency-dependent version of the Wright-Fisher process was introduced, and the large population regime was shown to be described by a generalised version of the Kimura Equation. More recently a number of works have explored these links providing various approaches to a unified view on these models, in the weak selection regime with infinite population limit and suitable scaling relations between the time step and population size [24, 19, 21]. In addition, the Kimura Equation and the Replicator Dynamics (RD) are connected for short times and strong selection. However, despite all these connections, there are also important differences — see [22] for results on the qualitative difference of fixation probability in large populations compared to infinite ones, and [23] for notable features of the fixation probability in finite populations that are not in the weak selection regime. These results suggest that the impact of all underlying assumptions in each of these three models is not yet fully understood.

The aim of this work is to investigate this classical triad from yet a different perspective which, as far as we are aware of, seems to have been overlooked: namely, the fact that these models can be formulated (or at least made compatible with) some sort of local maximisation principle. This approach has a long tradition in the biological literature, cf. [69, 42, 9]; for the relation between optimal principles in evolution and the Fundamental Theorem of Natural Selection, due to Fisher [49], see also [44, 45, 46].

We should point out that we are not attempting to obtain a global maximisation principle. The existence (or usefulness) of global principles is a quite controversial topic in evolution, and we refer the reader to [102] for a review on optimising techniques in evolution and to [85, 84] for a critique on this approach.

1.2 Gradient flow formulations of evolutionary models and main results

It will turn out that the fundamental tool that will allow us to accomplish the previous task is the concept of gradient flows. This is a rather classical topic in differential equations that has raised much interest recently after the water-shedding work [6] — see also [95] for a very gentle introduction.

Gradient flows are hardly new in evolutionary dynamics: under some hypotheses, the RD can be reinterpreted as a gradient flow with respect to a specific metric — known by now as the Shahshahani metric. In particular, for the one dimensional case, the RD is always a gradient flow in this metric [98, 1, 2].

Motivated by the positive results of a research program undertaken by two of the authors (FACCC and MOS) in studying this triad starting from the discrete processes [19, 21, 22], we will follow the same pattern with a slight modification: we will start from the continuous-time generalisation of the well-known Moran process [40]. This will allow us to adapt the framework recently developed by [79] to our case, and obtain a formulation of the transient part of the Moran process as a gradient flow. These adaptations turn out to be deeply connected with the so-called associated $Q$-process, which describes the probability law of eternal paths in this absorbing system. This formulation will also provide a “geometrisation” for finite populations, and answers a question raised in [1].

Finally, from well-known gradient flow formulations of Fokker-Planck type equations[66], and based on the $Q$-process in the continuous setting, we obtain a gradient flow formulation for the Kimura Partial Differential Equation.

Summing up, we were able to derive independent gradient flow formulations of the triad. Subsequently, we study the natural compatibility between these models, which is expected to hold due to some of the authors’ previous results.

More precisely, we will:

Reformulate all three models as gradient flows, i.e, in each case we find a free-energy $\mathsf{H}$ defined on a Riemannian manifold with distance $\mathsf{w}$ such that the model is equivalent, in a sense to be made precise, to the steepest descent flow of $\mathsf{H}$ as measured by $\mathsf{w}$. See [6, 94, 107, 5] for an overview on the topic. 2. 2.

Obtain these processes as a time-step minimisation of $\mathsf{H}$. Namely, consider that at time $t$ the system is at state $p$, where $p$ is a probability density distribution describing all possible states of the system. In the next time step, the system will be at state $p^{\prime}$ such that the infimum of $q\mapsto\frac{h}{2}\mathsf{w}^{2}(p,q)+\mathsf{H}(q)$ is attained, for a small and positive $h$. After the seminal paper by Jordan, Kinderlehrer and Otto [66], this second approach became known as JKO scheme.

Under appropriate convexity assumptions these two approaches are equivalent in a very general setting [6], but the construction can also be made rigorous in the absence of convexity on a case-by-case basis; see, e.g., [74, 36, 12, 76, 72]. The three models are linked by two limiting processes: the partial differential equation (PDE) is connected to the Markov Chain by an infinite population limit, and to the ODE by a vanishing viscosity limit. One may therefore wonder whether (or hope that) all these gradient flow formulations are compatible in some sense. An appropriate tool to investigate this turns out to be the $\Gamma$-convergence of gradient flows [97, 93] and we will discuss how this connection can be obtained.

A more thorough discussion of modelling implications to evolution will be postponed to a subsequent work. However, we should state that both the free energy and metric will be derived from a set of common assumptions in evolutionary dynamics and are not artificial quantities. The metric is the Wasserstein distance between two probability measures built upon a generalisation of the so-called Shahshahani metric [98], introduced in the framework of gradient systems in the Replicator Equation; cf. [58, 99, 100, 60], see also the discussion on Kimura’s Maximum Principle and the (Svirezhev-) Shahshahani metric in [13]. We note also that, for finite populations, short-term and long-term information on the dynamics has been recently obtained from the free energy [18].

1.3 Outline

The structure of this paper is as follows. In the remainder of this section, we first introduce and fix the basic notation that will be used throughout this work and we recall the theory previously established related to the connections between the models. In the sequel, we review the main results by FACCC, MOS, and collaborators related to the current work.

In Section 2, we introduce a class of matrices that encompass celebrated classes of matrices used in population genetics, and introduce a class of continuous-in-time discrete Markov chains. In order to reformulate this class of processes as gradient flows, we appeal to the so-called $Q$-processes and generalise recent results by [79] for irreducible Markov chains to the case with two or more absorbing states.

In Section 3, we digress and discuss the relation between the discrete and continuous time Markov chains. While the discrete ones are not amenable to the treatment developed in Section 2, we opt to include it here as it is the most traditional framework for finite-population models in population genetics. Furthermore, it has been the starting point of previous work from part of the authors, as described above. In particular, we show how the results from the previous section can be used to define entropies in population dynamics. The framework developed herein is sufficiently general to treat simultaneously the Boltzmann-Gibbs-Shannon (BGS) entropy and the Tsallis entropies (both in discrete and continuous time).

In Section 4, we describe the gradient flow formulation of the Kimura Equation, a degenerated PDE of drift-diffusion type. As already shown in [20, 21] the appropriate formulation satisfies two conservation laws and has measure-valued solutions. However, as in the discrete case, it turns out that only the interior dynamic can be formulated as a gradient flow and we focus on this part of the solutions.

In Section 5 we study the Replicator Equation, which is well known [98, 1] to have a variational structure. In fact, we shall rather study its fully equivalent PDE version. The latter is more appropriate for our framework, but all our results immediately apply to the former.

In Section 6, we discuss the compatibility between the variational structure (both gradient flow and JKO formalism) for all the models discussed in Sections 2, 4, and 5. As mentioned earlier the $\Gamma$-convergence of gradient flows will be the important tool in this section — see [30] for a general introduction on $\Gamma$-convergence and [97, 93] for $\Gamma$-convergence of gradient flows — see also [51].

We finish with some comments in Section 7.

A roadmap of the paper summarising the main relationships between the three models is given in Figure 1.

1.4 Notations

Let us consider a $N+1$ Markov chain on an abstract finite space $\mathcal{X}_{N}$. By abuse of notation, we write $\mathcal{X}_{N}\mathrel{\mathop{:=}}\left\{0,\dots,N\right\}$. We will use bold symbols for discrete vectors and matrices. Vectors are considered as column vectors by default. Given a matrix, and unless otherwise specified, the words stochastic and substochastic mean column-stochastic and column-substochastic, respectively.

By $\mathcal{P}(\mathcal{X}_{N})$ we denote the space of probability measures in $\mathcal{X}_{N}$, canonically identified with vectors in the $N$-dimensional simplex $\Delta^{N}\mathrel{\mathop{:=}}\{\mathbf{p}\in\mathbb{R}_{+}^{N+1}|\sum_{i}p_{i}=1\}$. We shall use lowercase to denote probability vectors, and uppercase for their densities with respect to some reference measure: Typically $\mathbf{p}=\left(p_{i}\right)_{i\in\mathcal{X}_{N}}$ will denote an arbitrary probability vector; if $\bm{\pi}=\left(\pi_{i}\right)_{i\in\mathcal{X}_{N}}$ is a particular reference probability measure, we write $\mathbf{P}\mathrel{\mathop{:=}}\frac{\mathrm{d}\mathbf{p}}{\mathrm{d}\bm{\pi}}=\left(\frac{p_{i}}{\pi_{i}}\right)_{i\in\mathcal{X}_{N}}$ for the density of the measure $\mathbf{p}$ with respect to the measure $\bm{\pi}$.

For Markov chains with more than one absorbing states, transient parts are particularly relevant. Therefore, we denote them by tilded quantities $\widetilde{\mathbf{p}}=(p_{i})_{i=0,\dots,N-k}$, where $k$ is the number of linearly independent absorbing states. In a slightly inconsistent notation, and whenever we focus on chains with only two absorbing states, these will sometimes be labelled instead as $i=0$ and $i=N$ for convenience, and tilded quantities will denote projections to the “interior”, e.g. $\tilde{\mathbf{p}}=(p_{i})_{i=1,\dots,N-1}$. In simple words, we discard the $i=0,N$ entries in the vector $\mathbf{p}=(p_{i})_{i=0,\dots,N}$.

The transition matrix of the Markov chain on the state space $\mathcal{X}_{N}$ is given by a $(N+1)\times(N+1)$ matrix $\mathbf{M}=\left(M_{ij}\right)_{i,j\in\mathcal{X}_{N}}$. The matrix $\widetilde{\mathbf{M}}\mathrel{\mathop{:=}}\left(M_{ij}\right)_{i,j=0,\dots,N-k}$ associated to the transient part of the process is called the core matrix associated to $\mathbf{M}$. (Depending on the particular labelling, we also consider $\widetilde{\mathbf{M}}=\left(M_{ij}\right)_{i,j=1,\dots,N-1}$.)

If $\mathbf{A},\mathbf{B}\in\mathbb{R}^{n\times m}$, then we will write $\mathbf{A}\circ\mathbf{B}\in\mathbb{R}^{n\times m}$ to denote their entry-wise (Hadamard) product — i.e. $\left(\mathbf{A}\circ\mathbf{B}\right)_{ij}=A_{ij}B_{ij}$. Note that $\left(\mathbf{A}\circ\mathbf{B}\right)^{\dagger}=\mathbf{A}^{\dagger}\circ\mathbf{B}^{\dagger}$. For $\mathbf{a}\in\mathbb{R}^{n}$, we will write $\mathop{\mathrm{Diag}}[\mathbf{a}]\in\mathbb{R}^{n\times n}$ to denote the $n\times n$ matrix with its main diagonal given by $\mathbf{a}$ and zero elsewhere. Recall that the usual matrix product and the Hadamard product agree for diagonal matrices — in particular, $\mathop{\mathrm{Diag}}[\mathbf{a}\circ\mathbf{b}]=\mathop{\mathrm{Diag}}[\mathbf{a}]\mathop{\mathrm{Diag}}[\mathbf{b}]$, and $\frac{\mathrm{d}\mathbf{p}}{\mathrm{d}\bm{\pi}}=\mathbf{p}\circ\bm{\pi}^{-1}$.

In the continuous setting, $x\in\Omega\mathrel{\mathop{:=}}[0,1]$ denotes the state of the system, i.e., $x$ is the fraction of individuals of the focal type. In particular, the absorbing states are indicated by $x\in\partial\Omega\mathrel{\mathop{:=}}\{0,1\}$. This explains the reordering in the indexes used in the discrete setting when there are exactly 2 absorbing states. We shall write $p\in\mathcal{P}(\Omega)$ for measures on the whole domain, while $\widetilde{p}$ will denote the restriction $p\mathop{\lfloor}{\Omega}$ to the interior, i.e. $p=a\delta_{0}+\widetilde{p}+b\delta_{1}$, for $a,b\geq 0$. We denote by $q\ll p$ and $q\perp p$ if $q$ is absolutely continuous with respect to $p$, and if $q$ is singular with respect to $p$, respectively.

In general, tilded quantities, both in the discrete and continuous settings, are not probability measures and we will often need to rescale and/or renormalise those interior projections so that they become again probability measures: the resulting scaled variables will be denoted by $\mathbf{q}=(q_{i})_{i=1\dots N-1}$ or $q\in\mathcal{P}(\Omega)$, while we keep the letter $\mathbf{p},p$ for the initial variables.

Unless otherwise specified dotted quantities will denote time derivatives $\dot{X}(t)=\frac{\mathrm{d}X}{\mathrm{d}t}(t)$, while primes will stand for spatial derivatives $V^{\prime}(x)=\frac{\mathrm{d}V}{\mathrm{d}x}(x)$.

A time step will be denoted $\Delta t$ if it is a parameter of the model, or $h$ if it corresponds to a discretization of the continuous time variable $t$.

1.5 State of the art

We finish this introduction with a more detailed explanation of the triad, and the links between these three different classes. Note that in part of the current work, however, we opt to describe more general processes than the one exemplified here: in Sections 2 and 3, we consider stochastic processes with an arbitrary number $k>1$ of absorbing states, while in Sections 4 and 5, we consider more general diffusion coefficients than in this subsection.

The Moran process was introduced in [89] as a mathematical simplification of the older Wright-Fisher process [110, 48] and it is a particular example of a birth and death process. In the setting we are interested in, we consider a population of fixed size $N$ divided into two groups of individuals, indistinguishable apart from the characteristic under study. Let us call these two types $\mathbb{A}$ and $\mathbb{B}$. Every $\Delta t$ seconds, one individual is chosen to die with uniform probability $1/N$, while a second one (possibly the same one) is chosen to reproduce according to a certain type selection probability vector $(\mathsf{s}_{0},\mathsf{s}_{1},\dots,\mathsf{s}_{N})$. Here $\mathsf{s}_{i}$ indicates the probability to select for reproduction an individual of type $\mathbb{A}$ in a population with $i$ individuals of type $\mathbb{A}$. In this case, we say that the population is at state $i$. Because the Moran process has no mutations (i.e., $\mathsf{s}_{0}=1-\mathsf{s}_{N}=0$) there are two absorbing states, namely $i=0$ and $i=N$.

At each time step the transition probability from state $j$ to $i$ is thus given by a $(N+1)\times(N+1)$ matrix $\mathbf{M}=\left(M_{ij}\right)_{i,j=0,\dots,N}$ with

[TABLE]

The evolution equation (also known as master equation) is given by

[TABLE]

where $\mathbf{p}(t)=(p_{0}(t),\dots,p_{N}(t))$, and $p_{i}(t)$ indicates the probability to find the system at state $i$ at time $t$.

An alternative view of the Moran process will be presented in Section 2, where time will be considered a continuous variable and therefore, the evolution equation will rather be given by

[TABLE]

being $\mathbf{I}$ the $(N+1)\times(N+1)$ identity matrix. The obvious link between equations (2) and (3) will be fully exploited in Section 3 and will be instrumental to build the link between finite and infinite population evolutionary models, and also in order to translate to the more usual setting (2) all the results found for (3). In fact, as we will see in a moment, infinite population in previous works of some of the authors is derived from the discrete-time evolution (2), while the gradient flow formulation, the main object of the present work, will require from the start a continuous time. By contrast with previous works, the infinite population model will be derived here from the continuous time (3).

In order to study the limit $N\to\infty$, $\Delta t\to 0$ it is necessary to assume a certain scaling relation between the population size and time-step, as well as the so-called weak selection principle: At leading order, the type selection probability must be of the specific form

[TABLE]

for a given potential $V(x)$ (the gradient $-V^{\prime}$ representing the fitness difference between the focal type, $\mathbb{A}$, and its opponent $\mathbb{B}$). The parameter $\kappa>0$ is the inverse of the selection strength, see [18] for a detailed analysis of each parameter in Equation (4). In fact, if $N\left(\Delta t\right)^{1/2}\to\kappa^{-1}$, then the vector $\mathbf{p}$ obtained from Equation (2), given a certain initial condition, converges in an appropriate sense to a measure $p$, where $p$ is the solution of a certain degenerate parabolic partial differential equation of drift-diffusion type known as the Kimura equation

[TABLE]

for $(x,t)\in[0,1]\times\mathbb{R}_{+}$. This equation must be supplemented by two conservation laws

[TABLE]

where $\varphi$ is the unique solution of $\varphi^{\prime\prime}-V^{\prime}(x)\varphi^{\prime}=0$ with boundary conditions $\varphi(0)=0$, $\varphi(1)=1$. The initial condition will be the limit, in the same sense, of the initial conditions of the discrete process. See [19] for the derivation of Equation (5) from the Moran process, and [21] for its generalisation to an arbitrary number of types derived from the Wright-Fisher process, i.e. a process such that the transition matrix probability from state $j$ to $i$ is given by $M_{ij}=\binom{N}{i}\mathsf{s}_{j}^{i}(1-\mathsf{s}_{j})^{N-i}$. Finally, see [20, 31] for the detailed study of the Kimura Equation (5) from the partial differential equation point of view.

As a last remark, we note that when $\kappa\ll 1$, the limit of the Kimura Equation is the transport equation $\partial_{t}p=\partial_{x}\left(x(1-x)V^{\prime}(x)p\right)$, which is a PDE version of the well-known Replicator Dynamics

[TABLE]

2 Continuous in time, discrete Markov chains

In this section, we consider continuous time Markov chains on $\mathcal{X}_{N}$, given by Equation (3) with $\mathbf{M}$ a given stochastic matrix, and initial condition $\mathbf{p}(0)=\mathbf{p}^{\mathrm{I}}$. When $\mathbf{M}$ is irreducible, $p(t)$ converges as $t\to\infty$ to a unique and strictly positive invariant probability measure [67]. Under the additional assumption that $\mathbf{M}$ is reversible (see Subsection 2.2), one can construct as in [79] a discrete Wasserstein distance $\mathcal{W}_{N}$ on the space of probabilities $\mathcal{P}(\mathcal{X}_{N})$ such that (3) is the gradient flow of the relative Boltzmann-Shannon-Gibbs entropy, to be defined at Subsection 4.4.

The same discrete distance was constructed independently in [86, 26, 79] and will be discussed later on in Subsection 2.4. Therefore, we refrain from giving the details and precise definitions at this early stage, but it is worth pointing out that this theory of discrete optimal transport crucially requires irreducibility and reversibility of the Markov kernel $\mathbf{M}$. The non-irreducibility of a Markov process is typically due to the existence of absorbing states.

Our goal here is to obtain the aforementioned variational framework for a class of reducible Markov chains. In particular we aim at providing a gradient flow structure for some models of population dynamics that are not covered by a straightforward application of the results in [79], yet include the Moran process.

With that goal in mind, we will first rephrase the aforementioned results to substochastic, irreducible and reversible chains, and subsequently apply the results to our chains of interest, introduced in Subsection 2.1.

Roughly speaking we shall focus on Markov processes for which a particular subdynamics can be identified and allows to reconstruct the whole dynamics, and such that the subdynamics can be recast into an irreducible, reversible Markov process. As alluded to in the introduction, this subdynamics is the core dynamics and corresponds to the evolution of the transient states only. We first make these structural assumptions on the Markov kernels $\mathbf{M}$ more precise and technically explicit in terms of linear algebra (Subsections 2.1 and 2.2). We discuss next the relation between those technical assumptions and probabilistic conditioning of the original processes (Subsection 2.3), we discuss the resulting variational framework (Subsection 2.4), and finally we apply this framework to a time continuous version of one of the models of triad (Subsection 2.5).

2.1 Admissible matrices

In the sequel we will consider Markov processes with $N+1$ states and $1\leq k<N-1$ absorbing states and such that the chain conditioned on non-absorption is irreducible. Moreover, we assume that all absorbing states are accessible. More explicitly,

Definition 1** (Admissible chains).**

Let $1\leq k<N-1$. We say that a $(N+1)\times(N+1)$ stochastic matrix $\mathbf{M}$ is $(N+1,k)$-admissible (admissible, in short), denoted $\mathbf{M}\in\mathcal{A}_{N+1,k}$, if there exists a permutation matrix $\mathbf{R}$ such that

[TABLE]

where $\mathbf{I}$ is the $k\times k$ identity, $\widetilde{\mathbf{M}}$ is an $(N+1-k)\times(N+1-k)$ irreducible matrix, and no row of $\mathbf{A}$ is identically zero. We will refer to $\widetilde{\mathbf{M}}$ as the core matrix associated to $\mathbf{M}$.

One should think here of $k$ absorbing states that do not interact with each other and ${(N+1)-k}$ transient states that possibly self-interact or lose information by getting absorbed. These features are encoded by the matrices $\mathbf{I}$, $\widetilde{\mathbf{M}}$, and $\mathbf{A}$, respectively.

Observe that our structural normalisation (7) is not unique, since one can always further permute any of the first $(N+1-k)$ columns and rows (corresponding to relabelling the transient states) while keeping a similar structure. By abuse of notations we will still talk of the core matrix $\widetilde{\mathbf{M}}$, which is thus defined only up to permutations. As a consequence we always think of the permutation matrix $\mathbf{R}$ as the identity matrix, and of the Markov kernel as already in the canonical form

[TABLE]

In what follows, we shall denote the dynamics of the $(N+1-k)$ transient states as

[TABLE]

By definition 1 of admissible chains the absorption matrix $\mathbf{A}$ has non-zero rows and $\widetilde{\mathbf{M}}$ is strictly substochastic: $\sum_{i}\widetilde{p}_{i}(t)$ is therefore non-increasing in time, and the transient states leak information to the absorbed states.

Remark 1** (Kimura matrices).**

The class of $(N+1,k)$-admissible matrices is an extension of a number of classes previously investigated. In particular, $\mathcal{A}_{N+1,2}$ denotes the so-called Kimura matrices [23], which is relevant to evolutionary dynamics. For this class, which includes the Moran and Wright-Fisher processes discussed in the introduction, a different presentation and notation were used: the matrix is naturally ordered with $i$ indicating the presence of a focal type, and the two absorbing states are labelled $i=0,N$ with

[TABLE]

Here $\bm{\alpha}$, $\bm{\beta}$ and $\mathbf{0}$ are $(N-1)\times 1$ vectors, with $\bm{\alpha}$ and $\bm{\beta}$ being nonnegative and nonzero.

By stochasticity of $\mathbf{M}$ we always have as an eigenvector $\mathbf{1}\mathrel{\mathop{:=}}(1,\dots,1)$, and we can choose a basis $\{\mathbf{F}_{1},\ldots,\mathbf{F}_{k}\}$ of the left-eigenspace of $\mathbf{M}$ (which is indeed $k$-dimensional from Definition 1) comprised of non-negative vectors such that $\sum_{j=1}^{k}\mathbf{F}_{j}=\mathbf{1}$. One readily checks that solutions to (3) automatically satisfy the conservation laws

[TABLE]

We refer to  [19, 21] for a discussion on how these conservation laws matter for the dynamics when considering the diffusive (continuous) approximation of Markov chains with absorbing states.

Remark 2**.**

In the particular case of Kimura matrices $\mathbf{M}\in\mathcal{A}_{N+1,2}$, the distinguished left eigenvectors are taken to be $\mathbf{1}$ and $\mathbf{F}$, with the first and the last entries of $\mathbf{F}$ being zero and one according to the representation of Equation (9). The vector $\mathbf{F}$ is the fixation probability of the focal type — see [23] and references therein and Remark 10 for the continuous version.

Remark 3** (Multi-type admissible processes).**

As we shall see below, the Moran process with $k$-types is not admissible because the matrix $\widetilde{M}$ is not irreducible. This is a consequence of non-interaction requirements that must be satisfied by the inner dynamics (when a type is extinct, it cannot reappear). On the other hand, it is possible to construct a birth-death process with $k$ absorbing states, such that $\widetilde{M}$ is irreducible. As an example, we mention a process akin to the Moran process with frequency dependent mutations. The easiest example of such a process can be briefly defined as follows: consider this process with a three-type population and set to zero all transition probabilities from the homogeneous states to any other state — this modification will still allow mutations from states with two types into states with three types. Such a process will have an irreducible $\widetilde{M}$ and three absorbing states. Naturally, this can be extended to $k$-type process, with $k$ absorbing states.

It will be also convenient to make explicit the structure of the semigroup associated to the forward Equation (3):

Lemma 1**.**

Let $\mathbf{M}$ be admissible. Then the fundamental solution to (3) is given by

[TABLE]

where we write indistinctly $\mathbf{I}$ for the identity matrix of dimension either $N+1$ or $N+1-k$.

Proof.

This follows from (7). ∎∎

As already discussed $\widetilde{\mathbf{M}}$ is strictly substochastic – the rows of $\mathbf{A}$ being non zero in (7) – and therefore its spectral radius

[TABLE]

The Perron-Frobenius theorem implies next that $\mu$ is the dominant eigenvalue of $\widetilde{\mathbf{M}}$, and both its associated left and right eigenvectors can be chosen positive, c.f. [64]. Following up on the rough idea that the transient dynamics determines the whole evolution, we define next

Definition 2** (Characteristic triple).**

Let $\mathbf{M}$ be admissible with core matrix $\widetilde{\mathbf{M}}$, and $\mu=\rho(\widetilde{\mathbf{M}})$. In addition, let $\widetilde{\mathbf{w}}$ and $\widetilde{\mathbf{z}}$ be the unique positive left and right eigenvectors associated to $\mu$, normalised as $\langle\widetilde{\mathbf{w}},\widetilde{\mathbf{1}}\rangle=\langle\widetilde{\mathbf{w}},\widetilde{\mathbf{z}}\rangle=1$. We will term $(\mu,\widetilde{\mathbf{w}},\widetilde{\mathbf{z}})$ the characteristic triple.

In general, it is difficult to derive explicit expressions for either $\mu$ or $\widetilde{\mathbf{w}},\widetilde{\mathbf{z}}$. Nevertheless, formulas can be obtained in certain particular cases or asymptotically in the limit of large populations $N\to\infty$, see Subsections 2.5 and 6.5.

2.2 Micro-reversible processes

In order to exploit the results from [79], we need to restrict ourselves to processes that satisfy some reversibility, at least to some extent. We introduce below a generalised notion of reversibility, adapted for the case of substochastic matrices, that we shall call micro-reversibility. Intuitively, micro-reversible processes should be time-reversible in the meta-stable regime. More precisely, micro-reversibility means that, when considering the difference between the mass flow from $i$ to $j$ and the mass flow from $j$ to $i$, the total loss of mass at site $i$ is independent of both $i$ and $j$, provided neither $i$ or $j$ are absorbing sates:

Definition 3** (Micro-reversible matrices).**

Let $\mathbf{M}$ be an admissible matrix with core $\widetilde{\mathbf{M}}$ in the sense of Definition 1 and with characteristic triple $(\mu,\widetilde{\mathbf{w}},\widetilde{\mathbf{z}})$ as in Definition 2. We say that $\mathbf{M}$ is micro-reversible if

[TABLE]

In matrix notation this amounts to requiring symmetry of $\mathop{\mathrm{Diag}}[\widetilde{\mathbf{w}}]\widetilde{\mathbf{M}}\mathop{\mathrm{Diag}}[\widetilde{\mathbf{z}}]$.

In a certain sense, this means that for these slow processes, a strong equilibrium relation is valid at each step, which resembles the quasi-stationary or ergodic processes in physics.

Note that this definition generalises the usual notion of reversibility for irreducible stochastic processes. Indeed for irreducible and column-stochastic Markov chains we have by definition $\widetilde{\mathbf{M}}=\mathbf{M}$, the left leading eigenvector is $\widetilde{\mathbf{w}}=\mathbf{w}=\mathbf{1}$, and our condition (11) reduces to the usual reversibility (detailed balance) $M_{ij}z_{j}=M_{ji}z_{i}$, $\forall i,j$. For irreducible row-stochastic matrices the right-eigenvector $\mathbf{z}=\mathbf{1}$, and reversibility reads instead $w_{i}M_{ij}=w_{j}M_{ji}$. In these cases, we say that $\mathbf{M}$ is column- or row-reversible, respectively. When no confusion arises we simply say that $\mathbf{M}$ is reversible, cf. [68, 86].

As we will see in a moment, micro-reversibility is satisfied at least for a particular class of processes, the birth-death processes with two absorbing states. This includes the Moran process (but not the Wright-Fisher one). For irreducible chains, birth-death processes with two absorbing states are among the simplest examples of reversibility, cf. [68]. For micro-reversibility, it is not difficult to prove that

Lemma 2**.**

Let $\mathbf{M}\in\mathcal{A}_{N+1,2}$ such that $\widetilde{\mathbf{M}}$ is tridiagonal. Then, $\mathbf{M}$ is micro-reversible.

Note that this is not true in general for admissible matrices with $k>2$ absorbing states, even if $\widetilde{\mathbf{M}}$ is tridiagonal.

Proof.

Let $(\mu,\widetilde{\mathbf{w}},\widetilde{\mathbf{z}})$ be the characteristic triple of $\mathbf{M}$. By assumption and up to permutation if needed, the core $\widetilde{\mathbf{M}}$ is irreducible and tridiagonal, hence from standard linear algebra [64] there exists a positive vector $\mathbf{d}$ such that $\mathbf{T}=\mathop{\mathrm{Diag}}[\mathbf{d}]\widetilde{\mathbf{M}}\mathop{\mathrm{Diag}}[\mathbf{d}^{-1}]$ is symmetric, i.e, $\widetilde{\mathbf{w}}\circ\mathbf{d}^{-1}=\mathbf{d}\circ\widetilde{\mathbf{z}}$. From the identity

[TABLE]

and the symmetry of $\mathbf{T}$, the micro-reversibility condition follows immediately. ∎∎

2.3 The associated $Q$-process

Given an absorbing Markov process $X_{t}$, the associated $Q$-process consists in conditioning it to non-absorption. Roughly speaking, the corresponding law at a fixed time $s\geq 0$ is given by $\lim_{T\to\infty}\mathbb{P}\left[X_{s}=x\,|\,\tau>T\right]$, where $\tau$ is the absorbing stopping time — cf. [75, 8, 15, 16, 28]. In our current setting, the importance of these processes is twofold: (i) there is a one-to-one correspondence between the transition matrix of the original and conditioned processes; (ii) the conditioned chain fits the framework of [79]. We refrain from going into the technical details and give instead a more direct definition in terms of linear algebra. The interested reader can check that the evolution of the above limiting process is indeed given by the transition matrix $\mathbf{K}$ below:

Lemma 3** ($Q$-process kernels).**

Let $\mathbf{M}\in\mathcal{A}_{N+1,k}$ with characteristic triple $(\mu,\widetilde{\mathbf{w}},\widetilde{\mathbf{z}})$. The associated $Q$-process is defined by its $(N+1-k)\times(N+1-k)$ transition matrix

[TABLE]

The kernel $\mathbf{K}$ is irreducible and row-stochastic, and its unique positive stationary probability distribution is given by

[TABLE]

Furthermore, if $\mathbf{M}$ is micro-reversible then $\mathbf{K}$ is row-reversible.

Proof.

Since $\widetilde{\mathbf{w}}$ is strictly positive and $\widetilde{\mathbf{M}}$ is irreducible, clearly $\mathbf{K}$ is irreducible. First, we check that $\bm{\pi}$ defined by (13) is indeed a left-eigenvector:

[TABLE]

By standard Perron-Frobenius theory $\bm{\pi}$ is thus the unique dominant eigenvector, and positive. From definition 2 we see that $\bm{\pi}$ is correctly normalised to be a probability vector. The fact that $\mathbf{K}$ is row-stochastic follows from $\mathbf{K}\mathbf{1}=\mathbf{1}$. The row-reversibility of $\mathbf{K}$, i.e. $\pi_{i}K_{ij}=\pi_{j}K_{ji}$ for all $i,j$ follows from the definition of $\bm{\pi}$ and $\mathbf{K}$, and from the micro-reversibility of $\mathbf{M}$, Equation (11). ∎∎

As already discussed, the transient dynamics leaks mass to the absorbed states. More explicitly, from (8) and because $\widetilde{\mathbf{w}}$ is a left-eigenvector, we have

[TABLE]

which implies

[TABLE]

In particular, since $\widetilde{\mathbf{M}}$ is substochastic with spectral radius $\mu\in(0,1)$ the $\widetilde{\mathbf{w}}$-weighted mass $\langle\widetilde{\mathbf{w}},\widetilde{\mathbf{p}}(t)\rangle$ of the transient states decays at an exponential rate $\mu-1<0$. Moreover, since $\widetilde{w}_{i}>0$ we see that either the initial data is completely absorbed $\widetilde{\mathbf{p}}(0)=\mathbf{0}$ and the dynamics trivially remains absorbed $\widetilde{\mathbf{p}}(t)\equiv\mathbf{0}$, or $\widetilde{\mathbf{p}}(0)\neq\mathbf{0}$ and therefore $\widetilde{\mathbf{p}}(t)\neq\mathbf{0}$ for all $t>0$. We therefore discard the trivial case $\widetilde{\mathbf{p}}(0)=\mathbf{0}$, and thus we can assume that $\langle\widetilde{\mathbf{w}},\widetilde{\mathbf{p}}(t)\rangle\neq 0$ for all times. We then define the two new $(N+1-k)$-dimensional variables

[TABLE]

By definition of $\mathbf{q}$ it is clear that $\sum_{i}q_{i}=1$. We think here of $\bm{\pi}=\widetilde{\mathbf{w}}\circ\widetilde{\mathbf{z}}$ as the new reference stationary measure, of $\mathbf{q}$ as a new probability evolving on a new time-scale $\nicefrac{{t}}{{\mu}}$, and $\mathbf{Q}=\frac{d\mathbf{q}}{d\bm{\pi}}$ as the density of $\mathbf{q}$ with respect to $\bm{\pi}$. J. Maas’ theory [79] of discrete Wasserstein distances rather take the density $\mathbf{Q}$ as a primary variable, while the $\mathbf{q}$ probability variable will be more convenient to address the diffusive limit of large populations later on. Direct substitution and elementary matrix algebra yields:

Lemma 4**.**

The following three dynamics are equivalent:

$\dfrac{\mathrm{d}\widetilde{\mathbf{p}}}{\mathrm{d}t}=\left(\widetilde{\mathbf{M}}-\mathbf{I}\right)\widetilde{\mathbf{p}}$; 2. 2.

$\dfrac{\mathrm{d}\mathbf{q}}{\mathrm{d}t}=\left(\mathbf{K}^{\dagger}-\mathbf{I}\right)\mathbf{q}$; 3. 3.

$\dfrac{\mathrm{d}\mathbf{Q}}{\mathrm{d}t}=\left(\mathbf{K}-\mathbf{I}\right)\mathbf{Q}$* .*

In addition, $\langle\mathbf{Q}(t),\bm{\pi}\rangle=\langle\mathbf{q}(t),\mathbf{1}\rangle=1$.

For the sake of brevity we omit the (elementary) proof.

It is worth pointing out that the change of timescale $\nicefrac{{t}}{{\mu}}$ in (15) is needed for notational convenience only, otherwise an additional factor $\mu$ would appear in the evolution laws below. Also, in the limit of large populations $N\to\infty$ considered later on, the subdominant eigenvalue $\mu=\mu_{N}\to 1$ so this rescaling becomes irrelevant.

2.4 Gradient flow formulation

In the previous section, we gave a canonical construction of the irreducible, reversible $Q$-process starting from the initial reducible, irreversible process. With irreducibility and reversibility newly satisfied by the transition matrix $\mathbf{K}$ of this $Q$-process (Lemma 3), we can now apply Maas’ theory [79] and identify the $\mathbf{q}$ evolution as a gradient flow for some discrete optimal transport structure (to be recalled in a moment).

Given an irreducible, reversible Markov kernel $\mathbf{K}$ (indexed as above by $i=0\dots N-k$) and its unique stationary distribution $\bm{\pi}$, the BGS entropy, also known as Kullback-Leibler divergence, of a probability $\mathbf{q}$ computed relatively to $\bm{\pi}$, is defined as

[TABLE]

Note that we use here the definition of entropy with reverted sign with respect to the historical definition – and also most common among physicists. Therefore, we expect its value to be non-increasing in time.

Let $\beta$ be the logarithmic mean

[TABLE]

J. Maas defined the following discrete optimal transport distance between probability densities:

Definition 4** (Discrete Wasserstein distance [79]).**

Let $\mathbf{K}$ be a stochastic, irreducible, and reversible transition kernel, and let $\bm{\pi}$ denote its unique stationary measure. Given two probability densities $\mathbf{Q}^{0}=\frac{\mathrm{d}\mathbf{q}^{0}}{\mathrm{d}\bm{\pi}},\mathbf{Q}^{1}=\frac{\mathrm{d}\mathbf{q}^{1}}{\mathrm{d}\bm{\pi}}$ with respect to $\bm{\pi}$, the discrete squared Wasserstein distance is

[TABLE]

where the infimum runs over all piecewise $\mathcal{C}^{1}$ curves of probability densities $\mathbf{Q}:[0,1]\ni\tau\mapsto\mathbf{Q}^{\tau}\in\mathbb{R}_{+}^{N}$ and all measurable functions $\bm{\psi}:[0,1]\ni\tau\mapsto\bm{\psi}^{\tau}\in\mathbb{R}^{N}$ satisfying the discrete continuity equation with endpoints $\mathbf{Q}^{0},\mathbf{Q}^{1}$

[TABLE]

This is a discrete counterpart to the celebrated Benamou-Brenier formula [11] for the continuous Wasserstein distance [107, 94], and we emphasise the dependence of $\mathcal{W}_{N}$ on $N$, since we will typically consider the large population limit $N\to\infty$ later on.

In what follows, we will slightly abuse notation by noticing that one can canonically define the discrete Wasserstein distance between probabilities in terms of their densities $\mathcal{W}_{N}(\mathbf{q}_{0},\mathbf{q}_{1})\mathrel{\mathop{:=}}\mathcal{W}_{N}\left(\frac{\mathrm{d}\mathbf{q}_{0}}{\mathrm{d}\bm{\pi}},\frac{\mathrm{d}\mathbf{q}_{1}}{\mathrm{d}\bm{\pi}}\right)$. In the sequel we will keep abusing the notations, and we shall simply speak of the discrete Wasserstein distance and Riemannian structure – see e.g. [39] for similar issues on the subtle distinction between measures and densities. Likewise, we will also write $H(\mathbf{q}|\bm{\pi})$ for $H\left(\mathbf{Q}\circ\bm{\pi}|\bm{\pi}\right)$.

Theorem 1** (Properties of the discrete Wasserstein distance [79]).**

With the same assumptions as before, denote by $\mathcal{P}_{\mathrm{dens}}$ the space of (finite) probability densities with respect to $\bm{\pi}$, and by $\mathcal{P}_{\mathrm{dens}}^{*}$ the subspace of everywhere strictly positive densities. Then

1. (i)

$\mathcal{W}_{N}$* defines a distance on $\mathcal{P}_{\mathrm{dens}}$;* 2. (ii)

The metric space $(\mathcal{P}_{\mathrm{dens}}^{*},\mathcal{W}_{N})$ is a Riemannian manifold; 3. (iii)

Given $f:\mathbb{R}^{+}\to\mathbb{R}$, the gradient of a functional $F(\mathbf{Q}):=\sum_{i}f\left(Q_{i}\right)\pi_{i}$ with respect to the Riemannian structure in (ii) reads, in local coordinates,

[TABLE] 4. (iv)

The three equations of Lemma 4 are equivalent to the gradient flow

[TABLE]

of the relative BGS entropy.

In (iii) and (iv) the intrinsic Riemannian gradient $\operatorname{grad}_{\mathcal{W}_{N}}F(\mathbf{Q})$ of a function $F:\mathcal{P}_{\mathrm{dens}}(\mathcal{X}_{N})\to\mathbb{R}$ is defined such that, along any differentiable curve $\tau\mapsto\mathbf{Q}^{\tau}$ with speed $\dot{\mathbf{Q}}^{\tau}\in T_{\mathbf{Q}_{t}}\mathcal{P}_{\mathrm{dens}}$, the chain rule holds as

[TABLE]

The precise definition of the scalar product $\langle.,.\rangle_{\mathbf{Q}}$ in the tangent plane $T_{\mathbf{Q}}\mathcal{P}_{\mathrm{dens}}$ at a point $\mathbf{Q}$ involves a certain weighted Onsager operator $L_{\mathbf{Q}}\mathrel{\mathop{:=}}-\mathbf{\operatorname{div}}_{\mathbf{Q}}(\bm{\nabla}\cdot)$, defined in terms of suitable discrete divergence and gradient operators $\mathbf{\operatorname{div}}_{\mathbf{Q}},\bm{\nabla}$. This allows to formally rewrite (20) in a more intrinsic fashion as

[TABLE]

This will have a clear counterpart later on in the continuous world, see in particular (54) and section 4.3. This Onsager operator precisely gives the one-to-one correspondence between $\psi^{\tau}$ and $\frac{\mathrm{d}\mathbf{Q}^{\tau}}{\mathrm{d}\tau}=L_{\mathbf{Q}^{\tau}}\psi^{\tau}$ implicitly appearing in the continuity equation (19); see [79, section 3]. The fact that (iv) holds is consequence of [79, Theorem 4.7]. This can be checked directly with (20): with $F(\mathbf{Q})=H(\mathbf{Q}\circ\bm{\pi}|\bm{\pi})=\sum_{i}\{Q_{i}\log Q_{i}\}\pi_{i}=\sum_{i}\{Q_{i}\log Q_{i}-Q_{i}+1\}\pi_{i}$ we have $f^{\prime}(Q_{i})=\log(Q_{i})$ in (20), hence

[TABLE]

and $\dfrac{\mathrm{d}\mathbf{Q}}{\mathrm{d}t}=-\operatorname{grad}_{\mathcal{W}_{N}}H(\mathbf{Q}\circ\bm{\pi}|\bm{\pi})$ reads indeed $\dfrac{\mathrm{d}\mathbf{Q}}{\mathrm{d}t}=\left(\mathbf{K}-\mathbf{I}\right)\mathbf{Q}$ as in Lemma 4.

As pointed out in [79], the restriction to positive densities in (ii) is not an issue: since the kernel $\mathbf{K}$ is irreducible any solution of the Heat Equation $\frac{\mathrm{d}\mathbf{Q}}{\mathrm{d}t}=(\mathbf{K}-\mathbf{I})\mathbf{Q}$ becomes instantaneously positive, $Q_{i}(t)>0$ for all $i$ and $t>0$.

We would like to stress that our main interest lies in (iv): although the original evolution $\frac{\mathrm{d}\mathbf{p}}{\mathrm{d}t}=(\mathbf{M}-\mathbf{I})\mathbf{p}$ is not truly speaking a gradient flow (due to absorbing states causing reducibility and irreversibility) one can in fact change the relevant variables so that the new effective (lower dimensional) $Q$-process kernel $\mathbf{K}$ becomes irreducible and reversible, and obtain a variational structure via discrete mass transport. Summarising the previous discussions, we have established in this section:

Theorem 2** (gradient flow structure for reducible irreversible kernels).**

Let $\mathbf{M}$ be admissible in the sense of Definition 1, let $\mathbf{K}$ be the associated $Q$-process with stationary measure $\bm{\pi}$, and let $\mathbf{q}(t)$ be defined by (15). Then the evolution of $\mathbf{p}$ is “variational” in the sense that $\mathbf{q}$ is driven by the gradient flow

[TABLE]

and that the complete dynamics of $\mathbf{p}$ can be uniquely recovered from $\mathbf{q}$.

Remark 4**.**

The fact that the evolution of $\mathbf{q}$ fully determines the dynamics of $\mathbf{p}$ can be seen by undoing the change of variables from $\widetilde{\mathbf{p}}$ to $\mathbf{q}$ and then using the conservation laws to recover $\mathbf{p}$, see Lemma 1 and Equation (15).

The above identification of the gradient flow structure only involved the BGS entropy $H(\mathbf{q}|\bm{\pi})$ as a driving functional. One can also define more general relative $\phi$-entropies of the form

[TABLE]

for a given convex function $\phi:\mathbb{R}^{+}\to\mathbb{R}$, still computed relatively to the reference measure $\bm{\pi}$. This covers the so-called Tsallis entropies, a generalisation of the Boltzmann-Gibbs entropy for non-additive systems with growing importance in biology (see [111] and references therein); see also [71] to its application in a model of prebiotic evolution akin to the replicator equation. Tsallis entropies will be further discussed in Section 3.

It turns out that the $G_{\phi}$ functional is a Lyapunov functional for the $\mathbf{q}$-evolution (which is rather driven by the $H$ entropy.)

Lemma 5**.**

Let $\mathbf{M}\in\mathcal{A}_{N+1,k}$ be an admissible and micro-reversible transition kernel, let $G_{\phi}$ be as in (23) for a given differentiable convex function $\phi:\mathbb{R}^{+}\to\mathbb{R}$, and consider a solution $\mathbf{Q}(t)$ of the previous gradient flow (Lemma 4 and Theorem 1). Then

[TABLE]

Moreover, if $\phi$ is locally strongly convex in the sense that $\phi^{\prime\prime}(x)\geq c_{M}>0$ in any bounded interval $x\in[0,M]$, then there exists $C>0$ depending only on $\phi,\bm{\pi}$ (but not explicitly on $\mathbf{K}$ or on the solution $\mathbf{Q}$) such that there holds the improved dissipation estimate

[TABLE]

Furthermore, if $\phi(x)\geq x-1$, then $G_{\phi}(\mathbf{Q}|\bm{\pi})\geq 0$.

We stress that the right-hand side in (24) is nothing but the expression in local coordinates of the Riemannian chain rule

[TABLE]

That this quantity is indeed nonnegative simply follows from our assumption that $\phi^{\prime}$ is nondecreasing, the summand in (24) being nonnegative for all $i,j$. In (25) we stress that $|.|^{2}_{\bm{\pi}}$ is the squared norm induced by the weighted scalar product $\langle a,b\rangle_{\bm{\pi}}=\sum a_{i}b_{i}\,{\pi}_{i}$ introduced earlier, which should not be confused with the Riemannian scalar product $\langle.,.\rangle_{\mathbf{Q}}$ at a point $\mathbf{Q}$. Note that the right-hand side vanishes if and only if $\mathbf{Q}$ is a stationary point, $\frac{\mathrm{d}\mathbf{Q}}{\mathrm{d}t}=0\Leftrightarrow\mathbf{Q}=\mathbf{K}\mathbf{Q}$.

Proof.

Let $\mathbf{K}$ be given by (12), with $\frac{\mathrm{d}\mathbf{Q}}{\mathrm{d}t}=(\mathbf{K}-\mathbf{I})\mathbf{Q}$ and $Q_{i}(t)=\frac{q_{i}(t)}{\pi_{i}}$. We compute

[TABLE]

where the last equality simply comes from the stochasticity $\sum_{j}K_{ij}=1$. Leveraging the microreversibility $K_{ij}\pi_{i}=K_{ji}\pi_{j}$, a straightforward summation by parts leads to

[TABLE]

Recalling the definition of the logarithmic mean, Equation (17), this proves Equation (24).

Assume now that $\phi$ satisfies the additional strong local convexity as in our statement. Just as before we can write

[TABLE]

Observe that due to $Q_{i}=\frac{q_{i}}{\pi_{i}}\leq\frac{\sum_{j}q_{j}}{\pi_{i}}=\frac{1}{\pi_{i}}$ we have $0\leq Q_{i}\leq M:=\frac{1}{\min_{j}\pi_{j}}<+\infty$. Recalling that $\mathbf{K}$ is row-stochastic we see that the convex combination $\sum_{j=0}^{N-k}K_{ij}Q_{j}\leq M$ as well, for any fixed $i$. Exploiting our strong convexity assumption on $\phi$, a straightforward application of Taylor’s theorem gives that

[TABLE]

where $c_{M}>0$ is a lower bound for $\phi^{\prime\prime}(x)$ in the interval $[0,M]$. Substituting in (26) gives then, by convexity of $\phi$ and reversibility $K_{ij}\pi_{i}=K_{ji}\pi_{j}$,

[TABLE]

as desired.

Finally, the fact that $G_{\phi}(\mathbf{Q}|\bm{\pi})\geq 0$ immediately follows from

[TABLE]

where we used $\sum_{i}Q_{i}\pi_{i}=\sum_{i}q_{i}=\sum_{i}\pi_{i}=1$. ∎∎

We will state later on in Section 3 a discrete-in time monotonicity of the entropy, Proposition 1. This discrete-time framework will be based on a “Euclidean” discretization (linear 1st order difference quotient), and we believe that in this context a discrete equivalent of (25) can be established. However, (24) strongly leverages the Riemannian structure from discrete optimal transport, and seems therefore incompatible with the Euclidean time discretization. An alternative and natural discretization in time would rather be the JKO scheme, more compatible with the discrete Wasserstein and gradient flow structure (see Section 6 for a short discussion). For such JKO discretization, it should be possible to write down a discrete version of (24). However, even if this turns out to be possible, the relationship of this new discrete process with the embedded chain is not clear, and we will not pursue this direction here.

We also point out that $G_{\phi}=H$ is admissible in Lemma 5, in which case one is actually computing the dissipation of $H$ along its own gradient flow and the right-hand side of (24) is a discrete version of the Fisher information $F(\rho|\pi)=\int\left|\nabla\log\left(\frac{\rho}{\pi}\right)\right|^{2}\rho$.

2.5 The continuous in time Moran process

In this subsection, we study the Moran process with the techniques developed above. As already explained, the only two absorbing states are labelled here $i=0$ and $i=N$. We start with explicit results for the so-called neutral evolution, in which both types $\mathbb{A}$ and $\mathbb{B}$ have the same reproductive viability, namely, $\mathsf{s}_{i}^{(\mathsf{n})}=i/N$ in the transition matrix given by Equation (1). We indicate the neutral evolution in this subsection by the superscript $\mathsf{n}$. We will obtain an explicit formula for the Wasserstein distance between adjacent sites in birth and death processes, and finish with some comments on more general evolutionary processes including the Wright-Fisher process.

For the neutral Moran process, a simple calculation shows that $z_{i}^{(\mathsf{n})}=\frac{1}{N-1}$, $w_{i}^{(\mathsf{n})}=\frac{6i(N-i)}{N(N+1)}$ and $\lambda^{(\mathsf{n})}=-\frac{2}{N^{2}}$. Introducing the auxiliary notation $x_{i}=i/N$ and considering the BGS entropy given by Equation (16) as a function of $\mathbf{p}$ (rather than $\mathbf{q}$) through the change of variables (15), we arrive at the entropy for the neutral evolution:

[TABLE]

For non-neutral Moran process, one should not expect analytical formulas for $\lambda$, $\mathbf{w}$ and $\mathbf{z}$. However, as a consequence of Lemma 2 and Equation (11), we find that

[TABLE]

Therefore

[TABLE]

where $C$ is a certain normalisation constant. This equation will be used in Subsection 6.5 to identify the correct macroscopic limit of the Moran process, i.e., the correct model when $N\to\infty$.

For general Moran processes, we are able to compute explicitly the Wasserstein distance between adjacent sites. Namely,

Lemma 6**.**

Let $\mathbf{K}$ be a tridiagonal, irreducible, and reversible transition matrix, let $\beta(x,y)$ be the logarithmic mean defined in (17), and let us write $\mathbf{e}_{i}$ for the discrete probability vector concentrated on the $i$-th state. For adjacent sites $j=i+1$ we have

[TABLE]

Proof.

Since $\mathbf{K}$ is tridiagonal, the minimising curve in Definition 4 should only involve the two neighbouring sites, i.e. the density $\mathbf{Q}^{\tau}=r(\tau)\frac{\mathbf{e}_{i}}{\pi_{i}}+(1-r(\tau))\frac{\mathbf{e}_{i+1}}{\pi_{i+1}}$ for some function $r$ with $r(0)=1$ and $r(1)=0$ to be determined. From (19), we conclude that

[TABLE]

where $\dot{}=\frac{\mathrm{d}\ }{\mathrm{d}\tau}$. Plugging this into the action (18) and leveraging the reversibility of $\mathbf{K}$, we find a functional of $r$ only. Namely

[TABLE]

Thus $\mathcal{W}_{N}^{2}(\mathbf{e}_{i},\mathbf{e}_{i+1})$ is the minimal value of this functional over all $r\in\mathcal{C}^{1}([0,1])$, such that $r(0)=1$, $r(1)=0$, the Euler-Lagrange Equation implies that

[TABLE]

i.e., $r^{\prime}(\tau)=c\sqrt{F(r(\tau))}$, where $F(r)=\beta(K_{i,i+1}r,K_{i+1,i}(1-r))$ and

[TABLE]

The result follows immediately. ∎∎

Note that the fact that the curve linking two neighbouring states is a linear combination of $\mathbf{e}_{i}$ and $\mathbf{e}_{i+1}$ only holds here because the matrix is tridiagonal. Otherwise the nonlocal effects start building up and many other states play a role too, and the problem becomes more complicated than the optimisation over scalar functions $r(\tau)$ above.

We finish with one more comment about general processes in the Kimura class:

Remark 5**.**

For general $\mathbf{M}\in\mathcal{A}_{N+1,2}$, the micro-reversibility does not necessarily hold. In particular, for the Wright-Fisher process this is not true even in the neutral case. However, for sufficiently large $n$

[TABLE]

where $\mathbf{F}$ is the unique solution of $\mathbf{F}^{\dagger}=\mathbf{F}^{\dagger}\mathbf{M}$ with $F_{0}=0$, $F_{N}=1$. Therefore ${(\mathbf{M}^{n})_{ij}\approx\lambda^{n}z_{i}w_{j}}$ for $i,j=1,\dots,N-1$, i.e., $w_{i}M_{ij}^{n}z_{j}\approx w_{j}M_{ji}^{n}z_{i}$. In a loose sense, we would say that any process in the Kimura class is asymptotically micro-reversible. However, we will not explore these ideas in the current work.

3 Interlude: discrete in time, discrete Markov chains

The notion of relative entropy, and even the entropy itself, is usually defined for irreducible Markov processes only, see e.g. [103] for BGS entropy, and [50] for Tsallis entropies. Due to the existence of absorbing states in the transition kernel, this notion of course does not apply here to our admissible matrices.

In this section, we study discrete in time Markov processes. Discrete in time models cannot be recast in the gradient flow formalism, however we opt to include this discussion in the present work because, as for the Moran process, these are frequently used in population genetics. Furthermore, the results of the previous section can shed light on the precise notion of entropy for processes with two or more absorbing states.

We shall see that the right notion of entropy only depends on the transient states, and we will therefore speak of substochastic entropies. We would like to point out that this issue is a priori non trivial: when $\mathbf{M}$ has absorbing states, usual entropy functionals often lead to quantities that increase along the evolution. (We recall that with our sign conventions entropies should rather be nonincreasing in time.)

We first show in Subsection 3.1 how to make discrete in time and continuous in time models compatible, given a certain transition matrix $\mathbf{M}$ for the continuous in time model. In the sequel, we show that the entropy functional is independent of the time step and show how a generalisation of the BGS and Tsallis entropies for substochastic models naturally arises from our analysis.

3.1 The embedded chain

Every stochastic matrix $\mathbf{M}$ can always be seen as a discrete Markov chain in either discrete or continuous time. Indeed, $\mathbf{M}$ can be seen as a transition matrix and the evolution of $\mathbf{p}$ is given by Equation (2); alternatively, it can be seen as the so-called embedded chain, whose evolution is given by Equation (3), see [91]. Let $h>0$ be a small time step and define the kernel

[TABLE]

Let also $\mathbf{p}^{(h)}(t)$ be the (piecewise constant) time-interpolation defined by the recursion

[TABLE]

with $\mathbf{p}^{(h)}(0)=\mathbf{p}^{\mathrm{I}}$. From the fact that

[TABLE]

we see that (31) is of course the explicit Euler discretization of (3). Thus, the following convergence is not surprising:

Lemma 7**.**

In the limit of small time steps we have

[TABLE]

where $\mathbf{p}(t)$ is the solution of Equation (3) with $\mathbf{p}(0)=\mathbf{p}^{\mathrm{I}}$ and the convergence is locally uniform in time.

Proof.

The convergence of the Euler scheme is a standard result, cf. [104]. In our particular context of linear ODEs, a simple proof consists in writing $n=\left\lfloor t/h\right\rfloor$, and noticing that $\lim_{h\to 0}\left(\mathbf{M}^{(h)}\right)^{n}=\exp\left(t\left(\mathbf{M}-\mathbf{I}\right)\right)$. ∎∎

3.2 Dynamics and generalised entropies

In this section we only consider admissible transition kernels, $\mathbf{M}\in\mathcal{A}_{N+1,k}$.

Note that the entropy $H(\mathbf{q}|\bm{\pi})$, given by (16) as a function of the $(N+1)-k$ dimensional $Q$-process $\mathbf{q}$, depends on $\mathbf{M}$ only through the characteristic eigenvectors $\widetilde{\mathbf{w}},\widetilde{\mathbf{z}}$ - since by definition the reference measure $\bm{\pi}=\widetilde{\mathbf{w}}\circ\widetilde{\mathbf{z}}$. With this fact in mind, it is clear that changing the transition kernel from $\mathbf{M}$ to $\mathbf{M}^{(h)}$, for any $h>0$ will not change the entropy. More precisely,

Lemma 8**.**

Let $(\mu,\widetilde{\mathbf{w}},\widetilde{\mathbf{z}})$ be the characteristic triple of $\mathbf{M}$ in the sense of Definition 2, and let $\mathbf{M}^{(h)}$ be the discrete-in-time transition kernel defined in (30). Then $\mathbf{M}^{(h)}$ has characteristic triple $(\mu_{h},\widetilde{\mathbf{w}},\widetilde{\mathbf{z}})$, with

[TABLE]

In particular $\mathbf{M}$ and $\mathbf{M}^{(h)}$ share their characteristic eigenvectors. Moreover, $\mathbf{M}^{(h)}$ is micro-reversible if and only if $\mathbf{M}$ is. Finally, if $\mathbf{M}\in\mathcal{A}_{N+1,2}$ is a Kimura matrix with fixation probability $\mathbf{F}$, then so is $\mathbf{M}^{(h)}$.

Proof.

This easily follows from the expression (30) for $\mathbf{M}^{(h)}$. ∎∎

For notational convenience, we extend the previous transient reference measure $\bm{\pi}=(\pi_{i})_{i=0,\dots,N-k}$ by $k$ zeros to form the corresponding full $(N+1)$-dimensional probability measure

[TABLE]

It is also natural to extend definition (23) to the full dynamics, i.e., by abuse of notation, using the variable $\mathbf{p}$ instead of $\mathbf{q}$, which are equivalent in view of Remark 4, we write

[TABLE]

This definition is completely independent of the time step $h$.

Remark 6**.**

The entropy is defined on $\mathbf{q}$, and therefore on the transient states $\widetilde{\mathbf{p}}$, and not on the whole probability $\mathbf{p}$. In terms of measure theoretic considerations, our definition can be summarised as follows: performing a Lebesgue decomposition with respect to the reference measure $\overline{\bm{\pi}}$, i.e. $\mathbf{p}=\widetilde{\mathbf{p}}+\mathbf{p}^{\mathrm{s}}$ with $\widetilde{\mathbf{p}}\ll\overline{\bm{\pi}}$ and $\mathbf{p}^{\mathrm{s}}\perp\overline{\bm{\pi}}$, our entropy is finite if and only if $\widetilde{\mathbf{p}}\neq 0$. This departs from the usual definition of relative entropies, where one usually sets the entropy to $+\infty$ if $\mathbf{p}$ is not absolutely continuous w.r.t. the reference measure $\overline{\bm{\pi}}$ (i.e. if the singular part $\mathbf{p}^{\mathrm{s}}\neq 0$). Here we have a whole freedom between those two scenarios, and our entropy takes finite values even for $\mathbf{p}$’s that are partially absorbed with both $\widetilde{\mathbf{p}}\neq 0$ and $\mathbf{p}^{\mathrm{s}}\neq 0$. Note also that our reference measure $\overline{\bm{\pi}}$ does not have full support.

For the reader aiming to apply the previous results in the traditional (i.e., discrete time) Moran process, we state the entropy inequality in the natural variables. Namely, taking as primary variables the probability distribution $\mathbf{p}$ and considering generalised entropies, we state the following result.

Proposition 1** (Discrete monotonicity of generalised entropies).**

Let $\mathbf{M}\in\mathcal{A}_{N+1,k}$ be an admissible and micro-reversible transition kernel, let $\phi:\mathbb{R}^{+}\to\mathbb{R}$ be a convex function, and let $G_{\phi}(\cdot|\overline{\bm{\pi}})$ be the associated generalised entropy from Definition 23. Then the one-step monotonicity

[TABLE]

holds for any sufficiently small time step $h$ (depending only on $\mathbf{M}$ but not on $\mathbf{p}$). Furthermore, if $\phi(x)\geq x-1$, then $G_{\phi}(\mathbf{p}|\overline{\bm{\pi}})\geq 0$ for any probability vector $\mathbf{p}$.

Proof.

Let $(\mu,\widetilde{\mathbf{w}},\widetilde{\mathbf{z}})$ be the characteristic triple of $\mathbf{M}$ and $h\in(0,1/(1-\mu))$. If $G_{\phi}(\mathbf{p}|\overline{\bm{\pi}})=+\infty$ there is nothing to prove, hence we only consider $\langle\widetilde{\mathbf{w}},\widetilde{\mathbf{p}}\rangle>0$. Let us first prove that $\mathbf{M}^{(h)}\mathbf{p}$ has finite entropy, i.e. $\langle\widetilde{\mathbf{w}},\widetilde{\mathbf{M}^{(h)}\mathbf{p}}\rangle>0$ as well. To this end we first extend $\widetilde{\mathbf{w}}$ to $\overline{\mathbf{w}}$ by zeros. Then, since $\langle\widetilde{\mathbf{w}},\widetilde{\mathbf{p}}\rangle=\langle\overline{\mathbf{w}},\mathbf{p}\rangle=\overline{\mathbf{w}}^{\dagger}\mathbf{p}$ for any $\mathbf{p}$, we get

[TABLE]

if the time step $h\in(0,1/(1-\mu))$ is small enough to guarantee $\mu_{h}=1-h(1-\mu)>0$. This proves that $\mathbf{M}^{(h)}\mathbf{p}$ cannot be completely absorbed if $\mathbf{p}$ is not, hence $G_{\phi}(\mathbf{M}^{(h)}\mathbf{p}|\overline{\bm{\pi}})<\infty$.

The remainder of the proof follows from Lemma 8, identity (33), and the convexity of $\phi$:

[TABLE]

∎∎

Remark 7**.**

Choosing $\phi(x)=x\log x\geq x-1$ we recover the previous expression (16), providing a generalisation of the relative Boltzmann-Gibbs-Shannon entropy for reducible processes. On the other hand, choosing $\phi(x)=\frac{x^{m}-x}{m-1}\geq x-1$ for exponents $m>1$ also provides a generalisation of the relative Tsallis entropies for reducible kernels, i.e.

[TABLE]

When $\mathbf{M}$ is irreducible we have $k=0$, $\widetilde{\mathbf{w}}=\mathbf{1}$, and the reference measure is $\overline{\bm{\pi}}=\widetilde{\mathbf{z}}$: in this case our generalised Tsallis entropies coincides with that in [50]. Let us also point out that, at the continuous level $x\in[0,1]$, both the relative Kullback-Leibler divergence $\mathcal{H}(\rho|\pi)=\int\rho\log(\rho/\pi)\mathrm{d}x=\int\rho\log\rho\,\mathrm{d}x+\int\rho V\,\mathrm{d}x$ (with Gibbs distribution $\pi(x)=e^{-V(x)}/Z$) and the Tsallis entropies $\mathcal{E}_{m}(\rho)=\int\frac{\rho^{m}-\rho}{m-1}\mathrm{d}x$ play a particular role in continuous optimal transport, since the Fokker-Planck Equation $\partial_{t}\rho=\Delta\rho+\operatorname{div}(\rho\nabla V)$ and the Porous Medium Equation $\partial_{t}\rho=\Delta\rho^{m}$ can be viewed as their respective Wasserstein gradient flows, see [66, 92, 94, 107] and the next section.

Remark 8**.**

Results from Lemma 5 apply for Tsallis entropies in the range $m\in(1,2]$. This particular range seems to be of relevance to applications in ecology, as one possible interpretation of $m$ is as an interpolation parameter between two well established diversity measures for populations: the Shannon-Wiener index (the BGS entropy, in our notation) in the limit $m\to 1$ and the Simpson index, defined as the limit $m\to 2$ of the Tsallis entropies; see [111].

4 Continuous time, continuous Markov chains

In this section we take interest in the gradient flow formulation of the Kimura Equation (5). The relation between continuous and discrete models will be investigated in Section 6, but the reader may want to keep in mind the following picture: The Kimura Equation is the continuous counterpart of the Moran process, where $x\in[0,1]$ plays the role of $i/N$, i.e., the fraction of individuals of a given type in the population. Although motivated by this large population limit, we focus in this section on a self-contained presentation of a slightly more general setting than the specific Kimura PDE (5).

We start by defining this generalised Kimura Equation in Subsection 4.1. We continue with the exploration of the two ingredients of the gradient flow: an entropy in Subsection 4.2 and a distance in Subsection 4.3. We then precisely define the gradient flow for the Kimura Equation in Subsection 4.4.

4.1 Background

The Kimura Equation is a particular example of a stochastic differential equation of the form

[TABLE]

$W_{t}$ denoting the standard one-dimensional Brownian, cf. [41, 24] and the study of this class is of paramount importance in population genetics, see e.g. [43].

One particular feature of the models that we are interested in is that the dispersion coefficient $\sigma(x)$ is positive in the interior $x\in(0,1)$, but with $\sigma(0)=\sigma(1)=0$, representing two absorbing states on the boundary.

For the resulting PDEs, the presence of those two absorbing states makes the notion of measure-valued solutions natural and actually necessary, and the probability laws typically take the form

[TABLE]

Remark 9**.**

In the same spirit as in the previous sections, the continuous part $\widetilde{p}(t,x)$ corresponds to the previous transient densities $N\widetilde{\mathbf{p}}=(Np_{i})_{i=1,\dots,N-1}$, while the boundary contributions $a\delta_{0},b\delta_{1}$ correspond to the previous absorbed states $p_{0},p_{N}$.

In this work, we call the Kimura Equation a generalised version of Equation (5), i.e, given by

[TABLE]

where $\kappa>0$ is a diffusion coefficient and $V:[0,1]\to\mathbb{R}$ is the gradient potential as defined in Subsection 1.5. Furthermore,

[TABLE]

The typical case arising e.g. in (5) is $\Theta(x)=x(1-x)$, see Subsection 1.5. The fact that the zeroes are simple is crucial: the random motion is strong enough to counteract the deterministic advection driven by the velocity field $-\Theta(x)\partial_{x}V$, and hence the trajectories are absorbed in finite time almost surely. In terms of PDEs, the diffusion is locally uniform in any compact set $K\subset(0,1)$, but degenerates at the boundaries.

As is common, we shall refer to the operator

[TABLE]

in (35) as the forward operator, while we speak of its formal adjoint

[TABLE]

as the backward operator. If not required from the context, the index $\kappa$ will be omitted from the operators $\mathcal{L}$ and $\mathcal{L}^{\dagger}$.

According to (35) it is clear that the interior density $\tilde{p}(t,x)$ in (34) evolves according to

[TABLE]

However, because the diffusion is degenerate $\Theta(0)=\Theta(1)=0$, the evolution problem (35) cannot be supplemented with standard boundary conditions as usual. Additional conservation laws must be used instead to make sense of the Cauchy problem, and those are reminiscent from the discrete conservation laws (10) in Section 2. In [19, 21] two of the authors studied the forward equation, and obtained a characterisation of the boundary measures in terms of those conservation laws. More precisely, let $F$ be the solution of

[TABLE]

Remark 10**.**

This eigenfunction $F(x)$ is known as the fixation probability, which encodes the probability of the population ending with a homogeneous population of type $\mathbb{A}$, when starting from an initial ratio $x$ (hence $F(0)=0$ and $F(1)=1$ for the two absorbing states) – see also remark 2 at the discrete level.

Imposing two additional conservation laws for the total mass and fixation

[TABLE]

makes Equation (35) well-posed [20], given an initial condition $p^{\mathrm{I}}\in\mathcal{BM}^{+}([0,1])$ – the space of positive Radon measures in $[0,1]$. We should stress that both integrals are computed in the closed interval $x\in[0,1]$, since the measure $p(t)$ may – and typically does – charge the boundaries. With those conservation laws newly enforced, the resulting Cauchy problem is studied in [20] by means of (weighted, singular) Sturm-Liouville theory. There, it was proved that the transient component ($\widetilde{p}$ in Equation (34)) is smooth up to the boundary, $\widetilde{p}\in C^{\infty}(0,T;C^{k+1}[0,1])$ if $\Theta,V\in C^{k}([0,1])$. More importantly, the two global conservation laws (40) are equivalent to the two local flux conditions

[TABLE]

driving the loss of mass from the continuous inner (transient) component $\tilde{p}(t,x)\mathrm{d}x$ towards the boundary (absorbing) points $a\delta_{0}+b\delta_{1}$. These flux conditions can be heuristically obtained by inserting the representation for $p$ given in Equation (34) into the conservation laws given by Equation (40):

[TABLE]

and

[TABLE]

where we used that the fixation $F(x)$ satisfies by definition $F(0)=0$ and $F(1)=1$. Using the evolution law $\partial_{t}\tilde{p}=\mathcal{L}\tilde{p}$ for the smooth interior part, the stationary elliptic equation (39) for the fixation $F$, and integrating by parts in the inner integrals then gives

[TABLE]

For more details, see [20].

4.2 Entropy

Observe that we introduced the discrete entropy (16) in order to cope with the irreducibility of the relevant Markov chain. Therefore, we expect that simple definitions of the entropy (such as, e.g., $\int p\log p$) should not work either in the continuous setting. In fact, our construction in the discrete case crucially relied upon the (sub)dominant characteristic triple $(\mu,\widetilde{\mathbf{w}},\widetilde{\mathbf{z}})$ of $\mathbf{M}$ (definition 2). Therefore, our first step will be to introduce the corresponding continuous objects:

Definition 5** (Continuous characteristic triple).**

The characteristic triple is the triple $(\lambda,w,z)$ defined as the principal eigenvalue/eigenfunctions $-\mathcal{L}z=\lambda z$ and $-\mathcal{L}^{\dagger}w=\lambda w$ with homogeneous Dirichlet boundary conditions, i-e

[TABLE]

and

[TABLE]

We always choose $w,z$ to be positive in $(0,1)$ and normalised as

[TABLE]

The fact that those principal eigenfunctions/values are well-defined follows from standard spectral theory, after observing that the above Sturm-Liouville problems are of limit-circle-non-oscillatory type, see e.g. [112]. The eigenvalue $\lambda>0$ will quantify the exponential decay of $\widetilde{p}(t)$. Just as in the discrete case, we define the new reference measure

[TABLE]

We slightly abuse the notations and identify the measure $\pi$ to its density $\pi(x)$ with respect to the Lebesgue measure. Equation (43) should be be compared to the discrete counterpart (13). Our normalisation of $z$ and $w$ yields $\int_{0}^{1}\mathrm{d}\pi(x)=1$, and we view $\pi\in\mathcal{P}(0,1)$ as a reference probability measure. The measure $\pi$ only charges $x\in(0,1)$, but should in fact be viewed as a probability measure in the whole underlying space $[0,1]$. A useful information on this reference measure will be

Lemma 9**.**

With the same notations, there holds

[TABLE]

where $C=\int_{0}^{1}\frac{w^{2}(x)\mathrm{e}^{-2V(x)/\kappa}}{\Theta(x)}\mathrm{d}x>0$ is a normalising constant such that $\int\mathrm{d}\pi(x)=1$.

Proof.

Defining $Z(x):=\frac{w(x)\mathrm{e}^{-2V(x)/\kappa}}{\Theta(x)}$ and exploiting $-\mathcal{L}^{\dagger}w=\lambda w$, straightforward algebra leads to $-\mathcal{L}Z=\lambda Z$. Since $Z>0$ satisfies the boundary conditions $\lim\limits_{x\to\{0,1\}}\Theta(x)Z(x)=0$, we see by uniqueness of the principal eigenfunction in Definition 5 that necessarily $Z(x)=Cz(x)$ for some constant $C>0$. Hence $wz=\frac{1}{C}wZ$ takes the desired form, and the value of $C>0$ follows from our normalisation $\int\mathrm{d}\pi=1$.∎∎

Standard results from measure theory allow to decompose any probability $p\in\mathcal{P}([0,1])$ as

[TABLE]

Whenever $\widetilde{p}=\widetilde{p}(x)\mathrm{d}x$ is absolutely continuous and the $w$-weighted mass

[TABLE]

is non-zero, we define the renormalised transient probability measure

[TABLE]

Mimicking our definition in the discrete case, we define the relative entropy by setting

[TABLE]

Whenever $q\ll\pi$ we have of course $\frac{\mathrm{d}q}{\mathrm{d}\pi}(x)=\frac{q(x)}{\pi(x)}$. For the dynamics, the renormalised $q(t,x)$ variable will be derived from the original $p(t,x)$ probability through (45), but we will in fact use $q$ as a primary variable/unknown as before in the discrete setting. Formula (46) accordingly defines the entropy on the whole space $q\in\mathcal{P}(0,1)$, whether $q$ actually arose from some $p$ or not. Moreover, the $\kappa/2$ factor appears in (46) due to the particular diffusion scaling in (35). The reference measure $\pi=\pi_{\kappa}$ and the entropy functional $\mathcal{H}=\mathcal{H}_{\kappa}$ both depend on $\kappa$. We shall in fact take the so-called deterministic limit $\kappa\to 0$ later on, but we dispense at this stage from tracking the $\kappa$-dependence.

Just as in the discrete case, one can define more general entropies in terms of the original $p$ measure, which we called previously the substochastic or reducible entropies. Similarly to (23) and (32), if $\phi:\mathbb{R}^{+}\to\mathbb{R}$ is a convex, superlinear function, those read

[TABLE]

The above expression should be understood to be $+\infty$ whenever $q\not\ll\pi$, and makes sense for general $q$. However, when $q=\frac{w\widetilde{p}}{\langle w,\widetilde{p}\rangle}$ is obtained as the $Q$-process corresponding to some $p$ with $\langle w,\widetilde{p}\rangle\neq 0$ (which propagates from $t=0$ to later times as in the discrete case, see also (55) below), and recalling that $\pi(x)=w(x)z(x)$, we abuse the notations and also express this same entropy in terms of the original $p$ measure as

[TABLE]

The continuous BGS entropy (46) corresponds of course to the particular choice $\phi(\eta)=\frac{\kappa}{2}\eta\log\eta$, and one can also consider Tsallis entropies $\phi(\eta)=\frac{\kappa}{2}\frac{\eta^{m}-\eta}{m-1}$.

4.3 The Wasserstein distance

In this section we introduce a suitable Wasserstein distance in the space of probability measures $\mathcal{P}([0,1])$ that will allow to write the Kimura Equation as a gradient flow. Due to the presence of the variable coefficient $\Theta(x)$ in (35), this quadratic Wasserstein distance will not be based on the usual Euclidean distance, but will rely instead upon viewing the underlying $\Omega=(0,1)$ as a suitable Riemannian manifold. More precisely, we consider the Riemannian metric with scalar product on the tangent plane $T_{x}\Omega$ at a point $x\in\Omega$ induced by $\frac{1}{\Theta(x)}$, namely the norm of a tangent vector $\zeta\in T_{x}\Omega$ is defined as

[TABLE]

The induced generalised Shahshahani distance is

[TABLE]

for $x,y\in\Omega$. As can be expected, this distance is well-behaved:

Lemma 10**.**

Assume that $\Theta$ satisfies (36). Then the infimum in (49) is achieved for a unique constant-speed geodesic, $\mathsf{d}$ can be uniformly extended to $\overline{\Omega}^{2}$, and

[TABLE]

Moreover, $\mathsf{d}$ defines a distance in $\overline{\Omega}$ and the metric space $(\overline{\Omega},\mathsf{d})$ is Polish.

Proof.

For $x,y\in\Omega$ in the interior, the existence of a unique minimising curve is a standard exercise in the calculus of variations and we omit the details. As in the proof of Lemma 6, we start by writing the Euler-Lagrange Equation $\frac{\mathrm{d}\ }{\mathrm{d}t}\left(\frac{\xi^{\prime}(t)^{2}}{\Theta(\xi(t))}\right)=0$ and conclude that the intrinsic speed is constant, i.e., $|\xi^{\prime}(t)|^{2}_{T_{\xi(t)}\Omega}=\mathsf{d}^{2}(x,y)$. In particular $\xi^{\prime}(t)$ never vanishes and (50) immediately follows from the change of variables $u=\xi(t)$ in (49). From the fact that $x=0,1$ are simple zeros of $\Theta$ as assumed in (36), the extension to $\overline{\Omega}$ follows; finally, completeness is an easy consequence of the explicit representation (50). ∎∎

It is worth pointing out that this Shahshahani distance is locally equivalent to the Euclidean one in the interior (i.e. $c_{K}|x-y|\leq\mathsf{d}(x,y)\leq C_{K}|x-y|$ in any compact set $K\subset\subset(0,1)$), but behaves differently close to the boundary (e.g. $\mathsf{d}(0,x)\sim\int_{0}^{x}\frac{\mathrm{d}u}{\sqrt{u}}\sim\sqrt{x}$ for small $x$). This is reflected in the behaviour of the Kimura Equation (35), which is locally uniformly parabolic in the interior, but degenerate at the boundaries. With the Polish space $(\overline{\Omega},\mathsf{d})$ at hand, one classically defines the corresponding Wasserstein distance on the space of probabilities $\mathcal{P}(\overline{\Omega})$ as

[TABLE]

Here $\Gamma(\mu,\nu)$ denotes the set of admissible transport plans, i.e., the set of probability measures $\gamma\in\mathcal{P}(\overline{\Omega}\times\overline{\Omega})$ with first marginal $\gamma_{1}=\mu$ and second marginal $\gamma_{2}=\nu$. The superposition principle

[TABLE]

gives the natural correspondence between the underlying Polish space $(\overline{\Omega},\mathsf{d})$ and the overlying Wasserstein space $(\mathcal{P}(\Omega),\mathcal{W})$, and we refer to [107, 108, 94] for an extended account on the optimal transport theory and bibliography.

As in the discrete case, we have the dynamical representation

Proposition 2** (Benamou-Brenier formula [11, 78]).**

For $q_{0},q_{1}\in\mathcal{P}(\overline{\Omega})$ there holds

[TABLE]

where the infimum runs over narrowly continuous curves $[0,1]\ni t\mapsto q_{t}\in\mathcal{P}(\overline{\Omega})$ with endpoints $q|_{t=0}=q_{0},\,q|_{t=1}=q_{1}$ and satisfying the continuity equation

[TABLE]

with zero-flux boundary conditions.

This is the exact equivalent of the dynamical definition of the discrete Wasserstein distance – Definition 4 – where the discrete continuity equation appear (see [79] for discussions). In the Lagrangian action (53) the velocity-field $\mathbf{v}=\nabla\psi$ is measured not with respect to the standard Euclidean norm, but rather with respect to the intrinsic Shahshahani metrics $|\nabla\psi|^{2}_{T_{x}\Omega}=\frac{|\nabla\psi|^{2}}{\Theta(x)}$. We refer to [78] for a discussion on Wasserstein distances with variable coefficients, and to [108] for optimal transport on abstract Riemannian manifolds. Since we called the underlying metrics $\mathsf{d}$ the generalised Shahshahani distance, we shall sometimes speak of the corresponding Wasserstein distance as the Wasserstein-Shahshahani distance.

From the works of Otto [92] it is known that the Wasserstein distance endows $\mathcal{P}(\Omega)$ with a (formal, infinite-dimensional) Riemannian structure, see also [107, 94] for a comprehensive introduction. In our setting with the $\Theta(x)$ intrinsic tensor, the (formal) gradient of a functional $\mathcal{F}(q)$ with respect to this Riemannian structure reads

[TABLE]

see [78]. Here $\frac{\delta\mathcal{F}}{\delta q}$ denotes the first variation computed in the usual Euclidean sense, e.g. if $\mathcal{F}(q)=\int E(q(x))\mathrm{d}x+\int V(x)q(x)\mathrm{d}x$ then $\frac{\delta\mathcal{F}}{\delta q}(x)=E^{\prime}(q(x))+V(x)$. This is the exact counterpart of the discrete formula (21), with the subtle difference that the geometry on the underlying space $\mathcal{X}=[0,1]$ is now encoded by the mobility $\Theta(x)$, while the “discrete geometry” on $\mathcal{X}_{N}=\{i/N,\,\,i=0\dots N\}$ was previously encoded directly by the kernel $\mathbf{K}$ in (21).

4.4 Gradient flow formulation

We want to identify now the Kimura Equation (35) as a gradient flow, based on the formula (54). To this end we first need to retrieve the evolution equation for the rescaled $Q$-process $q(t,x)$.

From [23] we know that the absolutely continuous part $\widetilde{p}(t,x)$ satisfies $\partial_{t}\widetilde{p}=\mathcal{L}\widetilde{p}$ in the classical sense and remains smooth up to the boundary. Since $\mathcal{L}^{\dagger}w=-\lambda w$ with zero boundary values we see that

[TABLE]

hence the weighted mass decays exponentially

[TABLE]

as in the discrete counterpart (14). Discarding the case of completely absorbed initial data (leading to a trivial stationary evolution $p(t)\equiv p^{\mathrm{I}}=a^{\mathrm{I}}\delta_{0}+0+b^{\mathrm{I}}\delta_{1}$), we can assume that the transient dynamics never gets absorbed, $\langle w,\widetilde{p}(t)\rangle\neq 0$ for all $t>0$, and thus renormalise

[TABLE]

Note that $q(t,\cdot)$ is a probability (density) by construction. As in the discrete case, $q(t,\cdot)$ can be obtained as the law of the natural $Q$-process, i.e. the original stochastic process $X_{t}$ conditioned to non-extinction in infinite time – see for example [80, 25]. Since $\widetilde{p}(t,x)$ is uniformly bounded, cf. [23], and the eigenfunction $w(x)$ satisfies Dirichlet boundary conditions, we see that $q$ automatically satisfies $q(t,0)=q(t,1)=0$ on the boundary, see Remark 11.

Remark 11**.**

In optimal transport and Wasserstein gradient flows, the evolution takes place by construction in the space of probabilities $\mathcal{P}([0,1])$, and one therefore usually enforces no-flux boundary conditions in the PDEs so as to comply with the conservation of mass—see, however, [47] for an application to gradient flows with Dirichlet boundary conditions. However in our framework, since $\widetilde{p}$ is bounded and $w$ vanishes on the boundaries, our new variable $q=\frac{w\widetilde{p}}{\langle w,\widetilde{p}\rangle}$ should also vanish and the evolution (57) is implicitly understood here with Dirichlet boundary conditions $q|_{\partial\Omega}=0$. This is of course not a contradiction: since $w$ and $\Theta$ vanish linearly we see from (44) that $\pi=wz$ does too, and the effective flux in (57) is

[TABLE]

on the boundaries $x=0,1$. Thus the usual no-flux condition is here equivalent to our Dirichlet condition for $q$.

Let us now identify the evolution law for our new variable $q$. From [19] we know that $\widetilde{p}$ satisfies $\partial_{t}\widetilde{p}=\mathcal{L}\widetilde{p}$ in the classical sense, whence

[TABLE]

Substituting the explicit expressions (37),(38) of $\mathcal{L},\mathcal{L}^{\dagger}$, respectively, we find after a straightforward calculation

[TABLE]

Using Equation (56) and Lemma 9 to identify $\pi=\frac{1}{C}\frac{w^{2}e^{-2V/\kappa}}{\Theta}$ inside the logarithm, we get

[TABLE]

Computing the first variation $\frac{\delta\mathcal{H}}{\delta q}=\frac{\kappa}{2}\left(\log\left(\frac{q}{\pi}\right)-1\right)$ of the entropy (46) and applying formula (54) for the Wasserstein gradient, we finally recognise the gradient flow structure

[TABLE]

for the Kimura Equation.

Remark 12**.**

For any solution $q$ of Equation (57), it is not difficult to check that $\mathcal{G}_{\phi}(q|\pi)$ given by Equation (47) is nonincreasing in time as in Lemma 5 and Proposition 1. In fact we can even compute the dissipation

[TABLE]

because $\phi$ is convex and because $q$ satisfies homogeneous Dirichlet boundary conditions at the end points — see Remark 11.

5 The replicator dynamics

The replicator dynamics was introduced in [105] and termed so in [96]. The model consists in an infinite population of $n$ possible types, with frequency $x_{i}$ of type $i$. This dynamics is based upon a simple postulate: the per-capita growth rate $\left(\frac{\dot{x}_{i}}{x_{i}}\right)$ of $x_{i}$ is given by the difference between the expected fitness of type $i$ and the population average fitness, i.e.

[TABLE]

where $\psi_{i}$ is the fitness of type $i$ when the populations is at state $\mathbf{x}$, and $\bar{\psi}(\mathbf{x})=\sum_{i=1}^{n}x_{i}\psi_{i}(\mathbf{x})$. More recently, [61] popularised the Replicator dynamics under two-player games, i.e., with $\psi_{i}(\mathbf{x})=(\mathbf{A}\mathbf{x})_{i}$ for a given $n\times n$ matrix $\mathbf{A}$ — typically $\mathbf{A}$ is associated to payoffs of a $n$-strategy, two-player game.

Equation (58) is a cornerstone of evolutionary game theory, and it has been discussed and reviewed in various works [62, 109, 101, 61]. Incidentally, as we shall review below, it is also associated with the vanishing viscosity limit $\kappa\to 0$ of the Kimura Equation (35). It is worth noticing that, while Kimura Equation arises in the infinite population limit as a delicate balance between selection effects and genetic drift, deviations from this balance lead to either a pure diffusive model or to a hyperbolic one — the latter arises from selection dominating the genetic drift in the large population limit, and as discussed in [21] it is equivalent to the Replicator Equation. We refer to [19, 21] for a discussion about the different scalings and corresponding limits, see also [22] for a discussion on the different regimes both in finite and infinite populations.

In what follows, we will consider the case of $n=2$ types only, and in this case we may write $x_{1}=X$ and $x_{2}=1-X$ and write the generalised one-dimensional Replicator Equation

[TABLE]

with same coefficient $\Theta(x)$ and potential $V(x)$ as in Section 4. Through natural embedding of point sets into empirical probability measures, any $L$-tuple of solutions $X_{1}(t)\dots X_{L}(t)$ to the ODE (59) immediately gives a (probability) measure-valued solution

[TABLE]

to the corresponding hyperbolic PDE

[TABLE]

which we also call the Replicator Equation with a slight abuse of notation. The characteristics ODE is the Lagrangian counterpart of the Eulerian hyperbolic PDE. Note that (60) is obtained formally by taking the diffusion $\kappa=0$ in the Kimura Equation (35).

Since (60) is written in terms of the $p$ variable, one might wonder as in Section 2 what the conditioned $Q$-process might be, and what the resulting dynamics would be. The main difference is that, since $\kappa=0$ here, no random fluctuation arises and the process is purely deterministic, $\mathrm{d}X_{t}=-\Theta(X_{t})\partial_{x}V(X_{t})\mathrm{d}t$. Therefore by Cauchy-Lipschitz uniqueness of trajectories for the ODE (59), any Lagrangian particle $X(t)$ initially on the boundaries remains absorbed, while a particle starting from the interior cannot reach the boundaries in finite time. For the PDE this implies that the mass of $p^{\mathrm{I}}=a^{\mathrm{I}}\delta_{0}+\widetilde{p}^{\mathrm{I}}(x)\mathrm{d}x+b^{\mathrm{I}}\delta_{1}$ initially on the boundaries remains absorbed, while only the transient mass can evolve in the interior: In other words the distribution remains of the form $p(t)=a^{\mathrm{I}}\delta_{0}+\widetilde{p}(t,x)\mathrm{d}x+b^{\mathrm{I}}\delta_{1}$ for $t>0$, with $\|\widetilde{p}(t)\|_{L^{1}}=\|\widetilde{p}^{I}\|_{L^{1}}$. As a consequence absorption never occurs, mass no longer leaks from the interior to the boundaries, and the previous transient rescaling from Section 4 now simply reads

[TABLE]

Up to the constant-in time scaling factor $\frac{1}{\|\widetilde{p}(t)\|_{L^{1}}}=\frac{1}{\|\widetilde{p}^{\mathrm{I}}\|_{L^{1}}}$ we have thus $q=\widetilde{p}$ for the replicator dynamics, and in fact one should think of the replicator Equation (60) as acting on the $q$ variable rather than on $p$.

That being said, we have two equivalent gradient flow formulations for the replicator dynamics:

It is well known that the Replicator ODE (59) is a gradient flow with respect to the (generalised) Shahshahani metric, see [7, 1, 2]. Indeed, choosing again to view $\Omega=(0,1)$ as a Riemannian manifold with the scalar product induced by $\frac{1}{\Theta}$ on $T_{x}\Omega$ (see Subsection 4.3), an immediate computation allows to obtain the intrinsic gradients as $\operatorname{grad}_{\mathsf{d}}=\Theta(x)\partial_{x}$, whence

[TABLE] 2. 2.

For the Replicator PDE (60), the Eulerian energy corresponding to the previous Lagrangian $V(X)$ for single particles is naturally

[TABLE]

Our previous formula (54) for Wasserstein gradients gives next, with the first variation $\frac{\delta\mathcal{V}}{\delta q}=V$, the gradient flow structure

[TABLE]

The convergence of the Kimura Equation (35) towards the Replicator Equation (60) in the deterministic limit $\kappa\to 0$ is well-known from a classical PDE point of view [19, 21], but we will show in the next section that, using the right $q$ variable dictated by the conditioning of the corresponding $Q$-process, the convergence is variational (in some precise sense to be discussed later). This is why we carefully and intentionally wrote the gradient flow (61) in terms of $q$ instead of using the original $p$ variable. As just discussed, this is completely equivalent for the Replicator Equation (up to multiplicative scaling $p=q$), but the situation is drastically different in the presence of diffusion.

6 Variational structures and their compatibility

So far, we have discussed three different gradient flow structures for models that are relevant to evolutionary biology:

1. (i)

the finite population dynamics discussed in Section 2, defined for population size $N<\infty$ 2. (ii)

the continuous population counterpart $x\in[0,1]$ discussed in Section 4, defined for $N=\infty$ and diffusion parameter $\kappa>0$ 3. (iii)

the replicator dynamics discussed in Section 5, which in its Eulerian formulation (60) is defined for $N=\infty$ and $\kappa=0$

The convergence of (i) towards (ii) in the large population limit $N\to\infty$ as well as that of (ii) towards (iii) in the deterministic limit $\kappa\to 0$ have been proved to hold in particular models and in some appropriate sense (e.g. weak convergence of measures, or uniform convergence), see [41, 24, 19, 21]. In this section we intend to convinde the reader that, under reasonable assumptions satisfied by many processes used in population genetics, our framework allows to further identify these convergences as natural variational $\Gamma$-convergence of gradient flows, in a sense to be discussed shortly. To the best of our knowledge this was never considered before for our classical triad: Not only do we provide a gradient flow structure for each of the three settings, but our structures are moreover energetically compatible with the relevant limits $N\to\infty,\kappa\to 0$. We should however stress that we do not aim at proving new convergence results here. In addition, writing down a full, rigorous proof for $\Gamma$-convergence of gradient flow is usually a nontrivial task involving significant technical work. Here we will only provide partial results in this direction, and we will be content with the convergence of the driving energy (entropy) functionals and of the metric structures.

Perhaps, one of the surprising features of the discussion in the previous sections was its reliance on the $Q$-process variable $q$ for the analysis, since the latter is related to eternal paths in a system where the dynamics is almost surely absorbed in finite time. Traditionally, when dealing with models where absorption is certain, one relies on quasi-stationary distributions (which happen to exist in most models of interest) in order to understand the fate of trajectories prior to absorption. However, when investigating a possible variational structure one should not expect this approach to be appropriate. Indeed, as already discussed, absorption is a non-reversible process, and reversibility was a key feature in obtaining a variational structure. Thus, for a generic trajectory that has not been absorbed at time $t$, the probability that it will remain non-absorbed decreases exponentially over time. When this probability decreases very slowly, meta-stable sates arise. However, even if these meta-stable states persist for very long times, the dynamics eventually becomes non-reversible in the long run. As a consequence one should not expect these trajectories to have a variational dynamics. These observations suggest that interesting trajectories should then be the immortal ones, i.e. those that never get absorbed. Two remarkable facts then happen: (i) this subset of trajectories is not empty, and one can indeed obtain a variational dynamics for these trajectories; (ii) the knowledge of the dynamics on this very restricted and small (zero-measure, negligible) set of trajectories is sufficient to recover the full transient dynamics (hence the whole dynamics, since the evolution of the absorbed states can be deduced from the transient dynamics with the help of the additional conservation laws). Roughly speaking, this is why one should rather consider the $Q$-process $q$ instead of the original $p$ distribution when seeking for a variational (gradient flow) structure, whether it be at the discrete or continuous level.

6.1 Gamma-convergence of gradient flows

We first discuss shortly the notion of variational convergence of gradient flows needed for our purpose, and follow closely the exposition in [93, 97]. Let us remind that the notion of $\Gamma$-convergence, introduced by E. De Giorgi in the 70’s, is a notion of convergence of functionals that essentially guarantees convergence of the minimisers – see the classical monograph [30] for a detailed introduction. In some sense, this is precisely the notion of convergence needed when handling minimisation problems, and $\Gamma$-convergence is ubiquitous nowadays in variational analysis and modelling.

Often times one deals in practice with sequences of functionals that are not necessarily defined on the same space. In the following one should roughly keep in mind the idea of a $\Gamma$-converging sequence of functionals defined on a sequence of converging spaces. In order to illustrate this general idea, assume for simplicity that we are given a sequence of Hilbert spaces $\{H_{\epsilon}\}_{\epsilon>0}$ and a “limit” Hilbert space $H$, with some “projections” $\Pi_{\epsilon}:H_{\epsilon}\to H$. We say that a sequence $u_{\epsilon}\in H_{\epsilon}$ converges to $u\in H$ as $\epsilon\to 0$, denoted $u_{\epsilon}\overset{S}{\rightharpoonup}u$, if $\Pi_{\epsilon}(u_{\epsilon})\xrightarrow[]{\sigma}u$ in $H$ for some topology $\sigma$. Both the projection $\Pi_{\epsilon}$ and topology $\sigma$ are crucial choices that one should make, depending on the model and applications under consideration. The Gamma-convergence of functionals on varying spaces is then defined as

Definition 6** ($\Gamma$-convergence).**

We say that a sequence of functionals $F_{\epsilon}:H_{\epsilon}\to\mathbb{R}\cup\{+\infty\}$ $\Gamma$-converges to $F:H\to\mathbb{R}\cup\{+\infty\}$ as $\epsilon\to 0$, denoted by $F=\mathop{\operatorname{\Gamma-\lim}}F_{\epsilon}$ or $F_{\epsilon}\overset{\Gamma}{\to}F$, if the following $\mathop{\operatorname{\Gamma-\liminf}}$ and $\mathop{\operatorname{\Gamma-\limsup}}$ conditions hold

1. (i)

for any sequence $u_{\epsilon}\overset{S}{\rightharpoonup}u$ there holds

[TABLE] 2. (ii)

for any $u\in H$ there exists a recovery sequence $u_{\epsilon}\overset{S}{\rightharpoonup}u$ such that

[TABLE]

Consider now the sequence of gradient flows $u_{\epsilon}(t):[0,T]\to H_{\epsilon}$ given by

[TABLE]

where we emphasise the fact that the gradient of $F_{\epsilon}$ is computed with respect to the $H_{\epsilon}$ structure. Then, since gradient flows tend to minimise the energy along the evolution, and because $\Gamma$-convergence guarantees convergence of minimisers towards minimisers, one expects that limits of $F_{\epsilon}$-gradient flows should be gradient flows for the limiting functional $F=\mathop{\operatorname{\Gamma-\lim}}F_{\epsilon}$. (We shall refer to any such convergence as $\Gamma$-convergence of gradient flows.) This was proved in [93, 97] under additional conditions:

Theorem 3** ([97, Theorem 1]).**

Assume that $F=\mathop{\operatorname{\Gamma-\lim}}F_{\epsilon}$ for $\epsilon\to 0$ as in Definition 6, and let $u_{\epsilon}(t)$ be a solution of (64) such that

[TABLE]

for some limiting curve $u:[0,T]\to H$. If, additionally, $F(u(0))=\lim F_{\epsilon}(u_{\epsilon}(0))$, and if the lower bounds on the velocity and slope from [97] hold, then $u$ is a solution of the limit gradient flow

[TABLE]

We deliberately remain formal at this level and refrain from further discussing the precise definition of the above velocity and slope lower bounds, see again [93, 97]. We should however stress that these two additional conditions are not a mere technical detail, and checking their validity is usually the most difficult part when trying to prove completely rigorous $\Gamma$-convergence of gradient flows. Since our main concern here is the new biological paradigm rather than the careful and rigorous mathematical analysis, we deliberately choose to omit this technical part in order not to obfuscate the exposition. However, we do not want to convey the wrong idea that $\Gamma$-convergence of the driving functional $F_{\epsilon}\xrightarrow[]{\Gamma}F$ and of the Hilbert spaces $H_{\epsilon}\to H$ are sufficient for the convergence of the associated gradient flows: There are of course counterexamples, but we believe that both the speed and slope conditions should hold in practice for our models in the large population and deterministic limits, $N\to\infty$ and $\kappa\to 0$.

Obtaining the convergence (65) of the sequence $u_{\epsilon}(t)$ towards some limit curve $u(t)$ is in general not involved (standard weak compactness arguments typically apply), the challenge is rather to conclude that this limit is in fact a gradient flow for the limiting functional. (In such nonlinear settings this usually requires strong convergence.)

One can actually build a theory of gradient flows in mere metric spaces (thus dispensing from any Hilbert or differential structures), as originally formulated by De Giorgi [32] in terms of curves of maximal slope. This is the notion we shall implicitly refer to in the sequel when we speak of variational evolution or metric gradient flow. For the sake of exposition we refrain from discussing this delicate definition and refer instead to the classical monograph [6] (see also Subsection 6.4 below). It was observed in [93] that the above scenario of $\Gamma$-convergence of gradient flows should hold in this very general metric setting, namely: if $(\mathcal{X}_{\epsilon},d_{\epsilon})$ is a sequence of metric spaces “converging” to a limit metric space $(\mathcal{X},d)$ – e.g. in the Gromov-Hausdorff sense – then the $\Gamma$-convergence of the driving functionals $F_{\epsilon}:\mathcal{X}_{\epsilon}\to\mathbb{R}$ towards $F:\mathcal{X}\to\mathbb{R}$ should reasonably suffice (under additional speed and slope conditions) to guarantee the convergence of the corresponding metric gradient flows. Again, we should stress that this is not completely rigorous: the additional speed and slope lower bounds from [97], actually required for a full convergence, correspond to suitable $\Gamma-\liminf$ estimates for the metric structure and metric derivative of the driving functional. The former can sometimes be related to Gromov-Hausdorff convergence of metric spaces through Benamou-Brenier formulations, and both together provide suitable energy dissipation in the limit [87]. Although the $\Gamma$-convergence of the functional and the Gromov-Hausdorff convergence of the underlying metric spaces alone do not suffice in general for the $\Gamma$-convergence of the gradient flows [87, 88], we claim that both the additional speed an slope convergences hold for our three specific evolutionary models, but we will not push further the rigorous mathematical analysis.

In the next section we shall exemplify this general scenario in two particular cases:

In the limit of large populations $N\to\infty$ the space of discrete probability measures $\mathcal{X}_{N}=\mathcal{P}(\Delta_{N+1})$, endowed with the discrete Wasserstein distance $\mathcal{W}_{N}$ from Definition 4, will Gromov-Hausdorff converge to the continuous Wasserstein space $\mathcal{X}=\mathcal{P}([0,1])$, endowed with the continuous Wasserstein-Shahshahani distance $\mathcal{W}$ from Section 4. For the driving functionals we shall consider the sequence of discrete relative BGS entropies $H(\mathbf{q}|\bm{\pi}_{N})$ that $\Gamma$-converge to the continuous counterpart $\mathcal{H}(q|\pi)$ (up to scaling factors). As a result and loosely speaking, the Moran process will $\Gamma$-converge to the Kimura model; more precisely, Equation (22) will $\Gamma$-converge to Equation (57). 2. 2.

In the deterministic limit of small diffusion $\kappa\to 0$, the metric space will be fixed to be the Wasserstein space $\mathcal{P}([0,1])$, endowed with the fixed Wasserstein-Shahshahani distance $\mathcal{W}$, and we will consider the sequence of functionals

[TABLE]

Here $V(x)$ is the same potential initially prescribed for the Kimura Equation (35), and the effective potential $V_{\kappa}=-\frac{\kappa}{2}\log\pi_{\kappa}$ (to be defined below in more details) will converge to $V$. As a result the Kimura gradient flow (57) will $\Gamma$-converge to the Replicator gradient flow (61).

In both cases the $\Gamma$-convergence will be taken relatively to the weak-$*$ convergence of measures $q_{n}\overset{*}{\rightharpoonup}q$, which we recall is defined by duality with bounded, continuous test-functions as $\int\varphi(x)\mathrm{d}q_{n}(x)\to\int\varphi(x)\mathrm{d}q(x)$ for all $\varphi\in\mathcal{C}_{b}$. This is a very reasonable choice, because the Wasserstein distance metrises the weak-$*$ convergence [108].

For each case we will try to justify below why the above scenario should hold, namely we discuss the relevant $\Gamma$-convergences and the Gromov-Hausdorff convergence of the metric spaces. However, we will neither address the speed and slope lower bounds, nor discuss the well-preparedness of the initial data in Theorem 3. We believe that those important assumptions hold true at least for the examples in Figure 1, but, as pointed out by Sandier and Serfaty in [93], this is a case-to-case issue to be addressed by hand based on the specific structure of the problem under consideration. We also refer to [37] for the rigorous derivation of such lower bounds in the particular setting of finite-volume discretization of the Fokker-Planck Equation. For most of the classical processes, the convergence of solutions was proved most of the time in very strong topologies in previous works, using various techniques essentially based on PDE methods [41, 24, 19, 21]. Again, our interest does not lie here in rigorous proofs of convergence, and we rather wish to illustrate the energetic compatibility between the gradient flow structures.

6.2 The large population limit $N\to\infty$

With the particular Moran process and Kimura Equation in mind, we restrict here the general statements from Section 2 to $k=2$ absorbing states. We chose to do so mainly for the ease of exposition, but the discussion below however extends to more general situations.

In this particular setup the original $\mathbf{p}$ variable is thus $(N+1)$ dimensional, and the rescaled $Q$-process $\mathbf{q}$ is $(N-1)$-dimensional. The $(N+1)$-dimensional kernel $\mathbf{M}$ gives rise to the rescaled $(N-1)$-dimensional kernel $\mathbf{K}$, and we often write $\mathbf{L}:=\mathbf{K}^{\dagger}-\mathbf{I}$ for the effective kernel driving the transient evolution $\frac{\mathrm{d}\,}{\mathrm{d}t}\mathbf{q}=(\mathbf{K}^{\dagger}-\mathbf{I})\mathbf{q}=\mathbf{L}\mathbf{q}$ from Lemma 4.

Writing again $x_{i}=\frac{i}{N}$ for the uniform partition of $[0,1]$, we map canonically any $(N-1)$ probability vector $\mathbf{q}$ to the associated empiric measure

[TABLE]

(and in fact $\widehat{q}\in\mathcal{P}(0,1)$ since the absorbing states $i=0,N$ were discarded in the construction of the $Q$-process). Note that $\mathbf{q}\mapsto\widehat{q}$ actually defines the “projections $\Pi_{\epsilon}$” from Subsection 6.1, allowing to embed a sequence of varying discrete (probability) spaces into the limit (probability) space $\mathcal{P}([0,1])$. Slightly abusing the notations, we say that a sequence $\mathbf{q}_{N}$ of probability vectors converges weakly-$*$ to the probability measure $q\in\mathcal{P}([0,1])$ as $N\to\infty$, denoted $\mathbf{q}_{N}\overset{*}{\rightharpoonup}q$, if $\widehat{q}_{N}\overset{*}{\rightharpoonup}q$ weakly-$*$ in the sense of measures. Note that this defines the abstract convergence “$\overset{S}{\rightharpoonup}$” from Subsection 6.1.

We have then

Lemma 11**.**

Assume that $\bm{\pi}_{N}\overset{*}{\rightharpoonup}\pi$, let $H(\mathbf{q}|\bm{\pi}_{N})$ be given by Equation (16) if $\mathbf{q}\ll\bm{\pi}_{N}$ and $+\infty$ otherwise, and $\mathcal{H}(q|\pi)$ be given by Equation (46). Then

[TABLE]

in the sense of Definition 6.

This only requires convergence of the reference measures $\bm{\pi}_{N}\overset{*}{\rightharpoonup}\pi$, which should be satisfied in practice for reasonably physical models. A typical situation that we shall consider below is when the effective kernel $\mathbf{L}_{N}=\mathbf{K}_{N}^{\dagger}-\mathbf{I}$ converges to the continuous operator $\mathcal{L}$, in which case one should also expect the eigenvectors $\mathbf{w}_{N},\mathbf{z}_{N}$ to converge to the continuous eigenfunctions $w(x),z(x)$ and thus the measures $\bm{\pi}_{N}=\mathbf{w}_{N}\circ\mathbf{z}_{N}\rightharpoonup w(x)z(x)\mathrm{d}x=\pi(x)\mathrm{d}x=\pi$ at least in some weak sense. One can check that this holds at least for the convergence of the neutral Moran process towards the Kimura Equation, see Subsection 6.5 below.

Proof.

We need to check the two conditions in Definition 6, and we begin with the $\mathop{\operatorname{\Gamma-\liminf}}$ part. Let $\mathbf{q}_{N}\overset{*}{\rightharpoonup}q$ be an arbitrary converging sequence, and notice that by definition we have $H(\mathbf{q}_{N}|\bm{\pi}_{N})=\mathcal{H}(\widehat{q}_{N}|\widehat{\pi}_{N})$ with $\widehat{q}_{N}\overset{*}{\rightharpoonup}q$ and $\widehat{\pi}_{N}\overset{*}{\rightharpoonup}\pi$ in the sense of measures. From the convexity and lower semi-continuity of $\eta\mapsto\eta\log\eta$ we can immediately apply [4, Theorem 2.34] to conclude that the $\mathop{\operatorname{\Gamma-\liminf}}$ condition (i) holds as

[TABLE]

For the $\mathop{\operatorname{\Gamma-\limsup}}$ part (ii), fix any $q\in\mathcal{P}([0,1])$. If $\mathcal{H}(q|\pi)=+\infty$ there is nothing to prove, hence we can assume that $\mathcal{H}(q|\pi)<+\infty$ and in particular $q\ll\pi$ and $f:=\frac{\mathrm{d}q}{\mathrm{d}\pi}\in L\log L(\mathrm{d}\pi)\subset L^{1}(\mathrm{d}\pi)$. By approximation it is enough to consider $f\in\mathcal{C}_{b}([0,1])$ positive, in particular $f\log f\in\mathcal{C}_{b}([0,1])$. Let $f_{N}(x):=\frac{1}{\sum_{j=1}^{N}f(x_{j})\pi_{j}}f(x)$ be the renormalisation of $f$, and observe that $f_{N}\to f$ uniformly because $\sum_{j}f(x_{j})\pi_{j}=\int f(x)\mathrm{d}\widehat{\pi}_{N}\to\int f(x)\mathrm{d}\pi(x)=\int\mathrm{d}q(x)=1$ (our main assumption is precisely that $\widehat{\pi}_{N}\overset{*}{\rightharpoonup}\pi$). Defining next the discrete probability vector $\mathbf{q}_{N}$ by $q_{i}:=f_{N}(x_{i})\pi_{i}$, we have of course $\mathbf{q}_{N}\ll\bm{\pi}_{N}$, $\frac{\mathrm{d}\mathbf{q}_{N}}{\mathrm{d}\bm{\pi}_{N}}=f_{N}$ and $\mathbf{q}_{N}\overset{*}{\rightharpoonup}f\pi=q$. Moreover

[TABLE]

as desired, where the convergence follows from our initial assumption that $\bm{\pi}_{N}\overset{*}{\rightharpoonup}\pi$ (i-e $\widehat{\pi_{N}}\overset{*}{\rightharpoonup}\pi$) and the strong uniform convergence $f_{N}\log f_{N}\to f\log f$ (the function $\eta\mapsto\eta\log\eta$ being uniformly continuous in any bounded interval $\eta\in[0,M]$). ∎∎

Turning now to the convergence of the metric spaces, N. Gigli and J. Maas proved in [51] that, if $\mathbf{L}_{N}=\mathbf{K}_{N}^{\dagger}-\mathbf{I}$ is the standard Laplacian matrix in the discrete $d$-dimensional torus $\mathbb{T}_{N}^{d}=(\mathbb{Z}/N\mathbb{Z})^{d}$ with mesh size $1/N$ (i-e if $\mathbf{K}_{N}$ is the uniform random walk), then the discrete Wasserstein space $(\mathcal{P}(\mathbb{T}_{N}^{d}),\mathcal{W}_{N})$ converges to the continuous Wasserstein space $(\mathcal{P}(\mathbb{T}^{d}),\mathcal{W})$ in the sense of Gromov-Hausdorff. This is a non trivial result per se, whose proof is the whole purpose of [51]. Let us mention that Gigli and Maas already pointed out potential applications to $\Gamma$-convergence of gradient flows, see also [3, 37] for a description of finite volume schemes for Fokker-Planck Equations as discrete Wasserstein gradient flows and their convergence towards continuous counterparts [66].

Here we do not pretend to prove any rigorous statement in this direction, and we shall be content with the following heuristics:

Claim 1**.**

Let $\Theta(x)$ satisfy our assumptions (36), and $U\in\mathcal{C}^{2}(0,1)$. If the effective kernel $\mathbf{L}_{N}=\mathbf{K}_{N}^{\dagger}-\mathbf{I}$ arises from any reasonable finite difference discretization of the operator

[TABLE]

on the domain $\Omega=(0,1)$ with no-flux boundary conditions and uniform mesh size $1/N$, then the discrete Wasserstein space $(\mathcal{P}_{N},\mathcal{W}_{N})$ Gromov-Hausdorff converges, as $N\to\infty$, to the continuous Wasserstein-Shahshahani space $(\mathcal{P},\mathcal{W})$.

Remark 13**.**

Similar questions were investigated in [52, 53], where it was shown that some homogeneity and uniformity of the space meshing is essential (as assumed here). Note that part of our statement in Claim 1 is that the limiting distance $\mathcal{W}$ does not “see” the potential $U$, while the discrete Wasserstein distance does depend on the whole kernel $\mathbf{K}_{N}$ (hence a priori on $U$). This means that, as the population size increases, the discrete distance only retains the purely diffusive part $\operatorname{div}(\Theta\nabla q)$ of the elliptic operator, while the influence of the drift $\operatorname{div}(\Theta q\nabla U)$ smears out and vanishes in the limit. We believe that this should be the case in higher dimensions as well, see again [51] for a particular and rigorous $d$-dimensional statement in the torus. Let us point out that, in general, the discrete distance $\mathcal{W}_{N}$ is highly non local. Our assumption that the kernel arises from a finite difference approximation essentially localises the distance as $N\to\infty$: for example if the discretization is given by a fixed $m$-points stencil, then each discrete state interacts with a neighbouring fraction $m/N\to 0$ of all possible states (the kernel $\mathbf{K}_{N}$ becomes sparser and sparser, and asymptotically concentrates on the diagonal). This dependency of the discrete metric on the potential can be seen as a drawback, and raises the question whether our processes possess more appropriate gradient structures, independently of the potential. For the sake of exposition we prefer keeping this question open for future work, and refer to [88] for recent results and related concepts of tilt and contact convergence of gradient flows.

Remark 14**.**

Here we slightly abuse the notations: the operator $\mathcal{L}$ in our claim is not the forward operator (37) acting on $p$ as before, but rather the divergence-form operator (57) acting on $q$, namely $\mathcal{L}q=\partial_{x}\left(\Theta q\partial_{x}\left(\log\left(\frac{q}{\pi}\right)\right)\right)=\partial_{x}(\Theta\nabla q)-\partial_{x}(\Theta q\partial_{x}\log\pi)$ with potential $U=-\log\pi$ and $\kappa=2$.

To give a hint of why the claim should hold, we remember Lemma 6 and state that

Corollary 1**.**

Let $\mathbf{L}_{N}=\mathbf{K}_{N}^{\dagger}-\mathbf{I}$ be the three-points stencil, forward finite difference discretization of $\mathcal{L}q=\partial_{x}(\Theta\partial_{x}q)+\partial_{x}(\Theta q\partial_{x}U)$ as in our Claim 1. For fixed $x\in(0,1)$ choose $i=i_{N}\sim\lfloor Nx\rfloor$ such that $x_{i}=\frac{i}{N}\to x$. Then

[TABLE]

This means that the discrete Wasserstein distance fully encodes, at least asymptotically, the local geometry of the continuous Shahshahani space $(\Omega,\mathsf{d})$ – i.e. the Riemannian tensor $\frac{1}{\Theta(x)}$ on $T_{x}\Omega$. Indeed formula (66) gives the $\frac{1}{N}$-infinitesimal variations of the discrete distance between point-measures located at $x_{i}\sim x$ and $x_{i+1}\sim x+1/N$. This is the exact counterpart of the continuous Wasserstein distance between point-masses

[TABLE]

which is easily checked using the representation formula (50) for the Shahshahani distance. Since the overlying Wasserstein space $(\mathcal{P}([0,1],\mathcal{W}))$ is built upon – and reflects the geometry of – the underlying Shahshahani space $(\Omega,\mathsf{d})$, this explains now our claim 1.

Proof of Corollary 1.

In order to keep the notations light we write $x_{i}=i/N$, $\Delta x=x_{i+1}-x_{i}=\frac{1}{N}$, $\Theta_{i}=\Theta(x_{i})$, and $U_{i}=U(x_{i})$. Our assumption that $\mathbf{L}$ is the forward finite difference approximation of $\mathcal{L}q=\partial_{x}(\Theta\partial_{x}q)+\partial_{x}(\Theta q\partial_{x}U)$ means here that

[TABLE]

The off-diagonal coefficients of $\mathbf{K}=\mathbf{L}-\mathbf{I}$ thus read

[TABLE]

and

[TABLE]

where the leading $\frac{1}{\Delta x^{2}}$ order only stems from the higher-order diffusive part of the operator $\mathcal{L}$. The lower $\mathcal{O}(1/\Delta x)$ order arises from the drift part exclusively, after absorbing one order of $1/\Delta x$ into $\left|\frac{U_{i+1}-U_{i}}{\Delta x}\right|\approx|\partial_{x}U|\lesssim C$ uniformly in $\Delta x$. This explains why the limit distance $\mathcal{W}$ does not see the potential $U$ at leading order, and the same argument would carry through with any consistent discretization (centered, backwards, five-point stencil…) Appealing next to Lemma 6, substituting the above expressions in (29), and by scaling properties $\beta(\lambda x,\lambda y)=\lambda\beta(x,y)$ of the logarithmic mean, we finally get

[TABLE]

and the proof is complete. ∎∎

It is worth noticing that all the computations in this proof are locally uniform in $i$ as long as $x_{i}=i/N$ remains bounded away from the boundaries, which in our statement was guaranteed since we considered $x_{i}\to x\in(0,1)$. Of course the behaviour close to the boundaries is drastically different because $\Theta(x)$ vanishes and the diffusion degenerates.

6.3 The deterministic limit $\kappa\to 0$

Here we aim at recovering the Replicator Equation (60) from the Kimura Equation (35). For both dynamics the metric space is fixed once and for all to be the Wasserstein-Shahshahani space $(\mathcal{P}([0,1]),\mathcal{W})$, hence the only delicate point is to check the $\Gamma$-convergence for the sequence of functionals from Section 4 when $\kappa\to 0$. Emphasising now the dependence on $\kappa$, our relative entropy (46) was

[TABLE]

Let us recall that we implicitly set $\mathcal{F}_{\kappa}(q)=+\infty$ whenever $q$ is not absolutely continuous with respect to $\pi_{\kappa}=\pi_{\kappa}(x)\mathrm{d}x$ (or equivalently with respect to $\mathrm{d}x$), and from Lemma 9 the reference measure can be written

[TABLE]

Here $C_{\kappa}$ is a normalising constant such that $\int\pi_{\kappa}=1$, and $w_{\kappa},z_{\kappa}$ are the principal eigenfunctions of the operators $\mathcal{L}^{\dagger}_{\kappa},\mathcal{L}_{\kappa}$ from Definition 5. Using the above expression of $\pi_{\kappa}$ in terms of the fixed potential $V(x)$ in order to expand $\log\pi_{\kappa}$, we rewrite the entropy in the more convenient form

[TABLE]

As anticipated, we have next

Proposition 3**.**

Assume that the sequence of probability measures $\{\pi_{\kappa}\}_{\kappa>0}$ satisfies a Large Deviation Principle (LDP) with speed $2\kappa^{-1}$ and rate function $V(x)$ in the deterministic limit $\kappa\to 0$, i.e.

[TABLE]

for all Borel sets $E\subset[0,1]$, where $\mathring{E}$ and $\overline{E}$ are the interior and the closure of $E$, respectively. For $q\in\mathcal{P}([0,1])$, define

[TABLE]

Then

[TABLE]

We refer to [35] for an introduction to Large Deviations. We stress that the problem strongly resembles the so-called entropic regularisation in optimal transport, see e.g. [14, 77, 29].

Although not surprising at least formally in view of the expansion (67), this convergence is not trivial in our particular setting: For any fixed $\kappa>0$ the eigenfunction $w_{\kappa}(x)$ vanishes at the boundary, hence the effective potential blows-up as $\lim_{x\to\{0,1\}}V_{\kappa}(x)=+\infty$. However, since the initial potential $V(x)$ is smooth up to the boundaries we see that the convergence $V_{\kappa}\to V$ cannot be uniform in $[0,1]$: the degenerate diffusion induces a thin boundary layer, whose length-scale should be accounted for by the $\Gamma$-convergence. Moreover, and even less evidently, there can exist internal transition layers near interior local maximum points of $V$ leading to metastable states.

Proof.

With our strong assumption (68) this is an immediate consequence of the abstract results in [81, §3]. ∎∎

In practice, one should check by hand the strong hypothesis (68) in each case of interest, and this is a nontrivial task that strongly depends on the structure of the potential $V$. Let us point out that this trivially holds at least for the neutral Kimura Equation, i.e., Equation (5) with constant potential $V$. In this case $\Theta(x)=x(1-x)$ and $z_{\kappa}(x)=1,w_{\kappa}(x)=6x(1-x)$, hence $\pi_{\kappa}(x)=w_{\kappa}(x)z_{\kappa}(x)=6x(1-x)$ is independent of $\kappa$ and trivially satisfies the LDP (68).

6.4 Minimising movements and JKO schemes

In Section 3, we briefly discussed some popular models in finite population evolutionary dynamics, namely the Moran and the Wright-Fisher processes. In these models, as time is a discrete variable, no gradient flow formulation is possible. However, the so-called JKO scheme provides a direct view of the functional that is minimised by evolution, even when considering discrete in time models. In fact, one possible way to make sense of gradient flows in mere metric space is De Giorgi’s minimising movement [33], which is roughly speaking an implicit Euler time-stepping. More precisely, given a driving functional $F$ in the abstract metric space $(\mathcal{X},d)$, one wishes to make sense of $\partial_{t}x=-\operatorname{grad}_{d}F(x)$. To this end, choose a small time step $h>0$, and define the minimising scheme by solving recursively

[TABLE]

This is sometimes called a discrete (in time) gradient flow: Indeed, for the particular instance of a Hilbert space $H$, the distance is $d^{2}(x,y)=\|x-y\|^{2}$, and the Euler-Lagrange Equation associated with $\min\limits_{x}\{\|x-x^{n}\|^{2}/2h+F(x)\}$ is nothing but the Euler implicit scheme $\frac{x^{n+1}-x^{n}}{h}=-\operatorname{grad}_{H}F(x^{n+1})$. Defining next the piecewise-constant interpolation in time $x_{h}(t)\mathrel{\mathop{:=}}x^{n+1}$ for $t\in(nh,(n+1)h]$ one should in general expect convergence $x_{h}(t)\to x(t)$ as the time step $h\to 0$, and obtain in the limit a solution to the abstract gradient flow. We refer to [6] for a detailed overview of this metric theory.

When the space under consideration is that of probability measures endowed with the quadratic Wasserstein distance, Jordan, Kinderlehrer and Otto identified the classical Fokker-Planck Equation as the Wasserstein gradient flow of the BGS entropy [66], precisely by proving the convergence of the minimising movement in the limit of small time steps. In this particular context this time discretization is often called the JKO scheme, after the three authors, and this has proved to be a powerful tool in order to prove existence of weak solutions of PDEs and construct numerical approximations.

Let us illustrate this general idea in the simplest possible case here, namely the one-dimensional replicator ODE (59). In this context, the minimising movement reads

[TABLE]

where $\mathsf{d}$ is the Shahshahani distance from (49) and we initialise $X^{0}=X^{I}$ for some given initial condition $X^{I}$. With our smoothness assumptions on the potential $V$, standard convexity arguments show that $X^{n+1}$ is uniquely well-defined, and as expected we have

Theorem 4**.**

Let $X_{h}(t):[0,\infty)\to[0,1]$ denote the piecewise-constant interpolation of the sequence $X^{n}$ defined in (69). Then $X_{h}(t)$ converges uniformly to $X(t)$ in any finite interval $t\in[0,T]$ as $h\to 0$, where $X(t)$ is the unique solution to the Replicator Equation (59) with initial condition $X(0)=X^{I}$.

Proof.

Testing $X=X^{n}$ as a competitor in the variational scheme (69) and summing over $n$, one gets the classical total square distance estimate

[TABLE]

This is a fairly general but crucial property of the abstract minimising movement, see [6]. One should think of this as an $H^{1}$ estimate $\int_{0}^{\infty}|\dot{X}_{h}(t)|^{2}_{X_{h}(t)}\mathrm{d}t\leq C$, where the metric speed $\dot{X}_{h}(t)\in T_{X_{h}(t)}\Omega$ is measured in the intrinsic Shahshahani metrics. This in turn gives equicontinuity in time and thus compactness of $\{X_{h}\}_{h>0}$ in the uniform topology. Up to extraction of a subsequence we can therefore assume that $X_{h}$ converges uniformly to some $X$, and we only have left to prove that this limit curve is the unique solution to the Replicator ODE.

To this end we use the definition (49) of the Shahshahani distance to differentiate $\mathsf{d}^{2}(\cdot,X^{n})$ and write the optimality condition for (69) as

[TABLE]

Exploiting formula (50) and applying the mean value theorem, there exists $X_{*}\in[X(t),X(t+h)]$ such that

[TABLE]

Taking $h\to 0$, and from the uniform convergence $X_{h}(t)\to X(t)$, we find the Replicator Equation (59), together with the initial condition $X_{h}(0)=X^{\mathrm{I}}\Rightarrow X(0)=X^{\mathrm{I}}$. ∎∎

Similarly, all models discussed so far can be obtained by recursive minimisation of the following JKO functionals (denoting the previous steps $\mathbf{q}_{0},q_{0}$ as parameters): The finite state model given by Equation (22) is associated to the minimisation of the functional

[TABLE]

The generalised Kimura Equation given by (57) is obtained via minimisation of the functional

[TABLE]

and, finally, the hyperbolic Replicator PDE (60) corresponds to

[TABLE]

From our previous discussions on the various $\Gamma$-limits and Gromov-Hausdorff convergence, the one-step operators converge accordingly: on the one hand $J_{N}(\cdot\,;\mathbf{q}_{0}^{N})\xrightarrow[]{\Gamma}\mathcal{J}_{\kappa}(\cdot\,;q_{0})$ in the large population limit $N\to\infty$ (presuming that the previous step $\mathbf{q}^{N}_{0}\overset{*}{\rightharpoonup}q_{0}$), and on the other hand $\mathcal{J}_{\kappa}(\cdot;q_{0})\xrightarrow[]{\Gamma}\mathcal{J}(\cdot;q_{0})$ in the deterministic limit $\kappa\to 0$. Since $\Gamma$-convergence guarantees convergence of the minimisers towards minimisers of the limit functional, this means roughly speaking that the discrete-in-time gradient flows (JKO schemes) also converge.

We finish this subsection with some words on the relevance of JKO schemes in biological modelling. On one hand, there is an old argument if there is a functional which is minimised by biological evolution, cf. e.g. [9, 42]; on the other hand, all models discussed here were introduced without any explicit reference to such functionals. We just showed that not only do such functionals appear in all the models discussed so far, but also that they naturally split into a “free energy” – which nature tries to minimise in time – and an inertial term making such changes more difficult. These ideas will be further explored in a subsequent work.

6.5 From Moran processes to the Kimura Equation

From the discussion in Subsection 1.5, it is clear that the Kimura Equation (5) can be derived from the Moran process 1. In this section, we will reproduce this known result (see [19]) with a different approach. In fact, we will show that the precise form of the degenerate diffusion in Equation (5) is a consequence of the limiting behaviour of the characteristic triple, both in the discrete and continuous cases. More precisely

Lemma 12**.**

Consider the Moran process and assume the weak selection principle (4). Let $(\mu,\widetilde{\mathbf{w}},\widetilde{\mathbf{z}})$ be its characteristic triple. Consider the generalised Kimura Equation (35) and its continuous characteristic triple $(\lambda,w,z)$. Assume further that $\frac{1}{N}\widetilde{\mathbf{w}}\circ\widetilde{\mathbf{z}}^{-1}$ converges pointwise uniformly to $\frac{w}{z}$. Then $\Theta(x)=x(1-x)$.

Proof.

From Equation (28) we first have

[TABLE]

Following [80] we get next

[TABLE]

On the other hand, from Lemma 9, it is clear that $\frac{w(x)}{z(x)}=C\Theta(x)\mathrm{e}^{2V(x)/\kappa}$, for a possibly different constant $C$, and imposing $x=i/N$ we finish the proof. ∎∎

Remark 15**.**

The convergence assumed in Lemma 12 is usually not straightforward to verify, but it is easily checked in the neutral case.

As a direct consequence of Lemmas 10 and 12, we recover that the distance between two populations at deterministic states is given by the so-called Shahshahani distance

[TABLE]

This distance turns out to be the same as the “genetic distance” introduced by Edwards and Cavalli-Sforza in the 1970’s [17, 38]— see also [2, 7].

Remark 16**.**

We hasten to point out that, while the Gromov-Hausdorff convergence of the metrics $\mathcal{W}_{N}\to\mathcal{W}$ is neither necessary nor sufficient for the $\Gamma$-convergence of gradient flows discussed above. However, the compatibility of such metrics suggests that there are generalisations of the metric displayed in Equation (73) that are both compatible and relevant to discrete and continuous non-deterministic models.

Finally, we point out that in most cases of biological interest one expects that the leading eigenvectors of the discrete process will converge pointwise uniformly to the eigenfuctions of the continuous operator, and this is sufficient to verify the assumptions in Lemma 11. In this situation, we will also have that $\mathcal{G}_{\phi}(p|\pi)=G_{\phi}(\mathbf{p}|\overline{\bm{\pi}})+\mathcal{O}\left(N^{-1}\right)$, provided $\phi$ is sufficiently regular — e.g. of locally bounded variation in $\mathbb{R}^{+}$, cf. [27].

7 Conclusion

This work was born from a crossbreeding between two unrelated research programs: (i) to clarify differences and similarities in the triad of evolutionary models, and hence to understand them in a unified way; (ii) to investigate the existence and relevance of local maximisation principles for evolutionary models and, by extension, in evolutionary biology.

From this viewpoint, the main consequence of this work is to show that local maximisation principles — and, consequently, variational structures — may be formulated for the main models in evolutionary dynamics and are compatible among themselves.

Perhaps not surprisingly, the correct approach to address this question was to understand the gradient flow formulation of all considered processes. This program has been already carry out for the Replicator Dynamics in [1, 98]. For finite populations, the gradient flow formulation introduced by Maas in [79] was the basic tool used to obtain a similar formulation for the continuous in time Moran process. The latter was then used as a starting point to obtain this formulation for the continuous process (Kimura Equation). In addition, following earlier work [21], we embed the replicator dynamics in the hyperbolic equation formally obtained from the Kimura Equation in the vanishing genetic drift limit, and verify that the gradient flow structure found in [98, 1] is preserved as expected.

Once these gradient structures are well established, a very natural question arises next: Since there is a limiting relation among the processes involved, are the gradient structures compatible with one another? A very natural framework to answer this question was for us the use of $\Gamma$-convergence of gradient flows.

The current work can also be seen as part of a long program on the geometrisation of evolution initiated by [1] (see also [7, 98, 2, 63]). The gradient structure for the continuous Moran process extends this geometric approach to stochastic models for finite populations without mutation.

As a byproduct, we introduced several possible free-energies (or entropies) for evolutionary processes. A priori, there is no reason to favour any particular one, and it seems fruitful to understand the dynamics induced by all these entropies in different scenarios (i.e., for several fitness potentials $V$), in particular the classical ones in evolutionary dynamics: dominance, convergence, and coexistence. This will require a detailed study of the entropic dynamics of discrete and reducible Markov chains, not limited to the standard BGS entropy [103].

An important question is how far the results obtained here can be extended. Ideally, there should be at least some relevant classes of more complex models (multi-type, structured) that should be amenable to a similar analysis. An initial attempt towards generalisations appears in [18], with promising results; in particular the study of the fitness potential $V$ seems to provide hints on 3-types dynamics without mutations and in the 2-types dynamics with mutations.

Finally, it is important to point out that the optimisation process addressed here are all local — i.e. they are myopic in economics parlance — and this does not guarantee the existence of a global optimisation process. This framework is compatible with the adaptive dynamics point of view of evolution; see [57, 55, 83, 85, 84, 90, 54, 73]; it is also compatible with Maximum Entropy Production Principles [82, 113, 65]. An extreme illustration of the difference between local and global maximisation in biological dynamics is the so-called evolutionary suicide [56].

Acknowledgements FACCC and AMR were partially supported by FCT/Portugal Strategic Project UID/MAT/00297/2019 (Centro de Matemática e Aplicações, Faculdade de Ciências e Tecnologia, Universidade Nova de Lisboa). FACCC also benefited from an “Investigador FCT” grant. LM was supported by FCT/Portugal projects PTDC/MAT-STA/0975/2014 and PTDC/MAT-STA/28812/2017. MOS was partially supported by CNPq under grants # 309079/2015-2, 310293/2018-9 and by CAPES — Finance Code 001.