\documentclass[10pt]{article}
\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}
\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
linkcolor = blue,
urlcolor = blue,
citecolor = blue,
anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
\errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
\usepackage{natbib}
\renewenvironment{abstract}
{{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
{\bigskip}
\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}
\usepackage{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}
\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}
\usepackage{float}
\begin{document}
\title{Modelling predation: Theoretical criteria and empirical evaluation of
functional form equations for predator-prey systems}
\author[1]{Julien Malard}%
\author[1]{Jan Adamowski}%
\author[1]{Jessica Bou Nassar}%
\author[2]{Nallusamy Anandaraja}%
\author[3]{Héctor Tuy}%
\author[1]{Hugo Melgar-Quiñonez}%
\affil[1]{McGill University - MacDonald Campus}%
\affil[2]{Tamil Nadu Agricultural University}%
\affil[3]{Rafael Landivar University}%
\vspace{-1em}
\date{\today}
\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup
\selectlanguage{english}
\begin{abstract}
Correct modelling of relationships between predators and prey is crucial
to ecological and population dynamics models. However, and despite a
long-standing competition between ratio and prey-dependent models (and a
few alternative intermediate forms) in the literature, most equations
currently used to represent such relationships do not meet theoretical
criteria for biological consistency. This research proposes a set of
universally applicable criteria for all predation equations and shows
that the most commonly used predation equations in the literature fail
to meet these same criteria. We follow with a proposal for a new
predation equation that does meet these criteria, which combines both
prey and ratio-dependent concepts while giving reasonable predictions in
the cases of both high predator or high prey densities. We show its
empirical performance by applying the new equation, along with existing
alternatives, to various experimental predation datasets from the
literature. Results show that the new equation is not only more
mathematically consistent than existing equations, but also performs
more consistently empirically across different datasets from various
ecological situations. This research is the first to propose a
systematic set of criteria to evaluate predation equations and then to
offer an equation that meets these criteria and also performs well both
theoretically and empirically across datasets from a wide range of
predation systems.%
\end{abstract}%
\sloppy
\section*{Introduction}
{\label{introduction}}
Correctly representing the relationship between predators and prey is a
key challenge to all attempts at modelling trophic interactions and
population dynamics (Kratina \emph{et al.} 2009). A wide range of
predation equations is currently available in the ecological literature
to calculate predation rates (prey consumed per predator per unit time)
based on an independent variable (prey density, predator density, or
both) and one or more empirical parameters (Abrams \& Ginzburg 2000).
These equations can be, very generally, classified according to the
independent variable and the mathematical form (functional response) of
the equation that links this variable to predation rates.
According to the independent variable, predation equations are often
divided into prey-dependent (where the prey population (N) is the only
independent variable) and ratio-dependent equations (where the ratio
between prey and predator (P) populations, \(NP\), is used
as an independent variable). A wide variety of intermediate forms have
also been proposed. Of these, we here mention some of the most commonly
used, the Hassell-Varley forms (Hassell \& Varley 1969) and the
Beddington-DeAngelis equation (Beddington 1975; DeAngelis \emph{et
al.}1975). The former represent a compromise between prey and ratio
dependence, where the \(\frac{N}{P}\) prey-predator ratio is replaced
by\(\frac{N}{P^{m}}\), where \emph{m} is a parameter in the
range\(\left[0,1\right]\) that allows the function to transition smoothly
between either prey or predator dependence. The Beddington-DeAngelis
equation, on the other hand, incorporates both N and P, but not in pure
ratio form.
Functional responses, on the other hand, are generally categorised as
Type I (linearly increasing, \(y=\alpha x\)), Type II (logistic-type
growth; \(y=\selectlanguage{greek}\frac{\text{αx}}\selectlanguage{english}{x+b}\)), Type III (sigmoidal type growth;
\(y=\frac{\alpha x^{2}}{x^{2}+b}\)), and Type IV (increasing to a peak, followed by a
decreasing response) functional responses (Holling 1959), where x is the
independent variable (N in the case of prey dependence
and\(\ NP\) in the case of ratio dependence). Type I rarely
occurs (with the notable exception of filter feeders) (Jeschke\emph{et
al.} 2004), Type II indicates predator saturation, and type III combines
predator saturation at high prey densities with sigmoidal growth in
predation at lower prey densities (which could be caused by the
saturation of limited prey refuges in the habitat, or be indicative of
preferential prey switching or learning by the predator). Type IV would
occur in the case of dangerous prey (L\selectlanguage{ngerman}íznarová \& Pekár 2013), where
higher prey densities actually diminish the predator's capacity to feed.
Some equations (Hassell-Varley, prey-dependent, and ratio-dependent
forms) can be combined with all response types, while others (such as
Beddington-DeAngelis) do not lend themselves to such modifications but
intrinsically generate functional forms similar to one or more of the
Holling forms.
None of these equations, however, is completely satisfactory from a
theoretical standpoint. Ratio dependence is unrealistic at low absolute
prey densities, when the physical scarcity of prey is likely to decrease
predation rates as predators spend most of their time searching, even if
the predator population is itself very low and the prey-predator ratio
correspondingly high (Abrams \& Ginzburg 2000), and also produces what
have been described as ``unusual'' results (namely, infinite prey
availability per predator) as predator density declines towards 0
(Gleeson 1994). On the other hand, prey dependence is unrealistic at low
prey-predator ratios, when predators must necessarily divide available
prey amongst themselves (Abrams \& Ginzburg 2000), and studies have
found experimental data partially between ratio and prey dependence
(Schenk \emph{et al.} 2005).
There is therefore a need for an equation that encompasses both
processes; while alternatives include equations of intermediate form
partially between ratio and prey dependence, their formulations, as will
be discussed later, do not necessarily guarantee mathematically sound
resolution of the different issues faced by each extreme. For instance,
the transition between prey and predator dependence in the
Hassell-Varley forms is dependent on a parameter, not on prey or
predator populations, while the theoretical considerations mentioned
above suggest that population levels are important in determining when
either relationship is more appropriate. Switching equations based on
the nature of predator interference in the system has also been proposed
(Skalski \& Gilliam 2001).
There is therefore a need for clear criteria and expected behaviours for
all useful predation equations, which may promote a more standardised
development of predation equations that are widely useful across various
population densities and ratios. Such a list of criteria has been
previously proposed by (Berryman \emph{et al.} 1995), but most relate to
complete predator-prey models and not to the predation function itself.
Of the criteria that do apply are that predators ``must have finite
appetites'' and that the potential for intraspecific competition should
be included, criteria which disqualify prey-dependent equations.
However, these criteria apply only to per predator feeding rates and do
not include important conditions regarding total predation rates under
extreme (high or low) prey or predator populations.
Despite the long-standing debate, many of these issues have still not
been satisfactorily resolved (Abrams 2015), and these equations are
widely used in theoretical and practical research, despite their
previously identified shortcomings. Even prey dependence and Type I
functional responses, which have been shown to be in most ecological
systems generally inferior to other (though imperfect) alternatives
(Abrams \& Ginzburg 2000; Jeschke \emph{et al.} 2004) have been recently
used as part of larger ecological and population dynamics models
(Chesson \& Kuang 2008; Sanchez \emph{et al.} 2018; Sanders \emph{et
al.} 2018; Imbert \emph{et al.} 2020).
These shortcomings occur because most of these equations are built upon
assumptions (mainly mostly stable prey and predator population
densities) that are unlikely to hold true under natural conditions.
These are structural equation (model) errors, which cannot be
circumvented through parametrisation and calibration. As such, the
predictions of these equations under extreme conditions (such as those
often seen in agriculture or other human-managed systems) may be
invalid. For instance, equations whose assumptions do not hold under
conditions of very high predator density are unlikely to give correct
predictions with regards to augmentative biocontrol efforts, where large
numbers of predators or parasitoids are released into the field. This is
a major impediment to the successful application of predator-prey
equations to ecological, and, especially, agroecological modelling. As
(Kratina \emph{et al.} 2009) have observed, ``choosing appropriate
functional responses is crucial for adequate predictions of food web
dynamics. When predator dependence is incorporated into predator--prey
models their stability is usually enhanced,'' highlighting the
importance of correct functional response equations to ecological model
behaviour and reliability.
This paper addresses this gap in the literature by 1) proposing a
standard set of criteria that all useful predation equations should
observe, 2) showing that the majority of equations available in the
literature fail to meet one or several of these conditions, and 3)
proposing a new equation (``Kovai equation'') for modelling predation
rates that does meet these criteria. The paper concludes with 4) a
comparative evaluation of the performance of all reviewed equations by
applying them to a range of predation datasets from the literature to
demonstrate the practical improvements offered by the improved predation
equation formulation. This is the first research article to
comprehensively address limitations in existing predation equations and
to offer a mathematically sound improvement to these equations that
performs well both theoretically and empirically.
\section*{Mathematical background}
{\label{mathematical-background}}
Predation equations may be defined as the predation rate per predator
per unit time, \(f\left(N,P\right)\), or as the total predation rate per
unit time, \(g\left(N,P\right)=P*f\left(N,P\right)\) where\emph{N} and \emph{P} are the prey
and predator densities, respectively. The following sections present a
proposed set of standard criteria for predation equations, followed by a
proposed equation that satisfies them. In the following sections,
\(\alpha\) denotes the maximum feeding rate for the predator, or
the amount of prey per unit time that an individual predator can consume
if no time is wasted searching for prey.
As comparison, commonly used equations from the literature were
evaluated according to the same criteria, including prey-dependent,
ratio-dependent, and the intermediate Hassel-Varley equation (Hassell \&
Varley 1969), each combined with Type I, II and III Holling functional
response forms (Holling 1959; Abrams \& Ginzburg 2000). A final
equation, the Beddington-DeAngelis equation (which does not lend itself
to combination with Holling functional response forms) was also tested
(Beddington 1975; DeAngelis \emph{et al.} 1975). See Table 2 for a full
definition of each equation.
\subsection*{Behaviour at high prey
density}
{\label{behaviour-at-high-prey-density}}
As prey density increases, predation per predator should asymptotically
approach the maximum feeding rate for the predator. In mathematical
terms, \(\operatorname{}{f\left(N,P\right)}=\alpha\), and, by extension, \(\operatorname{}{g\left(N,P\right)}=\alpha P\).
\subsection*{Behaviour at low prey
density}
{\label{behaviour-at-low-prey-density}}
As prey density approaches 0, both total as well as per predator
predation should tend towards 0, that is,\(\operatorname{}{f\left(N,P\right)}=\operatorname{}{g\left(N,P\right)}=0\).
\subsection*{Behaviour at high predator
density}
{\label{behaviour-at-high-predator-density}}
As predator density increases, overall predation should increase
asymptotically towards the total prey population available (or perhaps
slightly less, if it is assumed that some habitat niches offer complete
predation protection to a limited number of prey), while prey
availability per predator decreases as predators begin competing for
prey. In mathematical terms,\(\operatorname{}{g\left(N,P\right)}\leq N\),
though\(\operatorname{}{f\left(N,P\right)}=0\).
\subsection*{Behaviour at low predator
density}
{\label{behaviour-at-low-predator-density}}
On the other hand, as predator density decreases, the total predation
should also approach 0: \(\operatorname{}{g\left(N,P\right)}=0\).
\subsection*{Other considerations}
{\label{other-considerations}}
Finally, many equations consider the functional response of predation
rates to a measure of prey abundance (often prey density or else the
prey-predator ratio), based on the assumption that predation rates are
not a constant function of prey abundance and will instead likely change
nonlinearly along with prey abundance (for instance, due to predator
saturation, or else to the existence of niches where a subpopulation of
prey may be partially shielded from predation). Two of the three types
(Type II and Type III) of functional response proposed by Holling
(Holling 1959), and widely used since, were designed to model such
effects.
However, although one may reasonably expect that additional prey added
to a system may be more vulnerable to predation than already present
prey (due, for instance, to the saturation of protected niches), the
addition of new prey should in no circumstance increase the
instantaneous predation risk for prey that was already present. In other
words, adding 1 unit of prey cannot increase total predation by more
than 1 unit of prey, or, mathematically,\(0\leq\frac{\text{dg}}{\text{dN}}\leq 1\).
\subsection*{The Kovai equation}
{\label{the-kovai-equation}}
As can be later seen in Table 2, most equations currently used for
modelling predation fail in at least one of the above-mentioned
criteria. We here describe a new equation, which we propose to call by
the name of the city where research leading to its development was
conducted.
\begin{equation}
f\left(N,P\right)=\alpha\selectlanguage{greek}\left(1-e^{-\frac{\text{Nu}}{\left(\text{αP}\right)}}\selectlanguage{english}\right)\nonumber \\
\end{equation}\begin{equation}
u=1-\frac{b}{N}\left(1-e^{-\frac{N}{b}}\right)\nonumber \\
\end{equation}
The first part of the Kovai equation takes the ratio between predator
and prey into account, partitioning available prey amongst the predators
present. In this sense, the Kovai equation is similar to a
ratio-dependent equation, but instead of using the raw prey density to
calculate the prey-predator ratio, a measure of effective prey
availability, \emph{u} , is used that takes the scarcity of the prey
(and prey refuges) into account and ranges between 0 and 1. This
performs a similar purpose as the Holling Type III functional response,
but retains the notable advantage that\(\frac{d}{\text{dN}}uN=1-e^{\frac{-N}{b}}\), which is in
the range\(\left[\left.\ 0,\ 1\right)\right.\ \) for non-negative\emph{b} , thereby ensuring
that the effective prey availability will never increase faster than the
prey density itself does. \emph{b} may be interpreted as the maximum
capacity of prey refuges (\(\operatorname{}\left(N-Nu\right)=\lim_{N\rightarrow\infty}\left[N-N+b\left(1-e^{-\frac{N}{b}}\right)\right]=b\)), or, alternatively, as
the prey density at which\(e^{-1}\approx 36.8\%\) of the prey is accessible
to the predator.
At low values of \emph{b} , the equation presents a functional response
form similar to the Holling Type II form, allowing this parameter to
control the functional response form of the equation. In addition, the
fact that \emph{u} is dependent on N and not on the ratio N/P enables
the Kovai equation to account for prey scarcity in a similar manner to
prey-dependent equations.
\section*{Methods (empirical study)}
{\label{methods-empirical-study}}
12 experimental studies with predation (functional response) data were
taken from the literature, for a total of 29 distinct datasets (Table
1). Where data was only available in the form of figures instead of
tables, these were first digitised using the web version of
WebPlotDigitizer (Rohatgi 2019). All analyses were run on Python 3.7
(Python Software Foundation (PSF) 2019), using the libraries Numpy (van
der Walt \emph{et al.} 2011), Pandas (McKinney 2010) and PyMC3
(Salvatier \emph{et al.} 2016) for mathematical analyses; graphs were
plotted using Seaborn (Michael Waskom \emph{et al.} 2018) and Matplotlib
(Hunter 2007). All data and code used to run the analyses and plot the
figures present in this article is included in the supplementary
material.
Table : Studies and datasets used in the analyses. PD indicates prey
dependence, RD ratio dependence, HV Hassell-Varley, BD
Beddington-DeAngelis, and O other equations not reviewed in this study.
Roman numerals indicate Holling functional response form. Only equations
finally recommended or identified as giving good fit by the authors of
the corresponding study are mentioned in the last column.\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
Reference & No. datasets & System studied & Equations
used\tabularnewline
\midrule
\endhead
(Kratina \emph{et al.} 2009) & 1 & Benthic flatworm-paramecium & BD*, HV
III, RD III\tabularnewline
(Fernandez-Maldonado \emph{et al.} 2017) & 1 & Cannibalism in
\emph{Nabis pseudoferus} & O I\tabularnewline
(L\selectlanguage{ngerman}íznarová \& Pekár 2013) & 5 & Spider-ant (various species) & O
IV\tabularnewline
(Schenk \emph{et al.} 2005) & 1 & Paper wasp-shield beetle & HV
III\tabularnewline
(Schenk \& Bacher 2002) & 2 & Paper wasp-shield beetle & O
III\tabularnewline
(Long \emph{et al.} 2012) & 4 & Cannibalism in red king crab & PD
II\tabularnewline
(Kfir 1983) & 1 & Egg parasitoid of \emph{Phthorimaea operculella} & O
II\tabularnewline
(Edwards 1961) & 3 & Parasitism (various species) & None\tabularnewline
(Eveleigh \& Chant 1982) & 4 & Predatory mites & HV\tabularnewline
(Mertz \& Davies 1968) & 1 & Cannibalism in flour beetles &
O\tabularnewline
(Taylor 1988) & 4 & Ectoparasitism (various prey) & O II\tabularnewline
(Uttley 1980)** & 2 & Damselfly/black swimmer-cladoceran
&\tabularnewline
\bottomrule
\end{longtable}
* Slightly modified to give a Type III functional response.
** Data from figure in (Skalski \& Gilliam 2001)
Each of the equations in Table 2 was calibrated using Bayesian inference
(Hamiltonian Monte Carlo No U-Turn Sampler algorithm (Hoffman \& Gelman
2014)) with each of the datasets, using weakly informative priors.
In the case of experiments where prey densities were held constant
(theoretically infinite predation opportunities), the data was modelled
with a Poisson distribution as follows:
\begin{equation}
y\sim Poisson(\mu*n)\nonumber \\
\end{equation}
Where \emph{y} is the total number of predation events observed (across
all replicates), \emph{n} is the number of replicates, and \selectlanguage{greek}\emph{µ} \selectlanguage{english}is
calculated according to each equation's specification.
When (as was most often the case) prey were not replaced to maintain
constant density throughout the experiment, the data was modelled with a
Binomial distribution:
\begin{equation}
y\sim Binomial(n=\mu*n,\ p=\max\left(\frac{\mu}{N},\ 0.999\right))\nonumber \\
\end{equation}
With the same parameter definitions as above, and \emph{N} being the
initial number of prey. Note that \emph{p} had to be bounded to avoid
overflows with several equation types.
Following calibration, each equation's performance within a dataset was
ranked according to the Widely-applicable Information Criterion (WAIC)
(Watanabe 2013; Vehtari \emph{et al.} 2017), which takes the number of
parameters into account (and penalises equations with more parameters in
order to avoid overfitting); lower scores indicate better model fit.
\section*{Results}
{\label{results}}
\subsection*{Theoretical study}
{\label{theoretical-study}}
As indicated in Table 2, none of the equations currently available in
the literature fulfills these criteria, which casts doubt on their
ability to generate accurate predictions for systems where predator or
prey concentrations deviate significantly from static equilibrium
conditions. The new proposed equation, however, did yield satisfactory
results with regards to all conditions.
{[}Insert Table 2 here{]}
Only expectations with regards to equation behaviour under conditions
where prey populations approach 0 were uniformly respected, meaning that
all equations evaluated did accurately predict that per predator and
total predation rates fall to 0 when prey populations disappear.
However, each of the other conditions was violated by one or more
equation.
As could be expected, prey-dependent equations performed poorly in
conditions of high predator densities, while ratio-dependent and
Hassell-Varley equations fared better in these conditions. However, Type
I predation responses caused problematic behaviour under conditions of
high prey densities in all equation types.
Particularly problematic, however, were cases where total predation
rates under very high predator densities either exceeded the total prey
population outright (e.g., \(\operatorname{}\left(g\right)=\infty\)) or else could exceed it
under certain parameterisations of the equations (e.g.,
\(\operatorname{}\left(g\right)=\selectlanguage{greek}\frac{\text{αN}}\selectlanguage{english}{b}\)). In fact, only the new proposed Kovai equation
showed consistent behaviour with regards to total predation rates under
high predator densities. Additionally, none of the equations from the
literature correctly predicted the rate of total predation increase upon
the addition of an additional prey. In many circumstances, increasing
the prey density by 1 could lead to a much larger increase in total
predation rates, suggesting that the addition of one prey individual to
a system could greatly \emph{increase} the predation of the other prey
individuals already in the system. Such behaviour seems unrealistic in
the real world and could, in modelling studies, lead to the prediction
of higher predation rates than would be expected. Only the Kovai
equation showed consistent behaviour in this regard.
\subsection*{Empirical study}
{\label{empirical-study}}
As can be seen in Table 1, the datasets included in the empirical study
span a wide range of predator-prey relationships as well as equation
types and functional responses (though type II and III responses were
more common). Figure 1 shows the results of different models applied to
predation data from (Kratina \emph{et al.} 2009), including three
equations recommended by those authors as well as the Kovai equation.
At low values of P, where the prey-predator ratio is high enough to
avoid predator competition, the ratio-dependent model was unable to
model the impact of absolute prey scarcity on predation rates and
therefore confounded the effects of various prey densities (Figure 1,
a). On the other hand, intermediate form equations (Hassel-Varley Type
III, Beddington-DeAngelis, and Kovai; see Figure 1 b, c and d,
respectively), were successful in distinguishing the effects of prey and
predator densities. (Prey-dependent models were overall unsuccessful in
modelling the data, as these models cannot distinguish between different
values of P at all, and are not shown in the figure.)\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}
\end{center}
\end{figure}
Figure : Functional response (prey consumed per predator) observations
for various predator and prey densities of the Stenostomum
virginanum-Paramecium aurelia predator-prey system, along with fitted
model predictions (line represents mean; shaded area the standard
deviation of model posterior distribution). a) Ratio dependent Type III
equation, b) Hassell-Varley Type III equation, c) Beddington-DeAngelis
equation, d) Kovai equation. N indicates the prey density. Data from
(Kratina et al. 2009).
Figure 2 shows the applicability of various equation types across all
datasets used in this study. Given the Kovai equation's ability to
adjust to both Type II and Type III functional responses, we expected it
to show more consistent applicability across datasets with varying
functional response forms than pure Type II or Type III equations
(regardless of prey or ratio dependence). Datasets were classified by
their level of (dis)similarity to a Type II functional response, using
the average rank of Type II equations' fit (as compared to the fit of
equations with different functional responses) when fitted to the data.
Figure 2 compares the average rank (lower indicates better fit relative
to other equations) of different equation groups with the average rank
of Type II response equations for each of the different datasets.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}
\end{center}
\end{figure}
Figure : Rankings of various equations, where lower ranks indicate
better fit relative to the other equations, versus the Type II equation
rank for the same experimental dataset. Each point represents ranking of
one equation or group of equations (Type III, BD or KV) when fitted to
one experimental dataset. Lines indicate linear model regression for
each equation or group, and shaded areas represent 95\% confidence
intervals. Type II includes Hassell-Varley, prey-dependent and
ratio-dependent equations with Holling Type II response; Type III
includes the same equations with a Holling Type III response; BD
indicates the Beddington-DeAngelis equation; KV indicates the Kovai
equation.
As expected, the ranking of Type III equations improved as datasets
moved away from a Type II response, and the BD equation, whose shape
resembles a Holling Type II form, predictably improves in rank the
closer the dataset is to a Type II functional response form.
The Kovai equation, however, ranked consistently across datasets,
regardless of whether the data showed Type II functional response or
not. This is most likely due to its parametrisation, which allows it to
model both Type II and Type III functional responses. It is important to
note, however, that as the Kovai equation supposes that prey
availability will increase along with prey densities, it is not suitable
for representing Type IV functional responses where increasing densities
of dangerous prey reduce predators' ability to hunt effectively.
In conclusion, the Kovai equation provides several main improvements
over existing predation equations in the literature. 1) Its main
contribution is its improved theoretical consistency, in particular in
its predictions of total predation at high prey or predator densities
and the rate of change of total predation as more prey is added to the
system. 2) The Kovai equation also combines the best of prey and ratio
dependent concepts in a mathematically consistent manner, by adjusting
for predator competition (correctly partitioning available prey between
predators when prey become limiting) while still distinguishing between
low and high absolute prey density in the case of otherwise identical
prey-predator ratios. 3) The equation also allows for a (mathematically
consistent) smooth transition between Type II and Type III functional
responses based on its parametrisation.
\section*{Conclusions}
{\label{conclusions}}
While predator-prey relationships are central to all ecological
modelling of population dynamics, a mathematically consistent equation
to represent these relationships has proved elusive. This research
showed that the vast majority of the most common equations currently in
use in both theoretical and empirical (field) ecological research fail
to meet several basic criteria for mathematical consistency. We propose
a list of such criteria to facilitate the evaluation of predation
equations developed in the future and propose a new equation that does
meet these criteria. The empirical evaluation of the new proposed
equation also suggests that it is more robust to varying functional
response forms found in different ecological contexts.
\section*{Author contributions}
{\label{author-contributions}}
JM conducted the research and wrote the paper. JBN contributed to the
mathematical analysis. JA, HQM, NA, and HTY supervised the research
project. The authors declare no competing interests.
\section*{Acknowledgements}
{\label{acknowledgements}}
This research was supported by an~NSERC Discovery Grant~(Grant
number:~RGPIN-2015-05554) held by Dr. Jan Adamowski, as well as an~FQRNT
Bourse de 3e~cycle, a~Bourse du CRDI aux chercheurs candidats au
doctorat (Award number 107759-99906075-017), a~Bourse d'\selectlanguage{ngerman}études
supérieures du Canada Vanier and a Supplément pour études à l`étranger
Michael-Smith~(Programme de bourses d`études supérieures du Canada)
scholarship held by Julien Malard.
Table : Comparison of critical extreme condition properties of common
predatation equations.\(f\left(N,P\right)\) represents the predation rate
per predator function and \(\backslash ng\left(N,P\right)=P*f\left(N,P\right)\)represents the total
predation rate, where \selectlanguage{greek}α \selectlanguage{english}is (theoretically) the maximum predation rate
per predator, N and P are the prey and predator densities, respectively,
and other variables are empirical constants.\(a,\ b,\ c\ \in\left[\left.\ 0,\ \infty\right)\right.\ ;m\in\left[0,\ 1\right]\). Shaded
squares indicate undesirable behaviours. See Appendix A: Calculation of
non-trivial limits and derivatives for details on the calculation of the
results presented here.\selectlanguage{english}
\begin{longtable}[]{@{}lllllllllll@{}}
\toprule
\begin{minipage}[b]{0.07\columnwidth}\raggedright\strut
\textbf{Equation form}\strut
\end{minipage} & \begin{minipage}[b]{0.07\columnwidth}\raggedright\strut
\textbf{Equation form}\strut
\end{minipage} & \begin{minipage}[b]{0.07\columnwidth}\raggedright\strut
\(\mathbf{f}\left(\mathbf{N,P}\right)\)\strut
\end{minipage} & \begin{minipage}[b]{0.07\columnwidth}\raggedright\strut
\(\operatorname{}\left(\mathbf{f}\right)\)\strut
\end{minipage} & \begin{minipage}[b]{0.07\columnwidth}\raggedright\strut
\(\operatorname{}\left(\mathbf{f}\right)\)\strut
\end{minipage} & \begin{minipage}[b]{0.07\columnwidth}\raggedright\strut
\(\operatorname{}\left(\mathbf{f}\right)\)\strut
\end{minipage} & \begin{minipage}[b]{0.07\columnwidth}\raggedright\strut
\(\operatorname{}\left(\mathbf{g}\right)\)\strut
\end{minipage} & \begin{minipage}[b]{0.07\columnwidth}\raggedright\strut
\(\operatorname{}\left(\mathbf{g}\right)\)\strut
\end{minipage} & \begin{minipage}[b]{0.07\columnwidth}\raggedright\strut
\(\operatorname{}\left(\mathbf{g}\right)\)\strut
\end{minipage} & \begin{minipage}[b]{0.07\columnwidth}\raggedright\strut
\(\operatorname{}\left(\mathbf{g}\right)\)\strut
\end{minipage} & \begin{minipage}[b]{0.07\columnwidth}\raggedright\strut
\(\frac{\mathbf{\text{dg}}}{\mathbf{\text{dN}}}\)\strut
\end{minipage}\tabularnewline
\midrule
\endhead
\begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Ideal equation\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Ideal equation\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\alpha\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left[0,N\right]\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\selectlanguage{greek}α\selectlanguage{english}P\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left[0,1\right]\)\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Prey-dependent\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Type I\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\selectlanguage{greek}α\selectlanguage{english}N\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\selectlanguage{greek}α\selectlanguage{english}N\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\infty\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\infty\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\infty\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left[\left.\ 0,\ \infty\right)\right.\ \)\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Type II\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\selectlanguage{greek}\frac{\text{αN}}\selectlanguage{english}{N+b}\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\selectlanguage{greek}\frac{\text{αN}}\selectlanguage{english}{N+b}\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\alpha\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\infty\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\selectlanguage{greek}α\selectlanguage{english}P\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left[\left.\ 0,\ \infty\right)\right.\ \)\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Type III\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\frac{\alpha N^{2}}{N^{2}+b}\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\frac{\alpha N^{2}}{N^{2}+b}\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\alpha\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\infty\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\selectlanguage{greek}α\selectlanguage{english}P\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left[\left.\ 0,\ \infty\right)\right.\ \)\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Ratio-dependent\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Type I\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\alpha\frac{N}{P}\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\infty\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\selectlanguage{greek}α\selectlanguage{english}N\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\selectlanguage{greek}α\selectlanguage{english}N\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\infty\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left[\left.\ 0,\ \infty\right)\right.\ \)\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Type II\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\frac{\alpha\frac{N}{P}}{\frac{N}{P}+b}\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\alpha\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\selectlanguage{greek}\frac{\text{αN}}\selectlanguage{english}{b}\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\selectlanguage{greek}α\selectlanguage{english}P\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left[\left.\ 0,\ \infty\right)\right.\ \)\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Type III\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\frac{\alpha\left(\frac{N}{P}\right)^{2}}{\left(\frac{N}{P}\right)^{2}+b}\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\alpha\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\selectlanguage{greek}α\selectlanguage{english}P\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left[\left.\ 0,\ \infty\right)\right.\ \)\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Beddington-DeAngelis\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Beddington-DeAngelis\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\selectlanguage{greek}\frac{\text{αN}}\selectlanguage{english}{N+cP+b}\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\alpha\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\selectlanguage{greek}\frac{\text{αN}}\selectlanguage{english}{c}\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\selectlanguage{greek}α\selectlanguage{english}P\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left[\left.\ 0,\ \infty\right)\right.\ \)\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Hassell-Varley\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Type I\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\alpha\left(\frac{N}{P^{m}}\right)\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left\{\begin{matrix}\alpha N,\ \ \&m=0\\
0,\ \ \&m>0\\
\end{matrix}\right.\ \)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\infty\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left\{\begin{matrix}0,\ \ \&0\leq m<1\\
\alpha N,\ \ \&m=1\\
\end{matrix}\right.\ \)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left\{\begin{matrix}\infty,\ \ \&0\leq m<1\\
\alpha N,\ \ \&m=1\\
\end{matrix}\right.\ \)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\infty\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left[\left.\ 0,\ \infty\right)\right.\ \)\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Type II\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\frac{\alpha\left(\frac{N}{P^{m}}\right)}{\left(\frac{N}{P^{m}}\right)+b}\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left\{\begin{matrix}\selectlanguage{greek}\frac{\text{αN}}\selectlanguage{english}{N+b},\ \ \&m=0\\
0,\ \ \&m>0\\
\end{matrix}\right.\ \)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\alpha\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left\{\begin{matrix}\infty,\ \ \&0\leq m<1\\
\selectlanguage{greek}\frac{\text{αN}}\selectlanguage{english}{b},\ \ \&m=1\\
\end{matrix}\right.\ \)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\selectlanguage{greek}α\selectlanguage{english}P\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left[\left.\ 0,\ \infty\right)\right.\ \)\strut
\end{minipage}\tabularnewline
\begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
Type III\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\frac{\alpha\left(\frac{N}{P^{m}}\right)^{2}}{\left(\frac{N}{P^{m}}\right)^{2}+b}\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left\{\begin{matrix}\frac{\alpha N^{2}}{N^{2}+b},\ \ \&m=0\\
0,\ \ \&m>0\\
\end{matrix}\right.\ \)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\alpha\)\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\emph{0}\strut
\end{minipage} & \begin{minipage}[t]{0.07\columnwidth}\raggedright\strut
\(\left\{\begin{matrix}\infty,\ \ \&0\leq m<0.5\\
\frac{\alpha N^{2}}{b},\ \ \&m=0.5\\
0,\ \ \&0.50\\
\end{matrix}\right.\ \nonumber \\
\end{equation}\begin{equation}
\operatorname{}\left(f\right)=\operatorname{}\left[\frac{\alpha\frac{1}{P^{m}}\ }{\frac{1}{P^{m}}+\frac{b}{N}}\right]=\frac{\alpha\frac{1}{P^{m}}\ }{\frac{1}{P^{m}}}=\alpha\nonumber \\
\end{equation}\begin{equation}
\operatorname{}\left(f\right)=\frac{0}{0+b}=0\nonumber \\
\end{equation}\begin{equation}
\operatorname{}\left(g\right)=\operatorname{}\left[\selectlanguage{greek}\frac{\text{αN}P^{1-m}}\selectlanguage{english}{\left(\frac{N}{P^{m}}\right)+b}\right]=\left\{\begin{matrix}\selectlanguage{greek}\text{αN}\selectlanguage{english}\frac{\infty}{N+b}=\infty,\ \ \&m=0\\
\selectlanguage{greek}\text{αN}\selectlanguage{english}\frac{\infty}{\frac{N}{\infty}+b}=\infty,\ \ \&0