Geometric Numerical Integration in Ecological Modelling

A major neglected weakness of many ecological models is the numerical method used to solve the governing systems of differential equations. Indeed, the discrete dynamics described by numerical integrators can provide spurious solution of the corresponding continuous model. The approach represented by the geometric numerical integration, by preserving qualitative properties of the solution, leads to improved numerical behaviour expecially in the long-time integration. Positivity of the phase space, Poisson structure of the flows, conservation of invariants that characterize the continuous ecological models are some of the qualitative characteristics well reproduced by geometric numerical integrators. In this paper we review the benefits induced by the use of geometric numerical integrators for some ecological differential models.


Introduction
Ecological modelling based on non linear differential systems of equations can be divided into two main categories [1]. Predictive modelling aims to reproduce, based on the knowledge provided by historical time series, the real observed phenomena and to predict their state in the future. Due to the large amount of uncertainty contained in ecological data, however, the approach of conceptual modelling is an effective alternative aimed to understand the main features of the ecological dynamics by making scenario analysis [2]. In both approaches, the mathematical model is based on governing laws described by non linear differential systems of equations. Their correct solution requires robust and accurate numerical algorithms; nevertheless, numerical schemes for ecological models have received little attention in the literature as most descriptions of models outline the governing equations but do not discuss their numerical solution. Indeed, it seems reasonable to assume that, within a predictive approach, numerical errors are always overwhelmed by uncertainties in the data and governing equations and to pay attention to numerical robustness is, therefore, unwarranted. On the other side, for conceptual modelling, first-order explicit numerical procedures are generally used in the belief that they are sufficient to infer potential future scenarios.
However, in literature, a rigorous numerical analysis is a recognized need due to the often neglected fact that the dynamics of the numerical flow can differ significantly from that of the original differential system: • 'An inadequate choice of a numerical method may have a detrimental effect on the study by making simulations costly and even producing wrong results' [1]; • 'Erroneous or misleading conclusions of model analysis and prediction arising from numerical artifacts in hydrological models are intolerable' [3]; • 'The literature abounds with examples of spurious behavior of numerical solutions, e.g., spurious fixed points, numerical 'chaos', and spurious periodicity' [4].
In the last decades, the theory of numerical methods allowed excellent general-purpose codes, mainly based on Runge Kutta methods or linear multistep methods. The born of the geometric numerical integration operated in the literature a shift of view-point: a numerical method is viewed as a discrete dynamical system which approximates the flow of the differential equation while preserving some of its specific properties. These are crucial for a qualitatively correct simulation of the dynamics and always improve long-time numerical behaviour [5].
In ecological modelling, the approximated solutions should be able to reproduce the main physical qualitative characteristics of the observed quantities (e.g., positive concentrations, mass/energy conservation) in order to make accurate previsions or outline realistic scenarios. For this reason, geometric numerical integration plays a crucial role in the analysis of ecosystem models described by systems of differential equations. The principal flow structures that arise in ecological modelling and that have to be preserved are Poisson maps. By means of the Darboux-Lie theorem, it is possible to find a transformation such that in the novel coordinates, a Poisson problem assumes a canonical Hamiltonian form. Then, symplectic methods, as for example collocation methods [6], are convenient for the structure preserving long-time integration. In this regard, as the order of the error of the collocation method is basically the same as of the underlying quadrature rule, one can exploit the Gaussian splines rules [7][8][9], for numerical integration to keep the same order of error, while requiring fewer quadrature points [5]. Instead of trasforming Poisson maps in Hamiltonian maps, in this review we will describe discrete dynamical systems that directly preserve the Poisson structure. They will correspond to Poisson numerical integrators of Poisson ecological systems. We will show their properties and the gain in performance with respect to classical procedures.
Other aspects that, in the nonlinear context, are alternative to the preservation of the structure of the flow, but not less important, are the phase-space preservation [10] and the conservation of invariants (see, for example, [11,12]). Indeed, since in computational ecology all involved quantities (i.e., densities, masses or concentrations) should be nonnegative, is of outmost importance to preserve the positive phase-space of the dynamics in order to produce physically valid numerical approximations. Moreover, approximations that preserve linear invariants are able to provide more accurate description of biogeochemical processes in ecosystem models [13].
In this paper, we will describe the application and the results of the analysis of geometric numerical integration in the framework of ecological modeling. For both predictive and conceptual real world modelling, the use of geometric numerical integration was indeed the favourite tool for making accurate qualitative scenario analysis and saving computational time [14,15]. Limiting to first-order procedures, we will illustrate the most suitable geometric numerical integrators for some ecosystem dynamics models. For optimal control problems that model management policies, we will describe the advantage of the symplectic Runge-Kutta pairs in order to get efficient methods easy to be implemented within any numerical software [16]. Finally, as future perspective, we will present some examples of ecological models featured by both Poisson and biochemical structure [17] and we will discuss about some open questions related to their numerical approximation.

