remarkRemark hypothesisHypothesis claimClaim assumAssumption dfDefinition
Modeling the process of speciation using a multi-scale frameworkM. K. Brun, E.
Ahmed, J. M. Nordbotten, and N. C. Stenseth ex_supplement

# Modeling the process of speciation using a multi-scale framework including a
posteriori error estimates††thanks: Submitted to the editors DATE. This work
was funded in part by Norwegian Research Council project no. 263149.

Mats K. Brun CEES, Dept. of Biosciences, University of Oslo, NO-0316 Oslo,
Norway. Current address: Expert Analytics, NO-0179 Oslo, Norway (/).
<EMAIL_ADDRESS><EMAIL_ADDRESS>Elyes Ahmed SINTEF, NO-0314 Oslo, Norway ().
<EMAIL_ADDRESS>Jan M. Nordbotten Department of Mathematics, University
of Bergen, NO-5020 Bergen, Norway<EMAIL_ADDRESS>Nils Chr. Stenseth
CEES, Dept. of Biosciences, University of Oslo, NO-0316 Oslo, Norway ().
<EMAIL_ADDRESS>
###### Abstract

This paper concerns the modeling and numerical simulation of the process of
speciation. In particular, given conditions for which one or more speciation
events within an ecosystem occur, our aim is to develop the necessary modeling
and simulation tools. Care is also taken to establish a solid mathematical
foundation on which our modeling framework is built. This is the subject of
the first half of the paper. The second half is devoted to developing a multi-
scale framework for eco-evolutionary modeling, where the relevant scales are
that of species and individual/population, respectively. The species level
model we employ can be considered as an extension of the classical Lotka-
Volterra model, where in addition to the species abundance, the model also
governs the evolution of the species mean traits and species trait
covariances, and in this sense generalizes the purely ecological Lotka-
Volterra model to an eco-evolutionary model. Although the model thus allows
for evolving species, it does not (by construction) allow for the branching of
species, i.e., speciation events. The reason for this is related to that of
separate scales; the unit of species is too coarse to capture the fine-scale
dynamics of a speciation event. Instead, the branching species should be
regarded as a population of individuals moving along a selection of trait axes
(i.e., trait-space). For this, we employ a trait-specific population density
model governing the dynamics of the population density as a function of
evolutionary traits. At this scale there is no a priori definition of species,
but both species and speciation may be defined a posteriori as e.g., local
maxima and saddle points of the population density, respectively. Hence, a
system of interacting species can be described at the species level, while for
branching species a population level description is necessary. Our multi-scale
framework thus consists of coupling the species and population level models
where speciation events are detected in advance and then resolved at the
population scale until the branchin is complete. Moreover, since the
population level model is formulated as a PDE, we first establish the well-
posedness in the time-discrete setting, and then derive the a posteriori error
estimates which provides a fully computable upper bound on an energy-type
error, including also for the case of general smooth distributions (which will
be useful for the detection of speciation events). Several numerical tests
validate our framework in practice.

###### keywords:

Adaptive dynamics; speciation; multi-scale; a posteriori error estimates.

92-10, 92B05, 92B99

##  1 Introduction

Mathematical models have a long history in ecology and evolutionary biology,
with the most famous example being the Lotka-Volterra model [18], which
describes the interaction of two species (commonly referred to as predator and
prey) at a timescale where individual traits remain constant, i.e., no
evolution. Going beyond strictly ecological models, the interaction between
ecology and evolution has been studied extensively in the so-called “adaptive
dynamics” literature (see e.g., [6, 7, 16, 17], and the references therein).
In particular, the concept known as “adaptive speciation”, i.e., the idea that
a series of small adaptive changes in the traits of individuals over long
enough time leads to a diversity of species, has received considerable
attention [8, 9, 10, 11, 12, 24]. In the current work, we are studying the
process of adaptive speciation at the level of model formulation where the
_species_ is regarded as the fundamental unit of the eco-evolutionary system,
and where the species dynamics is emergent from the dynamics at the
population/individual level.

The process of _speciation_ is a biological phenomenon in which a population
within one species gradually evolves into two (or more) distinct species (for
classical literature on speciation, see e.g., [4, 16, 21, 23]). As such,
speciation is an emergent property of natural selection acting on a population
of individual organisms. If the individual is regarded as the fundamental unit
of the evolutionary process, the point at which the gradual evolutionary
changes has accumulated sufficiently to produce a new species becomes a matter
of definition (e.g., reproductive isolation). On the other hand, as species is
regarded as the fundamental unit, it becomes important to separate evolution
into two categories; _cladogenesis_ , which is the splitting of a parent
species into new distinct child species, and _anagenesis_ , which is the
gradual evolution of a species that continues to exists. Modeling the process
of speciation using mathematical models therefore presents a different set of
challenges whether one takes the individual or the species as the fundamental
unit of the eco-evolutionary system.

The difficulty in modeling the process of speciation using any species
interaction-type model is inherent in the model itself, i.e., it is assumed
a-priori that distinct species can be identified at all times within the
population. Since a speciation event will necessarily imply some ambiguity in
the identification of distinct species (at least for a period of time), a
breakdown in the underlying assumptions of the model is therefore unavoidable.
In mathematical terms; the _coarse-scale_ species level model is by
construction unable to capture the _fine-scale_ population dynamics of a
speciation event. The mathematical challenge in modeling speciation as an
emergent property is therefore related to the connection of these two scales.
Relevant studies concerning this difficulty has covered e.g., evolutionary
branching driven by stochastic mutations [28], derivation of an eco-
evolutionary model at the species level but not including the variability in
species abundance [5], and the derivation of a model describing the
interaction of different morphs within the same species [25].

A fundamental observation of real biological systems is that most of the time
individuals are clearly grouped into distinct species. Thus, it is natural to
take the species as the fundamental unit of any eco-evolutionary model.
However, another observation is that over evolutionary time speciation events
do occur. Hence, from the time there is only the parent species to the time
when increasingly diverging traits among the constituent individuals has
resulted in new child species, the natural unit to consider is instead the
individual (or population). These observations form the basis of our
developments; most of the time we assume the species level formulation is the
correct description of a given biological system (i.e., distinct species
evolve as species), and only for the (short) intermediate time interval
between the existence of parent and child species, we temporarily discard the
species-centric view and instead regard the relevant population as trait
specific distributions.

The species-interaction model we take as the starting point is formulated as a
system of ordinary differential equations (ODE) describing the dynamics of a
biological system in which any finite number of species interact and evolve
adaptively over evolutionary time (referred to in the sequel as the _species
level_ , SLM, or _macro-scale_ model). This model extends the Lotka-Volterra
system by representing each species not only by its abundance, but also by its
mean traits coordinate and trait covariance matrix. It originates in a recent
work by two of the authors [20], wherein the model equations were derived by
applying a moment closure averaging technique on a deterministic trait-
specific population density model (referred to in the sequel as the
_population level_ , PLM, or _micro-scale_ model), under the assumptions that
(A1) distinct species can be identified at all times, and (A2) that the
species distributions are of a known statistical quality (in this case, the
normal distribution). See also related derivations in e.g., [5, 25].

The population level model is formulated as a partial differential equation
(PDE), and governs the abundance density as a function of evolutionary traits.
Moreover, the authors showed in [20] that the species and population level
models are consistent, i.e., they describe the same system dynamics as long as
the system behavior is such that both models are valid. Still, these
descriptions are fundamentally different in that they operate on separate
scales; the species level model operates with the coarser unit of distinct
species, while the population level model operates with the finer unit of the
individual (i.e., the category of species is not imposed upon the individuals
who make up the population). Thus, the species level model described above
approximates the system dynamics governed by the population level model, under
the conditions that the population can be grouped into distinct species, and
that each species’ distribution conforms to a normal distribution. While this
requirement that trait distributions of a species are essentially normally
distributed may seem restrictive, we note that the shape of a distribution is
not an intrinsic quality of the observed system, but is equally dependent on
the scale of measurement. This assumption can thus be considered not as much
as a limitation on the biological system, as a constraint on the scale of
measurement. Note also that due to the deterministic nature of the population
level model, this property is also inherited at the species level.

The first part of our developments concerns monitoring the species in the
system and estimating the deviation from the true population density function.
We do this as follows: Since each species is described by an abundance-trait-
covariance tuple, we map these to time dependent distributions. Then, using
derived a posteriori error estimates of the population level model, we track
the modeling error of each species’ associated distribution function. The
second part concerns the case of residual blow-up, which implies a model
breakdown at the species level, here associated with a speciation event.
Concerning the actual speciation event we compare two different approaches
with regards to splitting the parent species into new child species; (1)
_Heuristic approach_ : When a residual blow-up is detected, we proceed to
split the relevant reconstructed distribution function along the trait
directions orthogonal to the direction of divergence. Each of these sub-
distributions are then mapped back to abundance-trait-variance tuples, thus
yielding the new child species to be incorporated at the species level. (2)
_Multi-scale approach_ : When a residual blow-up is detected, we delegate the
relevant species’ reconstructed density function to the population level model
to be solved in an appropriate local region of trait space while coupled to
the species level model which governs the interaction of the remaining species
in the system. Simultaneously, we map the local abundance density function to
abundance-trait-variance tuples and measure the distance between the mean
trait coordinates. At such time when the new child species are sufficiently
separated in trait space (e.g., by a multiple of the maximum trait standard
deviation), we incorporate the new child species at the species level and
decouple the multi-scale model. We assess the accuracy of each approach by
comparing the species parameters (pre- and post-speciation) to the
corresponding statistical moments of the reference (global) PLM solution.

Based on our results, we find that the multi-scale approach gives the better
approximation of the reference PLM solution post-speciation, and thus conclude
that a multi-scale approach is indeed necessary for robust eco-evolutionary
modeling at the level of species interaction, capable of handling speciation
events. Although we base our developments on the specific models described
above, our methodology is more general; it is applicable on any species
interaction model which, including the species abundance, also incorporates
the evolution of species mean traits and species trait covariance in some way
(thus allowing for the reconstruction of population density distributions),
whether the model is based on deterministic or stochastic dynamics. Figure 1
below shows a schematic representation of the micro/macro-scale coupling
strategy for speciation events.

Figure 1: Coupling of SLM and PLM during speciation events. Outside the green
area the macro-scale (SLM) model is solved (here, the curve represents the
reconstructed distribution function associated with each species’ abundance-
trait-covariance tuple), the green area indicates the micro-scale (PLM) model
is solved (here, the curve represents the abundance density, or, solution
function of the PLM).

The article is organized as follows: In Section 2 we present the relevant
population and species level models. In Section 3 we introduce the time
discrete population level model and proceed to analyze this. In Section 4 we
derive the a posteriori error bound of the population level model. In Section
5 we present in detail the methodology for species splitting and micro-macro
scale coupling. In Section 6 we refine the derived error bound for the case of
individual species. In Section 7 we provide numerical examples where we test
our multi-scale framework in detail. Finally, in Section 8 we provide some
concluding remarks.

##  2 Models

In this section we present both the (micro-scale) population level model (PLM)
and (macro-scale) species level model (SLM) as described in [19]. We also
briefly discuss the consistency of these model formulations.

###  2.1 The trait-specific population level model

The PLM describes the abundance of a population in terms of the population
density as a function of evolutionary traits over evolutionary time. The model
relates the rate of change in population density to birth/death rate of
individuals, interaction between individuals with different evolutionary
traits (cooperation/competition), and interaction between individuals with
equal evolutionary traits (self-limitation). The evolution of the population,
i.e., change in traits from one generation to the next, is incorporated as a
diffusion process.

Given an open, bounded domain $\Omega\subset{\mathbb{R}}^{d}$ in trait space
(i.e., the number of evolutionary traits considered is equal to $d\geq 1$),
and final time $T>0$, let $\Omega_{T}:=\Omega\times(0,T)$ be the space-time
domain, and let $n:\Omega_{T}\rightarrow{\mathbb{R}}$ be the trait-specific
population density, or, number of individuals per measure on
${\mathbb{R}}^{d}$. The domain $\Omega$ can then be regarded as all
evolutionary traits attainable by the relevant population over evolutionary
time, and the total number of individuals present at any time $t\in(0,T)$ is
then given by integration over all attainable traits, i.e.,
$\int_{\Omega}n(x,t)\textnormal{d}x$, where we require non-negativity of the
population density, i.e.,

(1) $n\geq 0,\textnormal{ for all }(x,t)\in\Omega_{T}.$

Furthermore, let $f:\Omega_{T}\rightarrow{\mathbb{R}}$ be the source term (net
rate of migration/immigration). The model then reads as the following
nonlinear and non-local equation:

(2) $\partial_{t}n-rn+\Phi(n)n-\nabla\cdot(g\nabla n)=f.$

The ecological processes consists of the growth term, $-rn$, where
$r:\Omega_{T}\rightarrow{\mathbb{R}}$ is the inherent per-capita growth-rate,
and the interaction term,

(3) $\Phi(n)n:=bn^{2}-n\int_{\Omega}\alpha(x,y,t)n(y,t)\textnormal{d}y,$

where $b:\Omega_{T}\rightarrow{\mathbb{R}}_{\geq 0}$ is the local individual
limitation coefficient, and
$\alpha:\Omega\times\Omega_{T}\rightarrow{\mathbb{R}}$ is the non-local
individual interaction coefficient, i.e., $\alpha(x,y,t)$ is the effect of an
individual with trait coordinate $y$ on an individual with trait coordinate
$x$. The diffusive term represents the evolutionary process, where
$g:\Omega_{T}\rightarrow{\mathbb{R}}^{d\times d}$ is the intergenerational
trait diffusion tensor. Generally speaking, the model equation (2) belongs to
the class of reaction-diffusion type parabolic conservation equations (in this
case, with a non-local reaction term). We remark that several generalizations
of (2) are possible. In particular, the diffusive term as stated models
incremental evolutionary processes: To allow for rare but possibly non-
incremental evolutionary processes more general non-local operators (such as
e.g., fractional derivatives), should be considered.

Since homogeneous boundary conditions of Dirichlet type is always justified
from a biological perspective given a large enough trait domain $\Omega$, we
shall only consider this situation for the present purposes. We let initial
data be specified by $n(x,0)=n_{0}(x)$. The PLM then reads as the following
initial/boundary value problem:

Find $n:\Omega_{T}\rightarrow{\mathbb{R}}$ such that $n\geq 0$ and

(4a) $\displaystyle\partial_{t}n-rn+\Phi(n)n-\nabla\cdot(g\nabla n)$
$\displaystyle=f,\quad$ in $\displaystyle\Omega_{T}$ (4b) $\displaystyle n$
$\displaystyle=0,$ on $\displaystyle\partial\Omega\times(0,T),$ (4c)
$\displaystyle n$ $\displaystyle=n_{0},$ in
$\displaystyle\Omega\times\\{0\\}.$

For additional details regarding the above model, see [20].

###  2.2 The species level model

The SLM describes the temporal evolution of the abundance, mean traits
coordinate, and trait covariance matrix of all species in an eco-evolutionary
system. Similarly to the PLM, the rate of change of the species abundance is
related to the growth/death rate, self-limitation and cooperation/competition,
but here among distinct species. Evolution is thus incorporated in the model
by the change in mean traits (i.e., location in trait space) and change in
trait covariance (i.e., spread in trait space).

Indexing the species present in the ecosystem by $i=1,\cdots,s$ (where $s\geq
1$ is the total number of species), we denote the abundance, mean traits
coordinate, and trait covariance matrix for species $i$ at time $t\in(0,T)$ by
$n_{i}(t)\in{\mathbb{R}}$, $x_{i}(t)\in{\mathbb{R}}^{d}$ and
$\upsilon_{i}(t)\in{\mathbb{R}}^{d\times d}$, respectively, and let
$N_{i}:=(n_{i},x_{i},\upsilon_{i})$ be the tuple representing the $i$’th
species, and similarly let $N:=(N_{1},\cdots,N_{s})$ represent the full
$s$-species ecosystem. Furthermore, let $F_{0,i}(t)\in{\mathbb{R}}$ be the
source of individuals for species $i$, with associated mean and covariance
$F_{1,i}(t)\in{\mathbb{R}}^{d}$ and $F_{2,i}(t)\in{\mathbb{R}}^{d\times d}$,
respectively. The model then reads as the following system of ODE’s (for
$i=1,\cdots,s$):

(5a) $\displaystyle\frac{\textnormal{d}n_{i}}{\textnormal{d}t}$
$\displaystyle=\mathbf{N}_{i}(N):=R_{i}n_{i}-B_{i}n_{i}^{2}+n_{i}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{s}A_{i,j}n_{j}+F_{0,i},$ $\displaystyle t\in(0,T),$
(5b) $\displaystyle\frac{\textnormal{d}x_{i}}{\textnormal{d}t}$
$\displaystyle=\mathbf{X}_{i}(N):=\upsilon_{i}\left(\nabla
r_{i}-\frac{n_{i}}{2\Lambda(\upsilon_{i})}\nabla
b_{i}+\sum_{j=1}^{s}n_{j}\nabla_{i}\alpha_{i,j}+F_{1,i}\right),\ $
$\displaystyle t\in(0,T),$ (5c)
$\displaystyle\frac{\textnormal{d}\upsilon_{i}}{\textnormal{d}t}$
$\displaystyle=\boldsymbol{\Sigma}_{i}(N):=V_{0,i}+V_{1,i}\upsilon_{i}+\upsilon_{i}V_{2,i}\upsilon_{i}+F_{2,i},$
$\displaystyle t\in(0,T),$

where $\Lambda(\upsilon_{i}):=\sqrt{(2\pi)^{d}}|\upsilon_{i}|^{\frac{1}{2}}$
and where (for $i,j=1,\cdots,s$)

(6a) $\displaystyle R_{i}$ $\displaystyle:=r_{i}+\frac{1}{2}\nabla\nabla
r_{i}:\upsilon_{i},$ (6b) $\displaystyle B_{i}$
$\displaystyle:=\frac{1}{\Lambda(\upsilon_{i})}\left(b_{i}+\frac{1}{4}\nabla\nabla
b_{i}:\upsilon_{i}\right)-\left(\alpha_{i,i}+\nabla_{i}\nabla_{i}\alpha_{i,i}:\upsilon_{i}\right),$
(6c) $\displaystyle A_{i,j}$
$\displaystyle:=\alpha_{i,j}+\frac{1}{2}\nabla_{i}\nabla_{i}\alpha_{i,j}:\upsilon_{i}+\frac{1}{2}\nabla_{j}\nabla_{j}\alpha_{i,j}:\upsilon_{j},$
(6d) $\displaystyle V_{0,i}$ $\displaystyle:=2g_{i},$ (6e) $\displaystyle
V_{1,i}$ $\displaystyle:=\frac{n_{i}b_{i}}{2\Lambda(\upsilon_{i})},$ (6f)
$\displaystyle V_{2,i}$ $\displaystyle:=\nabla\nabla
r_{i}+\frac{1}{4}\sum_{j=1}^{s}\nabla_{i}\nabla_{i}\alpha_{i,j}n_{j}+\frac{n_{i}}{4\Lambda(\upsilon_{i})}\left(\frac{1}{2}(\upsilon_{i}:\nabla\nabla
b_{i})\upsilon_{i}^{-1}-\nabla\nabla b_{i}\right).$

Here, we employed the following notational convention:
$r_{i}(t):=r(x_{i}(t),t)$, $\nabla r_{i}(t):=(\nabla r)(x_{i}(t),t)$ and
$\nabla\nabla r_{i}(t):=(\nabla\nabla r)(x_{i}(t),t)$, and similarly for the
other PLM coefficients. Subscripts on differential operators indicate on which
argument the operator is acting. For the full $s$-species ecosystem, we
abbreviate the above model as

(7)
$\mathcal{N}(N):=(\mathbf{N}_{1}(N),\mathbf{X}_{1}(N),\boldsymbol{\Sigma}_{1}(N),\cdots,\mathbf{N}_{s}(N),\mathbf{X}_{s}(N),\boldsymbol{\Sigma}_{s}(N)),$

and let initial data be given by

(8) $N(0)=N_{0}:=(N_{1,0},\cdots,N_{s,0}).$

The species level model then reads as the following initial value problem:

Find $N:(0,T)\rightarrow{\mathbb{R}}^{s(1+d+d^{2})}$ such that $n_{i}\geq 0$,
for $1\leq i\leq s$, and

(9a) $\displaystyle\dfrac{\textnormal{d}N}{\textnormal{d}t}$
$\displaystyle=\mathcal{N}(N),\quad$ $\displaystyle t\in(0,T),$ (9b)
$\displaystyle N$ $\displaystyle=N_{0},$ $\displaystyle t=0.$

In practice, the problem (9a)–(9b) can be reduced to $s(1+d+d(d+1)/2)$
equations due to the symmetry of the trait covariance matrix. For additional
details including a derivation of the above model, see [19]. Note also that
when indexing the components of a point $x\in\Omega$ we use subscripts, i.e.,
$x=(x_{1},\cdots,x_{d})$, not to be confused with the species mean trait
coordinate $x_{i}$, where we write $x_{i}(t)=(x_{i,1}(t),\cdots,x_{i,d}(t))$.

###  2.3 Relationship between population and species level models

If there is a time interval $J:=(t_{0},t_{1})\subset(0,T)$ where $s\geq 1$
species can unambiguously be identified, there exists a collection of non-
overlapping (possibly time dependent), open and bounded regions in trait
space, i.e., $\\{B_{i}\subset\Omega,i\leq s,B_{i}\cap B_{j}=\emptyset,i\neq
j\\}$, where the population density is compactly supported and contains only
one local maximum. Thus, if we let $\mu_{B_{i}}^{k}(n)$, be the $k$’th order
moment of $n$ (for $k\in\\{0,1,2\\}$, and assuming $n$ is smooth enough for
these moments to exist), within the region $B_{i}$, i.e.,

(10a) $\displaystyle\mu_{B_{i}}^{0}(n)$
$\displaystyle:=\int_{B_{i}}n(x)\textnormal{d}x,$ (10b)
$\displaystyle\mu_{B_{i}}^{1}(n)$
$\displaystyle:=\frac{1}{\mu_{B_{i}}^{0}(n)}\int_{B_{i}}n(x)x\textnormal{d}x,$
(10c) $\displaystyle\mu_{B_{i}}^{2}(n)$
$\displaystyle:=\frac{1}{\mu_{B_{i}}^{0}(n)}\int_{B_{i}}n(x)(x-\mu_{B_{i}}^{1})\otimes(x-\mu_{B_{i}}^{1})\textnormal{d}x,$

we can initialize the SLM at $t=t_{0}$ with initial data

(11)
$N_{0}=(\mu_{B_{1}}^{0},\mu_{B_{1}}^{1},\mu_{B_{1}}^{2},\cdots,\mu_{B_{s}}^{0},\mu_{B_{s}}^{1},\mu_{B_{s}}^{2})|_{t=t_{0}}.$

Solving (9a)–(9b) then amounts to approximating (10a)–(10c) for $t\in J$.

##  3 Time discrete PLM

In this section we introduce the time discrete PLM, which will serve as the
basis for the a posteriori error bounds:

For $k=0,\cdots,M$, we let $t_{k}$ denote the discrete times, such that
$0=t_{0}<t_{1}<\cdots<t_{M}:=T$, and where
$t_{k}=t_{k-1}+\tau_{k}=\sum_{j=1}^{k}\tau_{k}$, for some set of time
increments $\\{\tau_{k}\\}_{k=1}^{M}$ (in the following, we use superscripts
to indicate dependency on the discrete times, e.g., $u^{k}:=u(t_{k})$). The
time-discrete version of eq. 2 then reads as (for $k\geq 1$)

(12)
$n^{k}-\tau_{k}r^{k}n^{k}+\tau_{k}\Phi^{k}(n^{k})n^{k}-\tau_{k}\nabla\cdot(g^{k}\nabla
n^{k})=\tau_{k}f+n^{k-1},$

where we have discretized in time using a standard implicit first order method
(i.e., the backward Euler method). While we do not prove convergence of the
above discrete scheme to the original continuous problem, the local truncation
error for the backward Euler method is known to be of second order, and one
would expect first-order convergence of the time-discrete problem to the
continuous problem for sufficiently smooth problems.

We remark that the time-discrete problem has significant biological relevance
by itself, as many processes are naturally modeled on a discrete time-scale
due to the presence of strong temporal cycles (e.g., day, year, and
generation).

In the following, we first establish the existence of a solution to (12) in
the weak sense for every discrete time, then we derive an identity for the
residual which also provides uniqueness of the time-discrete solution.

###  3.1 Notation

We employ standard notation for function spaces. For a function space $V$ we
denote by $\lVert\cdot\rVert_{V}$ its energy norm, by $V^{\ast}$ its dual
space, and by $\langle\cdot,\cdot\rangle_{V^{\ast},V}$ the dual pairing. In
particular, for $1\leq p\leq\infty$, and a domain $D\subset\Omega$ we denote
by $L^{p}(D)$ the Lebesgue-spaces of integrable functions defined on $D$, with
associated norm $\lVert\cdot\rVert_{p,D}$ (domain subscript omitted if this is
clear from the context). Furthermore, let $L^{p}_{+}(D):=\\{u\in
L^{p}(D):u\geq 0\textnormal{ a.e.}\\}$. For the case $p=2$, we let $\lVert
u\rVert:=\lVert u\rVert_{2}=(u,u)^{\frac{1}{2}}$, where $(\cdot,\cdot)$ is the
standard $L^{2}$ inner product. Moreover, we denote by
$\mathbf{H}({\textnormal{div}},\Omega)$ the space of functions in
$[L^{2}(D)]^{d}$ admitting a weak divergence, and by $H^{1}(D)$ the space of
functions in $L^{2}(D)$ admitting weak gradients, with its vanishing trace
subspace denoted by $H_{0}^{1}(D)$, where $H^{-1}(D):=H_{0}^{1}(D)^{\ast}$.
Finally, denote by $H_{0,+}^{1}(D):=H_{0}^{1}(D)\cap L^{2}_{+}(D)$.

###  3.2 Preliminaries

Regarding the coefficients and source term of the PLM we introduce the
following assumption.

######  Assumption 1 (Coefficients and source term).

assume that $r,b,\alpha,g$ and $f$ are defined for all discrete times $t_{k}$
such that the following holds for all for $1\leq k\leq M$:

  1. 1.

$r^{k},\,b^{k}\in L^{\infty}(\Omega)$, and $\alpha^{k}\in
L^{\infty}(\Omega\times\Omega)$, such that

(13a) $|(r^{k}u,u)|\leq R^{k}\lVert u\rVert^{2},\quad\forall u\in
L^{2}(\Omega),$

for some set of positive constants $\\{R_{k}\\}_{k=1}^{M}$.

  2. 2.

The time increments $\tau_{k}$ are chosen such that

(13b) $\gamma_{k}:=1-\tau_{k}R_{k}>0.$

  3. 3.

$g^{k}\in[L^{\infty}(\Omega)]^{d\times d}$, symmetric and uniformly positive
definite, and such that

(13c)
$0<G_{k}|\zeta|^{2}\leq\zeta^{\textnormal{T}}g^{k}(x)\zeta,\quad\forall\zeta\in{\mathbb{R}}^{d}\setminus\\{0\\},$

for some set of positive constants $\\{G_{k}\\}_{k=1}^{M}$.

  4. 4.

$n_{0}\in H_{0,+}^{1}(\Omega)$.

  5. 5.

$f^{k}\in H^{-1}(\Omega)$.

Due to (13b) and (13c), we now equip the space $H_{0}^{1}(\Omega)$ with the
following equivalent inner product

(14) $(u,v)_{k}:=((1-\tau_{k}r^{k})u,v)+\tau_{k}(g^{k}\nabla u,\nabla v),$

and associated energy norm

(15)
${|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}^{2}_{k}:=(v,v)_{k},$

and define the Hilbert space $X_{k}:=(H_{0,+}^{1}(\Omega),(\cdot,\cdot)_{k})$,
and corresponding dual space $X_{k}^{\ast}$ with respect to $L^{2}(\Omega)$.
For $\xi\in X_{k}^{\ast}$, we denote the dual norm by

(16)
${|\kern-1.07639pt|\kern-1.07639pt|\xi|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}}:=\sup_{\underset{{|\kern-0.75346pt|\kern-0.75346pt|v|\kern-0.75346pt|\kern-0.75346pt|}_{k}=1}{v\in
X_{k},}}\langle\xi,v\rangle_{k},$

where
$\langle\cdot,\cdot\rangle_{k}:=\langle\cdot,\cdot\rangle_{X_{k}^{\ast},X_{k}}$
is the dual pairing. Furthermore, let $\mathcal{L}(X_{k},X_{k}^{\ast})$ be the
space of all linear maps from $X_{k}$ into $X_{k}^{\ast}$. For
$\Xi\in\mathcal{L}(X_{k},X_{k}^{\ast})$, we denote the operator norm by

(17)
${|\kern-1.07639pt|\kern-1.07639pt|\Xi|\kern-1.07639pt|\kern-1.07639pt|}_{k,k^{\ast}}:=\sup_{\underset{{|\kern-0.75346pt|\kern-0.75346pt|v|\kern-0.75346pt|\kern-0.75346pt|}_{k}=1}{v\in
X_{k},}}{|\kern-1.07639pt|\kern-1.07639pt|\Xi(v)|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}}.$

Due to Assumption 1 and the Poincaré inequality, there exists a constant
$c_{k,\Omega}>0$ (depending on $r^{k}$, $g^{k}$ and $\tau_{k}$ in addition to
the domain $\Omega$) such that

(18) $\lVert v\rVert\leq
c_{k,\Omega}{|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{k},\quad\forall
v\in X_{k}.$

For $1\leq k\leq M$ we now introduce the following assumption regarding the
nonlinearities of the PLM:

######  Assumption 2 (Nonlinear terms).

Assume that $\Phi^{k}(u)\in\mathcal{L}(X_{k},X_{k}^{\ast})$ for any $u\in
X_{k}$, such that

  1. 1.

(19a) $\langle\Phi^{k}(u)v,v\rangle_{k}\geq 0,\qquad\forall v\in X_{k}.$

  2. 2.

For every ball $\mathcal{B}_{a}\subset X_{k}$ of finite radius $a>0$, there
exists $L^{k}>0$ (depending on $k$ and $a$) such that there holds

(19b)
${|\kern-1.07639pt|\kern-1.07639pt|\Phi^{k}(u)-\Phi^{k}(v)|\kern-1.07639pt|\kern-1.07639pt|}_{k,k^{\ast}}\leq
L^{k}\lVert u-v\rVert_{1}\quad\forall u,v\in\mathcal{B}_{a}.$

  3. 3.

$\Phi^{k}(u)$ is a monotone operator, i.e., there holds

(19c) $\langle\Phi^{k}(u)u-\Phi^{k}(v)v,u-v\rangle_{k}\geq 0,\quad\forall
u,v\in X_{k}.$

With this, we introduce now the weak time-discrete formulation of the problem
(4a)–(4c): For $k\geq 1$ and given $n^{k-1}\in X_{k}$ find $n^{k}\in X_{k}$
such that the following integral equality holds

(20a)
$\displaystyle(n^{k},v)_{k}+\tau_{k}\langle\Phi^{k}(n^{k})n^{k},v\rangle_{k}=\tau_{k}\langle
f^{k},v\rangle_{k}+(n^{k-1},v),\qquad\forall v\in X_{k},$ and such that
initial condition (4c) is satisfied in the weak sense, i.e., (20b)
$(n^{0},v)=(n_{0},v),\qquad\forall v\in X_{k}.$

######  Remark 3.1 (Assumptions).

Assumption 1 is, to the best of the authors’ understanding, natural from a
biological point of view. On the other hand, Assumption 2 is introduced for
the analysis, and may be more restrictive than desired for certain biological
systems. We will return to this issue in the numerical examples (cf. Section
7).

###  3.3 Solvability

In this section we discuss the solvability of the weak time-discrete PLM (20).
We summarize the result in the following theorem.

######  Theorem 3.2 (Well posedness).

For $d\geq 3$ and given Assumptions 1 and 2, then for every $1\leq k\leq M$
there exists a unique solution $n^{k}\in X_{k}$ to problem (20) satisfying

(21)
${|\kern-1.07639pt|\kern-1.07639pt|n^{k}|\kern-1.07639pt|\kern-1.07639pt|}_{k}\leq\tau_{k}{|\kern-1.07639pt|\kern-1.07639pt|f^{k}|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}}+c_{k,\Omega}\lVert
n^{k-1}\rVert.$

######  Proof 3.3.

The proof follows by a series of calculations done in the next sections.
First, we show existence of a linearized problem in 3.3.1, which we then use
to infer existence to the nonlinear problem in 3.3.2. Then, we derive an
identity for the dual norm of the residual in 3.3.3, which we then apply to
obtain the uniqueness in 3.3.4. Our existence proof is based on techniques
from [1].

####  3.3.1 Linearization

We establish the well-posedness of a linearized problem. Choose
$\hat{\Phi}\in\mathcal{L}(X_{k},X_{k}^{\ast})$ such that
$\langle\hat{\Phi}u,u\rangle_{k}\geq 0$ for all $u\in X_{k}$, and define the
bilinear form $a^{k}:X_{k}\times X_{k}\rightarrow{\mathbb{R}}$ by

(22) $a^{k}(u,v):=(u,v)_{k}+\tau_{k}\langle\hat{\Phi}u,v\rangle_{k}.$

The linearized problem then reads: For $k\geq 1$ and given $n^{k-1}\in X_{k}$,
find $\hat{n}^{k}\in X_{k}$ such that

(23) $a^{k}(\hat{n}^{k},v)=\tau_{k}\langle
f^{k},v\rangle_{k}+(n^{k-1},v),\qquad\forall v\in X_{k}.$

The form $a^{k}$ is continuous and coercive on $X_{k}$:

(24a) $\displaystyle|a^{k}(u,v)|$
$\displaystyle\leq(1+\tau_{k}{|\kern-1.07639pt|\kern-1.07639pt|\hat{\Phi}|\kern-1.07639pt|\kern-1.07639pt|}_{k,k^{\ast}}){|\kern-1.07639pt|\kern-1.07639pt|u|\kern-1.07639pt|\kern-1.07639pt|}_{k}{|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{k},$
(24b) $\displaystyle a^{k}(u,u)$
$\displaystyle\geq{|\kern-1.07639pt|\kern-1.07639pt|u|\kern-1.07639pt|\kern-1.07639pt|}_{k}^{2}.$

Thus, by the Lax-Milgram Theorem [15] there exists a unique solution
$\hat{n}^{k}\in X_{k}$ to (23) satisfying

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\hat{n}^{k}|\kern-1.07639pt|\kern-1.07639pt|}_{k}$
$\displaystyle\leq{|\kern-1.07639pt|\kern-1.07639pt|\tau_{k}f^{k}+n^{k-1}|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}}$
(25)
$\displaystyle\leq\tau_{k}{|\kern-1.07639pt|\kern-1.07639pt|f^{k}|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}}+c_{k,\Omega}\lVert
n^{k-1}\rVert,$

since the action of $\tau_{k}f^{k}+n^{k-1}\in X_{k}^{\ast}$ on $X_{k}$ is
defined by the right hand side of (20a), and where we used the triangle
inequality and (18) in the second line.

####  3.3.2 Nonlinear problem

We establish existence of a solution to the nonlinear problem. For given $w\in
X_{k}$, let $\hat{n}^{k}_{w}\in X_{k}$ be the unique solution to the following
linear problem

(26)
$(\hat{n}^{k}_{w},v)_{k}+\tau_{k}\langle\Phi^{k}(w)\hat{n}^{k}_{w},v\rangle_{k}=\tau_{k}\langle
f^{k},v\rangle_{k}+(n^{k-1},v).$

Define the map $\mathcal{F}^{k}:X_{k}\rightarrow X_{k}$ by
$\mathcal{F}^{k}(w)=\hat{n}^{k}_{w}$. Due to the result of the previous
section this is well defined. Moreover, let

(27) $\mathcal{B}^{k}_{a}:=\\{w\in
X_{k}:{|\kern-1.07639pt|\kern-1.07639pt|w|\kern-1.07639pt|\kern-1.07639pt|}_{k}\leq{|\kern-1.07639pt|\kern-1.07639pt|\tau_{k}f^{k}+n^{k-1}|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}}=:a\\}.$

Due to (25), it follows that
$\mathcal{F}^{k}(\mathcal{B}_{a}^{k})\subset\mathcal{B}_{a}^{k}$.

Next step is to show that $\mathcal{F}^{k}$ is continuous and compact on
$\mathcal{B}^{k}_{a}$. Let
$\\{w_{j}\\}_{n=1}^{\infty}\subset\mathcal{B}^{k}_{a}$ be a sequence. Since
$X_{k}$ is a Hilbert space, and since $\\{w_{j}\\}_{n=1}^{\infty}$ is bounded,
the Eberlein-Šmulian Theorem [3] applies; there exists a subsequence (denoted
the same way) and $w\in X_{k}$ such that $w_{j}\rightharpoonup w$ weakly in
$X_{k}$. Moreover, since $d\geq 3$, by the Rellich-Kondrachov Theorem [15],
the embedding $X_{k}\subset L^{1}(\Omega)$ is compact. Thus, $w_{j}\rightarrow
w$ strongly in $L^{1}(\Omega)$. Now, let
$\hat{n}^{k}_{j}=\mathcal{F}^{k}(w_{j})$. Then $\hat{n}^{k}_{j}$ solves