Poisson Integrators for Poisson Systems in Ecology
Ecosystem population dynamics, including organism such as algae, invertebrates, fish, large herbivores and carnivores, are characterized by the interaction between predator and prey populations that controls and drives changes in populations over time. When resources are limited, individuals compete for access to resources, and populations decline. To survive and reproduce, they must obtain sufficient food while simultaneously avoiding becoming food for a predator. To address the trade-off between food intake and predator avoidance, ecologists have turned to mathematical models to better understand foraging behavior and predator-prey dynamics. Lotka-Volterra models provide the main tool to help population ecologists understand the factors that influence population dynamics. Although the models greatly simplify actual conditions, they demonstrate that under certain circumstances, predator and prey populations can oscillate over time in a manner similar to that observed in the real-word.
The Lotka-Volterra model falls in the more general class of ecological models called M-systems [18] for which population dynamics is described by Poisson systems analogous to the Hamiltonian formulation of classical mechanics. Other important examples of ecological differential models featured by a Poisson structure can be found among multispecies Lotka-Volterra systems [19]. The classical example is the Rock-Paper-scissor model based on relationships that characterize communities without strict competitive hierarchies [20].
We will see that the geometric numerical integration of ecological models featured by a Poisson dynamics, performed by means of Poisson integrators [5], especially over long time horizons, improves the qualitative results and, consequently the analysis of ecosystems scenarios.
Before entering into the details, we recall that, given a vector-valued scalar function H = H(y; t), called Hamiltonian function, a Poisson system is defined by the differential equation where B(y) is a matrix skew-symmetric operator and satisfies for all i, j, k ( the Jacobi identity). The flow ϕ t (y) of a Poisson-system is a Poisson map i.e., its Jacobian matrix φ t satisfies If the skew-symmetric operator B(y), which allows to define a Poisson bracket, is of constant rank n − q = 2r, then there exist q functions C 1 (y), C 2 (y), . . . C q (y) called Casimirs, which satisfy ∇C i (y) T B(y) = 0. Consequently, resulting Casimirs first integrals of Poisson system for all Hamiltonian function H.
In the next section, we provide two model examples representative of Poisson dynamics in ecology and we review the principal Poisson integrators developed in the recent literature for their numerical approximation. Moreover, we will illustrate the rock-scissor-paper model and the Hirota's algorithm for its numerical approximation.

Dynamical M-System
There exists a wide class of systems in ecology, that are called M-systems, for which population dynamics can be described by a phase-space theory analogous to the Hamiltonian formulation of classical mechanics [18]. An M-system is a dynamical system defined by a given resource function M (u, v; t) and equation of motions given bẏ By setting y = (u, v), it is easy to see that M-systems are special cases of Poisson systems dy dt = B(y)∇M with Hamiltonian function M (y; t) and skew-symmetric matrix operator which satisfies the Jacobi identity (1). By considering the change of variables:ũ = ln (u/σ),ṽ = ln (v/ξ) where σ, ξ > 0, the class of M-system on R 2 + has an associated Hamiltonian systeṁũ where H(ũ,ṽ; t) = M(u(ũ), v(ṽ); t) is the Hamiltonian function. Notice that, the relations comprise the positive Poisson phase plane R 2 + into the Hamiltonian phase plane R 2 ; conversely, the inverted relations comprise the Hamilton phase space R 2 into the positive Poisson phase plane R 2 + . The first example of a positive M-system is given by the Lotka-Volterra dynamics that, in adimensional variables is written aṡ In this case, the resource function M (u, v) is given by M(u, v) = −a ln v + v − b ln u + u (see also [5], Sec VII.2.2). It represents an example of separable M-system as the resource function is separable in two parts, one depending on u and the other on v i.e., The second example of positive M-system considers a two-species system described by the equationsu where the second equation defines the logistic evolution, γ and k are the growing rate for the population u and v respectively ; µ v represents the density dependent mortality rate, and u s is a saturation constant. The corresponding M -function is defined by Notice that, in this case, M is not separable in two parts each depending on only one variable.

Poisson Integrators for M-Systems
Given a dynamical M-system, a numerical one-step method y n+1 = Φ ∆t (y n ) is called Poisson integrator for a M -system if it is a Poisson map whenever applied to the M -system i.e., if the Jacobian Φ ∆t (y n ) satisfies Φ ∆t (y n ) B(y n )Φ ∆t (y n ) T = B(Φ ∆t (y n )).
where B (of full rank) is defined in Equation (2). The Symplectic Euler scheme is a Poisson integrator for separable M-system. It is implicit, it is not a Poisson integrator for general non-separable M-system and, in general, starting from positive values, it does not provide positive solutions without constraining the stepsize. This constraint, for Lotka-Volterra dynamics results stronger than the condition required by the linear stability analysis [21]. The explicit variant of symplectic Euler [22] u n+1 − u n ∆t = −u n v n ∂M ∂v is a Poisson integrator for a separable M-systems. As symplectic Euler is not a Poisson integrator for general non-separable M-systems and in general, starting from positive solution, it does not provide positive solutions. When applied to Lotka-Volterra dynamics, the condition for positivity for exponentially long-time intervals gained by means of backward error analysis, is stronger than conditions given by the linear stability analysis [21]. The main advantage is its explicit functional form. A Poisson integrator, suitable for positive M-systems, is Poisson Euler [5]: together with its adjoint. In the following, we recall the main properties proved in [23] Theorem 1. The first-order method (4) is a Poisson integrator for M-systems for any resource M-function such It provides positive solutions for all ∆t > 0 whenever v 0 , u 0 are positive. In case of separable M-systems it is also explicit.
, the numerical scheme (4), linearized around the equilibria, has the same stability properties of the Lotka-Volterra model i.e., it has a saddle point at the origin and a neutral center at the internal equilibrium (b, a).
Form the above properties, Poisson Euler provides positive solutions without constraints on the choice of the stepsize and has the same linear stability constraints of the Symplectic Euler method when applied to Lotka-Volterra dynamics. For separable M-systems it is explicit while it results implicit for non separable M-systems. This drawback is not overcome by the Explicit Variant of Poisson Euler scheme, given by as it results a discrete Poisson map only in case of separable M-system. To illustrate the performance of the methods, we approximate the solution of Equation (3) which represents a not separable M-system. In Figure 1 the gain in performance of the Poisson Euler schemes with respect to Symplectic Euler schemes is illustrated.