(28)
$(\hat{n}^{k}_{j},v)_{k}+\tau_{k}\langle\Phi^{k}(w_{j})\hat{n}^{k}_{j},v\rangle=\tau_{k}\langle
f^{k},v\rangle_{k}+(n^{k-1},v).$

Furthermore, $\hat{n}^{k}_{w}=\mathcal{F}^{k}(w)$ solves (26). This
establishes the continuity of $\mathcal{F}^{k}$.

Next step is to show $\hat{n}^{k}_{j}\rightarrow\hat{n}^{k}_{w}$ strongly in
$X_{k}$. To this end, subtract (26) from (28) to obtain

(29)
$(\hat{n}^{k}_{j}-\hat{n}^{k}_{w},v)_{k}+\tau_{k}\langle\Phi^{k}(w_{j})(\hat{n}^{k}_{j}-\hat{n}^{k}_{w}),v\rangle_{k}=-\tau_{k}\langle(\Phi^{k}(w_{j})-\Phi^{k}(w))\hat{n}^{k}_{j},v\rangle_{k}.$

Since we have

$\displaystyle|\langle(\Phi^{k}(w_{j})-\Phi^{k}(w))\hat{n}^{k}_{j},v\rangle_{k}|$
$\displaystyle\leq{|\kern-1.07639pt|\kern-1.07639pt|(\Phi^{k}(w_{j})-\Phi^{k}(w))\hat{n}^{k}_{j}|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}}{|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{k}$
$\displaystyle\leq{|\kern-1.07639pt|\kern-1.07639pt|\Phi^{k}(w_{j})-\Phi^{k}(w)|\kern-1.07639pt|\kern-1.07639pt|}_{k,k^{\ast}}{|\kern-1.07639pt|\kern-1.07639pt|\hat{n}^{k}_{j}|\kern-1.07639pt|\kern-1.07639pt|}_{k}{|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{k}$
(30) $\displaystyle\leq L^{k}\lVert
w_{j}-w\rVert_{1}{|\kern-1.07639pt|\kern-1.07639pt|\tau_{k}f^{k}+n^{k-1}|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}}{|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{k},$

it follows that

(31)
$|\langle(\Phi^{k}(w_{j})-\Phi^{k}(w))\hat{n}^{k}_{j},v\rangle_{k}|\rightarrow
0,\qquad\forall v\in X_{k},$

and consequently $\hat{n}^{k}_{j}\rightarrow\hat{n}^{k}_{w}$ in $X_{k}$ as
$j\rightarrow\infty$, since the equation (29) holds for any $v\in X_{k}$ and
is linear in the argument $\hat{n}^{k}_{j}-\hat{n}^{k}_{w}$. Thus,
$\mathcal{F}^{k}$ is a continuous and compact operator on
$\mathcal{B}^{k}_{a}$, and the Schauder Theorem [3] applies; there exists
$n^{k}\in\mathcal{B}^{k}_{a}$ such that $n^{k}=\mathcal{F}^{k}(n^{k})$, which
is the solution to the nonlinear problem (20).

####  3.3.3 Residual identity

Based on the previous existence result, we introduce now the residual operator
$\mathcal{R}^{k}(\psi)\in X_{k}^{\ast}$ for a given $\psi\in X_{k}$ and for
all $v\in X_{k}$ as

(32) $\displaystyle\langle\mathcal{R}^{k}(\psi),v\rangle_{k}:=\tau_{k}\langle
f^{k},v\rangle_{k}+(n^{k-1},v)-(\psi,v)_{k}-\tau_{k}\langle\Phi^{k}(\psi)\psi,v\rangle_{k}.$

We also introduce the semi-metric (cf. the monotonicity condition (19c))

$\displaystyle\mathcal{J}^{k}(n^{k},\psi)$
$\displaystyle:=\tau_{k}^{2}{|\kern-1.07639pt|\kern-1.07639pt|\Phi^{k}(n^{k})n^{k}-\Phi^{k}(\psi)\psi|\kern-1.07639pt|\kern-1.07639pt|}^{2}_{k^{\ast}}$
(33)
$\displaystyle\qquad+2\tau_{k}\langle\Phi^{k}(n^{k})n^{k}-\Phi^{k}(\psi)\psi,n^{k}-\psi\rangle_{k}.$

The following Lemma now makes precise the connection between the dual norm of
the residual and the resulting error.

######  Lemma 3.4 (Residual-error identity).

For all $\psi\in X_{k}$ there holds

(34)
${|\kern-1.07639pt|\kern-1.07639pt|\mathcal{R}^{k}(\psi)|\kern-1.07639pt|\kern-1.07639pt|}^{2}_{k^{\ast}}={|\kern-1.07639pt|\kern-1.07639pt|n^{k}-\psi|\kern-1.07639pt|\kern-1.07639pt|}^{2}_{k}+\mathcal{J}^{k}(n^{k},\psi).$

######  Proof 3.5.

Choose $\phi\in X_{k}$ such that

(35)
$(\phi,v)_{k}=\tau_{k}\langle\Phi^{k}(n^{k})n^{k}-\Phi^{k}(\psi)\psi,v\rangle_{k},\qquad\forall
v\in X_{k}.$

Thus, we have the identity

(36)
${|\kern-1.07639pt|\kern-1.07639pt|\phi|\kern-1.07639pt|\kern-1.07639pt|}_{k}=\tau_{k}{|\kern-1.07639pt|\kern-1.07639pt|\Phi^{k}(n^{k})n^{k}-\Phi^{k}(\psi)\psi|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}}.$

We can then write the dual norm of the residual as follows

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{R}^{k}(\psi)|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}}=$
$\displaystyle\sup_{\underset{{|\kern-0.75346pt|\kern-0.75346pt|v|\kern-0.75346pt|\kern-0.75346pt|}_{k}=1}{v\in
X_{k},}}\\{(n^{k}-\psi,v)_{k}+\tau_{k}\langle\Phi^{k}(n^{k})n^{k}-\Phi^{k}(\psi)\psi,v\rangle_{k}\\}$
$\displaystyle=$
$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\phi+n^{k}-\psi|\kern-1.07639pt|\kern-1.07639pt|}_{k}$
$\displaystyle=$
$\displaystyle\Big{(}{|\kern-1.07639pt|\kern-1.07639pt|\phi|\kern-1.07639pt|\kern-1.07639pt|}_{k}^{2}+2(\phi,n^{k}-\psi)_{k}+{|\kern-1.07639pt|\kern-1.07639pt|n^{k}-\psi|\kern-1.07639pt|\kern-1.07639pt|}_{k}^{2}\Big{)}^{\frac{1}{2}}$
$\displaystyle=$
$\displaystyle\Big{(}\tau_{k}^{2}{|\kern-1.07639pt|\kern-1.07639pt|\Phi^{k}(n^{k})n^{k}-\Phi^{k}(\psi)\psi|\kern-1.07639pt|\kern-1.07639pt|}^{2}_{k^{\ast}}$
(37) $\displaystyle\hskip
56.9055pt+2\tau_{k}\langle\Phi^{k}(n^{k})n^{k}-\Phi^{k}(\psi)\psi,n^{k}-\psi\rangle_{k}+{|\kern-1.07639pt|\kern-1.07639pt|n^{k}-\psi|\kern-1.07639pt|\kern-1.07639pt|}_{k}^{2}\Big{)}^{\frac{1}{2}}.$

The assertion follows.

######  Remark 3.6 (Semi-metric).

Due to the nonlinear nature of the model equation (20), the norm of the
residual can not be shown to be equal to some energy norm of the error.
Instead, it is a combination of a norm and a semi-metric. The semi-metric
satisfies non-negativity, symmetry, and identity of indiscernibles, but not
the triangle inequality. See e.g., [30], where a nonlinear parabolic equation
was analyzed and similar terms appear.

####  3.3.4 Uniqueness

It remains to show the uniqueness of $n^{k}$. To this end, we suppose that
$\psi$ satisfies (20), which implies that
${|\kern-1.07639pt|\kern-1.07639pt|\mathcal{R}^{k}(\psi)|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}}=0$.
Then, since the first term on the right hand side of (34) is a norm and since
the second term is convex due to the monotonicity condition (19c), we obtain
that $\psi=n^{k}$ in $X_{k}$. This section concludes the proof of Theorem 3.2.

##  4 A posteriori error bound

In this section we adopt energy-type a posteriori error bounds based on the
dual norm of the residual for the weak time-discrete PLM, i.e., the problem
(20) using similar techniques as in [13, 14, 22, 26]. We then provide a
guaranteed and fully computable upper bound on the energy–type error in terms
of various error estimators, using suitable flux and density reconstructions.

###  4.1 Space approximations

Let $\mathcal{Z}_{h}$ be a partition of the domain $\Omega$, consisting of
rectangular or simplicial elements, such that
$\bar{\Omega}=\cup_{K\in\mathcal{Z}_{h}}K$. We let $h_{K}$ denote the diameter
of the element $K\in\mathcal{Z}_{h}$, and define
$h:=\max_{K\in\mathcal{Z}_{h}}h_{K}$ as the maximum diameter over all elements
in $\mathcal{Z}_{h}$. For two elements $K,L\in\mathcal{Z}_{h}$, we require
their intersection to be either a common face, edge, vertex, or the empty set.
We use subscripts to indicate dependency on the discrete mesh
$\mathcal{Z}_{h}$, e.g., $u_{h}$. Our a posteriori error estimates given next
will involve approximations of coefficients and source term appearing in (20).
Thus, for any function
$\varphi^{k}\in\left\\{r^{k},\,b^{k},\alpha^{k},\,g^{k},\,f^{k}\right\\}$ (see
Assumption 1), we denote by $\varphi^{k}_{h}$ its approximation satisfying the
following orthogonality

(38) $(\varphi^{k}_{h}-\varphi^{k},v)_{K}=0,\quad\forall v\in L^{2}(\Omega),$

with the approximation of the interaction $\Phi_{h}^{k}$ corresponding to (3).
Furthermore, for a convex element $K\in\mathcal{Z}_{h}$, there holds the
Poincaré–Friedrichs inequality (see e.g., [27]), i.e.,

(39) $\lVert v-v_{k}\rVert_{K}\leq\frac{h_{K}}{\pi}\lVert\nabla
v\rVert_{K},\quad\forall v\in H^{1}(K),$

where $v_{k}$ is the mean value of $v$ over $K$.

###  4.2 First upper bound

To be agnostic with regards to the scheme (micro or/and macro levels) we use
to approximate the weak solution, we make the following definitions: [Density
reconstruction] We call a density reconstruction any function $s_{h}^{k}\in
X_{k}$. [Equilibrated flux and density reconstructions] We call an
equilibrated flux reconstruction any function
$\boldsymbol{\sigma}_{h}^{k}:\Omega\rightarrow{\mathbb{R}}^{d}$ which
satisfies:

(40a)
$\displaystyle\boldsymbol{\sigma}_{h}^{k}\in\mathbf{H}({\textnormal{div}},\Omega),$
$\displaystyle(\tau_{k}\nabla\cdot\boldsymbol{\sigma}^{k}_{h},1)_{K}$ (40b)
$\displaystyle\
=(\tau_{k}f^{k}_{h}+s_{h}^{k-1}-(1-\tau_{k}r^{k}_{h})s_{h}^{k}-\tau_{k}\Phi^{k}_{h}(s_{h}^{k})s_{h}^{k},1)_{K},\
$ $\displaystyle 1\leq k\leq M,\forall K\in\mathcal{Z}_{h}.$

For all $K\in\mathcal{Z}_{h}$, define the local residual, flux, r-data,
$\Phi$-data, g-data, and f-data oscillation estimators:

(41a) $\displaystyle\eta_{\textnormal{R},K}^{k}$
$\displaystyle:=\tilde{\omega}^{k}_{K}\lVert\tau_{k}f_{h}^{k}+n^{k-1}-(1-\tau_{k}r_{h}^{k})s_{h}^{k}-\tau_{k}\Phi_{h}^{k}(s_{h}^{k})s_{h}^{k}-\tau_{k}\nabla\cdot\boldsymbol{\sigma}^{k}_{h}\rVert_{K}$
(41b) $\displaystyle\eta_{\textnormal{DF},K}^{k}$
$\displaystyle:=\tau_{k}G_{k}^{-\frac{1}{2}}\lVert g_{h}^{k}\nabla
s_{h}^{k}+\boldsymbol{\sigma}^{k}_{h}\rVert_{K},$ (41c)
$\displaystyle\eta_{r,K}^{k}$
$\displaystyle:=\tau_{k}\tilde{\omega}^{k}_{K}\lVert(r^{k}-r_{h}^{k})s_{h}^{k}\rVert_{K},$
(41d) $\displaystyle\eta_{\Phi}^{k}$
$\displaystyle:=\tau_{k}{|\kern-1.07639pt|\kern-1.07639pt|(\Phi^{k}(s_{h}^{k})-\Phi_{h}^{k}(s_{h}^{k}))s_{h}^{k}|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}},$
(41e) $\displaystyle\eta_{g,K}^{k}$
$\displaystyle:=\tau_{k}G_{k}^{-\frac{1}{2}}\lVert(g^{k}-g_{h}^{k})\nabla
s_{h}^{k}\rVert_{K},$ (41f) $\displaystyle\eta_{f}^{k}$
$\displaystyle:=\tau_{k}{|\kern-1.07639pt|\kern-1.07639pt|f^{k}-f_{h}^{k}|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}}.$

where we introduced also the following weight

(42)
$\displaystyle\tilde{\omega}^{k}_{K}:=\min(\gamma_{k}^{-\frac{1}{2}},G_{k}^{-\frac{1}{2}}\frac{h_{K}}{\pi}).$

We now give the first (abstract) a posteriori error bound of the problem (20).

######  Theorem 4.1 (Energy-error bound).

For $1\leq k\leq M$, let $n^{k}$ be the (unknown) exact solution to the weak-
time discrete problem (20), $s_{h}^{k}$ the reconstructed density of
Definition 4.2, and $\boldsymbol{\sigma}_{h}^{k}$ is the reconstructed flux of
Definition 4.2. Then, the following estimate holds true

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|n^{k}-s_{h}^{k}|\kern-1.07639pt|\kern-1.07639pt|}^{2}_{k}+\mathcal{J}^{k}(n^{k},s_{h}^{k})$
(43)
$\displaystyle\qquad\quad\leq\left(\left\\{\sum_{K\in\mathcal{Z}_{h}}\left[\eta_{\textnormal{R},K}^{k}+\eta_{r,K}^{k}+\eta_{\textnormal{DF},K}^{k}+\eta_{g,K}^{k}\right]^{2}\right\\}^{\frac{1}{2}}+\eta_{\Phi}^{k}+\eta_{f}^{k}\right)^{2}.$

######  Proof 4.2.

We let $\psi=s_{h}^{k}$ in (32), then we subtract the terms
$\tau_{k}(\boldsymbol{\sigma}_{h}^{k},\nabla
v)+\tau_{k}(\nabla\cdot\boldsymbol{\sigma}^{k}_{h},v)=0$, followed by adding
and subtracting the terms
$\tau_{k}(r^{k}_{h}s_{h}^{k},v)$,
$\tau_{k}(\Phi_{h}^{k}(s_{h}^{k})s_{h}^{k},v)$, $\tau_{k}(g_{h}^{k}\nabla
s_{h}^{k},\nabla v)$, and $\tau_{k}(f_{h}^{k},v)$, which leads to

$\displaystyle\langle\mathcal{R}^{k}(s^{k}_{h}),v\rangle_{k}$
$\displaystyle=(\tau_{k}f_{h}^{k}+n^{k-1}-(1-\tau_{k}r_{h}^{k})s_{h}^{k}-\tau_{k}\Phi_{h}^{k}(s^{k}_{h})s_{h}^{k}-\tau_{k}\nabla\cdot\boldsymbol{\sigma}_{h}^{k},v).$
$\displaystyle\qquad-\tau_{k}(g_{h}^{k}\nabla
s_{h}^{k}+\boldsymbol{\sigma}_{h}^{k},\nabla
v)+\tau_{k}((r^{k}-r^{k}_{h})s_{h}^{k},v)$
$\displaystyle\qquad+\tau_{k}\\{(\Phi^{k}_{k}(s_{h}^{k})s_{h}^{k},v)-\langle\Phi_{h}^{k}(s_{h}^{k})s_{h}^{k},v\rangle_{k}\\}$
$\displaystyle\qquad+\tau_{k}((g_{h}^{k}-g^{k})\nabla s_{h}^{k},\nabla
v)+\tau_{k}\\{\langle f^{k},v\rangle_{k}-(f_{h}^{k},v)\\}$ (44)
$\displaystyle=:\sum_{i=1}^{6}T_{i}^{k},$

for all $v\in X_{k}$. Recalling the local estimators (41), we estimate each of
the terms $T_{1}^{k}$–$T_{6}^{k}$ on the right hand side of (44), i.e., for
the first term, we can use (13b) to get

$\displaystyle
T_{1}^{k}\leq\sum_{K\in\mathcal{Z}_{h}}\gamma^{-\frac{1}{2}}_{k}\lVert\tau_{k}f_{h}^{k}+n^{k-1}-(1-\tau_{k}r_{h}^{k})s_{h}^{k}-$
$\displaystyle\tau_{k}\Phi_{h}^{k}(s^{k}_{h})s_{h}^{k}-\tau_{k}\nabla\cdot\boldsymbol{\sigma}_{h}^{k}\rVert_{K}$
(45) $\displaystyle\cdot\lVert(1-\tau_{k}r^{k})^{\frac{1}{2}}v\rVert_{K}.$

Also, we can use the equilibration property (40b) together with the
orthogonality (38) to obtain

$\displaystyle T_{1}^{k}$
$\displaystyle=\sum_{K\in\mathcal{Z}_{h}}(\tau_{k}f_{h}^{k}+n^{k-1}-(1-\tau_{k}r_{h}^{k})s_{h}^{k}-\tau_{k}\Phi_{h}^{k}(s^{k}_{h})s_{h}^{k}-\tau_{k}\nabla\cdot\boldsymbol{\sigma}_{h}^{k},v)_{K},$
$\displaystyle=\sum_{K\in\mathcal{Z}_{h}}(\tau_{k}f_{h}^{k}+n^{k-1}-(1-\tau_{k}r_{h}^{k})s_{h}^{k}-\tau_{k}\Phi_{h}^{k}(s^{k}_{h})s_{h}^{k}$
(46) $\displaystyle\hskip
170.71652pt-\tau_{k}\nabla\cdot\boldsymbol{\sigma}_{h}^{k},v-v_{K})_{K}.$

By applying the Poincaré–Friedrichs inequality (39) we get

$\displaystyle T_{1}^{k}\leq
G_{k}^{-\frac{1}{2}}\frac{h_{K}}{\pi}\sum_{K\in\mathcal{Z}_{h}}\lVert\tau_{k}f_{h}^{k}+n^{k-1}-(1-\tau_{k}r_{h}^{k})s_{h}^{k}-$
$\displaystyle\tau_{k}\Phi_{h}^{k}(s^{k}_{h})s_{h}^{k}-\tau_{k}\nabla\cdot\boldsymbol{\sigma}_{h}^{k}\rVert_{K}$
(47) $\displaystyle\cdot\lVert(g^{k})^{\frac{1}{2}}\nabla v\rVert_{K}.$

Combining (4.2) and (47), we obtain