Positive Integrators for Generic Predator-Prey Dynamics
Consider a generic two-species predator-prey systeṁ The idea to perform the log transform, apply a symplectic Runge-Kutta scheme and then transform back by the exponential, allows to build numerical integrators able to provide positive approximations [23]. The first order variants of these schemes, that we continue to call Poisson Euler as it generalizes the Poisson Euler for M-dynamics, are given by: with its explicit variant its adjoint and its explicit variant Notice that, when f (u, v) = u F(v) then Poisson Euler method (5) results explicit; similarly, when g(u, v) = v G(u) the adjoint method (7) results explicit too. In the sequel, the performance of these positive schemes will be illustrated when applied to the approximation of the reaction semiflow of a spatially explicit dynamics.

Comparison Among Integrators
In Table 1 we summarize the properties of the main first-order geometric numerical integrators here proposed and we compare their computational complexity with respect to classical explicit and implicit Euler schemes. To generalize the proposed approaches to 2N dimensional interactions with N > 1, we refer to the solution of the ODE systeṁ where u, v, f, g ∈ R N . We leave out the costs due to the evaluation of N-dimensional vector functions and we account only for the solution of systems of equations. We distinguish the methods according to their functional form: The method considered are listed below.
• Explicit Euler (EE) The classic system, that involves a community of three competing species satisfying a relationship similar to the children's game Rock-Paper-scissors, has been considered in [20] to study how local dispersal promotes biodiversity for non-transitive communities, that is, those without strict competitive hierarchies. Such relationships, where rock crushes scissors, scissors cuts paper, and paper covers rock have been demonstrated in several natural systems (see for example [24,25]). The dynamical model is described by a three species Lotka-Volterra system: By setting y = (u, z, v), the three species model may be defined ad Poisson systems in two ways: and A numerical one-step method y n+1 = Φ ∆t (y n ) is a Poisson integrator for the above three species Lotka-Volterra system if the Jacobian Φ ∆t (y n ) is a Poisson map w.r.t. B 1 which preserves the Casimir H 1 or is a Poisson map w.r.t. B 2 which preserves the Casimir H 2 . Hirota's algorithm given in [26] for the above three-species example, is described by the following iterative scheme [19] : It results that it is a Poisson integrator for dy dt = B 2 ∇H 1 i.e., it is a Poisson map w.r-t B 2 which preserves the Casimir H 2 .
Another important geometric aspect of the Hirota's algorithm is the preservation of positivity. Indeed, starting from positive values, it provides positive solutions without constraints on the choice of stepsize. This can be shown by writing the Hirota's algorithm as:

Poisson Integrators for Spatially Extended Ecosystem Dynamics
In the previous section, we have recalled that the predator-prey interactions of different species, are extremely important in some temporal ecological problems. Although the space has been neglected so far, to describe the onset of important phenomena in some ecosystems dynamics, the spatial factors should be included in the model.

Poisson Euler for Turing Patterns Occurrence in a Ratio-Dependent Phytoplankton-Zooplankton Model
When space is included in the models, the formation of spatial patterns constitutes a very active research area in ecosystem modelling, based on the pioneering work of Turing [27]. There is a huge literature on the subject; however, to show the benefits of using Poisson integrators for predicting the onset of Turing patterns here we will focus on the specific model described in [28]. It is a three-chain coupled map lattice model built for exploring the spatiotemporal complexity of a predator-prey system with migration and diffusion. Based on Turing instability analysis, the authors derive pattern formation conditions for the predator-prey system. Via numerical simulation, Turing patterns can be found with subtle self organized structures under diffusion-driven and migration-driven mechanisms. It is evident that, if the aim of numerical simulation is to establish the formation of patterns, it is crucial to have numerical flows which share the same qualitative characteristics of the theoretical ones.
The model is a ratio-dependent phytoplankton-zooplankton model in a two dimensional space Ω = [0, L] 2 . Periodic boundary conditions are set and dynamics is described by Linear stability analysis established that (P 1 , Z 1 ) = (k, 0) is an unstable equilibrium of the local dynamics. Under the hypothesis that m z < 1 and r + k(m z − 1) > 0, the equilibrium (P 2 , i,j representing prey and predator density in the (i, j) spatial cell at the time t n = t 0 + n∆t, the discrete flow is split in advection-diffusion-reaction (ADR) semiflows: 3. discrete reaction semiflow: advances in time with the explicit variant of (the adjoint of) Poisson-Euler (8): In [28] the discrete reaction semiflow advances in time with explicit Euler. Here, we use Poisson Euler scheme to investigate the performance of a geometric numerical integrator.

Turing Instability Analysis of ADR Discrete Semiflows
We consider a spatially heterogeneous perturbation of the stable equilibrium (P 2 , Z 2 ) and we consider the linearized solution (P Given λ + , and λ − the eigenvalues ofJ, the occurrence condition for Turing instability can be written as is the Jacobian evaluated at (P 2 , Z 2 ) of the discrete reaction semiflows defined by the numerical methods. It is evident that, changing the numerical method used to discretize the reaction semiflow, changes the entries of the matrix A and, consequently, the conditions for Turing instability. Of course, letting both the stepsizes ∆t and ∆x go to zero, all the numerical methods should predict the occurrence or non-occurence of Turing instability for the continuous model (9). In Figure 2, in correspondence of parameters k P = 0.02, k Z = 0.2, k = 1.15, m z = 0.6, r = 0.5 and T = 2250, L = 100 we show the phytoplankton distribution P at T = 2250 of the discrete ADR model which proceeds with ∆t = 0.05 in the square [0, 100] 2 discretized with ∆x = 0.5: we can observe the transition from hot spots-stripes to banded patterns obtained by increasing the advective coefficient W g . In the same Figure we show that, in correspondence of that values of W g , there is always a non empty set of wavenumbers where the occurrence condition for Turing instability is satisfied.
In our simulations, we used Poisson Euler as discrete reaction semiflow; however, results remained unchanged if, in the third step of the algorithm described in Section 3.1, we use Explicit Euler in place of Poisson Euler to approximate the reaction step. Now, set the parameters k = 0.96 and W g = 0.01 and keep all the others unchanged. In Figure 3 we show the the Turing instability occurrence condition at T = 2250 of the discrete ADR model in the square [0, 100] 2 discretized with ∆x = 1 and Explicit Euler semiflow (left) that proceeds with three different values of temporal stepsize ∆t = 0.1, 0.5, 0.9. While with larger values of the time step the Turing instability occurrence condition is satisfied, with ∆t = 0.1 the condition is not satisfied. This means that, if we want to infer the behaviour of the continuous model (9) with the ADR model, the explicit Euler discrete semiflow can give a wrong prediction if the temporal stepsize is not chosen sufficiently small.  (d) With the same set of parameters, instead, Poisson Euler predicts the non occurrence of Turing instability even in correspondence of the largest stepsize ∆t = 0.9 (Figure 3b). This seems to indicate that using geometric numerical method to discretize a semiflow featured by a given structure (in this case the reaction semiflow featured by a Poisson structure) produces solutions qualitatively more similar to the solution of the theoretical flow.

Geometric Integration for the Spatiotemporal Dynamics of Aquatic Ecosystems
The spatiotemporal complexity of plankton and fish dynamics is explored in the seminal paper [29]. The authors show that the formation of a patchy spatial distribution of species can be described by a two-species prey-predator (i.e., Phytoplankton-Zooplankton) interaction modelled by the well-known Rosenzweig-MacArthur system [30], here given in non dimensional form: The populations are supposed distributed in an horizontal two-dimensional domain where zero-flux boundary conditions are set. The model (10) couples logistic prey growth with Holling II functional predator response and shows a very rich dynamics. Indeed, for a purely homogeneous initial distribution of species, the system stays homogeneous forever and no spatial pattern emerges; for a very weakly perturbed initial distribution, a smooth pattern arises that is not persistent and gradually evolves to the homogeneous distribution; for somehow stronger initial perturbations, the system evolves to the formation of a jagged spatial pattern which is persistent in time. The formation of the irregular patchy structure can be preceded by the evolution of a regular spiral spatial pattern.
Similarly to the well-known Rosenzweig-MacArthur model (10), the principal population dynamics models are based on logistic prey growth and 'Holling type' functional response of the predators [31]. They can exhibit spiral waves, target waves, and spatiotemporal chaos; for all those complex dynamics, the geometric numerical approach can represent a powerful tool to reproduce the correct behavior of the continuous solutions.

Positive Schemes for Spatially Extended Predator-Prey Dynamics
In this section, we focus on spatially-extended predator-prey models described by reaction-diffusion systems in the following general form where u(x, t) and v(x, t) represent population densities of prey and predators at time t and position x and D u and D v are positive constant diffusion coefficients. The equations evolve in Ω T := Ω × (0, T) where the domain Ω is a bounded and open subset of R d , d ≤ 3. The system is augmented with initial conditions and homogeneous zero-Neumann boundary conditions. To assure the non-negativity of solutions which correspond to biologically meaningful densities, the reaction kinetics have to satisfy [32] Consequently, if the initial data (u 0 (x), v 0 (x)) is chosen in [0, +∞) 2 for all x ∈ Ω, then by a maximum principle the solution (u(x, t), v(x, t)) also lies in [0, +∞) 2 , which is a positively invariant region for the system. For the numerical approximation, in [23] a positive first order accurate integrator is obtained by compositions of exact and numerical flow. Firstly, the spatial discretization of the Laplacian operators transforms the problem (11) and (12) into the system where The diffusive flow dy dt = L h y is supposed to be exactly computed being ϕ ∆t the exact flow.
The reaction flow dy dt = F h (y) is approximated by Poisson Euler schemes Φ ∆t in Equation (5), or the adjoint Φ [R] * ∆t in Equation (7). Then the methods: together with their adjoints ∆t .
are first order positive schemes, here denoted as EXponential-Positive Poisson (EXPP) schemes.
For example, EXPP schemes defined by Φ * ∆t = ϕ ∆t , advances in time according to the steps: for all x ∈ Ω h . The above method, when applied for solving Rosenzweig-MacArthur model (10) reads as follows: As we will see in the following, EXPP schemes outperform, in qualitative comparison, any existing scheme in detecting the onset of chaos predicted by model (10)  ( perfomed in [23] with polynomial Krylov approximation), together with a rigorous theoretical analysis, is still missing.

Implicit-Symplectic Schemes for Spatially Extended Predator-Prey Dynamics
Implicit-symplectic (IMSP) schemes are numerical integrators based on an implicit scheme for the stiff diffusive term and a geometric integrator for the reaction function. IMSP schemes were proposed in [33], as novel numerical schemes for the simulation of population and metapopulation predator-prey dynamics.
In order to illustrate the methods, consider the spatial semidiscretization of the Laplacian that yields to the system (15). The IMSP schemes are based on the expression of F h in (15) as the sum of two operators for all x ∈ Ω h . Then, for i = 1, 2, 3, define y i such that y = y 1 + y 2 + y 3 and satisfy A partitioned Runge-Kutta scheme is employed to solve the previous system: the dynamics of y 1 is approximated by a diagonally implicit method, then the evolution of y 2 and y 3 variables is approximated by a partitioned symplectic Runge-Kutta method defined as In [33] the authors set s = 1 in (17) and apply the following two-stage schemes in order to approximate y 1 , y 2 and y 3 , respectively. The resulting method for the approximation of y = [u h , v h ] T is given by where y n ∈ V h approximates y(t n ), for t n = t 0 + n ∆t and Y n,1 h , Y n,2 h , Y n,3 h ∈ V h are intermediate steps. For β = 1 and α = 0 this scheme is featured by first order approximation; when β = α = 1/2, then it is a second order accurate method. Moreover, under the assumption that diffusivity coefficients D u and D v nullify, the first order scheme represents a symplectic Euler scheme which handles the variable u h in an explicit way and the variable v h implicitly. The second order approximation reduces to the classical Störmer-Verlet scheme when it is written as a partitioned Runge-Kutta method exploiting the trapezoidal rule in order to approximate the stiff diffusive term.
In terms of u h and v h , the first order IMSP scheme, when applied to the Rosenzweigh-MacArthur model (10) reads

A Linear Stability Analysis
A stability analysis of IMSP schemes in terms of the diffusion and the reaction time-scales was recently developed in [34]. Their numerical simulations reveal that IMSP schemes provide the best choice for spatio-temporal dynamics of standing oscillations around an equilibrium of centre type. Their dissipation analysis is based on the application of IMSP schemes to the ODE test system: The first order IMSP scheme (19), introduced in [33], which approximate the diffusive term related to the coefficient λ with implicit Euler and the reaction term with symplectic Euler, can be written as By eliminating u n 1 and v n 1 and setting ξ = λ ∆t, ν = β ∆t the numerical solution can be written as By denoting with λ ± M = ρ M e ±i θ M the eigenvalues of M, and comparing the numerical with the exact solution written as = e ξ cos(ν n) −sin(ν n) sin(ν n) cos(ν n) where the quantities δ(ξ, ν) = e ξ − ρ M (ξ, ν) and φ(ξ, ν) = ν − θ M (ξ, ν) represent the dissipative and dispersion errors, respectively. Since ξ < 0, we require that the IMSP method is dissipative i.e., ρ M (ξ, ν) < 1. Since the eigenvalues of M are To study the dispersion properties of the IMSP method firstly we notice that λ ± . The dispersion error is calculated by It results that φ(ξ, ν) = −ν 3 /24 + O(ν 5 ) for ξ → 0, with a constant error smaller than that of the explicit ADI scheme and of IMplicit-EXplicit (IMEX)-Euler scheme [34].

Analysis of Semi-Discrete in Time IMSP First Order Scheme in Weak Form
The results of a more technical methodology based on the analysis of a semi-discrete in time formulation of IMSP schemes has been performed in [32]. Let D = diag( D u , D v ) be a linear matrix operator and consider the vector G(y) = [ f (u, v), g(u, v)] T . With the previous notations, the system (11)-(12) leads to the following continuous-in-time weak formulation: find y(·, t) = [u(·, t), v(·, t)] ∈ H 1 (Ω) × H 1 (Ω) such that y t , χ + (D∇y, ∇χ) = (G(y), χ) (20) for all χ ∈ H 1 (Ω) × H 1 (Ω) and for almost every t ∈ (0, T). Define the vectors G (u) = [ f (u, v), 0] T and G (v) = [0, g(u, v)] T ; then IMSP scheme (in weak form) is defined as follows: for n = 0, . . . , N − 1, find y n 1 , y n 2 , y n 3 , y n+1 ∈ H 1 (Ω) × H 1 (Ω) such that ∀χ ∈ H 1 (Ω) × H 1 (Ω): For β = 1 and α = 0 we recover the first-order accurate IMSP scheme. In terms of the variables u and v, it can be written as follows. For n = 0, . . . , N − 1, find v n 1 , u n 1 , v n+1 , u n+1 ∈ H 1 (Ω) such that for all χ ∈ H 1 (Ω) where we adopted a uniform mesh grid of N + 1 points t n = n∆t, n = 0, . . . , N with constant time step ∆t = T/N to discretize the temporal horizon (0, T). Moreover, v 0 = v 0 (·) and u 0 = u 0 (·) are the initial densities of predators and prey respectively. We summarize the main theoretical results proven in [32] obtained under the assumptions that f (u, v) has logistic dominated growth in the first variable, namely and the function g satisfies a sub-linear growth in the second variable, i.e., there exists C g > 0 such that Theorem 3 (Positivity). Assume the time step ∆t < 1/L and Ω is a domain of class C 1 . Provided that the initial conditions are positive, i.e., u 0 ( for all x ∈ Ω, then the solutions u n (x), v n (x) of the first-order scheme (21) are positive for all n ≥ 0.
and assume the nonlinearities f , g are globally Lipschitz, i.e., there exists L > 0 such that for all u i and v i in a compact subset of R + × R + . Then, there exists a constant C(u, v) > 0 such that max 0≤n≤N ( e n u + e n v ) + ∆t The above error estimate implies convergence and first-order accuracy in time for the solution u n , v n of Equation (21) with the assumption that e 0 u , e 0 v are also of order ∆t.

Numerical Comparison between IMEX and Geometric Integrators
Consider the dynamics of Equation (10) The five-point central difference approximation of the Laplacian in two dimensions on a grid Ω h with stepsize h = 1 was used. The spatially explicit dynamics of the phytoplankton-zooplankton model is shown in Figure 4 where the prey (u) approximations obtained with first order EXPP method given in Equation (16) are plotted. We see that spirals appear with their centers located in the vicinity of equilibria u * = 6/35 and v * = 116/245 of the dynamics. At t = 300 (see Figure 4) we observe the onset of destruction of the spirals which leads to the formation of two growing embryos of the patchy spatial pattern and finally, at t = 1200, the irregularities occupy the whole domain. In [29] the onset of the spatial patterns near the critical points is already evident at t = 120. To explain such a difference in the transient to the formation of irregular patchy spatial patterns the dynamics evolution in the temporal interval [0 120] in a smaller domain Ω = [0, 300] × [0, 300] was considered. The IMEX first order scheme was compared with EXPP in (16) and IMSP given in (19). The results are shown in Figures 5-7 and summarized below.
• For IMEX and IMSP schemes (Figures 5 and 6) we start with a temporal stepsize ∆t = 1/3 (in order to have convergence) and then we reduce the stepsize. As ∆t is reduced, the initial spiral patterns generated by the IMEX scheme with large temporal stepsize disappear, thus confirming that the spiral pattern in [29] is a numerical artifact. • This is not the case of EXPP method (Figure 7), which exhibits no spiral, even with the larger initial stepsize ∆t = 1. By reducing the stepsize the solution does not change significantly, with this meaning that convergence is reached and no artifacts arised.   As a conclusion, we have observed that a numerical method applied to the Rosenzweig-MacArthur model (10), may alter the solution in the short run. To make correct qualitative considerations on the transient dynamics, classical methods as IMEX scheme need a very small stepsize. Geometric IMSP and EXPP schemes both show a more stable behaviour; however, the positive EXPP scheme provides the best performance.

Symplectic-Exponential Lawson Integration for the Control of Invasive Population
Mathematical modeling and optimization play an important role in improving strategies for the control and eradication of invasive species [15,35,36], a crucial aspect in nature conservation [37].
In [16], a reaction-diffusion partial differential equation models the spatiotemporal dynamics of the invasion of alien fish population inside a very complex geometry representing a realistic lake and the optimal control theory has been applied to obtain the optimal management strategy with a given budget constraint. The model describes a complex and realistic situation by including a control term that has Holling-II type behavior, a budget constraint and the habitat suitability function which represents the suitability of habitat for a given species based on known affinities with environmental parameters.
The goal is to minimize the environmental damage over time at the minimum cost, in terms of the resources allocated to the species harvesting. The objective function is a sum of three terms. The first term takes into account of the environmental damage through the function ω(x, t) > 0 over Ω × [0, T], ω ∈ L ∞ (Ω × [0, T]), which is assumed to have a cost which increases with the presence of the invasive species: where u(x, t) represents the population density at time t ∈ [0, T] and (vector) position x in an open bounded ν ∈ L ∞ (Ω), ν(x) ≥ 0 is a weight for the final population density and δ ∈ (0, 1) is the discount factor. The second term is the cost due to the control effort: The third term accounts for the budget constraint 0 < E < B introduced as a penalty term with weight c ≥ 0; The optimal control model searches for a control such that J(E * ) = min E∈U J (E), subjects to the state equation where D > 0 represents the constant diffusion coefficient and n is the outward normal vector on ∂Ω.
The optimality conditions yield the following boundary value system As before, the approximation process of the whole dynamics assumes a semi-discretization in the space variable; the resulting procedure leads to a boundary value problem governed by the ordinary differential system (15) in the time variable, which can be numerically integrated by a partitioned method which consists of a Runge-Kutta scheme and the exponential Lawson algorithm [38] related to the symplectic counterpart [39].

Exponential Lawson Symplectic Integration of the Reaction Semiflow
We illustrate the advantage of using a symplectic numerical scheme when we approximate the boundary value problem described by the reaction semiflow. For simplicity, we assume ρ = 1 and ν > 0 and ω > 0 as constant functions. Consider the boundary value ODE system: with u(0) = u 0 and λ(T) = ν, equipped with and ϕ(u, λ) is defined in (27). We use the forward-backward sweep approach described in [40] joint with an exponential Lawson symplectic scheme [39]. This choice is motivated by the fact that, for all x ∈ Ω h the ODE system (28) and (29) is what is referred in literature as a nearly Hamiltonian system : where It has been proven that a partitioned method which consists of a Runge-Kutta scheme for the state variable and the exponential Lawson [38] symplectic counterpart for the current costate, is an effective alternative to the standard ODE solvers in the framework of optimal control [39].
In more details, given a Runge-Kutta scheme defined by coefficients a i,j , b i , with b i = 0 and its symplectic counterpart defined by coefficientsâ i,j , b i , with b k b j = b kâkj + b j a jk for each k, j and c k =ĉ k = ∑ s j=1âk,j with c s = 1, the scheme applied to (28) and (29) proceeds according to the following steps: • make an initial guess for λ n k on the mesh time t n k = t n + c k ∆t for n = 0, . . . N − 1 and k = 1, . . . s, with ∆t = T/N; • using initial conditions u 0 , final condition λ(T) = ν and the guess values λ n k , solve u forward in time according to the Runge-Kutta scheme: a kj F(u n j , λ n j ), k = 1, . . . , s, • using the transversality condition λ N = ν, solve λ(t) backward in time by means of the exponential Lawson symplectic counterpart b k e −δ∆tc k G(u n k , λ n k ).
Notice that, for k = s, since c s = 1, it results that λ n s = λ n+1 . Moreover, (b j −â kj ) = b j b k a jk . This implies that, if the Runge-Kutta scheme used on the state variable u is explicit (i.e., a jk = 0 for j ≤ k), the corresponding Lawson backward scheme for the costate variable λ(t) is explicit too: for n = N − 1, N − 2, . . . , 0 and k = s − 1, . . . , 1.

Implicit-Exponential Lawson Euler Integration
By using the previous notations, the spatial discretization of the Laplacian operator transforms the problem (26) into the BVP problem where ∆ h represents an approximation of ∆ on Ω h . The solution y := [u h , λ h ] T represents an approximate function of [u, λ] T such that, for all t ∈ [0, T], u h (t) and λ h (t) belongs to V h , the vector space of all the functions from Ω h to R. For t ∈ [0, T] and F h is the operator that associates to (31). Consider the splitting with f [D] (y) denoting the diffusion semiflow and f [R] (y) the reaction semiflow. We want to solve the diffusion semiflow with implicit Euler and the reaction semiflow with the exponential Lawson symplectic scheme [39] and apply the schemes within a forward-backward procedure [40] for solving the boundary value problem (32). In details, set λ N h = ν and make an initial guess for λ n h for n = 0, . . . , N − 1 on the mesh time t n = t n−1 + ∆t. Evaluate u n h with the forward formulā Then, the backward formula proceeds evaluating λ n−1 h starting from λ N h = ν, Finally, check for the convergence and evaluate E n h (x) = min max ϕ(u n h (x), λ n h (x)), 0 , B , for x ∈ Ω h and n = 0, . . . N.

Simulation of an Optimal Abatement Program
In [16] a control action over a time horizon with length T = 4 for the invasion of an alien fish population inside the hypothetical lake is simulated. A penalty term weighted by coefficient c = 0.22 for the budget constraint 0 ≤ E ≤ B = 0.5 was set. A spatial mesh Ω h consisting of 692 triangles and 427 nodes was used and a time-step ∆t = 1/36 was chosen for the temporal procedure. The parameters for fish population are set to r = k = 1, µ = 20, h = 1, with a diffusivity coefficient D = 0.05; the parameters related to control dynamics are set as ν = 5 e −0.4 , δ = 0.1, ω = 1. The dynamical evolution at t = 1, 2, 3, 4 , starting at t 0 with a uniform population distribution corresponding to the maximum density n 0 (x, t) = 1, is shown in Figure 8.
The optimal control tends to nullify the fish population almost everywhere except for some small areas where it gets its maximum allowed value. Fish population tends to a final distribution with mean value 0.0012, integrated over the whole spatial domain. The maximum allowed value E = B = 0.5 is achieved in the areas where the sensitivity function gets its largest values.

Future Challenge: Geometric Numerical Integrators for Multi-Structured Systems
A simplified spatiotemporal model developed in [41] to study the implication of the destabilization of phytoplankton dynamics due to rising temperatures, is motivating the future research on geometric numerical schemes for dynamical flows featured by multiple structures.
Here we will focus on systems which fall in the class of both Poisson and biochemical systems. We recall the definition of biochemical system and the main qualitative properties associated to their dynamics. A biochemical system consists of a set of M reactions of the type where X i are the chemical species in the system, σ − ij and σ + ij are the stoichiometric coefficients of the reactants and products, respectively. The dynamics is governed by the ODE systeṁ y = f(y) = S r(y), y(0) = y (0) (34) where y = (y 1 , . . . , y N ) T is the vector of the concentration of chemical species, y 0 > 0 is the fixed vector of initial concentrations, S − = σ − i,j and S + = σ + i,j are N × M negative and positive parts of the N × M stoichiometric matrix S = σ i,j i.e., S = S + − S − . The reaction function r(y) = (r 1 (y), . . . , r M (y)) T is defined elementwise as Well-posedness and positivity of the solution are guaranteed under the assumptions stated in [17]. All linear functions

Biochemical-Poisson System for Modelling the Oceanic Deep Chlorophyll Maximum (DCM)
To study the implication for oceanic primary production, species composition and carbon export of the destabilization of phytoplankton dynamics due to rising temperatures, a simplified version of the vertical water column model given in [41] is here described.
Let x indicate the depth in the water column. The population dynamics of the phytoplankton population P and nutrient N is described by a reaction-advection-diffusion system: where the specific growth rate of the phytoplankton is proportional through a costant g p to an increasing saturating function F(N, I) of nutrient availability N and light intensity I, m P is the specific loss rate of the phytoplankton, q is the nutrient content of the phytoplankton, W g is the phytoplankton sinking velocity and k P and k N are the vertical turbulent diffusivity of phytoplankton and nutrient, respectively. In this simplified model we assume that all the nutrient in dead phytoplankton is recycled.
In an operator splitting perspective, to construct approximated solutions of (35), which separate diffusive and advection semiflows from reaction, a step requires the numerical integration of the reaction semiflow. For for constant value of light intensity I, the reaction semiflow is given by the autonomous and spatially homogeneous system ∂N ∂t = − q g p F(N, I) P + q m p P ∂P ∂t = g p F(N, I) P − m p P.
It results that (36) is an example of biochemical-Poisson system. Indeed, by neglecting the dependence on the spatial coordinate and setting y = [N, P], system (36) possess two main structures:

•
Biochemical structure:ẏ = Sr(y) with r(y) = g p F(N, I) P, m p P T , and is the stoichiometric matrix with rank(S) = 1.

•
Poisson structure:ẏ = B(y)∇H with H (y) = N + q P and is a full rank matrix which satisfies the Jacobi identity (1).

Geometric Integrators for Biochemical-Poisson Systems
With the aim of constructing geometrical numerical integrators for a biochemical-Poisson systems, we need to enumerate all the desired characteristics. In general, an ideal geometric numerical integrator y n+1 = Φ h (y n ) for an (autonomous) biochemical-Poisson systemẏ = f(y), with y ∈ R n should: 1.
preserve positivity as negative concentrations are non-physical; 2.
preserve n − s linear invariants (as biochemical system) where s is the rank of the stoichiometrix matrix S ∈ R n×m , 3.
preserve the Hamiltonian and n − q Casimirs as Poisson system where q is the rank of B ∈ R n×n ; 4.
preserve the Poisson structure of the flow i.e., Φ ∆t (y n )B(y n )Φ ∆t (y n ) T = B(Φ ∆t (y n )); where Φ ∆t (y n ) denotes the Jacobian of the discrete flow.
For a generic biochemical-Poisson systems, a numerical integrator cannot be featured by all the above properties (it is enough to consider that, in general, for a Poisson system a numerical integrator cannot be both a Poisson and an Hamiltonian-preserving map); hence, an ideal geometric numerical integrator cannot be constructed.
For example, the application of Poisson Euler in the form (7) to the biochemical-Poisson system (36) leads to the numerical method P n+1 = P n e ∆t (g p F(N n , I) − m p ) , N n+1 = N n e ∆t −q g p F(N n , I) P n+1 + q m P P n+1 N n which is a positive explicit integrator for the system (36). However, it does not preserve linear invariants nor the Hamiltonian and the question if it is a Poisson integrator in correspondence of the above matrix B(y) is still an open question as it depends on the functional form of the function F. The most recent numerical integrators developed with the aim of preserving the qualitative features of a biochemical dynamics are gBBKS schemes [43]. They provide positive approximations and preserve all linear invariants of the continuous flow. The first-order variant applied toẏ = f(y) is given by y (n+1) = y (n) + ∆t   ∏ j∈J n y (n+1) j y (n) j   1 q f(y (n) ), J n = i : f i (y (n) ) < 0, i ∈ {1, . . . , N} , with q ≥ card(J n ). The application of the above method for solving (36) guarantees preservation of positivity, preservation of Hamiltonian which is a linear invariant, and again is still an open question if there exists a class of function F for which gBBKS scheme may result also a Poisson map. In the following, we suggest two test models, one linear and one nonlinear, that can be useful in order to compare the performance of existing numerical methods and can provide the basis for the construction of novel geometrical integrators for biochemical-Poisson dynamics [44].

Biochemical-Poisson Test Models
In order to compare the schemes when applied to biochemical-Poisson systems, a possible way is to evaluate their performances on solving the two test problems we are going to describe in the following section. We will see that in the general framework of non-standard finite difference approximation [45], their exact solutions can be reproduced by suitable non-standard numerical procedures. These have been the basis for the construction of novel geometrical schemes for biochemical systems which, despite their explicit fuctional form, are able to preserve both positivity and linear invariants [44].
As linear test problem we consider the two-dimensional systeṁ with a, b, q > 0. Having set y = [u, v] T , is is easy to see that the above test system is a biochemical-Poisson system since it can be written both asẏ = S r(y) with S = −q q 1 −1 r(y) = a u b v and asẏ = B(y)∇H with H (y) = u + q v and a skewsymmetric matrix which satisfies the Jacobi identity (1). The stoichiometric matrix S has rank one, hence there exists a unique linear invariant given by H (y). On the other side B has full rank, this implying that there are no Casimirs. The theoretical solution can be written as follows (see [45]): where u(0) = u 0 and v(0) = v 0 denote the initial values. Then, the exact scheme corresponds to the non-standard procedure u n+1 = u n + φ(∆t) (−q a u n + q b v n ), v n+1 = v n + φ(∆t) (a u n − b v n ), where φ(∆t) = 1 − e −∆t (q a+b) q a + b .
As nonlinear test system we consideru = − q u v v = u v with q > 0. Having set y = [u, v] T , is is easy to see that the above test system is a biochemical-Poisson system since it can be written both asẏ = Sr(y) with S = −q 1 with r(y) = u v and result an M-systemẏ = B(y)∇M with M(y) = u + q v and B given in (2). The stoichiometric matrix S has rank one, hence there exists a unique linear invariant given by M(y). On the other side B has full rank, this implying that there are no Casimirs. The unique invariant is given by the Hamiltonian M(y) which results to be linear. The theoretical solution can be explicitly written as Then, the exact scheme corresponds to the non-standard procedure u n+1 = u n − φ(∆t) q u n v n , v n+1 = v n + φ(∆t) u n v n , where φ(∆t) = e (u n + q v n ) ∆t − 1 u n + q v n e (u n + q v n ) ∆t .

Conclusions
In ecological modelling, the approximated solutions should be able to reproduce the main physical qualitative characteristics of the observed quantities (e.g., positive concentrations, mass/energy conservation) in order to make accurate previsions or outline realistic scenarios. For this reason, geometric numerical integration should play a crucial role in the analysis of ecosystem models described by systems of differential equations. Nevertheless, numerical schemes for ecological models have received little attention in the literature as most descriptions of models outline the governing equations but do not discuss their numerical solution. Unfortunately, an inadequate choice of a numerical method may have a detrimental effect especially within a conceptual modelling which is thought to make scenario analysis.
In this paper, we detect the more frequent flow structures that characterize some ecological models encountered in our research activity. Flow structures that more frequently arise in ecological modelling and that have to be preserved are Poisson maps, dynamics evolving in positive phase-space, preservation of invariants. We have described the application and the results of the analysis of geometric numerical integration in the framework of ecological modeling, focusing on first-order procedures. In general, what we experienced is that the mechanism that produces positive solutions, as in case of Poisson Euler scheme, is also able to produce more stable solutions. This produces some benefits not only in the long run but even in the transient phase. For optimal control problems that model management policies, we have described how to use a Lawson symplectic Runge-Kutta pairs in order to get an efficient method for solving the boundary value problem for the nearly Hamiltonian system that provides optimality conditions.
As future perspective, for dynamical flows which are featured by multiple structures, the identification of the main qualitative properties that are crucial to preserve in order to have the best performances, will deserve a deep investigation.