(48)
$T_{1}^{k}\leq\sum_{K\in\mathcal{Z}_{h}}\eta_{\textnormal{R},K}^{k}{|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{K,k}.$

The same procedure applied to the term $T^{k}_{3}$ yields

(49)
$T_{3}^{k}\leq\sum_{K\in\mathcal{Z}_{h}}\eta_{r,K}^{k}{|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{K,k}.$

On the terms $T_{2}$ and $T_{5}$ we use (13c) to obtain

(50a) $\displaystyle T_{2}^{k}$
$\displaystyle\leq\sum_{K\in\mathcal{Z}_{h}}\eta_{\textnormal{DF},K}^{k}\lVert(g^{k})^{\frac{1}{2}}\nabla
v\rVert_{K},$ (50b) $\displaystyle T_{5}^{k}$
$\displaystyle\leq\sum_{K\in\mathcal{Z}_{h}}\eta_{g,K}^{k}\lVert(g^{k})^{\frac{1}{2}}\nabla
v\rVert_{K}.$

Finally, on $T_{4}$ and $T_{6}$ we use the definition of the dual norm (16) to
obtain

(51a) $\displaystyle T_{4}^{k}$
$\displaystyle\leq\eta_{\Phi}^{k}{|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{k},$
(51b) $\displaystyle T_{6}^{k}$
$\displaystyle\leq\eta_{f}^{k}{|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{k}.$

Combining the above bounds leads to

$\displaystyle\langle\mathcal{R}^{k}(s_{h}^{k}),v\rangle_{k}$
$\displaystyle\quad\leq\sum_{K\in\mathcal{Z}_{h}}\left\\{\left[\eta_{\textnormal{R},K}^{k}+\eta_{r,K}^{k}+\eta_{\textnormal{DF},K}^{k}+\eta_{g,K}^{k}\right]{|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{K,k}\right\\}+(\eta_{f}^{k}+\eta_{\Phi}^{k}){|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{k}$
(52)
$\displaystyle\quad\leq\left(\left\\{\sum_{K\in\mathcal{Z}_{h}}\left[\eta_{\textnormal{R},K}^{k}+\eta_{r,K}^{k}+\eta_{\textnormal{DF},K}^{k}+\eta_{g,K}^{k}\right]^{2}\right\\}^{\frac{1}{2}}+\eta_{f}^{k}+\eta_{\Phi}^{k}\right){|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{k},$

since
${|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{k}^{2}=\sum_{K\in\mathcal{Z}_{h}}{|\kern-1.07639pt|\kern-1.07639pt|v|\kern-1.07639pt|\kern-1.07639pt|}_{K,k}^{2}$.
By (16) and (34) we prove the assertion.

##  5 Speciation

In this section we present our strategies for the actual speciation event. The
first is based on splitting the support of the relevant reconstructed density
function into subregions of trait space, and mapping back to new species
abundance-trait-variance tuples by calculating the statistical moments of the
reconstructed density in each subregion. The second approach is based on
coupling the PLM and SLM in a multi-scale framework for the duration of the
speciation event, where the diverging species is delegated to the PLM to be
solved locally in trait space.

###  5.1 Trait space density regions

For any integer $\nu=1,2,\cdots$, and time step $k\geq 1$, we associate with
species $i$ the ($d$-hyper-rectangular) density region,
$B_{i,\nu}^{k}\subseteq\Omega$ centered on $x_{i}^{k}$ and with orientation
according to the (orthonormal) spectrum of $\upsilon^{k}_{i}(t)$, and with
side lengths equal to
$2\nu(\lambda_{1}^{k})^{1/2},\cdots,2\nu(\lambda_{d}^{k})^{1/2}$, where
$\\{\lambda_{j}^{k},j\leq d\\}$ are the eigenvalues of $\upsilon_{i}^{k}$
(i.e., for any time discrete time $t_{k}$ the density region $B_{i,\nu}^{k}$
extends $\nu$ times the $j$’th standard deviation of the trait covariance
along the $j$’th trait dimension from the mean traits coordinate $x_{i}^{k}$).

###  5.2 Heuristic approach

If at time $t_{k}$ there is detected diverging traits within a species
$N_{i}^{k}$, we split this into two child species, $N_{i,1}^{k}$ and
$N_{i,2}^{k}$, as follows: The associated density region $B_{i,\nu}^{k}$ is
divided along the directions orthogonal to the largest eigenvector of the
trait covariance matrix, $\upsilon_{i}^{k}$, into two sub-regions
$B_{i,\nu}^{k,1}$ and $B_{i,\nu}^{k,2}$ such that $B_{i,\nu}^{k,1}\cup
B_{i,\nu}^{k,2}=B_{i,\nu}^{k}$, and the new species are initialized by the
moments of the reconstructed density function of the parent species in each
sub-region, i.e., for $j\in\\{1,2\\}$, we define

(53)
$N_{i,j}^{k}:=(\mu_{B_{i,\nu}^{k,j}}^{0}(\mathcal{D}N^{k}_{i}),\mu_{B_{i,\nu}^{k,j}}^{1}(\mathcal{D}N^{k}_{i}),\mu_{B_{i,\nu}^{k,j}}^{2}(\mathcal{D}N^{k}_{i})),$

where $\mathcal{D}:{\mathbb{R}}^{1+d+d^{2}}\rightarrow C^{\infty}(\Omega)\cap
X_{k}$ is the density reconstruction operator (i.e., macro-to-micro scale
mapping), such that $\mathcal{D}N^{k}_{i}$ is any smooth distribution
characterized by a mean and variance, and adhering to the relevant boundary
condition (in this case, homogenous Dirichlet). Figure 2 below shows the
splitting of an example density region in 2D, where $x_{i}^{k}=(0,0)$ and
$\upsilon_{i}^{k}=\operatorname{diag}(\upsilon_{i,11}^{k},\upsilon_{i,22}^{k})$.

(a) Density region, with spectrum of
covariance matrix scaled with $\nu$.

(b) Splitting of density region along direction orthogonal to direction of
trait divergence into two subregions.

Figure 2: Splitting of example density region in 2D. Dotted line represents
species distribution.

###  5.3 Multi-scale approach

In this section we couple the PLM, (4a)–(4c), and the SLM, (9a)–(9b), in a
multi-scale framework, following [29]. In particular, if a speciation event is
detected for species $i$ at time step $k\geq 1$, we delegate its reconstructed
density function to the PLM to be solved locally in the associated density
region of trait space, while coupled to the SLM which governs the remaining
$s-1$ species in the system. With this, we write the coupled multi-scale model
as

Micro-scale: (54a)
$\displaystyle\quad\begin{cases}\partial_{t}n^{\ell}-rn^{\ell}+\Phi(n^{\ell})n^{\ell}-\nabla\cdot(g\nabla
n^{\ell})+n^{\ell}{\mathbb{C}}^{\ell}(N^{\ell})=f^{\ell},&\textnormal{ in
}B_{i,\nu}^{k},\,\ell>k\\\ n^{k}=\mathcal{D}N_{i}^{k},&\textnormal{ in
}B_{i,\nu}^{k}.\end{cases}$ Macro-scale: (54b)
$\displaystyle\quad\begin{cases}\dfrac{\textnormal{d}N^{\ell}}{\textnormal{d}t}=\mathcal{N}(N^{\ell})+{\mathbb{D}}^{\ell}(n^{\ell})N^{\ell},&\ell>k\\\\[4.30554pt]
N^{k}=(N_{1}^{k},\cdots,N_{i-1}^{k},0,N_{i+1}^{k},\cdots,N_{s}^{k}),&\end{cases}$

where $\partial_{t}n^{\ell}$ and
$\dfrac{\textnormal{d}N^{\ell}}{\textnormal{d}t}$ denotes the discrete time
derivatives for the PLM and SLM, respectively, and where
${\mathbb{C}}^{\ell}:{\mathbb{R}}^{s(1+d+d^{2})}\rightarrow{\mathbb{R}}$ and
${\mathbb{D}}^{\ell}:H_{0,+}^{1}(B_{i,\nu}^{k})\rightarrow{\mathbb{R}}^{s(1+d+d^{2})}$
are the micro/macro couplings, approximating respectively the individual and
species interactions. These are defined as follows:

(55a) $\displaystyle({\mathbb{C}}^{\ell}(N^{\ell}))(x)$
$\displaystyle:=\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{s}\int_{B_{j,\nu}^{\ell}}\alpha^{\ell}(x,y)(\mathcal{D}N_{j}^{\ell})(y)\textnormal{d}y,$
$\displaystyle{\mathbb{D}}^{\ell}(n^{\ell}):=($ $\displaystyle
n_{1}^{\ell}A_{1,i}^{\ell},\ \nabla_{1}\alpha_{1,i}^{\ell},\
\upsilon_{1}^{\ell}\nabla_{1}\nabla_{1}\alpha_{1,i}^{\ell}\upsilon_{1}^{\ell}/4,$
$\displaystyle\cdots,$ $\displaystyle n_{i-1}A_{i-1,i},\
\nabla_{i-1}\alpha_{i-1,i},\
\upsilon_{i-1}^{2}\nabla_{i-1}\nabla_{i-1}\alpha_{i-1,i}\upsilon_{i-1}^{2}/4,$
$\displaystyle 0,\,$ $\displaystyle n_{i+1}^{\ell}A_{i+1,i}^{\ell},\
\nabla_{i+1}\alpha_{i+1,i}^{\ell},\
\upsilon_{i+1}^{\ell}\nabla_{i+1}\nabla_{i+1}\alpha_{i+1,i}^{\ell}\upsilon_{i+1}^{\ell}/4,$
(55b) $\displaystyle\cdots,$ $\displaystyle n_{s}^{\ell}A_{s,i}^{\ell},\
\nabla_{s}\alpha_{s,i}^{\ell},\
\upsilon_{s}^{\ell}\nabla_{s}\nabla_{s}\alpha_{s,i}^{\ell}\upsilon_{s}^{\ell}/4)|_{\mathcal{C}_{B^{k}_{i,\nu}}(n^{\ell})}\int_{B^{k}_{i,\nu}}n^{\ell}(x)\textnormal{d}x.$

During the speciation event, we map the local abundance density function to
$m$ ‘virtual’ species, i.e.,

(56)
$\hat{N}^{\ell}=(\hat{N}_{1}^{\ell},\cdots,\hat{N}^{\ell}_{m}):=\mathcal{C}^{m}_{B_{i,\nu}^{k}}(n^{\ell}),$

where
$\mathcal{C}^{m}_{B_{i,\nu}^{k}}(n^{\ell}):H_{0}^{1}(B_{i,\nu}^{k})\rightarrow{\mathbb{R}}^{m(1+d+d^{2})}$
is the $m$-species compression operator, defined by

(57)
$\mathcal{C}^{m}_{B_{i,\nu}^{k}}(n^{\ell}):=\operatorname*{arg\,min}_{\begin{subarray}{c}n_{j},x_{j},\upsilon_{j},\\\
1\leq i\leq
m\end{subarray}}\Big{\lVert}n^{\ell}-\sum_{j=1}^{m}\mathcal{D}(n_{j},x_{j},\upsilon_{j})\Big{\rVert}^{2}_{B_{i,\nu}^{k}}.$

The speciation event is then complete when the distances between the mean
trait coordinates of the virtual species are larger than some specified
tolerance (e.g., some multiple of the largest trait standard deviation). At
this point, we initialize the SLM with the new species configuration
consisting of $s-1+m$ species, i.e.,

(58)
$(N_{1}^{\ell},\cdots,N_{i-1}^{\ell},N_{i+1},\cdots,N_{s}^{\ell})\cup\hat{N}^{\ell},$

and decouple the micro and macro scales until the next speciation event is
detected.

##  6 Application: distinguishing the modeling error

In this section we apply the results from Section 4 to derive an energy-type a
posteriori error bound for the heuristic and multi-scale methods. The error
bound distinguishes in particular the micro-to-macro modeling error, which
will be employed for the purposes of detecting speciation events in the SLM.

###  6.1 Abundance density reconstruction and flux equilibration
with smooth sources

To apply the results of Theorem 43, we reconstruct the species density with
smooth distributions, and for the sake of simplicity we only consider smooth
data. At the time step $k\geq 0$, for the species $i\leq s$, we construct its
density function by setting

(59) $s_{h,i}^{k}:=\mathcal{D}N^{k}_{i}.$

Thus, the global ($s$-species) abundance density reconstruction is
$H^{1}$-conforming in space i.e.,

(60) $s_{h}^{k}:=\sum_{i=1}^{s}s^{k}_{h,i}\in C^{\infty}(\Omega)\cap X_{k}.$

What remains is to construct the equilibrated flux,
$\boldsymbol{\sigma}_{h}^{k}$ satisfying Definition 4.2. The idea is to
reconstruct this equilibrated flux as the sum of a discretization flux
$\boldsymbol{\sigma}_{h,\textnormal{disc}}^{k}$ and the remainder flux
$\boldsymbol{\sigma}_{h,\textnormal{rem}}^{k}$. In this section we assume the
source term $f$ is smooth enough to calculate the discrete residual for single
species, $\mathbf{r}^{k}_{h,i}$ explicitly, i.e., for $1\leq i\leq s$, $1\leq
k\leq M$, let

(61) $\displaystyle\mathbf{r}^{k}_{h,i}$
$\displaystyle:=\frac{1}{s}(\tau_{k}f^{k}_{h}+s_{h}^{k-1})-(1-\tau_{k}r_{h}^{k})s^{k}_{h,i}-\tau_{k}\Phi_{h}^{k}(s_{h}^{k})s^{k}_{h,i}+\tau_{k}\nabla\cdot(g_{h}^{k}\nabla
s^{k}_{h,i}),$

with the global residual defined by
$\mathbf{r}^{k}_{h}:=\sum_{i=1}^{s}\mathbf{r}^{k}_{h,i}$. The equilibrated
discretization flux for species $i$ is then given by

(62) $\boldsymbol{\sigma}^{k}_{h,\textnormal{disc},i}:=-g_{h}^{k}\nabla
s^{k}_{h,i}.$

The construction of the micro-macro misfit (remainder) flux for species $i$ is
given component-wise by

(63)
$\displaystyle\boldsymbol{\sigma}^{k}_{h,\textnormal{rem},i}:=\frac{1}{d\tau_{k}}\Bigg{[}\int_{x_{i,1}}^{x_{1}}\mathbf{r}^{k}_{h,i}|_{(\zeta,x_{2},\cdots,x_{d})}\textnormal{d}\zeta,\cdots,\int_{x_{i,d}}^{x_{d}}\mathbf{r}^{k}_{h,i}|_{(x_{1},\cdots,x_{d-1},\zeta)}\textnormal{d}\zeta\Bigg{]}^{\top},$

Then, for $1\leq k\leq M$, we let

(64)
$\boldsymbol{\sigma}^{k}_{h}:=\boldsymbol{\sigma}^{k}_{h,\textnormal{disc}}+\boldsymbol{\sigma}^{k}_{h,\textnormal{rem}}:=\sum_{i=1}^{s}\\{\boldsymbol{\sigma}^{k}_{h,\textnormal{disc},i}+\boldsymbol{\sigma}^{k}_{h,\textnormal{rem},i}\\}.$

The following key result shows that $\boldsymbol{\sigma}_{h}^{k}$ from the
above definition leads to an equilibrated flux in the sense of Definition 4.2.

######  Proposition 6.1.

(Flux equilibration) Let the flux reconstruction
$\boldsymbol{\sigma}_{h}^{k}$, be defined by (64), where
$\boldsymbol{\sigma}^{k}_{h,\textnormal{disc},i}$ is defined by (62) and
$\boldsymbol{\sigma}^{k}_{h,\textnormal{rem},i}$ by (63). Then
$\boldsymbol{\sigma}_{h}^{k}\in\mathbf{H}({\textnormal{div}},\Omega)$ and we
have the flux equilibration property (40b) satisfied in the strong sense,
i.e.,

(65)
$\displaystyle(1-\tau_{k}r^{k}_{h})s_{h}^{k}+\tau_{k}\Phi^{k}_{h}(s_{h}^{k})s_{h}^{k}+\tau_{k}\nabla\cdot\boldsymbol{\sigma}_{h}^{k}=\tau_{k}f_{h}^{k}+s_{h}^{k-1}.$

######  Proof 6.2.

By construction, from relations (60) and (61), we have

(66a)
$\displaystyle\tau_{k}\nabla\cdot\boldsymbol{\sigma}^{k}_{h,\textnormal{disc},i}$
$\displaystyle=\frac{1}{s}(\tau_{k}f^{k}_{h}+s_{h}^{k-1})-(1-\tau_{k}r_{h}^{k})s^{k}_{h,i}-\tau_{k}\Phi_{h}^{k}(s_{h}^{k})s^{k}_{h,i}-\mathbf{r}^{k}_{h,i},$
(66b)
$\displaystyle\tau_{k}\nabla\cdot\boldsymbol{\sigma}^{k}_{h,\textnormal{rem},i}$
$\displaystyle=\mathbf{r}^{k}_{h,i}.$

Thus, from (64) and (60) we obtain

$\displaystyle\tau_{k}\nabla\cdot\boldsymbol{\sigma}_{h}^{k}$
$\displaystyle=\tau_{k}f^{k}_{h}+s_{h}^{k-1}-\sum_{i=1}^{s}\left\\{(1-\tau_{k}r_{h}^{k})s^{k}_{h,i}+\tau_{k}\Phi_{h}^{k}(s_{h}^{k})s^{k}_{h,i}\right\\}$
(67)
$\displaystyle=\tau_{k}f^{k}_{h}+s_{h}^{k-1}-(1-\tau_{k}r_{h}^{k})s^{k}_{h}-\tau_{k}\Phi_{h}^{k}(s_{h}^{k})s^{k}_{h}\in
L^{2}(\Omega).$

######  Remark 6.3 (Equilibrated flux).

The choice of the equilibrated discretization flux in (62) is motivated by the
fact that we wish to be as general as possible with regards to from where
$s_{h,i}^{k}$ is obtained. In practice, this means we may overestimate the
error. In principle, one could solve an optimization problem to get a better
estimate (see e.g., [2]), but then additional requirements on $s_{h,i}^{k}$ is
needed.

###  6.2 Distinguishing the error components

The preceding developments lead to the following result.

######  Theorem 6.4 (Error components).

For the species $i\leq s$, let $s_{k,i}^{k}$ be the population density as
given in (59) and the equilibrated flux as characterized in (65). We have the
following a posteriori error bound distinguishing the error components For
$1\leq k\leq M$, the following estimate holds true

(68)
$\displaystyle\left\\{{|\kern-1.07639pt|\kern-1.07639pt|n^{k}-s_{h}^{k}|\kern-1.07639pt|\kern-1.07639pt|}^{2}_{k}+\mathcal{J}^{k}(n^{k},s_{h}^{k})\right\\}^{\frac{1}{2}}\leq\sum_{i=1}^{s}(\eta_{\textnormal{rem},i}^{k}+\eta_{r,i}^{k}+\eta_{g,i}^{k}+\eta_{\Phi,i}^{k}+\frac{1}{s}\eta_{f}^{k}),$

where

(69)
$\eta_{\star,i}^{k}=\left\\{\sum_{K\in\mathcal{Z}_{h}}\left(\eta_{\star,K,i}^{k}\right)^{2}\right\\}^{\frac{1}{2}},$

with the species-dependent counterparts of the estimators (41) defined by

(70a) $\displaystyle\eta_{r,K,i}^{k}$
$\displaystyle:=\tau_{k}\tilde{\omega}^{k}_{K}\lVert(r^{k}-r_{h}^{k})s_{h,i}^{k}\rVert_{K},$
(70b) $\displaystyle\eta_{g,K,i}^{k}$
$\displaystyle:=\tau_{k}G_{k}^{-\frac{1}{2}}\lVert(g^{k}-g_{h}^{k})\nabla
s_{h,i}^{k}\rVert_{K},$ (70c) $\displaystyle\eta_{\Phi,i}^{k}$
$\displaystyle:=\tau_{k}{|\kern-1.07639pt|\kern-1.07639pt|(\Phi^{k}(s_{h,i}^{k})-\Phi_{h}^{k}(s_{h,i}^{k}))s_{h,i}^{k}|\kern-1.07639pt|\kern-1.07639pt|}_{k^{\ast}},$
and the modeling-remainder estimator defined by (70d)
$\eta_{\textnormal{rem},K,i}^{k}:=\tau_{k}G_{k}^{-\frac{1}{2}}\lVert\boldsymbol{\sigma}^{k}_{h,\textnormal{rem},i}\rVert_{K}.$

######  Proof 6.5.

First, observe that $\eta_{\textnormal{R},K}^{k}=0$ since the flux satisfies
the strong equilibration property (65). Next, substitute (60) and (64) in
(43). Due to (62), this gives

$\eta_{\textnormal{DF},K}^{k}=\tau_{k}G_{k}^{-\frac{1}{2}}\lVert\sum_{i=1}^{s}\boldsymbol{\sigma}^{k}_{h,\textnormal{rem},i}\rVert_{K}.$

Finally, use the triangle inequality to separate error components for
individual species to arrive at (68).

###  6.3 Detection of speciation events

Using Theorem 6.4 we can now detect speciation events in the SLM as follows:
At each time step $k>1$, we reconstruct the density function of each species
in the system using the reconstruction operator, $\mathcal{D}$, and compute
the a posteriori error bound (68). Then, if the estimate exceeds some given
tolerance we infer that a speciation event is about to happen for species $i$.

##  7 Numerical examples

In this section we present two numerical examples where we employ both the
heuristic and multi-scale methods. The examples are chosen such that the first
example conforms Assumption 2, and thus all the theoretical results apply. The
first example thus allows us to validate the applicability of the methodology.
The second example is motivated by the biological setting of speciation in a
predator-prey setting, for which Assumption 2 do not hold. Most notably, $d<3$
in addition to the operator $\Phi$ not satisfying the monotonicity condition
(19c), thus invalidating the assertion of Lemma 3.4. This example provides
numerical evidence to the efficacy of the multi-scale algorithm outside of the
parameter space where we have been able to prove error bounds. In this
context, the bound (68) must be considered more loosely as an error indicator.

For both examples we assess the accuracy of the speciation methods by solving
simultaneously the PLM globally and calculating the statistical moments which
are then compared to the corresponding species parameters. Moments of the
reference PLM solution are always shown as dotted lines, and indicated by the
letter $\mu$ in the legends. For the multi-scale method, the start and end of
the speciation event is indicated by two vertical dotted lines. For the
heuristic method, speciation is indicated by a single vertical dotted line.

We assume the discretization errors associated with solving the models are
negligible compared to the modeling error which we are interested in, hence in
practice we calculate only the remainder estimator
$\eta_{\textnormal{rem},i}^{k}$ from (70d) when estimating (68). Moreover, the
density reconstruction operator $\mathcal{D}$ is implemented as a normal
distribution, i.e., for a species
$N_{i}^{k}=(n_{i}^{k},x_{i}^{k},\upsilon_{i}^{k})$ we have

(71)
$(\mathcal{D}N_{i}^{k})(x):=\frac{1}{\Lambda(\upsilon_{i}^{k})}\exp\left(-\frac{1}{2}(x-x_{i}^{k})^{\top}(\upsilon_{i}^{k})^{-1}(x-x_{i}^{k})\right).$

Note that this does not satisfy the homogenous Dirichlet boundary conditions.
However, in practice for a large enough domain, boundary conditions are still
satisfied within working precision. An alternative approach would be to
truncate the tails of the normal distribution, but for the problems considered
herein, we do not expect this to make any difference in the results. Finally,
the $m$-species compression operator, $\mathcal{C}^{m}_{B_{i,\nu}^{k}}$, is
implemented using a nonlinear least squares iteration (here for $m=2$), where
we iterate until convergence with a relative tolerance of $1e-3$.

###  7.1 Numerical approximations

We advance in time using a fourth order Runge-Kutta scheme (RK4). Numerical
integration is by the midpoint rule, and the differential term of the PLM is
approximated by the two-point flux approximation (TPFA) method on a regular
Cartesian grid. All numerical examples are implemented in MatLab v. R2019b.

###  7.2 Example 1

For the first numerical example we let $d=3$, and consider a system initially
consisting of $s=1$ species, and where one speciation event occurs. Here, we
let the micro and macro time increments be given by $\tau_{m}=\tau_{M}=5e-2$,
respectively, and set $T=600$ as the final time. Spatial grid size is $\Delta
x=(1,1,1)/20$. The residual relative tolerance is chosen as
$\textnormal{TOL}_{\textnormal{res}}=5e1$, and the number of standard
deviations for the trait space density regions (and tolerance distance between
mean trait coordinates during speciation) as $\nu=10$. When the computed error
bound exceeds the tolerance, we backtrack 100 time units before initiating the
multi-scale/heuristic algorithms. Furthermore, since the Assumption 2 is
fulfilled for this example, we calculate the difference between the global PLM
solution and the reconstructed SLM solution in the energy norm, i.e.,
${|\kern-1.07639pt|\kern-1.07639pt|n^{k}-s^{k}|\kern-1.07639pt|\kern-1.07639pt|}_{k}$.
The trait space domain is the unit cube, i.e., $\Omega=[0,1]^{3}$. Initial
data is given by $N_{0}=(n_{0},x_{0},\upsilon_{0})$, where

(72) $n_{0}=2e-1,\quad x_{0}=(0.2,0.2,0.2),\quad\textnormal{ and
}\quad\upsilon_{0}=5e-3\times{\mathbb{I}}.$

####  7.2.1 Parameters

We impose a speciation event on the system by having a time dependent growth-
rate, where a single attractor point in trait space gradually transitions into
two attractor points. In particular, for $\gamma>0$ the growth rate is defined
as

(73) $r(x,t):=1-\gamma(f_{0}(t)r_{0}(x)+f_{\infty}(t)r_{\infty}(x)),$

where

(74a) $\displaystyle r_{0}(x)$ $\displaystyle:=\lVert
x-(0.5,0.5,0.5)\rVert^{2},$ (74b) $\displaystyle r_{\infty}(x)$
$\displaystyle:=\lVert x-(0.2,0.8,0.8)\rVert^{2}\lVert
x-(0.8,0.2,0.2)\rVert^{2},$

and where (for $\theta>0$)

(75) $f_{0}(t):=\begin{cases}1-t/\theta,&t<\theta,\\\
0,&t\geq\theta,\end{cases}\quad\textnormal{ and }\quad
f_{\infty}(t):=\begin{cases}t/\theta,&t<\theta,\\\ 1,&t\geq\theta.\end{cases}$

Thus, $\gamma$ determines the speed at which the species moves towards the
attractor points set by (74a) or (74b), and $\theta$ is the transition time
between $r_{0}$ and $r_{\infty}$. For the present situation we choose the
following parameter values

(76) $\gamma=2,\quad\textnormal{ and }\quad\theta=300.$

The remaining coefficients are defined by constant values, i.e.,

(77a) $\displaystyle b(x,t)$ $\displaystyle=1e-3,$ (77b)
$\displaystyle\alpha(x,y,t)$ $\displaystyle=-1,$ (77c) $\displaystyle g(x,t)$
$\displaystyle=5e-6\times{\mathbb{I}}.$

####  7.2.2 Simulation

With the parameters given in the previous section, we employ both the
heuristic and multi-scale speciation algorithms. Figure 3 below shows the
species abundance, figure 4 the species mean trait coordinates, figure 5 the
largest eigenvalue of the trait covariance matrix, figure 6 the a posteriori
error bound, and figure 7 the error in energy norm, as functions of time. Note
that due to the presence of $G^{-\frac{1}{2}}$ in the definition of
$\eta_{rem,i}^{k}$ (eq. (70d)), the magnitudes in figure 6 can not be directly
compared to those of the preceding figures.

We observe from the results a very close match between the reference solution
and the multi-scale method throughout the simulation time. Indeed, the error
in the energy norm, as seen in Figure 7 is lower during the multi-scale window
than during the pure SLM simulation, indicating that the modeling error of the
SLM in terms of capturing species dynamics dominates over the error associated
with the speciation event. In terms of species-level parameters, as shown in
the remaining figures, the qualitative match is also quite satisfactory for
the multi-scale mathod outside the speciation event (these quantities are of
course not defined during the speciation event itself).

In contrast, while the heuristic method is somewhat acceptable in terms of
capturing the initial dynamics and final state, it does not capture the
speciation dynamics themselves as accurately as the multi-scale method, thus
emphasizing the value of the multi-scale simulation framework. The lack of
accuracy during the speciation event leads to later errors in the higher
moments of the solution (i.e., velocities and covariances in trait space), as
seen in particular in Figures 4 and 5.

(a) Multi-scale method.

(b) Heuristic method.

Figure 3: Species abundance as functions of time.

(a) Multi-scale method.

(b) Heuristic method.

Figure 4: Species mean traits as functions of time (first component).

(a) Multi-scale method.

(b) Heuristic method.

Figure 5: Maximum eigenvalue of trait covariance matrix as functions of time.

(a) Multi-scale method.

(b) Heuristic method.

Figure 6: A posteriori modeling-remainder estimator,
$\eta_{\textnormal{rem},i}^{k}$, as a function of time.

(a) Multi-scale method.

(b) Heuristic method.

Figure 7: Global error measured in energy norm as function of time.

###  7.3 Example 2

For the second numerical example we let $d=2$, and consider a predator-prey
system initially consisting of $s=2$ species, and where the ‘prey’-species
undergoes a speciation event by traveling along a Y-shaped ridge in trait
space. Hence, the divergence will happen at the branching point of the ridge.
Here, we let the micro and macro time increments be given by $\tau_{m}=1e-3$
and $\tau_{M}=1e-2$, respectively, and set $T=600$ as the final time. Spatial
grid size is $\Delta x=(1,1)/50$. The residual relative tolerance is chosen as
$\textnormal{TOL}_{\textnormal{res}}=1e1$, and the number of standard
deviations for the trait space density regions as $\nu=10$. When the computed
error bound exceeds the tolerance, we backtrack 250 time units before
initiating the multi-scale/heuristic algorithms. The trait-space domain
consists of two disjoint regions; $\Omega_{1}=[0,1]^{2}$ and
$\Omega_{2}=[2,3]^{2}$, i.e., $\Omega=\Omega_{1}\cup\Omega_{2}$, wherein the
two species are located, respectively. Initial data is given by
$N_{0}=(N_{\textnormal{prey},0},N_{\textnormal{pred},0})$, where

(78a) $\displaystyle n_{\textnormal{prey},0}=2e-1,\quad
x_{\textnormal{prey},0}=(0.5,0.3),\quad\textnormal{ and
}\quad\upsilon_{\textnormal{prey},0}=5e-3\times{\mathbb{I}},$ (78b)
$\displaystyle n_{\textnormal{pred},0}=2e-1,\quad
x_{\textnormal{pred},0}=(2.5,2.5),\quad\textnormal{ and
}\quad\upsilon_{\textnormal{pred},0}=5e-3\times{\mathbb{I}}.$

####  7.3.1 Parameters

We impose a speciation event upon the ‘prey’-species by initializing it at the
foot of a ridge in trait space, and as the species travels along this ridge,
the ridge splits into two branches. This branching ridge is incorporated in
the growth-rate coefficient, which for $c_{0},\delta,r_{0}>0$ and
$\mathbf{c}\in{\mathbb{R}}^{2}$, is defined as

(79) $r(x):=\begin{cases}c_{0}+\delta\mathbf{c}\cdot
x^{\top}-\varphi_{\epsilon}(x)\ast\textnormal{dist}(x,Y),\quad&x\in\Omega_{1},\\\
-r_{0},&x\in\Omega_{2},\end{cases}$

where the set of points $Y\subset\Omega_{1}$, is the three line segments
connecting the nodes $\\{(0.5,0.3),(0.5,0.5),(0.2,0.7),(0.8,0.7)\\}$ to form a
Y-shape, and where $\varphi_{\epsilon}$ is the mollifier function centered on
the origin with radius of support $\epsilon>0$, and where
$\textnormal{dist}(x,Y)$ is the distance from the point $x$ to the set $Y$.
Hence, $c_{0}$ is the growth-rate of the prey, $\delta$ is the speed at which
it travels along the ridge, $\mathbf{c}$ is a direction vector, $\epsilon$ is
the steepness of the ridge, and $r_{0}$ is the loss rate of the predator. The
interaction coefficient is defined by

(80) $\alpha(x,y):=\begin{cases}0,\quad&x\in\Omega_{1},y\in\Omega_{1},\\\
-\gamma,\quad&x\in\Omega_{1},y\in\Omega_{2},\\\
\beta,\quad&x\in\Omega_{2},y\in\Omega_{1},\\\
0,\quad&x\in\Omega_{2},y\in\Omega_{2},\end{cases}$

where $\gamma>0$ is the rate of predation upon the prey, and where $\beta>0$
is the growth rate of the predator. For the present situation we choose the
following parameter values

$\displaystyle c_{0}=0,\quad\delta$ $\displaystyle=0.8,\quad
r_{0}=0.5,\quad\epsilon=0.2,\quad\gamma=3.0,$ (81) $\displaystyle\beta$
$\displaystyle=8.0,\text{ and }\mathbf{c}=(0,1)^{\top}.$

The remaining coefficients are defined by constant values, i.e.,

(82a) $\displaystyle b(x,t)$ $\displaystyle:=0,$ (82b) $\displaystyle g(x,t)$
$\displaystyle:=2e-6\times{\mathbb{I}}.$

####  7.3.2 Simulation

With the parameters given in the previous section, we employ both the
heuristic and multi-scale algorithms. Figures 8–9 below shows the species
abundance, figures 10–11 the species mean trait coordinates, figure 12 the
largest eigenvalue of the trait covariance matrix, and figures 13–14 the a
posteriori error bounds, as functions of time.

We recall that for this example, Assumption 2 is not satisfied, and the
results are thus not expected to be as strong as for the previous example.
Indeed, we observe that while the multi-scale method performs fairly overall,
the a posteriori error bound, now being only an error indicator, is not
precise enough to identify the speciation event early enough, necessitating a
larger backtrack window than the previous example. Note also that the
abundance plots, figures 8–9, indicate that the predator-prey cycles have been
shifted out of phase, but the correct structure is still retained. However,
the multi-scale method still clearly outperforms the heuristic method, which
due to the reliance on the same error indicator also suffers from a somewhat
delayed speciation event. On the other hand, the heuristic method models the
speciation event less accurately than the multi-scale method, and thus the
errors after the speciation event are significantly larger when seen in terms
of the species-level parameters.

(a) Multi-scale method.

(b) Heuristic method.

Figure 8: Prey abundance as function of time. Curves for child-species are
overlapping.

(a) Multi-scale method.

(b) Heuristic method.

Figure 9: Predator abundance as function of time.

(a) Multi-scale method.

(b) Heuristic method.

Figure 10: Prey mean traits coordinate as functions of time (first component).

(a) Multi-scale method.

(b) Heuristic method.

Figure 11: Prey mean traits coordinate as functions of time (second
component). Curves for child-species are overlapping

(a) Multi-scale method.

(b) Heuristic method.

Figure 12: Prey maximum eigenvalue of trait covariance matrix as functions of
time. Curves for child-species are overlapping.

(a) Multi-scale method.

(b) Heuristic method.

Figure 13: A posteriori modeling-remainder estimator,
$\eta_{\textnormal{rem},i}^{k}$, for ‘prey’-species as functions of time.

(a) Multi-scale method.

(b) Heuristic method.

Figure 14: A posteriori modeling-remainder estimator,
$\eta_{\textnormal{rem},i}^{k}$, for ‘predator’-species as a function of time.

######  Remark 7.1 (Computation times).

The heuristic method will have a computation time approximately equal to that
of solving only the species level model for the same time interval. On the
other hand, the multi-scale method will have a computation time higher than
that of solving the population level model for a time interval equal to the
speciation interval (with the same grid size as the local density region)
since the compression operator is implemented as an iteration procedure. Thus,
the heuristic method will in general be orders of magnitude faster than the
multi-scale method.

##  8 Conclusions

We have developed two strategies for modeling speciation
events within the context of species interaction models. The first, _heuristic
approach_ , is based on splitting the diverging species according to the
spectrum of the trait covariance matrix. The second, _multi-scale approach_ ,
is based on resolving the diverging species as a population density
distribution using a fine-scale population level model for the duration of the
speciation event (i.e., until ‘child’-species are sufficiently separated in
trait space). Crucial to both these approaches is the connection between the
species scale and the population scale, i.e., the ability to view the species
either as an abundance-trait-covariance tuple, or as a population density
distribution. This allows for defining a multi-scale framework in which these
two scales are coupled, and to calculate the a posteriori error bound of the
reconstructed macro-scale solution which then indicates the modeling error. We
have also given conditions on the nonlinearities of the micro-scale model for
which well-posedness of the time-discrete problem is guaranteed. Using
explicit equilibrated flux and density reconstructions, we presented a
posteriori error estimates for an error measure composed of an energy
$H^{1}$-norm and a semi-metric in terms of a residual monotone operator. In
particular, the dual norm of the residual is found to be equal to the error
between the exact and approximate solutions, again given conditions on the
nonlinearities. Even when the theoretical conditions are not fulfilled, our
framework provides a working algorithm in practice, as our second numerical
example shows. Finally, regarding the heuristic and multi-scale methods, by
comparison of the species parameters with the corresponding moments from the
reference (global) PLM solution, it is clear that the multi-scale approach is
superior to the heuristic approach. In fact, from our experiments it appears
that a multi-scale approach to speciation is indeed required for eco-
evolutionary modeling at the species level, in which speciation events are
allowed. We propose that this multi-scale approach might serve as a productive
way of integrating ecological processes and evolutionary processes.

## Acknowledgments

The research was funded in part by Norwegian Research Council project no.
263149.

## References

  * [1] G. V. Alekseev, R. V. Brizitskii, and Z. Y. Saritskaya, Stability estimates of solutions to extremal problems for a nonlinear convection-diffusion-reaction equation, Journal of Applied and Industrial Mathematics, 10 (2016), pp. 155–167.
  * [2] R. Becker and R. Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta numerica, 10 (2001), pp. 1–102.
  * [3] W. Cheney, Analysis for applied mathematics, vol. 208 of Graduate Texts in Mathematics, Springer-Verlag, New York, 2001, https://doi.org/10.1007/978-1-4757-3559-8, https://doi.org/10.1007/978-1-4757-3559-8.
  * [4] J. A. Coyne, H. A. Orr, et al., Speciation, vol. 37, Sinauer Associates Sunderland, MA, 2004.
  * [5] F. Débarre, S. Nuismer, and M. Doebeli, Multidimensional (co) evolutionary stability, The American Naturalist, 184 (2014), pp. 158–171.
  * [6] U. Dieckmann, Can adaptive dynamics invade?, Trends in Ecology & Evolution, 12 (1997), pp. 128–131.
  * [7] U. Dieckmann, Å. Brännström, R. HilleRisLambers, and H. C. Ito, The adaptive dynamics of community structure, in Mathematics for ecology and environmental sciences, Springer, 2007, pp. 145–177.
  * [8] U. Dieckmann and M. Doebeli, On the origin of species by sympatric speciation, Nature, 400 (1999), pp. 354–357.
  * [9] U. Dieckmann, M. Doebeli, J. A. Metz, and D. Tautz, Adaptive speciation, Cambridge University Press, 2004.
  * [10] M. Doebeli, A quantitative genetic competition model for sympatric speciation, Journal of evolutionary biology, 9 (1996), pp. 893–909.
  * [11] M. Doebeli and U. Dieckmann, Evolutionary branching and sympatric speciation caused by different types of ecological interactions, The american naturalist, 156 (2000), pp. S77–S101.
  * [12] M. Doebeli and U. Dieckmann, Speciation along environmental gradients, Nature, 421 (2003), pp. 259–264.
  * [13] A. Ern, I. Smears, and M. Vohralík, Guaranteed, locally space-time efficient, and polynomial-degree robust a posteriori error estimates for high-order discretizations of parabolic problems, SIAM Journal on Numerical Analysis, 55 (2017), pp. 2811–2834.
  * [14] A. Ern and M. Vohralík, A posteriori error estimation based on potential and flux reconstruction for the heat equation, SIAM Journal on Numerical Analysis, 48 (2010), pp. 198–223.
  * [15] L. C. Evans, Partial differential equations, Graduate studies in mathematics, 19 (1998).
  * [16] S. Gavrilets, Fitness landscapes and the origin of species (MPB-41), Princeton University Press, 2004.
  * [17] S. A. Geritz and K. Éva, Adaptive dynamics in diploid, sexual populations and the evolution of reproductive isolation, Proceedings of the Royal Society of London. Series B: Biological Sciences, 267 (2000), pp. 1671–1678.
  * [18] E. R. Leigh, The ecological role of volterra’s equations, Some mathematical problems in biology, (1968).
  * [19] J. M. Nordbotten, F. Bokma, J. S. Hermansen, and N. C. Stenseth, The dynamics of trait variance in multi-species communities, Royal Society open science, 7 (2020), p. 200321.
  * [20] J. M. Nordbotten and N. C. Stenseth, Asymmetric ecological conditions favor red-queen type of continued evolution over stasis, Proceedings of the National Academy of Sciences, 113 (2016), pp. 1847–1852.
  * [21] P. Nosil, Ecological speciation, Oxford University Press, 2012.
  * [22] J. Papež, Z. Strakoš, and M. Vohralík, Estimating and localizing the algebraic and total numerical errors using flux reconstructions, Numerische Mathematik, 138 (2018), pp. 681–721.
  * [23] T. Price et al., Speciation in birds, Roberts and Co., 2008.
  * [24] M. L. Rosenzweig, Competitive speciation, Biological Journal of the Linnean Society, 10 (1978), pp. 275–289.
  * [25] A. Sasaki and U. Dieckmann, Oligomorphic dynamics for analyzing the quantitative genetics of adaptive speciation, Journal of mathematical biology, 63 (2011), pp. 601–635.
  * [26] I. Smears and M. Vohralík, Simple and robust equilibrated flux a posteriori estimates for singularly perturbed reaction–diffusion problems, ESAIM: Mathematical Modelling and Numerical Analysis, 54 (2020), pp. 1951–1973.
  * [27] M. Vohralík, A posteriori error estimates for lowest-order mixed finite element discretizations of convection-diffusion-reaction equations, SIAM J. Numer. Anal., 45 (2007), pp. 1570–1599, https://doi.org/10.1137/060653184, https://doi.org/10.1137/060653184.
  * [28] A. Vukics, J. Asbóth, and G. Meszéna, Speciation in multidimensional evolutionary space, Physical Review E, 68 (2003), p. 041903.
  * [29] E. Weinan, Principles of multiscale modeling, Cambridge University Press, 2011.
  * [30] C. S. Woodward and C. N. Dawson, Analysis of expanded mixed finite element methods for a nonlinear parabolic equation modeling flow into variably saturated porous media, SIAM Journal on Numerical Analysis, 37 (2000), pp. 701–724.